[gmx-users] Electrostatic field in atoms - GROMACS is not reading table interaction functions for bonded interactions

[Justin Lemkul]

On 11/29/12 8:18 AM, Ariana Torres Knoop wrote:
> Dear GROMACS users,
>
> We are trying to determine the electric fields at the C and N atoms in a protein. To do so, we thought of performing a one step energy minimization with the Lennard-Jones as well as the bonded interactions turned off, in order to get only the coulomb forces.
> For that, according to www.gromacs.org/Documentation/How-tos/Tabulated_Potentials , we made a table (table.xvg) in which the Lennard-Jones potential was set to zero and we also made three tables (table_b1.xvg, table_a5.xvg and table_d5.xvg) according to section 4.2.13 of the Manual (Tabulated interaction functions), in which the f(x)=0 for the whole range of distances, angles, dihedrals present in the simulation box. The numbers 1, 5, 2 correspond to the function number referred to in the topology file. For example for the bonds interaction, in the topology file we have:
> [ bonds ]
> ;  ai    aj funct            c0            c1            c2            c3
>      1     2     1
>      1     3     1
>      1     4     1
>      1     5     1
>      5     6     1
> .
> .
> .
> so we named the table "table_b1.xvg".
>
> We changed the run.mdp file as follows:
> vdw-type                 = User
> coulombtype              = PME-User
>
> And then we try to run it  using:
> mdrun  -v  -deffnm  run.tpr
>
> When checking the energy contribution using g_energy, Lennard-Jones contribution was indeed zero, but
> for the bonded interactions this is not  the case:
>
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> U-B                         18180.1         --    276.416   -540.638  (kJ/mol)
> Proper Dih.                 11732.3         --    27.6829   -54.0396  (kJ/mol)
> Improper Dih.               1027.26         --    47.1598   -91.7097  (kJ/mol)
> LJ-14                             0         --          0          0  (kJ/mol)
> LJ (SR)                           0         --          0          0  (kJ/mol)
>
> We though that maybe the tables for the bonded interactions were not taken, so we tried:
>   mdrun  -v -tableb table_b1.xvg  table_a5.xvg table_d2.xvg  -deffnm  run.tpr
>
> But it also did not work.
>
> We are not sure whether the problem is in the table definition, the way of running the program, the run type, the name of the tables, or if there should be some other modifications in the *mdp file.
>
> We would be very grateful if anyone could help us out here.
>

I think the problem comes from the fact that the function types you are 
specifying in the topology correspond to actual potentials and thus the tables 
are being ignored.  Per Table 5.5, it seems that any value for the table is 
acceptable (n >= 0) but using numbers corresponding to actual function types may 
not work.  The specified format also requires that both the table number and 
force constant (k) be specified in the topology, but the [bonds] section you 
show above has no value of k and the function type corresponds to a simple 
harmonic bond.  There is no "bond" energy in the output, which suggests you are 
using constraints - is that the case?

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================