Reaction field method

The reaction field method is used in molecular simulations to simulate the effect of long range dipole-dipole interactions for simulations with periodic boundary conditions. Around each molecule there is a 'cavity' or sphere within which the Coulomb interactions are treated explicitly. Outside of this cavity the medium is assumed to have a uniform dielectric constant. The molecule induces polarization in this media which in turn creates a reaction field, sometimes called the Onsager reaction field. Although Onsager's name is often attached to the technique, because he considered such a geometry in his theory of the dielectric constant,^[1] the method was first introduced by Barker and Watts in 1973.^[2][3]

The effective pairwise potential becomes:

${\displaystyle U_{AB}=q_{A}q_{B}\left[{\frac {1}{r_{AB}}}+{\frac {(\varepsilon _{RF}-1)r_{AB}^{2}}{(2\varepsilon _{RF}+1)r_{c}^{3}}}\right]}$

where ${\displaystyle r_{c}}$ is the cut-off radius.

The reaction field in the center of the cavity is given by :

${\displaystyle E_{RF}={\frac {2(\varepsilon _{RF}-1)}{2\varepsilon _{RF}+1}}{\frac {\vec {M}}{r_{c}^{3}}}}$

where ${\displaystyle {\vec {M}}=\sum \mu _{i}}$ is the total dipole moment of all the molecules in the cavity. The contribution to the potential energy of the molecule ${\displaystyle i}$ at the center of the cavity is ${\displaystyle -1/2\mu _{i}\cdot E_{RF}}$ and the torque on molecule ${\displaystyle i}$ is simply ${\displaystyle \mu _{i}\times E_{RF}}$.

When a molecule enters or leaves the sphere defined by the cut-off radius, there is a discontinuous jump in energy.^[4] When all of these jumps in energy are summed, they do not exactly cancel, leading to poor energy conservation, a deficiency found whenever a spherical cut-off is used. The situation can be improved by tapering the potential energy function to zero near the cut-off radius. Beyond a certain radius ${\displaystyle r_{t}}$ the potential is multiplied by a tapering function ${\displaystyle f(r)}$. A simple choice is linear tapering with ${\displaystyle r_{t}=.95r_{c}}$, although better results may be found with more sophisticated tapering functions.

Another potential difficulty of the reaction field method is that the dielectric constant must be known a priori. However, it turns out that in most cases dynamical properties are fairly insensitive to the choice of ${\displaystyle \varepsilon _{RF}}$. It can be put in by hand, or calculated approximately using any of a number of well-known relations between the dipole fluctuations inside the simulation box and the macroscopic dielectric constant.^[4]

Another possible modification is to take into account the finite time required for the reaction field to respond to changes in the cavity. This "delayed reaction field method" was investigated by van Gunsteren, Berendsen and Rullmann in 1978.^[5] It was found to give better results—this makes sense, as without taking into account the delay, the reaction field is overestimated. However, the delayed method has additional difficulties with energy conservation and thus is not suitable for simulating an NVE ensemble.