17.3 Frozen Density Embedding calculations

17.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 enviromental effects through an effective embedding potential. Here we refer in particular to the (fully-variational) Frozen Density Embedding (FDE) [196] with the Kohn-Sham Constrained Electron Density (KSCED) equations [197,198].

In the FDE/KSCED method the embedding potential required by an embedded subsystem with density ρA to account for the presence of another (frozen) subsystem with density ρB is:

                          δT nadd[ρA;ρB]  δEnadd[ρA;ρB]
vemb(r) = vBext(r)+ vJ[ρB](r)+--s-δρ-(r)----+ ---xcδρ-(r)----
                                A              A

where vextB(r) and vJ[ρB](r) are the electrostatic potentials generated by the nuclei and electron density of the subsystem B, respectively, and

Ts   [ρA;ρB ] = Ts[ρA + ρB]- Ts[ρA ]- Ts[ρB],

Enadd[ρ ;ρ ] = E  [ρ  + ρ ]- E  [ρ  ]- E [ρ ]
 xc   A  B     xc A    B    xc A     xc B

are the non-additive non-interacting kinetic energy and exchange-correlation energy functionals, respectively. In the expressions above Ts[ρ] is the (unknown) non-interacting kinetic energy density functional and Exc[ρ] is the exchange-correlation energy functional. Note that, while the first two terms in Eq. (17.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 [199] cycles, the role of the frozen and the embedded subsystem is iteratively exchanged, till convergence. If expressions (17.2) and (17.3) are computed exactly, then the density ρA + ρ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 Exc[ρ], the non-additive exchange-correlation potential (δExcnadd∕δρA(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 the corresponding GGA calculation of the total system), because the exact explicit density dependence of Ts from the density is not known. Using GGA approximations for the kinetic energy functional (Ts ˜TsGGA) we have:

 nadd          GGA            GGA        GGA
Ts   [ρA;ρB] ≈ T˜s [ρA + ρB]- ˜Ts   [ρA]- T˜s   [ρB]


δT-nadd[ρA;ρB]   GGA              GGA
    δρA(r)    ≈ ˜vT   [ρA + ρB](r) -v˜T   [ρA](r).

where TGGA(r) = δT˜sGGA∕δρ(r).

The FDE total energy of total system is:

EFDE [˜ρA, ˜ρB] =   Ts[˜ρA]+ Ts[˜ρB]+ Tnsadd[˜ρA;˜ρB ]

             +   Vext[A+ B ]+ J[˜ρA + ˜ρB]+ Exc[˜ρA + ρ˜B ].     (17.6)
Note that this energy differs from the KS total energy of the total system due to the approximation in Eq. (17.4) as well as the approximated kinetic potential (see Eq. 17.5) which lead to approximated embedded densities (˜ρ A ρA and ρ˜ B ρ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. [200]. 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

Enxacdd[ρA;ρB ] = Exc[ΦA+B [ρA + ρB ]]- Exc[ΦA[ρA]]- Exc[ΦB [ρB]]

where ΦA+B[ρA + ρB] denotes the Slater determinant which yields the total density ρA + ρ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 [201]

  nadd           GGA            GGA        GGA
E xc [ρA;ρB] ≈ E xc [ρA + ρB ]- E xc [ρA]- E xc [ρB].

17.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), runs the calculations (only dscf is supported in the present versiob), 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 [199] procedure until convergence.
The converged FDE calculations are store in the subdirectories STEPN/SUBSYSTEM_A and STEPN/SUBSYSTEM_B, where N is the number of the FDE iteration. The subdirectory ISOLATED_SUBSYSTEM_A and ISOLATED_SUBSYSTEM_B contain instead the calculations for isolated subsystems (see also Section 17.3.3).

Current functionalities and limitations of FDE are:

In order to perform a FDE calculation, the files coord and control for the total system are necessary to take informations 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.

Given a closed-shell supramolecular system with a GGA/LDA exchange-correlation functional, the command

                        FDE -p 3

invokes an iterative resolution of the KSCED equations with revAPBEk [202,203] as approximation of the non-additive kinetic potential (see Eq. 17.5) in the monomolecular basis set approach. The two subsystems are defined via an integer m = 3 in the example above which identifies the first atom of the subsystem B in the file coord of the supramolecular system with n atoms, where the atoms 1m - 1 belong to the subsystem A while the atoms mn to the B one. Thus the file coord must contains first all the atoms of the system A and then all the atoms of the system B.

As an example we report here the FDE -p 3 output for the HF dimer:

                           FDE  Version 1.02  
                 Frozen Density Embedding Main Driver  
                          Scf-like procedure for  
                closed-shell interacting systems (dimers)  
                 program development: Savio Laricchia  
                                      Eduardo Fabiano  
                                      Fabio Della Sala  
          S. Laricchia, E. Fabiano, L. A. Constantin, F. Della Sala,  
          J. Chem. Theory Comp. (2011)  
          S. Laricchia, E. Fabiano, F. Della Sala,  
          J. Chem. Phys. 133, 164111 (2010)  
          L. A. Constantin, E. Fabiano, S. Laricchia, F. Della Sala,  
          Phys. Rev. Lett. 106, 186406 (2011)  
          S. Laricchia, E. Fabiano, F. Della Sala,  
          Chem. Phys. Lett. 518, 114 (2011)  
                   Sun Mar 25 23:00:01 CEST 2012  
Monomolecular basis set approach...  
Serial calculation will be performed...  
running /home/fabiods/REDO/branch64/TURBOMOLE/bin/em64t-unknown-linux-gnu/dscf  
b-lyp exchange-correlation potential in KS supermolecular calculation...  
revapbek kinetic energy approximation will be used...  
Default convergence criterion on the system dipole: 0.005  
Default value of starting damping parameter is 0.45  
Default value of step damping parameter is 0.10  
Default value of maximum damping parameter is 0.90  
Default value of maximum fde iterations  is 20  
Saving options in fde.input  
  | Subsystem A atomic coordinates and basis set information  |  
  |  x        y       z      atom    basis set         ecp    |  
  2.5015  -0.1705  -0.0000     f    def2-TZVP         none  
  3.2889   1.3859   0.0000     h    def2-TZVP         none  
  | Subsystem B atomic coordinates and basis set information  |  
  |  x        y       z      atom    basis set         ecp    |  
 -2.7537   0.0364  -0.0000     f    def2-TZVP         none  
 -1.0191  -0.1789   0.0003     h    def2-TZVP         none  
Running Isolated subsystems:  
Saved isolated subsystems data in:  
*   FDE - step 1    *  
FDE ENERGY (TOTAL SYSTEM):            -200.96417090754 Ha  
FDE BINDING ENERGY:                      5.865327 mHa  
                                         3.680548 kcal/mol  
Dipole convergence: 0.138071,  Damping: 0.45  
*   FDE - step 2    *  
FDE ENERGY (TOTAL SYSTEM):            -200.96418098234 Ha  
FDE BINDING ENERGY:                      5.875401 mHa  
                                         3.686870 kcal/mol  
Dipole convergence: 0.009246,  Damping: 0.35  
*   FDE - step 3    *  
FDE ENERGY (TOTAL SYSTEM):            -200.96418289036 Ha  
FDE BINDING ENERGY:                      5.877309 mHa  
                                         3.688067 kcal/mol  
Dipole convergence: 0.004395,  Damping: 0.25  
 See embedded susbsystems calculations in:  
 See total system in:  
                   Sun Mar 25 23:00:21 CEST 2012  
                   Total time: 20 secs.

The final energies are stored in the file fde_energy. The directory STEPN/ENERGY_SYSTEM contains the total system with density ρA + ρB; this directory can (only) be used for density analysis.

17.3.3 Options

All the options for the FDE can be specified as commandlines, and are described below. The options can be 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.

Subsystem definition

The flag -p integer is required or it must be present in the fde.input file.

Equivalent command: --pos-cut integer

fde.input option: pos-cut= integer

Kinetic-energy functionals

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

For example, the command

                      FDE -p 3 -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:

vemb(r) = vext(r)+ vJ[ρB](r)

with vextB(r) and vJ[ρB](r) the electrostatic potentials generated respectively by the nuclei and electron density of the subsystem B. To perform an electrostatic embedding calculation use

                      FDE  -p 3 -k electro

and can be performed for both Kohn-Sham (only for LDA/GGA exchange-correlation functionals) and Hartree-Fock methods.

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 polarizabilies at all orders).

Equivalent command: --kin string

fde.input option: kin= string

FDE charged subsystems

FDE can perform calculations for charged closed-shell systems whose charge is localized on one or both subsystems. To localize the charge on a given subsystem, --chargeA= integer must be used for the subsystem A and --chargeB= integer for the B one. Here integer denotes the charge added to the neutral subsystem. For example, the command

                      FDE -p 3 --chargeA= 2

performs a FDE calculation for a negative charged closed shell system (for example Zn( H2O)22+) whose subsystem B has charge 2. Note that in this case the starting control file must have a charge +2.

fde.input option: chargeA= integer , chargeB= integer

FDE with subsystem B taken frozen

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

                         FDE -p 3 --frozen

The frozen embedding calculation is store in the subdirectory STEP1/SUBSYSTEM_A. The control file is modified with the following keywords:

$fde read  
$fde-input zj file=fde_ZJ.mat  
$fde-input kxc file=fde_KXC.mat

The program dscf will read the submatrices fde_ZJ.mat and fde_KXC.mat and add them to the Hamiltonian.

fde.input option: frozen=1

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 -p 3 -nth 4

will use 4 threads.

Equivalent command: --nthreads integer

fde.input option: nthreads= integer

Monomolecular and supermolecular basis set approach

The ρA and ρB densities can be expanded using the supermolecular or monomolecular basis set. In a supermolecular basis set expansion the basis functions {χ} of both subsystems are employed to expand the subsystem electron densities. In a monomolecular basis set expansion, instead, only basis functions {χ} centered on the atoms in the -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 εj for the j-th step in the freeze-and-thaw procedure is computed by means the following expression

     |Δ μj|+ |Δ μj |
εj = ---A------B--


|Δμji|  =  |μji|- |μji-1|  i = A, B
is the difference between the dipole moments of two consecutive steps for the i-th subsystem. Eq. (17.10) allows to consider changes in both subsystems or one of them because of the relaxation of their electron densities. By default, FDE stops when εj 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 η must be used for the matrix elements of the embedding potential (vemb)ij as perturbation to a given subsystem

d(vemb)j = (1 - η)(vemb)j + η(vemb)j-1
       ik              ik        ik

for the j-th iteration. Here d(vemb)ikj is the matrix element effectively used in the j-th iteration after the damping. In FDE the starting value of η 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

ΔE = EFDE [˜ρA;˜ρB]- EDFT [ρ]

where EDFT is the DFT total energy of total system with density ρ(r). In order to compute Δ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  
  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 (EFDE), the FDE binding energy, the embedding energy error (Δ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


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_A and ISOLATED_SUBSYSTEM_B will be overwritten by the converged calculations from previous run. The energy and the orbital from the isolated systems are saved in the current directory in the files: isolated_energy.ks, mos_A.ks and mos_B.ks.

Note that a restart is possible only if the same subsystem definition and the same basis set are used (e.g. the same -p flag and the -s or -m flag). Other flags, e.g. kinetic and xc-functionals and convergence paramters, 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 inkoved without any parameters.

To force a calculation from scratch use:

                     FDE -p 3 --scratch

Table 17.1: Other options in the shell script FDE

-d or --dipole dipole=1 dipole moment each step
-v or --verbose  shows more informations
--save-mos save-mos=1 save the MOs of both subsystems
for each step
--save-matrix save-matrix=1 save fde matrices
in each step
--help list all commands

17.3.4 FDE with hybrid and orbital-dependent functionals

In order to use local approximations (17.1) and (17.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 vemb(r). All LDA/GGA functionals in TURBOMOLE can be considered as approximations.

For example, the command

                         FDE -p 3 -f b-lyp

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

                         FDE -p 3 -f pbe

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 -p 3 -f becke-exchange

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

Equivalent command: --func string

fde.input option: func= string