Conventional methods for extracting rare earth metals (REMs) from mined mineral ores are inefficient, expensive, and environmentally damaging. Recent discovery of lanmodulin (LanM), a protein that coordinates REMs with high-affinity and selectivity over competing ions, provides inspiration for new REM refinement methods. Here, we used quantum mechanical (QM) methods to investigate trivalent lanthanide cation (Ln) interactions with coordination systems representing bulk solvent water and protein binding sites. Energy decomposition analysis (EDA) showed differences in the energetic components of Ln interaction with representatives of solvent (water, HO) and protein binding sites (acetate, CHCOO), highlighting the importance of accurate description of electrostatics and polarization in computational modeling of REM interactions with biological and bioinspired molecules. Relative binding free energies were obtained for Ln with coordination complexes originating from binding sites in PDB structures of a lanthanum binding peptide (PDB entry 7CCO) and LanM, with explicit consideration of the first hydration shell waters, according to quasi-chemical theory (QCT). Beyond the first shell, the bulk solvent environment was represented with an implicit continuum model. Ln interactions with (HO) and both binding site models became more favorable, moving down the periodic series. This trend was more pronounced with the protein binding site models than with water, resulting in affinity increasing with periodic number, except for the last REM, Lu, which bound less favorably than the preceding element, Yb. Using the truncated 7CCO binding site model, the magnitude and trend of the experimental Ln relative binding free energies for the whole 7CCO peptide were reproduced. Conversely, the previously reported experimental data for LanM show a preference for the earlier lanthanides; this is likely due to longer-range interactions and cooperative effects, which are not represented by the reduced models. Using the truncated 7CCO binding site model, the magnitude and trend of the experimental Ln relative binding free energies for the whole 7CCO peptide were reproduced. In contrast to the previously reported experimental data for LanM, the peptide preferentially binds the earlier lanthanides. This difference likely arises due to longer-range interactions and cooperative effects not represented by the peptide. Further investigation of Ln interactions with whole proteins using polarizable molecular mechanics models with explicit solvent is warranted to understand the influence of longer-ranged interactions, cooperativity, and bulk solvent. Nevertheless, the present work provides new insights into Ln interactions with biomolecules and presents an effective computational platform for designing specific single-site REM binding peptides more efficiently.