Water model


In computational chemistry, a water model is used to simulate and thermodynamically calculate water clusters, liquid water, and aqueous solutions with explicit solvent. The models are determined from quantum mechanics, molecular mechanics, experimental results, and these combinations. To imitate a specific nature of molecules, many types of models have been developed. In general, these can be classified by the following three points; the number of interaction points called site, whether the model is rigid or flexible, whether the model includes polarization effects.
An alternative to the explicit water models is to use an implicit solvation model, also termed a continuum model, an example of which would be the COSMO solvation model or the polarizable continuum model or a hybrid solvation model.

Simple water models

The rigid models are considered the simplest water models and rely on non-bonded interactions. In these models, bonding interactions are implicitly treated by holonomic constraints. The electrostatic interaction is modeled using Coulomb's law, and the dispersion and repulsion forces using the Lennard-Jones potential. The potential for models such as TIP3P and TIP4P is represented by
where kC, the electrostatic constant, has a value of 332.1 Å·kcal/ in the units commonly used in molecular modeling; qi and qj are the partial charges relative to the charge of the electron; rij is the distance between two atoms or charged sites; and A and B are the Lennard-Jones parameters. The charged sites may be on the atoms or on dummy sites. In most water models, the Lennard-Jones term applies only to the interaction between the oxygen atoms.
The figure below shows the general shape of the 3- to 6-site water models. The exact geometric parameters vary depending on the model.

2-site

A 2-site model of water based on the familiar three-site SPC model has been shown to predict the dielectric properties of water using site-renormalized molecular fluid theory.

3-site

Three-site models have three interaction points corresponding to the three atoms of the water molecule. Each site has a point charge, and the site corresponding to the oxygen atom also has the Lennard-Jones parameters. Since 3-site models achieve a high computational efficiency, these are widely used for many applications of molecular dynamics simulations. Most of the models use a rigid geometry matching that of actual water molecules. An exception is the SPC model, which assumes an ideal tetrahedral shape instead of the observed angle of 104.5°.
The table below lists the parameters for some 3-site models.
TIPSSPCTIP3PSPC/E
r, Å0.95721.00.95721.0
HOH, deg104.52109.47104.52109.47
A, 103 kcal Å12/mol580.0629.4582.0629.4
B, kcal Å6/mol525.0625.5595.0625.5
q−0.80−0.82−0.834−0.8476
q+0.40+0.41+0.417+0.4238

The SPC/E model adds an average polarization correction to the potential energy function:
where μ is the electric dipole moment of the effectively polarized water molecule, μ0 is the dipole moment of an isolated water molecule, and αi is an isotropic polarizability constant, with a value of. Since the charges in the model are constant, this correction just results in adding 1.25 kcal/mol to the total energy. The SPC/E model results in a better density and diffusion constant than the SPC model.
The TIP3P model implemented in the CHARMM force field is a slightly modified version of the original. The difference lies in the Lennard-Jones parameters: unlike TIP3P, the CHARMM version of the model places Lennard-Jones parameters on the hydrogen atoms too, in addition to the one on oxygen. The charges are not modified. Three-site model has better performance in calculating specific heats.

Flexible SPC water model

The flexible simple point-charge water model is a re-parametrization of the three-site SPC water model. The SPC model is rigid, whilst the flexible SPC model is flexible. In the model of Toukan and Rahman, the O–H stretching is made anharmonic, and thus the dynamical behavior is well described. This is one of the most accurate three-center water models without taking into account the polarization. In molecular dynamics simulations it gives the correct density and dielectric permittivity of water.
Flexible SPC is implemented in the programs MDynaMix and Abalone.

Other models

The four-site models have four interaction points by adding one dummy atom near of the oxygen along the bisector of the HOH angle of the three-site models. The dummy atom only has a negative charge. This model improves the electrostatic distribution around the water molecule. The first model to use this approach was the Bernal–Fowler model published in 1933, which may also be the earliest water model. However, the BF model doesn't reproduce well the bulk properties of water, such as density and heat of vaporization, and is thus of historical interest only. This is a consequence of the parameterization method; newer models, developed after modern computers became available, were parameterized by running Metropolis Monte Carlo or molecular dynamics simulations and adjusting the parameters until the bulk properties are reproduced well enough.
The TIP4P model, first published in 1983, is widely implemented in computational chemistry software packages and often used for the simulation of biomolecular systems. There have been subsequent reparameterizations of the TIP4P model for specific uses: the TIP4P-Ew model, for use with Ewald summation methods; the TIP4P/Ice, for simulation of solid water ice; and TIP4P/2005, a general parameterization for simulating the entire phase diagram of condensed water.
Most of four-site water models use OH distance and HOH angle matching that of the free water molecule. An exception is OPC model, on which no geometry constraints are imposed other than the fundamental C2v molecular symmetry of the water molecule. Instead, the point charges and their positions are optimized to best describe the electrostatics of the water molecule. OPC reproduces a comprehensive set of bulk properties more accurately than commonly used rigid n-site water models. OPC model is implemented in AMBER force field.
BFTIPS2TIP4PTIP4P-EwTIP4P/IceTIP4P/2005OPCTIP4P-D
r, Å0.960.95720.95720.95720.95720.95720.87240.9572
HOH, deg105.7104.52104.52104.52104.52104.52103.6104.52
r, Å0.150.150.150.1250.15770.15460.15940.1546
A, 103 kcal Å12/mol560.4695.0600.0656.1857.9731.3865.1904.7
B, kcal Å6/mol837.0600.0610.0653.5850.5736.0858.1900.0
q−0.98−1.07−1.04−1.04844−1.1794−1.1128−1.3582−1.16
q+0.49+0.535+0.52+0.52422+0.5897+0.5564+0.6791+0.58

Others:
The 5-site models place the negative charge on dummy atoms representing the lone pairs of the oxygen atom, with a tetrahedral-like geometry. An early model of these types was the BNS model of Ben-Naim and Stillinger, proposed in 1971, soon succeeded by the ST2 model of Stillinger and Rahman in 1974. Mainly due to their higher computational cost, five-site models were not developed much until 2000, when the TIP5P model of Mahoney and Jorgensen was published. When compared with earlier models, the TIP5P model results in improvements in the geometry for the water dimer, a more "tetrahedral" water structure that better reproduces the experimental radial distribution functions from neutron diffraction, and the temperature of maximal density of water. The TIP5P-E model is a reparameterization of TIP5P for use with Ewald sums.
BNSST2TIP5PTIP5P-E
r, Å1.01.00.95720.9572
HOH, deg109.47109.47104.52104.52
r, Å1.00.80.700.70
LOL, deg109.47109.47109.47109.47
A, 103 kcal Å12/mol77.4238.7544.5554.3
B, kcal Å6/mol153.8268.9590.3628.2
q−0.19562−0.2357−0.241−0.241
q+0.19562+0.2357+0.241+0.241
RL, Å2.03792.0160--
RU, Å3.18773.1287--

Note, however, that the BNS and ST2 models do not use Coulomb's law directly for the electrostatic terms, but a modified version that is scaled down at short distances by multiplying it by the switching function S:
Thus, the RL and RU parameters only apply to BNS and ST2.

6-site

Originally designed to study water/ice systems, a 6-site model that combines all the sites of the 4- and 5-site models was developed by Nada and van der Eerden. Since it had a very high melting temperature when employed under periodic electrostatic conditions, a modified version was published later optimized by using the Ewald method for estimating the Coulomb interaction.

Other

The computational cost of a water simulation increases with the number of interaction sites in the water model. The CPU time is approximately proportional to the number of interatomic distances that need to be computed. For the 3-site model, 9 distances are required for each pair of water molecules. For the 4-site model, 10 distances are required. For the 5-site model, 17 distances are required. Finally, for the 6-site model, 26 distances are required.
When using rigid water models in molecular dynamics, there is an additional cost associated with keeping the structure constrained, using constraint algorithms.