Hexagonal cell instability during lipid bilayer sim in NAMD 2.9

From: Michael Bokoch (mbokoch_at_gmail.com)
Date: Mon Jul 07 2014 - 21:42:54 CDT

Dear NAMD forum,

I’m writing to ask for help simulating a membrane protein embedded in a lipid bilayer using a hexagonal unit cell in NAMD 2.9 with CHARMM forcefields.

The basic problem is that the equilibrated unit cell slowly morphs shape from a hexagon (with flat sides along the y-axis) to a cube (rotated 45 deg from normal to x-y plane) over the course of a 100 ns simulation. Otherwise, the bilayer and protein seem to behave properly.

This problem does not become apparent during the 4 ns equilibration, during which time the hexagonal cell retains it’s original orientation.

Do you think this is a Pressure Control problem, and I need to constrain the area of the unit cell in some way? For example, should I use:

useConstantRatio yes
-or-
useConstantArea yes

Currently, I am not using either of these options.
 
I only use the following options that seem relevant to this discussion:
wrapNearest on #on for hexagon
useGroupPressure yes
useFlexibleCell yes

Are there other options that I could be missing that are crucial for hexagonal simulations?

More details and specific questions are below, including the steps I’ve taken to create this simulation, and troubleshoot it. Perhaps this will be helpful to others trying to set up similar simulations.

With sincere thanks,
Mike

-- 
Michael Bokoch, MD, PhD
University of California, San Francisco
---
1.) Here are the basis vectors I use for the initial minimization, and a description of the geometry that is correct for a hexagonal cell with flat sides oriented along y.
cellBasisVector1    90.2985  0       0
cellBasisVector2    45.1493 78.2008  0
cellBasisVector3     0       0      91.191
cellOrigin           0.3541  0.3237  1.4828
# Vector 1 is (r, 0, 0)
# Vector 2 is (r/2, sqrt(3)*r/2, 0)
# Vector 3 is (0, 0, h)
# Where r is distance between two adjacent hexagonal cells
# which also equals the inner diameter of the hexagon
# h is the height of the unit cell
 
2.)  Here are the periodic cell parameters from the .xst file over the course of the simulation, demonstrating the problem:
End of equilibration:
90.9184 0 0 45.4592 71.5339 0 0 0 93.4775
After 5 ns:
91.6406 0 0 45.8203 72.1299 0 0 0 91.8068
After 25 ns:
92.0509 0 0 46.0255 67.4482 0 0 0 97.4556
After 50 ns:
104.511 0 0 52.2562 60.3603 0 0 0 96.2181
After 100 ns:
105.794 0 0 52.8978 57.2271 0 0 0 100.165
3.) Is it true that NAMD prefers a hexagonal cell with flat-sides oriented along the y-axis?  This is suggested by an earlier post to this mailing list.
http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2005-2006/3006.html
However, this post suggests that “r” (in my vector definitions above) is the inner radius of the hexagon, not the inner diameter (which equals the distance between two adjacent cell centers).
I believe we have used trial-and-error to prove that the basis vector length “r” should be equal to the cell-to-cell distance, not the inner radius of the hexagon, nor the edge length of the hexagon.
We have tried minimizing our model with all of these options, and the result is overcrowding of waters and lipids around the cell periphery.
We also tried using basis vectors that define adjacent vertices on the Bravais lattice (see #4) for the edge length of our hexagon, but this resulted in a unit cell that was much too small with overcrowding of waters and lipids.
http://en.wikipedia.org/wiki/Bravais_lattice
4.) The original hexagonal model of lipid-embedded protein was generated in CHARMM-GUI Membrane builder.
http://www.charmm-gui.org/?doc=input/membrane
Interestingly, the output of CHARMM-GUI Step 5 yields a hexagonal cell that is rotated 15 degrees clockwise from the orientation with flat-sides along the y-axis. 
I had to use VMD to manually rotate the entire .pdb 15 degrees counterclockwise such that the flat-sides are again parallel to the y-axis, prior to minimization.

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:22:34 CST