# How do water models fit into force field methods?

What I Understand:

In molecular mechanics, it is recognised that molecules tend to be composed of units that are structurally similar in a vast array of environments so atom types are defined from which molecular models can be built.

A force field is defined that describes the energy required to distort the molecule from equilibrium.[1]

$$E_\mathrm{FF} = E_\mathrm{str} + E_\mathrm{bend} + E_\mathrm{tors} + E_\mathrm{vdw} + E_\mathrm{el} + E_\mathrm{cross}$$

Once a functional form has been decided on, the field is parameterised based on empirical results.

What I Don't:

A large array of force fields have been developed, each different in their functional form or parameterisation. Numerous "water models" exist, and these appear to be force fields in their own right, often with far fewer energy terms.

The SPC model, for example, appears to contain only Coloumbic and Lennard-Jones terms, as described in Berendsen.[2]

We developed a three—point charge model (on hydrogen and oxygen positions) with a Lennard— Jones 6—12 potential on the oxygen positions only.

What is the purpose of these water models?

Question:

Why are standalone water models required, and how are they incorporated into more general forcefields like MM2?

References

1. F. Jensen, Introduction to Computational Chemistry, John Wiley & Sons, Chichester, 2007.
2. Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, W. F., and Hermans, J., Interaction Models for Water in Relation to Protein Hydration, in: Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry Held in Jerusalem, Israel, April 13--16, 1981. DOI: 10.1007/978-94-015-7658-1_21; pdf mirror
• Are you familiar with continuum models and their shortcomings? Mar 1, 2019 at 16:37
• The reason I ask is because you wondered why standalone models were required, and I assumed this could mean, why not use something simpler, like a continuum. Perhaps you meant why are force field parameters from other molecules not used? Mar 1, 2019 at 16:50
• I suppose I'm looking for the latter. Why isn't water taken care of in, for example, the MM2 forcefield, like every other molecule in the simulation. Why does it require its own forcefield? And then, how does that run alongside the MM2 part? Mar 1, 2019 at 16:55
• I don't know the history of mm2 or force-field dev to give you a detailed answer, but I know water is hard to simulate, because the point is to replicate properties. Capturing properties over a wide range of T/P for water, and including also behavior in a range of environments, is hard because you want to minimize computational costs. TIP3P, SPC/E etc all vary in accuracy under different conditions. Parameters were optimized for selected conditions. Most MD packages allow choice of different water models. Mar 1, 2019 at 17:43
• Using a special parameterization optimized for a selected molecule can be done by providing a parameter labels/files that the program understands, and if the program has rules for mixing molecules with different parameters (particularly LJ). Mar 1, 2019 at 17:48

I guess for starters, try finding a torsion angle on a 3 point molecule.

## Is it fast enough, accurate enough?

The real issue is why do some force-fields have intra-molecular bonds and others do not.

Including bond and angle forces in water models make their computation more expensive. There are classes of flexible forcefields for water. There are even more indepth forcefields that are polarizable.

Why don't we use the most sophisticated forcefield models for water?

The more sophisticated the model, the longer our simulations take. Often, rigid water models are good enough for certain properties. So, we use the simplest model that gives us the right answer for the property we are after. This is just a general good rule.

## Parameters

Force-fields use parameters, and proteins use ALOT of parameters. These parameters are often fit for a specific water model. Hence, you are kind of stuck using that water model if you are using that Force-Field for the protein/solute molecule.

## Complexity

As our interests move towards more involved calculations that involve water it does become important to include more effects such as flexibility and polarizability. This is increasingly the case in calculating Free Energies, whether it is the Free Energy of hydration or the Free Energy of binding a ligand to a protein. Another example is modelling supercritical water - only polarizable models have been successful. Rigid models such as SPC/E and TIP3P simply are not good enough for many properties under many conditions. You can of course fudge parameters to make anything fit a particular problem, but we want to use generalizable force-fields as much as possible. In this case, models with more parameters (flexible, polarisable) give us additional accuracy for the difficult problems.

## Bandwagon

Another reason many people use one model more than another is because their colleagues use it, so they use it because other people use it and it likely won't be questioned since other people use it. This is not a good reason for using a certain model! usually the first people to use a certain model had a reason, it solved their problem. However often the following bandwagon users don't read the strength/weaknesses of the model but instead just use it because they see it used lots for seemingly similar areas of research. Welcome to molecular modelling.

## Mixing Rules

Generally atoms have dispersion and electrostatic terms. In the case of dispersion, when the force or energy between two atoms is being calculated a single parameter for each term is required for the calculation. Take the well known Lennard-Jones for instance

$$U(r_{ij}) = 4 \epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}}\right)^6 \right]$$

different forcefields may have different numeric values for $$\epsilon_i$$, $$\epsilon_j$$, $$\sigma_i$$,$$\sigma_j$$ but they will all have a value. You need to use mixing rules to take the respective $$i$$ and $$j$$ parameters and create an $$\epsilon_{ij}$$ and $$\sigma_{ij}$$. We call these mixing rules a.k.a. combining rules. Wikipedia lists several of the common ones. https://en.wikipedia.org/wiki/Combining_rules

For instance for the Lorentz-Berthelot $$\epsilon_{ij} = \sqrt{\epsilon_{i} \epsilon_{j} }$$ and $$\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}$$

Finally, in SPC and SPC/E the hydrogens do not have Lennard-Jones parameters. It is not uncommon for hydrogens to have no dispersion parameters. This can lead to problems though, since dispersion parameters keep other atoms/particles from getting too close.

So a specific example of an SPC water molecule interacting with a Sodium ion, the electrostatics just needs to know that the Oxygen has a negative charge of $$q_o =−0.82$$ and the sodium has a charge of $$q_{N_a} =+1$$. A naive approach would be to treat this as a bare coulomb interaction which is just

$$U(r_{ij})^{\rm electrostatic} = k \frac{q_o q_{N_a}} {r_{ij}}$$

where

$$k = \frac{e^2}{4 \pi \beta}$$

$$e$$ is the charge of an electron and $$\beta$$ is the permittivity and a global and constant value for the simulation.

A better approach is to use EWALD's summation for electrostatics. It is more complicated, but the only charge specific parameters it needs to know is still just $$q_o$$ and $$q_{N_a}$$. The other parameters it uses for damping can be considered global parameters. Note: also do this for each of the hydrogen charges interacting with the Sodium ion.

For Lennard-Jones you would need to use mixing rules for the oxygen (there is no LJ parameters for the hydrogens)

$$\sigma_{(o-N_a)} = \frac{\sigma_o + \sigma_{N_a}}{2}$$

and

$$\epsilon_{(o-N_a)} = \sqrt{\epsilon_o \epsilon_{Na}}$$

then do the standard LJ calculation for force and or potential energy.

• This is a good overview of why a particular forcefield might be used, but what I'd like to know next is how are they integrated? For example, if I were to simulate an organic molecule solvated in water, where would each force field apply. Presumably the water-water interactions would be handled by SCP, and solute-solute interactions would be handled by MM2/CHARMM/..., but how are water-solute interactions handled? Mar 2, 2019 at 1:28
• Ahh, very perceptive! For the dispersion forces (typically Lennard jones) you use mixing rules. For electrostatic charge is charge.... There is no mixing rules. I will add this to my response when I have a keyboard rather than my mobile :) Mar 2, 2019 at 1:58
• In short you either take the average of the two parameters or you multiply them and take the square root. Mar 2, 2019 at 2:06
• I would also recommend using modern forcefields, not MM2. Check out Allen & Tildesley or Frenkel & Smit for references. Mar 2, 2019 at 3:29