2. Electronic Structure Theories

2.1. Density-Functional Theory

Excited state calculations require a reference ground state calculation within density-functional theory. VOTCA-XTP provides both an automated interface to the ORCA package [Neese:2012] and a lightweight internal DFT engine based on atom-centered Gaussian-type orbitals for method developing and testing. It solves the Kohn-Sham Equations for the molecular orbitals \(\phi_n^\textrm{KS}(\mathbf{r})\) with orbital energies \(\varepsilon_n^\textrm{KS}\)

(2.1)\[\left\{ -\frac{\hbar^2}{2m}\nabla^2 + V_\mathrm{ext}(\mathbf{r}) + V_\textrm{H}(\mathbf{r}) +V_\textrm{xc}(\mathbf{r})\right\}\phi_n^\textrm{KS}(\mathbf{r}) =\varepsilon_n^\textrm{KS} \phi_n^\textrm{KS}(\mathbf{r}) ,\]

where \(V_\textrm{ext}\) is the external potential, \(V_\textrm{H}\) the Hartree potential, and \(V_\textrm{xc}\) the exchange-correlation potential. VOTCA-XTP also contains functionality for projector-based-embedding DFT-in-DFT ground state calculations [Manby:2012], in which a chosen active subregion of a molecular system is embedded into an inactive one, reproducing the total energy of the full system ground state exactly.

For closed-shell systems, Eq.2.1 is solved once for a common set of spatial orbitals doubly occupied by electrons of opposite spin (restricted Kohn-Sham, RKS). Open-shell systems – radicals, cations, anions, or any system with an odd number of electrons or a triplet or higher-multiplicity ground state – instead require separate equations for the \(\alpha\) and \(\beta\) spin channels, each with its own set of orbitals, orbital energies, and occupied-orbital count (unrestricted Kohn-Sham, UKS). The internal DFT engine automatically switches to the UKS formalism whenever the numbers of \(\alpha\) and \(\beta\) electrons differ, as is always the case, for example, for a molecule with an odd total charge or an unpaired electron.

2.1.1. Self-consistent field convergence

Because \(V_\textrm{H}\) and \(V_\textrm{xc}\) in Eq.2.1 themselves depend on the orbitals being solved for, the Kohn-Sham equations must be solved iteratively: a trial density produces a Fock (Kohn-Sham) matrix, which is diagonalized to give new orbitals and a new density, and this cycle repeats until the change in energy and density falls below a set of convergence thresholds. Convergence of this self-consistent field (SCF) procedure is normally accelerated using the direct inversion of the iterative subspace (DIIS) [Pulay:1980] and its energy-based variant (A-DIIS) [Hu:2010], which extrapolate a new Fock matrix guess from a short history of previous iterations rather than using the current iteration’s own Fock matrix directly. In addition, level shifting of the virtual orbitals and linear mixing of the density with that of the previous iteration are available to further stabilize convergence.

Open-shell systems with a small HOMO-LUMO gap, or in which the physically correct electronic state involves a specific localization of charge or spin density that is not obviously the lowest-energy solution to plain SCF, are systematically more difficult to converge with DIIS-based methods alone – a well-documented limitation of DIIS-family extrapolation rather than a peculiarity of any one implementation. For such cases, the internal DFT engine falls back to direct energy minimization in orbital-rotation space once DIIS-based convergence has repeatedly failed, or once the SCF error is observed to be improving too slowly on average over a trailing window of iterations. This fallback solves an augmented-Hessian eigenvalue problem [Helmich-Paris:2021] – the same family of algorithm used, for example, in ORCA’s Trust-Region Augmented-Hessian (TRAH) SCF converger – within a trust-region framework that only accepts a step once it is confirmed, after the fact, to actually lower the energy by an amount consistent with the quadratic model used to propose it. Further implementation details are given in Implementation Notes.

2.2. Constrained Density-Functional Theory

Ordinary Kohn-Sham DFT finds the ground-state density by minimizing the total energy functional \(E[n]\) over all densities integrating to the correct total number of electrons. For many systems of interest – localized charge or spin states relevant to electron- and hole-transfer reactions, diabatic states for computing electronic couplings, or simply a well-defined starting point for systems with a near-degenerate manifold of low-lying electronic configurations – the physically relevant state is not necessarily the unconstrained minimum. Self-interaction error in approximate exchange-correlation functionals in particular tends to systematically favor spurious delocalization of charge or spin density over two or more fragments, even when the physically correct picture is a localized one.

Constrained density-functional theory (CDFT) [Wu:2005] addresses this by minimizing \(E[n]\) subject to an additional constraint: that some population \(N_c\), obtained by integrating the density against a chosen weight function \(w_c(\mathbf{r})\), take on a specified target value,

(2.2)\[\int w_c(\mathbf{r})\, n(\mathbf{r})\, d\mathbf{r} = N_c .\]

The weight function \(w_c(\mathbf{r})\) defines which part of the system the constraint applies to. VOTCA-XTP uses Hirshfeld (stockholder) partitioning [Hirshfeld:1977], in which each atom \(A\) is assigned a weight

(2.3)\[w_A(\mathbf{r}) = \frac{n_A^{\textrm{at}}(\mathbf{r})}{\sum_B n_B^{\textrm{at}}(\mathbf{r})}\]

built from a superposition of free-atom reference densities \(n_A^{\textrm{at}}\), and the constraint weight for a group (fragment) of atoms is the sum of the weights of its constituent atoms. This has been found to give a more physically faithful charge partitioning than simpler, purely geometric schemes for the class of charge- and spin-localization problems CDFT is typically used for [Ahart:2022].

The constraint in Eq.2.2 is enforced via a Lagrange multiplier \(V_c\), leading to the constrained functional

(2.4)\[W[n, V_c] = E[n] + V_c\left(\int w_c(\mathbf{r})\, n(\mathbf{r})\, d\mathbf{r} - N_c\right) .\]

The constrained ground state is the density that minimizes \(W\) for the multiplier value \(V_c\) that itself maximizes \(W\) – a saddle-point problem, though one with a unique solution [Wu:2005]. In practice, this is solved as two nested loops: an inner Kohn-Sham SCF loop for fixed \(V_c\), in which the constraint enters as an additional, weight-function-dependent potential added directly to the Kohn-Sham Hamiltonian in Eq.2.1,

(2.5)\[V_\textrm{c}(\mathbf{r}) = V_c\, w_c(\mathbf{r}) ,\]

and an outer bisection loop over \(V_c\) itself, adjusting the multiplier until the population in Eq.2.2 matches the requested target. Since the inner SCF for a nearby value of \(V_c\) typically provides a good starting guess for the next, the outer loop warm-starts each new inner SCF from the orbitals of the previous one.

Because CDFT constrains a population defined via a weight function that depends on the nuclear positions (both directly, through the definition of the atomic fragments, and indirectly, through the free-atom reference densities entering Eq.2.3), the force acting on the nuclei is not given by the ordinary Hellmann-Feynman expression alone. Differentiating the constrained functional with respect to a nuclear coordinate yields additional terms describing how the weight function itself, the associated real-space integration grid, and the basis functions centered on the moving atoms all change with the nuclear positions [Oberhofer:2009]. VOTCA-XTP evaluates these terms explicitly, so that CDFT forces are available for geometry optimization and molecular dynamics on the constrained potential energy surface, not just single-point energies.

2.3. Many-Body Green’s Functions and the Bethe-Salpeter Equation

Using the ground-state reference, many-body Green’s functions theory with the \(GW\) approximation first calculates single-particle excitations (electron addition or removal) as solutions to the quasiparticle equations

(2.6)\[\left\{ -\frac{\hbar^2}{2m}\nabla^2 + V_\textrm{ext}(\mathbf{r}) + V_\textrm{H}(\mathbf{r})\right\}\phi_n^\textrm{QP}(\mathbf{r}) + \int{\Sigma(\mathbf{r},\mathbf{r}',\varepsilon_n^\textrm{QP})\phi_n^\textrm{QP}(\mathbf{r}')d\mathbf{r}'} = \varepsilon_n^\textrm{QP} \phi_n^\textrm{QP}(\mathbf{r}) .\]

In place of the exchange-correlation potential in Eq.2.1, the energy-dependent self-energy operator \(\Sigma(\mathbf{r},\mathbf{r}',E)\) occurs in the QP equations. This operator is evaluated using the one-body Green’s function in quasi-particle approximation

(2.7)\[G(\mathbf{r},\mathbf{r}',\omega) = \sum_n{\frac{\phi_n(\mathbf{r})\phi_n^*(\mathbf{r}')}{\omega-\varepsilon_n+i0^+\textrm{sgn}(\varepsilon_n -\mu)}}\]

as

(2.8)\[\Sigma(\mathbf{r},\mathbf{r}',E) = \frac{i}{2\pi} \int{e^{-i\omega 0^+}G(\mathbf{r},\mathbf{r}',E-\omega)W(\mathbf{r},\mathbf{r}',\omega)\,d\omega},\]

where \(W\) denotes the dynamically screened Coulomb interaction. Assuming that \(\phi^\textrm{QP}_n\approx \phi^\textrm{KS}_n\), the quasiparticle energies can be evaluated perturbatively according to

(2.9)\[\varepsilon_n^\textrm{QP}= \varepsilon_n^\textrm{KS} + \Delta \varepsilon_n^{GW} = \varepsilon_n^\textrm{KS} + \left\langle\phi^\textrm{KS}_n\left\vert \Sigma(\varepsilon_n^\textrm{QP})-V_\text{xc} \right\vert\phi^\textrm{KS}_n\right\rangle .\]

As the correction \(\Delta \varepsilon_n^{GW}\) itself depends on \(\varepsilon_n^\textrm{QP}\), Eq.2.9 needs to be solved self-consistently.

Neutral excitations with a conserved number of electrons can be obtained from the Bethe-Salpeter Equation (BSE) by expressing coupled electron-hole amplitudes of excitation \(S\) in a product basis of single-particle orbitals, i.e.,

(2.10)\[\chi_S(\mathbf{r}_\textrm{e},\mathbf{r}_\textrm{h})=\sum_{v}^{\mathrm{occ}}\sum_c^{\mathrm{unocc}}A_{vc}^S\phi_{c}(\mathbf{r}_\textrm{e})\phi^*_{v}(\mathbf{r}_\textrm{h})+B_{vc}^S\phi_{v}(\mathbf{r}_\textrm{e})\phi^{*}_{c}(\mathbf{r}_\textrm{h}),\]

where \(\mathbf{r}_\textrm{e}\) (\(\mathbf{r}_\textrm{h}\)) is for the electron (hole) coordinate and \(A_{vc}\) (\(B_{vc}\)) are the expansion coefficients of the excited state wave function in terms of resonant (anti-resonant) transitions between occupied \(v\) and unoccupied \(c\) states, respectively. In this basis, the BSE turns into an effective two-particle Hamiltonian problem of the form

(2.11)\[\begin{split}\begin{pmatrix} \underline{\mathbf{H}}^{\text{res}}&\underline{\mathbf{K}} \\ -\underline{\mathbf{K}} & -\underline{\mathbf{H}}^{\text{res}} \end{pmatrix} \begin{pmatrix} \mathbf{A}^S\\ \mathbf{B}^S \end{pmatrix} =\Omega_S \begin{pmatrix} \mathbf{A}^S\\ \mathbf{B}^S \end{pmatrix}.\end{split}\]

Specifically, the matrix elements of the blocks \(\underline{\mathbf{H}}^{\text{res}}\) and \(\underline{\mathbf{K}}\) are calculated as

(2.12)\[\begin{split}\begin{align} H^{\text{res}}_{vc,v'c'}&=D_{vc,v'c'}+\eta K^\mathrm{x}_{vc,v'c'}+K^\mathrm{d}_{vc,v'c'}\\ K_{cv,v'c'}&=\eta K^\mathrm{x}_{cv,v'c'}+K^\mathrm{d}_{cv,v'c'}\, , \end{align}\end{split}\]

with

(2.13)\[\begin{split}\begin{align} D_{vc,v'c'}&=(\varepsilon_c-\varepsilon_v)\delta_{vv'}\delta_{cc'},\\ K^\text{x}_{vc,v'c'}&=\iint \phi_c^*(\mathbf{r}_\textrm{e})\phi_v(\mathbf{r}_\textrm{e})v_{\mathrm{C}}(\mathbf{r}_\textrm{e},\mathbf{r}_\textrm{h}) \phi_{c'}(\mathbf{r}_\textrm{h})\phi_{v'}^*(\mathbf{r}_\textrm{h}) d^3\mathbf{r}_\textrm{e} d^3\mathbf{r}_\textrm{h}\\ K^\text{d}_{vc,v'c'}&=-\iint \phi_c^*(\mathbf{r}_\textrm{e})\phi_{c'}(\mathbf{r}_\textrm{e})W(\mathbf{r}_\textrm{e},\mathbf{r}_\textrm{h},\omega=0) \phi_v(\mathbf{r}_\textrm{h})\phi_{v'}^*(\mathbf{r}_\textrm{h})d^3\mathbf{r}_\textrm{e} d^3\mathbf{r}_\textrm{h} \, . \end{align}\end{split}\]

and \(\eta=2\) (\(\eta=0\)) for singlet (triplet) excitations. Here, \(K^\text{x}\) is the repulsive exchange interaction originating from the bare Coulomb term \(v_\mathrm{C}\), while the direct interaction \(K^\text{d}\) contains the attractive, but screened, interaction \(W\) between electron and hole, causing the binding of the electron-hole pair. In Eq.2.13 it is assumed that the dynamic properties of \(W(\omega)\) are negligible, and the computationally less demanding static approximation \(\omega=0\) is employed.