theory_energy_force_stress.rst

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 = 2Tr[KH] + ΔE_Ha + Δ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:

Ĥ = T̂ + V̂_L + 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 = ∫dr^′n(r^′)/ ∣ r − r^′∣, and the exchange-correlation potential. The alternative form, often known as the DFT energy, is:

E_DFT = 2Tr[K(T+V_L+V_NL)] + 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 = ∫drϵ_XC[n]n(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

ΔE_XC = ∫dr(ϵ_XC[n]−V_XC[n])n(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 <theory_energy_force_stress>.

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 δn(r) = n(r) − ∑_in_i(r) for atomic densities n_i(r). We write:

E_DFT, NA = 2Tr[K(T+V_NA+V_NL)] + E_δ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(r) is the Hartree potential from the atomic density n_i(r). We define the neutral atom potential for an atom as V_NA, i(r) = V_L, i(r) + V_Ha, i(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(r) = ∑_iV_NA, i(r). If we write the pseudo-atomic density as n_PAD(r) = ∑_in_i(r) then we can also write V_NA(r) = V_L(r) + V_Ha, PAD(r).

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

E_HF = 2Tr[KH] + ΔE_Ha + ΔE_XC + E_SII

with the Hamiltonian defined as:

Ĥ = T̂ + V̂_NA + V̂_NL + V_δHa + V_XC

where V_δHa(r) = ∫dr^′δn(r^′)/ ∣ r − r^′∣. Accordingly, the double counting Hartree correction term has to change:

ΔE_Ha =  − E_δHa − ∫drδn(r)∑_iV_Ha, i(r).

Go to top <theory_energy_force_stress>.

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(r) and write:

E_NSC = 2Tr[KH] + ΔE_Ha[n_in] + ΔE_XC[n_in] + 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 = ∑_ijϕ_iK_ijϕ_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 <th_efs_nap>, δn(r) = 0 by definition, which also means that E_δHa = 0 and ΔE_Ha = 0.

Go to top <theory_energy_force_stress>.

Partial core corrections

Also known as non-linear core corrections, partial core corrections (PCC) efs-Louie:1982aa 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(r) + n_c(r) where the valence charge is input or output charge density defined above: V_XC[n_v+n_c]. The exchange-correlation energy becomes:

E_XC = ∫dr(n_v(r)+n_c(r))V_XC[n_v+n_c].

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] = ∫drn_v(r)V_XC[n_v + n_c]. We write:

$$\begin{aligned} \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{aligned}$$

There is an extra factor of ∫drn_c(r)ϵ_XC[n_v + n_c] over and above the usual term.

Go to top <theory_energy_force_stress>.

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 <theory_energy_force_stress>.

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 efs-Miyazaki2004,efs-Torralba:2009nr though the changes needed to account for PCC <th_efs_pcc>, particularly in the non-self-consistent <th_efs_nsc> 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 <th_efs_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:

F_i^PCC =  − ∫dr∇_in_i^c(r)V_XC[n_v + n_c]

If the non-self-consistent <th_efs_nsc> 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:

F_i^NSC =  − ∫drV_δHa(r)∇_in_i^v(r) − ∫drδn(r)V_XC^′[n_vⁱⁿ+n_c](∇_in_i^v(r)+∇_in_i^c(r))

where V_XC^′ is the derivative of the exchange-correlation potential with respect to charge density.

Go to top <theory_energy_force_stress>.

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 (σ_αα) are calculated.

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

σ_αα = ∑_iF_iαR_iα

where R_iα 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 <advanced_general_tags> section of the manual, though we do not recommend changing it.)

Go to top <theory_energy_force_stress>.

references.bib

Go to top <theory_energy_force_stress>.