Re: question on membrane crushing

From: Arneh Babakhani (ababakha_at_mccammon.ucsd.edu)
Date: Tue Mar 14 2006 - 15:59:53 CST

Hello Himanshu,

I'm wouldn't say that I observe no systematic contraction. My simulation
thus far is only about 200 ps, not nearly enough time to say that
there's no systemic or gradual contraction.

What I did observe is this:
Without a "melting" step (in which the headgroups were fixed, and the
membrane ran for 20 ps at a high temp, to randomize the tails) and
without a gradual release of constraints, the membrane is "crushed".
Over a timespan of about 200 ps, the thickness contracts from 38
Angstroms to 32. A very dramatic effect . . .

With a "melting" step and a gradual release of constraints (as described
previously), over that same timespan of 200 ps, the contraction is less
dramatic. The membrane thickness contracts from 38 to 37 Angstroms. At
the end of this 200 ps run, the thickness seems to have steadied out at
37 Angstroms, but again, this isn't a large enough timescale to tell if
it really has equilibrated. The only thing I can tell is that there is
no dramatic "crushing", and I think this can be attributed to the two
steps mentioned above.

To answer your other questions:
Used periodic boundary conditions, with dimensions of about 100 x 100 x
60, to accommodate my entire system of 200 lipids (100 per leaflet) and
water. Constant pressure and temp applied (at 310K). I did not use
flexcell, constant area, nor did I use constant ratio. (I've read many
mixed opinions about these parameters, some say they're essential for
membrane simulations. Others say that they're not. Do you have any
thoughts?). Anyway, without those parameters, I get a decent thickness
and my lipid orderparameters seem right, but again, I need to
equilibrate longer before I can be sure.

I've pasted my config file below,

Thanks,

Arneh

# Heating of my solvated membrane

## Adjustable Parameters

structure ../newsolvated.psf
coordinates ../newsolvated.pdb

binvelocities ../minimizationandheating/lipwaterheatup-2i.vel
bincoordinates ../minimizationandheating/lipwaterheatup-2i.coor
extendedsystem ../minimizationandheating/lipwaterheatup-2i.xsc

set temperature 310
set outputname lipwaterconstpress-1

firsttimestep 0

#############################
##Simulation Parameters
##############################

paraTypeCharmm on
parameters /u1/ababakha/Common/toppar/par_all27_prot_lipid.inp

rigidbonds all

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 12.0
switching on
switchdist 10.0
pairlistdist 13.5

# Integrator Parameters
timestep 2
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10

# Constant Temperature Control
langevin on ;# do langevin dynamics
langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
langevinTemp $temperature
langevinHydrogen off ;# don't couple langevin bath to hydrogens

# Constant Pressure Control (variable volume)
useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no

langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
langevinPistonPeriod 100.
langevinPistonDecay 200.
langevinPistonTemp $temperature

wrapAll on

# PME (for full-system periodic electrostatics)
PME yes
PMEGridSizeX 150
PMEGridSizeY 150
PMEGridSizeZ 90

# Output
outputName $outputname

restartfreq 500 ;# 500steps = every 1ps
dcdfreq 500
outputEnergies 100
outputPressure 100

#########Execution
run 12500

Himanshu Khandelia wrote:

>Arneh,
>
>This is very very interesting. So you do not observe any systematic
>contraction in your membrane in the NPT ensemble ? I do not recall any
>simulation longer than 10 ns in CHARMM/NAMD where this has been achieved.
>
>Could you please answer a few questions, which will really help us out:
>
>1. how long are your simulations ?
>2. what kind of boundary conditions are you using ?
>3. what exactly do you mean by the NPT ensemble ? Are all the dimensions
>of the unit cell flexible ? Is there a constant ratio in the x-y plane ?
>Or is this an NPzAT ensemble, where the area is fixed, and the pressure
>can fluctuate in the z-direction (the direction of the bilayer normal). It
>would also help if you could show us a config file,
>
>Thanks very much !
>
>-Himanshu
>
>===================================================
>Himanshu Khandelia
>
>Doctoral Candidate,
>Kaznessis Research,
>Department of Chemical Engineering and Materials Science,
>University of Minnesota
>
>Mailing Address:
>
>499, Walter Library,
>117, Pleasant St. SE,
>Minneapolis, MN 55455
>
>Phone(o): 612-624-4945
>===================================================
>
>
>On Mon, 13 Mar 2006, Arneh Babakhani wrote:
>
>
>
>>Hi Longzhu,
>>
>>Ok, I think I know what's going on here, regarding this "crushing" effect.
>>
>>Through trial and error, I figured out that the most effect way to
>>equilibrate my membrane (in my case, a DMPC one) is the following:
>>
>>1. Build your membrane, using whatever template and script. I would
>>make the thickness and area per lipid a little bit larger than your
>>target thickness and area. So for instance, for DMPC (which normally has
>>a thickness of 36-37 Angstroms), I originally constructed my membrane to
>>have a thickness of 38.
>>
>>2. Fix the tails. Minimize the headgroups
>>3. Fix the headgroups. Minimize the tails.
>>4. "Melt the tails". With the headgroups fixed, heatup your membrane
>>to a high temp. I set mine to 450K, for about 20 ps. This step is
>>crucial to randomize your tails.
>>5. Solvate your system
>>6. Minimize.
>>
>>Then perform constrained heating, and release the constraints gently. I
>>think this is really crucial. There's several ways you can do this, I
>>did the following:
>>7. Ran 20 ps heating to 310K, with HG/Tails/Water constrained with
>>force constants of 10/5/2 kcal per mol, respectively.
>>8. Ran 20 ps heating to 310K, with HG/Tails/Water constrained with
>>force constants of 5/2/0 kcal per mol, respectively.
>>9. Ran 20 ps heating to 310K, with HG/Tails/Water constrained with
>>force constants of 2/0/0 kcal per mol, respectively.
>>10. Ran 20 ps heating to 310K, with HG/Tails/Water constrained with
>>force constants of 0/0/0 kcal per mol, respectively. (In other words,
>>all free, no constraints).
>>
>>Then I applied constant pressure and ran my production runs.
>>
>>Doing this, I found my membrane to equilibrate and retain more
>>recognizable features. No crushing.
>>
>>Hope this helps,
>>
>>Arneh
>>
>>Longzhu Shen wrote:
>>
>>
>>
>>>Dear All,
>>>
>>>I was trying simulating the POPE membrane with namd. I manually inserted a
>>>peptide into the POPE model, ran minimization, heated the system, performed
>>>position restrained md and free md. Everything went with along the
>>>simulation. However, when I checked the trajectories with VMD, I found the
>>>water on the two ends both moved toward the middle of the system. This made
>>>the two leaves of the POPE further inserted to each other and lipid bilayer
>>>severely crushed. The energy of the system looked stable during the whole
>>>process with VDW slightly less than zero. I'd attache the configuration file
>>>and wonder whether any kind guy could point out where I made mistakes. Many
>>>thanks.
>>>
>>>sincerely,
>>>
>>>Longzhu Shen
>>>
>>>
>>>------------------------------------------------------------------------
>>>
>>># Input Force-Field Parameters
>>>paraTypeCharmm on
>>>parameters par_all27_prot_lipid.inp
>>>
>>>#
>>>#molecules
>>>
>>>
>>>set import ions_added
>>>set export eqmd
>>>
>>>structure ${import}.psf
>>>coordinates ${import}.pdb
>>>#bincoordinates ${import}.coor
>>>#binvelocities ${import}.vel
>>>#extendedSystem ${import}.xsc
>>>
>>>set init_temp 0
>>>set desired_temp 310
>>>
>>>#
>>># Output files & writing frequency for DCD
>>># and restart files
>>>#
>>>
>>>outputname ${export}/eqmd
>>>binaryoutput no
>>>restartfreq 1000
>>>binaryrestart yes
>>>#dcdFile output/solvated_heat_out.dcd
>>>
>>>#
>>># Frequencies for logs and the xst file
>>>#
>>>outputEnergies 100
>>>outputTiming 100
>>>outputPressure 100
>>>xstFreq 1000
>>>dcdFreq 1000
>>>
>>>#
>>># Timestep & friends
>>>#
>>>
>>>timestep 1.0 ;# 1fs/step
>>>nonbondedFreq 1
>>>rigidBonds all
>>>rigidTolerance 0.00000001
>>>fullElectFrequency 4
>>>stepspercycle 20
>>>
>>>
>>>#
>>># Simulation space partitioning
>>>#
>>>switching on
>>>switchDist 10
>>>cutoff 12
>>>pairlistdist 16
>>>
>>>#
>>># Basic dynamics
>>>#
>>>temperature $init_temp
>>>COMmotion no
>>>dielectric 1.0
>>>exclude scaled1-4
>>>1-4scaling 1.0
>>>
>>>
>>>#
>>># Particle Mesh Ewald parameters.
>>>#
>>>
>>>PME yes
>>>PMEGridSizeX 54
>>>PMEGridSizeY 54
>>>PMEGridSizeZ 75
>>>
>>>#
>>># Periodic boundary things
>>>#
>>>wrapWater on
>>>wrapNearest on
>>>wrapAll on
>>>
>>>margin 3
>>>
>>>
>>>cellBasisVector1 52.45 0 0
>>>cellBasisVector2 0 51.35 0
>>>cellBasisVector3 0 0 74.67
>>>cellOrigin -22.69 -22.61 -0.76
>>>
>>>#
>>># Fixed atoms for initial heating-up steps
>>>#
>>>fixedAtoms on
>>>fixedAtomsForces on
>>>fixedAtomsFile fix_backbone.pdb
>>>fixedAtomsCol B
>>>
>>># Restrained atoms for initial heating-up steps
>>>#
>>>constraints on
>>>consRef restrain_ca.pdb
>>>consKFile restrain_ca.pdb
>>>consKCol B
>>>
>>>#
>>># Langevin dynamics parameters
>>>#
>>>langevin on
>>>langevinDamping 1
>>>langevinTemp $desired_temp #
>>>langevinHydrogen on
>>>
>>>langevinPiston on
>>>langevinPistonTarget 1.01325
>>>langevinPistonPeriod 200
>>>langevinPistonDecay 100
>>>langevinPistonTemp $desired_temp #
>>>
>>>useGroupPressure yes
>>>
>>># The actual minimisation and heating-up
>>>
>>># run one step to get into scripting mode
>>>minimize 0
>>>
>>># turn off pressure control until later
>>>langevinPiston off
>>>
>>># minimize nonbackbone atoms
>>>minimize 2500 ;#
>>>output ${export}/min_fix
>>>
>>>#
>>># min all atoms
>>>#
>>>fixedAtoms off
>>>minimize 2500 ;#
>>>output ${export}/min_all
>>>
>>>#
>>># heat with Ps restrained
>>>#
>>>set temp $init_temp;
>>>while { $temp <= $desired_temp } { ;#
>>>langevinTemp $temp
>>>run 10000 ;#
>>>output ${export}/heating_Pr
>>>set temp [expr $temp + 31]
>>>}
>>>
>>>#
>>># equilibrate volume with Ps restrained
>>>#
>>>langevinPiston on
>>>run 100000 ;#
>>>output ${export}/equil_Pr
>>>
>>>#
>>># equilibrate volume without restraints
>>>#
>>>constraintScaling 0
>>>run 500000 ;#
>>>output ${export}/equil
>>>
>>>
>>>
>>>
>>--
>>----------------------------
>>Arneh Babakhani
>>University of California at San Diego
>>Physical Chemistry / Department of Chemistry & Biochemistry
>>Laboratory of Prof. J. A. McCammon
>>9500 Gilman Drive MC 0365
>>La Jolla, CA 92093-0365
>>
>>(619)895-6540
>>(858)534-4974 (FAX)
>>
>>http://mccammon.ucsd.edu/~ababakha/
>>ababakha_at_mccammon.ucsd.edu
>>
>>
>>
>>
>>
>>
>
>
>

-- 
----------------------------
Arneh Babakhani
University of California at San Diego
Physical Chemistry / Department of Chemistry & Biochemistry
Laboratory of Prof. J. A. McCammon
9500 Gilman Drive MC 0365
La Jolla, CA 92093-0365
(619)895-6540
(858)534-4974 (FAX)
http://mccammon.ucsd.edu/~ababakha/
ababakha_at_mccammon.ucsd.edu

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:43:23 CST