[gmx-users] Problem with equilibrated lipid bilayer structure

[Peter C. Lai]

Is there a reason to switch water models? Gromacs supports both the 
traditional TIP3P and the CHARMM TIP3P (TIPS3P).

Since you need to reinitialize velocities, it may not be a good idea to use 
Nose-Hoover and Parinello-Rahman straight away. Furthermore, you need 
to gen_vel=yes since gromacs has no velocity from your previous charmm run 
and subjecting an otherwise frozen system to a parinello-rahman thermostat is 
asking for trouble. 

Also, pressure coupling for NPT should probably be semiisotropic for bilayers. 
Finally you may want to experiment with EnerPres; see also 
http://pubs.acs.org/doi/abs/10.1021/ct3003157 (although this isn't related to
your crash).

On 2012-10-15 01:48:47PM +0800, Jernej Zidar wrote:
> 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

-- 
==================================================================
Peter C. Lai			| University of Alabama-Birmingham
Programmer/Analyst		| KAUL 752A
Genetics, Div. of Research	| 705 South 20th Street
pcl at uab.edu			| Birmingham AL 35294-4461
(205) 690-0808			|
==================================================================