[gmx-users] LINCS WARNING - Protein-Membrane-system

[Justin Lemkul]

On 1/1/13 10:53 AM, Shima Arasteh wrote:
>
>
>   If I tell you that I changed the constraints and now the NVT is running, would that be ok? However I don't know the scientific reason of this incident!
>
> But I bring you the EM.mdp as follow:
> define         = -DSTRONG_POSRES
> integrator    = steep        ; Algorithm (steep = steepest descent minimization)
> emtol        = 100.0      ; Stop minimization when the maximum force < 100.0 kJ/mol/nm
> emstep            = 0.01      ; Energy step size
> nsteps        = 50000          ; Maximum number of (minimization) steps to perform
> ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
> nstlist        = 1            ; Frequency to update the neighbor list and long range forces
> ns_type        = grid        ; Method to determine neighbor list (simple, grid)
> rlist        = 1.2        ; Cut-off for making neighbor list (short range forces)
> coulombtype    = PME        ; Treatment of long range electrostatic interactions
> rcoulomb    = 1.2        ; Short-range electrostatic cut-off
> rvdw        = 1.2        ; Short-range Van der Waals cut-off
> pbc        = xyz         ; Periodic Boundary Conditions
>
>
> And the NVT.mdp as follow (after change of all-bond constraint to h-bond):
> title        = NVT equilibration for dimer-POPC   -  CHARMM36
> define        = -DPOSRES    ; Protein is position restrained (uses the posres.itp file information)
> ; Parameters describing the details of the NVT simulation protocol
> integrator    = md        ; Algorithm ("md" = molecular dynamics [leap-frog integrator]; "md-vv" = md using velocity verlet; sd = stochastic dynamics)
> dt        = 0.002        ; Time-step (ps)
> nsteps        = 50000    ; Number of steps to run (0.002 * 50000 = 100 ps)
>
> ; Parameters controlling output writing
> nstxout        = 1        ; Write coordinates to output .trr file every 2 ps
> nstvout        = 1000        ; Write velocities to output .trr file every 2 ps
> nstenergy    = 1000        ; Write energies to output .edr file every 2 ps
> nstlog        = 1000        ; 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     = 0.8
> ; Parameters for treating bonded interactions
> continuation    = no        ; Whether a fresh start or a continuation from a previous run (yes/no)
> constraint_algorithm = LINCS    ; Constraint algorithm (LINCS / SHAKE)
> constraints    = h-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        = v-rescale            ; Modified Berendsen thermostat using velocity rescaling
> tc-grps        = Protein POPC    SOL_CL    ; Define groups to be coupled separately to temperature bath
> tau_t        = 0.1    0.1    0.1        ; Group-wise coupling time constant (ps)
> ref_t        = 310     310    310        ; Group-wise reference temperature (K)
>
> ; Pressure coupling parameters
> pcoupl        = no         ; Under NVT conditions pressure coupling is not done
>
> ; Miscellaneous control parameters
> ; Dispersion correction
> DispCorr    = EnerPres    ; Dispersion corrections for Energy and Pressure for vdW cut-off
> ; Initial Velocity Generation
> gen_vel        = yes        ; Generate velocities from Maxwell distribution at given temperature
> gen_temp    = 310        ; Specific temperature for Maxwell distribution (K)
> gen_seed    = -1        ; Use random seed for velocity generation (integer; -1 means seed is calculated from the process ID number)
> ; 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 SOL_CL    ; COM removal relative to the specified groups
>
>
> Any suggestions please?
>

Yes, the problem is the constraints.  You minimized without any, and then 
constrained all bonds during NVT.  During EM, you can get bad geometry sometimes 
when atomic clashes force atoms to move in bad ways.  Without constraints, mdrun 
doesn't object to these deformations.  Then, as you see from step zero of the 
NVT, trying to apply constraints to these bad areas immediately fails.

-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

========================================