[gmx-users] Energy Minimisation and Equilibration Problems
Bruce D. Ray
brucedray at yahoo.com
Fri Aug 7 00:48:55 CEST 2009
On Thursday, August 6, 2009 at 4:08:18 PM, Nancy <nancy5villa at gmail.com> wrote:
> I have attempted to perform energy minimisation from scratch again, and these are the
> commands that I am using to do so:
>
> $ .../topolbuild1_2_1/src/topolbuild -n ethanediol -dir .../topolbuild1_2_1/dat/gromacs -ff gmx53a6
>
> which generates the files:
>
> ethanediol.gro
> ethanediol.log
> ethanediolMOL.mol2
> ethanediol.top
> ffethanediol.itp
> posreethanediol.itp
>
> in the "ethanediol.log" file I noticed the following lines with asterisks:
>
> ========================================
> Angles Force Field Results
> Angle Atoms force angle method measured
> 1 H12- C1- H11 ****** ****** ****** 109.403
> 2 C2- C1- H11 ****** ****** ****** 109.456
>
3 O1- C1- H11 ****** ****** ****** 109.438
> 4 C2- C1- H12 ****** ****** ****** 109.537
> 5 O1- C1- H12 ****** ****** ****** 109.468
> 6 O1- C1- C2 320.0 109.50 1 109.526
>
7 C1- C2- O2 320.0 109.50 1 109.526
> 8 C1- C2- H22 ****** ****** ****** 109.453
> 9 C1- C2- H21 ****** ****** ****** 109.537
> 10 H22- C2- O2 ****** ****** ****** 109.439
>
11 H21- C2- O2 ****** ****** ****** 109.468
> 12 C2- O2- HO2 450.0 109.50 1 106.864
> 13 H21- C2- H22 ****** ****** ****** 109.404
> 14 C1- O1- HO1 450.0 109.50 1 106.864
>
========================================
These are fine. Only the nonpolar hydrogens lack parameters and they are removed from
consideration in the final topology because this is a united atoms model force field.
> I then enlarge the box and solvate the molecule:
>
> $ editconf -f ethanediol.gro -o ethanediol_box.gro -box 3 3 3 -c
>
> $ genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -p ethanediol.top -box 3 3 3
>
> There are 2679 atoms (891 waters and 1 ethanediol). I then configure the minimisation with grompp:
>
> $ grompp -f minim.mdp -c ethanediol_solv.gro -p ethanediol.top -o em.tpr
>
> using the following .mdp file:
>
> ========================================
> define = -DFLEXIBLE
Just out of curiosity, why are you defining flexible waters at this point?
> integrator = steep
> emtol = 10.0
> emstep = 0.1
> nsteps = 15000
>
> nstlist = 1
> ns_type = grid
>
rlist = 1.0
> coulombtype = PME
> rcoulomb_switch = 1.0
> rvdw_switch = 1.3
> pbc = xyz
>
> pme_order = 4
> constraints = none
>
> nstxout = 1
> nstvout = 1
> nstenergy = 1
>
nstlog = 1
> nstcomm = 1
> Tcoupl = no
> Pcoupl = no
> gen_vel = no
> ========================================
>
> I then run the simulation for 15,000 steps:
>
> $ mdrun -v -deffnm em
>
> the last few lines of output are as follows:
>
> ========================================
> Step=14991, Dmax= 5.3e-06 nm, Epot= -6.24538e+04 Fmax= 2.50258e+04, atom= 2397
> Step=14992, Dmax= 6.4e-06 nm, Epot= -6.24539e+04 Fmax= 3.52982e+04, atom= 1171
>
Step=14993, Dmax= 7.6e-06 nm, Epot= -6.24540e+04 Fmax= 3.71731e+04, atom= 2397
> Step=14994, Dmax= 9.2e-06 nm, Epot= -6.24540e+04 Fmax= 4.93409e+04, atom= 1171
> Step=14995, Dmax= 1.1e-05 nm, Epot= -6.24541e+04 Fmax= 5.53109e+04, atom= 2397
>
Step=14997, Dmax= 6.6e-06 nm, Epot= -6.24546e+04 Fmax= 8.50832e+03, atom= 1171
> Step=14998, Dmax= 7.9e-06 nm, Epot= -6.24547e+04 Fmax= 6.37691e+04, atom= 2397
> Step=14999, Dmax= 9.5e-06 nm, Epot= -6.24553e+04 Fmax= 2.67488e+04, atom= 1171
>
Step=15000, Dmax= 1.1e-05 nm, Epot= -6.24548e+04 Fmax= 8.09451e+04, atom= 2397
> writing lowest energy coordinates.
>
> Steepest Descents did not converge to Fmax < 10 in 15001 steps.
> Potential Energy = -6.2455254e+04
>
Maximum force = 2.6748828e+04 on atom 1171
> Norm of force = 9.5851715e+02
> ========================================
>
> After
viewing the trajectory with ngmx, I noticed that most of the motion of
the water
> molecules occurs within the first 50 picoseconds.
>
> Please advise on how to further minimise this system.
I downloaded NSC93876 ethylene glycol as a Sybyl (.sy2) file from NCI, processed it with
dos2unix, added the missing@<TRIPOS> last line, and named it ethylene_glycol.mol2
Charge information is absent on this file as downloaded, and I did not do anything to
account for charges. This may make my results less than useful. I then ran the following commands:
/usr/bdray/bin/topolbuild1_3 -dir /usr/bdray/tables -ff gmx53a6 -n glycol/ethylene_glycol
editconf -f ethylene_glycol.gro -o ethanediol_box.gro -box 3 3 3 -c
genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -p ethanediol.top -box 3 3 3
grompp -f em.mdp -c ethanediol_solv.gro -p ethanediol.top -o em.tpr
File em.mdp includes:
constraints = none
integrator = steep
nsteps = 20000
;
; Energy minimizing stuff
;
emtol = 10.0
emstep = 0.01
nstlist = 1
coulombtype = PME
nstcomm = 1
ns_type = grid
rlist = 1
rcoulomb = 1.0
rvdw = 1.3
nstxout = 1
pbc = xyz
pme_order = 4
I deliberately did not use the flexible waters model with this energy minimization
I then ran
mdrun -v -deffnm em
I got the result:
Stepsize too small, or no change in energy.
Converged to machine precision,
but not to the requested precision Fmax < 10
Double precision normally gives you higher accuracy.
You might need to increase your constraint accuracy, or turn
off constraints alltogether (set constraints = none in mdp file)
Steepest Descents converged to machine precision in 6291 steps,
but did not reach the requested Fmax < 10.
Potential Energy = -4.9487391e+04
Maximum force = 9.0675652e+01 on atom 5
Norm of force = 2.5599347e+02
I note that atom 5 is part of the ethylene glycol molecule. Furthermore, with very few
exceptions, most of the activity in the energy minimization was devoted to the ethylene
glycol molecule. I presume that I would reach the desired Fmax if I were to use double
precision. However, an Fmax < 90.7 would seem to me to be usable.
Perhaps the problem stems from the use of flexible waters in the energy minimization?
--
Bruce D. Ray, Ph.D.
Associate Scientist
IUPUI
Physics Dept.
402 N. Blackford St.
Indianapolis, IN 46202-3273
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090806/9af27d21/attachment.html>
More information about the gromacs.org_gmx-users
mailing list