Lipid_bilayer_area_per_lipid_problem_using_NPT_CharmmForcefield

From: yousef nademi (yousef.namd_at_gmail.com)
Date: Thu Jul 06 2017 - 20:49:48 CDT

Dear all,

I build POPC lipid bilayer with the dimension of 18*18 using Charmm36
Forcefield and membrane builder plugin in VMD. the total number of lipid
molecules generated in this stage is 921 lipid molecules which mean each
leaflet has approximately 461 lipid molecules.
the experimental value for area per lipid for DOPC lipid bilayer is 65.8
A^2 at 271K. I used following procedure to equilibrate the membrane at 1
bar and 310k.
1- I run NVT simulation by keeping everything fixed except lipid tail to
create the fluid like structure as suggested in the membrane-protein
tutorial.
2- then, I run NPT simulation with following configuration file.
the temperature of my simulation was 310K that means I should have a less
packed structure in my system in comparison with mentioned experimental
value. to obtain the experimental value of 65.8 (A^2) for my lipid bilayer,
I need my area to be 3.03338 *10^4 (A^2). after 0.5 ns, I have 2.85*10^4 (A^2)
area. after 1ns, it will go below 2.8*10^4 (A^2) area. my area keeps
shrinking by time.
should I stop the shrinking at nearest possible experimental value and
continue my simulation at fixed area suing useConstantArea =yes or is it
something wrong with my configuration file.
the membrane tutorial stated that "Simulation of membrane bilayers have
shown that the current CHARMM forcefield parameters do no reproduce the
experimentally observed area per lipid over long MD trajectories. During
the previous simulation steps, we let the are in the xy-plane fluctuate to
permit packing of lipids against the protein. However, after good packing
has been observed, one should keep the area in the xy-plane constant."

Thank you in advance,

Yousef Nademi

PhD student, Chemical and Material Department, University of Alberta, Canada

set fs 748000
#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################
set X 183.0400047302246
set Y 183.06000518798828
set Z 79.99099731445313
set CX 41.34382247924805
set CY 42.519981384277344
set CZ 0.03781001642346382

set PX 192
set PY 192
set PZ 80

structure ./ionized.psf
coordinates ./ionized.pdb
extendedsystem ./popc_melting_lipid_tail-02-498000.restart.xsc
bincoordinates ./popc_melting_lipid_tail-02-498000.restart.coor
binvelocities ./popc_melting_lipid_tail-02-498000.restart.vel

#set consfileName ./restrain_ionized_exceptlipidtail.pdb
set temperature 310
set outputname popc_melting_lipid_tail-02-748000

firsttimestep $fs

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters ./par_all36_prot_na_lipid_pei_heparin-Jan2016.prm
#temperature 0

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 12
switching on ;# to be checked the result w or w/o switching
switchdist 10
pairlistdist 12

# Integrator Parameters
timestep 2.0 ;# 1fs/step
rigidBonds all ;# needed for 2fs steps "'shake bonh' in CHARMM?"
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10 ;# 'inbfrq 10 in CHARMM?'

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

# Periodic Boundary Conditions
cellBasisVector1 $X 0. 0.
cellBasisVector2 0. $Y 0.
cellBasisVector3 0. 0 $Z
cellOrigin $CX $CY $CZ

wrapAll on
wrapWater on
wrapNearest off

margin 3

# PME (for full-system periodic electrostatics)
if {1} {
PME yes
# PMEGridSpacing 1.0
PMEGridSizeX $PX
PMEGridSizeY $PY
PMEGridSizeZ $PZ
}

# Constant Pressure Control (variable volume)
if {1} {
useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell yes ;# no for water box, yes for membranees
useConstantArea no ;# no for water box, maybe for membrane

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

# Output
outputName $outputname
binaryoutput off
restartfreq 1000 ;#1000 steps= every 2 ps
dcdfreq 1000
xstFreq 1000
outputEnergies 50
outputPressure 50

#############################################################
## EXTRA PARAMETERS ##
#############################################################
if {0} {
fixedAtoms on
fixedAtomsForces on
fixedAtomsFile $consfileName
fixedAtomsCol B
}

tclforces on
set waterCheckFreq 100
set lipidCheckFreq 100
set allatompdb ./ionized.pdb
tclForcesScript keep_water_out.tcl

#############################################################
## EXECUTION SCRIPT ##
#############################################################

run 250000 ;# 0.5 ns

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:21:26 CST