The human telomeric G-quadruplex (G4) is an attractive target for developing anticancer drugs. Natural products protoberberine alkaloids are known to bind human telomeric G4 and inhibit telomerase. Among several structurally similar protoberberine alkaloids, epiberberine (EPI) shows the greatest specificity in recognizing the human telomeric G4 over duplex DNA and other G4s. Recently, NMR study revealed that EPI recognizes specifically the hybrid-2 form human telomeric G4 by inducing large rearrangements in the 50-flanking segment and loop regions to form a highly extensive four-layered binding pocket. Using the NMR structure of the EPI-human telomeric G4 complex, here we perform molecular dynamics free energy calculations to elucidate the ligand selectivity in the recognition of protoberberines by the human telomeric G4. The MM-PB(GB)SA (molecular mechanics-Poisson Boltzmann/Generalized Born) Surface Area) binding free energies calculated using the Amber force fields bsc0 and OL15 correlate well with the NMR titration and binding affinity measurements, with both calculations correctly identifying the EPI as the strongest binder to the hybrid-2 telomeric G4 wtTel26. The results demonstrated that accounting for the conformational flexibility of the DNA-ligand complexes is crucially important for explaining the ligand selectivity of the human telomeric G4. While the MD-simulated (molecular dynamics) structures of the G-quadruplex-alkaloid complexes help rationalize why the EPI-G4 interactions are optimal compared with the other protoberberines, structural deviations from the NMR structure near the binding site are observed in the MD simulations. We have also performed binding free energy calculation using the more rigorous double decoupling method (DDM); however, the results correlate less well with the experimental trend, likely due to the difficulty of adequately sampling the very large conformational reorganization in the G4 induced by the protoberberine binding.