Pablo Englebienne wrote:
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!
If not all your bonds are constraints, you may have to use dt = 0.001. Also,
for the initial equilibration, you might try the Berendsen weak coupling method,
then after the system stabilizes, switch to Parrinello-Rahman.
-Justin
Regards,
Pablo
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
_______________________________________________
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
