シミュレーション装置、シミュレーション方法、シミュレーションプログラムを格納した記録媒体
【課題】数千個以上の原子集団を扱う大規模系の分子動力学シミュレーションにおいて、多元素が混在する複雑な組成を有する物質の、共有結合の組み換えを伴う構造変化を、効率的に再現する。
【解決手段】共有結合部位を決定する役割を持つ内部自由度を持った原子モデルを考え、原子核座標と内部自由度を併せた拡張力学系の運動を古典分子動力学法で計算し、結合の滑らかな繋ぎ変えを表現する。共有結合のような方位依存性のある引力相互作用の発現を、上記の方法で所望の原子間に限定することにより、従来法で問題となっていた企図しない余計な相互作用を取り除き、多体相互作用項同士の干渉を大幅に減少させることができる。これにより、多体相互作用項のポテンシャル定数の決定が大幅に単純化され、多元素が混在する複雑な系のポテンシャル定数を用意することが可能となる。従来法では不可能であった、多元素混在系において、共有結合の組み換えを伴う化学変化の大規模分子シミュレーションが、本発明によって可能となる。
【解決手段】共有結合部位を決定する役割を持つ内部自由度を持った原子モデルを考え、原子核座標と内部自由度を併せた拡張力学系の運動を古典分子動力学法で計算し、結合の滑らかな繋ぎ変えを表現する。共有結合のような方位依存性のある引力相互作用の発現を、上記の方法で所望の原子間に限定することにより、従来法で問題となっていた企図しない余計な相互作用を取り除き、多体相互作用項同士の干渉を大幅に減少させることができる。これにより、多体相互作用項のポテンシャル定数の決定が大幅に単純化され、多元素が混在する複雑な系のポテンシャル定数を用意することが可能となる。従来法では不可能であった、多元素混在系において、共有結合の組み換えを伴う化学変化の大規模分子シミュレーションが、本発明によって可能となる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は半導体集積回路製作等に用いる材料設計支援シミュレーション技術、加工プロセスのシミュレーション技術、あるいは、蛋白質等の生体分子やその他有機化合物質の分子シミュレーション技術に係り、特にナノメートルスケールのサイズを有する数千原子規模以上の系を対象とする大規模分子シミュレーション方法、シミュレーション装置、シミュレーションプログラムを格納した記録媒体に関する。
【背景技術】
【０００２】
半導体集積回路の高集積化に伴い、ナノメートルオーダーの極めて高い精度の微細加工技術が要求されている。例えば、金属−絶縁膜−半導体型電界効果トランジスタの絶縁膜部分は僅か1ナノメートルまで薄層化されている。これはシリカの構成単位であるSiO4正四面体構造の４層程度の厚さであり、積層構造を１原子層レベルで再現性良く制御する技術が要求される。ここまで微細化が進行すると、連続的な拡散方程式を解く従来型のプロセス・シミュレーションでは、界面形状や組成分布の予測能力に限界があり、原子１個１個の挙動を視野に入れた分子シミュレーション技術の導入が必要となっている。
【０００３】
一方創薬分野においては、分子シミュレーションは既に開発現場で欠かせない技術となっている。薬剤候補となる膨大な数の低分子化合物と蛋白質の親和性を計算機で予測することで、候補薬剤の絞込み、化合物最適化の高効率化に重要な役割を果たしている。
【０００４】
一般に全ての物質の性質は、電子および原子核間のクーロン相互作用を含む相対論的なディラック方程式、もしくは非相対論的なシュレディンガー方程式をパウリの排他原理の下で解くことで、基本的な物理定数のみから非経験的に予測することができる。原子核を固定した場合の近似的な数値解法としてハートリー・フォック（Hartree-Fock）方程式を解く分子軌道法や、コーン・シャム（Kohn-Sham）方程式を解く密度汎関数法が確立している。また、密度汎関数法で得られる全エネルギーを断熱ポテンシャルとして動力学計算するカー・パリネロ（Car-Parrinello）法も広く普及している。しかし、これら第一原理計算手法が取り扱える系の規模は、その計算量の多さ故、スーパーコンピュータを使用しても数百原子程度に限定される。１辺にたった１０個の原子を並べた立方晶ですら１０００原子に達することを鑑みても分かるように、第一原理計算がナノスケールの構造体の加工プロセスのシミュレータの役割を果たすことは期待できない。
【０００５】
数千原子以上の原子団の挙動を計算するには、原子間相互作用を単純なポテンシャル関数で近似し、その原子間ポテンシャルを用いてニュートンの運動方程式を数値的に解く古典的分子動力学法が有効である。電子波動関数の固有値問題を解かないこの方法では、扱う元素の組み合わせに応じて原子間ポテンシャル関数を設計する必要があるが、一度ポテンシャル関数が設定された系に対しては、コンピュータシステムにより、数十原子レベルのミクロな系から数十万原子レベルのナノスケール領域に至る広い原子数規模の範囲の動的シミュレーションを、短時間で、効率的に実行できる。
【０００６】
古典的分子動力学計算で用いられる原子間ポテンシャル関数は、大きく２種類に分類される。１つは、原子座標と共に共有結合のトポロジーを入力情報とするポテンシャルで、MM2やAMBERなどが知られている。パラメータの決定が比較的容易で多元素混在系が扱えることから、蛋白質分子の立体構造を基にした計算機創薬で中心的役割を果たしている。ただし、共有結合のトポロジーは固定されるため、例えば加水分解や脱水縮合によるペプチド結合の組み換えような化学変化は再現できない。
【０００７】
これに対し、原子座標のみを入力情報とするポテンシャルも数多く存在する。希ガスやイオン結合性の無機化合物を表現する単純な２体相互作用項からなるポテンシャルを始め、金属や有機化合物、半導体など、強い結合を含む物質を表現できる複雑な多体項を有するポテンシャルまで、様々な関数が提案されている。これらは結合箇所を露に指定しないため、結合組み換えを伴う様々な構造変化も扱える。半導体装置のプロセス・シミュレーションにおいては、結晶化やアモルファス化、表面や界面での再構成、欠陥の挙動など、共有結合の組み換えを伴う構造変化の再現が重要となるため、導入されつつある分子シミュレーション技術では主に複雑な多体相互作用ポテンシャル関数が用いられる。しかし、多元組成系に適用できるポテンシャル定数を揃えることが極めて難しく、複雑な多体相互作用ポテンシャルの場合、同時に扱える元素は通常２種類程度に限定される。
【０００８】
【非特許文献１】R. Car and M. Parrinello, Physical Review Letters Vol. 55, (1985) p.2471.
【非特許文献２】N. L. Allinger, Journal of American Chemical Society Vol. 99, (1977) p.8127.
【非特許文献３】S. J. Weiner, P. A. Kollman, D. T. Nguyen and D. A. Case, Journal of Computational Chemistry Vol. 4, (1983) p.187.
【非特許文献４】F. H. Stillinger and T. A. Weber, Physical Review B Vol.31, (1985), p.5262.
【非特許文献５】J. Tersoff, Physical Review B Vol.39, (1989) p.5566.
【非特許文献６】T. Watanabe, H. Fujiwara, H. Noguchi, T. Hoshino, and I. Ohdomari, Japanese Journal of Applied Physics Vol. 38, (1999) p.L366.
【非特許文献７】J. Dalla Torre, J.-L. Bocquet, Y. Limoge, J.-P. Crocombette, E. Adam, G. Martin, T. Baron, P. Rivallin and P. Mur, Journal of Applied Physics Vol. 92, (2002) p.1084.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
半導体集積回路のゲート絶縁膜はHfSiO2など、複雑な組成の高誘電率材料に置き換えられようとしており、半導体プロセスで用いられる材料の多様化が急速に進んでいる。このため、従来の結合組み換え型ポテンシャルでも、より多くの元素を同時に扱える必要が生じている。
【００１０】
結合組み替え型の多体相互作用ポテンシャルには、ポテンシャル定数の決定が困難であること以外に、もう一つ深刻な欠点を抱えている。多体相互作用ポテンシャルでは、３体項以上の多体項が配位数を調節する役割を果たしているが、これが企図しないポテンシャル障壁を形成する。このため、圧縮ストレスを受けた材料のシミュレーションでは材料の内部応力が１〜２桁も過大評価される場合がある。また、一般に３体項を含むポテンシャルを用いた融解シミュレーションでは、融点が実験値より高くなる傾向があるが、これも余計なポテンシャル障壁が融解を阻止していることが原因と考えられる。
【００１１】
ターゾフ（Tersoff）ポテンシャルなど、２体相互作用項に周辺の原子配置に応じて変化する結合次数を導入するボンドオーダー型ポテンシャルでは、３体項以上の多体項を導入することなく配位数調節が実現されているが、第一原理計算では見られない準安定構造が再現されるなど、複雑な結合次数項が原因と見られる副作用が、安定構造から大きく外れた領域で顕在化する。また、関数の複雑さ故、ポテンシャル定数と力定数の関係が明瞭でなく、多元素混在系への拡張を一層困難にしている。
【００１２】
多元素混在系に適用でき、かつ結合の組み換えが許されるようなポテンシャル関数を実現するには、原子核配座のみで一意にエネルギーを決定する従来法では限界がある。共有結合は、原子間で電子軌道の重なりが生じている部位に発現するものであり、原子間距離で一意に決まるものではない。従って、結合箇所を指定するため、電子の空間分布に相当する付加的な情報も必要と考えられる。
【００１３】
上記問題点に鑑み、本発明は、多元素混在系においても結合組み換えを伴う構造変化を正確に再現する、大規模分子動力学シミュレーション装置を提供することを目的とする。
【００１４】
具体的には、本発明は、従来の多体ポテンシャルが有していた余計なポテンシャル障壁を含むことなく、圧縮ストレスを帯びた物質の内部応力も正しく表現する、結合組み換え型の原子間相互作用モデルとその運動方程式を提供し、数千原子以上の多元素混在系においても、結合組み換えを伴う化学変化を再現するシミュレーション装置を提供することである。
【００１５】
本発明の更に他の目的は、数千原子以上の多元素混在系においても、結合組み換えを伴う化学変化を高速で計算するシミュレーション方法を提供することである。
【００１６】
本発明の更に他の目的は、数千原子以上の多元素混在系においても、結合組み換えを伴う化学変化を高速で計算するシミュレーションプログラムを格納した記録媒体を提供することである。
【課題を解決するための手段】
【００１７】
上記目標を達成するため、本発明の第１の主要な観点によれば、原子核座標及び内部座標の情報を格納する座標情報格納部と、前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価部と、前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部と、前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算部と、前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション部と、前記シミュレーション結果を出力する出力部とを有するシミュレーション装置が提供される。
【００１８】
また、この発明の一の実施形態によれば、前記シミュレーション部は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理部と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理部と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価部と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものである。また、この場合、前記シミュレーション実行部は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算部とを有するものであることが好ましい。
【００１９】
さらに、この発明の一の実施形態によれば、この装置は、前記出力部で出力されたデータの解析を行うデータ解析部をさらに有する。この場合、データ解析部は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価部を有するものであることが好ましい。また、このデータ解析部は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う構造解析部を有するものであることが好ましい。さらに、前記データ解析部は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算部を有するものであることが好ましい。さらに、前記データ解析部は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算部を有するものであることが好ましい。また、前記データ解析部は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析部を有するものであることが好ましい。また、前記データ解析部は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化部を有するものであることが好ましい。この場合、前記構造可視化部は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものである。
【００２０】
この発明の別の一の実施形態によれば、前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものである。この場合、前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものである。
【００２１】
さらに、この発明の別の一の実施形態によれば、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものである。
【００２２】
本発明によるシミュレーション装置によれば、各原子間の結合次数を、原子間距離と各原子の内部自由度変数の値から簡単に決定できるため、共有結合エネルギーを表す２体引力相互作用を、結合が組まれている原子間のみに発現させることができる。結合を組んでいない原子間には２体引力相互作用をはじめ、３体以上の多体相互作用項が働かないため、企図しないポテンシャル障壁の形成を防ぐことができる。
【００２３】
また、本発明によるシミュレーション装置によれば、結合を組んでいない原子間に共有結合を表す２体引力相互作用が適用されないため、余計な２体引力項を打ち消すように３体相互作用項を調節して配位数の調節を行う必要がなく、共有結合手の方向を決定する内部自由度の数そのもので、適切な配位数を保つことができる。このため、２体ポテンシャル定数と３体ポテンシャル定数を独立に決定することができ、多元素混在系のポテンシャル定数を用意することが可能となる。
【００２４】
また、本発明によるシミュレーション装置によれば、結合が途切れた３原子間に３体相互作用項を働かせないことができる。このため、結合が途切れた３原子間に不要に働く３体項を打ち消すようにポテンシャル定数を調節する必要がなく、異なる３体相互作用項のポテンシャル定数を独立に決定することができ、多元素混在系のポテンシャル定数を用意することが可能となる。
【００２５】
また、本発明によるシミュレーション装置によれば、原子核同様、各原子の共有結合手の方向を決定する内部自由度変数にも力が働き、その力に基づいて運動することができるため、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができる。これにより、シミュレーション中に共有結合の組み換えを連続的に起こすことができ、共有結合の組み換えを伴う動的過程を再現することができる。
【００２６】
本発明によるシミュレーション方法によれば、共有結合の所在を明確に示しながら、その組み換えを連続的に起こすことができ、共有結合の組み換えを伴う化学変化を再現することができる。
【００２７】
本発明によるシミュレーション方法によれば、内部座標情報管理部を実現するための記憶容量と、動的結合評価部の計算手続きが増える他は、従来の古典分子動力学計算と同程度の計算量となるため、数千個以上の原子から構成される大規模系の分子動力学計算が、短時間で効率的に実行できる。
【００２８】
本発明によるシミュレーション方法において、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度の相互作用が通常弱く、熱エネルギーの交換が起こりにくいため、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御することが好ましい。
【００２９】
また、本発明によるシミュレーション方法の１階微分計算部において、１階微分を計算する際、動的結合次数評価部における計算結果を繰り返し参照する必要があるため、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑えることが好ましい。
【００３０】
さらに、本発明によるシミュレーション方法において、構成元素の数を増やしていくほど、必要となるポテンシャル定数が増えていくが、シミュレーションに前もって全ての相互作用項のポテンシャル定数を用意する必要はない。シミュレーションで扱う系毎に、必要となる２体相互作用項や３体相互作用項などの種類が限定されることが多いため、これら最小限必要な相互作用項のポテンシャル定数のみで計算を実行しても良い。その場合、シミュレーション実行中に結合の組み換えが起こり、最初にポテンシャル定数を用意していなかった相互作用項が新たに必要となる場面が生じうるが、その際は、シミュレーションを一時中断し、新たに必要となったポテンシャル定数の入力を待って、シミュレーションを再開するようにしてもよい。
【００３１】
この発明の第２の主要な観点によれば、コンピュータシステムにより実行されるコンピュータソフトウェアプログラムであって、記録媒体と、前記記録媒体に格納され、前記コンピュータシステムに原子核座標及び内部座標の情報を格納させる座標情報格納部と、前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標及び内部座標に基いて原子間の結合次数を計算させる動的結合次数評価部と、前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算させるポテンシャルエネルギー計算部と、前記記録媒体に格納され、前記コンピュータシステムに前記ポテンシャルエネルギーを１階微分処理することによって各原子核座標および内部座標に作用する力を計算させる１階微分計算部と、前記記録媒体に格納され、前記コンピュータシステムに前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算させる圧力テンソル計算部と、前記記録媒体に格納され、前記コンピュータシステムに、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるシミュレーション部と、前記記録媒体に格納され、前記コンピュータシステムに前記シミュレーション結果を出力させる出力部とを有することを特徴とするコンピュータソフトウエアプログラムが提供される。
【００３２】
このような構成によれば、コンピュータシステムにインストールされることで、前記第１の主要な観点による装置を構成することができ、同様の作用効果をえることができる。
【００３３】
また、この発明の第３の主要な観点によれば、原子核座標及び内部座標の情報を格納する座標情報格納工程と、前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価工程と、前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算工程と、前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算工程と、前記微分計算工程で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算工程と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション工程と、前記シミュレーション結果を出力する出力工程と
を有することを特徴とする分子動力学シミュレーション方法が提供される。
【００３４】
このような構成によれば、前記第１の主要な観点によるシミュレーション装置で実行可能なシミュレーション方法が提供される。
【００３５】
なお、ここで、説明していない本発明の他の特徴及び顕著な作用及び効果は、次に説明する発明を実施するための最良の形態の項及び図面から当業者が容易に理解することができる。
【発明を実施するための最良の形態】
【００３６】
（シミュレーション装置）
図１は、本発明のシミュレーション装置の機能的な構成を示すブロック図である。このシミュレーション装置は、図１に示すように、データ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０を少なくとも具備し、一連の多体分子シミュレーションの機能を有した処理制御部４１と、操作者からのデータや命令などの入力を受け付ける入力部４２と、シミュレーション結果を出力する出力部４５と、原子間相互作用ポテンシャルの定数、および原子核座標と各原子の内部座標などの入力データなどを格納したデータ記憶部４３と、シミュレーションプログラムなどを格納したプログラム記憶部４４とから少なくとも構成されている。
【００３７】
データ解析部１０は、構造評価部１１、統計情報計算部１２、構造可視化部１３等のような、原子座標およびその時系列データの解析機能を少なくとも備えている。例えば、構造評価部１１は、原子間相互作用計算部の原子座標情報管理部の情報から、動径分布関数等を計算する機能を有する。また、原子間相互作用計算部の動的結合次数評価部で得られる各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う機能を有する。統計情報計算部１２は、原子間相互作用計算部の動的結合次数評価部で得られる各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う機能を有する。更に、シミュレーション実行部の古典分子動力学計算部で得られる各原子核座標の時系列情報から、拡散係数などの輸送係数を計算することができる。構造可視化部１３は、原子間相互作用計算部の原子核座標情報管理部に格納されている原子核座標および動的結合次数評価部で得られる各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる機能を有する。結合次数に応じて原子間結合手の太さを変えて表示させることができ、各原子間の結合の強さを直観的に把握できる。
【００３８】
シミュレーション実行部２０は、古典分子動力学計算部２１と、構造最適化計算部２２を少なくとも具備している。古典分子動力学計算部２１は、原子間相互作用計算部３０で計算される原子核座標および各原子の内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する。構造最適化計算部２２は、原子間相互作用計算部３０で計算される原子核座標および各原子の内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する。
【００３９】
原子間相互作用計算部３０は、原子核座標情報管理部３１、内部座標情報管理部３２、動的結合次数評価部３３、ポテンシャルエネルギー計算部３４、１階微分計算部３５、圧力テンソル計算部３６を少なくとも有する。原子核座標情報管理部３１は、原子核座標およびその時間微分の数値を保持する。内部座標情報管理部３２は、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を保持する。動的結合次数評価部３３は、原子核座標と内部座標の値から各原子間の結合次数を評価する機能を有する。ポテンシャルエネルギー計算部３４は、原子核座標と内部座標の値、さらに動的結合次数評価部３３で評価された結合次数から、系の全ポテンシャルエネルギーを計算する。１階微分計算部３５は、原子核座標と内部座標の値、さらに動的結合次数評価部３３で評価された結合次数から、原子核および各原子の内部自由度に加わる力を計算する。圧力テンソル計算部３６は、１階微分計算部３５で計算される原子間力から、系の内部応力テンソルを計算する。
【００４０】
図１において、入力部４２はキーボード、マウス、ライトペン、フレキシブルディスク装置、またはＬＡＮ端子などで構成される。データ記憶部４３およびプログラム記憶部４４は、周知の磁気テープ、磁気ドラム、磁気ディスク、光ディスク、光磁気ディスク、あるいはＲＯＭ、ＲＡＭなどの半導体メモリー等を用いた記録部である。また出力部４５は、ディスプレイ装置やプリンタ装置などにより構成されている。図１に示した処理制御部４１、データ記憶部４３、およびプログラム記憶部４４等はＣＰＵ，およびこのＣＰＵに接続されたＲＯＭ、ＲＡＭ、磁気ディスクなどの記憶装置を含む通常のコンピュータシステムで構成すればよい。図１に示した処理制御部４１で実行される各処理の入力データはデータ記憶部４３に格納され、プログラム命令はプログラム記憶部４４に記憶される。そしてこれらの入力データやプログラム命令は必要に応じてＣＰＵに読み込まれ、演算処理が実行されるとともに、データ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０で発生した数値情報などのデータはデータ記憶部４３に格納される。
【００４１】
次に、原子間相互作用計算部３０の動作について説明する。原子間相互作用計算部３０においては、まず動的結合次数評価部３３により各原子間の結合次数が評価される。その結果を用いてポテンシャルエネルギー計算部３４において系の全ポテンシャルエネルギーが計算され、１階微分計算分３５において各原子核および内部自由度に加わる力が計算される。また、１階微分計算分３５で計算された原子間力から、圧力テンソル計算部３６において系の内部応力テンソルが計算される。
【００４２】
各原子間の結合次数を判定する動的結合次数評価部３３の機能的な構成を示すブロック図を図２に示す。図２に示すように、動的結合次数評価部３３は、結合手の重なり判定部３３１、結合次数判定部３３２、及びスイッチング関数部３３３を少なくとも有している。
【００４３】
図２において結合手の重なり判定部３３１は、原子核座標情報管理部３１及び内部座標情報管理部３２より座標情報を読取り、各原子間の結合手の重なり度を表すスカラ量Ψijを、図２の式（＊１）乃至式（＊３）で計算する。式（＊１）乃至式（＊３）でｉ，ｊは各原子を示す指数（インデックス）で、α，βはそれぞれ、原子ｉ，ｊに属する結合手を示す指数である。また式（＊３）で定義される関数ψの添え字Ｉ，Ｊは、原子i，ｊの種類を示す。ＣＩ，ｐＩ，ｎＩは、Iで指定される原子の種類に固有の定数である。ｓ，θ，φは、原子の内部自由度変数を示し、式（＊１）に見られるそれらの添え字は、各内部自由度変数が属する原子の指数と結合手の指数を示している。ｒijは原子iと原子ｊの間の距離を示し、ｅijは原子iの原子核位置を始点とし原子ｊの原子核位置を終点とするベクトルの単位ベクトルを示す。結合手の重なり度は、式（＊１）に示すように、原子iと原子ｊの間の各結合手の組み合わせについて式（＊２）で示される関数Ｓを計算し、その総和として定義される。式（＊２）において、関数ｆcはスイッチング関数部３３３で評価され、それは式（＊５）で定義される。近距離では１、遠距離では０となる連続関数である。式（＊５）Ｒ，Ｄは、スイッチングを起こす距離を規定するパラメータで、それらは式（＊２）において、２原子の組み合わせＩ，Ｊ毎に指定できる。
【００４４】
図２において結合次数判定部３３２は、結合手の重なり判定部３３１で決定されるスカラ量Ψijを引数としてスイッチング関数部を呼び出し、式（＊４）に従って計算する。式（＊４）におけるＲ，Ｄは、結合次数判定部３３２固有のパラメータである。
【００４５】
以上の動的結合次数評価部３３における結合次数の決定方法を直観的に把握するための模式図を図３に示す。原子ｉと原子ｊ間の結合次数ｇ(Ψij)は、原子ｉと原子ｊ間の結合手の重なり度を表すスカラ量Ψijに応じて０から１まで連続的に増加する。この結合手の重なり度Ψijは、図３のＡ及びＢに示すように、原子ｉと原子ｊ間の距離と、それらの結合手の相対的な向きによって定まる。図３のＡ及びＢに描かれている結合手の形状は、式（＊３）で定義される関数ψを正確に反映している。即ち、式（＊３）で定義された関数ψを、引数である単位ベクトルｅの方位を変え、対応する関数値を単位ベクトルｅの向きをもつ動径の長さで表すと、動径の先端をプロットして描かれる曲面は、Ａ及びＢに示したような偏った亜鈴形を成す。もし、２原子が相互作用できる距離に近づいても、Ａに示したように、互いの結合手に重なりが生じていない場合には、式（＊１）乃至式（＊３）より、Ψijは負の値となり、結合次数ｇ(Ψij)は０と判定される。一方、Ｂに示したように、互いの結合手に重なりが生じている場合には、式（＊１）乃至式（＊３）より、Ψijは正の値を持つと評価され、結合次数ｇ(Ψij)は１と判定される。尚、図３には描かれていないが、動的結合次数項ｇ(Ψij)はΨijの単調増加関数ではなく、Ψijが一定値を超えると再び０に減少する式（＊４）で定義されている。これは、相手の原子を完全に覆い尽くすほど異常に結合手が伸びた場合、相手の原子の更に先に存在する原子とも同時に結合して不自然な安定化を起こす可能性があるため、これを防ぐための処置として、Ψijが一定値を超えると再びｇ(Ψij)が０に減少するようになっている。
【００４６】
系の全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部３４の機能的な構成を示すブロック図を図４に示す。図４に示すように、ポテンシャルエネルギー計算部３４は、１体項計算部３４１、２体項計算部３４２、３体項計算部３４３、及び全エネルギー合算部３４４を少なくとも有している。
【００４７】
図４において１体項計算部３４１は、内部座標情報管理部３２より座標情報を読取り、各結合手の第一自由度変数sの値から、原子ｉの１体ポテンシャルエネルギーｆ１（ｉ）を図２の式（＊６）で計算する。式（＊６）でｉは原子を示す指数で、αは原子ｉに属する結合手を示す指数である。ωｉは、原子ｉに属する結合手の本数を示す。添え字Ｉは、原子ｉの種類を示す。ａＩ，ｍＩは、Iで指定される原子の種類に固有の定数である。
【００４８】
図４において２体項計算部３４２は、原子核座標情報管理部３１より座標情報を読取り、２原子間距離のみで決まる２体ポテンシャルエネルギーｆ２typeを計算する。ｆ２の上付き添え字typeは、２体ポテンシャルエネルギーの種類を指定する標識で、本シミュレーション装置では、等方的に働く反発相互作用を表す添え字repと、結合が組まれた原子間のみに作用する引力相互作用を表す添え字atrの２種類が用意されている。２体ポテンシャルエネルギーｆ２typeは、式（＊７）に示すように、式（＊８）で定義される関数fsw2の２項の和で定義される。式（＊８）で、ｒは２原子間距離を示し、A，ｐ，ａはパラメータを示す。これらのパラメータは式（＊７）に添え字付きで示されているように、２つの原子の種類Ｉ，Ｊと２体ポテンシャルエネルギーの種類によって固有の値を取り得る。特にａはカットオフ距離と呼び、距離の次元を有し、引力相互作用が完全に打ち切られて０になる距離を指定する。
【００４９】
図４において３体項計算部３４３は、原子核座標情報管理部３１より座標情報を読取り、３原子間の位置関係で決まる３体ポテンシャルエネルギーｆ3を計算する。３体ポテンシャルエネルギーｆ3は、式（＊９）に示すように、式（＊１０）で定義される３体関数Λと式（＊１１）で定義される3体関数Θの積で定義される。３体関数Λは、原子i−ｊ間の距離ｒｉｊと原子i−ｋ間の距離ｒｉｋの関数で、λ，γ，ａはパラメータである。各パラメータの下付き添え字は、３体項を構成する原子の種類を指定する。例えば、下付き添え字JIKは、Ｉで指定される原子の種類を頂点とし、それがJおよびＫで指定される原子と結合を組んで結合角を構成している場合を表す。また、γ，ａの上付き添え字は、上記結合角を張る２つの結合の一方を指定するための標識である。ａはカットオフ距離を表し、３体ポテンシャルエネルギーが作用する範囲を指定する。また、3体関数Θは、結合角の角度θjikの関数で、θ0，αはパラメータを示す。各パラメータの下付き添え字は、３体項を構成する原子の種類を指定する。特にパラメータθ0JIKは、Ｉで指定される原子の種類を頂点とし、それがJおよびＫで指定される原子と結合を組むことで構成される結合角の再安定角度を決定する。
【００５０】
図４において、全エネルギー合算部３４４は、１体項計算部３４１、２体項計算部３４２、及び３体項計算部３４３で計算されたｆ１，ｆ２rep，ｆ２atr，ｆ3を式（＊１２）のように合計し、系の全ポテンシャルエネルギーΦを計算する。式（＊１２）において、ｉ，ｊ，ｋは各原子を示す指数で、Ｎは系に含まれる全粒子数を表す。またｇ(Ψ)は、動的結合次数評価部３３で計算される動的結合次数項を示す。式（＊１２）に示されるように、共有結合に起因する引力相互作用ｆ２atrと３体項ｆ3が、動的結合次数(Ψ)が０でない有限の値となる原子間に限定して適用される。このため、企図しない原子間に不要な引力相互作用が発生することがなく、また企図しない原子間に３体項が働いて余計なポテンシャル障壁が発生するのを防ぐことができる。動的結合次数(Ψ)は、各原子の内部座標変数に応じて連続的に変化するため、共有結合の所在を明示しつつその繋ぎ変えを滑らかに起こすことを可能にしている。
【００５１】
各原子核および内部自由度に加わる力を計算する、１階微分計算部３５の機能的な構成を示すブロック図を図５に示す。図５に示すように、１階微分計算部３５は、原子核座標に関する微分実行部３５１、及び内部座標に関する微分実行部３５２を少なくとも有している。
【００５２】
図５において原子核座標に関する微分実行部３５１は、動的結合次数評価部３３及びポテンシャルエネルギー計算部３４で得られる計算結果を利用し、式（＊１３）より全ポテンシャルエネルギーの各原子核位置に関する１階偏微分係数を計算する。式（＊１３）において、ｉ，ｊ，ｋは、各原子を示す指数である。
【００５３】
図５において内部座標に関する微分実行部３５２は、動的結合次数評価部３３及びポテンシャルエネルギー計算部３４で得られる計算結果を利用し、式（＊１４）より全ポテンシャルエネルギーの各内部座標変数に関する１階偏微分係数を計算する。式（＊１４）において、ξiαは、式（＊１）で示した、原子ｉに属する結合手αの、３つの内部自由度変数siα，θiα，φiαをまとめたベクトル変数である。siα，θiα，φiαはこの順番で、ξiαの第一要素、第二要素、第三要素に対応する。
【００５４】
１階微分計算部３５で原子核座標及び内部座標に関するポテンシャルの１階微分が計算されると、それらを質量で除算することで、直ちに原子核及び内部座標変数に加わる力を求めることができる。ただし、ここで言う質量とは、原子核に加わる力を求める際はkgの次元を有する原子核の質量そのものであるが、内部座標変数に加わる力を求める際は、内部座標変数用に別途指定される、エネルギーの次元を有する仮想的な質量を用いる。質量は一般に、運動する自由度の変化の容易さを決めるパラメータであり、その次元はkgに限らない。本シミュレーション装置において内部座標変数を原子核に素早く追従させたい場合には、原子核の質量に比べて十分小さな質量を指定する。ただし、内部座標変数に対してあまり小さな質量を指定すると、計算の精度低下を招くため、扱う系に応じて適切な質量の値を使用者が指定することになる。
【００５５】
系の内部応力テンソルを計算する、圧力テンソル計算部３６の機能的な構成を示すブロック図を図６に示す。図６に示すように、圧力テンソル計算部３６は、各２原子間に働く力を計算する２原子間力計算部３６１、及び圧力テンソル合算部３６２を少なくとも有している。
【００５６】
図６において２原子間力計算部３６１は、１階微分計算部３５で得られる計算結果を利用し、式（＊１５）より、原子ｉが原子ｊより受ける力を、全ての原子ｉ，ｊの組み合わせについて計算する。
【００５７】
図６において圧力テンソル合算部３６２は、２原子間力計算部３６１で得られる計算結果を利用し、式（＊１６）より、系の全内部応力テンソル（圧力テンソル）を計算する。式（＊１６）において|vi>は、原子iの速度ベクトルを示す。各原子核の速度ベクトルは、シミュレーション実行部２０の古典分子動力学計算部２１において逐次生成され、原子核座標情報管理部３１に格納される。シミュレーション実行部２０の構造最適化計算部２２を実行する場合は、この速度ベクトルは生成されず、原子核座標情報管理部３１において速度ベクトル情報を格納するメモリには常に数値０がセットされる。従って、構造最適化計算部２２を実行する場合の系の全内部応力テンソルは、式（＊１６）の中括弧内のΣΣ|fij><rij|の項のみが有効となる。
【００５８】
図６において圧力テンソル合算部３６２では、内部自由度変数の指数は表れていない。これは、内部自由度変数が各原子核位置を中心とする極座標系を運動しており、圧力テンソルの起源である運動量輸送は、原子間相互作用と原子核の移動によって行われると考えられるためである。従って、圧力テンソルは従来法通り、原子核の運動エネルギーと原子間相互作用のみで定義し、内部自由度変数の運動エネルギー分は含まれない。
【００５９】
以上のように、原子間相互作用計算部３０によれば、原子核同様、各原子の共有結合手の方向を決定する内部自由度変数にも力が働くため、シミュレーション実行部２０によって、その力に基づいて結合手の方向を変化させることができる。これにより、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができ、共有結合の組み換えを伴う分子シミュレーションを行うことが可能となる。この拡張に伴って増える計算は、内部座標情報管理部３２のメモリ操作及び動的結合次数評価部３３に止まるため、従来の古典分子動力学計算と同様、数千個以上の原子から構成される大規模系のシミュレーションが、短時間で効率的に実行できる。
【００６０】
また、原子間相互作用計算部３０によれば、２体引力相互作用および３体相互作用の発現を結合が組まれている原子間に限定することができる。このため、企図しない原子間に余計な引力やポテンシャル障壁が形成される事を防ぐことができる。更に、結合を組んでいない原子間の２体引力相互作用を打ち消すように３体相互作用項を調節する必要がなくなるため、２体ポテンシャル定数と３体ポテンシャル定数を独立に決定することができ、多元素混在系用のポテンシャル定数を用意することが可能となる。
【００６１】
(シミュレーション方法)
図１に示した原子間相互作用計算部３０を用いて分子動力学シミュレーションを実行するには、所望の統計アンサンブルを生成するように設計されたオイラー−ラグランジュ方程式の数値積分を行えばよい。また、構造最適化計算を実行するには、原子間相互作用計算部３０で計算される力の情報を用いて、共役勾配法等、多次元空間上のポテンシャルの極小点探索アルゴリズムを実行すればよい。以下において、分子動力学計算及び構造最適化計算を行うための具体例、即ち、シミュレーション方法１及びシミュレーション方法２を示す。
【００６２】
シミュレーション方法１：図７のフローチャートを参照して、本発明のシミュレーション方法１を説明する。
（イ）まず、ステップＳ１１０で、分子動力学計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標を入力部４２から読み込む。ここで分子動力学計算を実行する上で必要となる制御データとは、数値積分を実行する際の1ステップ当たりの時間刻み、全ステップ数など、分子動力学計算を実行する上で必要となるデータである。また温度制御や圧力制御アルゴリズムを使用するかどうかの選択や、その際の設定温度及び設定圧力の数値なども含む。分子動力学計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標の情報は基本的に全て入力部４２より読み込むが、入力が省略された情報については、必要に応じて自動的に補完する。例えば、内部座標情報が省略された場合、原子核座標情報に基づいて、各原子の最も近い原子から順に結合を組むように、内部自由度の座標を自動的に決定する。
（ロ）次にステップＳ１２０で、原子核座標及び内部座標の情報から、各原子間の動的結合次数を評価し、原子間結合のトポロジーを定める。
（ハ）次にステップＳ１３０で、上記原子間結合のトポロジーに含まれる２体ポテンシャルおよび３体ポテンシャルの定数が完備されているか判定する。完備されていない時はステップＳ１４０に進み、不足しているポテンシャル定数を要求するエラーメッセージを出力し、計算の途中データを記録して終了する。必要なポテンシャル定数が揃っている時はステップＳ１５０に進む。
（ニ）ステップＳ１５０ではポテンシャルエネルギー、１階微分、圧力テンソルを計算する。
（ホ）ステップＳ１６０では原子座標、内部座標に関する１階微分の値から力を計算し、これをオイラー−ラグランジュ方程式に代入して各座標の加速度を計算し、分子動力学計算を１ステップ進める。
（ヘ）ステップＳ１７０では入力データで指定された全ステップ数に達したか判定する。達していない場合はＳ１２０に戻り、達している場合には計算結果を出力して終了する。
【００６３】
シミュレーション方法２：図８のフローチャートを参照して、本発明のシミュレーション方法２を説明する。
（イ）まず、ステップＳ２１０で、構造最適化計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標を入力部４２から読み込む。ここで構造最適化計算を実行する上で必要となる制御データとは、構造最適化計算の収束判定条件など、構造最適化計算を実行する上で必要となるデータである。また、周期的境界条件を設けて無限に大きな系の部分系を模す場合、部分系を構成する単位セルのサイズ及び形状も最適化できるが、この単位セルのサイズ及び形状の変化を許すか否かを示す制御情報も含む。構造最適化計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標の情報は基本的に全て入力部４２より読み込むが、入力が省略された情報については、必要に応じて自動的に補完する。例えば、内部座標情報が省略された場合、原子核座標情報に基づいて、各原子の最も近い原子から順に結合を組むように、内部自由度の座標を自動的に決定する。
（ロ）次にステップＳ２２０で、原子核座標及び内部座標の情報から、各原子間の動的結合次数を評価し、原子間結合のトポロジーを定める。
（ハ）次にステップＳ２３０で、上記原子間結合のトポロジーに含まれる２体ポテンシャルおよび３体ポテンシャルの定数が完備されているか判定する。完備されていない時はステップＳ２４０に進み、不足しているポテンシャル定数を要求するエラーメッセージを出力し、計算の途中データを記録して終了する。必要なポテンシャル定数が揃っている時はステップＳ２５０に進む。
（ニ）ステップＳ２５０ではポテンシャルエネルギー、１階微分、圧力テンソルを計算する。
（ホ）ステップＳ２６０では、次の１次元極小点探索の方向を決定する。ここで言う方向とは、原子核座標および各原子の内部自由度座標を併せた全自由度から構成される多次元空間内における方向ベクトルを指す。最初にこのステップに来た場合には１階微分から得られる力の方向を指定し、２回目以降にこのステップに来た場合は、既に探索した方向に共役な方向を、共役勾配法の手続きで計算する。
（ヘ）ステップＳ２７０では、前のステップで指定された、原子核座標および各原子の内部自由度座標を併せた全自由度から構成される多次元空間内における方向ベクトル上で、１次元の極小点探索を実行する。
（ト）ステップＳ２８０では、構造最適化計算の収束判定条件が満たされるか判定する。収束判定条件とは、具体的には、ステップＳ２７０の前後で、原子核座標および内部自由度座標の移動幅を計算し、その合計が入力データで指定した範囲に収まっていることを指す。収束判定条件が満たされない時は、ステップＳ２２０に戻る。収束判定条件が満たされた時は、計算結果を出力して終了する。
【００６４】
以上のシミュレーション方法１及びシミュレーション方法２によれば、共有結合の所在を明示しながら、その滑らかな繋ぎ変えを可能にする分子動力学計算及び構造最適化計算を行うことができる。共有結合の発現を明示された原子間に限定できるため、従来法で問題となっていた余計な相互作用が取り除かれ、多元素が混在する複雑な系のポテンシャル定数を用意することが可能となる。つまり、シミュレーション方法１及びシミュレーション方法２によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子シミュレーションを、数千原子以上の大きな規模で実行できる。
【００６５】
（シミュレーションプログラムを格納した記録媒体）
図７及び図８に示すシミュレーション方法を含め、図１に示したデータ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０におけるシミュレーション方法を実現するためのプログラムはコンピュータ読取り可能な記録媒体に保存しておいても良い。この記録媒体をコンピュータシステムによって読み込ませ、図１に示したプログラム記憶部４４に格納し、このプログラムを処理制御部４１で実行して本シミュレーション方法を実現することもできる。ここで、記録媒体とは、例えば磁気ディスク、光ディスク、光磁気ディスク、磁気テープなどのプログラムを記録することができるような媒体などが含まれる。コンピュータの外部メモリも、ここで言う記録媒体に含まれる。
【００６６】
（実施例）
本発明の実施例として、結合組み換えを伴う３元素混在系の分子動力学計算、原子核系と内部自由度の温度制御分子動力学計算、Ｓｉ結晶モデルの圧力制御分子動力学計算、シリコン酸化膜とシリコン界面のモデリングを示す。
【００６７】
図９はＳｉＯ４Ｈ４分子の分子動力学シミュレーションを実施した際の初期構造（0ステップ）、2000ステップ、4000ステップ、5000ステップ、6000ステップ、8000ステップ経過後の構造を示す。各原子を球で表示し、各原子の結合手に相当する内部自由度を円柱で表示している。結合手１本あたり３次元の自由度を有し、円柱の配向が２自由度、長さが１自由度に対応する。結合手の本数はＳｉ原子は４、Ｏ原子は２、Ｈ原子は１とした。この表示では、例えば0ステップ目の状態で結合手は重なっていないように見えるが、計算上は向かいあった円柱の間に結合が組まれている。適当な初速度を与えて分子動力学計算を開始すると、分子動力学計算のステップが進むにつれて、原子核及び結合手の向き、長さが変化していく。向かい合った円柱は、結合が切れないよう、向かい合った状態を維持するように変化していく。これは、結合が組まれているときのみ、２体引力相互作用による安定化が起こるためである。5000ステップ目ではＡで示したＯ原子の結合手がそれまで結合していたＨ原子と反対側に結合手を伸ばし、一時的にＯ−Ｈ結合の切断が起こったが、間もなくＨ原子と再び結合しなおす様子が再現された。初期温度をより高く設定すると、原子間結合が完全に切断され、原子核が遠くに離れていく様子も再現される。この小規模の分子動力学計算から、本発明により、２元素以上の多元素混在系において、結合の組み換えを伴う分子動力学計算が可能であることが示される。
【００６８】
図１０は、図９で示した分子動力学計算中の、ハミルトニアンの変動を示す。古典力学では、運動エネルギーとポテンシャルエネルギーの総和であるハミルトニアンは保存量であり、このハミルトニアンの保存を確かめることで、力の計算、即ち、ポテンシャル関数の１階微分が正しいか否かを判定できる。図１０では、時間刻み0.00008τ、0.00005τ、0.00001τで行った分子動力学計算におけるハミルトニアンの変動を示している。時間刻みが小さくなるほど、離散化誤差の影響が減少し、ハミルトニアンの変動が抑えられ、一定の値に保存していくことがわかる。
【００６９】
図１１は、温度制御分子動力学計算を実施した際の、原子核系の温度、内部自由度系の温度、及び原子核と内部自由度系を合わせた全系の温度の時間変化を表す。ここで温度とは、粒子の平均運動エネルギーから次の式(１)で定義される量である。
【００７０】
Ｔ ＝ （2／3kB）<1/2 m v2> (１)
【００７１】
ここでTは温度、kBはボルツマン定数、mは原子の質量、vは原子の速度を表す。内部自由度系の温度も同様に、次の式(２)で定義する。
【００７２】
Ｔ ＝ （2／3kB）<1/2 μ (dξ/dt)2> (２)
【００７３】
ここでμは内部自由度変数の仮想的な質量、ξは内部自由度変数ベクトルを表す。図１１(a)は、原子核系と内部自由度系を合わせた全体の温度が、能勢の温度制御アルゴリズムにより計算開始時から1000Kに保持されている。しかし、原子核系と内部自由度系のそれぞれの温度は、初期状態で原子核系は約2000K、内部自由度系は500K程度と大きな差があり、両者が熱平衡状態に達するまでに長い分子動力学計算を必要とする。これは、原子核座標系と内部自由度座標系の相互作用が、動的結合次数項を介してしか起こらないため、結合組み換え時の極限られた時にしか運動エネルギーを交換できないからである。
【００７４】
原子核座標系と内部自由度座標系の温度を所望の温度に制御するには、同時に複数の熱浴を設けて温度制御が行える、能勢-Hoover形式を使えばよい。図１１(b)は、能勢-Hoover形式を用いて、原子核座標系と内部自由度座標系の温度を1000Kに制御した例を示す。初期状態では、図１１aと同様、原子核座標系と内部自由度座標系の温度に大きな差があるが、この差は直ちに解消され、計算開始直後から平均温度が1000Kに制御されていることがわかる。この能勢-Hoover形式を用いれば、原子核座標系と内部自由度座標系を別々の温度で制御することも可能である。
【００７５】
圧力テンソルの起源である運動量輸送は、原子間相互作用と原子核の移動によって行われると考えられるため、本発明においても、圧力テンソルは従来法通り、原子核の運動エネルギーと原子間相互作用のみで定義し、内部自由度変数の運動エネルギー分は含めていない。この方法でも、正しく圧力制御の分子動力学計算を実施できる。
【００７６】
図１２は、圧力制御の分子動力学計算で用いたＳｉ結晶５７６原子モデルである。結晶モデルを囲む立方体は、本計算で適用した３次元周期的境界条件の単位セルを表す。パリネロ−ラーマン(Parrinello-Rahman)の圧力制御分子動力学法により、この単位法セルの大きさを調節して、内部の圧力を設定圧力5Gpaに制御した結果を図１３に示す。Aは、各瞬間で計算した圧力テンソルの対角成分の平均値をプロットしたものである。変動の幅が大きく、設定圧力に保たれている様子を確認することが難しいが、Ｂのように計算開始時からの平均値をプロットすると、設定圧力5Gpaに正確に保持されていることが確認できる。
【００７７】
従来の結合組み換え型多体ポテンシャルでは、余計なポテンシャル障壁により圧縮ストレスを過大評価してしまう問題が生じたが、共有結合による相互作用を結合手が組まれた原子間に限定できる本発明では、この問題が解決される。この効果を確認するため、圧縮ストレスが顕著に観察される例として、シリコン基板上のシリコン酸化膜構造を選び、従来の結合組み換え型多体ポテンシャル、及び本発明の動的結合手を含む原子間相互作用ポテンシャルで上記シリコン酸化膜構造モデルを最適化し、最適化後の構造、及び残留ストレスの値を比較した。
【００７８】
Ｓｉ酸化膜モデルは、まずＳｉ単結晶モデルを作成し、そこにＯ原子を徐々に導入して作成する。初期のＳｉ単結晶構造は4800個のSi原子から構成される、表面が(001)面を向いた単結晶モデルで、幅及び奥行が3.84nm、厚さが6.51nmに相当する。幅及び奥行方向に２次元の周期的境界条件を設けている。このＳｉ単結晶モデルの表面から、Ｏ原子を１個ずつ、Si-Si結合の中点に挿入し、これを繰り返してＳｉ酸化膜部分を形成する。Ｏ原子を挿入するサイトは面内で無作為に選ぶが、１原子層の酸化が完了するまで次の層へ導入はしない。この酸化方法は、Ｓｉ結晶表面の酸化が１原子層ずつ進行することを示す実験事実に基づいている。Ｏ原子を１個導入する度に、構造最適化計算を行う。また、Ｏ原子を２０個導入する度に、温度制御分子動力学計算を20000ステップ実行し、構造のアニール（焼き鈍し）を行う。これは、高温のアニールで結合の組み換えを促進し、酸化膜構造のアモルファス化及び欠陥の消滅を起こすためである。設定温度は800℃（1073K）とする。時間刻みは2.30×10-5ピコ秒とした。分子動力学計算の数値積分アルゴリズムには、5値のGearの予測子・修正子法を用いた。
【００７９】
図１４(a)に、Si,O混在系用に開発された従来の結合組み換え型多体ポテンシャルを用いて、上記の方法で作成したＳｉ酸化膜とＳｉ基板の界面モデルを示す。酸化膜は2800個のＯ原子を含み、膜厚は約3.0nmに相当する。この酸化膜モデルのストレスは、約13GPaと見積もられた。酸化膜部のストレスは、積層構造全体の内部応力テンソルの対角成分の、幅及び奥行方向の成分の平均値とする。内部応力テンソルを計算する際の体積は、酸化膜部分の体積 (3.84 x 3.84 x 3.0 nm3)を用いた。酸化膜の面内ストレスは、実験では0.4GPaから0.7GPaと報告されているが、これら実験値に比べて著しく大きい。これは、O-Si-O３体項等が余計なポテンシャル障壁を形成し、構造の圧縮を強く阻んでいるためと考えられる。
【００８０】
図１４(b)は、本発明におけるシミュレーション方法で、本発明におけるシミュレーション装置を用いて、図１４(a)で得られたＳｉ酸化膜とＳｉ基板の界面構造をアニール及び構造最適化した結果得られた構造を示す。アニールは、温度制御分子動力学計算を、上記の条件で100000ステップ実行することで行われた。図１４(a)及び図１４(b)の構造は、一見大きな違いが見られないが、余計なポテンシャル障壁の消滅により、本発明で作成した図１４(b)の構造の方が、極僅か酸化膜厚が減少している。図１４(b)の酸化膜部のストレスは、約2GPaであり、従来の結合組み換え型多体ポテンシャルに比べて大幅に減少し、実験値の数倍程度の大きさまで落ち着いた。依然実験値よりも大きい理由として、上記分子動力学計算ならびに構造最適化計算において、単位胞サイズを完全に固定している点が挙げられる。実際のＳｉ酸化膜では、酸化膜のストレスにより、下地のSi基板が横方向に反り、若干ではあるが横方向に膨張している。計算では横方向の膨張を完全に抑えていることを鑑みれば、2GPaという酸化膜部のストレスは、十分妥当な値であり、少なくとも従来の結合組み換え型多体ポテンシャルで得られる13GPaと比較すれば、著しい改善である。
【００８１】
図１５は、図１４(a)及び図１４(b)の構造を用いて計算した、２体相関関数を示す。２体相関関数は、
【００８２】
ｇ（ｒ） ＝ < n(r,r+Δr)> / 4πr2ρΔr (3)
で計算される、距離の関数である。< n(r,r+Δr)>は、ある原子を中心とした半径rと半径r+Δrの球殻の間に存在する原子の数を、系の全ての原子を中心に計算してその平均をとったもので、ρは系全体の原子数密度を表す。２体分布関数は位置の相関が無くなると１に近づく。
【００８３】
1.62Å（オングストローム）付近に立つピークは酸化膜中のSi-O結合を表し、2.35Å付近に立つピークは下地のSi基板中のSi-Si結合を表す。2.6Å付近に立つ緩やかなピークは、SiO4正四面体構造のO-O間距離に起因する。3.1Å付近に立つ小さなピークは、酸化膜中の第２近接のSi-Si対を表す。プロットAでは、3.3Åに鋭いピークが見られる。これは解析から、Si-O対に起因するピークであることがわかっている。この3.3Å付近の鋭いピークは、アモルファスSiO2のX線回折から得られる２体相関関数では観察されない不自然なピークである。この不自然なピークは、図14(b)に示した、本発明で作成した構造から計算した２体分布関数を示すプロットBには見られない。その分、3.1Åから3.3Åの成分が増えている。この結果は、従来型ポテンシャルでは余計なポテンシャル障壁によって接近が阻まれていたSi-O対が、本発明におけるシミュレーション方法では、正しく接近できるように改善されたことを示している。
【００８４】
本発明によれば、共有結合のような方位依存性のある引力相互作用の発現箇所を明示しつつ、その滑らかな繋ぎ変えを許す分子動力学シミュレーション及び構造最適化計算が可能となる。
【００８５】
つまり、本発明によれば、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができ、共有結合の組み換えを伴う分子シミュレーションを行うことが可能となる。
【００８６】
また、本発明によれば、従来の古典分子動力学計算と同様、数千個以上の原子から構成される大規模系のシミュレーションが、短時間で効率的に実行できる。
【００８７】
さらに、本発明によれば、企図しない原子間に余計な引力やポテンシャル障壁の形成を防ぐことができる。
【００８８】
さらに、本発明によれば、多体相互作用項同士の干渉を大幅に減少させることができる。これにより、多体相互作用項のポテンシャル定数の決定が大幅に単純化され、多元素が混在する複雑な系のポテンシャル定数を用意することが可能となる。
【００８９】
つまり、本発明によれば、従来法では不可能であった、多元素混在系において、共有結合の組み換えを伴う化学変化の大規模分子シミュレーションが可能となる。
【００９０】
さらに、本発明によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子動力学シミュレーション及び構造最適化計算が可能なシミュレーション方法を提供することができる。
【００９１】
さらに、本発明によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子動力学シミュレーション及び構造最適化計算が可能なシミュレーションプログラムを格納した記録媒体を提供することができる。
【図面の簡単な説明】
【００９２】
【図１】本発明の実施の形態に係るシミュレーション装置の機能的な構成を示すブロック図である。
【図２】図１に示した動的結合次数部３３の機能的な構成を示すブロック図である。
【図３】動的結合次数評価部３３における結合次数の決定方法を模式的に示した図である。
【図４】図１に示したポテンシャルエネルギー計算部３４の機能的な構成を示すブロック図である。
【図５】図１に示した１階微分計算部３５の機能的な構成を示すブロック図である。
【図６】図１に示した圧力テンソル計算部３６の機能的な構成を示すブロック図である。
【図７】本発明のシミュレーション方法１を示すフローチャートである。
【図８】本発明のシミュレーション方法２を示すフローチャートである。
【図９】図９は、３元素が混在するＳｉＯ４Ｈ４分子の分子動力学シミュレーションを実施した際の構造変化を示す。
【図１０】図１０は、図９で示した分子動力学計算中の、ハミルトニアンの変動を示す。
【図１１】図１１は、温度制御分子動力学計算を実施した際の、原子核系の温度、内部自由度系の温度、及び原子核と内部自由度系を合わせた全系の温度の時間変化を表す。
【図１２】Ｓｉ結晶５７６原子モデルを示す図である。
【図１３】図１２のＳｉ結晶モデルの単位法セルの大きさを調節して、内部の圧力を実施した際の圧力の変動を示した図である。
【図１４】図１４(a)は、従来の結合組み換え型多体ポテンシャルを用いて作成したＳｉ酸化膜とＳｉ基板の界面モデルを示す。図１４(b)は、本発明におけるシミュレーション方法で、図１４(a)の界面構造をアニール及び構造最適化した結果得られた構造を示す。
【図１５】図１４(a)及び図１４(b)の構造を用いて計算した、２体相関関数を示した図である。
【符号の説明】
【００９３】
１０ データ解析部
１１ 構造評価部
１２ 統計情報計算部
１３ 構造可視化部
２０ シミュレーション実行部
２１ 古典分子動力学計算部
２２ 構造最適化計算部
３０ 原子間相互作用計算部
３１ 原子核座標情報管理部
３２ 内部座標情報管理部
３３ 動的結合次数評価部
３４ ポテンシャルエネルギー計算部
３５ １階微分計算部
３６ 圧力テンソル計算部
４１ 処理制御部
４２ 入力部
４３ データ記憶部
４４ プログラム記憶部
４５ 出力部
３３１ 結合手の重なり判定部
３３２ 結合次数判定部
３３３ スイッチング関数部
３４１ １体項計算部
３４２ ２体項計算部
３４３ ３体項計算部
３４４ 全エネルギー合算部
３５１ 原子核座標に関する微分実行部
３５２ 内部座標に関する微分実行部
３６１ ２原子間力計算部
３６２ 圧力テンソル計算部
【技術分野】
【０００１】
本発明は半導体集積回路製作等に用いる材料設計支援シミュレーション技術、加工プロセスのシミュレーション技術、あるいは、蛋白質等の生体分子やその他有機化合物質の分子シミュレーション技術に係り、特にナノメートルスケールのサイズを有する数千原子規模以上の系を対象とする大規模分子シミュレーション方法、シミュレーション装置、シミュレーションプログラムを格納した記録媒体に関する。
【背景技術】
【０００２】
半導体集積回路の高集積化に伴い、ナノメートルオーダーの極めて高い精度の微細加工技術が要求されている。例えば、金属−絶縁膜−半導体型電界効果トランジスタの絶縁膜部分は僅か1ナノメートルまで薄層化されている。これはシリカの構成単位であるSiO4正四面体構造の４層程度の厚さであり、積層構造を１原子層レベルで再現性良く制御する技術が要求される。ここまで微細化が進行すると、連続的な拡散方程式を解く従来型のプロセス・シミュレーションでは、界面形状や組成分布の予測能力に限界があり、原子１個１個の挙動を視野に入れた分子シミュレーション技術の導入が必要となっている。
【０００３】
一方創薬分野においては、分子シミュレーションは既に開発現場で欠かせない技術となっている。薬剤候補となる膨大な数の低分子化合物と蛋白質の親和性を計算機で予測することで、候補薬剤の絞込み、化合物最適化の高効率化に重要な役割を果たしている。
【０００４】
一般に全ての物質の性質は、電子および原子核間のクーロン相互作用を含む相対論的なディラック方程式、もしくは非相対論的なシュレディンガー方程式をパウリの排他原理の下で解くことで、基本的な物理定数のみから非経験的に予測することができる。原子核を固定した場合の近似的な数値解法としてハートリー・フォック（Hartree-Fock）方程式を解く分子軌道法や、コーン・シャム（Kohn-Sham）方程式を解く密度汎関数法が確立している。また、密度汎関数法で得られる全エネルギーを断熱ポテンシャルとして動力学計算するカー・パリネロ（Car-Parrinello）法も広く普及している。しかし、これら第一原理計算手法が取り扱える系の規模は、その計算量の多さ故、スーパーコンピュータを使用しても数百原子程度に限定される。１辺にたった１０個の原子を並べた立方晶ですら１０００原子に達することを鑑みても分かるように、第一原理計算がナノスケールの構造体の加工プロセスのシミュレータの役割を果たすことは期待できない。
【０００５】
数千原子以上の原子団の挙動を計算するには、原子間相互作用を単純なポテンシャル関数で近似し、その原子間ポテンシャルを用いてニュートンの運動方程式を数値的に解く古典的分子動力学法が有効である。電子波動関数の固有値問題を解かないこの方法では、扱う元素の組み合わせに応じて原子間ポテンシャル関数を設計する必要があるが、一度ポテンシャル関数が設定された系に対しては、コンピュータシステムにより、数十原子レベルのミクロな系から数十万原子レベルのナノスケール領域に至る広い原子数規模の範囲の動的シミュレーションを、短時間で、効率的に実行できる。
【０００６】
古典的分子動力学計算で用いられる原子間ポテンシャル関数は、大きく２種類に分類される。１つは、原子座標と共に共有結合のトポロジーを入力情報とするポテンシャルで、MM2やAMBERなどが知られている。パラメータの決定が比較的容易で多元素混在系が扱えることから、蛋白質分子の立体構造を基にした計算機創薬で中心的役割を果たしている。ただし、共有結合のトポロジーは固定されるため、例えば加水分解や脱水縮合によるペプチド結合の組み換えような化学変化は再現できない。
【０００７】
これに対し、原子座標のみを入力情報とするポテンシャルも数多く存在する。希ガスやイオン結合性の無機化合物を表現する単純な２体相互作用項からなるポテンシャルを始め、金属や有機化合物、半導体など、強い結合を含む物質を表現できる複雑な多体項を有するポテンシャルまで、様々な関数が提案されている。これらは結合箇所を露に指定しないため、結合組み換えを伴う様々な構造変化も扱える。半導体装置のプロセス・シミュレーションにおいては、結晶化やアモルファス化、表面や界面での再構成、欠陥の挙動など、共有結合の組み換えを伴う構造変化の再現が重要となるため、導入されつつある分子シミュレーション技術では主に複雑な多体相互作用ポテンシャル関数が用いられる。しかし、多元組成系に適用できるポテンシャル定数を揃えることが極めて難しく、複雑な多体相互作用ポテンシャルの場合、同時に扱える元素は通常２種類程度に限定される。
【０００８】
【非特許文献１】R. Car and M. Parrinello, Physical Review Letters Vol. 55, (1985) p.2471.
【非特許文献２】N. L. Allinger, Journal of American Chemical Society Vol. 99, (1977) p.8127.
【非特許文献３】S. J. Weiner, P. A. Kollman, D. T. Nguyen and D. A. Case, Journal of Computational Chemistry Vol. 4, (1983) p.187.
【非特許文献４】F. H. Stillinger and T. A. Weber, Physical Review B Vol.31, (1985), p.5262.
【非特許文献５】J. Tersoff, Physical Review B Vol.39, (1989) p.5566.
【非特許文献６】T. Watanabe, H. Fujiwara, H. Noguchi, T. Hoshino, and I. Ohdomari, Japanese Journal of Applied Physics Vol. 38, (1999) p.L366.
【非特許文献７】J. Dalla Torre, J.-L. Bocquet, Y. Limoge, J.-P. Crocombette, E. Adam, G. Martin, T. Baron, P. Rivallin and P. Mur, Journal of Applied Physics Vol. 92, (2002) p.1084.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
半導体集積回路のゲート絶縁膜はHfSiO2など、複雑な組成の高誘電率材料に置き換えられようとしており、半導体プロセスで用いられる材料の多様化が急速に進んでいる。このため、従来の結合組み換え型ポテンシャルでも、より多くの元素を同時に扱える必要が生じている。
【００１０】
結合組み替え型の多体相互作用ポテンシャルには、ポテンシャル定数の決定が困難であること以外に、もう一つ深刻な欠点を抱えている。多体相互作用ポテンシャルでは、３体項以上の多体項が配位数を調節する役割を果たしているが、これが企図しないポテンシャル障壁を形成する。このため、圧縮ストレスを受けた材料のシミュレーションでは材料の内部応力が１〜２桁も過大評価される場合がある。また、一般に３体項を含むポテンシャルを用いた融解シミュレーションでは、融点が実験値より高くなる傾向があるが、これも余計なポテンシャル障壁が融解を阻止していることが原因と考えられる。
【００１１】
ターゾフ（Tersoff）ポテンシャルなど、２体相互作用項に周辺の原子配置に応じて変化する結合次数を導入するボンドオーダー型ポテンシャルでは、３体項以上の多体項を導入することなく配位数調節が実現されているが、第一原理計算では見られない準安定構造が再現されるなど、複雑な結合次数項が原因と見られる副作用が、安定構造から大きく外れた領域で顕在化する。また、関数の複雑さ故、ポテンシャル定数と力定数の関係が明瞭でなく、多元素混在系への拡張を一層困難にしている。
【００１２】
多元素混在系に適用でき、かつ結合の組み換えが許されるようなポテンシャル関数を実現するには、原子核配座のみで一意にエネルギーを決定する従来法では限界がある。共有結合は、原子間で電子軌道の重なりが生じている部位に発現するものであり、原子間距離で一意に決まるものではない。従って、結合箇所を指定するため、電子の空間分布に相当する付加的な情報も必要と考えられる。
【００１３】
上記問題点に鑑み、本発明は、多元素混在系においても結合組み換えを伴う構造変化を正確に再現する、大規模分子動力学シミュレーション装置を提供することを目的とする。
【００１４】
具体的には、本発明は、従来の多体ポテンシャルが有していた余計なポテンシャル障壁を含むことなく、圧縮ストレスを帯びた物質の内部応力も正しく表現する、結合組み換え型の原子間相互作用モデルとその運動方程式を提供し、数千原子以上の多元素混在系においても、結合組み換えを伴う化学変化を再現するシミュレーション装置を提供することである。
【００１５】
本発明の更に他の目的は、数千原子以上の多元素混在系においても、結合組み換えを伴う化学変化を高速で計算するシミュレーション方法を提供することである。
【００１６】
本発明の更に他の目的は、数千原子以上の多元素混在系においても、結合組み換えを伴う化学変化を高速で計算するシミュレーションプログラムを格納した記録媒体を提供することである。
【課題を解決するための手段】
【００１７】
上記目標を達成するため、本発明の第１の主要な観点によれば、原子核座標及び内部座標の情報を格納する座標情報格納部と、前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価部と、前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部と、前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算部と、前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション部と、前記シミュレーション結果を出力する出力部とを有するシミュレーション装置が提供される。
【００１８】
また、この発明の一の実施形態によれば、前記シミュレーション部は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理部と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理部と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価部と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものである。また、この場合、前記シミュレーション実行部は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算部とを有するものであることが好ましい。
【００１９】
さらに、この発明の一の実施形態によれば、この装置は、前記出力部で出力されたデータの解析を行うデータ解析部をさらに有する。この場合、データ解析部は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価部を有するものであることが好ましい。また、このデータ解析部は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う構造解析部を有するものであることが好ましい。さらに、前記データ解析部は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算部を有するものであることが好ましい。さらに、前記データ解析部は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算部を有するものであることが好ましい。また、前記データ解析部は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析部を有するものであることが好ましい。また、前記データ解析部は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化部を有するものであることが好ましい。この場合、前記構造可視化部は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものである。
【００２０】
この発明の別の一の実施形態によれば、前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものである。この場合、前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものである。
【００２１】
さらに、この発明の別の一の実施形態によれば、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものである。
【００２２】
本発明によるシミュレーション装置によれば、各原子間の結合次数を、原子間距離と各原子の内部自由度変数の値から簡単に決定できるため、共有結合エネルギーを表す２体引力相互作用を、結合が組まれている原子間のみに発現させることができる。結合を組んでいない原子間には２体引力相互作用をはじめ、３体以上の多体相互作用項が働かないため、企図しないポテンシャル障壁の形成を防ぐことができる。
【００２３】
また、本発明によるシミュレーション装置によれば、結合を組んでいない原子間に共有結合を表す２体引力相互作用が適用されないため、余計な２体引力項を打ち消すように３体相互作用項を調節して配位数の調節を行う必要がなく、共有結合手の方向を決定する内部自由度の数そのもので、適切な配位数を保つことができる。このため、２体ポテンシャル定数と３体ポテンシャル定数を独立に決定することができ、多元素混在系のポテンシャル定数を用意することが可能となる。
【００２４】
また、本発明によるシミュレーション装置によれば、結合が途切れた３原子間に３体相互作用項を働かせないことができる。このため、結合が途切れた３原子間に不要に働く３体項を打ち消すようにポテンシャル定数を調節する必要がなく、異なる３体相互作用項のポテンシャル定数を独立に決定することができ、多元素混在系のポテンシャル定数を用意することが可能となる。
【００２５】
また、本発明によるシミュレーション装置によれば、原子核同様、各原子の共有結合手の方向を決定する内部自由度変数にも力が働き、その力に基づいて運動することができるため、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができる。これにより、シミュレーション中に共有結合の組み換えを連続的に起こすことができ、共有結合の組み換えを伴う動的過程を再現することができる。
【００２６】
本発明によるシミュレーション方法によれば、共有結合の所在を明確に示しながら、その組み換えを連続的に起こすことができ、共有結合の組み換えを伴う化学変化を再現することができる。
【００２７】
本発明によるシミュレーション方法によれば、内部座標情報管理部を実現するための記憶容量と、動的結合評価部の計算手続きが増える他は、従来の古典分子動力学計算と同程度の計算量となるため、数千個以上の原子から構成される大規模系の分子動力学計算が、短時間で効率的に実行できる。
【００２８】
本発明によるシミュレーション方法において、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度の相互作用が通常弱く、熱エネルギーの交換が起こりにくいため、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御することが好ましい。
【００２９】
また、本発明によるシミュレーション方法の１階微分計算部において、１階微分を計算する際、動的結合次数評価部における計算結果を繰り返し参照する必要があるため、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑えることが好ましい。
【００３０】
さらに、本発明によるシミュレーション方法において、構成元素の数を増やしていくほど、必要となるポテンシャル定数が増えていくが、シミュレーションに前もって全ての相互作用項のポテンシャル定数を用意する必要はない。シミュレーションで扱う系毎に、必要となる２体相互作用項や３体相互作用項などの種類が限定されることが多いため、これら最小限必要な相互作用項のポテンシャル定数のみで計算を実行しても良い。その場合、シミュレーション実行中に結合の組み換えが起こり、最初にポテンシャル定数を用意していなかった相互作用項が新たに必要となる場面が生じうるが、その際は、シミュレーションを一時中断し、新たに必要となったポテンシャル定数の入力を待って、シミュレーションを再開するようにしてもよい。
【００３１】
この発明の第２の主要な観点によれば、コンピュータシステムにより実行されるコンピュータソフトウェアプログラムであって、記録媒体と、前記記録媒体に格納され、前記コンピュータシステムに原子核座標及び内部座標の情報を格納させる座標情報格納部と、前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標及び内部座標に基いて原子間の結合次数を計算させる動的結合次数評価部と、前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算させるポテンシャルエネルギー計算部と、前記記録媒体に格納され、前記コンピュータシステムに前記ポテンシャルエネルギーを１階微分処理することによって各原子核座標および内部座標に作用する力を計算させる１階微分計算部と、前記記録媒体に格納され、前記コンピュータシステムに前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算させる圧力テンソル計算部と、前記記録媒体に格納され、前記コンピュータシステムに、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるシミュレーション部と、前記記録媒体に格納され、前記コンピュータシステムに前記シミュレーション結果を出力させる出力部とを有することを特徴とするコンピュータソフトウエアプログラムが提供される。
【００３２】
このような構成によれば、コンピュータシステムにインストールされることで、前記第１の主要な観点による装置を構成することができ、同様の作用効果をえることができる。
【００３３】
また、この発明の第３の主要な観点によれば、原子核座標及び内部座標の情報を格納する座標情報格納工程と、前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価工程と、前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算工程と、前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算工程と、前記微分計算工程で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算工程と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション工程と、前記シミュレーション結果を出力する出力工程と
を有することを特徴とする分子動力学シミュレーション方法が提供される。
【００３４】
このような構成によれば、前記第１の主要な観点によるシミュレーション装置で実行可能なシミュレーション方法が提供される。
【００３５】
なお、ここで、説明していない本発明の他の特徴及び顕著な作用及び効果は、次に説明する発明を実施するための最良の形態の項及び図面から当業者が容易に理解することができる。
【発明を実施するための最良の形態】
【００３６】
（シミュレーション装置）
図１は、本発明のシミュレーション装置の機能的な構成を示すブロック図である。このシミュレーション装置は、図１に示すように、データ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０を少なくとも具備し、一連の多体分子シミュレーションの機能を有した処理制御部４１と、操作者からのデータや命令などの入力を受け付ける入力部４２と、シミュレーション結果を出力する出力部４５と、原子間相互作用ポテンシャルの定数、および原子核座標と各原子の内部座標などの入力データなどを格納したデータ記憶部４３と、シミュレーションプログラムなどを格納したプログラム記憶部４４とから少なくとも構成されている。
【００３７】
データ解析部１０は、構造評価部１１、統計情報計算部１２、構造可視化部１３等のような、原子座標およびその時系列データの解析機能を少なくとも備えている。例えば、構造評価部１１は、原子間相互作用計算部の原子座標情報管理部の情報から、動径分布関数等を計算する機能を有する。また、原子間相互作用計算部の動的結合次数評価部で得られる各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う機能を有する。統計情報計算部１２は、原子間相互作用計算部の動的結合次数評価部で得られる各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う機能を有する。更に、シミュレーション実行部の古典分子動力学計算部で得られる各原子核座標の時系列情報から、拡散係数などの輸送係数を計算することができる。構造可視化部１３は、原子間相互作用計算部の原子核座標情報管理部に格納されている原子核座標および動的結合次数評価部で得られる各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる機能を有する。結合次数に応じて原子間結合手の太さを変えて表示させることができ、各原子間の結合の強さを直観的に把握できる。
【００３８】
シミュレーション実行部２０は、古典分子動力学計算部２１と、構造最適化計算部２２を少なくとも具備している。古典分子動力学計算部２１は、原子間相互作用計算部３０で計算される原子核座標および各原子の内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する。構造最適化計算部２２は、原子間相互作用計算部３０で計算される原子核座標および各原子の内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する。
【００３９】
原子間相互作用計算部３０は、原子核座標情報管理部３１、内部座標情報管理部３２、動的結合次数評価部３３、ポテンシャルエネルギー計算部３４、１階微分計算部３５、圧力テンソル計算部３６を少なくとも有する。原子核座標情報管理部３１は、原子核座標およびその時間微分の数値を保持する。内部座標情報管理部３２は、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を保持する。動的結合次数評価部３３は、原子核座標と内部座標の値から各原子間の結合次数を評価する機能を有する。ポテンシャルエネルギー計算部３４は、原子核座標と内部座標の値、さらに動的結合次数評価部３３で評価された結合次数から、系の全ポテンシャルエネルギーを計算する。１階微分計算部３５は、原子核座標と内部座標の値、さらに動的結合次数評価部３３で評価された結合次数から、原子核および各原子の内部自由度に加わる力を計算する。圧力テンソル計算部３６は、１階微分計算部３５で計算される原子間力から、系の内部応力テンソルを計算する。
【００４０】
図１において、入力部４２はキーボード、マウス、ライトペン、フレキシブルディスク装置、またはＬＡＮ端子などで構成される。データ記憶部４３およびプログラム記憶部４４は、周知の磁気テープ、磁気ドラム、磁気ディスク、光ディスク、光磁気ディスク、あるいはＲＯＭ、ＲＡＭなどの半導体メモリー等を用いた記録部である。また出力部４５は、ディスプレイ装置やプリンタ装置などにより構成されている。図１に示した処理制御部４１、データ記憶部４３、およびプログラム記憶部４４等はＣＰＵ，およびこのＣＰＵに接続されたＲＯＭ、ＲＡＭ、磁気ディスクなどの記憶装置を含む通常のコンピュータシステムで構成すればよい。図１に示した処理制御部４１で実行される各処理の入力データはデータ記憶部４３に格納され、プログラム命令はプログラム記憶部４４に記憶される。そしてこれらの入力データやプログラム命令は必要に応じてＣＰＵに読み込まれ、演算処理が実行されるとともに、データ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０で発生した数値情報などのデータはデータ記憶部４３に格納される。
【００４１】
次に、原子間相互作用計算部３０の動作について説明する。原子間相互作用計算部３０においては、まず動的結合次数評価部３３により各原子間の結合次数が評価される。その結果を用いてポテンシャルエネルギー計算部３４において系の全ポテンシャルエネルギーが計算され、１階微分計算分３５において各原子核および内部自由度に加わる力が計算される。また、１階微分計算分３５で計算された原子間力から、圧力テンソル計算部３６において系の内部応力テンソルが計算される。
【００４２】
各原子間の結合次数を判定する動的結合次数評価部３３の機能的な構成を示すブロック図を図２に示す。図２に示すように、動的結合次数評価部３３は、結合手の重なり判定部３３１、結合次数判定部３３２、及びスイッチング関数部３３３を少なくとも有している。
【００４３】
図２において結合手の重なり判定部３３１は、原子核座標情報管理部３１及び内部座標情報管理部３２より座標情報を読取り、各原子間の結合手の重なり度を表すスカラ量Ψijを、図２の式（＊１）乃至式（＊３）で計算する。式（＊１）乃至式（＊３）でｉ，ｊは各原子を示す指数（インデックス）で、α，βはそれぞれ、原子ｉ，ｊに属する結合手を示す指数である。また式（＊３）で定義される関数ψの添え字Ｉ，Ｊは、原子i，ｊの種類を示す。ＣＩ，ｐＩ，ｎＩは、Iで指定される原子の種類に固有の定数である。ｓ，θ，φは、原子の内部自由度変数を示し、式（＊１）に見られるそれらの添え字は、各内部自由度変数が属する原子の指数と結合手の指数を示している。ｒijは原子iと原子ｊの間の距離を示し、ｅijは原子iの原子核位置を始点とし原子ｊの原子核位置を終点とするベクトルの単位ベクトルを示す。結合手の重なり度は、式（＊１）に示すように、原子iと原子ｊの間の各結合手の組み合わせについて式（＊２）で示される関数Ｓを計算し、その総和として定義される。式（＊２）において、関数ｆcはスイッチング関数部３３３で評価され、それは式（＊５）で定義される。近距離では１、遠距離では０となる連続関数である。式（＊５）Ｒ，Ｄは、スイッチングを起こす距離を規定するパラメータで、それらは式（＊２）において、２原子の組み合わせＩ，Ｊ毎に指定できる。
【００４４】
図２において結合次数判定部３３２は、結合手の重なり判定部３３１で決定されるスカラ量Ψijを引数としてスイッチング関数部を呼び出し、式（＊４）に従って計算する。式（＊４）におけるＲ，Ｄは、結合次数判定部３３２固有のパラメータである。
【００４５】
以上の動的結合次数評価部３３における結合次数の決定方法を直観的に把握するための模式図を図３に示す。原子ｉと原子ｊ間の結合次数ｇ(Ψij)は、原子ｉと原子ｊ間の結合手の重なり度を表すスカラ量Ψijに応じて０から１まで連続的に増加する。この結合手の重なり度Ψijは、図３のＡ及びＢに示すように、原子ｉと原子ｊ間の距離と、それらの結合手の相対的な向きによって定まる。図３のＡ及びＢに描かれている結合手の形状は、式（＊３）で定義される関数ψを正確に反映している。即ち、式（＊３）で定義された関数ψを、引数である単位ベクトルｅの方位を変え、対応する関数値を単位ベクトルｅの向きをもつ動径の長さで表すと、動径の先端をプロットして描かれる曲面は、Ａ及びＢに示したような偏った亜鈴形を成す。もし、２原子が相互作用できる距離に近づいても、Ａに示したように、互いの結合手に重なりが生じていない場合には、式（＊１）乃至式（＊３）より、Ψijは負の値となり、結合次数ｇ(Ψij)は０と判定される。一方、Ｂに示したように、互いの結合手に重なりが生じている場合には、式（＊１）乃至式（＊３）より、Ψijは正の値を持つと評価され、結合次数ｇ(Ψij)は１と判定される。尚、図３には描かれていないが、動的結合次数項ｇ(Ψij)はΨijの単調増加関数ではなく、Ψijが一定値を超えると再び０に減少する式（＊４）で定義されている。これは、相手の原子を完全に覆い尽くすほど異常に結合手が伸びた場合、相手の原子の更に先に存在する原子とも同時に結合して不自然な安定化を起こす可能性があるため、これを防ぐための処置として、Ψijが一定値を超えると再びｇ(Ψij)が０に減少するようになっている。
【００４６】
系の全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部３４の機能的な構成を示すブロック図を図４に示す。図４に示すように、ポテンシャルエネルギー計算部３４は、１体項計算部３４１、２体項計算部３４２、３体項計算部３４３、及び全エネルギー合算部３４４を少なくとも有している。
【００４７】
図４において１体項計算部３４１は、内部座標情報管理部３２より座標情報を読取り、各結合手の第一自由度変数sの値から、原子ｉの１体ポテンシャルエネルギーｆ１（ｉ）を図２の式（＊６）で計算する。式（＊６）でｉは原子を示す指数で、αは原子ｉに属する結合手を示す指数である。ωｉは、原子ｉに属する結合手の本数を示す。添え字Ｉは、原子ｉの種類を示す。ａＩ，ｍＩは、Iで指定される原子の種類に固有の定数である。
【００４８】
図４において２体項計算部３４２は、原子核座標情報管理部３１より座標情報を読取り、２原子間距離のみで決まる２体ポテンシャルエネルギーｆ２typeを計算する。ｆ２の上付き添え字typeは、２体ポテンシャルエネルギーの種類を指定する標識で、本シミュレーション装置では、等方的に働く反発相互作用を表す添え字repと、結合が組まれた原子間のみに作用する引力相互作用を表す添え字atrの２種類が用意されている。２体ポテンシャルエネルギーｆ２typeは、式（＊７）に示すように、式（＊８）で定義される関数fsw2の２項の和で定義される。式（＊８）で、ｒは２原子間距離を示し、A，ｐ，ａはパラメータを示す。これらのパラメータは式（＊７）に添え字付きで示されているように、２つの原子の種類Ｉ，Ｊと２体ポテンシャルエネルギーの種類によって固有の値を取り得る。特にａはカットオフ距離と呼び、距離の次元を有し、引力相互作用が完全に打ち切られて０になる距離を指定する。
【００４９】
図４において３体項計算部３４３は、原子核座標情報管理部３１より座標情報を読取り、３原子間の位置関係で決まる３体ポテンシャルエネルギーｆ3を計算する。３体ポテンシャルエネルギーｆ3は、式（＊９）に示すように、式（＊１０）で定義される３体関数Λと式（＊１１）で定義される3体関数Θの積で定義される。３体関数Λは、原子i−ｊ間の距離ｒｉｊと原子i−ｋ間の距離ｒｉｋの関数で、λ，γ，ａはパラメータである。各パラメータの下付き添え字は、３体項を構成する原子の種類を指定する。例えば、下付き添え字JIKは、Ｉで指定される原子の種類を頂点とし、それがJおよびＫで指定される原子と結合を組んで結合角を構成している場合を表す。また、γ，ａの上付き添え字は、上記結合角を張る２つの結合の一方を指定するための標識である。ａはカットオフ距離を表し、３体ポテンシャルエネルギーが作用する範囲を指定する。また、3体関数Θは、結合角の角度θjikの関数で、θ0，αはパラメータを示す。各パラメータの下付き添え字は、３体項を構成する原子の種類を指定する。特にパラメータθ0JIKは、Ｉで指定される原子の種類を頂点とし、それがJおよびＫで指定される原子と結合を組むことで構成される結合角の再安定角度を決定する。
【００５０】
図４において、全エネルギー合算部３４４は、１体項計算部３４１、２体項計算部３４２、及び３体項計算部３４３で計算されたｆ１，ｆ２rep，ｆ２atr，ｆ3を式（＊１２）のように合計し、系の全ポテンシャルエネルギーΦを計算する。式（＊１２）において、ｉ，ｊ，ｋは各原子を示す指数で、Ｎは系に含まれる全粒子数を表す。またｇ(Ψ)は、動的結合次数評価部３３で計算される動的結合次数項を示す。式（＊１２）に示されるように、共有結合に起因する引力相互作用ｆ２atrと３体項ｆ3が、動的結合次数(Ψ)が０でない有限の値となる原子間に限定して適用される。このため、企図しない原子間に不要な引力相互作用が発生することがなく、また企図しない原子間に３体項が働いて余計なポテンシャル障壁が発生するのを防ぐことができる。動的結合次数(Ψ)は、各原子の内部座標変数に応じて連続的に変化するため、共有結合の所在を明示しつつその繋ぎ変えを滑らかに起こすことを可能にしている。
【００５１】
各原子核および内部自由度に加わる力を計算する、１階微分計算部３５の機能的な構成を示すブロック図を図５に示す。図５に示すように、１階微分計算部３５は、原子核座標に関する微分実行部３５１、及び内部座標に関する微分実行部３５２を少なくとも有している。
【００５２】
図５において原子核座標に関する微分実行部３５１は、動的結合次数評価部３３及びポテンシャルエネルギー計算部３４で得られる計算結果を利用し、式（＊１３）より全ポテンシャルエネルギーの各原子核位置に関する１階偏微分係数を計算する。式（＊１３）において、ｉ，ｊ，ｋは、各原子を示す指数である。
【００５３】
図５において内部座標に関する微分実行部３５２は、動的結合次数評価部３３及びポテンシャルエネルギー計算部３４で得られる計算結果を利用し、式（＊１４）より全ポテンシャルエネルギーの各内部座標変数に関する１階偏微分係数を計算する。式（＊１４）において、ξiαは、式（＊１）で示した、原子ｉに属する結合手αの、３つの内部自由度変数siα，θiα，φiαをまとめたベクトル変数である。siα，θiα，φiαはこの順番で、ξiαの第一要素、第二要素、第三要素に対応する。
【００５４】
１階微分計算部３５で原子核座標及び内部座標に関するポテンシャルの１階微分が計算されると、それらを質量で除算することで、直ちに原子核及び内部座標変数に加わる力を求めることができる。ただし、ここで言う質量とは、原子核に加わる力を求める際はkgの次元を有する原子核の質量そのものであるが、内部座標変数に加わる力を求める際は、内部座標変数用に別途指定される、エネルギーの次元を有する仮想的な質量を用いる。質量は一般に、運動する自由度の変化の容易さを決めるパラメータであり、その次元はkgに限らない。本シミュレーション装置において内部座標変数を原子核に素早く追従させたい場合には、原子核の質量に比べて十分小さな質量を指定する。ただし、内部座標変数に対してあまり小さな質量を指定すると、計算の精度低下を招くため、扱う系に応じて適切な質量の値を使用者が指定することになる。
【００５５】
系の内部応力テンソルを計算する、圧力テンソル計算部３６の機能的な構成を示すブロック図を図６に示す。図６に示すように、圧力テンソル計算部３６は、各２原子間に働く力を計算する２原子間力計算部３６１、及び圧力テンソル合算部３６２を少なくとも有している。
【００５６】
図６において２原子間力計算部３６１は、１階微分計算部３５で得られる計算結果を利用し、式（＊１５）より、原子ｉが原子ｊより受ける力を、全ての原子ｉ，ｊの組み合わせについて計算する。
【００５７】
図６において圧力テンソル合算部３６２は、２原子間力計算部３６１で得られる計算結果を利用し、式（＊１６）より、系の全内部応力テンソル（圧力テンソル）を計算する。式（＊１６）において|vi>は、原子iの速度ベクトルを示す。各原子核の速度ベクトルは、シミュレーション実行部２０の古典分子動力学計算部２１において逐次生成され、原子核座標情報管理部３１に格納される。シミュレーション実行部２０の構造最適化計算部２２を実行する場合は、この速度ベクトルは生成されず、原子核座標情報管理部３１において速度ベクトル情報を格納するメモリには常に数値０がセットされる。従って、構造最適化計算部２２を実行する場合の系の全内部応力テンソルは、式（＊１６）の中括弧内のΣΣ|fij><rij|の項のみが有効となる。
【００５８】
図６において圧力テンソル合算部３６２では、内部自由度変数の指数は表れていない。これは、内部自由度変数が各原子核位置を中心とする極座標系を運動しており、圧力テンソルの起源である運動量輸送は、原子間相互作用と原子核の移動によって行われると考えられるためである。従って、圧力テンソルは従来法通り、原子核の運動エネルギーと原子間相互作用のみで定義し、内部自由度変数の運動エネルギー分は含まれない。
【００５９】
以上のように、原子間相互作用計算部３０によれば、原子核同様、各原子の共有結合手の方向を決定する内部自由度変数にも力が働くため、シミュレーション実行部２０によって、その力に基づいて結合手の方向を変化させることができる。これにより、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができ、共有結合の組み換えを伴う分子シミュレーションを行うことが可能となる。この拡張に伴って増える計算は、内部座標情報管理部３２のメモリ操作及び動的結合次数評価部３３に止まるため、従来の古典分子動力学計算と同様、数千個以上の原子から構成される大規模系のシミュレーションが、短時間で効率的に実行できる。
【００６０】
また、原子間相互作用計算部３０によれば、２体引力相互作用および３体相互作用の発現を結合が組まれている原子間に限定することができる。このため、企図しない原子間に余計な引力やポテンシャル障壁が形成される事を防ぐことができる。更に、結合を組んでいない原子間の２体引力相互作用を打ち消すように３体相互作用項を調節する必要がなくなるため、２体ポテンシャル定数と３体ポテンシャル定数を独立に決定することができ、多元素混在系用のポテンシャル定数を用意することが可能となる。
【００６１】
(シミュレーション方法)
図１に示した原子間相互作用計算部３０を用いて分子動力学シミュレーションを実行するには、所望の統計アンサンブルを生成するように設計されたオイラー−ラグランジュ方程式の数値積分を行えばよい。また、構造最適化計算を実行するには、原子間相互作用計算部３０で計算される力の情報を用いて、共役勾配法等、多次元空間上のポテンシャルの極小点探索アルゴリズムを実行すればよい。以下において、分子動力学計算及び構造最適化計算を行うための具体例、即ち、シミュレーション方法１及びシミュレーション方法２を示す。
【００６２】
シミュレーション方法１：図７のフローチャートを参照して、本発明のシミュレーション方法１を説明する。
（イ）まず、ステップＳ１１０で、分子動力学計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標を入力部４２から読み込む。ここで分子動力学計算を実行する上で必要となる制御データとは、数値積分を実行する際の1ステップ当たりの時間刻み、全ステップ数など、分子動力学計算を実行する上で必要となるデータである。また温度制御や圧力制御アルゴリズムを使用するかどうかの選択や、その際の設定温度及び設定圧力の数値なども含む。分子動力学計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標の情報は基本的に全て入力部４２より読み込むが、入力が省略された情報については、必要に応じて自動的に補完する。例えば、内部座標情報が省略された場合、原子核座標情報に基づいて、各原子の最も近い原子から順に結合を組むように、内部自由度の座標を自動的に決定する。
（ロ）次にステップＳ１２０で、原子核座標及び内部座標の情報から、各原子間の動的結合次数を評価し、原子間結合のトポロジーを定める。
（ハ）次にステップＳ１３０で、上記原子間結合のトポロジーに含まれる２体ポテンシャルおよび３体ポテンシャルの定数が完備されているか判定する。完備されていない時はステップＳ１４０に進み、不足しているポテンシャル定数を要求するエラーメッセージを出力し、計算の途中データを記録して終了する。必要なポテンシャル定数が揃っている時はステップＳ１５０に進む。
（ニ）ステップＳ１５０ではポテンシャルエネルギー、１階微分、圧力テンソルを計算する。
（ホ）ステップＳ１６０では原子座標、内部座標に関する１階微分の値から力を計算し、これをオイラー−ラグランジュ方程式に代入して各座標の加速度を計算し、分子動力学計算を１ステップ進める。
（ヘ）ステップＳ１７０では入力データで指定された全ステップ数に達したか判定する。達していない場合はＳ１２０に戻り、達している場合には計算結果を出力して終了する。
【００６３】
シミュレーション方法２：図８のフローチャートを参照して、本発明のシミュレーション方法２を説明する。
（イ）まず、ステップＳ２１０で、構造最適化計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標を入力部４２から読み込む。ここで構造最適化計算を実行する上で必要となる制御データとは、構造最適化計算の収束判定条件など、構造最適化計算を実行する上で必要となるデータである。また、周期的境界条件を設けて無限に大きな系の部分系を模す場合、部分系を構成する単位セルのサイズ及び形状も最適化できるが、この単位セルのサイズ及び形状の変化を許すか否かを示す制御情報も含む。構造最適化計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標の情報は基本的に全て入力部４２より読み込むが、入力が省略された情報については、必要に応じて自動的に補完する。例えば、内部座標情報が省略された場合、原子核座標情報に基づいて、各原子の最も近い原子から順に結合を組むように、内部自由度の座標を自動的に決定する。
（ロ）次にステップＳ２２０で、原子核座標及び内部座標の情報から、各原子間の動的結合次数を評価し、原子間結合のトポロジーを定める。
（ハ）次にステップＳ２３０で、上記原子間結合のトポロジーに含まれる２体ポテンシャルおよび３体ポテンシャルの定数が完備されているか判定する。完備されていない時はステップＳ２４０に進み、不足しているポテンシャル定数を要求するエラーメッセージを出力し、計算の途中データを記録して終了する。必要なポテンシャル定数が揃っている時はステップＳ２５０に進む。
（ニ）ステップＳ２５０ではポテンシャルエネルギー、１階微分、圧力テンソルを計算する。
（ホ）ステップＳ２６０では、次の１次元極小点探索の方向を決定する。ここで言う方向とは、原子核座標および各原子の内部自由度座標を併せた全自由度から構成される多次元空間内における方向ベクトルを指す。最初にこのステップに来た場合には１階微分から得られる力の方向を指定し、２回目以降にこのステップに来た場合は、既に探索した方向に共役な方向を、共役勾配法の手続きで計算する。
（ヘ）ステップＳ２７０では、前のステップで指定された、原子核座標および各原子の内部自由度座標を併せた全自由度から構成される多次元空間内における方向ベクトル上で、１次元の極小点探索を実行する。
（ト）ステップＳ２８０では、構造最適化計算の収束判定条件が満たされるか判定する。収束判定条件とは、具体的には、ステップＳ２７０の前後で、原子核座標および内部自由度座標の移動幅を計算し、その合計が入力データで指定した範囲に収まっていることを指す。収束判定条件が満たされない時は、ステップＳ２２０に戻る。収束判定条件が満たされた時は、計算結果を出力して終了する。
【００６４】
以上のシミュレーション方法１及びシミュレーション方法２によれば、共有結合の所在を明示しながら、その滑らかな繋ぎ変えを可能にする分子動力学計算及び構造最適化計算を行うことができる。共有結合の発現を明示された原子間に限定できるため、従来法で問題となっていた余計な相互作用が取り除かれ、多元素が混在する複雑な系のポテンシャル定数を用意することが可能となる。つまり、シミュレーション方法１及びシミュレーション方法２によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子シミュレーションを、数千原子以上の大きな規模で実行できる。
【００６５】
（シミュレーションプログラムを格納した記録媒体）
図７及び図８に示すシミュレーション方法を含め、図１に示したデータ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０におけるシミュレーション方法を実現するためのプログラムはコンピュータ読取り可能な記録媒体に保存しておいても良い。この記録媒体をコンピュータシステムによって読み込ませ、図１に示したプログラム記憶部４４に格納し、このプログラムを処理制御部４１で実行して本シミュレーション方法を実現することもできる。ここで、記録媒体とは、例えば磁気ディスク、光ディスク、光磁気ディスク、磁気テープなどのプログラムを記録することができるような媒体などが含まれる。コンピュータの外部メモリも、ここで言う記録媒体に含まれる。
【００６６】
（実施例）
本発明の実施例として、結合組み換えを伴う３元素混在系の分子動力学計算、原子核系と内部自由度の温度制御分子動力学計算、Ｓｉ結晶モデルの圧力制御分子動力学計算、シリコン酸化膜とシリコン界面のモデリングを示す。
【００６７】
図９はＳｉＯ４Ｈ４分子の分子動力学シミュレーションを実施した際の初期構造（0ステップ）、2000ステップ、4000ステップ、5000ステップ、6000ステップ、8000ステップ経過後の構造を示す。各原子を球で表示し、各原子の結合手に相当する内部自由度を円柱で表示している。結合手１本あたり３次元の自由度を有し、円柱の配向が２自由度、長さが１自由度に対応する。結合手の本数はＳｉ原子は４、Ｏ原子は２、Ｈ原子は１とした。この表示では、例えば0ステップ目の状態で結合手は重なっていないように見えるが、計算上は向かいあった円柱の間に結合が組まれている。適当な初速度を与えて分子動力学計算を開始すると、分子動力学計算のステップが進むにつれて、原子核及び結合手の向き、長さが変化していく。向かい合った円柱は、結合が切れないよう、向かい合った状態を維持するように変化していく。これは、結合が組まれているときのみ、２体引力相互作用による安定化が起こるためである。5000ステップ目ではＡで示したＯ原子の結合手がそれまで結合していたＨ原子と反対側に結合手を伸ばし、一時的にＯ−Ｈ結合の切断が起こったが、間もなくＨ原子と再び結合しなおす様子が再現された。初期温度をより高く設定すると、原子間結合が完全に切断され、原子核が遠くに離れていく様子も再現される。この小規模の分子動力学計算から、本発明により、２元素以上の多元素混在系において、結合の組み換えを伴う分子動力学計算が可能であることが示される。
【００６８】
図１０は、図９で示した分子動力学計算中の、ハミルトニアンの変動を示す。古典力学では、運動エネルギーとポテンシャルエネルギーの総和であるハミルトニアンは保存量であり、このハミルトニアンの保存を確かめることで、力の計算、即ち、ポテンシャル関数の１階微分が正しいか否かを判定できる。図１０では、時間刻み0.00008τ、0.00005τ、0.00001τで行った分子動力学計算におけるハミルトニアンの変動を示している。時間刻みが小さくなるほど、離散化誤差の影響が減少し、ハミルトニアンの変動が抑えられ、一定の値に保存していくことがわかる。
【００６９】
図１１は、温度制御分子動力学計算を実施した際の、原子核系の温度、内部自由度系の温度、及び原子核と内部自由度系を合わせた全系の温度の時間変化を表す。ここで温度とは、粒子の平均運動エネルギーから次の式(１)で定義される量である。
【００７０】
Ｔ ＝ （2／3kB）<1/2 m v2> (１)
【００７１】
ここでTは温度、kBはボルツマン定数、mは原子の質量、vは原子の速度を表す。内部自由度系の温度も同様に、次の式(２)で定義する。
【００７２】
Ｔ ＝ （2／3kB）<1/2 μ (dξ/dt)2> (２)
【００７３】
ここでμは内部自由度変数の仮想的な質量、ξは内部自由度変数ベクトルを表す。図１１(a)は、原子核系と内部自由度系を合わせた全体の温度が、能勢の温度制御アルゴリズムにより計算開始時から1000Kに保持されている。しかし、原子核系と内部自由度系のそれぞれの温度は、初期状態で原子核系は約2000K、内部自由度系は500K程度と大きな差があり、両者が熱平衡状態に達するまでに長い分子動力学計算を必要とする。これは、原子核座標系と内部自由度座標系の相互作用が、動的結合次数項を介してしか起こらないため、結合組み換え時の極限られた時にしか運動エネルギーを交換できないからである。
【００７４】
原子核座標系と内部自由度座標系の温度を所望の温度に制御するには、同時に複数の熱浴を設けて温度制御が行える、能勢-Hoover形式を使えばよい。図１１(b)は、能勢-Hoover形式を用いて、原子核座標系と内部自由度座標系の温度を1000Kに制御した例を示す。初期状態では、図１１aと同様、原子核座標系と内部自由度座標系の温度に大きな差があるが、この差は直ちに解消され、計算開始直後から平均温度が1000Kに制御されていることがわかる。この能勢-Hoover形式を用いれば、原子核座標系と内部自由度座標系を別々の温度で制御することも可能である。
【００７５】
圧力テンソルの起源である運動量輸送は、原子間相互作用と原子核の移動によって行われると考えられるため、本発明においても、圧力テンソルは従来法通り、原子核の運動エネルギーと原子間相互作用のみで定義し、内部自由度変数の運動エネルギー分は含めていない。この方法でも、正しく圧力制御の分子動力学計算を実施できる。
【００７６】
図１２は、圧力制御の分子動力学計算で用いたＳｉ結晶５７６原子モデルである。結晶モデルを囲む立方体は、本計算で適用した３次元周期的境界条件の単位セルを表す。パリネロ−ラーマン(Parrinello-Rahman)の圧力制御分子動力学法により、この単位法セルの大きさを調節して、内部の圧力を設定圧力5Gpaに制御した結果を図１３に示す。Aは、各瞬間で計算した圧力テンソルの対角成分の平均値をプロットしたものである。変動の幅が大きく、設定圧力に保たれている様子を確認することが難しいが、Ｂのように計算開始時からの平均値をプロットすると、設定圧力5Gpaに正確に保持されていることが確認できる。
【００７７】
従来の結合組み換え型多体ポテンシャルでは、余計なポテンシャル障壁により圧縮ストレスを過大評価してしまう問題が生じたが、共有結合による相互作用を結合手が組まれた原子間に限定できる本発明では、この問題が解決される。この効果を確認するため、圧縮ストレスが顕著に観察される例として、シリコン基板上のシリコン酸化膜構造を選び、従来の結合組み換え型多体ポテンシャル、及び本発明の動的結合手を含む原子間相互作用ポテンシャルで上記シリコン酸化膜構造モデルを最適化し、最適化後の構造、及び残留ストレスの値を比較した。
【００７８】
Ｓｉ酸化膜モデルは、まずＳｉ単結晶モデルを作成し、そこにＯ原子を徐々に導入して作成する。初期のＳｉ単結晶構造は4800個のSi原子から構成される、表面が(001)面を向いた単結晶モデルで、幅及び奥行が3.84nm、厚さが6.51nmに相当する。幅及び奥行方向に２次元の周期的境界条件を設けている。このＳｉ単結晶モデルの表面から、Ｏ原子を１個ずつ、Si-Si結合の中点に挿入し、これを繰り返してＳｉ酸化膜部分を形成する。Ｏ原子を挿入するサイトは面内で無作為に選ぶが、１原子層の酸化が完了するまで次の層へ導入はしない。この酸化方法は、Ｓｉ結晶表面の酸化が１原子層ずつ進行することを示す実験事実に基づいている。Ｏ原子を１個導入する度に、構造最適化計算を行う。また、Ｏ原子を２０個導入する度に、温度制御分子動力学計算を20000ステップ実行し、構造のアニール（焼き鈍し）を行う。これは、高温のアニールで結合の組み換えを促進し、酸化膜構造のアモルファス化及び欠陥の消滅を起こすためである。設定温度は800℃（1073K）とする。時間刻みは2.30×10-5ピコ秒とした。分子動力学計算の数値積分アルゴリズムには、5値のGearの予測子・修正子法を用いた。
【００７９】
図１４(a)に、Si,O混在系用に開発された従来の結合組み換え型多体ポテンシャルを用いて、上記の方法で作成したＳｉ酸化膜とＳｉ基板の界面モデルを示す。酸化膜は2800個のＯ原子を含み、膜厚は約3.0nmに相当する。この酸化膜モデルのストレスは、約13GPaと見積もられた。酸化膜部のストレスは、積層構造全体の内部応力テンソルの対角成分の、幅及び奥行方向の成分の平均値とする。内部応力テンソルを計算する際の体積は、酸化膜部分の体積 (3.84 x 3.84 x 3.0 nm3)を用いた。酸化膜の面内ストレスは、実験では0.4GPaから0.7GPaと報告されているが、これら実験値に比べて著しく大きい。これは、O-Si-O３体項等が余計なポテンシャル障壁を形成し、構造の圧縮を強く阻んでいるためと考えられる。
【００８０】
図１４(b)は、本発明におけるシミュレーション方法で、本発明におけるシミュレーション装置を用いて、図１４(a)で得られたＳｉ酸化膜とＳｉ基板の界面構造をアニール及び構造最適化した結果得られた構造を示す。アニールは、温度制御分子動力学計算を、上記の条件で100000ステップ実行することで行われた。図１４(a)及び図１４(b)の構造は、一見大きな違いが見られないが、余計なポテンシャル障壁の消滅により、本発明で作成した図１４(b)の構造の方が、極僅か酸化膜厚が減少している。図１４(b)の酸化膜部のストレスは、約2GPaであり、従来の結合組み換え型多体ポテンシャルに比べて大幅に減少し、実験値の数倍程度の大きさまで落ち着いた。依然実験値よりも大きい理由として、上記分子動力学計算ならびに構造最適化計算において、単位胞サイズを完全に固定している点が挙げられる。実際のＳｉ酸化膜では、酸化膜のストレスにより、下地のSi基板が横方向に反り、若干ではあるが横方向に膨張している。計算では横方向の膨張を完全に抑えていることを鑑みれば、2GPaという酸化膜部のストレスは、十分妥当な値であり、少なくとも従来の結合組み換え型多体ポテンシャルで得られる13GPaと比較すれば、著しい改善である。
【００８１】
図１５は、図１４(a)及び図１４(b)の構造を用いて計算した、２体相関関数を示す。２体相関関数は、
【００８２】
ｇ（ｒ） ＝ < n(r,r+Δr)> / 4πr2ρΔr (3)
で計算される、距離の関数である。< n(r,r+Δr)>は、ある原子を中心とした半径rと半径r+Δrの球殻の間に存在する原子の数を、系の全ての原子を中心に計算してその平均をとったもので、ρは系全体の原子数密度を表す。２体分布関数は位置の相関が無くなると１に近づく。
【００８３】
1.62Å（オングストローム）付近に立つピークは酸化膜中のSi-O結合を表し、2.35Å付近に立つピークは下地のSi基板中のSi-Si結合を表す。2.6Å付近に立つ緩やかなピークは、SiO4正四面体構造のO-O間距離に起因する。3.1Å付近に立つ小さなピークは、酸化膜中の第２近接のSi-Si対を表す。プロットAでは、3.3Åに鋭いピークが見られる。これは解析から、Si-O対に起因するピークであることがわかっている。この3.3Å付近の鋭いピークは、アモルファスSiO2のX線回折から得られる２体相関関数では観察されない不自然なピークである。この不自然なピークは、図14(b)に示した、本発明で作成した構造から計算した２体分布関数を示すプロットBには見られない。その分、3.1Åから3.3Åの成分が増えている。この結果は、従来型ポテンシャルでは余計なポテンシャル障壁によって接近が阻まれていたSi-O対が、本発明におけるシミュレーション方法では、正しく接近できるように改善されたことを示している。
【００８４】
本発明によれば、共有結合のような方位依存性のある引力相互作用の発現箇所を明示しつつ、その滑らかな繋ぎ変えを許す分子動力学シミュレーション及び構造最適化計算が可能となる。
【００８５】
つまり、本発明によれば、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができ、共有結合の組み換えを伴う分子シミュレーションを行うことが可能となる。
【００８６】
また、本発明によれば、従来の古典分子動力学計算と同様、数千個以上の原子から構成される大規模系のシミュレーションが、短時間で効率的に実行できる。
【００８７】
さらに、本発明によれば、企図しない原子間に余計な引力やポテンシャル障壁の形成を防ぐことができる。
【００８８】
さらに、本発明によれば、多体相互作用項同士の干渉を大幅に減少させることができる。これにより、多体相互作用項のポテンシャル定数の決定が大幅に単純化され、多元素が混在する複雑な系のポテンシャル定数を用意することが可能となる。
【００８９】
つまり、本発明によれば、従来法では不可能であった、多元素混在系において、共有結合の組み換えを伴う化学変化の大規模分子シミュレーションが可能となる。
【００９０】
さらに、本発明によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子動力学シミュレーション及び構造最適化計算が可能なシミュレーション方法を提供することができる。
【００９１】
さらに、本発明によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子動力学シミュレーション及び構造最適化計算が可能なシミュレーションプログラムを格納した記録媒体を提供することができる。
【図面の簡単な説明】
【００９２】
【図１】本発明の実施の形態に係るシミュレーション装置の機能的な構成を示すブロック図である。
【図２】図１に示した動的結合次数部３３の機能的な構成を示すブロック図である。
【図３】動的結合次数評価部３３における結合次数の決定方法を模式的に示した図である。
【図４】図１に示したポテンシャルエネルギー計算部３４の機能的な構成を示すブロック図である。
【図５】図１に示した１階微分計算部３５の機能的な構成を示すブロック図である。
【図６】図１に示した圧力テンソル計算部３６の機能的な構成を示すブロック図である。
【図７】本発明のシミュレーション方法１を示すフローチャートである。
【図８】本発明のシミュレーション方法２を示すフローチャートである。
【図９】図９は、３元素が混在するＳｉＯ４Ｈ４分子の分子動力学シミュレーションを実施した際の構造変化を示す。
【図１０】図１０は、図９で示した分子動力学計算中の、ハミルトニアンの変動を示す。
【図１１】図１１は、温度制御分子動力学計算を実施した際の、原子核系の温度、内部自由度系の温度、及び原子核と内部自由度系を合わせた全系の温度の時間変化を表す。
【図１２】Ｓｉ結晶５７６原子モデルを示す図である。
【図１３】図１２のＳｉ結晶モデルの単位法セルの大きさを調節して、内部の圧力を実施した際の圧力の変動を示した図である。
【図１４】図１４(a)は、従来の結合組み換え型多体ポテンシャルを用いて作成したＳｉ酸化膜とＳｉ基板の界面モデルを示す。図１４(b)は、本発明におけるシミュレーション方法で、図１４(a)の界面構造をアニール及び構造最適化した結果得られた構造を示す。
【図１５】図１４(a)及び図１４(b)の構造を用いて計算した、２体相関関数を示した図である。
【符号の説明】
【００９３】
１０ データ解析部
１１ 構造評価部
１２ 統計情報計算部
１３ 構造可視化部
２０ シミュレーション実行部
２１ 古典分子動力学計算部
２２ 構造最適化計算部
３０ 原子間相互作用計算部
３１ 原子核座標情報管理部
３２ 内部座標情報管理部
３３ 動的結合次数評価部
３４ ポテンシャルエネルギー計算部
３５ １階微分計算部
３６ 圧力テンソル計算部
４１ 処理制御部
４２ 入力部
４３ データ記憶部
４４ プログラム記憶部
４５ 出力部
３３１ 結合手の重なり判定部
３３２ 結合次数判定部
３３３ スイッチング関数部
３４１ １体項計算部
３４２ ２体項計算部
３４３ ３体項計算部
３４４ 全エネルギー合算部
３５１ 原子核座標に関する微分実行部
３５２ 内部座標に関する微分実行部
３６１ ２原子間力計算部
３６２ 圧力テンソル計算部
【特許請求の範囲】
【請求項１】
原子核座標及び内部座標の情報を格納する座標情報格納部と、
前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価部と、
前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部と、
前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算部と、
前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算部と、
原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション部と
前記シミュレーション結果を出力する出力部と
を有することを特徴とする分子動力学シミュレーション装置。
【請求項２】
請求項１記載の分子動力学シミュレーション装置において、
前記シミュレーション部は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理部と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理部と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価部と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものであることを特徴とする分子動力学シミュレーション装置。
【請求項３】
請求項１記載の分子動力学シミュレーション装置において、
前記シミュレーション実行部は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算部とを有するものであることを特徴とする分子動力シミュレーション装置。
【請求項４】
請求項１記載の分子動力学シミュレーション装置において、
前記出力部で出力されたデータの解析を行うデータ解析部をさらに有することを特徴とする分子動力学シミュレーション装置。
【請求項５】
請求項４記載の分子動力学シミュレーション装置において、
データ解析部は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項６】
請求項４記載の分子動力学シミュレーション装置において、
データ解析部は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う構造解析部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項７】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項８】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項９】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１０】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１１】
請求項１０記載の分子動力学シミュレーション装置において、
前記構造可視化部は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１２】
請求項１記載の分子動力学シミュレーション装置において、
前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１３】
請求項１２記載の分子動力学シミュレーション装置において、
前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１４】
請求項１記載の分子動力学シミュレーション装置において、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１５】
コンピュータシステムにより実行されるコンピュータソフトウェアプログラムであって、
記録媒体と、
前記記録媒体に格納され、前記コンピュータシステムに原子核座標及び内部座標の情報を格納させる座標情報格納部と、
前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標及び内部座標に基いて原子間の結合次数を計算させる動的結合次数評価部と、
前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算させるポテンシャルエネルギー計算部と、
前記記録媒体に格納され、前記コンピュータシステムに前記ポテンシャルエネルギーを１階微分処理することによって各原子核座標および内部座標に作用する力を計算させる１階微分計算部と、
前記記録媒体に格納され、前記コンピュータシステムに前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算させる圧力テンソル計算部と、
前記記録媒体に格納され、前記コンピュータシステムに、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるシミュレーション部と
前記記録媒体に格納され、前記コンピュータシステムに前記シミュレーション結果を出力させる出力部と
を有することを特徴とするコンピュータソフトウエアプログラム。
【請求項１６】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記シミュレーション部は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理部と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理部と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価部と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項１７】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記シミュレーション実行部は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算部とを有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項１８】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記出力部で出力されたデータの解析を行うデータ解析部をさらに有することを特徴とするコンピュータソフトウエアプログラム。
【請求項１９】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
データ解析部は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２０】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
データ解析部は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う構造解析部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２１】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２２】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２３】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２４】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２５】
請求項２４記載のコンピュータソフトウエアプログラムにおいて、
前記構造可視化部は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２６】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２７】
請求項２６記載のコンピュータソフトウエアプログラムにおいて、
前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものであることを特徴とする分子動力学コンピュータソフトウエアプログラム。
【請求項２８】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２９】
原子核座標及び内部座標の情報を格納する座標情報格納工程と、
前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価工程と、
前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算工程と、
前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算工程と、
前記微分計算工程で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算工程と、
原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション工程と
前記シミュレーション結果を出力する出力工程と
を有することを特徴とする分子動力学シミュレーション方法。
【請求項３０】
請求項２９記載の分子動力学シミュレーション方法において、
前記シミュレーション工程は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理工程と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理工程と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価工程と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものであることを特徴とする分子動力学シミュレーション方法。
【請求項３１】
請求項２９記載の分子動力学シミュレーション方法において、
前記シミュレーション実行工程は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算工程と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算工程とを有するものであることを特徴とする分子動力シミュレーション方法。
【請求項３２】
請求項２９記載の分子動力学シミュレーション方法において、
前記出力工程で出力されたデータの解析を行うデータ解析工程をさらに有することを特徴とする分子動力学シミュレーション方法。
【請求項３３】
請求項３２記載の分子動力学シミュレーション方法において、
データ解析工程は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３４】
請求項３２記載の分子動力学シミュレーション方法において、
データ解析工程は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥シミュレーション方法。
【請求項３５】
請求項３４記載の分子動力学シミュレーション方法において、
前記データ解析工程は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３６】
請求項３２記載の分子動力学シミュレーション方法において、
前記データ解析工程は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３７】
請求項３２記載の分子動力学シミュレーション方法において、
前記データ解析工程は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３８】
請求項３２記載の分子動力学シミュレーション方法において、
前記データ解析工程は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３９】
請求項３８記載の分子動力学シミュレーション方法において、
前記構造可視化工程は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものであることを特徴とする分子動力学方法。
【請求項４０】
請求項２９記載の分子動力学シミュレーション方法において、
前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項４１】
請求項４０記載の分子動力学シミュレーション方法において、
前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものであることを特徴とする分子動力学シミュレーション方法。
【請求項４２】
請求項２９記載の分子動力学シミュレーション方法において、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項１】
原子核座標及び内部座標の情報を格納する座標情報格納部と、
前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価部と、
前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部と、
前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算部と、
前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算部と、
原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション部と
前記シミュレーション結果を出力する出力部と
を有することを特徴とする分子動力学シミュレーション装置。
【請求項２】
請求項１記載の分子動力学シミュレーション装置において、
前記シミュレーション部は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理部と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理部と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価部と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものであることを特徴とする分子動力学シミュレーション装置。
【請求項３】
請求項１記載の分子動力学シミュレーション装置において、
前記シミュレーション実行部は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算部とを有するものであることを特徴とする分子動力シミュレーション装置。
【請求項４】
請求項１記載の分子動力学シミュレーション装置において、
前記出力部で出力されたデータの解析を行うデータ解析部をさらに有することを特徴とする分子動力学シミュレーション装置。
【請求項５】
請求項４記載の分子動力学シミュレーション装置において、
データ解析部は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項６】
請求項４記載の分子動力学シミュレーション装置において、
データ解析部は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う構造解析部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項７】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項８】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項９】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１０】
請求項４記載の分子動力学シミュレーション装置において、
前記データ解析部は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１１】
請求項１０記載の分子動力学シミュレーション装置において、
前記構造可視化部は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１２】
請求項１記載の分子動力学シミュレーション装置において、
前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１３】
請求項１２記載の分子動力学シミュレーション装置において、
前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１４】
請求項１記載の分子動力学シミュレーション装置において、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものであることを特徴とする分子動力学シミュレーション装置。
【請求項１５】
コンピュータシステムにより実行されるコンピュータソフトウェアプログラムであって、
記録媒体と、
前記記録媒体に格納され、前記コンピュータシステムに原子核座標及び内部座標の情報を格納させる座標情報格納部と、
前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標及び内部座標に基いて原子間の結合次数を計算させる動的結合次数評価部と、
前記記録媒体に格納され、前記コンピュータシステムに前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算させるポテンシャルエネルギー計算部と、
前記記録媒体に格納され、前記コンピュータシステムに前記ポテンシャルエネルギーを１階微分処理することによって各原子核座標および内部座標に作用する力を計算させる１階微分計算部と、
前記記録媒体に格納され、前記コンピュータシステムに前記微分計算部で計算された力に基づいて系の内部応力テンソルを計算させる圧力テンソル計算部と、
前記記録媒体に格納され、前記コンピュータシステムに、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるシミュレーション部と
前記記録媒体に格納され、前記コンピュータシステムに前記シミュレーション結果を出力させる出力部と
を有することを特徴とするコンピュータソフトウエアプログラム。
【請求項１６】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記シミュレーション部は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理部と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理部と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価部と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項１７】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記シミュレーション実行部は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算部とを有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項１８】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記出力部で出力されたデータの解析を行うデータ解析部をさらに有することを特徴とするコンピュータソフトウエアプログラム。
【請求項１９】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
データ解析部は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２０】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
データ解析部は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う構造解析部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２１】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２２】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２３】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２４】
請求項１８記載のコンピュータソフトウエアプログラムにおいて、
前記データ解析部は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２５】
請求項２４記載のコンピュータソフトウエアプログラムにおいて、
前記構造可視化部は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２６】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、
前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２７】
請求項２６記載のコンピュータソフトウエアプログラムにおいて、
前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものであることを特徴とする分子動力学コンピュータソフトウエアプログラム。
【請求項２８】
請求項１５記載のコンピュータソフトウエアプログラムにおいて、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものであることを特徴とするコンピュータソフトウエアプログラム。
【請求項２９】
原子核座標及び内部座標の情報を格納する座標情報格納工程と、
前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価工程と、
前記原子核座標、内部座標及び前記結合次数に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算工程と、
前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算工程と、
前記微分計算工程で計算された力に基づいて系の内部応力テンソルを計算する圧力テンソル計算工程と、
原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシミュレーション工程と
前記シミュレーション結果を出力する出力工程と
を有することを特徴とする分子動力学シミュレーション方法。
【請求項３０】
請求項２９記載の分子動力学シミュレーション方法において、
前記シミュレーション工程は、原子核座標およびその時間微分の数値を格納する原子核座標情報管理工程と、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を格納する内部座標情報管理工程と、原子核座標と内部座標の値から各原子間の結合次数を評価する動的結合次数評価工程と、を有するものであり、前記ポテンシャルエネルギー計算、微分計算、及び圧力テンソル計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行させるものであることを特徴とする分子動力学シミュレーション方法。
【請求項３１】
請求項２９記載の分子動力学シミュレーション方法において、
前記シミュレーション実行工程は、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する古典分子動力学計算工程と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する構造最適化計算工程とを有するものであることを特徴とする分子動力シミュレーション方法。
【請求項３２】
請求項２９記載の分子動力学シミュレーション方法において、
前記出力工程で出力されたデータの解析を行うデータ解析工程をさらに有することを特徴とする分子動力学シミュレーション方法。
【請求項３３】
請求項３２記載の分子動力学シミュレーション方法において、
データ解析工程は、前記原子核座標の数値から動径分布関数及び／若しくはＸ線回折像を計算する構造評価工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３４】
請求項３２記載の分子動力学シミュレーション方法において、
データ解析工程は、前記各原子間の結合次数情報から、結合欠損等の結晶欠陥シミュレーション方法。
【請求項３５】
請求項３４記載の分子動力学シミュレーション方法において、
前記データ解析工程は、前記各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う統計情報計算工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３６】
請求項３２記載の分子動力学シミュレーション方法において、
前記データ解析工程は、各原子核座標の時系列情報を記録することにより、拡散係数、粘性係数、熱伝導度などの輸送係数を計算する統計情報計算工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３７】
請求項３２記載の分子動力学シミュレーション方法において、
前記データ解析工程は、基準振動解析により、赤外吸収スペクトルやラマンスペクトルなどの分光学的性質を計算する基準振動解析工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３８】
請求項３２記載の分子動力学シミュレーション方法において、
前記データ解析工程は、原子核座標および各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる構造可視化工程を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項３９】
請求項３８記載の分子動力学シミュレーション方法において、
前記構造可視化工程は、結合次数に応じて結合手の太さを変えることにより、各原子間の結合の強さを表示するものであることを特徴とする分子動力学方法。
【請求項４０】
請求項２９記載の分子動力学シミュレーション方法において、
前記実行部は、中の平均温度を所望の設定温度に保つ温度制御部を有するものであることを特徴とする分子動力学シミュレーション方法。
【請求項４１】
請求項４０記載の分子動力学シミュレーション方法において、
前記シミュレーション部は、温度を制御した分子動力学シミュレーションを実行する場合、原子核座標と内部自由度を別々の熱浴変数と相互作用させ、原子核座標の温度と内部自由度変数の温度を別々に制御するものであることを特徴とする分子動力学シミュレーション方法。
【請求項４２】
請求項２９記載の分子動力学シミュレーション方法において、前記一階微分計算部は、動的結合次数評価部の結果を一時記憶領域に格納して計算回数を抑える計算回数制御部を有するものであることを特徴とする分子動力学シミュレーション方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【公開番号】特開２００６−１９０２３４（Ｐ２００６−１９０２３４Ａ）
【公開日】平成１８年７月２０日（２００６．７．２０）
【国際特許分類】
【出願番号】特願２００５−５５３９２（Ｐ２００５−５５３９２）
【出願日】平成１７年３月１日（２００５．３．１）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り ２００４年９月１日 社団法人応用物理学会発行の「２００４年（平成１６年）秋季 第６５回 応用物理学会学術講演会講演予稿集 第０分冊、第２分冊」に発表
【出願人】（８９９００００６８）学校法人早稲田大学 (602)
【Ｆターム（参考）】
【公開日】平成１８年７月２０日（２００６．７．２０）
【国際特許分類】
【出願日】平成１７年３月１日（２００５．３．１）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り ２００４年９月１日 社団法人応用物理学会発行の「２００４年（平成１６年）秋季 第６５回 応用物理学会学術講演会講演予稿集 第０分冊、第２分冊」に発表
【出願人】（８９９００００６８）学校法人早稲田大学 (602)
【Ｆターム（参考）】
[ Back to top ]