[gmx-users] unstable system

Justin Lemkul jalemkul at vt.edu
Wed May 8 14:37:25 CEST 2013



On 5/8/13 2:35 AM, Shima Arasteh wrote:
> OK.
>
> 1. Exact commands given in the preparation protocol (EM and equilibration)
>
> EM:
> # grompp -f minim.mdp -c input.gro -p topol.top -o minim.tpr
> #mdrun -deffnm minim
>
> NVT:
> #grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
> #mdrun -deffnm nvt -v
>
> NPT:
> # grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr
> # mdrun_mpi -deffnm npt
>
>
> *** Next removed some disturbing water molecules, so edited topol file and made a new index file index_mod_2.ndx:
> Ran a new EM:
> # grompp -f npt.mdp -c npt_mod.gro -p topol.top -o npt_minim.tpr
> # mdrun -deffnm npt_minim
>
>
> Ran a new NPT:
> # grompp -f npt.mdp -c npt_minim.gro -p topol.top -n index_mod_2.ndx -o npt.tpr
> # mdrun_mpi -deffnm npt
>
> I repeated the NPT the same as above in 3 steps:
> a) restraints on lipid phosphorous and protein for 1 ns
>
> b) restraints on protein for 1ns
>
> c) restraint on protein with less force for 2 ns
>

All 3 of these steps complete successfully?

>
> MD:
> # grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index_mod_2.ndx -o md_0_1.tpr
> # mdrun -deffnm npt
>

There are a few discrepancies in the commands you have shown here, in particular 
naming differences between input and output file names.  I ask for exact 
commands (not what you think you remember) for a very important reason - the 
diagnostic information you have provided (from the .log file) seems to indicate 
that the previous equilibrated state was not passed along to the MD step, so I 
suspect you've left out a checkpoint file somewhere or you have regenerated 
velocities by accident.

As a tip, write all of your commands in a shell script and execute the script. 
Then, if something goes wrong and someone asks for exact commands, just copy and 
paste.  I still can't tell if this is actually what you did.

>
> 3 particles communicated to PME node 2 are more than 2/3 times the cut -
> off out of the domain decomposition cell of their charge group in
> dimension y.
> This usually means that your system is not well equilibrated.
>
>
> 2. The .mdp files being used for all steps, most importantly NPT and MD
>
> npt.mdp:
>
> ;NPT equlibration Dimer-POPC-Water - CHARMM36
> define        = -DPOSRES
> integrator      = md            ; leap-frog integrator
> nsteps          =1000000         ; 1 * 1000000 = 1000 ps
> dt              = 0.001         ; 1 fs
> ; Output control
> nstxout         = 2000          ; save coordinates every 0.4 ps
> nstvout         = 2000           ; save velocities every 0.2 ps
> nstenergy       = 1000           ; save energies every 0.2 ps
> nstlog          = 1000           ; update log file every 0.2 ps
>
> continuation    = yes            ; first dynamics run
> constraint_algorithm = lincs    ; holonomic constraints
> constraints     = h-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 cells
> nstlist         = 5             ; 10 fs
> rlist           = 1.2           ; short-range neighborlist cutoff (in nm)
> rlistlong       = 1.4
> rcoulomb        = 1.2           ; short-range electrostatic cutoff (in nm)
> rvdw            = 1.2           ; short-range van der Waals cutoff (in nm)
> vdwtype         = switch
> rvdw_switch     = 1.0
> ; Electrostatics
> coulombtype     = PME           ; Particle Mesh Ewald for long-range electrostatics
> pme_order       = 4             ; cubic interpolation
> fourierspacing  = 0.16          ; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl          = Nose-Hoover     ; modified Berendsen thermostat
> tc-grps         = Protein_POPC    Water_and_ions    ; two coupling groups - more accurate
> tau_t           = 0.5   0.5       ; time constant, in ps
> ref_t           = 310   310     ; reference temperature, one for each group, in K
> pcoupl          = Berendsen            ; no pressure coupling in NVT
> pcoupltype      = semiisotropic
> tau_p           = 4
> ref_p           = 1.01325 1.01325
> compressibility = 4.5e-5 4.5e-5
>
> ; Periodic boundary conditions
> pbc             = xyz           ; 3-D PBC
> ; Dispersion correction
> DispCorr        = no    ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel         = no           ; assign velocities from Maxwell distribution
> ;gen_temp       = 310           ; temperature for Maxwell distribution
> ;gen_seed       = -1            ; generate a random seed
> nstcomm         = 1
> comm_mode       = Linear
> comm_grps       = Protein_POPC Water_and_ions
> refcoord_scaling = com
>
>
> md.mdp:
>
> title        = Production run forGLC-Water-POPC system
> define        = -DPOSRES
> ; Parameters describing the details of the NVT simulation protocol
> integrator    = md
> dt        = 0.001
> nsteps        = 2500000    ; Number of steps to run (0.002 * 2500000 = 5 ns)
>
> ; Parameters controlling output writing
> nstxout        = 1000        ; Write coordinates to output .trr file every 10 ps
> nstvout        = 2000    ; Write velocities to output .trr file every 10 ps
> nstenergy    = 1000        ; Write energies to output .edr file every 2 ps
> nstlog        = 5000        ; Write output to .log file every 2 ps
> ; Parameters describing neighbors searching and details about interaction calculations
> ns_type        = grid        ; Neighbor list search method (simple, grid)
> nstlist        = 5        ; Neighbor list update frequency (after every given number of steps)
> rlist        = 1.2        ; Neighbor list search cut-off distance (nm)
> rlistlong       = 1.4
> rcoulomb    = 1.2        ; Short-range Coulombic interactions cut-off distance (nm)
> rvdw        = 1.2        ; Short-range van der Waals cutoff distance (nm)
> pbc        = xyz        ; Direction in which to use Perodic Boundary Conditions (xyz, xy, no)
> vdwtype         = switch
> rvdw_switch     = 1.0
> ; Parameters for treating bonded interactions
> continuation    = yes        ; Whether a fresh start or a continuation from a previous run (yes/no)
> constraint_algorithm = LINCS    ; Constraint algorithm (LINCS / SHAKE)
> constraints    = all-bonds    ; Which bonds/angles to constrain (all-bonds / hbonds / none / all-angles / h-angles)
> lincs_iter    = 1        ; Number of iterations to correct for rotational lengthening in LINCS (related to accuracy)
> lincs_order    = 4        ; Highest order in the expansion of the constraint coupling matrix (related to accuracy)
>
> ; Parameters for treating electrostatic interactions
> coulombtype    = PME        ; Long range electrostatic interactions treatment (cut-off, Ewald, PME)
> pme_order    = 4        ; Interpolation order for PME (cubic interpolation is represented by 4)
> fourierspacing    = 0.16        ; Maximum grid spacing for FFT grid using PME (nm)
>
> ; Temperature coupling parameters
> tcoupl        = Nose-Hoover            ; Modified Berendsen thermostat using velocity rescaling
> tc-grps        = Protein_POPC Water_and_ions        ; Define groups to be coupled separately to temperature bath
> tau_t        = 0.5    0.5             ; Group-wise coupling time constant (ps)
> ref_t        = 310     310            ; Group-wise reference temperature (K)
>
> ; Pressure coupling parameters
> pcoupl        = Parrinello-Rahman        ; Pressure coupler used under NPT conditions
> pcoupltype    = semiisotropic            ; Isotropic scaling in the x-y direction, independent of the z direction
> tau_p        = 2.0                ; Coupling time constant (ps)
> ref_p        = 1.01325 1.01325        ; Reference pressure for coupling, x-y, z directions (bar)
> compressibility = 4.5e-5    4.5e-5        ; Isothermal compressibility (bar^-1)
>
> ; Miscellaneous control parameters
> ; Dispersion correction
> DispCorr    = EnerPres        ; Dispersion corrections for Energy and Pressure for vdW cut-off
> ; Initial Velocity Generation
> gen_vel        = no            ; Velocity is read from the previous run
> ; Centre of mass (COM) motion removal relative to the specified groups
> nstcomm        = 1            ; COM removal frequency (steps)
> comm_mode    = Linear        ; Remove COM translation (linear / angular / no)
> comm_grps    =Protein_POPC Water_and_ions    ; COM removal relative to the specified groups
>
>

I see no explanation as to why the run fails, aside from what I have said above.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
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