Ewald energy calculator

Molecules and solids, by definition, have more than one atom. Thus, the nuclear-nuclear interaction $U_{NN}$, which is zero for an isolated atom, will soon become for us an important term in the energy. When using periodic boundary conditions, we are effectively repeating each atom in our cell infinitely many times in a periodic arrangement throughout space. Summing the Coulomb energy,

$$U_{NN} \equiv \frac{1}{2} \sum_{I \ne J} \frac{Z_I Z_J}{\vert\vert \vec R_I - \vec R_J\vert\vert},$$

directly then becomes a very tricky task.

A very powerful technique for computing this term is the ``Ewald summation'' technique. This technique is based on the fact that, outside a spherically symmetric charge distribution, the distribution acts just like a point charge concentrated at the center. Thus, we can compute $U_{NN}$ with our Poisson solver by simply replacing each of the nuclear charges with a Gaussian containing the same net charge. The only subtlety is that the Gaussians will inevitably overlap somewhat and not be entirely ``outside'' of each other. The full Ewald algorithm handles this correction explicitly. For us, by choosing the Gaussians to narrow enough that the overlaps are insignificant, we can make a workable (but not optimal) implementation of the Ewald technique. Because your ``poisson.m'' script already gives the energy of continuous, periodic charge distributions, you are a few changes away from having a quick-and-dirty Ewald summation program.

From the above discussion, we see that the Ewald technique hinges on construction of a density of the form

$$n(\vec x) = \sum_I g(\vec x - \vec R_I),$$ (6)

where $g(\vec x)$ is a Gaussian normalized to the charge of each nucleus $Z$ (assumded here to be the same for all nuclei), and the $\vec R_I$ are the locations of the nuclei. Notice that the expansion coefficients $\hat{n_\alpha}$ are easy to compute once we have the expansion coefficients $\hat g_\alpha$ for $g(x)$,

$$g(\vec x) = \sum_\alpha \hat g_\alpha e^{i \vec G_\alpha \cdot \vec x}.$$ (7)

Inserting (9) into (8) gives,

$$\begin{eqnarray*} n(\vec x) & = & \sum_I \sum_\alpha \hat g_\alpha e^{i \vec G_\... ...ha \cdot \vec R_I} \right) e^{i \vec G_\alpha \cdot \vec x}, \\ \end{eqnarray*}$$

where we have used the product rule for addition of exponents and rearranged the terms in the sums. Comparing the result to the expansion for $n(\vec x)$,

$$n(\vec x) = \sum_\alpha \hat n_\alpha e^{i \vec G_\alpha \cdot \vec x},$$

we see immediately that the expansion coefficients which we need for $n(\vec x)$ are just

$$\hat n_\alpha = \hat g_\alpha S_{\vec G_\alpha},$$

where the so-called ``structure factor'' is defined as

$$S_{\vec G_\alpha} \equiv \sum_I e^{-\vec G_\alpha \cdot \vec R_I}.$$