Total Energy and Forces: Ver. 1.3

1 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) ${\varphi }_{i\alpha }$ centered on site ${\tau }_i$ by

${\psi }_{\sigma \mu }^{(𝐤)}(𝐫)={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{\displaystyle \sum _{i\alpha }}{c}_{\sigma \mu ,i\alpha }^{(𝐤)}{\varphi }_{i\alpha }(𝐫-{\tau }_i-{𝐑}_{\mathrm{n}}),$ (1)

where $c$ and $\varphi$ are an expansion coefficient and pseudo-atomic function, ${𝐑}_{\mathrm{n}}$ a lattice vector, $i$ a site index, $\sigma$ $(\uparrow$ $\mathrm{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 ${\widehat{n}}_{\sigma }$ for the spin index $\sigma$ is given by

${\widehat{n}}_{\sigma }={\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}|{\psi }_{\sigma \mu }^{(𝐤)}⟩⟨{\psi }_{\sigma \mu }^{(𝐤)}|,$ (2)

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

$n_{\sigma }(𝐫)$ $=$ $⟨𝐫|{\widehat{n}}_{\sigma }|𝐫⟩,$ (3)

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}⟨𝐫|{\psi }_{\sigma \mu }^{(𝐤)}⟩⟨{\psi }_{\sigma \mu }^{(𝐤)}|𝐫⟩,$

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\varphi }_{i\alpha }(𝐫-{\tau }_i){\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}}),$

$=$ ${\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{\varphi }_{i\alpha }(𝐫-{\tau }_i){\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})$

with a density matrix defined by

${\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}={\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\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(𝐫)=n_{\uparrow }(𝐫)+n_{\downarrow }(𝐫).$ (5)

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

$\delta n(𝐫)$ $=$ $n(𝐫)-n^{(\mathrm{a})}(𝐫),$ (6)

$=$ $n(𝐫)-{\displaystyle \sum _i}{n}_i^{(\mathrm{a})}(𝐫),$

where $n_i^{(\mathrm{a})}(𝐫)$ 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_{\mathrm{kin}}$, the electron-core Coulomb energy $E_{\mathrm{ec}}$, the electron-electron Coulomb energy $E_{\mathrm{ee}}$, the exchange-correlation energy $E_{\mathrm{xc}}$, and the core-core Coulomb energy $E_{\mathrm{cc}}$ as

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

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

$E_{\mathrm{kin}}$ $=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}⟨{\psi }_{\sigma \mu }^{(𝐤)}|\widehat{T}|{\psi }_{\sigma \mu }^{(𝐤)}⟩,$ (8)

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|\widehat{T}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$

$=$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{h}_{i\alpha j\beta ,\mathrm{kin}}^{({𝐑}_{\mathrm{n}})}.$

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

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

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}⟨{\psi }_{\sigma \mu }^{(𝐤)}|{\displaystyle \sum _I}\{V_{\mathrm{core},I}(𝐫-{\tau }_I)+V_{\mathrm{NL},I}(𝐫-{\tau }_I)\}|{\psi }_{\sigma \mu }^{(𝐤)}⟩,$

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\displaystyle \sum _I}{V}_{\mathrm{core},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\displaystyle \sum _I}{V}_{\mathrm{NL},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$

$=$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\displaystyle \sum _I}{V}_{\mathrm{core},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$

$+$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\displaystyle \sum _I}{V}_{\mathrm{NL},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$

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

$E_{\mathrm{ec}}^{(\mathrm{L})}$ $=$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\displaystyle \sum _I}{V}_{\mathrm{core},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$ (10)

$=$ ${\displaystyle \int 𝑑r^{3}n(𝐫)\sum _I{V}_{\mathrm{core},I}(𝐫-{\tau }_I)}.$

$E_{\mathrm{ec}}^{(\mathrm{NL})}$ $=$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\displaystyle \sum _I}{V}_{\mathrm{NL},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩.$ (11)

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

$E_{\mathrm{ee}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \int \int 𝑑r^3𝑑r^{\prime 3}\frac{n(𝐫)n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}},$ (12)

$=$ ${\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^{3}n(𝐫)V_{\mathrm{H}}(𝐫)},$

$=$ ${\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^{3}n(𝐫)\{V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)+\delta V_{\mathrm{H}}(𝐫)\}},$

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

$V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)$ $=$ ${\displaystyle \sum _I}{\displaystyle \int 𝑑r^{\prime 3}\frac{n_I^{(\mathrm{a})}(𝐫)}{|𝐫-{𝐫}^{\prime }|}},$ (13)

$=$ ${\displaystyle \sum _I}{V}_{\mathrm{H},I}^{(\mathrm{a})}(𝐫-{\tau }_I),$

$\delta V_{\mathrm{H}}(𝐫)$ $=$ ${\displaystyle \int 𝑑r^{\prime 3}\frac{\delta n(𝐫)}{|𝐫-{𝐫}^{\prime }|}}.$ (14)

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

$E_{\mathrm{xc}}={\displaystyle \int 𝑑r^3\{n_{\uparrow }(𝐫)+n_{\downarrow }(𝐫)+n_{\mathrm{pcc}}(𝐫)\}{ϵ}_{\mathrm{xc}}(n_{\uparrow }+\frac{1}{2}{n}_{\mathrm{pcc}},n_{\downarrow }+\frac{1}{2}{n}_{\mathrm{pcc}})},$ (15)

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

$E_{\mathrm{cc}}={\displaystyle \frac{1}{2}}{\displaystyle \sum _{I,J}}{\displaystyle \frac{Z_I{Z}_J}{|{\tau }_I-{\tau }_J|}}.$ (16)

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

$E_{\mathrm{ec}}^{(L)}+E_{\mathrm{ee}}+E_{\mathrm{cc}}$ $=$ ${\displaystyle \int 𝑑r^{3}n(𝐫)\sum _I{V}_{\mathrm{core},I}(𝐫-{\tau }_I)}+{\displaystyle \int 𝑑r^{3}n(𝐫)V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)}-{\displaystyle \int 𝑑r^{3}n(𝐫)V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)}$ (17)

$+{\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^{3}n(𝐫)\{V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)+\delta V_{\mathrm{H}}(𝐫)\}}+{\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^3{n}^{(\mathrm{a})}(𝐫)V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)}-{\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^3{n}^{(\mathrm{a})}(𝐫)V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)}$

$+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{I,J}}{\displaystyle \frac{Z_I{Z}_J}{|{\tau }_I-{\tau }_J|}},$

$=$ ${\displaystyle \int 𝑑r^{3}n(𝐫)\sum _I{V}_{\mathrm{na},I}(𝐫-{\tau }_I)}+{\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^{3}n(𝐫)\delta V_{\mathrm{H}}(𝐫)}-{\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^3\delta n(𝐫)V_{\mathrm{H}}^{(\mathrm{a})}(𝐫)}$

$+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{I,J}}\left[{\displaystyle \frac{Z_I{Z}_J}{|{\tau }_I-{\tau }_J|}}-{\displaystyle \int 𝑑r^3{n}_I^{(\mathrm{a})}(𝐫)V_{\mathrm{H},J}^{(\mathrm{a})}(𝐫)}\right],$

$=$ ${\displaystyle \int 𝑑r^{3}n(𝐫)\sum _I{V}_{\mathrm{na},I}(𝐫-{\tau }_I)}+{\displaystyle \frac{1}{2}}{\displaystyle \int \delta n(𝐫)\delta V_{\mathrm{H}}(𝐫)𝑑𝐫}$

$+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{I,J}}\left[{\displaystyle \frac{Z_I{Z}_J}{|{\tau }_I-{\tau }_J|}}-{\displaystyle \int n_I^{(\mathrm{a})}(𝐫)V_{\mathrm{H},J}^{(\mathrm{a})}(𝐫)𝑑𝐫}\right],$

where

$V_{\mathrm{na},I}(𝐫-{\tau }_I)$ $=$ $V_{\mathrm{core},I}(𝐫-{\tau }_I)+V_{\mathrm{H},I}^{(\mathrm{a})}(𝐫-{\tau }_I).$ (18)

Therefore, we can reorganize these three terms as follows:

$E_{\mathrm{ec}}^{(L)}+E_{\mathrm{ee}}+E_{\mathrm{cc}}=E_{\mathrm{na}}+E_{\delta \mathrm{ee}}+E_{\mathrm{scc}},$ (19)

$E_{\mathrm{na}}$ $=$ ${\displaystyle \int 𝑑r^{3}n(𝐫)\sum _I{V}_{\mathrm{na},I}(𝐫-{\tau }_I)},$ (20)

$=$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{\displaystyle \sum _I}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{na},I}(𝐫-{\tau }_I)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$

$E_{\delta \mathrm{ee}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑r^3\delta n(𝐫)\delta V_{\mathrm{H}}(𝐫)},$ (21)

$E_{\mathrm{scc}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \sum _{I,J}}\left[{\displaystyle \frac{Z_I{Z}_J}{|{\tau }_I-{\tau }_J|}}-{\displaystyle \int 𝑑r^3{n}_I^{(\mathrm{a})}(𝐫)V_{\mathrm{H},J}^{(\mathrm{a})}(𝐫)}\right].$ (22)

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

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

2 Kohn-Sham equation

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

$F=E_{\mathrm{tot}}+{\displaystyle \sum _{𝐤}}{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\nu {\nu }^{\prime }}}{\epsilon }_{\sigma \nu \sigma {\nu }^{\prime }}^{(𝐤)}\left({\delta }_{\sigma \nu \sigma {\nu }^{\prime }}^{(𝐤)}-⟨{\psi }_{\sigma \nu }^{(𝐤)}|{\psi }_{\sigma {\nu }^{\prime }}^{(𝐤)}⟩\right).$ (24)

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

$H_{\sigma }^{(𝐤)}{c}_{\sigma }^{(𝐤)}=S^{(𝐤)}{c}_{\sigma }^{(𝐤)}{\epsilon }_{\sigma \sigma }^{(𝐤)}.$ (25)

Moreover, noting that the matrix ${\epsilon }_{\sigma \sigma }^{(𝐤)}$ for Lagrange’s multipliers is Hermitian, and introducing $V$ diagonalizing the matrix, we can transform above equation as:

$H_{\sigma }^{(𝐤)}{c}_{\sigma }^{(𝐤)}V=S^{(𝐤)}{c}_{\sigma }^{(𝐤)}V{V}^{†}{\epsilon }_{\sigma \sigma }^{(𝐤)}V.$ (26)

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

$H_{\sigma }^{(𝐤)}{c}_{\sigma \mu }^{(𝐤)}={\epsilon }_{\sigma \mu }^{(𝐤)}{S}^{(𝐤)}{c}_{\sigma \mu }^{(𝐤)},$ (27)

where the Hamiltonian and overlap matrices are given by

$H_{\sigma ,i\alpha j\beta }^{(𝐤)}$ $=$ ${\displaystyle \sum _{\mathrm{n}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\widehat{H}}_{\sigma }|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$ (28)

$=$ ${\displaystyle \sum _{\mathrm{n}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{h}_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}.$

$S_{i\alpha j\beta }^{(𝐤)}$ $=$ ${\displaystyle \sum _{\mathrm{n}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$ (29)

$=$ ${\displaystyle \sum _{\mathrm{n}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{s}_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}.$

with a Kohn-Sham Hamiltonian operator

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

and the effective potential

$V_{\mathrm{eff},\sigma }={\displaystyle \sum _k}{V}_{\mathrm{NL},k}+{\displaystyle \sum _k}{V}_{\mathrm{na},k}+\delta V_{\mathrm{H}}+V_{\mathrm{xc},\sigma }.$ (31)

3 Projector expansion of $V_{\mathrm{na},I}$

The neutral atom potential $V_{\mathrm{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_{\mathrm{na},I}$ can be expressed by a projector expansion as follows:

${\widehat{V}}_{\mathrm{na},I}={\displaystyle \sum _{lm}^{L_{\max }}}{\displaystyle \sum _{\zeta }^{N_{\mathrm{rad}}}}|V_{\mathrm{na},I}{\overline{R}}_{l\zeta }{Y}_{lm}⟩{\displaystyle \frac{1}{c_{l\zeta }}}⟨Y_{lm}{\overline{R}}_{l\zeta }{V}_{\mathrm{na},I}|,$ (32)

where a set of radial functions $\{{\overline{R}}_{l\zeta }\}$ is an orthonormal set defined by a norm $\int r^2𝑑rR{V}_{\mathrm{na},k}{R}^{\prime }$ for radial functions $R$ and $R^{\prime }$, and is calculated by the following Gram-Schmidt orthogonalization [4]:

${\overline{R}}_{l\zeta }$ $=$ $R_{l\zeta }-{\displaystyle \sum _{\eta }^{\zeta -1}}{\overline{R}}_{l\eta }{\displaystyle \frac{1}{c_{l\eta }}}{\displaystyle \int r^2𝑑r{\overline{R}}_{l\eta }{V}_{\mathrm{na},k}{R}_{l\zeta }},$ (33)

$c_{l\zeta }$ $=$ ${\displaystyle \int r^2𝑑r{\overline{R}}_{l\zeta }{V}_{\mathrm{na},k}{\overline{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].

4 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 ${\varphi }_{i\alpha }(𝐫)$ in Eq. (1) is given by a product of a real spherical harmonic function ${\overline{Y}}_{lm}$ and a radial wave function $R_{pl}$:

${\varphi }_{i\alpha }(𝐫)={\overline{Y}}_{lm}(\widehat{𝐫})R_{pl}(r),$ (35)

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

${\varphi }_{i\alpha }(𝐫)=Y_{lm}(\widehat{𝐫})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 ${\overline{Y}}_{lm}(\widehat{𝐫})$ by matrix operations. Then, the pseudo-atomic function ${\varphi }_{i\alpha }(𝐫)$ given by Eq. (36) can be Fourier-transformed as follows:

${\tilde{\varphi }}_{i\alpha }(𝐤)={\left({\displaystyle \frac{1}{\sqrt{2\pi }}}\right)}^3{\displaystyle \int 𝑑r^3{\varphi }_{i\alpha }(𝐫){\mathrm{e}}^{-\mathrm{i}𝐤\cdot 𝐫}}.$ (37)

The back transformation is defined by

${\varphi }_{i\alpha }(𝐫)={\left({\displaystyle \frac{1}{\sqrt{2\pi }}}\right)}^3{\displaystyle \int 𝑑k^3{\tilde{\varphi }}_{i\alpha }(𝐤){\mathrm{e}}^{\mathrm{i}𝐤\cdot 𝐫}}.$ (38)

Noting the Rayleigh expansion

${\mathrm{e}}^{\mathrm{i}𝐤\cdot 𝐫}$ $=$ $4\pi {\displaystyle \sum _{L=0}^{\mathrm{\infty }}}{\displaystyle \sum _{M=-L}^L}{\mathrm{i}}^L{j}_L(kr)Y_{LM}^{*}(\widehat{𝐤})Y_{LM}(\widehat{𝐫}),$ (39)

${\mathrm{e}}^{-\mathrm{i}𝐤\cdot 𝐫}$ $=$ $4\pi {\displaystyle \sum _{L=0}^{\mathrm{\infty }}}{\displaystyle \sum _{M=-L}^L}{(-\mathrm{i})}^L{j}_L(kr)Y_{LM}(\widehat{𝐤})Y_{LM}^{*}(\widehat{𝐫}),$ (40)

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

${\tilde{\varphi }}_{i\alpha }(𝐤)$ $=$ ${\left({\displaystyle \frac{1}{\sqrt{2\pi }}}\right)}^3{\displaystyle \int 𝑑r^3{Y}_{lm}(\widehat{𝐫})R_{pl}(r)\left\{4\pi \sum _{L=0}^{\mathrm{\infty }}\sum _{M=-L}^L{(-\mathrm{i})}^L{j}_L(kr)Y_{LM}(\widehat{𝐤})Y_{LM}^{*}(\widehat{𝐫})\right\}},$ (41)

$=$ ${\left({\displaystyle \frac{1}{\sqrt{2\pi }}}\right)}^{3}4\pi {\displaystyle \sum _{L=0}^{\mathrm{\infty }}}{\displaystyle \sum _{M=-L}^L}{(-\mathrm{i})}^L{Y}_{LM}(\widehat{𝐤}){\displaystyle \int 𝑑r{r}^2{R}_{pl}(r)j_L(kr)\int 𝑑\theta 𝑑\varphi \sin (\theta )Y_{lm}(\widehat{𝐫})Y_{LM}^{*}(\widehat{𝐫})},$

$=$ $\left[{\left({\displaystyle \frac{1}{\sqrt{2\pi }}}\right)}^{3}4\pi {(-\mathrm{i})}^l{\displaystyle \int 𝑑r{r}^2{R}_{pl}(r)j_l(kr)}\right]Y_{lm}(\widehat{𝐤}),$

$=$ ${\tilde{R}}_{pl}(k)Y_{lm}(\widehat{𝐤}),$

where we defined

${\tilde{R}}_{pl}(k)=\left[{\left({\displaystyle \frac{1}{\sqrt{2\pi }}}\right)}^{3}4\pi {(-\mathrm{i})}^l{\displaystyle \int 𝑑r{r}^2{R}_{pl}(r)j_L(kr)}\right].$ (42)

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

$⟨{\varphi }_{i\alpha }(𝐫)|{\varphi }_{j\beta }(𝐫-\tau )⟩$ $=$ ${\displaystyle \int 𝑑r^3{\varphi }_{i\alpha }^{*}(𝐫){\varphi }_{j\beta }(𝐫-\tau )},$ (43)

$=$ ${\displaystyle \int 𝑑r^3{\left(\frac{1}{\sqrt{2\pi }}\right)}^3\int 𝑑k^3{\tilde{R}}_{pl}^{*}(k)Y_{lm}^{*}(\widehat{𝐤}){\mathrm{e}}^{-\mathrm{i}𝐤\cdot 𝐫}{\left(\frac{1}{\sqrt{2\pi }}\right)}^3\int 𝑑k^{\prime 3}{\tilde{R}}_{p^{\prime }{l}^{\prime }}(k^{\prime })Y_{l^{\prime }{m}^{\prime }}({\widehat{𝐤}}^{\prime }){\mathrm{e}}^{\mathrm{i}{𝐤}^{\prime }\cdot (𝐫-\tau )}},$

$=$ ${\left({\displaystyle \frac{1}{2\pi }}\right)}^3{\displaystyle \int 𝑑k^3\int 𝑑k^{\prime 3}{\mathrm{e}}^{-\mathrm{i}{𝐤}^{\prime }\cdot \tau }{\tilde{R}}_{pl}^{*}(k)Y_{lm}^{*}(\widehat{𝐤}){\tilde{R}}_{p^{\prime }{l}^{\prime }}(k^{\prime })Y_{l^{\prime }{m}^{\prime }}({\widehat{𝐤}}^{\prime })\int 𝑑r^3{\mathrm{e}}^{\mathrm{i}({𝐤}^{\prime }-𝐤)\cdot 𝐫}},$

$=$ ${\displaystyle \int 𝑑k^3{\mathrm{e}}^{-\mathrm{i}𝐤\cdot \tau }{\tilde{R}}_{pl}^{*}(k)Y_{lm}^{*}(\widehat{𝐤}){\tilde{R}}_{p^{\prime }{l}^{\prime }}(k)Y_{l^{\prime }{m}^{\prime }}(\widehat{𝐤})},$

$=$ ${\displaystyle \int 𝑑k^3\left[4\pi \sum _{L=0}^{\mathrm{\infty }}\sum _{M=-L}^L{(-\mathrm{i})}^L{j}_L(k\tau )Y_{LM}(\widehat{𝐤})Y_{LM}^{*}(\widehat{\tau })\right]{\tilde{R}}_{pl}^{*}(k)Y_{lm}^{*}(\widehat{𝐤}){\tilde{R}}_{p^{\prime }{l}^{\prime }}(k)Y_{l^{\prime }{m}^{\prime }}(\widehat{𝐤})},$

$=$ $4\pi {\displaystyle \sum _{L=0}^{\mathrm{\infty }}}{\displaystyle \sum _{M=-L}^L}{(-\mathrm{i})}^L{Y}_{LM}^{*}(\widehat{\tau })C_{l(-m),l^{\prime }{m}^{\prime },LM}{\displaystyle \int 𝑑k{k}^2{j}_L(k|\tau |){\tilde{R}}_{pl}^{*}(k){\tilde{R}}_{p^{\prime }{l}^{\prime }}(k)},$

where an integral in the third line was evaluated as

${\displaystyle \int 𝑑r^3{\mathrm{e}}^{\mathrm{i}({𝐤}^{\prime }-𝐤)\cdot 𝐫}}$ $=$ ${(2\pi )}^3\delta ({𝐤}^{\prime }-𝐤),$ (44)

and one may notice Gaunt coefficients defined by

$C_{l(-m),l^{\prime }{m}^{\prime },LM}={\displaystyle \int 𝑑\theta 𝑑\varphi \sin (\theta )Y_{lm}^{*}(\widehat{𝐤})Y_{l^{\prime }{m}^{\prime }}(\widehat{𝐤})Y_{LM}(\widehat{𝐤})}.$ (45)

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

$⟨{\varphi }_{i\alpha }(𝐫)|\widehat{T}|{\varphi }_{j\beta }(𝐫-\tau )⟩$ $=$ ${\displaystyle \int 𝑑r^3{\varphi }_{i\alpha }^{*}(𝐫)\left\{-\frac{1}{2}(\frac{d^2}{d{x}^2}+\frac{d^2}{d{y}^2}+\frac{d^2}{d{z}^2})\right\}{\varphi }_{j\beta }(𝐫-\tau )},$ (46)

$=$ $2\pi {\displaystyle \sum _{L=0}^{\mathrm{\infty }}}{\displaystyle \sum _{M=-L}^L}{(-\mathrm{i})}^L{C}_{l(-m),l^{\prime }{m}^{\prime },LM}{\displaystyle \int 𝑑k{k}^4{j}_L(kr){\tilde{R}}_{pl}^{*}(k){\tilde{R}}_{p^{\prime }{l}^{\prime }}(k)}.$

5 Discretization

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

$E_{\mathrm{cut}}^{(1)}$ $=$ ${\displaystyle \frac{1}{2}}{\mathrm{𝐠𝐛}}_1\cdot {\mathrm{𝐠𝐛}}_1,E_{\mathrm{cut}}^{(2)}={\displaystyle \frac{1}{2}}{\mathrm{𝐠𝐛}}_2\cdot {\mathrm{𝐠𝐛}}_2,E_{\mathrm{cut}}^{(3)}={\displaystyle \frac{1}{2}}{\mathrm{𝐠𝐛}}_3\cdot {\mathrm{𝐠𝐛}}_3,$ (47)

${\mathrm{𝐠𝐚}}_1$ $=$ ${\displaystyle \frac{{𝐚}_1}{N_1}},{\mathrm{𝐠𝐚}}_2={\displaystyle \frac{{𝐚}_2}{N_2}},{\mathrm{𝐠𝐚}}_3={\displaystyle \frac{{𝐚}_3}{N_3}},$ (48)

${\mathrm{𝐠𝐛}}_1$ $=$ $2\pi {\displaystyle \frac{{\mathrm{𝐠𝐚}}_2\times {\mathrm{𝐠𝐚}}_3}{\mathrm{\Delta }V}},{\mathrm{𝐠𝐛}}_2=2\pi {\displaystyle \frac{{\mathrm{𝐠𝐚}}_3\times {\mathrm{𝐠𝐚}}_1}{\mathrm{\Delta }V}},{\mathrm{𝐠𝐛}}_3=2\pi {\displaystyle \frac{{\mathrm{𝐠𝐚}}_1\times {\mathrm{𝐠𝐚}}_2}{\mathrm{\Delta }V}},$ (49)

$\mathrm{\Delta }V$ $=$ ${\mathrm{𝐠𝐚}}_1\cdot \left({\mathrm{𝐠𝐚}}_2\times {\mathrm{𝐠𝐚}}_3\right),$ (50)

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

$E_{\delta \mathrm{ee}}={\displaystyle \frac{1}{2}}\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩})\delta V_{\mathrm{H}}({𝐫}_{𝐩}).$ (51)

The same regular mesh $\{{𝐫}_{𝐩}\}$ is also applied to the solution of Poisson’s equation to find $\delta V_{\mathrm{H}}$. Then, the charge density is Fourier transformed by

$\delta \tilde{n}(𝐪)$ $=$ ${\displaystyle \frac{1}{N_1{N}_2{N}_3}}{\displaystyle \sum _{n_1}^{N_1-1}}{\displaystyle \sum _{n_2}^{N_2-1}}{\displaystyle \sum _{n_3}^{N_3-1}}\delta n({𝐫}_{n_1{n}_2{n}_3}){\mathrm{e}}^{-i𝐪\cdot {𝐫}_{n_1{n}_2{n}_3}},$ (52)

$=$ ${\displaystyle \frac{1}{N_1{N}_2{N}_3}}{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\mathrm{e}}^{-i𝐪\cdot {𝐫}_{𝐩}}.$

From $\{\delta \tilde{n}(𝐪)\}$, we can evaluate $\delta {\tilde{V}}_{\mathrm{H}}(𝐪)$ by

$\delta {\tilde{V}}_{\mathrm{H}}(𝐪)={\displaystyle \frac{4\pi }{{|𝐪|}^2}}\delta \tilde{n}(𝐪).$ (53)

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

$\delta V_{\mathrm{H}}(𝐫)$ $=$ ${\displaystyle \sum _{𝐪}}\delta {\tilde{V}}_{\mathrm{H}}(𝐪){\mathrm{e}}^{i𝐪\cdot 𝐫},$ (54)

$=$ ${\displaystyle \frac{4\pi }{N_1{N}_2{N}_3}}{\displaystyle \sum _{𝐪}}{\displaystyle \frac{1}{{|𝐪|}^2}}{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\mathrm{e}}^{i𝐪\cdot (𝐫-{𝐫}_{𝐩})}.$

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

$E_{\mathrm{xc}}$ $=$ $\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}\{n_{\uparrow }({𝐫}_{𝐩})+n_{\downarrow }({𝐫}_{𝐩})+n_{\mathrm{pcc}}({𝐫}_{𝐩})\}{ϵ}_{\mathrm{xc}}(n_{\uparrow }({𝐫}_{𝐩})+{\displaystyle \frac{1}{2}}{n}_{\mathrm{pcc}}({𝐫}_{𝐩}),n_{\downarrow }({𝐫}_{𝐩})+{\displaystyle \frac{1}{2}}{n}_{\mathrm{pcc}}({𝐫}_{𝐩})).$ (55)

For the GGA, $E_{\mathrm{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({𝐫}_{𝐩})$ with respect to $c_{\sigma \mu ,i\alpha }^{(𝐤)*}$ is given by

${\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}}={\displaystyle \sum _{j\beta }}{c}_{\sigma \mu ,j\beta }{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{\varphi }_{i\alpha }({𝐫}_{𝐩}){\varphi }_{j\beta }({𝐫}_{𝐩}),$ (56)

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

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}}$ $=$ ${\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}}{\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n({𝐫}_{𝐩})}},$ (57)

$=$ ${\displaystyle \sum _{j\beta }}{c}_{\sigma \mu ,j\beta }\left[{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}{\varphi }_{i\alpha }({𝐫}_{𝐩})\delta V_{\mathrm{H}}({𝐫}_{𝐩}){\varphi }_{j\beta }({𝐫}_{𝐩})\right]$

and

${\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}}$ $=$ ${\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n_{\sigma }({𝐫}_{𝐩})}{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}}{\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial n_{\sigma }({𝐫}_{𝐩})}},$ (58)

$=$ ${\displaystyle \sum _{j\beta }}{c}_{\sigma \mu ,j\beta }\left[{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}{\varphi }_{i\alpha }({𝐫}_{𝐩})v_{\mathrm{xc},\sigma }({𝐫}_{𝐩}){\varphi }_{j\beta }({𝐫}_{𝐩})\right],$

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

6 Force on atom

${\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial {\tau }_k}}={\displaystyle \frac{\partial E_{\mathrm{kin}}}{\partial {\tau }_k}}+{\displaystyle \frac{\partial E_{\mathrm{na}}}{\partial {\tau }_k}}+{\displaystyle \frac{\partial E_{\mathrm{ec}}^{(\mathrm{NL})}}{\partial {\tau }_k}}+{\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial {\tau }_k}}+{\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial {\tau }_k}}+{\displaystyle \frac{\partial E_{\mathrm{scc}}}{\partial {\tau }_k}}.$ (59)

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

${\displaystyle \frac{\partial E_{\mathrm{kin}}}{\partial {\tau }_k}}$ $=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{\tau }_{\mathrm{n}}\cdot 𝐤}{\displaystyle \frac{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}{\partial {\tau }_k}}{c}_{\sigma \mu ,j\beta }^{(𝐤)}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|\widehat{T}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$ (60)

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{\displaystyle \frac{\partial c_{\sigma \mu ,j\beta }^{(𝐤)}}{\partial {\tau }_k}}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|\widehat{T}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\displaystyle \frac{\partial ⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|\widehat{T}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩}{\partial {\tau }_k}}.$

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

${\displaystyle \frac{\partial E_{\mathrm{na}}}{\partial {\tau }_k}}$ $=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{\displaystyle \frac{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}{\partial {\tau }_k}}{c}_{\sigma \mu ,j\beta }^{(𝐤)}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{na}}(𝐫)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$ (61)

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{\displaystyle \frac{\partial c_{\sigma \mu ,j\beta }^{(𝐤)}}{\partial {\tau }_k}}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{na}}(𝐫)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\displaystyle \frac{\partial \left\{⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{na}}(𝐫)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩\right\}}{\partial {\tau }_k}}.$

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

${\displaystyle \frac{\partial E_{\mathrm{ec}}^{(\mathrm{NL})}}{\partial {\tau }_k}}$ $=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{\displaystyle \frac{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}{\partial {\tau }_k}}{c}_{\sigma \mu ,j\beta }^{(𝐤)}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{NL}}(𝐫)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$ (62)

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{\displaystyle \frac{\partial c_{\sigma \mu ,j\beta }^{(𝐤)}}{\partial {\tau }_k}}⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{NL}}(𝐫)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩$

$+$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\displaystyle \frac{\partial \left\{⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{NL}}(𝐫)|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩\right\}}{\partial {\tau }_k}}.$

The derivative of the Hartree energy for the difference charge density $\delta n(𝐫)$ with respect to ${\tau }_k$ is given by

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial {\tau }_k}}={\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial {\tau }_k}}{\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n({𝐫}_{𝐩})}}+{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n^{(\mathrm{a})}({𝐫}_{𝐩})}{\partial {\tau }_k}}{\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n^{(\mathrm{a})}({𝐫}_{𝐩})}}.$ (63)

Considering Eq. (54) and

${\displaystyle \frac{\partial \delta V_{\mathrm{H}}({𝐫}_{𝐩})}{\partial n({𝐫}_{𝐪})}}={\displaystyle \frac{4\pi }{N_1{N}_2{N}_3}}{\displaystyle \sum _{𝐪}}{\displaystyle \frac{1}{{|𝐪|}^2}}{\mathrm{e}}^{i𝐪\cdot ({𝐫}_{𝐩}-{𝐫}_{𝐪})},$ (64)

we have

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n({𝐫}_{𝐪})}}$ $=$ ${\displaystyle \frac{1}{2}}\mathrm{\Delta }V\{\delta V_{\mathrm{H}}({𝐫}_{𝐪})+{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\displaystyle \frac{\partial \delta V_{\mathrm{H}}({𝐫}_{𝐩})}{\partial n({𝐫}_{𝐪})}}\},$ (65)

$=$ ${\displaystyle \frac{1}{2}}\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐪})+{\displaystyle \frac{4\pi \mathrm{\Delta }V}{2N_1{N}_2{N}_3}}{\displaystyle \sum _{𝐪}}{\displaystyle \frac{1}{{|𝐪|}^2}}{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\mathrm{e}}^{i𝐪\cdot ({𝐫}_{𝐩}-{𝐫}_{𝐪})},$

$=$ $\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐪}),$

and

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n^{(\mathrm{a})}({𝐫}_{𝐪})}}$ $=$ $-{\displaystyle \frac{1}{2}}\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐪})+{\displaystyle \frac{1}{2}}\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\displaystyle \frac{\partial \delta V_{\mathrm{H}}({𝐫}_{𝐩})}{\partial n^{(\mathrm{a})}({𝐫}_{𝐪})}},$ (66)

$=$ $-{\displaystyle \frac{1}{2}}\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐪})-{\displaystyle \frac{1}{2}}\mathrm{\Delta }V{\displaystyle \frac{4\pi }{N_1{N}_2{N}_3}}{\displaystyle \sum _{𝐪}}{\displaystyle \frac{1}{{|𝐪|}^2}}{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\mathrm{e}}^{i𝐪\cdot ({𝐫}_{𝐩}-{𝐫}_{𝐪})},$

$=$ $-\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐪}).$

Moreover, the derivative of $n({𝐫}_{𝐩})$ with respect to ${\tau }_k$ is given by

${\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial {\tau }_k}}$ $=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\displaystyle \sum _{i\alpha ,j\beta }}{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}\left({\displaystyle \frac{\partial c_{\sigma \mu ,i\alpha }^{(𝐤)*}}{\partial {\tau }_k}}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\varphi }_{i\alpha }({𝐫}_{𝐩}){\varphi }_{j\beta }({𝐫}_{𝐩})+c_{\sigma \mu ,i\alpha }^{(𝐤)*}{\displaystyle \frac{\partial c_{\sigma \mu ,j\beta }^{(𝐤)}}{\partial {\tau }_k}}{\varphi }_{i\alpha }({𝐫}_{𝐩}){\varphi }_{j\beta }({𝐫}_{𝐩})\right)$ (67)

$+{\displaystyle \frac{2}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\alpha ,j\beta }}{\displaystyle \sum _{\mu }}{\displaystyle \sum _{\mathrm{n}}^N}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{c}_{\sigma \mu ,k\alpha }^{(𝐤)*}{c}_{\sigma \mu ,j\beta }^{(𝐤)}{\displaystyle \frac{\partial {\varphi }_{k\alpha }({𝐫}_{𝐩})}{\partial {\tau }_k}}{\varphi }_{j\beta }({𝐫}_{𝐩}).$

The derivative of $n^{(\mathrm{a})}({𝐫}_{𝐩})$ with respect to ${\tau }_k$ is simply given by

${\displaystyle \frac{\partial n^{(\mathrm{a})}({𝐫}_{𝐩})}{\partial {\tau }_k}}={\displaystyle \frac{\partial n_k^{(\mathrm{a})}({𝐫}_{𝐩})}{\partial {\tau }_k}}.$ (68)

The derivative of $E_{\mathrm{xc}}$ with respect to the atomic coordinate ${\tau }_k$ is easily evaluated by

${\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial {\tau }_k}}$ $=$ ${\displaystyle \sum _{\sigma }}{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n_{\sigma }({𝐫}_{𝐩})}{\partial {\tau }_k}}{\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial n_{\sigma }({𝐫}_{𝐩})}}+{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n_{\mathrm{pcc}}({𝐫}_{𝐩})}{\partial {\tau }_k}}{\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial n_{\mathrm{pcc}}({𝐫}_{𝐩})}},$ (69)

$=$ $\mathrm{\Delta }V{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n_{\sigma }({𝐫}_{𝐩})}{\partial {\tau }_k}}{V}_{\mathrm{xc},\sigma }(n({𝐫}_{𝐩},n_{\mathrm{pcc}}({𝐫}_{𝐩}))+{\displaystyle \frac{\mathrm{\Delta }V}{2}}{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n_{\mathrm{pcc}}({𝐫}_{𝐩})}{\partial {\tau }_k}}{V}_{\mathrm{xc},\sigma }(n({𝐫}_{𝐩}),n_{\mathrm{pcc}}({𝐫}_{𝐩})).$

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

${\displaystyle \frac{\partial E_{\mathrm{scc}}}{\partial {\tau }_k}}={\displaystyle \frac{1}{2}}{\displaystyle \sum _{I,J}}\left[{\displaystyle \frac{\partial \left\{\frac{Z_I{Z}_J}{|{\tau }_I-{\tau }_J|}\right\}}{\partial {\tau }_k}}-{\displaystyle \frac{\partial \left\{\int 𝑑r^3{n}_I^{(\mathrm{a})}(𝐫)V_{\mathrm{H},J}^{(\mathrm{a})}(𝐫)\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={\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}\left\{\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)}){\displaystyle \frac{\partial c_{\sigma }^{(𝐤)†}}{\partial {\tau }_k}}{H}_{\sigma }^{(𝐤)}{c}_{\sigma }^{(𝐤)}\right)+\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)})c_{\sigma }^{(𝐤)†}{H}_{\sigma }^{(𝐤)}{\displaystyle \frac{\partial c_{\sigma }^{(𝐤)}}{\partial {\tau }_k}}\right)\right\},$ (71)

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

$A$ $=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}\left\{\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)}){\displaystyle \frac{\partial c_{\sigma }^{(𝐤)†}}{\partial {\tau }_k}}{S}^{(𝐤)}{c}_{\sigma }^{(𝐤)}{\epsilon }_{\sigma }^{(𝐤)}\right)+\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)}){\epsilon }_{\sigma }^{(𝐤)}{c}_{\sigma }^{(𝐤)†}{S}^{(𝐤)}{\displaystyle \frac{\partial c_{\sigma }^{(𝐤)}}{\partial {\tau }_k}}\right)\right\},$ (72)

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}\left\{\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)}){\displaystyle \frac{\partial c_{\sigma }^{(𝐤)†}}{\partial {\tau }_k}}{S}^{(𝐤)}{c}_{\sigma }^{(𝐤)}{\epsilon }_{\sigma }^{(𝐤)}\right)+\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)})c_{\sigma }^{(𝐤)†}{S}^{(𝐤)}{\displaystyle \frac{\partial c_{\sigma }^{(𝐤)}}{\partial {\tau }_k}}{\epsilon }_{\sigma }^{(𝐤)}\right)\right\},$

$=$ ${\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\sigma }}\left\{\mathrm{Tr}\left(\mathrm{\Theta }(E_{\mathrm{F}}-{\epsilon }_{\sigma }^{(𝐤)})\left({\displaystyle \frac{\partial c_{\sigma }^{(𝐤)†}}{\partial {\tau }_k}}{S}^{(𝐤)}{c}_{\sigma }^{(𝐤)}+c_{\sigma }^{(𝐤)†}{S}^{(𝐤)}{\displaystyle \frac{\partial c_{\sigma }^{(𝐤)}}{\partial {\tau }_k}}\right){\epsilon }_{\sigma }^{(𝐤)}\right)\right\}.$

Moreover, taking account of the derivative of the orthogonality relation $c_{\sigma }^{(𝐤)†}{S}^{(𝐤)}{c}_{\sigma }^{(𝐤)}=I$ with respect to ${\tau }_k$, we have the following relation:

${\displaystyle \frac{\partial c_{\sigma }^{(𝐤)†}}{\partial {\tau }_k}}{S}^{(𝐤)}{c}_{\sigma }^{(𝐤)}+c_{\sigma }^{(𝐤)†}{S}^{(𝐤)}{\displaystyle \frac{\partial c_{\sigma }^{(𝐤)}}{\partial {\tau }_k}}=-c_{\sigma }^{(𝐤)†}{\displaystyle \frac{\partial S^{(𝐤)}}{\partial {\tau }_k}}{c}_{\sigma }^{(𝐤)}.$ (73)

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

$A=-{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{E}_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{\displaystyle \frac{\partial S_{i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}}{\partial {\tau }_k}},$

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

$E_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}={\displaystyle \frac{1}{V_{\mathrm{B}}}}{\displaystyle {\int }_{\mathrm{B}}}𝑑k^3{\displaystyle \sum _{\mu }^{\mathrm{occ}}}{\mathrm{e}}^{\mathrm{i}{𝐑}_{\mathrm{n}}\cdot 𝐤}{\epsilon }_{\sigma \mu }^{(𝐤)}{c}_{\sigma \mu ,i\alpha }^{(𝐤)*}{c}_{\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={\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}{\displaystyle \frac{\partial \left\{h_{i\alpha j\beta ,\mathrm{kin}}^{({𝐑}_{\mathrm{n}})}+h_{i\alpha j\beta ,V_{\mathrm{na}}}^{({𝐑}_{\mathrm{n}})}+h_{i\alpha j\beta ,V_{\mathrm{NL}}}^{({𝐑}_{\mathrm{n}})}\right\}}{\partial {\tau }_k}},$ (75)

where

$h_{i\alpha j\beta ,\mathrm{kin}}^{({𝐑}_{\mathrm{n}})}$ $=$ $⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|\widehat{T}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$ (76)

$h_{i\alpha j\beta ,V_{\mathrm{na}}}^{({𝐑}_{\mathrm{n}})}$ $=$ $⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{na}}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$ (77)

$h_{i\alpha j\beta ,V_{\mathrm{NL}}}^{({𝐑}_{\mathrm{n}})}$ $=$ $⟨{\varphi }_{i\alpha }(𝐫-{\tau }_i)|V_{\mathrm{NL}}|{\varphi }_{j\beta }(𝐫-{\tau }_j-{𝐑}_{\mathrm{n}})⟩,$ (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={\displaystyle \sum _{\sigma }}{\displaystyle \sum _{\mathrm{n}}^N}{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{\sigma ,i\alpha j\beta }^{({𝐑}_{\mathrm{n}})}2\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial {\varphi }_{k\alpha }({𝐫}_{𝐩})}{\partial {\tau }_k}}\{\delta V_{\mathrm{H}}({𝐫}_{𝐪})+V_{\mathrm{xc},\sigma }({𝐫}_{𝐪},n_{\mathrm{pcc}}({𝐫}_{𝐩})\}{\varphi }_{j\beta }({𝐫}_{𝐩}).$ (79)

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

$D=-\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}\delta V_{\mathrm{H}}({𝐫}_{𝐩}){\displaystyle \frac{\partial n_k^{(\mathrm{a})}({𝐫}_{𝐩})}{\partial {\tau }_k}}+{\displaystyle \frac{\mathrm{\Delta }V}{2}}{\displaystyle \sum _{\sigma }}{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n_{\mathrm{pcc},\mathrm{k}}({𝐫}_{𝐩})}{\partial {\tau }_k}}{V}_{\mathrm{xc},\sigma }(n({𝐫}_{𝐩}),n_{\mathrm{pcc}}({𝐫}_{𝐩})).$ (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

${\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial {\tau }_k}}=A+B+C+D+{\displaystyle \frac{\partial E_{\mathrm{scc}}}{\partial {\tau }_k}}.$ (81)