[gmx-users] Nucleic Acid in vacuo
matthias.ernst2 at student.kit.edu
Tue Sep 13 17:07:16 CEST 2011
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,
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 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
More information about the gromacs.org_gmx-users