Frozen-Density Embedding for including environmental effects in the Dirac-Kohn-Sham theory: an implementation based on density fitting and prototyping techniques

The Frozen Density Embedding scheme represents an embedding method in which environmental effects onto a given subsystem are included by representing the other subsystems making up the surroundings quantum mechanically, by means of their electron densities. In the present paper, we extend the full 4-component relativistic Dirac-Kohn-Sham method, as implemented in the BERTHA code, to include environmental and confinement effects with the FDE scheme. This implementation has been enormously facilitated by BERTHA's python API (PyBERTHA), which provides a flexible framework of development by using all Python advantages in terms of code re-usability, portability while facilitating the interoperability with other FDE implementations available through the PyADF framework. The computational performance has been evaluated on a series of gold clusters (Au$_n$, with n=2,4,8) embedded into an increasing number of water molecules (5, 10, 20, 40 and 80 water molecules). We found that the procedure scales approximately linearly both with the size of the frozen surrounding environment (in line with the underpinnings of the FDE approach) and with the size of the active system (in line with the use of density fitting). Finally, we applied the code to a series of Heavy (Rn) and Super-Heavy elements (Cn, Fl, Og) embedded in a C_60 cage to explore the confinement effect induced by C_60 on their electronic structure. We compare the results from our simulations with more approximate models employed in the atomic physics literature, in which confinement is represented by a radial potential slightly affected by the nature of the central atom. Our results indicate that the specific interactions described by FDE are able to improve upon the cruder approximations currently employed, and thus provide a basis from which to generate more realistic radial potentials for confined atoms.


Introduction
Molecular systems, clusters, and materials containing heavy atoms have drawn considerable recent attention because of their rich chemistry and physics. [1][2][3][4] In order to model computationally systems containing heavy elements, the methods of relativistic quantum mechanics 2 must be necessarily adopted to capture scalar and spin-orbit interactions that are neglected in the conventional non-relativistic formulation of quantum chemistry. Furthermore, most of the chemistry occurs in solution and the environment plays a key role in the determination of the properties and reactivity of substances in condensed phases. [5][6][7][8][9] Thus, the complexity of chemical phenomena in solution has made it necessary to develop a variety of models and computational techniques to be combined with (relativistic) quantum chemistry methods.
Among different approaches to include environmental effects we mention the quantum mechanics/molecular mechanics (QM/MM) approach, 10 which includes the molecular environment explicitly and at a reduced cost using classical mechanical description, or in polarizable continuous medium (PCM) 11 (i.e., where the solvent degree of freedom are replaced by an effective classical dielectric). Despite being widely and successfully applied these methods may have drawbacks. For instance methods based on PCM cannot describe specific interactions with the environments (e.g. hydrogen, halogen bonds), while the QM/MM approach, which is based on classical force fields, may be limited by the availability of accurate parameterizations which may reduce its predictive power, in particular when heavy elements are involved.
An alternative is to use quantum embedding theories (for an overview see Ref. 12-15 and references therein), in which a QM description for a subsystem of interest is combined with a QM description of the environment (QM/QM). A notable example of QM/QM methods is the frozen-density embedding (FDE) scheme introduced by Wesołowski and Warshel, 16,17 based on the approach originally proposed by Senatore and Subbaswamy, 18 and later Cortona,19 for solid-state calculations. The method has been further generalized 20,21 and directed to the simultaneous optimization of the subsystems electronic densities.
FDE is a DFT-in-DFT embedding method that allows to partition a larger Kohn-Sham system into a set of smaller, and coupled, Kohn-Sham subsystems. The coupling term is defined by a local embedding potential depending only on the electron densities of both the sole active subsystem and the environment (i.e. no orbital information is shared among subsystems). This feature gives to the FDE scheme an enormous flexibility, as indeed virtually arbitrary methods can be combined to treat different subsystems. For example, wavefunction theory (WFT) based methods can be used for the active system while one can take advantage of the efficiency of DFT to describe a large environment (WFT-in-DFT). 12,13,17,22-25 Also one can employ very different computational protocols for different subsystem including i) using Hamiltonian dealing with different relativistic approximations (from the full 4-component methods to the non-relativistic ones); [26][27][28][29] ii) different basis sets size and type (Gaussian and Slater type functions, relativistic four-component spinors) and even iii) different quantum chemical packages. 26,30,31 We mention that the FDE scheme has been extended both to the linear-response TDDFT, [32][33][34] including to account for charge-transfer excitations 35,36 and to real time TDDFT (rt-TDDFT). 31,37 FDE-based calculations are shown to be accurate in the case of weakly interacting systems including hydrogen bond systems, 38,39 whereas, their use for subsystems interacting with a larger covalent character is problematic (see Ref. 38 and references therein) due to the use of approximate kinetic energy functional (KEDF) in the non-additive contribution to the embedding potential. The research for more accurate KEDFs is a key aspect for the applicability of the FDE scheme as a general scheme, [40][41][42][43] including the partitioning of the system also breaking covalent bonds. 44 We mention here that alternative QM/QM approaches, avoiding the use of KEDFs and allowing also for fragmentation in subsystems through covalent bonds, have been recently proposed (see for instance Ref. 14,15,[45][46][47][48][49][50][51] ).
Thanks to its flexibility the FDE scheme has been implemented in different flavors into computational packages such as: embedded Quantum Espresso, 52 ADF, 21,53,54 Turbomole, 55,56 Dalton, 26,57 Koala, 58 Molpro, 45 Serenity, 59 and Q-Chem, 60 (the first two based on plane waves and Slater type functions respectively, the others on Gaussian type functions).
FDE has also been implemented to treat the subsystems at full-relativistic four-component level based on the Dirac equation within the DIRAC code, 61 and can be used with DFT and different wavefunction methods both for molecular properties and energies involving the ground or excited electronic states. 26,28,29,[62][63][64] Despite its conceptual simplicity, its actual implementations may lead to relatively complicated workflows. A simpler approach is therefore to integrate such legacy codes as computational engines to handle the different FDE steps, which are then glued together and their execution automatized using suitable frameworks such as for instance that implemented in PyADF, 30,65 that can be easily extensible due to its object-oriented implementation in the Python programming language. 66 Prototyping techniques also based on Python are very useful to build reference implementations, as for instance the Psi4-rt-PyEmbed code, 31,67 where the Python interface of Psi4Numpy and PyADF 30,68 (including its PyEmbed module 69,70 and XCFun library 71,72 to evaluate non-additive exchange-correlation and kinetic energy contributions) has been used by some of us to build real time non-relativistic TDDFTin-DFT FDE 31 and projection-based embedding 73 implementations.
In this work we extend the Dirac-Kohn-Sham (DKS) method implemented in the BERTHA code (with its new Python API, PyBERTHA) 74,75 to the FDE scheme to include environmental/confinement effects in the DKS calculations (DKS-in-DFT FDE). The implementation takes advantages of the DKS formulation implemented in BERTHA, including the density fitting algorithms at the core of the computation (i.e. in the evaluation of the embedding potential and of its matrix representation is relativistic G-spinor functions), and the FDE implementation already available in the PyEmbed module of the PyADF framework.
The outline of the paper is as follows. In Section 2 we present the basic theory of FDE and a brief description of the DKS method as implemented in BERTHA. In section 3, we then describe in detail our implementation. In section 4 we present some numerical results, including the computational burden and scalability of this new implementation with respect to the size of the active system as well as of the embedding one. We will also present an application to a series of Heavy (Rn) and Super-Heavy elements (Cn, Fl, Og) confined into a C 60 cage. Finally, concluding remarks are given in section 5.

Theory
In this section, we briefly review the basic formalism of the FDE scheme and its extension to use the DKS theory for the active system (DKS-in-DFT FDE). We will remark also some details of the DKS implementation in BERTHA, mainly focusing on those aspects (including density fitting techniques), which are relevant for an efficient implementation of the FDE scheme. Finally, we will illustrate the basic characteristics of the our recent BERTHA Python API, PyBERTHA (and the related pyberthamod module available under GPLv3 license at Ref. 67, for additional and technical details see Refs. [74][75][76] which is a key tool here to devise a simple work-flow for the DKS-in-DFT FDE scheme.

Subsystem DFT and Frozen Density Embedding formulation
In the subsystem formulation of DFT the entire system is partitioned into N subsystems, and the total density ρ tot (r) is represented as the sum of electron densities of the various subsystems [i.e., ρ a (r) (a = I, .., N )]. In the following we consider the total density as partitioned in only two contributions as ρ tot (r) = ρ I (r) + ρ II (r). (1) The total energy of the system can then be written as with the energy of each subsystem (E i [ρ i ], with i = I, II) given according to the usual definition in DFT as In the above expression, v i nuc (r) is the nuclear potential due to the set of atoms which defines the subsystem, and E i nuc is the related nuclear repulsion energy. T s [ρ i ] is the kinetic energy of the auxiliary non-interacting system, which is, within the Kohn-Sham (KS) approach, commonly evaluated using the KS orbitals. The interaction energy is given by the expression: with v I nuc and v II nuc being the nuclear potentials due to the set of atoms associated with the subsystem I and II, respectively. The repulsion energy for nuclei belonging to different The electron density of a given fragment (ρ I or ρ II in this case) can be determined by minimizing the total energy functional (Eq.2) with respect to the density of the fragment while keeping the density of the other subsystem frozen. This procedure is the essence of the FDE scheme and leads to a set of Kohn-Sham-like equations (one for each subsystem) which are coupled by the embedding potential term v I emb (r), which carries all dependence on the other fragment's density. Here T denotes the kinetic energy operator, which in a nonrelativistic framework has the form −∇ 2 /2, whereas for a relativistic framework is cα · p (see discussion below). We also note that in the relativistic framework, the FDE expressions above correspond to the case in which an external vector potential is absent. Further details for their generalization can be found in Ref. 28 where the non-additive exchange-correlation and kinetic energy contributions are defined as the difference between the associated exchange-correlation and the kinetic potentials defined using ρ tot (r) (i.e., ρ I (r) + ρ II (r)) and ρ I (r). For both potentials one needs to account for the fact that only the density is known for the total system, so that potentials that require input in the form of KS orbitals are prohibited. For the exchange-correlation potential, one may make use of accurate density functional approximations and its quality is therefore similar to that of ordinary KS. The potential for the non-additive kinetic term ( δT nadd δρ I (r) , in Eq.6) is more problematic as less accurate orbital-free kinetic energy density functionals (KEDFs) are available for this purpose. Examples of popular functional approximations applied in this context are the Thomas-Fermi (TF) kinetic energy functional 77 or the GGA functional PW91k. 78 As already mentioned in the introduction, the research for more accurate KEDFs is a key aspect for the applicability of the FDE scheme as a general scheme, including the partitioning of the system also breaking covalent bonds. 44 In general, the set of coupled equations that arise in the FDE scheme for the subsystems have to be solved iteratively and a freeze-and-thaw scheme, where one relaxes the electron density of one subsystem at a time keeping frozen the others, until electron densities of all subsystems reach a required convergence. In this work we limit to one subsystem (active) while keeping the density of the environment frozen to their ground state density. In this case, the implementation of FDE reduces to the evaluation of v I emb (r) potential (which is an one-electron operator) that needs to be added to the Hamiltonian of the active system.
The matrix representation of the embedding potential may be evaluated using numerical integration grids, similar to those used for the exchange-correlation term in the KS method.

8
This contribution is then added to the KS matrix and the eigenvalue problem is solved with the usual self-consistent field (SCF) procedure. We note here that two approaches can be taken: the first is the use of a pre-calculated embedding potential 26 (for instance, from a prior subsystem DFT calculation) that is used as a one-body operator (referred to in the literature as a "static" embedding potential) added to the one-body Fock matrix at the start of the (4-component) calculations. The second approach involves the regeneration of the embedding potential using the (4-component) actual electron density of the active system.
In this case, the matrix representation of the embedding potential is updated during SCF procedure, because of its dependence on the active subsystem density (see Eq.6) that itself changes during the SCF iterations. As discussed below, in this work we shall mostly make use of the latter approach.

Dirac-Kohn-Sham scheme in BERTHA and its extension to FDE based on density fitting
For the detailed theoretical basis of the Dirac-Kohn-Sham methodology we refer the reader to previous works [79][80][81][82][83][84][85] and references therein. Here, we summarize only the main aspects of the implementation of an all-electron DKS method based on the use of G-spinor basis sets and the density-fitting techniques as is implemented in BERTHA. 75,86 In atomic units, and including only the longitudinal electrostatic potential, the DKS equation reads where c is the speed of light in vacuum, p is the electron momentum, while: where σ = (σ x , σ y , σ z ), σ q is a 2 × 2 Pauli spin matrix and I is the 2 × 2 identity matrix.
The longitudinal interaction term is represented by a diagonal operator borrowed from nonrelativistic theory and made up of: a nuclear potential term v N (r), a Coulomb interaction We mention that the Breit interaction contributes to the transverse part of the Hartree interaction and is not considered here, as we restrict ourselves to using non-hybrid, non-relativistic functionals of the electron density.
In BERTHA, the spinor solution (Ψ i (r) in Eq. 7) is expressed as a linear combination of the G-spinor basis functions, 87 M T µ (r) (T = L, S with L and S refer to the so-called "large" and "small" component, respectively). The G-spinors do not suffer from the variational problems of kinetic balance (see Ref. 88 and references therein) and, regarding the evaluation of multicentre integrals, retain the advantages that have made Gaussian-type functions the most widely-used expansion set in non-relativistic quantum chemistry. The matrix representation of the DKS operator in the G-spinor basis is given by where The eigenvalue equation in the algebraic representation is given by where c (T ) are the spinor expansion vectors. The H DKS matrix is defined in terms of the v (T T ) , J (T T ) , K (T T ) , S (T T ) , and Π (T T ) matrices, being respectively the basis representation of the nuclear, Coulomb, and exchange-correlation potentials, the overlap matrix, and the matrix of the kinetic operator, respectively. The nuclear charges have been modeled by a finite Gaussian distribution. 89 The resulting matrix elements are defined by The terms ρ ν (r)), which can be exactly expressed as linear combination of standard Hermite Gaussian-type functions (HGTFs). 86,87,90 The H DKS matrix depends on ρ(r) in v νi , with T and T equal to both L and S), where the sum runs over the occupied positive-energy states. The total electron density is obtained according to The computation of the Coulomb and exchange-correlation contributions to the DKS matrix, that is Eq. 13 and 14, respectively, is the most demanding computational step in a DKS calculation involving a G-spinor basis set. The current version of BERTHA takes advantage of both density fitting 86,91-93 and of advanced parallelization techniques 74,75,[94][95][96] for the evaluation of these two contributions. The relativistic density (which is a real scalar function) is thereby expanded in a set of N aux auxiliary atom-centered functions.
In the Coulomb metric the expansion coefficients d t are defined as the solution of the linear system, given by where A is a real and symmetric matrix, representing the Coulomb interaction in the auxiliary basis, A st = f s || f t while the elements (v s ) of the vector v are the projection of the electrostatic potential on the fitting functions, f s || ρ . In our implementation the latter integrals are efficiently evaluated using the relativistic generalization of the scalar Hermite density matrix proposed by Almlöf. 97,98 For the exchange-correlation K matrix we adopt a similar strategy by solving for z in the linear system where the vector w is the projection of the exchange-correlation potential (ṽ The elements of the vector w, which involve integrals of the exchange-correlation potential, are computed numerically by the integration scheme already implemented in the code. 99 Once the vectors d and z have been worked out, the Coulomb and the exchange-correlation contributions to the DKS matrix can be evaluated in terms of 3-center two-electron integrals The extension of this scheme to include the contribution FDE potential, v I emb [ρ I , ρ II ](r), is straightforward. For the evaluation of the embedding potential matrix representation in G-spinors (Ṽ emb(T T ) µν ) we can strictly follow the procedure already employed above for the exchange-correlation term. We solve for c in the linear system where the vector g is the projection of the embedding potential, v I emb [ρ I , ρ II ](r), on the fitting functions The elements of the vector g are computed numerically on a suitable integration grid (further details of the implementation will be given in the next section). Once the vector g has been computed the embedding potential contribution can be evaluated in a single step together with the Coulomb and the exchange-correlation ones (see Eq.25), and finally added to the DKS matrix.J Note that, in the evaluation of the embedding potential on a numerical grid,ṽ emb (r) = v emb [ρ I (r), ρ II (r)], we use the fitted density of the active system,ρ I (see Eq.18). An important aspect for the computational efficiency arises from the use of an auxiliary fitting basis set of primitive HGTFs. Indeed, they are grouped together in sets sharing the same exponents. 91 In particular, each set is defined so that a specified auxiliary function of a given

The PyBerthaEmbed DKS-in-DFT FDE implementation
In this section we outline the computational strategy we adopted to implement the DKS-
Algorithm 1 reports the most important part of the pyberthaemb.py code, and it well illustrates how we can gain a relatively simple work-flow to implement FDE using the DKS level of theory for the active system using PyBERTHA and the new pyemb class for the environmental system. The pybertha class is instantiated (line 4) with the shared object bertha_wrapper.so specified as an input. The SO contains the cited c_wrapper and bertha_wrapper code which are based on the core FORTRAN libraries, namely: libbertha.so and libberthaserial.so (see Figure 1 above). After the initialization, the full DKS calculation is worked out (line 10) using the bertha.run() method. At line 14 the pyemb class is also instantiated specifying the files (specified in xyz format) for the geometries of both the active and embedding systems. The quantum chemistry software employed for the actual calculation of the environment system is specified at this stage (in the current example, and throughout this work, we used the ADF package 53 ).
All the details for the computation of the embedding system are set at line 15. This includes: the selection of the type of basis set functions, Hamiltonian, exchange-correlation functional and also the non-additive kinetic functional used to define the embedding potential. At this stage all the basis sets and exchange-correlation functionals available in the ADF library can be used. Similarly, the numerical integration grid used for the numerical representation of the embedding potential is set, both the type (global grid or active system) and the quality. By default the numerical grids internally defined by the ADF program are used, however other options are available, including the possibility to use an user-defined grid (see for instance line 16, commented).
The embfactory.initialize() method performs a stand-alone single point calculation on the embedding system. The method evaluates the nuclear and Coulomb potentials of the environment (see Eq.6) and its ground state electron density (ρ II ). All these quantities are mapped on the numerical grid. At this time, the numerical grid defined within PyADF is made available (as a numpy.array) using the get_grid() method (line 19) and used as an input for the get_density_on_grid() method of the pybertha class. This method allows to define the ground state densityρ I of the active system at DKS level of theory.ρ I is also available as a numpy.array that, after a reshape (line 22), can be used as an input of the get_potential() method of the pyemb class (line 25) to obtain the final embedding potential.
After this initial setup, we proceed to the actual FDE calculation. In this example, the embedding potential will be generated using the active subsystem density (loop structure, lines 27 to 40), using the split-SCF scheme. 104 Thus, in each of the spin-SCF iterations, the new set_embpot_on_grid() method of pybertha class makes both the numerical grid and the embedding potential available at the FORTRAN layer. Thus, the numerical integration of the v emb (r) on the fitting basis functions (Eq.24) and linear system solution (Eq.23) are efficiently evaluated in FORTRAN. The DKS matrix employed in the intervening BERTHA calculation (bertha.run(), line 31) is updated using the G-spinor representation of the embedding potential, and here the split-SCF scheme is interesting as it does not require the evaluation of the embedding potential at each SCF step taking place on the BERTHA side. The new fitted density (line 32) is used to compute a new embedding potential (line 36) which is used in for the next iteration of the split-SCF procedure. 104 This scheme is iterated till a convergence criteria is satisfied (line 40).
In Figure 2 we present the workflow which emphasizes the interoperability between different tasks and modules or programs involved including the layers where the actual computations are carried out. This schematic picture highlights also how the key quantities, required to implement the DKS-in-DFT FDE scheme, have different representations along the computation. As example, we focus on the electron density of the active system,ρ I (r).
This quantity is evaluated at the DKS level of theory activated by the bertha.run() method within the PyBERTHAembed program. The actual calculation is done within the FOR-TRAN layer. At this level,ρ I (r) is represented in terms of the expansion coefficients of auxiliary fitting functions (d, see Eq.19) and is stored as a FORTRAN array (of dimension N aux ). However, this representation is not useful itself for the evaluation of the embedding potential. Indeed, its evaluation and in particular the non-additive contribution requires thatρ I (r) is represented on a grid. Thus, within pyberthaemb.py, the numerical grid (GRID) evaluated within PyADF (see panel init) is made available as numpy.array and via the bertha.get_density_on_grid() method is made accessible to the FORTRAN layer (bertha_wrapper module), and stored as a FORTRAN array of dimension npoints.
The calculation of the numerical representation ofρ I (r) on the grid is done efficiently in FOR-TRAN (see in Figure 2, panel a) and the latter is accessible within pyberthaemb.py as a numpy.array. Analogously, different representations are used also for the embedding potential along the workflow. Note that all quantities accessible from pyberthaemb.py, namelỹ ρ I (r), v emb (r k ) and GRID (labelling the arrows in figure), are defined as numpy.array that can be easily manipulated within a Python source code. The computational steps which involve BERTHA are instead implemented in FORTRAN (panel a) and panel c) in Figure) 19 and have been efficiently parallelized using OpenMP.  In the present section we report a series of numerical results mainly devoted to assess the correctness of our new implementation of DKS-in-DFT FDE scheme. In addition we are also reporting the computational cost and scalability with respect to the size of both the active and the embedding system. Finally, we will present an application to a series of Heavy (Rn) and Super-Heavy elements (Cn, Fl, Og) confined into a C 60 cage. In the adduct the water molecule is the active system that is bound to an ammonia molecule, which instead plays the role of the embedding environment. In the Psi4-rt-PyEmbed case we use basis sets obtained by the decontraction of the Gaussian cc-pVDZ, cc-pVTZ and aug-cc-pVDZ basis sets 106,107 for the active system (these basis sets are referred as cc-pvdz-decon, cc-pvtz-decon and aug-cc-pvdz-decon, respectively). The same basis set has been used in the DKS calculation to define the large component of the G-spinor basis set.
The corresponding small component was generated using restricted kinetic balance relation. 88 Noteworthy, for these reference calculations, we have used an extremely large auxiliary fitting basis set (A4 spdf g ) which gives an error on the Coulomb energy even below 10 −6 Eh. The computational details for the definition of the environment, including parameters to define the embedding potential, are identical in both PyBERTHAembed and Psi4-rt-PyEmbed.
In particular, the basis set used in PyADF for the calculation of the environment frozen density (ammonia) and the embedding potential is the AUG-TZ2P Slater-type set from the ADF library. 108 As numerical grid we used the supermolecular Voronoi Polyhedra grid defined in ADF which is set defining an the integration parameter equal to 4 (this corresponds to a total number of 33280 grid points). The PBE 109 exchange-correlation functional has been used for the active system while the BLYP 110,111 exchange-correlation functional has been used for the ammonia molecule. The Thomas-Fermi and LDA functionals 112,113 has been employed for the non-additive kinetic and non-additive exchange-correlation potential, respectively. The molecular structure of the adduct is reported in SI (Table S.1). The effect of the environment (ammonia) on the active system (water) have been evaluated comparing the dipole moment components and diagonal elements of the polarizability tensor (α xx , α yy and α zz ) of the isolated (Free) respect to the embedded (Emb) water. We note here that the quoted polarizability values are not those for the supermolecular system, but only for the active subsystem.
The numerical results, reported in Table 1, show an evident quantitative agreement between the two implementations. Indeed, both the variations induced by the presence of the embedding system (∆ values) and the absolute values show a good agreement. Noteworthy, independently by the basis set used, the differences are below of 0.001 a.u. and 0.01 a.u. for the dipole moment components and for the polarizability tensor components, respectively.
We mention that we also performed the calculations increasing the speed of light by 1 order of magnitude (i.e., c = 1370.36 a.u.) to approximate the non-relativistic limit and, as expected, we obtain almost indistinguishable results (see Table S.2). All the above findings make us confident that our implementation is both numerically stable and correct.
As we have extensively described in the previous section, our implementation strongly 22 Table 1: Dipole moment (components µ x , µ y , µ z and module |µ|) and dipole polarizability (tensor diagonal components α xx , α yy , α zz and isotropic contribution α iso ) of both the isolated (Free) and embedded (Emb) water molecule. In the embedded water molecule, an ammonia molecule is used as environment. Data have been obtained using our new PyBERTHAembed implementation and the reference Psi4-rt-PyEmbed implementation (see text for details). The shift ∆ is also reported. All numerical data are reported in atomic units (a.u.). The diagonal components of the dipole polarizability tensor have been calculated using a finite field approach using an external electric of 0.001 a.u.

Psi4-rt-PyEmbed
PyBERTHAembed  Table 2: Dipole moment (components µ x , µ y , µ z and module |µ|) and dipole polarizability (tensor diagonal components α xx , α yy , α zz and isotropic contribution α iso ) of the embedded water molecule (water-ammonia system). Data have been obtained with our new PyBERTHAembed implementation (using a G-spinor basis functions derived from the ccpvtz-decon basis) and several auxiliary density fitting basis sets (A2 s , A2 sp , A2 spd , A2 spdf g , A3 spdf g and A4 spdf g ). The size of different fitting basis sets (N aux ) are reported in parenthesis. ∆E J is the absolute error on the Coulomb energy due to the density fitting. The diagonal components of the dipole polarizability tensor have been calculated using a finite field approach using an external electric of 0.001. All numerical data are reported in atomic units (a.u.). See text for the fitting basis set definition and further details.   Table 2. In the Table we also show the absolute error in the Coulomb energy (∆E J ), which is the quantity that is variationally optimized in the fitting procedure and typically regarded as its quality index. This numerical test shows that the use of density fitting does not introduce any significant instability in the DKS calculation of the active system, also in presence of the embedding potential. The Very accurate results can already be obtained starting from the A2 spd auxiliary basis set (i.e., 163 functions for the water molecule). It is interesting to note that the ∆E J associated with this basis set is of the same order of magnitude of that typically required (0.1 mEh per atom) in standard calculations based on density fitting without including FDE. Thus, these preliminary results suggest that the variational density fitting scheme can safely be applied in the implementation of DKS-in-DFT method without jeopardizing its accuracy.

Computational efficiency : gold clusters in water
It is interesting now to put forward some assessments on the computational efficiency of our DKS-in-DFT FDE implementation, together with its scaling properties in terms of time statistics and memory usage. This analysis will give us a detailed overview of the the computational burden, and possible bottlenecks, along the relatively complex workflow we implemented (using different quantum chemistry packages and programming languages).
Furthermore, it will be a solid starting point for future optimizations and developments (e.g. The results are reported in Table 3 and in Table 4, where, together with the total elapsed time (t d ) for each SCF iteration including the FDE contribution, we also partition between different tasks related with the FDE implementation, namely: a) numerical representation of active system fitted density on grid,ρ(r); b) calculation of the non-additive terms of embedding potential by PyADF (with the PyEmbed class); c) projection of the embedding potential onto fitting basis functions. In the Tables we also report the maximum memory usage for the SCF procedure ("Mem"), the number of points of grid and the timing for the "init" phase which involves: the evaluation of the ground state electronic density of the environmental together with the associated Coulomb potential, and their mapping on the numerical grid. We recall that the electron density of the environment is kept frozen, thus this initial step is done once at the beginning of the procedure. All tasks are also highlighted (using the same labeling: a, b, c and init) in Figure 2.

DKS-in-DKS
As general remark we may state that the FDE contribution to the total time is relatively small. By increasing the size of the active system (Au 2 , Au 4 and Au 8 ), and keeping fixed the environment (using ten water molecules), see Table 3, the relative impact of the FDE computational phase decreases. It passes from 13.3% for Au 2 (H 2 O) 10 to 0.9% for Au 8 This may be expected since tasks a), b) and c) have a more favorable scaling than the DKS calculation (i.e., O(N 3 )). The computational cost for the step a) and c) is proportional to the 26 product N aux · N gridpoints , where: N aux is the total number of the auxiliary fitting functions in the active system, and N gridpoints total number of grid points. Thus, the computational cost should scale as O(N 2 ) (being N the dimension of the active system). The actual scaling is much lower (i.e., slightly higher than O(N )) mainly because the total grid points are largely dominated by the environmental system (see number of grid points, N gridpoints , as reported in Table 3 ). Concerning the step b) and considering the fact that the environment is maintained fixed, it scales, as expected, linearly with number of points of the grid, N gridpoints .
The maximum use of memory, during the whole DKS-in-DFT FDE procedure, increases with respect to the number of Au atoms being N 1.7 , which is close to the theoretical value N 2 .
When we fix the active system (Au 4 ) increasing instead the size of the environment, see Table 4, the relative computational cost to include the embedding passes from 2.2%, in the case of Au 4 @(H 2 O) 5 , to 16.1% for Au 4 @(H 2 O) 80 . In this case, all tasks associated with the FDE procedure (a,b, and c) have a computational burden which increases linearly with the size of the environment (and the number of total grid points, see Figure S1 in SI), while the maximum memory usage during the SCF procedure is almost independent from the number of water molecules in environment, as only a slight increase can be observed.

The generation of atom-endohedral fullerenes model potentials
We conclude our work by showcasing how we can leverage our FDE implementation to determine fullerene-atom model potentials, that are applicable for species across the periodic where r 0 = 5.8 a.u. and U 0 = −0.30134 a.u. and ∆ = 1.9 a.u represent the finite thickness of the spherical potential. 123 Other model potentials have been proposed, 126 as well as other approaches to avoid numerical instability related to the sharp form of those potentials. 127,128 An alternative approach may be to start from the embedding potential generated in the FDE scheme, which is expected to be highly accurate and without artificial discontinuities. In Before proceeding in the comparison of different models, we have compared for the Rn atom the ability of FDE to capture environment effects on orbital energies, with respect to standard (supramolecular) DFT calculations. Our results, shown in Figure S2 in the supplementary information, show that there is a good agreement between the FDE and supramolecular calculations, with FDE yielding overall slightly larger orbital energy shifts, that are nevertheless very homogeneous across the different orbitals. From these results we conclude that FDE is, in effect, capable of correctly describing the fullerene cage's effect onto the atom's electronic structure.
All calculations, reported in the following, were carried out using a basis set for the  centered at nuclei positions, and negative values located in correspondence of the bonds. As one may expect the EMBP potential, while maintaining an overall spherical shape, is clearly different with respect to a simple short-range attractive V c (r) spherical potential. Indeed, if we consider the spherical average of the EMBP (see SI for details on the spherical average procedure employed) extracted for the various A@C 60 systems, as reported in Figure 4, while the EMBP seems to detect the same short-range attractive values surely it shows a more complex radial structure. The spherical average of the EMBP shows a positive repulsive value immediately before the inner C 60 surface and, maybe more importantly, never completely goes to zero, not even at the center of the fullerene where the atom A is placed. Not surprisingly, the final results (see Table 6) in terms of HOMO-LUMO gap for the spherical model potential are quite different respect to the results obtained using the FDE procedure. Indeed, as we mentioned, the overall shape of the EMBP is quite different respect to a simple spherical ones, see Figure 3. Nevertheless, is interesting to note as the spherical average seems to work well. Comparing the results obtained from the full FDE procedure respect to the ones computed using a model potential that is the spherical average of the EMBP, respectively columns 3 and 4 of Table 6, one can easily note that the spherical average is able to well reproduce the electronic structures of the active system (i.e., the central atom) including HOMO-LUMO gaps values with an error that is generally less the 1%. Similar conclusions can be drawn looking at Figure 5, where we report instead all the differences in orbital energies with respect to the isolated Rn atom for all the occupied orbitals. Once again both the EMBP (i.e., the FDE procedure) and its spherical average lead to similar results. Instead using the simple spherical model the energy shift is always the opposite. Finally, and maybe more interestingly, if we consider column 5 of Table 6 we see that we can obtain the HOMO-LUMO gaps for all the atoms with an error that is always less then 33 4% using a Rn based EMBP spherical average. Thus, we computed the spherical average of the EMBP for the Rn@C 60 system and using this single model potential we have been able to quantitatively well reproduce the HOMO-LUMO gap for all the other A@C 60 systems. This latter result let us clearly envision a practical approach to be used to build model potential as a result of the FDE procedure.

Conclusions and perspectives
Including environmental effects based on first-principles is of paramount importance in order to obtain an accurate description of molecular species in solution and in confined spaces.
Among others, the Frozen Density Embedding (FDE) density functional theory represents a embedding scheme in which environmental effects are included by considering explicitly the environmental system by means of its "frozen" electron density. In the present paper, (24) Huang, P.; Carter, E. A. Self-consistent embedding theory for locally correlated con-   Table S2: Dipole moment (components µ x , µ y , µ z and module |µ|) and dipole polarizability (tensor diagonal components α xx , α yy , α zz and isotropic contribution α iso ) of both the isolated (free) and embedded (emb.) water molecule. In the embedded water molecule, an ammonia molecule is used as environment. Data have been obtained using our new pyberthaemb implementation increasing the speed of light c = 10 · 137.036 au.

Fitting basis set
The exponents generation depends on the smallest and largest exponent in the primitive Gaussian exponents of the chosen basis set and by a parameter (N=2, 3,4) which together determine the total number of exponents. We recall that the auxiliary functions we employed are grouped in s, sp, spd, spdf and spdfg sets and that the exponents are shared within each of these sets. For instance, the A4 spdf g (generated with N=4) and associated with cc-pvtz basis for the oxygen atom have 17 different exponents with the angular part described with (5,7,5) auxiliary function notation adopted by Calaminici et al. S1 It describes five s sets together with 5 functions, seven spd sets together with 70 functions, and 5 spdfg sets together with functions 175. For the hydrogen atom the A4 spdf g fitting basis set correspond to a (2,4,3) fitting basis set. The total number of auxiliary function for the water molecule is 544. Using the same automatic generation scheme (even if may be not optimal) we generate different sets of exponents of reduced size depending on the teger parameter N. Thus, for the same cc-pvdz basis, we generate other five fitting basis sets (A2 s , A2 sp , A2 spd , A2 spdf g and A3 spdf g ).
See the definition in the following.    5 Evaluation of the spherical average of the EMBP As in the case of the the contour plot, reported in the main text, it is important to underlined as the spherical average of the EMBP for the neutral endohedral fullerenes A@C 60 (A=Rn,Og,Fl,Cn) is the result of a Nearest-neighbor interpolation performed starting from the potential represented on the original non homogeneous ADF grid. Initially we took advantage of the capability of PyADF to dump the Embedding Potential into a file. Subsequently the file is processed via a simple Python script (i.e. gridtodx.py S2 ) that, using a Nearest neighbor interpolation, as implemented in the SciPy, S3 transforms the potential, originally represented on a non homogeneous grid to a potential evaluated on a homogeneous grid, producing a DX file via the gridDataFormats package. S4 Once the EMBP has been dumped in a DX file, using an homogeneously spaced grid, quite easily the Spherical overage can be evaluated (i.e. avgalongr.py S2 ) simply considering the average of the Em-  Figure S2: Differences in orbital energies with respect to the isolated Rn atom for the Rn-C 60 system using the frozen density embedding (FDE) scheme (where Rn is used as active system while C 60 as frozen environment) and the supramolecular calculation (Supramolecular). All calculations have been carried out with the ADF code (2016.104 version) S5 using the BLYP functional together with the TZP basis set, the ZORA Hamiltonian for the inclusion of scalar relativistic effects, no frozen core approximation and the quality of both the density fitting and the numerical integration set to "Very good". The FDE calculations have been carried out using a full supramolecular grid for the frozen subsystem.