<p dir="ltr" style="text-align: left;">Over the past couple days I&#39;ve
attempted to implement the water model developed by Ishiyama &amp;
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&#39;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&#39;ll most likely post follow-ups with more questions to round
out the discussion.<br>
   <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&#39;*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 &#39;-tableb&#39; 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&#39;ve
calculated the 3 columns needed in a table_b0.xvg type of file: x,
f(x), f&#39;(x) -- corresponding to the bondlength, potential, and first
derivative of the potential (a.k.a. the force), respectively. In doing
so I&#39;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&#39;s rather straightforward, and gives a clean analytical solution for the first derivative:</p>
   <p dir="ltr" style="text-align: left;">f&#39;(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&#39;m doing things right? Furthermore,
why bother with this 3rd column if it is something already calculated
(and compared to) by gromacs?<br>
   <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&#39;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
&#39;mdrun -tableb table.xvg&#39; 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&#39;ll feel confident enough to
boldly stride towards finessing the water model .itp file for sharing,
and more coherent documentation. I&#39;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