最小自由エネルギー経路計算法
寺田透

平成 26 年 7 月 22 日

1

自由エネルギー地形

自由度 N の系のカノニカルアンサンブルを考える。系の座標を x = (x1 , . . . , xN ) とす
ると、確率密度関数は以下で与えられる。

ρ (x) = Z

−1

[

V (x)
exp −
kB T

]

(1)

ここで、V (x) はポテンシャルエネルギー関数、kB はボルツマン因子、T は温度、Z は
∫

Z=

[

]

V (x)
exp −
dx
kB T

(2)

である。N 次元の x から、n 次元の z = (z1 , . . . , zn ) への変換 z = θ (x) を考える。

θ (x) = (θ1 (x) , . . . , θn (x))(3)

z 空間における、確率密度分布 P (z) は、∫P (z) =

ρ (x) δ (z − θ (x)) dx(4)となる。従って、自由エネルギー地形 F (z) は以下で与えられる。

F (z) = −kB T ln P (z)= −kB T ln Z

2

−1

∫

[

]

V (x)
exp −
δ (z − θ (x)) dx
kB T

(5)

最小自由エネルギー経路

状態 A から状態 B への遷移（立体構造変化）を考える。２つの状態を繋ぐ経路のうち、実際の遷移で通過する確率が最も高い経路が最小自由エネルギー経路である。
最小自由エネルギー経路と類似した概念に、最小エネルギー経路がある。これは、遷移状態から出発して、ポテンシャルエネルギーの勾配の方向に進み、始状態と終状態に至る経路のことである。経路上の点に静止した状態から出発した微小時間 ∆t の間の運動を考える。運動方程式は、


 ẋk

= vk


 v̇k

= −

(6)

1 ∂V (x)mk ∂xk1と書ける。ここで vk は速度、mk は質量である。これを解くと、移動度 ∆x は、∆xk = −

1 ∂V (x) ∆t2
mk ∂xk2(7)で与えられる。従って、経路を x̂ (α) = (x̂1 (α) , . . . , x̂N (α)) と表すと、この経路上では、
経路とポテンシャルエネルギーの勾配を質量で割ったものが平行になる。

∂ x̂k (α)
1 ∂V (x̂ (α))
∥
∂α
mk
∂xk

(8)

ここで、α (0 ≤ α ≤ 1) は経路上の位置を表すパラメータであり、x̂ (0) は始状態、x̂ (1) は終状態を表す。このような運動を仮定して得られた経路は、固有反応座標（Intrinsic Reaction
Coordinate、IRC）と呼ばれる [1]。
これを z 空間における経路 ẑ (α) に変換する。
N∂ ẑi (α)∂θi (x̂ (α)) ∑∂θi (x̂ (α)) ∂ x̂k (α)
=
=
∂α
∂α
∂xk
∂α
k=1

(9)

であるから、式 (8) を代入すると、N∂ ẑi (α) ∑∂θi (x̂ (α)) 1 ∂V (x̂ (α))
∥
∂α
∂xk
mk
∂xk
k=1

(10)

となる。z 空間におけるポテンシャルエネルギー関数を U (z) とする（ただし z 空間で x
空間がすべて記述できるものと仮定する）と、n∂U (ẑ (α)) ∑∂U (ẑ (α)) ∂θj (x̂ (α))∂V (x̂ (α))
=
=
∂xk
∂xk
∂zj
∂xk
j=1

(11)

であるから、式 (10) は、
n ∑N∂ ẑi (α) ∑1 ∂θi (x̂ (α)) ∂θj (x̂ (α)) ∂U (ẑ (α))
∥
∂α
mk
∂xk
∂xk
∂zj
j=1 k=1

(12)

と書ける。
これを一般的な粗視化空間 z に拡張するには、自由エネルギー地形 F (z) を用いて、n∂ ẑi (α) ∑∂F (ẑ (α))∥
Mij (ẑ (α))
∂α
∂zj
j=1

(13)

のように書き換える [2]。この妥当性については??節で議論する。行列 M の要素は、−1

Mij (z) = P (z)∫ ∑N

1 ∂θi (x) ∂θj (x)ρ (x) δ (z − θ (x)) dx
mk ∂xk
∂xk
k=1

(14)

で与えられる [2]。従って、最小自由エネルギー経路を得るためには、状態 A と状態 B を繋ぎ、経路上の至る所で式 (13) を満たす経路を求めれば良い。

2

3

String 法

式 (13) を満たす経路 ẑ (α) を求めるために、経路を最適化する方法を考える。このために、経路が時間に依存して変化すると考え、変化は M ∂F/∂z のうち、経路に垂直な成分に比例するものとする。

γdẑ (α, t)∂F (ẑ (α, t))∂ ẑ (α, t)= −M (ẑ (α, t))+ λ (α, t)
dt
∂z
∂α

(15)

ここで γ は摩擦係数であり、最適化の速度を決める。また、λ (α, t) は、経路に平行な成分dẑ (α, t)を差し引くための関数である。経路が収束した段階では、= 0 であるから、dt

M (ẑ (α))∂F (ẑ (α))∂ ẑ (α)= λ (α)
∂z
∂α

(16)

となり、式 (13) が満たされていることがわかる。
M (ẑ (α, t)) を求めるためには、式 (14) より、ẑ (α, t) = θ (x) を満たす x について、アンサンブル平均を計算すればよい。この条件を満たす x をサンプルするために、束縛付き
定温分子動力学シミュレーションを行う。ポテンシャルエネルギーは、以下で与えられる。
nκ∑Uκ,z (x) = V (x) +(θj (x) − zj )2
2 j=1

(17)

ここで、κ は束縛の強さを表す定数である。このポテンシャルエネルギーの下で得られる
確率密度関数は、[−1
ρκ,z (x) = Zκ,z exp −Uκ,z (x)
kB T

]

(18)

となる。ここで、∫Zκ,z =

]

[

exp −

Uκ,z (x)
dx
kB T

(19)

である。デルタ関数は分散が 0 の正規分布で近似できることを利用すると、
[

]





n
V (x)
κ ∑
−1
−
lim ρκ,z (x) = lim Z −1 exp −(θj (x) − zj )2 
ZZκ,
exp
z
κ→∞
κ→∞
kB T
2kB T j=1

= P (z)−1 ρ (x) δ (z − θ (x))(20)

と書ける。従って十分大きい κ を用いてシミュレーションを行えば、Mij (z) ≃≃∫ ∑N

1 ∂θi (x) ∂θj (x)ρκ,z (x) dx
mk ∂xk
∂xk
k=1

Ns ∑
N
1 ∑
1 ∂θi (x (τ )) ∂θj (x (τ ))Ns τ =1 k=1 mk
∂xk
∂xk

(21)

のように近似的に求めることができる。ここで、τ はシミュレーションにおける時刻、Ns
はシミュレーションから得られたサンプル数である。

3

次に ∂F/∂z の計算法を考える。束縛付きのシミュレーションの結果得られる自由エネ
ルギーを Fκ (z) とすると、

Fκ (z) = −kB T ln Zκ,z(22)であるから、∂Fκ (z)
∂zj

[

∫

]

∂
Uκ,z (x)
=
exp −
dx
∂zj
kB T
[
]
∫
∂Uκ,z (x) −1
Uκ,z (x)=Zκ,z exp −
dx
∂zj
kB T
−1
−kB T Zκ,z∫κ (z j − θj (x)) ρκ,z (x) dx

=

(23)

を得る。一方、式 (22) は、κ が大きい時は、∫[]Uκ,z (x)
dx
kB T
[
]
∫
V (x)
= −kB T ln exp −kB TFκ (z) = −kB T ln

exp −





n
κ ∑
× exp −(θj (x) − zj )2  dx
2kB T j=1{

≃ −kB T ln Z −1(×Z

2πkB T
κ

∫

[

exp −

)n
2

]

}

δ (z − θ (x)) dx(= F (z) − kB T ln Z

V (x)
kB T

2πkB T
κ

)n
2

(24)

と書くことができる。従って、∂F (z)∂zj≃∂Fκ (z)
∂zj

∫

=
≃

κ (z j − θj (x)) ρκ,z (x) dx

Ns1 ∑κ [zj − θj (x (τ ))]Ns τ =1(25)

により計算することができる。
実際の経路の最適化計算は次のように行う。経路上の点を運動を計算するために、まず、
pとして、ẑ p = ẑ (αp ) とす
経路の離散化を行う。R 個の点で離散化する場合、αp =R−1る。次に、∆t 秒後の各点の位置を、Euler 法などを用いて式 (15) を解くことによって求めた後、各点の間隔が等しくなるように再配置する。これを、経路が収束するまで繰り返
すことによって、最適な経路を求める。この方法は、string 法 [2] と呼ばれる。

4

On-the-fly string 法

4String 法では、経路最適化の各ステップで、M と ∂F/∂z を求めるために、Ns 個のサンプルを得るシミュレーションが必要となる。この手続きは煩雑であり、時間もかかるため、より簡便な “on-the-fly” string 法が提案されている [3]。

4.1

方法

ここでは、式 (15) は、以下のように書き直される。

γdẑ (α, t)f (x (t)) κ [ẑ (α, t) − θ (x (t))] + λ (α, t) ∂ ẑ (α, t)
= −M
dt
∂α

(26)

f の成分は、式 (21) で Ns = 1 としたものである。
ここで、Mfij (z) =MN∑
1 ∂θi (x) ∂θj (x)k=1

mk ∂xk(27)∂xk

また、式 (15) における、∂F/∂z の置き換えには、式 (25) で Ns = 1 としたものを用いている。ここで、M 、∂F/∂z は、本来別々にアンサンブル平均をとったものであることに注意する必要がある。⟨X⟩ を X のアンサンブル平均と書くことにすると、一般に ⟨A⟩ ⟨B⟩と、⟨AB⟩ は異なる。
この問題を解決するために、x とは独立に運動する変数 y を導入し、式 (26) を以下のように書き換える。

γdẑ (α, t)f (x (t)) κ [ẑ (α, t) − θ (y (t))] + λ (α, t) ∂ ẑ (α, t)
= −M
dt
∂α

(28)

このようにすると、右辺第 1 項のアンサンブル平均は、∫

f (x) κ [z − θ (y)] ρκ,z (x) ρκ,z (y) dxdyM{∫=f (x) ρκ,z (x) dx
M= M (z)} {∫}κ [z − θ (y)] ρκ,z (y) dy∂F (z)∂z(29)となり、別々にアンサンブル平均をとったものの積と等しくなる。
この方法に基づく経路の最適化は、次のように行う。まず、前節と同様に、経路を R 個の点 ẑ p (p = 1, . . . , R) で離散化する。更に、束縛付きポテンシャルエネルギー関数

Uκ,ẑ p (x) = V (x) +nκ∑[θj (x) − ẑpj ]2
2 j=1

(30)

で平衡化した系を 2 コピーずつ、計 2R コピー用意する。次いで、コピーそれぞれについて、束縛付き分子動力学シミュレーションを行う。ここではまず、時刻 t における x (t)、y (t)、ẑ p (t) から x (t + ∆t)、y (t + ∆t) を求める。並行して、以下により経路上の点の位置を更新する。
f (x (t)) κ [ẑ p (t) − θ (y (t))] ∆tẑ ∗p (t + ∆t) = ẑ p (t) − γ −1 M

5

(31)

この段階では、まだ式 (26) の第 2 項を考慮していないので、ẑ p (t + ∆t) に ∗ を付した。式(26) の第 2 項は、経路に沿った移動をキャンセルするためのものである。従って、経路上に等間隔に点を配置している場合は、ẑ ∗p (t + ∆t) が等間隔になるように再配置すれば良い。これらの操作を経路が収束するまで繰り返すことにより、最適な経路を求めることができる。

Alanine dipeptide の α → β 構造転移の最小自由エネルギー経路

4.2

Alanine dipeptide (Ace-Ala-Nme) は水溶液中で、α ヘリックス様の構造と、polyproline
II ヘリックス様の構造をとる。ここでは、両安定構造の間の構造転移の最小自由エネル
ギー経路を求める。反応座標として Ala の二面角 ϕ、ψ を用いることにすると、式 (27) よf の成分は、
り、M


(
)
N
∑
1 ∂ϕ (x) 2



mk
∂xk
f =
M
 Nk=1
 ∑ 1 ∂ψ (x) ∂ϕ (x)
k=1

mk ∂xk∂xk

N∑
1 ∂ϕ (x) ∂ψ (x)mk ∂xk
k=1

∂xk

(
)
N
∑
1 ∂ψ (x) 2

k=1

mk









(32)

∂xk

と書ける。本シミュレーションを行う、µ2 lib のアプリケーションプログラムの概要は以
下の通り。

4.2.1

プログラムの概要

まず初めに、二面角を束縛する DiheRest クラスを作成する（説明に使用しない変数、関数は省略している）。

1 class DiheRest : public InteractionBase {2 protected:
3
int num;4
int *list[4];5double *z;
6double *THETA;
7
double *kappa;
8
double **dTHETA;
9 public:
10
int get_ene_num(void);11
void setup(Param& p, Conf& c, Parallel& para,Molecule &m);
12
void force(double *x,double **th_f,double *ene,MoleculeBase &m,
13
DynamicsContext &context);
14 };µ2 lib では、相互作用を計算するクラスは、InteractionBase クラスの派生クラスとして作成する。束縛する二面角の数を変数 num に格納する（ここでは num = 2）。list は 4×num
6

の 2 次元配列で、list[j][i] に i 番目の二面角の j 番目の原子のインデックスを格納する。z、THETA、kappa は大きさ num の配列で、それぞれ、そのコピーにおける ẑ p の座標、二面角の値、力の定数を格納する。dTHETA は num ×12 の 2 次元配列で、dTHETA[i][3*j]に i 番目の二面角の j 番目の原子の x 座標に関する微分を、dTHETA[i][3*j+1] に i 番目の二面角の j 番目の原子の y 座標に関する微分を、dTHETA[i][3*j+2] に i 番目の二面角の j 番目の原子の z 座標に関する微分を格納する（list とはインデックスの順番が逆になっていることに注意）。force() 関数は以下の通り。

1 void DiheRest::force(double *x,double **th_f,double *ene,2
MoleculeBase &m,DynamicsContext &context)3 {4
int i,i3;5
double d,fac;
6
double *cosTHETA,*sinTHETA;
7
double *f;89ene[0]=0.0;
10 // Only master process calculates force and energy.
11
if(master) {12
f=th_f[0];1314ALLOCATE_DOUBLE_ARRAY(cosTHETA,num);
15
ALLOCATE_DOUBLE_ARRAY(sinTHETA,num);16
Dihe::calc_torsion(x,list[0],list[1],list[2],list[3],170,num,dTHETA,cosTHETA,sinTHETA);
18
vdAcos(num,cosTHETA,THETA);
19
for(i=0;i<num;i++) {20
if(sinTHETA[i] < 0.0) {21
THETA[i]*=-1.0;
22
}
23
}
24
for(i=0;i<num;i++) {25
d=THETA[i]-z[i];26d-=2.0*PI*floor(d/(2.0*PI)+0.5);27
ene[0]+=0.5*kappa[i]*SQR(d);28fac=-kappa[i]*d;29
i3=list[0][i]*3;30
f[i3++]+=fac*dTHETA[i][0];31f[i3++]+=fac*dTHETA[i][1];32f[i3 ]+=fac*dTHETA[i][2];33/***** 中略 *****/
34
i3=list[3][i]*3;
7

35
f[i3++]+=fac*dTHETA[i][9];36f[i3++]+=fac*dTHETA[i][10];37f[i3 ]+=fac*dTHETA[i][11];38}39
FREE_ARRAY(cosTHETA);40
FREE_ARRAY(sinTHETA);
41
}
42 }
16 行目の Dihe::calc torsion() は Dihe クラスのユーティリティ関数で、座標 x と原子リスト（list[0]∼list[3]）、担当領域（0、num）から、二面角の座標に関する微分
dTHETA と二面角のコサイン、サインの値（cosTHETA、sinTHETA）を計算する。
次に、ẑ p の時間発展を担う Tstring クラスを作成する。
1 class Tstring {2 protected:
3
int num;
4
double *z;
5double **M;
6
double gamma,dt;
7 public:8
void setup(DiheRest& r,double g,Conf& c);
9void update(DiheRest& r,Param& p,Parallel& para);
10
void rearrange(DiheRest& r,Parallel& para,bool write_crd);11
void write_rst(DiheRest& r,char *fname);
12 };ここで、gamma は ẑ p の質量、dt は時間刻みである。setup()、write rst() は、それぞれ、Tstring クラスオブジェクトの初期化と、シミュレーションを再開する際に使用する、Amber フォーマットの束縛設定ファイルの書き出しを行う関数である。式 (27)、(31) の計算を行う update() 関数は以下の通りである。

1 void Tstring::update(DiheRest& r,Param& p,Parallel& para)
2 {3
int i,j,k,l,m,n;4
int **list;
5
double msum,fac,gamma_inv,dti;
6
double *THETAy,*kappa,**dTHETA,*mass,mass_inv,*z_new;7
int size,rank;
8
MPI_Comm comm;
9
MPI_Status status;
10
11
if(!para.get_local_master()) {12
return;
8

13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

}
size
rank
comm

= para.get_global_size();= para.get_global_rank();= para.get_global_comm();mass= p.get_molecule_param().get_mass();list= r.get_list();kappa= r.get_kappa();dTHETA= r.get_dTHETA();gamma_inv = 1.0/gamma;
if(rank % 2 == 0) {/* x */
ALLOCATE_DOUBLE_ARRAY(THETAy,num);
MPI_Recv(THETAy,num,MPI_DOUBLE,rank+1,rank+1,comm,&status);
ALLOCATE_DOUBLE_ARRAY(z_new,num);} else {/* y */
THETAy= r.get_THETA();MPI_Send(THETAy,num,MPI_DOUBLE,rank-1,rank,comm);}
if(rank % 2 == 0) {/* x */
for(i=0;i<num;i++) {fac=0.0;
for(j=0;j<num;j++) {msum=0.0;
for(k=0;k<4;k++) {m=list[k][i];n=-1;for(l=0;l<4;l++) {if(list[l][j] == m) {n=l;}}if(n >= 0) {mass_inv=1.0/mass[m];msum+=mass_inv*dTHETA[i][k*3 ]*dTHETA[j][n*3 ];msum+=mass_inv*dTHETA[i][k*3+1]*dTHETA[j][n*3+1];
msum+=mass_inv*dTHETA[i][k*3+2]*dTHETA[j][n*3+2];
}
}

9

54
M[i][j]=msum;
55
}
56
}
57
for(i=0;i<num;i++) {58
fac=0.0;59
for(j=0;j<num;j++) {60
fac+=M[i][j]*kappa[j]*(z[j]-THETAy[j]);61}62z_new[i]=z[i]-gamma_inv*fac*dt;
63
}
64
65
for(i=0;i<num;i++) {66
z[i]=z_new[i];67}68FREE_ARRAY(z_new);69
FREE_ARRAY(THETAy);
70
}
71 }
この関数は、第 3 引数で、プロセスをコピー（イメージの数 ×2）ごとに分割した Parallel
クラスオブジェクトを受け取る。この計算は、各コピーのローカルマスタープロセスのみが行う（11∼13 行目）。また、前節で述べた、独立に運動する 2 つのコピーのうち、偶数のグローバルランクを持つ x のコピーを担うプロセスが ẑ p の時間発展の計算を行うため、奇数のグローバルランクを持つ y のコピーを担うプロセスは二面角のデータを、x のコピーを担うプロセスに送る（25∼32 行目）。また、x のコピーを担うプロセスは、35∼56 行目
f の計算を行う。ここで、∂ϕ (x) /∂xk は k が Ace のカルボニル炭素、Ala ので、行列 M
アミド窒素、Cα、カルボニル炭素の各原子を指している場合のみ値を持ち、∂ψ (x) /∂xk
は k が Ala のアミド窒素、Cα、カルボニル炭素、Nme のアミド窒素の各原子を指している場合のみ値を持つことに注意すると、積が値を持つのは、2 つの二面角が共有している
Ala のアミド窒素、Cα、カルボニル炭素のみであることがわかる。従って、i 番目の二面角を構成する原子のリストを順に見ていき（39、40 行目）、その原子が j 番目の二面角を構成する原子のリストにも存在する（43 行目）時のみ、二面角の微分の積を計算すればよ
f と、y から ẑ ∗ (t + ∆t)
いことになる（47∼52 行目）。最後に、このようにして計算した M
pを計算する（57∼63 行目）。ẑ ∗p (t + ∆t) を等間隔になるように再配置する rearange() 関数は以下の通り。

1 void Tstring::rearrange(DiheRest& r,Parallel& para,bool write_crd)
2 {3
int i,j,q,total_steps,save_crd_freq;
4
double d,d2,d_inv,len;5
double *L,*z_all,*z_new,**dz;
6
int size,rank,master;
7
MPI_Comm comm;
10

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

double z_offset;
double *THETA;
if(para.get_local_master()) {
size
rank
mastercomm= para.get_global_size();= para.get_global_rank();= para.get_global_master();= para.get_global_comm();

if(master) {ALLOCATE_DOUBLE_ARRAY(z_all,(size*num));
ALLOCATE_DOUBLE_ARRAY(z_new,(size*num));
ALLOCATE_DOUBLE_MATRIX(dz,size,num);
ALLOCATE_DOUBLE_ARRAY(L,size);}
MPI_Gather(z,num,MPI_DOUBLE,z_all,num,MPI_DOUBLE,0,comm);
if(master) {L[0]=0.0;
for(i=1;i<size;i++) {d2=0.0;for(j=0;j<num;j++) {dz[i][j]=z_all[i*num+j]-z_all[(i-1)*num+j];z_offset=2.0*PI*floor(dz[i][j]/(2.0*PI)+0.5);dz[i][j]-=z_offset;
z_all[i*num+j]-=z_offset;
d2+=SQR(dz[i][j]);}d_inv=1.0/sqrt(d2);for(j=0;j<num;j++) {dz[i][j]*=d_inv;}
L[i]=L[i-1]+d2*d_inv;}
d=L[size-1]/(double)(size-1);
q=1;
for(i=1;i<size-1;i++) {len=d*(double)i;
while(L[q] < len && q < size-1) q++;for(j=0;j<num;j++) {z_new[i*num+j]=z_all[(q-1)*num+j]+(len-L[q-1])*dz[q][j];

11

49
}
50
}
51
for(j=0;j<num;j++) {52
z_new[j]=z_all[j];53z_new[(size-1)*num+j]=z_all[(size-1)*num+j];54}
55
for(i=0;i<size*num;i++) {56
z_new[i]-=2.0*PI*floor((z_new[i]+PI)/(2.0*PI));
57
}
58
}
59
MPI_Scatter(z_new,num,MPI_DOUBLE,z,num,MPI_DOUBLE,0,comm);
60
if(master) {61
FREE_ARRAY(z_all);62
FREE_ARRAY(z_new);63
FREE_MATRIX(dz,size);64
FREE_ARRAY(L);
65
}
66
}
67
68
MPI_Bcast(z,num,MPI_DOUBLE,0,para.get_local_comm());
69
r.set_z(z);70/***** z の書き出し（省略） *****/
71 }

q–2
dz[q]
0

q

len=d*iq–1
図 1: イメージの再配置

この関数は、第 2 引数で、プロセスをイメージごと（各イメージは x、y の 2 つのコピーを含む）に分割した Parallel クラスのオブジェクトを受け取る（update() 関数とは異なるので注意）。ẑ ∗p (t + ∆t) の再配置は、グローバルマスターのみが行う（18 行目、25 行12目）。各イメージのマスターから z の内容をグローバルマスターの z all に集めた後（24行目）、グローバルマスターはイメージ間の距離から string の全長を求め（40 行目）、新
たなイメージの間隔を求める（42 行目）。また、同時に隣接したイメージ間を結ぶ単位ベ
クトル dz を求める（38 行目）。これらを用いて、図 1 のように、イメージの新たな座標を求める（43∼54 行目）。最後に、これを各イメージのマスターに分配し（59 行目）、各
イメージのマスターは、そのイメージに属するローカルプロセスにブロードキャストする（68 行目）。各ローカルプロセスは、次のステップの計算に備えて新しい z を DiheRest ク
ラスのオブジェクトにセットする（69 行目）。main() 関数は以下の通り。

1 int main(int ac,char **av)
2 {3
AmberParam p;
4
AmberConf c;
5
Molecule m;
6
Dynamics a;
7
Parallel para,para2;
8bool restart_flag;
9
int i,ncopy,nsteps,thermostat,rearrange_freq;
10
double gamma;
11char *rst_out_fname;
12
DiheRest r;
13
Tstring z;
1415 // check args
16
ncopy=Conf::get_int_option(ac,av,"-ncopy");17
gamma=Conf::get_double_option(ac,av,"-gamma");18
rst_out_fname=Conf::get_string_option(ac,av,"-rst_out");
19
rearrange_freq=Conf::get_int_option(ac,av,"-rearrange_freq");20
21 // start parallel
22
para.init(&ac,&av);
23
para2=para;
24
para.multicopy_mode(ncopy);
25
para2.multicopy_mode(ncopy/2);
26
27 /***** ファイルの読み込み（省略） *****/
28
29 // start simulation
30
srand48((long)c.get_random_seed());
31
c.set_restint(false);
32
m.add_interaction(&r);
33a.setup(p,c,para,&m);13

34
thermostat=c.get_thermostat();35
nsteps=c.get_nsteps();36
restart_flag=c.get_restart();37
z.setup(r,gamma,c);38
for(i=0;i<nsteps;i++) {39
z.update(r,p,para);40
if(thermostat == Conf::THERMOSTAT_TYPE_NONE) {41
a.leapfrog(1,restart_flag);42} else if(thermostat == Conf::THERMOSTAT_TYPE_BERENDSEN) {43
a.berendsen(1,restart_flag);44} else if(thermostat == Conf::THERMOSTAT_TYPE_LANGEVIN) {45
a.langevin(1,restart_flag);46}
47
restart_flag=true;48
if((i+1)%rearrange_freq == 0) {49
z.rearrange(r,para2,((i+1)%c.get_save_crd_freq()==0));
50
}
51
}
52
if(para.get_local_master()) {53
z.write_rst(r,rst_out_fname);54}
55
para.finalize();56 }ここでは、MoleculeBase クラスに内蔵されている RestInt クラスによって束縛の計算が行われないように設定し（31 行目）、代わりに、DiheRest クラスのオブジェクトをMoleculeBase クラスのオブジェクトに追加している（32 行目）。Tstring クラスのオブジェクトを初期化し（37 行目）、nsteps のシミュレーションを開始する（38 行目）。こ
こでは、まず時刻 t における x もしくは y の座標からその時点における自由エネルギーの勾配を求め、ここから ẑ ∗p (t + ∆t) を求める（39 行目）。次いで、x もしくは y の座標を ∆t 進める（40∼46 行目）。ここでは、DiheRest クラスのオブジェクトが保持しているイメージの座標 z は更新されていないので、時刻 t における束縛力が計算される。その後、ẑ ∗p (t + ∆t) の再配置を行い、ẑ p (t + ∆t) を得る。更に、新しいイメージの座標を
DiheRest クラスのオブジェクトに登録する（49 行目）。再配置は、x のコピーを担う全てのプロセスの間の通信を伴うため、毎ステップ再配置するのではなく rearrange freq ス
テップごとに再配置するようにすることもできる。この間は、ẑ p は更新されないが、ẑ ∗pは更新されていくので、結局、この期間の自由エネルギーの勾配の平均を用いて、ẑ p を更新するのと同じことになる。シミュレーション終了後は、シミュレーションの再開に備え、現在の ẑ p を Amber の束縛設定ファイル形式で書き出す（53 行目）。

14

4.2.2シミュレーションの結果

まず、ϕ-ψ 空間上の初期パスに 8 つの点（イメージ）を配置した（R = 8）。各イメージの (ϕ, ψ) 座標は、それぞれ、(−40.0, 130.0)、(−40.0, 105.0)、(−40.0, 80.0)、(−40.0, 55.0)、(−40.0, 30.0)、(−40.0, 5.0)、(−40.0, −20.0)、(−40.0, −45.0) とした。それぞれのイメージ
は独立に運動する 2 つのコピーを持つため、全体で 16 コピーのシミュレーションを行った。Alanine dipeptide の主鎖の二面角 (ϕ, ψ) を各イメージの座標付近に束縛した 300 K、
1 bar の定温定圧分子動力学シミュレーションを行い、string 法のための初期構造とした。
これらを用いて、300 K、1 bar の定温定圧条件下で string 法の計算を 1 ns 行った。ここでは、パラメータ κ は 1000 kcal mol−1 rad−2 とし、γ は 500 とした。再配置の頻度は 10
ステップごと（20 fs ごと）とした。Intel Xeon E5-2670 CPU を 2 基備える計算ノード 2台（合計 32 コア）を用い、コピー当たり 1 プロセス 2 スレッドの条件でシミュレーショ
ンを行ったところ、1 ns のシミュレーションに約 1 時間 30 分かかった。図 2 は、初期パス（赤）、25 ps 後（緑）、50 ps 後（青）、100 ps 後（赤紫）、500 ps 後（水色）、1 ns 後（黄色）のパスを、別途求めた自由エネルギー地形の上に重ね合わせてプロットした図である。ここから、収束したパスは、polyproline II ヘリックス様構造の自由エネルギー極小状態から、α ヘリックス様構造の自由エネルギー極小状態をつなぐ、自由エネルギー障
壁が最も低いパスに一致していることがわかる。図 3 の赤線は、
d (t) = 

R
1 ∑

R p=1

1
2

2

|ẑ p (t) − ẑ p (0)|(33)

の時間変化を表している。ここから、0.4 ns 程度でパスが収束していることがわかる。緑線は、再配置の頻度を 1（毎ステップ）として同様なシミュレーションを行った結果を表しているが、収束の速度はほとんど変わらないことから、計算の効率を考慮すると、10 ス
テップに 1 回程度再配置すれば十分であることがわかる。図 4 において、赤線は後半 0.5
ns の ẑ p の平均位置 z p における自由エネルギー勾配から計算した potential of mean forceを表している。これは以下の式に従って計算した。
Ns1∑∑[z i − θ (y (τ ))] · (z i+1 − z i−1 )2 i=1 τ =1
p

Fp =

(34)

ここで、

z0 = z1(35)

z R+1 = z R(36)である。このシミュレーションは 1 ns 行い、5 ステップ（10 fs）毎に保存した、θ (y (τ ))のトラジェクトリから、後半 0.9 ns 分を用いた（Ns = 9 × 104 ）。また、緑線は、別途行った定温定圧シミュレーションにおいて、z p ± 5◦ の範囲の二面角をとる確率を求め、この値から自由エネルギーを求めたものである。String 法においては、離散化による誤差のた
めにピークの位置がずれているものの、自由エネルギー障壁の高さについては良く一致していることがわかる。

15

180
120

ψ [degree]

60
0
-60
-120
-180
-180

-120

-60

0

60

120

180

ϕ [degree]

図 2: パスの時間変化

50
45
40
d (t) [degree]

35
30
25
20
15
10
5
0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

t [ns]

図 3: 初期パスからのずれの時間変化

16

0.8

0.9

1

2.5

Free energy [kcal mol−1 ]

2
1.5
1
0.5
0
-0.5
1

2

3

4

5

6

7

8

Image index

図 4: Potential of mean force の比較

参考文献[1] Kenichi Fukui, Shigeki Kato, and Hiroshi Fujimoto “Constituent analysis of thepotential gradient along a reaction coordinate. Method and an application to CH4+ T reaction” J. Am. Chem. Soc. 97, 1–7 (1975).
[2] Luca Maragliano, Alexander Fischer, Eric Vanden-Eijnden, and Giovanni Ciccotti“String method in collective variables: Minimum free energy paths and isocommittor
surfaces” J. Chem. Phys. 125, 024106 (2006).
[3] Luca Maragliano and Eric Vanden-Eijnden “On-the-fly string method for minimum
free energy paths calculation” Chem. Phys. Lett. 446, 182–190 (2007).

17