JP2010015284A - シミュレーション方法及びプログラム - Google Patents

Abstract

【課題】 高精度なシミュレーションを実現する。
【解決手段】 （ａ）シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦを、ポテンシャルが最小となる距離に対応するパラメータｒ_０、ポテンシャルの形状を規定する少なくとも１つのパラメータＣ_０を用い、２つの粒子の中心間距離ｒの関数Φ（ｒ，ｒ_０，Ｃ_０）として決定する。（ｂ）パラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えたとき、シミュレーション対象の材料のヤング率の変化が抑制されるように、パラメータＣ_０をＣ_１に置き換える。（ｃ）新たなパラメータｒ_１とＣ_１とで定義される二体間ポテンシャルΦ（ｒ，ｒ_１，Ｃ_１）を用いてシミュレーションを行う。
【選択図】 図８

Description

本発明は、シミュレーション方法及びプログラムに関し、特に、分子動力学を用いたシミュレーション方法、及びそれをコンピュータに実行させるためのプログラムに関する。

コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。

分子動力学法に基づく分子シミュレーションに関して、種々の研究がなされている。たとえば少ない計算量で分子動力学計算を実行するシミュレーション方法の発明が開示されている（たとえば、特許文献１参照）。

しかしながらこれまでのところ、分子動力学シミュレーションを行う際のポテンシャルパラメータとヤング率との関係は明確にされておらず、たとえばヤング率を一定としたまま、ポテンシャルパラメータを変更してシミュレーションを行うことは困難であった。

特開２００６−２８５８６６号公報

本発明の目的は、高精度なシミュレーションを実現するシミュレーション方法を提供することである。

また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することである。

本発明の一観点によれば、（ａ）シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦを、ポテンシャルが最小となる距離に対応するパラメータｒ_０、ポテンシャルの形状を規定する少なくとも１つのパラメータＣ_０を用い、２つの粒子の中心間距離ｒの関数Φ（ｒ，ｒ_０，Ｃ_０）として決定する工程と、（ｂ）パラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えたとき、シミュレーション対象の材料のヤング率の変化が抑制されるように、パラメータＣ_０をＣ_１に置き換える工程と、（ｃ）新たなパラメータｒ_１とＣ_１とで定義される二体間ポテンシャルΦ（ｒ，ｒ_１，Ｃ_１）を用いてシミュレーションを行う工程とを有するシミュレーション方法が提供される。

また、本発明の他の観点によると、種々のシミュレーションを行うためにコンピュータを、シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦ（ｒ，ｒ_０，Ｃ_０）が最小となる距離に対応するパラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えるための情報を入力する手段、前記ｒ_０をｒ_１に置き換えたとき、前記シミュレーション対象の材料のヤング率の変化が抑制されるように、前記二体間ポテンシャルの形状を規定するパラメータＣ_０を置き換えるべきＣ_１を得る手段、新たなパラメータｒ_１とＣ_１とで定義される二体間ポテンシャルΦ（ｒ，ｒ_１，Ｃ_１）を用いてシミュレーションを行う手段として機能させるためのプログラムが提供される。

本発明によれば、高精度なシミュレーションを実現するシミュレーション方法を提供することができる。

また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することができる。

二体間ポテンシャルは、たとえば下式（１）で表される。

・・（１）
二体間ポテンシャルとは、ポテンシャルφが２個の原子の中心間の距離ｒのみを独立変数として定義される原子間ポテンシャルをいい、式（１）において、Ｄ、αは、ポテンシャルの形状を規定するパラメータであり、具体的には、それぞれ結合エネルギ、原子の大きさを示すポテンシャルパラメータである。また、ｒ_０は２個の原子が最安定となる原子間距離を示す。右辺のｆ（ｒ）という形は左辺のφがｒの関数であることを示す。

図１は、二体間ポテンシャルを示す概略的なグラフである。グラフの横軸は、２個の原子間の距離ｒを表し、グラフの縦軸はポテンシャルφを表す。最小（極小）となるポテンシャルφを与える原子間距離ｒがｒ_０である。

式（１）で表されるポテンシャルφを、ポテンシャルの底（ｒ＝ｒ_０）近傍でテイラー展開することで下の式（２）を得る。

・・（２）

ここでｆ”は、関数ｆを変数ｒで２回微分したものである。

式（２）の両辺をｒで微分することで、下の式（３）が得られる。式（３）の左辺のＦは、基点にある原子が距離ｒ_０近傍に位置する原子に及ぼす力を示す。

・・（３）

また、歪ε、応力σはそれぞれ下式（４）、（５）で表され、両者の間にはヤング率Ｅを介して、下式（６）の関係がある。なお、式（５）において、ｓは原子間のクロスセクションを示す。原子間のクロスセクションは、相隣る原子間において、力の働く断面積を意味し、たとえば図２に示す面積を表す。

・・（４）

・・（５）

・・（６）

式（４）、（５）を式（６）に代入することにより、下式（７）が得られる。

・・（７）

そして、式（３）と（７）とから、下式（８）を得る。

・・（８）

ここで、クロスセクションｓは、ｒ_０ ^２のオーダーであるため、

・・（９）

と書ける。

式（９）を式（８）に代入することにより、下式（１０）を得る。

・・（１０）

式（１０）は、ヤング率Ｅとポテンシャルパラメータとの関係を示している。式（１０）によると、ヤング率ＥはポテンシャルパラメータＤ、１／α^２、１／ｒ_０に依存（比例）する。

ところで式（１）は、Ｍｏｒｓｅポテンシャルの場合、

・・（１１）

と表すことができる。ここで、Ａ＝１／αの関係がある。この場合、式（１０）は、

・・（１２）

と表される。

なお、二体間ポテンシャルは、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルの場合、

・・（１３）

と表すことができる。この場合、式（１０）に相当する式は、

・・（１４）

と表される。

本願発明者らは、ヤング率Ｅとポテンシャルパラメータとの関係（（式１０）で示される関係）を確認するために、銅（Ｃｕ）のＭｏｒｓｅポテンシャルの場合（（式１２））を例にとり、以下の第１〜第３のシミュレーションを行った。

図３及び図４を参照して第１のシミュレーションについて説明する。第１のシミュレーションは、ポテンシャルパラメータＤと、ヤング率Ｅとの関係を確認するシミュレーションである。

銅のＭｏｒｓｅポテンシャルのポテンシャルパラメータは、ｒ_０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）である。第１のシミュレーションにおいては、ポテンシャルパラメータｒ_０及びＡを上記一定値とし、ポテンシャルパラメータＤを変化させて、ヤング率Ｅの変化を調べた。

図３に、第１のシミュレーションモデルを示す。第１のシミュレーションモデルは、規則正しくｆｃｃ構造をなすように配列した銅原子が、全体として直方体状の結晶を形成している。［１００］、［０１０］、［００１］方向をそれぞれｘ、ｙ、ｚ方向とし、結晶のサイズは、ｘ、ｙ、ｚ方向に対して、１．４４８（ｎｍ）、１．４４８（ｎｍ）、１．８１０（ｎｍ）である。［１００］、［０１０］、［００１］はミラー指数によって表された方向指数である。更に、ｘ方向、ｙ方向の境界条件は、周期境界条件とした。

ヤング率Ｅは、引っ張り試験シミュレーションにより、図３に示した直方体状の結晶を構成する銅原子の最上層をｚ方向に引っ張ったときの、直方体状の結晶全体のｚ方向への伸びから計算した。

図４は、第１のシミュレーションの結果を示すグラフである。グラフの横軸は、ポテンシャルパラメータＤを単位「ｅＶ」で表し、縦軸は、ヤング率Ｅを単位「ＧＰａ」で表す。

本図に示す結果から、ヤング率ＥがポテンシャルパラメータＤに比例することが確認された。

図５及び図６を参照して第２のシミュレーションについて説明する。第２のシミュレーションは、ポテンシャルパラメータｒ_０及びＡと、ヤング率Ｅとの関係を確認するシミュレーションである。

第２のシミュレーションにおいては、ポテンシャルパラメータＤを一定とし、ポテンシャルパラメータｒ_０及びＡを、ｒ_０とＡとの積を一定としつつ変化させて、ヤング率Ｅの変化を調べた。すなわち、ｒ_０をｘｒ_０、ＡをＡ／ｘとし、パラメータｘを変化させることでｘｒ_０及びＡ／ｘを変化させて、ヤング率Ｅとｘとの関係を調べた。

なお、銅のＭｏｒｓｅポテンシャルのポテンシャルパラメータ、ｒ_０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）の値はそのまま用いた。

図５は、ポテンシャルパラメータを変化させる前後のポテンシャルφを示すグラフである。グラフの横軸は、基点原子（ｒ＝０の位置に存在する原子）からの距離ｒを単位「ｍ」で表し、グラフの縦軸はポテンシャルφを単位「ｅＶ」で表す。

曲線ｐ（ｒ）はパラメータを変化させない場合（ｘ＝１の場合）の基点原子からの距離ｒとポテンシャルとの関係を示す。曲線ｑ（ｒ）はパラメータを変化させた場合の基点原子からの距離ｒとポテンシャルとの関係を示す。曲線ｑ（ｒ）は、ｘ＝２としたとき、すなわちｒ_０を２ｒ_０、ＡをＡ／２としたときの関係を示すグラフである。

Ｐ_１、Ｐ_２、Ｐ_３はそれぞれポテンシャルパラメータを変化させない場合（ｘ＝１の場合）における、基点原子に最も近い原子、２番めに近い原子、３番めに近い原子の原子位置（基点原子からの距離）を表す。

また、Ｑ_１、Ｑ_２、Ｑ_３はそれぞれｘ＝２としてポテンシャルパラメータを変化させた場合における、基点原子に最も近い原子、２番めに近い原子、３番めに近い原子の原子位置（基点原子からの距離）を表す。

パラメータｒ_０をｘｒ_０、ＡをＡ／ｘと変化させることは、１次元的に考えたとき、ｘの値にかかわらず、基点原子にｔ番め（ｔ＝１、２、３）に近い原子のポテンシャル値を、変化前後でそれぞれ常に同一に保つ、という意味を有する。たとえば、図５においては、ｐ（Ｐ_１）＝ｑ（Ｑ_１）、ｐ（Ｐ_２）＝ｑ（Ｑ_２）、ｐ（Ｐ_３）＝ｑ（Ｑ_３）の関係が満たされるように、ポテンシャルパラメータｒ_０及びＡは変化されることになる。また、ｘの値にかかわらず、変化前のｒにおけるポテンシャル値と、変化後のｘｒにおけるポテンシャル値とを同一とすることを意味する。

第２のシミュレーションにおいても、第１のシミュレーションと同様の引っ張り試験シミュレーションモデルを用いてシミュレーションを行った。

図６は、第２のシミュレーションの結果を示すグラフである。グラフの横軸は１／ｘ^３を表し、縦軸はヤング率Ｅを単位「ＧＰａ」で表す。

本図に示す結果から、ヤング率Ｅが１／ｘ^３に比例することが確認される。このことはヤング率Ｅが、Ａ^２／ｒ_０に比例することをも意味する。

最後に、第３のシミュレーションについて説明する。第３のシミュレーションにおいては、ｒ_０をｘｒ_０、ＡをＡ／ｘ、Ｄをｘ^３Ｄとして、パラメータｘを変化させ、ヤング率Ｅとｘとの関係を調べた。

第３のシミュレーションにおいても、第１及び第２のシミュレーションと同様の引っ張り試験シミュレーションモデルを用いた。また、第２のシミュレーション同様、銅のＭｏｒｓｅポテンシャルのポテンシャルパラメータ、ｒ_０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）の値はそのまま用いた。

図７は、第３のシミュレーションの結果を示すグラフである。グラフの横軸はｘを表し、縦軸はヤング率Ｅを単位「ＧＰａ」で表す。グラフには、Ｗ、Ｃｒ、Ｆｅ、Ｐｔ、Ｃｕ、Ａｌ、及びＰｂの各金属のヤング率を書き入れてある。

本図に示す結果から、ヤング率Ｅはｘの値に関らず、ほぼ一定に保たれることが確認される。第３のシミュレーションで計算されたヤング率Ｅの値は、実際の銅（Ｃｕ）のヤング率に接近し、ＰｔやＡｌ、その他の金属のヤング率と混同を生じさせるものではないこともわかる。また、第３のシミュレーションの結果は、ヤング率Ｅが、ＤＡ^２／ｒ_０に比例（比例定数は約１）することをも示している。

第３のシミュレーションの結果から、ヤング率を一定としたまま、ｒ_０、Ａ、Ｄのポテンシャルパラメータを変更して種々のシミュレーションを行うことができることが理解される。たとえばヤング率を一定に保ったままｒ_０を大きくし、種々のシミュレーションを行うことができる。

図８を参照して、実施例によるシミュレーション方法を説明する。

ステップＳ１０１において、シミュレーション対象の材料を構成する粒子の二体間ポテンシャルφを、ポテンシャルが最小となる距離に対応するパラメータｒ_０、ポテンシャルの形状を規定するパラメータＣ_０を用い、２つの粒子の中心間距離ｒの関数φ（ｒ，ｒ_０，Ｃ_０）として決定する。

たとえばシミュレーション対象の材料が銅である場合、二体間ポテンシャルφは、

・・（１１）

で表され、ｒ_０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）である。

次に、ステップＳ１０２において、パラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えたとき、シミュレーション対象の材料、ここでは銅のヤング率の変化が抑制されるように、パラメータＣ_０をＣ_１に置き換える。

たとえば最安定原子間距離ｒ_０を１０倍として、シミュレーションを行うのであれば、ｒ_０を１０＊ｒ_０に置き換える。この置き換えに伴って、ＡをＡ／１０、Ｄを１０００＊Ｄに置き換える。

そして、ステップＳ１０３において、新たなパラメータｒ_１（実施例においては、１０＊ｒ_０）とＣ_１（実施例においては、Ａ／１０、及び、１０００＊Ｄ）とで定義される二体間ポテンシャルφ（ｒ，ｒ_１，Ｃ_１）を用いてシミュレーションを行う。

実施例によるシミュレーション方法によれば、高精度のシミュレーションを実現することができる。任意のヤング率に対し、ヤング率を変化させずにシミュレーションを行うことも可能である。

また、最安定原子間距離を任意に設定することができるため、多様なシミュレーションを実現することができる。

シミュレーションに要する種々の計算は、プログラムにより、コンピュータで実行することができる。

図９（Ａ）は、実施例によるシミュレーション方法を実行するシミュレーション実行装置のシステム構成図であり、図９（Ｂ）は、当該装置により実行される処理のフローである。

まず、ステップＳ２０１において、キーボードなどの入力装置から、シミュレーション対象の材料を構成する粒子の二体間ポテンシャルが最小となる距離に対応するパラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えるための情報を入力する。

たとえば、シミュレーション対象の材料が銅であること、及び、銅の最安定原子間距離ｒ_０を１０倍にしてシミュレーションを行うことを入力する。

ここでファイル装置のデータファイルには、たとえば銅の二体間ポテンシャルを示すＭｏｒｓｅポテンシャルにおいては、ｒ_０＝０．２８６６（ｎｍ）、Ａ＝１３．５８８（１／ｎｍ）、Ｄ＝０．３４２９（ｅＶ）であることが記録されている。

続いて、ステップＳ２０２において、ｒ_０を１０＊ｒ_０に置き換えたとき、シミュレーション対象の材料（銅）のヤング率の変化が抑制されるように、Ｍｏｒｓｅポテンシャルの形状を規定するパラメータＡ及びＤを置き換える。

たとえば中央処理装置は、メインメモリ中の制御プログラムの指令を受け、データファイルから上記数値を読み出し、シミュレーションを行うに当たっては、ｒ_０＝２．８６６（ｎｍ）、Ａ＝１．３５８８（１／ｎｍ）、Ｄ＝３４２．９（ｅＶ）とすることを算出する。

そして、ステップＳ２０３において、中央処理装置は、メインメモリ中の制御プログラムの指令を受け、銅で形成された対象につき、ｒ_０＝２．８６６（ｎｍ）、Ｄ＝３４２．９（ｅＶ）、Ａ＝１．３５８８（１／ｎｍ）として、メインメモリに保存されているシミュレーション実行プログラムを実行し、シミュレーションを行う。シミュレーションの結果は、ディスプレイなどの出力装置を介して表示される。

上述の説明においては、ステップＳ２０１において、シミュレーション対象の材料が銅であること、及び、銅の最安定原子間距離ｒ_０を１０倍にしてシミュレーションを行うことを入力したが、たとえば直接、ｒ_０（ｒ_１）を２．８６６（ｎｍ）としてシミュレーションを行うことが入力されてもよい。

またステップＳ２０２において、Ａ＝１．３５８８（１／ｎｍ）、Ｄ＝３４２．９（ｅＶ）とすることを算出したが、たとえばファイル装置に変換テーブルを記憶させておき、演算ではなく、ｒ_０（ｒ_１）の値をもとに、参照してこれらの値を得ることも可能である。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。

たとえばＭｏｒｓｅポテンシャルの場合について、ヤング率Ｅとポテンシャルパラメータとの関係を検証したが、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルその他の場合についても、同様に確認することが可能であろう。

このため、Ｍｏｒｓｅポテンシャル以外のポテンシャルで表される材料についても、たとえばヤング率を一定としたまま、ポテンシャルパラメータを変更して種々のシミュレーションを行うことができる。

その他、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

様々なシミュレーション、殊に分子動力学シミュレーションに好適に利用することができる。

二体間ポテンシャルを示す概略的なグラフである。 原子間のクロスセクションを示す図である。 第１のシミュレーションモデルを示す。 第１のシミュレーションの結果を示すグラフである。 ポテンシャルパラメータを変化させる前後のポテンシャルφを示すグラフである。 第２のシミュレーションの結果を示すグラフである。 第３のシミュレーションの結果を示すグラフである。 実施例によるシミュレーション方法を説明するためのフローチャートである。 （Ａ）は、シミュレーション実行装置のシステム構成図であり、（Ｂ）は、当該装置により実行される処理のフローである。

符号の説明

Ｓ１０１、Ｓ１０２、Ｓ１０３、Ｓ２０１、Ｓ２０２、Ｓ２０３ ステップ

Claims (8)

1. （ａ）シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦを、ポテンシャルが最小となる距離に対応するパラメータｒ_０、ポテンシャルの形状を規定する少なくとも１つのパラメータＣ_０を用い、２つの粒子の中心間距離ｒの関数Φ（ｒ，ｒ_０，Ｃ_０）として決定する工程と、
   （ｂ）パラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えたとき、シミュレーション対象の材料のヤング率の変化が抑制されるように、パラメータＣ_０をＣ_１に置き換える工程と、
   （ｃ）新たなパラメータｒ_１とＣ_１とで定義される二体間ポテンシャルΦ（ｒ，ｒ_１，Ｃ_１）を用いてシミュレーションを行う工程と
   を有するシミュレーション方法。
2. 前記工程（ｂ）において、二体間ポテンシャルΦを与える基点となる粒子から、１番め、２番め、及び３番めに近い粒子のポテンシャル値が、それぞれ置き換え前後で同一となるように、パラメータＣ_０をＣ_１に置き換える請求項１に記載のシミュレーション方法。
3. 前記工程（ａ）において、二体間ポテンシャルを、
   
   
   
   
   
   として決定する請求項１または２に記載のシミュレーション方法。
4. 前記工程（ｂ）において、ｘを定めて、ｒ_０をｘｒ_０、ＡをＡ／ｘ、Ｄをｘ^３Ｄに置き換える請求項３に記載のシミュレーション方法。
5. 種々のシミュレーションを行うためにコンピュータを、
   シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦ（ｒ，ｒ_０，Ｃ_０）が最小となる距離に対応するパラメータｒ_０を、ｒ_０とは異なるｒ_１に置き換えるための情報を入力する手段、
   前記ｒ_０をｒ_１に置き換えたとき、前記シミュレーション対象の材料のヤング率の変化が抑制されるように、前記二体間ポテンシャルの形状を規定するパラメータＣ_０を置き換えるべきＣ_１を得る手段、
   新たなパラメータｒ_１とＣ_１とで定義される二体間ポテンシャルΦ（ｒ，ｒ_１，Ｃ_１）を用いてシミュレーションを行う手段
   として機能させるためのプログラム。
6. 前記Ｃ_１を得る手段は、
   前記二体間ポテンシャルΦを与える基点となる粒子から、１番め、２番め、及び３番めに近い粒子のポテンシャル値が、それぞれ置き換え前後で同一となるパラメータＣ_１を得る請求項５に記載のプログラム。
7. 前記二体間ポテンシャルが、
   
   
   
   
   
   で表現される請求項５または６に記載のプログラム。
8. 前記Ｃ_１を得る手段は、ｒ_０がｘｒ_０となるように置き換えられたとき、ＡをＡ／ｘ、Ｄをｘ^３Ｄに置き換えて、新たなパラメータＡ_１とＤ_１とを得る請求項７に記載のプログラム。