ABSTRACT: The ability to accurately predict the binding site, binding pose, and binding affinity for ligand-RNA binding is important for RNA-targeted drug design. Here, we describe a new computational method, RLDOCK, for predicting the binding site and binding pose for ligand-RNA binding. By developing an energy-based scoring function, we sample exhaustively all of the possible binding sites with flexible ligand conformations for a ligand-RNA pair based on the geometric and energetic scores. The model distinguishes from other approaches in three notable features. First, the model enables exhaustive scanning of all of the possible binding sites, including multiple alternative or coexisting binding sites, for a given ligand-RNA pair. Second, the model is based on a new energy-based scoring function developed here. Third, the model employs a novel multistep screening algorithm to improve computational efficiency. Specifically, first, for each binding site, we use a gird-based energy map to rank the binding sites according to the minimum Lennard-Jones potential energy for the different ligand poses. Second, for a given selected binding site, we predict the ligand pose using a two-step algorithm. In the first step, we quickly identify the probable ligand poses using a coarse-grained simplified energy function. In the second step, for each of the probable ligand poses, we predict the ligand poses using a refined energy function. Tests of the RLDOCK for a set of 230 RNA-ligand-bound structures indicate that RLDOCK can successfully predict ligand poses for 27.8, 58.3, and 69.6% of all of the test cases with the root-mean-square deviation within 1.0, 2.0, and 3.0 Å, respectively, for the top three predicted docking poses. The computational method presented here may enable the development of a new, more comprehensive framework for the prediction of ligand-RNA binding with an ensemble of RNA conformations and the metal-ion effects.