(57)【要約】【課題】粒子数が大きい系に対しても適切な計算時間で計算することが可能な分子動力学シミュレーション方法およびプログラムを提供する。【解決手段】各粒子に対して粒子間の相互作用力が働く探索範囲内の粒子を用いて各粒子との粒子ペアをリスト化するリスト化ステップと、リスト化した粒子ペアについて、少なくともポテンシャル関数に含まれるペア配列関数の値を計算し、計算結果を２次元配列変数とするペア配列関数計算ステップと、ペア配列関数計算ステップの計算結果を用いて系内の粒子間の相互作用力を計算する相互作用力計算ステップと、相互作用力計算ステップで計算した結果を用いて最新の時刻における各粒子の速度を計算する速度計算ステップと、相互作用力計算ステップで計算した結果を用いて最新の時刻における位置座標を計算する座標計算ステップと、を含む。【選択図】図４
2013211006.tif 000002
【特許請求の範囲】
【請求項１】
複数の粒子からなる系における各粒子の速度、位置座標およびポテンシャル関数が有するポテンシャルパラメータを含む複数のパラメータを記憶する記憶部を備えたコンピュータが実行する分子動力学シミュレーション方法であって、
各粒子に対して粒子間の相互作用力が働く探索範囲内の粒子を用いて各粒子との粒子ペアをリスト化し、該リスト化した粒子ペアを前記記憶部へ格納するリスト化ステップと、
前記リスト化ステップでリスト化した粒子ペアについて、少なくとも前記ポテンシャル関数に含まれるペア配列関数の値を計算し、計算結果を２次元配列変数として前記記憶部へ格納するペア配列関数計算ステップと、
前記ペア配列関数計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて前記系内の粒子間の相互作用力を計算して前記記憶部へ格納する相互作用力計算ステップと、
前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における各粒子の速度を計算して前記記憶部へ格納する速度計算ステップと、
前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における位置座標を計算して前記記憶部へ格納する座標計算ステップと、
を含むことを特徴とする分子動力学シミュレーション方法。
【請求項２】
前記ペア配列関数は、２つの粒子間の距離に応じて該２つの粒子間の相互作用力を打ち切る効果を表すカットオフ関数および２つの粒子の近傍に第３の粒子が存在する場合に該第３の粒子が２つの粒子間の相互作用を遮蔽する効果を表すスクリーン関数の少なくともいずれか一方を含むことを特徴とする請求項１に記載の分子動力学シミュレーション方法。
【請求項３】
前記２次元配列変数は、前記系の全粒子数をＮとし、該Ｎより小さい自然数をＸとするとき、Ｘ×Ｎのマトリックスの各要素に変数が格納されたものであることを特徴とする請求項１または２に記載の分子動力学シミュレーション方法。
【請求項４】
前記自然数Ｘは、前記リスト化ステップでリスト化された任意の粒子とペアを構成する粒子数の前記系内の全ての粒子に関する最大値であることを特徴とする請求項３に記載の分子動力学シミュレーション方法。
【請求項５】
前記リスト化ステップは、
前記探索範囲の外に位置する粒子が前記ペア配列関数がゼロでない領域に入ってくるまでに要する時間よりも短い周期で行うことを特徴とする請求項１〜４のいずれか一項に記載の分子動力学シミュレーション方法。
【請求項６】
前記リスト化ステップから前記座標計算ステップに至る処理を規定回数繰り返すことを特徴とする請求項１〜５のいずれか一項に記載の分子動力学シミュレーション方法。
【請求項７】
前記複数のパラメータは、前記系内の粒子の種類、該粒子の質量、前記系の設定温度、およびシミュレーションの時間刻みを含み、
前記リスト化ステップは、リスト化を行う際に前記記憶部から前記粒子の位置座標を呼び出し、
前記ペア配列関数計算ステップは、計算の際に前記記憶部から前記粒子の種類および前記位置座標を呼び出し、
前記相互作用力計算ステップは、計算の際に前記記憶部から前記粒子の種類、前記位置座標および前記ポテンシャルパラメータを呼び出し、
前記速度計算ステップは、計算の際に前記記憶部から前記粒子の種類、前記粒子の質量、前記速度および前記時間刻みを呼び出し、
前記座標計算ステップは、計算の際に前記記憶部から前記粒子の種類、前記粒子の質量、前記位置座標、前記速度および前記時間刻みを呼び出すことを特徴とする請求項１〜６のいずれか一項に記載の分子動力学シミュレーション方法。
【請求項８】
複数の粒子からなる系における各粒子の速度、位置座標およびポテンシャル関数が有するポテンシャルパラメータを含む複数のパラメータを記憶する記憶部を備えたコンピュータに、
各粒子に対して粒子間の相互作用力が働く探索範囲内の粒子を用いて各粒子との粒子ペアをリスト化し、該リスト化したデータを前記記憶部へ格納するリスト化ステップと、
前記リスト化ステップでリスト化した粒子ペアについて、少なくとも前記ポテンシャル関数に含まれるペア配列関数の値を計算し、計算結果を２次元配列変数として前記記憶部へ格納するペア配列関数計算ステップと、
前記ペア配列関数計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて前記系内の粒子間の相互作用力を計算して前記記憶部へ格納する相互作用力計算ステップと、
前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における各粒子の速度を計算して前記記憶部へ格納する速度計算ステップと、
前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における位置座標を計算して前記記憶部へ格納する座標計算ステップと、
を実行させることを特徴とする分子動力学シミュレーションプログラム。
【請求項９】
前記ペア配列関数は、２つの粒子間の距離に応じて該２つの粒子間の相互作用力を打ち切る効果を表すカットオフ関数および２つの粒子の近傍に第３の粒子が存在する場合に該第３の粒子が２つの粒子間の相互作用を遮蔽する効果を表すスクリーン関数の少なくともいずれか一方を含むことを特徴とする請求項８に記載の分子動力学シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学法に基づいた計算を行う分子動力学シミュレーション方法およびプログラムに関する。
【背景技術】
【０００２】
分子動力学法は、電子計算機（コンピュータ）に原子、分子などの粒子の座標と速度を入力し、粒子間の相互作用力を所定のポテンシャル関数を用いて算出し、ニュートンの運動法則に基づいて粒子の座標と速度の時間的な変化を求める方法である。本方法の適用対象として、例えばグラフェンやダイヤモンドなどのカーボン材料、シリコンなどの半導体材料のように共有結合を含む材料を挙げることができる。これらの材料に適したポテンシャル関数として、粒子間の結合角および周囲の粒子の配位に依存するボンドオーダー（結合次数）を含むＴｅｒｓｏｆｆポテンシャル、Ｂｒｅｎｎｅｒポテンシャル、ＲＥＢＯ（Reactive  Empirical  Bond-Order）ポテンシャルなどが知られている。
【０００３】
通常、ポテンシャル関数は粒子間距離の関数であるが、現実の系では粒子間距離がある値以上になると粒子間に相互作用力（ポテンシャル関数の一次微分）が働かなくなる。このため、分子動力学法においては、ある値以上の粒子間距離におけるポテンシャルをうち切って０とするカットオフを施すのが一般的である。この場合、カットオフの対象となる粒子間距離の最小値であるカットオフ距離以上では、粒子間に働く相互作用力を考慮する必要がない。したがって、粒子間の相互作用力をすべての粒子間で計算しなくて済む。具体的には、粒子毎にカットオフ距離より少し大きい領域内（以下、探索範囲と称する）にある粒子の番号をリスト化してコンピュータのメモリに予め格納しておき、相互作用力を計算する際にメモリからリストを呼び出して、リストに含まれる粒子との間の相互作用力のみを計算する。このようにリスト化した要素を用いて所定の計算を行う方法は、帳簿法またはブックキーピング法と呼ばれている。
【０００４】
カットオフ距離以上では、ポテンシャル関数およびその一次微分である相互作用力はいずれも０となるが、そのままではカットオフ距離においてポテンシャル関数や相互作用力が不連続になる場合がある。この場合、ポテンシャル関数および相互作用力がカットオフ距離において連続になるように、ポテンシャル関数を定義する必要がある。例えば、上述したＴｅｒｓｏｆｆポテンシャル、Ｂｒｅｎｎｅｒポテンシャル、ＲＥＢＯポテンシャル等の共有結合型のポテンシャル関数は、２つの粒子の粒子間距離の関数であるカットオフ関数ｆC（ｒij）を陽に含む。ここで、変数ｒijは、粒子ｉと粒子ｊの粒子間距離である。カットオフ関数ｆC（ｒij）は、粒子間距離ｒijが所定の範囲で大きくなるにつれて１から０に連続的に変化する関数である。このようなカットオフ関数ｆC（ｒij）は、粒子数がＮの計算モデルにおいては、粒子ｉと粒子ｊの組み合わせに応じてＮ（Ｎ−１）通りの値を有する。
【０００５】
カットオフ関数ｆC（ｒij）を陽に含む共有結合型ポテンシャルには、３体力あるいは４体力効果を含むものがある。このため、例えば３体力効果を含むポテンシャルを用いて粒子ｉと粒子ｊとの間に働く相互作用力を計算する場合に、粒子ｉと粒子ｊ以外の粒子ｋ（≠ｉ，ｊ）による効果を考慮する必要がある。粒子ｉと粒子ｊの相互作用力を計算するには、粒子ｉと粒子ｋの粒子間距離ｒikのカットオフ関数ｆC（ｒik）およびその一次微分の値が必要になる。そこで、通常は相互作用力の計算を行う前に、予めカットオフ関数ｆC（ｒij）およびその一次微分を計算してコンピュータのメモリに格納し、相互作用力を計算する際に必要な値をメモリから呼び出す方法が採られている。
【０００６】
分子動力学法では、時間とともに粒子の座標が変化するため、各粒子についてリスト化された粒子の番号も一定であるとは限らない。そのため、従来は、粒子数Ｎの系における粒子ｉと粒子ｊの全ての組み合わせに対してカットオフ関数ｆC（ｒij）およびその一次微分の値を計算し、この計算した値を格納するためにＮ×Ｎ＝Ｎ2のサイズの２次元配列を用いてメモリに格納していた（例えば、非特許文献１参照）。
【０００７】
しかしながら、上述した従来の方法では、粒子数が多くなると、カットオフ関数の値を格納するために必要とするメモリ容量がＮ2に比例して増加するため、比較的少ない粒子数でもメモリ不足となってしまい、粒子数の大きい大規模計算を行うことは困難であった。特に、ファンデルワールス力のような長距離力を扱う場合には、探索範囲が大きくなってさらに多くの計算時間を必要とするため、計算可能な粒子数はたかだか１万程度であった（例えば、非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】泉聡志、分子動力学冬季講習、ポテンシャル３  共有結合性材料、[online]、２００３年１１月２日、[平成２４年２月２４日検索]、インターネット<URL:  http://www.fml.t.u-tokyo.ac.jp/~izumi/sotsuron/MD＿Potential3.ppt>
【非特許文献２】Bohayra  Mortazavi,  Yves  Remond,  Said  Ahzi,  and  Valerie  Toniazzo,  “Thickness  and  chirality  effects  on  tensile  behavior  of  few-layer  graphene  by  molecular  dynamics  simulations”,  Computational  Materials  Science,  53,  298〜302  (2012).
【発明の開示】
【発明が解決しようとする課題】
【０００９】
一般に、シミュレーションの計算結果の妥当性を検証するために、問題としている特性について実験結果と比較することが行われることがある。例えば、材料の強度を実験的に測定する場合、試験片の把持や応力の検出限界等の測定上の制約から、通常はミクロン／サブミクロンオーダのサイズを有する試験片を扱うことが多い。実験結果との比較を厳密に行うためには、試験片を構成する粒子数と同程度の粒子数の計算モデルが必要となる。このため、ミクロン／サブミクロンオーダのサイズを有する試験片の場合、粒子数は億単位の数となる。しかしながら、上述した従来技術で扱える粒子数は、たかだかその１万分の１（１万）程度に過ぎなかった。このため、分子動力学法において、適切な計算時間で計算可能である粒子数の拡大が待望されているのが現状である。
【００１０】
本発明は、上記に鑑みてなされたものであって、粒子数が大きい系に対しても適切な計算時間で計算することが可能な分子動力学シミュレーション方法およびプログラムを提供することにある。
【課題を解決するための手段】
【００１１】
上述した課題を解決し、目的を達成するために、本発明に係る分子動力学シミュレーション方法は、複数の粒子からなる系における各粒子の速度、位置座標およびポテンシャル関数が有するポテンシャルパラメータを含む複数のパラメータを記憶する記憶部を備えたコンピュータが実行する分子動力学シミュレーション方法であって、各粒子に対して粒子間の相互作用力が働く探索範囲内の粒子を用いて各粒子との粒子ペアをリスト化し、該リスト化した粒子ペアを前記記憶部へ格納するリスト化ステップと、前記リスト化ステップでリスト化した粒子ペアについて、少なくとも前記ポテンシャル関数に含まれるペア配列関数の値を計算し、計算結果を２次元配列変数として前記記憶部へ格納するペア配列関数計算ステップと、前記ペア配列関数計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて前記系内の粒子間の相互作用力を計算して前記記憶部へ格納する相互作用力計算ステップと、前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における各粒子の速度を計算して前記記憶部へ格納する速度計算ステップと、前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における位置座標を計算して前記記憶部へ格納する座標計算ステップと、を含むことを特徴とする。
【００１２】
本発明に係る分子動力学シミュレーション方法は、上記発明において、前記ペア配列関数は、２つの粒子間の距離に応じて該２つの粒子間の相互作用力を打ち切る効果を表すカットオフ関数および２つの粒子の近傍に第３の粒子が存在する場合に該第３の粒子が２つの粒子間の相互作用を遮蔽する効果を表すスクリーン関数の少なくともいずれか一方を含むことを特徴とする。
【００１３】
本発明に係る分子動力学シミュレーション方法は、上記発明において、前記２次元配列変数は、前記系の全粒子数をＮとし、該Ｎより小さい自然数をＸとするとき、Ｘ×Ｎのマトリックスの各要素に変数が格納されたものであることを特徴とする。
【００１４】
本発明に係る分子動力学シミュレーション方法は、上記発明において、前記自然数Ｘは、前記リスト化ステップでリスト化された任意の粒子とペアを構成する粒子数の前記系内の全ての粒子に関する最大値であることを特徴とする。
【００１５】
本発明に係る分子動力学シミュレーション方法は、上記発明において、前記リスト化ステップは、前記探索範囲の外に位置する粒子が前記ペア配列関数がゼロでない領域に入ってくるまでに要する時間よりも短い周期で行うことを特徴とする。
【００１６】
本発明に係る分子動力学シミュレーション方法は、上記発明において、前記リスト化ステップから前記座標計算ステップに至る処理を規定回数繰り返すことを特徴とする。
【００１７】
本発明に係る分子動力学シミュレーション方法は、上記発明において、前記複数のパラメータは、前記系内の粒子の種類、該粒子の質量、前記系の設定温度、およびシミュレーションの時間刻みを含み、前記リスト化ステップは、リスト化を行う際に前記記憶部から前記粒子の位置座標を呼び出し、前記ペア配列関数計算ステップは、計算の際に前記記憶部から前記粒子の種類および前記位置座標を呼び出し、前記相互作用力計算ステップは、計算の際に前記記憶部から前記粒子の種類、前記位置座標および前記ポテンシャルパラメータを呼び出し、前記速度計算ステップは、計算の際に前記記憶部から前記粒子の種類、前記粒子の質量、前記速度および前記時間刻みを呼び出し、前記座標計算ステップは、計算の際に前記記憶部から前記粒子の種類、前記粒子の質量、前記位置座標、前記速度および前記時間刻みを呼び出すことを特徴とする。
【００１８】
本発明に係る分子動力学シミュレーションプログラムは、複数の粒子からなる系における各粒子の速度、位置座標およびポテンシャル関数が有するポテンシャルパラメータを含む複数のパラメータを記憶する記憶部を備えたコンピュータに、各粒子に対して粒子間の相互作用力が働く探索範囲内の粒子を用いて各粒子との粒子ペアをリスト化し、該リスト化したデータを前記記憶部へ格納するリスト化ステップと、前記リスト化ステップでリスト化した粒子ペアについて、少なくとも前記ポテンシャル関数に含まれるペア配列関数の値を計算し、計算結果を２次元配列変数として前記記憶部へ格納するペア配列関数計算ステップと、前記ペア配列関数計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて前記系内の粒子間の相互作用力を計算して前記記憶部へ格納する相互作用力計算ステップと、前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における各粒子の速度を計算して前記記憶部へ格納する速度計算ステップと、前記相互作用力計算ステップで計算した結果を前記記憶部から呼び出し、この計算結果を用いて最新の時刻における位置座標を計算して前記記憶部へ格納する座標計算ステップと、を実行させることを特徴とする。
【００１９】
本発明に係る分子動力学シミュレーションプログラムは、上記発明において、前記ペア配列関数は、２つの粒子間の距離に応じて該２つの粒子間の相互作用力を打ち切る効果を表すカットオフ関数および２つの粒子の近傍に第３の粒子が存在する場合に該第３の粒子が２つの粒子間の相互作用を遮蔽する効果を表すスクリーン関数の少なくともいずれか一方を含むことを特徴とする。
【発明の効果】
【００２０】
本発明によれば、ペア配列関数の計算を予め登録した粒子ペアに対してのみ実行するため、粒子数の大きい系に対しても適切な計算時間で大規模計算を実行することが可能となる。
【図面の簡単な説明】
【００２１】
【図１】図１は、本発明の実施の形態に係る分子動力学シミュレーション方法を実行する装置の機能構成を示すブロック図である。
【図２】図２は、本発明の実施の形態１に係るリスト化処理における粒子の探索範囲を模式的に示す図である。
【図３】図３は、カットオフ関数の一例を示す図である。
【図４】図４は、本発明の実施の形態に係る分子動力学シミュレーション方法の概要を示すフローチャートである。
【図５】図５は、本発明の実施の形態に係る分子動力学シミュレーション方法で使用する２次元配列変数の構成例を示す図である。
【図６】図６は、従来の分子動力学シミュレーション方法で使用する２次元配列変数の構成例を示す図である。
【図７】図７は、本発明の実施の形態２に係るリスト化処理における粒子の探索範囲を模式的に示す図である。
【図８】図８は、本発明の実施例において、ペア配列関数としてカットオフ関数を適用し、計算が完了した場合の全粒子数と計算時間との関係を両対数グラフで示した図である。
【図９】図９は、本発明の実施例において、ペア配列関数としてスクリーン関数を適用し、計算が完了した場合の全粒子数と計算時間との関係を両対数グラフで示した図である。
【発明を実施するための形態】
【００２２】
以下、添付図面を参照して、本発明を実施するための形態（以下、「実施の形態」という）を詳細に説明する。
【００２３】
（実施の形態１）
図１は、本発明の実施の形態１に係る分子動力学シミュレーション方法（以下、ＭＤシミュレーション方法という）を実行可能な装置（以下、ＭＤシミュレーション装置という）の機能構成を示すブロック図である。同図に示すＭＤシミュレーション装置１は、キーボード、マウス、タッチパネル等のインタフェースを用いて実現される入力部２と、液晶または有機ＥＬ等からなる表示パネルを有し、各種情報を表示する表示部３と、ＭＤシミュレーションを行うための情報を含む各種情報を記憶する記憶部４と、ＭＤシミュレーション装置１の動作を制御する制御部５と、を備える。
【００２４】
記憶部４は、ＭＤシミュレーションを行うためのパラメータを記憶するパラメータ記憶部４１と、ＭＤシミュレーションにおける各種計算結果を記憶する計算データ記憶部４２と、本実施の形態１に係る分子動力学シミュレーションプログラムや所定のＯＳを起動するプログラム等を記憶するプログラム記憶部４３とを有する。記憶部４は、ＲＡＭ（メモリ）、ＲＯＭ等を用いて構成される。
【００２５】
パラメータ記憶部４１は、ＭＤシミュレーションの計算対象である系に含まれる粒子の種類、質量、位置座標、速度、系の設定温度、系のポテンシャル関数を与えるポテンシャルパラメータ、シミュレーションの時間刻み等を記憶する。
【００２６】
以下、パラメータについて具体的に説明する。
粒子の種類は系に応じて定まるが、１種類とは限らない。例えば、炭素原子と窒素原子、水分子と酸素分子、というように、異なる種類の粒子が混合する場合も可能である。
粒子の質量は、系に複数の粒子が存在する場合、各粒子の質量である。
粒子の位置座標は、ＭＤシミュレーション装置１が読み取り可能な形式で記述されていればよい。
系の設定温度は、実験結果との比較に際して温度一定の条件を課す場合に用いられる。
ポテンシャルパラメータは、各粒子間に働く相互作用力を求めるためのポテンシャル関数に含まれる。系に異なる種類の粒子が混在する場合には、それぞれのポテンシャルパラメータが記憶される。
シミュレーションの時間刻みは、結晶を構成している粒子が熱エネルギーによって振動している周期（例えばピコ秒以下）の１/１００〜１/２００に設定するのが好ましい。
【００２７】
制御部５は、ＣＰＵ等を用いて実現され、制御対象であるＭＤシミュレーション装置１の各構成部位とバスラインを介して接続されている。制御部５は、リスト化部５１と、ペア配列関数計算部５２と、相互作用力計算部５３と、速度計算部５４と、座標計算部５５と、を有する。
【００２８】
リスト化部５１は、パラメータ記憶部４１から位置座標を呼び出し、この位置座標を用いて探索範囲内の粒子ペアをリスト化し、該リスト化したデータを計算データ記憶部４２へ格納する。図２は、リスト化処理における粒子の探索範囲を模式的に示す図である。図２に示すように、探索範囲Ｄsは、注目粒子Ｐiを中心とする半径（以下、探索半径という）Ｒsの球形領域である。探索範囲Ｄsは、ポテンシャル関数のカットオフ領域Ｄcよりも少し広い範囲とするのが好ましい。すなわち、カットオフ領域Ｄcを、注目粒子Ｐiを中心とする半径（以下、カットオフ半径という）Ｒcの球形領域とするとき、Ｒs＞Ｒcである。なお、図２では、探索範囲Ｄsに含まれる粒子Ｐjを他の粒子と区別するために斜線を用いて記載している。
【００２９】
ペア配列関数計算部５２は、パラメータ記憶部４１から粒子の種類、位置座標およびポテンシャルパラメータを呼び出し、このパラメータを用いて、リスト化された粒子ペアについて、ポテンシャル関数に含まれるペア配列関数として、２つの粒子間の距離に応じて相互作用力を打ち切る効果を表すカットオフ関数を適用する。具体的には、ペア配列関数計算部５２は、カットオフ関数および該カットオフ関数の一次微分の値を計算し、計算結果を２次元配列変数として計算データ記憶部４２へ格納する。
【００３０】
図３は、カットオフ関数の一例を示す図である。カットオフ関数１０１は、粒子間距離ｒの関数である。カットオフ関数１０１は、ｒ＜Ｒminの範囲で一定値１をとり、Ｒmin＜ｒ＜Ｒmaxの範囲でｒが大きくなるとともに滑らかに減少し、ｒ＜Ｒmaxで０となる。したがって、図３におけるＲmaxは、図２におけるカットオフ半径Ｒcに等しい。このような関数として、例えば三角関数、指数関数、多項式関数、対数関数あるいはこれらの有理関数を適用することができる。
【００３１】
カットオフ関数１０１を含むポテンシャル関数を用いて材料の特性を評価する場合、カットオフ関数１０１の値が１から減少し始める粒子間距離Ｒmin近傍でのポテンシャルが重要になる場合が多い。このため、ポテンシャル関数は、関数に含まれる係数の変更によってポテンシャル形状の調整が容易な多項式関数であればより好ましい。また、自然現象を模擬するための計算式は、連続関数によって記述されることが好ましい。このため、カットオフ関数１０１は、相互作用力がｒ＝Ｒmin，Ｒmaxで滑らかになるように、Ｒmin，Ｒmaxにおける連続性も含めて、Ｎを２以上の整数としてＮ回微分可能であり、かつ０次微分からＮ次微分まで全て連続な関数であることが好ましい。
【００３２】
相互作用力計算部５３は、パラメータ記憶部４１から粒子の種類、位置座標およびポテンシャルパラメータを呼び出し、このパラメータを用いて系内の粒子間の相互作用力を計算して計算データ記憶部４２へ格納する。
【００３３】
速度計算部５４は、パラメータ記憶部４１から粒子の種類、質量、位置座標、設定温度、速度、時間刻みを呼び出すとともに、計算データ記憶部４２から相互作用力を呼び出し、このパラメータおよび相互作用力を用いて、時間刻み分だけ進めた時刻における速度を計算してパラメータ記憶部４１が記憶する速度を更新する。
【００３４】
座標計算部５５は、パラメータ記憶部４１から粒子の種類、質量、位置座標、設定温度、速度、時間刻みを呼び出すとともに、計算データ記憶部４２から相互作用力を呼び出し、このパラメータおよび相互作用力を用いて、時間刻み分だけ進めた時刻における座標を計算してパラメータ記憶部４１が記憶する座標を更新する。
【００３５】
以上の構成を有するＭＤシミュレーション装置１は、１つのコンピュータを用いて実現してもよいし、複数のコンピュータを用いて実現してもよい。後者の場合には、インターネット等の通信ネットワークを介してデータの送受信を行いながら、互いに連携して処理を行うようにすることも可能である。
【００３６】
本実施の形態１に係る分子動力学シミュレーションプログラムは、ハードディスク、フラッシュメモリ、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体に記録して広く流通させることも可能である。
【００３７】
図４は、本実施の形態１に係る分子動力学シミュレーション方法の概要を示すフローチャートである。まず、入力部２によって入力された粒子の種類、質量、位置座標、設定温度、速度、時間刻み、ポテンシャルパラメータ等のパラメータをパラメータ記憶部４１に格納する（ステップＳ１）。このステップＳ１は、後述するステップＳ２〜Ｓ５とは別に、予め行っておけばよい。
【００３８】
続くステップＳ２〜Ｓ６の処理は、ＭＤシミュレーション装置１が規定回数だけ繰り返し実行する処理である。以下、ステップＳ２〜Ｓ６の処理を具体的に説明する。
【００３９】
まず、リスト化部５１は、パラメータ記憶部４１が記憶する粒子の位置座標を用いて粒子間の距離を計算し、探索範囲内の粒子ペアをリスト化して計算データ記憶部４２へ格納する（ステップＳ２）。
【００４０】
リスト化の頻度すなわちリストの更新周期は、ある時点で探索範囲Ｄsの外に位置する粒子がカットオフ領域Ｄcの中に入ってくるまでに要する時間よりも短い時間で設定される。この意味で、ステップＳ２のリスト化処理は、ループごとに実行されるとは限らない。ここで、具体的なリストの更新頻度（リストの有効期限）の決定方法を説明する。前提として、注目粒子Ｐiから見て探索半径Ｒsとカットオフ半径Ｒcの差Ｒs−Ｒcの部分の球殻状領域Ｄs-cを粒子が横切る間、リストは有効であり、リストに登録されていない粒子ペアが相互作用する可能性はないものとする。この前提に立った上で、リスト化した時刻をｔとし、系の中で最速の粒子の速度をｖmaxとすると、リストの有効期限ｔeは、
【数１】
2013211006.tif 000003
で与えられる。式（１）の右辺の第１項の分母が２ｖmaxとなっているのは、同程度の速度をもつ粒子が互いに近づいていく状況を考慮に入れているためである。
【００４１】
リストの有効期限ｔe内ではリストを再構築する必要はないが、系の最高速度ｖmaxが変わるたびに有効期限ｔeを修正する必要がある。ある時点までの最高速度をｖoldとし、最新の最高速度をｖnewとすると、新しい有効期限ｔe’は、
【数２】
2013211006.tif 000004
で表される。式（２）から容易に分かるように、最新の最高速度ｖnewがそれ以前の最高速度ｖoldよりも速ければ（ｖnew＞ｖold）、有効期限は短くなる（ｔe’＜ｔe）一方、最新の最高速度ｖnewがそれ以前の最高速度ｖoldよりも遅ければ（ｖnew＜ｖold）、有効期限は長くなる（ｔe’＞ｔe）。また、探索半径Ｒsを大きくするほどリストの有効期限は長くなるため、登録する粒子が多くなり、相互作用力の計算において、実際には相互作用していない粒子間の距離を計算することが多くなる。リストの有効期限は探索半径Ｒsに比例する一方、探索に要する時間は探索半径Ｒsの３乗（Ｒs3）に比例するため、探索半径Ｒsには適切な値が存在することになる。
【００４２】
ステップＳ３において、ペア配列関数計算部５２は、パラメータ記憶部４１から粒子の種類、位置座標、ポテンシャルパラメータを呼び出すとともに、計算データ記憶部４２からリスト化した粒子ペアを呼び出し、これらを用いてペア配列関数であるカットオフ関数およびその一次微分の値を計算し、計算結果を２次元配列変数として計算データ記憶部４２へ格納する（ステップＳ３）。
【００４３】
図５は、ステップＳ３で用いる２次元配列変数の構成例を示す図である。同図に示す２次元配列変数２０１は、サイズがＸ×Ｎのマトリックスの各要素に変数が格納されたものである。マトリックスの各要素に格納される変数Ａ（ｉ，ｊ）（ｉ＝１〜Ｘ、ｊ＝１〜Ｎ）には、カットオフ関数およびその一次微分の値のいずれかが対応付けられる。２次元配列変数２０１の列数Ｎは系の全粒子数であり、行数ＸはＮより小さい自然数である。Ｘの値は、例えば任意の粒子ｉと相互作用する他の粒子の数（以下、Ｎiとする）の全てのｉに関する最大値とすることができる。このようなＸをとる場合には、ＭＤシミュレーション装置１の計算量を低減することができる。なお、Ｘの値は定数としてもよいし、時間刻みごとに変化する変数としてもよい。
【００４４】
ステップＳ４において、相互作用力計算部５３は、パラメータ記憶部４１から粒子の種類、位置座標、ポテンシャルパラメータを呼び出し、ポテンシャル関数によって各粒子に働く相互作用力を計算して計算データ記憶部４２に格納する（ステップＳ４）。このステップＳ４で適用するポテンシャル関数はカットオフ関数を含むポテンシャルであればよく、例えばＴｅｒｓｏｆｆポテンシャル、Ｂｒｅｎｎｅｒポテンシャル、ＲＥＢＯポテンシャル、ＡＩＲＥＢＯ（Adaptive  Intermolecular  REBO）ポテンシャル、ＬＣＢＯＰ（Long-range  Carbon  Bond-Order  Potential）、ＥＤＩＰ（Environment-Dependent  Interaction  Potential）などを挙げることができる。
【００４５】
続いて、速度計算部５４は、パラメータ記憶部４１から粒子の種類、質量、位置座標、設定温度、速度、時間刻みを呼び出すとともに、計算データ記憶部４２から相互作用力を呼び出し、時間刻み分だけ進めた新しい時刻における速度を計算する（ステップＳ５）。このステップＳ５において、速度計算部５４は、数値積分法および温度制御法を適用して速度の計算を行う。
【００４６】
ステップＳ５で適用される数値積分法としては、差分近似法や予測子−修正子法などの公知の方法を適用することができる。このうち、差分近似法として、オイラー法、修正オイラー法、ベルレ法、速度ベルレ法、かえる跳び法などが適用可能であるが、計算精度の観点からは速度ベルレ法またはかえる跳び法が好適である。
【００４７】
ステップＳ５で適用される温度制御法としては、速度スケーリング法、Ｎｏｓｅ−Ｈｏｏｖｅｒ法、Ｎｏｓｅ−Ｐｏｉｎｃａｒｅ法、Ｂｅｒｅｎｄｓｅｎ法などの公知の方法を適用することができる。このうち、計算効率化の観点からは速度スケーリング法が好適である。
【００４８】
速度スケーリング法では、粒子の質量ｍ、現在の時刻における粒子の平均自乗速度Ｖ2、ボルツマン定数ｋBを用いて、現在の時刻における系の温度を
【数３】
2013211006.tif 000005
によって求める。また、粒子ｉの速度Ｖiを、以下の式（４）にしたがって変換する。
【数４】
2013211006.tif 000006
式（４）の右辺のＴSETは、設定温度を表す。
【００４９】
なお、実験結果と比較する際に温度一定の条件を課す必要がない場合、速度計算部５４は速度の計算に温度制御法を適用しない。
【００５０】
ステップＳ５の後、座標計算部５５は、パラメータ記憶部４１から粒子の種類、質量、位置座標、速度、時間刻みを呼び出すとともに、計算データ記憶部４２から相互作用力を呼び出し、時間刻み分だけ進めた新しい時刻の座標を数値積分によって計算する（ステップＳ６）。このステップＳ６においても、上述したステップＳ５と同様の数値積分法を適用することができる。
【００５１】
ＭＤシミュレーション装置１は、以上説明したステップＳ２〜Ｓ６の処理を規定回数だけ繰り返して実行した後、一連の処理を終了する。
【００５２】
ここで、従来のＭＤシミュレーション方法との差異を説明する。従来法において、ペア配列関数計算部５２は、系内の粒子の全ての組み合わせについて、カットオフ関数およびその一次微分の値を計算し、２次元配列変数３０１として計算データ記憶部４２へ格納する。図６は、従来のＭＤシミュレーション方法で用いる２次元配列変数の構成例を示す図である。図６に示す２次元配列変数３０１は、サイズがＮ×Ｎのマトリックスの各要素に変数が格納されたものである。この場合、計算におけるメモリの使用量は粒子数の自乗（Ｎ2）に比例しており、粒子数Ｎが大きくなると加速度的にメモリ使用量が大きくなり、ＭＤシミュレーション装置１の計算能力の限界を超えてメモリ不足となってしまう。
【００５３】
これに対し、本実施の形態１において、２次元配列変数２０１の行数Ｘを上述したＮiの最大値とすると、行数Ｘは粒子数Ｎに依存せずほぼ一定となる。このため、メモリ使用量は、粒子数Ｎのみに比例することとなる。また、本実施の形態１において、相互作用していない粒子ペアについてのカットオフ関数およびその一次微分の値はどちらも０になるため、相互作用力の計算には影響しない。したがって、本実施の形態１では、従来法でメモリ不足となるような粒子数であっても、ＭＤシミュレーション装置１の計算能力の限界を超えることがない。
【００５４】
以上説明した本発明の実施の形態１によれば、カットオフ領域よりもわずかに広い探索範囲内に含まれる粒子に対してリスト化を行うため、カットオフ関数や相互作用力の計算を行う際、系内の全ての粒子ペアについて計算を行う必要がない。したがって、本実施の形態１によれば、粒子数が大きい系に対しても適切な計算時間で計算することができ、例えば計算モデルの大きさを現実の材料の大きさに近づけることが可能となる。その結果、本実施の形態１によれば、計算によって得られた材料物性や構造に関して、実験や測定による結果との比較が可能となり、物性、構造予測の信頼性を向上させることができる。
【００５５】
（実施の形態２）
本発明の実施の形態２は、ペア配列関数として、２つの粒子の近傍に第３の粒子が存在する場合、この第３の粒子が２つの粒子間の相互作用を遮蔽する効果を表すスクリーン関数を適用することを特徴とする。
【００５６】
本実施の形態２に係る分子動力学シミュレーション装置の構成は、実施の形態１で説明したＭＤシミュレーション装置１の構成と同様である。このため、以下の説明では、ＭＤシミュレーション装置１と同じ符号を用いて説明を行う。また、本実施の形態２に係る分子動力学シミュレーション方法は、ペア配列関数としてスクリーン関数を適用する点を除いて、上述した実施の形態１に係る分子動力学シミュレーション方法と同様である。
【００５７】
ここで、スクリーン関数について説明する。スクリーン関数は、カットオフ関数と同様にポテンシャル関数全体に陽に含まれる関数であり、０から１までの値をとる関数である。スクリーン関数は、カットオフ関数とは異なり、第３の粒子との位置関係によって変化する。スクリーン関数の値について、より具体的に説明する。以下、２つの粒子を結ぶ線分を長径として、短径が異なる２つの楕円体を設定し、小さい短径を有する楕円体を第１の楕円体、大きい短径を有する楕円体を第２の楕円体と称する。第１の楕円体の内部に第３の粒子が存在する場合、スクリーン関数の値は０となり、この第３の粒子が前記２つの粒子間の相互作用を生じさせないように結合を遮蔽する。これに対し、第２の楕円体の内部に第３の粒子が存在しない場合、スクリーン関数の値は１、すなわち、遮蔽効果は０となり、前記２つの粒子間に相互作用力を生じさせる。上記以外の場合、スクリーン関数は０より大きく、かつ１より小さい値をとり、その値に応じて前記２つの粒子間の相互作用を弱める効果を表す。ここで、「楕円体の内部に第３の粒子が存在する」とは、楕円体の内部に第３の粒子の座標点が少なくとも一つ入っていることを意味する。また、「楕円体の内部に第３の粒子が存在しない」とは、楕円体の内部に第３の粒子の座標点が一つも入っていないことを意味する。このような座標点として、例えば粒子を球とみなしたときのその球の中心を取ることができる。以下、第１の楕円体の内部に第３の粒子が存在する場合、その２つの粒子の結合を、第３の粒子によって遮蔽された結合という。
【００５８】
上述した特性を有するスクリーン関数を導入した場合には、アモルファスカーボン中に含まれるｓｐ3炭素濃度を精度良く再現できることが知られている（例えば、T.  Kumagai,  S.  Hara,  J.  Choi,  S.  Izumi,  and  T.  Kato,  “Development  of  empirical  bond-order-type  interatomic  potential  for  amorphous  carbon  structures”,  Journal  of  Applied  Physics,  105,  064310-1〜10  (2009)参照）。
【００５９】
次に、リスト化部５１が行うリスト化処理における粒子の探索範囲について説明する。なお、リスト化部５１が行う処理は、実施の形態１と同じである。図７は、リスト化処理における粒子の探索範囲を模式的に示す図である。図７において、注目粒子Ｐiと別の粒子Ｐjとを結ぶ線分を長径とする第１の楕円体Ｅjの内部には、第３の粒子が存在していない。したがって、注目粒子Ｐiと粒子Ｐjとの間には相互作用力が働く。一方、注目粒子Ｐiと粒子Ｐj’とを結ぶ線分を長径とする第１の楕円体Ｅj’の内部には第３の粒子Ｐkが存在している。したがって、注目粒子Ｐiと粒子Ｐj’との間には相互作用力が働かない。このようにして、第３の粒子によって遮蔽された結合に関しては、相互作用力を考慮する必要がない。図７に示す場合、注目粒子Ｐiとの間に相互作用力が働く粒子は、斜線で示した粒子に限られる。
【００６０】
探索範囲に含まれる粒子は、以下の２つの規定によって決定される。
＜第１の規定＞
注目粒子Ｐiとの間で形成される第１の楕円体（図７の楕円体Ｅj、Ｅj’等）よりも短径が小さい第３の楕円体（図７の楕円体Ｅs、Ｅs’等）を設定したとき、全ての第３の粒子がこの第３の楕円体の外部に位置するような粒子。ここで、第３の楕円体の短径を小さく取るほど探索範囲を大きくすることができ、リスト化の頻度は少なくなる。
＜第２の規定＞
注目粒子Ｐiを中心として限界距離Ｒtを半径とする球形領域Ｄtよりも少しだけ大きい半径Ｒssを有する球形領域Ｄssに含まれる粒子。ここで、「限界距離」とは、２つの粒子の相互作用を考慮する必要がない程度に十分大きい距離のことをいう。また、半径Ｒssの値は、粒子の速度を考慮して定められる。
【００６１】
リスト化の頻度については、実施の形態１で説明したカットオフ関数の場合と同様に設定される。第１の規定に関しては、式（１）のカットオフ半径ＲSを、対象とする第１の楕円体の短径の１／２倍に置き換えるとともに、式（１）の半径ＲCを、その第１の楕円体に対応する第３の楕円体の短径の１／２倍に置き換えたものを有効期限ｔe(1)とする。この有効期限は、第３の楕円体の内部にある粒子が第１の楕円体の外部に出て行くまでの時間に相当する。また、第２の規定に関しては、式（１）のカットオフ半径ＲSを球Ｄsの半径ＲSSに置き換えるとともに、式（１）のカットオフ半径Ｒcを限界距離Ｒtに置き換えたものを有効期限ｔe(2)とする。この有効期限は、球形領域Ｄssの外部にある粒子が球殻状領域Ｄss-tを通過して球形領域Ｄtの内部に入ってくる時間に相当する。リスト化部５１は、これら２つの有効期限ｔe(1)、ｔe(2)のうち小さい方の値を実際の有効期限ｔeとして選択する。
【００６２】
ペア配列関数計算部５２は、パラメータ記憶部４１から粒子の種類、位置座標およびポテンシャルパラメータを呼び出し、このパラメータを用いて、リスト化された粒子ペアについて、ポテンシャル関数に含まれるペア配列関数として、スクリーン関数を適用する。具体的には、ペア配列関数計算部５２は、スクリーン関数の値を計算し、計算結果を２次元配列変数として計算データ記憶部４２へ格納する（図４のステップＳ３を参照）。ペア配列関数計算部５２が計算データ記憶部４２へ格納する２次元配列変数は、図５に示す２次元配列変数２０１と同様の構成を有する。したがって、本実施の形態２においても、実施の形態１と同様、ＭＤシミュレーション装置１の計算量を低減することができる。
【００６３】
以上説明した本発明の実施の形態２によれば、スクリーン関数や相互作用力の計算を行う際、系内の全ての粒子ペアについて計算を行う必要がないため、粒子数が大きい系に対しても適切な計算時間で計算することができ、例えば計算モデルの大きさを現実の材料の大きさに近づけることが可能となる。これにより、本実施の形態２によれば、計算によって得られた材料物性や構造に関して、実験や測定による結果との比較が可能となり、物性、構造予測の信頼性を向上させることができる。
【００６４】
ここまで、本発明の２つの実施の形態を説明してきたが、本発明は、それらの実施の形態によってのみ限定されるべきものではない。例えば、本発明では、ペア配列関数として、カットオフ関数およびスクリーン関数を含む関数を適用することも可能である。
【実施例】
【００６５】
以下、本発明の実施例を詳細に説明する。なお、本発明は、後述する実施例によって限定されるものではない。
【００６６】
（実施例１）
実施例１では、系の全粒子数Ｎ＝８４６７２のグラフェンの系について圧力、温度一定の条件下におけるＭＤシミュレーションを実行した。ＭＤシミュレーションの計算は、本出願人が独自に開発したＦｏｒｔｒａｎプログラムを用いて行った。計算には、メモリ１０２４ＭＢ、クロック周波数２．５０ＧＨｚのインテル（登録商標）製Ｃｅｌｅｒｏｎ（登録商標）を搭載したＨＰ（登録商標）製パソコンＣＯＭＰＡＱ（登録商標）  ６０００  Ｐｒｏを用いた。以下、上述した実施の形態１におけるＭＤシミュレーション処理のステップごとに実施例を詳細に説明する。
【００６７】
＜ステップＳ１＞
Ｆｏｒｔｒａｎプログラム中にグラフェンの原子座標を計算するコードを記載して初期座標を決定した。設定温度を３００Ｋとし、初期速度として３００ＫのＭａｘｗｅｌｌ分布に従う速度を、乱数を用いて各原子にランダムに与えた。時間刻みを０．２ｆｓとした。緩和時間を４ｆｓ（ループ換算で２０回繰り返しに相当）とした。ポテンシャル関数はＴｅｒｓｏｆｆポテンシャルとし、Ｔｅｒｓｏｆｆが提唱したオリジナルのパラメータを用いた。ポテンシャル関数の具体的な数式は、次式（５）で表される。
【数５】
2013211006.tif 000007
式（５）の右辺の第１項は斥力項であり、第２項は引力項である。式（５）の右辺のｒijは、原子ｉと原子ｊの原子間距離を表す。また、カットオフ関数ｆC（ｒ）は、
【数６】
2013211006.tif 000008
で与えられる（図３を参照）。
【００６８】
式（５）の右辺第２項の係数ｂijは、
【数７】
2013211006.tif 000009
で与えられる。ここで、ζijは、
【数８】
2013211006.tif 000010
であり、この式（８）の右辺のｇ（θijk）は、
【数９】
2013211006.tif 000011
である。変数θijkは、ｉ−ｊ結合とｉ−ｋ結合の間の結合角である。
【００６９】
本実施例１において、式（５）〜（９）に含まれるパラメータＡ、Ｂ、λ、μ、β、ｎ、ｃ、ｄ、ｈ、Ｒmin、Ｒmax（＝Ｒc）については、以下の表１に示す値を用いた。
【表１】
2013211006.tif 000012
また、ステップＳ２においてリスト化を行う際に必要な探索半径Ｒsは、２．６Åとした。
【００７０】
＜ステップＳ２＞
各原子の座標を用いて原子間距離ｒijを計算し、探索範囲内の原子ペアをリスト化した。このとき、任意の原子ｉとペアをなす原子の数Ｎiの最大値Ｘを同時に求めた。その結果、Ｘの値は全ての時間ステップを通じて１０程度であった。
【００７１】
＜ステップＳ３＞
探索範囲内の原子ペアについてのみ、カットオフ関数およびカットオフ関数の一次微分の値を計算した。計算した値を入力するための変数として、図５に示す２次元配列変数２０１を用いた。
【００７２】
＜ステップＳ４＞
探索範囲内の原子ペア、原子の座標、カットオフ関数の値、およびカットオフ関数の一次微分の値を用いて、原子間に働く相互作用力を計算した。
【００７３】
＜ステップＳ５＞
相互作用力および現在の時刻における原子の速度を用いて、速度ベルレ法および速度スケーリング法により新しい時刻における原子の速度を計算した。
【００７４】
＜ステップＳ６＞
現在の時刻における原子の座標、速度および原子に働く相互作用力の値を用いて、速度ベルレ法により新しい時刻における原子の座標を計算した。
【００７５】
上記ステップＳ２〜Ｓ６を２０回繰り返したところ、問題なく計算が完了した。計算の所要時間は９４．５秒であった。
【００７６】
（実施例２）
実施例２では、全粒子数Ｎを２１１６８として、実施例１と同様に計算を実行した。その結果、所要時間７．１秒で問題なく計算が完了した。
【００７７】
（実施例３）
実施例３では、全粒子数Ｎを１２１６０として、実施例１と同様に計算を実行した。その結果、所要時間２．７秒で問題なく計算が完了した。
【００７８】
（実施例４）
実施例４では、全粒子数Ｎを５２９２として、実施例１と同様に計算を実行した。その結果、所要時間０．７秒で問題なく計算が完了した。
【００７９】
（実施例５）
実施例５では、全粒子数Ｎを１２１６０とし、実施例１の計算におけるステップＳ２、Ｓ３の内容を以下のように変更して計算を実行した。
＜ステップＳ２＞
各原子の座標を用いて、各原子ｉとｊの間の原子間距離を計算し、探索範囲内の原子ペアをリスト化する。
＜ステップＳ３＞
探索範囲内の原子ペアおよび原子の座標を用いて、探索範囲内の原子ペアについてのみ、カットオフ関数およびカットオフ関数の一次微分の値を計算する。計算した値を入力するための変数として、図５に示す２次元配列変数２０１を用いる。ここで、Ｘ＝Ｎ／２＝６０８０とした。
【００８０】
以上の条件で計算を実行したところ、所要時間２．７秒で問題なく計算が完了した。
【００８１】
（実施例６）
実施例６では、全粒子数Ｎを１２１６０とし、実施例１の計算におけるステップＳ３でカットオフ関数およびカットオフ関数の一次微分の値を入力する変数を、図６に示す２次元配列変数３０１として、探索範囲内の原子ペアについてのみ計算を実行した。その結果、所要時間３．７秒で問題なく計算が完了した。
【００８２】
（実施例７）
実施例７では、全粒子数Ｎを５２９２とし、実施例６と同様に計算を実行した。その結果、所要時間０．８秒で問題なく計算が完了した。
【００８３】
（参考例１）
参考例１では、全粒子数Ｎを２１１６８とし、実施例５と同様の計算を実行しようとした（Ｘ＝Ｎ／２＝１０５８４）。しかしながら、メモリ不足となって計算を実行することができなかった。
【００８４】
（参考例２）
参考例２では、実施例１（Ｎ＝８４６７２）において、ステップＳ３におけるカットオフ関数およびカットオフ関数の一次微分の計算を、探索範囲内の原子ペアについてのみ実行し、図６に示す２次元配列変数３０１に値を入力してステップＳ４以降の計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【００８５】
（参考例３）
参考例３では、全粒子数Ｎを２１１６８とし、参考例２と同様にて計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【００８６】
（比較例１）
比較例１では、実施例１（Ｎ＝８４６７２）において、ステップＳ３におけるカットオフ関数およびカットオフ関数の一次微分の値の計算を、２つの原子の間の全ての組み合わせであるＮ×（Ｎ−1）通りについて実行し、図６に示す２次元配列変数３０１に値を入力してステップＳ４以降の計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【００８７】
（比較例２）
比較例２では、全粒子数Ｎを２１１６８とし、比較例１と同様の計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【００８８】
（比較例３）
比較例３では、全粒子数Ｎを１２１６０とし、比較例１と同様の計算を実行した。その結果、計算は完了したものの、所要時間は１０分以上であった。
【００８９】
（比較例４）
比較例４では、全粒子数Ｎを５２９２とし、比較例１と同様の計算を実行した。その結果、計算は完了したものの、所要時間は４５秒であった。
【００９０】
以下に示す表２は、実施例１〜７、参考例１〜３および比較例１〜４の全粒子数Ｎ、計算の所要時間（計算時間）、２次元配列変数、２次元配列変数３０１を採用した場合の計算方法、および２次元配列変数２０１を採用した場合の行数Ｘを示している。なお、表２では、２次元配列変数２０１、３０１をそれぞれＡ、Ｂと記載している。また、２次元配列変数３０１を使用する例では、ステップＳ３におけるカットオフ関数およびカットオフ関数の一次微分の値の計算方法について、探索範囲内の原子ペアについてのみ行う場合をＩ、全ての組み合わせであるＮ×（Ｎ−1）通りについて行う場合をＩＩと記載している。
【表２】
2013211006.tif 000013
【００９１】
図８は、実施例１〜７、参考例１〜３および比較例１〜４のうち計算を完了した例における全粒子数（Ｎ）と計算時間（秒）との関係を両対数グラフで示した図である。図８において、点Ｅ１〜Ｅ７は、実施例１〜７に順次対応している。また、点Ｃ３、Ｃ４は、比較例３、４にそれぞれ対応している。
【００９２】
全粒子数Ｎ＝８４６７２で計算した実施例１および比較例１、ならびに全粒子数Ｎ＝２１１６８で計算した実施例２および比較例２とを比較すると、本発明においては、全粒子数Ｎが、従来法で適用不可能な原子数の４倍以上であっても、メモリ不足とならずに実行できることが分かる。
【００９３】
また、同じ全粒子数Ｎ＝１２１６０で計算した実施例３と比較例３とを比較することにより、本発明は従来法に対して最大２００分の１以下の時間短縮となっていることも分かる。
【００９４】
（実施例８）
実施例８では、系の全粒子数Ｎ＝８４６７２のグラフェンの系について体積、温度一定の条件下におけるＭＤシミュレーションを実行した。ＭＤシミュレーションの計算は、本出願人が独自に開発したＦｏｒｔｒａｎプログラムを用いて行った。計算に使用したパソコンは、実施例１と同じである。以下、上述した実施の形態２におけるＭＤシミュレーション処理のステップごとに実施例を詳細に説明する。
＜ステップＳ１＞
Ｆｏｒｔｒａｎプログラム中にグラフェンの原子座標を計算するコードを記載して初期座標を決定した。設定温度を３００Ｋとし、初期速度として３００ＫのＭａｘｗｅｌｌ分布に従う速度を、乱数を用いて各原子にランダムに与えた。時間刻みを０．２ｆｓとした。緩和時間を４ｆｓ（ループ換算で２０回繰り返しに相当）とした。ポテンシャル関数は、Ｔｅｒｓｏｆｆポテンシャルを参考として、カットオフ関数をスクリーン関数に置き換えた関数を適用した。スクリーン関数以外のパラメータは、Ｔｅｒｓｏｆｆが提唱したオリジナルのパラメータを用いた。ポテンシャル関数の具体的な数式は、次式（１０）で表される。
【数１０】
2013211006.tif 000014
式（１０）の右辺の第１項は斥力項であり、第２項は引力項である。式（１０）の右辺のｒijは、原子ｉと原子ｊの原子間距離を表す。式（１０）の右辺第２項の係数ｂijは、式（７）で与えられる。
【００９５】
次に、式（１０）の右辺にあるスクリーン関数ｆSC（ｒij）について説明する。スクリーン関数ｆSC（ｒij）は、
【数１１】
2013211006.tif 000015
で与えられる。ここで、式（１１）の右辺のＳijkは、
【数１２】
2013211006.tif 000016
である。第１および第２の楕円体の短径はぞれぞれ、ｒij（Ｃmin）1/2、ｒij（Ｃmax）1/2で表される。また、Ｃminより小さい定数をＣsとして、第３の楕円体の短径はｒij（Ｃs）1/2で表される。
また、式（１２）の右辺のＣijkは、
【数１３】
2013211006.tif 000017
で与えられる。ここで、式（１３）の右辺のＸik、Ｘjkは、
【数１４】
2013211006.tif 000018
である。
また、式（１２）の右辺の関数ｆs（ｒ）は、
【数１５】
2013211006.tif 000019
で与えられる。
【００９６】
本実施例８において、式（１０）〜（１５）に含まれるパラメータＡ、Ｂ、λ、μ、β、ｎ、ｃ、ｄ、ｈについては、表１に示す値とし、式（１２）のパラメータＣmin、Ｃmaxについては、それぞれ１．０、２．０とした。また、ステップＳ２においてリスト化を行う際に必要な探索半径Ｒsを４．０Åとし、Ｃsを０．２とした。このとき、楕円体Ｅsの短径ｄsは、ｄs＝ｒij（Ｃｓ）1/2となる。
【００９７】
＜ステップＳ２＞
各原子の座標を用いて原子間距離ｒijを計算し、探索範囲内の原子ペアをリスト化した。このとき、任意の原子ｉとペアをなす原子の数Ｎiの最大値Ｘを同時に求めた。その結果、Ｘの値は全ての時間ステップを通じて１０程度であった。
＜ステップＳ３＞
探索範囲内の原子ペアについてのみ、スクリーン関数の値を計算した。計算した値を入力するための変数として、図５に示す２次元配列変数２０１を用いた。
＜ステップＳ４＞
探索範囲内の原子ペア、原子の座標、スクリーン関数の値を用いて、原子間に働く相互作用力を計算した。
＜ステップＳ５＞
相互作用力および現在の時刻における原子の速度を用いて、速度ベルレ法および速度スケーリング法により新しい時刻における原子の速度を計算した。
＜ステップＳ６＞
現在の時刻における原子の座標、速度および原子に働く相互作用力の値を用いて、速度ベルレ法により新しい時刻における原子の座標を計算した。
【００９８】
上記ステップＳ２〜Ｓ６を２０回繰り返したところ、問題なく計算が完了した。計算の所要時間は１８６秒であった。
【００９９】
（実施例９）
実施例９では、全粒子数Ｎを２１１６８として、実施例８と同様に計算を実行した。その結果、所要時間２８．８秒で問題なく計算が完了した。
【０１００】
（実施例１０）
実施例１０では、全粒子数Ｎを１２１６０として、実施例８と同様に計算を実行した。その結果、所要時間１５．８秒で問題なく計算が完了した。
【０１０１】
（実施例１１）
実施例１１では、全粒子数Ｎを５２９２として、実施例８と同様に計算を実行した。その結果、所要時間６．７２秒で問題なく計算が完了した。
【０１０２】
（実施例１２）
実施例１２では、全粒子数Ｎを１２１６０とし、実施例８の計算におけるステップＳ２、Ｓ３の内容を以下のように変更して計算を実行した。
＜ステップＳ２＞
各原子の座標を用いて、各原子ｉとｊの間の原子間距離を計算し、探索範囲内の原子ペアをリスト化する。
＜ステップＳ３＞
探索範囲内の原子ペアおよび原子の座標を用いて、探索範囲内の原子ペアについてのみ、スクリーン関数の値を計算する。計算した値を入力するための変数として、２次元配列変数２０１を用いる。ここで、Ｘ＝Ｎ／２＝６０８０とした。
以上の条件で計算を実行したところ、所要時間１９．８秒で問題なく計算が完了した。
【０１０３】
（実施例１３）
実施例１３では、全粒子数Ｎを１２１６０とし、実施例８の計算におけるステップＳ３でスクリーン関数の値を入力する変数を、２次元配列変数３０１として、探索範囲内の原子ペアについてのみ計算を実行した。その結果、所要時間５７．５秒で問題なく計算が完了した。
【０１０４】
（実施例１４）
実施例１４では、全粒子数Ｎを５２９２とし、実施例１３と同様に計算を実行した。その結果、所要時間７．８４秒で問題なく計算が完了した。
【０１０５】
（参考例４）
参考例４では、全粒子数Ｎを２１１６８とし、実施例１２と同様の計算を実行しようとした（Ｘ＝Ｎ／２＝１０５８４）。しかしながら、メモリ不足となって計算を実行することができなかった。
【０１０６】
（参考例５）
参考例５では、実施例８（Ｎ＝８４６７２）において、ステップＳ３におけるスクリーン関数の値の計算を、探索範囲内の原子ペアについてのみ実行し、２次元配列変数３０１に値を入力してステップＳ４以降の計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【０１０７】
（参考例６）
参考例６では、全粒子数Ｎを２１１６８とし、参考例５と同様にて計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【０１０８】
（比較例５）
比較例５では、実施例８（Ｎ＝８４６７２）において、ステップＳ３におけるスクリーン関数の値の計算を、２つの原子の間の全ての組み合わせであるＮ×（Ｎ−1）通りについて実行し、図６に示す２次元配列変数３０１に値を入力してステップＳ４以降の計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【０１０９】
（比較例６）
比較例６では、全粒子数Ｎを２１１６８とし、比較例５と同様の計算を実行しようとした。しかしながら、メモリ不足となって計算を実行することができなかった。
【０１１０】
（比較例７）
比較例７では、全粒子数Ｎを１２１６０とし、比較例５と同様の計算を実行した。その結果、計算は完了したものの、所要時間は１０分以上であった。
【０１１１】
（比較例８）
比較例８では、全粒子数Ｎを５２９２とし、比較例５と同様の計算を実行した。その結果、計算は完了したものの、所要時間は１０分以上であった。
【０１１２】
以下に示す表３は、実施例８〜１４、参考例４〜６および比較例５〜８の全粒子数Ｎ、計算の所要時間（計算時間）、２次元配列変数、２次元配列変数３０１を採用した場合の計算方法、および２次元配列変数２０１を採用した場合の行数Ｘを示している。なお、表３では、２次元配列変数２０１、３０１をそれぞれＡ、Ｂと記載している。また、２次元配列変数３０１を使用する例では、ステップＳ３におけるスクリーン関数の値の計算方法について、探索範囲内の原子ペアについてのみ行う場合をＩ、全ての組み合わせであるＮ×（Ｎ−1）通りについて行う場合をＩＩと記載している。
【表３】
2013211006.tif 000020
【０１１３】
図９は、実施例８〜１４、参考例４〜６および比較例５〜８のうち計算を完了した例における全粒子数（Ｎ）と計算時間（秒）との関係を両対数グラフで示した図である。図９において、点Ｅ８〜Ｅ１４は、実施例８〜１４に順次対応している。また、点Ｃ７、Ｃ８は、比較例７、８にそれぞれ対応している。
【０１１４】
このように、本発明によれば、従来法よりも適用可能な粒子数を大幅に増加させることができ、なおかつ計算時間を短縮することができる。
【符号の説明】
【０１１５】
１  ＭＤシミュレーション装置
２  入力部
３  表示部
４  記憶部
５  制御部
４１  パラメータ記憶部
４２  計算データ記憶部
４３  プログラム記憶部
５１  リスト化部
５２  ペア配列関数計算部
５３  相互作用力計算部
５４  速度計算部
５５  座標計算部
１０１  カットオフ関数
２０１、３０１  ２次元配列変数
【図１】
2013211006.tif 000021
【図２】
2013211006.tif 000022
【図３】
2013211006.tif 000023
【図４】
2013211006.tif 000024
【図５】
2013211006.tif 000025
【図６】
2013211006.tif 000026
【図７】
2013211006.tif 000027
【図８】
2013211006.tif 000028
【図９】
2013211006.tif 000029