(57)【要約】【課題】離散値のサンプル点に依存せずに動径分布関数を取得可能にして、信頼性を高めた粗視化モデルの非結合ポテンシャルの算出方法を提供する。【解決手段】算出方法は、原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの動径分布関数の基礎となる離散データを算出し、当該離散データを近似して目標となる動径分布関数Ｇ（ｒ）を取得するステップＳＴ２〜５と、粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの動径分布関数の基礎となる離散データを算出し、当該離散データを近似して動径分布関数ｇ（ｉ）（ｒ）を取得し、当該動径分布関数ｇ（ｉ）（ｒ）と上記目標となる動径分布関数Ｇ（ｒ）を含む評価関数が小さくなるように、非結合ポテンシャルを補正及び粗視化分子動力学計算を繰り返すステップＳＴ６〜１２とを含む。離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得する。【選択図】図１
2014102554.tif 000002
【特許請求の範囲】
【請求項１】
原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを、コンピュータが算出する方法であって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して目標となる動径分布関数を取得するステップと、
前記目標となる動径分布関数に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルにパラメータとして設定し、粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して動径分布関数を取得し、粗視化分子動力学に基づく動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下となるか否かを判定し、前記評価関数が前記所定閾値よりも大きいときは、前記評価関数が小さくなる方向に補正した前記非結合ポテンシャルを新たなパラメータとして前記評価関数が前記所定閾値以下となるまで粗視化分子動力学計算を繰り返すステップと、を含み、
前記粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得することを特徴とする粗視化モデルの非結合ポテンシャルの算出方法。
【請求項２】
前記窓関数は、次の式（１）で表されるガウス関数である請求項１に記載の粗視化モデルの非結合ポテンシャルの算出方法。
ｈ（ｕ）＝（１／√（２πσ２））ｅｘｐ（−ｕ２／２σ２）  …（１）
σは標準偏差である。
【請求項３】
動径分布関数の基礎となる離散データに基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定するステップと、
上記ガウス関数のσの値を種々異ならせて平滑化した動径分布関数を算出するステップと、
各々の算出結果について、平滑化後の動径分布関数をｇ（ｒ）とし、平滑化前の動径分布関数をｇ’（ｒ）とする場合に、０．９×ｒmax≦ｒ≦１．１×ｒmaxの範囲において、下記式（６）に示すＲを算出するステップと、
上記複数の算出結果から、Ｒ≦１０−２を満たすσの最大値σmaxを抽出すると共に、Ｒ≧１０−４を満たすσの最小値σminを抽出し、σmin以上且つσmax以下となるσを決定するステップと、を更に含む請求項２に記載の粗視化モデルの非結合ポテンシャルの算出方法。
2014102554.tif 000006
（ただし、Ｎ＝（１．１×ｒmax−０．９×ｒmax）／ｄｒ、Ｎは四捨五入した整数）
【請求項４】
前記離散データに基づき波形のピーク又はボトムを特定するステップと、
特定した波形のピーク又はボトムに基づきノイズの波長又は周波数を特定するステップと、
特定したノイズの波長又は周波数に応じてガウス関数の標準偏差σを決定するステップと、を更に含む請求項２に記載の粗視化モデルの非結合ポテンシャルの算出方法。
【請求項５】
原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを算出する算出装置であって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出する原子モデル分子動力学算出部と、
前記原子モデル分子動力学算出部が算出した離散データを近似して目標となる動径分布関数を取得する第一動径分布関数取得部と、
前記第一動径分布関数取得部が取得した目標となる動径分布関数に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルにパラメータとして設定する非結合ポテンシャル設定部と、
粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出する粗視化モデル分子動力学算出部と、
前記粗視化モデル分子動力学算出部が算出した離散データを近似して動径分布関数を取得する第二動径分布関数取得部と、
前記第二動径分布関数取得部が取得した粗視化分子動力学に基づく動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下となるか否かを判定する判定部と、を備え、
前記評価関数が前記所定閾値よりも大きいと前記判定部が判定した場合には、前記評価関数が小さくなる方向に補正した前記非結合ポテンシャルを新たなパラメータとし、前記評価関数が前記所定閾値以下となるまで前記粗視化モデル分子動力学算出部による粗視化分子動力学計算を繰り返すように構成されており、
前記第一及び第二の動径分布関数取得部は、前記粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得するように構成されていることを特徴とする粗視化モデルの非結合ポテンシャルの算出装置。
【請求項６】
前記窓関数は、次の式（１）で表されるガウス関数である請求項５に記載の粗視化モデルの非結合ポテンシャルの算出装置。
ｈ（ｕ）＝（１／√（２πσ２））ｅｘｐ（−ｕ２／２σ２）  …（１）
σは標準偏差である。
【請求項７】
動径分布関数の基礎となる離散データに基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定し、上記ガウス関数のσの値を種々異ならせて平滑化した動径分布関数を算出し、各々の算出結果について、平滑化後の動径分布関数をｇ（ｒ）とし、平滑化前の動径分布関数をｇ’（ｒ）とする場合に、０．９×ｒmax≦ｒ≦１．１×ｒmaxの範囲において、下記式（６）に示すＲを算出し、上記複数の算出結果から、Ｒ≦１０−２を満たすσの最大値σmaxを抽出すると共に、Ｒ≧１０−４を満たすσの最小値σminを抽出し、σmin以上且つσmax以下となるσを決定する標準偏差値決定部を更に備える請求項６に記載の粗視化モデルの非結合ポテンシャルの算出装置。
2014102554.tif 000007
（ただし、Ｎ＝（１．１×ｒmax−０．９×ｒmax）／ｄｒ、Ｎは四捨五入した整数）
【請求項８】
離散データに基づき波形のピーク又はボトムを特定し、特定した波形のピーク又はボトムに基づきノイズの波長又は周波数を特定し、特定したノイズの波長又は周波数に応じてガウス関数の標準偏差σを決定する標準偏差値決定部を更に備える請求項６に記載の粗視化モデルの非結合ポテンシャルの算出装置。
【請求項９】
原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを算出するプログラムであって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して目標となる動径分布関数を取得するステップと、
前記目標となる動径分布関数に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルにパラメータとして設定し、粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して動径分布関数を取得し、粗視化分子動力学に基づく動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下となるか否かを判定し、前記評価関数が前記所定閾値よりも大きいときは、前記評価関数が小さくなる方向に補正した前記非結合ポテンシャルを新たなパラメータとして前記評価関数が前記所定閾値以下となるまで粗視化分子動力学計算を繰り返すステップと、をコンピュータに実行させ、
前記粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得することを特徴とする粗視化モデルの非結合ポテンシャルの算出プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粗視化分子動力学によるシミュレーションに用いる粗視化モデルの非結合ポテンシャルの算出方法、粗視化モデルの非結合ポテンシャルの算出装置、及び粗視化モデルの非結合ポテンシャルの算出プログラムに関する。
【背景技術】
【０００２】
従来から、評価対象となる材料の物性の予測及び評価を行うためのシミュレーションの手法として、原子を最小単位とする原子モデルを用いた分子動力学が知られている。しかし、原子モデルを用いた分子動力学では、評価対象となる材料が原子レベルで再現されるので、大規模及び長時間のシミュレーションを実施するためには、莫大な計算コストが必要となり、現実的ではないという問題がある。
【０００３】
この問題を解決する有効な手段の一つとして、粗視化分子動力学が知られている。粗視化分子動力学は、複数の原子を一括りとした粗視化粒子を最小単位のモデルとして取り扱うことにより、計算コストを低減し、原子モデルを用いた分子動力学よりも大規模及び長時間のシミュレーションを取り扱うことができる手法である。粗視化分子動力学を用いるためには、粗視化モデル間に作用する非結合ポテンシャルを原子レベルの分子動力学計算に基づき算出する必要があり、この算出方法が例えば非特許文献１に開示されている。
【０００４】
具体的には、原子モデルの分子動力学計算に基づき、平衡状態の粗視化粒子の重心位置の動径分布関数の基礎となる離散データ（図５Ａ参照）を算出する。次に、動径分布関数の基礎となる離散データを、スプライン関数などの補間法を用いて平滑且つ近似し、目標となる動径分布関数を得る。次に、Boltzmann Inversionを用いて目標となる動径分布関数に基づき平均力ポテンシャルを得る。次に、平均力ポテンシャルを粗視化モデルの非結合ポテンシャルとし、粗視化モデルの分子動力学計算を行い、平衡状態の粗視化粒子の重心位置の動径分布関数の基礎となる離散データを算出する。次に、当該離散データをスプライン関数などの補間法を用いて平滑処理することで、離散データに近似する粗視化モデルの動径分布関数を得る。次に、粗視化モデルの動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下であるか否かを判定し、前記評価関数が所定閾値よりも大きいときは、評価関数が小さくなる方向に前記パラメータとしての非結合ポテンシャルを補正し、前記評価関数が所定閾値よりも小さくなるまで、非結合ポテンシャルの補正及び補正後の非結合ポテンシャルに用いた粗視化モデルの分子動力学計算を繰り返す。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】"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
【発明の概要】
【発明が解決しようとする課題】
【０００６】
すなわち、上記従来技術として、分子動力学計算により得られる動径分布関数の基礎となる離散データから動径分布関数を得るために、スプライン補間などの補間法を利用して平滑化処理及び近似処理を実施している。
【０００７】
しかしながら、スプライン補間などの補間法により得られる動径分布関数は、サンプル点の取り方によって異なってしまう。例えば図５Ｂ及び図５Ｃに例示する。図中では、丸印が離散データのサンプル点であり、サンプル点をスプライン関数で補間して得た関数を実線で示す。図５Ｃでは、比較のために図５Ｂの関数を破線で示している。図５Ｂに示すサンプル点が多い場合の動径分布関数と、図５Ｃに示すサンプル点の少ない場合の動径分布関数とを見れば明らかであり、サンプル点の取り方に応じて得られる動径分布関数が変動してしまう。すると、動径分布関数に基づき算出される非結合ポテンシャルも変動することとなり、粗視化分子動力学計算によって得られる物性が真値から変化してしまう問題が生じる。
【０００８】
本発明は、このような課題に着目してなされたものであって、その目的は、離散値のサンプル点に依存せずに動径分布関数を取得可能にして、信頼性を高めた粗視化モデルの非結合ポテンシャルの算出方法、算出装置及び算出プログラムを提供することである。
【課題を解決するための手段】
【０００９】
本発明は、上記目的を達成するために、次のような手段を講じている。
【００１０】
すなわち、本発明の粗視化モデルの非結合ポテンシャルの算出方法は、原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを、コンピュータが算出する方法であって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して目標となる動径分布関数を取得するステップと、
前記目標となる動径分布関数に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルにパラメータとして設定し、粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して動径分布関数を取得し、粗視化分子動力学に基づく動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下となるか否かを判定し、前記評価関数が前記所定閾値よりも大きいときは、前記評価関数が小さくなる方向に補正した前記非結合ポテンシャルを新たなパラメータとして前記評価関数が前記所定閾値以下となるまで粗視化分子動力学計算を繰り返すステップと、を含み、
前記粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得することを特徴とする。
【００１１】
本発明は、装置としても特定可能である。
すなわち、本発明の粗視化モデルの非結合ポテンシャルの算出装置は、原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを算出する算出装置であって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出する原子モデル分子動力学算出部と、
前記原子モデル分子動力学算出部が算出した離散データを近似して目標となる動径分布関数を取得する第一動径分布関数取得部と、
前記第一動径分布関数取得部が取得した目標となる動径分布関数に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルにパラメータとして設定する非結合ポテンシャル設定部と、
粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出する粗視化モデル分子動力学算出部と、
前記粗視化モデル分子動力学算出部が算出した離散データを近似して動径分布関数を取得する第二動径分布関数取得部と、
前記第二動径分布関数取得部が取得した粗視化分子動力学に基づく動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下となるか否かを判定する判定部と、を備え、
前記評価関数が前記所定閾値よりも大きいと前記判定部が判定した場合には、前記評価関数が小さくなる方向に補正した前記非結合ポテンシャルを新たなパラメータとし、前記評価関数が前記所定閾値以下となるまで前記粗視化モデル分子動力学算出部による粗視化分子動力学計算を繰り返すように構成されており、
前記第一及び第二の動径分布関数取得部は、前記粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得するように構成されていることを特徴とする。
【００１２】
窓関数は、微分可能な関数であり、且つ、定義区間の中心が最大値となり、定義区間の中心に対して対称となる形状をなし、積分可能な関数で、積分値が１となる関数である。
【００１３】
このように、粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算（積分）することにより、動径分布関数を取得するので、従来のようにスプライン補間などの補間法により平滑化処理及び近似処理を行う場合に比べて、離散データのサンプリングの取り方に依存せずに一様に動径分布関数を取得でき、動径分布関数に基づく非結合ポテンシャルの算出精度を向上させることが可能となる。したがって、物性を真値に近づけることができ、より適切なシミュレーションが実施可能になる。
【００１４】
畳み込み演算による近似処理及び平滑化処理の効果を適切に得るためには、前記窓関数は、次の式（１）で表されるガウス関数であることが好ましい。σは標準偏差である。
ｈ（ｕ）＝（１／√（２πσ２））ｅｘｐ（−ｕ２／２σ２）  …（１）
【００１５】
ノイズの残留や関数の絶対値が小さくなる弊害を回避して、適切な動径分布関数の取得を可能にするためには、動径分布関数の基礎となる離散データに基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定するステップと、
上記ガウス関数の標準偏差σの値を種々異ならせて平滑化した動径分布関数を算出するステップと、
各々の算出結果について、平滑化後の動径分布関数をｇ（ｒ）とし、平滑化前の動径分布関数をｇ’（ｒ）とする場合に、０．９×ｒmax≦ｒ≦１．１×ｒmaxの範囲において、下記式（６）に示すＲを算出するステップと、
上記複数の算出結果から、Ｒ≦１０−２を満たすσの最大値σmaxを抽出すると共に、Ｒ≧１０−４を満たすσの最小値σminを抽出し、σmin以上且つσmax以下となるσを決定するステップと、を更に含むことが好ましい。
装置の場合には、上記ステップを実行する標準偏差値決定部を設けるのが好ましい。
【００１６】
適切な動径分布関数の取得を可能にする他の一つの手段としては、前記離散データに基づき波形のピーク又はボトムを特定するステップと、
特定した波形のピーク又はボトムに基づきノイズの波長又は周波数を特定するステップと、
特定したノイズの波長又は周波数に応じてガウス関数のσを決定するステップと、を更に含むことが好ましい。
装置の場合には、上記ステップを実行する標準偏差値決定部を設けるのが好ましい。
【００１７】
本発明は、上記算出方法を構成するステップを、プログラムの観点から特定することも可能である。
【００１８】
すなわち、本発明の粗視化モデルの非結合ポテンシャルの算出プログラムは、原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを算出するプログラムであって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して目標となる動径分布関数を取得するステップと、
前記目標となる動径分布関数に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルにパラメータとして設定し、粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して動径分布関数を取得し、粗視化分子動力学に基づく動径分布関数と上記目標となる動径分布関数との差値に対応する評価関数が所定閾値以下となるか否かを判定し、前記評価関数が前記所定閾値よりも大きいときは、前記評価関数が小さくなる方向に補正した前記非結合ポテンシャルを新たなパラメータとして前記評価関数が前記所定閾値以下となるまで粗視化分子動力学計算を繰り返すステップと、をコンピュータに実行させ、
前記粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得することを特徴とする。
このプログラムを実行することによっても、上記算出方法が奏する作用効果を得ることができる。
【図面の簡単な説明】
【００１９】
【図１】本発明に係る粗視化モデルの非結合ポテンシャルの算出方法を示すフローチャート。
【図２】図１における標準偏差値決定処理ルーチンを示すフローチャート。
【図３】本発明に係る非結合ポテンシャルの算出装置を模式的に示すブロック図。
【図４】本実施形態で用いる原子モデル及び粗視化モデルの一例を模式的に示す説明図。
【図５Ａ】分子動力学により算出される離散データを模式的に示す説明図。
【図５Ｂ】従来の方法により離散データから取得した動径分布関数を示す模式図。
【図５Ｃ】従来の方法により離散データから取得した動径分布関数を示す模式図。
【図６Ａ】本発明の方法により離散データから取得した動径分布関数を示す模式図。
【図６Ｂ】本発明の方法により離散データから取得した動径分布関数を示す模式図。
【図６Ｃ】本発明の方法により離散データから取得した動径分布関数を示す模式図。
【図６Ｄ】本発明の方法により離散データから取得した動径分布関数を示す模式図。
【図６Ｅ】本発明の方法により離散データから取得した動径分布関数を示す模式図。
【図７】本発明の他の実施形態に係る標準偏差値決定処理を示すフローチャート。
【図８】本発明の他の実施形態に係る標準偏差値決定処理を説明する図。
【発明を実施するための形態】
【００２０】
以下、本発明の一実施形態を、図面を参照して説明する。
【００２１】
［粗視化モデルの非結合ポテンシャルの算出装置］
本実施形態の粗視化モデルの非結合ポテンシャルの算出装置は、粗視化分子動力学計算に必要となる粗視化モデルの非結合ポテンシャルを算出する装置である。具体的には、算出装置は、原子モデルの分子動力学に基づき目標となる動径分布関数を取得し、目標となる動径分布関数を粗視化モデルにおいて再現する非結合ポテンシャルを算出する。本実施形態では、図４に例示するように、原子レベルの分子動力学では、５つの原子Ｍ１〜Ｍ５の各々を最小単位として取り扱うのに対し、粗視化モデルでは、例えば５つ原子からなる分子を一括りとした粗視化粒子Ｍ６を最小単位のモデルとして取り扱う。なお、粗視化粒子を構成する原子の数は目的に応じて任意に変更可能である。
【００２２】
図３に示すように、算出装置１は、初期設定部１０と、原子モデル分子動力学算出部１１と、第一動径分布関数取得部１２と、非結合ポテンシャル設定部１３と、粗視化モデル分子動力学算出部１４と、第二動径分布関数取得部１５と、標準偏差値決定部１６と、判定部１７とを有する。これら各部１０〜１７は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図示しない処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。
【００２３】
図３に示す初期設定部１０は、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、解析対象となる原子モデルや粗視化モデルに関する情報の設定、分子動力学計算に必要な設定等の各種設定を実行し、これら設定値を図示しないメモリに記憶する。原子モデルや粗視化モデルは、図４に例示するものが挙げられる。
【００２４】
図３に示す原子モデル分子動力学算出部１１は、予め設定されている原子モデル、その結合ポテンシャル、非結合ポテンシャル等の各種パラメータを用い、原子モデルの分子動力学計算に基づき平衡状態における各原子モデルの重心位置を算出する。これに伴い、粗視化モデルを構成する原子モデルの重心位置に基づき粗視化モデルの重心位置を算出し、粗視化モデルの動径分布関数の基礎となる離散データを算出する。離散データは、図５Ａに示すように、距離ｒ及び当該距離ｒを入力値とする動径分布関数ｇ（ｒ）の値で構成されるデータ（図中にて丸印で示す）が離散的に多数存在し、動径分布関数の基礎となるデータである。なお、図５Ａに示す離散データは、説明のために間引いて表している。
【００２５】
図３に示す第一動径分布関数取得部１２は、原子モデル分子動力学算出部１１が算出した離散データ（図５Ａ参照）に近似し且つ目標となる動径分布関数Ｇ（ｒ）を取得する。ここでは、離散データと予め定めた窓関数ｈ（ｕ）とを畳み込み演算（積分）することにより、目標となる動径分布関数Ｇ（ｒ）を得る。窓関数ｈ（ｕ）は、ガウス関数であり、次の式（１）で表される。
ｈ（ｕ）＝（１／√（２πσ２））ｅｘｐ（−ｕ２／２σ２）  …（１）
本実施形態では、窓関数としてガウス関数を用いているが、これに限定されない。窓関数として、微分可能な関数であり、且つ、定義区間の中心が最大値となり、定義区間の中心に対して対称となる形状をなし、積分可能な関数で、積分値が１となる関数であれば、任意の釣り鐘型などの関数を用いることができる。例えば、ハニング関数やハミング関数などが挙げられる。σは関数の標準偏差を表す。
【００２６】
図３に示す非結合ポテンシャル設定部１３は、第一動径分布関数取得部１２で得られた目標となる動径分布関数Ｇ（ｒ）に基づき下記式（２）を用いて平均力ポテンシャルＵＰＭＦを算出し、これを下記式（３）に示す粗視化モデルの非結合ポテンシャルＵ（０）としてメモリにパラメータとしてセットする。
ＵＰＭＦ（ｒ）＝−ｋＢＴ ln Ｇ（ｒ）  …（２）
Ｕ（０）＝−ｋＢＴ ln Ｇ（ｒ）  …（３）
式（２）のように動径分布関数から平均力ポテンシャルエネルギー関数に変換する処理はBoltzmann Inversionと呼ばれる。Ｕ（０）は初期値を意味する。ｋＢはボルツマン定数を表し、Ｔは温度を表す。
【００２７】
図３に示す粗視化モデル分子動力学算出部１４は、予め設定されている粗視化モデル、パラメータとしての非結合ポテンシャルＵ（ｉ）を用いて粗視化モデルの分子動力学計算に基づき平衡状態における粗視化モデルの重心位置を算出する。これに伴い、粗視化モデルの動径分布関数の基礎となる離散データ（図５Ａ参照）を算出する。ｉは、粗視化モデルの分子動力学計算の繰り返し回数である。
【００２８】
図３に示す第二動径分布関数取得部１５は、粗視化モデル分子動力学算出部１４が算出した離散データ（図５Ａ参照）に近似する粗視化モデルの動径分布関数ｇ（ｉ）（ｒ）を取得する。処理内容は、上記第一動径分布関数取得部１２が行う処理と同じである。
【００２９】
図３に示す標準偏差値決定部１６は、第一及び第二の動径分布関数取得部１２、１５で用いるガウス関数の標準偏差値σを決定する。ガウス関数などの窓関数は、平滑化フィルタの機能を有するが、σの値が小さいほど平滑化の効果が少なく、σの値が大きいほど平滑化の効果が大きくなる。また、動径分布関数には、図６Ａに示すように、高周波数成分のノイズが乗っているが、ガウス関数などの窓関数はいわゆるローパスフィルタの機能を有するため、高周波数成分を除去するための適切な標準偏差値σを設定することが好ましい。図６Ａ〜図６Ｅは、平滑化前の関数をＢｆとし、平滑化後の関数をＡｆとて表している。σは、図６Ａの例が最も小さく、図６Ｂから図６Ｅの順にσが大きくなる。仮に、σが過大であると、図６Ｅに示すように、平滑化後の動径分布関数ｇ（ｒ）の絶対値が小さくなり、上限値σmaxが必要となる。逆に、σが過小であると、図６Ａに示すように、高周波成分のノイズが残るため、下限値σminが必要となる。高周波成分のノイズを除去でき、且つ、関数の絶対値が損なわれないようにするためには、図６Ｂ〜Ｄの範囲が好ましく、１０−４≦Ｒ≦１０−２を満たすσが好ましいことが分かる。
【００３０】
図３に示す標準偏差値決定部１６は、動径分布関数の基礎となる離散データ（図５Ａ参照）に基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定する（図６Ａ参照）。次に、標準偏差値決定部１６は、σの値を種々異ならせて離散データとガウス関数を用いて平滑化後の動径分布関数ｇ（ｒ）を算出する。次に、標準偏差値決定部１６は、各々の算出結果について、平滑化前の動径分布関数をｇ’（ｒ）とした場合に、０．９×ｒmax≦ｒ≦１．１×ｒmaxの範囲において、次の式（６）を算出する。
2014102554.tif 000003
（ただし、Ｎ＝（１．１×ｒmax−０．９×ｒmax）／ｄｒ、Ｎは、四捨五入した整数）
次に、標準偏差値決定部１６は、上記の複数の算出結果から、Ｒ≦１０−２を満たすσの最大値σmaxを抽出すると共に、Ｒ≧１０−４を満たすσの最小値σminを抽出し、σmin以上且つσmax以下となるσを決定する。決定方法は、最小値σminと最大値σmaxの中間値に決定する方法や、最大値σmaxと最小値σminの範囲内においてランダムに決定する方法が挙げられる。
【００３１】
ここで、範囲を０．９×ｒmax≦ｒ≦１．１×ｒmaxにしているのは、動径分布関数が最大となる領域では、粒子が多数存在するために、ノイズが他の領域よりも顕著になるため、当該領域のみに限定した方が効果的だからである。勿論、領域を限定しなくてもよい。なお、本実施形態では、標準偏差値σを逐次決定しているが、予めメモリに設定した所定値をσに用いる場合には、標準偏差値決定部１６を設けなくてもよい。
【００３２】
図３に示す判定部１７は、粗視化モデルの動径分布関数ｇ（ｉ）（ｒ）と、上記目標となる動径分布関数Ｇ（ｒ）との差値に対応する評価関数が所定閾値以下であるか否かを判定する。評価関数は、粗視化モデルと原子モデルの動径分布関数の平均二乗偏差である。判定部は、式（４）を満足するか否かを判定する。Δは所定閾値であり、Ｒはカットオフ長さである。これらは予めメモリに設定されている。
2014102554.tif 000004
【００３３】
図３に示す判定部１７は、評価関数が所定閾値以下のときは、パラメータとして設定されている非結合ポテンシャルＵ（ｉ）が解であるとする。一方、判定部１７は、評価関数が所定閾値よりも大きいときには、評価関数が小さくなる方向に非結合ポテンシャルを補正し、補正した非結合ポテンシャルを新たなパラメータとして粗視化モデルの分子動力学計算を繰り返す。下記式（５）に示すように、補正した非結合ポテンシャルはＵ（ｉ＋１）であり、式（５）の右辺の第２項が補正値となる。ｋＢはボルツマン定数を表し、Ｔは温度を表す。
2014102554.tif 000005
【００３４】
［粗視化モデルの非結合ポテンシャルの算出方法］
上記算出装置１を用いた粗視化モデルの非結合ポテンシャルの算出方法について説明する。
【００３５】
まず、ステップＳＴ１（図１参照）において、図３に示す初期設定部１０は、解析対象となる原子モデルや粗視化モデルに関する情報の設定、分子動力学計算に必要な設定等の各種設定を実行し、これら設定値を図示しないメモリに記憶する。
【００３６】
次のステップＳＴ２（図１参照）において、図３に示す原子モデル分子動力学算出部１１は、予め設定された各種パラメータを用い、原子モデルの分子動力学計算に基づき平衡状態にある各原子モデルの重心位置を算出する。次のステップＳＴ３（図１参照）において、原子モデル分子動力学算出部１１は、粗視化モデルの重心位置を算出し、粗視化モデル間の重心位置に関する動径分布関数の基礎となる離散データを算出する。
【００３７】
次のステップＳＴ４（図１参照）において、図３に示す標準偏差値決定部１６は、図２に示す標準偏差値決定処理ルーチンを実行する。なお、標準偏差値σが予め設定されているときには、この処理を省略できる。
【００３８】
次のステップＳＴ５（図１参照）において、図３に示す第一動径分布関数取得部１２は、原子モデル分子動力学算出部１１が算出した離散データ（図５Ａ参照）に平滑化且つ近似し、図６Ｂ〜Ｄにて符号Ａｆに示すように、目標となる動径分布関数Ｇ（ｒ）を取得する。
【００３９】
次のステップＳＴ６（図１参照）において、図３に示す非結合ポテンシャル設定部１３は、第一動径分布関数取得部１２で得られた目標となる動径分布関数Ｇ（ｒ）に基づき平均力ポテンシャルＵＰＭＦを算出し、これを粗視化モデルの非結合ポテンシャルＵ（０）としてメモリにパラメータとしてセットする。
【００４０】
次のステップＳＴ７（図１参照）において、図３に示す粗視化モデル分子動力学算出部１４は、パラメータとしての非結合ポテンシャルＵ（ｉ）を用いて粗視化モデルの分子動力学計算に基づき平衡状態における粗視化モデルの重心位置を算出する。次のステップＳＴ８（図１参照）において、図３に示す粗視化モデル分子動力学算出部１４は、粗視化モデル間の重心位置に関する動径分布関数の基礎となる離散データ（図５Ａ参照）を算出する。
【００４１】
次のステップＳＴ９（図１参照）において、図３に示す標準偏差値決定部１６は、図２に示す標準偏差値決定処理ルーチンを実行する。なお、標準偏差値σが予め設定されているときには、この処理を省略できる。次のステップＳＴ１０（図１参照）において、図３に示す第二動径分布関数取得部１５は、粗視化モデル分子動力学算出部１４が算出した離散データ（図５Ａ参照）に平滑化且つ近似し、図６Ｂ〜Ｄにて符号Ａｆに示すように、動径分布関数ｇ（ｉ）（ｒ）を取得する。
【００４２】
次のステップＳＴ１１（図１参照）において、図３に示す判定部１７は、粗視化モデルの動径分布関数ｇ（ｉ）（ｒ）と、ステップＳＴ５で算出した目標となる動径分布関数Ｇ（ｒ）との差値に対応する評価関数が所定閾値以下であるか否かを判定する。評価関数が所定閾値以下でないと判定された場合（ＳＴ１１：ＮＯ参照）には、判定部１７は、パラメータである非結合ポテンシャルを補正し（ＳＴ１２参照）、ステップＳＴ７に戻る。一方、評価関数が所定閾値以下であると判定された場合（ＳＴ１１：ＹＥＳ参照）には、パラメータの非結合ポテンシャルＵ（ｉ）を解として特定する。
【００４３】
図２に示すように、標準偏差値決定処理ルーチンが実行されると、標準偏差値決定部１６は、ステップＳＴ２０において、離散データ（図５Ａ参照）に基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定する。次に、標準偏差値決定部１６は、ステップＳＴ２１において、σの値を種々異ならせて離散データとガウス関数を用いて平滑化後の動径分布関数ｇ（ｒ）を複数算出する。次に、標準偏差値決定部１６は、ステップＳＴ２２において、各々の動径分布関数に対してＲ値を算出する。次に、標準偏差値決定部１６は、ステップＳＴ２３において、σの最大値と最小値を算出し、使用するσを特定する。
【００４４】
以上のように、本実施形態の粗視化モデルの非結合ポテンシャルの算出方法は、原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数Ｇ（ｒ）を粗視化モデルにおいて再現する非結合ポテンシャルＵを、コンピュータが算出する方法であって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し、当該離散データを近似して目標となる動径分布関数を取得Ｇ（ｒ）するステップ（ＳＴ２〜５）と、
目標となる動径分布関数Ｇ（ｒ）に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルＵ（０）にパラメータとして設定し（ＳＴ６）、粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出し（ＳＴ７）、当該離散データを近似して動径分布関数ｇ（ｉ）（ｒ）を取得し（ＳＴ１０）、粗視化分子動力学に基づく動径分布関数ｇ（ｉ）（ｒ）と目標となる動径分布関数Ｇ（ｒ）との差値に対応する評価関数が所定閾値以下となるか否かを判定し（ＳＴ１１）、評価関数が所定閾値よりも大きいときは、評価関数が小さくなる方向に補正した非結合ポテンシャルＵ（ｉ＋１）を新たなパラメータとして評価関数が所定閾値以下となるまで粗視化分子動力学計算を繰り返すステップ（ＳＴ６〜ＳＴ１２）と、を含み、
粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数を取得する（ＳＴ５、ＳＴ１０）。
【００４５】
本実施形態の粗視化モデルの非結合ポテンシャルの算出装置１は、原子モデルの分子動力学計算に基づき得られる目標となる動径分布関数Ｇ（ｒ）を粗視化モデルにおいて再現する非結合ポテンシャルＵを算出する算出装置であって、
原子モデルの分子動力学計算に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出する原子モデル分子動力学算出部１１と、
原子モデル分子動力学算出部１１が算出した離散データを近似して目標となる動径分布関数Ｇ（ｒ）を取得する第一動径分布関数取得部１２と、
第一動径分布関数取得部１２が取得した目標となる動径分布関数Ｇ（ｒ）に基づき算出される平均力ポテンシャルを粗視化モデルの非結合ポテンシャルＵ（０）にパラメータとして設定する非結合ポテンシャル設定部１３と、
粗視化モデルの分子動力学に基づき平衡状態にある粗視化モデルの重心位置の動径分布関数の基礎となる離散データを算出する粗視化モデル分子動力学算出部１４と、
粗視化モデル分子動力学算出部１４が算出した離散データを近似して動径分布関数ｇ（ｉ）（ｒ）を取得する第二動径分布関数取得部１５と、
第二動径分布関数取得部１５が取得した粗視化分子動力学に基づく動径分布関数ｇ（ｉ）（ｒ）と目標となる動径分布関数Ｇ（ｒ）との差値に対応する評価関数が所定閾値以下となるか否かを判定する判定部１７と、を備え、
評価関数が所定閾値よりも大きいと判定部１７が判定した場合には、評価関数が小さくなる方向に補正した非結合ポテンシャルＵ（ｉ＋１）を新たなパラメータとし、評価関数が所定閾値以下となるまで粗視化モデル分子動力学算出部１４による粗視化分子動力学計算を繰り返すように構成されており、
第一及び第二の動径分布関数取得部１２、１５は、粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算することにより、平滑化した動径分布関数ｇ（ｉ）（ｒ）を取得するように構成されている。
【００４６】
このように、粗視化モデルの重心位置の動径分布関数の基礎となる離散データと予め定めた窓関数とを畳み込み演算（積分）することにより、動径分布関数ｇ（ｉ）（ｒ）を取得するので、図５Ｂ及び図５Ｃに示す従来のスプライン補間などの補間法により平滑化処理及び近似処理を行う場合に比べて、離散データのサンプリングの取り方に依存せずに一様に動径分布関数を取得でき、動径分布関数ｇ（ｉ）（ｒ）に基づく非結合ポテンシャルＵの算出精度を向上させることが可能となる。したがって、信頼性の高い非結合ポテンシャルを得ることができ、物性を真値に近づけることができ、より適切なシミュレーションが実施可能になる。
【００４７】
本実施形態では、窓関数は、式（１）で表されるガウス関数であるので、他の窓関数に比して平滑化処理及び近似処理の効果を得るうえで、好ましい。
【００４８】
さらに、本実施形態の方法では、動径分布関数の基礎となる離散データに基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定するステップＳＴ２０と、
上記ガウス関数のσの値を種々異ならせて平滑化した動径分布関数を算出するステップＳＴ２１と、
各々の算出結果について、平滑化後の動径分布関数をｇ（ｒ）とし、平滑化前の動径分布関数をｇ’（ｒ）とする場合に、０．９×ｒmax≦ｒ≦１．１×ｒmaxの範囲において、式（６）に示すＲを算出するステップＳＴ２２と、
上記複数の算出結果から、Ｒ≦１０−２を満たすσの最大値σmaxを抽出すると共に、Ｒ≧１０−４を満たすσの最小値σminを抽出し、σmin以上且つσmax以下となるσを決定するステップＳＴ２３と、を更に含む。
【００４９】
また、本実施形態の装置では、動径分布関数の基礎となる離散データに基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定し、上記ガウス関数のσの値を種々異ならせて平滑化した動径分布関数を算出し、各々の算出結果について、平滑化後の動径分布関数をｇ（ｒ）とし、平滑化前の動径分布関数をｇ’（ｒ）とする場合に、０．９×ｒmax≦ｒ≦１．１×ｒmaxの範囲において、式（６）に示すＲを算出し、上記複数の算出結果から、Ｒ≦１０−２を満たすσの最大値σmaxを抽出すると共に、Ｒ≧１０−４を満たすσの最小値σminを抽出し、σmin以上且つσmax以下となるσを決定する標準偏差値決定部を更に備える。
【００５０】
このように、σの最大値及び最小値を特定し、両者の間にσの値を設定しているので、高周波成分のノイズが残ったり、動径分布関数の絶対値が小さくなりすぎたりする弊害を回避して、適切に動径分布関数を取得することが可能となる。したがって、信頼性の高い非結合ポテンシャルを得ることが可能となる。
【００５１】
本実施形態に係る粗視化モデルの非結合ポテンシャルの算出プログラムは、上記算出方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記算出方法の奏する作用効果を得ることが可能となる。言い換えると、上記算出方法を使用しているとも言える。
【００５２】
以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【００５３】
例えば、図３に示す各部１０〜１７は、所定プログラムをコンピュータのＣＰＵで実行することで実現しているが、各部を専用メモリや専用回路で構成してもよい。
【００５４】
［他の実施形態］
また、標準偏差値決定部１６を、離散データに基づき高周波成分のノイズの波長又は周波数を特定し、特定した波長又は周波数に応じてガウス関数のσの値を決定するように構成してもよい。具体的には、図７及び図８に示すように、動径分布関数の基礎となる離散データに基づき当該動径分布関数が最大となる入力パラメータである最大距離値ｒmaxを特定する（ＳＴ３０参照）。ｒ＝ｒmaxの所定周辺（関数値の最大値付近）において、波形のピークｐｔ又はボトムｂｔを特定する（ＳＴ３１参照）。次に、特定した波形のピークｐｔ又はボトムｂｔに基づきノイズの波長又は周波数を特定する（ＳＴ３２参照）。具体的には、図８に示すように所定数（複数）のノイズ波を特定し、これらの波長λ１〜λ３、‥の算術平均をノイズの波長λavgとする。次に、ノイズの波長λavg又は周波数に応じてガウス関数のσを決定する（ＳＴ３３参照）。具体的には、σ＝ｆ（λavg）のように、波長又は周波数をσに適切に変換するための予め設定した関数を用いることが一例として挙げられる。上記において、波形のピークｐｔ又はボトムｂｔを特定する際に、ｒmaxの周辺に限定しているが、これは、動径分布関数が最大となる領域では、粒子が多数存在するために、ノイズが他の領域よりも顕著になるため、当該領域のみに限定した方が効果的だからである。勿論、領域を限定しなくてもよい。この場合、最大距離値ｒmaxを特定するステップを省略可能である。
【００５５】
このように構成すれば、動径分布関数の基礎となる離散データに含まれるノイズ成分の波長又は周波数を特定して、ノイズを適切に除去するとともに、動径分布関数の値が損なわれることを抑制でき、適切に動径分布関数を取得することが可能となる。したがって、信頼性の高い非結合ポテンシャルを得ることが可能となる。
【００５６】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【００５７】
Ｇ（ｒ）…目標となる動径分布関数（原子モデルに基づく動径分布関数）
Ｕ（ｉ）…非結合ポテンシャル
ｇ（ｉ）（ｒ）…粗視化モデルに基づく動径分布関数
【図１】
2014102554.tif 000008
【図２】
2014102554.tif 000009
【図３】
2014102554.tif 000010
【図４】
2014102554.tif 000011
【図５Ａ】
2014102554.tif 000012
【図５Ｂ】
2014102554.tif 000013
【図５Ｃ】
2014102554.tif 000014
【図６Ａ】
2014102554.tif 000015
【図６Ｂ】
2014102554.tif 000016
【図６Ｃ】
2014102554.tif 000017
【図６Ｄ】
2014102554.tif 000018
【図６Ｅ】
2014102554.tif 000019
【図７】
2014102554.tif 000020
【図８】
2014102554.tif 000021