[gmx-users] how to use the new potential

[Matthew Hoopes]

>He, Yang wrote:
>> Hi all users,
>> 
>> I am engaged in deal with the course grain model,in which the 
>> potential functions are not the same as in gromacs.Hence, If I want to 
>> use the new potentials , I guess I can change the source code or use 
>> the table. But I never have this experience about  that.Can anyone of 
>> you give me some suggestions or examples?
>
>Start with the manual. :-) Hopefully the version 4 manual will be out 
>soon, but even if not, the 3.3 manual will be fine for nonbonded 
>tabulated potentials.
>
>Mark

Here is an example:

#
# Tabulated WCA potentials, dr=0.002, rc=3.5, no coulomb
#
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
2.0000000000e-03  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
...
3.6000000000e-02  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
3.8000000000e-02  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
4.0000000000e-02  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  2.3841857812e+17  7.1525573584e+19
4.2000000000e-02  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  1.3276038530e+17  3.7931538762e+19
4.4000000000e-02  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  7.5967506265e+16  2.0718410875e+19
4.6000000000e-02  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  4.4562136749e+16  1.1624905294e+19
...
3.4980000000e+00  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
3.5000000000e+00  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
0.0000000000e+00  0.0000000000e+00  0.0000000000e+00

Column 1) is r (distance between atoms) 6) h(r) the repulsive term of the
potential and 7) -h'(r) as defined in section 6.7.2

1) Just use your favorite script (I used octave) to generate the following
text file with printf statements or something similar.
2) If some of the early values are too large, it may cause an error so I
just zeroed them out.
3) The information you need is in a few places in the manual. The table
above is defined in section 6.7.2 on page 131 of the GROMACS 4.0 manual.
4) You will also want to read section 7.3.12 on page 148 about entries in
your parameter file (e.g. energy groups, table file naming conventions)
5) You will need an index file that has particle groups for your new
potentials. The group names in the index file are what you use for energy
groups.
6) Pay attention to which combination rules you use in the topology file
7) Note that in the LJC(12-6-1) table in the top directory of GROMACS, the
minimum of the LJ is -0.25 and not -1.0 so I guess there is a factor of
4*epsilon multiplied to the table values.

What I don't know, and perhaps someone else can answer this is what should
be done if you are unsure of how your potential is separated in to
attractive and repulsive terms (e.g. LJ with COS tail). First, it's
piecewise, and second the COS section is not the sum of two terms so how
should it be broken across g(r) and h(r) or is it always OK to use the g(r)
column for the whole potential?

Hope this helps.

-Matt