Total Energy Minimization

Total Energy Minimization#

In this document, we start from the objective function of solid state DFT, starting from the abstract form and derive it until it is implementable. Solid state DFT solves the following energy minimization problem (we leave why the objective looks like this in the preliminaries of DFT and crystalgraphy)

\[\begin{split}\begin{align} &\min_{f,\psi_{ik}} E_\text{total}[\{\psi_{ik}\}, f] \\ s.t.\quad & \psi_{ik}(r)=e^{ikr}u_{ik}(r); \\ & \langle u_{ik}|u_{jk}\rangle=\delta_{ij}; \\ & \sum_{k=1}^K\sum_{i=1}^I f_{ik} = N; \; f_{ik}\in [0, 1] \end{align}\end{split}\]

The total energy functional is a functional of a set of wave functions \(\{\psi_{ik}\}\) and the occupation vectors \(f\). \(k\in K\) and \(i\in[1..I]\), we call \(I\) the number of bands and \(K\) the set of k points, and \(N\) is the number of electrons in the system. \(f\) specifies how the \(N\) electrons are distributed over the k points and bands. With the wave function and occupation vectors, we can calculate the electron density function as \(\rho(r)=\sum_i^I\sum_{k\in K} f_{ik}\psi^2_{ik}(r)\). As for all KS DFT calculations, the energy functional contains four terms, kinetic energy, external energy, hartree energy, and exchange-correlation energy.

\[E_\text{total}[\{\psi_{ik}\},f] = E_\text{kin}[\{\psi_{ik}\},f] + E_\text{ext}[\rho] + E_\text{har}[\rho] + E_\text{xc}[\rho]\]

And notice that except that the kinetic energy term is a direct functional of the wave functions and occupation, all other energies are direct functionals of the electron density (which in term depends on the wave function and occupation). We write below the definition of each energy term except the exchange-correlation functional. As the crystal system is periodic over the entire \(R^3\) space, we only calculate a fraction of the energy within a single unit cell \(\Omega\).

\[E_\text{kin}[\{\psi_{ik}\},f]=-\frac{1}{2}\sum_i\sum_k f_{ik}\int_{\Omega} \psi^*_{ik}(r)\nabla^2_r\psi_{ik}(r) dr\]
\[E_\text{ext}[\rho] = -\sum_a \int_{\Omega} \rho(r) \frac{Z_a}{r-R_a}dr\]
\[E_\text{har}[\rho] = \frac{1}{2}\int_\Omega dr \int dr'\rho(r)\frac{1}{r-r'}\rho(r')\]

To this point, we have introduced DFT as an optimization problem with constraints in the function space.

  • The objective functional: \(E_\text{total}[\{\psi_{ik}\},f]\), which can be expanded into the above terms.

  • The parameter: \(f\) and \(\psi_{ik}\).

  • The constraints:

    • \(\psi_{ik}(r)=e^{ikr}u_{ik}(r)\), where \(u_{ik}\) is periodic over the unit cell, and \(\langle u_{ik}|u_{jk}\rangle=\delta_{ij}\).

    • \(f_{ik}\in [0,1]\) and \(\sum_k \sum_i f_{ik}=N\).

To make the problem computable, we just need to parameterize \(u_{ik}\) and \(f\) in a way that satisfy the constraints, plugging them back into the objective function and then perform the optimization in the parameter space. This is what we will do in the rest of this document.

Parameterizing \(u_{ik}(r)\) and \(f\)#

In planewave calculations, \(u_{ik}\) is parameterized as a linear combination over fourier components of different frequencies

\[u_{ik}(r)=\frac{1}{\sqrt{\Omega_\text{cell}}}\sum_{G} c_{ikG} e^{iGr}\]

where we limit the \(G\) in \(e^{iGr}\) to be frequency components that is periodic over the unit cell, making \(u_{ik}(r)\) satisfy the periodic constraint. At the same time, plugging the planewave back to \(\langle u_{ik}|u_{jk}\rangle=\delta_{ij}\), we translate the orthogonality constraint into \(\sum_G c^*_{ikG}c_{jkG}=I_{ij}\). A orthogonal matrix can be easily generated via reparameterization, for example, with the QR decomposition

key = jax.random.PRNGKey(0)
w = jax.random.normal(key, (num_K, num_G, I))
c, _ = jnp.linalg.qr(w)

The other parameter \(f\) needs to satisfy \(f_{ik}\in [0,1]\) and \(\sum_k \sum_i f_{ik}=N\). Similarly, we can use reparameterize it as

key = jax.random.PRNGKey(0)
v = jax.random.normal(key, (I*num_K, N*num_K))
Q, _ = jnp.linalg.qr(v)  # Q has shape (I*num_K, N*num_K) and Q.T @ Q = I
f = jnp.diag(Q @ Q.T).reshape(I, num_K)  # f has shape (I*num_K)

It is easily verifiable that \(\sum_{ik}f_{ik}=\Tr(QQ^\top)=\Tr(Q^\top Q)=N\) and \(f_{ik}=\|Q_{i*|K|+k}\|^2\in[0,1]\)

Casting into the parameter space#

Now that \(u_{ik}\) becomes a function parameterized by \(c\), we can substitute it back to the energy terms to cast the energy into a function of the finite dimensional parameters \(c\) and \(f\).

The Kinetic Energy#

Firstly, we apply the kinetic operator on the parameterized wave function

\[\nabla^2_r\psi_{ik}(r) = \nabla^2\left[\frac{1}{\sqrt{\Omega_\text{cell}}}e^{ikr}\sum_{G} c_{ikG}e^{iGr}\right] = \nabla^2\left[\frac{1}{\sqrt{\Omega_\text{cell}}}\sum_{G} c_{ikG}e^{i(k+G)r}\right] = -\|k+G\|^2\psi_{ik}(r)\]

The kinetic energy is then reduced to the following using the property that \(\int_{\Omega} \psi^*_{ik}(r)\psi_{jk}(r) dr=\delta_{ij}\).

\[\begin{split}\begin{align} E_\text{kin}[\{\psi_{ik}\},f]=&\frac{1}{2}\sum_i\sum_k f_{ik}\int_{\Omega} \psi^*_{ik}(r)\nabla^2_r\psi_{ik}(r) dr \\ =& \frac{1}{2}\sum_i\sum_{k}\sum_G f_{ik} c_{ikG}^2\|k+G\|^2 \int_{\Omega} \psi^*_{ik}(r)\psi_{ik}(r) dr \\ =& \frac{1}{2}\sum_i\sum_{k}\sum_G f_{ik} c_{ikG}^2\|k+G\|^2 \end{align}\end{split}\]

Reciprocal representation of the Coulombic potential#

The Coulombic potential generated from a charge \(\rho\) is

\[\begin{split}\begin{equation} \begin{split} V(\vb{r}) =& \rho \star \frac{1}{r} = \int_{\Omega + \vb{R} } \dd{\vb{r}'} \frac{1}{\norm{\vb{r} - \vb{r}'}} \rho(\vb{r}') \\ \end{split} \end{equation}\end{split}\]

where \(\vb{R}\) is a Bravais lattice vector. Its reciprocal representation is given by (see Ewald Summation for derivation)

\[\begin{equation} \tilde{V} (\vb{G}) = \frac{4\pi \tilde{\rho}(\vb{G})}{\norm{\vb{G}}^{2}}. \end{equation}\]

The External Energy#

The atomic point charge within the unit cell is

\[\begin{equation} \rho ^{\text{atom}}(\vb{r})=-\sum_{\ell} Z_{\ell }\delta (\vb{r}-\vb*{\tau }_{\ell }) \end{equation}\]

Therefore its reciprocal representation is given by

\[\begin{equation} \tilde{V}_{\text{ext}}(\vb{G}) = -\sum_{\ell} \frac{4\pi Z_{\ell } e^{-\text{i} \vb{G} ^{\top} \vb*{\tau }_{\ell}}}{\norm{\vb{G}}^{2}}. \end{equation}\]

Now by Parseval’s theorem we have

\[\begin{split}\begin{align} E_\text{ext}[\rho] &= \int_{\Omega} \rho(\vb{r}) V_{\text{ext}}(\vb{r}) \dd{\vb{r}} \\ &= \sum_{\vb{G}\neq \vb{0}} \tilde{V} _{\text{ext}}(\vb{G})^{*} \tilde{\rho} (\vb{G}) \\ &= - 4\pi \sum_{\vb{G} \neq \vb{0}} \tilde{\rho} (\vb{G}) \sum_\ell e^{ -\text{i}\vb{G}^\top \vb*{\tau}_\ell} \dfrac{Z_\ell}{ \Vert \vb{G} \Vert^2} \end{align}\end{split}\]

where the \(\vb{G}=0\) term is removed due to neutral charge requirement (TODO: add doc on this). The reciprocal density \(\tilde{\rho} (\vb{G})\) can be calculated effciently using FFT (see jrystal.pw.density_grid_reciprocal()).

The Hartree Energy#

Again using Parseval’s theorem we have

\[\begin{equation} E_\text{har}[\rho] = \frac{1}{2} \int_\Omega \rho(\vb{r}) \left( \rho \star \frac{1}{r} \right) \dd{\vb{r}} = \frac{1}{2} \sum_{\vb{G}\neq 0} \tilde{\rho}(\vb{G}) \frac{4\pi \tilde{\rho}(\vb{G})}{\norm{\vb{G}}^{2}} \end{equation}\]
Contents