[gmx-users] Problems with NVT simulation with CHARMM27
sa
sagmx.mail at gmail.com
Thu Oct 21 20:00:12 CEST 2010
Dear All,
I am trying to equilibrate in NVT ensemble of a peptide with glycolipid (GL)
molecules in a cubic box filled with TIP3P water (from Mackerell et al.).
The force field for GL were converted to GROMACS format from CHARMM27 force
field with additional parameters non present in the forcefield added in
lipids.rtp, ffbonded.itp and ffnonboded.itp files.
I use the git version of GMX4.5.1 downloaded yesterday.
I can minimize successfully the system with the following em.mdp before to
perform the nvt run
------- em.mdp
title = Glyco + peptide in water
; Preprocessor - specify a full path if necessary.
cpp = cpp
include = -I../top
define = -DFLEXIBLE
integrator = steep
;nstcgsteep = 10000
emstep = 0.01
emtol = 400.0
;dt = 0.002
pbc = xyz
nsteps = 10000
nstlist = 1
ns_type = grid
vdw-type = Cut-off ; twin range cut-off’s with neighbor
list
rlistlong = 1.0
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.2
;fourierspacing = 0.12
;pme_order = 4
;ewald_rtol = 1e-05
;optimize_fft = yes
I obtained the final energies :
Energies (kJ/mol)
Bond Angle U-B Proper Dih. Improper Dih.
3.65918e+04 1.86096e+04 6.02333e+03 3.68911e+04 8.61501e+00
CMAP Dih. LJ-14 Coulomb-14 LJ (SR) LJ (LR)
-1.31952e+02 1.17641e+04 1.76586e+05 1.84759e+05 -1.92541e+03
Coulomb (SR) Coul. recip. Potential Pressure (bar)
-1.26284e+06 -1.61597e+05 -9.55265e+05 0.00000e+00
Steepest Descents converged to Fmax < 400 in 518 steps
Potential Energy = -9.5526519e+05
Maximum force = 3.7991144e+02 on atom 313
Norm of force = 9.2632399e+00
When I try to do the NVT equilibration step at 300 K with the following
nvt.mdp, the simulation crashed quickly (app. after 20 ps of run) with
several LINCS warnings for TIP3 water atoms
------- nvt.mdp
title = Glyco + peptide in water
define = -DPOSRES ; position restrain for the peptide
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 24000 ; 2 * 50000 = 100 ps
; to test
dt = 0.001 ; 2 fs
; Output control
nstxout = 1000 ; save coordinates every 0.2 ps
nstvout = 10000 ; save velocities every 0.2 ps
nstenergy = 10000 ; save energies every 0.2 ps
nstlog = 100 ; update log file every 0.2 ps
energygrps = Protein Non-Protein
; Bond parameters
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 cels
nstlist = 5 ; 10 fs
vdw-type = Cut-off ; twin range cut-off’s with neighbor list
rlistlong = 1.0 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; 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.12 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; one 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
----- Here the message obtained
step 1499: Water molecule starting at atom 47947 can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
Step 1500:
The charge group starting at atom 36 moved than the distance allowed by the
domain decomposition (1.200000) in direction Y
distance out of cell -14795603.000000
Old coordinates: 5.159 5.115 3.453
New coordinates: 8325234.000 -14795599.000 3.894
Old cell boundaries in direction Y: 4.371 9.400
New cell boundaries in direction Y: 4.370 9.400
Step 1500:
The charge group starting at atom 47947 moved than the distance allowed by
the domain decomposition (1.200000) in direction Y
distance out of cell 19727472.000000
Old coordinates: 5.530 4.419 2.912
New coordinates: -11100299.000 19727476.000 5.506
Old cell boundaries in direction Y: 0.000 4.579
New cell boundaries in direction Y: 0.000 4.585
---------
I aware that these errors are probably caused by a system not well
equilibrated but since the Steepest Descents minimization are well
converged, i am puzzled for the NVT simulation. I have checked the
parameters of the force field and for me they are corrects, tried several
protocols (reduce the time step (2 ->1 fs, above, lower the temperature, or
increase the number of minimization steps, etc.) with no success.
Any advices will be greatly appreciated.
Thank you for your help.
sa
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20101021/ec831d3b/attachment.html>
More information about the gromacs.org_gmx-users
mailing list