[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field
Pablo Englebienne
p.englebienne at tue.nl
Tue Oct 27 12:14:23 CET 2009
Thanks for the suggestions, Justin.
I'm still having issues with instabilities and large fluctuations in
unconstrained NPT simulations of the CHCl3 box, so I'll appreciate
comments from the list on my run parameters. I think there is a problem
with the pressure coupling, but I'm not sure what to change to fix it.
I modified the topology, using the H-Cl and Cl-Cl as bonds, and the H-C
and C-Cl bond distances as constraints:
---[chcl3.itp]---
[...]
[ bonds ]
1 3 2 gb_47
1 4 2 gb_47
1 5 2 gb_47
3 4 2 gb_48
3 5 2 gb_48
4 5 2 gb_48
[ constraints ]
1 2 1 0.1100
2 3 1 0.1758
2 4 1 0.1758
2 5 1 0.1758
[...]
---[chcl3.itp]---
This prevented the grompp warning message about too many constraints to
appear.
I first minimized the 216-molecule box with the following mdp:
---[em.mdp]---
integrator = steep
nsteps = 50000
emtol = 100
emstep = 0.01
;
; Electrostatics
;
coulombtype = PME
rlist = 1
rcoulomb = 1.0
ns_type = grid
nstlist = 1
;
; vdW
;
rvdw = 1.0
;
; PBC
;
pbc = xyz
;
; constraints
;
constraint_algorithm = lincs
constraints = none
lincs_iter = 1
lincs_order = 4
---[em.mdp]---
This converged well:
---[em.log]---
Energies (kJ/mol)
G96Bond G96Angle LJ (SR) Coulomb (SR) Coul. recip.
3.05135e+00 4.35449e-01 -5.41548e+03 -2.09196e+01 -2.57100e+01
Potential Pressure (bar) Cons. rmsd ()
-5.45863e+03 1.17031e+03 1.50663e-07
Steepest Descents converged to Fmax < 100 in 297 steps
Potential Energy = -5.45862698325888e+03
Maximum force = 9.95795243267423e+01 on atom 798
Norm of force = 3.22657158607021e+01
---[em.log]---
Then, I heated the system to 300K in an NVT simulation:
---[nvt.mdp]---
;
;NVT equilibration
;
; Run parameters
integrator = md
nsteps = 25000
dt = 0.002
; Output control
nstxout = 100
nstvout = 100
nstenergy = 100
nstlog = 100
; Bond parameters
continuation = no
constraint_algorithm = lincs
constraints = none
lincs_iter = 1
lincs_order = 4
; Neighborsearching
ns_type = grid
nstlist = 5
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; Temperature coupling is on
tcoupl = V-rescale
tau_t = 0.1
tc_grps = SYSTEM
ref_t = 300
; Pressure coupling is off
pcoupl = no
; PBC
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = yes
gen_temp = 300
gen_seed = -1
---[nvt.mdp]---
This simulation yielded a temperature with large fluctuations:
Energy Average RMSD Fluct. Drift
Tot-Drift
-------------------------------------------------------------------------------
Temperature 298.899 10.1857 10.1724 -0.0360483
-1.80249
The fluctuations look fairly large, but I'm not sure if these are
reasonable for a 216-molecule system.
I then applied pressure to the system:
---[npt.mdp]---
; NPT equilibration
; Run parameters
integrator = md
nsteps = 50000
dt = 0.002
; Output control
nstxout = 100
nstvout = 100
nstenergy = 100
nstlog = 100
; Bond parameters
continuation = yes
constraint_algorithm = lincs
constraints = none
lincs_iter = 1
lincs_order = 4
; Neighborsearching
ns_type = grid
nstlist = 5
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; Temperature coupling is on
tcoupl = V-rescale
tc-grps = SYSTEM
tau_t = 0.1
ref_t = 300
; Pressure coupling is on
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 1e-4
; Periodic boundary conditions
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = no
---[npt.mdp]---
This simulation yields large fluctuations of temperature as well:
I tried playing around with the value of tau_p:
- tau_p = 5.0 ==> slow increase of density, not stable after 100 ps;
temperature with similar fluctuations as NVT (300K, RMSD 9); pressure
starts oscillating wildly after ~10 ps
- tau_p = 2.0 ==> increase of density, not stable at end of simulation;
huge fluctuations in potential energy (-4000 kJ; RMSD 1000)
- tau_p = 1.0 ==> increase of Density to 1850 in the first 10ps, then
stable with fluctuations of ~34 in the remaining 90ps; potential energy
with huge fluctuations
- tau_p = 0.2 (as used in Mol. Phys. 1994, 83, 381) ==> LINCS warnings
and posterior crash:
Step 23, time 0.046 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.079726, max 0.094860 (between atoms 672 and 673)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
242 243 36.2 0.2504 0.1919 0.1758
[...]
-------------------------------------------------------
Program mdrun_mpi_d, VERSION 4.0.5
Source code file: nsgrid.c, line: 357
Range checking error:
Explanation: During neighborsearching, we assign each particle to a grid
based on its coordinates. If your system contains collisions or parameter
errors that give particles very high velocities you might end up with some
coordinates being +-Infinity or NaN (not-a-number). Obviously, we cannot
put these on a grid, so this is usually where we detect those errors.
Make sure your system is properly energy-minimized and that the potential
energy seems reasonable before trying again.
Variable ci has value -2147483631. It should have been within [ 0 .. 84 ]
-------------------------------------------------------
I'll appreciate pointers as to what to try next!
Regards,
Pablo
--
Pablo Englebienne, PhD
Institute of Complex Molecular Systems (ICMS)
Eindhoven University of Technology, TU/e
PO Box 513, HG -1.26
5600 MB Eindhoven, The Netherlands
Tel +31 40 247 5349
More information about the gromacs.org_gmx-users
mailing list