[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field
Justin A. Lemkul
jalemkul at vt.edu
Tue Oct 27 13:03:01 CET 2009
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
========================================
More information about the gromacs.org_gmx-users
mailing list