2
$\begingroup$

I generally work on fluid flow, so some basic concepts in chemistry are new to me. Work I'm doing requires use of molecular dynamics (MD) to do classical potential simulations of water–NaCl systems but this question is about water models more than simulation and software details.

The simulations require definition of a "water model." There's apparently one already implemented in an example/tutorial I'm using (below). I want to compare the results of simulations based on different water models. My problem is that I'm having a lot of trouble understanding how to change my water model (which I think is SPC/E) because I don't really understand what it is.

I have found tables that provide parameters for different water models, like this table at Wikipedia: https://en.wikipedia.org/wiki/Water_model

However, when I look at my software, I can't relate what I'm seeing to the table. I can't find an example of where someone has taken information from this sort of table and explained how to implement it in a simulation. So my first request is, please point me to anything out there! Here's the part of my software that I think defines my water model:

&FORCEFIELD
  &BEND
    ATOMS H O H
    K 0.
    THETA0 1.8
  &END BEND
  &BEND
    ATOMS O H H
    K 0.
    THETA0 1.8
  &END BEND
  &BOND
    ATOMS O H
    K 0.
    R0 1.8
  &END BOND
  &BOND
    ATOMS H H
    K 0.
    R0 1.8
  &END BOND
  &CHARGE
    ATOM O
    CHARGE -0.8476
  &END CHARGE
  &CHARGE
    ATOM H
    CHARGE 0.4238
  &END CHARGE
  &NONBONDED
    &LENNARD-JONES
      ATOMS O O
      EPSILON 78.198 ! this is K, = 0.155 kcal/mol = 0.650 kJ/mol
      SIGMA 3.166
      RCUT 11.4
    &END LENNARD-JONES
  &LENNARD-JONES
      ATOMS O H
      EPSILON 0.0
      SIGMA 3.6705
      RCUT 11.4
    &END LENNARD-JONES
    &LENNARD-JONES
      ATOMS H H
      EPSILON 0.0
      SIGMA 3.30523
      RCUT 11.4
    &END LENNARD-JONES
  &END NONBONDED
&END FORCEFIELD

I can generally see that everything up to &NONBONDED relates to water somehow because it's about charge, H, and O (I have no idea what "K" is). But what about the "LENNARD-JONES" stuff? I know in general what a Lennard Jones model is, but cannot see how that would be related to water models. If it is, how? Thank you for any advice on how to understand this.

$\endgroup$
  • $\begingroup$ jonsca, sorry for leaving another question, but my reputation doesn't allow me to comment on your comment. I don't see a checkmark icon. I know what it should look like because I had an account here in the past. If you could possibly mark the answer as correct, that would be great. I would like to make sure that the answerer gets credit for his great explanation. Thank you. $\endgroup$ – Antst Sep 26 '18 at 18:18
3
$\begingroup$

I can explain some of this, and actually some of it is explained on the Wikipedia page you referenced.

Most water models represent the oxygen using a Lennard-Jones potential, which encapsulates the repulsive interaction that stops the molecules getting too close (this is physically because the electron shells don't want to overlap each other) and a long-range attraction due to Van der Waals forces. The Lennard-Jones potential can be parameterized in several ways, two of which are: $$ v_{\text{LJ}}(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] = \frac{A}{r^{12}} - \frac{B}{r^6} $$ The first of these appears in your force field specification file. The $\epsilon$ parameter is actually given in Kelvin, which is $\epsilon/k_{\text{B}}$, where $k_{\text{B}}$ is Boltzmann's constant. Hence the note "this is K". This may seem an odd choice, but there are some historical reasons. However, you'll see on the same line, in a comment for the user, that this number is translated into SI units and also into old-time cgs units of kcal/mol. The $\sigma$ parameter is also given in your file, in old-style Angstroms. Just to be clear, this is the interaction between oxygen atoms, hence it is described as the O-O Lennard-Jones term.

The second form appears on the Wikipedia page. The values of $A$ and $B$ appear in the table, in old-time units. If you do the necessary multiplications you can check that (to within a few significant figures, anyway), $4\epsilon\sigma^{12}$ and $4\epsilon\sigma^{6}$ do indeed agree with the tabulated $A$ and $B$, taking $\epsilon$ in kcal/mol and $\sigma$ in Angstroms, for the SPC/E model.

There are some redundant entries in your file: BEND and BOND. For a molecule like this, one could allow the OH bonds to vibrate, and the HOH bond angle to vary, and these terms would need some force constants and equilibrium bond lengths and angles. But SPC/E is a rigid model: the bond lengths and angles are fixed, so it looks like the force constants ($K$) have been set to zero.

One might imagine that the H atoms deserve also to be represented as Lennard-Jones atoms, but they are quite small, and in this model this aspect is neglected. Hence we see that $\epsilon$ is zero for the H-O and H-H interactions.

So then, all that is needed to complete the model is the charge distribution, and it looks like you've identified the three charges, one on each atom, characterizing SPC/E.

This isn't the most sophisticated model, as you can see if you scroll down that page.

For more details concerning the input file, I would hope that your software package and/or exercise instructions, give enough information. Hope this helps at least to understand the correspondence between the numbers in your file, and the numbers on the Wikipedia page.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.