[gmx-users] Nucleic Acid in vacuo

[Justin A. Lemkul]

Matthias Ernst wrote:
> Dear Gromacs users,
> 
> I am currently working on porting forcefields from Gromacs to another
> program with the aim of simulating DNA. Therefore, I wrote a little
> script which reads in a PDB file, evaluates the Lennard-Jones and the
> Coulomb energy in vacuo with the forcefield data from Gromacs.  I tried
> to implement the search for the neighbours by comparing inter-atomic
> distances to the van-der-Waals-radii and a flexible list-based factor
> weighting of the various contributions (where -1 is a default value, 0
> is the atom itsself, 1 a bond to a next neighbour and so on, so my
> factor "3" should be equivalent to the "1-4-interaction" in Gromacs).
> 
> To verify what I have done, I tried to do a calculation of a single "DA"
> nucleic acid with Gromacs to compare the LJ and Coulomb energy.
> Unfortunately, the results do not match and I do not understand why.
> 
> Please find the PDB code I used here http://pastebin.com/zB3sQHdZ and
> the python code I used for evaluating the potentials here
> http://pastebin.com/2v4nVdg3
> 
> The results I got with my code are 
> -> Total Lennard-Jones Energy [kJ/mol]:   0.712752493372
> -> Total 1-4-Lennard-Jones Energy [kJ/mol]:   37.1121824156
> -> Total Coulomb Energy [kJ/mol]:   -4354.55160305
> -> Total Coulomb 1-4 Energy [kJ/mol]:   -207.712332451
> 
> Using Gromacs, I got
>    Energies (kJ/mol)
> Angle    Proper Dih.  Improper Dih.          LJ-14     Coulomb-14
> 2.00975e+01    6.05485e+01    4.28869e-04    3.69509e+01   -2.06010e+02
> LJ (SR)   Coulomb (SR)      Potential Pressure (bar)   Constr. rmsd
> -1.45022e+00   -5.67518e+01   -1.46615e+02    0.00000e+00    1.80690e-04
> 
> The option in my mdp are
> title           =       INITIAL_EM0
> cpp             =       /lib/cpp
> constraints     =       all-bonds
> integrator      =       steep
> nsteps          =       50
> nstcomm         =       1
> ns_type         =       grid
> rlist           =       0.0
> rcoulomb        =       0.0
> rvdw            =       0.0
> 
> ;
> ; Energy minimizing stuff
> ;
> emtol           =       1500.0
> emstep          =       0.01
> 
> Tcoupl          =       no
> Pcoupl          =       no
> gen_vel         =       no
> pbc = no
> 
> I have not yet done a lot of calculations with Gromacs so maybe I did
> something wrong. If so, please correct me.
> 

You're better off with a single-point energy calculation, rather than a 50-step 
EM.  Use the MD integrator with nsteps = 0 to compare configurations.  If I 
recall, steepest descents makes a move before step zero, so you're not comparing 
the same quantities.

> What I wonder is that the 1-4-energies seem quite reasonable (LJ: 37.11
> vs. 36.95, Coulomb -207.7 vs. -206.0) but the rest is rather different,
> especially the sign of the LJ term and the order of magnitude of the
> Coulomb term. Why is that? What does the "SR" in LJ and Coulomb in the
> Gromacs output mean?
> 

SR means short-range, i.e. within the cutoff.  Since you've set all cutoffs to 
zero, this effectively means all interactions are calculated in this case.

I can't offer an explanation for the different values observed, but I'll suggest 
you look into the Gromacs code to see how it calculates the energies and compare 
your program to see if you are, in fact, doing things the same way.

-Justin

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

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

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