[gmx-users] Nucleic Acid in vacuo

[Matthias Ernst]

Thank you very much for your help, Justin. 
I tried it with those mdp options:

title           =       INITIAL_EM0
cpp             =       /lib/cpp
constraints     =       all-bonds
integrator      =       md
rlist           =       0.0
rcoulomb        =       0.0
rvdw            =       0.0
pbc 		= 	no

This resulted in
Angle    Proper Dih.  Improper Dih.          LJ-14     Coulomb-14
1.88961e+01    6.04179e+01    5.60480e-04    3.74069e+01   -2.06394e+02
LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
6.94305e-03   -5.66241e+01   -1.46290e+02    2.88457e-01   -1.46002e+02
Temperature Pressure (bar)   Constr. rmsd
1.21730e+00    0.00000e+00    0.00000e+00

So, now the 1-4-values are slightly better (37.11 <-> 37.04 LJ, -207.7
<-> -206.3 Coulomb), while the SR-LJ and Coulomb energys are at two
orders of magnitude too small (LJ: 0.7127 <-> 0.0069430, Coulomb:
-4354.5 <-> -56.62). 

In fact, I considered looking through the code but wanted to avoid it
because I fear that there are so much optimisations implemented in
Gromacs that I probably will not understand whats going on and why.
Maybe I will try it, but first I will play with the forcefield
parameters to see what happens if I switch everything off except certain
pairs, then I can compare them.

Besides, I'm wondering, is there a way to disable all optimisations and
just do a "plain" Lennard-Jones and Coulomb evaluation in vacuo without
any box, solvent, boundary conditions and without cutoffs in Gromacs?

Thanks and best regards,
Matthias


Am Montag, den 12.09.2011, 14:37 +0200 schrieb 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
> 
> ========================================