# Extending Libraries of Extremely Localized Molecular Orbitals to Metal Organic Frameworks: A Preliminary Investigation

Erna K Wieduwilt, Giovanni Macetti, Rebecca Scatena, Piero Macchi, Alessandro Genoni

<figure><img src="image_1.png" /></figure>

## To cite this version:

Erna K Wieduwilt, Giovanni Macetti, Rebecca Scatena, Piero Macchi, Alessandro Genoni. Extending Libraries of Extremely Localized Molecular Orbitals to Metal Organic Frameworks: A Preliminary Investigation. Crystals, 2021, 11 (2), pp.207. 10.3390/cryst11020207. hal-03154255

HAL Id: hal-03154255

https://hal.univ-lorraine.fr/hal-03154255

Submitted on 27 Feb 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Extending Libraries of Extremely Localized Molecular Orbitals to Metal Organic Frameworks: A Preliminary Investigation

Erna K. Wieduwilt <fnref n="1" />,*, Giovanni Macetti <fnref n="1" />,*, Rebecca Scatena <fnref n="2" />, Piero Macchi <fnref n="3" />, and Alessandro Genoni <fnref n="1" />,†

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Wieduwilt, E.K.; Macetti, G.; Scatena, R.; Macchi, P.; Genoni, A. Extending Libraries of Extremely Localized Molecular Orbitals to Metal Organic Frameworks: A Preliminary Investigation. *Crystals* **2021**, *11*, 207. https://doi.org/10.3390/cryst11020207

Academic Editor: Linda J. W. Shimon

Received: 20 January 2021

Accepted: 18 February 2021

Published: 20 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Libraries of extremely localized molecular orbitals (ELMOs) have been recently assembled to reconstruct approximate wavefunctions of very large biological systems, such as polypeptides and proteins. In this paper, we investigate for the first time the possibility of using ELMO transferability to also quickly obtain wavefunctions, electron densities, and electrostatic potentials of three-dimensional coordination polymers such as metal organic frameworks (MOFs). To accomplish this task, we propose a protocol that, in addition to exploiting the usual exportability of extremely localized molecular orbitals, also takes advantage of the novel QM/ELMO (quantum mechanics/extremely localized molecular orbital) approach to properly describe the secondary building units of MOFs. As a benchmark test, our technique has been applied to the well-known metal organic framework HKUST-1 ([(Cu₃(BTC)₂)ₙ]ₙ, with BTC=1,3,5-benzenetricarboxylate) to quickly calculate electrostatic potential maps in the small and large cavities inside the network. On the basis of the obtained results, we envisage further improvements and applications of this strategy, which can be also seen as a starting point to perform less computationally expensive quantum mechanical calculations on metal organic frameworks with the goal of investigating transformation phenomena such as chemisorption.

**Keywords:** extremely localized molecular orbitals (ELMOs); ELMO transferability; ELMO libraries; metal organic frameworks (MOFs); electron density; electrostatic potential

## 1. Introduction

The investigation of large systems by means of reliable and computationally affordable techniques has been one of the main topics in theoretical chemistry for a long time. Several strategies have been introduced over the years, most of them consisting of subdividing the (macro)molecule/(macro)system under exam into different subunits that are treated separately, sometimes even at different levels of theory [1,2].

In this area of research, prominent examples are the so-called embedding methods, according to which the chemically crucial region of the system is treated at a higher quantum chemical level, while the remaining part is described through a more approximate strategy. Techniques belonging to this category are the well-known QM/MM (quantum mechanics/molecular mechanics) approaches [3–6] originally introduced by Warshel [7], Levitt [8], and Karplus [9]; the ONIOM (Our own N-layer Integrated molecular Orbital molecular Mechanics) method proposed by Morokuma and collaborators [10,11]; and the more recent and promising fully quantum mechanical embedding approaches. Within the last group, archetypal examples are the frozen density embedding theory (FDET)

<footnote n="1">¹ Université de Lorraine & CNRS, Laboratoire de Physique et Chimie Théoriques (LPCT), UMR CNRS 7019, 1 Boulevard Arago, F-57078 Metz, France; erna-katharina.wieduwilt@univ-lorraine.fr (E.K.W.); giovanni.macetti@univ-lorraine.fr (G.M.)</footnote>

<footnote n="2">² Clarendon Laboratory, Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, UK; rebecca.scatena@physics.ox.ac.uk</footnote>

<footnote n="3">³ Department of Chemistry, Materials and Chemical Engineering, Politecnico di Milano, via Mancinelli 7, 20131 Milano, Italy</footnote>

<footnote>* Correspondence: piero.macchi@polimi.it (P.M.); Alessandro.Genoni@univ-lorraine.fr (A.G.)</footnote>
developed by Wesolowski and coworkers [12–14] or the projection-based embedding (PbE) approach jointly devised by the Manby and Miller research groups [15–24].

Fragmentation methods are also worth mentioning in this context. The idea is to parti-
tion the large system under investigation into smaller and tractable subunits on which fully
quantum mechanical calculations are carried out. The results of these computations are af-
terwards combined to obtain the global wavefunction and/or electron density of the target
macrosystem. Examples are the pioneering Divide and Conquer approach, initially pro-
posed by Yang in the framework of Density Functional Theory (DFT) [25,26] and later
extended by Merz and coworkers to Hartree-Fock (HF) and semiempirical schemes [27–29],
and the molecular tailoring approach devised by Gadre and collaborators at different
quantum chemical levels [30,31]. Special cases of fragmentation strategies are the “frag-
ment interaction methods”. In these techniques, the total energies/electron densities
of the examined macromolecules are the sums of the energies/electron densities of the
different subunits, corrected for the interactions occurring in the dimers, trimers, and
tetramers of the considered subsystems. Notable approaches falling into this category are
the popular fragment molecular orbital (FMO) strategy [32–35], the kernel energy method
(KEM) [36–38], and the molecular fractionation with conjugate caps (MFCC) technique in
its different variants [39–42].

Finally, another important family of computational approaches that allow one to
rapidly obtain reliable (although approximate) wavefunctions and electron densities of
very large molecular systems are those that exploit the concept of transferability, which de-
rives from the fundamental observation that molecules are generally composed of units that
preserve their main features in different environments (see, for instance, the well-known
concept of functional group). On the basis of this observation, different methods have
been introduced since the 1990s, such as the techniques devised by Mezey and coworkers,
who initially proposed the molecular electron density LEGO assembler (MEDLA) strat-
egy [43] and afterwards the adjustable density matrix assembler (ADMA) approach [44,45],
respectively consisting of libraries of fuzzy electron densities and density matrices as-
sociated with molecular fragments. Similar methods include also the transferable atom
equivalent (TAE) technique [46] developed by the Breneman group or other real-space
strategies proposed by Bader [47,48] and Matta [49] within the quantum theory of atoms in
molecules (QTAIM) [50].

In this context, an approach that uses extremely localized molecular orbitals (ELMOs)
has also been recently proposed. ELMOs are molecular orbitals strictly localized on small
molecular subunits (e.g., atoms, bonds, and functional groups) [51–53]. Thanks to their
extreme localization, they can be easily exported from small model molecules to a macro-
system, such as a protein, with the purpose of reconstructing approximate wavefunctions,
electron densities, and electrostatic potentials. The positive results obtained from many
preliminary assessments of the ELMO transferability [53–60] stimulated the compilation of
ELMO libraries [59–61]. They currently contain ELMOs corresponding to all the possible
elementary fragments of the twenty natural amino acids.

These libraries not only allow one to reconstruct the electron density properties of
polypeptides and proteins but are also starting points for advanced structure refinement
techniques. For instance, within quantum crystallography [62–67] the ELMO databanks
have been coupled with the emerging Hirshfeld atom refinement (HAR) approach [68–73],
giving rise to the HAR-ELMO method [74] which allowed successful refinements of crystal
structures of relatively large polypeptides and of the small protein crambin. Furthermore,
the database has been also used to improve the non-covalent interaction (NCI) [75,76] and
independent gradient model (IGM) [77–80] strategies for the detection of non-covalent
interactions in large molecular systems [81,82]. Finally, the ELMO libraries have been the
basis to develop the novel QM/ELMO (quantum mechanics/extremely localized molecular
orbital) approach [83–86]—namely, a new fully quantum mechanical embedding strategy in
which the chemically relevant region of the system is treated at a high quantum mechanical
level, while the rest is described by means of transferred and frozen extremely localized
molecular orbitals. Interestingly, the QM/ELMO technique has also been recently exploited
in the framework of Hirshfeld atom refinement to properly account for crystal environment
effects, thus improving the performances of traditional HAR [87].

The current ELMO databanks were initially conceived to investigate extended molecules
of biological interest [61]. However, given the flexibility of the *ELMOdb* software [61]
interfaced with the constructed ELMO libraries (see Section 2.1), the ELMO transferability
holds for almost any kind of chemical system. For example, HAR-ELMO refinements of
crystal structures of organometallic complexes were recently reported [74].

Having all of this in mind, in this paper we present the first preliminary study on
the ELMO transferability to investigate metal organic frameworks (MOFs) [88–91], which,
differently from chain-like polymers such as polypeptides and proteins, are extended
three-dimensional coordination polymers and need a slightly different treatment (see
Section 2.2). ELMO databanks would also enable instantaneous calculations of quantum
mechanically rigorous electron densities and electrostatic potentials for MOFs, whereas QM
calculations with periodic boundary conditions are extremely demanding. The ELMO
strategy would allow the mapping of the most likely binding sites available for extra-
framework guest molecules trapped and adsorbed inside the MOFs cavities. Moreover,
a database of ELMOs for MOFs will also be the starting point to perform QM/ELMO
calculations for metal organic frameworks, which could be of invaluable help to investigate
chemisorption/physisorption phenomena.

The paper is organized as follows. In the next section, in addition to briefly reviewing
the concept of ELMOs and the recently constructed ELMO libraries, we will provide the
computational details of the present investigation, mainly focusing on the protocol that
we have fine-tuned to determine and transfer extremely localized molecular orbitals for
MOFs. In Section 3, the obtained results will be shown and discussed. Finally, in the last
part of the paper (Section 4) conclusions will be drawn and possible future perspectives
will be presented.

**2. Theoretical and Computational Details**

2.1. *Extremely Localized Molecular Orbitals and ELMO Libraries*

Localized molecular orbitals have always been used by theoretical chemists as useful
tools to connect quantum chemistry calculations with traditional chemical concepts that
derive from the localized picture of electronic structure based on Lewis diagrams.

Since the canonical molecular orbitals obtained through Hartree–Fock computations
are completely delocalized on all atoms of the examined systems, several localization
methods have been proposed. Some of them consist in unitary transformations of the
canonical Hartree–Fock molecular orbitals and aim at the minimization/maximization
of physically sound functionals [92–97]. The molecular orbitals obtained through these
techniques are mainly localized on small molecular fragments, but they present tiny
orthogonalization tails that prevent a straightforward transferability. If one really wants to
have strictly localized molecular orbitals easily transferable from molecule to molecule, it is
necessary to resort to localization approaches that constrain a priori (i.e., before starting the
calculations) the molecular orbital expansions over local sets of basis functions [51,98–106].

The extremely localized molecular orbitals considered in the present paper adhere to
the technique proposed by Stoll and coworkers [51]. Following this strategy, the system is
partitioned into small molecular fragments, such as atoms or bonds, according to the Lewis
structure of the molecule. To be more precise, the examined system is subdivided into
groups of electrons localized on atoms (core and lone pair electrons) and groups of electrons
localized on bonds (bond electrons). This subdivision automatically leads to the definition
of local basis sets, which consist of atomic orbitals centered on the atoms constituting the
fragments and are used to expand the ELMOs that describe the electrons of the different
subunits. These extremely localized molecular orbitals are afterwards determined by
variationally minimizing the energy of the single Slater determinant constructed with
them—namely, by variationally minimizing the energy associated with the global ELMO
wavefunction (see Supporting Information for more theoretical details about the Stoll method). Proceeding in this way, it is possible to determine molecular orbitals easily associable with elementary molecular subunits and really transferable from one molecule to another. As an example, in Figure 1 we have depicted the ELMOs that one can obtain through a simple ELMO calculation on acetamide by adopting a localization scheme corresponding to the Lewis diagram of the molecule—namely, ELMOs directly associable with the C-H, C-N, C=O, and N-H bonds and to the lone pairs of nitrogen and oxygen atoms (ELMOs corresponding to core electrons are omitted in Figure 1). As an example, the carbonyl C=O fragment corresponds to two ELMOs (see Figure 1G,H) that overall describe four bond electrons. At the same time, the oxygen atom fragment is associated with three ELMOs, two of which (see Figure 1I,J) describe the four lone pair electrons, while the remaining one (not depicted in Figure 1) describes the two core electrons.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Examples of extremely localized molecular orbitals computed for acetamide (see panel (A)) using a localization scheme corresponding to the Lewis diagram of the molecule: (B) ELMO describing the lone pair of the nitrogen atom, (C) ELMO associated with one of the two N-H bonds, (D) ELMO corresponding to the N-C bond, (E) ELMO describing the C-C bond, (F) ELMO associated with one of the three C-H bonds, (G,H) ELMOs corresponding to the C=O double bond ($\sigma$ and $\pi$ orbitals, respectively), (I,J) ELMOs for the lone pairs of the oxygen atom. All the orbitals were computed with the cc-pVDZ basis set and were plotted considering a 0.15 a.u. isosurface.</figcaption></figure>

As already mentioned in the Introduction, due to the absolute localization of the extremely localized molecular orbitals, the Stoll technique was exploited to construct ELMO libraries that cover all the possible elementary subunits of the twenty natural amino acids in all their possible protonation states and forms (i.e., N-terminal, C-terminal, and non-terminal) [61]. The current databases contain (i) ELMOs localized on one-atom fragments for the description of core and lone pair electrons, (ii) ELMOs localized on two-atom subunits for the description of ordinary bonds, and (iii) ELMOs localized on three-atom fragments for the description of a particular bonding situation in which the electronic structure is more delocalized (i.e., peptide bonds, carboxylic/carboxylate groups, aromatic rings). The databases are available for five standard basis sets of quantum chemistry: 6-31G, 6-31G(d,p), 6-311G, 6-311G(d,p), and cc-pVDZ.
The extremely localized molecular orbitals stored in the ELMO libraries can be transferred to the examined target systems by exploiting the associated *ELMOdb* program [61], which implements the Philipp and Friesner technique for the rotation/transfer of strictly localized bond orbitals (see again Supporting Information for more theoretical details) [59,107]. It is worth noting that *ELMOdb* has an open interface towards ELMOs that are not currently included in the databases, but that are user-calculated ad hoc (with any basis set of choice) for the target system under exam. This is in fact the strategy adopted in this study to read and transfer the extremely localized molecular orbitals of each subunit of HKUST-1 (see Section 2.2). After transferring the databank or user-defined ELMOs, *ELMOdb* computes the one-particle density matrix of the whole system and eventually provides both a Gaussian formatted checkpoint file and a wavefunction file (wfx format), useful for subsequent analyses and calculations, such as topological and/or NCI analyses or electron density/electrostatic potential mappings.

## 2.2. ELMO-Protocol for MOFs

The metal organic framework chosen for our preliminary study was {[Cu₃(BTC)₂]n (where BTC stands for 1,3,5-benzenetricarboxylate), which is also commonly known as HKUST-1 (see Figure 2) [108,109]. Since the current ELMO libraries cover only the basic fragments of the twenty natural amino acids, it was primarily necessary to obtain tailor-made ELMOs computed on suitable model molecules of the MOF subunits.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Crystal structure of {[Cu₃(BTC)₂]n (HKUST-1): (left panel) unit-cell content with large and small cavities highlighted by ochre and blue spheres, respectively; (right panel) paddlewheel model system for the QM/ELMO calculations.</figcaption></figure>

To accomplish this task, we have initially determined the ELMOs for the elementary fragments of the 1,3,5-benzenetricarboxylate linker using the 6-311G(2d,2p) basis set on the geometry of the BTC molecule optimized at B3LYP/6-311G(2d,2p) level.

Afterwards, we took into account the secondary building units (SBUs), which consist of two Cu(II) atoms. The secondary building units represent an increased complexity in the hierarchy of ELMO libraries and, in general, the ELMO description is not the most suitable one for them. This was especially the case for HKUST-1, where the two Cu(II) atoms of the secondary building units can be coupled in three different ways: (i) diamagnetically, (ii) anti-ferromagnetically, and (iii) ferromagnetically.

Therefore, to account for the previous three possibilities and to obtain suitable molecular orbitals for the description of the SBUs, we decided to perform QM/ELMO calculations on paddlewheel model systems consisting of the Cu-Cu atom pair surrounded by four BTC linkers (see the right panel of Figure 2). In all cases, the QM region coincided with the Cu-
Cu atom pair and was treated at DFT-B3LYP level, while the ELMO subsystem consisted of the four surrounding 1,3,5-benzenetricarboxylate units and was described with the frozen ELMOs pre-calculated for the BTC linker (see above). The geometries of the model systems for the QM/ELMO computations were extracted from the experimental crystal structures of HKUST-1 containing different amounts of water molecules [109]. For the model system of each crystal structure, three different types of B3LYP/ELMO calculations were performed, with a QM region in: (i) singlet-restricted closed-shell state (diamagnetic case), (ii) singlet-broken symmetry state (anti-ferromagnetic case), and (iii) triplet-unrestricted state (ferromagnetic case). The standard triple-zeta basis set 6-311G(2d,2p) was used for both the quantum mechanical and the ELMO subsystem.

For all the performed QM/ELMO computations, the molecular orbitals of the QM region were mainly localized on the subunit constituted by the Cu pair, thus allowing us to extract localized molecular orbitals for the secondary building units that also properly take into account the spin state of the system. The localized molecular orbitals corresponding to the Cu-Cu atom pair were then transferred to the other symmetry-related secondary building units of the examined crystal structure by simply exploiting the crystal symmetry operations.

The results that we obtained indicate a suitable strategy/protocol for the future extension of the ELMO libraries to MOFs:

1. Construct a library of ELMOs describing the elementary fragments of the most common linkers employed for MOFs design.

2. For each MOF crystal structure, perform ad hoc QM/ELMO calculations on model systems consisting of a SBU (at QM level) and the connected linkers (using ELMOs computed at step 1).

3. The ELMOs for the linkers (see point 1) and the localized molecular orbitals for the SBU (obtained from point 2) will be finally transferred to the symmetry-related positions of the crystal structure using the *ELMOdb* program [61], thus quickly obtaining an approximate wavefunction/electron density for the periodic system.

In this work, all the ELMO calculations were performed using a modified version of the *GAMESS-UK* quantum chemistry package [110], where the Stoll equations are implemented [52]. The QM/ELMO computations were carried out by exploiting a modified version of the *Gaussian09* quantum chemical software [111].

**3. Results and Discussion**

As the first step of this investigation, we focused only on the model for the QM/ELMO calculations necessary for the determination of the molecular orbitals corresponding to the secondary building units (see Section 2.2 and the right panel of Figure 2). The geometry was taken from the HKUST-1 experimental crystal structure with a water chemisorption degree equal to 14.3(6)% [109], the experimental geometry which is very close to an exactly empty structure.

We performed standard DFT-B3LYP and B3LYP/ELMO calculations with basis set 6-311G(2d,2p) for each of the three possible spin states of the model system: restricted closed-shell singlet (diamagnetic case), broken-symmetry singlet (anti-ferromagnetic case), and open-shell triplet (ferromagnetic case). Since the primary goal of the present study is to evaluate the reliability of the ELMO and QM/ELMO approximations in providing electron densities and electrostatic potentials, we compared the NCI plots (see Figures 3 and 4) and the electrostatic potential maps (see Figure 5) resulting from all the above-mentioned computations, using the standard quantum mechanical results as benchmarks.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** NCI plots at different levels of theory with basis set 6-311G(2d,2p) for the model system extracted from the HKUST-1 experimental crystal structure with water chemisorption degree of 14.3(6)%. All the reduced-density gradient isosurfaces correspond to the 0.4 a.u. isovalue and are colored according to the BGR (blue-green-red) scheme over the range -5.0 a.u < sign($\lambda_2$)$\rho$ < 5.0 a.u.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Zoom on the Cu-Cu subunit of the NCI plots depicted in Figure 3. All the reduced-density gradient isosurfaces correspond to the 0.4 a.u. isovalue and are colored according to the BGR (blue-green-red) scheme over the range -5.0 a.u < sign($\lambda_2$)$\rho$ < 5.0 a.u.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Electrostatic potential maps at different levels of theory with basis set 6-311G(2d,2p) for the model system extracted from the HKUST-1 experimental crystal structure with a water chemisorption degree of 14.3(6)%. The electrostatic potentials were plotted on the corresponding 0.05 a.u. electron density isosurfaces according to the RWB (red-white-blue) scheme over the range [−0.1 a.u., 0.1 a.u.].</figcaption></figure>

The NCI plots are practically identical in the linker regions, regardless of the chosen level of theory (B3LYP or B3LYP/ELMO) and spin state (ferromagnetic, anti-ferromagnetic or diamagnetic); see Figure 3. On the contrary, in the quantum mechanical subsystem of the B3LYP/ELMO calculations (i.e., the Cu-Cu subunit), a strong attractive interaction occurs between the copper atoms in all cases, although the blue isosurface is larger in the plots associated with the standard B3LYP calculations (see Figure 4). The most important differences are in the regions of the Cu-linker coordination (see again Figure 4). In fact, the B3LYP/ELMO calculations seem to introduce some artifacts. This is not surprising, because those regions are the frontiers between the fully quantum mechanical subunit and the ELMO subsystem in the B3LYP/ELMO computations.

For a more quantitative comparison of the electron densities, the values of different topological properties were also examined (see Tables 1 and 2). In particular, we analyzed the Cu–Cu and Cu–O interactions. A bond critical point (bcp) was always detected for both the examined interactions, independently of the chosen level of theory and of the spin state considered in the calculations.

Let us first consider the Cu-Cu case. At variance from many other M-M bonds linked through genuine covalent interactions [112], the Cu-Cu bond in this sub-unit is anomalous because the bond direction coincides with the overlap between fully occupied 3d<sub>z<sup>2</sup></sub> orbitals, which are also responsible for the typical Jahn-Teller distortion [113] along this direction for ligand molecules occupying the apical site (such as water molecules). Nevertheless, the bond presents similar topological indicators as typical metal-metal bonds—namely, a small electron density at the critical point associated with a local charge depletion (∇<sup>2</sup>ρ<sub>bcp</sub> > 0), but with a dominance of the potential energy density over the kinetic energy density (|V<sub>bcp</sub>|/G<sub>bcp</sub> > 1), hence with a negative total energy density (H<sub>bcp</sub> = −K<sub>bcp</sub> = V<sub>bcp</sub> + G<sub>bcp</sub> < 0). This reflects the dominant contribution of the 4s orbitals at the bcp, in agreement with their diffuse nature and despite the nominal oxidation state of the copper atoms. For the same reason, the topological indices are substantially invariant on changing the spin states, because this mainly affects the d-orbitals of the metals.
On the other hand, the delocalization index DI (i.e., the number of electron pairs shared between the two atoms) is severely affected either by the type of calculation and by the forced spin state. The QM/ELMO approximation overestimates the role of the metal-metal sharing in the diamagnetic state, whereas it slightly underestimates it in the ferro- and anti-ferromagnetic states. The overestimation in the B3LYP/ELMO diamagnetic state can be explained with the fact that the electrons are practically constrained to be localized in the QM and ELMO regions, while in the standard QM calculations the electrons are free to spread over the whole molecule. This can be also seen in Table S1 of the Supporting Information, where we reported how the delocalization index changes as a function of the size of the QM region. If the four carboxylate groups are included in the QM subunit, the DI for the Cu-Cu interaction approaches the result of the standard QM computation on the paddlewheel model system.

<table><caption>**Table 1.** Properties at the bond critical point between the Cu atoms and between the Cu and O atoms, as obtained from the calculations performed at different levels of theory with basis set 6-311G(2d,2p) on the model system extracted from the HKUST-1 experimental crystal structure with a water chemisorption degree of 14.3(6)%. (a)</caption><thead><tr><th>Interaction</th><th>Spin State</th><th>Calculation</th><th>ρ<sub>bcp</sub></th><th>∇<sup>2</sup>ρ<sub>bcp</sub></th><th>|V<sub>bcp</sub>|/G<sub>bcp</sub></th><th>-K<sub>bcp</sub></th><th>DI</th></tr></thead><tbody><tr><td rowspan="4">Cu-Cu</td><td rowspan="2">Diamagnetic</td><td>B3LYP/ELMO</td><td>0.029</td><td>0.075</td><td>1.344</td><td>-0.0098</td><td>1.031</td></tr><tr><td>Standard B3LYP</td><td>0.036</td><td>0.069</td><td>1.423</td><td>-0.0127</td><td>0.426</td></tr><tr><td rowspan="2">Anti-Ferromagnetic</td><td>B3LYP/ELMO</td><td>0.029</td><td>0.074</td><td>1.348</td><td>-0.0099</td><td>0.100</td></tr><tr><td>Standard B3LYP</td><td>0.035</td><td>0.070</td><td>1.410</td><td>-0.0122</td><td>0.146</td></tr><tr><td rowspan="4">Cu-O</td><td rowspan="2">Ferromagnetic</td><td>B3LYP/ELMO</td><td>0.029</td><td>0.074</td><td>1.348</td><td>-0.0099</td><td>0.100</td></tr><tr><td>Standard B3LYP</td><td>0.035</td><td>0.070</td><td>1.409</td><td>-0.0121</td><td>0.144</td></tr><tr><td rowspan="2">Diamagnetic</td><td>B3LYP/ELMO</td><td>0.066</td><td>0.614</td><td>0.956</td><td>0.0065</td><td>0.227</td></tr><tr><td>Standard B3LYP</td><td>0.095</td><td>0.462</td><td>1.140</td><td>-0.0189</td><td>0.429</td></tr><tr><td rowspan="4">Anti-Ferromagnetic</td><td rowspan="2">Ferromagnetic</td><td>B3LYP/ELMO</td><td>0.065</td><td>0.626</td><td>0.943</td><td>0.0084</td><td>0.233</td></tr><tr><td>Standard B3LYP</td><td>0.093</td><td>0.477</td><td>1.127</td><td>-0.0174</td><td>0.466</td></tr><tr><td rowspan="2">Ferromagnetic</td><td>B3LYP/ELMO</td><td>0.065</td><td>0.626</td><td>0.943</td><td>0.0084</td><td>0.233</td></tr><tr><td>Standard B3LYP</td><td>0.093</td><td>0.478</td><td>1.127</td><td>-0.0173</td><td>0.466</td></tr></tbody></table>

<table><caption>**Table 2.** Bader charges and volumes (0.001 e/bohr³ isosurface) for the Cu and O atoms, as obtained from calculations performed at different levels of theory with basis set 6-311G(2d,2p) on the model system extracted from the HKUST-1 experimental crystal structure with a water chemisorption degree of 14.3(6)%. (a)</caption><thead><tr><th>Spin State</th><th>Calculation</th><th>q<sub>Cu1</sub></th><th>q<sub>Cu2</sub></th><th>q<sub>O</sub></th><th>V<sub>Cu1</sub></th><th>V<sub>Cu2</sub></th><th>V<sub>O</sub></th></tr></thead><tbody><tr><td rowspan="2">Diamagnetic</td><td>B3LYP/ELMO</td><td>1.79</td><td>1.79</td><td>-1.42</td><td>61.50</td><td>61.50</td><td>109.45</td></tr><tr><td>Standard B3LYP</td><td>1.10</td><td>1.10</td><td>-1.13</td><td>75.03</td><td>75.03</td><td>104.62</td></tr><tr><td rowspan="2">Anti-Ferromagnetic</td><td>B3LYP/ELMO</td><td>1.79</td><td>1.79</td><td>-1.42</td><td>61.52</td><td>61.52</td><td>109.46</td></tr><tr><td>Standard B3LYP</td><td>1.20</td><td>1.20</td><td>-1.16</td><td>73.41</td><td>73.41</td><td>104.95</td></tr><tr><td rowspan="2">Ferromagnetic</td><td>B3LYP/ELMO</td><td>1.79</td><td>1.79</td><td>-1.42</td><td>61.52</td><td>61.52</td><td>109.46</td></tr><tr><td>Standard B3LYP</td><td>1.21</td><td>1.21</td><td>-1.16</td><td>73.35</td><td>73.35</td><td>104.97</td></tr></tbody></table>

The partial lack of charge transfer in the QM/ELMO computations also affects the Bader charges and volumes reported in Table 2. In fact, the B3LYP/ELMO computations for all spin states returned more positive charges and smaller volumes for the copper atoms, and, consequently, more negative charges and larger volumes for the oxygen atoms of the surrounding carboxylate groups.

<footnote>(a) $ρ_{bcp}$ and $∇^2ρ_{bcp}$ are, respectively, the electron density (e/bohr³) and the Laplacian of the electron density (e/bohr⁵) at the bond critical point; $|V_{bcp}|/G_{bcp}$ is the ratio between the potential and the kinetic energy density at the bond critical point; $-K_{bcp}$ is the total bond energy density at the bond critical point (hartree/bohr³); and DI is the delocalization index (electron pairs shared between two atoms).</footnote>

<footnote>(a) Charges in electrons (e) and volumes in bohr³.</footnote>
Concerning the Cu-O interaction, the results of the B3LYP/ELMO and standard
calculations differ, but they are independent from the spin states. At the bond critical
point, the B3LYP/ELMO calculations always returned smaller electron density values and
a more positive Laplacian. Moreover, the ratio $|V_{bcp}|/G_{bcp}$ is always slightly lower than
1.0 for the QM/ELMO computations, while it is slightly greater than 1.0 for the standard
quantum mechanical calculations. Deviations are also observed for the DI index, which is
always larger in the standard B3LYP cases. However, also in this case the QM/ELMO
results converge towards the fully quantum mechanical ones as the size of the QM region
becomes larger (see again Table S1 in the Supporting Information). This can be interpreted
again with the fact that in standard quantum mechanical computations the electrons
are free to delocalize all over the molecule. As already seen for the Cu-Cu interaction,
this discrepancy in the values of delocalization index is consistent with the trend observed
for the Bader charges and volumes and, particularly for this case, with the more negative
charges observed for the oxygen atoms of the carboxylate groups when the B3LYP/ELMO
calculations are taken into account.

The deviations observed for the topological properties of the Cu–O interactions and
the artifacts in the NCI plots are certainly due to the imposed frontier between the SBU and
the surrounding linkers in B3LYP/ELMO calculations. However, we want to point out that,
in this preliminary investigation, we only chose the smallest and simplest QM region for all
the QM/ELMO computations, with the goal of obtaining transferable molecular orbitals for
the secondary building units. As we have seen from the results obtained with a larger QM
region (see Table S1 in the Supporting Information), the definition of a more extended and
suitable quantum mechanical subsystem improves the quality of the results and reduces
the deviations from the outcomes of the corresponding standard QM computations.

In Figure 5, we reported the maps of the electrostatic potentials at different levels of
theory and spin-states for the model system of HKUST-1. All the maps are very similar,
with very small discrepancies only in the frontier regions, in agreement with the results of
the NCI and topological analyses shown and discussed in the previous paragraphs.

One could be tempted to use the energies computed with the QM/ELMO method for
the paddlewheel and predict also the coupling constants. However, as already pointed
out above, the QM region was very small, with the sole purpose of determining molecular
orbitals as much localized as possible on the Cu-Cu secondary building unit for subsequent
reconstructions of the electron density/electrostatic potential of larger portions of HKUST-1.
For reliable energetic values, one should enlarge the QM subsystem and use a higher
quantum mechanical level of theory because the energy differences associated with different
spin configurations, which are needed for the calculation of magnetic exchange coupling
constants, are small.

On the basis of the results obtained from the preliminary test-calculations, by follow-
ing the procedure described in Section 2.2 we afterwards exploited the transferability of the
ELMOs for the linkers and the exportability of the localized molecular orbitals for the sec-
ondary building units to reconstruct the wavefunctions, electron densities, and electrostatic
potentials of large portions of the metal organic framework HKUST-1. In particular, as
target systems we considered both a portion of HKUST-1 that forms one of its small cavities,
and one portion of HKUST-1 that forms one of its large pores (see again left panel of Figure
2 for the graphical representation of the small and large cavities of the considered MOF).

For these target systems, we studied how the presence of water molecules affects
the electrostatic potential inside the pores, which are relevant sites to study the water
sorption ability of HKUST-1. For the small and large cavity portions of the investigated
MOF, we considered the experimental geometries of the filled HKUST-1 [109]. Then,
we compared the electrostatic potential maps calculated with and without the water
molecules keeping the MOF geometry frozen. For this calculation we transferred the
localized molecular orbitals for the secondary building units (anti-ferromagnetic case),
the ELMOs for the linkers and, when necessary, the ELMOs for the water molecules. In the
large cavity (see Figure 6), we see that in the absence of water the potential is highly
symmetric, and a minimum (i.e., most negative) site occurs in the center of the cavity, whereas strongly positive regions coincide with the apical sites of the paddlewheel where in fact water molecules bind the framework. Upon the water saturation of the Cu sites, the high symmetry is broken. Indeed, in the experimental structure a disorder is supposed to take place and the ordered simulation here simply reflects one of the many possible conformations. The potential is radically modified and the negative area is reduced in size. This situation reflects the field experienced by additional water molecules hosted in the MOF upon the saturation of all Cu sites. For the small cavities (see Figure 7), the potential is stronger, but it is less perturbed by the water coordination as no Cu unsaturated site points toward the cavity center, at variance from the large cavity.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Electrostatic potential maps for the large cavity of HKUST-1 (A) without water molecules and (B) in the presence of water molecules, plotted both on the 0.05 a.u. isosurface of the electron density and as isocontour lines in the middle of the pore. The color scale for the maps is given on the right-hand side of the figure.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Electrostatic potential maps for the small cavity of HKUST-1 (A) without water molecules and (B) in the presence of water molecules, plotted both on the 0.05 a.u. isosurface of the electron density and as isocontour lines in the middle of the pore. The color scale for the maps is given on the right-hand side of the figure.</figcaption></figure>

**4. Conclusions and Perspectives**

In this paper, we have investigated for the first time the possibility of exploiting the transferability of extremely localized molecular orbitals to quickly and reliably reconstruct wavefunctions, electron densities, and electrostatic potentials of large portions of metal organic frameworks. To accomplish this task, we have proposed a protocol that, in addition to exploiting the ELMOs exportability to describe the MOFs linkers, also takes advantage of the recent QM/ELMO embedding approach to compute transferable molecular orbitals suitable for the MOFs secondary building units. After preliminary benchmark calcula-
tions, the new proposed strategy has been profitably applied to the well-known MOF HKUST-1, for which electrostatic potential maps in its small and large cavities have been easily reconstructed.

In future works, we plan to generalize the protocol proposed in this paper. The first step will consist in extending the current ELMO libraries by including extremely localized molecular orbitals to describe the elementary fragments of most linkers occurring in three-dimensional coordination polymers. Afterwards, it will be also necessary to fine tune the QM/ELMO-based procedure to obtain the molecular orbitals for secondary building units.

The availability of ELMO libraries for metal organic frameworks will allow not only a rapid identification of suitable binding sites of guest molecules but also more adequate QM/ELMO calculations extending the QM region to one or more paddlewheel structures, while using frozen ELMOs for the rest of the MOF. This will be important for investigating chemical transformations involving MOFs (e.g., chemisorption processes or post-synthetic modifications) and also to obtain reliable energies for specific spin states, which are of fundamental importance to determine magnetic and spectroscopic properties.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2073-4352/11/2/207/s1: Theoretical details about the Stoll method for the determination of extremely localized molecular orbitals; theoretical details on the Philipp and Friesner strategy for the ELMOs rotation. Table S1 showing the values of the electron density at the bond critical point and of the delocalization index for the Cu-Cu and Cu-O interactions at different spin states, as obtained from standard B3LYP and B3LYP/ELMO calculations with different sizes of the QM region on the paddlewheel model structure.

**Author Contributions:** Conceptualization, E.K.W., G.M., R.S., P.M., and A.G.; methodology, E.K.W., G.M., P.M., and A.G.; software, E.K.W., G.M., and A.G.; validation, E.K.W.; formal analysis, E.K.W.; investigation, E.K.W.; resources, R.S., P.M., and A.G.; data curation, G.M., E.K.W., and A.G.; writing—original draft preparation, A.G. and P.M.; writing—review and editing, E.K.W., G.M., R.S., P.M., and A.G.; visualization, E.K.W., R.S.; supervision, A.G. and P.M.; project administration, A.G. and P.M.; funding acquisition, A.G. and P.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** A.G., E.K.W. and G.M. gratefully acknowledge the French Research Agency (ANR) for financial support of this research through the Young Investigator Project QuMacroRef (Grant No. ANR-17-CE29-0005-01). P.M. and R.S. thank the Swiss National Science Foundation (project nr. 160157).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding authors.

**Acknowledgments:** The High-Performance Computing Center EXPLOR of the University of Lorraine is thanked for providing computing time through the projects 2019CPMXX0966, 2019CPMXX0886 and 2019CPMXX1332. Fabien Pascale is also acknowledged for the set-up and maintenance of our local cluster, which was used to perform most of the calculations reported in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Gordon, M.S.; Slipchenko, L.V. Introduction: Calculations on Large Systems. *Chem. Rev.* **2015**, *115*, 5605–5606. [CrossRef]

2. Jones, L.O.; Mosquera, M.A.; Schatz, G.C.; Ratner, M.A. Embedding Methods for Quantum Chemistry: Applications from Materials to Life Sciences. *J. Am. Chem. Soc.* **2020**, *142*, 3281–3295. [CrossRef]

3. Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Car-bonium ion in the Reaction of Lysozyme. *J. Mol. Biol.* **1976**, *103*, 227–249. [CrossRef]

4. Field, M.J.; Bash, P.A.; Karplus, M. A Combined Quantum Mechanical and Molecular Mechanical Potential for Molecular Dynamics Simulations. *J. Comput. Chem.* **1990**, *11*, 700–733. [CrossRef]

5. Gao, J. Methods and Applications of Combined Quantum Mechanical and Molecular Mechanical Potentials. In *Reviews in Computational Chemistry*; Lipkowitz, K.B., Boyd, D.B., Eds.; VCH Publishers, Inc.: Weinheim, Germany, 1996; Volume 7, pp. 119–186. [CrossRef]
6. Senn, H.M.; Thiel, W. QM/MM Methods for Biomolecular Systems. *Angew. Chem. Int. Ed.* **2009**, *48*, 1198–1229. [CrossRef]

7. Warshel, A. Multiscale Modeling of Biological Functions: From Enzymes to Molecular Machines (Nobel Lecture). *Angew. Chem. Int. Ed.* **2014**, *53*, 10020–10031. [CrossRef]

8. Levitt, M. Birth and Future of Multiscale Modeling for Macromolecular Systems (Nobel Lecture). *Angew. Chem. Int. Ed.* **2014**, *53*, 10006–10018. [CrossRef]

9. Karplus, M. Development of Multiscale Models for Complex Chemical Systems: From H+H₂ to Biomolecules (Nobel Lecture). *Angew. Chem. Int. Ed.* **2014**, *53*, 9992–10005. [CrossRef] [PubMed]

10. Svensson, M.; Humbel, S.; Froese, R.D.J.; Matsubara, T.; Sieber, S.; Morokuma, K. ONIOM: A Multilayered Integrated MO + MM Method for Geometry Optimizations and Single Point Energy Predictions. A Test for Diels–Alder Reactions and Pt(P(t-Bu)₃)₂+ H₂ Oxidative Addition. *J. Phys. Chem.* **1996**, *100*, 19357–19363. [CrossRef]

11. Chung, L.W.; Sameera, W.M.C.; Ramozzi, R.; Page, A.J.; Hatanaka, M.; Petrova, G.P.; Harris, T.V.; Li, X.; Ke, Z.; Liu, F.; et al. The ONIOM Method and Its Applications. *Chem. Rev.* **2015**, *115*, 5678–5796. [CrossRef] [PubMed]

12. Wesolowski, T.A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. *J. Phys. Chem.* **1993**, *97*, 8050–8053. [CrossRef]

13. Wesolowski, T.A. Embedding a multideterminantal wave function in an orbital-free environment. *Phys. Rev. A* **2008**, *77*, 012504. [CrossRef]

14. Wesolowski, T.A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. *Chem. Rev.* **2015**, *115*, 5891–5928. [CrossRef]

15. Manby, F.R.; Stella, M.; Goodpaster, J.D.; Miller, T.F., III. A Simple, Exact Density-Functional-Theory Embedding Scheme. *J. Chem. Theory Comput.* **2012**, *8*, 2564–2568. [CrossRef]

16. Barnes, T.A.; Goodpaster, J.D.; Manby, F.R.; Miller, T.F., III. Accurate basis set truncation for wavefunction embedding. *J. Chem. Phys.* **2013**, *139*, 24103. [CrossRef]

17. Goodpaster, J.D.; Barnes, T.A.; Manby, F.R.; Miller, T.F., III. Accurate and systematically improvable density functional theory embedding for correlated wave functions. *J. Chem. Phys.* **2014**, *140*, 18A507. [CrossRef]

18. Bennie, S.J.; Stella, M.; Miller, T.F., III; Manby, F.R. Accelerating wavefunction in density-functional-theory embedding by truncating the active basis set. *J. Chem. Phys.* **2015**, *143*, 024105. [CrossRef]

19. Pennifold, R.C.R.; Bennie, S.J.; Miller, T.F., III; Manby, F.R. Correcting density-driven errors in projection-based embedding. *J. Chem. Phys.* **2017**, *146*, 084113. [CrossRef]

20. Welborn, M.; Manby, F.R.; Miller, T.F., IIII. Even-handed subsystem selection in projection-based embedding. *J. Chem. Phys.* **2018**, *149*, 144101. [CrossRef] [PubMed]

21. Chulhai, D.V.; Goodpaster, J.D. Improved Accuracy and Efficiency in Quantum Embedding through Absolute Localization. *J. Chem. Theory Comput.* **2017**, *13*, 1503–1508. [CrossRef]

22. Chulhai, D.V.; Goodpaster, J.D. Projection-Based Correlated Wave Function in Density Functional Theory Embedding for Periodic Systems. *J. Chem. Theory Comput.* **2018**, *14*, 1928–1942. [CrossRef]

23. Bennie, S.J.; van der Kamp, M.W.; Pennifold, R.C.R.; Stella, M.; Manby, F.R.; Mulholland, A.J. A Projector-Embedding Approach for Multiscale Coupled-Cluster Calculations Applied to Citrate Synthase. *J. Chem. Theory Comput.* **2016**, *12*, 2689–2697. [CrossRef] [PubMed]

24. Lee, S.J.R.; Welborn, M.; Manby, F.R.; Miller, T.F., III. Projection-Based Wavefunction-in-DFT Embedding. *Acc. Chem. Res.* **2019**, *52*, 1359–1368. [CrossRef] [PubMed]

25. Yang, W. Direct calculation of electron density in density-functional theory. *Phys. Rev. Lett.* **1991**, *66*, 1438–1441. [CrossRef] [PubMed]

26. Yang, W. Direct calculation of electron density in density-functional theory: Implementation for benzene and a tetrapeptide. *Phys. Rev. A* **1991**, *44*, 7823–7826. [CrossRef]

27. Dixon, S.L.; Merz, K.M., Jr. Semiempirical molecular orbital calculations with linear system size scaling. *J. Chem. Phys.* **1996**, *104*, 6643–6649. [CrossRef]

28. Dixon, S.L.; Merz, K.M., Jr. Fast, accurate semiempirical molecular orbital calculations for macromolecules. *J. Chem. Phys.* **1997**, *107*, 879–893. [CrossRef]

29. He, X.; Merz, K.M., Jr. Divide and Conquer Hartree–Fock Calculations on Proteins. *J. Chem. Theory Comput.* **2010**, *6*, 405–411. [CrossRef] [PubMed]

30. Gadre, S.R.; Shirsat, R.N.; Limaye, A.C. Molecular Tailoring Approach for Simulation of Electrostatic Properties. *J. Phys. Chem.* **1994**, *98*, 9165–9169. [CrossRef]

31. Sahu, N.; Gadre, S.R. Molecular Tailoring Approach: A Route for ab Initio Treatment of Large Clusters. *Acc. Chem. Res.* **2014**, *47*, 2739–2747. [CrossRef]

32. Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment molecular orbital method: An approximate computational method for large molecules. *Chem. Phys. Lett.* **1999**, *313*, 701–706. [CrossRef]

33. Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment molecular orbital method: Application to polypeptides. *Chem. Phys. Lett.* **2000**, *318*, 614–618. [CrossRef]
34. Fedorov, D.G.; Kitaura, K. Theoretical development of the fragment molecular orbital (FMO) method. In *Modern Methods for Theoretical Physical Chemistry and Biopolymers*; Starikov, E.B., Lewis, J.P., Tanaka, S., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; Chapter 1; pp. 3–38. [CrossRef]

35. Fedorov, D.G.; Kitaura, K. Theoretical Background of the Fragment Molecular Orbital (FMO) Method and Its Implementation in GAMESS. In *The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems*; Fedorov, D.G., Kitaura, K., Eds.; CRC Press-Taylor & Francis Group: Boca Raton, FL, USA, 2009; Chapter 2; pp. 5–36. [CrossRef]

36. Huang, L.; Massa, L.; Karle, J. Kernel energy method illustrated with peptides. *Int. J. Quantum Chem.* **2005**, *103*, 808–817. [CrossRef]

37. Huang, L.; Massa, L.; Karle, J. Kernel energy method applied to vesicular stomatitis virus nucleoprotein. *Proc. Natl. Acad. Sci. USA* **2009**, *106*, 1731–1736. [CrossRef] [PubMed]

38. Huang, L.; Bohorquez, H.; Matta, C.F.; Massa, L. The Kernel Energy Method: Application to Graphene and Extended Aromatics. *Int. J. Quantum Chem.* **2011**, *111*, 4150–4157. [CrossRef]

39. Zhang, D.W.; Zhang, J.Z.H. Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein-molecule interaction energy. *J. Chem. Phys.* **2003**, *119*, 3599–3605. [CrossRef]

40. Gao, A.M.; Zhang, D.W.; Zhang, J.Z.H.; Zhang, Y. An efficient linear scaling method for ab initio calculation of electron density of proteins. *Chem. Phys. Lett.* **2004**, *394*, 293–297. [CrossRef]

41. He, X.; Zhang, J.Z.H. The generalized molecular fractionation with conjugate caps/molecular mechanics method for direct calculation of protein energy. *J. Chem. Phys.* **2006**, *124*, 184703. [CrossRef]

42. Li, S.; Li, W.; Fang, T. An Efficient Fragment-Based Approach for Predicting the Ground-State Energies and Structures of Large Molecules. *J. Am. Chem. Soc.* **2005**, *127*, 7215–7226. [CrossRef]

43. Walker, P.D.; Mezey, P.G. Ab Initio Quality Electron Densities for Proteins: A MEDLA Approach. *J. Am. Chem. Soc.* **1994**, *116*, 12022–12032. [CrossRef]

44. Exner, T.E.; Mezey, P.G. Ab Initio-Quality Electrostatic Potentials for Proteins: An Application of the ADMA Approach. *J. Phys. Chem. A* **2002**, *106*, 11791–11800. [CrossRef]

45. Exner, T.E.; Mezey, P.G. Ab initio quality properties for macromolecules using the ADMA approach. *J. Comput. Chem.* **2003**, *24*, 1980–1986. [CrossRef] [PubMed]

46. Breneman, C.M.; Thompson, T.R.; Rhem, M.; Dung, M. Electron density modeling of large systems using the transferable atom equivalent method. *Comput. Chem.* **1995**, *19*, 161–179. [CrossRef]

47. Chang, C.; Bader, R.F.W. Theoretical construction of a polypeptide. *J. Phys. Chem.* **1992**, *96*, 1654–1662. [CrossRef]

48. Bader, R.F.W.; Martín, F.J. Interdeterminancy of basin and surface properties of an open system. *Can. J. Chem.* **1998**, *76*, 284–291. [CrossRef]

49. Matta, C.F. Theoretical Reconstruction of the Electron Density of Large Molecules from Fragments Determined as Proper Open Quantum Systems: The Properties of the Oripavine PEO, Enkephalins, and Morphine. *J. Phys. Chem. A* **2001**, *105*, 11088–11101. [CrossRef]

50. Bader, R.F.W. *Atoms in Molecules: A Quantum Theory*; Oxford University Press: Oxford, UK, 1990.

51. Stoll, H.; Wagenblast, G.; Preuβ, H. On the use of local basis sets for localized molecular orbitals. *Theor. Chim. Acta* **1980**, *57*, 169–178. [CrossRef]

52. Fornili, A.; Sironi, M.; Raimondi, M. Determination of Extremely Localized Molecular Orbitals and Their Application to Quantum Mechanics/Molecular Mechanics Methods and to the Study of Intramolecular Hydrogen Bonding. *J. Mol. Struct.* **2003**, *632*, 157–172. [CrossRef]

53. Sironi, M.; Genoni, A.; Civera, M.; Pieraccini, S.; Ghitti, M. Extremely localized molecular orbitals: Theory and applications. *Theor. Chem. Acc.* **2007**, *117*, 685–698. [CrossRef]

54. Genoni, A.; Sironi, M. A novel approach to relax extremely localized molecular orbitals: The extremely localized molecular orbital-valence bond method. *Theor. Chem. Acc.* **2004**, *112*, 254–262. [CrossRef]

55. Genoni, A.; Fornili, A.; Sironi, M. Optimal Virtual Orbitals to Relax Wave Functions Built Up with Transferred Extremely Localized Molecular Orbitals. *J. Comput. Chem.* **2005**, *26*, 827–835. [CrossRef]

56. Genoni, A.; Ghitti, M.; Pieraccini, S.; Sironi, M. A novel extremely localized molecular orbitals based technique for the one-electron density matrix computation. *Chem. Phys. Lett.* **2005**, *415*, 256–260. [CrossRef]

57. Genoni, A.; Merz, K.M., Jr.; Sironi, M. A Hylleras functional based perturbative technique to relax extremely localized molecular orbitals. *J. Chem. Phys.* **2008**, *129*, 054101. [CrossRef]

58. Sironi, M.; Ghitti, M.; Genoni, A.; Saladino, G.; Pieraccini, S. DENPOL: A new program to determine electron densities of poly-peptides using extremely localized molecular orbitals. *J. Mol. Struct.* **2009**, *898*, 8–16. [CrossRef]

59. Meyer, B.; Guillot, B.; Ruiz-Lopez, M.F.; Genoni, A. Libraries of Extremely Localized Molecular Orbitals. 1. Model Molecules Approximation and Molecular Orbitals Transferability. *J. Chem. Theory Comput.* **2016**, *12*, 1052–1067. [CrossRef]

60. Meyer, B.; Guillot, B.; Ruiz-Lopez, M.F.; Jelsch, C.; Genoni, A. Libraries of Extremely Localized Molecular Orbitals. 2. Comparison with the Pseudoatoms Transferability. *J. Chem. Theory Comput.* **2016**, *12*, 1068–1081. [CrossRef] [PubMed]

61. Meyer, B.; Genoni, A. Libraries of Extremely Localized Molecular Orbitals. 3. Construction and Preliminary Assessment of the New Databanks. *J. Phys. Chem. A* **2018**, *122*, 8965–8981. [CrossRef]
62. Genoni, A.; Bučinský, L.; Claiser, N.; Contreras-García, J.; Dittrich, B.; Dominiak, P.M.; Espinosa, E.; Gatti, C.; Giannozzi, P.; Gillet, J.-M.; et al. Quantum Crystallography: Current Developments and Future Perspectives. *Chem. Eur. J.* **2018**, *24*, 10881–10905. [CrossRef]

63. Grabowsky, S.; Genoni, A.; Bürgi, H.-B. Quantum crystallography. *Chem. Sci.* **2017**, *8*, 4159–4176. [CrossRef] [PubMed]

64. Genoni, A.; Macchi, P. Quantum Crystallography in the Last Decade: Developments and Outlooks. *Crystals* **2020**, *10*, 473. [CrossRef]

65. Grabowsky, S.; Genoni, A.; Thomas, S.P.; Jayatilaka, D. The Advent of Quantum Crystallography: Form and Structure Factors from Quantum Mechanics for Advanced Structure Refinement and Wavefunction Fitting. In *21st Century Challenges in Chemical Crystallography II. Structure and Bonding*; Springer International Publishing: Berlin/Heidelberg, Germany, 2020; pp. 65–144. [CrossRef]

66. Macchi, P. The connubium between crystallography and quantum mechanics. *Crystallogr. Rev.* **2020**, *26*, 209–268. [CrossRef]

67. Massa, L.; Matta, C.F. Quantum crystallography: A perspective. *J. Comput. Chem.* **2018**, *39*, 1021–1028. [CrossRef]

68. Jayatilaka, D.; Dittrich, B. X-ray structure refinement using aspherical atomic density functions obtained from quantum mechanical calculations. *Acta Crystallogr. Sect. A* **2008**, *64*, 383–393. [CrossRef] [PubMed]

69. Capelli, S.C.; Bürgi, H.-B.; Dittrich, B.; Grabowsky, S.; Jayatilaka, D. Hirshfeld atom refinement. *IUCrJ* **2014**, *1*, 361–379. [CrossRef]

70. Woińska, M.; Grabowsky, S.; Dominiak, P.M.; Woźniak, K.; Jayatilaka, D. Hydrogen atoms can be located accurately and pre-cisely by x-ray crystallography. *Sci. Adv.* **2016**, *2*, e1600192. [CrossRef]

71. Fugel, M.; Jayatilaka, D.; Hupf, E.; Overgaard, J.; Hathwar, V.R.; Macchi, P.; Turner, M.J.; Howard, J.A.K.; Dolomanov, O.V.; Puschmann, H.; et al. Probing the accuracy and precision of Hirshfeld atom refinement with HARt interfaced with Olex2. *IUCrJ* **2018**, *5*, 32–44. [CrossRef] [PubMed]

72. Wieduwilt, E.K.; Macetti, G.; Malaspina, L.A.; Jayatilaka, D.; Grabowsky, S.; Genoni, A. Post-Hartree-Fock methods for Hirshfeld atom refinement: Are they necessary? Investigation of a strongly hydrogen-bonded molecular crystal. *J. Mol. Struct.* **2020**, *1209*, 127934. [CrossRef]

73. Kleemiss, F.; Dolomanov, O.V.; Bodensteiner, M.; Peyerimhoff, N.; Midgley, L.; Borhis, L.J.; Genoni, A.; Malaspina, L.A.; Jayatilaka, D.; Spencer, J.L.; et al. Accurate Crystal Structures and Chemical Properties from NoSpherA2. *Chem. Sci.* **2021**, *12*, 1675–1692. [CrossRef]

74. Malaspina, L.A.; Wieduwilt, E.K.; Bergmann, J.; Kleemiss, F.; Meyer, B.; Ruiz-López, M.F.; Pal, R.; Hupf, E.; Beckmann, J.; Piltz, R.O.; et al. Fast and Accurate Quantum Crystallography: From Small to Large, from Light to Heavy. *J. Phys. Chem. Lett.* **2019**, *10*, 6973–6982. [CrossRef] [PubMed]

75. Johnson, E.R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A.J.; Yang, W. Revealing Noncovalent Interactions. *J. Am. Chem. Soc.* **2010**, *132*, 6498–6506. [CrossRef]

76. Contreras-García, J.; Johnson, E.R.; Keinan, S.; Chaudret, R.; Piquemal, J.-P.; Beratan, D.N.; Yang, W. NCIPLOT: A Program for Plotting Noncovalent Interaction Regions. *J. Chem. Theory Comput.* **2011**, *7*, 625–632. [CrossRef]

77. Lefebvre, C.; Rubez, G.; Khartabil, H.; Boisson, J.-C.; Contreras-García, J.; Hénon, E. Accurately Extracting the Signature of In-termolecular Interactions Present in the NCI Plot of the Reduced Density Gradient versus Electron Density. *Phys. Chem. Chem. Phys.* **2017**, *19*, 17928–17936. [CrossRef] [PubMed]

78. Lefebvre, C.; Khartabil, H.; Boisson, J.-C.; Contreras-García, J.; Piquemal, J.-P.; Hénon, E. The Independent Gradient Model: A New Approach for Probing Strong and Weak Interactions in Molecules from Wave Function Calculations. *ChemPhysChem* **2018**, *19*, 724–735. [CrossRef]

79. Ponce-Vargas, M.; Lefebvre, C.; Boisson, J.-C.; Hénon, E. Atomic Decomposition Scheme of Noncovalent Interactions Applied to Host–Guest Assemblies. *J. Chem. Inf. Model.* **2020**, *60*, 268–278. [CrossRef]

80. Klein, J.; Khartabil, H.; Boisson, J.-C.; Contreras-García, J.; Piquemal, J.-P.; Hénon, E. New Way for Probing Bond Strength. *J. Phys. Chem. A* **2020**, *124*, 1850–1860. [CrossRef] [PubMed]

81. Arias-Olivares, D.; Wieduwilt, E.K.; Contreras-García, J.; Genoni, A. NCI-ELMO: A New Method to Quickly and Accurately Detect Noncovalent Interactions in Biosystems. *J. Chem. Theory Comput.* **2019**, *15*, 6456–6470. [CrossRef]

82. Wieduwilt, E.K.; Boisson, J.-C.; Terraneo, G.; Hénon, E.; Genoni, A. A Step toward the Quantification of Noncovalent Interactions in Large Biological Systems: The Independent Gradient Model-Extremely Localized Molecular Orbital Approach. *J. Chem. Inf. Model.* **2021**. [CrossRef] [PubMed]

83. Macetti, G.; Genoni, A. Quantum Mechanics/Extremely Localized Molecular Orbital Method: A Fully Quantum Mechanical Embedding Approach for Macromolecules. *J. Phys. Chem. A* **2019**, *123*, 9420–9428. [CrossRef] [PubMed]

84. Macetti, G.; Wieduwilt, E.K.; Assfeld, X.; Genoni, A. Localized Molecular Orbital-Based Embedding Scheme for Correlated Methods. *J. Chem. Theory Comput.* **2020**, *16*, 3578–3596. [CrossRef]

85. Macetti, G.; Genoni, A. Quantum Mechanics/Extremely Localized Molecular Orbital Embedding Strategy for Excited States: Coupling to Time-Dependent Density Functional Theory and Equation-of-Motion Coupled Cluster. *J. Chem. Theory Comput.* **2020**, *16*, 7490–7506. [CrossRef]

86. Macetti, G.; Wieduwilt, E.K.; Genoni, A. QM/ELMO: A Multi-Purpose Fully Quantum Mechanical Embedding Scheme Based on Extremely Localized Molecular Orbitals. *J. Phys. Chem. A*.
87. Wieduwilt, E.K.; Macetti, G.; Genoni, A. Climbing Jacob's Ladder of Structural Refinement: Introduction of a Localized Molecular Orbital-Based Embedding for Accurate X-ray Determinations of Hydrogen Atom Positions. *J. Phys. Chem. Lett.* **2021**, *12*, 463-471. [CrossRef] [PubMed]

88. Furukawa, H.; Cordova, K.E.; O'Keeffe, M.; Yaghi, O.M. The Chemistry and Applications of Metal-Organic Frameworks. *Science* **2013**, *341*, 1230444. [CrossRef] [PubMed]

89. Ockwig, N.W.; Delgado-Friedrichs, O.; O'Keeffe, M.; Yaghi, O.M. Reticular Chemistry: Occurrence and Taxonomy of Nets and Grammar for the Design of Frameworks. *Acc. Chem. Res.* **2005**, *38*, 176-182. [CrossRef] [PubMed]

90. Yaghi, O.M.; O'Keeffe, M.; Ockwig, N.W.; Chae, H.K.; Eddaoudi, M.; Kim, J. Reticular Synthesis and the Design of New Materials. *Nature* **2003**, *423*, 705-714. [CrossRef]

91. Zagorodniy, K.; Seifert, G.; Hermann, H. Metal-organic frameworks as promising candidates for future ultralow-k dielectrics. *Appl. Phys. Lett.* **2010**, *97*, 251905. [CrossRef]

92. Boys, S.F. Construction of Some Molecular Orbitals to Be Approximately Invariant for Changes from One Molecule to Another. *Rev. Mod. Phys.* **1960**, *32*, 296-299. [CrossRef]

93. Foster, J.M.; Boys, S.F. Canonical Configurational Interaction Procedure. *Rev. Mod. Phys.* **1960**, *32*, 300-302. [CrossRef]

94. Edmiston, C.; Ruedenberg, K. Localized Atomic and Molecular Orbitals. *Rev. Mod. Phys.* **1963**, *35*, 457-464. [CrossRef]

95. Edmiston, C.; Ruedenberg, K. Localized Atomic and Molecular Orbitals. II. *J. Chem. Phys.* **1965**, *43*, S97-S116. [CrossRef]

96. Von Niessen, W. Density Localization of Atomic and Molecular Orbitals. I. *J. Chem. Phys.* **1972**, *56*, 4290. [CrossRef]

97. Pipek, J.; Mezey, P.G. A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions. *J. Chem. Phys.* **1989**, *90*, 4916-4926. [CrossRef]

98. McWeeny, R. The density matrix in many-electron quantum mechanics. I. Generalized product functions. Factorization and physical interpretation of the density matrices. *Proc. R. Soc. London. Ser. A Math. Phys. Sci.* **1959**, *253*, 242-259. [CrossRef]

99. Adams, W.H. On the Solution of the Hartree-Fock Equation in Terms of Localized Orbitals. *J. Chem. Phys.* **1961**, *34*, 89. [CrossRef]

100. Huzinaga, S.; Cantu, A. A Theory of Separability of Many-Electron Systems. *J. Chem. Phys.* **1971**, *55*, 5543-5549. [CrossRef]

101. Gilbert, T.L. Multiconfiguration Self-Consistent-Field Theory for Localized Orbitals. II. Overlap constraints, Lagrangian multipliers, and the screened interaction field. *J. Chem. Phys.* **1974**, *60*, 3835-3844. [CrossRef]

102. Matsuoka, O. Expansion methods for Adams-Gilbert equations. I. Modified Adams-Gilbert equation and common and fluctuating basis sets. *J. Chem. Phys.* **1977**, *66*, 1245-1254. [CrossRef]

103. Smits, G.F.; Altona, C. Calculation and properties of non-orthogonal, strictly local molecular orbitals. *Theor. Chem. Acc.* **1985**, *67*, 461-475. [CrossRef]

104. Francisco, E.; Martín Pendás, A.; Adams, W.H. Generalized Huzinaga building-block equations for nonorthogonal electronic groups: Relation to the Adams-Gilbert theory. *J. Chem. Phys.* **1992**, *97*, 6504-6508. [CrossRef]

105. Ordejón, P.; Drabold, D.A.; Grumbach, M.P.; Martin, R.M. Unconstrained minimization approach for electronic computations that scales linearly with system size. *Phys. Rev. B* **1993**, *48*, 14646-14649. [CrossRef]

106. Couty, M.; Bayse, C.A.; Hall, M.B. Extremely localized molecular orbitals (ELMO): A non-orthogonal Hartree-Fock method. *Theor. Chem. Acc.* **1997**, *97*, 96-109. [CrossRef]

107. Philipp, D.M.; Friesner, R.A. Mixed Ab Initio QM/MM Modeling Using Frozen Orbitals and Tests with Alanine Dipeptide and Tetrapeptide. *J. Comput. Chem.* **1999**, *20*, 1468-1494. [CrossRef]

108. Chui, S.S.-Y.; Lo, S.M.-F.; Charmant, J.P.H.; Orpen, A.G.; Williams, I.D. Chemically Functionalizable Nanoporous Material [Cu₃(TMA)₂(H₂O)₃]ₘ. *Science* **1999**, *283*, 1148-1150. [CrossRef] [PubMed]

109. Scatena, R.; Guntern, Y.T.; Macchi, P. Electron Density and Dielectric Properties of Highly Porous MOFs: Binding and Mobility of Guest Molecules in Cu₃(BTC)₂ and Zn₃(BTC)₂. *J. Am. Chem. Soc.* **2019**, *141*, 9382-9390. [CrossRef]

110. Guest, M.F.; Bush, I.J.; van Dam, H.J.J.; Sherwood, P.; Thomas, J.M.H.; van Lenthe, J.H.; Havenith, R.W.A.; Kendrick, J. The GAMESS-UK electronic structure package: Algorithms, developments and applications. *Mol. Phys.* **2005**, *103*, 719-747. [CrossRef]

111. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, USA, 2009.

112. Macchi, P.; Sironi, A. Chemical bonding in transition metal carbonyl clusters: Complementary analysis of theoretical and ex-perimental electron densities. *Coord. Chem. Rev.* **2003**, 238-239, 383-412. [CrossRef]

113. Jahn, H.A.; Teller, E. Stability of polyatomic molecules in degenerate electronic states-I—Orbital degeneracy. *Proc. R. Soc. London. Ser. A Math. Phys. Sci.* **1937**, *161*, 220-235. [CrossRef]