JP2018049608A - 自由エネルギーの計算方法 - Google Patents

Abstract

【課題】分子集合体に溶質分子が溶解している系における、スケールされた非結合性相互作用を用いた分子動力学シミュレーションにおいて、カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を考慮した、高精度な自由エネルギー計算方法を提供する。【解決手段】コンピュータに、原子の座標、スケール値など初期条件を入力し、（１）分子内力・分子間力を計算し、（２）カットオフ半径内のファンデルワールス相互作用から生じる圧力ＰvdW0→rcを計算し、（３）カットオフ半径から無限遠のファンデルワールス相互作用から生じる圧力ＰvdWrc→∞を計算し、（４）ＰvdW0→rcとＰvdWrc→∞との和からＰvdW0→∞計算し、（５）微小時間Δt後の原子の座標、速度を数値積分によって更新し、（５）の更新値を初期値として（１）〜（５）を再度実行させ、これを規定回数繰り返すことを含む。【選択図】図１

Description

本発明は、コンピュータを利用した分子動力学シミュレーションにおける、分子集合体に溶質分子が溶解している系の自由エネルギーの計算方法に関する。

分子動力学シミュレーションは、分子や分子集合体の高次構造や物性を詳細に解析するための手法として広く用いられる。

分子動力学シミュレーションより得られる、分子集合体に関する重要な物性の一つとして、分子集合体へ溶質が溶解するときの自由エネルギーが挙げられる。

分子集合体が低分子液体である場合、溶質が低分子液体へ溶解するかの議論をする上で、自由エネルギーは非常に重要な指標となる。また低分子液体として水と１−オクタノールを選択し、これらへの溶質の溶解自由エネルギーを計算すれば、溶質の水／１−オクタノール分配係数をシミュレーションより求めることができる（非特許文献１）。水／１−オクタノール分配係数は、薬物動態や環境動態の指標として用いられ産業応用上も非常に重要な物理量である。溶質として実験が困難な分子や、新規候補分子を用い、水／１−オクタノール分配係数を予測したいという需要が大きいため、シミュレーションによる予測も活発に行われている（非特許文献２）。

また分子集合体が高分子分離膜、溶質分子が分離物質の場合、分離物質の高分子膜への溶解自由エネルギーを表し、分離物質の膜透過性を議論するに当たり重要な物理量となる。この自由エネルギーから溶解度係数を求めることにより、高分子分離膜における物質透過性をシミュレーションより予測することができる（非特許文献３）。なお溶解度係数には、自由エネルギーの誤差が自然対数の底のべき乗で影響を及ぼすため、高分子分離膜系で膜透過性を議論する場合は非常に高精度な溶解自由エネルギー計算が必要とされる。

分子動力学シミュレーションを用いた代表的な自由エネルギー計算方法は、熱力学的積分法、自由エネルギー摂動法など、分子間相互作用をスケールした中間状態を考える錬金術的自由エネルギー計算法が挙げられる。（非特許文献４、５）これらの手法を溶質の分子集合体への溶解の問題に適用する場合、溶質‐分子集合体に属する分子の分子間非結合相互作用を、通常の相互作用から完全に相互作用が消滅した状態まで、徐々に非結合相互作用をスケールした中間状態を計算していかなければならない。

ここで分子動力学シミュレーションにおいて得られる、自由エネルギーを含む物理量は、理論上は取りうる分子配置を全てサンプリングした場合に真の値が得られる。そのため物理量の精度は分子の配置のサンプリング数に依存する。

前記の通り、高分子分離膜の物質透過性予測には非常に高精度な溶解自由エネルギー計算が必要とされる。しかし高分子の遅い分子運動のため、高分子膜中の分離物質は通常の分子動力学シミュレーションの計算時間ではほとんど移動しない。そのため、高分子は、通常の非結合相互作用に近い大きさの相互作用を持つ中間状態では、分離物質の運動性が低く分子配置を効率的にサンプリングできない。そこで、非結合相互作用のスケール値を時間と共に増減させる方法（特許文献１）や、レプリカと呼ばれる全く同じ系を発生させ、レプリカ毎に異なる非結合相互作用のスケール値を設定し、レプリカ毎に独立に分子動力学シミュレーションを行い、一定時間ごとにある確率でスケール値のみ交換するというハミルトニアンレプリカ交換分子動力学法（非特許文献６）といったスケール値を変化させるサンプリング方法を組み合わせる必要がある。

他にもタンパク質といった生体高分子系でも盛んに自由エネルギー計算が行われている。例えば分子集合体をタンパク質と水、溶質を薬剤分子とすると、薬剤のタンパク質レセプターへのドッキング時の自由エネルギーを求めることができ、薬剤設計において重要な指標となる。

特開２０１１−１２３８７４号公報

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における全ポテンシャルエネルギーV^total(λ_i)は、下記数式（１）のように表される。

（ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
ファンデルワールス（vdW）相互作用エネルギーV^vdW(0→∞)は、カットオフ半径内のvdW相互作用エネルギーV^vdW(0→rc)とカットオフ半径から無限遠のvdW相互作用エネルギーV^{vdW(rc→∞)}の和で表される。

ここで中間状態の分子配置を求めるため、スケールした非結合性相互作用を使用した場合、カットオフ半径から無限遠のvdW相互作用エネルギーV^{vdW(rc→∞)}を計算できず、中間状態λ_iにおける全ポテンシャルエネルギーV^total(λ_i)も正確に求めることができなかった。V^{vdW(rc→∞)}の寄与は小さいため、考慮されない場合も多いが、高分子分離膜系の溶解度係数のように非常に高精度な自由エネルギー計算が求められる場合は必要である。

カットオフ半径から無限遠のvdW相互作用を無視した場合、本当に問題になるのはポテンシャルエネルギーではなく、圧力を正しく求められなくなることである。ここで、分子動力学シミュレーションにおいて、系の全圧力P_totalは、下記数式（２）のように表される。

（ただし、P_intraは分子内相互作用から生じる圧力、P_vdW ^0→∞はvdW相互作用から生じる圧力、P_coulomb ^0→∞はクーロン相互作用から生じる圧力、P_Kは原子の運動エネルギーより生じる圧力である。）
vdW相互作用から生じる圧力Ｐ_vdW ^0→∞は、カットオフ半径内のvdW相互作用から生じる圧力Ｐ_vdW ^0→rcとカットオフ半径から無限遠のvdW相互作用から生じる圧力Ｐ_vdW ^rc→∞の和で表される。

中間状態の分子配置を分子動力学シミュレーションより求めるため、スケールした非結合性相互作用を使用した場合、ポテンシャルエネルギーと同様に、正確な圧力を計算することができない。

正確な圧力が計算できない場合、分子動力学シミュレーションにおいて標準状態を再現するために行う、圧力を一定に制御した状態のシミュレーションを正しく行えず、シミュレーションセルの大きさを過大評価あるいは過小評価してしまう。結果として、分子間の距離が通常よりも近くまたは遠くなるため、正確な相互作用エネルギーが得られず、自由エネルギーを過大または過小評価することになる。このようにスケールされたvdW相互作用自体のカットオフ半径から無限遠の寄与よりも、圧力を正確に計算できなくなり、密度が変化し、正確でない分子間距離で相互作用を計算してしまうことに問題があった。

特に、高分子膜中の分離物質系のように高分子と低分子の混合系の場合、圧力を正確に計算できないという影響が大きく、自由エネルギーを高精度に計算できないという問題があった。

以上述べたように、分子集合体への溶質分子の溶解自由エネルギーを求めるためには、スケールされた非結合性相互作用を用いた分子動力学シミュレーションを実施する必要がある。しかし、スケールされた非結合相互作用を用いたときのカットオフ半径から無限遠のvdW相互作用から生じるポテンシャルエネルギーV^{vdW(rc→∞)}および圧力Ｐ_vdW ^rc→∞を求める方法が未だ提唱されていないため、高精度な自由エネルギー計算が行えなかった。そこで本発明では、上記問題を解決した、高精度な自由エネルギー計算方法を提供することを目的とする。

本発明は上記課題を解決するために、次の手段を採用するものである。すなわち、本発明は以下のとおりである。

コンピュータを利用した、分子集合体に溶質分子が溶解している系の分子動力学シミュレーションにおいて、下記数式（３）のように、
（i）溶質分子の分子内vdW相互作用V^vdw _sol-intra(r, λ_sol-intra)と、
（ii）分子集合体に属する分子の分子内vdW相互作用V^vdw _mol-intra(r, λ_mol-intra)と、
（iii）溶質分子‐溶質分子間vdW相互作用V^vdw _sol-sol(r, λ_sol-sol)と、
（iv）分子集合体に属する分子の分子間vdW相互作用V^vdw _mol-mol(r, λ_mol-mol)と、
（v）溶質分子‐分子集合体に属する分子の分子間vdW相互作用V^vdw _sol-mol(r, λ_sol-mol)
との和（ただし0≦λ≦1）で表される、非結合性相互作用のスケール値λを含むスケールvdW相互作用ポテンシャルV^vdW(r, λ)を用いた場合の自由エネルギー計算法であって、

コンピュータに、
・前記溶質分子と分子集合体に含まれる原子の座標r
および、必要に応じ
・原子の速度v
の初期値、
ならびに
・結合情報
・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
・原子の質量m
・相互作用計算のカットオフ半径r_c
・中間状態を定義する非結合性相互作用のスケール値λ_i（λ₁、…λ_i…、λ_nint、0〜1を中間状態の数nintで分割した値）
を入力した後、
非結合性相互作用のスケール値λ_iをもつ中間状態iそれぞれにおいて、コンピュータ上で少なくとも、
（ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、
前記数式（１）で示した全ポテンシャルエネルギーV^total(λ_i) を原子の座標rで微分した式を用いて、前記分子集合体および前記溶質の分子内力(数式（１）第１項から第４項の微分の和)および分子間力（数式（１）第５項と第６項の微分の和）を計算するステップ；
（ステップ２）ステップ１で計算した分子間力の一つであるvdW力（数式（１）第５項の微分）と、
原子の座標rからビリアル定理を用いて圧力Ｐ_vdW ^0→ｒｃを計算するステップ；
（ステップ３）原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λ_iおよび相互作用計算のカットオフ半径ｒ_ｃから、
下記数式(４Ａ)および（４Ｂ）を用いて圧力Ｐ_vdW ^rc→∞を計算するステップ；

（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎ_α、Ｎ_βはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
（ステップ４）下記数式（５）によりvdW相互作用から生じる圧力Ｐ_vdW ^0→∞を計算するステップ；

（ステップ５）原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
・原子の座標r
・原子の速度ｖ
を数値積分によって更新するステップ：
を実行させ、
さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含む、自由エネルギーの計算方法。

本発明における自由エネルギーの計算方法を示すフローチャートである。 図１の計算に、カットオフ半径から無限遠のvdW相互作用のポテンシャルエネルギーへの寄与を考慮するステップＡ〜Ｄを加えた自由エネルギーの計算方法を示すフローチャートである。 図２の計算に、ハミルトニアンレプリカ交換を行うステップ６を加えた自由エネルギーの計算方法を示すフローチャートである。

図１、図２、図３に本発明の計算方法のフローを示す。以下、本発明の好ましい実施形態を、図１、図２、図３のフローに従って説明する。なお図１はカットオフ半径から無限遠のvdW相互作用の圧力への寄与を考慮する自由エネルギー計算であり、図２は図１の計算に加えカットオフ半径から無限遠のvdW相互作用のポテンシャルエネルギーへの寄与を考慮する自由エネルギー計算、図３は図２の計算に加え、ハミルトニアンレプリカ交換法を組み合わせた自由エネルギー計算法である。

本発明は分子動力学シミュレーションにおける、分子集合体に溶質分子が溶解している系の自由エネルギーの計算方法である。ここで分子集合体とは、有機溶媒や水といった溶液、合成高分子、生体高分子またはそれらの混合系など分子が凝集している系ならばどのようなものでも良い。凝集しているとは、気体状態でないことを示し、より具体的には、分子動力学シミュレーションにおいて、密度が0.3g/cm³以上で分子が基本セル内に存在する状態を意味するものとする。

合成高分子とは、ポリエチレン、ポリエチレンテレフタラートなど人工的にモノマーを合成して得られる高分子を指し、重合度が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を用いて、前記数式（１）における全ポテンシャルエネルギーV^total(λ_i)の中で、V^vdW(r,λ_i)以外のポテンシャルエネルギー{V^total (λ_i)‐V^vdW(r,λ_i)}を計算し、メモリに格納する。

ここで、数式（１）の右辺第１項V^bond(R)は結合長、第２項V^angle(θ)は結合角、第３項V^torsion(φ)は二面角、第４項V^inversion(ω)は反転、第５項V^vdw(r,λ)はｖｄＷ、第６項V^coulomb(r,λ)はクーロンのそれぞれのポテンシャルエネルギーを定義するポテンシャル関数を表わす。なお、数式（１）の右辺第１項から第４項を総称して分子内ポテンシャル関数と呼び、右辺第５項、第６項をまとめて分子間ポテンシャル関数と呼ぶ。また、分子内ポテンシャル関数、分子間ポテンシャル関数を原子の位置で微分した式を用いて計算した力をそれぞれ分子内力、分子間力と呼ぶ。

数式（１）における各ポテンシャル関数の具体的な表式は下記数式（６）で表される。

数式（６）の右辺第５項、第６項はスケールされた非結合性相互作用を表し、非結合性相互作用のスケール値λ_i=1のときスケールされていない一般的な分子動力学シミュレーションで使用するポテンシャル関数に等しい。

前記数式（６）中の結合長、結合角、二面角、反転おけるそれぞれの力の定数K^bond、K^angle、K^torsion、K^inversionおよびそれぞれの平衡位置R⁰、θ⁰、φ⁰、ω⁰には、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.〕が好適に用いられる。

前記数式（１）および数式（６）の第６項V^coulomb(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）溶質分子の分子内非結合性相互作用V_sol-intra(r, λ_sol-intra)と、
（ii）分子集合体に属する分子の分子内非結合性相互作用V_mol-intra(r, λ_mol-intra)と、
（iii）溶質分子‐溶質分子間非結合性相互作用V_sol-sol(r, λ_sol-sol)と、
（iv）分子集合体に属する分子の分子間非結合性相互作用V_mol-mol(r, λ_mol-mol)と、
（v）溶質分子‐分子集合体に属する分子の分子間非結合性相互作用V_sol-mol(r, λ_sol-mol)
との和で表される。

非結合性相互作用のスケール値λの範囲として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とすることによって、分子内水素結合を持つ低分子の配座を維持しながら低分子の運動性を上昇させることができる。エチレングリコールのような分子内水素結合を持つ低分子は、分子内非結合性をスケールすると水素結合が弱まるため本来とるべき配座を再現しない。このように分子内非結合性相互作用が低分子のコンフォメーションに大きく寄与する場合、低分子内の非結合性相互作用をスケールしなければ低分子内の配座を維持しながら運動性を上昇させることができる。

さらにカットオフ半径r_cは、計算時間と計算精度の観点から8Å≦r_c≦12Åが好ましく、r_c＝10Åが特に好ましい。

〔ステップＢ〕
スケールされたvdW相互作用である前記数式（６）第５項V^vdw(r,λ_i)は、下記数式（８）のようにカットオフ半径以内のvdW相互作用エネルギーV^vdw(0→rc)(r,λ_i)とカットオフ半径から無限遠のvdW相互作用エネルギーV^{vdw(rc→∞)}(r,λ_i)の和で表される。

数式（８）で示したカットオフ半径以内のvdW相互作用ポテンシャルV^vdW(0→rc)を計算するステップＢにおいては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、相互作用計算のカットオフ半径ｒ_ｃを用いてV^vdW(0→rc)を計算し、メモリに格納する。

カットオフ半径以内のvdW相互作用エネルギーV^vdw(0→rc)(r,λ_i)の具体的な関数形としては、下記数式（９）に示したレナードジョーンズポテンシャルやモースポテンシャル、バッキンガムポテンシャルが好ましく用いられる。

数式（９）において、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型でも好適に用いることができる。

前記数式（９）中のvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０には、前述のDREIDING、AMBER、hybrid/UA、桑島らのパラメータ等の公知のパラメータを好適に用いることができる。計算精度について考慮した場合、D^０、R^０については、密度と蒸発熱の実験値を再現するように決定することがより好ましい。

前記数式（９）、数式（１０）におけるエネルギーと力の距離依存性を調節するパラメータδについては、δ＝0ｍ²とし、かつ非結合性相互作用のスケール値λが非常に小さい値をとる場合に、前記数式（１）および数式（６）の第５項V^vdw(r,λ)で表わされるvdW相互作用が急激に変化し、数値積分の安定性が低下する。この傾向は原子間の距離が近距離となった場合に特に顕著になる。この問題を解決するために、非結合性相互作用のスケール値λの最小値を非常に小さな値に設定する場合には、δ＞0ｍ²とすることが好ましい。さらに好ましくは、δ≧5.0×10^-21ｍ²である。一方、δの値が大きすぎると、スケール値が変化した際、vdW力の変化が大きくなりすぎるため、δ≦5.0×10^-19ｍ²とすることが好ましい。このδの値としては、δ＝1.0×10^-20ｍ²が好適に用いられる。

〔ステップＣ〕
前記数式（８）で示したカットオフ半径から無限遠のvdW相互作用ポテンシャルV^{vdW(rc→∞)}を計算するステップＣにおいては、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、相互作用計算のカットオフ半径ｒ_ｃおよび下記数式（１１）を用いてV^{vdW(rc→∞)}を計算し、メモリに格納する。

（ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎ_α、Ｎ_βはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数を表す。）
数式（１１）におけるスケールされたvdW相互作用として前述の様々な関数形を用いることができる。これまで、非結合性相互作用のスケール値λ_iを含むスケールされたvdW相互作用ポテンシャルにおける、カットオフ半径から無限遠のvdW相互作用ポテンシャルV^{vdW(rc→∞)}の計算方法が確立されていなかったが、本発明者らの検討の結果、数式（１０）のDreiding型レナードジョーンズポテンシャルを数式（１１）に代入して得られた解析解である下記数式（１２Ａ）および（１２Ｂ）からカットオフ半径から無限遠のvdW相互作用ポテンシャルV^{vdW(rc→∞)}を計算する方法が好ましいことが見出された。

〔ステップＤ〕
ステップＤは、ステップＡ、ステップＢおよびステップＣのすべてのステップの実行後に、前記数式（１）で示した全ポテンシャルエネルギーV^totalを計算するステップである。ステップＤにおいては、前記数式（８）により、ステップＢおよびステップＣで計算したカットオフ半径内のvdW相互作用ポテンシャルV^vdW(0→rc)とカットオフ半径から無限遠のvdW相互作用ポテンシャルV^{vdW(rc→∞)}の和から全vdW相互作用ポテンシャルV^vdWを求め、ステップＡで計算した{V^total (λ_i)‐V^vdW(r,λ_i)}とV^vdWとの和より全ポテンシャルエネルギーV^totalを計算し、メモリに格納する。

なお本発明において、数式（６）の右辺第６項であるクーロン相互作用の計算には、右辺第６項を展開し、無限遠のクーロン相互作用の寄与まで考慮できる方法であるEwald法やParticle Mesh Ewald(PME)法、Particle-Particle Particle-Mesh Ewald(PPPME)法、Fast Multipole Method(FMM)などを用いることが好ましい。

ステップＡ〜Ｄを実施した場合、カットオフ半径から無限遠のvdW相互作用ポテンシャルの寄与を考慮した全ポテンシャルエネルギーV^totalを計算することができる。後述の自由エネルギー計算では、全ポテンシャルエネルギーV^totalを用いるため、自由エネルギーにもカットオフ半径から無限遠のvdW相互作用ポテンシャルの寄与が考慮された、より高精度な値を求めることできる。

〔ステップ１〕
分子集合体および溶質分子の分子内力(数式（１）第１項から第４項の微分の和)および分子間力（数式（１）第５項と第６項の微分の和）を計算するステップ１においては、メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、ならびに数式（１）で示した全ポテンシャルエネルギーV^totalを原子の座標rで微分した式を用いて分子内力・分子外力を計算し、メモリに格納する。なお数式（１）における具体的な関数の形やパラメータは、ステップＡの項に記したものが好適に用いることができ、ステップＡを計算する場合は、ステップ１においても同じ関数形と同じパラメータを用いなければならない。

〔ステップ２〕
カットオフ半径内のvdW相互作用から生じる圧力Ｐ_vdW ^0→ｒｃを計算するステップ２では、ステップ１で計算したvdW力（数式（１）第５項の微分）と、原子の座標rとから、下記数式（１３Ａ）および（１３Ｂ）で表されるビリアル定理を用いて圧力Ｐ_vdW ^0→ｒｃを計算する。

（ここでF^vdWは注目する原子i、j間に働くvdW力、l, m(=1、2、3)はベクトルの成分を表す。）
〔ステップ３〕
カットオフ半径から無限遠のvdW相互作用から生じる圧力Ｐ_vdW ^rc→∞を計算するステップ３では、原子の座標r、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λ_iおよび相互作用計算のカットオフ半径ｒ_ｃから前記数式(４Ａ)および（４Ｂ）を用いて圧力Ｐ_vdW ^rc→∞計算する。

数式（４Ａ）の通り圧力Ｐ_vdW ^rc→∞は対角項のみで、非対角項は０となる。なお数式（４Ｂ）におけるスケールされたvdW相互作用として前述の様々な関数形を用いることができる。

前述のように、これまで、非結合性相互作用のスケール値λを含むポテンシャルでの圧力Ｐ_vdW ^rc→∞の計算方法が確立されていなかったが、本発明者らの検討の結果、数式（１０）のDreiding型レナードジョーンズポテンシャルを数式（４Ｂ）に代入して得られた解析解である下記数式（１４Ａ）および（１４Ｂ）と数式（４Ａ）から圧力Ｐ^rc→∞を計算する方法が好ましいことが見出された。

〔ステップ４〕
ステップ４は、ステップ２およびステップ３の実行後に、vdW相互作用から生じる圧力Ｐ_vdW ^0→∞を計算するステップである。ステップ４においては、数式（５）により、ステップ２およびステップ３で計算したカットオフ半径内のvdW相互作用から生じる圧力Ｐ_vdW ^0→ｒｃとカットオフ半径から無限遠のvdW相互作用から生じる圧力Ｐ^rc→∞の和から圧力Ｐ_vdW ^0→∞を計算する。

また数式（１）に、数式（４Ａ）および数式（４Ｂ）、数式（５）、数式（１３Ａ）および数式（１３Ｂ）より求めたＰ_vdW ^0→∞を代入することにより全圧力Ｐ_totalを求めることができる。なおvdW相互作用から生じる圧力Ｐ_vdW ^0→∞以外の圧力の詳細な表式もビリアル定理から導くことができ、参考文献［Glenn J. Martyna, Mark E. Tuckerman, Douglas J. Tobias, Michael L. Klein, Mol. Phys., 1996, 87, p.1117.］に記述されている。この文献の全圧力の表式はクーロン相互作用から生じる圧力Ｐ_coulomb ^0→∞を、無限遠のクーロン相互作用の寄与を考慮できる方法である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アンサンブルでシミュレーションを行う場合は、正確な密度における分子配置が得られるため、その全ポテンシャルエネルギーV^totalから高精度な自由エネルギーを求めることができる。

〔ステップ６〕
ステップ６はハミルトニアンレプリカ交換を行う場合のみ、レプリカ交換を実施する間隔lrepに１回実施する。レプリカの交換方法は、一般的に次に示すメトロポリスの方法が好適に用いられる。

i番目のレプリカとスケール値λ_i+1を持つｉ＋１番目のレプリカにおける、全ポテンシャルエネルギーV^total(λ_{i or i+1})と、下記数式（１５）よりi番目とi+1番目のレプリカのスケール値を交換する確率αを計算する。０〜１までの疑似乱数をコンピュータ上に発生させ、交換確率以下の値である場合にλ_iとλ_i+1との交換を行う。

（ここでk_Bはボルツマン定数、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ナノ秒）以上がより好ましい。

〔分子動力学シミュレーション結果の出力〕
本発明の計算方法においては、ステップ１〜ステップ５の実行、またはステップＡ〜Ｄおよびステップ１〜５の実行、またはステップＡ〜Ｄおよびステップ１〜６を所定回数繰り返し実行した後に、解析のために必要な様々な物理量を出力できる。この物理量は所定回数の繰り返し時点における瞬間値、あるいは一定回数繰り返して得られた平均値のどちらか一方または両方を出力することができる。

出力する物理量は後述する様々な自由エネルギー計算に必要なポテンシャルエネルギー以外にも、原子の座標、圧力などが挙げられる。ポテンシャルエネルギーについては自由エネルギー計算のために必須であり、原子の座標については、分子動力学シミュレーション最大の特徴である原子の時間経過による運動をとらえることができるという性質から、ある繰り返し回数における瞬間値を出力することが好ましい。

〔自由エネルギー計算〕
本発明における自由エネルギー計算法は、分子動力学シミュレーションで用いられる一般的な方法であれば特に限定されない。しかし、スケールされた非結合性相互作用を用いたときのシミュレーションであるという観点から、非現実的な中間状態を定義する錬金術的自由エネルギー計算法を用いるとき、本発明の効果が発揮される。

錬金術的自由エネルギー計算法として熱力学的積分法、自由エネルギー摂動法、ベネット受容比法またはマルチステイトベネット受容比法などが挙げられる。これらの手法は高精度に自由エネルギーを計算できることが知られており、どの手法も好適に用いられる。

自由エネルギー摂動法は、各中間状態におけるポテンシャルエネルギーV^total(λ_i) およびi番目の中間状態の分子配置におけるi+1番目の中間状態のスケール値λ_i+1を用いた場合の全ポテンシャルエネルギー<V^total(λ_i+1)>_iと、下記数式（１６Ａ）および数式（１６Ｂ）を用いて自由エネルギーを求める手法である。

（ここで、<…>はアンサンブル平均を、<…>_iはi番目のレプリカの分子配置を用いた物理量であることを表す。）
ベネット受容比法は、各中間状態における全ポテンシャルエネルギーV^total(λ_i)およびi番目の中間状態の分子配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iに加え、i+1番目の中間状態の分子配置とi番目またはi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i)>_i+1、V^total(λ_i+1)と、下記数式（１７Ａ）および数式（１７Ｂ）を用いて自由エネルギーを求める手法である。

（ここで、L_F、L_Rはそれぞれ、i番目の中間状態の分子配置とi+1番目の中間状態のスケール値を用いてサンプリングした全ポテンシャルエネルギーと、i+1番目の中間状態の分子配置とi番目の中間状態のスケール値を用いてサンプリングした全ポテンシャルエネルギーのデータ数を表す。さらにβ=1/k_BT、C= k_BTln(L_F/L_R)を表す。）
マルチステイトベネット受容比法は、i番目の中間状態の分子配置においてk番目(k=1,…, i,…, nint)の中間状態のスケール値λ_kを用いた時の全ポテンシャルエネルギー<V^total(λ_k)>_iと、下記数式（１８Ａ）および数式（１８Ｂ）を用いて自由エネルギーを求める手法である。

（ここで、x_nはN_k回分子配置をサンプリングしたときのn番目の分子配置を表す。）
前述の通り、中間状態を定義する錬金術的自由エネルギー計算法を用いるとき、本発明の効果が最大限に発揮されるが、特定の計算条件では中間状態を定義しない自由エネルギー計算法でも大きな効果がある。

例えば中間状態を定義せず、始状態と終状態のみ定義し、それらの相互作用エネルギーからなるエネルギー分布関数を求め近似汎関数を構築し、これより自由エネルギー計算するエネルギー表示法を、高分子分離膜中の分離物質のような分子配置のサンプリングが困難な系に適用する場合が挙げられる。高分子分離膜中の分離物質のような系では、高精度な自由エネルギーを得るため、分離物質の運動を促す目的でスケールされた非結合性相互作用を用いることがある。スケールされた非結合性相互作用を用いる場合でも、自由エネルギー計算に必要な分子配置のサンプリングは、始状態と終状態のみであり中間状態は関係ない。しかしNPTアンサンブルの計算を行った場合、中間状態の構造がシミュレーション的に不正確な構造となるため、始状態と終状態の構造にも影響を及ぼすこととなる。

〔プログラムおよび記録媒体〕
本発明のプログラムは、コンピュータの構成要素に上記の自由エネルギーの計算方法の各ステップにおける計算手段として動作させるものである。このプログラムを用いれば、上述の計算方法を実行することができる。

また、本発明の記録媒体は、当該プログラムをコンピュータ読み取り可能に記録したものである。この記録媒体のプログラムをコンピュータにて読み取って実行すれば、上述の計算方法を実行することができる。

以上の手順を実行すれば、 分子集合体に溶質分子が溶解している系における、スケールされた非結合性相互作用を用いた分子動力学シミュレーションにおいて、カットオフ半径から無限遠のvdW相互作用から生じる圧力Ｐ_vdW ^rc→∞を考慮した自由エネルギー計算を実行できないという課題を解決することができる。

以下に実施例を挙げて本発明を具体的に説明するが、これらは例示的なものであって限定的なものではない。

［実施例１］
ポリジメチルシロキサン（PDMS）鎖４本＋水１分子系に対し、分子動力学シミュレーションを実施した後、水のPDMSへの溶解自由エネルギー計算を行った。

分子動力学シミュレーションはFORTRAN90プログラムを用いて行い、分子軌道計算は汎用プログラムGaussian（登録商標）09を用いて行った。また、計算にはクロック周波数2.93GHzのCore i7（登録商標）870を搭載したコンピュータ、およびクロック周波数2.7HzのXeon（登録商標） E5-2680を搭載したコンピュータを用いた。

バイオビア社製の分子設計システムMaterials Studio（登録商標）8.0上で下記構造式（１）に示すポリマーを作成した。

構造式(1)において、n=40（分子量＝2996）とした。次いで、Materials Studio（登録商標）8.0 Amorphous Cell Moduleを用いて、表１に示した組成でPDMS分子と水分子を三次元周期境界条件下でランダムに配置し系を作成した。Materials Studio（登録商標）8.0の外部出力機能を使用し、原子の座標r 、結合情報、原子のポテンシャルパラメータタイプ、原子の質量ｍ、基本セルベクトルhが記載された分子構造ファイルを作成した。

次に、出力した分子構造ファイルに記載されたポテンシャルパラメータタイプに対応したポテンシャルパラメータを記述したパラメータファイルを作成した。

さらにテキストエディタを用いて中間状態数であるnint、相互作用計算のカットオフ半径ｒ_ｃ、図２計算フローにおけるループの繰り返しを規定する回数nstep、出力する物理量とその出力間隔、設定温度T₀、設定圧力P₀、熱浴の慣性量Q、barostatの慣性量Wを記述した計算条件ファイルを作成した。

まず作成した分子構造ファイル、パラメータファイル、計算条件ファイルを分子動力学シミュレーションプログラムで読み取らせ、データをメモリに格納し、中間状態数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）アンサンブルを構成することにより、系を設定温度T₀、設定圧力P₀に制御した状態で下記ステップＡ〜Ｄ、ステップ１〜５を実施した。なお、設定温度T₀=20℃、設定圧力P₀=1atm、熱浴の慣性量Q=120.95(kJ・ps²)/mol、barostatの慣性量W=483.79(kJ・ps²)/molとした。

ＮＰＴアンサンブルにおける分子動力学シミュレーションを実施するためには、下記ステップ５において原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新することが必要である。そのため分子構造ファイルおよび計算条件ファイルで入力していない、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvに初期値を与える必要がある。そこで原子の速度vの初期値としてT₀=20℃におけるマクスウェル分布に従うよう計算された速度を、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルの速度hvの初期値として全ての成分に0を与えメモリに格納した。

各中間状態における非結合性相互作用のスケール値λ_iをコンピュータ上で計算し、各中間状態のメモリ領域に格納した。λ_iは１をnintで等分した値を用い、nintが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ポテンシャルを含む、数式（１）で示したポテンシャル関数V^totalを用いて、V^vdWを除くポテンシャルエネルギー{V^total (λ_i)‐V^vdW(r,λ_i)}およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの{<V^total (λ_i+1)>_i‐<V^vdW(r,λ_i+1)>_i}を計算しメモリに格納した。また、数式（１０）中のδ＝1.0×10^-20ｍ^２とし、数式（１）第６項のクーロンポテンシャルおよびその力の計算にはEwald法を用い、Ewald法のパラメータとして逆空間のクーロンポテンシャルにα＝0.21オングストローム^-1、｜ｎ｜^２ _ｍａｘ＝50を用いて計算した。

ポテンシャルパラメータの設定は以下の通りである。まず、PDMS分子については、分子内ポテンシャルパラメータとして、結合長・結合角の平衡位置R⁰、θ⁰にはHF/6-31G(d,p)レベルの分子軌道計算によって決定した値を、力の定数K^bond、K^angle、K^torsion、K^inversionおよび平衡位置φ⁰、ω⁰にはAMBERパラメータ、hybrid/UAパラメータおよび桑島ポテンシャルパラメータを用いた。vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０にはAMBERパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。

水分子については、結合長・結合角の平衡位置R⁰、θ⁰にはTIP3Pパラメータ[W. L. Jorgensen 、R. W. Impey 、J. Chandrasekhar 、 J. D. Madura、M. L. Klein 、J. Chem. Phys.、1983、79、p.926.〕を用い、結合長・結合角の力の定数K^bond、K^angleはＩＲの実測値を再現するように決定した値を用い、vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０と電荷にはTIP3Pパラメータを用いた。

〔ステップＢ〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目およびi+1番目の中間状態における非結合性相互作用のスケール値λ_iおよびλ_i+1、ならびに相互作用計算のカットオフ半径ｒ_ｃを用いて、前記数式（１０）で示したDreiding型のvdWポテンシャルを用いて、原子間距離が0からカットオフ半径r_c以内のvdW相互作用のポテンシャルエネルギーV^vdW(0→rc)(r,λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときのポテンシャルエネルギー<V^vdW(0→rc)(r,λ_i+1)>_iを計算しメモリに格納した。

〔ステップＣ〕
メモリに格納された結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目およびi+1番目の中間状態における非結合性相互作用のスケール値λ_iおよびλ_i+1、相互作用計算のカットオフ半径rc、前記数式（１２Ａ）、（１２Ｂ）を用いて、カットオフ半径r_cから無限遠までのvdW相互作用のポテンシャルエネルギーV^{vdW(rc→∞)}(r,λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときのポテンシャルエネルギー<V^{vdW(rc→∞)}(λ_i+1)>_iを計算しメモリに格納した。

〔ステップＤ〕
ステップＡ〜Ｃでメモリに格納した各成分のポテンシャルエネルギーと前記数式（１）より全ポテンシャルエネルギーV^total(λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iを計算しメモリに格納した。

〔ステップ１〕
メモリに格納された原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、i番目の中間状態における非結合性相互作用のスケール値λ_i、ならびに、前記数式（１）における全ポテンシャル関数V^total(λ_i)を原子の座標rで微分した式を用いて分子内力・分子外力を計算しメモリに格納した。

〔ステップ２〕
ステップ１でメモリに格納した分子間力の一つであるvdW力（数式（１０）の微分）と、原子の座標rからビリアル定理（数式（１３Ａ）および（１３Ｂ））を用いて圧力Ｐ_vdW ^0→ｒｃを計算しメモリに格納した。

〔ステップ３〕
結合情報、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、i番目の中間状態における非結合性相互作用のスケール値λ_iおよび相互作用計算のカットオフ半径ｒ_ｃと数式（４Ａ）、数式（１４Ａ）および数式（１４Ｂ）から圧力Ｐ_vdW ^rc→∞を計算しメモリに格納した。

〔ステップ４〕
ステップ２および３でメモリに格納した、カットオフ半径内のvdW相互作用から生じる圧力Ｐ_vdW ^0→ｒｃと、カットオフ半径から無限遠のvdW相互作用から生じる圧力Ｐ_vdW ^rc→∞と、前記数式（５）により圧力Ｐ_vdW ^0→∞を計算しメモリに格納した。

〔ステップ５〕
メモリに格納された原子の質量ｍ、原子の座標r、原子の速度v、熱浴の慣性量Q、熱浴の位置ζ、熱浴の速度ζv、barostatの慣性量W、基本セルベクトルh、基本セルベクトルの速度hv、分子内力、分子間力を用いて、微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを数値積分によって計算し、更新した。数値積分には多重時間刻み数値積分法（RESPA法）を用い、時間刻みΔtは表２に示した値を用いた。

ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを初期値としてステップＡ〜Ｄ、ステップ１〜５を再度実行させることをnstep＝16000000回繰り返した。RESPA法における最小の時間刻みが0.25フェムト秒であるため、この繰り返し回数は4ナノ秒に相当する。またステップ７を実施する間隔である8000回は２ピコ秒に相当する。

なお本実施例ではNose−Hoover法により温度、Parrinello−Rahman法により圧力を制御しているためスッテプ６において熱浴の位置ζ、熱浴の速度ζv、基本セルベクトルh、基本セルベクトルの速度hvを更新した。

ステップＡ〜Ｄ、ステップ１〜５を4000回（１ピコ秒に相当）繰り返すごとに、原子の座標r、各レプリカにおけるポテンシャルエネルギーV^total(λ_i) およびi番目のレプリカの分子配置におけるi+1番目のレプリカのスケール値λ_i+1を用いた場合の全ポテンシャルエネルギー<V^total(λ_i+1)>_iを出力した。

分子動力学シミュレーション終了後、V^total(λ_i)および<V^total(λ_i+1)>_iの時系列データと前記数式（１６Ａ）、（１６Ｂ）に示した自由エネルギー摂動法を用いて、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを計算した。

なお、本実施例の計算では構造緩和のため、上述の計算の前に20℃の原子数・温度・体積一定（ＮＶＴ）アンサンブルの分子動力学シミュレーションを10ピコ秒、1atm／20℃のＮＰＴアンサンブルの分子動力学シミュレーションを90ピコ秒、1atmの分子動力学シミュレーションにおいて温度を20℃〜500℃まで200ピコ秒周期で上昇下降を繰り返し緩和させるシミュレーテッドアニーリングを14サイクル（2800ピコ秒に相当）、1atm／20℃のＮＰＴアンサンブルの分子動力学シミュレーションを300ピコ秒実施した。さらにステップ１を始める前に、構造緩和計算後の構造を用いて各中間状態に対応するλ_iを用いて分子動力学シミュレーションを300ピコ秒実施し、中間状態での構造緩和を行った。

［実施例２］
次に示す条件以外は実施例１と同じ分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを計算した。

計算条件ファイルにハミルトニアンレプリカ交換を行う中間状態数nint＝21（＝レプリカ数）、交換を実施する間隔lrep＝8000回に１回と記述しメモリに格納した。実施例１における中間状態をレプリカと読み替え、各レプリカにおいて独立に実施例１のステップＡ〜Ｄおよびステップ１〜５を実施した。

すべてのレプリカにおいてステップ５までの計算が終わった時点で、これまでステップ１〜５まで8000回繰り返してきたならば、1番目のレプリカからnint番目のレプリカまで順に次のステップ６を実施する。それ以外の場合はステップ６を実施しない。

〔ステップ６〕
i番目のレプリカとスケール値λ_i+1を持つｉ＋１番目のレプリカにおける、ポテンシャルエネルギーV^total(λ_{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％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^solの計算を行った。なお計算条件は、下記の条件以外は実施例１と同様である。

Materials Studio（登録商標）8.0 Amorphous Cell Moduleを用いて、表３に示した組成でメタノール分子と水分子を三次元周期境界条件下でランダムに配置し系を作成した。

計算条件ファイルにおいて設定したパラメータとして、ループの繰り返しを規定する回数nstep=4000000回とし、この繰り返し数は1ナノ秒に相当する。

メタノール分子のポテンシャルパラメータの設定は以下の通りである。分子内ポテンシャルパラメータとして、結合長・結合角の平衡位置R⁰、θ⁰にはHF/6-31G(d,p)レベルの分子軌道計算によって決定した値を、力の定数K^bond、K^angle、K^torsion、K^inversionおよび平衡位置φ⁰、ω⁰にはAMBERパラメータおよび桑島ポテンシャルパラメータを用いた。vdWポテンシャルの深さと斥力が生じる位置を決めるパラメータD^０、R^０にはAMBERパラメータ、および桑島ポテンシャルパラメータを用い、電荷はB3LYP/6-31G(d,p)レベルの分子軌道計算によって決定した値を用いた。

［実施例４］
次に示す条件以外は実施例２と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを計算した。

ステップＡ〜Ｄにおいて全ポテンシャルエネルギーV^total(λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i-1)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V^total(λ_i)、<V^total(λ_i+1)>_i、<V^total(r,λ_i-1)>_iの時系列データと数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^polyを計算した。

［実施例５］
次に示す条件以外は実施例３と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^solを計算した。

スステップＡ〜Ｄにおいて全ポテンシャルエネルギーV^total(λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i-1)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V^total(λ_i)、<V^total(λ_i+1)>_i、<V^total(r,λ_i-1)>_iの時系列データと数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて溶解自由エネルギーΔG_solv ^solを計算した。

［実施例６］
次に示す条件以外は実施例２と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを計算した。

ステップＡ〜Ｄにおいて、i番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λ_kを用いたときの全ポテンシャルエネルギー<V^total(λ_k)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V<V^total(λ_k)>_iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^polyを計算した。

［実施例７］
次に示す条件以外は実施例３と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^solを計算した。

ステップＡ〜Ｄにおいて、i番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λ_kを用いたときの全ポテンシャルエネルギー<V^total(λ_k)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V<V^total(λ_k)>_iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^solを計算した。

［比較例１］
次に示す条件以外は実施例２と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを自由エネルギー摂動法により計算した。

ステップＡ〜Ｂを実施した後、ステップＣを実施せずにステップＤを実施した。またステップ１およびステップ２を実施した後、ステップ３を実施せず、ステップ４およびステップ５を実施した。すべてのレプリカにおいてステップ５までの計算が終わった時点で、これまでステップ１〜５まで8000回繰り返してきたならば、1番目のレプリカからnint番目のレプリカまで順に次のステップ６を実施した。以後ステップＡ〜Ｂ、ステップＤ、ステップ１〜２、ステップ４〜５を16000000回繰り返し実施し、繰り返しを実施した数が8000で割り切れるときはステップ５の後にステップ６を実施した。

ステップＣおよびステップ３を実施しないということは、カットオフ半径から無限遠までのvdW相互作用から生じるポテンシャルエネルギーと圧力を無視することに相当する。

［比較例２］
次に示す条件以外は実施例３と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^solを自由エネルギー摂動法により計算した。

ステップＡ〜Ｂを実施した後、ステップＣを実施せずにステップＤを実施した。またステップ１およびステップ２を実施した後、ステップ３を実施せず、ステップ４およびステップ５を実施し、すべてのレプリカにおいてステップ５までの計算が終わった時点で、これまでステップ１〜５まで8000回繰り返してきたならば、1番目のレプリカからnint番目のレプリカまで順に次のステップ６を実施した。以後ステップＡ〜Ｂ、ステップＤ、ステップ１〜２、ステップ４〜５を4000000回繰り返し実施し、繰り返しを実施した数が8000で割り切れるときはステップ５の後にステップ６を実施した。

［比較例３］
次に示す条件以外は比較例１と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを計算した。

ステップＡ〜Ｂ、ステップＤにおいて全ポテンシャルエネルギーV^total(λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i-1)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V^total(λ_i)、<V^total(λ_i+1)>_i、<V^total(λ_i-1)>_iの時系列データと前記数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^polyを計算した。

［比較例４］
次に示す条件以外は比較例２と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^solを計算した。

ステップＡ〜Ｂ、ステップＤにおいて全ポテンシャルエネルギーV^total(λ_i)およびi番目の中間状態の配置とi+1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iに加え、i番目の中間状態の配置とi-1番目の中間状態のスケール値を用いたときの全ポテンシャルエネルギー<V^total(λ_i-1)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V^total(λ_i)、<V^total(λ_i+1)>_i、<V^total(λ_i-1)>_iの時系列データと前記数式（１７Ａ）、（１７Ｂ）に示したベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^solを計算した。

［比較例５］
次に示す条件以外は比較例１と同じハミルトニアンレプリカ交換分子動力学シミュレーションを実施し、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyを計算した。

ステップＡ〜Ｂ、ステップＤにおいてi番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λ_kを用いたときの全ポテンシャルエネルギー<V^total(λ_k)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V<V^total(λ_k)>_iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^polyを計算した。

［比較例６］
次に示す条件以外は比較例２と同じ分子動力学シミュレーションを実施し、95wt％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^solを計算した。

ステップＡ〜Ｂ、ステップＤにおいてi番目の中間状態の配置とk番目(k=1,…, i,…, nint)の中間状態のスケール値λ_kを用いたときの全ポテンシャルエネルギー<V^total(λ_k)>_iを計算し、メモリに格納した。

分子動力学シミュレーション終了後、V<V^total(λ_k)>_iの時系列データと数式（１８Ａ）、（１８Ｂ）に示したマルチステイトベネット受容比法を用いて、溶解自由エネルギーΔG_solv ^solを計算した。 実施例１〜７および比較例１〜６で得られた自由エネルギーを表４にまとめた。なお表４中の自由エネルギー計算法では、自由エネルギー摂動法をFEP、ベネット受容比法をBAR、マルチステイトベネット受容比法をMBARと表した。

ここで、PDMSに水が１分子溶解するときの溶解自由エネルギーΔG_solv ^polyと95wt％メタノール水溶液への水分子の溶解自由エネルギーΔG_solv ^sol、下記数式（１９）より95wt％メタノール水溶液からPDMS膜に水が１分子溶解したときの自由エネルギーΔG_sorpを求めることができる。

このΔG_sorpは高分子膜を水溶液に含浸させるという簡便な実験より求めることができるため、シミュレーション結果と直接比較でき、シミュレーションの精度確認の指標に適している。非特許文献〔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).〕より吸水率と下記数式（２０）を用いることによりΔG_sorpの実験値を求めることができる。

（ここでWは吸水率、ρは密度、k_Bはボルツマン定数、Tは絶対温度をそれぞれ表す。）
数式（２０）よりΔG_sorpの実験値は3.76kcal/molと求められる。なおポリマーの密度は0.971g/cm3、９５ｗｔ％メタノール水溶液の密度は、水の密度1.0g/cm3およびメタノールの密度0.792g/cm3と混合比率から求めた0.802g/cm3を用いた。計算値および実験値のΔG_sorpも併せて表４に示した。なお実施例１〜３におけるΔG_sorpの（）外の値は、実施例２と実施例３から求めた値で、（）内の値は実施例１と実施例３から求めた値である。すなわち（）内の値はPDMS・水系においてハミルトンレプリカ交換を行わなかった場合のΔG_sorpを表す。

実施例１および実施例３より求めたΔG_sorpと、実施例２および実施例３より求めたΔG_sorpを比較すると、実施例２および実施例３より求めた方が実験値に近いことがわかる。これは、実施例１においてはハミルトニアンレプリカ交換を行わなかったことにより、スケール値λ_iが大きいときにPDMS中の水の運動性が悪くなり、局所的な分子配置しかサンプリングできなかった影響であると考える。

また、ベネット受容比法を用いた実施例４および実施例５より求めたΔG_sorpとマルチステイトベネット受容比法を用いた実施例６および実施例７より求めたΔG_sorpは、実施例２および実施例３のエネルギー摂動法より求めた値よりも実験値に近く、非常に精度良く溶解自由エネルギーを計算できているといえる。一般的にエネルギー摂動法よりもベネット受容比法およびマルチステイトベネット受容比法の方が高い計算精度を持つといわれており、今回の結果もそれを支持する。

カットオフ半径から無限遠のvdW相互作用を無視した比較例１〜６までの結果を用いた場合のΔG_sorpと実施例の結果を比較する。比較例１および比較例２または比較例３および比較例４または比較例５および比較例６から求めたΔG_sorpは、実験値から約0.6kcal/mol以上かい離し、実施例の結果よりも計算精度が低い。したがって本発明である自由エネルギー計算手法は、分子集合体に溶質分子が溶解している系において非常に高精度な自由エネルギーを与えるといえる。

Claims (11)

1. コンピュータを利用した、分子集合体に溶質分子が溶解している系の分子動力学シミュレーションにおいて、下記数式（３）のように、
   （i）溶質分子の分子内vdW相互作用V^vdw _sol-intra(r, λ_sol-intra)と、
   （ii）分子集合体に属する分子の分子内vdW相互作用V^vdw _mol-intra(r, λ_mol-intra)と、
   （iii）溶質分子‐溶質分子間vdW相互作用V^vdw _sol-sol(r, λ_sol-sol)と、
   （iv）分子集合体に属する分子の分子間vdW相互作用V^vdw _mol-mol(r, λ_mol-mol)と、
   （v）溶質分子‐分子集合体に属する分子の分子間vdW相互作用V^vdw _sol-mol(r, λ_sol-mol)
   との和（ただし0≦λ≦1）で表され、非結合性相互作用のスケール値λを含むスケールファンデルワールス（vdW）相互作用ポテンシャルV^vdW(r, λ)を用いた場合の自由エネルギー計算法であって、
   

   

   コンピュータに、
   ・前記溶質分子と分子集合体に含まれる原子の座標r
   および、必要に応じ
   ・原子の速度v
   の初期値、
   ならびに
   ・結合情報
   ・原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ
   ・原子の質量m
   ・相互作用計算のカットオフ半径r_c
   ・中間状態を定義する非結合性相互作用のスケール値λ_i（λ₁、…λ_i…、λ_nint、0〜1を中間状態の数nintで分割した値）
   を入力した後、
   非結合性相互作用のスケール値λ_iをもつ中間状態iそれぞれにおいて、コンピュータ上で少なくとも、
   （ステップ１）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、
   下記数式（１）で示した全ポテンシャルエネルギーV^total(λ_i) を原子の座標rで微分した式を用いて、前記分子集合体および前記溶質の分子内力(数式（１）第１項から第４項の微分の和)および分子間力（数式（１）第５項と第６項の微分の和）を計算するステップ；
   

   

   （ただし、R、θ、φ、ωはそれぞれ注目する結合の距離、変角の角度、二面角の角度、反転の角度を表す。）
   （ステップ２）ステップ１で計算した分子間力の一つであるvdW力（前記数式（１）第５項の微分）と、
   原子の座標rからビリアル定理を用いて圧力Ｐ_vdW ^0→ｒｃを計算するステップ；
   （ステップ３）結合情報、原子のポテンシャルパラメータタイプ、ポテンシャルパラメータ、非結合性相互作用のスケール値λ_iおよび相互作用計算のカットオフ半径ｒ_ｃから、
   下記数式(４Ａ)および（４Ｂ）を用いて圧力Ｐ_vdW ^rc→∞を計算するステップ；
   

   

   （ただし、Ｍはポテンシャルパラメータタイプの総数、Ｎ_α、Ｎ_βはそれぞれポテンシャルパラメータタイプα、βを持つ原子の数、diag[・・・]は対角行列を表す。）
   （ステップ４）ステップ２およびステップ３の実行後に、下記数式（５）によりvdW相互作用から生じる圧力Ｐ_vdW ^0→∞を計算するステップ；
   

   

   （ステップ５）ステップ４の実行後に、原子の質量ｍ、分子内力、分子間力を用いて、微小時間Δt後の
   ・原子の座標r
   ・原子の速度ｖ
   を数値積分によって更新するステップ：
   を実行させ、
   さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度ｖを初期値としてステップ１〜ステップ５を再度実行させ、これを繰り返すことを含む、自由エネルギーの計算方法。
2. ステップ１〜５のいずれかのステップの前または後に、
   （ステップＡ）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_iを用いて、
   前記数式（１）における全ポテンシャルエネルギーV^total(λ_i)の中で、V^vdW(r,λ_i)以外のポテンシャルエネルギーを計算するステップ：
   （ステップＢ）原子の座標r、結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、相互作用計算のカットオフ半径ｒ_ｃを用いて、V^vdW(r,λ_i)における、原子間距離が0からカットオフ半径r_c以内のvdW相互作用ポテンシャルV^vdW(0→rc)を計算するステップ：
   （ステップＣ）結合情報、原子のポテンシャルパラメータタイプおよびポテンシャルパラメータ、非結合性相互作用のスケール値λ_i、相互作用計算のカットオフ半径rc、下記数式（１１）を用いて、V^vdW(r,λ_i)における、カットオフ半径r_cから無限遠までのvdW相互作用ポテンシャルV^{vdW(rc→∞)}を計算するステップ：
   

   

   （ステップＤ）ステップＡ、ステップＢおよびステップＣのすべてのステップの実行後に、ステップＢおよびステップＣで計算したV^vdW(0→rc)とV^{vdW(rc→∞)}と下記数式（８）よりV^vdW(r,λ_i)を求め、V^vdW(r,λ_i)とステップＡで計算した各エネルギー成分のポテンシャルエネルギーと前記数式（２）より全ポテンシャルエネルギーV^totalを計算するステップ：
   

   

   を実行することを含む請求項１に記載の自由エネルギー計算方法。
3. 中間状態を定義する非結合性相互作用のスケール値λ_iを、スケール値を除き全く同じ入力情報を持つnint個のレプリカに割り当て、
   それぞれのレプリカにおいて独立にステップＡ〜Ｄおよびステップ１〜５を実行後、
   （ステップ６）ｉ番目のレプリカと隣り合う、スケール値λ_i+1を持つｉ＋１番目のレプリカに、
   ステップＤで求めた全ポテンシャルエネルギーV^totalとメトロポリスの方法による判定基準を用い、
   基準を満たす場合はλ_iとλ_i+1との交換を行うステップ：
   を実行させ、
   さらに、ステップ５で更新した微小時間Δt後の原子の座標r、原子の速度v、
   ステップ６でメトロポリスの判定基準を満たした場合は交換されたλを初期値にステップＡ〜Ｄおよびステップ１〜５を再度実行させ、これを繰り返すことを含む、請求項２に記載の自由エネルギーの計算方法。
4. 数式（１）におけるスケールされたvdW相互作用ポテンシャルV^vdW(r, λ)が数式（１０）に示したDreiding型レナードジョーンズポテンシャルであり、数式（１２Ａ）および（１２Ｂ）を用いてカットオフ半径rcから無限遠までのvdW相互作用ポテンシャルV^{vdW(rc→∞)}を、数式（４Ａ）、数式（１４Ａ）および（１４Ｂ）を用いてＰ_vdW ^rc→∞を計算することを特徴とする、請求項２または３に記載の自由エネルギーの計算方法。
   

   

   （ここで注目する二つの原子をj,kとし、δはエネルギーと力の距離依存性を調節するパラメータ、D^０、R^０はvdWポテンシャルの深さと斥力が生じる位置を決めるパラメータ、r² _jkはij間の距離とし、0≦λ_i≦1、δ≧0ｍ^２とする。）
   

   

   
5. 非結合性相互作用のスケール値λ_iをもつ中間状態iの計算のステップＡ〜Ｄにおいて、
   スケール値λ_iのときの全ポテンシャルエネルギーV^total(λ_i)に加え、同じ分子配置を用いたスケール値λ_i+1のときの全ポテンシャルエネルギー<V^total(λ_i+1)>_iを求め、
   自由エネルギー摂動法による自由エネルギー計算を実施することを特徴とする、
   請求項２〜４のいずれか記載の自由エネルギーの計算方法。
6. 非結合性相互作用のスケール値λ_iをもつ中間状態iの計算のステップＡ〜Ｄにおいて、
   スケール値λ_iのときの全ポテンシャルエネルギーV^total(λ_i)、スケール値λ_iのときの分子配置を用いたスケール値λ_i+1における全ポテンシャルエネルギー<V^total(λ_i+1)>_i、スケール値λ_iのときの分子配置を用いたスケール値λ_i-1における全ポテンシャルエネルギー<V^total(λ_i-1)>_iを求め、
   ベネット受容比法による自由エネルギー計算を実施することを特徴とする、
   請求項２〜４のいずれか記載の自由エネルギーの計算方法。
7. 非結合性相互作用のスケール値λ_iをもつ中間状態iの計算のステップＡ〜Ｄにおいて、スケール値λ_iのときの分子配置を用いたスケール値λ_k(k=1,…, i,…, nint)における全ポテンシャルエネルギー<V^total(λ_k)>_iを求め、
   マルチステイトベネット受容比法による自由エネルギー計算を実施することを特徴とする、
   請求項２〜４のいずれか記載の自由エネルギーの計算方法。
8. 分子集合体として高分子、溶質として低分子を用いることを特徴とする請求項１〜７のいずれか記載の自由エネルギーの計算方法。
9. 非結合性相互作用のスケール値をλ_sol-intra＝λ_mol-intra＝λ_sol-sol＝λ_mol-mol＝１、λ_sol-mol=λ_iとすることを特徴とする請求項１〜８のいずれか記載の自由エネルギーの計算方法。
10. 請求項１〜９のいずれか記載の自由エネルギーの計算方法をコンピュータに実行させるためのプログラム。
11. 請求項１０に記載のプログラムをコンピュータに読み取り可能に記録した記録媒体。
    

Images