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

[shahab shariati]

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?

Should I use new parameters in mdp file?

I am beginner in gromacs. Please help me to fix this issue.

Any help will highly appreciated

Best wishes