In this work, we introduce a new method for sampling multiple binding modes of flexible ligands with rotatable bonds. We use a hybrid method which combines molecular dynamics (MD) and nonequilibrium candidate Monte Carlo (NCMC) for this purpose. Previously it has been shown than MD/NCMC can help to accelerate sampling of multiple binding modes of rigid fragments and side-chain rotameric states of proteins. We validated our MD/NCMC protocol using ligands binding to the protein target c-Jun Nterminal kinase-1 (JNK1). These ligands have different substitutions in a terminal phenyl ring which revolves around a rotatable bond, resulting in two distinct binding modes. Kaus et al. performed extensive alchemical free energy calculations to determine the free energy difference between the two binding modes for these ligands, which can be converted to population estimates. Additionally, we also validated our method using a different ligand binding to the cyclin-dependent kinase 2 (CDK2) protein with experimental population estimates of the two binding modes based on X-ray electron density maps of crystal structures. Using MD/NCMC, we computed population estimates close to those obtained using experimental or other computational techniques. MD/NCMC also required much less computation than standard MD and umbrella sampling, thus suggesting that MD/NCMC is an efficient computational tool for determining multiple binding modes of flexible ligands.