\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}







\def\ISO{[A|0|0]}
\def\DIRAC{[A|11B_f|9B_f]^{216}_{\mathcal{D}irac}}
\def\SUP{[A+20B|0|0]}
\def\FnT{[A|11B_f|9B_f]}
\def\FDE{[A|0|20B_f]}

\iflatexml
\usepackage[separate-uncertainty=true]{siunitx}
\usepackage[space]{grffile}
\usepackage{longtable}
\usepackage{booktabs,array,multirow}

\else

\corraddress{Andr\'e Severo Pereira Gomes, Universit\'e de Lille, CNRS, UMR 8523 -- PhLAM -- Physique des Lasers, Atomes et Mol\'ecules, F-59000 Lille, France}
\corremail{andre.gomes@univ-lille.fr}
\fundinginfo{
Polish National Science Centre (2016/23/D/ST4/03217), 
Labex CaPPA (ANR-11-LABX-0005-01) and I-SITE ULNE project OVERSEE (ANR-16-IDEX-0004), 
CPER CLIMIBIO (European Regional Development Fund, Hauts de France council, French  Ministry of Higher
Education and Research), 
CNRS Institute of Physics (INP) via the PICS program (grant 6386),
French national supercomputing facilities (DARI A0030801859, A0050801859)
}
\fi


\begin{document}

\title{Investigating solvent effects on the magnetic properties of molybdate
ions (MoO\textsubscript{4}\textsuperscript{2-}) with relativistic
embedding}



\author[1]{Loïc Halbert}%
\author[2]{Malgorzata Olejniczak}%
\author[1]{Valerie Vallet}%
\author[1]{Andre Severo Pereira Gomes}%
\affil[1]{Université de Lille}%
\affil[2]{University of Warsaw}%


\vspace{-1em}



  
  \date{May 19, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
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.%
\end{abstract}%



\sloppy


\section*{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 so-called NMR shielding constant, $\sigma^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 ($\vec{m}_K$) and an external magnetic field ($\vec{B}$) at the zero-perturbation limit:

\begin{equation}
\sigma^K_{\alpha\beta} =  \frac{d^2E}{dB_\alpha dm_{K;\beta}} \biggr|_{\vec{B},\{\vec{m}_A\} = 0}.
\label{eq:shield1} 

\end{equation}
In most cases, what is measured in experiments is not $\sigma^K$ but rather a signal relative to a chosen reference species, the chemical shift $\delta^K$, though in recent years experimental advances have revived the interest in determining absolute shielding scales~\hyperref[csl:1]{(Aucar \& Aucar, 2019}; \hyperref[csl:2]{Jackowski \& Makulski, 2011}; \hyperref[csl:3]{\textit{{The Determination of Accurate Nuclear Magnetic Dipole Moments and Direct Measurement of Nmr Shielding Constants}}, 2012)}.

NMR spectra can be quite difficult to interpret. The difficulties may come from the complexity due to the size 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.~\hyperref[csl:4]{(\textit{{Solid-State Nmr of Quadrupolar Nuclei for Investigations Into Supported Organometallic Catalysts: Scope and Frontiers}}, 2572)} Because of these difficulties, theoretical modeling has become indispensable for the interpretation of experiments~(missing citation); \hyperref[csl:6]{Charpentier, 2011}; (missing citation). 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~\hyperref[csl:8]{(Gomes \& Jacob, 2012}; \hyperref[csl:9]{Jacob \& Neugebauer, 2014}; \hyperref[csl:10]{Tomasz~Adam Wesolowski \& Zhou, 5891}; \hyperref[csl:11]{Sun \& Chan, 2705)} 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$~\hyperref[csl:12]{(\textit{{Evaluation of $\mathrm ^ 95 Mo Nuclear Shielding and Chemical Shift of \mathrm [mo_6x_ 14 ]^ 2 $ Clusters in the Liquid Phase}}, 7673)}. 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~\hyperref[csl:13]{(Schrock, 5512}; \hyperref[csl:14]{Jang \& Goddard~III, 2001}; (missing citation), as fission products~\hyperref[csl:16]{(\textit{{High-Temperature Behavior of Dicesium Molybdate Cs2moo4: Implications for Fast Neutron Reactors}}, 2014)}, and also make up materials that can be used for applications ranging from nonlinear optical materials~\hyperref[csl:17]{(\textit{{Synthesis and Tunable Nonlinear Absorption Properties of Zn3mo2o9 Nanosheet Ceramic Material}}, 1095)} to glasses for vitrification of nuclear waste~\hyperref[csl:18]{(\textit{{Molybdenum in Glasses Containing Vitrified Nuclear Waste}}, 2005)}. 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)~\hyperref[csl:6]{(Charpentier, 2011}; (missing citation) 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~(missing citation); \hyperref[csl:20]{Cui \& Karplus, 3721}; \hyperref[csl:21]{Andrea Frank \& Exner, 1480)}, 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)~\hyperref[csl:22]{(\textit{{Structural Interpretation of the 31 P Nmr Chemical Shifts in Thiophosphate and Phosphate: Key Effects Due to Spin-Orbit and Explicit Solvent}}, 9924)}, 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~\hyperref[csl:23]{(Jacob \& Visscher, 1941}; \hyperref[csl:24]{Rosa~E. Bulo \& Visscher, 2640}; \hyperref[csl:25]{Andreas~W G{\"{o}}tz \& Visscher, 1041}; \hyperref[csl:26]{Ma{\l}gorzata Olejniczak \& Gomes, 8400)}. 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~\hyperref[csl:27]{(Pyykk{\"o}, 1988}; \hyperref[csl:28]{Dyall \& Knut~Faegri, 2007}; \hyperref[csl:29]{Reiher \& Wolf, 2009}; \hyperref[csl:30]{Autschbach, 1509)} (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~\hyperref[csl:31]{(Martin Kaupp \& Malkin, 2004}; \hyperref[csl:32]{Autschbach, 2011)}.

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~\hyperref[csl:33]{(Bartlett \& Musia{\l}, 2007}; \hyperref[csl:34]{\textit{{Recent Advances in Wave Function-Based Methods of Molecular-Property Calculations}}, 2012)}, thus paving the way for a quantitative leap in the interpretation of NMR spectra. As the hydrated MoO$_4^{2-}$ molybdate dianion appears to have well-structured hydrations shells~\hyperref[csl:12]{(\textit{{Evaluation of $\mathrm ^ 95 Mo Nuclear Shielding and Chemical Shift of \mathrm [mo_6x_ 14 ]^ 2 $ Clusters in the Liquid Phase}}, 7673)}, 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~\ref{sec:theory} we outline the key theoretical aspects underpinning the work. This is followed in section~\ref{computational-details} 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~\ref{results} 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~\ref{conclusions}.

\section*{Theory}
\label{sec:theory}

In what follows, we employ the standard notation for molecular orbital indices (with $i, j, \ldots$ corresponding to occupied, $a, b, \ldots$ - to virtual and $p, q, \ldots$ - to general molecular orbitals) and the summation convention over repeated indices. We use the SI-based atomic units (~$\hbar = m_e = e =1/(4\pi\epsilon_0) = 1$)~\hyperref[csl:35]{(Whiffen, 1978)}. Employing the second-quantization formalism and an exponential parametrization of the closed-shell ground state with the orbital-rotation operator, $\hat{\kappa} = \sum_{pq} \kappa_{pq}p^\dagger q$,

\begin{equation}
\psi_p = \psi_q \exp(-\kappa_{pq}),
\label{eq:exppar} 

\end{equation}
the derivative in Eq.~\ref{eq:shield1} can be further expressed as

\begin{equation}
\frac{d^2E}{dB_\alpha dm_{K;\beta}} \biggr|_{\varepsilon = 0} = 
\frac{\partial^2 E }{ \partial \kappa_{pq} \partial m_{K;\beta}}  \frac{\partial \kappa_{pq}}{\partial B_\alpha} \biggr|_{\varepsilon = 0}
+
\frac{\partial^2E}{\partial B_\alpha \partial m_{K;\beta}} \biggr|_{\varepsilon = 0},
\label{eq:shield2} 

\end{equation}
where all perturbing fields are collected in $\varepsilon$.
Partial derivatives of the orbital rotation parameters ($\{\kappa_{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:

\begin{equation}
E_0^{[2]} X_{\varepsilon_1} = - E_{\varepsilon_1}^{[1]},
\label{eq:lr} 

\end{equation}
with $E_0^{[2]}$, $E_{\varepsilon_1}^{[1]}$ and $X_{\varepsilon_1}$ corresponding to the electronic Hessian, the property gradient and the solution vector ($X_{\varepsilon_1} = ( \partial \kappa_{pq} / \partial {\varepsilon_1})_{\varepsilon = 0}$), respectively~\hyperref[csl:36]{(Saue \& Jensen, 2003)}.
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
~\hyperref[csl:37]{(Saue, 3077)} and on the spin-density functional theory (SDFT)~\hyperref[csl:38]{(Barth \& L.~Hedin, 1629}; \hyperref[csl:39]{Rajagopal \& J.~Callaway, 1912}; \hyperref[csl:40]{Jacob \& Reiher, 3661)}. 
The closed-shell system is described by the electron density ($\rho_0$) and the spin-density calculated in a non-collinear fashion as a norm of the spin-magnetization vector ($s = |\rho_\mu|$, $\mu \in \{x,y,z\}$). Following previous works~\hyperref[csl:41]{(Radovan Bast \& Saue, 2091}; \hyperref[csl:42]{\textit{{A Fully Relativistic Method for Calculation of Nuclear Magnetic Shielding Tensors with a Restricted Magnetically Balanced Basis in the Framework of the Matrix Dirac--Kohn--Sham Equation}}, 1041}; \hyperref[csl:43]{Ma\l{}gorzata Olejniczak \& Pecul, 141AD}; \hyperref[csl:26]{Ma{\l}gorzata Olejniczak \& Gomes, 8400)}, we collect the electron density and the spin-magnetization vector into one variable - a general density component, $\rho_k$ ($k \in \{0,x,y,z\}$).

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~\hyperref[csl:8]{(Gomes \& Jacob, 2012}; \hyperref[csl:44]{Goez \& Neugebauer, 2018}; \hyperref[csl:9]{Jacob \& Neugebauer, 2014)},
which is realized by expressing the electron density of the full system as a sum of electron densities of subsystems~\hyperref[csl:45]{(Wesolowski \& Warshel, 8050}; \hyperref[csl:46]{Wesolowski \& Weber, 1996}; \hyperref[csl:47]{Wesolowski \& Weber, 2002}; \hyperref[csl:10]{Tomasz~Adam Wesolowski \& Zhou, 5891)} (the same applies to the spin density if it is included in the formalism),

\begin{equation}
\rho_k^\text{tot} (\vec{r}) = \rho_k^I (\vec{r}) +\rho_k^{II} (\vec{r}) 
\label{eq:FDE_sum_density} 

\end{equation}
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,

\begin{equation}
E_\text{tot}[\rho_k^\text{tot}] = E_I[\rho_k^I] + E_{II}[\rho_k^{II}] + E_\text{int}[\rho^{I}_k, \rho^{II}_k]
\label{eq:FDE_sum_energy} 

\end{equation}
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:

\begin{align}
    E_\text{int}[\rho^{I}_k, \rho^{II}_k]
      &=  \int \left[ \rho_0^{I}(\vec{r}) v_\text{nuc}^{II}(\vec{r}) + \rho_0^{II}(\vec{r}) v_\text{nuc}^{I}(\vec{r}) \right]d\vec{r}  
      + \iint \frac{\rho_0^{I}(\vec{r}_1) \rho_0^{II}(\vec{r}_2)}{|\vec{r}_1 - \vec{r}_2|} d\vec{r}_1d\vec{r}_2
      + E_\text{nuc}^{I,II}      
      + E_\text{xc}^\text{nadd} + T_{s}^\text{nadd},
      \label{eq:Eint} 

\end{align}
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.~\ref{eq:Eint}), the non-additive exchange-correlation ($E_\text{xc}^\text{nadd}$) and kinetic energy ($T_{s}^\text{nadd}$) contributions emerge. These non-additive contributions~\hyperref[csl:48]{(Sebastian H{\"o}fener \& Visscher, 441AD)}  are defined as:

\begin{equation}
    X^\text{nadd} \equiv X^\text{nadd}[\rho_k^{I}, \rho_k^{II}] = X[\rho_k^\text{tot}] - X[\rho_k^{I}] - X[\rho_k^{II}]
    \label{eq:nonadd} 
  
\end{equation}
and they depend on the general density components of all subsystems.  
In typical FDE calculations, one subsystem is chosen as \emph{active} , while another subsystem constitute its \emph{environment} . In this setting, the \emph{active} density (e.g. $\rho_k^I$) is optimized by solving the so-called Kohn-Sham (KS) equations for a constrained electron density (KSCED)~\hyperref[csl:49]{(Wesolowski \& Warshel, 8050)}, 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 \emph{active} subsystem with the \emph{environment} (here: $\rho_k^{II}$):

\begin{equation}
\label{eq:embpot} 
    v_\text{emb;k}^I(\vec{r}) = \frac{\delta E_\text{int}}{\delta \rho_k^I(\vec{r})} 
                   = \frac{\delta E_\text{xc}^\text{nadd}}{\delta \rho_k^I} +  \frac{\delta T_{s}^\text{nadd}}{\delta \rho_k^I}  + v_\text{nuc}^{II}(\vec{r}) + \int \frac{\rho_0^{II}(\vec{r}')}{|\vec{r} - \vec{r}'|}d\vec{r}' \nonumber.
\end{equation}
It is also possible to relax the density of the \emph{environment} by solving the analogous KSCED equations set up for $\rho_k^{II}$. Repeating this procedure in an iterative manner, known as the \emph{freeze-thaw} ~\hyperref[csl:50]{(Christoph~R. Jacob \& Visscher, 1011)} (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~\hyperref[csl:8]{(Gomes \& Jacob, 2012}; \hyperref[csl:51]{Jacob \& Neugebauer, 2014}; \hyperref[csl:52]{\textit{{Subsystem-Dft Potential-Energy Curves for Weakly Interacting Systems}}, 1432}; \hyperref[csl:53]{\textit{{Accurate Frozen-Density Embedding Potentials As a First Step Towards a Subsystem Description of Covalent Bonds}}, 1641}; \hyperref[csl:54]{Denis~G. Artiukhin \& Neugebauer, 2341)}.
From the presence of the interaction energy term in the total energy expression (Eq.~\ref{eq:FDE_sum_density}) 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_0^{[2]}$, $E_{\varepsilon_1}^{[1]}$ in Eq.~\ref{eq:lr}) 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.~\ref{eq:embpot}) and the embedding kernel (see~\hyperref[csl:26]{(Ma{\l}gorzata Olejniczak \& Gomes, 8400)} for details). In effect, the property gradient and the electronic Hessian (Eq.~\ref{eq:lr}) are elegantly subdivided into contributions from isolated subsystems and terms related to the interaction input~\hyperref[csl:48]{(Sebastian H{\"o}fener \& Visscher, 441AD)}. In practice, these derivatives can be coded analogously to the exchange-correlation 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 system - comprising of all subsystems - to certain external perturbations. In particular, for the NMR shielding tensor of a nucleus K this can be symbolically expressed as

\begin{equation}
\label{eq:shield_fde} 
\sigma_{\alpha\beta}^K = \sigma_{\alpha\beta}^{K, I} + \sigma_{\alpha\beta}^{K, II}.
\end{equation}

Assuming that the nucleus K belongs to subsystem I, the term $\sigma_{\alpha\beta}^{K, II}$ can be interpreted as an effect that the perturbation due to $B_\alpha$ 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)~\hyperref[csl:55]{(\textit{{Nucleus-Independent Chemical Shifts: A Simple and Efficient Aromaticity Probe}}, 6317)}.
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)~\hyperref[csl:56]{(F.~London, 1937}; \hyperref[csl:57]{Trygve Helgaker \& Ruud, 1999)} instead of atomic orbitals, however with the price of the need to calculate additional terms in the property gradient~\hyperref[csl:58]{(Miroslav Ilia\v{s} \& Jensen, 1241}; \hyperref[csl:43]{Ma\l{}gorzata Olejniczak \& Pecul, 141AD}; \hyperref[csl:26]{Ma{\l}gorzata Olejniczak \& Gomes, 8400)}.

The NMR shielding tensor (Eqs.~\ref{eq:shield1},~\ref{eq:shield2}) can also be determined from numerical integration of the NMR shielding density~\hyperref[csl:59]{(Jameson \& Buckingham, 5684}; \hyperref[csl:60]{Jameson \& AD~Buckingham, 3366)}. 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~\hyperref[csl:61]{(Radovan Bast \& Saue, 2009}; \hyperref[csl:62]{David Sulzer \& Saue, 2068)} reads

\begin{equation}
\label{eq:jBj} 
j^B = \left. \frac{d\vec{j}}{d\vec{B}} \right|_{\vec{B} = 0}; 
\quad
\vec{j} = -e \psi_i^\dagger  c \vec{\alpha} \psi_i.
\end{equation}

The NMR shielding density is an integrand in the following expression:

\begin{equation}
\sigma^K_{\alpha\beta} = -\frac{1}{c^2} \int \frac{1}{r_K^3} \left( \vec{r}_K \times \vec{j}^{\,B_\alpha} \right)_\beta d\vec{r}.
\label{eq:shield3} 

\end{equation}
This reformulation of $\sigma^K_{\alpha\beta}$ 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~\hyperref[csl:26]{(Ma{\l}gorzata Olejniczak \& Gomes, 8400)}. 
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.~\ref{eq:shield_fde}~\hyperref[csl:23]{(Jacob \& Visscher, 1941)}.

\section*{Computational details\label{computational-details}}

The solvated molybdate structures used in this work are taken from the Car-Parrinello molecular dynamics (CPMD~\hyperref[csl:63]{(Hutter, 2011)}) calculations of Nguyen and coworkers~\hyperref[csl:12]{(\textit{{Evaluation of $\mathrm ^ 95 Mo Nuclear Shielding and Chemical Shift of \mathrm [mo_6x_ 14 ]^ 2 $ Clusters in the Liquid Phase}}, 7673)}. 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:

\begin{itemize}
\item a supermolecular system, containing the molybdate ion and all water molecules;

\item 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 \emph{freeze-thaw} 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);

\item 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;

\item 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''.

\end{itemize}

Using the labels ``A'' for the molybdate ion and ``B'' for the water molecules, and indices \textbf{g} and \textbf{f} to denote the ``\textbf{g} rouped'' or ``\textbf{f} ragmented'' environment, we shall use the following shorthand notation to represent these four cases: $\left[Active\,|\, Relaxed  \, |\, Frozen \right]$. This means, for example, that the supermolecular model is denoted by $\left[A + 20 B | 0 | 0 \right]$ (all species are in the active fragment), the mechanical embedding by $\left[A  | 0 | 0 \right]$, 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 $\left[A | 11 B_f | 9 B_f \right]$.

The execution of all calculations, including the preparation of the different fragments, has been handled with the PyADF~\hyperref[csl:64]{(\textit{{Pyadf  a Scripting Framework for Multiscale Quantum Chemistry}}, 2328)} scripting framework.
In all calculations, we employ Becke integration grids, and a gaussian nuclear charge distribution~\hyperref[csl:65]{(Visscher \& Dyall, 1997)} to describe the finite nuclear volume~\hyperref[csl:66]{(Andrae, 2000)}, PBE~\hyperref[csl:67]{(John~P. Perdew \& Ernzerhof, 3865}; \hyperref[csl:68]{John~P. Perdew \& Ernzerhof, 1396)} exchange-correlation functional. For subsystem embedding and FDE embedding calculations (see below), we have employed the PW91k~\hyperref[csl:69]{(A.~Lembarki \& H.~Chermette, 5328)} and PBE~\hyperref[csl:67]{(John~P. Perdew \& Ernzerhof, 3865}; \hyperref[csl:68]{John~P. Perdew \& Ernzerhof, 1396)} 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 \emph{freeze-thaw} cycles. In the NMR calculations we employ London orbitals.

\subsection*{One- and two-component calculations}

The scalar relativistic (SR) and spin-orbit (SO) ZORA~\hyperref[csl:70]{(E.~{Van Lenthe} \& Baerends, 6505}; \hyperref[csl:71]{E.~{Van Lenthe} \& Snijders, 9783}; \hyperref[csl:72]{E.~{Van Lenthe} \& Snijders, 1996}; \hyperref[csl:73]{{Van Lenthe}, 8943)} calculations have been performed with the ADF~\hyperref[csl:74]{(software for chemistry \& materials, n.d.)} code, using the TZ2P basis sets~\hyperref[csl:75]{(E.~{Van Lenthe} \& Baerends, 1142)}. The NMR calculations have been performed with the NMR program.

\subsection*{Four-component calculations}

The calculations employing the Dirac-Coulomb (DC) Hamiltonian have been performed with the Dirac code~\hyperref[csl:76]{({DIRAC}, 2018)}, revision \texttt{54ab939}. For the Mo atom the dyall.cv3z basis~\hyperref[csl:77]{(Dyall, 2007)} is used, and for all lighter atoms the aug-cc-pVTZ~\hyperref[csl:78]{(Jr, 1007)} 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.~\hyperref[csl:26]{(Ma{\l}gorzata Olejniczak \& Gomes, 8400)}. 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~\ref{analysis-shielding}) have been obtained with the analysis module of the Dirac code, and are based on four-component calculations with the DC Hamiltonian and equation~\ref{eq:shield3}. In the case of embedded calculations, following equation~\ref{eq:shield_fde}, it was necessary to also perform calculations for the environment using the position of the Mo atom as the $\vec{r}_K$ position.

\subsection*{Plotting}

All graphs have been prepared with Matplotlib~\hyperref[csl:79]{(Hunter, 2007)} python library, with the exception of the volumetric density and NMR shielding density plots, for which the  Mayavi~\hyperref[csl:80]{(Ramachandran \& Varoquaux, 2011)} python library has been used.

\section*{Results and discussion\label{results}}

In what follows we are interested in the isotropic part of the NMR shielding tensor of the $^{95}$Mo nucleus, $\sigma_{iso}$, which is 
defined as:

\begin{eqnarray}
\sigma_{iso}=1/3(\sigma_{11}+\sigma_{22}+\sigma_ {33})

\end{eqnarray}
where the indices refer to the principal axes of the tensor according to the Mason's notation~\hyperref[csl:81]{(J.~Mason, 1993)}, with $\sigma_{11}\leq \sigma_{22} \leq \sigma_{33} $.

\subsection*{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$_4^{2-}$ ion in the active subsystem.  In this subsection we explore this model.

\subsubsection*{Structural and electronic effects on NMR shieldings}

We start the discussion with the analysis of the influence of the geometry of molecular complexes on $\sigma_{iso}$ based on the SR-ZORA calculations. This is summarized in Figure~\ref{fig:Evolution_Shielding_snap400_500}, which shows the values of $\sigma_{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 $\sigma_{iso}$ values, for others these contributions are negligible and the mechanical embedding is sufficient to reproduce the supermolecular baseline.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-pyadf-all-ZORA-snap-400-500/fig-pyadf-all-ZORA-snap-400-500}
\caption{{\label{fig:Evolution_Shielding_snap400_500} $\sigma_{iso}$ (in ppm) for a subset of the snapshots considered (snapshots 400 to 500), obtained from SR-ZORA calculations for the
   supermolecule ($\SUP$), electronic (FnT:$\FnT$ and FDE:$\FDE$) and mechanical ($\ISO$) embedding.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{table}[h]
\centering
\begin{tabular}{
        ll
}
\hline
{Model} & {$\bar{\sigma}_{iso} \pm std$} \\
\hline
{$ \SUP ^ {216} $} & $-845.7 \pm 45.4$ \\ 
{$ \FnT ^ {216} $} & $-839.1 \pm 45.2$ \\
{$ \DIRAC $} & $-428.3 \pm 45.5$ \\
\hline
\end{tabular}
\caption{{Mean values for the isotropic NMR shielding constant ($\overline{\sigma}_{iso}$, in ppm) for the different computational models, and the
 respective standard deviations ($std$), obtained from a subset of 216 CPMD snapshots chosen following the analysis of the normal probability plot (see supplemental information for details). The 4-component embedding calculations ({$\DIRAC$}) are performed using embedding potentials and frozen densities obtained from the SR-ZORA $\FnT^{216}$ model.}}
\label{tab:table_all_results_Dirac}

\end{table}

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 $\bar{\sigma}_{iso}$ in Table~\ref{tab:table_all_results_ZORASR}. However, the choice of the computational model does affect the standard deviation of the $\sigma_{iso}$ values distribution - in this case the variance of the NMR shielding values over all the snapshots is significantly smaller 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  $\sigma_{iso}$.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-ZORA-SR-vs-meanDistance-dataZORASRiso/fig-ZORA-SR-vs-meanDistance-dataZORASRiso}
\caption{{\label{fig:ZORASR_vs_meandist} Correlations between $\sigma_{iso}$ (in ppm) and the mean distance between Mo and the oxygen atoms in the Molybdate ion ($\bar{r}$[Mo-O], in \r{A}), obtained from SR-ZORA calculations for various models: (a) isolated $\ISO$; (b) FDE $\FDE$; (c) FnT $\FnT$; and (d) Supermolecule $\SUP$ for all snapshots.%
}}
\end{center}
\end{figure}

In order to understand these structural dependencies of $\sigma_{iso}$ values, we decided to define a measure, the mean distance between the centers of Mo atom and O atoms in MoO$_4^{2-}$, hereafter referred to as $\bar{r}$[Mo-O], and to plot the $\sigma_{iso}$ vs $\bar{r}$[Mo-O] in Figure~\ref{fig:ZORASR_vs_meandist}.

The first conclusion of this figure is that there is a significant correlation between the $\sigma_{iso}$  and  $\bar{r}$[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 $\sigma_{iso}$ values cannot be simply attributed to the structural variations, or at least the ones represented by the $\bar{r}$[Mo-O] measure. The FDE and FnT methods recover some of the $\bar{r}$[Mo-O]-dependent contributions, as illustrated by the correlation coefficients ($R^2 = 0.962$ and $0.976$ for FnT and FDE, respectively). The slope of the linear regression also shows a better agreement between supermolecule and electronic embedding calculations than between supermolecule and mechanical embedding.

However, the linearity of the relationship between   $\sigma_{iso}$  and  $\bar{r}$[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 $\bar{r}$[Mo-O] and the values of the differences of $\sigma_{iso}$ obtained from embedding calculations and the supermolecular ones, denoted by $\Delta$ in Figure~\ref{fig:Error_ZORASR_vs_meandist}. From these results, we see that there is in fact no correlation between the two variables ($R^2 \leq 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 $\Delta$ 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 the $\bar{r}$[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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-Error-SR-vs-meanDistance-dataZORASRiso/fig-Error-SR-vs-meanDistance-dataZORASRiso}
\caption{{\label{fig:Error_ZORASR_vs_meandist} Difference between the embedded models (a) isolated $\ISO$; (b) FDE $\FDE$; and (c) FnT $\FnT$, and the supermolecule SR-ZORA results ($\Delta$, in ppm) as a function of the mean distance between Mo and the oxygen atoms in the Molybdate ion ($\bar{r}$[Mo-O]).%
}}
\end{center}
\end{figure}

\subsubsection*{Assessment of spin-orbit Hamiltonians}

In order to assess the importance of spin-orbit interactions on the the solvent effects, we begin by considering the SO-ZORA Hamiltonian. The 
$\bar{\sigma}_{iso}$ values and standard deviations for supermolecular and embedded calculations are shown in Table~\ref{tab:table_all_results_ZORASO}.

From these we observe that there are no major differences with respect to the SR-ZORA results (see  Table~\ref{tab:table_all_results_ZORASR}) 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 $\bar{\sigma}_{iso}$ values than the mechanical embedding case, though once more the FDE $\bar{\sigma}_{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.

Comparing the results in Tables~\ref{tab:table_all_results_ZORASR} and~\ref{tab:table_all_results_ZORASO} we also observe significant changes on the absolute $\bar{\sigma}_{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~\hyperref[csl:82]{(\textit{{Relativistic Heavy-Atom Effects on Heavy-Atom Nuclear Shieldings}}, 2006}; \hyperref[csl:83]{Micha{\l} Jaszu{\'{n}}ski \& Ruud, 2017)}, and therefore we refer the readers to the review~\hyperref[csl:83]{(Micha{\l} Jaszu{\'{n}}ski \& Ruud, 2017)} and recent examples.~\hyperref[csl:82]{(\textit{{Relativistic Heavy-Atom Effects on Heavy-Atom Nuclear Shieldings}}, 2006}; \hyperref[csl:84]{Alkan \& Dybowski, 2018}; \hyperref[csl:85]{Y.~Ootani \& H.~Fukui, 2006}; \hyperref[csl:86]{Valentin~A. Semenov \& Krivdin, 4908}; \hyperref[csl:87]{Holmes \& Schurko, 1785)}

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$_4^{2-}$ 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.

The exploration of a large number of snapshots also provided data with which to guide us in reducing the number of CPMD snapshots used for the more expensive calculations employing the DC Hamiltonian. With the help of normal probability plots (see supplemental information), we have been able to use the SR-ZORA data to reduce the number of snapshots to 216, roughly half of the original sample space, while at the same time assuring to keep a normal distribution of $\sigma_{iso}$ values for the SR-ZORA (supermolecule and embedded) calculations. We subsequently used this reduced sample to perform the DC calculations, for which results are shown in Table~\ref{tab:table_all_results_Dirac}.

While we cannot directly compare the standard deviations from the DC results to those of tables~\ref{tab:table_all_results_ZORASR} and~\ref{tab:table_all_results_ZORASO} 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 $\bar{\sigma}_{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.

\subsection*{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 explicitly 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 $\overline{\sigma}_{iso}$ value and towards the  tails of the distribution. Figure~\ref{graph:snap10_85_influence_N_eau} 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 (${\sigma}_{iso}$ closer to $\overline{\sigma}_{iso}$) we observe significant differences between embedding methods, with the subsystem embedding model ([$A +p B | (11-p)B_f | 9 B_f$ ]) yielding results which are only a few ppm below the reference supermolecular calculations. The FDE model ([$A +p B | 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 ${\sigma}_{iso}$. At around six explicit water molecules, which roughly corresponds to the first solvation shell around the molybdate ion~\hyperref[csl:12]{(\textit{{Evaluation of $\mathrm ^ 95 Mo Nuclear Shielding and Chemical Shift of \mathrm [mo_6x_ 14 ]^ 2 $ Clusters in the Liquid Phase}}, 7673)}, 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-ZORA-SR-activeH2O-snap10-snap85/fig-ZORA-SR-activeH2O-snap10-snap85}
\caption{{\label{graph:snap10_85_influence_N_eau} Variations of $\sigma_{iso}$ (in ppm) with increasing number of water molecules in the active subsystem for selected structures (top: snapshot 10; bottom: snapshot 85). All values are computed with the SR-ZORA Hamiltonian.%
}}
\end{center}
\end{figure}

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 (${\sigma}_{iso}$ away from $\overline{\sigma}_{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 ${\sigma}_{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~\ref{graph:snap10_85_influence_N_eau} 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$_4^{2-}$, 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 ${\sigma}_{iso}$ values.

\subsection*{Real-space analysis of embedding\label{analysis-shielding}}

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~\hyperref[csl:24]{(Rosa~E. Bulo \& Visscher, 2640)} for NMR shieldings.

Figure~\ref{graph:densite_2_B} presents the difference in SR-ZORA density between the reference and the electronic embedding models ($\FnT$ and $\FDE$). 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-Snap10-Diff-dens-FnT-Supermol-cubegrid-0p01/fig-Snap10-Diff-dens-FnT-Supermol-cubegrid-0p01}
\caption{{\label{graph:densite_2_B} Volumetric plot of the difference of SR-ZORA electron densities between (a) $[A|11B_g|9B_f]$ and $\SUP$; and (b) $\FnT$ and $[A|11B_g|9B_f]$ for snapshot 10. Displayed are difference density values at each grid point between 0.00 and 0.01 a.u.%
}}
\end{center}
\end{figure}

From Figure~\ref{graph:densite_2_B} 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$_4^{2-}$ unit, and for the $\FnT$ 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$_4^{2-}$. 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 ${\sigma}_{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~\hyperref[csl:26]{(Ma{\l}gorzata Olejniczak \& Gomes, 8400)}. 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 \emph{freeze-thaw} calculations using DIRAC. The results of these calculations are found in  Figure~\ref{graph:densite_shielding_SUP}, 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 location 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$_4^{2-}$ 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 ${\sigma}_{iso}$ (described as the cross product of the hyperfine and the magnetically induced current density) have a much less localized nature, and with that non-negligible 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 ${\sigma}_{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~\ref{graph:densite_shielding_SUP} 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-MoO4m2-2H2O-snap10-Net-vis-Shield-MoO4-2H2O-iso-oMoH-0p1/fig-MoO4m2-2H2O-snap10-Net-vis-Shield-MoO4-2H2O-iso-oMoH-0p1}
\caption{{\label{graph:densite_shielding_SUP} Volumetric plots of shielding densities and differences in shielding densities between 0.0 and 0.1 ppm, computed with the DC Hamiltonian, for [MoO$_4$(H$_2$O)$_p$]$^{2-}$ with increasing number of $p$ water molecules. Left column: Supermolecule ($[A+pB|0|0]$); Middle column: difference between supermolecular and electronic embedding ($[A|6B_g|0]$) calculations; Right column:difference between supermolecular and mechanical embedding ($[A|0|0]+[0|6B_g|0]$) calculations. For embedded models, water molecule contributions to $^{95}$Mo shielding densities are calculated with the NICS method.%
}}
\end{center}
\end{figure}

\section*{Conclusions and perspectives\label{conclusions}}

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$_4^{2-}$),  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 compare the scalar relativistic (SR) and 2-component spin-orbit (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 non-negligible 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~\hyperref[csl:22]{(\textit{{Structural Interpretation of the 31 P Nmr Chemical Shifts in Thiophosphate and Phosphate: Key Effects Due to Spin-Orbit and Explicit Solvent}}, 9924)}.

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, 4-component 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 errors in embedding calculations comes from the regions around the Mo nucleus, with less important discrepancies coming from all around the MoO$_4^{2-}$ 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 contributions, this suggests that for magnetic properties it would be of interest to try to improve the performance of such functionals.

\section*{Acknowledgements}

The authors thank Dr.\ J\selectlanguage{ngerman}érôme Cuny (Laboratoire de Chimie et Physique Quantiques, UMR 5626, Toulouse) for kindly providing us with the snapshots from the CPMD simulations on the molybdate ion.

\section*{Conflict of interest}
There are no conflicts to declare.

\section*{Supporting Information Available:}

The supporting information is available at the publisher website.

\appendix

\clearpage

\section*{Supplemental information}

\subsection*{Influence of the embedding model on $^{95}$Mo shielding values}\selectlanguage{english}
\begin{table}[htbp]
\centering
\begin{tabular}{
        lccccc
}
\toprule
& \multicolumn{2}{c}{snapshot 10}&& \multicolumn{2}{c}{snapshot 85}\\
\cline{2-3}\cline{5-6}
{Model} & {$\sigma_{iso}$} & {$\Delta\sigma_{iso}^{ref}$} && {$\sigma_{iso}$} & {$\Delta\sigma_{iso}^{ref}$} \\ 
\midrule
{$\SUP$ [$ref$]} & -805.141 & && -1054.585 & \\
{$\FnT$}& -808.887 & -3.746 && -1034.507 & 20.078\\
{$[A|11B_g|9B_f]$}& -808.879 & -3.738&& -1035.921 & 18.664\\
{$\FDE$}& -822.851 & -17.710 && -1028.881 & 25.704\\
{$[A|0|11B_g+9B_f]$}& -822.446 & -17.305&& -1032.424 & 22.161\\  
\bottomrule
\end{tabular}
\caption{{SR-ZORA values of $\sigma_{iso}$ as a function of the fragmentation of the first solvation layer of MoO$_4^{2-}$, and $\Delta\sigma_{iso}^{ref}$, with respect to the supermolecule $\SUP$ reference, for snapshots 10 and 85.}}
\label{tab:val_snap10_11grouped}

\end{table}

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~\ref{tab:val_snap10_11grouped} reveals that the $\sigma_{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.

\clearpage

\subsection*{Influence of the relativistic Hamiltonian on $\sigma_{iso}$}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-Evol-shield-pH2O-snap10-essai1/fig-Evol-shield-pH2O-snap10-essai1}
\caption{{\label{graph:All_snap10_bis} Evolution of the isotropic shielding $ \sigma_{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.%
}}
\end{center}
\end{figure}

In the manuscript, we have noted that the change of relativistic Hamiltonian from SR-ZORA (ADF) to DC (DIRAC) translates into a change \SI{410}{ppm} of the absolute shielding constant for the $\FnT$ model (See Table~3). To verify whether this shift is independent of the chemical model, we have performed a systematic comparison of $\sigma_{iso}$ values obtained with SR-ZORA and DC Hamiltonians for molybdate [MoO$_4$(H2O)$_p$]$^{2-}$ hydrated by $p=1,6$ water molecules either treated explicitly ($[A+pB|0|0]_{ADF}$ and $[A+pB|0|0]_{Dirac}$), or by an embedding potential treating the (H2O)$_p$ water molecules as a ``grouped'' water cluster  ($[A|pB_g|0]_{ADF}$ and $[A|pB_g|0]_{Dirac}$. The values drawn on Figure~\ref{graph:All_snap10_bis} reveal very similar trends between the SR-ZORA and DC values. For both Hamiltonians, the difference between FnT ($[A|pB_g|0]$) values and the supermolecular $[A+pB|0|0]$), vanishes at $p=6$ water molecules.

\subsection*{Statistical analysis of $\sigma_{iso}$ values; selection of a snapshot subset}
Figure~\ref{graph:All_ADF_DIRAC_b} plots the distribution of the $\sigma_{iso}$ values obtained from the 517~snapshots (red dotted line). With the aim of performing DC calculations, which are computationally more expensive that the ADF ones, it was necessary to reduce the number of snapshots, but keeping a normal distribution. The Normal probability plots of Figures~\ref{graph:Henry_line_ADF_FnT_216} and~\ref{graph:Henry_line_Dirac_et_ADF_216}, reveal that the ADF values, both for supermolecular calculations (dashed red line) and FnT calculations (black line) obtained for a subset of 216~snapshots follow a normal law, though with a smaller variance. The DC values (red lines) have the same statistical distribution, shifted by \SI{411}{ppm}.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-All-results-Dirac-with-ADF-embedding-216/fig-All-results-Dirac-with-ADF-embedding-216}
\caption{{\label{graph:All_ADF_DIRAC_b} Normalized Gaussians representing $\sigma_{iso} $ (in ppm).
The red dotted line refers to supermolecular ADF for all snapshots ($ [A+20B|0|0]_{ADF} $), while the red dashed line refers to the subset of 216 selected snapshots ($ [A+20B|0|0]^{216}_{ADF}$). The FnT ADF calculations ($ [A|11B_f|9B_f]^{216}_{ADF} $) for this 216 snapshot subset are drawn with a black line, while the DIRAC calculations using SR-ZORA (ADF) embedding potentials and electrostatic potentials are drawn in red.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-Henry-line-ADF-216-FnT/fig-Henry-line-ADF-216-FnT}
\caption{{\label{graph:Henry_line_ADF_FnT_216} Normal Probability Plot of $ \sigma_{iso} $ for the 216 
selected snapshots with ADF (FnT)~:~$ [A|11B_f|9B_f]^{216}_{ADF} $.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fig-Henry-line-Dirac-et-ADF-216/fig-Henry-line-Dirac-et-ADF-216}
\caption{{\label{graph:Henry_line_Dirac_et_ADF_216} Normal Probability Plot of $ \sigma_{iso} $ for the 216 
selected snapshots with $ \DIRAC $~:~$ [A|0|0]^{216-ADF}_{Dirac} $.%
}}
\end{center}
\end{figure}

\clearpage

\selectlanguage{english}
\FloatBarrier
\section*{References}\sloppy
\phantomsection
\label{csl:1}\textit{{Recent Developments in Absolute Shielding Scales for Nmr Spectroscopy}}. (2019).

\phantomsection
\label{csl:2}\textit{{Mathrm ^ 13 C $ Shielding Scale for Mas Nmr Spectroscopy}}. (2011).

\phantomsection
\label{csl:3}(2012).

\phantomsection
\label{csl:4}(2572).

\phantomsection
\label{csl:6}\textit{{The Paw/gipaw Approach for Computing Nmr Parameters: A New Dimension Added to Nmr Study of Solids}}. (2011).

\phantomsection
\label{csl:8}\textit{{Quantum-Chemical Embedding Methods for Treating Local Electronic Excitations in Complex Chemical Systems}}. (2012).

\phantomsection
\label{csl:9}\textit{{Subsystem Density-Functional Theory}}. (2014).

\phantomsection
\label{csl:10}\textit{{Frozen--Density Embedding Strategy for Multilevel Simulations of Electronic Structure}}. (5891).

\phantomsection
\label{csl:11}\textit{{Quantum Embedding Theories}}. (2705).

\phantomsection
\label{csl:12}(7673).

\phantomsection
\label{csl:13}\textit{{C Atalytic R Eduction of D Initrogen to A Mmonia by M Olybdenum T Heory Versus E Xperiment}}. (5512).

\phantomsection
\label{csl:14}\textit{{Selective Oxidation and Ammoxidation of Propene on Bismuth Molybdates, Ab Initio Calculations}}. (2001).

\phantomsection
\label{csl:16}(2014).

\phantomsection
\label{csl:17}(1095).

\phantomsection
\label{csl:18}(2005).

\phantomsection
\label{csl:20}\textit{{Molecular Properties from Combined Qm/mm Methods. 2. Chemical Shifts in Large Molecules}}. (3721).

\phantomsection
\label{csl:21}\textit{{Toward the Quantum Chemical Calculation of Nmr Chemical Shifts of Proteins. 2. Level of Theory, Basis Set, and Solvents Model Dependence}}. (1480).

\phantomsection
\label{csl:22}(9924).

\phantomsection
\label{csl:23}\textit{{Calculation of Nuclear Magnetic Resonance Shieldings Using Frozen-Density Embedding}}. (1941).

\phantomsection
\label{csl:24}\textit{{Nmr Solvent Shifts of Acetonitrile from Frozen Density Embedding Calculations}}. (2640).

\phantomsection
\label{csl:25}\textit{{Calculation of Nuclear Spin-Spin Coupling Constants Using Frozen Density Embedding}}. (1041).

\phantomsection
\label{csl:26}\textit{{On the Calculation of Second-Order Magnetic Properties Using Subsystem Approaches in a Relativistic Framework}}. (8400).

\phantomsection
\label{csl:27}\textit{{Relativistic Effects in Structural Chemistry}}. (1988).

\phantomsection
\label{csl:28}\textit{{Em Introduction to Relativistic Quantum Chemistry}}. (2007).

\phantomsection
\label{csl:29}\textit{{Em Relativistic Quantum Chemistry. T He Fundamental Theory of Molecular Science}}. (2009).

\phantomsection
\label{csl:30}\textit{{Perspective: Relativistic Effects}}. (1509).

\phantomsection
\label{csl:31}\textit{{Em Calculation of Nmr and Epr Parameters: Theory and Applications}}. (2004).

\phantomsection
\label{csl:32}\textit{{Relativistic Calculations of Magnetic Resonance Parameters: Background and Some Recent Developments}}. (2011).

\phantomsection
\label{csl:33}\textit{{Coupled-Cluster Theory in Quantum Chemistry}}. (2007).

\phantomsection
\label{csl:34}(2012).

\phantomsection
\label{csl:35}\textit{{Expression of Results in Quantum Chemistry}}. (1978).

\phantomsection
\label{csl:36}\textit{{Linear Response At the 4-Component Relativistic Level: Application to the Frequency-Dependent Dipole Polarizabilities of the Coinage Metal Dimers}}. (2003).

\phantomsection
\label{csl:37}\textit{{Relativistic Hamiltonians for Chemistry: a Primer}}. (3077).

\phantomsection
\label{csl:38}\textit{{A Local Exchange-Correlation Potential for the Spin Polarized Case. I}}. (1629).

\phantomsection
\label{csl:39}\textit{{Inhomogeneous Electron Gas}}. (1912).

\phantomsection
\label{csl:40}\textit{{Spin in Density--Functional Theory}}. (3661).

\phantomsection
\label{csl:41}\textit{{Relativistic Adiabatic Time-Dependent Density Functional Theory Using Hybrid Functionals and Noncollinear Spin Magnetization}}. (2091).

\phantomsection
\label{csl:42}(1041).

\phantomsection
\label{csl:43}\textit{{A Simple Scheme for Magnetic Balance in Four-Component Relativistic Kohn--Sham Calculations of Nuclear Magnetic Resonance Shielding Constants in a Gaussian Basis}}. (141AD).

\phantomsection
\label{csl:44}\textit{{Embedding Methods in Quantum Chemistry}}. (2018).

\phantomsection
\label{csl:45}\textit{{Frozen Density Functional Approach for Ab Initio Calculations of Solvated Molecules}}. (8050).

\phantomsection
\label{csl:46}\textit{{Kohn-Sham Equations with Constrained Electron Density: an Iterative Evaluation of the Ground-State Electron Density of Interacting Molecules}}. (1996).

\phantomsection
\label{csl:47}\textit{{Em Recent Advances in Density Functional Methods , Chapter~24, Page 371}}. (2002).

\phantomsection
\label{csl:48}\textit{{Molecular Properties Via a Subsystem Density Functional Theory Formulation: a Common Framework for Electronic Embedding}}. (441AD).

\phantomsection
\label{csl:49}\textit{{Frozen Density Functional Approach for Ab Initio Calculations of Solvated Molecules}}. (8050).

\phantomsection
\label{csl:50}\textit{{A Flexible Implementation of Frozen-Density Embedding for Use in Multilevel Simulations}}. (1011).

\phantomsection
\label{csl:51}\textit{{Subsystem Density-Functional Theory}}. (2014).

\phantomsection
\label{csl:52}(1432).

\phantomsection
\label{csl:53}(1641).

\phantomsection
\label{csl:54}\textit{{Excitation Energies from Frozen-Density Embedding with Accurate Embedding Potentials}}. (2341).

\phantomsection
\label{csl:55}(6317).

\phantomsection
\label{csl:56}\textit{{Th\' E Orie Quantique Des Courants Interatomiques Dans Les Combinaisons Aromatiques}}. (1937).

\phantomsection
\label{csl:57}\textit{{Ab Initio Methods for the Calculation of Nmr Shielding and Indirect Spin--Spin Coupling Constants}}. (1999).

\phantomsection
\label{csl:58}\textit{{Gauge Origin Independent Calculations of Nuclear Magnetic Shieldings in Relativistic Four-Component Theory}}. (1241).

\phantomsection
\label{csl:59}\textit{{Molecular Electronic Property Density Functions: The Nuclear Magnetic Shielding Density}}. (5684).

\phantomsection
\label{csl:60}\textit{{Nuclear Magnetic Shielding Density}}. (3366).

\phantomsection
\label{csl:61}\textit{{4-Component Relativistic Calculation of the Magnetically Induced Current Density in the Group 15 Heteroaromatic Compounds}}. (2009).

\phantomsection
\label{csl:62}\textit{{4-Component Relativistic Magnetically Induced Current Density Using London Atomic Orbitals}}. (2068).

\phantomsection
\label{csl:63}\textit{{Car-Parrinello Molecular Dynamics}}. (2011).

\phantomsection
\label{csl:64}(2328).

\phantomsection
\label{csl:65}\textit{{Dirac--Fock Atomic Electronic Structure Calculations Using Different Nuclear Charge Distributions}}. (1997).

\phantomsection
\label{csl:66}\textit{{Finite Nuclear Charge Density Distributions in Electronic Structure Calculations for Atoms and Molecules}}. (2000).

\phantomsection
\label{csl:67}\textit{{Generalized Gradient Approximation Made Simple}}. (3865).

\phantomsection
\label{csl:68}\textit{{Generalized Gradient Approximation Made Simple [phys. Rev. Lett. 77, 3865}}. (1396).

\phantomsection
\label{csl:69}\textit{{Obtaining a Gradient-Corrected Kinetic-Energy Functional from the Perdew-Wang Exchange Functional}}. (5328).

\phantomsection
\label{csl:70}\textit{{The Zero-Order Regular Approximation for Relativistic Effects: The Effect of Spin-Orbit Coupling in Closed Shell Molecules}}. (6505).

\phantomsection
\label{csl:71}\textit{{Relativistic Total Energy Using Regular Approximations}}. (9783).

\phantomsection
\label{csl:72}\textit{{Relativistic Regular Two-Component Hamiltonians}}. (1996).

\phantomsection
\label{csl:73}\textit{{Geometry Optimizations in the Zero Order Regular Approximation for Relativistic Effects}}. (8943).

\phantomsection
\label{csl:74}\textit{{Https://www.scm.com}}.

\phantomsection
\label{csl:75}\textit{{Optimized Slater-Type Basis Sets for the Elements 1-118}}. (1142).

\phantomsection
\label{csl:76}\textit{{Written by T.~saue, L.~visscher, H.~j.~ \relax Aa .~jensen and R.~bast, with Contributions from V.~bakken, K.~g.~dyall, S.~dubillard U.~ekstr \"o M, E.~eliav, T.~enevoldsen, E.~fa \ss Hauer, T.~fleig O.~fossgaard, A.~s.~p.~gomes, E.~d.~hedeg \aa Rd, T.~helgaker, J.~henriksson M.~ilia \v S , Ch.~r.~jacob, S.~knecht, S.~komorovsk \'y , O.~kullie J.~k.~l \ae Rdahl, C.~v.~larsen, Y.~s.~lee, H.~s.~nataraj, M.~k.~nayak P.~norman, G.~olejniczak, J.~olsen, J.~m.~h.~olsen, Y.~c.~park J.~k.~pedersen, M.~pernpointner, R.~di~remigio, K.~ruud, P.~sa \l Ek B.~schimmelpfennig, A.~shee, J.~sikkema, A.~j.~thorvaldsen, J.~thyssen J.~van~stralen, S.~villaume, O.~visser, T.~winther, and S.~yamamoto Available At \url Https://doi.org/10.5281/zenodo.2253986 , See Also Url Http://www.diracprogram.org}}. (2018).

\phantomsection
\label{csl:77}\textit{{Relativistic Double-Zeta, Triple-Zeta, and Quadruple-Zeta Basis Sets for the 4d Elements Y  Cd}}. (2007).

\phantomsection
\label{csl:78}\textit{{Gaussian Basis Sets for Use in Correlated Molecular Calculations. I the Atoms Boron Through Neon and Hydrogen}}. (1007).

\phantomsection
\label{csl:79}\textit{{Matplotlib: A 2d Graphics Environment}}. (2007).

\phantomsection
\label{csl:80}\textit{{Mayavi: 3d Visualization of Scientific Data}}. (2011).

\phantomsection
\label{csl:81}\textit{{C Onventions for the Reporting of Nuclear Magnetic Shielding (or Shift) Tensors Suggested by Participants in the N A T O A R W on N M R Shielding Constants At the U Niversity of M Aryland, C Ollege P Ark, J Uly 1992}}. (1993).

\phantomsection
\label{csl:82}(2006).

\phantomsection
\label{csl:83}\textit{{Molecular Electric, Magnetic, and Optical Properties}}. (2017).

\phantomsection
\label{csl:84}\textit{{Spin-Orbit Effects on the $\mathrm ^ 125 Te $ Magnetic-Shielding Tensor: A Cluster-Based Zora/dft Investigation}}. (2018).

\phantomsection
\label{csl:85}\textit{{Relativistic Calculation of Nuclear Magnetic Shielding Tensor Including Two-Electron Spin-Orbit Interactions}}. (2006).

\phantomsection
\label{csl:86}\textit{{Computational Multinuclear Nmr of Platinum Complexes: A Relativistic Four-Component Study}}. (4908).

\phantomsection
\label{csl:87}\textit{{A Dft/zora Study of Cadmium Magnetic Shielding Tensors: Analysis of Relativistic Effects and Electronic-State Approximations}}. (1785).

\end{document}

