6.1 POPC lipid bilayers

One of the widely used model for biological membranes is a bilayer of POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) lipids. The chemical formula of POPC is shown below. It consists of hydrophobic tail groups and hydrophilic head groups. In POPC, one of the fatty acid tail has a double bond with a cis conformation, and there is a stereocenter in the middle (both indicated in red in the figure). Because of its amphiphilic character, POPC lipids self-assemble to form a bilayer structure.

The aim of this tutorial is to learn,

  1. How to prepare an all-atom model for a POPC bilayer
  2. How to equilibrate the POPC bilayer system.
  3. How to analyze the trajectory and what we learn about the membrane.

In this tutorial, we carry out all-atom MD simulations of a lipid bilayer with GENESIS and CHARMM-GUI using the CHARMM36 force fields [1].


Let’s download the tutorial file (tutorial-6.1.tar.gz ). This tutorial consists of four steps: 1) system setup using CHARMM-GUI, 2) energy minimization and equilibration, 3) production run, and 4) trajectory analysis. Since we use the CHARMM force field, we make a symbolic link to the CHARMM toppar directory (see Tutorial 2).

# Download the tutorial file
$ cd /home/user/GENESIS/Tutorials
$ mv ~/Downloads/tutorial-6.1.tar.gz ./
$ tar -xzvf tutorial-6.1.tar.gz
$ cd tutorial-6.1

# Make a symbolic link to the CHARMM toppar directory
$ ln -s ../../Others/CHARMM/toppar ./
$ ls
1_charmm-gui  2_min_equil  3_production  4_analysis  toppar

1. Setup initial structure with CHARMM-GUI

One of the difficulties is a setup of a bilayer structure. Because each lipid molecule has a random structure, it takes much time to equilibrate the whole system, if one starts the simulation from an ordered structure. CHARMM-GUI is one of the useful tools to generate initial structures of a lipid bilayer for MD simulations [2]. It provides an automatic membrane builder not only for a lipid bilayer but also for protein-membrane systems.

In case the server is down, you can skip this section and use the files included in a folder 1_charmm-gui, which we obtained from CHARMM-GUI.

Click [here] to access the Bilayer builder to build a POPC lipid bilayer. (It can be accessed from the top page, Input Generator, Membrane Builder, and then Bilayer Builder.) At the first step of the bilayer builder, we choose “Membrane Only System”:

tutorial1.2-step1In the next step, we decide the number of lipids and water molecules to be added in a system. Here, we put 40 POPC and 40 POPC lipids in the upper and lower leaflets, respectively (Note that these numbers might be a lower limit for stable MD simulation with spdyn). Because each lipid molecule has a certain area per lipid (Aexp = 68.3Å2 for POPC), the initial box size (X, Y) = (52.27Å, 52.27Å) can be estimated from sqrt{(number of lipids in one leaflet)*(area per lipid)}. tutorial1.2-step2

Note that one may alternatively choose “Ratios of lipid components”, set  “Length of X and Y”, and set “Upperleaflet Ratio” and “Lowerleaflet Ratio” for each lipids. This scheme is useful for building a mixed lipid system.

In the next step, we add 0.15 M KCl in the system, mimicking the cellular environment, with the Monte Carlo method for ion placing. Then, we click “Next Step” and continue to Step5. We check “GENESIS” in “Input Generation Options” and select the NPT ensemble at T = 303.15 K in “Equilibration Options”. At this temperature, the POPC bilayer is in the liquid phase rather than the gel phase (Tm = 271 K). After clicking the “Next Step”, we can download a set of input files for minimization, equilibration, and production runs at the final step of the membrane builder. The download begins when you click a “download.tgz” button in the upper right corner.


After the download, let us check the structure of the downloaded PDB file using VMD:

# move the downloaded tar file and untar
$ mv ~/Download/charmm-gui.tgz ./
$ tar -xzf charmm-gui.tgz
$ ls
1_charmm-gui/ 3_production/ charmm-gui-XXXXXXXXXXXX/
2_min_equil/ 4_analysis/ charmm-gui.tgz

# Check the final pdb file
$ cd charmm-gui-XXXXXXXXXXXX
$ vmd step5_assembly.pdb


“XXXXXXXXXXXX” is a 12-digit ID number given by the CHARMM-GUI sever. Visually check that the system is constructed as you intended (the type of lipids, number of lipids, etc.). Note that the membrane normal is taken to be the z-axis, which is the convention in the field. You may also check step5_assembly.str, where the box size, the number of water molecules, and other information of the system are found. If no serious problem is found, we use the PDB file and an associated PSF file (step5_assembly.psf) as the initial coordinates of the simulation. Rename the downloaded folder to 1_charmm-gui,

# Save the original folder, and rename the new one 
$ cd ..
$ mv 1_charmm-gui 1_charmm-gui.org
$ mv charmm-gui-XXXXXXXXXXXX 1_charmm-gui
$ ls
1_charmm-gui/ 1_charmm-gui.org/ 3_production/ 
2_min_equil/ 4_analysis/ charmm-gui.tgz

Control files for GENESIS are contained in 1_charmm-gui/genesis:

$ ls 1_charmm-gui/genesis
membrane_restraint.genesis.str  step6.0_minimization.inp  step6.5_equilibration.inp
membrane_restraint2.genesis.str step6.1_equilibration.inp step6.6_equilibration.inp
restraints/                     step6.2_equilibration.inp step7_production.inp
step5_charmm2genesis.inp        step6.3_equilibration.inp
step5_charmm2genesis.str        step6.4_equilibration.inp

These files (but slightly modified) are used in the following sections.

2. Minimization and Equilibration

Let’s move to 2_min_equil and find control files for minimization (step6.0) and equilibration (step6.1 – 6.6),

$ cd ../../2_min_equil
$ ls
equil.vmd        restraints/                step6.1_equilibration.inp  step6.4_equilibration.inp
equil_area.gpi   run.sh                     step6.2_equilibration.inp  step6.5_equilibration.inp
equil_temp.gpi   step6.0_minimization.inp   step6.3_equilibration.inp  step6.6_equilibration.inp

Note that psffile, pdbfile, and reffile are set to the files in 1_charmm-gui,

psffile = ../1_charmm-gui/step5_assembly.psf  # protein structure file
pdbfile = ../1_charmm-gui/step5_assembly.pdb  # PDB file
reffile = ../1_charmm-gui/step5_assembly.pdb  # reference PDB file for restraint

The control files for minimization and equilibration are almost the same as before [see Sec. 3.2]. The difference is in the way restraints are applied.

2.1. The restraints for POPC

The dihedral angles, C3-C1-C2-O21 and C28-C29-C210-C211 (shown in yellow in the figure below), are restrained to 120 and 0 degree, respectively, to keep the stereoisomer and the cis-isomer. Furthermore, the position of a phosphate atom, (shown in purple), is also restrained in the Z direction to keep the thickness of the bilayer.

The dihedral angles are restrained using localresfile in the [INPUT] section. For example, step6.0_minimization.inp is as follow,

topfile = ../toppar/top_all36_lipid.rtf             # topology file
parfile = ../toppar/par_all36_lipid.prm             # parameter file
strfile = ../toppar/toppar_water_ions.str           # stream file
psffile = ../1_charmm-gui/step5_assembly.psf        # protein structure file
pdbfile = ../1_charmm-gui/step5_assembly.pdb        # PDB file
reffile = ../1_charmm-gui/step5_assembly.pdb        # reference PDB file for restraint
localresfile = restraints/step6.0_minimization.dihe # local restraint file

localresfile lists the dihedral angles in the following format,

DIHEDRAL 60 63 65 67 250 0.0
DIHEDRAL 36 25 28 30 250 120.0
DIHEDRAL 194 197 199 201 250 0.0

The first line means that a dihedral angle of atom index, 60-63-65-67, (i.e. C28-C29-C210-C211 of the first POPC) is restrained to 0.0 degree with a force constant of 250 kcal/mol. Similarly, the second line restrains C3-C1-C2-O21 of the first POPC to 120.0 degree with the same force constant. The restraints for other POPCs follow in the file.

The positional restraint of phosphate atoms is specified in the control file by [RESTRAINTS] and [SELECTION] sections,

group1           = sid:MEMB and ((rnam:POPC and (an:P)))  # restriant group 1

nfunctions       = 1       # number of functions
function1        = POSI    # restraint function type
direction1       = Z       # direction [ALL,X,Y,Z]
constant1        = 2.5     # force constant
select_index1    = 1       # restrained groups

Note that the direction is set to the Z direction. Thus, POPC can move freely in the membrane plane, and restrained in the direction of the membrane normal.

2.2. Run GENESIS

The setting of each step is summarized in the following Table.

Step Ensemble Integrator Step dt / fs Time / ps k [dihed] k [Pz]
0 Min 2,000 250 2.5
1 NVT VVER 25,000 1 25 250 2.5
2 NVT VVER 25,000 1 25 100 2.5
3 NPT VVER 25,000 1 25 50 1.0
4 NPT VVER 50,000 2 100 50 0.5
5 NPT VVER 50,000 2 100 25 0.1
6 NPT VRES 40,000 2.5 100 0 0

Note that the force constants (k [dihed] and k[Pz]) are gradually reduced. Also, the timestep, which is initially set to 1 fs, is increased to 2 fs and 2.5 fs in the last step. Step 6 is carried out with the RESPA integrator.

The simulation with NPT is carried out with semi-iso option, which keeps the ratio of the box size in the x and y directions,

ensemble         = NPT      # [NVE,NVT,NPT]
tpcontrol        = BUSSI    # thermostat and barostat
temperature      = 303.15   # initial and target temperature (K)
pressure         = 1.0      # atm
isotropy         = SEMI-ISO # keep the ratio of X and Y

Let’s execute GENESIS. run.sh is a script to run SPDYN for all steps:

$ cat run.sh


# 2000 step energy minimization
# with dihedral angle restraint (k = 250 kcal/mol)
#      z-positional restraint on P atoms (k = 2.5 kcal/mol/Angs^2)
mpirun -np 16 $spdyn step6.0_minimization.inp >& log6.0

# 25ps MD in the NVT ensemble
# with dihedral angle restraint (k = 250 kcal/mol)
#      z-positional restraint on P atoms (k = 2.5 kcal/mol/Angs^2)
mpirun -np 16 $spdyn step6.1_equilibration.inp >& log6.1

# 25ps MD in the NVT ensemble
# with dihedral angle restraint (k = 100 kcal/mol)
#      z-positional restraint on P atoms (k = 2.5 kcal/mol/Angs^2)
mpirun -np 16 $spdyn step6.2_equilibration.inp >& log6.2

# 25ps MD in the NPT ensemble
# with dihedral angle restraint (k = 50 kcal/mol)
#      z-positional restraint on P atoms (k = 1.0 kcal/mol/Angs^2)
mpirun -np 16 $spdyn step6.3_equilibration.inp >& log6.3

# 100ps MD in the NPT ensemble
# with dihedral angle restraint (k = 50 kcal/mol)
#      z-positional restraint on P atoms (k = 0.5 kcal/mol/Angs^2)
mpirun -np 16 $spdyn step6.4_equilibration.inp >& log6.4

# 100ps MD in the NPT ensemble
# with dihedral angle restraint (k = 25 kcal/mol)
#      z-positional restraint on P atoms (k = 0.1 kcal/mol/Angs^2)
mpirun -np 16 $spdyn step6.5_equilibration.inp >& log6.5

# 100ps MD in the NPT ensemble without any restraints using VRES
mpirun -np 16 $spdyn step6.6_equilibration.inp >& log6.6

Execute this script. This may take about an hour or so.

# Run all steps
$ ./run.sh

When the calculation is done, it is a good practice to check the status before moving on to the production step.

2.3. Check the status: Visualization
$ vmd -e equil.vmd

The molecules can be wrapped into a simluation box by typing,

pbc wrap -center com -compound residue -all

in the console of VMD. The console of VMD opens by choosing “Extensions -> Tk Console” in the VMD Main window.

2.4. Check the status: Temperature and Area per lipid (AL)

Extract the “INFO” lines from output files to separate files,

$ for i in 1 2 3 4 5 6; do
grep INFO log6.${i} > energy6.${i}.log

Then, plot the data of temperature using gnuplot,

$ cat equil_temp.gpi
set xrange [0:375]
set yrange [285:325]

plot 'energy6.1.log' using 3:16 w l t 'equil6.1', \
     'energy6.2.log' using ($3+25):16 w l t 'equil6.2', \
     'energy6.3.log' using ($3+50):16 w l t 'equil6.3', \
     'energy6.4.log' using ($3+75):16 w l t 'equil6.4', \
     'energy6.5.log' using ($3+175):16 w l t 'equil6.5', \
     'energy6.6.log' using ($3+275):16 w l t 'equil6.6'

$ gnuplot equil_temp.gpi

AL is obtained by Lx * Ly / nL, where Lx and Ly are the simulation box size in x and y dimension, respectively, and nL is the number of lipids in one of the leaflet.

$cat equil_area.gpi
set xrange [50:375]
set yrange [60:70]

plot 'energy6.3.log' using ($3+50):($18*$19/40) w l t 'equil6.3', \
'energy6.4.log' using ($3+75):($18*$19/40) w l t 'equil6.4',  \
'energy6.5.log' using ($3+175):($18*$19/40) w l t 'equil6.5', \
'energy6.6.log' using ($3+275):($18*$19/40) w l t 'equil6.6'

$gnuplot equil_area.gpi

If no serious problem is found, we now continue to a production step.

3. Production

Let’s move to 3_production. In this section, we carry out a NPT simulation for 1 ns. run.sh is a script to run SPDYN:

$cd ../3_production
$ cat run.sh


# 1000ps MD in the NPT ensemble using VRES
mpirun -np 16 $spdyn step7_production.inp >& log7

Although the script uses 1 OpenMP thread x 16 MPI processes, it may not be optimal for your computing system. It is highly recommended to benchmark the cost before starting the production run. See Sec. 3.3 for details.

Execute this script. This may take several hours.

# Run all steps
$ ./run.sh

The trajectory can be visualized by,

$ vmd -e prod.vmd

4. Analysis

Let’s analyze the resulting trajectory. Move to 4_analysis folder,

$ cd ../4_analysis
$ ls
prod_area.gpi  prod_temp.gpi  thickness.gpi  thickness.inp
4.1. Area per lipid (AL)

AL is obtained as before by first extracting the “INFO” lines from the output, and calculating and plotting Lx * Ly / nL with gnuplot,

$ grep INFO ../3_production/log7 > energy7.log
$ cat prod_area.gpi
set xlabel 'Time (ps)'
set ylabel 'Area per lipid (\305^2)'

set xrange [0:1000]
set yrange [60:70]

plot 'energy7.log' using 3:($17*$18/40) w l t 'equil7'

$ gnuplot prod_area.gpi

Although the calculated AL is slightly smaller than the experimental value (Aexp = 68.3Å2 for POPC), GENESIS reproduces the averaged area per lipid of 64.7Å2 reported in the CHARMM36 paper [1].

4.2. Bilayer thickness

lipidthick_analysis calculates the average Z-position of specified atoms in the upper and lower leaflets, and gives the membrane thickness as a different between the two. The input file, thickness.inp, looks as follow,

psffile = ../1_charmm-gui/step5_assembly.psf
reffile = ../1_charmm-gui/step5_assembly.pdb

outfile = thickness.log

trjfile1      = ../3_production/step7_production.dcd
md_step1      = 400000
mdout_period1 = 800
ana_period1   = 800
trj_format    = DCD
trj_type      = COOR+BOX

group1        = resname:POPC and an:P

check_only    = NO
membrane_atom = 1     # atom group representing lipid bilayer
  • output file is set to thickness.log
  • membrane_atom = 1 of [OPTION] and group1 of [SELECTION] specify phosphate atoms of POPC lipid molecules to be analyzed.

The analysis program is invoked by the following command.

$ ../../bin/lipidthick_analysis thickness.inp >& thickness.out
$ cat thickness.log
     1    38.718    19.377   -19.341
     2    38.708    19.453   -19.255
     3    38.815    19.454   -19.361
     500  39.758    19.947   -19.811

The numbers in each column are the step number, the membrane thickness, the average Z-postion of the upper and lower leaflets, respectively. These data are plotted using gnuplot by the following script.

$ cat thickness.gpi
set xlabel 'Time (ps)'
set ylabel 'P-P distance (Angs)'

set xrange [0:1000]
set yrange [38.2:40.5]

plot "thickness.log" u ($1*2):2 w l t "POPC"

$ gnuplot prod_area.gpi


  1. J. B. Klauda et al., J. Phys. Chem., 114, 7830-7843 (2010).
  2. S. Jo, T. Kim, and W. Im, PLoS ONE, 2, e880 (2007).

Written by Takaharu Mori@RIKEN Theoretical molecular science laboratory
June, 30, 2016

Written by Kiyoshi Yagi @ RIKEN, CPR, Theoretical molecular science laboratory
October, 7, 2019