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

[Berk Hess]

----------------------------------------
> Date: Thu, 29 Nov 2012 08:58:50 -0500
> From: jalemkul at vt.edu
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Electrostatic field in atoms - GROMACS is not reading table interaction functions for bonded interactions
>
>
>
> 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?

No, the issue is that you are using function type 1. Tabulated bonds are function type 8 or 9,
that is what should be in the 3rd column.
The actual number for the table file should appear after that in the 4th column.

Cheers,

Berk

>
> -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
>
> ========================================
> --
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists