Selective separation of rare earth elements (REEs) from solutions of mixed heavy and light metals by solid adsorbents is an important challenge in the fields of water treatment and metal recovery. The main challenge is water instability of many adsorbents, specifically metal–organic frameworks (MOFs), and their low selectivity. Grafting particular organophosphorus compounds (OPCs) on the MIL-101(Cr) MOF can provide both stability and selectivity. When the tributyl phosphate (TBP), bis(2-ethylhexyl) hydrogen phosphate (D2EHPA or HDEHP) and bis(2,4,4-trimethylpentyl) phosphinic acid (Cyanex®-272) OPCs are grafted and applied to mixed-metal aqueous solutions containing Co²⁺, Ni²⁺, Cu²⁺, Zn²⁺, Nd³⁺, Gd³⁺ and Er³⁺, MIL-101(Cr) offers high selectivity towards the Nd³⁺, Gd³⁺ and Er³⁺ REEs (with stronger affinity towards Er³⁺). However, the underlying chemistry is unknown and the factors leading to the selectivity remain poorly understood. To uncover the key molecular-level factors, we performed state-of-the-art computational simulations using a combination of high-level density functional theory (DFT), semi-empirical calculations, and configurational sampling of the metal ion-MOF binding modes in aqueous solutions. Our simulation study reproduced the available experimental results, in addition to determining the contributing intermolecular interactions, uptake modes and the most significant structural features for improving selectivity towards the REEs. Therefore, our most important result is rationalization of the mechanism of REE separation by OPC-grafted MOFs using quantum mechanical and electrostatic principles. The results provide guidelines for synthesis of OPC-grafted MIL-101(Cr) structures with enhanced selectivity and stability. Moreover, an efficient computational framework is proposed to facilitate comprehensive modeling of similar systems.