21.3 Frozen Density Embedding calculations

21.3.1 Background Theory

In the subsystem formulation of the density-functional theory a large system is decomposed into several constituting fragments that are treated individually. This approach offers the advantage of focusing the attention and computational cost on a limited portion of the whole system while including all the remaining environmental effects through an effective embedding potential. Here we refer in particular to the (fully-variational) Frozen Density Embedding (FDE)[367] with the Kohn-Sham Constrained Electron Density (KSCED) equations[368, 369].

In the FDE/KSCED method the embedding potential required by an embedded subsystem with density \(\rho_A\) to account for the presence of another (frozen) subsystem with density \(\rho_B\) is: \[\begin{equation} v_{\text{emb}}(\mathbf{r})=v^B_{\text{ext}}(\mathbf{r})+v_J[\rho_B](\mathbf{r})+ \frac{\delta {T}_s^{\text{nadd}}[\rho_A;\rho_B]}{\delta\rho_A(\mathbf{r})}+ \frac{\delta E_{xc}^{\text{nadd}}[\rho_A;\rho_B]}{\delta\rho_A(\mathbf{r})} \end{equation}\](21.1) where \(v^B_{\text{ext}}(\mathbf{r})\) and \(v_J[\rho_B](\mathbf{r})\) are the electrostatic potentials generated by the nuclei and electron density of the subsystem B, respectively, and \[\begin{equation} {T}_s^{\text{nadd}}[\rho_A;\rho_B]={T}_s[\rho_A+\rho_B]- {T}_s[\rho_A]-{T}_s[\rho_B] \, , \end{equation}\](21.2) \[\begin{equation} E_{xc}^{\text{nadd}}[\rho_A;\rho_B]=E_{xc}[\rho_A+\rho_B]-E_{xc}[\rho_A]- E_{xc}[\rho_B] \end{equation}\](21.3) are the non-additive non-interacting kinetic energy and exchange-correlation energy functionals, respectively. In the expressions above \({T}_s[\rho]\) is the (unknown) non-interacting kinetic energy density functional and \(E_{xc}[\rho]\) is the exchange-correlation energy functional. Note that, while the first two terms in Eq. (21.1) refer to classical electrostatics (and could be described by e.g. external point-charges), the last two terms are related to quantum-mechanical effects.

Using freeze-and-thaw [370] cycles, the role of the frozen and the embedded subsystem is iteratively exchanged, till convergence. If expressions (21.2) and (21.3) are computed exactly, then the density \(\rho_A+\rho_B\) will coincide with the exact density of the total system.

Because the FDE/KSCED was originally developed in the Kohn-Sham framework, using standard GGA approximations for \(E_{xc}[\rho]\), the non-additive exchange-correlation potential (\(\delta E_{xc}^{\text{nadd}}/\delta\rho_A(\mathbf{r})\)) can be computed exactly as a functional of the density, leaving the expression of the non-additive kinetic energy term as the only approximation (with respect to the corresponding GGA calculation of the total system), because the exact explicit density dependence of \(T_s\) from the density is not known. Using GGA approximations for the kinetic energy functional (\({T}_s\approx \tilde{T}_s^{GGA}\)) we have: \[\begin{equation} {T}_s^{\text{nadd}}[\rho_A;\rho_B] \approx \tilde{T}^{GGA}_s[\rho_A+\rho_B] -\tilde{T}^{GGA}_s[\rho_A] -\tilde{T}^{GGA}_s[\rho_B] \end{equation}\](21.4) and \[\begin{equation} \frac{\delta T^{\text{nadd}}[\rho_A;\rho_B]}{\delta\rho_A(\mathbf{r})}\approx \tilde{v}_T^{\text{GGA}}[\rho_A+\rho_B](\mathbf{r})-\tilde{v}_T^{\text{GGA}}[\rho_A](\mathbf{r}). \end{equation}\](21.5) where \(\tilde{v}_T^{\text{GGA}}(\mathbf{r})=\delta \tilde{T}_s^{\text{GGA}}/\delta\rho(\mathbf{r})\).

The FDE total energy of the total system is: \[\begin{align} E^{FDE}[\tilde{\rho}_A,\tilde{\rho}_B] &= T_s[\tilde{\rho}_A] + T_s[\tilde{\rho}_B] + {T}_s^{\text{nadd}}[\tilde{\rho}_A;\tilde{\rho}_B] \\ &+ V_{ext}[A+B] + J[\tilde{\rho}_A+\tilde{\rho}_B] + E_{xc}[\tilde{\rho}_A+\tilde{\rho}_B] \, . \end{align}\](21.6) Note that this energy differs from the KS total energy of the total system due to the approximation in Eq. (21.4) as well as the approximated kinetic potential (see Eq. 21.5) which lead to approximated embedded densities (\(\tilde{\rho}_A\approx\rho_A\) and \(\tilde{\rho}_B\approx\rho_B\)). With the current state-of-the-art GGA kinetic approximations, the error in the binding energy for weakly interacting systems is close to chemical accuracy.

Using the Generalized Kohn-Sham (GKS) theory, also hybrid exchange-correlation functionals can be used in embedding calculations. To obtain a practical computational method, the obtained embedding potential must be approximated by a local expression as shown in Ref. [371]. This corresponds to performing for each subsystem hybrid calculations including the interaction with other subsystems through an embedding potential derived at a semilocal level of theory. When orbital dependent exchange-correlation functionals (e.g. hybrid functional and LHF) are considered within the FDE method, the embedding potential includes a non-additive exchange-correlation term of the form \[\begin{equation} E^{\text{nadd}}_{xc}[\rho_A;\rho_B]=E_{xc}[\Phi^{A+B}[\rho_A+\rho_B]]- E_{xc}[\Phi^{A}[\rho_A]]-E_{xc}[\Phi^{B}[\rho_B]] \end{equation}\](21.7) where \(\Phi^{A+B}[\rho_A+\rho_B]\) denotes the Slater determinant which yields the total density \(\rho_A+\rho_B\). Since such a determinant is not easily available, the non-additive exchange-correlation contribution cannot be determined directly and the non-additive exchange-correlation term can be approximated as [372] \[\begin{equation} E^{\text{nadd}}_{xc}[\rho_A;\rho_B]\approx E^{\text{GGA}}_{xc}[\rho_A+\rho_B]- E^{\text{GGA}}_{xc}[\rho_A]-E^{\text{GGA}}_{xc}[\rho_B] \,. \end{equation}\](21.8)

The approximation of the kinetic potential within the embedding potential can lead to electron spill-out problems, especially for excitation energy calculations using diffuse basis functions. Here the energy of virtual orbitals is unphysically lowered, and therefore also the corresponding excitation energy, because of a poor description of the Pauli repulsion through the embedding potential. In Ref. [373] all-electron pseudopotentials were introduced at the environment atoms to model the missing Pauli repulsion contributions and additional to the embedding potential the potentials from the all-electron pseudopotentials are added to the vacuum Fock operator(\(F^{vac}\)): \[\begin{equation} F^{FDE(ECP)}_{\alpha\beta} = F^{vac}_{\alpha\beta} + \ensuremath{\langle \alpha | \nu_{emb}(\mathbf{r}) | \beta \rangle} + \nu_{ECP}~. \end{equation}\](21.9) Here the ECP potential \(\nu_{ECP}\) is nonlocal as it employs projection operators to act on specific parts of the one-electron Hilbert space.

21.3.2 Frozen Density Embedding calculations using the FDE script

The shell script FDE controls and executes automatically FDE calculations. The script FDE prepares the input files (running define/fdetools), runs the calculations and combines the results (running fdetools). Because the FDE equations are coupled sets of one-electron equations (one for each subsystem), full relaxation of the electron densities of both subsystems is obtained by using a freeze-and-thaw [370] procedure until convergence.

The converged FDE calculations are stored in the subdirectories STEPN/SUBSYSTEM_AAA, STEPN/SUBSYSTEM_AAB and so on, where N is the number of the FDE iterations. The subdirectories ISOLATED_SUBSYSTEM_AAA, ISOLATED_SUBSYSTEM_AAB and so on contain instead the calculations of the isolated subsystems (see also Section 21.3).

Currently two FDE implementations are available, which differ in the way the embedding potential is computed. One is the original implementation by Laricchia et. al. [371], which needs supermolecular steps in order to compute the embedding potential from matrix elements of the Fock operators. The other one was added by Treß et. al. [373] and it computes the embedding potential on a supermolecular integration grid and therefore avoids any supermolecular steps to compute the embedding potential. This implementation is the current default since it offers applicability for much bigger systems.

Current functionalities and limitations of the default FDE implementation are:

• only \(C_1\) point group;

• only for closed-shell systems that consist of closed-shell subsystems

• only total and binding energy calculations (no gradients);

• LDA/GGA kinetic energy functionals (for weakly interacting systems);

• full or pure electrostatic embedding;

• LDA/GGA, hybrid or orbital-dependent exchange-correlation potentials;

• multilevel FDE calculation;

• FDE calculation with frozen subsystems.

• FDE(ECP) method

• serial, OMP and MPI

• works with dscf,ridft,ricc2 and pnoccsd

In order to perform a FDE calculation, the files coord and control for the total system are necessary to take information on atomic coordinates and basis sets. The input file for the total system can be generated, as usual, with define but no calculation on the total system is required. $denconv 1.d-7 option should be defined in file control in order to better converge the embedded densities and better describe the dipole moment. The subsystems have to be defined with the $frag keyword, which can either be done by hand or with define in case of a HF dimer it may look like this:

 $frag
    fragment 1 atoms 1-2
    fragment 2 atoms 3-4

In addition to the keywords above also the size of the grid and the exchange-correlation functional for the calculation of the embedding potential have to be defined within the keyword $fde. It uses the same syntax as the dft keyword and can also be defined with define.

After these keywords are defined a FDE calculation can be started by simply invoking the command

                        FDE 

and an iterative resolution of the KSCED equations with revAPBEk [374, 375] as approximation of the non-additive kinetic potential (see Eq. 21.5) in the monomolecular basis is performed.

21.3.3 Options

Options for FDE calculations can either be specified as command lines or within the $fde keyword. Options specified with the command line can also be specified in file fde.input, which is read by the FDE script. If fde.input is not present it is created by the FDE script. Command lines options overwrites options found in the fde.input file.

Kinetic-energy functionals

In order to use different GGA approximations of the non-additive kinetic potential, the flag -k string can be used. Here string is the acronym used to identify a given GGA kinetic energy approximation, that can be selected among the following functionals:

• string=revapbek: generalized gradient approximation with a PBE-like enhancement factor, obtained using the asymptotic expansions of the semiclassical neutral atom as reference [374, 375] (revAPBEk). This is the default choice;

• string=lc94: Perdew-Wang (PW91) exchange functional reparametrized for kinetic energy by Lembarki and Chermette [376] (LC94);

• string=t-f: gradient expansion truncated at the zeroth order (GEA0), corresponding to the Thomas-Fermi functional.

For example, the command

                      FDE -k lc94 

approximates the non-additive kinetic contribution to the embedding potential through the functional derivative of LC94 kinetic energy functional.

A pure electrostatic embedding can be also performed with FDE script, where the embedding potential required by a subsystem A to account for the presence of the B one will be merely: \[\begin{equation} v_{\text{emb}}(\mathbf{r})=v^B_{\text{ext}}(\mathbf{r})+v_J[\rho_B](\mathbf{r}) \end{equation}\](21.10) with \(v^B_{\text{ext}}(\mathbf{r})\) and \(v_J[\rho_B](\mathbf{r})\) the electrostatic potentials generated respectively by the nuclei and electron density of the subsystem B. To perform an electrostatic embedding calculation use

                      FDE  -k electro   .

The electrostatic embedding is implemented only for testing purpose. It resembles an electrostatic embedding with external point-charges and/or point-dipoles, but it is “exact” as it is based on the whole densities (i.e. it considers all multipole moments of the density and the polarizabilities at all orders).

Equivalent command: --kin string

fde.input option: kin= string

The kinetic density functional can also be defined within the $fde keyword, for example the LC94 kinetic energy functional:

 $fde
  gridsize 3
  functional pbe
  kin lc94

The $fde keyword uses the same abbreviations for the kinetic energy functional as described above for the command line option. In addition to these kinetic functionals also kinetic functionals within the libxc or xcfun library can be specified together with the exchange-correlation functional:

 $fde
  gridsize 3
  functional libxc 101             #exchange pbe
  functional libxc add 2 130 521   #correlation pbe + LC94

FDE charged subsystems

FDE can perform calculations for charged closed-shell systems whose charge is localized on one or more subsystems. The charge of a subsystem is defined within the $frag keywords. For a sodium cation (atom 1) and water (atoms 2-4) system the $frag keyword looks like this:

 $frag
    fragment 1 charge 1 atoms 1
    fragment 2 charge 0 atoms 2-4

FDE(ECP) method

In order to prevent electron spill-out all-electron ECPs can be added to the embedding potential. This can be done by defining the ECPs used for the embedding with the keyword $ecpemb and adding fdeecp to the keyword $fde. Both can be done with define in the fde sub menu.

Post Hartree-Fock methods

Post Hartree Fock methods can be used in three different ways for all of them the keyword for the post Hartree Fock (HF) method has to be added first. The first option to use post HF methods is the perturbation to the energy (PE) approach, where contributions in the post HF method only result from the embedded Orbitals obtained from the underlying FDE-HF calculation. This approach can be applied by adding fdehf and kernel 0 to the $fde keyword:

 $fde
   gridsize 3
   functional pbe
   kin lc94
   fdehf
   kernel 0

The second option to use post HF methods is the post SCF reaction field method. Here correlation contributions are added to the reaction field after the FaT procedure is converged. This approach can be applied by adding fdehf and kernel 1 to the $fde keyword:

 $fde
   gridsize 3
   functional libxc 101             #exchange pbe
   functional libxc add 2 130 521   #correlation pbe + LC94
   fdehf
   kernel 1

As soon as the kernel contributions are added (kernel 1) within the post HF methods the kinetic energy functional has to be defined within the libxc or xcfun library. The third option is the perturbation to the energy and density (PTED) approach. Here the post HF method is included within the FaT cycles and the density from the post HF calculation is used to compute the embedding potential. This approach can be applied by only adding kernel 1 to the $fde keyword:

 $fde
   gridsize 3
   functional libxc 101             #exchange pbe
   functional libxc add 2 130 521   #correlation pbe + LC94
   kernel 1

FDE with frozen environment

FDE can perform embedding calculations where the environment of a subsystem is taken frozen, i.e. without scf calculation on it using an embedding potential. Therefore only one step will be performed if the flag --frozenenv will be used

                         FDE --frozenenv 1

Here the environment of subsystem 1 is taken frozen and the frozen embedding calculation is stored in the subdirectory STEP1/SUBSYSTEM_AAA.

Multilevel ansatz within FDE

A multilevel ansatz is defined such that all keywords within the control file of the supermolecular system will also be used for each subsystem and specific keywords for specific subsystems have to be defined in a multilevel reference file. An information about the file has to be added to the $fde keyword.

 $fde
   gridsize 3
   functional pbe
   multilevel file=multi.inp

Such a multilevel reference file (here: multi.inp) has to start with the keyword $multilevel and end with $end. Between those two keywords lines defining the fragments followed by their specific keywords have to be added. For a system containing of three subsystems where subsystem one and two shall be calculated with the DFT method and subsystem three with the CC2 model together with the first two excitation energies the multilevel reference file would look like this:

 
 $multilevel
    fragment 1-2
       $dft 
          gridsize 3
          functional pbe
    fragment 3
       $ricc2
         cc2
       $excitations
         irrep=a nexc=2 npre=4 nstart=6
 $end

It is suggested to define all necessary basis sets already in the control file of the supermolecular system to ensure a correct execution of the calculation.

Parallel calculations

If PARA_ARCH=SMP and OMP calculation will be performed. The flag -nth nthreads can be used to specify the number of threads. For example, with the following command

                     FDE -nth 4 

will use 4 threads.

Equivalent command: --nthreads integer

fde.input option: nthreads= integer

Restarting

The script FDE checks in the current directory for previous FDE calculations. If these are present, then the FDE calculation will be restarted from the last iteration found. The directories ISOLATED_SUBSYSTEM_AAA, ISOLATED_SUBSYSTEM_AAB and so on will be overwritten by the converged calculations from the previous run. The energy and the orbital from the isolated systems are saved in the current directory in the files: isolated_energy.ks, mos_AAA.ks, mos_AAB.ks and so on.

Note that a restart is possible only if the same subsystem definition is used. Other flags, e.g. kinetic and xc-functionals and convergence parameters, can be instead modified.

As all the options are saved in the fde.input file, to restart a FDE calculation the FDE script can be invoked without any parameters.

To force a calculation from scratch use:

                     FDE --scratch 

21.3.4 Compatibility Mode

The original implementation of Laricchia et. al. can be used by either adding comp to the $fde keyword or by adding –comp to the command line when invoking the FDE script. This implementation has some additional limitations and features:

• serial and OMP dscf runs (no MPI);

• monomolecular and supermolecular basis set approach;

• energy-decomposition;

• FDE(ECP) method is not available

Monomolecular and supermolecular basis set approach

The \(\rho_A\) and \(\rho_B\) densities can be expanded using the supermolecular or monomolecular basis set. In a supermolecular basis set expansion the basis functions \(\{\chi\}\) of both subsystems are employed to expand the subsystem electron densities. In a monomolecular basis set expansion, instead, only basis functions \(\{\chi^\ell\}\) centered on the atoms in the \(\ell\)-th subsystem are used to expand the corresponding density.

Both monomolecular and supermolecular basis set expansion of the electron densities are implemented in FDE: with the flag -m a monomolecular expansion is performed, while for a supermolecular one -s is used. In the absence of both flags a monomolecular expansion is performed by default.

For an accurate calculation of binding-energies of weakly interacting molecular systems a supermolecular basis set is required (to avoid the basis-set superposition error). Otherwise a very large monomolecular basis set is necessary.

NOTE: The FDE script supports only basis-set in the TURBOMOLE library.

Equivalent command: --mono or --super

fde.input option: method=mono or method=super

Convergence of the freeze-and-thaw cycles

The script FDE runs a self-consistent calculation when a convergence criterion is fulfilled. The convergence criterion is the change in the total dipole moment. This is a tight convergence criterion, as the dipole moment is highly sensitive to small changes in electron density. The convergence parameter \(\varepsilon^j\) for the \(j\)-th step in the freeze-and-thaw procedure is computed by means the following expression \[\begin{equation} \varepsilon^j=\frac{\lvert\Delta\mu^j_A\rvert+\lvert\Delta\mu^j_B\rvert}{2} \end{equation}\](21.11) where \[\begin{align} \lvert\Delta\mu^j_i\rvert&= \lvert\mu_i^j\rvert - \lvert\mu_i^{j-1}\rvert \qquad i=A,B \end{align}\] is the difference between the dipole moments of two consecutive steps for the \(i\)-th subsystem. Eq. (21.11) allows to consider changes in both subsystems or one of them because of the relaxation of their electron densities. By default, FDE stops when \(\varepsilon^j \le 0.005\) a.u.. The default value for the convergence criteria can be changed using the flag --epsilon= real where real is a decimal number.

The maximum number of freeze-and-thaw cycles can be specified by --max-iter= integer, and the default value is 20.

In order to make easy the convergence of the iterative solution of the KSCED coupled equations, a damping factor \(\eta\) must be used for the matrix elements of the embedding potential \(\bigl(v_{\text{emb}}\bigr)_{ij}\) as perturbation to a given subsystem \[\begin{equation} {}_d\bigl(v_{\text{emb}}\bigr)_{ik}^j= (1-\eta) \bigl(v_{\text{emb}}\bigr)_{ik}^j + \eta \bigl(v_{\text{emb}}\bigr)_{ik}^{j-1} \end{equation}\](21.12) for the \(j\)-th iteration. Here \({}_d\bigl(v_{\text{emb}}\bigr)_{ik}^j\) is the matrix element effectively used in the \(j\)-th iteration after the damping. In FDE the starting value of \(\eta\) can be changed using --start-damp= real (default value is \(0.45\)) where real is a decimal number. The damping parameter can also dynamically change at each iterative step (according to the convergence process) of a quantity set by --step-damp= real (default value is \(0.10\)). The minimum value set by --max-damp= real (default value is \(0.90\)).

fde.input options:
epsilon= real
max-iter= integer
start-damp= real
max-damp= real
step-damp= real

Embedding energy error

The embedding error in the total energy is computed as \[\begin{equation} \Delta E=E^{\text{FDE}}[\tilde{\rho}_A;\tilde{\rho}_B]-E^{\text{DFT}}[\rho] \end{equation}\](21.13) where \(E^{\text{DFT}}\) is the DFT total energy of total system with density \(\rho(\mathbf{r})\). In order to compute \(\Delta E\) as well as its components, the flag --err-energy must be used. This flag will required also the DFT calculation on the total system. In this case the converged SCF output file must be named output.dscf.

An example of session output for the computation of embedding energy and energy error decomposition, when --err-energy flag is present, is the following:

FDE ENERGY (TOTAL SYSTEM):   -200.99720391651 Ha
FDE BINDING ENERGY:                      4.960885 mHa
                                         3.113002 kcal/mol
FDE ENERGY ERROR:                        2.003352 mHa
ERROR ENERGY DECOMPOSITION
  coulomb contribution:                 -0.693026 mHa
  nuclear contribution:                 -3.136544 mHa
  exchange-correlation contribution:    -1.156390 mHa
  kinetic contribution:                  6.989320 mHa

where the FDE energy (\(E^{\text{FDE}}\)), the FDE binding energy, the embedding energy error (\(\Delta E\)) and the error energy decomposition in its coulomb, nuclear, exchange-correlation and kinetic contributions are reported. This output is present at each FDE iteration.

fde.input option: err-energy=1

FDE with hybrid and orbital-dependent functionals

In order to use local approximations (21.1) and (21.8) with FDE, the flag -f string must be add to the options of the script. Here string denotes the local/semilocal approximation to hybrid or orbital-dependent exchange-correlation potentials in \(v_{\text{emb}}(\mathbf{r})\). All LDA/GGA functionals in TURBOMOLE can be considered as approximations.

For example, the command

                         FDE -f b-lyp --comp

can be used to approximate bh-lyp or b3-lyp hybrid non-additive potentials, while the command

                         FDE -f pbe --comp

approximates the pbe0 hybrid non-additive potentials. Other combinations of functionals are not recommended (meta-GGA are not supported).

Finally, also calculations with the Local Hartree-Fock (LHF) potential can be performed. In this case the command

                     FDE -f becke-exchange --comp 

can be used to approximate the LHF non-additive potential[372].

Equivalent command: --func string

fde.input option: func= string