ゆらぎの定理を利用した自由エネルギー計算に関する理論的背景と非平衡分子動力学シミュレーションへの応用

束村 晴; 菅波 祐介; 赤津 裕哉; 髙橋 颯人; 西村 好史; 中井 浩巳

doi:10.2477/jccj.2024-0033

Abstract

Enhanced sampling methods have been applied to handle chemical reactions rarely observed in limited-scale molecular dynamics simulations. Some of these methods intentionally put the system in a nonequilibrium state. We focus on free energy calculation methods based on the fluctuation theorem to directly estimate free energy from simulations of nonequilibrium systems. In this paper, we derive Jarzynski's equality and Bennett acceptance ratio method, and present their illustrative applications to free energy calculation.

Translated Abstract

1 はじめに

分子動力学 (MD) シミュレーションでは熱力学的なアンサンブルによって複雑系の現象を取り扱う．MDシミュレーションにおいて第一原理計算や反応力場によってあらわに化学反応を取り扱う場合，長時間のシミュレーションが必要となる．そのため，メタダイナミクス法 [1]，レプリカ交換MD法 [2]などの拡張アンサンブル法が用いられる．またこれらの手法を用いることで一定の条件のもとで自由エネルギーを計算することができる．例えばメタダイナミクス法の場合，特定の集団座標 (collective variable, CV) に対する自由エネルギー計算が可能である [3,4,5,6,7]．現在，我々は事前にCVを指定する必要がない反応解析手法としてナノリアクター分子動力学 (NMD) 法 [8]に注目している．NMD法では系をソフトポテンシャルにより圧縮し，系中の温度や圧力，粒子同士の衝突頻度を増加させることで，化学反応のサンプリングを加速させる．しかし見出された反応経路に沿った自由エネルギーを求める手法は確立されていない．NMD法において系は平衡状態にはないため，非平衡を取り扱う理論に基づいて自由エネルギーを得る手法が必要となる．本稿では平衡系の熱力学と統計力学から出発し，非平衡系を議論するにあたって必要となる設定や知識を要約する [9]．これらをもとに詳細ゆらぎの定理を経てJarzynski等式 [10]を導出し，非平衡系における自由エネルギー計算の理論を紹介する．最後にこれらの理論を利用したMDシミュレーションの具体例を挙げる．

2 非平衡系の熱力学・統計力学

まず従来の平衡系における熱力学・統計力学について簡単にまとめる．平衡系の熱力学は式(1)，(2)に表される熱力学第一法則，熱力学第二法則に基づいて記述される．

ΔU=δQ+δW

(1)

ΔS−δQT≥0

(2)

∆Uは内部エネルギー変化，δQは系への吸熱，δWは系へなされた仕事，∆Sはエントロピー変化，Tは系の絶対温度である．式(2)の左辺はエントロピー生成σと定義され，これを用いると熱力学第二法則は次式で表される．

σ≥0

(3)

熱力学においてHelmholtzの自由エネルギーFは次式で表される．

F=U−TS

(4)

Uは内部エネルギー，Sはエントロピーである．

平衡系の統計力学は，ミクロな状態の確率分布とマクロな熱力学量の関係を与える．平衡状態において，系が状態xを取る確率は次式に示すカノニカル分布に従う．

Pcan(x)=e−βExZ=eβ(F−Ex)

(5)

ただし，βは熱浴の逆温度，Zは分配関数，E_xは状態xにおけるエネルギーである．統計力学におけるFは次式で表される．

F=−1βlnZ

(6)

従来の熱力学・統計力学は高い普遍性を保つ一方で，非平衡系や系の時間発展を記述することができない．

この課題を解決するため，非平衡系における熱力学量では，時刻tにおける系の内部エネルギーU(t)は，状態xが現れる確率P(x,t)と状態xにおけるエネルギーE_x(t)により次式で定義される．

U(t)=∑xPx,tEx(t)

(7)

式(7)を時間微分することで次式が得られる．

dU(t)dt=∑x∂P(x,t)∂tEx(t)+∑xP(x,t)∂Ex(t)∂t

(8)

熱の時間微分をQ˙(t)，仕事の時間微分をW˙(t)とする．箱の中の理想気体の系から，熱と仕事に関するミクロな表式を考察する．理想気体に対する熱の出入りは温度の変化で表現される．温度は粒子のエネルギー状態の確率分布を特徴づける．よって，熱に関する操作はP(x,t)を変更する操作である．

気体に対する力学的な仕事は，ピストンの操作などの体積の変化に対応する．粒子のエネルギーは系の体積に依存する量なので，仕事に関する操作はエネルギーE_xを変える操作である．以上の考察から，式(8)より，次式の熱と仕事の時間微分の表式が得られる．

Q˙(t)=∑x∂P(x,t)∂tEx(t)

(9)

W˙(t)=∑xP(x,t)∂Ex(t)∂t

(10)

非平衡系に対して，エントロピーを定義する．式(5)より，分配関数Zは次式で表すことができる．

Z=e−βExPcan(x)

(11)

Σ_xP_can(x)=1に注意して，式(6)に式(11)を代入すると次式が成り立つ．

F=−1β∑xPcan(x)lne−βExPcan(x)=∑xPcan(x)Ex+1β∑xPcan(x)lnPcan(x)=U+1β∑xPcan(x)lnPcan(x)

(12)

式(4)に式(12)を代入することで，エントロピーSについて次式が得られる．

S=U−FT=−kB∑xPcan(x)lnPcan(x)

(13)

ここで，k_BはBoltzmann定数である．式(13)で定義される量はGibbsエントロピーと呼ばれ，温度に依存しない量である．Gibbsエントロピーを一般の時間変化する確率分布に拡張し，非平衡系におけるエントロピーとして定義する．

S(t)=−∑xP(x,t)lnP(x,t)

(14)

ただし，式(13)からk_Bを落とした．また，非平衡系におけるエントロピー生成は次式で定義される．

σ˙(t)=S˙(t)−βQ˙(t)

(15)

ただし，式(14) と同様にk_Bを落としている．

非平衡系は，時間変化するミクロな系の状態と，系の操作に対応するパラメータによって記述される．系のミクロな状態をx(t)，操作パラメータをλ(t)で表す．時刻t=0からt=τまでのx(t)をx_τと表し，これを経路と呼ぶ．同様に，t=0からt=τまでのλ(t)をλ_τで表し，これを操作プロトコルと呼ぶ．

操作プロトコルに対してどのような経路が実現するかは確率的に決定される．あるλ_τのもとでx_τが実現する確率密度は経路確率と呼ばれ，P[x_τ, λ_τ]で表される．x_τ，λ_τに依存するミクロな物理量をA^xτで表す．経路に依存するマクロな物理量Axτは，A^xτの経路積分によるアンサンブル平均で与えられる．

Axτ=∫DxτPxτ,λτA^xτ,λτ

(16)

x_τ，λ_τに対する確率的な熱と仕事をそれぞれQ^[xτ,λτ]，W^[xτ,λτ]と表す．Q^は系が経路を通る過程で熱浴から吸収する熱，W^は操作プロトコルによる仕事である．ある経路に沿った熱力学第一法則は次式で表される．

ΔE=Ex(τ)(τ)−Ex(0)(0)=Q^[xτ,λτ]+W^[xτ,λτ]

(17)

ここで，W^はλを介したエネルギー変化として次式で表すことができる．

W^[xτ,λτ]=∫0τdtdλ(t)dt∂Ex(t)(λ(t))∂λ

(18)

ここで，エネルギーE_x_(t)(t)は操作に依存して決定されることから，λの関数E_x_(t)(λ(t))として記した．

あるλ_τのもと初期状態がx(0)であるという条件で，x_τが実現する過程は順過程と呼ばれる．初期分布がP(x(0),0)であるときに順過程が実現する確率密度をP[x_τ|λ_τ]で表す．これに対し時間反転させたプロトコルをλ_τ^†のもと初期状態がx^†(0)であるという条件で，時間反転した経路x_τ^†が実現する過程は逆過程と呼ばれる．初期分布がP_R(x^†(0),0)=P(x(τ), τ)であるときに逆過程が実現する確率密度はP[x_τ^†|λ_τ^†]で表される．

Figure 1に順過程と逆過程を模式的に示した．ピストンの位置が操作パラメータλ(t)，ある時刻における系中のミクロな状態がx(t)に相当する．ピストンを押し込むときに系中で実現する経路が順過程，引き抜くときに実現する経路が逆過程に相当する．

Figure 1.

Schematic diagram of forward (upper panel) and reverse (lower panel) processes.

3 詳細ゆらぎの定理

平衡系での熱力学第二法則は，孤立系の状態変化はエントロピーが増大する方向に進行することを主張する．これは，エントロピーが増大する過程は起こりやすく，エントロピーが減少する過程は発生しにくいため，全体としてエントロピーが増大する方向に状態が変化することを意味する [11]．すなわち，ある過程の実現確率は熱力学量に依存していることを示唆している．

後述の詳細ゆらぎの定理は，順過程と逆過程の実現確率の比と熱力学量を結びつける関係式である．詳細ゆらぎの定理は非平衡を含む広い過程において成り立つ等式であり，非平衡系を記述する非常に強力な手段となる．特に，各過程の始状態がカノニカル分布を示す場合，仕事と自由エネルギーの関係で表されるゆらぎの定理はJarzynski等式 [10]に相当する．Jarzynski等式は非平衡過程の自由エネルギー計算に用いることができる．詳細ゆらぎの定理はHamiltonian系一般で成立する関係である [12]が，本節ではMarkovジャンプ過程に限定して定理の証明を行う．

Markov過程では系の状態が時刻tにおいて完全に決定され，状態の変化は直前の状態にのみ依存する．Markov過程の確率P(x,t)の時間発展は次式のマスター方程式によって記述される．

∂P(x′,t)∂t=∑x(≠x′)R(x′|x;t)P(x,t)−R(x|x′;t)P(x′,t)

(19)

ここでRは単位時間あたりの遷移確率であり，遷移速度と呼ばれる．時刻tの状態がxであるとき，時刻t+Δtの状態がx′(≠x)になる確率をR(x'|x;t)dtで与える．状態x'から別の状態へと遷移する確率は脱出率γ (x';t)として次式で定義される．

γ(x′;t):=∑x(≠x′)R(x|x′;t)

(20)

Markovジャンプ過程は状態が離散的で時間は連続なMarkov過程である．時刻0からτまでの時刻t_k (k=1,2,…,K), (t₀=0, t_K₊₁=τ)にK回の遷移が起きるMarkovジャンプ過程を考える．時刻t_k ≤ t ≤ t_k₊₁ (k=0,1,…,K) における状態をx(t)=x_kと表す．この過程における経路x_τは，初期状態x₀，遷移が起きた時刻t_k (k=1,2,…,K)，その遷移後の状態x_kにより一意に定められる．Figure 2 にMarkovジャンプ過程の経路を模式的に示す．

Figure 2.

Schematic diagram of Markov jump process. Horizontal axis represents time, vertical axis represents state. States are discrete in Markov jump process.

Markovジャンプ過程における経路に沿ったミクロな熱と仕事を求める．熱は状態xの遷移に伴うエネルギーの変化であるから，次式で与えられる．

Q^[xτ,λτ]=∑k=1K{Exk(tk)−Exk−1(tk)}

(21)

仕事は式(18)に示した一般的な定義により得られる．

W^[xτ,λτ]=∫0τ∂Ex(t)∂tx=x(t)dt

(22)

また，各時刻で系と熱浴との間に次式で表される詳細つり合いが成り立つと仮定する．

R(x′|x;t)exp−βEx(t)=R(x|x′;t)exp−βEx′(t)

(23)

以上の定義，仮定のもと，ある経路について順過程と逆過程の実現確率を考える．簡単のため，遷移速度Rが時刻tに依存しないとする．時刻t=0に初期状態x₀から出発し，t=t₁までx₀にとどまり続ける確率を計算する．微小時間Δtの間に状態が変化しない確率は脱出率γ(x₀)を用いて1−γ(x₀)Δtで表される．N時間ステップ後も状態が変化していない確率は次式で表される．

1−γ(x0)ΔtN

(24)

N=t₁/Δtを代入してΔt→0の極限をとる．

limΔt→01−γ(x0)Δtt1/Δt=exp−γ(x0)t1

(25)

x₀からx₁への遷移が起こる確率はR(x₁|x₀)dtで与えられる．その後，時刻t₂まで状態x₁にとどまり続ける確率はexp{−γ(x₁)(t₂−t₁)}である．これらを繰り返し，個々の遷移確率の積をとることで経路確率が得られる．

P[xτ]∝∏k=1KRxk|xk−1∏k=0Ke−γ(xk)(tk+1−tk)P(x0,0)

(26)

ただし，右辺からdt^Kの項を落とした．一般に遷移速度が時刻に依存する場合も同様に考える．遷移速度の時間依存性を特徴づけるパラメータをλ_τすると，遷移速度の時間変化を考慮した経路確率は次式で与えられる．

P[xτ|λτ]∝∏k=1KRxk|xk−1;t∏k=0Ke−∫tktk+1γ(xk;t)dtP(x0,0)

(27)

遷移確率の時間依存性を考慮した逆過程の経路確率は次式で与えられる．

P[xτ†|λτ†]∝∏k=1KRxk−1|xk;t∏k=0Ke−∫tktk+1γ(xk;t)dtP(xτ,τ)

(28)

順過程と逆過程の経路確率の比をとる．

P[xτ†|λτ†]P[xτ|λτ]=∏k=1KRxk−1|xk;tkRxk|xk−1;tkP(xτ,τ)P(x0,0)

(29)

式(23)の詳細つり合いの仮定より，式(29)の括弧内はある経路に沿った熱Q^[xτ,λτ]を用いて次式で表される．

∏k=1KRxk−1|xk;tkRxk|xk−1;tk=exp−∑k=1Kβ{Exk−1(tk)−Exk(tk)}=expβQ^[xτ,λτ]

(30)

次に，σ^[xτ,λτ]を次式で定義する．

σ^[xτ,λτ]=−lnP(xτ,τ)+lnP(x0,0)−βQ^[xτ,λτ]

(31)

これより，式(16)に従ってσ^[xτ,λτ]のアンサンブル平均を取ると，式(15)に示したエントロピー生成と一致することがわかる．σ^[xτ,λτ]は経路x_τにおけるエントロピー生成に相当する量であり，確率的なエントロピー生成と呼ばれる．

σ^[xτ,λτ]を用いて式(29)の右辺を表すことで詳細ゆらぎの定理に対する表式が得られる．

P[xτ†|λτ†]P[xτ|λτ]=e−σ^[xτ,λτ]

(32)

式(32)は，順過程と逆過程の起こる確率の比が，確率的なエントロピー生成によって定量的に記述できることを意味する．また，式(32)のアンサンブル平均を取ることにより，次式の積分型ゆらぎの定理が得られる [13]．

e−σ^[xτ,λτ]=∫DxτP[xτ†|λτ†]P[xτ|λτ]P[xτ|λτ]=1

(33)

次に，順過程および逆過程の始状態をカノニカル分布P_canと仮定し，仕事とHelmholtzの自由エネルギーを用いた積分型ゆらぎの定理の表式を求める．始状態の分布のエネルギーは時間反転対称性を満たすから，カノニカル分布は時間反転対称である．逆過程の始状態PR(x†(0),0)は順過程の終状態のカノニカル分布であるから，次式が成り立つ．

PR(x†(0),0)=Pcan(x(τ),τ)

(34)

この場合のエントロピー生成σ^w[xτ,λτ]について，式(5)，(17)，(34)を用いると次式が得られる．

σ^w[xτ,λτ]=−lnPR(x†(0),0)+lnP(x(0),0)−βQ^[xτ,λτ]=−lnPcan(x(τ),τ)Pcan(x(0),0)−βQ^[xτ,λτ]=−β(ΔF−ΔE)−β(ΔE−W^[xτ,λτ])=β(W^[xτ,λτ]−ΔF)

(35)

ただし，∆Fは順過程の始状態と逆過程の始状態の間のHelmholtzの自由エネルギー変化である．このような順過程および逆過程の取り方に対して，積分型ゆらぎの定理は次式のようになる．

∫DxτP[xτ|λτ]e−β(W^[xτ,λτ]−ΔF)=1

(36)

ΔFは状態量であるから右辺に移項することができる．両辺にexp(−βΔF)をかけることで，次式のJarzynski等式 [10]が得られる．

e−βW^[xτ,λτ]=e−βΔF

(37)

4 自由エネルギー計算

式(32)に示した詳細ゆらぎの定理は，非平衡過程における自由エネルギー変化を計算する手法に使用できる．ゆらぎの定理に基づく自由エネルギー計算手法について以下に述べる．また，以降の議論においてW^のハット記号 (^) は省略する．

式(37)に示したJarzynski等式を変形することで次式が得られる．

ΔF=−1βlne−βW

(38)

式(38)は，Helmholtzの自由エネルギー変化を系に行った仕事のアンサンブル平均として計算できることを意味する．式(38)の右辺をキュムラント展開すると，自由エネルギー変化を近似的に算出することができる [14]．

ΔF≈W−β2!W2−W2+β23!W3−3W2W+2W3−⋯

(39)

Bennett受容比 (Bennett acceptance ratio, BAR) 法 [15]は，順過程と逆過程，両方向についてなされた仕事を用いて自由エネルギーを計算する手法である．式(32)に式(35)を代入することで次式が得られる．

lnP[xτ|λτ]P[xτ†|λτ†]=βW−ΔF

(40)

ここで，P[x_τ|λ_τ]=P(W|F)は，順過程を用いたときの順過程での仕事Wが得られる条件付き確率，P[x_τ^†|λ_τ^†]=P(W|R)は逆過程を用いたときの順過程の仕事Wが得られる条件付き確率に相当する．これらの比は次式で表される．

P(W|F)P(W|R)=P(F|W)P(R|W)P(R)P(F)=P(F|W)1−P(F|W)nRnF

(41)

ここで，n_F, n_Rは順過程と逆過程の試行回数である．式(41)を式(40)に代入することで次式が得られる．

lnP(F|W)1−P(F|W)=βW−ΔF+1βlnnFnR=βW−ΔF+M

(42)

したがって，仕事Wが得られたときに，順過程から得られた確率と逆過程から得られた確率はそれぞれ式(43)，(44)で与えられる．

P(F|W)=11+e−β(W−ΔF+M)=exp−ln1+e−β(W−ΔF+M)

(43)

P(R|W)=11+eβ(W−ΔF+M)=exp−ln1+eβ(W−ΔF+M)

(44)

(n_F+n_R)回の試行から得られた仕事Wのうち，W_i (i=1,…,n_F) が順過程から，W_j (j=1,…,n_R) が逆過程から得られた確率の積を，ΔFをパラメータとする尤度関数Lとして次式で定義する．

L(ΔF)=∏i=1nFP(F|Wi)∏j=1nRP(R|Wj)

(45)

最尤推定法に基づき，式(45)のLを最大化することで尤もらしいΔF を推定する．ln Lを最大化する停留値の条件より，次式に示すBAR法の基本式が得られる．

∑i=1nF11+e−β(Wi−ΔF+M)=∑j=1nR11+eβ(Wj−ΔF+M)

(46)

BAR法は過程の始状態と終状態間のHelmholtzの自由エネルギー変化を推定する．式(46)では過程の中間地点までの自由エネルギーの推定はできない [16]．

詳細ゆらぎの定理に基づく自由エネルギー計算では，準静的でないプロトコルに対して系の緩和を待つことなく計算を行うことができる．一方で，多数の試行をしてアンサンブル平均をとる必要がある．サンプリングは平衡状態から出発しなくてはならないという制約もある [9]．

5 Steered MD法への応用

詳細ゆらぎの定理に基づく自由エネルギー計算の具体例としてsteered MD (SMD) 法 [17]による自由エネルギー計算を示す．

SMD法は系に時間依存のバイアスポテンシャルをかけるMDシミュレーション手法である．このバイアスポテンシャルをかけるという操作により算出される系への仕事を文献 [14]に従い導出する．バイアスポテンシャルVを反応座標ξにかける場合，Vは次式で与えられる．

V(r;λ)=k2ξ(r)−λ2

(47)

ただし，rは位置，kは力の定数，λはλ(t)=λ₀+vtで表される操作パラメータであり，vはポテンシャルの速度である．全HamiltonianH˜は，オリジナルのHamiltonian Hと運動量pにより次式で与えられる．

H˜(r,p;λ)=H(r,p)+V(r;λ)

(48)

時刻aからbまでにH˜系へされた仕事Wabは式(18)に基づいて次式で計算される．

Wab=∫abdtdλ(t)dt∂H˜(r,p;λ)∂λ=∫abdtdλ(t)dt∂∂λk2ξ(r)−λ2=−kv∫abdtξ(r)−λ0−vt

(49)

λをλ₀=λ(0)からλ_τ=λ(τ)まで変化させるものとし，その過程を順過程，反対方向の過程を逆過程とする．λ=λ(t)において，順過程で系にされた仕事をWF(λ)=W0t，逆過程で系にされた仕事をWR(λ)=Wτ−tτとそれぞれ定義する．Jarzynski等式に基づき，平均力ポテンシャル (potential of mean force, PMF) をλの関数として順過程と逆過程それぞれについて次式のように定める．

ΦF(λ)=−1βlne−βWF(λ)

(50)

ΦR(λ)=−1βlne−βWR(λ)

(51)

両方向の過程を用いたPMFは式(46)のBAR法により算出した自由エネルギー変化ΔF_Bを用いて次式により算出する [18]．

ΦB(λ)=−1β lnnFe−βWF(λ)nF+ nRe−βWF(λτ)−ΔFBF+nReβWR(λ)nF+ nReβWR(λ0)+ΔFBR

(52)

なお，ここで示したPMFはkが十分に大きければ自由エネルギーとして近似可能である [14]．

詳細ゆらぎの定理に基づく自由エネルギー計算がMDシミュレーションに利用されている例としては，生体分子に対するSMD法の適用 [14, 19,20,21,22,23]や，溶媒中の溶質の会合や溶媒和のシミュレーション [24, 25]がある．本稿では適用例として真空中のグラフェンシート (C₅₄H₁₈) のひずみ (Figure 3) に対する自由エネルギー計算の結果を示す．計算プログラムはDCDFTBMD [3]を用い，計算はDFTB3-D3(BJ)/3ob [26,27,28]レベルで行った．NVTアンサンブルでAndersen熱浴により300 Kに保ち，刻み幅は0.5 fsとした．SMD法のパラメータはk=2400 kcal mol⁻¹ Å⁻²，|v|=10⁻³, 10⁻⁵ Å fs⁻¹とした．順過程および逆過程のトラジェクトリ数はn_F=n_R=10とした．中央の炭素原子集団と外縁の炭素原子集団の重心間距離をλとした．

Figure 3.

Structure of (a) planar and (b) strained graphene sheet. The central and outer carbon atom groups for reaction coordinate are represented by black and tan ball models, respectively.

構造最適化後のグラフェンシートの構造に基づきλ₀=0.0 Åとした．200 psの平衡化の後，λ_τ=2.2 Åとして順過程のSMD計算を行った．再び200 psの平衡化の後，λ_τからλ₀の逆過程のSMD計算を行った．各過程のシミュレーション時間は|v|=10⁻³ Å fs⁻¹の場合は2.20 ps, 10⁻⁵ Å fs⁻¹の場合は220 psとした．

Figure 4 (a), (b) に各|v|に対するWを示す．W_FとW_Rは増加の傾向やλ₀とλ_τ間の絶対値において類似したλ依存性を示した．|v|=10⁻³ Å fs⁻¹の場合，トラジェクトリごとに異なるWが得られた．一方で，|v|=10⁻⁵ Å fs⁻¹の場合Wはすべてのトラジェクトリでよく一致した．したがって，前者は不可逆な過程，後者は限りなく可逆に近い過程である．

Figure 4.

((a) and (b)) Work obtained from forward process (red dotted lines) and reverse process (blue dashed lines) by straining graphene sheet with (a) |v|=10⁻³ Å fs⁻¹ and(b) |v|=10⁻⁵ Å fs⁻¹. ((c) and (d)) PMF calculated for forward process (red dotted line), reverse process (blue dashed line), and bidirectional method (purple solid line) with (c) |v|=10⁻³ Å fs⁻¹ and (d) |v|=10⁻⁵ Å fs⁻¹.

Table 1に式(38)のJarzynski等式および式(46)のBAR法により得られた始状態と終状態間の自由エネルギー変化ΔFを示す．|v|=10⁻⁵ Å fs⁻¹の場合，各手法から算出されたΔFは0.4 kcal mol⁻¹以内の範囲で一致した．|v|=10⁻³ Å fs⁻¹の場合，それぞれの手法について|v|=10⁻⁵ Å fs⁻¹の結果と比較すると，順過程は13.0 kcal mol⁻¹の差，逆過程は10.3 kcal mol⁻¹の差があったが，両方向の過程は1.3 kcal mol⁻¹の差であった．

Table 1. Free energy change of planar and strained graphene sheet obtained from SMD simulations. F and R represent the results of Jarzynski's equality for forward and reverse processes. B represents the results of the BAR method using bidirectional processes.

Calculation method	ΔF / kcal mol⁻¹
Calculation method	\|v\|=10⁻⁵ Å fs⁻¹	\|v\|=10⁻³ Å fs⁻¹
F	169.1	182.1
R	168.7	158.4
B	168.9	170.2

Figure 4 (c), (d) にFigure 4 (a), (b) で示したWおよびΔF_Bから式(50)-(52)により算出したPMFを示す．ΔF_BはTable 1に示したBAR法により算出した自由エネルギー変化であり，Ф_B(λ_τ)に対応する．Ф_F，Ф_R，Ф_Bは|v|=10⁻³ Å fs⁻¹において概ね一致し，|v|=10⁻⁵ Å fs⁻¹においてさらによく一致した．また|v|=10⁻³ Å fs⁻¹におけるФ_Bは|v|=10⁻⁵ Å fs⁻¹におけるФ_F，Ф_R，Ф_Bとよく一致した．

6 おわりに

本稿では非平衡を記述する理論として詳細ゆらぎの定理から出発し，非平衡過程から自由エネルギーを計算する手法としてJarzynski等式とBAR法を紹介した．詳細ゆらぎの定理はHamiltonian系において広く成り立つ定理であり，Jarzynski等式およびBAR法が適用可能な範囲もこれに準ずる．適用例では，ゆらぎの定理を用いることで不可逆な過程のシミュレーションから自由エネルギーを計算した．また，過程が不可逆な場合でも，過程が可逆に近い場合と同程度の精度の自由エネルギー計算を達成した．

参考文献

Corresponding author

Register with J-STAGE for free!