From ee45a5ba2e770229ce15aaf38b6e891b936bf959 Mon Sep 17 00:00:00 2001 From: Serguiei Vassiliev Date: Mon, 6 Nov 2023 17:27:13 -0400 Subject: [PATCH] bilayer simulation --- .../17-Preparing_and_Simulating_Bilayer.md | 52 +++++++++++++++---- 1 file changed, 42 insertions(+), 10 deletions(-) diff --git a/_episodes/17-Preparing_and_Simulating_Bilayer.md b/_episodes/17-Preparing_and_Simulating_Bilayer.md index 607bb19..c723ea3 100644 --- a/_episodes/17-Preparing_and_Simulating_Bilayer.md +++ b/_episodes/17-Preparing_and_Simulating_Bilayer.md @@ -18,7 +18,7 @@ keypoints: This in-depth tutorial guides you through the process of building a bilayer, preparing simulation input files, minimizing energy, equilibrating the system, and running an equilibrium molecular dynamics simulation. You will need to follow a number of steps to complete the tutorial. If you are already familiar with some of these topics, you can skip them and focus on the ones you don't know about. -### Step 1: A bilayer of lipids is formed by packing them together. +### 1. A bilayer of lipids is formed by packing them together. It has been shown that the composition of the phospholipids in the erythrocyte membrane is asymmetric. Here is an article that reviews this topic in detail: [Interleaflet Coupling, Pinning, and Leaflet Asymmetry—Major Players in Plasma Membrane Nanodomain Formation](https://www.frontiersin.org/articles/10.3389/fcell.2016.00155/full). @@ -43,14 +43,46 @@ packmol-memgen \ ~~~ {: .language-bash} + The job takes about 30 minutes. It generates several output files: - AMBER topology bilayer_only_lipids.top - AMBER coordinates bilayer_only_lipids.crd - PDB-formatted system bilayer_only_lipids.pdb -### Step 2: Creating a system with different force filed parameters. +By default parametrization uses LIPID17, ff14SB and TIP3P force fields. + +Figure: +Neutral: POPC green, POPE blue, PSM white +Charged: POPS red + +### 2. Preparing the bilayer for simulation with a different force field. + +The file bilayer_only_lipids.pdb can be used to create topology files using any AMBER force field distributed with AmberTools. This allows to prepare a simulation system using a force field that is not provided with packmol-memgen. + +However, if your system contains SPM this file will not work because packmol-memgen erroneoiusly adds TER record between the lipid headgroup (SPM) and second acyl tail (SA). This splits PSM molecule in 2 parts: PA+SPM and SA. + +For this reason we will, create a pdb file from AMBER topology using parmed and then use tleap to reparameterize it using a polarizable water model OPC3-Pol which is not yet available in packmol-memgen natively. + +~~~ +module load StdEnv/2020 gcc/9.3.0 openmpi/4.0.3 ambertools/23 +python +~~~ +{: .language-bash} + +~~~ +import parmed +amber = parmed.load_file("bilayer_only_lipid.top", "bilayer_only_lipid.crd") +amber.save("bilayer_only_lipid_fixed.pdb") +quit() +~~~ + +Alternatively you can fix the file bilayer_only_lipid.pdb manually using standard unix commands grep and awk. The following command removes these TER records: +~~~ +grep -n "H22 SPM " bilayer_only_lipid.pdb | cut -f1 -d: | awk '{print $0 + 1}' | awk -F, 'NR==FNR { nums[$0]; next } !(FNR in nums)' - bilayer_only_lipid.pdb > bilayer_only_lipid_fixed.pdb +~~~ +{: .language-bash} -There is a file named bilayer_only_lipids.pdb that can be used to create topology files using AMBER force fields that have not been included in packmol-memgen yet. We only need to know one more thing before we can proceed: the size of the periodic box. It can be extracted from the file bilayer_only_lipid.top: +We only need to know one more thing before we can proceed: the size of the periodic box. It can be extracted from the file bilayer_only_lipid.top: ~~~ grep -A2 BOX bilayer_only_lipid.top @@ -60,7 +92,7 @@ grep -A2 BOX bilayer_only_lipid.top ~~~ ... %FORMAT (5E168) - 9.00000000E+01 1.03537600E+02 1.03536600E+02 9.99353000E+01 + 9.00000000E+01 1.03536600E+02 1.03744300E+02 9.97203000E+01 ~~~ {: .output} @@ -71,8 +103,8 @@ module load StdEnv/2020 gcc/9.3.0 openmpi/4.0.3 ambertools/23 tleap -f <(cat << EOF source leaprc.water.opc3pol source leaprc.lipid21 -sys=loadpdb bilayer_only_lipid.pdb -set sys box {103.5 103.5 99.9} +sys=loadpdb bilayer_only_lipid_fixed.pdb +set sys box {103.5 103.7 99.7} saveamberparm sys bilayer.parm7 bilayer.rst7 quit EOF) @@ -81,7 +113,7 @@ EOF) Don't forget to replace {103.5 103.5 99.9} with your own numbers. -### Step 3: Energy minimization +### 3. Energy minimization A system must be optimized in order to eliminate clashes and prepare it for molecular dynamics. To run energy minimization we need three files: @@ -121,7 +153,7 @@ srun sander.MPI -O -i min.in \ When the optimization job is completed, the file minimized.rst7 will contain the optimized coordinates that were produced. -### Step 4: Heating +### 4. Heating ~~~ Heating &cntrl @@ -165,7 +197,7 @@ pmemd.cuda -O -i heating.in \ ~~~ {: .language-bash} -### Step 5: The first phase of the equilibration process. +### 5. The first phase of the equilibration process. ~~~ Equilibration 1 &cntrl @@ -209,7 +241,7 @@ pmemd.cuda -O -i equilibration.in \ ~~~ {: .language-bash} -### Step 6: The second phase of the equilibration process. +### 6. The second phase of the equilibration process. ~~~ Equilibration 2 &cntrl