[gmx-users] Very High Forces During Minimisation

[Nancy]

Hello,

I am trying to equilibrate and run MD on a system of solvated ethylene
glycol (ethanediol).  However I am running into numerous problems.

First, I try to minimise the system.  I start with an ethanediol.mol2 file
which contains the structure of the molecule.

$ .../topolbuild1_2_1/src/topolbuild -n ethanediol -dir
.../topolbuild1_2_1/dat/gromacs -ff gmx53a6

the above outputs the following files:

ethanediol.gro
ethanediol.log
ethanediolMOL.mol2
ethanediol.top
ffethanediol.itp
posreethanediol.itp

Note: the ethanediol.log file contains a section that has several 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
============================

$ editconf -f ethanediol.gro -o ethanediol_box.gro -c -d 1.0 -bt cubic

$ genbox -cp ethanediol_box.gro -cs spc216.gro -o ethanediol_solv.gro -p
ethanediol.top

$ grompp  -f minim.mdp -c ethanediol_solv.gro -p ethanediol.top -o em.tpr

the "minim.mdp" file is:

============================
define        = -DFLEXIBLE
integrator    = steep
emtol        = 10.0
emstep          = 0.01
nsteps        = 2000

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
============================

$ mdrun -v -deffnm em

the last few lines of the output of the minimisation are:

============================
Step= 1990, Dmax= 1.1e-05 nm, Epot= -3.38051e+04 Fmax= 4.80304e+04, atom=
559
Step= 1991, Dmax= 1.3e-05 nm, Epot= -3.38053e+04 Fmax= 7.73723e+04, atom=
795
Step= 1992, Dmax= 1.6e-05 nm, Epot= -3.38059e+04 Fmax= 7.31046e+04, atom=
559
Step= 1994, Dmax= 9.5e-06 nm, Epot= -3.38069e+04 Fmax= 2.67253e+04, atom=
739
Step= 1995, Dmax= 1.1e-05 nm, Epot= -3.38076e+04 Fmax= 4.96799e+04, atom=
559
Step= 1996, Dmax= 1.4e-05 nm, Epot= -3.38078e+04 Fmax= 8.04458e+04, atom=
795
Step= 1997, Dmax= 1.6e-05 nm, Epot= -3.38084e+04 Fmax= 7.56146e+04, atom=
559
Step= 1999, Dmax= 9.9e-06 nm, Epot= -3.38094e+04 Fmax= 2.68226e+04, atom=
739
Step= 2000, Dmax= 1.2e-05 nm, Epot= -3.38102e+04 Fmax= 5.49575e+04, atom=
559

writing lowest energy coordinates.

Steepest Descents did not converge to Fmax < 10 in 2001 steps.
Potential Energy  = -3.3810168e+04
Maximum force     =  5.4957539e+04 on atom 559
Norm of force     =  2.7119446e+03
============================

I have already tried increasing the number of steps, but I found no way of
lowering the forces.  I have found that by limiting the number of solvent
molecules to ~5, I can manage to acheive lower energies, however, the amount
of time it takes for a larger number of water molecules increases
exponentially.

So, I am trying to proceed to equilibration:

$ grompp -f nvt.mdp -c em.gro -p ethanediol.top -o nvt.tpr

"nvt.mdp" is:

============================
title        = Ethanediol NVT equilibration
define        = -DFLEXIBLE

integrator    = md
nsteps        = 2000
dt        = 0.002

nstxout        = 10
nstvout        = 10
nstenergy    = 10
nstlog        = 10

continuation    = no
constraint_algorithm = lincs
constraints    = all-angles
lincs_iter    = 1
lincs_order    = 4

ns_type        = grid
nstlist        = 5
rlist        = 1.0
rcoulomb    = 1.0
rvdw        = 1.2

coulombtype    = PME
pme_order    = 4
fourierspacing    = 0.16

tcoupl        = V-rescale
tc-grps        = EDO SOL
tau_t        = 0.1 0.1
ref_t        = 300 300

pcoupl        = no

pbc        = xyz

DispCorr    = EnerPres

gen_vel        = yes
gen_temp    = 300
gen_seed    = -1
============================

Unfortunately, I receive numerous errors such as this one:

============================
Step 3, time 0.006 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.067249, max 2.300835 (between atoms 691 and 692)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
    406    407   33.5    0.1005   0.0956      0.1000
    692    693   91.2    0.1631   0.3470      0.1633
    691    692   88.9    0.0999   0.3301      0.1000
    691    693   86.0    0.0999   0.0881      0.1000
    794    795   35.4    0.1707   0.1587      0.1633
    793    795   59.5    0.1063   0.0995      0.1000
Wrote pdb files with previous and current coordinates
Segmentation fault
============================

I believe that equilibration is failing because of the large forces on the
molecules during minimisation (although I am not sure, as I was able to run
a simulation without proper minimisation before).

If at all possible, please advise on how to reduce the energies and forces
during minimisation.


Thank you.

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