(57)【要約】【課題】全粒子モデルの動径分布関数との一致を損なうことなく、全粒子モデルの圧力値を再現する粗視化モデルの非結合ポテンシャルを算出する装置を提供する。【解決手段】装置１は、全粒子モデルの圧力値Ｐ全粒子を算出する全粒子モデル圧力値算出部１２と、全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルＵＩＢＩ（ｒ）を算出する非結合ポテンシャル算出部１３と、粗視化モデルの圧力値Ｐ粗視化が全粒子モデルの圧力値Ｐ全粒子と一致するように、粗視化モデルの補正後の非結合ポテンシャルを算出する補正後非結合ポテンシャル算出部１４と、を有する。算出部１４は、補正関数Ｕ補正（ｒ）を補正前の非結合ポテンシャルＵＩＢＩ（ｒ）に加えることで補正後の非結合ポテンシャルを算出する。補正関数Ｕ補正（ｒ）は、ｒで微分した関数のｒにカットオフ距離ｒcutを代入した値が０となる関数である。【選択図】図２
2015103030.tif 000002
【特許請求の範囲】
【請求項１】
全粒子モデルの分子動力学計算に基づき算出される動径分布関数及び非結合ポテンシャルを用い、当該動径分布関数及び非結合ポテンシャルのパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値を算出するステップと、
全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルを算出するステップと、
粗視化モデルの動径分布関数及び補正後の非結合ポテンシャルを用いて前記ｒの積分により算出される粗視化モデルの圧力値が前記全粒子モデルの圧力値と一致するように、粗視化モデルの補正後の非結合ポテンシャルを算出するステップと、を含み、
前記補正後の非結合ポテンシャルを算出するステップでは、パラメータの一つとして前記ｒを有する補正関数を補正前の非結合ポテンシャルに加えることで補正後の非結合ポテンシャルを算出し、前記補正関数は、前記ｒで微分した関数の前記ｒにカットオフ距離を代入した値が０となる関数であることを特徴とする粗視化モデルの非結合ポテンシャルを算出する方法。
【請求項２】
前記補正関数は、前記ｒをパラメータの一つとする二次関数である、請求項１に記載の粗視化モデルの非結合ポテンシャルを算出する方法。
【請求項３】
前記補正関数は、次の式（７）で表される関数である請求項１又は２に記載の粗視化モデルの非結合ポテンシャルを算出する方法。
2015103030.tif 000009
Ｎ粗視化は粗視化モデルの粒子数を表し、ｇ粗視化（ｒ）は粗視化モデルの動径分布関数を表し、ＵＩＢＩ（ｒ）はＩＢＩによって算出された粗視化モデルの非結合ポテンシャルを表し、Ｖは体積を表し、Ｔは温度を表し、Ｐ全粒子は全粒子モデルの圧力を表し、ｋＢはボルツマン定数を表す。
【請求項４】
請求項１〜３のいずれかに記載の方法をコンピュータに実行させるコンピュータプログラム。
【請求項５】
全粒子モデルの分子動力学計算に基づき算出される動径分布関数及び非結合ポテンシャルを用い、当該動径分布関数及び非結合ポテンシャルのパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値を算出する全粒子モデル圧力値算出部と、
全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルを算出するポテンシャル算出部と、
粗視化モデルの動径分布関数及び補正後の非結合ポテンシャルを用いて前記ｒの積分により算出される粗視化モデルの圧力値が前記全粒子モデルの圧力値と一致するように、粗視化モデルの補正後の非結合ポテンシャルを算出する補正後ポテンシャル算出部と、を備え、
前記補正後ポテンシャル算出部は、パラメータの一つとして前記ｒを有する補正関数を補正前の非結合ポテンシャルに加えることで補正後の非結合ポテンシャルを算出し、前記補正関数は、前記ｒで微分した関数の前記ｒにカットオフ距離を代入した値が０となる関数であることを特徴とする、粗視化モデルの非結合ポテンシャルを算出する装置。
【請求項６】
前記補正関数は、前記ｒをパラメータの一つとする二次関数である、請求項５に記載の粗視化モデルの非結合ポテンシャルを算出する装置。
【請求項７】
前記補正関数は、次の式（７）で表される関数である請求項５又は６に記載の粗視化モデルの非結合ポテンシャルを算出する装置。
2015103030.tif 000010
Ｎ粗視化は粗視化モデルの粒子数を表し、ｇ粗視化（ｒ）は粗視化モデルの動径分布関数を表し、ＵＩＢＩ（ｒ）はＩＢＩによって算出された粗視化モデルの非結合ポテンシャルを表し、Ｖは体積を表し、Ｔは温度を表し、Ｐ全粒子は全粒子モデルの圧力を表し、ｋＢはボルツマン定数を表す。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粗視化分子動力学によるシミュレーションに用いる粗視化モデルの非結合ポテンシャルを算出する方法、粗視化モデルの非結合ポテンシャルを算出する装置、及びコンピュータプログラムに関する。
【背景技術】
【０００２】
評価対象となる材料の物性の予測及び評価を行うためのシミュレーションの手法として、原子を最小単位とする原子モデルを用いた分子動力学が知られている。しかし、原子モデルを用いた分子動力学では、評価対象となる材料が原子レベルで再現されるので、大規模及び長時間のシミュレーションを実施するためには、莫大な計算コストが必要となり、現実的ではないという問題がある。
【０００３】
この問題を解決する有効な手段の一つとして、粗視化分子動力学が知られている。粗視化分子動力学は、複数の原子を一括りとした粗視化粒子を最小単位のモデルとして取り扱うことにより、計算コストを低減し、原子モデルを用いた分子動力学よりも大規模及び長時間のシミュレーションを取り扱うことができる手法である。粗視化分子動力学を用いるためには、粗視化モデル間に作用する非結合ポテンシャルを原子レベルの分子動力学計算に基づき算出する必要があり、この算出方法が例えば非特許文献１に開示されている。
【０００４】
具体的には、原子モデル（全粒子モデル）の分子動力学計算に基づき、平衡状態の粗視化粒子の重心位置の動径分布関数を得る。これが目標となる動径分布関数である。次に、Boltzmann Inversionを用いて目標となる動径分布関数に基づき平均力ポテンシャルを得る。次に、平均力ポテンシャルを粗視化モデルの非結合ポテンシャルとし、粗視化モデルの分子動力学計算を行い、平衡状態の粗視化粒子の重心位置の動径分布関数を得る。次に、粗視化モデルの動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下であるか否かを判定し、前記評価関数が所定閾値よりも大きいときは、評価関数が小さくなる方向に前記パラメータとしての非結合ポテンシャルを補正し、前記評価関数が所定閾値よりも小さくなるまで、非結合ポテンシャルの補正及び補正後の非結合ポテンシャルに用いた粗視化モデルの分子動力学計算を繰り返す。
【０００５】
このように、ＢＩ（Boltzmann Inversion）の繰り返し（ＩＢＩ；Iterative Boltzmann Inversion）により得られた粗視化モデルの動径分布関数及び非結合ポテンシャルを用いて、所定体積及び所定温度の条件下における圧力値を計算することができる。圧力値は、動径分布関数及び非結合ポテンシャルのパラメータの一つである粒子間距離ｒの積分によって得られる。しかしながら、粗視化モデルの圧力値が全粒子モデル（原子モデル）の圧力値と大きく異なってしまうことが知られている。非特許文献１では、補正前の非結合ポテンシャルＵIBI（ｒ）に、補正関数［Ｕ補正（ｒ）＝Ａ（１−（ｒ／ｒcut））］を加えることで補正後の非結合ポテンシャルを得て、この補正後の非結合ポテンシャルを用いて算出した粗視化モデルの圧力値が全粒子モデルの圧力値と一致することが開示されている。Ａは定数である。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】"Victor Ruhle, Christoph Junghans, Alexander Lukyanov, Kurt Kremer, and Denis Andrienko"著、「Versatile Object-Oriented Toolkit for Coarse-Graining Applications」、Journal of Chemical Theory and Conputation、Vol. 5, No. 12 ,2009、p.3211-3223
【発明の概要】
【発明が解決しようとする課題】
【０００７】
上記補正を実施すれば、全粒子モデルと粗視化モデルの圧力値とを一致させることができる。しかしながら、図４Ａに示すように、カットオフ距離ｒcut付近で力が不連続となり、その結果、図４Ｂに示すように、補正後の非結合ポテンシャルを用いた粗視化モデルの分子動力学計算により得られる動径分布関数が全粒子モデルの動径分布関数と大きく異なってしまう不具合が生じる場合がある。この不具合は、ポリマー等、粒子間に作用する力及びポテンシャルがソフトな場合に、力の不連続性が顕著に現れることがわかった。
【０００８】
本発明は、このような課題に着目してなされたものであって、その目的は、全粒子モデルの動径分布関数との一致を損なうことなく、全粒子モデルの圧力値を再現する粗視化モデルの非結合ポテンシャルを算出する方法、装置及びコンピュータプログラムを提供することである。
【課題を解決するための手段】
【０００９】
本発明は、上記目的を達成するために、次のような手段を講じている。
【００１０】
すなわち、本発明の粗視化モデルの非結合ポテンシャルを算出する方法は、全粒子モデルの分子動力学計算に基づき算出される動径分布関数及び非結合ポテンシャルを用い、当該動径分布関数及び非結合ポテンシャルのパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値を算出するステップと、全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルを算出するステップと、粗視化モデルの動径分布関数及び補正後の非結合ポテンシャルを用いて前記ｒの積分により算出される粗視化モデルの圧力値が前記全粒子モデルの圧力値と一致するように、粗視化モデルの補正後の非結合ポテンシャルを算出するステップと、を含み、前記補正後の非結合ポテンシャルを算出するステップでは、パラメータの一つとして前記ｒを有する補正関数を補正前の非結合ポテンシャルに加えることで補正後の非結合ポテンシャルを算出し、前記補正関数は、前記ｒで微分した関数の前記ｒにカットオフ距離を代入した値が０となる関数であることを特徴とする。
【００１１】
本発明の粗視化モデルの非結合ポテンシャルを算出する装置は、全粒子モデルの分子動力学計算に基づき算出される動径分布関数及び非結合ポテンシャルを用い、当該動径分布関数及び非結合ポテンシャルのパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値を算出する全粒子モデル圧力値算出部と、全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルを算出するポテンシャル算出部と、粗視化モデルの動径分布関数及び補正後の非結合ポテンシャルを用いて前記ｒの積分により算出される粗視化モデルの圧力値が前記全粒子モデルの圧力値と一致するように、粗視化モデルの補正後の非結合ポテンシャルを算出する補正後ポテンシャル算出部と、を備え、前記補正後ポテンシャル算出部は、パラメータの一つとして前記ｒを有する補正関数を補正前の非結合ポテンシャルに加えることで補正後の非結合ポテンシャルを算出し、前記補正関数は、前記ｒで微分した関数の前記ｒにカットオフ距離を代入した値が０となる関数であることを特徴とする。
【００１２】
このように、粗視化モデルの圧力値を全粒子モデルの圧力値に一致させるための補正関数に、粒子間距離ｒで微分した関数のｒにカットオフ距離を代入した値が０となる関数を用いているので、カットオフ距離での力の不連続性をなくすことができ、全粒子モデルの動径分布関数との一致を損なうことなく、全粒子モデルの圧力値を再現する粗視化モデルの非結合ポテンシャルを算出可能となる。
【００１３】
適用例としては、前記補正関数は、前記ｒをパラメータの一つとする二次関数であることが好ましい。
【００１４】
適用例としては、前記補正関数は、式（７）で表される関数であることが好ましい。
【００１５】
本発明は、上記方法を構成するステップをコンピュータに実行させるプログラムとして特定可能である。
【図面の簡単な説明】
【００１６】
【図１】本発明に係る粗視化モデルの非結合ポテンシャルを算出する方法を示すフローチャート。
【図２】本発明に係る非結合ポテンシャルを算出する装置を模式的に示すブロック図。
【図３】本実施形態で用いる原子モデル及び粗視化モデルの一例を模式的に示す説明図。
【図４Ａ】従来技術の不具合に関する説明図。
【図４Ｂ】従来技術の不具合に関する説明図。
【発明を実施するための形態】
【００１７】
以下、本発明の一実施形態を、図面を参照して説明する。
【００１８】
［粗視化モデルの非結合ポテンシャルを算出する装置］
本実施形態の装置は、粗視化分子動力学計算に必要となる粗視化モデルの非結合ポテンシャルを算出する装置である。具体的には、装置は、原子モデル（全粒子モデル）の分子動力学に基づき目標となる動径分布関数を取得し、目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを算出する。算出した非結合ポテンシャルに基づき導出された粗視化モデルの圧力値は、全粒子モデルの圧力値と大きく乖離するために、双方の圧力値が一致するように補正した非結合ポテンシャルを算出する必要がある。本装置は、適切に補正がなされた補正後の非結合ポテンシャルを算出する。
【００１９】
本実施形態では、図３に例示するように、全粒子モデルの分子動力学では、５つの粒子Ｍ１〜Ｍ５の各々を最小単位として取り扱うのに対し、粗視化モデルでは、例えば５つ粒子を一括りとした粗視化粒子Ｍ６を最小単位のモデルとして取り扱う。なお、粗視化粒子を構成する粒子の数及びその構造は目的に応じて任意に変更可能である。
【００２０】
図２に示すように、装置１は、初期設定部１０と、全粒子モデル分子動力学算出部１１と、全粒子モデル圧力値算出部１２と、粗視化モデルの非結合ポテンシャル算出部１３と、補正後非結合ポテンシャル算出部１４と、を有する。これら各部１０〜１４は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図示しない処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。
【００２１】
図２に示す初期設定部１０は、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、解析対象となる全粒子モデル（原子モデル、分子モデル）や粗視化モデルに関する情報の設定、分子動力学計算に必要な設定、圧力値の算出に必要な所定体積Ｖや所定温度Ｔ等の各種設定を実行し、これら設定値を図示しないメモリに記憶する。全粒子モデルや粗視化モデルは、図３に例示するものが挙げられる。
【００２２】
図２に示す全粒子モデル分子動力学算出部１１は、予め設定されている全粒子モデル、その粒子間相互作用（非結合ポテンシャル等）Ｕ全粒子（ｒ）等の各種パラメータを用い、全粒子モデルの分子動力学計算に基づき平衡状態における各粒子モデルの重心位置を算出する。これに伴い、粗視化モデルを構成する全粒子モデルの重心位置に基づき粗視化モデルの重心位置を算出し、粗視化モデルの目標となる動径分布関数Ｇ（ｒ）を算出する。また、全粒子モデルの動径分布関数ｇ全粒子（ｒ）を算出する。
【００２３】
図２に示す全粒子モデル圧力値算出部１２は、全粒子モデルの分子動力学計算に基づき算出される動径分布関数ｇ全粒子（ｒ）及び非結合ポテンシャルＵ全粒子（ｒ）を用い、所定体積Ｖ及び所定温度Ｔを条件として、動径分布関数ｇ全粒子（ｒ）及び非結合ポテンシャルＵ全粒子（ｒ）のパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値Ｐ全粒子を算出する。具体的に、全粒子モデル圧力値算出部１２は、次の式（１）を演算する。
2015103030.tif 000003
ここで、Ｐ全粒子は全粒子モデルの圧力を表し、Ｎ全粒子は全粒子モデルの粒子数を表し、ｇ全粒子（ｒ）は全粒子モデルの動径分布関数を表し、Ｕ全粒子（ｒ）は全粒子モデルの粒子間相互作用を表す。
【００２４】
図２に示す粗視化モデルの非結合ポテンシャル算出部１３は、目標となる動径分布関数Ｇ（ｒ）を粗視化モデルにおいて再現する非結合ポテンシャルＵＩＢＩ（ｒ）を算出する。具体的には、まず、目標となる動径分布関数Ｇ（ｒ）に基づき下記式（２）を用いて平均力ポテンシャルＵＰＭＦを算出し、これを下記式（３）に示す粗視化モデルの非結合ポテンシャルＵ（０）とする。
ＵＰＭＦ（ｒ）＝−ｋＢＴ ln Ｇ（ｒ）  …（２）
Ｕ（０）＝−ｋＢＴ ln Ｇ（ｒ）  …（３）
式（２）のように動径分布関数から平均力ポテンシャルエネルギー関数に変換する処理はBoltzmann Inversionと呼ばれる。Ｕ（０）は初期値を意味する。ｋＢはボルツマン定数を表し、Ｔは温度を表す。
【００２５】
次に、非結合ポテンシャル算出部１３は、予め設定されている粗視化モデル、パラメータとしての非結合ポテンシャルＵ（ｉ）を用いて粗視化モデルの分子動力学計算に基づき平衡状態における粗視化モデルの重心位置を算出し、粗視化モデルの動径分布関数ｇ（ｉ）（ｒ）を取得する。非結合ポテンシャル算出部１３は、粗視化モデルの動径分布関数ｇ（ｉ）（ｒ）と、上記目標となる動径分布関数Ｇ（ｒ）との差値に対応する評価関数が所定閾値以下であるか否かを判定する。評価関数は、粗視化モデルと原子モデルの動径分布関数の平均二乗偏差である。判定部は、式（４）を満足するか否かを判定する。Δは所定閾値であり、Ｒはカットオフ長さである。これらは予めメモリに設定されている。
2015103030.tif 000004
【００２６】
非結合ポテンシャル算出部１３は、評価関数が所定閾値以下のときは、パラメータとして設定されている非結合ポテンシャルＵ（ｉ）が解であるとして、解である非結合ポテンシャルＵＩＢＩ（ｒ）を出力する。一方、評価関数が所定閾値よりも大きいときには、評価関数が小さくなる方向に非結合ポテンシャルを補正し、補正した非結合ポテンシャルを新たなパラメータとして粗視化モデルの分子動力学計算を繰り返す。下記式（５）に示すように、補正した非結合ポテンシャルはＵ（ｉ＋１）であり、式（５）の右辺の第２項が補正値となる。ｋＢはボルツマン定数を表し、Ｔは温度を表す。
2015103030.tif 000005
【００２７】
補正後非結合ポテンシャル算出部１４は、粗視化モデルの動径分布関数ｇ粗視化（ｒ）及び補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を用いてｒの積分により算出される粗視化モデルの圧力値Ｐ粗視化が全粒子モデルの圧力値Ｐ全粒子と一致するように、粗視化モデルの補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を算出する。
【００２８】
非結合ポテンシャルを補正するのは、粗視化モデルにおける動径分布関数ｇ粗視化（ｒ）及び非結合ポテンシャルＵＩＢＩ（ｒ）に基づき次の式（６）を用いて算出される粗視化モデルの圧力値Ｐ粗視化が全粒子モデルの圧力値Ｐ全粒子に一致していないといけないところ、異なる（小さくなってしまうことが多い）という欠点を補償するためである。
2015103030.tif 000006
ここで、Ｐ粗視化は粗視化モデルの圧力を表し、Ｎ粗視化は粗視化モデルの粒子数を表し、ｇ粗視化（ｒ）は粗視化モデルの動径分布関数を表し、ＵＩＢＩ（ｒ）はＩＢＩによって算出された粗視化モデルの非結合ポテンシャルを表す。
【００２９】
このとき、補正後の非結合ポテンシャルを［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］とし、補正関数［Ｕ補正（ｒ）＝Ａ（１−（ｒ／ｒcut））］とした場合には、図４Ａに示すように、カットオフ距離ｒcut付近で補正関数の力が不連続となり、その結果、図４Ｂに示すように、補正後の非結合ポテンシャルを用いた粗視化モデルの分子動力学計算により得られる動径分布関数ｇ’（ｒ）が全粒子モデルの動径分布関数Ｇ（ｒ）と大きく異なってしまう不具合が生じる。
【００３０】
この原因について説明する。分子動力学計算では、粒子間相互作用は粒子間距離が離れるほど急激に小さくなる。そのため、ある適切な距離を設け、それよりも離れた粒子間の相互作用を無視する（力及びポテンシャルが共に０）ことで計算精度を保ったまま計算コストの低減を図っている。このような相互作用の打ち切り距離をカットオフ距離ｒcutとしている。補正関数をｒで微分してマイナスをとれば力（フォース）となるが、補正関数［Ａ（１−ｒ／ｒcut）］をｒで微分した関数は、Ａ／ｒcutとなり、ｒにカットオフ距離ｒcutを代入しても０にならず、定数となる。すなわち、力が不連続となっているのは、図４Ａにおけるカットオフ距離ｒcutの左側では力が０でないので、或る値が足し込まれるのに対し、カットオフ距離ｒcutの右側では力が０と見なされるので、カットオフ距離ｒcutにおいて力が不連続となることがわかった。この力の不連続性が、図４Ｂに示すように動径分布関数を狂わす要因となる。
【００３１】
そこで、本実施形態では、補正関数Ｕ補正（ｒ）には、ｒで微分した関数のｒにカットオフ距離ｒcutを代入した値が０となる関数を用いる。具体的には、ｒをパラメータの一つとする二次関数であり、次の式（７）で表される。
2015103030.tif 000007
Ｎ粗視化は粗視化モデルの粒子数を表し、ｇ粗視化（ｒ）は粗視化モデルの動径分布関数を表し、ＵＩＢＩ（ｒ）はＩＢＩによって算出された粗視化モデルの非結合ポテンシャルを表し、Ｖは体積を表し、Ｔは温度を表し、Ｐ全粒子は全粒子モデルの圧力を表し、ｋＢはボルツマン定数を表す。
【００３２】
この補正関数Ｕ補正（ｒ）をｒで微分してマイナスをとると、次の式（８）となる。式（８）のｒにカットオフ距離ｒcutを代入すれば、値が０となる。
2015103030.tif 000008
【００３３】
このように補正関数Ｕ補正（ｒ）を算出し、補正後の非結合ポテンシャルを算出すれば、動径分布関数のズレを防止すると共に、全粒子モデルの圧力値を再現できた。
【００３４】
なお、本実施形態では、補正関数Ｕ補正（ｒ）に、式（７）に示される二次関数を利用しているが、ｒで微分してｒにカットオフ距離ｒcutを代入した値が０になれば、これに特に限定されない。例えば、二次関数でなくとも三次関数、四次関数でもよい。また、その他、Ｕ補正（ｒ）＝Ａ（ｅｘｐ（ｒcut−ｒ）−１）２
Ｕ補正（ｒ）＝Ａ ｓｉｎ（（ｒcut−ｒ）２）
を用いてもよい。Ａは比例定数である。
【００３５】
［粗視化モデルの非結合ポテンシャルを算出する方法］
上記装置１を用いて粗視化モデルの非結合ポテンシャルを算出する方法について説明する。
【００３６】
まず、ステップＳＴ１（図１参照）において、図２に示す初期設定部１０は、解析対象となる全粒子モデル（原子モデル）や粗視化モデルに関する情報の設定、分子動力学計算に必要な設定等の各種設定を実行し、これら設定値を図示しないメモリに記憶する。
【００３７】
次のステップＳＴ２（図１参照）において、図２に示す全粒子モデル分子動力学算出部１１は、予め設定された各種パラメータを用い、分子動力学計算に基づき平衡状態にある全粒子モデルの重心位置を算出し、動径分布関数ｇ全粒子（ｒ）を算出する。
【００３８】
次のステップＳＴ３（図１参照）において、全粒子モデル圧力値算出部１２は、全粒子モデル分子動力学算出部１１が算出した分布関数ｇ全粒子（ｒ）と、全粒子モデルの非結合ポテンシャルＵ全粒子（ｒ）、所定体積Ｖと所定温度Ｔに基づき全粒子モデルでの圧力値を算出する。
【００３９】
次のステップＳＴ４（図１参照）において、非結合ポテンシャル算出部１３は、全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルを算出する。具体的には、全粒子モデル分子動力学算出部１１の算出結果から得られた目標となる動径分布関数Ｇ（ｒ）に基づき平均力ポテンシャルＵＰＭＦを算出し、これを粗視化モデルの非結合ポテンシャルＵ（０）としてメモリにパラメータとしてセットする。次に、非結合ポテンシャルＵ（ｉ）を用いて粗視化モデルの分子動力学計算に基づき平衡状態における粗視化モデルの重心位置を算出し、動径分布関数ｇ（ｉ）（ｒ）を取得する。次に、粗視化モデルの動径分布関数ｇ（ｉ）（ｒ）と、目標となる動径分布関数Ｇ（ｒ）との差値に対応する評価関数が所定閾値以下であるか否かを判定する。評価関数が所定閾値以下でないと判定された場合には、判定部１７は、パラメータである非結合ポテンシャルを補正し、上記粗視化モデルの分子動力学計算の実行に戻る。一方、評価関数が所定閾値以下であると判定された場合には、パラメータの非結合ポテンシャルＵ（ｉ）を解として特定する。
【００４０】
次のステップＳＴ５において、補正後非結合ポテンシャル算出部１４は、粗視化モデルの動径分布関数ｇ粗視化（ｒ）及び補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を用いて算出される粗視化モデルの圧力値Ｐ粗視化が全粒子モデルの圧力値Ｐ全粒子と一致するように、粗視化モデルの補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を算出する。具体的には、式（７）を用いて、補正関数Ｕ補正（ｒ）を算出して、ＵＩＢＩ（ｒ）に加算して補正し、補正後の非結合ポテンシャルを得る。
【００４１】
以上のように、本実施形態の粗視化モデルの非結合ポテンシャルを算出する方法は、全粒子モデルの分子動力学計算に基づき算出される動径分布関数ｇ全粒子（ｒ）及び非結合ポテンシャルＵ全粒子（ｒ）を用い、動径分布関数ｇ全粒子（ｒ）及び非結合ポテンシャルＵ全粒子（ｒ）のパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値Ｐ全粒子を算出するステップ（ＳＴ３）と、全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルＵＩＢＩ（ｒ）を算出するステップ（ＳＴ４）と、粗視化モデルの動径分布関数ｇ粗視化（ｒ）及び補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を用いてｒの積分により算出される粗視化モデルの圧力値Ｐ粗視化が全粒子モデルの圧力値Ｐ全粒子と一致するように、粗視化モデルの補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を算出するステップ（ＳＴ５）と、を含む。補正後の非結合ポテンシャルを算出するステップでは、パラメータの一つとしてｒを有する補正関数Ｕ補正（ｒ）を補正前の非結合ポテンシャルＵＩＢＩ（ｒ）に加えることで補正後の非結合ポテンシャルを算出する。補正関数Ｕ補正（ｒ）は、ｒで微分した関数のｒにカットオフ距離ｒcutを代入した値が０となる関数である。
【００４２】
本実施形態の粗視化モデルの非結合ポテンシャルを算出する装置１は、全粒子モデルの分子動力学計算に基づき算出される動径分布関数ｇ全粒子（ｒ）及び非結合ポテンシャルＵ全粒子（ｒ）を用い、動径分布関数ｇ全粒子（ｒ）及び非結合ポテンシャルＵ全粒子（ｒ）のパラメータの一つである粒子間距離ｒの積分により全粒子モデルの圧力値Ｐ全粒子を算出する全粒子モデル圧力値算出部１２と、全粒子モデルの動径分布関数を再現する粗視化モデルの非結合ポテンシャルＵＩＢＩ（ｒ）を算出する非結合ポテンシャル算出部１３と、粗視化モデルの動径分布関数ｇ粗視化（ｒ）及び補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を用いてｒの積分により算出される粗視化モデルの圧力値Ｐ粗視化が全粒子モデルの圧力値Ｐ全粒子と一致するように、粗視化モデルの補正後の非結合ポテンシャル［ＵＩＢＩ（ｒ）＋Ｕ補正（ｒ）］を算出する補正後非結合ポテンシャル算出部１４と、を備える。補正後非結合ポテンシャル算出部１４は、パラメータの一つとしてｒを有する補正関数Ｕ補正（ｒ）を補正前の非結合ポテンシャルＵＩＢＩ（ｒ）に加えることで補正後の非結合ポテンシャルを算出する。補正関数Ｕ補正（ｒ）は、ｒで微分した関数のｒにカットオフ距離ｒcutを代入した値が０となる関数である。
【００４３】
このように、粗視化モデルの圧力値を全粒子モデルの圧力値に一致させるための補正関数Ｕ補正（ｒ）に、粒子間距離ｒで微分した関数のｒにカットオフ距離ｒcutを代入した値が０となる関数を用いているので、カットオフ距離ｒcutでの力の不連続性をなくすことができ、全粒子モデルの動径分布関数との一致を損なうことなく、全粒子モデルの圧力値を再現する粗視化モデルの非結合ポテンシャルを算出可能となる。
【００４４】
本実施形態では、補正関数Ｕ補正（ｒ）は、ｒをパラメータの一つとする二次関数である。補正関数は、上記条件を満たせば、三次関数、四次関数を用いることができるが、比較的式が簡素となる二次関数が好ましい。
【００４５】
本実施形態では、補正関数は、式（７）で表される関数であり、本発明の好適な例である。
【００４６】
本実施形態に係るコンピュータプログラムは、上記方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。言い換えると、上記方法を使用しているとも言える。
【００４７】
以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【００４８】
例えば、図２に示す各部１０〜１７は、所定プログラムをコンピュータのＣＰＵで実行することで実現しているが、各部を専用メモリや専用回路で構成してもよい。
【００４９】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【００５０】
１２…全粒子モデル圧力値算出部
１３…非結合ポテンシャル算出部
１４…補正後非結合ポテンシャル算出部
【図１】
2015103030.tif 000011
【図２】
2015103030.tif 000012
【図３】
2015103030.tif 000013
【図４Ａ】
2015103030.tif 000014
【図４Ｂ】
2015103030.tif 000015