[gmx-users] Problem with equilibrated lipid bilayer structure

[Sebastien Cote]

Hello,
You might want to use TIPS3P (special CHARMM TIP3P) instead of regular TIP3P because some researchers observed that regular TIP3P can impact the area per lipid significantly when simulating CHARMM lipid bilayer (see reference below).  
Concerning the CHARMM36 lipid force field in Gromacs, there is also an issue about cutoffs. The switch cutoff scheme in the charmm software is different than in the gromacs software. It might impact your results for some bilayers, while it does not for others. 
You might find that previous post on the mailing list enlightening (and the cited reference) :
" Hi everyone,This message is just to let people know that there is finally a 
reference that we ask you to cite if you publish work using the CHARMM36 
lipid force field contribution available from the GROMACS website. The 
reference is provided in a new version of the forcefield.doc file, and 
is also given below:

Thomas J. Piggot, Ángel Piñeiro and Syma Khalid
Molecular Dynamics Simulations of Phosphatidylcholine Membranes: A 
Comparative Force Field Study
Journal of Chemical Theory and Computation
DOI: 10.1021/ct3003157
http://pubs.acs.org/doi/abs/10.1021/ct3003157

A validation of the conversion for DPPC and POPC membrane systems (in 
terms of single point energies calculated in both GROMACS and NAMD) is 
provided in the supporting information of this work.

Cheers

Tom"
Hope that helps,
Sebastien 
> Date: Mon, 15 Oct 2012 13:48:47 +0800
> From: jernej.zidar at gmail.com
> To: gmx-users at gromacs.org
> Subject: [gmx-users] Problem with equilibrated lipid bilayer structure
> 
> 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
> -- 
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists