I can help with only some of this...<br><br>----- Original Message -----<br>From: Eric Shamay &lt;eshamay@uoregon.edu&gt;<br>Date: Wednesday, October 6, 2010 18:35<br>Subject: [gmx-users] New shell water model, tabulated bonding interactions, and the documentation.<br>To: gmx-users@gromacs.org<br><br><p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>Over the past couple days I'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'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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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. <br></p><p dir="ltr" style="text-align: left;"><br></p><p dir="ltr" style="text-align: left;">Yes, per manual 4.2.13</p><p dir="ltr" style="text-align: left;"><br></p><p dir="ltr" style="text-align: left;">&gt; The function, f(x), is as follows:</p>         <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>It's rather straightforward, and gives a clean analytical solution for the first derivative:</p>        <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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?</p><p dir="ltr" style="text-align: left;"><br></p><p dir="ltr" style="text-align: left;">I'd expect you using analytic is best, and that GROMACS is doing an internal finite difference calculation as a cross-check. I'm not sure why they might disagree...<br></p><p dir="ltr" style="text-align: left;"><br> </p>        <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>[ bonds ]<br><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>; pair-potential harmonic bonds O-H terms<br><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>1 2 8 0 1.<br><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>1 3 8 0 1.<br></p>    <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>(where 1 == Oxygen; 2,3 == hydrogens)</p>        <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>Question: Does this all sound sane so far? <br></p><p dir="ltr" style="text-align: left;">Yes<br></p><p dir="ltr" style="text-align: left;">&gt; Is that value of 1.0 for the force constant actually keeping the potential term unaltered?&nbsp;</p><p dir="ltr" style="text-align: left;">Unaltered from the table, yes.<br></p><p dir="ltr" style="text-align: left;">&gt;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;">Various nonbonded terms might get excluded, but not in your case. See manual 4.6.1 and 5.4<br></p>         <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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;">I would expect so. This will be straightforward to test. Run one with the old non-table .itp, and another with the table .itp. If you're not convinced something different is happening, adjust k for the tabulated interactions to 100 or something :-)<br></p>        <p dir="ltr" style="text-align: left;"><font style="font-style: normal; font-weight: normal; background-color: rgb(245, 248, 240); font-size: 14px;">&gt; </font>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>That would be good. I've created a stub section for you here http://www.gromacs.org/Documentation/How-tos/Tabulated_Potentials, though you may need to email Rossen to get access to edit.<br><br>Mark