[gmx-users] Problem with equilibrated lipid bilayer structure

[Jernej Zidar]

Hi.
  I used CHARMM to prepare a cholesterol/POPC lipid bilayer. I
equilibrated it in CHARMM and the imported the resulting PDB to
GROMACS. In CHARMM I was the TIP3P water model, so before importing
the PDB I changed the atom types from the TIP3P's ones to SPCE ones.
Then I used this command to import  the PDB:
pdb2gmx -f lipid-wat.pdb -o bilayer.gro -p bilayer.top -i
bilayer-posres.itp -ff charmm36cgenff -water spce -noter -v -renum

  After the import I center the system in the unit cell with: editconf
-f bilayer.gro -o bilayer.gro -c

  Problem 1: The system size (unit cell) is not properly
computed/detected. The following numbers are for the same structure
after the CHARMM equilibration.
                   CHARMM reports the size as: 10.45820  x 10.45825  x
6.864382 (in nm; converted from Angstroms)
                   GROMACS states the size is: 12.30940  x 11.81980  x
7.138100 (in nm)

  Problem 2: While I'm able to run MD in CHARMM if I start from the
equilibrated structure, I am unable to do so in GROMACS. As the system
apparently blows up with the cryptic message:
Program mdrun, VERSION 4.5.5
Source code file: /build/buildd/gromacs-4.5.5/src/mdlib/pme.c, line: 538

Fatal error:
5 particles communicated to PME node 0 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.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

  I can minimize the system, but any attempt to run MD results in the
aforementioned message. Here's the MDP file I would like to use:
;define		= -DPOSRES	; position restrain the protein
; Run parameters
integrator	= md		; leap-frog integrator
nsteps		= 500000	; 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
nstenergy	= 100		; save energies every 0.2 ps
nstlog		= 100		; update log file every 0.2 ps
; 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.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		= Nose-Hoover		    ; More accurate thermostat
tc-grps		= LIPID SOL		; three coupling groups - more accurate
tau_t		= 0.5	0.5	        ; time constant, in ps
ref_t		= 300 	300	        ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl		= Parrinello-Rahman	    ; Pressure coupling on in NPT
pcoupltype	= isotropic		    ; uniform scaling of x-y box vectors, independent z
tau_p		= 5.0			        ; time constant, in ps
ref_p		= 1.0			        ; reference pressure, x-y, z (in bar)
compressibility = 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		= no		; Velocity generation is off
; COM motion removal
; These options remove motion of the protein/bilayer relative to the
solvent/ions
nstcomm         = 1
comm-mode       = Linear
comm-grps       = LIPID SOL
refcoord_scaling = all
- - - -

  I imagine I'm doing something wrong but I'm unable to be able to
pinpoint the error. I have also tried the NPT-simulated annealing path
suggested in the GROMACS' protein-membrane tutorial but to no avail.
I'm using the GROMACS version of the CHARMM36 lipid forcefield.

Thanks in advance for any advice,
Jernej Zidar