(57)【要約】      （修正有）【課題】分子集合体に溶質分子が溶解している系における、スケールされた非結合性相互作用を用いた分子動力学シミュレーションにおいて、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を考慮した、高精度な自由エネルギー計算方法を提供する。【解決手段】コンピュータに、原子の座標、スケール値など初期条件を入力し、（１）分子内力・分子間力を計算し、（２）カットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→rcを計算し、（３）カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を計算し、（４）ＰvdW0→rcとＰvdWrc→∞との和からＰvdW0→∞計算し、（５）微小時間Δt後の原子の座標、速度を数値積分によって更新し、（５）の更新値を初期値として（１）〜（５）を再度実行させ、これを規定回数繰り返すことを含む。【選択図】図１
2018049608.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを利用した、分子集合体に溶質分子が溶解している系の分子動力学シミュレーションにおいて、下記数式（３）のように、
（i）溶質分子の分子内vdW相互作用Vvdwsol-intra(r, λsol-intra)と、
（ii）分子集合体に属する分子の分子内vdW相互作用Vvdwmol-intra(r, λmol-intra)と、
（iii）溶質分子‐溶質分子間vdW相互作用Vvdwsol-sol(r, λsol-sol)と、
（iv）分子集合体に属する分子の分子間vdW相互作用Vvdwmol-mol(r, λmol-mol)と、
（v）溶質分子‐分子集合体に属する分子の分子間vdW相互作用Vvdwsol-mol(r, λsol-mol)
との和（ただし0≦λ≦1）で表され、非結合性相互作用のスケール値λを含むスケールファンデルワールス（vdW）相互作用ポテンシャルVvdW(r, λ)を用いた場合の自由エネルギー計算法であって、
【数１】
2018049608.tif 000028
コンピュータに、
・前記溶質分子と分子集合体に含まれる原子の座標r 
および、必要に応じ
・原子の速度v
の初期値、
ならびに
・結合情報
・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
・原子の質量m
・相互作用計算のカットオフ半径rc
・中間状態を定義する非結合性相互作用のスケール値λi（λ1、…λi…、λnint、0〜1を中間状態の数nintで分割した値）
を入力した後、
非結合性相互作用のスケール値λiをもつ中間状態iそれぞれにおいて、コンピュータ上で少なくとも、
（ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、
下記数式（１）で示した全ポテンシャルエネルギーVtotal(λi) を原子の座標rで微分した式を用いて、前記分子集合体および前記溶質の分子内力(数式（１）第１項から第４項の微分の和)および分子間力（数式（１）第５項と第６項の微分の和）を計算するステップ；
【数２】
2018049608.tif 000029
（ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
（ステップ２）ステップ１で計算した分子間力の一つであるvdW力（前記数式（１）第５項の微分）と、
原子の座標rからビリアル定理を用いて圧力ＰvdW0→ｒｃを計算するステップ；
（ステップ３）結合情報、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λiおよび相互作用計算のカットオフ半径ｒｃから、
下記数式(４Ａ)および（４Ｂ）を用いて圧力ＰvdWrc→∞を計算するステップ；
【数３】
2018049608.tif 000030
（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎα、Ｎβはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
（ステップ４）ステップ２およびステップ３の実行後に、下記数式（５）によりvdW相互作用から生じる圧力ＰvdW0→∞を計算するステップ；
【数４】
2018049608.tif 000031
（ステップ５）ステップ４の実行後に、原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
・原子の座標r
・原子の速度ｖ
を数値積分によって更新するステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含む、自由エネルギーの計算方法。
【請求項２】
ステップ１〜５のいずれかのステップの前または後に、
（ステップＡ）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λiを用いて、
前記数式（１）における全ポテンシャルエネルギーVtotal(λi)の中で、VvdW(r,λi)以外のポテンシャルエネルギーを計算するステップ：
（ステップＢ）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、相互作用計算のカットオフ半径ｒｃを用いて、VvdW(r,λi)における、原子間距離が0からカットオフ半径rc以内のvdW相互作用ポテンシャルVvdW(0→rc)を計算するステップ：
（ステップＣ）結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、相互作用計算のカットオフ半径rc、下記数式（１１）を用いて、VvdW(r,λi)における、カットオフ半径rcから無限遠までのvdW相互作用ポテンシャルVvdW(rc→∞)を計算するステップ：
【数５】
2018049608.tif 000032
（ステップＤ）ステップＡ、ステップＢおよびステップＣのすべてのステップの実行後に、ステップＢおよびステップＣで計算したVvdW(0→rc)とVvdW(rc→∞)と下記数式（８）よりVvdW(r,λi)を求め、VvdW(r,λi)とステップＡで計算した各エネルギー成分のポテンシャルエネルギーと前記数式（２）より全ポテンシャルエネルギーVtotalを計算するステップ：
【数６】
2018049608.tif 000033
を実行することを含む請求項１に記載の自由エネルギー計算方法。
【請求項３】
中間状態を定義する非結合性相互作用のスケール値λiを、スケール値を除き全く同じ入力情報を持つnint個のレプリカに割り当て、
それぞれのレプリカにおいて独立にステップＡ〜Ｄおよびステップ１〜５を実行後、
（ステップ６）ｉ番目のレプリカと隣り合う、スケール値λi+1を持つｉ＋１番目のレプリカに、
ステップＤで求めた全ポテンシャルエネルギーVtotalとメトロポリスの方法による判定基準を用い、
基準を満たす場合はλiとλi+1との交換を行うステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、
ステップ６でメトロポリスの判定基準を満たした場合は交換されたλを初期値にステップＡ〜Ｄおよびステップ１〜５を再度実行させ、これを繰り返すことを含む、請求項２に記載の自由エネルギーの計算方法。
【請求項４】
数式（１）におけるスケールされたvdW相互作用ポテンシャルVvdW(r, λ)が数式（１０）に示したDreiding型レナードジョーンズポテンシャルであり、数式（１２Ａ）および（１２Ｂ）を用いてカットオフ半径rcから無限遠までのvdW相互作用ポテンシャルVvdW(rc→∞)を、数式（４Ａ）、数式（１４Ａ）および（１４Ｂ）を用いてＰvdWrc→∞を計算することを特徴とする、請求項２または３に記載の自由エネルギーの計算方法。
【数７】
2018049608.tif 000034
（ここで注目する二つの原子をj,kとし、δはエネルギーと力の距離依存性を調節するパラメータ、D０、R０はvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータ、r2jkはij間の距離とし、0≦λi≦1、δ≧0ｍ２とする。）
【数８】
2018049608.tif 000035
【数９】
2018049608.tif 000036
【請求項５】
非結合性相互作用のスケール値λiをもつ中間状態iの計算のステップＡ〜Ｄにおいて、
スケール値λiのときの全ポテンシャルエネルギーVtotal(λi)に加え、同じ分子配置を用いたスケール値λi+1のときの全ポテンシャルエネルギー<Vtotal(λi+1)>iを求め、
自由エネルギー摂動法による自由エネルギー計算を実施することを特徴とする、
請求項２〜４のいずれか記載の自由エネルギーの計算方法。
【請求項６】
非結合性相互作用のスケール値λiをもつ中間状態iの計算のステップＡ〜Ｄにおいて、
スケール値λiのときの全ポテンシャルエネルギーVtotal(λi)、スケール値λiのときの分子配置を用いたスケール値λi+1における全ポテンシャルエネルギー<Vtotal(λi+1)>i、スケール値λiのときの分子配置を用いたスケール値λi-1における全ポテンシャルエネルギー<Vtotal(λi-1)>iを求め、
ベネット受容比法による自由エネルギー計算を実施することを特徴とする、
請求項２〜４のいずれか記載の自由エネルギーの計算方法。
【請求項７】
非結合性相互作用のスケール値λiをもつ中間状態iの計算のステップＡ〜Ｄにおいて、スケール値λiのときの分子配置を用いたスケール値λk(k=1,…, i,…, nint)における全ポテンシャルエネルギー<Vtotal(λk)>iを求め、
マルチステイトベネット受容比法による自由エネルギー計算を実施することを特徴とする、
請求項２〜４のいずれか記載の自由エネルギーの計算方法。
【請求項８】
分子集合体として高分子、溶質として低分子を用いることを特徴とする請求項１〜７のいずれか記載の自由エネルギーの計算方法。
【請求項９】
非結合性相互作用のスケール値をλsol-intra＝λmol-intra＝λsol-sol＝λmol-mol＝１、λsol-mol=λiとすることを特徴とする請求項１〜８のいずれか記載の自由エネルギーの計算方法。
【請求項１０】
請求項１〜９のいずれか記載の自由エネルギーの計算方法をコンピュータに実行させるためのプログラム。
【請求項１１】
請求項１０に記載のプログラムをコンピュータに読み取り可能に記録した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コンピュータを利用した分子動力学シミュレーションにおける、分子集合体に溶質分子が溶解している系の自由エネルギーの計算方法に関する。
【背景技術】
【０００２】
分子動力学シミュレーションは、分子や分子集合体の高次構造や物性を詳細に解析するための手法として広く用いられる。
【０００３】
分子動力学シミュレーションより得られる、分子集合体に関する重要な物性の一つとして、分子集合体へ溶質が溶解するときの自由エネルギーが挙げられる。
【０００４】
分子集合体が低分子液体である場合、溶質が低分子液体へ溶解するかの議論をする上で、自由エネルギーは非常に重要な指標となる。また低分子液体として水と１−オクタノールを選択し、これらへの溶質の溶解自由エネルギーを計算すれば、溶質の水／１−オクタノール分配係数をシミュレーションより求めることができる（非特許文献１）。水／１−オクタノール分配係数は、薬物動態や環境動態の指標として用いられ産業応用上も非常に重要な物理量である。溶質として実験が困難な分子や、新規候補分子を用い、水／１−オクタノール分配係数を予測したいという需要が大きいため、シミュレーションによる予測も活発に行われている（非特許文献２）。
【０００５】
また分子集合体が高分子分離膜、溶質分子が分離物質の場合、分離物質の高分子膜への溶解自由エネルギーを表し、分離物質の膜透過性を議論するに当たり重要な物理量となる。この自由エネルギーから溶解度係数を求めることにより、高分子分離膜における物質透過性をシミュレーションより予測することができる（非特許文献３）。なお溶解度係数には、自由エネルギーの誤差が自然対数の底のべき乗で影響を及ぼすため、高分子分離膜系で膜透過性を議論する場合は非常に高精度な溶解自由エネルギー計算が必要とされる。
【０００６】
分子動力学シミュレーションを用いた代表的な自由エネルギー計算方法は、熱力学的積分法、自由エネルギー摂動法など、分子間相互作用をスケールした中間状態を考える錬金術的自由エネルギー計算法が挙げられる。（非特許文献４、５）これらの手法を溶質の分子集合体への溶解の問題に適用する場合、溶質‐分子集合体に属する分子の分子間非結合相互作用を、通常の相互作用から完全に相互作用が消滅した状態まで、徐々に非結合相互作用をスケールした中間状態を計算していかなければならない。
【０００７】
ここで分子動力学シミュレーションにおいて得られる、自由エネルギーを含む物理量は、理論上は取りうる分子配置を全てサンプリングした場合に真の値が得られる。そのため物理量の精度は分子の配置のサンプリング数に依存する。
【０００８】
前記の通り、高分子分離膜の物質透過性予測には非常に高精度な溶解自由エネルギー計算が必要とされる。しかし高分子の遅い分子運動のため、高分子膜中の分離物質は通常の分子動力学シミュレーションの計算時間ではほとんど移動しない。そのため、高分子は、通常の非結合相互作用に近い大きさの相互作用を持つ中間状態では、分離物質の運動性が低く分子配置を効率的にサンプリングできない。そこで、非結合相互作用のスケール値を時間と共に増減させる方法（特許文献１）や、レプリカと呼ばれる全く同じ系を発生させ、レプリカ毎に異なる非結合相互作用のスケール値を設定し、レプリカ毎に独立に分子動力学シミュレーションを行い、一定時間ごとにある確率でスケール値のみ交換するというハミルトニアンレプリカ交換分子動力学法（非特許文献６）といったスケール値を変化させるサンプリング方法を組み合わせる必要がある。
【０００９】
他にもタンパク質といった生体高分子系でも盛んに自由エネルギー計算が行われている。例えば分子集合体をタンパク質と水、溶質を薬剤分子とすると、薬剤のタンパク質レセプターへのドッキング時の自由エネルギーを求めることができ、薬剤設計において重要な指標となる。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特開２０１１−１２３８７４号公報
【非特許文献】
【００１１】
【非特許文献１】N. M. Garrido et. al, AIChE Journal, 58, 6, 1929(2012)
【非特許文献２】C. Nieto-Draghi et. al, Chem. Rev. 115, 13093(2015)
【非特許文献３】Y. Tamai, H. Tanaka, and K. Nakanishi, Macromolecules, 28, 2544(1995)
【非特許文献４】岡崎進、コンピュータシミュレーションの基礎、  化学同人（２０００）
【非特許文献５】山下雄史、分子シミュレーション研究会会誌“アンサンブル”、Ｖｏｌ. １７、Ｎｏ. ２(２０１５)
【非特許文献６】G. Bussi, Mol. Phys. 112, 379(2014)
【発明の概要】
【発明が解決しようとする課題】
【００１２】
分子動力学シミュレーションにおいて、中間状態λiにおける全ポテンシャルエネルギーVtotal(λi)は、下記数式（１）のように表される。
【００１３】
【数１】
2018049608.tif 000003
【００１４】
（ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
ファンデルワールス（vdW）相互作用エネルギーVvdW(0→∞)は、カットオフ半径内のvdW相互作用エネルギーVvdW(0→rc)とカットオフ半径から無限遠のvdW相互作用エネルギーVvdW(rc→∞)の和で表される。
【００１５】
ここで中間状態の分子配置を求めるため、スケールした非結合性相互作用を使用した場合、カットオフ半径から無限遠のvdW相互作用エネルギーVvdW(rc→∞)を計算できず、中間状態λiにおける全ポテンシャルエネルギーVtotal(λi)も正確に求めることができなかった。VvdW(rc→∞)の寄与は小さいため、考慮されない場合も多いが、高分子分離膜系の溶解度係数のように非常に高精度な自由エネルギー計算が求められる場合は必要である。
【００１６】
カットオフ半径から無限遠のvdW相互作用を無視した場合、本当に問題になるのはポテンシャルエネルギーではなく、圧力を正しく求められなくなることである。ここで、分子動力学シミュレーションにおいて、系の全圧力Ptotalは、下記数式（２）のように表される。
【００１７】
【数２】
2018049608.tif 000004
【００１８】
（ただし、Pintraは分子内相互作用から生じる圧力、PvdW0→∞はvdW相互作用から生じる圧力、Pcoulomb0→∞はクーロン相互作用から生じる圧力、PKは原子の運動エネルギーより生じる圧力である。）
vdW相互作用から生じる圧力ＰvdW0→∞は、カットオフ半径内のvdW相互作用から生じる圧力ＰvdW0→rcとカットオフ半径から無限遠のvdW相互作用から生じる圧力ＰvdWrc→∞の和で表される。
【００１９】
中間状態の分子配置を分子動力学シミュレーションより求めるため、スケールした非結合性相互作用を使用した場合、ポテンシャルエネルギーと同様に、正確な圧力を計算することができない。
【００２０】
正確な圧力が計算できない場合、分子動力学シミュレーションにおいて標準状態を再現するために行う、圧力を一定に制御した状態のシミュレーションを正しく行えず、シミュレーションセルの大きさを過大評価あるいは過小評価してしまう。結果として、分子間の距離が通常よりも近くまたは遠くなるため、正確な相互作用エネルギーが得られず、自由エネルギーを過大または過小評価することになる。このようにスケールされたvdW相互作用自体のカットオフ半径から無限遠の寄与よりも、圧力を正確に計算できなくなり、密度が変化し、正確でない分子間距離で相互作用を計算してしまうことに問題があった。
【００２１】
特に、高分子膜中の分離物質系のように高分子と低分子の混合系の場合、圧力を正確に計算できないという影響が大きく、自由エネルギーを高精度に計算できないという問題があった。
【００２２】
以上述べたように、分子集合体への溶質分子の溶解自由エネルギーを求めるためには、スケールされた非結合性相互作用を用いた分子動力学シミュレーションを実施する必要がある。しかし、スケールされた非結合相互作用を用いたときのカットオフ半径から無限遠のvdW相互作用から生じるポテンシャルエネルギーVvdW(rc→∞)および圧力ＰvdWrc→∞を求める方法が未だ提唱されていないため、高精度な自由エネルギー計算が行えなかった。そこで本発明では、上記問題を解決した、高精度な自由エネルギー計算方法を提供することを目的とする。
【課題を解決するための手段】
【００２３】
本発明は上記課題を解決するために、次の手段を採用するものである。すなわち、本発明は以下のとおりである。
【００２４】
コンピュータを利用した、分子集合体に溶質分子が溶解している系の分子動力学シミュレーションにおいて、下記数式（３）のように、
（i）溶質分子の分子内vdW相互作用Vvdwsol-intra(r, λsol-intra)と、
（ii）分子集合体に属する分子の分子内vdW相互作用Vvdwmol-intra(r, λmol-intra)と、
（iii）溶質分子‐溶質分子間vdW相互作用Vvdwsol-sol(r, λsol-sol)と、
（iv）分子集合体に属する分子の分子間vdW相互作用Vvdwmol-mol(r, λmol-mol)と、
（v）溶質分子‐分子集合体に属する分子の分子間vdW相互作用Vvdwsol-mol(r, λsol-mol)
との和（ただし0≦λ≦1）で表される、非結合性相互作用のスケール値λを含むスケールvdW相互作用ポテンシャルVvdW(r, λ)を用いた場合の自由エネルギー計算法であって、
【００２５】
【数３】
2018049608.tif 000005
【００２６】
コンピュータに、
・前記溶質分子と分子集合体に含まれる原子の座標r 
および、必要に応じ
・原子の速度v
の初期値、
ならびに
・結合情報
・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
・原子の質量m
・相互作用計算のカットオフ半径rc
・中間状態を定義する非結合性相互作用のスケール値λi（λ1、…λi…、λnint、0〜1を中間状態の数nintで分割した値）
を入力した後、
非結合性相互作用のスケール値λiをもつ中間状態iそれぞれにおいて、コンピュータ上で少なくとも、
（ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、
前記数式（１）で示した全ポテンシャルエネルギーVtotal(λi) を原子の座標rで微分した式を用いて、前記分子集合体および前記溶質の分子内力(数式（１）第１項から第４項の微分の和)および分子間力（数式（１）第５項と第６項の微分の和）を計算するステップ；
（ステップ２）ステップ１で計算した分子間力の一つであるvdW力（数式（１）第５項の微分）と、
原子の座標rからビリアル定理を用いて圧力ＰvdW0→ｒｃを計算するステップ；
（ステップ３）原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λiおよび相互作用計算のカットオフ半径ｒｃから、
下記数式(４Ａ)および（４Ｂ）を用いて圧力ＰvdWrc→∞を計算するステップ；
【００２７】
【数４】
2018049608.tif 000006
【００２８】
（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎα、Ｎβはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
（ステップ４）下記数式（５）によりvdW相互作用から生じる圧力ＰvdW0→∞を計算するステップ；
【００２９】
【数５】
2018049608.tif 000007
【００３０】
（ステップ５）原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
・原子の座標r
・原子の速度ｖ
を数値積分によって更新するステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含む、自由エネルギーの計算方法。
【図面の簡単な説明】
【００３１】
【図１】本発明における自由エネルギーの計算方法を示すフローチャートである。
【図２】図１の計算に、カットオフ半径から無限遠のvdW相互作用のポテンシャルエネルギーへの寄与を考慮するステップＡ〜Ｄを加えた自由エネルギーの計算方法を示すフローチャートである。
【図３】図２の計算に、ハミルトニアンレプリカ交換を行うステップ６を加えた自由エネルギーの計算方法を示すフローチャートである。
【発明を実施するための形態】
【００３２】
図１、図２、図３に本発明の計算方法のフローを示す。以下、本発明の好ましい実施形態を、図１、図２、図３のフローに従って説明する。なお図１はカットオフ半径から無限遠のvdW相互作用の圧力への寄与を考慮する自由エネルギー計算であり、図２は図１の計算に加えカットオフ半径から無限遠のvdW相互作用のポテンシャルエネルギーへの寄与を考慮する自由エネルギー計算、図３は図２の計算に加え、ハミルトニアンレプリカ交換法を組み合わせた自由エネルギー計算法である。
【００３３】
本発明は分子動力学シミュレーションにおける、分子集合体に溶質分子が溶解している系の自由エネルギーの計算方法である。ここで分子集合体とは、有機溶媒や水といった溶液、合成高分子、生体高分子またはそれらの混合系など分子が凝集している系ならばどのようなものでも良い。凝集しているとは、気体状態でないことを示し、より具体的には、分子動力学シミュレーションにおいて、密度が0.3g/cm3以上で分子が基本セル内に存在する状態を意味するものとする。
【００３４】
合成高分子とは、ポリエチレン、ポリエチレンテレフタラートなど人工的にモノマーを合成して得られる高分子を指し、重合度が10以上、分子量が1000以上のものとする。重合度が10程度、あるいは分子量が1000程度の分子は、一般的には高分子ではなくオリゴマーと呼ばれる大きさであるが、本明細書においては合成高分子に含めて記載するものとする。生体高分子とは、タンパク質や核酸など生体由来の高分子を指す。
【００３５】
本発明の計算方法においては、まず、コンピュータに、計算に必要な諸条件を入力する。計算に必要な諸条件を入力するに際しては、コンピュータに読み取り可能な形態で、分子集合体に属する分子と溶質分子に含まれる原子の座標rの初期値、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、原子の質量ｍ、相互作用計算のカットオフ半径ｒｃ、中間状態数nintを入力する。ここで、後述するハミルトニアンレプリカ交換を実施する場合、中間状態数nintはレプリカ数に等しい。また、必要に応じて原子の速度ｖの初期値や、i番目の中間状態またはi番目のレプリカにおける非結合相互作用のスケール値λiを入力することもできる。しかし入力しなければ、速度はマクスウェル分布に従うよう計算された原子の速度が初期速度として与えられるよう構成すればよい。非結合相互作用のスケール値は０〜１の範囲を取るため、簡便に中間状態数nintで１を等分した値やプログラム内でnintの数に応じてあらかじめ設定してある値を用いても良い。
【００３６】
入力方法は特に限定されないが、一般に入手可能な分子モデル作成ソフトウェアやテキストエディタなどが利用でき、当該情報が記述されたファイルをコンピュータに読み取らせ、メモリに格納することが好ましい。
【００３７】
次に中間状態数nint分のメモリ領域を確保し、読み取った情報をコピーする。このときハミルトニアンレプリカ交換を行わない場合は、nint分のメモリ領域を確保し、同時に計算を行う必要は必ずしもない。中間状態は独立に計算可能であるため、図１または図２におけるフローチャートの中間状態を一つずつ計算していっても構わない。このとき中間状態1から中間状態nintまで、前の中間状態の最終座標を初期配置として順に計算していくことも、同じ初期配置を用いばらばらに計算していくことも、一つの中間状態で十分なサンプリングが行える時間の計算を行っているならばどちらでも構わない。十分なサンプリングが行える時間は系の構成によって異なるが、分子集合体に溶質分子が溶解している系が溶液である場合は、計算時間は100ピコ秒以上が好ましく、300ピコ秒以上がより好ましく、1000ピコ（1ナノ）秒以上がさらに好ましい。
【００３８】
中間状態（レプリカ）数nintは、錬金術的自由エネルギー計算法における計算精度に大きく影響する。スケール値λiが0〜1まで変化する中で、中間状態はなるべく細かく定義すると自由エネルギーの計算精度が上昇することが知られているが、計算コストも中間状態数に比例する。そのため好ましい精度かつ好適な計算時間で自由エネルギーを得るには、nintで１を等分したλiを割り当てる場合、6以上50以下が好ましく、10以上35以下がより好ましく、15以上25以下がさらに好ましいことを本発明者らは検討の結果見出した。
【００３９】
本発明では、i番目の中間状態において、図１、図２、図３に示したステップＡ〜Ｄ、ステップ１〜６を必要に応じて実施する。なおステップＡ〜Ｄは、ステップ１〜５のどのステップの前後においても実行することができる。またステップＡ〜Ｄに関しては、ステップＡ、ステップＢ、ステップＣの実行順序は問わず、ステップＡ、ステップＢおよびステップＣのすべてのステップの実行後にステップＤを実行すればよい。ステップ１〜５に関しても、ステップ２とステップ３の実行順序は問わず、ステップ２およびステップ３の実行後にステップ４を実行し、最後にステップ４の実行後にステップ５を実行すればよい。ステップ６に関してはステップＡ〜Ｄおよびステップ１〜５の実行後に実行する。なお、以下、本明細書においては、ステップＡ〜Ｄをこの順に実行した後、ステップ１〜６をこの順に実行した場合を例を前提に説明するものとする。
【００４０】
〔ステップＡ〕
ステップＡにおいては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λiを用いて、前記数式（１）における全ポテンシャルエネルギーVtotal(λi)の中で、VvdW(r,λi)以外のポテンシャルエネルギー{Vtotal (λi)‐VvdW(r,λi)}を計算し、メモリに格納する。
【００４１】
ここで、数式（１）の右辺第１項Vbond(R)は結合長、第２項Vangle(θ)は結合角、第３項Vtorsion(φ)は二面角、第４項Vinversion(ω)は反転、第５項Vvdw(r,λ)はｖｄＷ、第６項Vcoulomb(r,λ)はクーロンのそれぞれのポテンシャルエネルギーを定義するポテンシャル関数を表わす。なお、数式（１）の右辺第１項から第４項を総称して分子内ポテンシャル関数と呼び、右辺第５項、第６項をまとめて分子間ポテンシャル関数と呼ぶ。また、分子内ポテンシャル関数、分子間ポテンシャル関数を原子の位置で微分した式を用いて計算した力をそれぞれ分子内力、分子間力と呼ぶ。
【００４２】
数式（１）における各ポテンシャル関数の具体的な表式は下記数式（６）で表される。
【００４３】
【数６】
2018049608.tif 000008
【００４４】
数式（６）の右辺第５項、第６項はスケールされた非結合性相互作用を表し、非結合性相互作用のスケール値λi=1のときスケールされていない一般的な分子動力学シミュレーションで使用するポテンシャル関数に等しい。
【００４５】
前記数式（６）中の結合長、結合角、二面角、反転おけるそれぞれの力の定数Kbond、Kangle、Ktorsion、Kinversionおよびそれぞれの平衡位置R0、θ0、φ0、ω0には、DREIDING、AMBER、hybrid/UA[A. L. Frischknecht、J. G. Curro、Macromolecules、2003、36、p.2122]、桑島らのパラメータ〔桑島聖、野間祐之、逢坂俊郎、“非晶高分子のＭＤ計算用力場パラメータについてｉｉ：ポリオレフィン”、第４回計算化学シンポジウム&ＣＣＳ展  講演要旨集、日本化学プログラム交換機構（１９９４）〕等の公知のパラメータを用いることができる。
【００４６】
計算精度について考慮した場合、以下に述べる方法によってポテンシャルパラメータを決定することがさらに好ましい。
【００４７】
前記数式（１）および数式（６）の第１項から第４項で表わされる分子内ポテンシャル関数のパラメータは量子化学計算によって決定することが好ましい。量子化学計算としては、HF法〔Ａ．ザボ，Ｎ．Ｓ．オストランド，「新しい量子化学」，東京大学出版会（１９８７）．〕を用いることが好ましく、B3LYP法〔A. D. Becke, J. Chem. Phys., 1993, 98, p.5648／C. Lee, W. Yang, R. G. Parr, Phys. Rev. B, 1988, 37, p.785／B. Miehlich, A. Savin, H. Stoll,H. Preuss, Chem. Phys. Lett., 1989, 157,p.200.〕を用いることがさらに好ましい。
【００４８】
HF法やB3LYP法を用いる場合、１電子軌道を展開する基底関数を指定する必要がある。計算時間と計算精度を考慮すると、6-31G(d)基底関数、6-31G(d，p)基底関数、6-31+G(d，p)基底関数〔W. J.Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 1972, 56, p.2257， P. C. Hariharan, J.A. Pople, Theoret. Chimica Acta, 1973, 28, p.213.〕が好適に用いられる。
【００４９】
前記数式（１）および数式（６）の第６項Vcoulomb(r,λ)のクーロンポテンシャルのパラメータ、すなわち電荷については、QEq法〔A.K.Rappe,W.A.Goddard III, J. Phys. Chem. 95, 3358-3363(1991).〕等の学術的に信頼性の高い手法を用いて決定することが好ましく、量子化学計算から静電ポテンシャル（ESP）を求め、ESPフィッティングして電荷を決定することがさらに好ましい。静電ポテンシャルを求めるための量子化学計算としては、HF法やB3LYP法を好適に用いることができ、基底関数としては6-31G(d)基底関数、6-31G(d,p)基底関数、6-31G+(d,p)基底関数が好適に用いられ、さらに好ましくは、aug-cc-pVDZ基底関数、aug-cc-pVTZ基底関数、aug-cc-pVQZ基底関数〔T.H. Dunning, Jr.， J. Chem. Phys.，90，1007 (1989).〕が用いられる。ESPフィッティングして電荷を求める手法としては、MK法〔B. H. Besler, K. M. Merz Jr, P. A. Kollman, J. Comput. Chem., 1990 11, p.431.〕やCHelpG法〔C. M. Breneman, K. B. Wiberg, J. Comput. Chem., 1990, 11, p.361.〕といった方法を好適に用いることができる。
【００５０】
以上のような量子化学計算を実行するために、有償または無償で入手可能な種々の量子化学計算プログラムパッケージを利用することができる。例えば、Gaussian （登録商標）09 Revision C.01, M. J. Frisch他, ガウシアン社（Gaussian, Inc.）, Pittsburgh PA, 2002.〕、GAMESS〔M. W. Schmidt 他, J. Comput. Chem., 1993, 14, p.1347.〕といったプログラムを好適に用いることができる。
【００５１】
分子集合体に溶質分子が溶解している系の分子動力学シミュレーションにおいて、前記数式（１）および数式（６）の第５項、第６項における非結合性相互作用は、下記数式（７）のように
（i）溶質分子の分子内非結合性相互作用Vsol-intra(r, λsol-intra)と、
（ii）分子集合体に属する分子の分子内非結合性相互作用Vmol-intra(r, λmol-intra)と、
（iii）溶質分子‐溶質分子間非結合性相互作用Vsol-sol(r, λsol-sol)と、
（iv）分子集合体に属する分子の分子間非結合性相互作用Vmol-mol(r, λmol-mol)と、
（v）溶質分子‐分子集合体に属する分子の分子間非結合性相互作用Vsol-mol(r, λsol-mol)
との和で表される。
【００５２】
【数７】
2018049608.tif 000009
【００５３】
非結合性相互作用のスケール値λの範囲として0≦λ≦1を用いることができ、(i）から（v）それぞれのλにおいて別々のスケール値を用いることができる。しかし溶質分子が分子集合体に溶解するときの自由エネルギー計算には、中間状態を規定するスケール値をλiとすると、λsol-mol＝λi、λmol-intra＝λsol-sol＝λmol-mol＝１とする必要がある。なお、溶質分子を１分子のみとして計算を行う場合には、λsol-solは設定不要であり、水やメタンなど分子内非結合相互作用を持たない溶質分子の計算を行う場合には、λsol-intraの設定は不要である。一般的に分子動力学シミュレーションでは、分子内非結合相互作用は１−４相互作用以上を扱うため、１−２、１−３相互作用までしか存在しない水やメタンは、分子内非結合相互作用を計算しない。なお１−Ｘとは基準となる原子を１とし、共有結合でつながっている隣の原子を２、さらに２と共有結合でつながっている原子を３と順に番号付けしたときの相互作用ペアを表す。
【００５４】
また自由エネルギー計算においては、λmol-intra＝λsol-sol＝λmol-mol＝λsol-intra＝1とすることが好ましい。λsol-intra＝1とすることによって、分子内水素結合を持つ低分子の配座を維持しながら低分子の運動性を上昇させることができる。エチレングリコールのような分子内水素結合を持つ低分子は、分子内非結合性をスケールすると水素結合が弱まるため本来とるべき配座を再現しない。このように分子内非結合性相互作用が低分子のコンフォメーションに大きく寄与する場合、低分子内の非結合性相互作用をスケールしなければ低分子内の配座を維持しながら運動性を上昇させることができる。
【００５５】
さらにカットオフ半径rcは、計算時間と計算精度の観点から8Å≦rc≦12Åが好ましく、rc＝10Åが特に好ましい。
【００５６】
〔ステップＢ〕
スケールされたvdW相互作用である前記数式（６）第５項Vvdw(r,λi)は、下記数式（８）のようにカットオフ半径以内のvdW相互作用エネルギーVvdw(0→rc)(r,λi)とカットオフ半径から無限遠のvdW相互作用エネルギーVvdw(rc→∞)(r,λi)の和で表される。
【００５７】
【数８】
2018049608.tif 000010
【００５８】
数式（８）で示したカットオフ半径以内のvdW相互作用ポテンシャルVvdW(0→rc)を計算するステップＢにおいては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、相互作用計算のカットオフ半径ｒｃを用いてVvdW(0→rc)を計算し、メモリに格納する。
【００５９】
カットオフ半径以内のvdW相互作用エネルギーVvdw(0→rc)(r,λi)の具体的な関数形としては、下記数式（９）に示したレナードジョーンズポテンシャルやモースポテンシャル、バッキンガムポテンシャルが好ましく用いられる。
【００６０】
【数９】
2018049608.tif 000011
【００６１】
数式（９）において、m1、m2は型を決定する乗数であり、LJm1/m2型と表記するとき、LJ12/6型、LJ9/6型、LJ12/4型、LJ6/4型、LJ12/10型が好適に用いられる。さらにLJ12/6型には定数A,BやvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０の置き方により下記数式（１０）で示したDreiding型〔 S. L. Mayo 、B. D. Olafson 、W. A. Goddard III 、J. Phys. Chem.、1990、94、p.8897.〕やAmber型〔W. D. Cornell、P. Cieplak 、C. I. Bayly 、I. R. Gould 、K. M. Merz jr 、D. M. Ferguson 、D. C. Spellmeyer 、T. Fox 、J. W. Caldwell、P. A. Kollman 、J. Am. Chem. Soc.、1995、117、p.5179〕、OPLS型［W. L. Jorgensen、D. S. Maxwell、Julian Tirado-Rives、J. Am. Chem. Soc.、1996、118、p.11225］、CHARMM型［B. R. Brooks、R. E. Bruccoleri、B. D. Olafson、D. J. States、S. Swaminathan、M. Karplus、J. Comput. Chem.、1983、4、p.187］などいくつかの形式が存在し、どのLJ12/6型でも好適に用いることができる。
【００６２】
【数１０】
2018049608.tif 000012
【００６３】
前記数式（９）中のvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０には、前述のDREIDING、AMBER、hybrid/UA、桑島らのパラメータ等の公知のパラメータを好適に用いることができる。計算精度について考慮した場合、D０、R０については、密度と蒸発熱の実験値を再現するように決定することがより好ましい。
【００６４】
前記数式（９）、数式（１０）におけるエネルギーと力の距離依存性を調節するパラメータδについては、δ＝0ｍ2とし、かつ非結合性相互作用のスケール値λが非常に小さい値をとる場合に、前記数式（１）および数式（６）の第５項Vvdw(r,λ)で表わされるvdW相互作用が急激に変化し、数値積分の安定性が低下する。この傾向は原子間の距離が近距離となった場合に特に顕著になる。この問題を解決するために、非結合性相互作用のスケール値λの最小値を非常に小さな値に設定する場合には、δ＞0ｍ2とすることが好ましい。さらに好ましくは、δ≧5.0×10-21ｍ2である。一方、δの値が大きすぎると、スケール値が変化した際、vdW力の変化が大きくなりすぎるため、δ≦5.0×10-19ｍ2とすることが好ましい。このδの値としては、δ＝1.0×10-20ｍ2が好適に用いられる。
【００６５】
〔ステップＣ〕
前記数式（８）で示したカットオフ半径から無限遠のvdW相互作用ポテンシャルVvdW(rc→∞)を計算するステップＣにおいては、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、相互作用計算のカットオフ半径ｒｃおよび下記数式（１１）を用いてVvdW(rc→∞)を計算し、メモリに格納する。
【００６６】
【数１１】
2018049608.tif 000013
【００６７】
（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎα、Ｎβはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数を表す。）
数式（１１）におけるスケールされたvdW相互作用として前述の様々な関数形を用いることができる。これまで、非結合性相互作用のスケール値λiを含むスケールされたvdW相互作用ポテンシャルにおける、カットオフ半径から無限遠のvdW相互作用ポテンシャルVvdW(rc→∞)の計算方法が確立されていなかったが、本発明者らの検討の結果、数式（１０）のDreiding型レナードジョーンズポテンシャルを数式（１１）に代入して得られた解析解である下記数式（１２Ａ）および（１２Ｂ）からカットオフ半径から無限遠のvdW相互作用ポテンシャルVvdW(rc→∞)を計算する方法が好ましいことが見出された。
【００６８】
【数１２】
2018049608.tif 000014
【００６９】
〔ステップＤ〕
ステップＤは、ステップＡ、ステップＢおよびステップＣのすべてのステップの実行後に、前記数式（１）で示した全ポテンシャルエネルギーVtotalを計算するステップである。ステップＤにおいては、前記数式（８）により、ステップＢおよびステップＣで計算したカットオフ半径内のvdW相互作用ポテンシャルVvdW(0→rc)とカットオフ半径から無限遠のvdW相互作用ポテンシャルVvdW(rc→∞)の和から全vdW相互作用ポテンシャルVvdWを求め、ステップＡで計算した{Vtotal (λi)‐VvdW(r,λi)}とVvdWとの和より全ポテンシャルエネルギーVtotalを計算し、メモリに格納する。
【００７０】
なお本発明において、数式（６）の右辺第６項であるクーロン相互作用の計算には、右辺第６項を展開し、無限遠のクーロン相互作用の寄与まで考慮できる方法であるEwald法やParticle Mesh Ewald(PME)法、Particle-Particle Particle-Mesh Ewald(PPPME)法、Fast Multipole Method(FMM)などを用いることが好ましい。
【００７１】
ステップＡ〜Ｄを実施した場合、カットオフ半径から無限遠のvdW相互作用ポテンシャルの寄与を考慮した全ポテンシャルエネルギーVtotalを計算することができる。後述の自由エネルギー計算では、全ポテンシャルエネルギーVtotalを用いるため、自由エネルギーにもカットオフ半径から無限遠のvdW相互作用ポテンシャルの寄与が考慮された、より高精度な値を求めることできる。
【００７２】
〔ステップ１〕
分子集合体および溶質分子の分子内力(数式（１）第１項から第４項の微分の和)および分子間力（数式（１）第５項と第６項の微分の和）を計算するステップ１においては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λi、ならびに数式（１）で示した全ポテンシャルエネルギーVtotalを原子の座標rで微分した式を用いて分子内力・分子外力を計算し、メモリに格納する。なお数式（１）における具体的な関数の形やパラメータは、ステップＡの項に記したものが好適に用いることができ、ステップＡを計算する場合は、ステップ１においても同じ関数形と同じパラメータを用いなければならない。
【００７３】
〔ステップ２〕
カットオフ半径内のvdW相互作用から生じる圧力ＰvdW0→ｒｃを計算するステップ２では、ステップ１で計算したvdW力（数式（１）第５項の微分）と、原子の座標rとから、下記数式（１３Ａ）および（１３Ｂ）で表されるビリアル定理を用いて圧力ＰvdW0→ｒｃを計算する。
【００７４】
【数１３】
2018049608.tif 000015
【００７５】
（ここでFvdWは注目する原子i、j間に働くvdW力、l, m(=1、2、3)はベクトルの成分を表す。）
〔ステップ３〕
カットオフ半径から無限遠のvdW相互作用から生じる圧力ＰvdWrc→∞を計算するステップ３では、原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λiおよび相互作用計算のカットオフ半径ｒｃから前記数式(４Ａ)および（４Ｂ）を用いて圧力ＰvdWrc→∞計算する。
【００７６】
数式（４Ａ）の通り圧力ＰvdWrc→∞は対角項のみで、非対角項は０となる。なお数式（４Ｂ）におけるスケールされたvdW相互作用として前述の様々な関数形を用いることができる。
【００７７】
前述のように、これまで、非結合性相互作用のスケール値λを含むポテンシャルでの圧力ＰvdWrc→∞の計算方法が確立されていなかったが、本発明者らの検討の結果、数式（１０）のDreiding型レナードジョーンズポテンシャルを数式（４Ｂ）に代入して得られた解析解である下記数式（１４Ａ）および（１４Ｂ）と数式（４Ａ）から圧力Ｐrc→∞を計算する方法が好ましいことが見出された。
【００７８】
【数１４】
2018049608.tif 000016
【００７９】
〔ステップ４〕
ステップ４は、ステップ２およびステップ３の実行後に、vdW相互作用から生じる圧力ＰvdW0→∞を計算するステップである。ステップ４においては、数式（５）により、ステップ２およびステップ３で計算したカットオフ半径内のvdW相互作用から生じる圧力ＰvdW0→ｒｃとカットオフ半径から無限遠のvdW相互作用から生じる圧力Ｐrc→∞の和から圧力ＰvdW0→∞を計算する。
【００８０】
また数式（１）に、数式（４Ａ）および数式（４Ｂ）、数式（５）、数式（１３Ａ）および数式（１３Ｂ）より求めたＰvdW0→∞を代入することにより全圧力Ｐtotalを求めることができる。なおvdW相互作用から生じる圧力ＰvdW0→∞以外の圧力の詳細な表式もビリアル定理から導くことができ、参考文献［Glenn J. Martyna, Mark E. Tuckerman, Douglas J. Tobias, Michael L. Klein, Mol. Phys., 1996, 87, p.1117.］に記述されている。この文献の全圧力の表式はクーロン相互作用から生じる圧力Ｐcoulomb0→∞を、無限遠のクーロン相互作用の寄与を考慮できる方法であるEwald法やPME法を用いた場合の表式である。本発明においても無限遠のクーロン相互作用の寄与を考慮できる方法を用いることが好ましく、このような取り扱いをすることで無限遠までの相互作用を考慮した正確な全圧力Ｐtotalを得ることができる。
【００８１】
〔ステップ５〕
微小時間Δt後の原子の座標r、原子の速度ｖを数値積分によって更新するステップ５においては、原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の原子の座標と速度を数値積分によって計算し、更新する。圧力と温度を制御するために行うNPTアンサンブルを行う場合は、さらに熱浴の慣性量Q、熱浴の位置ζ、熱浴の速度ζv、barostatの慣性量W、基本セルベクトルh、基本セルベクトルの速度hv、分子内力、分子間力を用いて、微小時間Δt後の熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを数値積分によって計算し、更新する必要がある。
【００８２】
数値積分には、差分近似法や予測子−修正子法〔非特許文献４〕などの一般的な解法を好適に用いることができる。異なる時定数の力ごとに大きさの異なる時間刻みΔtを設定することのできる多重時間刻み数値積分（RESPA法）法〔G. J. Martyna、Mol. Phys.、1996、87、p.1117.〕を用いると、シミュレーション時間を短縮できるためさらに好ましい。
【００８３】
ステップ１〜５を実施すると、カットオフ半径から無限遠のvdW相互作用から生じる寄与が考慮された正確な圧力のもとでのシミュレーションを行うことができる。特にNPTアンサンブルでシミュレーションを行う場合は、正確な密度における分子配置が得られるため、その全ポテンシャルエネルギーVtotalから高精度な自由エネルギーを求めることができる。
【００８４】
〔ステップ６〕
ステップ６はハミルトニアンレプリカ交換を行う場合のみ、レプリカ交換を実施する間隔lrepに１回実施する。レプリカの交換方法は、一般的に次に示すメトロポリスの方法が好適に用いられる。
【００８５】
i番目のレプリカとスケール値λi+1を持つｉ＋１番目のレプリカにおける、全ポテンシャルエネルギーVtotal(λi or i+1)と、下記数式（１５）よりi番目とi+1番目のレプリカのスケール値を交換する確率αを計算する。０〜１までの疑似乱数をコンピュータ上に発生させ、交換確率以下の値である場合にλiとλi+1との交換を行う。
【００８６】
【数１５】
2018049608.tif 000017
【００８７】
（ここでkBはボルツマン定数、Tは系の絶対温度を表す。）
レプリカ交換を実施する間隔lrepは、本発明の最初に行われる諸条件の入力に際し初期条件として入力するか、またはあらかじめプログラム内に記述しておくことで規定する。
【００８８】
ここで微小時間Δtは、最も速い時定数の運動をとらえることができるよう十分小さな値でなければならない。一般的に最も速い時定数を持つ運動は水素原子を含む結合の伸縮運動である。この伸縮運動を精度良くとらえることができ、かつ計算時間の観点から微小時間Δtは、0.1フェムト秒以上0.3フェムト秒以下であることが好ましく、0.25フェムト秒がより好ましい。
【００８９】
lrepの設定に関しては、頻繁にレプリカ交換を実施すると、レプリカの構造が平衡にならず不安定な構造における分子の配置を多くサンプリングしてしまう。またレプリカ交換の頻度が少ないと、分子の配置のサンプリング効率が悪くなり、レプリカ交換を行っている意義が薄れる。そのためlrepは、1000回以上100000回以下（微小時間Δtが0.25フェムト秒のとき0.25ピコ秒以上25ピコ秒以下）が好ましく、2000回以上50000回以下（微小時間Δtが0.25フェムト秒のとき0.5ピコ秒以上12.5ピコ秒以下）がより好ましく、4000回以上20000回以下（微小時間Δtが0.25フェムト秒のとき2ピコ秒以上5ピコ秒以下）が特に好ましい。
【００９０】
ステップ６を実施する、すなわちハミルトニアンレプリカ交換を実施すると、ポリマーのように運動性の遅い分子集合体に溶解した溶質が、シミュレーション時間内に移動する空間が劇的に広くなる。そのため様々な分子配置がサンプリング可能になり、運動性の遅い分集合体に溶質が溶解した系において、高精度な自由エネルギーを求めることができる。  〔ステップＡ〜ステップＤおよびステップ１〜ステップ６の繰り返し〕
カットオフから無限遠のスケールされたvdW相互作用から生じる圧力のみを考慮したシミュレーションを実施する場合は、ステップ１〜５を、加えてカットオフから無限遠のスケールされたvdW相互作用ポテンシャルを考慮する場合は、ステップＡ〜Ｄおよびステップ１〜５を、さらにハミルトニアンレプリカ交換を実施する場合は、ステップＡ〜Ｄおよびステップ１〜６を繰り返し、あらかじめ規定した回数nstep実行する。
【００９１】
繰り返し回数nstepは、本発明の最初に行われる諸条件の入力に際し初期条件として入力するか、またはあらかじめプログラム内に記述しておくことで規定する。
【００９２】
繰り返し回数nstepは、系の平衡構造における構造の揺らぎが十分にシミュレーション時間中にあらわれる程度の回数が必要である。系の構造によって必要な繰り返し回数ｎは異なるが、分子集合体に低分子溶液、溶質分子に低分子を用いた場合であれば400000回（微小時間Δtが0.25フェムト秒のとき100ピコ秒）以上が好ましく、2000000回（微小時間Δtが0.25フェムト秒のとき500ピコ秒）以上がより好ましい。また分子集合体に高分子、溶質分子に低分子を用いた場合であれば8000000回（微小時間Δtが0.25フェムト秒のとき2ナノ秒）以上が好ましく、16000000回（微小時間Δtが0.25フェムト秒のとき4ナノ秒）以上がより好ましい。
【００９３】
〔分子動力学シミュレーション結果の出力〕
本発明の計算方法においては、ステップ１〜ステップ５の実行、またはステップＡ〜Ｄおよびステップ１〜５の実行、またはステップＡ〜Ｄおよびステップ１〜６を所定回数繰り返し実行した後に、解析のために必要な様々な物理量を出力できる。この物理量は所定回数の繰り返し時点における瞬間値、あるいは一定回数繰り返して得られた平均値のどちらか一方または両方を出力することができる。
【００９４】
出力する物理量は後述する様々な自由エネルギー計算に必要なポテンシャルエネルギー以外にも、原子の座標、圧力などが挙げられる。ポテンシャルエネルギーについては自由エネルギー計算のために必須であり、原子の座標については、分子動力学シミュレーション最大の特徴である原子の時間経過による運動をとらえることができるという性質から、ある繰り返し回数における瞬間値を出力することが好ましい。
【００９５】
〔自由エネルギー計算〕
本発明における自由エネルギー計算法は、分子動力学シミュレーションで用いられる一般的な方法であれば特に限定されない。しかし、スケールされた非結合性相互作用を用いたときのシミュレーションであるという観点から、非現実的な中間状態を定義する錬金術的自由エネルギー計算法を用いるとき、本発明の効果が発揮される。
【００９６】
錬金術的自由エネルギー計算法として熱力学的積分法、自由エネルギー摂動法、ベネット受容比法またはマルチステイトベネット受容比法などが挙げられる。これらの手法は高精度に自由エネルギーを計算できることが知られており、どの手法も好適に用いられる。
【００９７】
自由エネルギー摂動法は、各中間状態におけるポテンシャルエネルギーVtotal(λi) およびi番目の中間状態の分子配置におけるi+1番目の中間状態のスケール値λi+1を用いた場合の全ポテンシャルエネルギー<Vtotal(λi+1)>iと、下記数式（１６Ａ）および数式（１６Ｂ）を用いて自由エネルギーを求める手法である。
【００９８】
【数１６】
2018049608.tif 000018
【００９９】
（ここで、<…>はアンサンブル平均を、<…>iはi番目のレプリカの分子配置を用いた物理量であることを表す。）
ベネット受容比法は、各中間状態における全ポテンシャルエネルギーVtotal(λi)およびi番目の中間状態の分子配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi+1)>iに加え、i+1番目の中間状態の分子配置とi番目またはi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi)>i+1、Vtotal(λi+1)と、下記数式（１７Ａ）および数式（１７Ｂ）を用いて自由エネルギーを求める手法である。
【０１００】
【数１７】
2018049608.tif 000019
【０１０１】
（ここで、LF、LRはそれぞれ、i番目の中間状態の分子配置とi+1番目の中間状態のスケール値を用いてサンプリングした全ポテンシャルエネルギーと、i+1番目の中間状態の分子配置とi番目の中間状態のスケール値を用いてサンプリングした全ポテンシャルエネルギーのデータ数を表す。さらにβ=1/kBT、C= kBTln(LF/LR)を表す。）
マルチステイトベネット受容比法は、i番目の中間状態の分子配置においてk番目(k=1,…, i,…, nint)の中間状態のスケール値λkを用いた時の全ポテンシャルエネルギー<Vtotal(λk)>iと、下記数式（１８Ａ）および数式（１８Ｂ）を用いて自由エネルギーを求める手法である。
【０１０２】
【数１８】
2018049608.tif 000020
【０１０３】
（ここで、xnはNk回分子配置をサンプリングしたときのn番目の分子配置を表す。）
前述の通り、中間状態を定義する錬金術的自由エネルギー計算法を用いるとき、本発明の効果が最大限に発揮されるが、特定の計算条件では中間状態を定義しない自由エネルギー計算法でも大きな効果がある。
【０１０４】
例えば中間状態を定義せず、始状態と終状態のみ定義し、それらの相互作用エネルギーからなるエネルギー分布関数を求め近似汎関数を構築し、これより自由エネルギー計算するエネルギー表示法を、高分子分離膜中の分離物質のような分子配置のサンプリングが困難な系に適用する場合が挙げられる。高分子分離膜中の分離物質のような系では、高精度な自由エネルギーを得るため、分離物質の運動を促す目的でスケールされた非結合性相互作用を用いることがある。スケールされた非結合性相互作用を用いる場合でも、自由エネルギー計算に必要な分子配置のサンプリングは、始状態と終状態のみであり中間状態は関係ない。しかしNPTアンサンブルの計算を行った場合、中間状態の構造がシミュレーション的に不正確な構造となるため、始状態と終状態の構造にも影響を及ぼすこととなる。
【０１０５】
〔プログラムおよび記録媒体〕
本発明のプログラムは、コンピュータの構成要素に上記の自由エネルギーの計算方法の各ステップにおける計算手段として動作させるものである。このプログラムを用いれば、上述の計算方法を実行することができる。
【０１０６】
また、本発明の記録媒体は、当該プログラムをコンピュータ読み取り可能に記録したものである。この記録媒体のプログラムをコンピュータにて読み取って実行すれば、上述の計算方法を実行することができる。
【０１０７】
以上の手順を実行すれば、  分子集合体に溶質分子が溶解している系における、スケールされた非結合性相互作用を用いた分子動力学シミュレーションにおいて、カットオフ半径から無限遠のvdW相互作用から生じる圧力ＰvdWrc→∞を考慮した自由エネルギー計算を実行できないという課題を解決することができる。
【実施例】
【０１０８】
以下に実施例を挙げて本発明を具体的に説明するが、これらは例示的なものであって限定的なものではない。
【０１０９】
［実施例１］
ポリジメチルシロキサン（PDMS）鎖４本＋水１分子系に対し、分子動力学シミュレーションを実施した後、水のPDMSへの溶解自由エネルギー計算を行った。
【０１１０】
分子動力学シミュレーションはFORTRAN90プログラムを用いて行い、分子軌道計算は汎用プログラムGaussian（登録商標）09を用いて行った。また、計算にはクロック周波数2.93GHzのCore i7（登録商標）870を搭載したコンピュータ、およびクロック周波数2.7HzのXeon（登録商標） E5-2680を搭載したコンピュータを用いた。
【０１１１】
バイオビア社製の分子設計システムMaterials Studio（登録商標）8.0上で下記構造式（１）に示すポリマーを作成した。
【０１１２】
【化１】
2018049608.tif 000021
【０１１３】
構造式(1)において、n=40（分子量＝2996）とした。次いで、Materials Studio（登録商標）8.0 Amorphous Cell Moduleを用いて、表１に示した組成でPDMS分子と水分子を三次元周期境界条件下でランダムに配置し系を作成した。Materials Studio（登録商標）8.0の外部出力機能を使用し、原子の座標r 、結合情報、原子のポテンシャルパラメータタイプ、原子の質量ｍ、基本セルベクトルhが記載された分子構造ファイルを作成した。
【０１１４】
次に、出力した分子構造ファイルに記載されたポテンシャルパラメータタイプに対応したポテンシャルパラメータを記述したパラメータファイルを作成した。
【０１１５】
さらにテキストエディタを用いて中間状態数であるnint、相互作用計算のカットオフ半径ｒｃ、図２計算フローにおけるループの繰り返しを規定する回数nstep、出力する物理量とその出力間隔、設定温度T0、設定圧力P0、熱浴の慣性量Q、barostatの慣性量Wを記述した計算条件ファイルを作成した。
【０１１６】
【表１】
2018049608.tif 000022
【０１１７】
まず作成した分子構造ファイル、パラメータファイル、計算条件ファイルを分子動力学シミュレーションプログラムで読み取らせ、データをメモリに格納し、中間状態数nint分のコピーを作成した。
【０１１８】
計算条件ファイルにおいて設定したパラメータとして、中間状態数nint=21、カットオフ半径rc=10オングストローム、図２計算フローにおけるループの繰り返しを規定する回数nstep=16000000回、出力する物理量を原子の座標r、全ポテンシャルエネルギーとし、その出力間隔を下記ステップＡ〜Ｄ、ステップ１〜５のループを4000回繰り返す毎に１回とした。
【０１１９】
温度制御にNose−Hoover法〔M.Tuckerman,B.J.Berne,G.J.Martyna,J.Chem.Phys.,1992, 97,p.1990.〕、圧力制御にParrinello−Rahman法〔M. Parrinello、A. Rahman、J. Appl. Phys.、1981、52、p.7182〕を用いて、原子数・温度・圧力一定（NPT）アンサンブルを構成することにより、系を設定温度T0、設定圧力P0に制御した状態で下記ステップＡ〜Ｄ、ステップ１〜５を実施した。なお、設定温度T0=20℃、設定圧力P0=1atm、熱浴の慣性量Q=120.95(kJ・ps2)/mol、barostatの慣性量W=483.79(kJ・ps2)/molとした。
【０１２０】
ＮＰＴアンサンブルにおける分子動力学シミュレーションを実施するためには、下記ステップ５において原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新することが必要である。そのため分子構造ファイルおよび計算条件ファイルで入力していない、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvに初期値を与える必要がある。そこで原子の速度vの初期値としてT0=20℃におけるマクスウェル分布に従うよう計算された速度を、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvの初期値として全ての成分に0を与えメモリに格納した。
【０１２１】
各中間状態における非結合性相互作用のスケール値λiをコンピュータ上で計算し、各中間状態のメモリ領域に格納した。λiは１をnintで等分した値を用い、nintが21であるため(λ1, λ2, …, λ20, λ21)=(0.00, 0.05, …, 0.95, 1.00)となる。またλsol-mol＝λi、すなわち溶質（水）‐分子集合体（PDMS）間非結合性相互作用にのみ適用し、PDMSの分子内非結合性相互作用と、PDMS‐PDMS間非結合性相互作用のそれぞれのスケール値はλmol-intra＝λmol-mol＝１とし、スケールしなかった。本実施例における溶質は水１分子であるため、溶質の分子内非結合性相互作用のスケール値λsol-intra および溶質‐溶質間非結合性相互作用のスケール値λsol-solは設定しない。
【０１２２】
各中間状態において同時に次のステップＡ〜Ｄおよびステップ１〜５を実行した。
【０１２３】
〔ステップＡ〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目およびi+1番目の中間状態における非結合性相互作用のスケール値λiおよびλi+1、ならびに数式（１０）で示したDreiding型のvdWポテンシャルを含む、数式（１）で示したポテンシャル関数Vtotalを用いて、VvdWを除くポテンシャルエネルギー{Vtotal (λi)‐VvdW(r,λi)}およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの{<Vtotal (λi+1)>i‐<VvdW(r,λi+1)>i}を計算しメモリに格納した。また、数式（１０）中のδ＝1.0×10-20ｍ２とし、数式（１）第６項のクーロンポテンシャルおよびその力の計算にはEwald法を用い、Ewald法のパラメータとして逆空間のクーロンポテンシャルにα＝0.21オングストローム-1、｜ｎ｜２ｍａｘ＝50を用いて計算した。
【０１２４】
ポテンシャルパラメータの設定は以下の通りである。まず、PDMS分子については、分子内ポテンシャルパラメータとして、結合長・結合角の平衡位置R0、θ0にはHF/6-31G(d,p)レベルの分子軌道計算によって決定した値を、力の定数Kbond、Kangle、Ktorsion、Kinversionおよび平衡位置φ0、ω0にはAMBERパラメータ、hybrid/UAパラメータおよび桑島ポテンシャルパラメータを用いた。vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０にはAMBERパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。
【０１２５】
水分子については、結合長・結合角の平衡位置R0、θ0にはTIP3Pパラメータ[W. L. Jorgensen 、R. W. Impey 、J. Chandrasekhar 、 J. D. Madura、M. L. Klein 、J. Chem. Phys.、1983、79、p.926.〕を用い、結合長・結合角の力の定数Kbond、KangleはＩＲの実測値を再現するように決定した値を用い、vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０と電荷にはTIP3Pパラメータを用いた。
【０１２６】
〔ステップＢ〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目およびi+1番目の中間状態における非結合性相互作用のスケール値λiおよびλi+1、ならびに相互作用計算のカットオフ半径ｒｃを用いて、前記数式（１０）で示したDreiding型のvdWポテンシャルを用いて、原子間距離が0からカットオフ半径rc以内のvdW相互作用のポテンシャルエネルギーVvdW(0→rc)(r,λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときのポテンシャルエネルギー<VvdW(0→rc)(r,λi+1)>iを計算しメモリに格納した。
【０１２７】
〔ステップＣ〕
メモリに格納された結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目およびi+1番目の中間状態における非結合性相互作用のスケール値λiおよびλi+1、相互作用計算のカットオフ半径rc、前記数式（１２Ａ）、（１２Ｂ）を用いて、カットオフ半径rcから無限遠までのvdW相互作用のポテンシャルエネルギーVvdW(rc→∞)(r,λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときのポテンシャルエネルギー<VvdW(rc→∞)(λi+1)>iを計算しメモリに格納した。
【０１２８】
〔ステップＤ〕
ステップＡ〜Ｃでメモリに格納した各成分のポテンシャルエネルギーと前記数式（１）より全ポテンシャルエネルギーVtotal(λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi+1)>iを計算しメモリに格納した。
【０１２９】
〔ステップ１〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目の中間状態における非結合性相互作用のスケール値λi、ならびに、前記数式（１）における全ポテンシャル関数Vtotal(λi)を原子の座標rで微分した式を用いて分子内力・分子外力を計算しメモリに格納した。
【０１３０】
〔ステップ２〕
ステップ１でメモリに格納した分子間力の一つであるvdW力（数式（１０）の微分）と、原子の座標rからビリアル定理（数式（１３Ａ）および（１３Ｂ））を用いて圧力ＰvdW0→ｒｃを計算しメモリに格納した。
【０１３１】
〔ステップ３〕
結合情報、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、i番目の中間状態における非結合性相互作用のスケール値λiおよび相互作用計算のカットオフ半径ｒｃと数式（４Ａ）、数式（１４Ａ）および数式（１４Ｂ）から圧力ＰvdWrc→∞を計算しメモリに格納した。
【０１３２】
〔ステップ４〕
ステップ２および３でメモリに格納した、カットオフ半径内のvdW相互作用から生じる圧力ＰvdW0→ｒｃと、カットオフ半径から無限遠のvdW相互作用から生じる圧力ＰvdWrc→∞と、前記数式（５）により圧力ＰvdW0→∞を計算しメモリに格納した。
【０１３３】
〔ステップ５〕
メモリに格納された原子の質量ｍ、原子の座標r、原子の速度v、熱浴の慣性量Q、熱浴の位置ζ、熱浴の速度ζv、barostatの慣性量W、基本セルベクトルh、基本セルベクトルの速度hv、分子内力、分子間力を用いて、微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを数値積分によって計算し、更新した。数値積分には多重時間刻み数値積分法（RESPA法）を用い、時間刻みΔtは表２に示した値を用いた。
【０１３４】
【表２】
2018049608.tif 000023
【０１３５】
ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを初期値としてステップＡ〜Ｄ、ステップ１〜５を再度実行させることをnstep＝16000000回繰り返した。RESPA法における最小の時間刻みが0.25フェムト秒であるため、この繰り返し回数は4ナノ秒に相当する。またステップ７を実施する間隔である8000回は２ピコ秒に相当する。
【０１３６】
なお本実施例ではNose−Hoover法により温度、Parrinello−Rahman法により圧力を制御しているためスッテプ６において熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新した。
【０１３７】
ステップＡ〜Ｄ、ステップ１〜５を4000回（１ピコ秒に相当）繰り返すごとに、原子の座標r、各レプリカにおけるポテンシャルエネルギーVtotal(λi) およびi番目のレプリカの分子配置におけるi+1番目のレプリカのスケール値λi+1を用いた場合の全ポテンシャルエネルギー<Vtotal(λi+1)>iを出力した。
【０１３８】
分子動力学シミュレーション終了後、Vtotal(λi)および<Vtotal(λi+1)>iの時系列データと前記数式（１６Ａ）、（１６Ｂ）に示した自由エネルギー摂動法を用いて、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを計算した。
【０１３９】
なお、本実施例の計算では構造緩和のため、上述の計算の前に20℃の原子数・温度・体積一定（ＮＶＴ）アンサンブルの分子動力学シミュレーションを10ピコ秒、1atm／20℃のＮＰＴアンサンブルの分子動力学シミュレーションを90ピコ秒、1atmの分子動力学シミュレーションにおいて温度を20℃〜500℃まで200ピコ秒周期で上昇下降を繰り返し緩和させるシミュレーテッドアニーリングを14サイクル（2800ピコ秒に相当）、1atm／20℃のＮＰＴアンサンブルの分子動力学シミュレーションを300ピコ秒実施した。さらにステップ１を始める前に、構造緩和計算後の構造を用いて各中間状態に対応するλiを用いて分子動力学シミュレーションを300ピコ秒実施し、中間状態での構造緩和を行った。
【０１４０】
［実施例２］
次に示す条件以外は実施例１と同じ分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを計算した。
【０１４１】
計算条件ファイルにハミルトニアンレプリカ交換を行う中間状態数nint＝21（＝レプリカ数）、交換を実施する間隔lrep＝8000回に１回と記述しメモリに格納した。実施例１における中間状態をレプリカと読み替え、各レプリカにおいて独立に実施例１のステップＡ〜Ｄおよびステップ１〜５を実施した。
【０１４２】
すべてのレプリカにおいてステップ５までの計算が終わった時点で、これまでステップ１〜５まで8000回繰り返してきたならば、1番目のレプリカからnint番目のレプリカまで順に次のステップ６を実施する。それ以外の場合はステップ６を実施しない。
【０１４３】
〔ステップ６〕
i番目のレプリカとスケール値λi+1を持つｉ＋１番目のレプリカにおける、ポテンシャルエネルギーVtotal(λi or i+1)と前記数式（１５）に示したメトロポリスの方法による判定基準を用い、i番目とi+1番目のレプリカのスケール値を交換する確率αを計算する。Fortran90の組み込みサブルーチンにより0〜1までの乱数を発生させ、交換確率以下の値だったらλiとλi+1との交換を行い、新たなλをメモリに格納する。
【０１４４】
ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hv、ステップ６を実施した場合は更新したスケール値λiを初期値としてステップ１〜ステップ６を再度実行させることをnstep＝16000000回繰り返した。RESPA法における最小の時間刻みが0.25フェムト秒であるため、この繰り返し回数は4ナノ秒に相当する。またステップ６を実施する間隔である8000回は２ピコ秒に相当する。
【０１４５】
［実施例３］
95wt％メタノール水溶液系＋水１分子系に対し、分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsolの計算を行った。なお計算条件は、下記の条件以外は実施例１と同様である。
【０１４６】
Materials Studio（登録商標）8.0 Amorphous Cell Moduleを用いて、表３に示した組成でメタノール分子と水分子を三次元周期境界条件下でランダムに配置し系を作成した。
【０１４７】
【表３】
2018049608.tif 000024
【０１４８】
計算条件ファイルにおいて設定したパラメータとして、ループの繰り返しを規定する回数nstep=4000000回とし、この繰り返し数は1ナノ秒に相当する。
【０１４９】
メタノール分子のポテンシャルパラメータの設定は以下の通りである。分子内ポテンシャルパラメータとして、結合長・結合角の平衡位置R0、θ0にはHF/6-31G(d,p)レベルの分子軌道計算によって決定した値を、力の定数Kbond、Kangle、Ktorsion、Kinversionおよび平衡位置φ0、ω0にはAMBERパラメータおよび桑島ポテンシャルパラメータを用いた。vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD０、R０にはAMBERパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。
【０１５０】
［実施例４］
次に示す条件以外は実施例２と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを計算した。
【０１５１】
ステップＡ〜Ｄにおいて全ポテンシャルエネルギーVtotal(λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi+1)>iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi-1)>iを計算し、メモリに格納した。
【０１５２】
分子動力学シミュレーション終了後、Vtotal(λi)、<Vtotal(λi+1)>i、<Vtotal(r,λi-1)>iの時系列データと数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて、溶解自由エネルギーΔGsolvpolyを計算した。
【０１５３】
［実施例５］
次に示す条件以外は実施例３と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsolを計算した。
【０１５４】
スステップＡ〜Ｄにおいて全ポテンシャルエネルギーVtotal(λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi+1)>iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi-1)>iを計算し、メモリに格納した。
【０１５５】
分子動力学シミュレーション終了後、Vtotal(λi)、<Vtotal(λi+1)>i、<Vtotal(r,λi-1)>iの時系列データと数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて溶解自由エネルギーΔGsolvsolを計算した。
【０１５６】
［実施例６］
次に示す条件以外は実施例２と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを計算した。
【０１５７】
ステップＡ〜Ｄにおいて、i番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λkを用いたときの全ポテンシャルエネルギー<Vtotal(λk)>iを計算し、メモリに格納した。
【０１５８】
分子動力学シミュレーション終了後、V<Vtotal(λk)>iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔGsolvpolyを計算した。
【０１５９】
［実施例７］
次に示す条件以外は実施例３と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsolを計算した。
【０１６０】
ステップＡ〜Ｄにおいて、i番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λkを用いたときの全ポテンシャルエネルギー<Vtotal(λk)>iを計算し、メモリに格納した。
【０１６１】
分子動力学シミュレーション終了後、V<Vtotal(λk)>iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔGsolvsolを計算した。
【０１６２】
［比較例１］
次に示す条件以外は実施例２と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを自由エネルギー摂動法により計算した。
【０１６３】
ステップＡ〜Ｂを実施した後、ステップＣを実施せずにステップＤを実施した。またステップ１およびステップ２を実施した後、ステップ３を実施せず、ステップ４およびステップ５を実施した。すべてのレプリカにおいてステップ５までの計算が終わった時点で、これまでステップ１〜５まで8000回繰り返してきたならば、1番目のレプリカからnint番目のレプリカまで順に次のステップ６を実施した。以後ステップＡ〜Ｂ、ステップＤ、ステップ１〜２、ステップ４〜５を16000000回繰り返し実施し、繰り返しを実施した数が8000で割り切れるときはステップ５の後にステップ６を実施した。
【０１６４】
ステップＣおよびステップ３を実施しないということは、カットオフ半径から無限遠までのvdW相互作用から生じるポテンシャルエネルギーと圧力を無視することに相当する。
【０１６５】
［比較例２］
次に示す条件以外は実施例３と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsolを自由エネルギー摂動法により計算した。
【０１６６】
ステップＡ〜Ｂを実施した後、ステップＣを実施せずにステップＤを実施した。またステップ１およびステップ２を実施した後、ステップ３を実施せず、ステップ４およびステップ５を実施し、すべてのレプリカにおいてステップ５までの計算が終わった時点で、これまでステップ１〜５まで8000回繰り返してきたならば、1番目のレプリカからnint番目のレプリカまで順に次のステップ６を実施した。以後ステップＡ〜Ｂ、ステップＤ、ステップ１〜２、ステップ４〜５を4000000回繰り返し実施し、繰り返しを実施した数が8000で割り切れるときはステップ５の後にステップ６を実施した。
【０１６７】
［比較例３］
次に示す条件以外は比較例１と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを計算した。
【０１６８】
ステップＡ〜Ｂ、ステップＤにおいて全ポテンシャルエネルギーVtotal(λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi+1)>iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi-1)>iを計算し、メモリに格納した。
【０１６９】
分子動力学シミュレーション終了後、Vtotal(λi)、<Vtotal(λi+1)>i、<Vtotal(λi-1)>iの時系列データと前記数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて、溶解自由エネルギーΔGsolvpolyを計算した。
【０１７０】
［比較例４］
次に示す条件以外は比較例２と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsolを計算した。
【０１７１】
ステップＡ〜Ｂ、ステップＤにおいて全ポテンシャルエネルギーVtotal(λi)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi+1)>iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<Vtotal(λi-1)>iを計算し、メモリに格納した。
【０１７２】
分子動力学シミュレーション終了後、Vtotal(λi)、<Vtotal(λi+1)>i、<Vtotal(λi-1)>iの時系列データと前記数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて、溶解自由エネルギーΔGsolvsolを計算した。
【０１７３】
［比較例５］
次に示す条件以外は比較例１と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyを計算した。
【０１７４】
ステップＡ〜Ｂ、ステップＤにおいてi番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λkを用いたときの全ポテンシャルエネルギー<Vtotal(λk)>iを計算し、メモリに格納した。
【０１７５】
分子動力学シミュレーション終了後、V<Vtotal(λk)>iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔGsolvpolyを計算した。
【０１７６】
［比較例６］
次に示す条件以外は比較例２と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsolを計算した。
【０１７７】
ステップＡ〜Ｂ、ステップＤにおいてi番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λkを用いたときの全ポテンシャルエネルギー<Vtotal(λk)>iを計算し、メモリに格納した。
【０１７８】
分子動力学シミュレーション終了後、V<Vtotal(λk)>iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔGsolvsolを計算した。  実施例１〜７および比較例１〜６で得られた自由エネルギーを表４にまとめた。なお表４中の自由エネルギー計算法では、自由エネルギー摂動法をFEP、ベネット受容比法をBAR、マルチステイトベネット受容比法をMBARと表した。
【０１７９】
ここで、PDMSに水が１分子溶解するときの溶解自由エネルギーΔGsolvpolyと95wt％メタノール水溶液への水分子の溶解自由エネルギーΔGsolvsol、下記数式（１９）より95wt％メタノール水溶液からPDMS膜に水が１分子溶解したときの自由エネルギーΔGsorpを求めることができる。
【０１８０】
【数１９】
2018049608.tif 000025
【０１８１】
このΔGsorpは高分子膜を水溶液に含浸させるという簡便な実験より求めることができるため、シミュレーション結果と直接比較でき、シミュレーションの精度確認の指標に適している。非特許文献〔M. H. V. Mulder and C. A. Smolders, Separation Science and technology, 26, 85(1991). 〕によると、20℃の95wt％メタノール水溶液にPDMS膜を含浸させたとき、吸水率は0.13wt％であった。非特許文献〔T. Kawakami, I. Shigemoto, N. Matubayasi, J. Chem. Phys. 137, 234903 (2012).〕より吸水率と下記数式（２０）を用いることによりΔGsorpの実験値を求めることができる。
【０１８２】
【数２０】
2018049608.tif 000026
【０１８３】
（ここでWは吸水率、ρは密度、kBはボルツマン定数、Tは絶対温度をそれぞれ表す。）
数式（２０）よりΔGsorpの実験値は3.76kcal/molと求められる。なおポリマーの密度は0.971g/cm3、９５ｗｔ％メタノール水溶液の密度は、水の密度1.0g/cm3およびメタノールの密度0.792g/cm3と混合比率から求めた0.802g/cm3を用いた。計算値および実験値のΔGsorpも併せて表４に示した。なお実施例１〜３におけるΔGsorpの（）外の値は、実施例２と実施例３から求めた値で、（）内の値は実施例１と実施例３から求めた値である。すなわち（）内の値はPDMS・水系においてハミルトンレプリカ交換を行わなかった場合のΔGsorpを表す。
【０１８４】
【表４】
2018049608.tif 000027
【０１８５】
実施例１および実施例３より求めたΔGsorpと、実施例２および実施例３より求めたΔGsorpを比較すると、実施例２および実施例３より求めた方が実験値に近いことがわかる。これは、実施例１においてはハミルトニアンレプリカ交換を行わなかったことにより、スケール値λiが大きいときにPDMS中の水の運動性が悪くなり、局所的な分子配置しかサンプリングできなかった影響であると考える。
【０１８６】
また、ベネット受容比法を用いた実施例４および実施例５より求めたΔGsorpとマルチステイトベネット受容比法を用いた実施例６および実施例７より求めたΔGsorpは、実施例２および実施例３のエネルギー摂動法より求めた値よりも実験値に近く、非常に精度良く溶解自由エネルギーを計算できているといえる。一般的にエネルギー摂動法よりもベネット受容比法およびマルチステイトベネット受容比法の方が高い計算精度を持つといわれており、今回の結果もそれを支持する。
【０１８７】
カットオフ半径から無限遠のvdW相互作用を無視した比較例１〜６までの結果を用いた場合のΔGsorpと実施例の結果を比較する。比較例１および比較例２または比較例３および比較例４または比較例５および比較例６から求めたΔGsorpは、実験値から約0.6kcal/mol以上かい離し、実施例の結果よりも計算精度が低い。したがって本発明である自由エネルギー計算手法は、分子集合体に溶質分子が溶解している系において非常に高精度な自由エネルギーを与えるといえる。
【図１】
2018049608.tif 000037
【図２】
2018049608.tif 000038
【図３】
2018049608.tif 000039