JP2009080803A - 親和性計算方法 - Google Patents

Abstract

【課題】分子集合体と溶質分子との親和性を高速に計算するための方法を提供する。
【解決手段】分子集合体と溶質との親和性を求めるためのシミュレーションにおいて、分子集合体の初期座標データを計算機に格納するステップ、分子集合体の初期座標データに対して計算を実行し、少なくとも１つのスナップショットを出力するステップ、少なくとも１つ選択したスナップショットに対してランダムな初期位置に溶質分子を挿入するステップ、溶質分子が初期位置にあるときのポテンシャルエネルギー差を計算するステップ、計算機に保存するステップ、溶質分子位置を変化させるステップ、溶質分子のポテンシャルエネルギー差を計算するステップ、ΔＶを計算するステップ、変化させた溶質分子の位置を判定するステップ、計算機に保存するステップ、化学ポテンシャルを計算するステップ、化学ポテンシャルの値を出力するステップとを有する。
【選択図】なし

Description

本発明は、電子計算機を利用して分子集合体と溶質との親和性を高速に計算するための方法に関する。

分子動力学法やモンテカルロ法といった分子シミュレーション手法は、分子あるいは分子集合体について、実験では困難な原子レベルでの詳細な構造や運動に関する知見を得ることができる方法として広く用いられており、近年の電子計算機の飛躍的向上と方法論の発展に伴ってその適用分野を広めつつある。とりわけ、ポリオレフィンやポリエステルをはじめとする合成高分子、タンパク質やＤＮＡといった巨大分子への適用が盛んに試みられている。本明細書においては、分子集合体モデルに対して分子シミュレーション手法を適用することを分子集合体計算と呼ぶ。

一般に、分子シミュレーション手法で算出可能な物性として、（ｉ）分子あるいは分子集合体の高次構造、（ｉｉ）分子あるいは分子集合体の電子状態、（ｉｉｉ）分子集合体中での分子の拡散性、（ｉｖ）分子集合体と分子との親和性などを挙げることが出来る。（ｉ）、（ｉｉ）および（ｉｉｉ）については、以前より盛んに研究が行われており、信頼性のある結果を導き出せるようになりつつある。

これらに対して、（ｉｖ）は、分子集合体中での分子の透過性や分子の溶媒への溶解性などを計算し、フィルムや樹脂といった高分子成形体の溶媒への溶解性・耐薬品性・水蒸気透過性・ガス透過性を評価するためには欠かせない重要な量でありながら、現状ではシミュレーション方法の開発が十分に進んでいるとは言えない状況である。

分子集合体と分子との親和性を求めるためには、自由エネルギー計算を必要とする。一般に、分子シミュレーションにおいて、自由エネルギー計算は多大な計算時間を必要とすることが問題となっている。そのため、分子集合体中における分子の親和性について、分子シミュレーションを適用した例は極めて少ないのが現状であり、自由エネルギーを高速に計算することのできる方法が求められていた。

代表的な自由エネルギー計算方法として、熱力学的積分法（非特許文献１）、自由エネルギー摂動法（非特許文献２）、粒子挿入法（非特許文献３）を挙げることができる。

熱力学的積分法は、溶質分子を分子集合体中に徐々に浸透させながら、溶質分子に働くポテンシャルエネルギー差を求め、ポテンシャルエネルギー差と溶質分子の変位との積から自由エネルギーを求める方法である。

自由エネルギー摂動法は、分子集合体と溶質分子が隣接した状態と、分子集合体中に溶質分子が存在している状態とのポテンシャルエネルギー差から自由エネルギーを求める方法である。

これらの方法はいずれも非常に多くの計算時間を必要とするだけでなく、計算対象となる系の性質に応じて、ポテンシャルエネルギー差の計算方法や計算モデルの作成方法に種々のノウハウを必要とし、容易に実行できる方法ではない。

粒子挿入法は上記の方法に比べて簡便に自由エネルギーを求めることのできる方法として知られている。この方法は、温度Ｔ_１における溶質分子１個あたりの自由エネルギー、つまり化学ポテンシャルを、数式（４）に示す計算式から求める方法である。

ここで、Ｋ_Ｂはボルツマン定数、Ｎは溶質分子の位置を変化させる回数（サンプリング数という）、ｉはステップ番号（ｉ＝１、２、・・・Ｎ）を表す。溶質分子が位置ｒ_ｉにあるときのポテンシャルエネルギー差Ｖ（ｒ_ｉ）は、数式（５）に示す計算式から求める。

ここで、ｖ（ｒ_ｉ）は位置ｒ_ｉにおける溶質分子と分子集合体を合わせた全ポテンシャルエネルギー、ｖ_０分子集合体のみのポテンシャルエネルギーを表す。

前記粒子挿入法の具体的な計算手順を図１に示す。

分子集合体の初期座標データを電子計算機に格納し（ステップ１）、ステップ１で与えた分子集合体の初期座標データに対して分子集合体計算を実行し、少なくとも１つ以上のスナップショットを出力する（ステップ２）。次に、ステップ２で出力したスナップショットから少なくとも１つを選択し、選択したスナップショットに対してランダムな初期位置ｖ_０に溶質分子を挿入する（ステップ３）。その後、溶質分子が初期位置ｖ_０にあるときのポテンシャルエネルギー差Ｖ（ｒ_０）を計算し（ステップ４）、ｒ_０、Ｖ（ｒ_０）を電子計算機に保存する（ステップ５）。次いで、溶質分子を元の位置ｒ_ｉ−１から次の位置ｒ_ｉにランダムに変化させ（ステップ６）、溶質分子が位置ｒ_ｉにあるときのポテンシャルエネルギー差Ｖ（ｒ_ｉ）を計算する（ステップ７）。

その後、ΔV =Ｖ（ｒ_ｉ） -Ｖ（ｒ_ｉ−１）を計算し（ステップ８）、ｒ_ｉ、Ｖ（ｒ_ｉ）を電子計算機に保存する（ステップ９）。その後、化学ポテンシャルを計算し（ステップ１０）、値を出力する（ステップ１１）する。なお、分子集合体計算の終了後にステップ３を行い、ステップ６〜１０を規定回数に到達するまで繰り返し行う。

この方法では、溶質分子の位置ｒ_ｉをランダムに選ぶため、分子集合体と溶質分子が反発しあう位置を高頻度で発生させる。実際の計算において、サンプリングできる数は有限である。そのため、本来サンプリングされるはずの、溶質分子が分子集合体中に安定に存在する位置をほとんど発生させることができず、化学ポテンシャルが一定の値に収束しない。化学ポテンシャルを一定の値に収束させるためには、サンプリング数を増加させ、溶質分子が分子集合体中に安定に存在する位置を多く発生させる必要がある。ところが、この場合、非常に多くのサンプリング数を必要とするため計算時間が非常に大きくなるという問題がある。
エム・ピー・アレン （Ｍ. Ｐ. Ａｌｌｅｎ）、 ディー・ジェイ・ティルデスリー (Ｄ. Ｊ. Ｔｉｌｄｅｓｌｅｙ)著、コンピューター シミュレーション オブ リキッズ （Ｃｏｍｐｕｔｅｒ Ｓｉｍｕｌａｔｉｏｎ ｏｆ Ｌｉｑｕｉｄｓ）、１９８７、ｐ.２１９. アール・ダブリュー・ツヴァンツィッヒ （Ｒ. Ｗ. Ｚｗａｎｚｉｇ）、 ジャーナル オブ ケミカル フィジクス （Ｊ. Ｃｈｅｍ. Ｐｈｙｓ）、１９５４、２２、ｐ.１４２０. ビー・ワイダム （Ｂ. Ｗｉｄｏｍ）、 ジャーナル オブ ケミカル フィジクス （Ｊ. Ｃｈｅｍ. Ｐｈｙｓ）、１９６３、３９、ｐ.２０８０.

このように、従来の自由エネルギー計算方法は多大な計算時間を必要とする一方、計算時間を短縮し、効率的に自由エネルギーひいては親和性を計算するための方法は未だ提唱されていない。そこで本発明は、上記の問題を解決し、分子集合体と溶質分子との親和性を高速に計算するための方法を提供することを目的とする。

上記問題を解決するため、本発明は、分子集合体と溶質との親和性を求めるための分子シミュレーションにおいて、（ａ）分子集合体の初期座標データを電子計算機に格納するステップ（ステップ１）と、（ｂ）ステップ１で与えた分子集合体の初期座標データに対して分子集合体計算を実行し、少なくとも１つ以上のスナップショットを出力するステップ（ステップ２）と、（ｃ）ステップ２で出力したスナップショットから少なくとも１つを選択し、選択したスナップショットに対してランダムな初期位置ｒ_０に溶質分子を挿入するステップ（ステップ３）と、（ｄ）溶質分子が初期位置ｒ_０にあるときのポテンシャルエネルギー差Ｖ（ｒ_０）を計算するステップ（ステップ４）と、（ｅ）ｒ_０、Ｖ（ｒ_０）を電子計算機に保存するステップ（ステップ５）と、（ｆ）溶質分子を元の位置ｒ_ｉ−１から次の位置ｒ_ｉに変化させるステップ（ステップ６）と、（ｇ）溶質分子が位置ｒ_ｉにあるときのポテンシャルエネルギー差Ｖ（ｒ_ｉ）を計算するステップ（ステップ７）と、（ｈ）Ｖ（ｒ_ｉ） とＶ（ｒ_ｉ−１）の差であるΔＶを計算するステップ（ステップ８）と、（ｉ） ステップ６で変化させた溶質分子の位置ｒ_ｉを、判定基準を用いて採択または棄却を判定するステップ（ステップ９）と、（ｊ）ｒ_ｉ、Ｖ（ｒ_ｉ）を電子計算機に保存するステップ（ステップ１０）と、（ｋ）化学ポテンシャルを計算するステップ（ステップ１１）と、（ｌ）化学ポテンシャルの値を出力するステップ（ステップ１２）とを有し、分子集合体計算の終了後にステップ３を行い、ステップ６〜１０を規定回数に到達するまで繰り返し行うことを特徴とする計算方法であることを本旨とし、また、種々の改良された態様を提案する。

本発明の親和性計算方法により、従来は多大な計算時間を必要とした自由エネルギーを高速に計算することが可能となり、新規のフィルムや樹脂などの高分子成形体、新規の薬剤などの開発においてスピードアップ、ブレイクスルーが期待できる。

従来の粒子挿入法を用いて分子集合体と溶質分子との化学ポテンシャルを計算する場合、分子集合体中における溶質分子の位置をランダムに変化させるため、化学ポテンシャルの値が収束するには非常に多くのサンプリング数を必要とし、長大な計算時間を必要とする。

本発明者らは、分子集合体中に挿入した溶質分子の位置をメトロポリスのアルゴリズム（エヌ・メトロポリス （Ｎ.Ｍｅｔｒｏｐｏｌｉｓ）、 エー・ダブリュー・ローゼンブルース （Ａ. Ｗ. Ｒｏｓｅｎｂｌｕｔｈ）、 エム・エヌ・ローゼンブルース （Ｍ. Ｎ. Ｒｏｓｅｎｂｌｕｔｈ）、 エー・エイチ・テラー （Ａ. Ｈ. Ｔｅｌｌｅｒ）、 イー・ジェイ・テラー （Ｅ. Ｊ. Ｔｅｌｌｅｒ）、 ジャーナル オブ ケミカル フィジクス （Ｊ. Ｃｈｅｍ. Ｐｈｙｓ）、１９５３、２１、ｐ.１０８７.）に従うように変化させ、溶質分子が分子集合体中に安定に存在する位置を高頻度に発生させることで化学ポテンシャルの値が速やかに一定の値に収束することに気づき、また、メトロポリスのアルゴリズムにおける最適な温度範囲を見出した。さらに、分子集合体中で溶質分子の存在しうる位置をスナップショット内の自由体積に限定した状態で、メトロポリスのアルゴリズムを実行することにより、化学ポテンシャルの値がより速く収束することを見出し、本発明に到達したものである。

本発明は以下の構成からなる。すなわち、
（１）分子集合体と溶質との親和性を求めるための分子シミュレーションにおける、自由エネルギーの計算方法であって、
（ａ）分子集合体の初期座標データを電子計算機に格納するステップ（ステップ１）と、
（ｂ）ステップ１で与えた分子集合体の初期座標データに対して分子集合体計算を実行し、少なくとも１つ以上のスナップショットを出力するステップ（ステップ２）と、
（ｃ）ステップ２で出力したスナップショットから少なくとも１つを選択し、選択したスナップショットに対してランダムな初期位置ｒ_０に溶質分子を挿入するステップ（ステップ３）と、
（ｄ）溶質分子が初期位置ｒ_０にあるときのポテンシャルエネルギー差Ｖ（ｒ_０）を計算するステップ（ステップ４）と、
（ｅ）ｒ_０、Ｖ（ｒ_０）を電子計算機に保存するステップ（ステップ５）と、
（ｆ）溶質分子を元の位置ｒ_ｉ−１から次の位置ｒ_ｉに変化させるステップ（ステップ６）と、
（ｇ）溶質分子が位置ｒ_ｉにあるときのポテンシャルエネルギー差Ｖ（ｒ_ｉ）を計算するステップ（ステップ７）と、
（ｈ）Ｖ（ｒ_ｉ） とＶ（ｒ_ｉ−１）の差であるΔＶを計算するステップ（ステップ８）と、
（ｉ） ステップ６で変化させた溶質分子の位置ｒ_ｉを、判定基準を用いて採択または棄却を判定するステップ（ステップ９）と、
（ｊ）ｒ_ｉ、Ｖ（ｒ_ｉ）を電子計算機に保存するステップ（ステップ１０）と、
（ｋ）化学ポテンシャルを計算するステップ（ステップ１１）と、
（ｌ）化学ポテンシャルの値を出力するステップ（ステップ１２）とを有し、
分子集合体計算の終了後にステップ３を行い、ステップ６〜１０を規定回数に到達するまで繰り返し行うことを特徴とする計算方法。

（２）ステップ（９）における判定基準としてメトロポリスのアルゴリズムを用い、温度Ｔ_１における化学ポテンシャルの値を、ステップ１０で保存されたポテンシャルエネルギー差Ｖ（ｒ_ｉ）を用いて、数式（１）から計算する（１）に記載の計算方法。

ここで、μは化学ポテンシャル、ｋ_Ｂはボルツマン定数、Ｔ_２はメトロポリスのアルゴリズムを実行するときの温度、Ｎはサンプリング数、ｉはステップ番号を表す。

（３）Ｔ_２の値を４００Ｋより大きく、８０００Ｋ未満とすることを特徴とする（２）に記載の計算方法。

（４）ステップ３でスナップショットを２個以上選択し、統計処理することによって化学ポテンシャルの値を計算することを特徴とする（１）に記載の計算方法。

（５）ステップ（３）でスナップショットを２個以上選択し、さらに、ステップ（９）における判定基準としてメトロポリスのアルゴリズムを用い、温度Ｔ_１における化学ポテンシャルの値を、ステップ（１０）で保存されたポテンシャルエネルギー差Ｖ（ｒ_ｉ）を用いて、数式（２）を用いて計算することを特徴とする（２）または（３）のいずれかに記載の計算方法。

ここで、Ｍは選択したスナップショット数、ｍはスナップショット番号、Ｖ_ｍｉはｍ番目のスナップショットでのステップ番号ｉにおけるポテンシャルエネルギー差を表す。

（６）溶質分子が存在しうる位置を、スナップショット内の自由体積に限定し、ステップ１１の後で、数式（３）で表わされる補正項μ^ｃｏｒｒの値を化学ポテンシャルに加えることを特徴とする（２）、（３）または（５）のいずれかに記載の計算方法。

Ｍは選択したスナップショット数、C_fmはｍ番目のスナップショット内の特定の範囲が占める体積、C_smはｍ番目のスナップショットの基本セルが占める体積を表わす。

（７）自由体積は、ある半径ａを有するプローブによりスナップショット内を走査し、スナップショット内の分子集合体とプローブが接触しない範囲を求めることによって算出することを特徴とする（６）に記載の計算方法。

（８）半径ａの値を０．５オングストローム以上３オングストローム以下とすることを特徴とする（７）に記載の計算方法。

本発明は、分子集合体中で溶質分子の存在しうる位置を、スナップショット内の自由体積に限定しない場合は、図２に示す計算フローに従う。また、スナップショット内の自由体積に限定する場合は、図３および４に示す計算フローに従う。

まず、図２の計算フローに従う場合の計算手順について詳細を示す。

ステップ１においては、電子計算機に読み取り可能な形態で分子集合体の初期座標データを格納する。格納には、一般に入手可能な分子構造式作成ソフトウェアやテキストエディタなどが利用できる。

ステップ２においては、ステップ１において電子計算機に格納された分子集合体に対して分子集合体計算を実行し、設定したステップ数ごとにスナップショットを出力する。分子集合体計算終了後ステップ３に進む。

ステップ３においては、ステップ２で出力したスナップショットから少なくとも１つを選択し、選択したスナップショットに対して溶質分子をランダムな初期位置ｒ_０に挿入する。このとき、スナップショットを２個以上選択し、化学ポテンシャルの値を数式（２）から計算することが好ましい。スナップショットを複数個選択することで、統計的な誤差を少なくし、計算精度を向上させることができる。

ステップ４においては、ステップ３で選択したスナップショットに対して、溶質分子を分子集合体中のランダムな初期位置ｒ_０に挿入した後、数式（５）を用いて溶質分子が初期位置ｒ_０にあるときのポテンシャルエネルギー差Ｖ（ｒ_０）を計算する。

ステップ６においては、溶質分子を元の位置ｒ_ｉ−１とは異なる位置ｒ_ｉに変化させる。

ステップ８においては、ステップ６で変化させた溶質分子の位置ｒ_ｉでのポテンシャルエネルギー差Ｖ（ｒ_ｉ）と、変化させる前の元の位置ｒ_ｉ−１でのポテンシャルエネルギー差Ｖ（ｒ_ｉ−１）を用いてΔＶ = Ｖ（ｒ_ｉ） - Ｖ（ｒ_ｉ−１）を計算する。

ステップ９においては、ステップ６で変化させた溶質分子の位置ｒ_ｉを、判定基準を用いて採択または棄却を判定する。判定基準としては、例えば、ポテンシャルエネルギー差や溶質分子の位置などにある閾値を設け、その値に基づいてサンプリングするといった任意の方法を用いることができる。

しかし、これらの方法は必ずしも統計力学的に意味のあるアンサンブルを発生しないため、カノニカルアンサンブルを発生することが統計力学的に証明されている方法としてメトロポリスのアルゴリズムが好ましく用いられる。この方法は、数式（６）に示す重み関数を用いて、溶質分子が分子集合対中に安定に存在する位置を効率的に発生することができる方法であるため、本発明の目的のため好適に用いることができる。

ここで、Ｔ_２はメトロポリスのアルゴリズムを実行するときの温度である。さらに好ましくは、Ｔ_２を最適な温度範囲に設定し、メトロポリスのアルゴリズムを実行するのが良い。

メトロポリスのアルゴリズムの実行手順を以下に示す。まず、ステップ８で計算したΔＶの値を用いて重み関数Ｐの値を計算する。Ｐの値が１以上であれば位置ｒ_ｉを棄却しステップ６へと進む。一方、Ｐの値が１未満であれば、０から１の一様乱数Ｘを発生させ、Ｘの値と比較する。Ｘ＜Ｐであれば、位置ｒ_ｉを採択しステップ１０へと進む。Ｘ＞Ｐであれば、位置ｒ_ｉを棄却しステップ６へと進む。

数式（６）で示した重み関数に用いる温度Ｔ_２の値によって、位置ｒ_ｉの採択確率が変化する。Ｔ_２の値は４００Ｋより大きく８０００Ｋ未満の範囲に設定するのが好ましく、さらに好ましくは５５０Ｋより大きく４５００Ｋ未満の範囲に設定するのが良い。Ｔ_２が４００Ｋ以下の場合、小さなΔＶの値をもつ位置を主に採択するため、とりうるＶ（ｒ_ｉ）の値の範囲が狭くなる。また、準安定な配置を取ると、他の位置に移行しにくくなるため、溶質分子が分子集合体中に安定に存在する位置を効率的に発生させることができず、サンプリングの効率が悪くなる。

一方、Ｔ_２が高温の場合、様々なΔＶの値をもつ位置を採択することでき、溶質分子が分子集合体中に安定に存在する位置を高頻度で発生できるようになる。しかし、Ｔ_２を８０００Ｋ以上の高温に設定すると、分子集合体と溶質分子とが反発しあう位置を高頻度で発生してしまうため、サンプリングの効率が悪くなる。

ステップ１１においては、ステップ１０で保存されたポテンシャルエネルギー差Ｖ（ｒ_ｉ）を用いて化学ポテンシャルを計算する。一般には、判断基準に応じた重み関数の効果を取り入れた数式を用いて化学ポテンシャルを計算することができるが、ステップ９においてメトロポリスのアルゴリズムを用いた場合は、数式（１）に従って化学ポテンシャルを計算する。

次に、図３および４の計算フローに従う場合の計算手順について詳細を示す。ステップ１、２においては、図２のステップ１、２と同じ計算手順を用いる。ここで、溶質分子の位置が分子集合体に近接しすぎている配置は、立体反発相互作用のため非常に高いエネルギーを持ち、棄却される確率が高い。そこで本発明者らは、このような高エネルギーの配置を予め除外しておくことによってサンプリングの効率を向上できることに気づき、以下に示すような自由体積の考え方を導入するに至った。

すなわち、ステップ３においては、ステップ２で出力したスナップショットから少なくとも１つを選択し、スナップショット内に存在する自由体積を計算した後、自由体積の位置データを電子計算機に保存する。この自由体積の位置データを用いて、ステップ４、ステップ７において溶質分子が存在しうる位置を自由体積内に限定するのである。自由体積の計算方法しては、例えば、スナップショット全体の体積と、分子集合体のファンデルワールス体積との差から計算する方法や、公知の文献に記載の方法、あるいは、Ａｃｃｅｌｒｙｓ社製のＣｅｒｉｕｓ^２やＭａｔｅｒｉａｌｓ Ｓｔｕｄｉｏのような一般的な分子設計ソフトウェアを用いて計算する方法を用いることができる。

Ｃｅｒｉｕｓ^２では、ある半径ａを有するプローブによりスナップショット内を走査し、スナップショット内の分子集合体とプローブが接触しない範囲を求めることによって自由体積を算出する。この方法は、プローブ半径ａの値を調整することによって、溶質分子が分子集合体と近距離にある不安定な領域を自由体積から除外できるという利点をもつ。そのため、本発明における自由体積の計算方法として好適に用いることができる。

プローブ半径ａは、０．５オングストローム以上３．０オングストローム以下に設定するのが好ましく、さらに好ましくは１．０オングストローム以上２．５オングストローム以下に設定するのが良い。ａの値が０．５オングストローム未満の場合、溶質分子が分子集合体に近接しすぎる領域まで自由体積に含まれるため、サンプリングの効率が悪くなる。ａを３．０オングストロームより大きな値に設定すると、溶質分子と分子集合体が互いに近接して強い分子間引力を持つ領域を自由体積から除外してしまうため、溶質分子が分子集合体中に安定に存在する位置を十分にサンプリングできなくなり、計算精度が悪くなる。

ステップ４においては、ステップ３で計算した自由体積内のランダムな初期位置ｒ_０に溶質分子を挿入する。このとき、ステップ３において、スナップショットを２個以上選択し、それぞれのスナップショットに対して自由体積を計算し、化学ポテンシャルを数式（２）から計算した後、補正項μ^ｃｏｒｒを加えることが好ましい。スナップショットを複数個選択することで、統計的な誤差を少なくし、計算精度を向上させることができる。

ステップ５においては、図２のステップ４と同じ計算手順を用いる。

ステップ７においては、溶質分子を元の位置ｒ_ｉ−１とは異なる位置ｒ_ｉに変化させる。このとき、溶質分子の変化させる位置ｒ_ｉは、スナップショットの自由体積内に限定する。

ステップ８〜１２においては、図２のステップ７〜１１と同じ計算手順を用いる。

ステップ１３においては、ステップ１２で計算した化学ポテンシャルの値に数式（３）で示される補正項μ^ｃｏｒｒを加える。

本発明の具体的実施態様を以下に実施例をもって述べるが、本発明はこれに限定されるものではない。

実施例1
以下のステップ１〜１２の手順により、ポリエチレンテレフタレート（ＰＥＴ）中の水分子の化学ポテンシャルを計算した。

ステップ１：分子集合体の初期座標データを入力
まず、Ａｃｃｅｌｒｙｓ社製の分子設計システムＣｅｒｉｕｓ^２上で３Ｄ−Ｓｋｅｔｃｈｅｒを用いて構造式（１）に示すモノマー単位を作成した。次いで、量子化学計算ソフトウェア・ガウシアン９８を実行し、モノマーの電荷を決定した。作成したモノマーモデルに対してＢ３ＬＹＰ／６−３１Ｇ（ｄ、ｐ）レベルの構造最適化を行った後、ＣｈｅｌｐＧの方法を用いて電荷を決定した。計算に必要なパラメータについては、全てソフトウェアの初期設定値を用いた。以上の量子化学計算は、クロック周波数３.１９ＧＨｚのＰｅｎｔｉｕｍ（登録商標）４を搭載した電子計算機上で実行した。その後、炭素原子に結合した水素原子の電荷、質量、ファンデルワールス半径を炭素原子に組み入れ、水素原子を削除することで融合原子モデル（ｕｎｉｔｅｄ ａｔｏｍ ｍｏｄｅｌ）を作成した。

次に、Ｃｅｒｉｕｓ^２−Ｐｏｌｙｍｅｒ Ｂｕｉｌｄｅｒを用い、分子量約３０００のポリマー分子鎖を４本作成した。その後、Ｃｅｒｉｕｓ^２−Ａｍｏｒｐｈｏｕｓ Ｂｕｉｌｄｅｒを用い、周期境界条件を課して低密度（０．２ｇ／ｃｍ^３）の非晶構造を作成した。
最後に、作成した非晶構造の単位セルに対して、Ｃｅｒｉｕｓ^２−ＯＦＦ Ｍｉｎｉｍｉｚｅｒを用い、セルを固定した状態で構造最適化計算を行い、初期座標データとして電子計算機に格納した。

ステップ２：分子動力学計算を実行し、スナップショットを出力
ステップ１で作成したポリマーモデルについて定温定圧の分子動力学計算を実行した。分子動力学計算に用いた力場パラメータは、ドライディング（ＤＲＥＩＤＩＮＧ） ［ エス・エル・メイオー （Ｓ. Ｌ. Ｍａｙｏ） 、ビー・ディー・オラフソン （Ｂ. Ｄ. Ｏｌａｆｓｏｎ） 、ダブリュー・エー・ゴダード３世 （Ｗ. Ａ. Ｇｏｄｄａｒｄ III） 、ジャーナル オブ フィジカル ケミストリー （Ｊ. Ｐｈｙｓ. Ｃｈｅｍ.） 、１９９０、９４、ｐ.８８９７.］ 、アンバー （ＡＭＢＥＲ） [ダブリュー・ディー・コーネル （Ｗ. Ｄ. Ｃｏｒｎｅｌｌ） 、ピー・シープラック （Ｐ. Ｃｉｅｐｌａｋ） 、シー・アイ・ベイリー （Ｃ. Ｉ. Ｂａｙｌｙ） 、アイ・アール・グルード （Ｉ. Ｒ. Ｇｏｕｌｄ） 、ケー・エム・メルツ・ジュニア （Ｋ. Ｍ. Ｍｅｒｚ ｊｒ） 、ディー・エム・ファーガソン （Ｄ. Ｍ. Ｆｅｒｇｕｓｏｎ） 、ディー・シー・スペルマイアー （Ｄ. Ｃ. Ｓｐｅｌｌｍｅｙｅｒ） 、ティー・フォックス （Ｔ. Ｆｏｘ） 、ジェイ・ダブリュー・コールドウェル （Ｊ. Ｗ. Ｃａｌｄｗｅｌｌ）、ピー・エー・コールマン （Ｐ. Ａ. Ｋｏｌｌｍａｎ） 、ジャーナル オブ アメリカン ケミカル ソサイアティー （Ｊ. Ａｍ. Ｃｈｅｍ. Ｓｏｃ.）、１９９５、１１７、ｐ.５１７９]、桑島等のパラメータ [桑島聖、野間祐之、逢坂俊郎、“非晶高分子のＭＤ計算用力場パラメータについてｉｉ：ポリオレフィン”、第４回計算化学シンポジウム&ＣＣＳ展 講演要旨集、日本化学プログラム交換機構（１９９４）]である。

温度は能勢・フーバー （Ｈｏｏｖｅｒ）法[能勢修一 （Ｓ. Ｎｏｓｅ） 、エム・エル・クライン （Ｍ. Ｌ. Ｋｌｅｉｎ） 、モレキュラーフィジクス （Ｍｏｌ. Ｐｈｙｓ.）、１９８３、５０、ｐ.１０５５／能勢修一 （Ｓ. Ｎｏｓｅ）、モレキュラーフィジクス （Ｍｏｌ. Ｐｈｙｓ.）、１９８４、５２、ｐ.２５５／ 能勢修一 （Ｓ. Ｎｏｓｅ） 、ジャーナル オブ ケミカル フィジクス （Ｊ. Ｐｈｙｓ. Ｃｈｅｍ.）、１９８４、８１、ｐ.５１１.]により２５℃に制御し、圧力は斜交系セルを用いるパリネロ・ラーマン（Ｐａｒｒｉｎｅｌｌｏ−Ｒａｈｍａｎ）法 [エム・パリネロ （Ｍ. Ｐａｒｒｉｎｅｌｌｏ）、エー・ラーマン （Ａ. Ｒａｈｍａｎ）、ジャーナル オブ アプライド フィジクス （Ｊ. Ａｐｐｌ. Ｐｈｙｓ.）、１９８１、５２、ｐ.７１８２]に従い１気圧に制御した。クーロン相互作用はエワルド（Ｅｗａｌｄ）の方法[エム・ピー・トシ （Ｍ. Ｐ. Ｔｏｓｉ） 著、エフ・ザイツ （Ｆ. Ｓｅｉｔｚ） 、ディー・ターンブル （Ｄ. Ｔｕｒｎｂｕｌｌ）編、ソリッドステイト・フィジクス・アドバンシィイズ・イン・リサーチ・アンド・アプリケーションズ （Ｓｏｌｉｄ Ｓｔ. Ｐｈｙｓ. Ａｄｖａｎｃｅｓ ｉｎ Ｒｅｓｅａｒｃｈ ａｎｄ Ａｐｐｌｉｃａｔｉｏｎ）、１９６４、１６、ｐ.１−１２０.]を用いた。数値積分には異なる時定数の力ごとに大きさの異なる時間刻みを設定することのできる多重時間刻み数値積分ＲＥＳＰＡ法 [ジー・ジェイ・マルチナ （Ｇ. Ｊ. Ｍａｒｔｙｎａ）、モルキュラー フィジクス （Ｍｏｌ. Ｐｈｙｓ.）、１９９６、８７、ｐ.１１１７.]を用いた。時間刻みの設定を表１に示す。

シミュレーション時間は２００ｐｓとし、クロック周波数３．１９ＧＨｚのＰｅｎｔｉｕｍ（登録商標）４を搭載した電子計算機上で自作のプログラムを用いて実行した。その後、２００ｐｓ時のスナップショットを出力した。

次に、定温定圧の分子動力学計算から得られたスナップショットを初期構造として、定温・定積の分子動力学計算を行った。温度は能勢・フーバー（Ｈｏｏｖｅｒ）法を用いて２５℃に制御した。力場パラメータは定温定圧の分子動力学計算に用いたものと同様のパラメータを使用した。シミュレーション時間は５００ｐｓとし、１００ｐｓ刻みで計５本のスナップショットを出力した。

ステップ３〜１０：サンプリング
ステップ３においては、５本のスナップショットを選択した。その後、Ｃｅｒｉｕｓ^２−Ｓｏｒｐｔｉｏｎを用いて、選択したスナップショットそれぞれに対して、サンプリングを行った。水分子の力場パラメータにはＴＩＰ３Ｐモデル[ダブリュー・エル・ヨルゲンセン （Ｗ. Ｌ. Ｊｏｒｇｅｎｓｅｎ） 、アール・ダブリュー・インペイ （Ｒ. Ｗ. Ｉｍｐｅｙ） 、ジェイ・チャンドラセカール （Ｊ. Ｃｈａｎｄｒａｓｅｋｈａｒ） 、ジェイ・ディー・マドゥラ （Ｊ. Ｄ. Ｍａｄｕｒａ） 、エム・エル・クライン （Ｍ. Ｌ. Ｋｌｅｉｎ） 、ジャーナル オブ ケミカル フィジクス （Ｊ. Ｃｈｅｍ. Ｐｈｙｓ.）、１９８３、７９、ｐ.９２６.）を用いた。実行手順を以下に示す。

最初に、Ｖ（ｒ_ｉ）の値を１ステップごとに保存するように設定した。ステップ３で選択したスナップショットに対して、水分子１個を挿入し、サンプリングを行った。サンプリング数Ｎは２００００００とし、メトロポリスのアルゴリズムを実行する温度Ｔ_２は５００Ｋに設定した。

上記に記載の操作を各スナップショットに対してそれぞれ行った。

ステップ１１：化学ポテンシャルの計算
本実施例ではスナップショットを５本選択しているため、数式（２）においてＭ＝５とし、化学ポテンシャルを計算した。

ステップ１２：化学ポテンシャルの出力
本発明は化学ポテンシャルを高速に計算する方法を提供するものである。計算時間を比較するために、設定したサンプリング数の範囲内で、化学ポテンシャルの値がどれだけ早く収束しているかを確認することで評価することとした。化学ポテンシャルの収束は、数式（７）に示すように、あるステップ番号ＸからＹまでの区間の化学ポテンシャルの平均値を計算し、その変化を観察することで確認した。なお、本発明では、計算開始点Ｘを８０００００ステップとし、Ｙを８１００００、８２００００と、１００００ステップずつ変えて計算を行った。

図５に化学ポテンシャルの収束性の計算結果を示す。縦軸は化学ポテンシャルの値、横軸にステップ番号をとった。

Ｔ_２ ＝ ５００Ｋでサンプリングを実行した場合、化学ポテンシャルは、１８０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。

ここで、収束性の判断基準を以下に述べる。図６、図７に化学ポテンシャルの収束性を計算したものを示す。

まず、あるステップ数で化学ポテンシャルが一定の値近くで揺らぎ始めたら、その点を収束開始点とする。収束開始点からサンプリング終了時まで収束開始点の値±０．５（ｋＪ／ｍｏｌ）未満の値であれば収束しているとみなす。また、収束開始点以降に、収束開始点からサンプリング終了時までに±０．５（ｋＪ／ｍｏｌ）以上の値であれば収束していないとみなす。

図６では、化学ポテンシャルが、サンプリング終了時まで一定の値に収束していないことが分かる。一方、図７では、１００００００ステップ目から一定の値に収束し始め、計算終了時まで収束開始点の値±０．５（ｋＪ／ｍｏｌ）未満の値を保持し続けている。よって、図７では１００００００ステップ目を収束点として一定の値に収束したと判断する。

上記の判断基準に従うと、Ｔ_２ ＝ ５００Ｋでサンプリングを実行した場合では、化学ポテンシャルは収束しており、収束点は１８０００００ステップと評価できる。
以下実施例２〜１１および比較例１においても、上記の判断基準に従い、化学ポテンシャルの収束性を評価した。

実施例２
Ｔ_２を６００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１３０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１３０００００ステップであった。

実施例3
Ｔ_２を８００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１３０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１３０００００ステップであった。

実施例４
Ｔ_２を１０００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１２０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１２０００００ステップであった。

実施例５
Ｔ_２を２０００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１２０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１２０００００ステップであった。

実施例６
Ｔ_２を３０００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１２０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１２０００００ステップであった。

実施例７
Ｔ_２を４０００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１３０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１３０００００ステップであった。

実施例８
Ｔ_２を５０００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１９０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１９０００００ステップであった。

実施例９
Ｔ_２を７０００Ｋとした他は、実施例１と同じスナップショットを用い、同様の手順でサンプリングを実行した。化学ポテンシャルは、１８０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１８０００００ステップであった。

実施例１０
実施例１で使用したスナップショットを用い、Ｃｅｒｉｕｓ^２−Ｆｒｅｅ Ｖｏｌｕｍｅを用いてプローブ半径ａを１．０オングストロームに設定し、スナップショット内の自由体積を計算した。その後、自由体積の値と位置データを電子計算機に格納した。次いで、自由体積の位置データを読み込み、溶質分子の存在しうる領域を自由体積内に限定した後、Ｃｅｒｉｕｓ^２−Ｓｏｒｐｔｉｏｎを用いて水分子１個を挿入し、サンプリングを行った。サンプリング数Ｎは２００００００とし、メトロポリスのアルゴリズムの実行温度Ｔ_２を６００Ｋに設定した。

本実施例ではスナップショットを５本選択しているため、数式（８）においてＭ＝５とし、化学ポテンシャルの収束性を評価した。

化学ポテンシャルは、１００００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１００００００ステップであった。

実施例１１
Ｔ_２を５０００Ｋとした他は、実施例１と同じスナップショットを用い、実施例１０と同様の手順でサンプリングを実行した。化学ポテンシャルは、１００００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１００００００ステップであった。

比較例１
Ｔ_２を９９９９９９９９Ｋとした他は、実施例１と同じスナップショットを用い、実施例１と同様の手順でサンプリングを実行した。

ここで、数式（２）においてＴ_２→∞ の極限では exp[-V/kT₂] →１となるため、数式（２）は数式（９）に帰着する。そのため、Ｔ_２を９９９９９９９９Ｋとした場合、実質的に粒子挿入法と同等の計算を実行できる。

上記の手順で粒子挿入法を実行した場合、化学ポテンシャルは、サンプリング終了時まで収束しなかった。

比較例２
Ｔ_２を９９９９９９９９Ｋとした他は、実施例１と同じスナップショットを用い、実施例１０と同様の手順でサンプリングを実行した。化学ポテンシャルは、比較例１と同様に数式（９）から計算した。その後、化学ポテンシャルに補正項μ^ｃｏｒｒを加えた。

化学ポテンシャルは、１８０００００ステップ目から一定の値に収束し始め、サンプリング終了時までその値近くを保持し続けていた。よって、化学ポテンシャルは収束しており、収束点は１８０００００ステップであった。

実施例１〜１１および比較例１、２の化学ポテンシャルの収束性を表２に示す。表のＡは１５０００００ステップ以内に収束点があるもの、Ｂは１５０００００〜２００００００ステップ以内に収束点があるもの、Ｃは計算終了時まで収束しなかったものを表す。

この表から、粒子挿入法を用いた比較例１では、化学ポテンシャルの収束に少なくとも２００００００ステップのサンプリング数が必要であることがわかる。これに対して、溶質分子の位置をメトロポリスのアルゴリズムに従うように変化させた場合、実施例１、８、９では、１８０００００〜１９０００００ステップ目で化学ポテンシャルの値が収束していた。さらに、実施例２〜７では、１２０００００〜１３０００００ステップ目で化学ポテンシャルが収束していた。

このことから、分子集合体中に挿入した溶質分子の位置を、メトロポリスのアルゴリズムに従うように変化させた場合、Ｔ_２の値をおよそ５００Ｋから７０００Ｋの範囲に設定することで、最大１９０００００以下のサンプリング数で化学ポテンシャルを計算でき、また、およそ６００Ｋから４０００Ｋの範囲に設定することで、最大１３０００００以下のサンプリング数で化学ポテンシャルを計算できた。

また、溶質分子の存在しうる位置を自由体積に限定し、さらに溶質分子の位置をメトロポリスのアルゴリズムに従うように変化させた場合、実施例１０および１１では、１００００００ステップ目で化学ポテンシャルが収束していた。ここで、実施例４〜６と実施例１０、１１を比較すると、実施例１０および１１は、実施例４〜６よりも少ないステップ数で化学ポテンシャルが収束しており、より速く化学ポテンシャルを計算できた。

粒子挿入法の計算フローを説明するための図である。 本発明の方法の計算フローを説明するための図である。 本発明の方法の計算フローを説明するための図である。 本発明の方法の計算フローを説明するための図である。 化学ポテンシャルの収束を説明するための図である。 化学ポテンシャルの収束を説明するための図である。 化学ポテンシャルの収束を説明するための図である。

Claims (8)

1. 分子集合体と溶質との親和性を求めるための分子シミュレーションにおける、自由エネルギーの計算方法であって、
   （ａ）分子集合体の初期座標データを電子計算機に格納するステップ（ステップ１）と、
   （ｂ）ステップ１で与えた分子集合体の初期座標データに対して分子集合体計算を実行し、少なくとも１つ以上のスナップショットを出力するステップ（ステップ２）と、
   （ｃ）ステップ２で出力したスナップショットから少なくとも１つを選択し、選択したスナップショットに対してランダムな初期位置ｒ_０に溶質分子を挿入するステップ（ステップ３）と、
   （ｄ）溶質分子が初期位置ｒ_０にあるときのポテンシャルエネルギー差Ｖ（ｒ_０）を計算するステップ（ステップ４）と、
   （ｅ）ｒ_０、Ｖ（ｒ_０）を電子計算機に保存するステップ（ステップ５）と、
   （ｆ）溶質分子を元の位置ｒ_ｉ−１から次の位置ｒ_ｉに変化させるステップ（ステップ６）と、
   （ｇ）溶質分子が位置ｒ_ｉにあるときのポテンシャルエネルギー差Ｖ（ｒ_ｉ）を計算するステップ（ステップ７）と、
   （ｈ）Ｖ（ｒ_ｉ） とＶ（ｒ_ｉ−１）の差であるΔＶを計算するステップ（ステップ８）と、
   （ｉ） ステップ６で変化させた溶質分子の位置ｒ_ｉを、判定基準を用いて採択または棄却を判定するステップ（ステップ９）と、
   （ｊ）ｒ_ｉ、Ｖ（ｒ_ｉ）を電子計算機に保存するステップ（ステップ１０）と、
   （ｋ）化学ポテンシャルを計算するステップ（ステップ１１）と、
   （ｌ）化学ポテンシャルの値を出力するステップ（ステップ１２）とを有し、
   分子集合体計算の終了後にステップ３を行い、ステップ６〜１０を規定回数に到達するまで繰り返し行うことを特徴とする計算方法。
2. ステップ（９）における判定基準としてメトロポリスのアルゴリズムを用い、温度Ｔ_１における化学ポテンシャルの値を、ステップ１０で保存されたポテンシャルエネルギー差Ｖ（ｒ_ｉ）を用いて、数式（１）から計算する請求項１に記載の計算方法。
   

   

   ここで、μは化学ポテンシャル、ｋ_Ｂはボルツマン定数、Ｔ_２はメトロポリスのアルゴリズムを実行するときの温度、Ｎはサンプリング数、ｉはステップ番号を表す。
3. Ｔ_２の値を４００Ｋより大きく、８０００Ｋ未満とすることを特徴とする請求項２に記載の計算方法。
4. ステップ３でスナップショットを２個以上選択し、統計処理することによって化学ポテンシャルの値を計算することを特徴とする請求項１に記載の計算方法。
5. ステップ（３）でスナップショットを２個以上選択し、さらに、ステップ（９）における判定基準としてメトロポリスのアルゴリズムを用い、温度Ｔ_１における化学ポテンシャルの値を、ステップ（１０）で保存されたポテンシャルエネルギー差Ｖ（ｒ_ｉ）を用いて、数式（２）を用いて計算することを特徴とする請求項２または３のいずれかに記載の計算方法。
   

   

   ここで、Ｍは選択したスナップショット数、ｍはスナップショット番号、Ｖ_ｍｉはｍ番目のスナップショットでのステップ番号ｉにおけるポテンシャルエネルギー差を表す。
6. 溶質分子が存在しうる位置を、スナップショット内の自由体積に限定し、ステップ１１の後で、数式（３）で表わされる補正項μ^ｃｏｒｒの値を化学ポテンシャルに加えることを特徴とする請求項２、３または５のいずれかに記載の計算方法。
   

   

   Ｍは選択したスナップショット数、C_fmはｍ番目のスナップショット内の自由体積が占める体積、C_smはｍ番目のスナップショットの基本セルが占める体積を表わす。
7. 自由体積は、ある半径ａを有するプローブによりスナップショット内を走査し、スナップショット内の分子集合体とプローブが接触しない範囲を求めることによって算出することを特徴とする請求項６に記載の計算方法。
8. 半径ａの値を０．５オングストローム以上３オングストローム以下とすることを特徴とする請求項７に記載の計算方法。

Images