Molecular Simulation/Treatment of Long Range Forces

Treatment of long-range interactions in molecular simulations poses a serious challenge in terms of computational requirements. By employing periodic boundary conditions, bulk properties of a sample may be calculated efficiently using a limited number of molecules within a finite cell. Each molecule in the cell interacts with every other molecule, both within the same cell and within the other infinite periodic cell images. These interactions include dispersion, repulsion, and any other electrostatic interactions. Essentially, this means that there are an infinite number of intermolecular interactions that must be calculated during the simulation, which is a problem in practical computation.

As is common in molecular simulation, the Lennard-Jones potential is used to describe both dispersion and repulsion interactions. These interactions, relative to other electrostatic interactions, decay fairly quickly, and can thus be dealt with fairly easily in computation. Electrostatic interactions (such as charge-dipole, dipole-dipole, quadrupole-quadrupole, etc.) are often significantly stronger at a given distance $r_{}$, and decay at a much slower rate as $r_{}$ increases, with the rate of decay increasing as the exponent of $r_{}$ increases. This makes their treatment somewhat more complicated.

Cutoffs

The unmodified Lennard-Jones potential (blue) and the Lennard-Jones potential with a smooth potential-based cutoff using a switching function (red). The switching function is active in the grey region. The switched-potential is zero beginning in the blue region.

Non-bonded cutoffs are employed so as to save computational resources without sacrificing validity of results. These cutoffs are used to model long-range interactions between atoms or molecules in a simulation, and essentially truncate the interactions at a specific distance. Because intermolecular interaction potentials converge to a particular value (namely zero) as $r\to \infty$, the interaction potential can effectively be tuned to this value after a certain distance. In the case of dispersion and repulsion interactions, which are modelled by the Lennard-Jones potential, a switching function, $S(r)$, is employed in the simulation. This switching function causes the decay of the interaction to reach zero by a desired cutoff distance, at which point the two molecules or atoms involved in the intermolecular interaction do not interact through the Lennard-Jones potential at all,

$\nu (r)={\begin{cases}\nu _{LJ}(r)&{\text{ for }}r\leq r_{switch}\\\nu _{LJ}(r)S(r)&{\text{ for }}r_{switch}\leq r\leq r_{cutoff}\\0&{\text{ for }}r>r_{cutoff}\end{cases}}$

The intermolecular distance at which the switching function is turned on is the switch distance, which is generally about 2 Å less than the cutoff distance. In between the switch distance and the cutoff distance, the interaction potential is gradually scaled to zero by the switching function. This gradual change in interaction potential maintains the physical basis of the simulation without using an abrupt, non-physical cutoff of interactions (as can be seen in the hard-sphere model for repulsive interactions). This method avoids a discontinuity in the interaction potential that would be caused by an abrupt truncation of the interaction potential function.

Such a convenient treatment of long-range interactions is easily applied to dispersion and repulsion interactions modelled by the Lennard-Jones potential, but is not applicable in all cases. For electrostatic interactions, other methods (such as Ewald summation, see below) must be employed. This is because electrostatic interactions are typically much longer-range than dispersion and repulsion interactions, and truncation through the use of simple switching functions could result in a serious decrease in accuracy of results. These interactions are stronger than dispersion or repulsion, and must be treated more carefully.

Ewald Decomposition

While truncation of the Lennard-Jones potential for dispersion and repulsion interactions is made relatively simple through the use of switching functions and non-bonded cutoffs, dealing with long-range electrostatic interactions is less straightforward, and requires a different approach. Evaluation of these interactions is commonly performed using the Ewald method, which divides electrostatic interactions into short-range and long-range components. The short-range component quickly decays to zero in real space, thus it can be truncated using a switching function, as was the case for the Lennard-Jones potential. The long-range component is more difficult to deal with, though it can be calculated efficiently in reciprocal space. In the Ewald summation method, the Gauss error function (erf) and complementary error function (erfc) are employed. These functions can be seen as

The Ewald decomposition divides the Coulombic interaction into a short-range real space component and a long range reciprocal space component by factoring it with the erfc(r) and erf(r) identity. The real-space erfc(r)/r component converges to zero at short distances, avoiding truncation error.

${\text{erf}}(x)={\frac {2}{\sqrt {\pi }}}\int \limits _{0}^{x}e^{-t^{2}}dt$

and

${\text{erfc}}(x)=1-{\frac {2}{\sqrt {\pi }}}\int \limits _{0}^{x}e^{-t^{2}}dt$

with ${\text{erf}}(0)=0$ and ${\text{erfc}}(0)=1$. These functions are thus complementary. Also, note that

$\lim _{x\to \infty }{\text{erf}}(x)=1$

and

$\lim _{x\to \infty }{\text{erfc}}(x)=0$

thus, by definition,

${\textrm {erfc}}(x)+{\textrm {erf}}(x)=1$.

In using these Gauss functions, the ${\frac {1}{r}}$ term of the Coulombic interaction potential is broken down into the aforementioned short- and long-range terms,

${\frac {1}{r}}={\frac {1}{r}}\left[{\textrm {erfc}}(\alpha r)+{\textrm {erf}}(\alpha r)\right]={\frac {{\textrm {erfc}}(\alpha r)}{r}}+{\frac {{\textrm {erf}}(\alpha r)}{r}}$

where $\alpha$ is a tunable constant chosen to ensure the electrostatic interaction decays to zero at the desired $r$ value. In this case, the ${\frac {{\textrm {erfc}}(\alpha r)}{r}}$ term corresponds to the short-range interaction potential, which decays quickly to zero through evaluation in real space, while the ${\frac {{\textrm {erf}}(\alpha r)}{r}}$ term corresponds to the long-range interaction potential, which is more easily evaluated in reciprocal space. This potential is evaluated as a sum of all possible charge-charge interactions across the infinite periodic images of the simulation cell. In reciprocal space, these interactions converge quickly, allowing efficient computation through application of fast Fourier transforms. The sum of all charge-charge interactions at long range throughout all periodic images of the simulation system can be represented as

$\sum _{\mathbf {S} }{\frac {q_{1}q_{2}{\text{erf}}(\alpha |\mathbf {r_{ij}} +\mathbf {S} |)}{|\mathbf {r_{ij}} +\mathbf {S} |}}$

where $\mathbf {r_{ij}}$ is the distance vector between interacting charges $q_{1}$ and $q_{2}$ and $\mathbf {S}$ is the translational vector between different periodic cell images. Equivalently, this sum can be written as a Fourier series,

${\frac {1}{L^{3}}}\sum _{\mathbf {g} }C_{\mathbf {g} }e^{i\mathbf {g} \cdot \mathbf {r} }$

where $C_{\mathbf {g} }$ is a reciprocal space expansion coefficient and $\mathbf {g}$ is a reciprocal space vector. Note that while this sum is still infinite, $C_{\mathbf {g} }\to 0$ quickly as $\mathbf {g} \to \infty$, allowing efficient computation of long-range electrostatic interactions through application of fast Fourier transforms.

Particle Mesh Ewald

The particle mesh Ewald method is a numerical approximation to Ewald summation. It involves a "smearing" of charges across a finite set of lattice points, located on a grid (or "mesh") that is placed over the simulation cell. By localizing charge density on this grid, the overall charge distribution across the simulation cell and its infinite periodic images is reduced to a more manageable form. This makes calculation of long-range electrostatic interactions significantly easier and less expensive, particularly in systems with a large simulation cell length, and thus, a large number of $\mathbf {g}$ vectors. From the lattice of "smeared" charge densities, long-range interaction potentials can be calculated using a fast Fourier transform. While this approximation decreases the computational cost significantly, accuracy of collected data is still high. The particle mesh Ewald method is also relatively simple to integrate into a simulation. These factors combine to make it an extremely attractive method to calculate long-range electrostatic interaction potentials.

See Also

Wikipedia article on Ewald summation