Investigating solvent effects on the magnetic properties of molybdate ions (MoO$_{4}^{2-}$) with relativistic embedding

We investigate the ability of mechanical and electronic density functional theory (DFT)-based embedding approaches to describe the solvent effects on nuclear magnetic resonance (NMR) shielding constants of the $^{95}$Mo nucleus in the molybdate ion in aqueous solution. From the description obtained from calculations with two- and four-component relativistic Hamiltonians, we find that for such systems spin-orbit coupling effects are clearly important for absolute shielding values, but for relative quantities a scalar relativistic treatment provides a sufficient estimation of the solvent effects. We find that the electronic contributions to the solvent effects are relatively modest yet decisive to provide a more accurate magnetic response of the system, when compared to reference supermolecular calculations. We analyze the errors in the embedding calculations by statistical methods as well as through a real-space representation of NMR shielding densities, which are shown to provide a clear picture of the physical processes at play.


I. INTRODUCTION
Nuclear magnetic resonance (NMR) spectroscopy measures the interaction between the magnetic moments of nuclei and an applied external field, screened by the electrons of the system. NMR is extremely useful for characterizing molecules and materials since it provides detailed information of the local chemical environment around the responding nuclei, and does so in a non-destructive manner. In addition to that, it probes species in their electronic ground state and, in doing so, provides more direct information on the interactions in that state compared to other techniques.
The screening by the electron cloud of the magnetic interactions between a given atomic isotope K and the applied external magnetic field, is represented by the socalled NMR shielding constant, σ K , which can be calculated as a second-order derivative of the energy of the molecular system with respect to the magnetic dipole moment of that nucleus ( m K ) and an external magnetic field ( B) at the zero-perturbation limit: In most cases, what is measured in experiments is not σ K but rather a signal relative to a chosen reference species, the chemical shift δ K , though in recent years experimental advances have revived the interest in determining absolute shielding scales [1][2][3]. NMR spectra can be quite difficult to interpret. The difficulties may come from the complexity due to the size * andre.gomes@univ-lille.fr of composition of the system (species in solution, disordered materials, etc.), the broadening of signals due to quadrupole interactions (for nuclei with total angular momentum larger than I = 1/2) or a combination of all of these. [4] Because of these difficulties, theoretical modeling has become indispensable for the interpretation of experiments [5][6][7]. Furthermore, for certain systems exhibiting a wide range of chemical shifts, or long acquisition times, theory is crucial to guide experiments by providing the approximate regions of the spectra in which to search for signals.
In this paper we are interested in exploring the use of quantum embedding [8][9][10][11] to describe the solvent effects on the NMR shielding constants of the 95 Mo on the molybdate dianion in aqueous solution. The molybdate moiety is arguably the simplest experimentally relevant Mo-containing system that can be studied, and it is often used as an experimental reference system in the determination of 95 Mo chemical shifts (in the form of an aqueous solution of sodium molybdate, Na 2 MoO 4 [12]. As such, it can be seen as ideal test system to evaluate theoretical approaches that can be applied to the modeling of other, more complicated species. Examples of more complex systems can be found in different classes of molybdenum oxides containing the MoO 3 and MoO 4 moieties. These oxides are found as component of catalysts [13][14][15], as fission products [16], and also make up materials that can be used for applications ranging from nonlinear optical materials [17] to glasses for vitrification of nuclear waste [18]. In most cases, the interest lies in understanding how the Mo atom modifies the properties (optical, mechanical etc) of the starting material, and that is where the characterization by NMR becomes interesting.
For materials and complex systems, often treated with band structure methods, the gauge-including projector augmented wave (GIPAW) [6,7] approach, coupled to density functional theory (DFT) calculations, is the de facto standard for NMR calculations useful to experimentalists. Quantum embedding, on the other hand, represent a class of theoretical methods, often used in combination to approaches developed for discrete systems (based on DFT but also on more accurate treatments based upon wavefunction theory, WFT). In it a given system is treated quantum mechanically but on the basis of a set of interacting subsystems, for which one can tune the accuracy of the description of each subsystem (by choosing for instance a DFT or WFT treatment), according to a desired compromise between accuracy and computational cost. As such, one can describe the interactions between the species of interest and its surroundings at a fraction of the cost of applying standard DFT or WFT approaches to the whole system. Contrary to other embedding schemes [19][20][21], such as those in which only the changes in structure due to the environment are considered (mechanical embedding), or where the environment is represented classically (QM/MM) [22], quantum embedding has not yet been extensively used for NMR properties, and this is particularly the case for the frozen density embedding or subsystem DFT methods [23][24][25][26]. This is unfortunate since these approaches have the advantage of allowing for a seamless combination of different electronic structure approaches (DFT and WFT), but also the combination of different Hamiltonians, such as those taking into account relativistic effects.
Relativistic effects [27][28][29][30] (scalar relativistic effects and spin-orbit coupling) are now recognized as important (if not essential) to approach experimental accuracy in electronic structure calculations. While this is undisputed in the case of heavy elements (Z > 31), it is increasingly the case that the effects of relativity are recognized and accounted for even for light systems (first and second-row atoms) when the molecular properties of interest involve a description of the electronic structure near the nucleus. This is the case of NMR properties [31,32].
We expect that quantum embedding, if shown to be reliable, can provide a flexible framework and enable more sophisticated simulations, either by complementing DFT-only approaches such as GIPAW (by easily letting one explore the use of different relativistic Hamiltonians) or by allowing the use of more accurate methods, such as those based on the coupled cluster ansatz [33,34], thus paving the way for a quantitative leap in the interpretation of NMR spectra. As the hydrated MoO 2− 4 molybdate dianion appears to have well-structured hydrations shells [12], we are particularly interested in assessing how well quantum embedding can effectively reduce the size of the explicit quantum model, by accounting for hydration effects.
The paper is organized as follows: in section II we outline the key theoretical aspects underpinning the work. This is followed in section III by a summary of the computational details of our calculations, with a particular emphasis in detailing which structural models are considered and introducing a shorthand notation to refer to these. In section IV we discuss the ability of the different embedding approaches to represent the molybdate ion in solution, the influence of the different relativistic Hamiltonians on the NMR shielding constants, and characterize (including in real space) the errors of the different embedding methods. We conclude and offer perspectives for future work in section V.

II. THEORY
In what follows, we employ the standard notation for molecular orbital indices (with i, j, . . . corresponding to occupied, a, b, . . . -to virtual and p, q, . . . -to general molecular orbitals) and the summation convention over repeated indices. We use the SI-based atomic units ( = m e = e = 1/(4π 0 ) = 1) [35]. Employing the second-quantization formalism and an exponential parametrization of the closed-shell ground state with the orbital-rotation operator,κ = pq κ pq p † q, the derivative in Eq. 1 can be further expressed as where all perturbing fields are collected in ε. Partial derivatives of the orbital rotation parameters ({κ pq }) with respect to the perturbing field are optimized in the linear response equations, which in the static-field regime can be compactly presented as: with E [2] 0 , E [1] ε1 and X ε1 corresponding to the electronic Hessian, the property gradient and the solution vector (X ε1 = (∂κ pq /∂ε 1 ) ε=0 ), respectively [36]. The two former quantities require the formulation of various derivatives of Fock or Kohn-Sham matrices, therefore they are strictly connected to the formalism employed in the calculations.
In this work, we exploit the computational model based on relativistic Hamiltonians [37] and on the spindensity functional theory (SDFT) [38][39][40]. The closedshell system is described by the electron density (ρ 0 ) and the spin-density calculated in a non-collinear fashion as a norm of the spin-magnetization vector (s = |ρ µ |, µ ∈ {x, y, z}). Following previous works [26,[41][42][43], we collect the electron density and the spin-magnetization vector into one variable -a general density component, This framework becomes slightly modified in the embedding situation. The DFT-subsystem-based approaches rely on the partitioning of the total system into (interacting) subsystems [8,9,44], which is realized by expressing the electron density of the full system as a sum of electron densities of subsystems [10,[45][46][47] (the same applies to the spin density if it is included in the formalism), Consequently, the total energy may be decomposed into contributions corresponding to energies of subsystems and the interaction energy term which depends on the (spin) densities of all subsystems, The analytical formula for the interaction energy term can be derived from the (S)DFT expression for the energy of the full system and reads as: where, in addition to terms describing the interaction of the electron density of a subsystem with the electron density and the nuclear potential of another subsystem and a term related to the nuclear repulsion energy between subsystems (first four terms in Eq. 7), the nonadditive exchange-correlation (E nadd xc ) and kinetic energy (T nadd s ) contributions emerge. These non-additive contributions [48] are defined as: (8) and they depend on the general density components of all subsystems. In typical FDE calculations, one subsystem is chosen as active, while another subsystem constitute its environment. In this setting, the active density (e.g. ρ I k ) is optimized by solving the so-called Kohn-Sham (KS) equations for a constrained electron density (KSCED) [49], in which an effective KS potential known from a standard (S)DFT is augmented by the embedding potential responsible for describing the interaction of an active subsystem with the environment (here: It is also possible to relax the density of the environment by solving the analogous KSCED equations set up for ρ II k . Repeating this procedure in an iterative manner, known as the freeze-thaw [50] (FnT) cycle allows to optimize the densities of all subsystems.
At this point it should be noted that although FDE is formally exact, its practical realization requires approximation to the non-additive exchange-correlation and kinetic energies. This is typically done with approximate density functionals, which -especially for the kinetic energy part -have limited accuracy [8,[51][52][53][54]. From the presence of the interaction energy term in the total energy expression (Eq. 5) and from its functional dependence on the densities of all subsystems, it becomes evident that the quantities which in the linear response theory are formulated as various energy derivatives (E [2] 0 , E [1] ε1 in Eq. 4) now have to be augmented with analogous derivatives of E int .
These additional contributions, derived by using the chain rule, involve the first-and second-order derivatives of E int with respect to the densities of subsystems -the embedding potential (Eq. 9) and the embedding kernel (see [26] for details). In effect, the property gradient and the electronic Hessian (Eq. 4) are elegantly subdivided into contributions from isolated subsystems and terms related to the interaction input [48]. In practice, these derivatives can be coded analogously to the exchangecorrelation terms in standard (S)DFT technique. There is, however, an additional difficulty in the molecular property calculations, specific to the embedding setting. Namely, one needs to account for the fact that the property of interest arises as a response of the full systemcomprising of all subsystems -to certain external perturbations. In particular, for the NMR shielding tensor of a nucleus K this can be symbolically expressed as Assuming that the nucleus K belongs to subsystem I, the term σ K,II αβ can be interpreted as an effect that the perturbation due to B α in subsystem II has in the point of the center of nucleus K -an idea analogous to the concept of the nucleus independent chemical shift (NICS) [55]. A final complication in the formulation of the NMR shielding tensor arises as a consequence of the dependence of this property operator (more specifically -the Zeeman term) on the choice of the gauge origin, which leads to unphysical results in incomplete basis set regime. The typical remedy is to use the London atomic orbitals (LAOs) [56,57] instead of atomic orbitals, however with the price of the need to calculate additional terms in the property gradient [26,43,58].
The NMR shielding tensor (Eqs. 1, 3) can also be determined from numerical integration of the NMR shielding density [59,60]. This property density can be calculated (analytically) in real space by using the relation to the magnetically-induced current density, j B -also a second-order tensor quantity, which in the relativistic framework [61,62] reads The NMR shielding density is an integrand in the following expression: This reformulation of σ K αβ has multiple advantages. First, it allows to recalculate the NMR shielding values (therefore serves as a test of a fully-analytical approach), secondly it opens the possibility to visualize the NMR shielding density on meshes, which has already proved to be an invaluable asset in the post-production analysis process [26]. Additional advantage of this formulation is that the NMR shielding density of a nucleus K can easily be calculated in any point in space, even if that point and the center of nucleus K belong to different subsystems, what permits to evaluate the second term of Eq. 10 [23].

III. COMPUTATIONAL DETAILS
The solvated molybdate structures used in this work are taken from the Car-Parrinello molecular dynamics (CPMD [63]) calculations of Nguyen and coworkers [12]. We considered 517 snapshots, each containing the molybdate ion and 20 water molecules. For additional details, the reader is referred to the aforementioned paper. From these we defined four basic structural models: • a supermolecular system, containing the molybdate ion and all water molecules; • a subsystem embedding model, where the molybdate ion (and depending on the case, selected water molecules) is the active subsystem and the (remaining) water molecules make up the environment, but in which a given number of water molecules nearest to the active subsystem are relaxed through freezethaw cycles while the others remain frozen. We will refer to this model as FnT.
The electron densities and electrostatic potential for the environment are obtained either from calculations on individual water molecules (and we refer to a "fragmented" environment) or by these forming a single subsystem (and we refer to a "grouped" environment); • a frozen density embedding model, in which all water molecules in the environment are kept frozen.
We will refer to this model as FDE, and the same notations for the environment ("fragmented" or "grouped") as for FnT is used; • an isolated molybdate ion. As the structures of the latter are not obtained in the gas phase but rather through the CPMD calculations, this model is equivalent to that of mechanical embedding, and as such we will refer to it as "mechanical embedding".
Using the labels "A" for the molybdate ion and "B" for the water molecules, and indices g and f to denote the "grouped" or "f ragmented" environment, we shall use the following shorthand notation to represent these four cases: [Active | Relaxed | F rozen]. This means, for example, that the supermolecular model is denoted by [A + 20B|0|0] (all species are in the active fragment), the mechanical embedding by [A|0|0], and a FnT model in which 20 molecules make up the environment but in which 11 are relaxed (and all are taken as individual fragments), is denoted by The execution of all calculations, including the preparation of the different fragments, has been handled with the PyADF [64] scripting framework. In all calculations, we employ Becke integration grids, and a gaussian nuclear charge distribution [65] to describe the finite nuclear volume [66], PBE [67,68] exchange-correlation functional. For subsystem embedding and FDE embedding calculations (see below), we have employed the PW91k [69] and PBE [67,68] functionals for calculating the non-additive kinetic and exchange-correlation contributions, respectively, and used a monomer expansion for the basis sets. In the case of subsystem embedding, we have performed 5 freeze-thaw cycles. In the NMR calculations we employ London orbitals.

A. One-and two-component calculations
The scalar relativistic (SR) and spin-orbit (SO) ZORA [70][71][72][73] calculations have been performed with the ADF [74] code, using the TZ2P basis sets [75]. The NMR calculations have been performed with the NMR program.

B. Four-component calculations
The calculations employing the Dirac-Coulomb (DC) Hamiltonian have been performed with the Dirac code [76], revision 54ab939. For the Mo atom the dyall.cv3z basis [77] is used, and for all lighter atoms the aug-cc-pVTZ [78] basis is used. To aid in the convergence of the SCF procedure, in all calculations we have used the atomic start procedure, described in the DIRAC documentation.
For the subsystem embedding and FDE embedding calculations NMR calculations, we remain in the uncoupled response approximation, but take into account contributions from the spin density according to Ref. [26]. We note that unless otherwise noted, the embedding calculations make use of results of embedded ZORA calculations with ADF (embedding potentials, electrostatic potentials and frozen electron densities and gradients of the density for the environment).
The NMR shielding densities used in the analysis of the embedding methods (see section IV C) have been obtained with the analysis module of the Dirac code, and are based on four-component calculations with the DC Hamiltonian and equation 12. In the case of embedded calculations, following equation 10, it was necessary to also perform calculations for the environment using the position of the Mo atom as the r K position.

C. Plotting
All graphs have been prepared with Matplotlib [79] python library, with the exception of the volumetric density and NMR shielding density plots, for which the Mayavi [80] python library has been used.

IV. RESULTS AND DISCUSSION
In what follows we are interested in the isotropic part of the NMR shielding tensor of the 95 Mo nucleus, σ iso , which is defined as: where the indices refer to the principal axes of the tensor according to the Mason's notation [81], with σ 11 ≤ σ 22 ≤ σ 33 .

A. Minimal structural model
An embedding-based computational strategy takes its most effective form, from a computational cost perspective, when the smallest possible molecular moiety can be considered as the active subsystem, and the rest of the system as the environment. In this work such a situation would correspond, as discussed above, to having only the MoO 2− 4 ion in the active subsystem. In this subsection we explore this model.

Structural and electronic effects on NMR shieldings
We start the discussion with the analysis of the influence of the geometry of molecular complexes on σ iso based on the SR-ZORA calculations. This is summarized in Figure 1, which shows the values of σ iso obtained for a subset of 100 geometries out of the 517 considered MD simulations.
An interesting observation that emerges from this figure is that the description offered by different approximations (FnT, FDE, mechanical embedding) is not consistent for all the snapshots. While for certain geometries the electronic part of embedding is an essential component of the model in order to approach the supermolecular σ iso values, for others these contributions are negligible and the mechanical embedding is sufficient to reproduce the supermolecular baseline. This is also reflected in the similarity of the mean values of the 95 Mo NMR shielding constants calculated over all the snapshots for all the considered approximations (FnT, FDE, mechanical embedding), featured asσ iso in Table I  if electronic embedding is included in the calculations. This signifies that the data from electronic embedding calculation is statistically less disperse, therefore in this case these can be considered as more accurate for the modeling of σ iso . In order to understand these structural dependencies of σ iso values, we decided to define a measure, the mean distance between the centers of Mo atom and O atoms in MoO 2− 4 , hereafter referred to asr[Mo-O], and to plot the σ iso vsr[Mo-O] in Figure 2.
The first conclusion of this figure is that there is a significant correlation between the σ iso andr[Mo-O] for all approximations, as seen from the large values of R 2 . The largest correlation between these two quantities was found for the mechanical coupling (R 2 = 0.992) and the lowest for the supermolecular results (R 2 = 0.945), which indicates the possibility that the effect of the interaction between two subsystem on σ iso values cannot be simply attributed to the structural variations, or at least the ones represented by ther ding calculations than between supermolecule and mechanical embedding. However, the linearity of the relationship between σ iso andr[Mo-O] is interesting since the NMR shielding tensor is considered to be a local property. Yet, as we shall discuss below, there might be other long-range effects that dampen the r −2 dependence of the corresponding property operator.
Since there is a non-negligible difference between the description of the supermolecular and embedded schemes, we have investigated the correlation between r[Mo-O] and the values of the differences of σ iso obtained from embedding calculations and the supermolecular ones, denoted by ∆ in Figure 3. From these results, we see that there is in fact no correlation between the two variables (R 2 ≤ 0.214), even though we nevertheless observe that there is a much larger dispersion of the data for the mechanical embedding than for the electronic embedding, further indicating that the latter better approximates the supermolecular calculations.
These ∆ values are a measure of the error of a given embedding approximation with respect to the supermolecular standard, which therefore cannot be attributed to the molecular geometry represented by ther [Mo-O] descriptor. It remains to be tested whether this error can be associated with other structural descriptors (for instance, ones taking into account the asymmetry around the Mo nucleus in a given environment) or whether it should be attributed to more subtle effects, for instance due to spin-orbit coupling, which we discuss in the following.   In order to assess the importance of spin-orbit interactions on the the solvent effects, we begin by considering the SO-ZORA Hamiltonian. Theσ iso values and standard deviations for supermolecular and embedded calculations are shown in Table II.

Assessment of spin-orbit Hamiltonians
From these we observe that there are no major differences with respect to the SR-ZORA results (see Table I) if it comes to the magnitude of the solvent effects for all approaches. We do observe some differences in the error patterns, with both electronic embedded models now being closer to the supermolecularσ iso values than the mechanical embedding case, though once more the FDĒ σ iso is the closest to the supermolecular one. The standard deviations in each embedded model show the same pattern as in the SR-ZORA, with FnT yielding the smallest and mechanical embedding the largest std values. Tables I and II we also observe significant changes on the absoluteσ iso values, which for SO-ZORA are upshifted by around 260 ppm from the SR-ZORA values. Such a spin-orbit effect on the absolute NMR shielding values is well-known and extensively discussed in the literature [82,83], and therefore we refer the readers to the review [83] and recent examples. [82,[84][85][86][87] The nearly identical behavior of embedded SR-ZORA and SO-ZORA can be seen, first, as a manifestation of the fact that, since there are no significant spin-orbit coupling effects on the electronic structure of the solvent, the environment electron densities and electrostatic potentials will be rather similar in the two sets of calculations (SR and SO), as we can see from the rather similar standard deviations for all models considered. Second, since the Mo atom is not very heavy, spin-orbit coupling is not expected to bring about qualitative differences on the electron density of the MoO 2− 4 ion. An appealing feature of the SR-ZORA and SO-ZORA calculations is their relatively modest computational cost, which has allowed us to consider a rather large number of snapshots (including for the supermolecular system). This makes them ideal tools for exploring the importance of relativistic effects across the periodic table while considering a great number of configuration from molecular dynamics or the effect of increasing the size of of the active space. As discussed above, this allowed us to verify that SR-ZORA calculations already capture the solvent effects for the molybdate case.

Comparing the results in
The exploration of a large number of snapshots also provided data with which to guide us in reducing the   Table III. While we cannot directly compare the standard deviations from the DC results to those of tables I and II we see from the table that the DC standard deviation values are close to the SR-ZORA ones, calculated with the same reduced sample of snapshots, and therefore we expect to see a similar error behavior for DC as for ZORA, had we used the original sample size. Furthermore, as was the case when improving the ZORA Hamiltonian with the addition of spin-orbit coupling, the more accurate treatment of relativistic effects afforded by the DC Hamiltonian introduces a further upshift of about 150 ppm with respect to SO-ZORA, on the values ofσ iso , so that the difference between the SR-ZORA and DC results is now a little over 410 ppm. This is a significant difference, given the growing interest in recent years in determining absolute shielding scales. Taking into account the steep computational cost of performing 4-component supermolecular calculations over hundreds of structures from MD simulations, we consider electronic embedding to be a promising way to calculate absolute NMR shielding constants in complex environments.

B. Larger active subsystems
The results above suggest that the minimal structural model for solvated molybdate (comprising only the anion) is a rather good one for obtaining NMR shieldings in solution. In this section, we explore the effect of ex-plicitly including nearest neighbor water molecules in the calculation.
As we have found little quantitative differences in the behavior of the embedding approaches for the different Hamiltonians, and that it is computationally very expensive to perform a systematic expansion of the active subsystem while averaging over hundreds of snapshot, here we have decided to focus on selected SR-ZORA calculations.
We therefore take two snapshots (numbered 10 and 85), representing respectively structures yielding NMR shieldings close to the σ iso value and towards the tails of the distribution. Figure 4 summarizes our results. We only show results for up to 11 water molecules added to the active subsystem, due to the fact that no significant qualitative changes occur beyond this number.
For snapshot 10 (σ iso closer to σ iso ) we observe significant differences between embedding methods, with the subsystem embedding model ([A + pB|(11 − p)B f |9B f ]) yielding results which are only a few ppm below the reference supermolecular calculations. The FDE model ([A + pB|0|(20 − p)B f ])) introduces more significant errors (over 10 ppm underestimation with respect to supermolecule), which from prior findings we can can attribute to the importance of relaxing the environment species around a highly charged active subsystem. The mechanical embedding model fares the worst (over 40 ppm underestimation to supermolecule).
Adding one water molecule to the active subsystem greatly improves the results for all embedding models, and gets the electronic embedding ones in very close agreement to the reference. From two to four explicit waters, however, we see a degradation of the electronic embedding results, which now overestimate σ iso . At around six explicit water molecules, which roughly corresponds to the first solvation shell around the molybdate ion [12], the electronic embedding methods have converged to the reference result, and show no further significant variations. The mechanical embedding model, on the other hand, shows small errors but does not show converged results even after 11 explicit water molecules.
We note that a similar analysis has been carried out for 4-component calculations for snapshot 10, but restricting the reference calculations to 6 water molecules due to constraints in our computational resources (see computational details and supplemental information). As we observe the same trends in the DC calculations as in the SR-ZORA ones, we do not discuss the former explicitly.
For snapshot 85 (σ iso away from σ iso ) we also observe the embedding approaches are pretty much converged to the reference after the first solvation shell is explicitly included. However, the behavior is quite different from snapshot 10 for the smaller number of explicitly included waters. Without any water molecule in the active subsystem, all methods underestimate σ iso , and the mechanical embedding model is the worst among the three, but as the number of explicitly included water molecules increases, the electronic embedding methods do not show any improvement until that number reaches five. The analysis of the snapshots in Figure 4 provides a first, but somewhat indirect, indication that electronic embedding approaches are more reliable than mechanical embedding for two key structural models: the embedded MoO 2− 4 , and the [MoO 4 (H 2 O) 6 ] 2− species. For a partial first solvation shell, it is nevertheless difficult to comprehend what is actually taking place from just following the σ iso values.

C. Real-space analysis of embedding
One can try to identify the underlying differences between embedding and supermolecular calculations by investigating the systems in real space. One option is to follow the differences in the electron density for the two treatments, as often done for other properties and previously done by Bulo and coworkers [24] for NMR shieldings. Figure 5 presents the difference in SR-ZORA density between the reference and the electronic embedding models ([A|11B f |9B f ] and [A|0|20B f ]). It is important to note here that this figure does not present density isosurfaces, but rather a volumetric description (i.e. the accumulated values ranging from lower (0.00) to the upper (0.01) bound considered) of the electron density. This representation has the advantage of enabling the visualization of high and low density regions in the same plot.
From Figure 5 we clearly see that inaccuracies in the embedding calculations, due to (for instance) the approximate treatment of the non-additive kinetic energy term, introduce small errors in the densities throughout the whole space. These errors tend to build up (green to red colors) in the regions away from the core of the MoO 2− 4 unit, and for the [A|11B f |9B f ] model they tend to be more significant around the second solvation shell, given that first solvation shell is fully relaxed in the presence of MoO 2− 4 . We also see from that figure that there is a small but non-negligible effect on the densities of the environment when we consider the molecules as a single subsystem or as a collection of fragments-as, in the latter, we create additional frontiers between subsystems, due to the monomer expansion and the eventual buildup of errors due to the approximate kinetic energy functionals.
Whatever the case, the analysis of the electron density is at odds with our findings for the mean values of shielding, since there appears to be nothing that points to significant errors in the region of the Mo-O bonds. This suggests that trying to understand the behavior σ iso values from changes in the electron density is not a suitable strategy.
An alternative to the visualization of the electron density is the analysis of the property densities. We have provided a first example for such an analysis for the analysis of embedding in the description of our 4-component implementation of subsystem embedding theory for magnetic properties, but for rather small systems [26]. Here we provide the analysis of a more complex example, while at the same time forsaking the use of isosurfaces in favor of volumetric plots.
Because the analysis of shielding densities is only implemented in the DIRAC code, we have restricted ourselves to the [MoO 4 (H 2 O) p ] 2− species (p = 1 − 6). These species are then considered as our references, and both mechanical and electronic embedding calculations considering only the molybdate ion in the active subsystem have been performed for each value of p. All fragments have therefore been treated with the DC Hamiltonian, with the embedding potentials being determined with freeze-thaw calculations using DIRAC. The results of these calculations are found in Figure 6, where we present the shielding density for the 95 Mo atom on the references, along with the differences in shielding density between the references and the electronic and mechanical embedding results. For the embedding calculations, contributions from the solvent water molecules at the lo-cation of the Mo atom are obtained with the NICS procedure.
The first striking feature of the embedding results is that, unlike for the electron densities, the error in the embedding calculations is indeed localized within the MoO 2− 4 species, in particular around the Mo atom. In addition to that, there are smaller errors also at the positions of the molybdate oxygen atoms.
We believe this points to a physical process that, while involving a rather local operator (the hyperfine operator has an effective r −2 dependence), makes the resulting property σ iso (described as the cross product of the hyperfine and the magnetically induced current density) have a much less localized nature, and with that nonnegligible contributions rather far from the responding atom (in the vicinity of the oxygen atoms) arise.
This goes to explain why, in both types of embedding, the dependence of σ iso on the number of water molecules is so significant: while the water molecules themselves do not contribute significantly to the shielding density, their absence (in the case of mechanical embedding) or the relatively inaccurate representation of the active and environment subsystems (in the case of electronic embedding) is enough to perturb the contribution to the shielding density around the molybdate oxygen atoms.
Furthermore, we see that as the first solvation shell is built up by including the nearest water molecules to the active subsystem (we show in Figure 6 only the even p values, see supplemental information for the odd values), the perturbation on the molybdate oxygen atoms is not accounted for in a systematic manner by the embedding methods, so that errors may build up even in regions in which water molecules had already been added. We do not yet possess the analysis tools to fully understand the interplay between these different effects, and are currently pursuing the development of new analysis approaches to address the issue.
Beyond these similarities, we nevertheless see that the errors for the mechanical embedding calculations, though of similar magnitude than those for the electronic embedding, extend much farther than for the latter, and explain why electronic embedding is more reliable.

V. CONCLUSIONS AND PERSPECTIVES
In this manuscript we provide an assessment of approaches based upon the frozen density embedding (FDE) framework for the creation of computationally efficient models capable of capturing solvent effects on NMR shieldings of heavy elements. We have investigated how water influences the shielding of 95 Mo in the molybdate dianion (MoO 2− 4 ), which can be at times a reference for 95 Mo NMR experiments, or a precursor for building Mo-containing materials.
A particular strength of FDE in this context is the ease with which different relativistic Hamiltonians can be used for different parts of the system. Here it allowed to com-pare the scalar relativistic (SR) and 2-component spinorbit (SO) ZORA Hamiltonians, and the 4-component Dirac-Coulomb (DC) Hamiltonian to describe the active subsystem while using the SR ZORA Hamiltonian for the aqueous environment.
From our calculations on different structural models, employing a large set of structures obtained from CPMD trajectories from work previously described in the literature, we have established that for the 95 Mo mean isotropic shielding value, there are weak but nonnegligible solvent shifts, due to the modification of the response of the electronic wavefunction of solute by the solvent. We find that even an approach along the lines of mechanical embedding yields a mean isotropic shielding value in good agreement with reference results obtained with standard DFT calculations on the full system.
Such mechanical embedding calculations, however, show a much larger spread of shielding values (for the different CPMD snapshots) with respect to the reference calculations than either of the electronic embedding variants considered. This leads us to conclude that electronic embedding is in fact an essential component for computational models.
By investigating the dependence of the isotropic shielding value from embedding calculations on the size of the active subsystems for selected CPMD snapshots, we have observed that convergence to the reference calculations it not monotonic and large discrepancies may still occur until roughly a full first solvation shell is explicitly included in the active subsystem [22].
We have found our results to be stable with respect to changes in the way the environment was described by FDE, that is, irrespective of whether we used embedding potentials obtained with the same (SR-ZORA or SO-ZORA) Hamiltonian throughout, or by using a combination of SR-ZORA based environment electron densities and electrostatic potentials together with a DC based calculation for the active subsystem.
This points to the possibility of devising efficient computational models for Mo-containing compounds, in which SR-ZORA DFT calculations are used to prepare embedding potentials for more sophisticated, 4component based calculations on the Mo-containing regions of interest. The latter approach is particularly interesting to investigate the absolute shieldings for molecules in complex environments.
We have found a significant shift (around 260 ppm) when changing from the SR-ZORA to the SO-ZORA Hamiltonian, and another significant shift (around 150 ppm) when changing from the SO-ZORA to the DC Hamiltonians. However, given the rather systematic nature of the differences observed, for chemical shifts of Mo complexes it is likely that the SR-ZORA Hamiltonian provides sufficient accuracy.
Finally, we have explored the visualization of the electron and shielding densities, as a means to provide further insight on the mechanisms behind the solvent effects, in addition to the performance (and shortcomings) of the embedding approaches. Our results underscore the point we have made previously, in that changes in shielding densities are better understood from the analysis of the differences of the corresponding shielding densities rather than from the changes in electron density.
Furthermore, we observe that the major source of er-rors in embedding calculations comes from the regions around the Mo nucleus, with less important discrepancies coming from all around the MoO 2− 4 species. As the major difference between the FDE and reference calculations lies in the use of approximate kinetic energy functionals to describe the respective non-additive contribu-tions, this suggests that for magnetic properties it would be of interest to try to improve the performance of such functionals. CPMD simulations on the molybdate ion.
The members of the PhLAM laboratory acknowledge support from the CaPPA project (

CONFLICT OF INTEREST
There are no conflicts to declare.

SUPPORTING INFORMATION AVAILABLE:
The supporting information is available at the publisher website.
Influence of the embedding model on 95 Mo shielding values The electron densities and electrostatic potential for the environment are obtained either from calculations on individual water molecules (and we refer to a "fragmented" environment, noted f ) or by these forming a single subsystem (and we refer to a "grouped" environment, noted g). We have looked at the influence of these two wait of treating the water molecules for two snapshots of the CPMD trajectories, snapshot 10 yielding results close to the average and snapshot 85 standing in the tail of the result distribution.
Table S1 reveals that the σ iso values obtained with FDE and FnT embeddings are marginally sensitive to the way the embedding density is calculated. For the sake of minimizing the computational cost, we will rely on "fragmented" water clusters for the generation of the FnT and FDE embedding potentials.

S2
Influence of the relativistic Hamiltonian on σ iso Figure S1: Evolution of the isotropic shielding σ iso (in ppm) as a function of the number of water molecules p in the active system ([A + pB|0|0]) or in the FnT process ([A|pB g |0]), with Dirac (top) and with ADF (bottom), for snapshot 10.
In the manuscript, we have noted that the change of relativistic Hamiltonian from SR-ZORA (ADF) to DC (DIRAC) translates into a change 410 ppm of the absolute shielding constant for the [A|11B f |9B f ] model (See Table 3