Polarizable force fields

OpenMM 6.1 (just released!) added two new force fields: the AMOEBA 2013 force field and the CHARMM 2013 polarizable force field.  Both of them try in different ways to solve the same problem: providing a more realistic description of electrostatics than conventional force fields do.

Most force fields model electrostatics by assigning a "partial charge" to every atom.  This is supposed to represent the average amount of charge found near each atom.  The force field then treats atoms as point charges, using Coulomb's law to compute the interactions between them.

There are two main problems with this.  First, real atoms are not point charges.  The electrons are spread out in a diffuse cloud around the nuclei.  In many cases (covalent bonds, lone pairs, etc.) there are concentrations of charge density that aren't centered at the location of any atom.  Atom-centered partial charges only give a poor approximation to the electric field created by the real charge distribution.

Second, real systems are polarizable.  The electrons move around, responding to the presence of other nearby charges.  Every time the atoms move, the electrons rearrange themselves and the charge distribution changes.  In conventional force fields the partial charges are all fixed, so they cannot accurately reproduce polarization effects.

AMOEBA improves on this in two ways.  First, instead of just assigning a partial charge to each atom, it also assigns a dipole moment and a quadrupole moment.  This gives it much more power to describe the electric field around a molecule.  Second, it also assigns an "induced dipole" to every atom.  These are recomputed at each time step based on the electric field from all other atoms in the system.  This allows it to model polarization effects.

These improvements make AMOEBA much more accurate, but they also make it much more expensive to compute!  The equations for the interaction between two quadrupoles are enormously more complicated than those for the interaction between two point charges.  Also, the induced dipoles depend not just on fixed multipoles but also on other induced dipoles.  So every induced dipole affects every other induced dipole, and we need to use an iterative process to solve for them.  All of this makes AMOEBA very slow to work with.

The CHARMM polarizable force field takes a less expensive approach.  It sticks with point charges, but adds a small number of "virtual interaction sites" to represent lone pairs.  This lets it do a somewhat better job of describing the actual electron density, though not as accurate as AMOEBA.  It also models polarization by adding "Drude particles" to a subset of atoms.  These are extra point charges, tightly bound to the atoms by harmonic springs.  They never move far from the atoms' positions, but they can still move a little bit in response to the local electric field, thus creating dipole moments.

In principle one should solve for the Drude particle positions at each time step by minimizing their energy, a process that is about as expensive as computing the induced dipoles for AMOEBA.  OpenMM does offer this option, but it also offers a faster option: an "extended Lagrangian".  We treat the Drude particles as having mass, then let them move around under the influence of the force field, just like any other particles.  To make sure they stay near their minimum energy positions, we use a special integrator with two different Langevin thermostats.  It couples the center of mass of each atom/Drude-particle pair to the regular thermostat (at, say, 300K), but it couples their relative motion to a separate thermostat at very low temperature (say, 1K).  This keeps the Drude particles close to their minimum energy positions while adding very little cost to the computation.  It does, however, require us to use a somewhat smaller time step: usually no more than about 1 fs.