Hamiltonian

Octopus is based upon Density Functional Theory in the Kohn-Sham formulation. The Kohn-Sham Hamiltonian is the main part of this formulation; in this section we describe how the Hamiltonian is treated in Octopus and what are the variables that control it.

Introduction

Although the Hohenberg-Kohn theorem states that DFT is exact, the Kohn-Sham method of reducing an interacting many-particle problem to a non-interacting single-particle problem introduces an approximation: the exchange-correlation term.

Ground-State DFT

The Kohn-Sham method of DFT assumes that, for each interacting ground-state density $n(r)$, there exists a non-interacting electron system with the same ground-state density. The interacting ground state is obtainable through the solution of the Kohn-Sham equations

$$ \left[-\frac{\nabla^2}{2}+v_{\rm KS}[n](r)\right]\varphi_i(r) = \varepsilon_i \varphi_i(r) $$

The notation $v_{\rm KS}[n]$ means that the Kohn-Sham potential, $v_{\rm KS}$, has a functional dependence on $n$, the electronic density, which is defined in terms of the Kohn-Sham wave-functions by

$$ n(r) = \sum_i^{\rm occ} | \varphi_i(r) |^2\;. $$

The potential $v_{\rm KS}$ is defined as the sum of the external potential (normally the potential generated by the nuclei), the Hartree term, and the exchange-correlation (xc) potential

$$ v_{\rm KS}[n](r) = v_{\rm ext}(r) + v_{\rm Hartree}[n](r) + v_{\rm xc}[n](r)\;. $$

Due to the functional dependence on the density, these equations form a set of nonlinear coupled equations. The standard procedure to solve it is iterating until self-consistency is achieved.

The total energy of the electronic system is given by

$$ E_{\rm elec}[n] = T_s[n] + \int\!\!d^3r\;n(r)v_{\rm ext}(r) + E_{\rm Hartree}[n] + E_{\rm x}[n] + E_{\rm c}[n] $$

where $T_s[n]$ is the non-interacting kinetic energy, $v_{\rm ext}$ is the external potential, and $E_{\rm x,c}[n]$ are the exchange (x) and correlation (c) energies. The second term is called the “external energy” in the code. In practice, from the solution of the Kohn-Sham equations, we can evaluate the energy via the eigenvalues as

$$ E_{\rm elec}[n] = \sum \epsilon_i - \int\!\!d^3r\;n(r)v_{\rm xc}(r) - E_{\rm Hartree}[n] + E_{\rm x}[n] + E_{\rm c}[n] $$

where $v_{\rm xc}[n]$ is the exchange-correlation potential. To find the total energy of the entire system, we additionally include ion-ion interaction and ionic kinetic energy.

$E_{\rm xc} = E_{\rm x} + E_{\rm c}$ is an unknown object and includes all the non-trivial many-body effects required to make KS theory exact. Several approximations to $E_{\rm xc}$ have been proposed. The most used is the local density approximation (LDA). In this approximation $E_{\rm xc}[n(r)]$ is taken to be the exchange and correlation energy of a homogeneous electron gas with density $n = n(r)$. Although there exists an exact expression for the exchange energy in this model, the exact value of the correlation energy is known only in the limit of very high densities. Ceperley and Alder did a Monte Carlo simulation of the homogeneous electron gas at several densities. Several parameterizations of the correlation energy for any density were then obtained interpolating the Monte Carlo results. One particularly simple parameterization was proposed by Perdew and Zunger, and this option may be used by Octopus. You can, of course, choose other xc functionals (XCFunctional), via the extensive library provided by libxc.

Time-dependent DFT

Time-dependent density-functional theory (TDDFT) extends the basic ideas of ground-state density-functional theory (DFT) to the treatment of excitations and of more general time-dependent phenomena. TDDFT can be viewed as an alternative formulation of time-dependent quantum mechanics but, in contrast to the normal approach that relies on wave-functions and on the many-body Schrödinger equation, its basic variable is the one-body electron density, $n(r,t)$. The advantages are clear: The many-body wave-function, a function in a $3N$-dimensional space (where $N$ is the number of electrons in the system), is a very complex mathematical object, while the density is a simple function that depends solely on the 3-dimensional vector $r$. The standard way to obtain $n(r,t)$ is with the help of a fictitious system of non-interacting electrons, the Kohn-Sham system. The final equations are simple to tackle numerically, and are routinely solved for systems with a large number of atoms. These electrons feel an effective potential, the time-dependent Kohn-Sham potential. The exact form of this potential is unknown, and has therefore to be approximated.

The time-dependent Kohn-Sham equations are

$$ i \frac{\partial}{\partial t}\varphi_i(r,t) = \left[ -\frac{\nabla^2}{2} + v_{\rm KS}(r,t) \right]\varphi_i(r,t) $$

The density of the interacting system can be obtained from the time-dependent Kohn-Sham orbitals

$$ n(r, t) = \sum_i^{\rm occ} \left|\varphi_i(r,t)\right|^2 \;. $$

The time-dependent Kohn-Sham equations, having the form of a one-particle equation, is fairly easy to solve numerically. We stress, however, that the Kohn-Sham equation ‘‘is not’’ a mean-field approximation: If we knew the exact Kohn-Sham potential, $v_{\rm KS}$, we would obtain the exact Kohn-Sham orbitals, and from these the correct density of the system. The Kohn-Sham potential is conventionally separated in the following way

$$ v_{\rm KS}(r,t) = v_{\rm ext}(r,t) + v_{\rm Hartree}(r,t) + v_{\rm xc}(r,t) \,. $$

The first term is again the external potential. The Hartree potential accounts for the classical electrostatic interaction between the electrons

$$ v_{\rm Hartree}(r,t) = \int\!\!d^3r^{\prime} \frac{n(r,t)}{\left|r-r^{\prime}\right|} \;. $$

The time-dependence of the exchange and correlation potential introduces the need for an approximation beyond the one made in the time-independent case. The simplest method of obtaining a time-dependent xc potential consists in assuming that the potential is the time-independent xc potential evaluated at the time-dependent density, i.e.,

$$ v^{\rm adiabatic}_{\rm xc}(r, t) = \left.\tilde v_{\rm xc}[n](r)\right|_{n=n(t)} \;, $$

This is called the adiabatic approximation. If the time-independent xc potential chosen is the LDA, then we obtain the so-called adiabatic local density approximation (ALDA). This approximation gives remarkably good excitation energies but suffers from the same problems as the LDA, most notably the exponential fall-off of the xc potential. If a strong laser pushes the electrons to regions far from the nucleus, ALDA should not be expected to give an accurate description of the system. Other options for the time-dependent xc potential are orbital-dependent potentials like the exact exchange functional (EXX) (usually in the Krieger-Li-Iafrate (KLI) approximation).

Occupations

You can set the occupations of the orbitals by hand, using the variable Occupations. Alternatively, one can specify a smearing of the occupations using the variables Smearing , and SmearingFunction .

External Potential

You can add an external uniform and constant electric or magnetic field, set by the StaticElectricField or StaticMagneticField. If you want to add a more general potential, you can do it using a user-defined Species.

Electron-Electron Interaction

You can neglect this term by setting TheoryLevel = independent_particles variable. This implies that both the Hartree and exchange-correlation terms will not be calculated.

Exchange and correlation potential

Octopus does not distinguish, for the moment, between ground-state DFT xc functionals, and time-dependent DFT functionals. In other words, in all cases the adiabatic approximation is assumed. In mathematical terms, we may formalize this in the following way: let $n$ be the time-dependent density, a function living in the four-dimensional space-time world. We will call $n_t(\vec{r}) = n(\vec{r}, t),$, the electronic density at time $t,$, a function in the three-dimensional space. An exchange and/or correlation energy functional $E^{\alpha},$ in ground state DFT may then be used to build an adiabatic exchange and/or correlation action functional in the following way:

$$ A^{\alpha}[n] = \int_{t_0}^{t_f}{\rm d}\tau E^{\alpha}[n_\tau]\,, $$

The time-dependent potential functional is then:

$$ v^{\alpha}[n](\vec{r},t) \equiv { \delta A^{\alpha} \over \delta n(\vec{r}, t)} = { \delta E^{\alpha} \over \delta n_{t}(\vec{r}) } = v^{\alpha \rm (GS)}[n_t](\vec{r})\,. $$

We use the distinct notation $v^{\alpha}[n](\vec{r},t)\,$ and $v^{\alpha \rm (GS)}[n_t](\vec{r}),$ to stress that the exchange and correlation potential in TDDFT – the former – and the exchange and correlation potential in GS-DFT – the latter – are in principle two conceptually different objects, which coincide only thanks to the adiabatic approximation. This is the reason why we may actually only refer to the functionals in the GS-DFT context.

We may classify the xc functionals contained in the Octopus code following John Perdew’s Jacob’s Ladder scheme:

• LDA rung: A functional $\alpha,$ belonging to the LDA rung depends only on the electronic density (on the spin density in spin-polarized or spinors cases). Moreover, it has a local dependency on the density, i.e.:

$$ E^{\alpha}_{\rm LDA} = E^{\alpha}_{\rm LDA}[n] = \int {\rm d}^3r e^{\alpha}_{\rm LDA}(n(\vec{r}))\,. $$

The potential may be then derived by functional derivation:

$$ v^{\alpha}_{LDA}[n](\vec{r}) = {\delta {E^{\alpha}_{\rm LDA}} \over {\delta n(\vec{r})}} = {{{\rm d}e^{\alpha}_{\rm LDA} } \over {{\rm d}n }}(n(\vec{r}))\,. $$

• GGA rung: A functional $\alpha,$ belonging to the GGA rung depends on the electronic density, and also on its gradient. Moreover, it also has a local dependency (actually, the GGA is very often called a ‘‘semi-local’’ functional due to this).

$$ E^{\rm \alpha (GGA)} = E^{\rm \alpha (GGA)}[n, \vec{\nabla}n] = \int {\rm d}^3r e^{\alpha}_{\rm GGA}(n(\vec{r}), \vec{\nabla}n(\vec{r}))\,. $$

• meta-GGA rung: The next rung, the of of meta-GGAs, contains semi-local functionals that not only depends on the electronic density, its gradient or its Laplacian, but only on the kinetic energy density ($\tau$).

$$ E^{\rm \alpha (mGGA)} = E^{\rm \alpha (mGGA)}[n, \vec{\nabla}n, \Delta n, \tau] = \int {\rm d}^3r e^{\alpha}_{\rm mGGA}(n(\vec{r}), \vec{\nabla}n(\vec{r}), \Delta n(\vec{r}), \tau(\vec{r}))\,. $$

Note that in order to control if one uses the U(1) gauge-invariant kinetic energy-density or not, one case use the variable XCUseGaugeIndependentKED. By default, the mGGAs are treated within the generalized Kohn-Sham approach. However, it is also possible to use the OEP approach, see below.

• OEP rung: This is the family of functionals which are defined in terms of the occupied Kohn-Sham orbitals, $\lbrace \varphi_i \rbrace_{i=1}^{N},$. These are in fact the only non-local functionals. The name of the rung, OEP, stands for “optimized-effective-potential”, the reason being that in general this is the method used to derive the potential from the energy functional (direct functional derivation is in this case not possible). A more suitable name would be orbital-dependent functionals. In order to use mGGAs using the OEP framework, plus set TheoryLevel=kohn_sham.

• Hybrid rung: The rung of hybrid functionals refer to the family of functionals that include a fraction of the Fock operator, or a fraction of a range-separated part of it. Note that in the case of solids, a special treatment of the Coulomb singularity is needed. This is activated by default and is controled by the variables HFSingularity, HFSingularityNk, and HFSingularityNsteps.

Octopus comes with several Exchange and Correlation potentials, including several flavours of LDA, GGA and OEP. You can choose your favorites by setting the variable XCFunctional. (Note that until Octopus <= 2.1 the exchange-correlation functional was chosen with the variables XFunctional and CFunctional.)

When using OEP Exchange, the variable OEPLevel controls the level of approximation required.

You can also include Self-Interaction Correction (SIC), controlled by the variable SICCorrection. Note that SIC is only treated using the OEP approach.

Octopus also supports DFT+U and some extension of it. See the corresponding tutorials for me details.

van der Walls interaction

Common semi-local density functional approximations do not capture the long-range tail of the dispersion energy, requiring one to use more advanced functionals, or an additive corrective scheme. Octopus supports different flavors of additive corrections for van der Walls interactions, using DFT-D3, libvdwxc, and a internal implementation of the Tkatchenko-Scheffler scheme. See the variable VDWCorrection, and the related variables VDW_TS_cutoff, VDW_TS_damping, VDW_TS_sr, VDWD3Functional, and VDWSelfConsistent.

DFT-D2

DFT-D2 (also supplied by the DFT-D3 library) approximates London-dispersion as a pairwise function of the atomic separation, $|\mathbf{R}{AB}|$, accounting for the attractive dispersion forces between atoms. The term uses fixed dispersion coefficients, $C_6^{AB}$, and a damping function, $f{\text {damp }}^{\mathrm{DFT} 2}$, to interpolate between the long- and short-range limits^vdW-1:

$$ E_{\text {disp }}^{\mathrm{DFT}-\mathrm{D} 2}=-\frac{1}{2} s_6 \sum_{A \neq B} \frac{C_6^{A B}}{R_{A B}^6} f_{\text {damp }}^{\mathrm{DFT} 2}\left(R_{AB}\right), $$

where $s_6$ is a global scaling parameter, and the prefactor of $1/2$ prevents double-counting. DFT-D2 is well-suited to the first 34 main-group elements but has limited applicability to 3d and 4d metals^vdW-2.

DFT-D3

DFT-D3 improves on D2 by making the dispersion coefficients dependent on the local atomic environment, via the atomic coordination number. It is applicable to the first 94 elements of the periodic table, with main-group elements and transition metals treated on an equal footing. The most important two-body term is given at long range by:

$$ E_{disp} = - \frac{1}{2} \sum_{A \neq B} \sum_{n=6,8} s_n \frac{C_n^{AB}}{R^n_{AB}}, $$

where $C^{AB}n$ denotes the averaged (isotropic) nth-order dispersion coefficient for atom pair $AB$, and $s_n$ is a functional-dependent scaling factor. In order to avoid near singularities for small distances ($|R{AB}| \rightarrow 0$), the dispersion contribution needs to be damped at short distances. DFT-D3 offers a couple of options for damping.

One can damp the dispersion contribution to zero for short ranges:

$$ E_{disp} = - \frac{1}{2} \sum_{A \neq B} \sum_{n=6,8} s_n \frac{C_n^{AB}}{R^n_{AB}} f_{d,n}(R^n_{AB}). $$

Choosing VDWCorrection $=$ VDW_D3_no_damping uses the damping function^vdW-3:

$$ f_{d, n}=\frac{1}{1+6\left(R_{AB} /\left(s_{r, n} R^0_{A B}\right)\right)^{-\alpha_n}}, $$

and choosing VDWCorrection $=$ VDW_D3M_no_damping uses the modified damping function of Sherrill and coworkers^vdW-4:

$$ f_{d, n}=\frac{1}{1+6\left(R_{AB} /\left(s_{r, n} R^0_{A B}\right) + R^0_{A B} \beta \right)^{-\alpha_n}}, $$

the latter of which improves the short-range behavior of the method. Rational damping was proposed by Becke and Johnson^vdW-5:

$$ E_{disp} = - \frac{1}{2} \sum {A \neq B} \sum_{n=6,8} s_n \frac{C_n^{AB}}{R^n_{AB} + f(R^0_{AB})^n}, $$

with:

$$ R^0_{AB} =\sqrt{\frac{C_8^{A B}}{C_6^{A B}}}, \ f\left(R^0_{AB}\right) =a_1 R^0_{AB} + a_2. $$

This can be used by choosing VDWCorrection $=$ VDW_D3_BJ or VDW_D3M_BJ. The modified BJ damping of Sherrill and coworkers uses the same functional form but refits parameters $s_8$, $a_1$ and $a_2$. In general, the library authors recommend the use of Becke-Johnson (BJ)-damping, because it does not lead to artificial repulsive forces. Further details of the DFT-D3 program can be found here.

Relativistic Corrections

The variable RelativisticCorrection allows one to choose the relativistic correction to be used. Several options are possible. One option is spin-orbit coupling (RelativisticCorrection = spin_orbit), based on relativistic pseudopotentials, see below. The other approach, that can be used for all-electron calculations as well as pseudopotential ones, is based on the zero-th order regular approximation (ZORA). It comes in two flavors, one is the scalar relativitic ZORA (RelativisticCorrection = scalar_relativistic_zora), and the second one is the fully-relativistic ZORA (RelativisticCorrection = fully_relativistic_zora), that includes both scalar relativistic and spin-orbit terms.

Spin-orbit Coupling with fully-relativistic pseudopotentials

The spin-orbit coupling, as it is implemented in Octopus, is included in the relativistic pseudo-potentials. These can either be HGH pseudo-potentials or Troullier-Martins-like pseudo-potentials. In the latter case the pseudo-potentials need to be generated from fully relativistic calculations and their fully separable form is given by: $$ \hat{v}_{KB}=v_{local}+\sum_{l,j,m_j} \frac{|\varphi^{PP}_{ljm_j}\delta v_{lj}^{PP}><\varphi^{PP}_{ljm_j}\delta v_{lj}^{PP}|}{<\varphi^{PP}_{ljm_j}| \delta v_{lj}^{PP}|\varphi^{PP}_{ljm_j}>} $$

Since the angular part of the pseudo wave-functions $\varphi^{PP}_{ljm_j} $ are spherical spinors the wave-functions should be complex spinors and so the SpinComponents needs to be set to non_collinear. This is also true for HGH pseudo-potentials.

Note that currently Octopus is only able to read j-dependent Troullier-Martins-like pseudo-potentials that are provided in the UPF1 and UPF2 file format.