Lesson 4.5 · 14 min

Water models

In most biomolecular and electrolyte simulations, the majority of atoms are water. Simple rigid models with 3, 4 or 5 interaction sites reproduce different properties well or badly. Choosing (and not mixing) the right one matters.

3-, 4- and 5-site rigid water models, and the virtual site M
Why model dipoles are larger than the gas-phase value (effective polarisation)
How models compare for density, diffusion, dielectric constant and melting point, and how to choose one

The idea

Water is the solvent of life and of most electrolytes, and in a typical simulation box 80–95% of the atoms belong to water. A water model must be cheap (it is evaluated billions of times) yet reproduce the liquid's structure, dynamics and thermodynamics.

Rigid point-charge models are the workhorses. The geometry is fixed (Lesson 4.4, SETTLE), one Lennard-Jones site sits on the oxygen, and partial charges model the electrostatics:

3-site (TIP3P, SPC, SPC/E): charges on O and on both H atoms. Simple and cheap.
4-site (TIP4P, TIP4P/2005, TIP4P-Ew, OPC): the negative charge moves from O to a massless virtual site M on the H–O–H bisector, about 0.15 Å from O. This gives a better charge distribution (quadrupole). Forces on M are passed back to O and H.
5-site (TIP5P): two negative sites at the "lone pair" positions, for a tetrahedral charge distribution.

Effective polarisation. A water molecule in the liquid is polarised by its neighbours; its dipole grows from 1.85 D (gas) to an estimated 2.6–3.0 D. Fixed-charge models cannot respond, so their charges are boosted to give about 2.3–2.5 D: a compromise that works in the liquid but overbinds the gas-phase dimer. Polarizable models (SWM4, AMOEBA, Drude) and machine-learned potentials treat this explicitly.

No model gets everything right. TIP3P diffuses more than twice too fast and melts at 146 K, but it is the water that AMBER and CHARMM were parameterised with. SPC/E is a good all-rounder. TIP4P/2005 reproduces densities, the temperature of maximum density and the phase diagram well, but underestimates the dielectric constant. OPC was optimised for its charge distribution and gives a good dielectric constant.

Key idea. Use the water model your force field was parameterised with (TIP3P for most AMBER/CHARMM work, SPC for GROMOS, TIP4P-type for many OPLS and phase-behaviour studies), and know which properties it gets wrong.

The algorithm

Handling a 4-site model in an MD code:

every step, after moving O, H1, H2:
    M = O + a * ((H1 - O) + (H2 - O))          # rebuild the virtual site
compute forces on all sites (LJ on O, Coulomb on H1, H2, M)
spread the force on M back to the real atoms:
    F_O  += (1 - 2a) * F_M
    F_H1 += a * F_M;   F_H2 += a * F_M
integrate only O, H1, H2 (M has no mass), then SETTLE

a is fixed by the geometry: a = d_OM / |(H1 − O) + (H2 − O)|. Spreading this way conserves the total force and torque (checked in the code tab).
Models are fitted with particular simulation settings (cutoff, long-range electrostatics, LJ tail correction). TIP4P-Ew and OPC assume Ewald/PME.
Use the same model for every water in the system, including ion and ligand parameterisation: ion parameters (e.g. Joung–Cheatham) are water-model specific.
GROMACS: #include "amber99sb.ff/tip3p.itp" or tip4p.itp; build with gmx solvate -cs tip4p.gro. LAMMPS: pair_style lj/cut/tip4p/long for 4-site models.

See it inside LAMMPS

The same idea in the real code, run live in your browser and checked by hand.

A.8 · read_data files