Ewald sum and SPME

This section outlines how to optimise the accuracy of the Ewald sum parameters for a given simulation. In what follows the directive spme may be used anywhere in place of the directive ewald if the user wishes to use the Smoothed Particle Mesh Ewald method.

As a guide to beginners DL_POLY_2 will calculate reasonable parameters if the ewald precision directive is used in the CONTROL file (see section 4.1.1). A relative error (see below) of 10$^{-6}$ is normally sufficient so the directive

ewald precision 1d-6

will cause DL_POLY_2 to evaluate its best guess at the Ewald parameters $\alpha$, kmax1, kmax2 and kmax3. (The user should note that this represents an estimate, and there are sometimes circumstances where the estimate can be improved upon. This is especially the case when the system contains a strong directional anisotropy, such as a surface.) These four parameters may also be set explicitly by the ewald sum directive in the CONTROL file. For example the directive

ewald sum 0.35 6 6 8

would set $\alpha= 0.35$ Å$^{-1}$, kmax1 = 6, kmax2 = 6 and kmax3 = 8. The quickest check on the accuracy of the Ewald sum is to compare the Coulombic energy ($U$) and the coulombic virial ($\cal W$) in a short simulation. Adherence to the relationship $U = -{\cal W}$ shows the extent to which the Ewald sum is correctly converged. These variables can be found under the columns headed eng_cou and vir_cou in the OUTPUT file (see section 4.2.2).

The remainder of this section explains the meanings of these parameters and how they can be chosen. The Ewald sum can only be used in a three dimensional periodic system. There are three variables that control the accuracy: $\alpha$, the Ewald convergence parameter; $r_{\rm cut}$ the real space forces cutoff; and the kmax1,2,3 integers ^3.1 that effectively define the range of the reciprocal space sum (one integer for each of the three axis directions). These variables are not independent, and it is usual to regard one of them as pre-determined and adjust the other two accordingly. In this treatment we assume that $r_{\rm cut}$ (defined by the cutoff directive in the CONTROL file) is fixed for the given system.

The Ewald sum splits the (electrostatic) sum for the infinite, periodic, system into a damped real space sum and a reciprocal space sum. The rate of convergence of both sums is governed by $\alpha$. Evaluation of the real space sum is truncated at $r=r_{\rm cut}$ so it is important that $\alpha$ be chosen so that contributions to the real space sum are negligible for terms with $r>r_{\rm cut}$. The relative error ($\epsilon$) in the real space sum truncated at $r_{\rm cut}$ is given approximately by

$$\epsilon \approx {\rm erfc}(\alpha r_{\rm cut})/r_{\rm cut} \approx \exp[-(\alpha.r_{\rm cut})^2]/r_{\rm cut}$$ (3.1)

The recommended value for $\alpha$ is 3.2/$r_{\rm cut}$ or greater (too large a value will make the reciprocal space sum very slowly convergent). This gives a relative error in the energy of no greater than $\epsilon = 4\times 10^{-5}$ in the real space sum. When using the directive ewald precision DL_POLY_2 makes use of a more sophisticated approximation:

$${\rm erfc}(x) \approx 0.56 \exp(-x^2)/x$$ (3.2)

to solve recursively for $\alpha$, using equation 3.1 to give the first guess.

The relative error in the reciprocal space term is approximately

$$\epsilon \approx \exp(- k_{max}^2/4\alpha^2)/k_{max}^2$$ (3.3)

where

$$k_{max} = \frac{2\pi}{L}~{\tt kmax}$$ (3.4)

is the largest $k$-vector considered in reciprocal space, $L$ is the width of the cell in the specified direction and kmax is an integer.

For a relative error of $4\times 10^{-5}$ this means using $k_{max} \approx 6.2 \alpha$. kmax is then

$${\tt kmax} > 3.2~L/r_{\rm cut}$$ (3.5)

In a cubic system, $r_{\rm cut}~=~L/2$ implies ${\tt kmax}~=~7$. In practice the above equation slightly over estimates the value of kmax required, so optimal values need to be found experimentally. In the above example kmax = 5 or 6 would be adequate.

If your simulation cell is a truncated octahedron or a rhombic dodecahedron then the estimates for the kmax need to be multiplied by $2^{1/3}$. This arises because twice the normal number of $k$-vectors are required (half of which are redundant by symmetry) for these boundary contributions [30].

If you wish to set the Ewald parameters manually (via the ewald sum or spme sum directives) the recommended approach is as follows. Preselect the value of $r_{\rm cut}$, choose a working a value of $\alpha$ of about $3.2/r_{\rm cut}$ and a large value for the kmax (say 10 10 10 or more). Then do a series of ten or so single step simulations with your initial configuration and with $\alpha$ ranging over the value you have chosen plus and minus 20%. Plot the Coulombic energy (and $-{\cal W}$) versus $\alpha$. If the Ewald sum is correctly converged you will see a plateau in the plot. Divergence from the plateau at small $\alpha$ is due to non-convergence in the real space sum. Divergence from the plateau at large $\alpha$ is due to non-convergence of the reciprocal space sum. Redo the series of calculations using smaller kmax values. The optimum values for kmax are the smallest values that reproduce the correct Coulombic energy (the plateau value) and virial at the value of $\alpha$ to be used in the simulation.

Note that one needs to specify the three integers (kmax1, kmax2, kmax3) referring to the three spatial directions, to ensure the reciprocal space sum is equally accurate in all directions. The values of kmax1, kmax2 and kmax3 must be commensurate with the cell geometry to ensure the same minimum wavelength is used in all directions. For a cubic cell set kmax1 = kmax2 = kmax3. However, for example, in a cell with dimensions $2A = 2B = C$ (ie. a tetragonal cell, longer in the c direction than the a and b directions) use 2kmax1 = 2kmax2 = (kmax3).

If the values for the kmax used are too small, the Ewald sum will produce spurious results. If values that are too large are used, the results will be correct but the calculation will consume unnecessary amounts of cpu time. The amount of cpu time increases with ${\tt kmax1}\times{\tt kmax2} \times {\tt kmax3}$.