From: Himanshu Khandelia (hkhandel_at_dtc.umn.edu)
Date: Tue Mar 14 2006 - 16:07:32 CST
Hi Arneh,
Thanks for the details.
My guess is that you bilayer will start shrinking (in area) in less than
10 ns. But then again, I would be very glad if I am wrong.
In CHARMM, NPAT and NP-gamma-T (constant surface tension) have usually
been recommended for planar interfaces, because of the anisotropy in the
pressure tensor. There is some debate about whether this is the right
thing to do, but there is no choice if the NPT simulation ends up with an
indefinetely shrinking cell, to reduce the surface tension.
Anyhow, it will be nice if you can post your findings, when you have run
for a few nanoseconds. All the best !
-Himanshu
On Tue, 14 Mar 2006, Arneh Babakhani wrote:
> 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:41:44 CST