[gmx-users] g_energy inconsistent results

[Ehud Schreiber]

Dear Gromacs users,

 

I am working with version 4.5.3, using the opls-aa forcefield in an
implicit solvent, all-vs-all setting:

 

pdb2gmx -ter -ff oplsaa -water none -f file.pdb

 

I am energy-minimizing structures in 3 stages (steep, cg and l-bfgs).
The last stage is the following:

 

grompp -f em3.mdp -p topol.top -c em2.gro -t em2.trr -o em3.tpr -po
em3.mdout.mdp

mdrun -nice 0 -v -pd -deffnm em3

g_energy -s em3.tpr -f em3.edr -o em3.potential_energy.xvg 

 

where the mdp file is:

 

;;;;;;;;;;;;;;;;;;; em3.mdp ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

integrator       = l-bfgs

nsteps           = 50000

implicit_solvent = GBSA

gb_algorithm     = Still 

sa_algorithm     = Ace-approximation 

pbc              = no

rgbradii         = 0 

ns_type          = simple

nstlist          = 0

rlist            = 0

coulombtype      = cut-off

rcoulomb         = 0

vdwtype          = cut-off

rvdw             = 0

nstcalcenergy    = 1

nstenergy        = 1000

emtol            = 0 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 

The last line in the em3.potential_energy.xvg file should give the
(potential) energy of the minimized structure em3.gro .

 

I wish also to compute the potential energy of .gro files in general,
not necessarily obtained from a simulation. For that, I prepared a .mdp
file for a degenerate energy minimization, having 0 steps, designed just
to give the status of the file:

 

;;;;;;;;;;;;;;;;;;; status.mdp ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

integrator       = l-bfgs

nsteps           = 0

implicit_solvent = GBSA

gb_algorithm     = Still

sa_algorithm     = Ace-approximation

pbc              = no

rgbradii         = 0

ns_type          = simple

nstlist          = 0

rlist            = 0

coulombtype      = cut-off

rcoulomb         = 0

vdwtype          = cut-off

rvdw             = 0

nstcalcenergy    = 1

nstenergy        = 1

emtol            = 0

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 

The only changes from the former .mdp file are in nsteps and nstenergy.

 

However, when I run this potential energy status run on em3.gro itself,

 

grompp -f status.mdp -p topol.top -c em3.gro -o status.tpr -po
status.mdout.mdp

mdrun -nice 0 -v -pd -deffnm status

g_energy -s status.tpr -f status.edr -o status.potential_energy.xvg

 

and look at the (single) energy line in status.potential_energy.xvg I
find that the energy does not agree with the one obtained during
minimization (it's higher by some tens of kJ/mol).

 

What am I doing wrong? How should one reliably find the energy of a
given .gro file?

 

Moreover, when changing in status.mdp to integrator = steep, the results
also change dramatically - why should the algorithm matter if no steps
are performed and the initial structure is explored?

 

Thanks,

Ehud.

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110308/048b4386/attachment.html>