Background on energy, forces and stress

A number of different ways of formulating the energy exist in Conquest at the moment, involving both self-consistent and non-self-consistent densities and potentials, both with and without the neutral atom potential. With self-consistency, all formulations should give the same result, though numerical issues may give small differences; without self-consistency the Harris-Foulkes functional is more accurate.

Self-consistent calculations

We define the energy in Conquest in two ways that are equivalent at the self-consistent ground state. The Harris-Foulkes energy is given as:

\[E_{HF} = 2\mathrm{Tr}\left[KH\right] + \Delta E_{Ha} + \Delta E_{XC} + E_{II}\]

where the first term is the band structure energy, equivalent to the sum over the energies of the occupied states, the second two terms compensate for double counting and the final term gives the ion-ion interaction:

\[E_{II} = \frac{1}{2}\left( \sum_{ij} \frac{Z_i Z_j}{\mid \mathbf{R}_i - \mathbf{R}_j \mid} \right)\]

The Hamiltonian is defined as:

\[\hat{H} = \hat{T} + \hat{V}_{L} + \hat{V}_{NL} + V_{Ha} + V_{XC}\]

where the operators are the kinetic energy, the local and non-local pseudopotentials, the Hartree potential, defined as \(V_{Ha} = \int d\mathbf{r}^\prime n(\mathbf{r}^\prime)/\mid \mathbf{r} - \mathbf{r}^\prime\mid\), and the exchange-correlation potential. The alternative form, often known as the DFT energy, is:

\[E_{DFT} = 2\mathrm{Tr}\left[K(T + V_{L} + V_{NL})\right] + E_{Ha} + E_{XC} + E_{II}\]

with the Hartree energy defined as usual:

\[E_{Ha} = \frac{1}{2}\int\int d\mathbf{r}d\mathbf{r}^{\prime} \frac{n(\mathbf{r})n(\mathbf{r}^\prime)}{\mid \mathbf{r} - \mathbf{r}^\prime\mid}\]

along with the exchange-correlation energy:

\[E_{XC} = \int d\mathbf{r} \epsilon_{XC}\left[n\right] n(\mathbf{r})\]

For the Harris-Foulkes and DFT energies to be equal, it is easy to see that the double counting correction terms in the Harris-Foulkes formalism must be:

\[\Delta E_{Ha} = -E_{Ha} = -\frac{1}{2}\int\int d\mathbf{r}d\mathbf{r}^{\prime} \frac{n(\mathbf{r})n(\mathbf{r}^\prime)}{\mid \mathbf{r} - \mathbf{r}^\prime\mid}\]

and

\[\Delta E_{XC} = \int d\mathbf{r} \left(\epsilon_{XC}[n] - V_{XC}[n]\right)n(\mathbf{r})\]

When calculating forces and stress with self-consistency, we generally use the differentials of the DFT energy rather than the Harris-Foulkes energy; this enables us to separate contributions that are calculated in different ways (in particular on those that are calculated on the integration grid from those that are not).

Go to top.

Neutral atom potential

In a DFT code using local orbitals as basis functions, the total energy is most conveniently written in terms of the interaction of neutral atoms: this is simply a reformulation of the total energy which, in particular, reduces the ion-ion interaction to a sum over short-range pair-wise interactions. The charge density of interest is now the difference between the total charge density and a superposition of atomic densities, notated as \(\delta n(\mathbf{r}) = n(\mathbf{r}) - \sum_i n_i(\mathbf{r})\) for atomic densities \(n_i(\mathbf{r})\). We write:

\[E_{DFT, NA} = 2\mathrm{Tr}\left[K(T + V_{NA} + V_{NL})\right] + E_{\delta Ha} + E_{XC} + E_{SII}\]

where the second term is defined as:

\[E_{\delta Ha} = \frac{1}{2}\int\int d\mathbf{r}d\mathbf{r}^{\prime} \frac{\delta n(\mathbf{r})\delta n(\mathbf{r}^\prime)}{\mid \mathbf{r} - \mathbf{r}^\prime\mid}\]

The final term, the screened ion-ion interaction, is short-ranged, and written as:

\[E_{SII} = \frac{1}{2}\left( \sum_{ij} \frac{Z_i Z_j}{\mid \mathbf{R}_i - \mathbf{R}_j \mid} - \int d\mathbf{r} n_i(\mathbf{r})V_{Ha,j}(\mathbf{r}) \right)\]

where \(V_{Ha,i}(\mathbf{r})\) is the Hartree potential from the atomic density \(n_i(\mathbf{r})\). We define the neutral atom potential for an atom as \(V_{NA,i}(\mathbf{r}) = V_{L,i}(\mathbf{r}) + V_{Ha,i}(\mathbf{r})\), combining the local potential and the Hartree potential for the atomic density; the overall neutral atom potential is given as the sum over the atomic densities, \(V_{NA}(\mathbf{r}) = \sum_i V_{NA,i}(\mathbf{r})\). If we write the pseudo-atomic density as \(n_{PAD}(\mathbf{r}) = \sum_i n_i(\mathbf{r})\) then we can also write \(V_{NA}(\mathbf{r}) = V_L(\mathbf{r}) + V_{Ha, PAD}(\mathbf{r})\).

In this case, we can write the Harris-Foulkes energy as:

\[E_{HF} = 2\mathrm{Tr}\left[KH\right] + \Delta E_{Ha} + \Delta E_{XC} + E_{SII}\]

with the Hamiltonian defined as:

\[\hat{H} = \hat{T} + \hat{V}_{NA} + \hat{V}_{NL} + V_{\delta Ha} + V_{XC}\]

where \(V_{\delta Ha}(\mathbf{r}) = \int d\mathbf{r^\prime} \delta n(\mathbf{r^\prime})/\mid \mathbf{r} - \mathbf{r}^\prime\mid\). Accordingly, the double counting Hartree correction term has to change:

\[\Delta E_{Ha} = -E_{\delta Ha} - \int d\mathbf{r} \delta n(\mathbf{r})\sum_i V_{Ha,i}(\mathbf{r}).\]

Go to top.

Non-self-consistent calculations

In non-self-consistent calculations, we use the Harris-Foulkes functional, along with a reasonable guess for the input density, which is normally taken as the superposition of atomic densities, \(n_{in}(\mathbf{r})\) and write:

\[E_{NSC} = 2\mathrm{Tr}\left[KH\right] + \Delta E_{Ha}\left[n_{in}\right] + \Delta E_{XC}\left[n_{in}\right] + E_{II}\]

Notice that we effectively have two densities being used here: \(n_{in}\) (which is normally the superposition of atomic densities used in the neutral atom case) and and effective output density, \(n_{out} = \sum_{ij} \phi_i K_{ij} \phi_j\) which comes from the band energy (first term); this complicates the calculation of forces and stress compared to the self-consistent case, as we have to consider contributions from both densities.

For the neutral atom potential, \(\delta n(\mathbf{r}) = 0\) by definition, which also means that \(E_{\delta Ha} = 0\) and \(\Delta E_{Ha}=0\).

Go to top.

Partial core corrections

Also known as non-linear core corrections, partial core corrections (PCC) [EFS1] add a model core charge to the pseudopotential to allow for the non-linear exchange-correlation interation between core and valence charge (which is linearised in standard pseudopotentials); this generally improves the accuracy of the pseudopotential. The exchange-correlation potential is evaluated in terms of the combined charge density, \(n_v(\mathbf{r}) + n_c(\mathbf{r})\) where the valence charge is input or output charge density defined above: \(V_{XC}\left[ n_v + n_c \right]\). The exchange-correlation energy becomes:

\[E_{XC} = \int d\mathbf{r} \left(n_v(\mathbf{r}) + n_c(\mathbf{r})\right) V_{XC}\left[ n_v + n_c \right] .\]

Once this change to the charge density has been made, there is no change to the DFT energy. However, the double counting term for Harris-Foulkes needs redefining, since XC contribution to the band energy is \(2Tr[KV_{XC}] = \int d\mathbf{r} n_v(\mathbf{r}) V_{XC}[n_v + n_c]\). We write:

\[\begin{split}\Delta E_{XC} &=& \int d\mathbf{r} \left(n_v(\mathbf{r}) + n_c(\mathbf{r})\right)\epsilon_{XC}[n_v + n_c] - \int d\mathbf{r} n_v(\mathbf{r})V_{XC}[n_v + n_c]\\ &=& \int d\mathbf{r} n_c(\mathbf{r})\epsilon_{XC}[n_v + n_c] + \int d\mathbf{r} \left(\epsilon_{XC}[n_v + n_c] - V_{XC}[n_v + n_c]\right)n_v(\mathbf{r})\end{split}\]

There is an extra factor of \(\int d\mathbf{r} n_c(\mathbf{r})\epsilon_{XC}[n_v + n_c]\) over and above the usual term.

Go to top.

Forces and Stresses

It is important that the forces and stresses be the exact derivatives of the energy, for consistency. In particular, this means that as the energy is calculated in different ways for different contributions, the force or stress contribution must be calculated in the same way.

Go to top.

Forces

Forces are defined as the change in energy with respect to atomic positions; as the basis functions move with the atoms, these changes will also include Pulay terms. The forces found in Conquest are documented extensively elsewhere [EFS2, EFS3] though the changes needed to account for PCC, particularly in the non-self-consistent case, have not been published and are given here for completeness. As well as the Hellmann-Feynman forces (which come from the movement of the local and non-local pseudopotentials with the atoms) we define Pulay forces (divided into two parts, known as phi-Pulay which come from changes in the Hamiltonian matrix, and S-Pulay, which come from changes in the overlap matrix; the phi-Pulay forces are calculated in three contributions, which depend on how the respective parts of the Hamiltonian matrix are calculated: the kinetic energy; the non-local pseudopotential; and the remaining terms which are all found on the integration grid). The ion-ion interactions also contribute forces.

The inclusion of PCC adds an extra term to the forces in all calculations, which comes from the change of the core density as the atoms move; the force on atom \(i\) is given as:

\[\mathbf{F}^{PCC}_i = -\int d\mathbf{r} \nabla_i n^c_i(\mathbf{r}) V_{XC}[n_v + n_c]\]

If the non-self-consistent formalism is used, then a further term is added (the non-self-consistent force changes) to include the gradient of the core charge. The non-self-consistent force is now written as:

\[\mathbf{F}^{NSC}_i = -\int d\mathbf{r} V_{\delta Ha}(\mathbf{r}) \nabla_i n^v_i(\mathbf{r}) - \int d\mathbf{r} \delta n(\mathbf{r}) V_{XC}^\prime\left[n^{in}_v + n_c\right] \left( \nabla_i n^v_i(\mathbf{r}) + \nabla_i n^c_i(\mathbf{r}) \right)\]

where \(V^\prime_{XC}\) is the derivative of the exchange-correlation potential with respect to charge density.

Go to top.

Stress

The stress includes all contributions to the change of energy with the lattice constants; the calculation of stress in Conquest is documented in a paper being prepared for publication, but we give a brief overview here. As Conquest uses orthorhombic cells, only the diagonal stress components (\(\sigma_{\alpha\alpha}\)) are calculated.

In most cases, forces also contribute to the stress; it is easy to show that the stress contribution is given by:

\[\sigma_{\alpha\alpha} = \sum_i F_{i\alpha}R_{i\alpha}\]

where \(R_{i\alpha}\) is the position of the atom. As well as these contributions, there are more subtle terms. Any energies calculated on the grid will contribute to the stress as the integration grid changes with cell size (the stress is simply the energy calculated), and the Hartree potential contributes a term related to the change in the reciprocal lattice vectors (as it is calculated by Fourier transforming the charge density). If the exchange-correlation functional is a GGA functional, then a further term coming from the change of the gradient of the density with the cell size arises. (For non-self-consistent calculations this leads to some complications, as this term technically requires both input and output densities; at present, we approximate this as a mixture of the term calculated with input density and the term calculated with output density; the proportion can be adjusted using the parameter General.MixXCGGAInOut documented in the Advanced and obscure tags section of the manual, though we do not recommend changing it.)

Go to top.

[EFS1]

Steven G. Louie, Sverre Froyen, and Marvin L. Cohen. Nonlinear ionic pseudopotentials in spin-density-functional calculations. Phys. Rev. B, 26:1738–1742, Aug 1982. doi:10.1103/PhysRevB.26.1738.

[EFS2]

T. Miyazaki, D. R. Bowler, R. Choudhury, and M. J. Gillan. Atomic force algorithms in density functional theory electronic-structure techniques based on local orbitals. J. Chem. Phys., 121:6186, 2004. doi:10.1063/1.1787832.

[EFS3]

Antonio S Torralba, David R Bowler, Tsuyoshi Miyazaki, and Michael J Gillan. Non-self-consistent Density-Functional Theory Exchange-Correlation Forces for GGA Functionals. J. Chem. Theory Comput., 5:1499, 2009. doi:10.1021/ct8005425.

Go to top.