# FDE: Frozen Density Embedding¶

The Frozen-Density-Embedding (FDE) option invokes calculation of the effective embedding potential introduced by Wesolowski and Warshel [1] in order to take into account the effect of the environment on the electronic structure of an embedded system. The embedding potential (Eq. 3 in Ref. [5]) depends explicitly on electron densities corresponding to the embedded subsystem (e.g. a solvated molecule) and its environment (e.g. solvent). For a detailed review, see Ref. [3]. The ADF implementation of the method is described in detail in Ref. [2] [4].

A time-dependent linear-response generalization of this embedding scheme was derived in Ref. [6]. Its implementation in an approximate form, which assumes a localized response of the embedded system only (uncoupled FDE), is described in the supplementary material to Ref. [7]. For possible drawbacks and pitfalls in connection with this approximation, see Refs. [2] [8] [9].

The theory of coupled excited states for subsystems is described in Refs. [35] [34], and extended for general response properties in Ref. [10]. This theory (subsystem TDDFT, coupled FDE) allows to treat the mutual response of several subsystems, including the ones that are considered environment.

A generalization of the FDE scheme to the calculation of NMR shieldings has been given in Ref. [11], where also the approximations involved and possible problems are discussed.

With the exception of interaction energies, the current implementation in ADF only allows the calculation of molecular properties that only depend on the electron density and of response properties using TDDFT. For an application to the calculation of several molecular properties in solution and a comparison to the DRF model also available in ADF, see Ref. [8]. For further applications of the ADF implementation, see Ref. [12] (weakly interacting complexes), Refs. [2] [8] [37] [38] (solvent effects), and Refs. [39] [40] (other environment effects).

To invoke a frozen-density embedding calculation, two additional specifications in the input are required. First, one or more frozen fragments have to be included in the FRAGMENTS block, and second, the block key FDE has to be included. In the simplest case, this input should look like this:

```
FRAGMENTS
...
FragType FragFile type=FDE
...
END
FDE
PW91K
end
```

In the FRAGMENTS block, for any fragment it is possible to specify the option type=FDE to indicate that the density of this fragment is kept frozen. This density is imported from the file FragFile. The frozen fragments have to be included in addition to the usual, non-frozen fragments. The atoms of the frozen fragments have to be included in the ATOMS block. As with normal fragments, the fragment found in the file will be rotated and translated to its position specified in the ATOMS block. For more details on specifying fragments, see the section ‘fragment files’. In the FDE input block, the recommended PW91k (also known as GGA97) approximant is recommended for the non-additive kinetic energy (the default is the local density approximant). A recommended alternative is NDSD. For all other options the defaults will be used.

Please note that throughout the FDE part of the documentation, the word “approximant” is used instead of the more usual “functional” to emphasize that the exact functional is not known, also in the case of the kinetic energy functional. In the literature one may encounter both words used interchangeably.

By including more than one frozen fragment, it is possible to use a frozen fragment that is a superposition of the densities of isolated molecules (this was possible in the previous version of ADF using the DENSPREP option). For a discussion and tests of the use of such approximate environment densities, see Ref. [2].

There is no restriction on the use of symmetry in FDE calculations, and usually the correct symmetry will be detected automatically. However, in the preparation of frozen fragments that will be rotated and/or translated in the FDE calculation, one has to include the keyword NOSYMFIT for technical reasons.

In the current implementation, only the electron density of the embedded (non-frozen) 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 particular, the calculation of energy gradients is not implemented yet. All quantities given in the output refer (unless explicitly specified otherwise) to the non-frozen system only.

To employ the extension of FDE to the calculation of NMR shieldings, the file TAPE10 has to be used in the FDE calculation (by including the option SAVE TAPE10), and subsequently the NMR shielding has to be calculated using the program NMR (not with EPR). The TDDFT extension of the FDE formalism allows the calculation of electronic excitation energies and polarizabilities. This extension is automatically activated if FDE is used in combination with the EXCITATIONS or the RESPONSE key. To allow the mutual response of several subsystems, see the section on [subsystem TDDFT].

## 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 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].
`FDEOPTIONS OPTIMIZE`

If the key

`FDEOPTIONS OPTIMIZE`

is not specified, then only the active fragment is optimized in a Geometry Optimization run (the FDE fragments are not optimized). If you include`FDEOPTIONS OPTIMIZE`

, then the corresponding FDE fragment will also be optimized. The`OPTIMIZE`

key implicitly activates the option`RELAX`

, which is necessary to evaluate the nuclear gradients of the respective fragment [Ref. [14]. The`OPTIMIZE`

key only supports mono-molecular basis set expansions and can therefore NOT be combined with the`USEBASIS`

option. For Geometry Optimizations the usage of the`FULLGRID`

option is recommended. See the example Example: Geometry optimization ICW FDE/sSDFT.**Note**:: When performing geometry optimization ICW FDE one must:- Disable symmetry (
`Symmetry NOSYM`

) - Use the old optimizer branch (
`Branch Old`

) - Use Cartesian coordinates in the optimization procedure (
`Optim Cartesian`

)

- Disable symmetry (

`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 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 TAPE21 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 EXCITATION 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.

## 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 TAPE21 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 TAPE21 file of a previous uncoupled FDE calculation, which contains all subsystems.

Although it is technically possible to use TAPE21 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 TAPE21 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 EXCITATION block or the keyword DIFFUSE. Both options lead to a slight adaption of the integration grid. Apart from this, the EXCITATION block will be ignored.

The key ALLOW PARTIALSUPERFRAGS is currently necessary to be able to use subsystem information for only one subsystem from a TAPE21 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.

Warning

EVERYTHING ELSE IS NOT YET IMPLEMENTED. THE RESULTS OBTAINED FOR OTHER PROPERTIES MIGHT BE MEANINGLESS.

**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

[1] | (1, 2, 3) T.A. Wesolowski and A. Warshel, Frozen Density Functional Approach for ab-initio Calculations of Solvated Molecules, Journal of Physical Chemistry 97, 8050 (1993) |

[2] | (1, 2, 3, 4, 5) J. Neugebauer, C.R. Jacob, T.A. Wesolowski and E.J. Baerends, An Explicit Quantum Chemical Method for Modeling Large Solvation Shells Applied to Aminocoumarin C151, Journal of Physical Chemistry A 109, 7805 (2005) |

[3] | (1, 2) T.A. Wesolowski, in: Computational Chemistry: Reviews of Current Trends - Vol. 10, World Scientific, 2006. |

[4] | C.R. Jacob, J. Neugebauer and L. Visscher, A flexible implementation of frozen-density embedding for use in multilevel simulations, Journal of Computational Chemistry 29, 1011 (2008) |

[5] | (1, 2, 3, 4) T.A. Wesolowski and J. Weber, Kohn-Sham equations with constrained electron density: an iterative evaluation of the ground-state electron density of interacting molecules, Chemical Physics Letters 248, 71 (1996) |

[6] | (1, 2, 3) M.E. Casida and T.A. Wesolowski, Generalization of the Kohn-Sham equations with constrained electron density formalism and its time-dependent response theory formulation, International Journal of Quantum Chemistry 96, 577 (2004) |

[7] | (1, 2) T.A. Wesolowski, Hydrogen-Bonding-Induced Shifts of the Excitation Energies in Nucleic Acid Bases: An Interplay between Electrostatic and Electron Density Overlap Effects, Journal of the American Chemical Society 126, 11444 (2004) |

[8] | (1, 2, 3) C.R. Jacob, J. Neugebauer, L. Jensen and L. Visscher, Comparison of frozen-density embedding and discrete reaction field solvent models for molecular properties, Physical Chemistry Chemical Physics 8, 2349 (2006) |

[9] | J. Neugebauer and E.J. Baerends, Exploring the Ability of Frozen-Density Embedding to Model Induced Circular Dichroism, Journal of Physical Chemistry A 110, 8786 (2006) |

[10] | (1, 2) J. Neugebauer, On the calculation of general response properties in subsystem density functional theory, Journal of Chemical Physics 131, 084104 (2009). |

[11] | C.R. Jacob and L. Visscher, Calculation of nuclear magnetic resonance shieldings using frozen-density embedding, Journal of Chemical Physics 125, 194104 (2006) |

[12] | (1, 2, 3) C.R. Jacob, T.A. Wesolowski and L. Visscher, Orbital-free embedding applied to the calculation of induced dipole moments in CO2^{…} X (X=He, Ne, Ar, Kr, Xe, Hg) van der Waals complexes, Journal of Chemical Physics 123, 174104 (2005) |

[13] | (1, 2, 3, 4, 5, 6, 7, 8, 9, 10) T.A. Wesolowski and J. Weber, Kohn-Sham equations with constrained electron density: The effect of various kinetic energy functional parametrizations on the ground-state molecular properties, International Journal of Quantum Chemistry 61, 303 (1997) |

[14] | D.Schlüns, M. Franchini, A. W. Götz, J. Neugebauer, C. R. Jacob ,L. Visscher, Analytical gradients for subsystem density functional theory within the slater-function-based amsterdam density functional program, J. Comput. Chem. 38, 4 (2017) |

[15] | (1, 2, 3, 4) T.A. Wesolowski, Density functional theory with approximate kinetic energy functionals applied to hydrogen bonds, Journal of Chemical Physics 106, 8516 (1997) |

[16] | (1, 2, 3, 4, 5, 6, 7) T.A. Wesolowski, H. Chermette and J. Weber, Accuracy of Approximate Kinetic Energy Functionals in the Model of Kohn-Sham Equations with Constrained Electron Density: the FH^{…} NCH complex as a Test Case, Journal of Chemical Physics 105, 9182 (1996) |

[17] | (1, 2, 3, 4) Y.A. Bernard, M. Dulak, J.W. Kaminski and T.A. Wesolowski, The energy-differences based exact criterion for testing approximations to the functional for the kinetic energy of non-interacting electrons, Journal of Physics A 41, 55302 (2008) |

[18] | L.H. Thomas, The calculation of atomic fields, Mathematical Proceedings of the Cambridge Philosophical Society 23, 542 (1927) |

[19] | E. Fermi, Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente, Zeitschrift für Physik 48, 73 (1928) |

[20] | P. Cortona, Self-consistently determined properties of solids without band-structure calculations, Physical Review B 44, 8454 (1991) |

[21] | A. Lembarki and H. Chermette, Obtaining a gradient-corrected kinetic-energy functional from the Perdew-Wang exchange functional, Physical Review A 50, 5328 (1994) |

[22] | (1, 2, 3) J.M. Garcia Lastra, J.W. Kaminski and T.A. Wesolowski, Orbital-free effective embedding potential at nuclear cusps, Journal of Chemical Physics 129, 074107 (2008) |

[23] | H. Lee, C. Lee and R.G. Parr, Conjoint gradient correction to the Hartree-Fock kinetic- and exchange-energy density functionals, Physical Review A 44, 768 (1991) |

[24] | (1, 2) D.A. Kirzhnits, Soviet Physics JETP-USSR 5, 64 (1957) |

[25] | A.J. Thakkar, Comparison of kinetic-energy density functionals, Physical Review A 46, 6920 (1992) |

[26] | (1, 2) H. Ou-Yang and M. Levy, Approximate noninteracting kinetic energy functionals from a nonuniform scaling requirement, International Journal of Quantum Chemistry 40, 379 (1991) |

[27] | (1, 2, 3) A.W. Götz, S.M. Beyhan and L. Visscher, Performance of Kinetic Energy Functionals for Interaction Energies in a Subsystem Formulation of Density Functional Theory, Journal of Chemical Theory and Computation 5, 3161 (2009) |

[28] | M. Ernzerhof, The role of the kinetic energy density in approximations to the exchange energy, Journal of Molecular Structure: THEOCHEM 501-502, 59 (2000) |

[29] | J.P. Perdew, Generalized gradient approximation for the fermion kinetic energy as a functional of the density, Physics Letters A 165, 79 (1992) |

[30] | (1, 2) C.R. Jacob, S.M., Beyhan and L. Visscher, Exact functional derivative of the nonadditive kinetic-energy bifunctional in the long-distance limit, Journal of Chemical Physics 126, 234116 (2007) |

[31] | D.V. Chulhai and L. Jensen, Frozen Density Embedding with External Orthogonality in Delocalized Covalent Systems, Journal of Chemical Theory and Computation 11, 3080 (2015). |

[32] | A. Goez, and J. Neugebauer, A Local Variant of the Conductor-Like Screening Model for Fragment-Based Electronic-Structure Method, Journal of Chemical Theory and Computation 11, 5277 (2015) |

[33] | A. Goez, and J. Neugebauer, Including protein density relaxation effects in first-principles embedding calculations of cofactor excitation energies, Molecular Physics 115, 526 (2016) |

[34] | (1, 2, 3, 4) J. Neugebauer, Photophysical Properties of Natural Light-Harvesting Complexes Studied by Subsystem Density Functional Theory, Journal of Physical Chemistry B 112, 2207 (2008) |

[35] | (1, 2, 3) J. Neugebauer, Couplings between electronic transitions in a subsystem formulation of time-dependent density functional theory, Journal of Chemical Physics 126, 134116 (2007). |

[36] | (1, 2) C. König, N. Schlüter. J. Neugebauer, Direct Determination of Exciton Couplings from Subsystem TDDFT within the Tamm-Dancoff Approximation, Journal of Chemical Physics 138, 034104 (2013) |

[37] | J. Neugebauer, M.J. Louwerse, E.J. Baerends and T.A. Wesolowski, The merits of the frozen-density embedding scheme to model solvatochromic shifts, Journal of Chemical Physics 122, 94115 (2005) |

[38] | J. Neugebauer, M.J. Louwerse, P. Belanzoni, T.A. Wesolowski and E.J. Baerends, Modeling solvent effects on electron-spin-resonance hyperfine couplings by frozen-density embedding, Journal of Chemical Physics 123, 114101 (2005) |

[39] | M. Zbiri, M. Atanasov, C. Daul, J.-M. Garcia Lastra and T.A. Wesolowski, Application of the density functional theory derived orbital-free embedding potential to calculate the splitting energies of lanthanide cations in chloroelpasolite crystals, Chemical Physics Letters 397, 441 (2004) |

[40] | M. Zbiri, C.A. Daul and T.A. Wesolowski, Effect of the f-Orbital Delocalization on the Ligand-Field Splitting Energies in Lanthanide-Containing Elpasolites, Journal of Chemical Theory and Computation 2, 1106 (2006) |