[gmx-users] Not able to continue with Equilibration

chris.neale at utoronto.ca chris.neale at utoronto.ca
Thu Mar 29 17:19:04 CEST 2012


I didn't follow this whole thread, but I sometimes need to turn off  
all constraints when doing minimization. In fact, for that reason I  
entirely stopped ever using restraints during energy minimization. In  
extreem cases, I have had success also by forcing the water to be  
flexible with a -DFLEXIBLE (or whatever is appropriate for your water  
model; check the .itp) statement in the .mdp file.

Once EM is done, use rigid water and restraints and everything works out ok.

Chris.

-- original message --

Hallo Felix,
thank you for your answer. I tried the constraints = h-bonds but no
change in the output. If I look at the step.pdb file that is produced
after the running I have some strange outcome. For example some of my
atoms are not recognized as part of my protein any more and my
structure is destroied. Am I using some strange parameters for nvt
without realizing it? I've started from the mdp file provided by the
lysozyme tutorial non the Gromacs webpage.
if anyone has any ideas it is more than welcome.
Francesca

integrator      = md            ; leap-frog integrator
nsteps          = 10000         ; 2 * 10000 = 20 ps
dt              = 0.002         ; 2 fs
; Output control
nstxout         = 1             ; save coordinates every 0.2 ps
nstvout         = 1             ; save velocities every 0.2 ps
nstenergy       = 50            ; save energies every 0.2 ps
nstlog          = 50            ; update log file every 0.2 ps
; Bond parameters
unconstrained_start = no
;continuation   = no            ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints     = h-bonds       ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
ns_type         = grid          ; search neighboring grid cells
nstlist         = 5             ; 10 fs
rlist           = 3.0           ; short-range neighborlist cutoff (in nm)
rcoulomb        = 3.0           ; short-range electrostatic cutoff (in nm)
rvdw            = 3.0           ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = PME           ; Particle Mesh Ewald for long-range
electrostatics
pme_order       = 4             ; cubic interpolation
fourierspacing  = 0.16          ; grid spacing for FFT
; Temperature coupling is on
tcoupl          = berendsen     ; modified Berendsen thermostat
tc-grps         = Protein Non-Protein   ; two coupling groups - more accurate
tau_t           = 0.1 0.1       ; time constant, in ps
ref_t           = 300 300               ; reference temperature, one
for each group, in K
; Pressure coupling is off
pcoupl          = no            ; no pressure coupling in NVT
; Periodic boundary conditions
pbc             = xyz           ; 3-D PBC
; Dispersion correction
DispCorr        = EnerPres      ; account for cut-off vdW scheme
; Velocity generation
gen_vel         = yes           ; assign velocities from Maxwell distribution
gen_temp        = 300           ; temperature for Maxwell distribution
gen_seed        = -1            ; generate a random seed



2012/3/29 Rausch, Felix <frausch at ipb-halle.de>:

[Hide Quoted Text]
Hello again,

Well, I once had problems with simulations crashing randomly during  
production runs (sometimes after tens of nanoseconds) with the LINCS  
warnings you described. Switching LINCS from "all-bonds" to only  
"h-bonds" did the trick for me, although I never exactly figured out  
why.
Maybe it's worth a try for you, too?
Cheers,
Felix

-----Ursprüngliche Nachricht-----
Von: gmx-users-bounces at gromacs.org  
[mailto:gmx-users-bounces at gromacs.org] Im Auftrag von francesca vitalini
Gesendet: Donnerstag, 29. März 2012 15:15
An: Discussion list for GROMACS users
Betreff: Re: [gmx-users] Not able to continue with Equilibration

Hi!
I'm having a similar problem. I have a dimer solvated in a big box of  
water plus ions that I have managed to minimize correctly (see output  
of minimization at the end) but when I try to run NVT equilibration  
(see later) I get LINCS warnings(see below) refearred to atoms which  
are not in a cluster or in a strange position. I have added dihedral  
restraints on them but still the same type of error. I'm using GROMACS  
3.3.1. I have tried to switch to a newer version of GROMACS but still  
the same error.
Does anyone have any suggestions?
Thanks
Francesca


MINIMIZATION OUTPUT

Steepest Descents converged to Fmax < 1000 in 681 steps Potential  
Energy  = -1.9597512e+07
Maximum force     =  2.4159846e+02 on atom 13087
Norm of force     =  2.1520395e+04


MINIMIZATION.MDP

define         = -DEFLEXIBLE ; flexible water
integrator      = steep         ; Algorithm (steep = steepest descent
minimization)
emtol           = 1000.0         ; Stop minimization when the maximum
force < 1000.0 kJ/mol/nm
emstep          = 0.001          ; Energy step size
nsteps          = 50000         ; Maximum number of (minimization)
steps to perform

; Parameters describing how to find the neighbors of each atom and how  
to calculate the interactions
nstlist         = 1             ; Frequency to update the neighbor
list and long range forces
ns_type         = grid          ; Method to determine neighbor list
(simple, grid)
rlist           = 1.0           ; Cut-off for making neighbor list
(short range forces)
coulombtype     = PME           ; Treatment of long range
electrostatic interactions
rcoulomb        = 1.0           ; Short-range electrostatic cut-off
rvdw            = 1.0           ; Short-range Van der Waals cut-off
pbc             = xyz           ; Periodic Boundary Conditions (yes/no)


NVT.MDP

define          = -DPOSRES      ; position restrain the protein
; Run parameters
integrator      = md            ; leap-frog integrator
nsteps          = 10000         ; 2 * 10000 = 20 ps
dt              = 0.002         ; 2 fs
; Output control
nstxout         = 50            ; save coordinates every 0.2 ps
nstvout         = 50            ; save velocities every 0.2 ps
nstenergy       = 50            ; save energies every 0.2 ps
nstlog          = 50            ; update log file every 0.2 ps
; Bond parameters
unconstrained_start = no
;continuation   = no            ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints     = all-bonds     ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
ns_type         = grid          ; search neighboring grid cells
nstlist         = 5             ; 10 fs
rlist           = 3.0           ; short-range neighborlist cutoff (in nm)
rcoulomb        = 3.0           ; short-range electrostatic cutoff (in nm)
rvdw            = 3.0           ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = PME           ; Particle Mesh Ewald for long-range
electrostatics
pme_order       = 4             ; cubic interpolation
fourierspacing  = 0.16          ; grid spacing for FFT
; Temperature coupling is on
tcoupl          = berendsen     ; modified Berendsen thermostat
tc-grps         = Protein Non-Protein   ; two coupling groups - more accurate
tau_t           = 0.1 0.1       ; time constant, in ps
ref_t           = 300 300               ; reference temperature, one
for each group, in K
; Pressure coupling is off
pcoupl          = no            ; no pressure coupling in NVT
; Periodic boundary conditions
pbc             = xyz           ; 3-D PBC
; Dispersion correction
DispCorr        = EnerPres      ; account for cut-off vdW scheme
; Velocity generation
gen_vel         = yes           ; assign velocities from Maxwell distribution
gen_temp        = 300           ; temperature for Maxwell distribution
gen_seed        = -1            ; generate a random seed

; Dihedral restraints

dihre               =  simple ; Some dihedrals are restrained  for
instance peptide bonds are set to trans conformation.
dihre_fc            =  500
dihre_tau           =  0.0
nstdihreout         =  50

NVT LINCS ERROR

Step 0, time 0 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 0.002851 (between atoms 10721 and 10723) rms 0.000161 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    672    673   33.7    0.1000   0.1000      0.1000
   1081   1082   33.8    0.1000   0.1000      0.1000
   1647   1648   34.4    0.1000   0.1000      0.1000
   2819   2820   37.1    0.1000   0.1000      0.1000
   2920   2921   39.6    0.1000   0.1000      0.1000
   4090   4091   35.2    0.1000   0.1000      0.1000
   4727   4728   37.2    0.1000   0.1000      0.1000
   5160   5161   31.7    0.1000   0.1000      0.1000
   6824   6825   33.3    0.1000   0.1000      0.1000
...
step 0
Step 1, time 0.002 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 0.017976 (between atoms 14692 and 14693) rms 0.000520 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    672    673   67.8    0.1000   0.0999      0.1000
....

Step 2, time 0.004 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 0.060681 (between atoms 6906 and 6907) rms 0.001874 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    672    673   80.8    0.0999   0.0999      0.1000
   1081   1082   90.0    0.1000   0.1008      0.1000
   ...

Step 3, time 0.006 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 0.104238 (between atoms 10603 and 10604) rms 0.003281 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    672    673   67.1    0.0999   0.1000      0.1000
   1081   1082   69.6    0.1008   0.1000      0.1000
...

Step 4, time 0.008 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 0.569872 (between atoms 8709 and 8710) rms 0.013041 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    672    673   40.3    0.1000   0.1000      0.1000
   1081   1082   39.6    0.1000   0.1000      0.1000
   1249   1250   37.6    0.1000   0.1000      0.1000
...

Step 5, time 0.01 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 2.761976 (between atoms 8202 and 8204) rms 0.056483 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
   3320   3321   43.0    0.1090   0.1090      0.1090
   3875   3876   32.0    0.1090   0.1090      0.1090
   4013   4014   34.3    0.1090   0.1082      0.1090
   4070   4071   42.1    0.1092   0.1092      0.1090
   4072   4073   45.2    0.1090   0.1086      0.1090
   4727   4729   34.2    0.1007   0.1002      0.1000
   4961   4962   31.4    0.1091   0.1089      0.1090
   4992   4993   46.0    0.1091   0.1089      0.1090
   8159   8160   30.9    0.1251   0.1465      0.1230
  ...

Step 6, time 0.012 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 22.009092 (between atoms 8202 and 8204) rms 0.488487 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
   2920   2921   44.7    0.1000   0.1001      0.1000
   3320   3321   36.6    0.1090   0.1087      0.1090
   3592   3593   32.1    0.1090   0.1089      0.1090
  ....
Warning: 1-4 interaction between 8159 and 8168 at distance 1.580 which  
is larger than the 1-4 table size 1.000 nm These are ignored for the  
rest of the simulation This usually means your system is exploding, if  
not, you should increase table-extension in your mdp file

Step 7, time 0.014 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 403.393890 (between atoms 8202 and 8204) rms 8.430617 bonds that  
rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    302    303   31.7    0.1003   0.1001      0.1000
    672    673   66.4    0.1001   0.1001      0.1000
   1050   1051   33.1    0.1003   0.1003      0.1000
   1081   1082   46.7    0.1000   0.1000      0.1000
   1225   1226   35.0    0.1003   0.1003      0.1000
   ....

Wrote pdb files with previous and current coordinates

Step 8, time 0.016 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 7568.298828 (between atoms 8202 and 8204) rms 155.435715 bonds  
that rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length
    205    206   31.6    0.1003   0.1002      0.1000
    302    303   36.9    0.1001   0.1002      0.1000
    672    673   76.0    0.1001   0.0997      0.1000
   1050   1051   34.3    0.1003   0.1003      0.1000
   1081   1082   63.0    0.1000   0.0999      0.1000
   1225   1226   36.6    0.1003   0.1003      0.1000
   1249   1250   38.7    0.1001   0.1001      0.1000
   1647   1648   62.8    0.1000   0.0996      0.1000
   1679   1680   42.5    0.1000   0.1000      0.1000
   2214   2215   31.6    0.1002   0.1001      0.1000
   2577   2578   31.1    0.1003   0.1002      0.1000
   2819   2820   68.0    0.1000   0.0997      0.1000
   2920   2921   71.5    0.1023   0.1000      0.1000
   3468   3469   34.7    0.1092   0.1090      0.1090
   3676   3677   34.7    0.1090   0.1090      0.1090
   3714   3717   34.7    0.0998   0.0999      0.1000
   3746   3747   70.8    0.1086   0.1069      0.1090
....

Back Off! I just backed up step7.pdb to ./#step7.pdb.1# Wrote pdb  
files with previous and current coordinates

Step 9, time 0.018 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 148599.687500 (between atoms 8202 and 8204) rms 3002.328125 bonds  
that rotated more than 30 degrees:
  atom 1 atom 2  angle  previous, current, constraint length ....

Back Off! I just backed up step8.pdb to ./#step8.pdb.1# Wrote pdb  
files with previous and current coordinates Segmentation fault


2012/3/29 Mark Abraham <Mark.Abraham at anu.edu.au>:
On 29/03/2012 8:22 PM, Hendry wrote:

Hi,



I am using Gromacs 4.5.4. After successful minimization by SD, I
continued with equilibration step but I got the below errors. I tried
many times with different parameters but the problem still persists. I
have given errors and md parameters of equilibration step below. I
have also provided my minimization output at the end. Could you
provide some suggestions what went wrong?.


You are using PR-pressure coupling for equilibration, which is an
unstable combination. You are coupling ions to their own thermostat,
which is unstable. Do check out the GROMACS manual for discussion of
how these algorithms work, and also the website for some practical  
observations.

Mark

-





More information about the gromacs.org_gmx-users mailing list