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

[Ariana Torres Knoop]

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.

Thanks a lot in advance and with kind regards,

Ariana