Optimizing Kinetic Energy Functionals
for Application in Subsystem Density
Functional Theory
Stephanie Alexa Grimmel
Supervisors: Prof. Dr. Lucas Visscher Dr. Tiago Quevedo Teodoro
Vrije Universiteit Amsterdam
In this graduate thesis, the performance of different nonadditive kinetic energy func-tionals for the calculation of excitation energies via time-dependent frozen density em-bedding, commonly referred to as “uncoupled FDE”, is studied. Regular Kohn-Sham TD-DFT results were taken as reference. Re-parametrization of a PBE-like nonadditive kinetic energy functional was also assessed as an attempt of systematically improving results. Calculations were carried out on a set of molecular dimers. For each system the first singlet-triplet and one singlet-singlet excitation energy were assessed.
Results showed a limited dependence of the obtained accuracy with regard to the chosen nonadditive kinetic energy functional. With the exception of a few outliers (es-pecially Karasiev’s PBEn functionals), the resulting quality differences are negligible. Most notably, even a simple LDA functional yields deviations similar to those obtained with more complex GEA- and GGA-type approximants. Re-optimization of the param-eters of a PBE-like functional did not yield meaningful improvements either.
Overall, no particular nonadditive kinetic energy was found to be preferred over the others for computing excitation energies via uncoupled FDE.
I am very grateful to my supervisor, Professor Lucas Visscher, for his guidance and many fruitful discussions during the realization of this research project. Thanks are also due to Dr. Tiago Teodoro for his co-supervision, assistance and advice.
Furthermore, I highly appreciate the effort undertaken by the SCM support team in order to quickly debug some the features of ADF that were necessary to conduct the research presented here. The advice given by Micha l Handzlik upon the use of the PLAMS package has been a great help.
The whole Theoretical Chemistry group at the Vrije Universiteit Amsterdam is thanked for the friendly working atmosphere.
This research project was carried out in the framework of the European joint Master’s degree “AtoSiM”. I therefore want to express my gratitude towards everyone involved in the design and execution of this study program. Special thanks go to Alexandra Guilleminot at ENS Lyon for her continuous support throughout the last two years. I am also indebted to my fellow AtoSiM students for making the program such a remark-able experience. The financial support granted by the European Commission is highly acknowledged.
I am very grateful for having been a fellow of the Studienstiftung des deutschen Volkes throughout my bachelor’s and master’s studies.
1 Introduction 1
2 Theoretical Background 3
2.1 Kohn-Sham Density Functional Theory . . . 3
2.2 Subsystem Density Functional Theory . . . 4
2.2.1 Interaction Energy . . . 7
2.2.2 Nonadditive Kinetic Potential . . . 8
2.3 Time-dependent Density Functional Theory . . . 11
2.3.1 Time-dependent Kohn-Sham Density Functional Theory . . . 11
2.3.2 Time-dependent Subsystem Density Functional Theory . . . 13
2.4 Design and Optimization of Density Functionals . . . 14
2.4.1 Functional Form and Training Set . . . 15
2.4.2 Cost Functions . . . 15
2.4.3 The Nelder-Mead Optimization Method . . . 16
2.5 Overview of Previously Suggested ˜TnadS . . . 18
2.6 Density Functional Based Tight Binding . . . 21
3 Methodology 23 4 Results and Discussion 25 4.1 Selection of Reference Systems and Excitations . . . 25
4.2 Performance of different ˜TnadS on the Training Set . . . 27
4.3 Re-parametrization of a PBE-like ˜TnadS . . . 32
4.3.1 Re-parametrization for the Entire Training Set . . . 33
4.3.2 Re-parametrization for a Training Set without 3rd Row Elements . 35 4.4 Performance of Different ˜TnadS on the Test Set . . . 37
4.4.1 Comparison of the Performance of Different ˜TnadS on the Test Set . 37 4.4.2 Transferability of the excPBE-3 Parametrization . . . 40
5 Conclusions and Outlook 41
Bibliography 43
List of Figures and Tables 49
The efficient use of computational quantum chemistry for detailed modeling of extended systems (e.g., biomolecules, condensed phase species, etc.) has been a recurring is-sue [1, 2]. This is due to the fact that conventional quantum mechanical methods suffer from the so-called “curse-of-dimensionality” problem, i.e., they scale unfavorably with the system size [1]: While one would hope for a linear increase of computational effort with the number of basis functions M , the full configuration interaction (full CI) method shows a factorial scale and even the effort necessary to perform conventional Kohn-Sham DFT calculations grows as M3 [3, 4].
Nevertheless, procedures that achieve a more favorable scaling with the system size exist. These techniques usually rely on dividing the original, large system, also referred to as the supersystem, into smaller subsystems and performing calculations on the in-dividual subunits. Then, approximations are used in order to obtain information on the supersystem, or on the effect of the environment, e.g., a solvent, on one selected subsystem [1].
Subsystem density functional theory is one example of these techniques. It partitions the entire system based on the electronic density ρ(r) [1]. In the limit of exact energy functionals, subsystem density functional theory would be exact. However, these exact functionals are not known and approximations have to be used. Besides the exchange-correlation functionals, which are also employed in regular Kohn-Sham DFT and have therefore been well-developed, also the so-called “nonadditive kinetic energy functional” needs to be approximated. The latter arises as a correction term when trying to re-construct the total Kohn-Sham kinetic energy from the contributions of the individual subsystems.
The standard class of functionals used for approximating the kinetic energy is based on the general gradient approximation (GGA) [1]. However, compared to their exchange-correlation counterparts, the amount of work that has been done to systematically op-timize their parameters is rather limited [5]. Furthermore, the vast majority of these attempts focuses on obtaining accurate estimates for the total non-interacting kinetic energy as it would be necessary for orbital free density functional theory, and not only its nonadditive part that is required in subsystem DFT. Firstly, a re-parametrization with the aim of optimizing this contribution only seems promising as it is significantly smaller than the entire contribution of the kinetic energy functional. Secondly, and even more importantly, it has been reported that there is not necessarily a correlation between the performance of a GGA-type nonadditive kinetic energy functional and that of its parent approximative kinetic energy functional applied to model the entire Kohn-Sham kinetic energy [6, 7].1
A common application of time-dependent subsystem DFT is the computation of
ex-1
Obviously, this statement does not apply to a hypothetic, exact total non-interacting kinetic energy functional, for which also the nonadditive counterpart derived according to Section 2.2.2 would be exact.
citation energies [1, 8, 9, 10]. However, again, functionals provided in the literature are usually not constructed for yielding accurate results for this property. A systematic assessment of their capabilities with regard to spectroscopic properties has not yet been carried out either [5].
The research presented here tries to address this problem by comparing excitation energies obtained via uncoupled frozen density embedding (FDEu), a popular variant of time-dependent subsystem DFT, with a range of previously suggested nonadditive kinetic energy functionals to the corresponding results from standard, supermolecular time-dependent response DFT. Additionally, a re-optimization of the parameters of a given GGA-type PBE-like nonadditive kinetic energy functional with the aim of a better reproduction of supermolecular results by FDEu calculations was tried.
In the following sections, firstly, the background necessary to understand the issue described above is reviewed. Afterwards, the obtained results are presented and inter-preted. Finally, an outlook on possible future research is given.
Atomic units are used throughout the theoretical introduction of this report. If appli-cable, all formulas are given for the spin-compensated case to maintain good readability. A generalization of the kinetic energy functional for spin-polarized systems is presented in Reference 11.
2.1 Kohn-Sham Density Functional Theory
According to the Hohenberg-Kohn theorems the electronic energy can be expressed as a functional of the electron density, ρ(r) [12]. However, the exact form of this functional is unknown and its evaluation is expected to be as complicated as solving the Schr¨odinger equation directly [3]. Thus, approximations have to be applied. The energy functional resulting from the popular Kohn-Sham approach is the following [1, 3, 13]:
E[{ψi}] = Vnn+ Vnuc[ρ] + J [ρ] + TS[{ψi}] + Exc[ρ] (2.1)
Vnn arises from the repulsion between nuclei thus being constant for fixed nuclear
posi-tions: Vnn= X I,J ZIZJ |RI − RJ| (2.2)
I and J are indices of nuclei and Z and R their charge number and position. The interaction energy of nuclei and electrons, Vnuc[ρ], can be expressed easily as a functional
of the electron density:
Vnuc[ρ] = Z ρ(r)vnuc(r)dr (2.3) vnuc(r) = − X I ZI |RI − r| (2.4)
The exact form of the density functional accounting for the electron-electron interaction energy is unknown. Thus it is represented by the sum of the so-called Hartree term, J [ρ], that yields the interaction energy of the electron density with itself,
J [ρ] = 1 2 Z ρ(r)ρ(r0) |r − r0| drdr 0 (2.5)
and a correction that accounts for the quantum mechanical effects of correlation and exchange. The correction term is part of the exchange-correlation energy functional, that will be defined explicitly in Equation 2.9.
The kinetic energy of the electrons is not known explicitly either. The key point of the Kohn-Sham technique is to approximate it as the kinetic energy of a system of n non-interacting electrons, TS[ρ]: TS[ρ] := min Φ→ρ * Φ −1 2 n X i ∇2i Φ + (2.6)
Φ stands for a Slater determinant composed of the so-called Kohn-Sham orbitals ψi that
are orthonormal and related to the electron density as follows:
ρ(r) =
n
X
i
|ψi(r)|2 (2.7)
Consequently, the Kohn-Sham kinetic energy can be written explicitly as a functional of these orbitals: TS[{ψi}] = − 1 2 n X i ψi ∇2 ψi (2.8)
The correction terms required to obtain the correct kinetic and electron-electron inter-action energy, T [ρ] and Vee[ρ] respectively, from the approximations introduced above
are combined in the exchange-correlation energy functional:
Exc[ρ] = T [ρ] + Vee[ρ] − (TS[ρ] + J [ρ]) (2.9)
By minimizing the energy functional given in Equation 2.1 with respect to the Kohn-Sham orbitals, under the constraint that these remain orthonormal, the Kohn-Kohn-Sham equations can be derived:
−1 2∇ 2+ v S[ρ](r) ψi(r) = iψi(r) (2.10)
i is referred to as the Kohn-Sham orbital energy of state i and vS[ρ] is the Kohn-Sham
potential:
vS[ρ] = vnuc(r) + vH[ρ](r) + vxc[ρ](r) (2.11)
with vH[ρ](r) being the Hartree and vxc[ρ](r) the exchange-correlation potential:
vH[ρ](r) = Z ρ(r0) |r − r0|dr 0 (2.12) vxc[ρ](r) = δExc[ρ] δρ(r) (2.13)
To obtain the ground state energy and density of a system the Kohn-Sham equations are solved self-consistently for a given nuclear geometry. In the limit of exact functionals the result would be exact [14, 15].
2.2 Subsystem Density Functional Theory
The basic idea underlying all fragmentation methods is the nearsightedness principle, which states that local electronic properties only depend significantly on the effective ex-ternal potential at nearby points [16]. Subsystem density functional theory relies on the fact that the electron density is additive for subsystems and partitions the supersystem’s density into subsystem contributions, ρI(r) [17, 18, 19]:
ρ(r) =X
I
Consequently, the computation of wave functions for the entire system can be avoided. This approach to divide the system in subsystems benefits from the fact that ρ is an observable real-space quantity, so that a very chemically intuitive partitioning can be achieved [1].
Similarly to what has been outlined in the context of Kohn-Sham DFT in Section 2.1, the electron densities of the subsystems are expressed via the nI occupied orbitals within
subsystem I corresponding to a system of non-interacting particles with the target den-sity ρI(r): ρI(r) = nI X iI |ψiI(r)|2 (2.15)
In analogy to Equation 2.8 the non-interacting kinetic energy of a subsystem I can be defined: TS[{ψiI}] = − 1 2 nI X iI ψiI ∇2 ψiI (2.16)
Although the supersystem’s (Kohn-Sham) kinetic energy is not equal to the sum of the subsystem contributions (and, since we do not know the orbitals corresponding to the supersystem, cannot be evaluated directly), one can still use it as an approximation. Then, in analogy to the reasoning applied when introducing the exchange-correlation energy functional in Section 2.1, a correction is introduced to render the expression formally exact. The correction term is called the nonadditive kinetic energy TSnad.
TSnad[{ρI}] := TS[ρ] − X I TS[ρI] (2.17) ⇒ TS[ρ] = X I TS[ρI] + TSnad[{ρI}] (2.18)
Inserting Equation 2.18 in the Kohn-Sham energy expression given in Equation 2.1, yields the subsystem DFT energy functional:
E[{ψiI}] = Vnn+ Vnuc[ρ] + J [ρ] + Exc[ρ] + X I TS[{ψiI}] + T nad S [{ρI}] (2.19)
From this expression, the “Kohn-Sham equations with constrained electron density” (KSCED) can be obtained by performing a minimization of the energy with respect to the orbitals of one subsystem K under the constraint that these remain orthonormal and while enforcing that the densities of the other subsystems remain unchanged [20].
−1 2∇ 2+ v(K) ef f[ρK](r) + v (K) emb[ρK, ρ](r) ψiK = iKψiK (2.20)
The effective potential v(K)ef f[ρK](r) is defined as follows:
The nuclear potential v(K)nuc(r) takes into account exactly those nuclei that are part of
the subsystem under consideration. Thus, as can be seen from comparison with Equa-tion 2.11, the effective potential corresponds to the Kohn-Sham potential of the isolated subsystem.
The influence of the environment is included via the embedding potential v(K)emb[ρK, ρ](r).
It comprises the effect of all nuclei of the supersystem, that are not part of the subsystem under consideration, as well as the classical Coulomb interaction of the other electrons. Furthermore, the correction terms for the nonadditive parts of the kinetic and exchange-correlation energy are added.
vemb(K)[ρK, ρ](r) =
X
I6=K
v(I)nuc(r) + vH[ρ − ρK](r) + vxcnad[ρK, ρ](r) + vnadkin[ρK, ρ](r) (2.22)
The nonadditive kinetic potential vkinnad[ρK, ρ](r) is the functional derivative of the
nonad-ditive kinetic energy functional, as defined in Equation 2.17, with respect to the density of the subsystem under consideration. Analogously to TSnad[{ρI}], one can also define
the nonadditive contribution of the exchange-correlation energy Excnad[{ρI}] and the
non-additive exchange-correlation potential vnadxc [ρK, ρ](r).
vkinnad[ρK, ρ](r) = δTSnad[{ρI}] δρK(r) (2.23) = δTS[ρ] δρ(r) − δTS[ρK] δρK(r) (2.24) vxcnad[ρK, ρ](r) = δEnad xc [{ρI}] δρK(r) = δExc[ρ] δρ(r) − δExc[ρK] δρK(r) (2.25)
The nonadditive exchange-correlation potential can be calculated directly from the avail-able Exc[ρ]-approximations developed for conventional Kohn-Sham DFT, if their explicit
density dependence is known, i.e., for LDA and GGA-type functionals [21]. If one wants to move higher up on Jacob’s ladder of density functionals to employ, e.g., meta-GGA or hybrid functionals, the subsystem DFT scheme can be adapted [21, 22].
In contrast to that, for the nonadditive kinetic potential no counterpart is available in regular Kohn-Sham DFT. Constructing appropriate approximations can be regarded as the key challenge in applying subsystem DFT [21]. This issue will be discussed more thoroughly in Section 2.2.2.
In the limit of exact functionals, the embedding potential, as defined in Equation 2.22, is exact. Differences between subsystem DFT and Kohn-Sham DFT results for energies and densities obtained for the same system, basis set, exchange-correlation functional etc. can only arise due to approximations made in the nonadditive kinetic energy func-tional [1]. Consequently, under the assumption of an exact nonadditive kinetic energy functional, subsystem DFT is equivalent to standard Kohn-Sham DFT. It is also worth to note that the electrostatic effects of the nuclei and electrons of the environment of the
active subsystem, which are expected to give the main contribution to the embedding potential, are treated exactly in subsystem DFT contrary to, e.g., QM/MM schemes [1]. A popular variant of subsystem DFT, which is also used in the research here, is the frozen density embedding theory (FDE), in which the system is divided into two subsystems, one of which is assumed to be “frozen”, i.e., its electronic density remains unchanged. Only the Kohn-Sham orbitals corresponding to the second fragment, the so-called active subsystem, are optimized so that the energy of the entire system is minimized under the given constraints [20].
2.2.1 Interaction Energy
The interaction energy is the difference between the total energy of the entire system and the sum of the subsystem contributions [1]:
Eint[{ρI}] := E[ρ] −
X
I
EIsub[ρI] (2.26)
The subsystem energies are evaluated from EIsub[{ψiI}] = TS[{ψiI}] + V
(I)
nuc[ρI] + J [ρI] + Exc[ρI] + VnnI (2.27)
i.e., by an expression similar to that of the Kohn-Sham energy functional as given in Equation 2.1 applied to the subsystem density and orbitals. Inserting Equation 2.27 in 2.26 yields the following formula for obtaining the subsystem interaction energy [1]:
Eint[{ρI}] = TSnad[{ρI}] + Excnad[{ρI}] +
X I6=J VnucI [ρJ] + X I<J J [ρI, ρJ] + X I<J VnnIJ (2.28)
The Coulomb interaction energy of ρI and ρJ is
J [ρI, ρJ] =
Z ρ
I(r)ρJ(r0)
|r − r0| drdr
0 (2.29)
VnnIJ is the repulsion between nuclei in the subsystems I and J
VnnIJ =X A∈I X B∈J ZAZB |RA− RB| (2.30)
and Vnuc(I)[ρJ] represents the Coulomb attraction between nuclei in subsystem I and the
electron density of system J :
Vnuc(I)[ρJ] = − X A∈I Z ρJ(r) ZA |RA− r| (2.31)
The fact that the interaction energy is a scalar renders it an attractive measure to com-pare the quality of subsystem DFT results to Kohn-Sham DFT references. Nevertheless, for a more detailed investigation, the accuracy of other quantities, such as the electron density, should also be considered. However, if these functions live in three-dimensional space, a systematic accuracy rating necessarily depends on the construction of a repre-sentative single-valued variable.
2.2.2 Nonadditive Kinetic Potential
Before the most popular techniques to approximate the nonadditive kinetic energy shall be discussed, it is illustrative to explain in more detail why this term arises.
The general reason was already given above: The motivation to use subsystem DFT calculations is to reduce the computational effort by not computing the Kohn-Sham orbitals for the whole system, but only for individual subsystems, so that the Kohn-Sham kinetic energy cannot be evaluated directly from Equation 2.8.
One could imagine a case, though, where the densities of the subsystems arise from subsets of the Kohn-Sham orbitals representing the ground state of the supersystem [20]. If, for simplicity, two subsystems A and B are assumed, and the orbitals giving rise to ρA are chosen to be labelled as orbitals ψ1 to ψnA and the remaining nB= n − nA ones
yield ρB, the overall kinetic energy could indeed be written as the sum of subsystem
contributions: TS[ρ] = TS[{ψi}] = − 1 2 n X i=1 ψi ∇2 ψi = −1 2 nA X i=1 ψi ∇2 ψi + − 1 2 n X i=nA+1 ψi ∇2 ψi (2.32)
However, when taking a close look at the definition of TS in Equation 2.6, it becomes
clear that the two parts of the sum in Equation 2.32 are not necessarily equal to TS[ρA]
and TS[ρB]. The Kohn-Sham orbitals of the supersystem were only assumed to give rise
to the correct subsystem densities, but do not necessarily result in the minimum possible value of the given sums, as required for TS. Thus, even for this case of vAB-representable
densities, it only follows that the summands are equal or greater than TS[ρA] and TS[ρB],
respectively. Substracting TS[ρA] + TS[ρB] from Equation 2.32 yields [23]:
TSnad[{ρI}] = Ts[ρ] −
X
I=A,B
TS[ρI] ≥ 0 (2.33)
Equality to zero is reached for nonoverlapping densities, i.e., if ρA(r)ρB(r) = 0 ∀ r as
demonstrated in Reference 24.
2.2.2.1 Decomposable Approximations for the Nonadditive Kinetic Potential In most practical applications the nonadditive energy is approximated in a decomposable manner, i.e., by constructing approximations for the total non-interacting kinetic energy TS as a functional of ρ, ˜TS[ρ], and then applying Equation 2.24 to obtain an approximate
nonadditive kinetic potential ˜vkinnad [24]. In principle, the construction of a high-quality kinetic energy functional that depends directly on the electron density ρ(r) would allow us to no longer reintroduce orbitals as done in the Kohn-Sham approach, but to di-rectly use the density functional to evaluate the kinetic energy contribution instead. In practice, so-called orbital-free DFT calculations are actively investigated and significant progress has been achieved in recent years, but their applicability is still limited by the
low quality of available approximative kinetic energy functionals ˜TS[ρ] [25, 26].
In subsystem DFT only the nonadditive kinetic potential is obtained from ˜TS[ρ].
Con-sequently, since usually the nonadditive part of the kinetic energy is comparably small, the requirements are lower than for orbital free methods [1].
Exact, analytical forms of the kinetic energy density functional are only available for a few simple systems: For the homogeneous electron gas Thomas and Fermi derived the following expression [27, 28]: TST F = CT Fρ 5 3 CT F = 3 10 3π 223 (2.34)
The local density approximation (LDA) follows immediately by generalizing Equation 2.34 for systems with an inhomogeneous electron density:
˜
TSLDA= CT F
Z
ρ53(r)dr (2.35)
The corresponding nonadditive kinetic potential can be calculated easily by using Equa-tion 2.24. For systems with only one Kohn-Sham orbital the von Weizs¨acker functional is exact [29, 1]:
TSvW = 1 8
Z |∇ρ(r)|2
ρ(r) dr (2.36)
The so-called “regular gradient expansion approximation” (GEA) route uses the Thomas-Fermi functional as the zeroth order term and yields a linear combination of TST F and TSvW when expanded to the second order [15, 30]. However, in prior investigations the GEA approach did not prove to be very successful for higher orders and there is no systematic improvement arising when including higher order terms [1, 31, 32].
In contrast to that, generalized gradient approximation (GGA) functionals are widely used to approximate the kinetic energy as ˜TSGGA [1]:
˜
TSGGA = CT F
Z
ρ53(r)F (s(r))dr (2.37)
F (s) is called the enhancement factor and depends on the reduced density gradient s(r):
s(r) := |∇ρ(r)| 2(3π2ρ4(r))13
(2.38)
The enhancement factor F (s) can be of many different forms, from unity, leading to the reproduction of the LDA functional, to fairly complex expressions. An overview can be found in Reference 1. The functionals compared in this study are furthermore briefly described in Section 2.5.
A popular approach to construct F (s) is to make use of the so-called “conjointness hy-pothesis”, which states that for a given GGA exchange functional ˜Ex[ρ] one can construct
a kinetic energy functional ˜TS[ρ] with an enhancement factor of the same analytical form
and both will perform comparably well [33]. It was even observed that the optimal co-efficients employed within F (s) can be chosen similarly [33]. Although this statement is
attractive, since the research on exchange energy functionals is very advanced compared to that on kinetic energy functionals and a simple transfer of knowledge would thus be desirable, it has to be underlined that there is no theoretical proof for the conjointness conjecture, but only empirical evidence [33, 34]. Indeed it has been shown that the conjointness hypothesis is not strictly true, but still helpful [25, 35, 36].
Currently available kinetic energy functionals fail for systems composed of subsystems connected by covalent bonds, as one would like to treat, for example, when investigat-ing biomolecules, such as proteins with the amino acid buildinvestigat-ing blocks connected by peptide bonds chosen as subsystems [37]. Notably, this limitation is also present in the case of very weak covalent bonds, which indicates that not primarily the bond strength but rather the type of the inter-subsystem bonds determines whether the technique is successful [35].
Reasonable accuracies can be reached for van der Waals and hydrogen-bonded com-plexes, including systems with strong hydrogen bonds [37].
Nevertheless, even for cases where the densities of subsystems only overlap weakly, comparisons to exact potentials generated by numerical inversion of the Kohn-Sham equation (see Section 2.2.2.2) indicate that there is still room for improvements [38].
2.2.2.2 Nondecomposable Approximations for the Nonadditive Kinetic Potential Approximation schemes in which ˜TSnad and ˜vkinnad are not derived from the parent func-tional ˜TS, as described in the previous section, but instead are constructed directly,
are referred to as “nondecomposable” [1, 37]. Typically, they rely on first approximat-ing ˜vnadkin and then constructing the corresponding nonadditive kinetic energy functional,
˜ Tnad
S [37].
In order to solve the KSCED equations, given in Equation 2.20, one does not need to know the nonadditive kinetic energy functional ˜TSnad explicitly, but only ˜vnadkin [37]. Consequently, if the total energy of the system is not of interest, but only quantities that arise from the KSCED solutions directly, e.g., electron densities or excitation energies, then in principle, the construction of ˜TSnad is not even necessary [1, 37].
Starting from the Euler-Lagrange equations of DFT, the following relation between KS-potentials vS, defined in Equation 2.11, and the nonadditive kinetic potential can be
shown [37]:
vkinnad[ρK, ρ](r) = vS[ρK](r) − vS[ρ](r) + ∆µ (2.39)
where ∆µ is a constant shift.
For some special cases, e.g., two subsystems situated at large distances from each other, the above approach allows one to analytically derive the exact vkinnad [1, 39]. If the behavior of the correct vnad
kin differs from that of approximations for the case under
investigation, corrections that enforce the desired behavior may be introduced.
A similar approach revealed possible reasons for the inability of conventional LDAs and GGAs to correctly describe covalently bound subsystems, leading to the suggestion that nondecomposable approximations might fix this problem [40].
In the general case, the Kohn-Sham potential, vS, as a functional of the electron
however, be constructed by numerically inverting the Kohn-Sham equations [1]: Starting from an initial guess for vS, the resulting density is evaluated, which is then compared
to the target density. The potential is adapted accordingly by making it more repul-sive in regions where the obtained density is larger than the wanted one and vice versa. This procedure is repeated iteratively until convergence is reached. Examples of con-crete schemes to implement this procedure and further details can be obtained from Reference 1.
The described techniques, in principle, allow one to obtain exact nonadditive kinetic and therefore embedding potentials. However, these procedures require the performance of supermolecular calculations, rendering this approach useless for standard applications of subsystem DFT in terms of computational efficiency [37]. Nevertheless, the procedure can still be helpful to gain further insight on the correct behavior of nonadditive kinetic potentials and to adapt approximations accordingly [37, 38].
2.3 Time-dependent Density Functional Theory
The methodology described so far aimed at the investigation of ground states. However, to calculate electronic spectra, information on excited states are required. Within the DFT framework, excitation energies are usually obtained via time-dependent response DFT [1]. Besides excitation energies, transition dipole moments µf i and the derived
oscillator strengths (sometimes also referred to as “f values”) can be of interest when describing optical transitions. For an electronic transition from the initial state Ψi to
the final state Ψf the transition dipole moment is defined as follows:
µf i:= hΨf|ˆp|Ψii (2.40)
ˆ
p denotes the electric dipole moment operator. Physically speaking, the transition dipole moment is the amplitude of the oscillating electric moment in a molecule induced by the interaction with an electromagnetic field [41]. The direction of µf i is referred to as its
polarization, while its square relates to the probability of the transition to occur. If the modulus of the dipole moment of a given transition is zero, the transition is forbidden. The oscillator strength ff i of a transition from the initial state i to the final state f is
given by
ff i:=
2
3(Ef − Ei) |µf i|
2 (2.41)
where Ei denotes the electronic energy of state i. By definition, the dimensionless
quantity ff i is positive for absorptions, negative for emissions and zero for forbidden
transitions.
2.3.1 Time-dependent Kohn-Sham Density Functional Theory
The central principle underlying time-dependent DFT (TD-DFT) is the Runge-Gross theorem and its generalization provided by van Leeuwen [42, 43]. Being the time-dependent analog to the Hohenberg-Kohn theorems it states that the mapping between
an external, time-dependent potential vext(r, t) and the electronic density ρ(r, t)
evolv-ing from a fixed initial state is invertible, i.e., all observable quantities of a (no longer necessarily static) many-electron system can be obtained from the time-dependent den-sity ρ(r, t) [42, 44]. Based on this statement, excitation energies are computed via a linear response formalism, which arises from the investigation of the density changes, δρ, caused by a time-dependent perturbation, δvef f [44]. δvef f comprises both the
change of the external potential δvext itself and the potential change induced by the
changes in density, δvind [1]:
δvef f(r, t) = δvext(r, t) + δvind(r, t) (2.42)
δvindcan be expressed in terms of density changes up to the first order using the
Hartree-exchange-correlation kernel fHxc. In the following, equations are written for a description
in Fourier space, with ω denoting the angular frequency:
δvind(r, ω) = Z fHxc(r, r0, ω)δρ(r, r0, ω)dr0 (2.43) fHxc(r, r0, ω) = δvef f(r, ω) δρ(r0, ω) ≈ 1 |r − r0|+ fxc(r, r 0) (2.44) fxc(r, r0) = δ2Exc(r) δρ(r)δρ(r0) (2.45)
The common procedure of evaluating the exchange-correlation kernel fxcas a
frequency-independent functional derivative is referred to as the adiabatic approximation [44]. The density response is described by the following linear response equation:
δρ(r, ω) = Z
χS(r, r0, ω)δvef f(r0, ω)dr0 (2.46)
χS(r, r0, ω) is the independent-particle response function:
χS(r, r0, ω) =
δρ(r, ω) δvef f(r0, ω)
(2.47)
Inserting Equations 2.42 and 2.43 in Equation 2.46 yields
δρ(r, ω) = Z χS(r, r0, ω) δvext(r0, ω) + Z fHxc(r0, r00, ω)δρ(r0, r00, ω)dr00 dr0 (2.48)
or, written in a symbolic matrix notation [45]:
δρ = χS(δvext+ fHxcδρ) (2.49)
⇔ δvext= δρ χ−1S − fHxc
(2.50)
The linear density response, δρ, has poles at the frequencies corresponding to the exci-tation energies [46]. According to Equation 2.50, the exciexci-tation energies can therefore be obtained from the resonance frequencies of χ−1S − fHxc.
(a) Acetone spectrum.
(b) Acetone−DMSO spectrum obtained with TD-DFT.
(c) Acetone−DMSO spectrum obtained with FDEu.
Figure 2.1: Section of the spectra of an isolated acetone molecule (a), the acetone−DMSO dimer as obtained with supermolecular TD-DFT (b) and the acetone−DMSO dimer as obtained with the FDEu method (c). ff i
de-notes the oscillator strength.
2.3.2 Time-dependent Subsystem Density Functional Theory
In the corresponding subsystem TD-DFT variant, introduced by Casida et al. and gen-eralized by Neugebauer, the response of the electronic density to a change in the external potential, in analogy to the ground state subsystem DFT ansatz given in Equation 2.14, is expressed as a sum of subsystem contributions [45, 47]. Written in terms of the Fourier components of the density responses of the subsystems δρI this translates to
δρ(r, ω) =X
I
δρI(r, ω) (2.51)
In line with the FDE method, mentioned in the first part of Section 2.2, in the so-called “uncoupled FDE” (FDEu) approach, the supersystem is divided in one active subsystem and a “frozen” environment [47, 48]. It is then assumed that the response of the density of the environment to a change in the external potential, δρemb, is negligible and therefore
set to zero:
Only the response of the active subsystem, δρact, is treated explicitly [48]. In comparison
to the analysis of the isolated active system with standard TD-DFT, as described in Sec-tion 2.3.1, here the embedding potential originating from the environment influences the calculation of orbitals and their energies [1]. Instead of the Hartree-exchange-correlation kernel provided in Equation 2.44, an effective kernel fef fact arises [45]:
fef fact(r, r0, ω) = 1 |r − r0|+ δ2Exc[ρact+ ρemb] δρact(r0, ω)δρact(r, ω) + δ 2Tnad S [ρact, ρemb] δρact(r0, ω)δρact(r, ω) (2.53)
ρact and ρemb stand for the density of the active system and the embedding/frozen
en-vironment, respectively. Naturally, this procedure is only applicable for the calculation of excitations that are fully localized on the active subsystem [1]. In terms of interpre-tation this is often regarded as advantageous because in most applications the interest lies in exactly these excitations and not in the more complicated spectra that might also (correctly) contain solvent-solvent and intermolecular excitations. Example spectra illustrating this feature are presented in Figure 2.1: The spectrum calculated with FDEu looks qualitatively alike the one obtained for the isolated active system. Quantitative differences arise from the environment effects. The TD-DFT spectrum, however, shows several additional peaks.
However, in principle, the responses of the different subsystems are always coupled and for the investigation of systems with interacting chromophores it is crucial to take this coupling into account [1]. Therefore, the induced potential for each subsystem K should include the density changes within all subsystems explicitly, so that it takes the following form [1, 47]: δv(K)ind(r, ω) =X I Z δvK ef f(r, ω) δρI(r0, ω) +δv K emb(r, ω) δρI(r0, ω) ! δρI(r0, ω)dr0 (2.54)
2.4 Design and Optimization of Density Functionals
With regard to the design of density functionals, kinetic energy density functionals just as well as exchange-correlation functionals, two main approaches can be distinguished: non-empirical and (semi-)non-empirical developments [49]. While the non-non-empirical route relies on constructing functionals and determining all their coefficients from known physical constraints, (semi-)empirical functionals arise from fitting a more flexible functional form, for example a power series with adjustable coefficients enhancing different variables related to physical quantities such as the density, ρ, its gradient, ∇ρ, etc., to reference values. Furthermore, a wide range of combinations of these ideas exist.
As the founding father of the development of non-empirical exchange-correlation func-tionals, Perdew, who for example developed the GGA functional PBE and the meta-GGA TPSS, has to be mentioned [50, 51]. Ground-breaking ideas on the semi-empirical design of exchange-correlation functionals were established by Becke, as presented for instance in Reference 52. Also in the domain of kinetic energy functionals both, the non-empirical and the empirical approach, are employed, as can be seen for example
from References 53 and 34, respectively. However, the majority of the more complex popular density functionals used in DFT today includes some empirical parameters [49].
2.4.1 Functional Form and Training Set
The first elementary step in empirically designing a density functional is the decision upon an analytic form of the functional and the so-called training set. The latter is the collection of reference values used during the fitting procedure. They can be obtained from experiment or from superior computational techniques.
A critical feature of (semi-)empirical density functionals is the number of adjustable parameters they contain. On the one hand, a flexible functional form allows to better reproduce the data provided in the training set. On the other hand, with many pa-rameters comes the risk of losing transferability, since papa-rameters that are only relevant for the chosen reference set might be introduced, and of generating a less smooth func-tional [52]. In any case, the transferability of an optimized funcfunc-tional has to be tested by applying it on a second set of accurately known data points referred to as the test set and investigating to which extent the newly designed functional is able to reproduce these. The number of adjustable parameters can be chosen in advance or by investigating in how far the increase of the quality of the fit with the rising complexity of the functional justifies the introduction of further coefficients. Recently an alternative technique has been suggested [54]: The GGA enhancement factor was expressed as a power series. Then, instead of simply truncating it, it was aimed at finding the best possible number and combination of coefficients to be included.
2.4.2 Cost Functions
Before setting up an optimization routine, a decision upon a proper cost function has to be made, i.e., a function that computes a real number depending on the quality of the chosen parameters and that shall be minimized during the optimization procedure. Common choices are, for example, the mean absolute error (MAE) and the root mean square deviation (RMSD), which are defined as follows:
MAE = 1 N X i |ei| (2.55) RMSD = s 1 N X i e2i (2.56) ei = xi− yi (2.57)
with xi standing for one obtained result, yithe corresponding reference and N the size of
the sample. As a simple adaption to these measures, depending on the application, also non-uniform weights can be considered. Both measures are scale-dependent and return a value in the unit of the compared quantities. They are obviously both negatively oriented, i.e., lower MAEs/RMSDs stand for better results.
Still, it is crucial to also point out the differences between the two cost functions introduced in Equations 2.55 and 2.56: Most importantly, the fact that the errors are squared in the RMSD before taking the average leads to a higher relative importance of large deviations. Thus the RMSD is the measure of choice if outliers with extraordinarily high errors for some data points are especially unwanted. From the perspective of interpretation the MAE is usually easier to access. A popular scale-independent error measure is the mean absolute percentage error (MAPE):
MAPE = 1 N X i |pi| (2.58) pi= 100 · xi− yi yi (2.59)
2.4.3 The Nelder-Mead Optimization Method
Concerning the work that will be presented in the following, the main challenge with regard to the choice of an appropriate algorithm to minimize the cost function is that the gradient of the errors with regard to the parameter changes is unknown. Therefore, it either has to be approximated numerically or optimization methods that do not require its evaluation have to be employed.
The Nelder-Mead method is one of the most used algorithms from the latter cate-gory [55, 56]: It aims at optimizing a nonlinear function f : Rn→ R only via the eval-uation of function values f (x) and without trying to approximate its derivative. The basic idea of this approach, which is also referred to as the downhill-simplex method, is to construct a series of simplexes of reducing size and with decreasing function values at its vertices that should finally gather around the wanted minimum. A simplex of di-mension n is a special n-didi-mensional polygon, that is a convex hull of its n + 1 vertices. For example, a simplex in R2 is a triangle. In further detail, the algorithm, that is also illustrated in Figure 2.2, consists of the following steps [57]:
1. Generation of an initial, non-degenerate working simplex with the vertices
X = {x0, x1, ..., xn} (2.60)
2. Identification of the vertex xi yielding the lowest (best), xl, highest, xh, and
second highest value, xs, of f :
xl = argmin xi∈X f (xi) (2.61) xh= argmax xi∈X f (xi) (2.62) xs= argmax xi∈X\{xh} f (xi) (2.63)
3. Calculation of the centroid c of X\{xh}:
c = 1 n
X
i6=h
Identify best, xl, worst, xhand second worst vertex, xs Initialize simplex
Compute centroid c
Compute reflection point xr
f (xr) < f (xl)
f (xl) ≤ f (xr) ≤ f (xs) f (xr) > f (xs)
Replace xhwith xnew
xnew= xr Get expansion point xe Get contraction point xc
xnew= argmin xi∈{xr,xe} f (xi) f (xc) < min xi∈{xr,xh} f (xi) f (xc) > min xi∈{xr,xh} f (xi)
xnew= xc Shrinking transformation Termination?
No
Figure 2.2: Basic steps performed within the Nelder-Mead minimization algorithm.
4. Transformation of the simplex:
a) Computation of the reflection point xr of the worst vertex xh:
xr:= c + α (c − xh) (2.65)
The parameter α has to be greater than zero and is typically set to one. If f (xl) ≤ f (xr) ≤ f (xs), the vertex xh is replaced by xr and the
transforma-tion is finished. Otherwise the algorithm proceeds according to the following criteria:
If f(xr) < f (xl), that means, the reflection point is the best point
ob-served so far, proceed to b). If f(xr) > f (xs), proceed to c).
b) Determine the expansion point xe:
xe:= c + γ (xr− c) (2.66)
γ has to be greater than one. A common choice is to set it equal to two. If f (xe) < f (xr), xe is accepted as a new vertex instead of xh, otherwise xr is
c) Computation of the contraction point xc starting from xh or xr, depending
on which of the two performs better: xstart= argmin
xi∈{xh,xr}
f (xi) (2.67)
xc:= c + β (xstart− c) (2.68)
where β is a parameter between zero and one and typically set to 0.5. If f (xc) < f (xstart) the contraction point is accepted as a new vertex instead
of the worst one. Otherwise Step d) is carried out.
d) Performing a shrinking transformation: All but the best vertex are recom-puted according to the following protocol:
xj := xj+ δ (xj− xl) (2.69)
δ is to be chosen between zero and one and typically set to 0.5.
5. Decision upon termination: In order to choose upon finishing, it is tested whether the working simplex became sufficiently small (domain convergence test) and/or whether the corresponding function values are acceptably close to each other (function-value convergence test). Otherwise the optimization continues by re-peating the outlined actions from Step 2 onwards. In practical implementations there is an additional termination criterion in the form of a maximum number of iterations.
Its intuitive set-up is one of the reasons that render the Nelder-Mead algorithm a pop-ular choice for unconstrained optimizations. Especially for applications where function evaluations are expensive their limited number per iteration – usually one or two, unless a shrinking step is performed – is advantageous. However, there is also serious criticism of the method, including that convergence to non-stationary points is possible [58]. In general, a proper mathematical investigation of the algorithm has proved to be surpris-ingly difficult [56]. As for most optimization routines, there is no guarantee that the minimum found by the algorithm corresponds to a global minimum.
2.5 Overview of Previously Suggested Nonadditive Kinetic
Energy Functionals
The enhancement factors of various nonadditive kinetic energy functionals proposed in the literature are displayed in Table 2.3.
Following the well-known LDA functional, labeled as “THOMASFERMI”, the second order GEA derived by Kirzhnits, TF9W, is listed [30].
The PW91k functional is a GGA with six parameters [59]. Its form corresponds to the exchange functional suggested by Perdew and Wang in 1991 [64]. During its parametrization the difference between c4 and c3 was held constant at 275 2 and also c5
Functional Enhancement Factor F (s) References THOMASFERMI3 1 27, 28 TF9W 1 + 5 27s 2 30 PW91k4 1+c2s sinh −1(c 1s)+ c3−c4e−c5s 2 s2 1+c2s sinh−1(c1s)+c6s4 59 NDSD – (nondecomposable) 7 THAKKAR92 1 + 0.0055 (βs) 2 1 + 0.0253βs sinh−1(βs)− 0.072βs 1 + 25/3βs 60 PW86k 1 + 1.296s2+ 14s4+ 0.2s61/15 61, 62 LLP91 1 + 0.0044188 (βs) 2 1 + 0.0253βs sinh−1(βs) 33 OL91A 1 + 5 27s 2+ 0.0187βs 63 OL91B 1 + 5 27s 2+ 0.0245βs 1 + 25/3βs 63 P92 1 + 88.3960s 2+ 16.3683s4 1 + 88.2108s2 64 E00 135 + 28s 2+ 5s4 135 + 3s2 32 PBE2 1 + 2.0309s 2 1 + 0.2942s2 25 PBE3 1 − 3.7425s 2 1 + 4.1355s2 + 50.258s4 (1 + 4.1355s2)2 25 PBE4 1 − 7.2333s 2 1 + 1.7107s2 + 61.645s4 (1 + 1.7107s2)2 − 93.683s6 (1 + 1.7107s2)3 25 TW02 1 + 0.2319s 2 1 +0.23190.8438s2 34 APBEK 1 + 0.23889s 2 1 +0.238890.804 s2 53 revAPBEK 1 + 0.23889s 2 1 +0.238891.245 s2 53 3
Note that this expression does in general not yield the original Thomas-Fermi kinetic energy functional given in Equation 2.34, but the one provided in Equation 2.35, i.e., the label “THOMASFERMI”, that
remained fixed. The other four parameters were optimized in order to properly reproduce the total Kohn-Sham kinetic energies of a set of atoms.
The NDSD functional was constructed in a nondecomposable manner [7]: The cor-rect behavior of the nonadditive kinetic potential for a supersystem composed of two subsystems in the limit of non-overlapping subsystem densities and the uniform elec-tron gas limit was imposed. Additionally, also the limit of zero elecelec-tron density in the first subsystem with the density of the second one integrating to two was enforced to be described correctly. While there is no analytic form for the enhancement factor, an expression for the resulting nonadditive kinetic energy functional and potential can be found in Reference 7.
The conjointness conjecture was applied in order to construct the PW86k functional from the equally named GGA exchange functional [61, 62]. The form of the LLP91 enhancement factor was obtained from the exchange functional suggested by Becke in 1988 [33, 65]. Its parameters were fitted in order to correctly reproduce the Hartree-Fock kinetic energies of rare gas atoms in an orbital-free approach.
The functionals proposed by Ou-Yang and Levy, OL91A and OL91B, are modifications of the second order GEA, TF9W, with an additional term added such that a physically required scaling law, the so-called non-uniform scaling property, is still fulfilled [63]. For both functionals, the contained empirical parameter was obtained by setting the yielded kinetic energy equal to the Hartree-Fock kinetic energy of the xenon atom.
The form of the enhancement factor proposed by Thakkar is achieved by combining those of the OL91B and LLP91 functionals, while its parameters were obtained by a least-square fit to the kinetic energy of 77 molecules [60].
The kinetic energy enhancement factor derived by Ernzerhof, here denoted as E00, and Perdew’s P92 functional share a similar form. The E00 functional is a GGA that was constructed to reproduce the fourth order GEA while at the same time yielding the van Weizs¨acker expression given in Equation 2.36 in the limit of an infinite reduced density gradient s. The P92 GGA functional was obtained by resummation of the sixth order GEA of the kinetic energy density functional that suffers from divergence near nuclei and at large distances for atomic densities [31, 64].
Of special interest for the research project presented here are functionals with an en-hancement factor of the form of the one present in the exchange part of the famous PBE exchange-correlation functional developed by Perdew, Burke and Ernzerhof [50]. Representatives of this family of kinetic energy functionals, which will be referred to as “PBE-like” in the following, are PBE2, TW02, APBEK and revAPBEK. The enhance-ment factor they all share, with the only difference being the parametrization, is the following:
FP BE(s) = 1 +
c1s2
1 + c2s2
(2.70)
An overview of the parameters c1 and c2 for the different functionals when written
in the form of Equation 2.70 is provided in Table 2.4. The PBE2 parameters were fitted in order to reproduce Born-Oppenheimer forces along the silicon-oxygen bond for various bond distances via orbital-free DFT [25]. Thus, only the gradient of the
Table 2.4: Comparison of parameters employed in different PBE-like kinetic energy func-tionals. Functional c1 c2 PBE2 2.0309 0.2942 TW02 0.2319 0.2748 APBEK 0.23889 0.297 revAPBEK 0.23889 0.192
energy was taken into account in this parametrization scheme. Tran and Wesolowski determined the parameters in their TW02-functional by requiring that it would yield the correct overall non-interacting kinetic energy for the helium and xenon atoms [34]. The APBEK kinetic energy functional was constructed non-empirically starting from an expansion of the non-interacting kinetic energy functional for the non-relativistic semiclassical neutral atom [53, 66]. The nonadditive kinetic energy functional arising from APBEK was applied in subsystem DFT calculations of eight weakly interacting systems, during which it was observed that with the same parameter c1 as derived
earlier, a lower parameter c2 yields more accurate total energies. The optimal c2 that
was found during these investigations is included in the revised APBEK functional, revAPBEK [53]. The revAPBEK and the NDSD enhancement factors are, among those included in the present examination, the only ones designed specifically to yield accurate results in subsystem DFT opposed to entirely orbital-free approaches.
Finally, also the PBE3 and PBE4 functionals were included in the analysis [25]. These functionals contain three or four empirical parameters, respectively, which were deter-mined anagously to those of PBE2, but with a larger training set that includes silicon-oxygen bonds of three different molecules.
2.6 Density Functional Based Tight Binding
Density Functional Based Tight Binding (DFTB) can be understood as a simplification of KS-DFT that allows, even though at the cost of reduced accuracy, to speed up calcu-lations by two to three orders of magnitude in comparison to standard DFT [67, 68]. It is based on writing the Kohn-Sham energy functional, given in Equation 2.1, in terms of a density that is expressed as the sum of a reference density ρ0 and a
perturba-tion δρ. The exchange-correlaperturba-tion energy Exc[ρ0 + δρ] is then expanded as a Taylor
series. Depending on the range of the expansion one arrives at a pure tight binding scheme (TB, 0th order expansion), standard DFTB (1st order expansion), DFTB with self-consistent charge correction (SCC-DFTB, 2ndorder expansion) or DFTB3 (3rdorder expansion) [68, 69, 70, 71, 72]. For all schemes the reference density is usually expressed
as a superposition of the atomic densities ρα 0: ρ0(r) = X α ρα0(rα) (2.71) rα= r − Rα (2.72)
Rα denotes the position of the nucleus of atom α. The first contribution incorporated
in the energy obtained from any of the listed techniques is resulting from evaluating the Hamiltonian in the atomic basis set. The diagonalization of this Hamilton matrix constitutes the dominating part of the calculation in terms of computational cost [68]. Other terms include the so-called repulsive contribution and, for all but the simple TB scheme, corrections due to density fluctuations. Numerous different approximative schemes are applied in this regime in order to properly account for energy changes with regard to the reference system. The parameter sets that are employed in the approximations are in general obtained via DFT calculations. Further details on DFTB methods can be found, e.g., in References 68 and 73.
Throughout this research excitation energies obtained via uncoupled FDE are compared to reference values calculated with supermolecular TD-DFT. The difference between these two values is subsequently often referred to as “the error” of the FDEu excitation energies.
The geometries of the molecular dimers that were equally used for supermolecular and FDE calculations were optimized with the DFTB engine of the ADF 2017 Modeling Suite [74]. The DFTB3 model was applied along with the Slater-Koster files for organic and biomolecules, 3ob-3-1 [75, 68, 76, 77]. Furthermore, Grimme’s dispersion correction D3-BJ was turned on [78].
All TD-DFT and FDEu calculations were performed with ADF 2017.212, release 64756 [79, 80, 81, 82, 83]. A double zeta plus polarization basis set (DZP) was em-ployed [84]. The application of Slater-type basis functions, which is a general feature of ADF, requires density fitting in order to allow for an efficient evaluation of the Coulomb integral. Throughout this work the STOFIT density fitting approach was employed [85]. A “very good” Becke-type integration grid was utilized [86, 87]. The integration grid during FDE calculations was constructed according to the ADF default option outlined in Reference 88, i.e., by taking into account only atoms of the active system and atoms of the embedding environment that are close to the active subsystem. No frozen core was applied (i.e., all electrons were treated explicitly).
To approximate the exchange-correlation contributions, the range-separated hybrid functional CAMY-B3LYP from the XCFun library was used [89, 90, 91, 92]. The choice of the functional is motivated by the fact that range-separated hybrids allow for a better treatment of charge-transfer excitations [93]. Including exact exchange reduces the risk of observing spurious charge-transfer artifacts in the supermolecular reference calcula-tions [94]. The new Hartree-Fock RI scheme was not employed.
However, due to their orbital dependence, (range-separated) hybrid corre-lation functionals cannot be applied to model the non-additive part of the exchange-correlation potential in ADF. Consequently, this contribution was accounted for by the parent GGA BLYP [65, 95].
Excitation energies and oscillator strengths were calculated via the iterative David-son procedure considering the ten lowest orbital energy differences [96]. Dipole-allowed and dipole-forbidden excitation energies were computed. The tolerance key was set to 10−8 and the orthonormality key to 10−10. It has to be noted that the given tolerance specifies the maximum error between two subsequent iterations in atomic units of the squared excitation energies. Consequently, these settings correspond to a tolerance of 10−4Hartree in the obtained excitation energies.
The keywords used for testing different nonadditive kinetic energy functionals corre-spond to the abbreviations given in Table 2.3, except for the case of the APBEK and revAPBEK functionals. Since these are currently not implemented in ADF explicitly, to use them the PBE2 functional was called and its parameters adapted accordingly.
Symmetry detection was disabled in all calculations. The densities of the frozen envi-ronments were obtained at the same level of theory as described above, but considering only the isolated embedding molecule without the active system being present.
The large number of individual calculations in this study required a high degree of au-tomation. Hence, several small programs in the programming language Python, version 3.5, were written [97]. Extensive use was made of the Python Library for Automat-ing Molecular Simulations (PLAMS) [98]. Furthermore, optimization routines from the “SciPy” library were applied and the “NumPy” package was used [99, 100]. Plots were created with the help of the “matplotlib” library while making occasional use of the package “Pandas” [101, 102].
4.1 Selection of Reference Systems and Excitations
The systems under investigation are molecular dimers consisting of one “active”, organic solute, that is to be excited, and one solvent molecule, that constitutes the embedding environment. It is important to note that the calculations presented in the following do not aim at reproducing experiment, but purely on the comparison of supermolecular DFT and FDE outcomes. To accurately achieve the former, it would be necessary to account for more solvent molecules to properly model the microsolvation environment. Especially for small numbers of solvent molecules a pronounced relation between their quantity and the obtained solvatochromic shift has to be expected [103].
A straightforward approach with regard to choosing the reference set would be to apply one of the commonly used benchmark sets of weakly interacting systems, such as the S22 set of molecular dimers [104]. However, there are two additional requirements with regard to the chosen reference sets:
1. The excitations that shall be reproduced via FDEu have to be accessible by the method, which means, as outlined in Section 2.3.2, that they have to be localized on the active system.5 The contour plots in Figure 4.1 further illustrate this condition.
2. Besides the fact that some of the excitations that are accessible to supermolec-ular TD-DFT calculations cannot be found via FDEu, it is also not guaranteed that even the obtainable excitations remain in the same order with respect to the magnitude of the corresponding excitation energies. Consequently, the simple com-parison of the nthexcitation found via standard TD-DFT and the nthone obtained via FDEu is likely to be misleading. At the same time a human investigation and selection of excitation energies for every spectrum is unfeasible due to the very large number of calculations. That is why additional criteria, which are suitable to be implemented in an algorithm that performs the matching of supermolecular and FDEu excitation energies automatically, have to be found.
Due to the above mentioned reasons, no standard benchmark set was found for the purpose of this project. Instead, the training and test set systems were handpicked from a larger number of trial systems under the condition that the outlined demands can be met. The motivation to investigate a trial system with regard to its suitability for this research was often drawn from previous studies on similar solute-solvent com-binations [8, 104, 105, 106, 107, 108]. With regard to the identification of matching excitations, from a theoretical point of view, one might want to consider a classification in terms of symmetry. However, it is not expected to detect the symmetry groups that
5
In practice, regarding to the actual selection procedure, “localized on the active system” has to be understood as “mainly localized on the active system”.
(a) Charge-transfer transition.
(b) Excitation that is mainly localized on the acetone subsystem.
Figure 4.1: Contour plots of the states composing the main contributions to two different excitations of the acetone−DMSO system calculated with supermolecular TD-DFT. The states of lower energy are depicted on the left. While the transition illustrated in (a) is not accessible to FDEu, the one shown in (b) composes a reasonable reference excitation.
are predicted in FDE results and investigations of the isolated active system in super-molecular calculations. This is also why symmetry was disabled in the computations, as outlined in Section 3. Consequently, this concept cannot be applied here. A scheme to quantitatively assess the similarity between two excitations via the computation of the overlap of the corresponding transition densities was proposed by Neugebauer et al. [8]. While this procedure is attractive, due to the time limitation of the present project, yet a simpler approach is preferable. In a significant number of investigated systems the first (i.e., lowest in energy) singlet-triplet excitation appeared to be a safe match between supermolecular and FDE calculations. Therefore, it was chosen to only include systems, for which the first singlet-triplet excitation appeared to be directly comparable, in the reference sets. During the subsequent calculations the first singlet-triplet excitations were compared without further checks.
In experimental spectra, primarily “allowed”, i.e., singlet-singlet excitations with non-zero oscillator strengths, are observed. This renders them especially interesting also for computational investigations. That is why it was aimed at also including one representa-tive of this group for every system in the set of reference excitation energies. However, for these, no such easy and general correspondence, as for the first singlet-triplet transition, could be observed for all but a few systems. Instead, for each active system, one of the first excitations with non-negligible oscillator strength that was also found in the super-molecular reference calculations was selected. Then criteria with regard to the oscillator strength and the size of the components of the transition dipole moments (in practice that means minimum and/or maximum allowed value) were retrieved that should allow to safely identify it: The code that carries out and analyzes FDEu calculations loops through the computed singlet-singlet excitations, moving from low to high energies, and
Table 4.2: List of the systems included in the training set and the test set.
(a) Training Set
Active System Environment Acetone DMSO Acetone H2O Acetone Cyclohexane Pyridine CCl4 Pyridine H2O Pyridine Acetonitrile Uracil H2O Uracil Acetonitrile Uracil Cyclohexane (b) Test Set
Active System Environment
Acrolein H2O
Acrolein Acetonitrile Acrolein n-Hexane 1-Bromo-4-chlorobenzene Cyclohexane Chlorobenzene Diethyl ether Benzenesulfonic acid H2O Ethene Ethyne Formaldehyde Acetonitrile 8-Hydroxyquinoline H2O Indole Benzene Pyrazine Acetonitrile Pyrazine CCl4 Pyrazine CHCl3 Pyrazine H2O Pyrazine iso-Octane
checks whether they meet the selected conditions. The first one that does so is assumed to be the excitation of interest that is compared to the supermolecular reference. These reference excitation energies were extracted by exactly the same procedure. While this selection algorithm does not provide any guarantee for extracting the correct excitation energy, sample controls indicated that it is likely to be robust.
A list of the training and test set systems6 that were selected manually as described here is provided in Table 4.2. Hereinafter, a supersystem will be referred to with the notation “active subsystem−frozen subsystem”.
4.2 Performance of Known Nonadditive Kinetic Energy
Approximants on the Training Set
Figure 4.3 shows the mean absolute deviations of the excitation energies of the training systems described in Section 4.1 obtained with FDEu and different nonadditive kinetic energy functionals compared to the supermolecular reference. It is striking that the majority of the functionals yield results of very similar quality. Interestingly, this also includes the MAE obtained with the GGA-type functional PW91k by Lembarki and Chermette, which is frequently recommended for subsystem DFT and FDE calculations in the literature, as it has proven to yield high-quality results for other quantities, such
6
Due to the effort in identifying suitable reference systems, the test set was chosen later during the project, to address the points outlined in the introduction of Section 4.4.
Figure 4.3: Mean absolute error obained for FDEu calculations of excitation energies when taking supermolecular TD-DFT results as a reference. The data pre-sented here arises from calculations on the so-called training set given in Table 4.2.
as ground state energies and electronic densities [109, 110, 111]. The functional con-tains comparably many, namely five, empirical parameters.7 It is not surprising that
the PBEn functionals proposed by Karasiev, PBE2, PBE3 and PBE4 (and among these especially PBE2) yield considerably higher errors than the majority of the other func-tionals taking into account the underlying, force-based parametrization scheme outlined in Section 2.5. Also Ernzerhof’s E00 functional performs rather poorly, although still considerably better than PBE2 and PBE4. The most impressive feature of this data is that even a simple LDA functional yields results of an accuracy similar to those achieved with more complex functional forms. The nondecomposable functional NDSD was pre-viously reported to be of an equally good quality as PW91k and is therefore said to be a promising alternative [7, 111]. However, it does not provide any advantage over a simple LDA approach on the data set under investigation.
In addition to the average values discussed so far, also the range of the errors through-out the training set was investigated. Indeed, the standard deviations of the MAEs depicted in Figure 4.3 are of a size similar to the MAEs themselves, indicating a large variance across the training set. To illustrate this feature, the MAEs per excitation en-ergy, with the average taken over all functionals listed in Table 2.3, is shown in Graph 4.4:
7
As already mentioned in Section 2.5, two of the six parameters listed in Table 2.3, following a physical argument, have to yield a fixed sum, so that in fact there is one degree of freedom less than one might assume when investigating the functional form.
Figure 4.4: Absolute deviations of the first singlet-triplet and the selected singlet-singlet excitation energies averaged across the set of investigated known nonadditive kinetic energy approximants. The error bars indicate the standard deviation among the absolute errors. The choice to depict only the positive error bars is solely made to maintain good readability.
While, for example, the first singlet-triplet excitation of the pyridine−H2O dimer can be captured extremely well, the singlet-singlet excitations of the acetone−cyclohexane and pyridine−CCl4complexes proved to be much more challenging. Large standard de-viations of the MAEs depicted in Graph 4.4 indicate that the quality of the calculated excitation energies strongly varies with respect to the chosen nonadditive functional. It should however be noted that this result firmly relies on the functionals included in the test. As can be seen when investigating the individual results provided in Appendix A.1 large standard deviations do not necessarily imply that larger subgroups of the function-als cannot still perform comparably well. With regard to this statement, it might be illustrative to compare Figure 4.4 to its counterpart Graph B.1 in the appendices that depicts exactly the same data but after excluding the results obtained with Karasiev’s PBEn functionals from the analysis: Not considering these functional yields significantly smaller MAEs and standard deviations for many of the investigated excitation energies. For all training set systems but the acetone−H2O and uracil−H2O dimers, the mean
absolute error of the first singlet-triplet excitation is lower than that obtained for the selected singlet-singlet transition. It could be advantageous if one was able to conclude from a proper representation of one excitation energy on the quality of the description of other transitions of the same system. This would allow one to rate the reliability of
Figure 4.5: Absolute error of the FDEu singlet-singlet excitation energies drawn against the error of the singlet-triplet excitation energies of the same system. Each point corresponds to the result obtained with a different nonadditive kinetic energy functional.
the whole spectra obtained with a specific ˜TSnad after controlling only the quality of one excitation energy explicitly. In order to assess possible correlations in this regime, Fig-ure 4.5 displays the absolute deviation in the investigated singlet-singlet excitation drawn against the error observed for the first singlet-triplet transition. One might see such a correlation for some of the dimers, such as the pyridine−CCl4 and the pyridine−H2O systems, but it does not hold for the remaining cases. Furthermore, it should be noted that for the dimers where pyridine constitutes the active system, the orbitals that ac-count for the main part of the involved states are alike, which explains the comparably strong correlation for these systems. Plot 4.5 also gives an insight into how the quality of the computed excitations differs with the functional choice. This is achieved by ob-serving how widely (or narrowly) the data points for one system are distributed within either direction.
So far, the deviations from the FDEu calculation to the supermolecular TD-DFT reference were discussed in terms of absolute errors. In order to rate the overall quality of the results an analysis in terms of relative errors is also desirable. These are given in the third column of Table 4.6. While the relative errors do not provide any new insight in terms of the comparative performance of the different functionals, the obtained values trigger an important question: With average deviations between the supermolecular
