GENESIS Tutorial 4
String Method
Lecturer: Yasuhiro Matsunaga
TA: Jaewoon Jung, Chigusa Kobayashi,
Koichi Tamura, Motoshi Kamiya
2017/03/01
1

Contents
I. Most probable pathway with the string method
ü Equilibration along an initial pathwayü String method
ü Visualization of the pathways
II. Free energy profile along the pathwayü Equilibration along the pathway
ü Umbrella sampling along the pathway
ü Analysis of free energy profile2

Target system
Conformational change of alanine-tripeptide (same as REMD’s)

How can we characterize this conformational change
by MD simulation?
3

Two main results obtained by this tutorial
Pathways (=most probable pathway)
in the dihedral angle space

Gets insights into causalities

Free energy profile along the pathway

Identifies transition state4

Basic usage# enter the login node
$ ssh –Y -l userXX XXXXXX# enter a computational node
$ ssh –Y XXX# change the current directory to Tutorial_4$ cd ~/Tutorial_4/5

Flow of string method simulations
Converged
pathway

Initial pathway
Single structure

Converged
pathway

Image 1

Image 2
Pulling thesystem to each
image

VisualizeUmbrellasampling…

…

…

Image 8String method
Move images

Free energy
6Preparation$ cd ~/Tutorial_4/$ ls
0_setup
1_minimize

2_pre-equilibration
3_equilibration

4_string_method
5_visualize

6_umbrella
7_analysis

We will skip step 0 to 3 because these are same as REMD’s
We will change the directory from 3_initial_path to 6_visualization step by step7

Step3. Pulling the system toward images
along the initial pathway
Single structure Pulling

Initial pathway
Image 1
Image 2
Image 3
Image 4
Image 5
Image 6
Image 7
Image 88Control file for GENESIS spdyn
The RPATH module in GENESIS can simulate replicas
without string method performed (rpath_period = 0).
[INPUT]
rstfile = ../2_minimize_pre-equi/step2.4.rst[OUTPUT]
logfiledcdfilerstfilerpathfile===={}.log{}.dcd{}.rst{}.rpath[ENERGY]
forcefield
electrostatic
switchdist
cutoffdist
pairlistdist= CHARMM= PME= 10.0= 12.0= 13.5[DYNAMICS]
integrator
nstepstimestepeneout_period
crdout_period
rstout_period= LEAP= 2000= 0.002= 1000= 1000= 2000[CONSTRAINTS]rigid_bond= YES[ENSEMBLE]
ensemble
tpcontrol
temperaturepressure= NPT= LANGEVIN= 300.00= 1.0

Simulate replicas
without string method[RPATH]
nreplica= 8
rpath_period = 0
rest_function = 1 2

Replica index is
automaticallyinserted into {}of the output
filename.

[SELECTION]
group1group2group3group4group5= atomindex:15= atomindex:17= atomindex:19= atomindex:25= atomindex:27[RESTRAINTS]
nfunctions= 2

function1= DIHEDconstant1= 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
reference1= -40 -40 -40 -40 -40 -40 -40 -40
select_index1 = 1 2 3 4 # PHIImages along the initial# use SHAKE/SETTLE

pahtwayfunction2= DIHED
constant2= 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
reference2= -45.0 -19.3 6.4 32.1 57.9 83.6 109.3 135.0
select_index2 = 2 3 4 5 # PSI

9Running a job$ cd 3_initial_path/$ lsrun.inp
run.sh3 threads × 8 MPI procs = 24 CPU cores$ cat run.sh#!/bin/bash
export OMP_NUM_THREADS=3
source /home2/data/env/setpath.CPU.DPmpirun -np 8 spdyn run.inp >run.out$ ./run.shThis will finish in about 2 minutes
10Step4. String method

Convergedpathway(most proablepathway)

Initialpathway

11

Control file for GENESIS spdyn
The RPATH module updates images every rpath_period = 500[INPUT]
rstfile = ../3_initial_path/{}.rst[OUTPUT]
logfiledcdfilerstfilerpathfile===={}.log{}.dcd{}.rst{}.rpath[ENERGY]
forcefield
electrostatic
switchdist
cutoffdist
pairlistdist= CHARMM= PME= 10.0
= 12.0= 13.5[DYNAMICS]integrator
nstepstimestepeneout_period
crdout_period
rstout_period= LEAP= 10000= 0.002= 250= 250= 10000[CONSTRAINTS]
rigid_bond= YES# use SHAKE/SETTLE[ENSEMBLE]
ensemble
tpcontrol
temperaturepressure= NPT= LANGEVIN= 300.00= 1.0[RPATH]
nreplica
rpath_period
delta
rest_function= 8= 500= 0.035= 1 2[SELECTION]
group1group2group3group4group5= atomindex:15= atomindex:17= atomindex:19= atomindex:25= atomindex:27[RESTRAINTS]
nfunctions= 2Stronger force constant
Images are automaticallyread from rst files

function1= DIHEDconstant1= 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0
reference1= 0 0 0 0 0 0 0 0
select_index1 = 1 2 3 4 # PHIfunction2= DIHED
constant2= 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0
reference2= 0 0 0 0 0 0 0 0
select_index2 = 2 3 4 5 # PSI

12

Running a job$ cd 4_string_method/$ lsrun.inp
run.sh$ cat run.sh#!/bin/bash
export OMP_NUM_THREADS=3
source /home2/data/env/setpath.CPU.DP
mpirun -np 8 spdyn run.inp >run.out$ ./run.shThis will finish in about 13 minutes
13

Step5. Visualization of pathways using aPython script$ cd 5_visualize/$ lsplot.py
pmf.matxi.matyi.mat# plot.py visualize ../4_string_method/*.rpath$ python3 plot.py# display plot.png or please scp it to your local machine$ display plot.png

14

plot.png
Converged
pathway

Initialpathway15

Step6. Umbrella sampling around theconverged pathway
Umbrellasampling…
16

Control file for GENESIS spdyn
The RPATH module in GENESIS can simulate replicas
without string method performed (rpath_period = 0).
[INPUT]
rstfile = ../4_string_method/{}.rst[OUTPUT]
logfiledcdfilerstfilerpathfile===={}.log{}.dcd{}.rst{}.rpath[ENSEMBLE]
ensemble
tpcontrol
temperaturepressure= NPT= LANGEVIN= 300.00= 1.0[RPATH]
nreplica
rpath_period
rest_function

= 8
= 0
= 1 2

[ENERGY]
forcefield
electrostatic
switchdist
cutoffdist
pairlistdist= CHARMM= PME= 10.0= 12.0= 13.5[DYNAMICS]
integrator
nstepstimestepeneout_periodcrdout_period
rstout_period[SELECTION]
group1group2group3group4group5= atomindex:15= atomindex:17= atomindex:19= atomindex:25= atomindex:27= LEAP= 10000= 0.002= 50= 50= 10000[RESTRAINTS]
nfunctions= 2[CONSTRAINTS]
rigid_bond= YESStronger force constant
Images are automaticallyread from rst files

function1= DIHEDconstant1= 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0
reference1= 0 0 0 0 0 0 0 0
select_index1 = 1 2 3 4 # PHI# use SHAKE/SETTLE

function2= DIHEDconstant2= 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0
reference2= 0 0 0 0 0 0 0 0
select_index2 = 2 3 4 5 # PSI

17

Running a job$ cd 6_umbrella/$ lsrun.inp
run.sh$ cat run.sh#!/bin/bash
export OMP_NUM_THREADS=3
source /home2/data/env/setpath.CPU.DP
mpirun -np 8 spdyn run.inp >run.out$ ./run.shThis will finish in about 13 minutes
18

Step7. Analysis of free energy profile
along the converged pathway$ cd 7_analysis/$ lsmbar_analysis.inp
remd_convert.inppathcv_analysis.inp
plot.pypmf_analysis.inp
run.shtrj_analysis.inp# run the analysis$ ./run.sh# display plot.png or please scp it to your local machine$ display plot.png

19What run.sh does
• Calls trj_analysis to calculate dihedral angles
from umbrella sampling data
• Calls mbar_analysis to perform MBAR andobtain unbiasing weights from dihedral angle
data• Calls pathcv_analysis to get tangential andperpendicular coordinates to the pathway• Calls pmf_analysis to evaluate free energy
profile (pmf.dat) from the MBAR weights andtangential coordinates to the pathway
20

plot.png
Barrier suggests a transition state21

Advanced topics• REUS instead of Umbrella sampling• Cartesian coordinates instead of dihedral
angles• Targeted MD for generating an initial
pathway• Tool for preparing rst files along an initial
pathway (rpath_generator)

Please check GENESIS website (some topics will be soo
illustrated in updated tutorials)
22

References
http://www.aics.riken.jp/labs/cbrt/tutorial/

23