[gmx-users] pbc problem

[shahab shariati]

Dear gromacs users

My system contains DOPC + CHOLESTEROLO + WATER + drug molecules in a
rectangular box.

I put drug molecule in 2 position: a) drug in the center of bilayer
membrane, b) drug inside water molecules in top leaflet.

For both positions, 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.

For both positions, I did equilibration in NPT ensemble with following mdp
file.
---------------------------------------------------------------------------------------
; Run parameters
integrator    = md        ; leap-frog integrator
nsteps        = 250000    ; 2 * 500000 = 1000 ps (1 ns)
dt            = 0.002        ; 2 fs
; Output control
nstxout        = 100        ; save coordinates every 0.2 ps
nstvout        = 100        ; save velocities every 0.2 ps
nstxtcout   = 100        ; xtc compressed trajectory output every 2 ps
nstenergy    = 100        ; save energies every 0.2 ps
nstlog        = 100        ; update log file every 0.2 ps
energygrps  = CHOL DOPC drg SOL
; Bond parameters
continuation    = no            ; Restarting after NVT
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.0        ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.0        ; short-range electrostatic cutoff (in nm)
rvdw        = 1.0        ; 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            ; More accurate thermostat
tc-grps        = CHOL_DOPC    drg SOL    ; three coupling groups - more
accurate
tau_t        = 0.5    0.5       0.5       ; time constant, in ps
ref_t        = 323     323       323     ; reference temperature, one for
each group, in K
; Pressure coupling is on
pcoupl        = Parrinello-Rahman        ; Pressure coupling on in NPT
pcoupltype    = semiisotropic            ; uniform scaling of x-y box
vectors, independent z
tau_p        = 5.0                    ; time constant, in ps
ref_p        = 1.0    1.0                ; reference pressure, x-y, z (in
bar)
compressibility = 4.5e-5    4.5e-5    ; isothermal compressibility, bar^-1
; 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         = 1
comm-mode       = Linear
comm-grps       = CHOL_DOPC_drg  SOL
; Scale COM of reference coordinates
refcoord_scaling = com

---------------------------------------------------------------------------------------
For 2 positions, I chechked tempreture and pressure fluctuation and box
dimention during equilibration. All things were good. When I saw trajectory
by VMD (npt.gro and npt xtc), I had pbc problem (some atoms leave box and
enter the box in opposit direction).

For position (a): I corrected pbc problem by

trjconv -f npt.xtc -s npt.tpr -n index.ndx -o 2npt.xtc -pbc mol -center

I selected CHOL_DOPC-drg group for centering. So problem was solved,
approximately.

For position (b) in which drug molecule is not in the center of lipid, I
can not use -center or I can not use group containing drug molecule for
centering.

How to fix this problem.

Any help will highly appreciated.