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.abra...@anu.edu.au> >Subject: Re: [gmx-users] g_energy inconsistent results >To: Discussion list for GROMACS users <gmx-users@gromacs.org> >Message-ID: <4d7622f8.7050...@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. >> -- gmx-users mailing list gmx-users@gromacs.org 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 gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists