Implementation and Evaluation

We will consider an Ewald implementation which is a modified version of the ewald code written for Berend Smit's Molecular Simulation course webpage. This code simply computes the Ewald energy for a cubic lattice, given an appropriate number of particles, and a value for $\sqrt{\alpha}$ (which is called $\alpha$ in the code), and a value for $k_{\rm max}$, the maximum integer index for enumerating $k$-vectors. ⁵

The units used in a system with electrostatics differ depending on community. So far, we have assumed that the units of electrostatic potential are charge $\mathscr{C}$, divided by length $\mathscr{L}$, because we write potential as $\phi = q/\vert{\bf r}\vert$, where $q$ is measured in units of $\mathscr{C}$ and distance in units of $\mathscr{L}$. Energy is therefore written in units of $\mathscr{C}^2$ over $\mathscr{L}$, and force in units of $\mathscr{C}^2$ over $\mathscr{L}^2$. If we want the final energy in more familiar units, we can choose $\mathscr{C}$ and $\mathscr{L}$, and use the standard prefactor $1/4\pi\epsilon_0$ to convert from ``charge squared per length'' to ``energy''. For example, in SI units, $\epsilon_0 = 8.85419\times10^{-12}$ (C$^2$/m)/J. In this implementation, we use a length of $\mathscr{L} = 1$ and $\mathscr{C} = 1$ and measure energy such that $1/4\pi\epsilon_0 = 1$.

We will examine two configurations, both with $N$ = 8$^3$ = 512 particles, with alternating + and - charges. One configuration has the particle on a cubic lattice with lattice spacing $r_0 = 1$, which is the standard NaCl crystal structure. We will call this the ``crystal'' configuration. The other is like the crystal, only each particle is displaced by a random amount from its lattice position with a maximum displacement of 0.3. We will call this the ``liquid'' configuration. We compute the total electrostatic energy via the Ewald sum technique for various values of $\sqrt{\alpha}$ and maximum $k$-vector index. As we increase the number of $k$-vectors taken in the sum, we would like to show that the total energy converges to a certain value. We will measure this in terms of the Madelung constant, $M$:

$$\mathscr{U}_{\rm Coul} = -\frac{N q^2 M}{4\pi\epsilon_0r_0}$$ (341)

The following results were obtained for the crystal when $\alpha = 0.8$:

 alpha kmax u_real    u_fourier     u_self    u_total     M
 1.20  2    -.206695  0.297020E-33  -0.67703  -0.88372    1.7674    
 1.20  3    -.206695  0.509667E-33  -0.67703  -0.88372    1.7674    
 1.20  4    -.206695  0.994016E-02  -0.67703  -0.87378    1.7476    
 1.20  5    -.206695  0.994016E-02  -0.67703  -0.87378    1.7476    
 1.20  6    -.206695  0.994016E-02  -0.67703  -0.87378    1.7476

This shows that, for the perfectly periodic crystal, very few $k$-vectors are needed to reach a converged energy.

For the liquid, the results are somewhat different, again for $\alpha = 1.2$:

 alpha kmax u_real    u_fourier     u_self    u_total     M
 1.20  2   -.244210   0.210535E-01  -0.67703  -0.90018    1.8004    
 1.20  3   -.244210   0.349887E-01  -0.67703  -0.88625    1.7725    
 1.20  4   -.244210   0.482931E-01  -0.67703  -0.87294    1.7459    
 1.20  5   -.244210   0.510164E-01  -0.67703  -0.87022    1.7404    
 1.20  6   -.244210   0.519170E-01  -0.67703  -0.86932    1.7386    
 1.20  7   -.244210   0.521388E-01  -0.67703  -0.86910    1.7382    
 1.20  8   -.244210   0.521741E-01  -0.67703  -0.86906    1.7381    
 1.20  9   -.244210   0.521802E-01  -0.67703  -0.86906    1.7381    
 1.20 10   -.244210   0.521808E-01  -0.67703  -0.86906    1.7381

Considering the converged results for various values of $\alpha$, we see that $M$ is not too sensitive to $\alpha$, once $\alpha > 0.6$:

Madelung constant $M$ vs. Ewald parameter $\alpha$ for a system of 512 ions.

Note that, in principle, $\mathscr{U}_{\rm Coul}$ doesn't depend on $\alpha$. However, errors due to too low a choice for $\alpha$ indicate that the Fourier sum is poorly converged, because of the factor $\exp[-k^2/4\alpha]$.