Total Energy and Forces: Ver. 1.2

Taisuke Ozaki, RCIS, JAIST

Total energy for the collinear case

The OpenMX is based on density functional theories (DFT), the norm-conserving pseudopotentials, and local pseudo-atomic basis functions. The Kohn-Sham (KS) Bloch functions $\psi_{\mu}$ are expanded in a form of linear combination of pseudo-atomic basis functions (LCPAO) $\phi_{i\alpha}$ centered on site $\tau_i$ by

$$\psi_{\sigma\mu}^{(\bf k)}({\bf r}) = \frac{1}{\sqrt{N}} \sum_{{\... ...ma\mu,i\alpha}^{(\bf k)} \phi_{i\alpha}({\bf r}-{\bf\tau}_{i}-{\bf R}_{\rm n}),$$ (1)

where $c$ and $\phi$ are an expansion coefficient and pseudo-atomic function, ${\bf R}_{\rm n}$ a lattice vector, $i$ a site index, $\sigma~(\uparrow~{\rm or}~\downarrow)$ spin index, $\alpha\equiv (plm)$ an organized orbital index with a multiplicity index $p$, an angular momentum quantum number $l$, and a magnetic quantum number $m$. The charge density operator ${\hat n}_{\sigma}$ for the spin index $\sigma$ is given by

$${\hat n}_{\sigma} = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum... ...vert \psi_{\sigma\mu}^{(\bf k)}\rangle \langle \psi_{\sigma\mu}^{(\bf k)}\vert,$$ (2)

where $\int_{\rm B}$ means the integration over the first Brillouin zone of which volume is $V_{\rm B}$, and $\sum^{\rm occ}$ means the summation over occupied states. The charge density $n_{\sigma}(\bf r)$ with the spin index $\sigma$ is found as

$$n_{\sigma}(\bf r) = \langle {\bf r} \vert {\hat n}_{\sigma} \vert {\bf r} \rangle,$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\mu}^{\rm occ} \la... ...a\mu}^{(\bf k)}\rangle \langle \psi_{\sigma\mu}^{(\bf k)}\vert {\bf r} \rangle,$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{{\rm n}}^{N} \sum_... ...a}({\bf r}-{\bf\tau}_{i}) \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}),$$

$$= \sum_{{\rm n}}^{N} \sum_{i\alpha, j\beta} \rho_{\sigma, i\alpha j... ...a}({\bf r}-{\bf\tau}_{i}) \phi_{j\beta} ({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n})$$ (3)

with a density matrix defined by

$$\rho_{\sigma, i\alpha j\beta}^{({\bf R}_{\rm n})} = \frac{1}{V_{\... ...\cdot {\bf k}} c^{(\bf k)*}_{\sigma\mu,i\alpha} c^{(\bf k)}_{\sigma\mu,j\beta}.$$ (4)

Although it is assumed that the electronic temperature is zero in this notes, OpenMX uses the Fermi-Dirac function with a finite temperature in the practical implementation. Therefore, the force on atom becomes inaccurate for metallic systems or very high temperature. The total charge density $n$ is the sum of $n_{\uparrow}$ and $n_{\downarrow}$ as follows:

$$n({\bf r}) = n_{\uparrow}({\bf r}) + n_{\downarrow}({\bf r}).$$ (5)

Also, it is convenient to define a difference charge density $\delta n({\bf r})$ for later discussion as

$$\delta n({\bf r}) = n({\bf r}) - n^{\rm (a)}({\bf r}),$$

$$= n({\bf r}) - \sum_{i}n_i^{\rm (a)}({\bf r}),$$ (6)

where $n_i^{\rm (a)}({\bf r})$ is an atomic charge density evaluated by a confinement atomic calculations associated with the site $i$. Within the local density approximation (LDA) and generalized gradient approximation (GGA), the total energy of the collinear case is given by the sum of the kinetic energy $E_{\rm kin}$, the electron-core Coulomb energy $E_{\rm ec}$, the electron-electron Coulomb energy $E_{\rm ee}$, the exchange-correlation energy $E_{\rm xc}$, and the core-core Coulomb energy $E_{\rm cc}$ as

$$E_{\rm tot} = E_{\rm kin} + E_{\rm ec} + E_{\rm ee} + E_{\rm xc} + E_{\rm cc}.$$ (7)

The kinetic energy $E_{\rm kin}$ is given by

$$E_{\rm kin} = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...i_{\sigma\mu}^{(\bf k)} \vert \hat{T} \vert \psi_{\sigma\mu}^{(\bf k)} \rangle,$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...ert \hat{T} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$

$$= \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{\sig... ...a j\beta}^{({\bf R}_{\rm n})} h_{i\alpha j\beta,{\rm kin}}^{({\bf R}_{\rm n})}.$$ (8)

The electron-core Coulomb energy $E_{\rm ec}$ is given by two contributions $E_{\rm ec}^{(\rm L)}$ and $E_{\rm ec}^{(\rm NL)}$ related to the local and non-local parts of pseudopotentials:

$$E_{\rm ec} = E_{\rm ec}^{(\rm L)} + E_{\rm ec}^{(\rm NL)},$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...+ V_{{\rm NL},I}({\bf r}-\tau_{I}) \} \vert \psi_{\sigma\mu}^{(\bf k)} \rangle,$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...r}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...r}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$= \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{\sig... ...r}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{\sig... ...}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$ (9)

where $V_{{\rm core},I}$ and $V_{{\rm NL},I}$ are the local and non-local parts of pseudopotential located on a site $I$. Thus, we have

$$E_{\rm ec}^{(\rm L)} = \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{\sig... ...}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$

$$= \int dr^3 n({\bf r}) \sum_{I}V_{{\rm core},I}({\bf r}-\tau_{I}).$$ (10)

$$E_{\rm ec}^{(\rm NL)} = \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{\sig... ...}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle.$$ (11)

The electron-electron Coulomb energy $E_{\rm ee}$ is given by

$$E_{\rm ee} = \frac{1}{2} \int \int dr^3 dr'^3 \frac{n({\bf r})n({\bf r'})}{\vert {\bf r}-{\bf r'}\vert},$$

$$= \frac{1}{2} \int dr^3 n({\bf r}) V_{\rm H}({\bf r}),$$

$$= \frac{1}{2} \int dr^3 n({\bf r}) \{V^{\rm (a)}_{\rm H}({\bf r}) + \delta V_{\rm H}({\bf r})\},$$ (12)

where $V_{\rm H}$ is decomposed into two contributions $V^{\rm (a)}_{\rm H}$ and $\delta V_{\rm H}({\bf r})$ coming from the superposition of atomic charge densities and the difference charge density $\delta n({\bf r})$ defined by

$$V^{\rm (a)}_{\rm H}({\bf r}) = \sum_{I} \int dr'^3 \frac{n_I^{\rm (a)}({\bf r})}{\vert {\bf r}-{\bf r'}\vert},$$

$$= \sum_{I} V^{\rm (a)}_{{\rm H},I}({\bf r}-\tau_{I}),$$ (13)

$$\delta V_{\rm H}({\bf r}) = \int dr'^3 \frac{\delta n({\bf r})}{\vert {\bf r}-{\bf r'}\vert}.$$ (14)

Within the LDA, the exchange-correlation energy $E_{\rm xc}$ is given by

$$E_{\rm xc} = \int dr^3 \{n_{\uparrow}({\bf r})+ n_{\downarrow}({\... ...c}(n_{\uparrow}+\frac{1}{2}n_{\rm pcc}, n_{\downarrow}+\frac{1}{2}n_{\rm pcc}),$$ (15)

where $n_{\rm pcc}$ is a charge density used for a partial core correction (PCC). The core-core Coulomb energy $E_{\rm cc}$ is given as repulsive Coulomb interactions among effective core charge $Z_{I}$ considered in the generation of pseudopotentials by

$$E_{\rm cc} = \frac{1}{2}\sum_{I,J} \frac{Z_{I}Z_{J}} {\vert \tau_{I}-\tau_{J} \vert}.$$ (16)

As discussed in the paper [1], for numerical accuracy and efficiency it is important to reorganize the sum of three terms $E^{(L)}_{\rm ec}$, $E_{\rm ee}$, and $E_{\rm cc}$, as follows:

$$E^{(L)}_{\rm ec} + E_{\rm ee} + E_{\rm cc} = \int dr^3 n({\bf r}) \sum_{I}V_{{\rm core},I}({\bf r}-\tau_{I}) +... ...^{\rm (a)}_{\rm H}({\bf r}) - \int dr^3 n({\bf r}) V^{\rm (a)}_{\rm H}({\bf r})$$

$$+ \frac{1}{2} \int dr^3 n({\bf r}) \{V^{\rm (a)}_{\rm H}({\bf r})... ...f r}) - \frac{1}{2} \int dr^3 n^{\rm (a)}({\bf r}) V^{\rm (a)}_{\rm H}({\bf r})$$

$$+ \frac{1}{2}\sum_{I,J} \frac{Z_{I}Z_{J}} {\vert \tau_{I}-\tau_{J} \vert},$$

$$= \int dr^3 n({\bf r}) \sum_{I} V_{{\rm na}, I}({\bf r}-\tau_{I}) +... ...{\bf r}) - \frac{1}{2} \int dr^3 \delta n({\bf r}) V^{\rm (a)}_{\rm H}({\bf r})$$

$$+ \frac{1}{2}\sum_{I,J} \left[ \frac{Z_{I}Z_{J}} {\vert \tau_{I}-... ... - \int dr^3 n^{\rm (a)}_{I}({\bf r}) V^{\rm (a)}_{{\rm H},J}({\bf r}) \right],$$

$$= \int dr^3 n({\bf r}) \sum_{I} V_{{\rm na}, I}({\bf r}-\tau_{I}) + \frac{1}{2} \int \delta n({\bf r}) \delta V_{\rm H}({\bf r}) d{\bf r}$$

$$+ \frac{1}{2}\sum_{I,J} \left[ \frac{Z_{I}Z_{J}} {\vert \tau_{I}-... ...int n^{\rm (a)}_{I}({\bf r}) V^{\rm (a)}_{{\rm H},J}({\bf r}) d{\bf r} \right],$$ (17)

where

$$V_{{\rm na}, I}({\bf r}-\tau_{I}) = V_{{\rm core},I}({\bf r}-\tau_{I}) + V^{\rm (a)}_{{\rm H},I}({\bf r}-\tau_{I}).$$ (18)

Therefore, we can reorganize these three terms as follows:

$$E^{(L)}_{\rm ec} + E_{\rm ee} + E_{\rm cc} = E_{\rm na} + E_{\rm\delta ee} + E_{\rm ecc},$$ (19)

$$E_{\rm na} = \int dr^3 n({\bf r}) \sum_{I} V_{{\rm na}, I}({\bf r}-\tau_{I}),$$

$$= \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{\sig... ...}-\tau_{I}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$ (20)

$$E_{\rm\delta ee} = \frac{1}{2} \int dr^3 \delta n({\bf r}) \delta V_{\rm H}({\bf r}),$$ (21)

$$E_{\rm scc} = \frac{1}{2}\sum_{I,J} \left[ \frac{Z_{I}Z_{J}} {\vert \tau_{I}-\t... ... - \int dr^3 n^{\rm (a)}_{I}({\bf r}) V^{\rm (a)}_{{\rm H},J}({\bf r}) \right].$$ (22)

Following the reorganization of energy terms, the total energy can be given by

$$E_{\rm tot} = E_{\rm kin} + E_{\rm na} + E_{\rm ec}^{(\rm NL)} + E_{\rm\delta ee} + E_{\rm xc} + E_{\rm scc}.$$ (23)

Kohn-Sham equation

Considering the orthogonality relation among one-particle wave functions, let us introduce a functional $F$ with Lagrange's multipliers $\varepsilon_{\sigma\nu \sigma ' \nu'}$:

$$F = E_{\rm tot} + \sum_{\bf k}\sum_{\sigma} \sum_{\nu\nu'} \varep... ...e \psi_{\sigma\nu}^{(\bf k)} \vert \psi_{\sigma\nu'}^{(\bf k)} \rangle \right).$$ (24)

The variation of the functional $F$ with respect to the LCPAO coefficients $c_{\sigma\mu,i\alpha}^{(\bf k)*}$ yields the following matrix equation:

$$H_{\sigma}^{(\bf k)} c_{\sigma}^{(\bf k)} = S^{(\bf k)} c_{\sigma}^{(\bf k)} \varepsilon_{\sigma \sigma}^{(\bf k)}.$$ (25)

Moreover, noting that the matrix $\varepsilon_{\sigma \sigma}^{(\bf k)}$ for Lagrange's multipliers is Hermitian, and introducing $V$ diagonalizing the matrix, we can transform above equation as:

$$H_{\sigma}^{(\bf k)} c_{\sigma}^{(\bf k)}V = S^{(\bf k)} c_{\sigma}^{(\bf k)} V V^{\dag } \varepsilon_{\sigma \sigma}^{(\bf k)}V.$$ (26)

By renaming $c_{\sigma}^{(\bf k)} V$ by $c_{\sigma}^{(\bf k)}$, and defining the diagonal element of $V^{\dag } \varepsilon_{\sigma \sigma}^{(\bf k)}V$ to be $\varepsilon_{\sigma\mu}^{(\bf k)}$, we have a well known Kohn-Sham equation in a matrix form as a generalized eigenvalue problem:

$$H_{\sigma}^{(\bf k)} c_{\sigma\mu}^{(\bf k)} = \varepsilon_{\sigma\mu}^{(\bf k)} S^{(\bf k)} c_{\sigma\mu}^{(\bf k)},$$ (27)

where the Hamiltonian and overlap matrices are given by

$$H_{\sigma,i\alpha j\beta}^{(\bf k)} = \sum_{\rm n} {\rm e}^{{\rm i}{\bf R}_{\rm n}\cdot {\bf k}} \langl... ...H}_{\sigma} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$

$$= \sum_{\rm n} {\rm e}^{{\rm i}{\bf R}_{\rm n}\cdot {\bf k}} h_{\sigma,i\alpha j\beta}^{({\bf R}_{\rm n})}.$$ (28)

$$S_{i\alpha j\beta}^{(\bf k)} = \sum_{\rm n} {\rm e}^{{\rm i}{\bf R}_{\rm n}\cdot {\bf k}} \langl... ...f\tau}_{i}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$

$$= \sum_{\rm n} {\rm e}^{{\rm i}{\bf R}_{\rm n}\cdot {\bf k}} s_{\sigma,i\alpha j\beta}^{({\bf R}_{\rm n})}.$$ (29)

with a Kohn-Sham Hamiltonian operator

$$\hat{H}_{\sigma} = \hat{T} + V_{\rm eff,\sigma},$$ (30)

and the effective potential

$$V_{\rm eff,\sigma} = \sum_{k} V_{{\rm NL},k} + \sum_{k} V_{{\rm na}, k} + \delta V_{\rm H} + V_{\rm xc,\sigma}.$$ (31)

Projector expansion of $V_{{\rm na}, I}$

The neutral atom potential $V_{{\rm na}, I}$ is spherical and defined within the finite range determined by the cutoff radius $r_{c,I}$ of the confinement potential. Therefore, the potential $V_{{\rm na}, I}$ can be expressed by a projector expansion as follows:

$$\hat{V}_{{\rm na}, I} = \sum_{lm}^{L_{\rm max}} \sum_{\zeta}^{N_{... ...\rangle \frac{1}{c_{l\zeta}} \langle Y_{lm}\bar{R}_{l\zeta}V_{{\rm na},I}\vert,$$ (32)

where a set of radial functions $\{\bar{R}_{l\zeta}\}$ is an orthonormal set defined by a norm $\int r^2 dr RV_{{\rm na},k}R'$ for radial functions $R$ and $R'$, and is calculated by the following Gram-Schmidt orthogonalization [4]:

$$\bar{R}_{l\zeta} = R_{l\zeta} -\sum_{\eta}^{\zeta-1} \bar{R}_{l\eta} \frac{1}{c_{l\eta}} \int r^2 dr \bar{R}_{l\eta}V_{{\rm na},k}R_{l\zeta},$$ (33)

$$c_{l\zeta} = \int r^2 dr \bar{R}_{l\zeta}V_{{\rm na},k}\bar{R}_{l\zeta}.$$ (34)

The radial function $R_{l\zeta}$ used is pseudo wave functions for both the ground and excited states under the same confinement potential as used in the calculation of the atomic electron density $n_i^{(a)}$ [2,3]. The most important feature in the projector expansion is that the deep neutral atom potential is expressed by a separable form, and thereby we only have to evaluate the two-center integrals to construct the matrix elements for the neutral atom potential. As discussed later, the two-center integrals can be accurately evaluated in momentum space. For details of the projector expansion method, see also the reference [1].

Two-center integrals

The overlap integrals, matrix elements for the non-local potentials, for the neutral atom potentials in a separable form discussed in the Sec. 3, and for the kinetic operator consist of two-center integrals. In this section, the evaluation of the two-center integrals is discussed. The pseudo-atomic function $\phi_{i\alpha}({\bf r})$ in Eq. (1) is given by a product of a real spherical harmonic function $\bar{Y}_{lm}$ and a radial wave function $R_{pl}$:

$$\phi_{i\alpha}({\bf r}) = \bar{Y}_{lm}(\hat{\bf r}) R_{pl}(r),$$ (35)

where $\hat{\bf r}$ means Euler angles, $\theta$ and $\phi$, for ${\bf r}$, and $r$ radial coordinate. Although the real function is used for the spherical harmonic function in OpenMX instead of the complex function $Y_{lm}(\hat{\bf r})$, firstly we consider the case with the complex function $Y_{lm}(\hat{\bf r})$ as

$$\phi_{i\alpha}({\bf r}) = Y_{lm}(\hat{\bf r}) R_{pl}(r).$$ (36)

After the evaluation of two-center integrals related to the complex function, they are transformed into two-center integrals for the real function $\bar{Y}_{lm}(\hat{\bf r})$ by matrix operations. Then, the pseudo-atomic function $\phi_{i\alpha}({\bf r})$ given by Eq. (36) can be Fourier-transformed as follows:

$$\tilde{\phi}_{i\alpha}({\bf k}) = \left( \frac{1}{\sqrt{2\pi}} \right )^3 \int dr^3 \phi_{i\alpha}({\bf r}) {\rm e}^{-{\rm i}{\bf k}\cdot {\bf r}}.$$ (37)

The back transformation is defined by

$$\phi_{i\alpha}({\bf r}) = \left( \frac{1}{\sqrt{2\pi}} \right )^3 \int dk^3 \tilde{\phi}_{i\alpha}({\bf k}) {\rm e}^{{\rm i}{\bf k}\cdot {\bf r}}.$$ (38)

Noting the Rayleigh expansion

$${\rm e}^{{\rm i}{\bf k}\cdot {\bf r}} = 4\pi \sum_{L=0}^{\infty} \sum_{M=-L}^{L} {\rm i}^{L} j_{L}(kr) Y_{LM}^{*}(\hat{\bf k})Y_{LM}(\hat{\bf r}),$$ (39)

$${\rm e}^{-{\rm i}{\bf k}\cdot {\bf r}} = 4\pi \sum_{L=0}^{\infty} \sum_{M=-L}^{L} (-{\rm i})^{L} j_{L}(kr) Y_{LM}(\hat{\bf k})Y_{LM}^{*}(\hat{\bf r}),$$ (40)

and putting Eqs. (36) and (40) into Eq. (37), we have

$$\tilde{\phi}_{i\alpha}({\bf k}) = \left( \frac{1}{\sqrt{2\pi}} \right )^3 \int dr^3 Y_{lm}(\hat{\bf... ...} (-{\rm i})^{L} j_{L}(kr) Y_{LM}(\hat{\bf k})Y_{LM}^{*}(\hat{\bf r}) \right\},$$

$$= \left( \frac{1}{\sqrt{2\pi}} \right )^3 4\pi \sum_{L=0}^{\infty} ... ...kr) \int d\theta d\phi \sin(\theta) Y_{lm}(\hat{\bf r})Y_{LM}^{*}(\hat{\bf r}),$$

$$= \left[ \left( \frac{1}{\sqrt{2\pi}} \right )^3 4\pi (-{\rm i})^{l} \int dr r^2 R_{pl}(r) j_{l}(kr) \right] Y_{lm}(\hat{\bf k}),$$

$$= \tilde{R}_{pl}(k)Y_{lm}(\hat{\bf k}),$$ (41)

where we defined

$$\tilde{R}_{pl}(k)= \left[ \left( \frac{1}{\sqrt{2\pi}} \right )^3 4\pi (-{\rm i})^{l} \int dr r^2 R_{pl}(r) j_{L}(kr) \right].$$ (42)

Using Eqs. (38), (40), and (41), the overlap integral is evaluated as follows:

$$\langle \phi_{i\alpha}({\bf r}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}) \rangle = \int dr^3 \phi^*_{i\alpha}({\bf r}) \phi_{j\beta}({\bf r}-{\bf\tau}),$$

$$= \int dr^3 \left( \frac{1}{\sqrt{2\pi}} \right )^3 \int dk^3 \tild... ...(k') Y_{l'm'}(\hat{\bf k}') {\rm e}^{{\rm i}{\bf k}'\cdot ({\bf r}-{\bf\tau})},$$

$$= \left( \frac{1}{2\pi} \right )^3 \int dk^3 \int dk'^3 {\rm e}^{-{... ...l'm'}(\hat{\bf k}') \int dr^3 {\rm e}^{{\rm i}({\bf k}'-{\bf k})\cdot {\bf r}},$$

$$= \int dk^3 {\rm e}^{-{\rm i}{\bf k}\cdot \bf\tau} \tilde{R}^*_{pl}(k) Y^*_{lm}(\hat{\bf k}) \tilde{R}_{p'l'}(k) Y_{l'm'}(\hat{\bf k}),$$

$$= \int dk^3 \left[ 4\pi \sum_{L=0}^{\infty} \sum_{M=-L}^{L} (-{\rm ... ...e{R}^*_{pl}(k) Y^*_{lm}(\hat{\bf k}) \tilde{R}_{p'l'}(k) Y_{l'm'}(\hat{\bf k}),$$

$$= 4\pi \sum_{L=0}^{\infty} \sum_{M=-L}^{L} (-{\rm i})^{L} Y_{LM}^{*... ...LM} \int dk k^2 j_{L}(k\vert\tau\vert) \tilde{R}^*_{pl}(k) \tilde{R}_{p'l'}(k),$$ (43)

where an integral in the third line was evaluated as

$$\int dr^3 {\rm e}^{{\rm i}({\bf k}'-{\bf k})\cdot {\bf r}} = (2\pi)^3\delta({\bf k}'-{\bf k}),$$ (44)

and one may notice Gaunt coefficients defined by

$$C_{l(-m),l'm',LM} = \int d\theta d\phi \sin(\theta) Y^*_{lm}(\hat{\bf k}) Y_{l'm'}(\hat{\bf k}) Y_{LM}(\hat{\bf k}).$$ (45)

The matrix element for the kinetic operator can be easily found in the same way as for the overlap integral as follows:

$$\langle \phi_{i\alpha}({\bf r}) \vert \hat{T} \vert \phi_{j\beta}({\bf r}-{\bf\tau}) \rangle = \int dr^3 \phi^*_{i\alpha}({\bf r}) \left\{ -\frac{1}{2}(\frac{d^... ...+\frac{d^2}{dy^2}+\frac{d^2}{dz^2}) \right \} \phi_{j\beta}({\bf r}-{\bf\tau}),$$

$$= 2\pi \sum_{L=0}^{\infty} \sum_{M=-L}^{L} (-{\rm i})^{L} C_{l(-m),l'm',LM} \int dk k^4 j_{L}(kr) \tilde{R}^*_{pl}(k) \tilde{R}_{p'l'}(k).$$ (46)

Discretization

The two energy terms $E_{\rm\delta ee}$ and $E_{\rm xc}$ are discretized by a regular mesh $\{{\bf r}_{\bf p}\}$ in real space [5], while the integrations in $E_{\rm kin}$, $E_{\rm na}$, $E_{\rm scc}$, and $E_{\rm ec}^{(\rm NL)}$ can be reduced to two center integrals which can be evaluated in momentum space. The regular mesh $\{{\bf r}_{\bf p}\}$ in real space is generated by dividing the unit cell with a same interval which is characterized by the cutoff energies $E_{\rm cut}^{(1)}$, $E_{\rm cut}^{(2)}$, and $E_{\rm cut}^{(3)}$:

$$E_{\rm cut}^{(1)} = \frac{1}{2} {\bf gb}_{1}\cdot {\bf gb}_{1}, \quad E_{\rm cut}^{(2... ... gb}_{2}, \quad E_{\rm cut}^{(3)} = \frac{1}{2} {\bf gb}_{3}\cdot {\bf gb}_{3},$$ (47)

$${\bf ga}_{1} = \frac{{\bf a}_{1}}{N_{1}}, \quad {\bf ga}_{2} = \frac{{\bf a}_{2}}{N_{2}}, \quad {\bf ga}_{3} = \frac{{\bf a}_{3}}{N_{3}},$$ (48)

$${\bf gb}_{1} = 2\pi\frac{{\bf ga}_{2}\times {\bf ga}_{3}} {\Delta V}, \quad {\bf... ...V}, \quad {\bf gb}_{3} = 2\pi\frac{{\bf ga}_{1}\times {\bf ga}_{2}} {\Delta V},$$ (49)

$$\Delta V = {\bf ga}_{1}\cdot \left({\bf ga}_{2}\times {\bf ga}_{3} \right ),$$ (50)

where ${\bf a}_{1}$, ${\bf a}_{2}$, and ${\bf a}_{3}$ are unit cell vectors of the whole system, and $N_1$, $N_2$, and $N_3$ are the number of division for ${\bf a}_{1}$, ${\bf a}_{2}$, and ${\bf a}_{3}$-axes, respectively. $N_1$, $N_2$, and $N_3$ are determined so that the differences among $E_{\rm cut}^{(1)}$, $E_{\rm cut}^{(2)}$, and $E_{\rm cut}^{(3)}$ can be minimized starting from the given cutoff energy. Using the regular mesh $\{{\bf r}_{\bf p}\}$, the Hartree energy $E_{\rm\delta ee}$ associated with the difference charge density $\delta n({\bf r})$ is discretized as

$$E_{\rm\delta ee} = \frac{1}{2} \Delta V \sum_{\bf p} \delta n({\bf r}_{\bf p}) \delta V_{\rm H}({\bf r}_{\bf p}).$$ (51)

The same regular mesh $\{{\bf r}_{\bf p}\}$ is also applied to the solution of Poisson's equation to find $\delta V_{\rm H}$. Then, the charge density is Fourier transformed by

$$\delta \tilde{n}({\bf q}) = \frac{1}{N_1 N_2 N_3} \sum_{n_1}^{N_1-1} \sum_{n_2}^{N_2-1} \sum_... ...ta \tilde{n}({\bf r}_{n_1n_2n_3}) {\rm e}^{-i{\bf q}\cdot {\bf r}_{n_1n_2n_3}},$$

$$= \frac{1}{N_1 N_2 N_3} \sum_{\bf p} \delta n({\bf r}_{\bf p}) {\rm e}^{-i{\bf q}\cdot {\bf r}_{\bf p} }.$$ (52)

From $\{\delta \tilde{n}({\bf q})\}$, we can evaluate $\delta \tilde{V}_{\rm H}({\bf q})$ by

$$\delta \tilde{V}_{\rm H}({\bf q}) = \frac{4\pi}{\vert{\bf q} \vert^2} \delta \tilde{n}({\bf q}).$$ (53)

Thus, we see that $\delta V_{\rm H}({\bf r})$ is also written in a discretized form with the regular mesh as follows:

$$\delta V_{\rm H}({\bf r}) = \sum_{\bf q} \delta \tilde{V}_{\rm H}({\bf q}) {\rm e}^{i{\bf q}\cdot {\bf r}},$$

$$= \frac{4\pi}{N_1 N_2 N_3} \sum_{\bf q} \frac{1}{\vert{\bf q} \vert... ...\delta n({\bf r}_{\bf p}) {\rm e}^{i{\bf q}\cdot ({\bf r} - {\bf r}_{\bf p}) }.$$ (54)

Within the LDA, $E_{\rm xc}$ can be easily discretized as well as $E_{\rm\delta ee}$ by

$$E_{\rm xc} = \Delta V \sum_{\bf p}\{n_{\uparrow} ({\bf r}_{\bf p}) + n_{\downa... ...{\bf p}), n_{\downarrow}({\bf r_p}) +\frac{1}{2} n_{\rm pcc}({\bf r}_{\bf p})).$$ (55)

For the GGA, $E_{\rm xc}$ is discretized with the gradient of charge density evaluated with a finite difference scheme in the same way in the LDA. Since the derivative of the charge density $n({\bf r_p})$ with respect to $c^{({\bf k})*}_{\sigma\mu,i\alpha}$ is given by

$$\frac{\partial n({\bf r_p})} {\partial c^{(\bf k)*}_{\sigma\mu,i\... ...\bf R}_{\rm n}\cdot {\bf k}} \phi_{i\alpha}({\bf r_p})\phi_{j\beta}({\bf r_p}),$$ (56)

the matrix elements for $V_{\rm\delta ee}$ and $V_{\rm xc}$ in the Kohn-Sham equation Eq. (27) are found by differentiating the energies $E_{\rm\delta ee}$ and $E_{\rm xc}$ with respect to $c^{({\bf k})*}_{\sigma\mu,i\alpha}$ as follows:

$$\frac{\partial E_{\rm\delta ee}} {\partial c^{({\bf k})*}_{\sigma\mu,i\alpha}} = \sum_{\bf p} \frac{\partial n({\bf r}_{\bf p})} {\partial c^{({\b... ...ma\mu,i\alpha}} \frac{\partial E_{\rm\delta ee}} {\partial n({\bf r}_{\bf p})},$$

$$= \sum_{j\beta} c_{\sigma\mu,j\beta} \left[ \sum_{{\rm n}}^{N} {\rm... ...}({\bf r_p}) \delta V_{\rm H}({\bf r}_{\bf p}) \phi_{j\beta}({\bf r_p}) \right]$$ (57)

and

$$\frac{\partial E_{\rm xc}} {\partial c^{({\bf k})*}_{\sigma\mu,i\alpha}} = \sum_{\bf p} \frac{\partial n_{\sigma}({\bf r}_{\bf p})} {\partia... ...mu,i\alpha}} \frac{\partial E_{\rm xc}} {\partial n_{\sigma}({\bf r}_{\bf p})},$$

$$= \sum_{j\beta} c_{\sigma\mu,j\beta} \left[ \sum_{{\rm n}}^{N} {\rm... ...}({\bf r_p}) v_{{\rm xc},\sigma} ({\bf r_p} ) \phi_{j\beta}({\bf r_p}) \right],$$ (58)

where the quantities in the parenthesis $[]$ correspond to the matrix elements.

Force on atom

$$\frac{\partial E_{\rm tot}}{\partial \tau_{k}} = \frac{\partial E... ..._{\rm xc}}{\partial \tau_{k}} + \frac{\partial E_{\rm scc}}{\partial \tau_{k}}.$$ (59)

The derivative of the kinetic energy with respect to $\tau_{k}$ is given by

$$\frac{\partial E_{\rm kin}}{\partial \tau_{k}} = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...vert \hat{T} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...vert \hat{T} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle}{\partial \tau_{k}}. \quad$$ (60)

The derivative of the neutral atom potential energy with respect to $\tau_{k}$ is given by

$$\frac{\partial E_{\rm na}}{\partial \tau_{k}} = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...a}}({\bf r}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...a}}({\bf r}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle \right\} } {\partial \tau_{k}}.\quad\quad$$ (61)

The derivative of the non-local potential energy with respect to $\tau_{k}$ is given by

$$\frac{\partial E_{\rm ec}^{(\rm NL)}}{\partial \tau_{k}} = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...L}}({\bf r}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...L}}({\bf r}) \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle$$

$$+ \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...u}_{j}-{\bf R}_{\rm n}) \rangle \right\} } {\partial \tau_{k}}. \quad\quad\quad$$ (62)

The derivative of the Hartree energy energy for the difference charge density $\delta n({\bf r})$ with respect to $\tau_{k}$ is given by

$$\frac{\partial E_{\rm\delta ee}} {\partial \tau_k} = \sum_{\bf p}... ...au_k} \frac{\partial E_{\rm\delta ee}} {\partial n^{\rm (a)}({\bf r}_{\bf p})}.$$ (63)

Considering Eq. (54) and

$$\frac{\partial \delta V_{\rm H}({\bf r}_{\bf p})} {\partial n({\b... ...frac{1}{\vert{\bf q} \vert^2} {\rm e}^{i{\bf q}\cdot ({\bf r_p} - {\bf r_q}) },$$ (64)

we have

$$\frac{\partial E_{\rm\delta ee}} {\partial n({\bf r}_{\bf q})} = \frac{1}{2} \Delta V \{ \delta V_{\rm H}({\bf r}_{\bf q}) + \sum_... ...c{\partial \delta V_{\rm H}({\bf r}_{\bf p})} {\partial n({\bf r}_{\bf q})} \},$$

$$= \frac{1}{2} \Delta V \delta V_{\rm H}({\bf r}_{\bf q}) + \frac{4\... ... p} \delta n({\bf r}_{\bf p}) {\rm e}^{i{\bf q}\cdot ({\bf r_p} - {\bf r_q}) },$$

$$= \Delta V \delta V_{\rm H}({\bf r}_{\bf q}),$$ (65)

and

$$\frac{\partial E_{\rm\delta ee}} {\partial n^{\rm (a)}({\bf r}_{\bf q})} = -\frac{1}{2} \Delta V \delta V_{\rm H}({\bf r}_{\bf q}) + \frac{1... ...ial \delta V_{\rm H}({\bf r}_{\bf p})} {\partial n^{\rm (a)}({\bf r}_{\bf q})},$$

$$= -\frac{1}{2} \Delta V \delta V_{\rm H}({\bf r}_{\bf q}) - \frac{1... ... p} \delta n({\bf r}_{\bf p}) {\rm e}^{i{\bf q}\cdot ({\bf r_p} - {\bf r_q}) },$$

$$= -\Delta V \delta V_{\rm H}({\bf r}_{\bf q}).$$ (66)

Moreover, the derivative of $n({\bf r_p})$ with respect to $\tau_k$ is given by

$$\frac{\partial n({\bf r_p})} {\partial \tau_k} = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \sum_{\mu}... ...ta}} {\partial \tau_k} \phi_{i\alpha}({\bf r_p})\phi_{j\beta}({\bf r_p})\right)$$

$$+ \frac{2}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma}\sum_{\al... ...\partial \phi_{k\alpha}({\bf r_p})} {\partial \tau_k} \phi_{j\beta}({\bf r_p}).$$ (67)

The derivative of $n^{\rm (a)}({\bf r_p})$ with respect to $\tau_k$ is simply given by

$$\frac{\partial n^{\rm (a)}({\bf r_p})} {\partial \tau_k} = \frac{\partial n^{\rm (a)}_{k}({\bf r_p})} {\partial \tau_k}.$$ (68)

The derivative of $E_{\rm xc}$ with respect to the atomic coordinate $\tau_k$ is easily evaluated by

$$\frac{\partial E_{\rm xc}} {\partial \tau_k} = \sum_{\sigma} \sum_{\bf p} \frac{\partial n_{\sigma}({\bf r}_{\bf... ...ial \tau_k} \frac{\partial E_{\rm xc}} {\partial n_{\rm pcc}({\bf r}_{\bf p})},$$

$$= \Delta V \sum_{\sigma} \sum_{\bf p} \frac{\partial n_{\sigma}({\b... ..._k} V_{{\rm xc},\sigma} ( n({\bf r_p}),n_{\rm pcc}({\bf r_p}) ).\quad\quad\quad$$ (69)

The derivative of the screened core-core Coulomb energy $E_{\rm scc}$ with respect to $\tau_k$ is given by

$$\frac{\partial E_{\rm scc}}{\partial \tau_{k}} = \frac{1}{2}\sum_... ...{\bf r}) V^{\rm (a)}_{{\rm H},J}({\bf r}) \right\} } {\partial \tau_k} \right].$$ (70)

Since the second term is tabulated in a numerical table as a function of distance due to the spherical symmetry of integrands, the derivative can be evaluated analytically by employing an interpolation scheme. The derivatives given by Eqs. (60), (61), (62), (63), and (69) contain the derivative of LCPAO coefficient $c$. The derivative of $c$ can be transformed to the derivative of the overlap matrix with respect to $\tau_k$ as shown below. By summing up all the terms including the derivatives of $c$ in Eqs. (60), (61), (62), (63), and (69), we have

$$A = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \left\... ...ac{\partial c^{(\bf k)}_{\sigma} }{\partial \tau_k} \right) \right\},\quad\quad$$ (71)

where $\Theta$ is a diagonal matrix consisting of Heaviside step functions. Noting that Eq. (27) can be written by $H_{\sigma}^{(\bf k)} c_{\sigma}^{(\bf k)}=S^{(\bf k)}c_{\sigma}^{(\bf k)}\varepsilon_{\sigma}^{(\bf k)}$ and $c_{\sigma}^{(\bf k)\dag }H_{\sigma}^{(\bf k)}=\varepsilon_{\sigma}^{(\bf k)}c_{\sigma}^{(\bf k)\dag } S^{(\bf k)}$ in a matrix form, the product of two diagonal matrices is commutable, and ${\rm Tr}(XY)={\rm Tr}(YX)$ for any square matrices, Eq. (71) can be written as

$$A = \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \left\{ {\... ...bf k)} \frac{\partial c^{(\bf k)}_{\sigma} }{\partial \tau_k} \right) \right\},$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \left\{ {\... ...)}_{\sigma} }{\partial \tau_k} \varepsilon_{\sigma}^{(\bf k)} \right) \right\},$$

$$= \frac{1}{V_{\rm B}} \int_{{\rm B}}dk^{3} \sum_{\sigma} \left\{ {\... ...ma} }{\partial \tau_k} \right) \varepsilon_{\sigma}^{(\bf k)} \right) \right\}.$$ (72)

Moreover, taking account of the derivative of the orthogonality relation $c^{(\bf k)\dag }_{\sigma}S^{(\bf k)}c_{\sigma}^{(\bf k)}=I$ with respect to $\tau_k$, we have the following relation:

$$\frac{\partial c^{(\bf k)\dag }_{\sigma} }{\partial \tau_k} S^{(\... ...g }_{\sigma} \frac{\partial S^{(\bf k)}}{\partial \tau_k} c_{\sigma}^{(\bf k)}.$$ (73)

Putting Eq. (73) into Eq. (72), we have

$$A = -\sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} E_{\s... ... n})} \frac{\partial S_{i\alpha j\beta}^{({\bf R}_{\rm n})}} {\partial \tau_k},$$

where the energy density matrix $E_{\sigma, i\alpha j\beta}^{({\bf R}_{\rm n})}$ is given by

$$E_{\sigma, i\alpha j\beta}^{({\bf R}_{\rm n})} = \frac{1}{V_{\rm ... ...\mu}^{(\bf k)} c^{(\bf k)*}_{\sigma\mu,i\alpha} c^{(\bf k)}_{\sigma\mu,j\beta}.$$ (74)

The terms including the derivative of matrix elements in Eqs. (60), (61), and (62) can be easily evaluated by

$$B = \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{... ...h_{i\alpha j\beta,V_{\rm NL}}^{({\bf R}_{\rm n})} \right\} } {\partial \tau_k},$$ (75)

where

$$h_{i\alpha j\beta,{\rm kin}}^{({\bf R}_{\rm n})} = \langle \phi_{i\alpha}({\bf r}-{\bf\tau}_{i}) \vert \hat{T} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$ (76)

$$h_{i\alpha j\beta,V_{\rm na}}^{({\bf R}_{\rm n})} = \langle \phi_{i\alpha}({\bf r}-{\bf\tau}_{i}) \vert V_{\rm na} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$ (77)

$$h_{i\alpha j\beta,V_{\rm NL}}^{({\bf R}_{\rm n})} = \langle \phi_{i\alpha}({\bf r}-{\bf\tau}_{i}) \vert V_{\rm NL} \vert \phi_{j\beta}({\bf r}-{\bf\tau}_{j}-{\bf R}_{\rm n}) \rangle,$$ (78)

The derivatives of these elements are evaluated analytically from the analytic derivatives of Eqs. (43) and (46). The remaining contributions in first terms of Eqs. (63) and (69) are given by

$$C = \sum_{\sigma} \sum_{{\rm n}}^{N} \sum_{i\alpha,j\beta} \rho_{... ...igma}({\bf r}_{\bf q},n_{\rm pcc}({\bf r_p}) \right\} \phi_{j\beta}({\bf r_p}).$$ (79)

The second terms in Eqs. (63) and (69) becomes

$$D = -\Delta V \sum_{\bf p} \delta V_{\rm H}({\bf r_p}) \frac{\par... ... {\partial \tau_k} V_{{\rm xc},\sigma} ( n({\bf r_p}),n_{\rm pcc}({\bf r_p}) ).$$ (80)

Thus, we see that the derivative of the total energy with respect to the atomic coordinate $\tau_k$ is analytically evaluated for any grid fineness as

$$\frac{\partial E_{\rm tot}}{\partial \tau_{k}} = A + B + C + D + \frac{\partial E_{\rm scc}}{\partial \tau_{k}}.$$ (81)

Bibliography

1

T. Ozaki and H. Kino, Phys. Rev. B 72, 045121 (2005).

2

T. Ozaki, Phys. Rev. B 67, 155108 (2003).

3

T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004).

4

P. E. Blochl, Phys. Rev. B 41, R5414 (1990).

5

J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002).