FDE: Frozen Density Embedding¶

Fragment-specific FDE options¶

For each frozen fragment, several additional options can be applied. To do this, the fragment specification is used as a subblock key by appending a & sign. The subblock is terminated with SubEnd. This subblock key looks, in the most general form, as follows:

FRAGMENTS
  ...
  FragType FragFile
  ...
END
FDEFRAGMENTS
  ...
  FragType type=FDE &
     {FDEOPTIONS [USEBASIS] [RELAX or FREEZEANDTHAW] [OPTIMIZE]}
     {FDEDENSTYPE [SCF | SCFexact | SCFfitted ]}
     {RELAXCYCLES n or FREEZEANDTHAWCYCLES n}
     {XC [LDA | GGA ggapotx ggapotc | MODEL SAOP]}
  SubEnd
  ...
END

FDEOPTIONS

FDEOPTIONS USEBASIS

If the USEBASIS option is specified, the basis functions of this frozen fragment will be included in the calculation of the embedded subsystem. This allows to expand the density of the embedded subsystem using not only atom-centered basis sets localized in the embedded subsystem but also the ones in the environment Ref. [13]. In large-scale simulations using the embedding potential, this option is recommended to be used in the preparation stage to investigate the basis set dependence of the results (chapter 5.3 in Ref. [3]). This option is also an indispensable element in the procedure introduced in Ref. [13] to test approximants to the kinetic-energy component of the embedding potential introduced by Wesolowski and Warshel.

FDEOPTIONS RELAX or FREEZEANDTHAW

If the RELAX option (or equivalent FREEZEANDTHAW option) is specified, the density of this frozen fragment will be relaxed in freeze-and-thaw cycles [5], i.e., the embedded subsystem is frozen, while this fragment is thawed. This is repeated, until convergence is reached or until the maximum number of iterations has been performed. By relaxing frozen fragments, it is possible to improve a given approximate environment density by including the polarization of the environment due to the embedded system. This option is recommended to be used in the preparation stage of a large-scale numerical simulation. The freeze-and-thaw calculations lead to a pair of electron densities (embedded system and environment) that minimizes the total energy. As a consequence, the electron density of the environment derived from the freeze-and-thaw calculations can be used as a reference to verify the adequacy of the assumed electron density for the environment in a large-scale simulation. Due to technical restrictions, freeze-and-thaw is not possible if an open-shell (unrestricted) fragment is present.

FDEOPTIONS USEBASIS RELAX or FDEOPTIONS USEBASIS FREEZEANDTHAW

It is further possible to combine USEBASIS and RELAX or FREEZEANDTHAW. In this case, the basis functions of the non-frozen fragment will be included when the density of the fragment is relaxed. This allows fully relaxed calculations with supermolecular expansion of the electron density of each subsystem. This option is to be used to test approximants to the kinetic-energy component of the embedding potential introduced by Wesolowski and Warshel by means of the procedure introduced in ref. [13].

FDEDENSTYPE

The FDEDENSTYPE option can be used to specify which density is read from the fragment file. The possible options are:

FDEDENSTYPE SCF (or FDEDENSTYPE SCFexact)

The exact density (not calculated using the fit functions) is used. This is the default.

FDEDENSTYPE SCFfitted

The fitted density is used. This is less accurate but can be significantly faster.

RELAXCYCLES n or FREEZEANDTHAWCYCLES n

This gives the maximum number of freeze-and-thaw cycles that are performed for this fragment. If the maximum number given in the FDE block is smaller, or if convergence is reached earlier, then fewer cycles are performed. For historical reasons, two equivalent keywords are available.

XC

The XC option can be used to select the exchange-correlation potential that is used for this fragment when it is relaxed. By default, the same potential as for the non-frozen system is used, but in some cases it might be preferable to use another approximation for certain fragments. An example is given in Ref. [12].

XC LDA

This option selects LDA as exchange-correlation potential for relaxing this fragment.

XC GGA ggapotx ggapotc

This selects a GGA potential for relaxing this fragment. The GGA potential is specified by giving the name of the exchange potential, followed by the name of the correlation potential. The available potentials are listed in the documentation for the XC key.

XC MODEL SAOP

This selects the model potential SAOP for relaxing this fragment.

Kinetic energy approximants¶

The approximants to the kinetic energy dependent component of the embedding potential are described here.

FDE
  {approximants to the kinetic energy dependent
                component of the embedding potential}
  {CJCORR [rho_cutoff]}
  {GGAPOTXFD exchange approximant}
  {GGAPOTCFD correlation approximant}
end

approximants to the kinetic energy dependent component of the embedding potential

Several approximants to the kinetic-energy-dependent component of the effective potential given in Eq. (21) of ref. [1] are available. None of them is applicable if the embedded system is covalently bound to its environment. The user is recommended to look at the numerical value of the TSNAD(LDA) parameter which is given in the units of energy and can be considered as a measure of the overlap. The following rule of thumb should be applied: if this parameter is smaller than the estimated interaction energy between the embedded subsystem and the environment, then the available approximants are most probably adequate. If it exceeds this limit, the results can be less reliable. Printing TSNAD(LDA) is not done by default, as it can be quite time-consuming. Its printing is switched on by including “EXTPRINTENERGY”, and “PRINTRHO2”, and “FULLGRID” in the FDE input block. If no kinetic energy approximant is specified, by default the local-density approximation (Thomas-Fermi approximant) is used. For an assessment of approximants for weakly overlapping pairs of densities see Refs. [13] [16] [15] [17]. Based on these studies, the use of PW91k (= GGA97) is recommended.

APPROXIMANTS TO BE USED IN NORMAL APPLICATIONS

THOMASFERMI (default)

Local-density-approximation form of vt[rhoA,rhoB] [20] derived from Thomas-Fermi expression for Ts[rho] [18] [19].

GGA97 (or PW91K)

Generalized-gradient-approximation form of vt[rhoA,rhoB] [16] derived from the Lembarki-Chermette [21] approximant to Ts[rho]. This approximant is currently the recommended one based on the numerical analysis of its accuracy [15] [16] and the fact that the used enhancement factor disappears at large reduced density gradients, i.e. where the second-order gradient-expansion approximation fails [13] [17].

NDSD

Similarly to GGA97, the NDSD approximant is constructed by taking into account the asymptotic behavior of the functional vt[rhoA,rhoB] at small density gradients. In the construction of NDSD, the exact property of vt[rhoA,rhoB] at rho_A \(\rightarrow\) 0 and for \(\int\) rhoB = 2 given in Eq. A6 of Ref. [22] is also taken into account. The analysis of the accuracy of this potential [22] shows that NDSD is of the same or superior quality as GGA97. NDSD is, therefore, recommended as the successor of GGA97 to be used anywhere where the quality of the results depends directly on the accuracy of the potential vt[rhoA,rhoB], i.e., for obtaining electronic-structure-dependent properties. The analytical form of the corresponding approximant to the functional Ts^nad [rho_A,rho_B]$ exists (Eq. 23 in Ref. [22]). It is not possible, however, to obtain the analytical form of the corresponding parent functional for the kinetic energy Ts[rho]. To reflect this and the fact that, similarly to the GGA approximants to vt[rhoA,rhoB], the numerical values of only first- and second derivatives of density are needed, the label NDSD (Non-Decomposable Second Derivatives) is used.

OBSOLETE APPROXIMANTS (can be used but GGA97 leads usually to a better embedding potential see refs. [13] [16])

LLP91

Generalized-gradient-approximation form of vt[rhoA,rhoB] [13] derived from Lee-Lee-Parr [23] approximant to Ts[rho].

PW86k

Generalized-gradient-approximation form of vt[rhoA,rhoB] [13] derived from the Fuentealba-Reyes approximant to Ts[rho] [24].

THAKKAR92

Generalized-gradient-approximation form of vt[rhoA,rhoB] [16] derived from the Thakkar approximant to Ts[rho] [25].

APPROXIMANTS WHICH MIGHT BE USEFUL ONLY FOR THEORY DEVELOPMENT

The accuracy of some of these approximants was investigated in detail [16] [13] [15] [17]. Each of them was shown to lead to a qualitatively incorrect embedding potential. They shouldn’t be used in practical applications.

COULOMB

Neglecting completely vt[rhoA,rhoB] (vt[rhoA,rhoB] equals zero) together with the exchange-correlation component of the embedding potential introduced by Wesolowski and Warshel.

TF9W

The approximant to vt[rhoA,rhoB] [1] derived from the second-order gradient expansion [24] for Ts[rho].

WEIZ

The approximant to vt[rhoA,rhoB] [17] derived from the von Weizsäcker approximant to Ts[rho] [6].

OL91A

Generalized-gradient-approximation form of vt[rhoA,rhoB] [13] derived from the first Ou-Yang and Levy approximant to Ts[rho] [26].

OL91B

Generalized-gradient-approximation form of vt[rhoA,rhoB] [16] derived from the second Ou-Yang and Levy approximant to Ts[rho] [26].

E00

Generalized-gradient-approximation form of vt[rhoA,rhoB] [27] derived from a kinetic energy functional by Ernzerhof [28] which represents the gradient expansion approximation up to the fourth order.

P92

Generalized-gradient-approximation form of vt[rhoA,rhoB] [27] derived from a kinetic energy functional by Perdew [29] which represents the gradient expansion approximation up to the sixth order.

LONG DISTANCE CORRECTIONS TO THE EFFECTIVE POTENTIAL

CJCORR

Option to switch on a long-distance correction. By default this option is not used. As was shown in Ref. [30], with the available approximate kinetic-energy approximants, the embedding potential has the wrong form in the limit of a large separation of the subsystems. In particular, it was shown that this can have serious consequences in the case of “supermolecular expansion of electron density of each subsystem” calculations (USEBASIS option). In Ref. [30], a correction is proposed that enforces the correct long-distance limit. (See also this reference for limitations of this correction.)

CJCORR [rho_cutoff]

This option switches on the long-distance correction. This option has to be used in combination with one of the above kinetic-energy approximants. By default, a density cut-off of 0.1 is employed.

NONADDITIVE EXCHANGE-CORRELATION APPROXIMANT

GGAPOTXFDGGAPOTCFD

Option to specify the non-additive exchange-correlation approximant. By default, in the construction of the effective embedding potential the exchange-correlation approximant that was specified in the XC block is used. It is possible to specify a different approximant with the GGAPOTXFD and GGAPOTCFD options. This is particularly useful in combination with the use of model potentials like SAOP, that can not be used in the embedding potential because of their orbital dependence. (For a discussion, see Ref. [12].)

GGAPOTXFD exchange approximant

The exchange approximant is used in the construction of the embedding potential. The same exchange approximants as in the XC key are available.

GGAPOTCFD correlation approximant

The correlation approximant is used in the construction of the embedding potential. The same correlation approximants as in the XC key are available.

General FDE options¶

In addition to the fragment-specific options and the kinetic energy approximants, there are also a number of options available in FDE calculations that will be described in the following.

FDE
  {FULLGRID}
  {RELAXCYCLES n or FREEZEANDTHAWCYCLES n}
  {RELAXPOSTSCF or FREEZEANDTHAWPOSTSCF}
  {EXTPRINTENERGY}
  {PRINTRHO2}
  {ENERGY}
  {SDFTENERGY}
  {DIPOLE}
end

FULLGRID

By default, FULLGRID is not used, and in FDE calculations the integration grid is generated as described in Ref. [2] by including only atoms of the frozen subsystem that are close to the embedded subsystem in the generation of the integration grid. The distance cutoff used is chosen automatically, based on the extent of the basis functions of the embedded subsystem. (It can also be chosen manually, see the option qpnear in the INTEGRATION key) This scheme results in an efficient and accurate integration grid. However, it is possible that the default integration scheme is not accurate enough. This can be the case for weakly interacting systems and when the distance between the frozen and the embedded system is large. It is therefore recommended to check the quality of the default integration grid by comparing to results obtained using the full supermolecular grid (FULLGRID option).

If the subkey FULLGRID is included, all atoms of the frozen system are included in the generation of the integration grid. This results in the same grid that would be used in a supermolecular calculation of the combined frozen and embedded system. The integration grid generated by this option might be much larger than the default grid. This option should be used to check the quality of the default integration grid.

RELAXCYCLES n or FREEZEANDTHAWCYCLES n

Specifies the maximum number n of freeze-and-thaw iterations [5] that are performed (for frozen fragments with the RELAX) option. If a smaller number of iterations is specified as a fragment-specific option, for this fragment this smaller number is used. Furthermore, if convergence is reached earlier, no more iterations will be performed.

RELAXPOSTSCF or FREEZEANDTHAWPOSTSCF

If this option is included, several post-SCF properties will be calculated after each freeze-and-thaw cycle [Ref. [5]. These are otherwise only calculated in the last cycle.

EXTPRINTENERGYPRINTRHO2

If the options EXTPRINTENERGY and PRINTRHO2 are included (both are needed and should be listed on separate lines), several additional quantities will be printed, including TSNAD(LDA). In order to obtain meaningful numbers, also the FULLGRID keyword (see above) has to be used.

ENERGY

Option to switch on the calculation of the FDE energy as the sum of the energy E[rhoA] of the active, embedded system and the interaction energy Eint[rhoA,rhoB] of the embedded system with the frozen environment. This relies on the calculation of the total energy for the embedded system and all caveats and restrictions for total energy evaluations apply (see keyword TOTALENERGY). All energy contributions are evaluated on the grid of the active subsystem. Some contributions to the interaction energy Eint[rhoA,rhoB] require an accurate integration grid in the region of the environment. Thus, in pure embedding calculations (without fragment-specific option RELAX), an accurate calculation of the FDE energy requires a full supermolecular integration grid (FULLGRID option). Details on the implementation and the performance of kinetic energy functionals for interaction energies are documented in Ref. [27]

The calculation of the full, variationally minimized subsystem DFT energy, that is, the sum of the energy of two subsystems E[rhoA] and E[rhoB] and their interaction energy Eint[rhoA,rhoB] in the framework of FDE, is invoked if then the fragment densities are relaxed in freeze-and-thaw cycles (option RELAXCYCLES and fragment-specific FDE option RELAX). In this case the supermolecular integration grid is not required. Instead, in each step of the freeze-and-thaw cycle, the critical energy terms are taken from the previous freeze-and-thaw step of the presently frozen fragment. The convergence of the energy contributions with the number of freeze-and-thaw iterations should be carefully monitored. Due to conceptual problems for the evaluation of the non-additive kinetic energy contribution, only two subsystems, that is, one frozen fragment, is supported for FDE energy calculations with freeze-and-thaw.

SDFTENERGY

This is a generalization of the original ENERGY implementation that allows for the evaluation of the FDE/sDFT total energy. Usage of a supermolecular integration grid is recommended (FULLGRID option), because of the non-additive kinetic and XC contributions. EXTPRINTENERGY can be combined with this option to print the subsystem DFT energy after each freeze-and-thaw cycle. For pure embedding calculations, the SDFTENERGY implementation reads the total energy from each frozen fragment’s adf.rkf (TAPE21) file. The total energy of frozen fragments is not (re-)evaluated during FDE calculations. It is therefore mandatory to add the TOTALENERGY keyword to the preparatory isolated fragment calculations.

DIPOLE

The dipole of the supersystem is calculated as the sum of analytically integrated fragment dipole moments. Note that the dipole moment becomes origin dependent for charged (sub-)fragments. The current implementation does NOT take care of this.

Frozen Density Embedding with External Orthogonality¶

An implementation of External Orthogonality (EO) into the FDE framework in ADF can be found in Ref. [31]. In this method in the fragment calculation one need to specify ghost atoms at the positions of the atoms in all the other fragments. Thus each fragment is calculated in the supermolecular basis.

AOMat2File
IgnoreOverlap
FDE
  EXTERNALORTHO {factor}
  ...
END

EXTERNALORTHO

Used to specify the use of external orthogonality (EO) in the FDE block, with an optional factor (default factor is 1e6). Note that the overlap region is defined by ghost atoms, and the general keyword IGNOREOVERLAP is necessary. Additionally, the general keyword AOMAT2FILE is required to save some important fragment AO matrices to adf.rkf for use in an EO calculation. it is recommended that one uses the STOFIT keyword for the STOFIT density fitting method.

FDE and (localized) COSMO¶

COSMO solvation can be included in combination with FDE. This means that a COSMO cavity will be created that holds both the active and frozen fragments.

For very large systems the solution of the COSMO equations can become the bottleneck of the calculation. A local COSMO variant, which exploits the subsystem nature of the underlying electronic description, is implemented in ADF, see Ref. [32]. This method, called LoCOSMO, is an approximation to regular COSMO for fragment-based electronic structure methods. If a given fragment is active (i.e. only its own density changes over the SCF cycles), the surface charges outside a given radius around this fragment will be kept fixed in magnitude (excluded from the COSMO optimization). However, no interactions are ignored by adding a corresponding constant term to the COSMO solution vector.

SOLVATION
  ..
  CHARGED LoCosmo LoCosmoDist
END

LoCosmo LoCosmoDist

If LoCosmo is included the local COSMO will be used in the calculation. LoCosmoDist is a cutoff radius (Angstrom), which must be specified by the user, if LoCOSMO is included. All surface charges within this radius from any atom of the active fragment are included in the active charge space. A sensible value for LoCosmoDist is 5.0 Angstrom.

To be really useful, Freeze-and-Thaw cycles should be carried out with LoCOSMO. After each fragment calculation, the KF file LOCSURCH is written to disk. If present, this file will be automatically read in during the next (LoCOSMO) fragment calculation. This means that if the internal FT cycles are used by specifying RELAXCYCLES, the code creates and reads in the file in a fully automatic fashion. LOCSURCH is neither created nor read during a regular COSMO calculation.

In case one calculates excitation energies with the EXCITATIONS block key, and LoCOSMO is specified, only the solvent response due to surface charges within the specified cutoff radius will be taken into account, see also Ref. [33]. It should be noted that in a regular (uncoupled) FDE TDDFT calculation, the response due to the frozen density is neglected.

FDE and QM/FQ(Fμ)¶

FDE can also be combined with the QM/FQ(Fμ) method which can be used to embed a quantum system in a polarizable classical environment. The fluctuating classical layer can be included in all stages of an FDE calculation, i.e. when performing the separate calculations for the two isolated fragments as well as the full FDE system. There is also the possibility of using the QM/FQ(Fμ) method to account for the mutual response of the subsystems by endowing the frozen layer’s atoms with fluctuating charges (and dipoles), provided a suitable parametrization is available. See the QM/FQ(Fμ) manual page for more details.

Subsystem TDDFT, coupled FDE¶

The linear-response subsystem TDDFT code implements the theory of coupled excited states for subsystems as described in Refs. [35] [34]. This theory is based on the FDE extension to excited states [6], which is implemented in ADF in a local response approximation, i.e., neglecting the dynamic response of the environment [7].

The subsystem TDDFT code allows to treat the mutual response of several subsystems, including the ones that are considered environment. A more typical situation would be a system composed of several equivalent chromophores treated as individual subsystems. In this case, the local response approximation leads to uncoupled excited states of the subsystems (hence the acronym FDEu is employed often), while the subsystem TDDFT code couples the monomer excitations to obtain the excited states of the total system (often denoted as coupled frozen density embedding, FDEc). This can be related to excitonic couplings between the monomers [34].

The current implementation is restricted to NOSYM calculations and Singlet-Singlet excitations without frozen core approximation. It makes use of the ALDA kernel (including a Thomas-Fermi part for the contribution arising from the non-additive kinetic energy) for consistency with the uncoupled FDE implementation for excited states. Some features have not or not extensively been tested and should be used with great care, e.g., linear dependencies in the basis set. Details on the calculation of transition moments, oscillator and rotational strengths are described in Ref. [10].

Subsystem TDDFT (FDEc) calculations can be invoked with the SUBEXCI key.

SUBEXCI input

FDEc calculations on coupled excited states first require that an uncoupled FDE-TDDFT calculation has been performed for every subsystem that should be included in the coupled calculation, and that the corresponding adf.rkf files, in which the considered subsystems are “active”, have been saved (see the separate FDE input description). This means that it is not possible to use the information on frozen/inactive fragments from a adf.rkf file of a previous uncoupled FDE calculation, which contains all subsystems.

Although it is technically possible to use adf.rkf files from non-FDE calculations on the separate subsystems, this would lead to results that are inconsistent with the subsystem TDDFT methodology from Ref. [35]. In any case, a previous TDDFT calculation for each subsystem that should be included in the coupling procedure is necessary. If that is not the case, the subsystem will still be considered in the calculation of the total electron density (needed in the setup of the exchange-correlation kernel), but will not be included in the coupling procedure.

The first subsystem should always be one of the coupled subsystems. The input will then look like the corresponding input for an uncoupled FDE-TDDFT calculation, but in addition should contain the following block:

SUBEXCI
  {LOWEST nlowest}
  {OPTSTATES list_of_optstates}
  {CTHRES coupling_threshold}
  {SFTHRES solutionfactor_threshold}
  {COUPLBLOCK}
  {TDA}
  {CICOUPL}
END

LOWEST nlowest

The selection of the excited states to be coupled consists of two steps. First, a number of reference states are selected. As a default, the nlowest (default: 10) lowest excited states present on the fragment file for the first subsystem are considered.

OPTSTATES list_of_optstates

If the keyword OPTSTATES is given, only those excited states of the first subsystem are considered as reference states that are given in the list_of_optstates (numbers of states separated by blanks).

CTHRES coupling_threshold

Second, all excitations of all subsystems (present on the fragment adf.rkf files) with an excitation energy that differs by less than coupling_threshold (to be given in units of eV; default: 30 eV) from one of the reference states are selected to be included in the coupling. Note that additional excited states of system 1 may be included here.

COUPLBLOCK

If COUPLBLOCK is specified in the input, all couplings between all of these local excited states are included. Otherwise (default), the coupling_threshold is also applied to select pairs of states for which couplings are calculated. I.e., couplings are not calculated if the two particular states to be coupled differ in energy by more than coupling_threshold.

SFTHRES solutionfactor_threshold

To reduce the computational effort, it is possible to ignore the effect of orbital pairs with coefficients less than solutionfactor_threshold in the solution factors (TDDFT eigenvectors) of the underlying uncoupled calculation in the construction of the exact trial densities during the calculation of the coupling matrix elements. These orbital pair contributions are not ignored in the subsequent calculation of transition moments, oscillator, and rotational strengths. The default value of 0.00001 typically leads to a precision of the coupled excitation energies of about 0.0001 eV.

TDA

TDA specifies the use of the Tamm-Dancoff-Approximation (Tamm-Dancoff approximation) in the underlying uncoupled FDE-TDDFT calculations (Ref. [36]). Contrary to the full SUBEXCI-TDDFT variant, SUBEXCI-TDA allows for the usage of hybrid functionals in the underlying uncoupled FDE-TDDFT calculations.

CICOUPL

Within the Tamm-Dancoff Approximation, the couplings between localized excited states on different subsystems correspond directly to so-called exciton couplings (see Ref. [36]). The CICOUPL keyword, in conjunction with TDA, prints these exciton couplings. It is also possible to use CICOUPL with full FDEc-TDDFT. In that case, the excitonic couplings between monomers are reconstructed from an effective 2x2 CIS-like eigenvalue problem, as e.g. done in Ref. [34].

In addition, the input file may contain either an EXCITATIONS block or the keyword DIFFUSE. Both options lead to a slight adaption of the integration grid. Apart from this, the EXCITATIONS block will be ignored.

The key ALLOW PARTIALSUPERFRAGS is currently necessary to be able to use subsystem information for only one subsystem from a adf.rkf file of a previous FDE calculation:

ALLOW PARTIALSUPERFRAGS

Restrictions and pitfalls¶

In the current implementation, only the electron density of the embedded system is calculated. Therefore, with the exception of interaction energies, only properties that depend directly on the electron density (e.g., dipole moments) are available. In addition, the TDDFT extension allows the calculation of electronic excitation energies and polarizabilities, and NMR shieldings can be calculated.

Kinetic energy approximant: Although the effective embedding potential is derived from first principles using universal density approximants, the ADF implementation relies on approximations. Currently, two implemented approximations are recommended [15]: PW91k (also known as GGA97) which uses electron densities and the corresponding gradients to express the non-additive kinetic energy component of the embedding potential, or TF (Thomas-Fermi LDA approximant), which does not use gradients at all. Either approximation is applicable only in cases where the overlap between electron densities of the corresponding interactions is small. Note: so far, no approximation has been developed for the strong-overlap case - two subsystem linked by covalent bonds for instance.

References