(57)【要約】【課題】  高精度なシミュレーションを実現する。【解決手段】  （ａ）シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦを、ポテンシャルが最小となる距離に対応するパラメータｒ０、ポテンシャルの形状を規定する少なくとも１つのパラメータＣ０を用い、２つの粒子の中心間距離ｒの関数Φ（ｒ，ｒ０，Ｃ０）として決定する。（ｂ）パラメータｒ０を、ｒ０とは異なるｒ１に置き換えたとき、シミュレーション対象の材料のヤング率の変化が抑制されるように、パラメータＣ０をＣ１に置き換える。（ｃ）新たなパラメータｒ１とＣ１とで定義される二体間ポテンシャルΦ（ｒ，ｒ１，Ｃ１）を用いてシミュレーションを行う。【選択図】  図８
2010015284.tif 000002
【特許請求の範囲】
【請求項１】
（ａ）シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦを、ポテンシャルが最小となる距離に対応するパラメータｒ０、ポテンシャルの形状を規定する少なくとも１つのパラメータＣ０を用い、２つの粒子の中心間距離ｒの関数Φ（ｒ，ｒ０，Ｃ０）として決定する工程と、
（ｂ）パラメータｒ０を、ｒ０とは異なるｒ１に置き換えたとき、シミュレーション対象の材料のヤング率の変化が抑制されるように、パラメータＣ０をＣ１に置き換える工程と、
（ｃ）新たなパラメータｒ１とＣ１とで定義される二体間ポテンシャルΦ（ｒ，ｒ１，Ｃ１）を用いてシミュレーションを行う工程と
を有するシミュレーション方法。
【請求項２】
前記工程（ｂ）において、二体間ポテンシャルΦを与える基点となる粒子から、１番め、２番め、及び３番めに近い粒子のポテンシャル値が、それぞれ置き換え前後で同一となるように、パラメータＣ０をＣ１に置き換える請求項１に記載のシミュレーション方法。
【請求項３】
前記工程（ａ）において、二体間ポテンシャルを、
2010015284.tif 000018
として決定する請求項１または２に記載のシミュレーション方法。
【請求項４】
前記工程（ｂ）において、ｘを定めて、ｒ０をｘｒ０、ＡをＡ／ｘ、Ｄをｘ３Ｄに置き換える請求項３に記載のシミュレーション方法。
【請求項５】
種々のシミュレーションを行うためにコンピュータを、
シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦ（ｒ，ｒ０，Ｃ０）が最小となる距離に対応するパラメータｒ０を、ｒ０とは異なるｒ１に置き換えるための情報を入力する手段、
前記ｒ０をｒ１に置き換えたとき、前記シミュレーション対象の材料のヤング率の変化が抑制されるように、前記二体間ポテンシャルの形状を規定するパラメータＣ０を置き換えるべきＣ１を得る手段、
新たなパラメータｒ１とＣ１とで定義される二体間ポテンシャルΦ（ｒ，ｒ１，Ｃ１）を用いてシミュレーションを行う手段
として機能させるためのプログラム。
【請求項６】
前記Ｃ１を得る手段は、
前記二体間ポテンシャルΦを与える基点となる粒子から、１番め、２番め、及び３番めに近い粒子のポテンシャル値が、それぞれ置き換え前後で同一となるパラメータＣ１を得る請求項５に記載のプログラム。
【請求項７】
前記二体間ポテンシャルが、
2010015284.tif 000019
で表現される請求項５または６に記載のプログラム。
【請求項８】
前記Ｃ１を得る手段は、ｒ０がｘｒ０となるように置き換えられたとき、ＡをＡ／ｘ、Ｄをｘ３Ｄに置き換えて、新たなパラメータＡ１とＤ１とを得る請求項７に記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びプログラムに関し、特に、分子動力学を用いたシミュレーション方法、及びそれをコンピュータに実行させるためのプログラムに関する。
【背景技術】
【０００２】
コンピュータを用いて物質科学全般の現象を探求する方法として、分子動力学法に基づく分子シミュレーションが知られている。分子シミュレーションによって、分子のポテンシャルエネルギや最安定構造等物質の特性を、分子レベルで解明することが可能である。
【０００３】
分子動力学法に基づく分子シミュレーションに関して、種々の研究がなされている。たとえば少ない計算量で分子動力学計算を実行するシミュレーション方法の発明が開示されている（たとえば、特許文献１参照）。
【０００４】
しかしながらこれまでのところ、分子動力学シミュレーションを行う際のポテンシャルパラメータとヤング率との関係は明確にされておらず、たとえばヤング率を一定としたまま、ポテンシャルパラメータを変更してシミュレーションを行うことは困難であった。
【０００５】
【特許文献１】特開２００６−２８５８６６号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本発明の目的は、高精度なシミュレーションを実現するシミュレーション方法を提供することである。
【０００７】
また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することである。
【課題を解決するための手段】
【０００８】
本発明の一観点によれば、（ａ）シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦを、ポテンシャルが最小となる距離に対応するパラメータｒ０、ポテンシャルの形状を規定する少なくとも１つのパラメータＣ０を用い、２つの粒子の中心間距離ｒの関数Φ（ｒ，ｒ０，Ｃ０）として決定する工程と、（ｂ）パラメータｒ０を、ｒ０とは異なるｒ１に置き換えたとき、シミュレーション対象の材料のヤング率の変化が抑制されるように、パラメータＣ０をＣ１に置き換える工程と、（ｃ）新たなパラメータｒ１とＣ１とで定義される二体間ポテンシャルΦ（ｒ，ｒ１，Ｃ１）を用いてシミュレーションを行う工程とを有するシミュレーション方法が提供される。
【０００９】
また、本発明の他の観点によると、種々のシミュレーションを行うためにコンピュータを、シミュレーション対象の材料を構成する粒子の二体間ポテンシャルΦ（ｒ，ｒ０，Ｃ０）が最小となる距離に対応するパラメータｒ０を、ｒ０とは異なるｒ１に置き換えるための情報を入力する手段、前記ｒ０をｒ１に置き換えたとき、前記シミュレーション対象の材料のヤング率の変化が抑制されるように、前記二体間ポテンシャルの形状を規定するパラメータＣ０を置き換えるべきＣ１を得る手段、新たなパラメータｒ１とＣ１とで定義される二体間ポテンシャルΦ（ｒ，ｒ１，Ｃ１）を用いてシミュレーションを行う手段として機能させるためのプログラムが提供される。
【発明の効果】
【００１０】
本発明によれば、高精度なシミュレーションを実現するシミュレーション方法を提供することができる。
【００１１】
また、高精度なシミュレーションを実現するシミュレーション方法を実行させるプログラムを提供することができる。
【発明を実施するための最良の形態】
【００１２】
二体間ポテンシャルは、たとえば下式（１）で表される。
2010015284.tif 000003
【００１３】
・・（１）
二体間ポテンシャルとは、ポテンシャルφが２個の原子の中心間の距離ｒのみを独立変数として定義される原子間ポテンシャルをいい、式（１）において、Ｄ、αは、ポテンシャルの形状を規定するパラメータであり、具体的には、それぞれ結合エネルギ、原子の大きさを示すポテンシャルパラメータである。また、ｒ０は２個の原子が最安定となる原子間距離を示す。右辺のｆ（ｒ）という形は左辺のφがｒの関数であることを示す。
【００１４】
図１は、二体間ポテンシャルを示す概略的なグラフである。グラフの横軸は、２個の原子間の距離ｒを表し、グラフの縦軸はポテンシャルφを表す。最小（極小）となるポテンシャルφを与える原子間距離ｒがｒ０である。
【００１５】
式（１）で表されるポテンシャルφを、ポテンシャルの底（ｒ＝ｒ０）近傍でテイラー展開することで下の式（２）を得る。
2010015284.tif 000004
【００１６】
・・（２）
ここでｆ”は、関数ｆを変数ｒで２回微分したものである。
【００１７】
式（２）の両辺をｒで微分することで、下の式（３）が得られる。式（３）の左辺のＦは、基点にある原子が距離ｒ０近傍に位置する原子に及ぼす力を示す。
2010015284.tif 000005
【００１８】
・・（３）
また、歪ε、応力σはそれぞれ下式（４）、（５）で表され、両者の間にはヤング率Ｅを介して、下式（６）の関係がある。なお、式（５）において、ｓは原子間のクロスセクションを示す。原子間のクロスセクションは、相隣る原子間において、力の働く断面積を意味し、たとえば図２に示す面積を表す。
2010015284.tif 000006
【００１９】
・・（４）
2010015284.tif 000007
【００２０】
・・（５）
2010015284.tif 000008
【００２１】
・・（６）
式（４）、（５）を式（６）に代入することにより、下式（７）が得られる。
2010015284.tif 000009
【００２２】
・・（７）
そして、式（３）と（７）とから、下式（８）を得る。
2010015284.tif 000010
【００２３】
・・（８）
ここで、クロスセクションｓは、ｒ０２のオーダーであるため、
2010015284.tif 000011
【００２４】
・・（９）
と書ける。
【００２５】
式（９）を式（８）に代入することにより、下式（１０）を得る。
2010015284.tif 000012
【００２６】
・・（１０）
式（１０）は、ヤング率Ｅとポテンシャルパラメータとの関係を示している。式（１０）によると、ヤング率ＥはポテンシャルパラメータＤ、１／α２、１／ｒ０に依存（比例）する。
【００２７】
ところで式（１）は、Ｍｏｒｓｅポテンシャルの場合、
2010015284.tif 000013
【００２８】
・・（１１）
と表すことができる。ここで、Ａ＝１／αの関係がある。この場合、式（１０）は、
2010015284.tif 000014
【００２９】
・・（１２）
と表される。
【００３０】
なお、二体間ポテンシャルは、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルの場合、
2010015284.tif 000015
【００３１】
・・（１３）
と表すことができる。この場合、式（１０）に相当する式は、
2010015284.tif 000016
【００３２】
・・（１４）
と表される。
【００３３】
本願発明者らは、ヤング率Ｅとポテンシャルパラメータとの関係（（式１０）で示される関係）を確認するために、銅（Ｃｕ）のＭｏｒｓｅポテンシャルの場合（（式１２））を例にとり、以下の第１〜第３のシミュレーションを行った。
【００３４】
図３及び図４を参照して第１のシミュレーションについて説明する。第１のシミュレーションは、ポテンシャルパラメータＤと、ヤング率Ｅとの関係を確認するシミュレーションである。
【００３５】
銅のＭｏｒｓｅポテンシャルのポテンシャルパラメータは、ｒ０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）である。第１のシミュレーションにおいては、ポテンシャルパラメータｒ０及びＡを上記一定値とし、ポテンシャルパラメータＤを変化させて、ヤング率Ｅの変化を調べた。
【００３６】
図３に、第１のシミュレーションモデルを示す。第１のシミュレーションモデルは、規則正しくｆｃｃ構造をなすように配列した銅原子が、全体として直方体状の結晶を形成している。［１００］、［０１０］、［００１］方向をそれぞれｘ、ｙ、ｚ方向とし、結晶のサイズは、ｘ、ｙ、ｚ方向に対して、１．４４８（ｎｍ）、１．４４８（ｎｍ）、１．８１０（ｎｍ）である。［１００］、［０１０］、［００１］はミラー指数によって表された方向指数である。更に、ｘ方向、ｙ方向の境界条件は、周期境界条件とした。
【００３７】
ヤング率Ｅは、引っ張り試験シミュレーションにより、図３に示した直方体状の結晶を構成する銅原子の最上層をｚ方向に引っ張ったときの、直方体状の結晶全体のｚ方向への伸びから計算した。
【００３８】
図４は、第１のシミュレーションの結果を示すグラフである。グラフの横軸は、ポテンシャルパラメータＤを単位「ｅＶ」で表し、縦軸は、ヤング率Ｅを単位「ＧＰａ」で表す。
【００３９】
本図に示す結果から、ヤング率ＥがポテンシャルパラメータＤに比例することが確認された。
【００４０】
図５及び図６を参照して第２のシミュレーションについて説明する。第２のシミュレーションは、ポテンシャルパラメータｒ０及びＡと、ヤング率Ｅとの関係を確認するシミュレーションである。
【００４１】
第２のシミュレーションにおいては、ポテンシャルパラメータＤを一定とし、ポテンシャルパラメータｒ０及びＡを、ｒ０とＡとの積を一定としつつ変化させて、ヤング率Ｅの変化を調べた。すなわち、ｒ０をｘｒ０、ＡをＡ／ｘとし、パラメータｘを変化させることでｘｒ０及びＡ／ｘを変化させて、ヤング率Ｅとｘとの関係を調べた。
【００４２】
なお、銅のＭｏｒｓｅポテンシャルのポテンシャルパラメータ、ｒ０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）の値はそのまま用いた。
【００４３】
図５は、ポテンシャルパラメータを変化させる前後のポテンシャルφを示すグラフである。グラフの横軸は、基点原子（ｒ＝０の位置に存在する原子）からの距離ｒを単位「ｍ」で表し、グラフの縦軸はポテンシャルφを単位「ｅＶ」で表す。
【００４４】
曲線ｐ（ｒ）はパラメータを変化させない場合（ｘ＝１の場合）の基点原子からの距離ｒとポテンシャルとの関係を示す。曲線ｑ（ｒ）はパラメータを変化させた場合の基点原子からの距離ｒとポテンシャルとの関係を示す。曲線ｑ（ｒ）は、ｘ＝２としたとき、すなわちｒ０を２ｒ０、ＡをＡ／２としたときの関係を示すグラフである。
【００４５】
Ｐ１、Ｐ２、Ｐ３はそれぞれポテンシャルパラメータを変化させない場合（ｘ＝１の場合）における、基点原子に最も近い原子、２番めに近い原子、３番めに近い原子の原子位置（基点原子からの距離）を表す。
【００４６】
また、Ｑ１、Ｑ２、Ｑ３はそれぞれｘ＝２としてポテンシャルパラメータを変化させた場合における、基点原子に最も近い原子、２番めに近い原子、３番めに近い原子の原子位置（基点原子からの距離）を表す。
【００４７】
パラメータｒ０をｘｒ０、ＡをＡ／ｘと変化させることは、１次元的に考えたとき、ｘの値にかかわらず、基点原子にｔ番め（ｔ＝１、２、３）に近い原子のポテンシャル値を、変化前後でそれぞれ常に同一に保つ、という意味を有する。たとえば、図５においては、ｐ（Ｐ１）＝ｑ（Ｑ１）、ｐ（Ｐ２）＝ｑ（Ｑ２）、ｐ（Ｐ３）＝ｑ（Ｑ３）の関係が満たされるように、ポテンシャルパラメータｒ０及びＡは変化されることになる。また、ｘの値にかかわらず、変化前のｒにおけるポテンシャル値と、変化後のｘｒにおけるポテンシャル値とを同一とすることを意味する。
【００４８】
第２のシミュレーションにおいても、第１のシミュレーションと同様の引っ張り試験シミュレーションモデルを用いてシミュレーションを行った。
【００４９】
図６は、第２のシミュレーションの結果を示すグラフである。グラフの横軸は１／ｘ３を表し、縦軸はヤング率Ｅを単位「ＧＰａ」で表す。
【００５０】
本図に示す結果から、ヤング率Ｅが１／ｘ３に比例することが確認される。このことはヤング率Ｅが、Ａ２／ｒ０に比例することをも意味する。
【００５１】
最後に、第３のシミュレーションについて説明する。第３のシミュレーションにおいては、ｒ０をｘｒ０、ＡをＡ／ｘ、Ｄをｘ３Ｄとして、パラメータｘを変化させ、ヤング率Ｅとｘとの関係を調べた。
【００５２】
第３のシミュレーションにおいても、第１及び第２のシミュレーションと同様の引っ張り試験シミュレーションモデルを用いた。また、第２のシミュレーション同様、銅のＭｏｒｓｅポテンシャルのポテンシャルパラメータ、ｒ０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）の値はそのまま用いた。
【００５３】
図７は、第３のシミュレーションの結果を示すグラフである。グラフの横軸はｘを表し、縦軸はヤング率Ｅを単位「ＧＰａ」で表す。グラフには、Ｗ、Ｃｒ、Ｆｅ、Ｐｔ、Ｃｕ、Ａｌ、及びＰｂの各金属のヤング率を書き入れてある。
【００５４】
本図に示す結果から、ヤング率Ｅはｘの値に関らず、ほぼ一定に保たれることが確認される。第３のシミュレーションで計算されたヤング率Ｅの値は、実際の銅（Ｃｕ）のヤング率に接近し、ＰｔやＡｌ、その他の金属のヤング率と混同を生じさせるものではないこともわかる。また、第３のシミュレーションの結果は、ヤング率Ｅが、ＤＡ２／ｒ０に比例（比例定数は約１）することをも示している。
【００５５】
第３のシミュレーションの結果から、ヤング率を一定としたまま、ｒ０、Ａ、Ｄのポテンシャルパラメータを変更して種々のシミュレーションを行うことができることが理解される。たとえばヤング率を一定に保ったままｒ０を大きくし、種々のシミュレーションを行うことができる。
【００５６】
図８を参照して、実施例によるシミュレーション方法を説明する。
【００５７】
ステップＳ１０１において、シミュレーション対象の材料を構成する粒子の二体間ポテンシャルφを、ポテンシャルが最小となる距離に対応するパラメータｒ０、ポテンシャルの形状を規定するパラメータＣ０を用い、２つの粒子の中心間距離ｒの関数φ（ｒ，ｒ０，Ｃ０）として決定する。
【００５８】
たとえばシミュレーション対象の材料が銅である場合、二体間ポテンシャルφは、
2010015284.tif 000017
【００５９】
・・（１１）
で表され、ｒ０＝０．２８６６（ｎｍ）、Ｄ＝０．３４２９（ｅＶ）、Ａ＝１３．５８８（１／ｎｍ）である。
【００６０】
次に、ステップＳ１０２において、パラメータｒ０を、ｒ０とは異なるｒ１に置き換えたとき、シミュレーション対象の材料、ここでは銅のヤング率の変化が抑制されるように、パラメータＣ０をＣ１に置き換える。
【００６１】
たとえば最安定原子間距離ｒ０を１０倍として、シミュレーションを行うのであれば、ｒ０を１０＊ｒ０に置き換える。この置き換えに伴って、ＡをＡ／１０、Ｄを１０００＊Ｄに置き換える。
【００６２】
そして、ステップＳ１０３において、新たなパラメータｒ１（実施例においては、１０＊ｒ０）とＣ１（実施例においては、Ａ／１０、及び、１０００＊Ｄ）とで定義される二体間ポテンシャルφ（ｒ，ｒ１，Ｃ１）を用いてシミュレーションを行う。
【００６３】
実施例によるシミュレーション方法によれば、高精度のシミュレーションを実現することができる。任意のヤング率に対し、ヤング率を変化させずにシミュレーションを行うことも可能である。
【００６４】
また、最安定原子間距離を任意に設定することができるため、多様なシミュレーションを実現することができる。
【００６５】
シミュレーションに要する種々の計算は、プログラムにより、コンピュータで実行することができる。
【００６６】
図９（Ａ）は、実施例によるシミュレーション方法を実行するシミュレーション実行装置のシステム構成図であり、図９（Ｂ）は、当該装置により実行される処理のフローである。
【００６７】
まず、ステップＳ２０１において、キーボードなどの入力装置から、シミュレーション対象の材料を構成する粒子の二体間ポテンシャルが最小となる距離に対応するパラメータｒ０を、ｒ０とは異なるｒ１に置き換えるための情報を入力する。
【００６８】
たとえば、シミュレーション対象の材料が銅であること、及び、銅の最安定原子間距離ｒ０を１０倍にしてシミュレーションを行うことを入力する。
【００６９】
ここでファイル装置のデータファイルには、たとえば銅の二体間ポテンシャルを示すＭｏｒｓｅポテンシャルにおいては、ｒ０＝０．２８６６（ｎｍ）、Ａ＝１３．５８８（１／ｎｍ）、Ｄ＝０．３４２９（ｅＶ）であることが記録されている。
【００７０】
続いて、ステップＳ２０２において、ｒ０を１０＊ｒ０に置き換えたとき、シミュレーション対象の材料（銅）のヤング率の変化が抑制されるように、Ｍｏｒｓｅポテンシャルの形状を規定するパラメータＡ及びＤを置き換える。
【００７１】
たとえば中央処理装置は、メインメモリ中の制御プログラムの指令を受け、データファイルから上記数値を読み出し、シミュレーションを行うに当たっては、ｒ０＝２．８６６（ｎｍ）、Ａ＝１．３５８８（１／ｎｍ）、Ｄ＝３４２．９（ｅＶ）とすることを算出する。
【００７２】
そして、ステップＳ２０３において、中央処理装置は、メインメモリ中の制御プログラムの指令を受け、銅で形成された対象につき、ｒ０＝２．８６６（ｎｍ）、Ｄ＝３４２．９（ｅＶ）、Ａ＝１．３５８８（１／ｎｍ）として、メインメモリに保存されているシミュレーション実行プログラムを実行し、シミュレーションを行う。シミュレーションの結果は、ディスプレイなどの出力装置を介して表示される。
【００７３】
上述の説明においては、ステップＳ２０１において、シミュレーション対象の材料が銅であること、及び、銅の最安定原子間距離ｒ０を１０倍にしてシミュレーションを行うことを入力したが、たとえば直接、ｒ０（ｒ１）を２．８６６（ｎｍ）としてシミュレーションを行うことが入力されてもよい。
【００７４】
またステップＳ２０２において、Ａ＝１．３５８８（１／ｎｍ）、Ｄ＝３４２．９（ｅＶ）とすることを算出したが、たとえばファイル装置に変換テーブルを記憶させておき、演算ではなく、ｒ０（ｒ１）の値をもとに、参照してこれらの値を得ることも可能である。
【００７５】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。
【００７６】
たとえばＭｏｒｓｅポテンシャルの場合について、ヤング率Ｅとポテンシャルパラメータとの関係を検証したが、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャルその他の場合についても、同様に確認することが可能であろう。
【００７７】
このため、Ｍｏｒｓｅポテンシャル以外のポテンシャルで表される材料についても、たとえばヤング率を一定としたまま、ポテンシャルパラメータを変更して種々のシミュレーションを行うことができる。
【００７８】
その他、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【産業上の利用可能性】
【００７９】
様々なシミュレーション、殊に分子動力学シミュレーションに好適に利用することができる。
【図面の簡単な説明】
【００８０】
【図１】二体間ポテンシャルを示す概略的なグラフである。
【図２】原子間のクロスセクションを示す図である。
【図３】第１のシミュレーションモデルを示す。
【図４】第１のシミュレーションの結果を示すグラフである。
【図５】ポテンシャルパラメータを変化させる前後のポテンシャルφを示すグラフである。
【図６】第２のシミュレーションの結果を示すグラフである。
【図７】第３のシミュレーションの結果を示すグラフである。
【図８】実施例によるシミュレーション方法を説明するためのフローチャートである。
【図９】（Ａ）は、シミュレーション実行装置のシステム構成図であり、（Ｂ）は、当該装置により実行される処理のフローである。
【符号の説明】
【００８１】
Ｓ１０１、Ｓ１０２、Ｓ１０３、Ｓ２０１、Ｓ２０２、Ｓ２０３  ステップ
【図１】
2010015284.tif 000020
【図３】
2010015284.tif 000021
【図５】
2010015284.tif 000022
【図９】
2010015284.tif 000023
【図２】
2010015284.tif 000024
【図４】
2010015284.tif 000025
【図６】
2010015284.tif 000026
【図７】
2010015284.tif 000027
【図８】
2010015284.tif 000028