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

: 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 ﬁrst 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 mechan-ics/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 3 (BTC) 2 } n , 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.


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][4][5][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) developed by Wesolowski and coworkers [12][13][14] or the projection-based embedding (PbE) approach jointly devised by the Manby and Miller research groups [15][16][17][18][19][20][21][22][23][24].
Fragmentation methods are also worth mentioning in this context. The idea is to partition 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 afterwards combined to obtain the global wavefunction and/or electron density of the target macrosystem. Examples are the pioneering Divide and Conquer approach, initially proposed 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][28][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 "fragment 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][33][34][35], the kernel energy method (KEM) [36][37][38], and the molecular fractionation with conjugate caps (MFCC) technique in its different variants [39][40][41][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 derives 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) strategy [43] and afterwards the adjustable density matrix assembler (ADMA) approach [44,45], respectively consisting of libraries of fuzzy electron densities and density matrices associated 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][52][53]. Thanks to their extreme localization, they can be easily exported from small model molecules to a macrosystem, 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][54][55][56][57][58][59][60] stimulated the compilation of ELMO libraries [59][60][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][63][64][65][66][67] the ELMO databanks have been coupled with the emerging Hirshfeld atom refinement (HAR) approach [68][69][70][71][72][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][78][79][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][84][85][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][89][90][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 extraframework 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.

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][93][94][95][96][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][99][100][101][102][103][104][105][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. 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 Figures 1G and 1H) that overall describe four bond electrons. At the same time, the oxygen atom fragment is associated with three ELMOs, two of which (see Figures 1I and 1J) describe the four lone pair electrons, while the remaining one (not depicted in Figure 1) describes the two core electrons. 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 nonterminal) [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 frag- 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 databanks 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.

ELMO-Protocol for MOFs
The metal organic framework chosen for our preliminary study was {Cu 3 (BTC) 2 } 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 tailormade ELMOs computed on suitable model molecules of the MOF subunits. ments 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 databanks 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.

ELMO-Protocol for MOFs
The metal organic framework chosen for our preliminary study was {Cu3(BTC)2}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 tailormade ELMOs computed on suitable model molecules of the MOF subunits. 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. 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].

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.  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 3 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-  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(λ 2 )ρ < 5.0 a.u. 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 following 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 secondary 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 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 z 2 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 (∇ 2 ρ bcp > 0), but with a dominance of the potential energy density over the kinetic energy density ( V bcp /G bcp > 1), hence with a negative total energy density (H bcp = −K bcp = V bcp + G bcp < 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. (a) ρ bcp and ∇ 2 ρ bcp are, respectively, the electron density (e/bohr 3 ) and the Laplacian of the electron density (e/bohr 5 ) 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 3 ); and DI is the delocalization index (electron pairs shared between two atoms). 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.
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 following 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 secondary 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.
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.  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.

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 threedimensional 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-4 352/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.

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