[gmx-users] Error in MD simulation

teklebrh at ualberta.ca teklebrh at ualberta.ca
Sun Jun 27 21:59:46 CEST 2010


Dear Gromacs users,

First I have to mention few points. I read all the LINCS error and  
tried a number of options.

1, Minimized well the system. and the Epot and Max Forces are very resonable.

Steepest Descents converged to Fmax < 100 in 379 steps
Potential Energy  = -1.0573750e+05
Maximum force     =  8.0876366e+01 on atom 1117
Norm of force     =  8.7360611e+00

2, Equlibrate well the system. It works for 100ps

Once I run NPT simulation using the following .mdp parameter set it  
shows the following error.

; The etire lines starting with ';' ARE COMMENTS

title		= NPT equilibration using the berendsen thermostat and  
barostat	; Title of run
cpp             = /usr/bin/cpp ; location of cpp on linux

; The following lines tell the program the standard locations where to  
find certain files

; VARIOUS PREPROCESSING OPTIONS
; Preprocessor information: use cpp syntax.
; e.g.: -I/home/joe/doe -I/home/mary/hoe
include                  =

; e.g.: -DI_Want_Cookies -DMe_Too

; RUN CONTROL PARAMETERS

integrator               = md ; leap-frog integrator
; Start time and timestep in ns
tinit                    = 0
dt                       = 0.002     ; 2fs
nsteps                   = 500000   ; 2*500 = 1ns

; THOSE OPTIONS REMOVE MOTION OF THE ASPHALTENE MODEL COMPUNDS  
RELATIVE TO THE SOLVENT/IONS

; mode for center of mass motion removal
comm-mode                = Linear
; number of steps for center of mass motion removal
nstcomm                  = 1
; group(s) for center of mass motion removal
comm-grps                = System

; LANGEVIN DYNAMICS OPTIONS
; Friction coefficient (amu/ps) and random seed
bd-fric                  = 0
ld-seed                  = 1993

; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol                    = 200
emstep                   = 0.01

; Frequency of steepest descents steps when doing CG
nstcgsteep               = 1000
nbfgscorr                = 10

; OUTPUT CONTROL OPTIONS

; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 2000
nstvout                  = 2000
nstfout                  = 0
; Output frequency for energies to log file and energy file
nstlog                   = 1000
nstenergy                = 1000
; Output frequency and precision for xtc file
nstxtcout                = 1000
xtc-precision            = 1000


; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps                 = BIS  HEP
; Selection of energy groups
energygrps               = BIS  HEP

; NEIGHBORSEARCHING PARAMETERS

; nblist update frequency
nstlist                  = 5
; ns algorithm (simple or grid)
ns-type                  = grid
; Periodic boundary conditions: xyz, no, xy
pbc                      = xyz
; nblist cut-off
rlist                    = 1.2

; OPTIONS FOR ELECTROSTATICS AND VDW

; Method for doing electrostatics
coulombtype              = PME
rcoulomb                 = 1.2


; Method for doing Van der Waals
vdw-type                 = Cut-off
; cut-off lengths
rvdw                     = 1.2

; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = EnerPres ; account for the cut-off vdw scheme

; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.16

; EWALD/PME/PPPM parameters
pme_order                = 4
ewald_rtol               = 1e-05
ewald_geometry           = 3d

; IMPLICIT SOLVENT ALGORITHM
implicit_solvent         = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling is on
tcoupl                   = V-rescale; modified berendsen thermostat
; Groups to couple separately
tc-grps                  = BIS  HEP ; three coupling gropus
; Time constant (ps) and reference temperature (K)
tau-t                    = 0.1  0.1   ; time constant
ref-t                    = 298  298   ; reference temprature
; Pressure coupling is on
Pcoupl                   = berendsen; Pressure coupling is on
Pcoupltype               = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p                    = 3  3
compressibility          = 1.47e-4  1.47e-4  ; isothermal  
compressibility, bar^-1
ref-p                    = 1  1
; Random seed for Andersen thermostat
andersen_seed            = 815131

; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel                  = yes ; assign velocities from Maxwell  
distribution (MXD)
gen-temp                 = 298 ; temperature for MXD
gen-seed                 = 173529    ; used to initialize random  
generator for random velocities

; OPTIONS FOR BONDS
constraints              = all-bonds
; Type of constraint algorithm
constraint-algorithm     = Lincs
; Do not constrain the start configuration
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter               = 2
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle          = 30
; Convert harmonic bonds to morse potentials
morse                    = no

==========

It run the simulation for quite sometime but it gives me an error of  
LINCS some how in the middle.

I try changing some parameter sets like step size,tau-p , tau-t,  
lincs-warnangle but did not overcame the problem.

Step 197685, LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.000000, max 0.000001 (between atoms 1856 and 15982)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>     97     98   35.3    0.1000   0.1000      0.1000

can some one give me an advice or help?

For my other 10 MD set up it works fine but with diffrent solvents.

rob



More information about the gromacs.org_gmx-users mailing list