[gmx-users] change of bilayer structure during NVT equilibration

[Justin Lemkul]

On 10/6/13 10:03 AM, shahab shariati wrote:
> Dear gromacs users
>
> My system contains DOPC + CHOLESTEROLO + WATER in a rectangular box.
>
> I did energy minimization successfully with following mdp file.
> --------------------------------------------------------------------------------------
> ; Parameters describing what to do, when to stop and what to save
> integrator    = steep        ; Algorithm (steep = steepest descent
> minimization)
> emtol        = 1000.0      ; Stop minimization when the maximum force <
> 1000.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
> 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
> ---------------------------------------------------------------------------------------
>
> After energy minimization, I saw obtained file (em.gro) by VMD. All things
> were true
> and intact.
>
> I did equilibration in NVT ensemble with following mdp file.
> --------------------------------------------------------------------------------------
> title        = NVT equilibration
> ; Run parameters
> integrator    = md                   ; leap-frog integrator
> nsteps        = 7500000                  ; 2 * 7500000 = 15 ns
> dt            = 0.002                   ; 2 fs
> ; Output control
> nstxout        = 1000                   ; save coordinates every 0.2 ps
> nstvout        = 1000                   ; save velocities every 0.2 ps
> nstxtcout   = 1000
> nstenergy    = 1000                   ; save energies every 0.2 ps
> nstlog        = 1000                   ; update log file every 0.2 ps
> ; Bond parameters
> continuation    = no            ; first dynamics run
> constraint_algorithm = lincs    ; holonomic constraints
> constraints    = all-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 cels
> nstlist        = 5                    ; 10 fs
> rlist        = 1.2                  ; short-range neighborlist cutoff (in
> nm)
> rcoulomb    = 1.2                  ; short-range electrostatic cutoff (in
> nm)
> rvdw        = 1.2                  ; short-range van der Waals cutoff (in
> nm)
> ; 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        = V-rescale            ; modified Berendsen thermostat
> tc-grps        = CHOL_DOPC SOL      ; three coupling groups - more accurate
> tau_t        = 0.1   0.1            ; time constant, in ps
> ref_t        = 323   323            ; reference temperature, one for each
> group, in K
> ; Pressure coupling is off
> pcoupl        = no                 ; no pressure coupling in NVT
> ; Periodic boundary conditions
> pbc            = xyz                 ; 3-D PBC
> ; Dispersion correction
> DispCorr    = EnerPres             ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel        = yes                   ; assign velocities from Maxwell
> distribution
> gen_temp    = 323                   ; temperature for Maxwell distribution
> gen_seed    = -1                   ; generate a random seed
> ; COM motion removal
> ; These options remove motion of the protein/bilayer relative to the
> solvent/ions
> nstcomm        = 5
> comm-mode    = Linear
> comm-grps    = CHOL_DOPC SOL
> ---------------------------------------------------------------------------------------
>
> After 15 ns NVT equilibration, I saw obtained file (nvt.gro) by VMD.
>
> Unfortunately, rectangular shape of box was converted to cylinder shape.
> DOPC and CHOL
> molecules moved from center of box to environs of box.
>
> What is reason of this issue?
>

Depends on how you prepared the system.  Probably there were voids in the 
solvent or lipid headgroups that caused distortion in the system.  Restraints on 
the lipid headgroups along z could help.  I also see no reason for 15 ns of NVT. 
  A short NVT is all that should be necessary, on the order of 100-500 ps, 
especially when using V-rescale, which relaxes quickly to the target temperature.

-Justin

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

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441

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