[gmx-users] RE: g_energy inconsistent results

[Ehud Schreiber]

Dear Mark Abraham (and anyone else interested),

I have implemented your suggestion, changing in the status.mdp file to
integrator = md and adding continuation = yes (this is the new name of
the unconstrained_start parameter); however, this did not have any
effect - the energy remained as before, different from the one obtained
during the minimization.

I have then encountered another (perhaps related) issue.
I thought that the problem may arise from the combination of the
generalized Born and all-vs-all settings.
I have therefore changed in the minimization to rgbradii = rlist =
rcoulomb = rvdw = 100 (from = 0).
As my system is far smaller than 100 nm, I expected these parameters to
provide also an all-vs-all setting (even if non-optimized).
Nevertheless, the resulting structure differed from the previous
minimized one (rmsd ~ 0.005 nm, delta energy ~ a few kJ/mol). Can this
arise from rounding effects only or does this signify a problem? I'm not
sure what rgbradii does, but the manual states that it must equal rlist.
In addition, changing also status.mdp to have all radii = 100 and
running on the new optimized output again does not yield the same energy
as during the optimization. 

Does all of this make any sense to you?

Thank,
Ehud.


>Message: 6
>Date: Tue, 08 Mar 2011 23:37:12 +1100
>From: Mark Abraham <Mark.Abraham at anu.edu.au>
>Subject: Re: [gmx-users] g_energy inconsistent results
>To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>Message-ID: <4D7622F8.7050407 at anu.edu.au>
>Content-Type: text/plain; charset="iso-8859-1"
>
>On 8/03/2011 9:44 PM, Ehud Schreiber wrote:
>>
>> 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:

> Zero-step EM does not calculate the initial energy because it is not
useful for gradient-based energy minimization.
> I don't recall the details, but perhaps the first EM step is reported
as step zero.
> Instead, you should use zero-step MD (with unconstrained_start = yes),
or (for multiple single points) mdrun -rerun.

> You will not necessarily reproduce the g_energy energies with
anything, because the energy is dependent on the state of the neighbour
lists.
> If nstenergy is a multiple of nstlist, then those energies should be
fairly reproducible.

> I have updated the grompp source to provide a note to the user to warn
against zero-step EM.

> Mark

>> ;;;;;;;;;;;;;;;;;;; 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.
>>