Skip to main content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

SPC/E Water Reference Calculations - 10Å cutoff

In this section, we provide sample configurations of SPC/E Water molecules[1] and report the various energy and force calculations for those configurations, where the coulombic contributions are computed using the traditional Ewald Summation Method[2]. These sample configurations and reference calculations can be used to validate the energy and force routines for an existing or new molecular simulation code. The discussion that follows places particular emphasis on the Ewald sum method for computing the Coulomb potential. The cutoff radius for real-space terms in this data set is 10Å.

1. SPC/E Model

The SPC/E Model represents water as a triatomic molecule with rigid bonds, i.e., both the length of the oxygen-hydrogen bonds and the angle between those bonds are fixed. The triatomic molecule has a single Lennard-Jones site centered on the oxygen atom and point charges centered on each atom:

Schematic of SPC/E Water Molecule

Fig. 1: Schematic of the SPC/E Molecule.

Red and white balls indicate oxygen and hydrogen atoms, respectively. The charges, indicated in the figure, are +q for both hydrogen atoms and -2q for the oxygen atom. The dashed line is a circle of diameter σ centered on the oxygen atom, signifying the effective size of the molecule. The hydrogen-oxygen bondlength is 1 Å and the hydrogen-oxygen-hydrogen bondangle is 109.47°. The remaining model parameters are [1]:

$$q = 0.42380~e $$

$$\sigma = 0.316555789~\textrm{nm}$$

$$\epsilon / k_B=78.19743111~K$$

2. SPC/E Potential Energy

The potential energy for the SPC/E model of water is:

$$ \Large E_{total}\left( \mathbf{r}^N \right) = E_{disp}\left( \mathbf{r}^N \right) + E_{LRC} + E_{coulomb}\left( \mathbf{r}^N \right) $$

The terms on the right hand side of the equality are the pair dispersion energy, the long-range correction to the pair dispersion energy, and the Coulombic potential energy, respectively. In the SPC/E Model, the pair dispersion energy is a Lennard-Jones intermolecular potential energy, which is discussed elsewhere in the NIST Standard Reference Simulation Website. The Coulombic energy is computed using the Ewald Summation Method, which is discussed below.

3. Ewald Summation Method

The Coulombic energy, and long-range energies, may be computed using the Ewald Summation Technique, which replaces the true Coulomb potential with the following [2,3]:

$$\Large \begin{eqnarray} E_{coulomb}\left(\mathbf{r}^N\right) & = & \dfrac{1}{2}  {\sum\limits_{\mathbf{n}}}^{\dagger}   \sum\limits_{j=1}^N \sum\limits_{l=1}^N
\dfrac{q_j q_l}{ 4 \pi \epsilon_0} \dfrac{\text{erfc}\left(\alpha \cdot \left| \mathbf{r}_{jl} + \mathbf{n}\right| \right)}{\left| \mathbf{r}_{jl} + \mathbf{n}\right|}\\
&& +\dfrac{1}{2 \pi V} \sum\limits_{\mathbf{k} \neq \mathbf{0}} \dfrac{1}{k^2} \exp \left[-\left( \dfrac{\pi k}{\alpha} \right)^2 \right] \dfrac{1}{4 \pi \epsilon_0} \cdot \left|\sum\limits_{j=1}^N q_j \exp \left(2\pi i \mathbf{k} \cdot \mathbf{r}_j \right) \right|^2 \\
&& - \dfrac{\alpha}{\sqrt{\pi}} \sum\limits_{j=1}^N \dfrac{q_j^2}{4 \pi \epsilon_0} \\
&& -\dfrac{1}{2}  {\sum\limits_{j=1}^M}^{\dagger^{-1}}   \sum\limits_{\kappa=1}^{N_j} \sum\limits_{\lambda=1}^{N_j} \dfrac{q_{j_\kappa} q_{j_\lambda}}{ 4 \pi \epsilon_0} \dfrac{\text{erf}\left(\alpha \cdot \left| \mathbf{r}_{j_\kappa j_\lambda} \right| \right)}{\left| \mathbf{r}_{j_\kappa j_\lambda} \right|}
\end{eqnarray}$$

The terms on the right hand side of the equality are 1) the real-space term Ereal, 2) the Fourier-space term, Efourier, 3) the self correction term Eself, and 4) the intramolecular term Eintra. We note that this form of the Ewald Summation 1) requires total charge neutrality for the configuration and 2) neglects the surface dipole term (equivalent to using the "tin-foil" or conducting surface boundary condition). The meaning of symbols in this nist-equation are:

nLattice vector of periodic cell images
kFourier-space vector of periodic cell images
kmodulus of k ; k2 = |k|2
qjValue of charge at site j
αEwald damping parameter
NTotal number of charged sites
MTotal number of molecules
NjTotal number of charged sites in molecule j
κ, λIndices of sites in a single molecule
VVolume of the simulation cell, LxLyLz
ε0Permittivity of vacuum (see below)
iImaginary unit, (-1)1/2
rjCartesian vector coordinate of site j
rjlrj -rl
erf(x)Error Function computed for abscissa x
erfc(x)Complimentary Error Function computed for abscissa x

In this form, the superscript "†" (dagger) in Ereal indicates that the sum skips all pairs i=j inside the original simulation cell (n = 0). The superscript "†-1" in Eintra indicates that the sum is over site pairs within molecules in the original simulation cell. Additionally, the Fourier vectors (k) in this nist-equation are composed of integer elements, e.g. k = 2ex+ey+4ez where ei is the unit vector for Cartesian direction i. The Fourier space term can alternatively be written using k vectors with elements proportional to 2∏. In practice, the above nist-equation is not how the Ewald Summation is actually implemented. Typically, one makes the following assumptions/reductions to simplify the summation:

  1. The real-space sum is done for the original simulation cell only, i.e. n=0.
  2. Site-site interactions in both Edisp and Ereal are truncated at a pre-defined cutoff rcut, where rcut ≤ min(Lx,Ly,Lz)/2. In practice, the damping parameter, α, is chosen so that contributions to the real-space term are negligible for rij > rcut.
  3. The Fourier space summation is truncated at a pre-defined maximum k or maximum value of k2.

Thus, the practical implementation of the Ewald Summation is [3]:

$$\Large \begin{eqnarray} E_{coulomb}\left(\mathbf{r}^N\right) & = &
\sum\limits_{j}  \sum\limits_{l>j}
\dfrac{q_j q_l}{ 4 \pi \epsilon_0} \dfrac{\text{erfc}\left(\alpha \cdot  \left| \mathbf{r}_{jl} \right| \right)}{\left| \mathbf{r}_{jl} \right|} \Theta\left( r_{cut} - \left|\mathbf{r}_{jl}\right| \right) \\
&& +\dfrac{1}{2 \pi V} \sum\limits_{\mathbf{k} \neq \mathbf{0}} \dfrac{1}{k^2} \exp \left[-\left( \dfrac{\pi k}{\alpha} \right)^2 \right] \dfrac{1}{4 \pi \epsilon_0} \cdot \left|\sum\limits_{j=1}^N q_j \exp \left(2\pi i \mathbf{k} \cdot \mathbf{r}_j \right) \right|^2 \\
&& - \dfrac{\alpha}{\sqrt{\pi}} \sum\limits_j \dfrac{q_j^2}{4 \pi \epsilon_0} \\
&& - \sum\limits_{j=1}^M  \sum\limits_{\kappa}  \sum\limits_{\lambda>\kappa} \dfrac{q_{j_\kappa} q_{j_\lambda}}{ 4 \pi \epsilon_0} \dfrac{\text{erf}\left(\alpha \cdot  \left| \mathbf{r}_{j_\kappa j_\lambda} \right| \right)}{\left| \mathbf{r}_{j_\kappa j_\lambda} \right|} 
\end{eqnarray}$$

We note that the real-space term now includes multiplication by the Heaviside Step Function, Θ(rcut - rij), which functionally truncates that term at rij = rcut.

4. Parameters and Physical Constants for Reference Calculations herein

For the reference calculations given below, we use the following parameters and apply certain conditions to the calculation of both the dispersion interactions and the Ewald Summation:

α5.6 / min(Lx,Ly,Lz)
kmax5 ; also only include k for which k2 < kmax2 +2, i.e. k2 < 27.
rcut10 Å
Truncation 
 DispersionTruncate at rcut, apply analytic long-range corrections
 CoulombTruncate real-space term at rcut
Boundary ConditionsPeriodic and tin-foil (conducting) boundary conditions for all Cartesian Directions
erfc(x)Implementation of Numerical Recipes ERFCC function; Ref. 4, page 164.

The reference calculations given below were done using fundamental constants of physics and chemistry recommended by CODATA in 2010 [5,6]. We report these constants because the calculation of each contribution to the intermolecular energy will depend, ever so slightly, on the choice of fundamental physical constants and, in particular, the number of digits in those constants that are carried in the simulation. We use the full constants (untruncated) given in the CODATA 2010 recommendation:

NameSymbolValueUnits
Boltzmann ConstantkB1.3806488E-23J/K
Avogadro ConstantNa6.02214129E+23mol-1
Planck constanth6.62606957E-34J s
Elementary chargee1.602176565E-19C
Permittivity of Vacuumε08.854187817E-12C2/(J m)

5. Sample Configurations of SPC/E Water Molecules

Four sample configurations of SPC/E molecules are available for download as a gzipped tarball archive. This archive contains five files: the four sample configuration files and one metadata file that explains the format of the sample configurations. These configurations should be converted to the configuration file format native to a user's simulation software.

6. Reference Calculations of Intermolecular Energy for SPC/E

Configuration ->1234
M (number of SPC/E molecules)100200300750
Lx=Ly=Lz (Å)2.00000E+012.00000E+012.00000E+013.00000E+01
Edisp/kB (K)9.95387E+041.93712E+053.54344E+054.48593E+05
ELRC/kB (K)-8.23715E+02-3.29486E+03-7.41343E+03-1.37286E+04
Ereal/kB (K)-5.58889E+05-1.19295E+06-1.96297E+06-3.57226E+06
Efourier/kB (K)6.27009E+036.03495E+035.24461E+037.58785E+03
Eself/kB (K)-2.84469E+06-5.68938E+06-8.53407E+06-1.42235E+07
Eintra/kB (K)2.80999E+065.61998E+068.42998E+061.41483E+07
Etotal/kB (K)-4.88604E+05-1.06590E+06-1.71488E+06-3.20501E+06


Lastly, we note that we have reported six significant figures in each contribution to the intermolecular energy, while the software used for these calculations (a code written in FORTRAN) carried as many digits as possible within the specified machine precision. Thus, the reported Etotal may differ in the final digits from a user's calculation.

References

  1. H. J. C. Berendsen, J. R. Griger, and T. P. Straatsma, J. Phys. Chem., 91, 6269 (1987)
  2. P. Ewald, Ann. Phys., 369, 253 (1921)
  3. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, New York, 1989).
  4.  W. H. Press, et al., Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, New York, 1986).
  5. P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys., 84, 1527 (2012)
  6. P. J. Mohr, B. N. Taylor, and D. B. Newell, J. Phys. Chem. Ref. Data, 41, 043109 (2012)
Created January 9, 2013, Updated August 25, 2016