JP2007004530A - 分子シミュレーション方法 - Google Patents

Abstract

【課題】 大規模な体系に対して、高精度で高速に現象解析を行うことが可能な分子シミュレーション方法を提供する。
【解決手段】 所定数の原子によって構成される第１原子集団と、第１原子集団の周囲に位置する第２原子集団と、第２原子集団の周囲に位置する第３原子集団と、から構成されるモデルに対して電子状態計算を行うことにより、原子の位置とエネルギー又は力との関係を示すポテンシャルマップを作成する。そして、現象解析手順では、ポテンシャルマップを参照して現象解析を行う。上記の分子シミュレーション方法によれば、ポテンシャルマップを用いることにより原子のエネルギーなどの算出に要する計算を削減することができるため、計算時間を短縮することができる。更に、ポテンシャルマップは、遠方に位置する原子の効果を平均的に考慮しつつ、電子状態計算を適切に用いることによって作成されているため、高精度に現象解析を行うことができる。
【選択図】 図５

Description

本発明は、分子シミュレーション方法に関する。

従来より、複数の原子から構成される体系の相転移や吸脱着などの現象解析を行う分子シミュレーション方法が知られている。例えば、特許文献１には、半導体装置の結晶シリコン中の点欠陥の拡散経路及び障壁エネルギーを算出する分子構造遷移シミュレーション方法が記載されている。更に、特許文献２乃至６には、分子動力学法などを用いた分子シミュレーション方法が記載されている。

特開２００２−２６０９７５号公報 特許２９５０２５３号公報 特許３１１０３１０号公報 特許２７９２４８７号公報 特開平５−２７４２７７号公報 特開平５−１５１１９０号公報

しかしながら、上記した特許文献１乃至６に記載された技術では、原子間ポテンシャルの存在や精度の面から応用範囲が制限されてしまう場合や、計算負荷が膨大となるためにモデルの規模を小さくした短時間の局所的な現象解析しか実行できない場合があった。

本発明は、このような問題点を解決するためになされたもので、その目的とするところは、大規模な体系に対して、高精度で高速に現象解析を行うことが可能な分子シミュレーション方法を提供することにある。

本発明の１つの観点では、所定数の原子によって構成される第１原子集団と、前記第１原子集団が存在する領域の周囲に位置する第２原子集団と、前記第１原子集団及び前記第２原子集団が存在する領域の周囲に位置する第３原子集団と、から構成されるモデルに対して電子状態計算を行うことにより、原子の位置とエネルギー又は力との関係を示すポテンシャルマップを作成するポテンシャルマップ作成手順と、作成された前記ポテンシャルマップを参照して、現象解析を行う現象解析手順と、を備える。

上記の分子シミュレーション方法は、複数の原子から構成される体系の相転移や吸脱着などの現象解析を行う場合などに好適に使用される。ポテンシャルマップ作成手順では、第１原子集団と第２原子集団と第３原子集団から構成されるモデルに対して電子状態計算を行うことにより、原子の位置とエネルギー又は力との関係を示すポテンシャルマップを作成する。そして、現象解析手順では、作成されたポテンシャルマップを参照して現象解析を行う。具体的には、第１原子集団は、所定数の原子数で構成され、モデルにおける中心部分に概ね位置する。第２原子集団は、第１原子集団が存在する領域の周囲に位置する。詳しくは、第２原子集団は、第１原子集団の中心原子から少なくとも所定距離離れた場所に位置する。一方、第３原子集団は、第１原子集団及び第２原子集団の存在する領域の周囲に位置する。以上より、上記の分子シミュレーション方法によれば、ポテンシャルマップを用いることにより原子のエネルギーなどに要する計算を削減することができるため、計算時間を短縮することができる。更に、ポテンシャルマップは、遠方に位置する原子の効果を平均的に考慮しつつ、電子状態計算を適切に用いることによって作成されているため、高精度に現象解析を行うことができる。したがって、上記の分子シミュレーション方法によれば、大規模な体系に対して、高精度で高速に現象解析を行うことが可能となる。

好適な１つの実施例では、前記第１原子集団及び前記第３原子集団は固定され、前記第２原子集団は、前記第１原子集団及び前記第３原子集団の影響によって移動させる。

この場合、第１原子集団を構成する原子は、ポテンシャルマップの各離散点に位置させるため、配置位置は固定する。一方、第２原子集団を構成する原子は、ポテンシャルマップにおける第１原子集団などの配置位置に応じて移動させる。即ち、第２原子集団を構成する原子は、ポテンシャルマップ作成の手順において、配置位置の最適化が行われる。一方、第３原子集団を構成する原子は、周囲の原子に対して影響を与えるが、配置位置は移動させない。

上記の分子シミュレーション方法の一態様では、前記電子状態計算は、バンド計算又は分子軌道計算によって行う。

この態様では、ポテンシャルマップ作成手順では、バンド計算又は分子軌道計算を行うことによってポテンシャルマップを作成する。これにより、反応などによって原子配置の組み換えによって電子状態が変化しても、古典的ポテンシャルとして扱うことが可能となる。

上記の分子シミュレーション方法の他の一態様では、前記現象解析手順は、前記ポテンシャルマップをスプライン補間することによって原子のエネルギー又は力を求める手順を有している。

この態様では、前記現象解析手順は、ポテンシャルマップから得られるエネルギー（ポテンシャルエネルギー）や力をスプライン補間して用いる。ポテンシャルマップを離散点によって作成した場合には、ポテンシャルエネルギーなどを求めるべき対象となっている原子は、必ずしも、この離散点上には位置しない。そのため、この原子の近傍の位置におけるポテンシャルエネルギーをポテンシャルマップから読み出し、読み出したポテンシャルエネルギーをスプライン補間した値を、求めるべきポテンシャルエネルギーとする。これにより、離散点によりポテンシャルマップを表すことによってデータ数の削減を図っても、現象解析の精度を確保することができる。

以下、図面を参照して本発明の好適な実施の形態について説明する。

［ポテンシャルマップ作成方法］
まず、本発明の実施形態に係るポテンシャルマップの作成方法について説明する。

図１は、ポテンシャルマップ作成に用いる電子状態計算モデルの基本概念を説明するための図である。

このモデルは、図１中の立方体の概ね中心部分に位置する原子集団（以下、「第１原子集団１」と呼ぶ。）と、第１原子集団１の存在する領域の周辺に位置する原子集団（以下、「第２原子集団２」と呼ぶ。）と、第２原子集団２の存在する領域の周辺に位置する原子集団（以下、「第３原子集団３」と呼ぶ。）によって構成される。第１原子集団１は、結晶の配位数程度の原子数によって構成され（図１の例では、４個）、後述するポテンシャルマップにおける各離散点に位置する。言い換えると、第１原子集団１を中心にしてポテンシャルマップを作成する。また、第２原子集団２は、第１原子集団１の中心原子から、結晶の第２又は第３近接程度離れた距離に位置する原子集団によって構成される。一方、第３原子集団３は、第１原子集団１の中心原子から第２又は第３近接程度の距離よりも更に離れた距離に位置する原子集団によって構成される。

本実施形態では、上記のような原子集団によって構成されるモデルに対して電子状態計算を行うことによりポテンシャルマップを作成する。具体的には、第１原子集団１を構成する原子は、ポテンシャルマップの各離散点に位置させるため、配置位置は固定する。即ち、第１原子集団１は、周囲の原子の力などによって移動させない。一方、第２原子集団２を構成する原子は、ポテンシャルマップにおける第１原子集団１などの配置位置に応じて移動させる。即ち、第２原子集団２を構成する原子は、周囲の原子の力などを受けて移動し、ポテンシャルマップ作成の手順において、配置位置の最適化が行われる。一方、第３原子集団３を構成する原子は、周囲の原子に対して影響を与えるが、配置位置は移動させない。

ここで、上記した電子状態計算の具体例について、図２を用いて説明する。

図２（ａ）は、バンド計算を説明するための図である。破線１１で示す領域内（以下、この領域を「クラスター」とも呼ぶ。）には第１原子集団１が位置する。また、破線１１で示す領域を除く、破線１２で示す領域内には第２原子集団２が位置する。この場合、第２原子集団２は、破線１５で示すような結晶格子を形成する。一方、破線１１及び破線１２で示す領域を除く、実線１３で示す領域内には第３原子集団３が位置する。バンド計算によれば、種々の原子配置について周期境界条件を課して力やエネルギーなどを計算することにより全エネルギーを求めることができる。この場合、実線１３で囲まれた領域内の全ての原子が、バンド計算におけるスーパーセルを構成する。即ち、結晶格子１５を適当に積層したスーパーセル１３の中心部分をクラスター１１に置き換えたユニットに対して周期境界条件を課して、バンド計算を行う。そして、求められた全エネルギーを第１原子集団１の原子配置に対してマップ化することにより、ポテンシャルマップが得られる。

図２（ｂ）は、分子軌道計算（ＭＯ；ＭＯＰＡＣ）を説明するための図である。図２（ｂ）においても、図２（ａ）と同様な配置位置に、第１原子集団１、第２原子集団２、及び第３原子集団３は位置している。分子軌道計算においては、周期境界条件を用いずに、上記したユニットを孤立分子と見なして全エネルギーを求める点でバンド計算とは異なる。この場合も、求められた全エネルギーを第１原子集団１の原子配置に対してマップ化することにより、ポテンシャルマップが得られる。

次に、ポテンシャルマップ作成方法について、図３に示すフローチャートを用いて説明する。ここでは、一例として、昇華法によるＳｉＣの単結晶成長についてのポテンシャルマップの作成方法について説明する。なお、ポテンシャルマップ作成方法は、パーソナルコンピュータなどによって実行される。具体的には、パーソナルコンピュータは、所定のプログラムを実行することによってポテンシャルマップを作成する。

まず、ステップＳ１０１では、クラスターに含まれる原子数Ｎ（即ち、第１原子集団１の数）と、クラスターサイズ（力の仕切り距離）を設定する。例えば、ガスと結晶を有するモデルが採用された場合、ガス相の原子数Ｎは「３」以下であり、結晶層の原子数Ｎは「３」より大きく「６」以下である。ステップＳ１０１の設定が終了すると、ステップＳ１０２の手順に進む。ステップＳ１０２では、シミュレーションを行う系に含まれる元素種を設定して、ステップＳ１０３に進む。

ステップＳ１０３では、原子数Ｎが「３」より大きいか否かの判定を行う。このような判定を行うのは、作成するポテンシャルマップの対象がガス相である場合（Ｎ≦３）と、結晶相である場合（Ｎ＞３）とを区別するためである。原子数Ｎが「３」より大きい場合（ステップＳ１０３；Ｙｅｓ）には、ステップＳ１２０に進み、原子数Ｎが「３」以下である場合（ステップＳ１０３；Ｎｏ）には、ステップＳ１０４に進む。

ステップＳ１０４では、変数ｉに「０」を代入し、ステップＳ１０５に進む。ステップＳ１０５では、原子数Ｎが「２」以上である原子集団において規定される原子間距離ｒを、演算式「ｒ＝ｒ０＋ｒｄ×ＲＥＡＬ(ｉ)」によって得る。この場合、「ｒ０」は最小原子間距離を示し、「ｒｄ」は更新する際の間隔を示している。また、「ＲＥＡＬ(ｉ)」は、変数ｉを実数に変換する演算式を示している。ステップＳ１０５の手順が終了すると、ステップＳ１０６に進む。なお、ステップＳ１０５の手順は、原子間距離ｒが最大原子間距離となるまで繰り返し行われる。

ステップＳ１０６では、原子数Ｎが「３」であるか否かの判定を行う。原子数Ｎが「３」である場合（ステップＳ１０６；Ｙｅｓ）には、ステップＳ１０８に進む。この場合には、３つの原子間に働く３体力を求める必要がある。一方、原子数Ｎが「３」ではない場合（ステップＳ１０６；Ｎｏ）、即ち原子数Ｎが「２」である場合には、ステップＳ１０７に進む。この場合には、２つの原子間に働く２体力を求める必要がある。

ステップＳ１０７では、分子軌道計算（ＭＯ）を用いて、原子数Ｎが「２」である場合の原子間相互作用エネルギーＶ_２（ｉ）を算出する（以下では、原子間相互作用エネルギーを、「ポテンシャルエネルギー」、又は単に「ポテンシャル」若しくは「エネルギー」とも呼ぶ。）。この場合、原子間距離ｒが最小原子間距離ｒ０から最大原子間距離となるまで変数ｉを更新している際に、その変数ｉに対応する原子間相互作用エネルギーＶ_２（ｉ）を繰り返し算出する。ステップＳ１０７の手順が終了すると、ステップＳ１１７に進む。

一方、ステップＳ１０８では、変数ｊを「０」に設定し、ステップＳ１０９に進む。ステップＳ１０９では、原子数Ｎが「３」以上である原子集団において規定される第３原子間距離ｓを、演算式「ｓ＝ｓ０＋ｓｄ×ＲＥＡＬ(ｊ)」によって得る。この場合、「ｓ０」は最小第３原子間距離を示し、「ｓｄ」は更新する際の間隔を示し、「ＲＥＡＬ(ｊ)」は変数ｊを実数に変換する演算式を示している。ステップＳ１０９の手順が終了すると、ステップＳ１１０に進む。なお、ステップＳ１０９の手順は、第３原子間距離ｒが最大第３原子間距離となるまで繰り返し行われる。

ステップＳ１１０では、変数ｋを「０」に設定し、ステップＳ１１１に進む。ステップＳ１１１では、原子数Ｎが「３」以上である原子集団において規定される、中心原子を挟む角度ｔ（以下、「原子間角度」と呼ぶ。）を演算式「ｔ＝ｔ０＋ｔｄ×ＲＥＡＬ(ｋ)」によって得る。この場合、「ｔ０」は最小原子間角度を示し、「ｔｄ」は更新する際の間隔を示し、「ＲＥＡＬ(ｋ)」は変数ｋを実数に変換する演算式を示している。ステップＳ１１１の手順が終了すると、ステップＳ１１２に進む。なお、ステップＳ１１１の手順は、原子間角度ｔが最大原子間角度となるまで繰り返し行われる。

ステップＳ１１２では、分子軌道計算（ＭＯ）を用いて、原子数Ｎが「３」である場合の原子間相互作用エネルギーＶ_３（ｉ、ｊ、ｋ）を算出する。なお、原子間距離ｒが最小原子間距離ｒ０から最大原子間距離となるまで、第３原子間距離ｓが最小第３原子間距離ｓ０から最大第３原子間距離となるまで、及び原子間角度ｔが最小原子間角度ｔ０から最大原子間角度となるまで、変数ｉ、ｊ、ｋを種々の値に更新した際に、これらの変数ｉ、ｊ、ｋに対応する原子間相互作用エネルギーＶ_３（ｉ、ｊ、ｋ）を繰り返し算出する。以上のステップＳ１１２の手順が終了すると、ステップＳ１１３に進む。

ステップＳ１１３では、変数ｋを「ｋ＋１」に更新し、ステップＳ１１４に進む。ステップＳ１１４では、変数ｋの更新が終了したか否かを判定する。具体的には、原子間角度ｔが最大原子間角度になったか否かを判定する。原子間角度ｔが最大原子間角度になっている場合（ステップＳ１１４；Ｙｅｓ）には、ステップＳ１１５に進む。一方、原子間角度ｔが最大原子間角度になっていない場合（ステップＳ１１４；Ｎｏ）には、ステップＳ１１１に戻る。

ステップＳ１１５では、変数ｊを「ｊ＋１」に更新し、ステップＳ１１６に進む。ステップＳ１１６では、変数ｊの更新が終了したか否かを判定する。具体的には、第３原子間距離ｓが最大第３原子間距離になったか否かを判定する。第３原子間距離ｓが最大第３原子間距離になっている場合（ステップＳ１１６；Ｙｅｓ）には、ステップＳ１１７に進む。一方、第３原子間距離ｓが最大第３原子間距離になっていない場合（ステップＳ１１６；Ｎｏ）には、ステップＳ１０９に戻る。

ステップＳ１１７では、変数ｉを「ｉ＋１」に更新し、ステップＳ１１８に進む。ステップＳ１１８では、変数ｉの更新が終了したか否かを判定する。具体的には、原子間距離ｒが最大原子間距離になったか否かを判定する。原子間距離ｒが最大原子間距離になっている場合（ステップＳ１１８；Ｙｅｓ）には、ステップＳ１２３に進む。一方、原子間距離ｒが最大原子間距離になっていない場合（ステップＳ１１８；Ｎｏ）には、ステップＳ１０５に戻る。

ステップＳ１２３では、得られたポテンシャルマップを保存する。具体的には、原子数Ｎが「２」である場合に種々の変数ｉを代入することによって得られた原子間相互作用エネルギーＶ_２（ｉ）、又は原子数Ｎが「３」である場合に種々の変数ｉ、ｊ、ｋを代入することによって得られた原子間相互作用エネルギーＶ_３（ｉ、ｊ、ｋ）がポテンシャルマップに対応する。以上の手順が終了すると、当該フローを抜ける。

一方、原子数Ｎが「３」より大きい場合、即ちポテンシャルマップの作成対象が結晶である場合には、ステップＳ１２０〜ステップＳ１２２の手順を行う。ステップＳ１２０では、前述したステップＳ１０５、Ｓ１０９、Ｓ１１１と同様に、原子間距離ｒや第３原子間距離ｓや原子間角度ｔなどを更新する手順を行う。そして、ステップＳ１２１の手順に進む。

ステップＳ１２１では、前述した電子状態計算のいずれかを選択し、原子数Ｎにおける原子間相互作用エネルギーＶ_Ｎ（ｉ、ｊ、ｋ、…）を算出する。具体的には、クラスターに含まれる複数の原子に対して中心原子を原点として極座標表示して、原子配置をあらゆる位置に移動させることによって、原子間相互作用エネルギーＶ_Ｎ（ｉ、ｊ、ｋ、…）を算出する。例えば、クラスターに含まれる原子数は、最大で５〜６個程度である。また、原子間相互作用エネルギーＶ_Ｎ（ｉ、ｊ、ｋ、…）は、前述したバンド計算や分子軌道計算などの電子状態計算によって計算することができる。ステップＳ１２１の手順が終了すると、ステップＳ１２２の手順に進む。

ステップＳ１２２では、原子間距離ｒや第３原子間距離ｓなどのパラメータが最大値に達するまで変数ｉ、ｊ、ｋ、…が更新されたか否かを判定している。更新が終了している場合（ステップＳ１２２；Ｙｅｓ）には、ステップＳ１２３の手順に進み、更新が終了していない場合（ステップＳ１２２；Ｎｏ）には、ステップＳ１２０の手順に戻る。

ステップＳ１２３では、得られたポテンシャルマップを保存する。具体的には、種々の変数ｉ、ｊ、ｋ、…を代入することによって得られた原子間相互作用エネルギーＶ_Ｎ（ｉ、ｊ、ｋ、…）をポテンシャルマップとして記憶する。以上の手順が終了すると、当該フローを抜ける。

ここで、作成されたポテンシャルマップの例を図４に示す。図４は、ＳｉＣ種結晶表面上のポテンシャルエネルギー場を示している。また、図４中の「Ｖ」はポテンシャルエネルギーを示している。具体的には、図４（ａ）は、３Ｃ−ＳｉＣの種結晶を示し、図４（ｂ）は、２Ｈ−ＳｉＣの種結晶を示している。これより、ポテンシャルマップ上に予想安定点（格子位置）などが現れていることがわかる。

このように、本実施形態に係るポテンシャルマップは、遠方に位置する原子の効果を平均的に考慮しつつ、電子状態計算を適切に用いることによって作成される。したがって、このポテンシャルマップによって決定されるポテンシャルは、成長させる結晶場で得られるポテンシャルに近い値となるとと共に、反応などによって原子配置の組み換えによって電子状態が変化しても、古典的ポテンシャルとして扱うことができる。

［現象解析方法］
次に、前述したポテンシャルマップを利用した現象解析方法について、図５に示すフローチャートを用いて説明する。ここでは、一例とし、昇華法によるＳｉＣの単結晶成長をシミュレーションした場合について説明する。なお、現象解析方法は、パーソナルコンピュータなどによって実行される。具体的には、パーソナルコンピュータは所定のプログラムを実行することにより、現象解析を行う。

まず、ステップＳ２０１では、解析する現象の条件（シミュレーション条件）を設定する。具体的には、ガス種の種類や、温度や、密度（分圧）や、供給条件や、種結晶の大きさや、方位や、温度条件などを設定する。ステップＳ２０１の設定が終了すると、ステップＳ２０２の手順に進む。

ステップＳ２０２では、現象解析の対象となっている原子の位置や速度などの初期状態を設定する。例えば、シミュレーションボックス（シミュレーションで設定しているモデルの単位）内の下方部分に、結晶系と表面の面方位を考慮してＳｉ原子、Ｃ原子を配置する。そして、シミュレーションボックスの上方部分に、所定の分圧に比例した密度の原料ガス分子を配置する。この場合、結晶相の原子とガス相の原子は、温度に応じた平均速度を有するように速度を設定する。また、シミュレーションボックスのサイズや、結晶相の温度や、ガス相の温度や、ガス分圧や、番号が付された各原子に対する元素種と座標と速度なども初期状態として設定する。以上のステップＳ２０２の手順が終了すると、ステップＳ２０３の手順に進む。

ステップＳ２０３では、現象解析の対象となっている全ての原子に対して、近傍に存在する原子の数、元素種、位置、速度をテーブルデータ（以下、「クラスタテーブル」と呼ぶ。）として記録する。このクラスタテーブルには、中心原子を原点に、最も近い原子をｚ軸に位置させ、次に近い原子をｘｚ面内に位置させた極座標によって表示した原子配置データが記録されている。更に、クラスタテーブルには、上記の極座標系と空間に固定した座標系との関係をオイラー角で示した変換行列によって表されたクラスター方位も記録されている。以上のステップＳ２０３の手順が終了すると、ステップＳ２０４の手順に進む。

ステップＳ２０４では、記憶しているポテンシャルマップを参照してポテンシャルエネルギーを読み出す。詳しくは、まず、ポテンシャルエネルギーの算出の対象となっている中心原子（以下、「対象原子」と呼ぶ。）の原子配置データをクラスタテーブルから読み出し、この原子配置データに最も近いポテンシャルマップにおける点（以下、「最近接点」と呼ぶ。）と、この最近接点に隣接する点（以下、「隣接点」と呼ぶ。）を決定する。そして、最近接点のポテンシャルエネルギーと、隣接点のポテンシャルエネルギーを読み出す。このような手順を行うのは、ポテンシャルマップは離散点によって規定されているため、対象原子が、必ずしも離散点上に位置しないためである。以上のステップＳ２０４の手順が終了すると、ステップＳ２０５の手順に進む。

ステップＳ２０５では、以下の手順を行う。まず、対象原子のポテンシャルエネルギーを、対象原子の原子配置データとポテンシャルマップ内の最近接点との距離に基づいて、ステップＳ２０４で決定した最近接点のポテンシャルエネルギーと隣接点のポテンシャルエネルギーとをスプライン補間することによって求める。更に、対象原子に働く力を、最近接点のポテンシャルエネルギーと隣接点のポテンシャルエネルギーとの差分を求め、求めた差分を各成分に分解した値に対してスプライン補間することによって求める。このように、スプライン補間によりポテンシャルエネルギーを求めることにより、ポテンシャルマップのデータ数の削減を図っても、現象解析の精度を確保することが可能となる。なお、ポテンシャルエネルギーを求める際にポテンシャルマップを補間する方法は、スプライン補間に限定はされない。また、計算速度の向上のために、補間を行わずに、最近接点のポテンシャルエネルギーをそのまま用いてもよい。

次に、対象原子の位置及び速度の変化を、運動方程式を差分化した式（例えばベルレの運動方程式）に、上記のようにして得られた対象原子に働く力を代入することによって求める。そして、前回ステップＳ２０５の手順において得られた値に対して、求められた対象原子の位置及び速度の変化を加算することによって、対象原子の位置と速度を更新する。なお、上記したエネルギー、力、位置、及び速度の計算は、系に含まれる全原子に対して行う。

更に、ステップＳ２０５では、ガス相に存在する全原子についての運動エネルギーの平均値及びビリアルを計算することによって、ガス相の温度と圧力を求める。これと同時に、結晶相に存在する全原子についての運動エネルギーの平均値を計算することによって、結晶相の温度を求める。以上のステップＳ２０５の手順が終了すると、ステップＳ２０６の手順を行う。

ステップＳ２０６の手順では、ステップＳ２０３〜Ｓ２０５までの手順を繰り返し行ったステップ数が所定のサイクル（補正サイクル）に達したか否かを判定する。例えば、補正サイクルは、温度や圧力などの熱力学データの補正を行うまでに繰り返し実行する回数に対応する。ステップ数が補正サイクルに達している場合（ステップＳ２０６；Ｙｅｓ）には、ステップＳ２０７の手順に進み、ステップ数が補正サイクルに達していない場合（ステップＳ２０６；Ｎｏ）には、ステップＳ２０３の手順に戻る。

ステップＳ２０７では、温度・圧力の平均値を更新すると共に、温度・圧力を補正する。具体的には、系の温度・圧力は所定の期間で更新し、この更新を行う間においては、温度・圧力の平均値を算出する。また、更新する前に、系の温度及び圧力と、ステップＳ２０１、Ｓ２０２などで設定した設定条件とを比較し、差がある場合には適切な方法によって補正を行う。以上のステップＳ２０７の手順が終了すると、ステップＳ２０８に進む。

ステップＳ２０８の手順では、ステップＳ２０３〜Ｓ２０７までの手順を繰り返し行ったステップ数が所定のステップ（最終ステップ）に達したか否かを判定する。この最終ステップは、シミュレーションを行うべき時間（即ち、繰り返し実行回数）などに対応する。ステップ数が最終ステップに達している場合（ステップＳ２０８；Ｙｅｓ）にはステップＳ２０９の手順に進み、ステップ数が最終ステップに達していない場合（ステップＳ２０８；Ｎｏ）には、ステップＳ２０３の手順に戻る。

ステップＳ２０９では、上記のシミュレーションを行うことによって得られた現象解析ファイルを作成する。具体的には、シミュレーションの実行時に随時記録された原子の位置・速度などのデータに基づいて、シミュレーションボックスのスナップショットの経時変化などの現象解析ファイルを作成する。以上の手順が終了すると、当該フローを抜ける。

このように、本実施形態に係る現象解析方法では、予め作成されたポテンシャルマップを用いて現象解析を行う。これにより、系の状態が変化するごとに原子のポテンシャルエネルギーなどをモンテカルロ法（ＭＣ）や分子動力学法などによって計算する必要がないため、計算時間を短縮することができる。更に、ポテンシャルマップは、遠方に位置する原子の効果を平均的に考慮しつつ、電子状態計算を適切に用いることによって作成されているため、ポテンシャルマップを用いることにより、高精度に現象解析を行うことができる。以上より、本実施形態に係る現象解析方法によれば、大規模な体系に対して、高精度で高速に現象解析を行うことが可能となる。

ポテンシャルマップ作成に用いる電子状態計算モデルの基本概念を説明するための図である。 電子状態計算の具体例を説明するための図である。 本実施形態に係るポテンシャルマップの作成方法を示すフローチャートである。 ポテンシャルマップの具体例を示す図である。 本実施形態に係る現象解析方法を示すフローチャートである。

符号の説明

１ 第１原子集団
２ 第２原子集団
３ 第３原子集団
１１ クラスター
１３ スーパーセル
１５ 結晶格子

Claims (4)

1. 所定数の原子によって構成される第１原子集団と、前記第１原子集団が存在する領域の周囲に位置する第２原子集団と、前記第１原子集団及び前記第２原子集団が存在する領域の周囲に位置する第３原子集団と、から構成されるモデルに対して電子状態計算を行うことにより、原子の位置とエネルギー又は力との関係を示すポテンシャルマップを作成するポテンシャルマップ作成手順と、
   作成された前記ポテンシャルマップを参照して、現象解析を行う現象解析手順と、を備えることを特徴とする分子シミュレーション方法。
2. 前記第１原子集団及び前記第３原子集団は固定され、
   前記第２原子集団は、前記第１原子集団及び前記第３原子集団の影響によって移動させることを特徴とする請求項１に記載の分子シミュレーション方法。
3. 前記電子状態計算は、バンド計算又は分子軌道計算によって行うことを特徴とする請求項１又は２に記載の分子シミュレーション方法。
4. 前記現象解析手順は、前記ポテンシャルマップをスプライン補間することによって原子のエネルギー又は力を求める手順を有していることを特徴とする請求項１乃至３のいずれか一項に記載の分子シミュレーション方法。
   

Images