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.
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
