Skip to content

Commit

Permalink
membrane simulation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ssvassiliev
ssvassiliev committed Nov 10, 2023
1 parent 01dce80 commit 195aed7
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 512 deletions.
326 changes: 0 additions & 326 deletions _episodes/14-Preparing_Protein-Membrane_System.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,332 +75,6 @@ packmol-memgen \
~~~
{: .language-bash}

This generates bilayer_only_lipids.pdb which can be used to create topology files with AMBER force fields not available in packmol-memgen. Example using OPC3-Pol polarizable water:

~~~
leap_input=$(cat << EOF
source leaprc.water.opc3pol
source leaprc.lipid21
sys=loadpdb bilayer_only_lipid.pdb
set sys box {111.3 111.3 81.5}
saveamberparm sys bilayer.parm7 bilayer.rst7
quit
EOF)
tleap -f <(echo "$leap_input")
~~~
{: .language-bash}

Box size can be extracted from the file bilayer_only_lipid.top:

~~~
grep -A2 BOX bilayer_only_lipid.top
~~~
{: .language-bash}


### Simulating a packed membrane system.

#### Energy minimization
Minimization input file (minimization.in).
~~~
Minimization
&cntrl
imin=1, ! Minimization
ntmin=3, ! Steepest Descent
maxcyc=1000, ! Maximum number of cycles for minimization
ntpr=10, ! Print to mdout every ntpr steps
ntwr=1000, ! Write a restart file every ntwr steps
cut=10.0, ! Nonbonded cutoff
/
~~~
{: .file-content}

Submission script
~~~
#!/bin/bash
#SBATCH --mem-per-cpu=4000 --ntasks=10 --time=6:0:0
module purge
module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
srun sander.MPI -O -i min.in \
-o minimization.out \
-p ../bilayer.parm7 \
-c ../bilayer.rst7 \
-r minimized.rst7
~~~
{: .language-bash}

#### Heating
~~~
Heating
&cntrl
imin=0, ! Molecular dynamics
ntx=1, ! Read coordinates with no initial velocities
ntc=2, ! SHAKE on for bonds with hydrogen
ntf=2, ! No force evaluation for bonds with hydrogen
tol=0.0000001, ! SHAKE tolerance
nstlim=20000, ! Number of MD steps
ntt=1, ! Berendsen thermostat
tempi=100, ! Initial temperature
temp0=300, ! Target temperature
tautp=5.0, ! Time constant for heat bath coupling
ntpr=100, ! Write energy info every ntwr steps
ntwr=10000, ! Write restart file every ntwr steps
ntwx=2000, ! Write to trajectory file every ntwx steps
dt=0.001, ! Timestep
ntb=2, ! Constant pressure periodic bc.
barostat=1, ! Berendsen barostat
taup=5, ! Pressure relaxation time
ntp=3, ! Semiisotropic pressure scaling
csurften=3, ! Constant surface tension with interfaces in the xy plane
cut=10.0, ! Nonbonded cutoff
/
~~~
{: .file-content}

Submission script
~~~
#!/bin/bash
#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=6:0:0
module purge
module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
pmemd.cuda -O -i heating.in \
-o heating.out \
-p ../bilayer.parm7 \
-c minimized.rst7 \
-r heated.rst7 \
~~~
{: .language-bash}

#### Equilibration - phase 1
~~~
Equilibration 1
&cntrl
imin=0, ! Molecular dynamics
irest=1, ! Restart
ntx=5, ! Read positions and velocities
ntc=2, ! SHAKE on for bonds with hydrogen
ntf=2, ! No force evaluation for bonds with hydrogen
tol=0.0000001, ! SHAKE tolerance
nstlim=100000, ! Number of MD steps
ntt=1, ! Berendsen thermostat
temp0=300, ! Set temperature
tautp=5.0, ! Time constant for heat bath coupling
ntpr=100, ! Write energy info every ntwr steps
ntwr=10000, ! Write restart file every ntwr steps
ntwx=2000, ! Write to trajectory file every ntwx steps
dt=0.001, ! Timestep (ps)
ntb=2, ! Constant pressure periodic bc.
barostat=1, ! Berendsen barostat
taup=5, ! Pressure relaxation time
ntp=3, ! Semiisotropic pressure scaling
csurften=3, ! Constant surface tension with interfaces in the xy plane
cut=10.0, ! Nonbonded cutoff
/
~~~
{: .file-content}

Submission script
~~~
#!/bin/bash
#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=6:0:0
module purge
module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
pmemd.cuda -O -i equilibration.in \
-o eq-1.out \
-p ../bilayer.parm7 \
-c heated.rst7 \
-r eq-1.rst7 \
~~~
{: .language-bash}

#### Equilibration - phase 2
~~~
Equilibration 2
&cntrl
imin=0, ! Molecular dynamics
irest=1 ! Restart
ntx=5, ! Read positions and velocities
ntc=2, ! SHAKE on for bonds with hydrogen
ntf=2, ! No force evaluation for bonds with hydrogen
tol=0.0000001, ! SHAKE tolerance
nstlim=2500000, ! Number of MD steps
ntt=3, ! Langevin thermostat
temp0=300, ! Target temperature
gamma_ln=2.0, ! Collision frquency
ntpr=100, ! Write energy info every ntwr steps
ntwr=10000, ! Write restart file every ntwr steps
ntwx=10000, ! Write to trajectory file every ntwx steps
dt=0.002, ! Timestep (ps)
ntb=2, ! Constant pressure periodic bc.
barostat=1, ! Berendsen barostat
taup=1.0, ! Pressure relaxation time
ntp=3, ! Semiisotropic pressure scaling
csurften=3, ! Constant surface tension with interfaces in the xy plane
cut=10.0, ! Nonbonded cutoff
/
~~~
{: .file-content}

Submission script
~~~
#!/bin/bash
#SBATCH --mem-per-cpu=4000 --gpus-per-node=v100:1 --time=6:0:0
module purge
module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 amber/20.12-20.15
pmemd.cuda -O -i equilibration-2.in \
-o eq-2.out \
-p ../bilayer.parm7 \
-c eq-1.rst7 \
-r eq-2.rst7 \
~~~
{: .language-bash}

AMBER offers only two barostats: Monte Carlo and Berendsen. As it was found that Monte Carlo barostat results in the depression of areas per lipid, it is not recommended for lipid bilayer simulations with AMBER [[Lipid21: Complex Lipid Membrane Simulations with AMBER]](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9007451/). Berendsen barostat does not strictly sample from the isothermal-isobaric ensemble and typically yields volume fluctuations that are too low.
If you need correct thermodynamocs ensemble you can continue simulating this system in another MD package that offers more pressure coupling algorithms. In the next section we will demonstrate transferring equilibrated MD system from AMBER to GROMACS

#### Moving to GROMACS

- Change directory to 4-equilibration.
- Create a directory for GROMACS simulation
- Load ambertools and gromacs modules and start python

~~~
mkdir ../5-amb2gro
cp eq-2.rst7 ../5-amb2gro
cd ../5-amb2gro
ln -s ../bilayer.parm7 bilayer.parm7
module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 ambertools/23 gromacs/2023.2
python
~~~
{: .language-bash}


#### Converting amber restart and topology to gromacs.
Use parmed AMBER utility to convert AMBER topology and restart files to GROMACS.

~~~
import parmed as pmd
amber = pmd.load_file("bilayer.parm7", "eq-2.rst7")
ncrest=pmd.amber.Rst7("eq-2.rst7")
amber.velocities=ncrest.vels
amber.save("restart.gro")
amber.save("bilayer.top")
quit()
~~~
{: .language-python}

These parmed commands will create two GROMACS files: restart.gro and bilayer.top

#### Creating and editing index file
Next we create and edit index file.
~~~
gmx make_ndx -f restart.gro << EOF
del 12
del 12
1&!rNa+&!rCl-
name 12 Lipids
11|rNa+|rCl-
name 13 WaterIons
q
EOF
~~~
{: .language-bash}

Check created index file:
~~~
gmx make_ndx -f restart.gro -n index.ndx
~~~
{: .language-bash}

You should have two new groups:
12 Lipids : 53850 atoms
13 WaterIons : 66792 atoms

#### Creating GROMACS restart file

Gromacs simulation input file

tcoupl = nose-hoover
tau-t = 1
tc-grps = 11 12
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau-p = 5 5
compressibility = 4.5e–5 4.5e–5


~~~
Title = bilayer
; Run parameters
integrator = md
nsteps = 400000
dt = 0.001
; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstenergy = 100
nstlog = 10000
nstxout-compressed = 5000
compressed-x-grps = System
; Bond parameters
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid
nstlist = 10
rcoulomb = 1.0
rvdw = 1.0
DispCorr = Ener ; anaytic VDW correction
; Electrostatics
coulombtype = PME
pme_order = 4
; Temperature coupling is on
tcoupl = V-rescale
tc-grps = WaterIons Lipids
tau_t = 0.1 0.1
ref_t = 300 300
; Pressure coupling is on
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 5.0
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
; Periodic boundary conditions
pbc = xyz
; Velocity generation
gen_vel = no
~~~

~~~
gmx trjconv -f restart.gro -o restart.trr
gmx grompp -p bilayer.top -c restart.gro -t restart.trr -f production.mdp -maxwarn 1
~~~
{: .language-bash}

#### Running simulation
~~~
#!/bin/bash
#SBATCH --mem-per-cpu 4000 --time 1:0:0
#SBATCH -c10 --gpus-per-node=v100:1
module load StdEnv/2020 gcc/9.3.0 cuda/11.4 openmpi/4.0.3 gromacs/2023.2
srun gmx mdrun -ntomp ${SLURM_CPUS_PER_TASK:-1} \
-nb gpu -pme gpu -update gpu -bonded cpu -s topol.tpr
~~~
{: .language-bash}


### Embedding a protein into a bilayer
We will use PDB file 6U9P (wild-type MthK pore in ~150 mM K+) for this exercise.
Expand Down
Loading

0 comments on commit 195aed7

Please sign in to comment.