Background

Your program ``dft.m'' now solves the Kohn-Sham equations in the potential coming from a simple harmonic oscillator. For most electronic structure calculations, however, the potential comes from the interaction between electrons and ions with a form like

$$V(\vec r) = - \sum_I \frac{Z_I}{\vert\vec r-\vec X_I\vert},$$ (1)

where $Z_I$ and $\vec X_I$ are the charge and location of each nucleus. We now proceed to determine the value Vdual which corresponds to this potential, including all relevant signs and normalizations.

Our approach to finding the potential Vdual from the nuclei for our periodic calculations is based on our standard approach to solving Poisson's equation,

$$\vec \phi_{nuc} = {\bf I}L^{-1} \left(-4\pi {\bf O}{\bf J}\vec \rho_{nuc} \right)$$

and the definition of Vdual from the previous exercise on the Harmonic oscillator,

$$\tilde V={\bf J}^\dagger {\bf O}{\bf J}\left(-\vec \phi_{nuc} \right),$$

where the extra `-' in front of $\phi_{nuc}$ corresponds to the minus sign in Eq. (1). Combining the above two results gives

$$\tilde V={\bf J}^\dagger {\bf O}{\bf L}^{-1} (4 \pi {\bf O}{\bf J}\vec \rho_{nuc})$$ (2)

where we have used the fact that ${\bf J}({\bf I}(\vec x))=\vec x$.

Now, the expansion coefficients $[{\bf J}\vec \rho_{nuc})]_{\vec G} \equiv \hat \rho_{\vec G}$ are normalized so that

$$\rho(\vec r) \equiv \sum_{\vec G} \hat \rho_{\vec G} e^{i\vec G \cdot \vec r}.$$

Therefore,

$$\begin{eqnarray*} \int_{\mbox{cell}} e^{-i\vec G \cdot \vec r} \rho(\vec r)\,dV ... ...' \cdot \vec r}\,dV\\ & = & \det({\bf R})\, \hat \rho_{\vec G}, \end{eqnarray*}$$

where $\det({\bf R})$ is the volume of the cell and the integral extends over that volume. This means

$$\hat \rho_{\vec G} = \frac{1}{\det({\bf R})} \int_{\mbox{cell}} e^{-i\vec G \cdot \vec r} \rho(\vec r)\,dV.$$

Finally, we have one copy of each nucleus in the cell, so that (over the range of the integral) $\rho(\vec r)=\sum_I Z \delta^{(3)}\left(\vec r-\vec X_I\right)$, where $\delta^{(3)}(\vec x)$ is the Dirac-delta function in three-dimensions, and $Z$ is the charge of each nucleus. Thus,

$$\hat \rho_{\vec G} = \frac{Z}{\det({\bf R})} \sum_I e^{-i\vec G \cdot \vec X_I} = \frac{Z}{\det({\bf R})} S_f(\vec G),$$

where $S_f(\vec G)$ is the same structure factor

$$S_f(\vec G) \equiv \sum_{I} e^{-i \vec G \cdot \vec X_I},$$ (3)

from the previous problem on computing Ewald energies. Inserting this into Eq. (2) and using the normalizations associated with the various operators gives

$$\tilde V={\bf J}^\dagger \hat v,$$ (4)

where

$$\hat v_{\vec G} = \left( \frac{-4 \pi Z}{G^2} \right) S_f(\vec G).$$ (5)

Here, the term $-4 \pi Z/G^2$ is just the Fourier-space version of the ionic potential $-Z/\vert\vec r\vert$, the structure factor $S_f(\vec G)$ replicates each of the ions in the appropriate locations, and ${\bf J}$ Fourier transforms back to real space with the appropriate normalization.