From: Allan Haldane (allan.haldane_at_mail.mcgill.ca)
Date: Tue Oct 10 2006 - 11:40:43 CDT
Hi,
I am quite new to MD, so I thought I would try to reproduce the result
Gullingsrud/Schulten (Biophysical Journal 86:3496-3509 (2004)) found
for the lateral pressure profile of a POPC membrane. It doesn't look
right, though. I have ideas of what might be wrong, but maybe someone
can quickly put me in the right direction.
Here is a plot of the pressure tensor components I get averaged over 5
ns (there is a bit of the water regions removed from the plot):
http://www.ugrad.physics.mcgill.ca/~ahalda/upload/protein/NAMD/press26b1.png
Here are the bonded + nonbonded components: (compare to simulation D4
in figure 6 of the paper)
http://www.ugrad.physics.mcgill.ca/~ahalda/upload/protein/NAMD/profComponents.png
I start with the popc membrane template (popc_box) that comes with vmd,
and after adding some extra water (10 A above and below) I heat over 10
ps (by langevin dynamics) and then equilibrate it for 1 ns at NPT. N,
P, T and V are indeed constant for most of this. Actually I calculate
that popc_box has an area per lipid of 76 A^2 By the end of the
equilibration which is high. (I calculate area/lipid as cell
area/lipids in monolayer (30). Is it truly that simple?) I see one of
Gullingsrud/Schulten's POPC runs was at around this area/lipid, though
(simulation D4).
I then do a pressure profile for 5 ns at NVT with nonbonded part done
seperately (with NAMD 2.6.b1). (I also tried with 2.6.2 + Ewald and got
quite different results...)
I see two things happen during the pressure profile run: First, the
global pressure increases by 53 bar over the 5 ns. It is initially -50
bar (unlike at the end of the NPT equilibration), but from a post on
this list this is apparently normal for NVT simulations with liquid.
Second, the bilayer moves quite a lot in the z (normal) direction, by a
total shift of 4 A. (For the NAMD 2.6.2 run I did, it shifted by 10 A)
This happens gradually and mostly continuously through the run.
Here is my config file for the profile calculation. The equilibration
uses the same parameters except that: langevinPiston on,
langevinDamping 10, pressureProfile off and I output more frequently.
My system is 50 by 46 by 69 A. I do nonbonded seperately, as shown in
the manual.
Thanks!
Allan
-----------------------
structure start.psf
coordinates requil.pdb
parameters par_all27_prot_lipid.inp
parameters par-extraterms.inp
paraTypeCharmm on
extendedSystem requil.xsc
velocities requil.vel
outputEnergies 1000
outputTiming 1000
xstFreq 100000
dcdFreq 10000
wrapAll on
restartname output/restart
restartfreq 100000
timestep 1
nonBondedFreq 2
fullElectFrequency 4
stepsPerCycle 20
switching on
switchDist 8.5
cutoff 10
pairlistdist 11.5
Pme on
PmeGridsizeX 64
PmeGridsizeY 64
PmeGridsizeZ 128
exclude scaled1-4
1-4scaling 1.0
langevin on
langevinDamping 1
langevinTemp 310
langevinHydrogen off
langevinPiston off
langevinPistonTarget 1.01325
langevinPistonPeriod 200
langevinPistonDecay 100
langevinPistonTemp 310
#these two only should affect equilibration:
useGroupPressure yes #hydrogen excluded
useFlexibleCell yes
binaryoutput off
outputname output/equil_out
pressureProfile on
pressureProfileSlabs 120
pressureProfileFreq 10000
run 5000000
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:42:40 CST
