7.3 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. (All of Prof. Smit’s codes are available in the FrenkelSmitCodes directory of the instructional-codes respository.) 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/|{\bf r}|\), 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\times 10^{-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 Self Part 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 \(1/\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\): \begin{equation} \mathscr {U}_{\rm Coul} = -\frac {N q^2 M}{4\pi \epsilon _0r_0} \end{equation}

Table 2 shows results of Ewald summation for the perfect lattice, and Table 3 shows resuts for the “liquid”. We see several interesting things from these example calculations:

1.

The self-energy is the dominant contributor from the long-range compensating charge smears.

2.

The Fourier-space contribution is relatively small, and is much smaller for the perfect lattice than for the liquid. This means the precision of the Fourier transform of the charge distribution (which is larger for larger \(k_{\rm max}\) is relatively unimportant to the calculation.

3.

The overall Coulomb energy is insensitive to the real-space cutoff for the perfect lattice, while it is weakly decreasing in magnitude (it is negative overall) for increasing real-space cutoff in the liquid.