JP6368637B2 - 分子間ポテンシャルの算出装置、分子間ポテンシャルの算出方法、及びコンピュータプログラム - Google Patents

Description

本発明は、分子間ポテンシャルの算出装置、分子間ポテンシャルの算出方法、及びコンピュータプログラムに関する。

分子の特性は、分子間に働く相互作用（ポテンシャル）によって定まるため、分子間ポテンシャルを知ることが重要である。アルゴンなどの単原子分子の場合には分子間ポテンシャルが既知であるが、高分子などの多原子分子の場合には分子間ポテンシャルが未知である。分子間ポテンシャルを導出する方法として、例えば非特許文献１には、積分方程式を用いて液体モデルのポテンシャルを算出することについての記載がある。分子動力学計算に基づきコンピュータでシミュレーションを行い、その結果と非特許文献１に記載の方法を用いれば、液体モデルの未知の分子間ポテンシャルを算出可能となる。

具体的に、積分方程式を用いたポテンシャルの算出は、一般的に、次の（１）〜（５）のステップを実施することで実現される。
（１）液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出
（２）算出した座標位置に基づき動径分布関数を取得
（３）動径分布関数に基づき全相関関数を算出してフーリエ変換
（４）フーリエ変換した全相関関数に基づき直接相関関数を算出
（５）動径分布関数、全相関関数及び直接相関関数に基づき分子間ポテンシャルを算出

"Y. Rosenfeld and G. Kahl"著，「The inverse problem for simple classical liquids: a density functional approach」, J. Phys.: Condens. Matter 9, 1997, L89-L98. Printed in the UK. IOP Publishing Ltd.

発明者は、液体分子モデルを単原子分子とし、分子間ポテンシャルにレナードジョーンズポテンシャル（以下、ＬＪポテンシャルとも表記する）を与えて、上記ステップ（１）において液体分子モデルの分子動力学計算を行い、上記ステップ（５）においてＨＮＣ（Hypernetted-chain）近似と呼ばれる積分方程式の近似式を用いて分子間ポテンシャルを算出した。しかしながら、本来であれば、算出したポテンシャルの形状がＬＪポテンシャルの形状に一致又は近似していないといけないところ、図４Ａに示すように、算出したポテンシャル（図中ではＨＮＣと表記）の形状がＬＪポテンシャルの形状から大きくずれてしまい、ポテンシャルの形状を適切に導出できていないことが判明した。この問題は、ＨＮＣ近似に限らず、他の近似式でも生じ得る。

本発明は、このような課題に着目してなされたものであって、その目的は、算出精度を向上させた分子間ポテンシャルの算出装置、分子間ポテンシャルの算出方法、及びコンピュータプログラムを提供することである。

本発明は、上記目的を達成するために、次のような手段を講じている。

すなわち、本発明の分子間ポテンシャルの算出装置は、液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出する分子動力学算出部と、前記分子動力学算出部が算出した座標位置に基づき動径分布関数を算出する動径分布関数算出部と、前記動径分布関数算出部が算出した前記動径分布関数に基づき全相関関数を算出し、算出した全相関関数をフーリエ変換により波数空間の全相関関数に変換する波数空間全相関関数算出部と、前記波数空間の全相関関数のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数に置き換えることで、波数空間の全相関関数を補正する補正部と、補正後の全相関関数に基づき直接相関関数を算出する直接相関関数算出部と、動径分布関数、補正後の全相関関数及び直接相関関数に基づき分子間ポテンシャルを算出する分子間ポテンシャル算出部と、を備える。

本発明の分子間ポテンシャルの算出方法は、コンピュータが実行する方法であって、（ａ）液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出するステップと、（ｂ）算出した座標位置に基づき動径分布関数を算出するステップと、（ｃ）算出した前記動径分布関数に基づき全相関関数を算出し、算出した全相関関数をフーリエ変換により波数空間の全相関関数に変換するステップと、（ｄ）波数空間の全相関関数のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数に置き換えることで、波数空間の全相関関数を補正するステップと、（ｅ）補正後の全相関関数に基づき直接相関関数を算出するステップと、（ｆ）動径分布関数、補正後の全相関関数及び直接相関関数に基づき分子間ポテンシャルを算出するステップと、を含む。

本発明のコンピュータプログラムは、上記方法をコンピュータに実行させる。

このように、波数空間において、全相関関数のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数に置き換えることで、波数空間の全相関関数を補正しているので、波数０付近にある値のバラツキを除去でき、ポテンシャル形状の算出精度を向上させることが可能となる。

本発明の第１実施形態の分子間ポテンシャルの算出装置を模式的に示すブロック図。 本発明の第２実施形態の分子間ポテンシャルの算出装置を模式的に示すブロック図。 第２実施形態の算出装置が実行する算出処理ルーチンを示すフローチャート。 従来方法により算出したポテンシャルと、ＬＪポテンシャルを示す図。 本発明の方法で算出したＨＮＣ近似のポテンシャル及びＭＨＮＣのポテンシャルと、ＬＪポテンシャルを示す図。 フーリエ変換した全相関関数を示す図。 フーリエ変換した全相関関数の補正処理に関する説明図。（ａ）漸近領域の特定に関する説明図。（ｂ）置換関数のフィッティング（係数決定）に関する説明図。（ｃ）置換処理に関する説明図。 （ａ）剛体球のポテンシャルを示す図。（ｂ）Ｓｔｉｃｋｙ剛体球のポテンシャルを示す図。 液体モデルの全相関関数と剛体球の全相関関数とのフィッティングに関する説明図。

以下、本発明の一実施形態を、図面を参照して説明する。

［分子間ポテンシャルの算出装置］
本実施形態の分子間ポテンシャルの算出装置は、液体分子モデルの動径分布関数に基づき分子間ポテンシャルを算出する装置である。分子モデルは液体分子モデルであればよい。

図１に示すように、算出装置１は、初期設定部１０と、分子動力学算出部１１と、動径分布関数算出部１２と、波数空間全相関関数算出部１３と、補正部１４と、直接相関関数算出部１５と、分子間ポテンシャル算出部１６と、を有する。これら各部１０〜１６は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図示しない処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。

図１に示す初期設定部１０は、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、解析対象となる液体分子モデルに関するデータの設定、分子動力学計算に必要な設定などの各種設定を実行し、これら設定値を図示しないメモリに記憶する。

図１に示す分子動力学算出部１１は、予め設定されている液体分子モデル、その分子を構成する原子の結合ポテンシャルなどの各種パラメータを用い、分子モデルの分子動力学計算に基づき平衡状態における各分子モデルの座標位置（重心位置など）を算出する。

図１に示す動径分布関数算出部１２は、分子動力学算出部１１により得られる各モデルの座標位置に基づいて動径分布関数ｇ（ｒ）を算出する。なお、動径分布関数は、ある粒子から距離ｒに他の粒子が存在する確率を表す関数である。

図１に示す波数空間全相関関数算出部１３は、全相関関数算出部１３ａと、フーリエ変換部１３ｂとを有する。全相関関数算出部１３ａは、動径分布関数算出部１２で得られる動径分布関数ｇ（ｒ）に基づき全相関関数ｈ（ｒ）を算出する。フーリエ変換部１３ｂは、全相関関数算出部１３ａが算出した全相関関数ｈ（ｒ）をフーリエ変換により波数空間の全相関関数Ｈ（ｋ）に変換する。なお、全相関関数は、２粒子間の相関を表したものであり、２粒子間に作用する直接的な影響と、２粒子以外の他の粒子の間接的な影響と、が含まれる。２粒子間に作用する直接的な影響は、直接相関関数として表される。

全相関関数ｈ（ｒ）は、式（１）により動径分布関数ｇ（ｒ）に基づき算出する。
ｈ（ｒ）＝ｇ（ｒ）−１ …（１）
ｒは距離を示す。

フーリエ変換した全相関関数Ｈ（ｋ）は、式（２）で表現でき、図５のように波数空間で表現される。

ｋは波数を示す。

図５に示すように、波数空間の全相関関数Ｈ（ｋ）は、波数０の周辺領域において値のバラツキを有する。本発明者は、この値のバラツキがポテンシャルの算出の精度を損なう原因であると考えた。

図１に示す補正部１４は、上記バラツキを是正するために、波数空間の全相関関数のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数に置き換えることで、波数空間の全相関関数を補正する。具体的に、補正部１４は、漸近領域特定部１４ａと、係数決定部１４ｂと、置換部１４ｃと、を有する。

図１に示す漸近領域特定部１４ａは、図６（ａ）に示すように、波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域において、傾きが０に漸近している漸近領域Ａｒを特定する。特定方法は、微分値が０又は０を中心とする予め定めた許容値以内であること、又は、連続する三点間において最大値と最小値の差が予め定めた許容値以内であることを条件として特定する方法が挙げられる。

図１に示す係数決定部１４ｂは、図６（ｂ）に示すように、波数空間において漸近領域Ａｒにおける全相関関数と置換関数との差に対応する評価値が最小となるように、置換関数を構成する係数を決定する。本実施形態では、置換関数として「Ａ・ｋ^２＋Ｂ」を用いている。Ａ，Ｂは係数である。この置換関数は、ｋ＝０において傾きが０となる関数である。この置換関数「Ａ・ｋ^２＋Ｂ」を最小二乗法を用いてフィティングしている。具体的には、漸近領域Ａｒにおける全相関関数と置換関数との差の二乗和を評価値とし、この評価値が最小となる係数Ａ，Ｂを決定する。

図１に示す置換部１４ｃは、波数空間において漸近領域Ａｒ及び漸近領域Ａｒよりも波数が０に近い領域を置換関数「Ａ・ｋ^２＋Ｂ」に置き換えることで、波数空間の全相関関数Ｈ（ｋ）を補正し、補正後の全相関関数Ｈ’（ｋ）を得る。具体的には、図６（ｃ）に示すように、漸近領域Ａｒの最大波数をｋ_１とした場合に、補正後の全相関関数Ｈ’（ｋ）は式（３）で表現される。

本実施形態では、置換関数として「Ａ・ｋ^２＋Ｂ」を用いているが、これに限定されない。例えば、「Ａ・ｋ^３＋Ｂ」、「Ａ・ｋ^４＋Ｂ」、「Ａ・ｅｘｐ（ｋ^２）＋Ｂ」、「Ａ・ｅｘｐ（ｋ^３）＋Ｂ」なども利用可能である。

図１に示す直接相関関数算出部１５は、補正後の全相関関数Ｈ’（ｋ）に基づき直接相関関数ｃ（ｒ）を算出する。まず、式（４）により波数空間において、補正後の全相関関数Ｈ’（ｋ）から直接相関関数Ｃ（ｋ）を算出する。

ρは平均密度である。
式（４）を得られる直接相関関数Ｃ（ｋ）は波数空間の関数であるので、実空間の関数に戻すために、式（５）を用いて逆フーリエ変換を行い、直接相関関数ｃ（ｒ）を算出する。

図１に示す分子間ポテンシャル算出部１６は、動径分布関数ｇ（ｒ）、補正後の全相関関数ｈ’（ｒ）、及び直接相関関数ｃ（ｒ）に基づき分子間ポテンシャルＵ（ｒ）を算出する。本実施形態では、式（６）で表されるＨＮＣ（Hypernetted-chain）近似を用いている。波数空間における補正後の全相関関数Ｈ’（ｋ）を実空間に戻すには、式（７）を用いて逆フーリエ変換すれば、補正後の全相関関数ｈ’（ｒ）を算出できる。全相関関数Ｈ’（ｋ）からｈ’（ｒ）への逆フーリエ変換は、分子間ポテンシャル算出部１６で行ってもよいし、図１に示す逆フーリエ変換部１７で行うようにしてもよい。

Ｔは温度、ｋ_Ｂはボルツマン定数を示す。

図４Ｂは、本実施形態の装置及び方法により算出したポテンシャルを示す。図４Ｂに示すように、分子間ポテンシャル算出部１６が算出したポテンシャルＵ_ＨＮＣ（ｒ）は、図４Ａに比べてＬＪポテンシャルに近い形状であることが分かる。したがって、上記装置１は、従来に比べて分子間ポテンシャルの算出精度が向上していることがわかる。

図２に示すように、分子間ポテンシャル算出部１６においてＨＮＣ近似を用いて算出した分子間ポテンシャルＵ_ＨＮＣ（ｒ）を、剛体球のブリッジ関数で修正するＭＨＮＣ（Modified Hypernetted-chain）近似を用いるようにすれば、更に精度を向上させることができる。図２は、ＭＨＮＣ近似を行う装置のブロック図である。

図２に示すように、算出装置１’の分子間ポテンシャル算出部１６は、式（６）によるＨＮＣ近似を用いて分子間ポテンシャルＵ_ＨＮＣ（ｒ）を算出する。算出装置１’は、更に、ＭＨＮＣ処理部１８を有する。ＭＨＮＣ処理部１８は、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルを算出する。ＭＨＮＣ処理部１８は、パラメータ決定部１８ａと、ブリッジ関数算出部１８ｂと、修正部１８ｃと、を有する。

図２に示すパラメータ決定部１８ａは、波数空間全相関関数算出部１３で得られる波数空間の全相関関数と剛体球モデルの全相関関数の差に対応する評価値が最小となるように、剛体球モデルの全相関関数を構成するパラメータを決定する。本実施形態では、ｓｔｉｃｋｙ剛体球に対するＰＹ近似の解析解を用いて計算を容易にしている。以下に、詳細について説明する。

剛体球モデルのポテンシャルを図７（ａ）に示す。数式では次のように表される。
Ｕ（ｒ）＝∞ ｒ＜σ
Ｕ（ｒ）＝０ ｒ＞σ
σは剛体球径を示す。
ｓｔｉｃｋｙ剛体球のポテンシャルを図７（ｂ）に示す。数式では次のように表される。
Ｕ（ｒ）＝∞ ｒ＜σ
Ｕ（ｒ）＝−ｋ_Ｂ Ｔ ｌｎ［Ｒ／１２τ（Ｒ−σ）］ σ＜ｒ＜Ｒ
Ｕ（ｒ）＝０ ｒ＞Ｒ
Ｒは粘着表面を持った剛体球径、Ｒ−σは無限小、τは温度依存パラメータである。
ｓｔｉｃｋｙ剛体球の解析解として、フーリエ変換した波数空間の直接相関関数Ｃ^ＳＨＳ（ｋ）は式（８）で表される。

ここで、密度ρは液体分子モデルで設定した値である。τについては、σ＜ｒ＜ＲにおいてＵ（ｒ）＝０となるようにＲとσとで決定されるとすれば、剛体球モデルのポテンシャルと同じになるので、式（８）に示すＣ^ＳＨＳ（ｋ）は、剛体球モデルの直接相関関数Ｃ^ＨＳ（ｋ）として取り扱うことができる。式（８）の可変パラメータは、Ｒとσとなる。

パラメータ決定部１８ａは、上記Ｒとσを決定することになる。具体的には、まず、パラメータＲ，σに初期値を与えて、式（８）により剛体球モデルの波数空間における直接相関関数Ｃ^ＨＳ（ｋ）を算出する。次に、式（９）を用いて剛体球モデルの直接相関関数Ｃ^ＨＳ（ｋ）に基づき全相関関数Ｈ^ＨＳ（ｋ）を算出する。

図８の上部に示すように、パラメータＲ，σに初期値を与えて算出した全相関関数Ｈ^ＨＳ（ｋ）_initは、波数空間全相関関数算出部１３で得られる流体モデルの全相関関数Ｈ（ｋ）からずれている。そこで、パラメータ決定部１８ａは、両関数が合うように、最小二乗法によってパラメータＲ，σを決定する。手法は、係数決定部１４ｂが行う最小二乗法と同じであり、パラメータ決定部１８ａは、液体分子モデルの全相関関数Ｈ（ｋ）と剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）の差に対応する評価値（実施形態では、差の二乗和）が最小となるように、剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）を構成するパラメータＲ，σを決定する。

図８の下部は、決定したパラメータＲ，σを用いて算出した全相関関数Ｈ^ＨＳ（ｋ）_fitと、液体分子モデルの全相関関数Ｈ（ｋ）とを示す。両者が合っていることが分かる。

図２に示すブリッジ関数算出部１８ｂは、パラメータ決定部１８ａが決定したパラメータＲ，σに基づき剛体球モデルの全相関関数ｈ^ＨＳ（ｒ）および直接相関関数ｃ^ＨＳ（ｒ）を算出し、剛体球モデルのブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）を式（１０）により算出する。なお、下付添え字ＰＹは、ＰＹ近似によるものであることを示している。
剛体球モデルの直接相関関数ｃ^ＨＳ（ｒ）は、パラメータＲ，σを式（８）に代入して得られる波数空間の直接相関関数Ｃ^ＨＳ（ｋ）を逆フーリエ変換すれば得ることができる。
剛体球モデルの全相関関数ｈ^ＨＳ（ｒ）は、波数空間の直接相関関数Ｃ^ＨＳ（ｋ）を式（９）に代入して全相関関数Ｈ^ＨＳ（ｋ）に変換し、逆フーリエ変換すれば得ることができる。
ｂ_ＰＹ ^ＨＳ（ｒ）＝ｈ^ＨＳ（ｒ）−ｃ^ＨＳ（ｒ）−ｌｎ［１＋ｈ^ＨＳ（ｒ）−ｃ^ＨＳ（ｒ）］ …（１０）

図２に示す修正部１８ｃは、分子間ポテンシャル算出部１６がＨＮＣ近似を用いて算出したポテンシャルＵ_ＨＮＣ（ｒ）を、ブリッジ関数算出部１８ｂが算出したブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）で修正することで、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルＵ_ＭＨＮＣ（ｒ）を算出する。具体的には、式（１１）を用いる。
Ｕ_ＭＨＮＣ（ｒ）＝Ｕ_ＨＮＣ（ｒ）−ｋ_Ｂ Ｔｂ_ＰＹ ^ＨＳ（ｒ） …（１１）

図４Ｂに示すように、修正部１８ｃで修正したポテンシャルＵ_ＭＨＮＣ（ｒ）は、図４Ａに比べてＬＪポテンシャルの形状に近いことが分かる。Ｕ_ＨＮＣ（ｒ）よりもＵ_ＭＨＮＣ（ｒ）の方がより算出精度が近いことが分かる。

［分子間ポテンシャルの算出方法］
図２に示す装置１’を用いた分子間ポテンシャルの算出方法について、図３を用いて説明する。

まず、ステップＳＴ１において、初期設定部１０は、解析対象となる液体分子モデルに関する情報の設定、分子動力学計算に必要な設定等の各種設定を実行し、これら設定値を図示しないメモリに記憶する。

次のステップＳＴ２において、分子動力学算出部１１は、予め設定された各種パラメータを用い、液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出する。

次のステップＳＴ３において、動径分布関数算出部１２は、分子動力学算出部１１が算出した座標位置に基づき動径分布関数ｇ（ｒ）を算出する。

次のステップＳＴ４において、全相関関数算出部１３ａが、動径分布関数算出部１２が算出した動径分布関数ｇ（ｒ）に基づき全相関関数ｈ（ｒ）を算出し、フーリエ変換部１３ｂが、全相関関数ｈ（ｒ）をフーリエ変換により波数空間の全相関関数Ｈ（ｋ）に変換する。

次のステップＳＴ５において、漸近領域特定部１４ａが、波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域において、傾きが０に漸近している漸近領域Ａｒを特定する。

次のステップＳＴ６において、係数決定部１４ｂが、波数空間において漸近領域Ａｒにおける全相関関数Ｈ（ｋ）と置換関数［Ａ・ｋ^２＋Ｂ］との差に対応する評価値が最小となるように、置換関数［Ａ・ｋ^２＋Ｂ］を構成する係数Ａ，Ｂを決定する。

次のステップＳＴ７において、置換部１４ｃが、波数空間において漸近領域Ａｒ及び漸近領域Ａｒよりも波数が０に近い領域を置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数Ｈ（ｋ）をＨ’（ｋ）に補正する。

すなわち、ステップＳＴ５〜７を実行することによって、波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数Ｈ（ｋ）をＨ’（ｋ）に補正することになる。

次のステップＳＴ８において、直接相関関数算出部１５が、補正後の全相関関数Ｈ’（ｋ）に基づき直接相関関数ｃ（ｒ）を算出する。

次のステップＳＴ９において、分子間ポテンシャル算出部１６が、動径分布関数ｇ（ｒ）、補正後の全相関関数ｈ’（ｒ）及び直接相関関数ｃ（ｒ）に基づき分子間ポテンシャルＵ（ｒ）を算出する。本実施形態では、ＨＮＣ近似による分子間ポテンシャルＵ_ＨＮＣ（ｒ）を算出している。図１に示す装置１の動作は以上であるが、図２に示す装置１’の動作は、以下に続く。

次のステップＳＴ１０において、パラメータ決定部１８ａが、波数空間において、ステップＳＴ４で得られる全相関関数Ｈ（ｋ）と剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）の差に対応する評価値が最小になるように、剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）を構成するパラメータＲ，σを決定する。

次のステップＳＴ１１において、ブリッジ関数算出部１８ｂは、パラメータ決定部１８ａが決定したパラメータＲ，σに基づき剛体球モデルの全相関関数ｈ^ＨＳ（ｒ）および直接相関関数ｃ^ＨＳ（ｒ）を算出し、剛体球モデルのブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）を算出する。

次のステップＳＴ１２において、修正部１８ｃは、分子間ポテンシャル算出部１６がＨＮＣ近似を用いて算出したポテンシャルＵ_ＨＮＣ（ｒ）を、ブリッジ関数算出部１８ｂが算出したブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）で修正することで、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルＵ_ＭＨＮＣ（ｒ）を算出する。

以上のように、本実施形態の分子間ポテンシャルの算出装置は、
液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出する分子動力学算出部１１と、
分子動力学算出部１１が算出した座標位置に基づき動径分布関数ｇ（ｒ）を算出する動径分布関数算出部１２と、
動径分布関数算出部１２が算出した動径分布関数ｇ（ｒ）に基づき全相関関数ｈ（ｒ）を算出し、算出した全相関関数ｈ（ｒ）をフーリエ変換により波数空間の全相関関数Ｈ（ｋ）に変換する波数空間全相関関数算出部１３と、
波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数Ｈ（ｋ）をＨ’（ｋ）に補正する補正部１４と、
補正後の全相関関数Ｈ’（ｋ）に基づき直接相関関数ｃ（ｒ）を算出する直接相関関数算出部１５と、
動径分布関数ｇ（ｒ）、補正後の全相関関数Ｈ’（ｋ）及び直接相関関数ｃ（ｒ）に基づき分子間ポテンシャルＵ（ｒ）を算出する分子間ポテンシャル算出部１６と、
を備える。

本実施形態の分子間ポテンシャルの算出方法は、
コンピュータが実行する方法であって、
（ａ）液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出するステップ（ＳＴ２）と、
（ｂ）算出した座標位置に基づき動径分布関数ｇ（ｒ）を算出するステップ（ＳＴ３）と、
（ｃ）算出した動径分布関数ｇ（ｒ）に基づき全相関関数ｈ（ｒ）を算出し、算出した全相関関数ｈ（ｒ）をフーリエ変換により波数空間の全相関関数Ｈ（ｋ）に変換するステップ（ＳＴ４）と、
（ｄ）波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数Ｈ（ｋ）をＨ’（ｋ）に補正するステップ（ＳＴ５〜７）と、
（ｅ）補正後の全相関関数Ｈ’（ｋ）に基づき直接相関関数ｃ（ｒ）を算出するステップ（ＳＴ８）と、
（ｆ）動径分布関数ｇ（ｒ）、補正後の全相関関数ｈ’（ｒ）及び直接相関関数に基づき分子間ポテンシャルＵ（ｒ）を算出するステップ（ＳＴ９）と、
を含む。

このように、波数空間において、全相関関数Ｈ（ｋ）のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数を補正しているので、波数０付近にある値のバラツキを除去でき、ポテンシャル算出の精度を向上させることが可能となる。

本実施形態の装置において、上記補正を実現する好適な例として、
補正部１４は、
波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域において、傾きが０に漸近している漸近領域Ａｒを特定する漸近領域特定部１４ａと、
波数空間において漸近領域Ａｒにおける全相関関数Ｈ（ｋ）と置換関数［Ａ・ｋ^２＋Ｂ］との差に対応する評価値が最小となるように、置換関数［Ａ・ｋ^２＋Ｂ］を構成する係数Ａ，Ｂを決定する係数決定部１４ｂと、
波数空間において漸近領域Ａｒ及び漸近領域Ａｒよりも波数が０に近い領域を置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数Ｈ（ｋ）をＨ’（ｋ）に補正する置換部１４ｃと、を含むことが挙げられる。

本実施形態の方法において、上記補正を実現する好適な例として、
前記波数空間の全相関関数Ｈ（ｋ）を補正するステップ（ｄ）は、
（ｄ１）波数空間の全相関関数Ｈ（ｋ）のうち波数０の周辺領域において、傾きが０に漸近している漸近領域Ａｒを特定するステップ（ＳＴ５）と、
（ｄ２）波数空間において漸近領域Ａｒにおける全相関関数Ｈ（ｋ）と置換関数［Ａ・ｋ^２＋Ｂ］との差に対応する評価値が最小となるように、置換関数［Ａ・ｋ^２＋Ｂ］を構成する係数Ａ，Ｂを決定するステップ（ＳＴ６）と、
（ｄ３）波数空間において漸近領域Ａｒ及び漸近領域Ａｒよりも波数が０に近い領域を置換関数［Ａ・ｋ^２＋Ｂ］に置き換えることで、波数空間の全相関関数Ｈ（ｋ）をＨ’（ｋ）に補正するステップ（ＳＴ７）と、
を含むことが挙げられる。

本実施形態の装置において、分子間ポテンシャル算出部１６は、ＨＮＣ（Hypernetted-chain）近似を用いており、
波数空間において、波数空間全相関関数算出部１３により得られる全相関関数Ｈ（ｋ）と剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）の差に対応する評価値が最小になるように、剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）を構成するパラメータＲ，σを決定するパラメータ決定部１８ａと、
パラメータ決定部１８ａが決定したパラメータＲ，σに基づき剛体球モデルの全相関関数ｈ^ＨＳ（ｒ）および直接相関関数ｃ^ＨＳ（ｒ）を算出し、剛体球モデルのブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）を算出するブリッジ関数算出部１８ｂと、
分子間ポテンシャル算出部１６がＨＮＣ近似を用いて算出したポテンシャルＵ_ＨＮＣ（ｒ）を、ブリッジ関数算出部１８ｂが算出したブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）で修正することで、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルＵ_ＭＨＮＣ（ｒ）を算出する修正部１８ｃと、
を更に備える。

本実施形態の方法において、分子間ポテンシャルを算出するステップ（ＳＴ９）は、ＨＮＣ（Hypernetted-chain）近似を用いており、
波数空間において、ステップＳＴ４で得られる全相関関数Ｈ（ｋ）と剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）の差に対応する評価値が最小になるように、剛体球モデルの全相関関数Ｈ^ＨＳ（ｋ）を構成するパラメータＲ，σを決定するステップ（ＳＴ１０）と、
ステップＳＴ１０で決定したパラメータＲ，σに基づき剛体球モデルの全相関関数ｈ^ＨＳ（ｒ）および直接相関関数ｃ^ＨＳ（ｒ）を算出し、剛体球モデルのブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）を算出するステップ（ＳＴ１１）と、
ステップＳＴ９でＨＮＣ近似を用いて算出したポテンシャルＵ_ＨＮＣ（ｒ）を、ステップＳＴ１１で算出したブリッジ関数ｂ_ＰＹ ^ＨＳ（ｒ）で修正することで、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルＵ_ＭＨＮＣ（ｒ）を算出するステップ（ＳＴ１２）と、
を含む。

このように、ＨＮＣ近似で導出したポテンシャルを更に修正することで、更に精度を向上させることが可能となる。

本実施形態に係るコンピュータプログラムは、上記算出方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記算出方法の奏する作用効果を得ることが可能となる。

以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。

例えば、図１、２に示す各部１０〜１８は、所定プログラムをコンピュータのＣＰＵで実行することで実現しているが、各部を専用回路で構成してもよい。

上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。

１、１’…分子間ポテンシャルの算出装置
１１…分子動力学算出部
１２…動径分布関数算出部
１３…波数空間全相関関数算出部
１４…補正部
１４ａ…漸近領域特定部
１４ｂ…係数決定部
１４ｃ…置換部
１５…直接相関関数算出部
１６…分子間ポテンシャル算出部
１８ａ…パラメータ決定部
１８ｂ…ブリッジ関数算出部
１８ｃ…修正部

Claims (7)

1. 液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出する分子動力学算出部と、
   前記分子動力学算出部が算出した座標位置に基づき動径分布関数を算出する動径分布関数算出部と、
   前記動径分布関数算出部が算出した前記動径分布関数に基づき全相関関数を算出し、算出した全相関関数をフーリエ変換により波数空間の全相関関数に変換する波数空間全相関関数算出部と、
   前記波数空間の全相関関数のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数に置き換えることで、波数空間の全相関関数を補正する補正部と、
   補正後の全相関関数に基づき直接相関関数を算出する直接相関関数算出部と、
   動径分布関数、補正後の全相関関数及び直接相関関数に基づき分子間ポテンシャルを算出する分子間ポテンシャル算出部と、
   を備える分子間ポテンシャルの算出装置。
2. 前記補正部は、
   波数空間の全相関関数のうち波数０の周辺領域において、傾きが０に漸近している漸近領域を特定する漸近領域特定部と、
   波数空間において前記漸近領域における全相関関数と前記置換関数との差に対応する評価値が最小となるように、前記置換関数を構成する係数を決定する係数決定部と、
   波数空間において前記漸近領域及び当該漸近領域よりも波数が０に近い領域を前記置換関数に置き換えることで、波数空間の全相関関数を補正する置換部と、
   を含む請求項１に記載の分子間ポテンシャルの算出装置。
3. 前記分子間ポテンシャル算出部は、ＨＮＣ（Hypernetted-chain）近似を用いており、
   波数空間において、前記波数空間全相関関数算出部により得られる全相関関数と剛体球モデルの全相関関数の差に対応する評価値が最小になるように、前記剛体球モデルの全相関関数を構成するパラメータを決定するパラメータ決定部と、
   前記パラメータ決定部が決定したパラメータに基づき前記剛体球モデルの全相関関数および直接相関関数を算出し、剛体球モデルのブリッジ関数を算出するブリッジ関数算出部と、
   前記分子間ポテンシャル算出部がＨＮＣ近似を用いて算出したポテンシャルを、前記ブリッジ関数算出部が算出したブリッジ関数で修正することで、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルを算出する修正部と、
   を更に備える請求項１又は２に記載の分子間ポテンシャルの算出装置。
4. コンピュータが実行する方法であって、
   （ａ）液体分子モデルの分子動力学計算に基づき平衡状態にある分子モデルの座標位置を算出するステップと、
   （ｂ）算出した座標位置に基づき動径分布関数を算出するステップと、
   （ｃ）算出した前記動径分布関数に基づき全相関関数を算出し、算出した全相関関数をフーリエ変換により波数空間の全相関関数に変換するステップと、
   （ｄ）波数空間の全相関関数のうち波数０の周辺領域を、波数０にて傾きが０となる置換関数に置き換えることで、波数空間の全相関関数を補正するステップと、
   （ｅ）補正後の全相関関数に基づき直接相関関数を算出するステップと、
   （ｆ）動径分布関数、補正後の全相関関数及び直接相関関数に基づき分子間ポテンシャルを算出するステップと、
   を含む分子間ポテンシャルの算出方法。
5. 前記波数空間の全相関関数を補正するステップ（ｄ）は、
   （ｄ１）波数空間の全相関関数のうち波数０の周辺領域において、傾きが０に漸近している漸近領域を特定するステップと、
   （ｄ２）波数空間において前記漸近領域における全相関関数と前記置換関数との差に対応する評価値が最小となるように、前記置換関数を構成する係数を決定するステップと、
   （ｄ３）波数空間において前記漸近領域及び当該漸近領域よりも波数が０に近い領域を前記置換関数に置き換えることで、波数空間の全相関関数を補正するステップと、
   を含む請求項４に記載の分子間ポテンシャルの算出方法。
6. 前記分子間ポテンシャルを算出するステップ（ｆ）は、ＨＮＣ（Hypernetted-chain）近似を用いており、
   （ｇ）波数空間において、ステップ（ｃ）で得られる全相関関数と剛体球モデルの全相関関数の差に対応する評価値が最小になるように、前記剛体球モデルの全相関関数を構成するパラメータを決定するステップと、
   （ｈ）ステップ（ｇ）で決定したパラメータに基づき前記剛体球モデルの全相関関数および直接相関関数を算出し、剛体球モデルのブリッジ関数を算出するステップと、
   （ｉ）ステップ（ｆ）でＨＮＣ近似を用いて算出したポテンシャルを、ステップ（ｈ）で算出したブリッジ関数で修正することで、ＭＨＮＣ（Modified Hypernetted-chain）近似による分子間ポテンシャルを算出するステップと、
   を更に含む請求項４又は５に記載の分子間ポテンシャルの算出方法。
7. 請求項４〜６のいずれかに記載の方法をコンピュータに実行させるコンピュータプログラム。