<p dir="ltr" style="text-align: left;">Over the past couple days I've
attempted to implement the water model developed by Ishiyama &
Morita (J. Phys. Chem. C 2007, 111, 721-737) using a combination of the
shell polarization technique and tabulated potentials for the bonded
terms, among other things. Before elaborating below, the documentation
I found was pretty weak for the tabulated bond potentials, and forum
posts were limited in scope, so I'd like to gather some information
from those in the know such that the wiki and documentation may be
updated on this front. Furthermore, as I slowly understand and make
progress, I'll most likely post follow-ups with more questions to round
out the discussion.<br>
</p>
<p dir="ltr" style="text-align: left;">From the
paper above, the water model is composed of 5 sites - Oxygen (O),
Hydrogen (2x H), center of mass (G), and an auxiliary site (M) along
the bisector a distance 0.021 nm from the O. All the nonbonding terms
are given, and the intramolecular potential is as follows:</p>
<p dir="ltr" style="text-align: left;">U_intra_i = Sum(n=2,6: K_n*(r1^n + r2^n)) + 1/2*K_theta*r3^2 + K_r_theta*r3*(r1+r2) + K_r_r'*r1*r2</p>
<p dir="ltr" style="text-align: left;">Where
r1 and r2 are the O-H bond displacements from the equilibrium value,
and r3 is the H-H displacement. The first term is a sum with force
constants K_2, K_3,...,K_6 multiplied by the corresponding bond
displacement exponentials, and is where I planned to use a tabulated
bond potential a la the '-tableb' option for mdrun. The higher-order
terms have been found necessary for reproducing the OH vibrational
red-shifting in calculated IR, Raman, and SFG spectra. The 2nd term is
the standard harmonic term but for the H-H bond. The 3rd and 4th terms
are bond-angle and bond-bond cross terms, and are readily available in
the gromacs suite already.</p>
<p dir="ltr" style="text-align: left;">Turning
to the tabulated potentials, this is where my questions begin. I've
calculated the 3 columns needed in a table_b0.xvg type of file: x,
f(x), f'(x) -- corresponding to the bondlength, potential, and first
derivative of the potential (a.k.a. the force), respectively. In doing
so I've assumed units of nm, kJ/mol/nm, and kJ/mol/nm^2, respectively.
The function, f(x), is as follows:</p>
<p dir="ltr" style="text-align: left;">f(x) = K_2*x^2 + K_3*x^3 + K_4*x^4 + K_5*x^5 + K_6*x^6</p>
<p dir="ltr" style="text-align: left;">It's rather straightforward, and gives a clean analytical solution for the first derivative:</p>
<p dir="ltr" style="text-align: left;">f'(x) = 2*K_2*x + 3*K_3*x^2 + 4*K_4*x^3 + 5*K_5*x^4 + 6*K_6*x^5</p>
<p dir="ltr" style="text-align: left;">Question:
Is the analytical version what I sohuld put in the 3rd column, or
should I use a finite difference method? Using the analytical version
produced a warning that my potential derivative values are an average
of 164% off from the internally calculated ones - does that message
mean anything? Does it mean that I'm doing things right? Furthermore,
why bother with this 3rd column if it is something already calculated
(and compared to) by gromacs?<br>
</p>
<p dir="ltr" style="text-align: left;">Using a
tabspace (resolution) of 1000 datapoints per nm, I constructed a
3-column table, table_b0.xvg, and added the following to the .itp
file's [ bond ] section:</p>
<p>[ bonds ]<br>; pair-potential harmonic bonds O-H terms<br>1 2 8 0 1.<br>1 3 8 0 1.<br></p>
<p dir="ltr" style="text-align: left;">(where 1 == Oxygen; 2,3 == hydrogens)</p>
<p dir="ltr" style="text-align: left;">This
signifies that a lookup table will be used (with exclusions) and that
the potential will be multiplied by a force constant of 1.0 kJ/nm/mol
-- effectively not changing the tabulated value.</p>
<p dir="ltr" style="text-align: left;">Question: Does this all sound
sane so far? Is that value of 1.0 for the force constant actually
keeping the potential term unaltered? Furthermore, what changes if I
use the function type 8(w/exclusions) vs. 9 (w/o exclusions)?? What is
being excluded if I use 8, and what does it mean if something is
excluded from bonding terms? Where are the intra-molecular bonding
exclusions defined?</p>
<p dir="ltr" style="text-align: left;">Next
step is to run the simulation. All is similar to previous successful
simulations in terms of my mdp input file, but I add the option of
'mdrun -tableb table.xvg' in order to use the tabulated bond potentials
for my O-H bonds. </p>
<p dir="ltr" style="text-align: left;">Question:
is -tableb the only parameter needed to pass information about the
tabulated data for mdrun and simulations? Are any options necessary in
the .mdp run input file?</p>
<p dir="ltr" style="text-align: left;">Once
I get some confirmation about the above I'll feel confident enough to
boldly stride towards finessing the water model .itp file for sharing,
and more coherent documentation. I'm sure others have come across this
issue before, and it would be great to have some type of example fully
illustrating use of the tabulated potential data in a simulation.</p>
<p dir="ltr" style="text-align: left;"><br></p><font color="#888888">~Eric
</font>