(57)【要約】      （修正有）【課題】数千個以上の原子集団を扱う大規模系の分子動力学シミュレーションにおいて、多元素が混在する複雑な組成を有する物質の、共有結合の組み換えを伴う構造変化を、効率的に再現しつつ、各原子の配位状態に依存した構造エネルギー変化を忠実に再現する。【解決手段】共有結合部位を決定する役割を持つ内部自由度を持った原子モデルを用いて、原子核座標と内部自由度を併せた拡張力学系の運動を古典分子動力学法で計算する。これにより、共有結合のような方位依存性のある引力相互作用の発現を明示しつつ、その滑らかな繋ぎ換えを表現する。この方法で明示される原子間結合の所在から各原子の配位状態を判定し、その配位状態に応じて異なるポテンシャル定数を採用することにより、従来法では不可能であった、動的に変化する配位状態に依存する構造エネルギー変化を忠実に再現する。【選択図】なし【特許請求の範囲】
【請求項１】
原子核座標及び内部座標の情報を格納する座標情報格納部と、各原子の配位状態に依存して異なる値をとるポテンシャル定数情報の格納部と、前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価部と、
前記結合次数に基づいて原子の配位状態を判定する配位状態判定部と、前記原子核座標、内部座標、前記結合次数及び配位状態に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部と、
前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算部と、
原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、及び微分計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシュミレーション部と、
前記シュミレーション結果を出力する出力部とを有することを特徴とする分子動力学シミュレーション装置。
【請求項２】
共有結合のような方位依存性を持つ相互作用の発現部位を決定する役割を持つ内部自由度を持った原子モデルを用いて、各原子の配位状態を識別し、その配位状態に応じて異なるポテンシャル定数を使い分けられる原子間相互作用モデルを用いる分子動力学シミュレーション手法。
【請求項３】
共有結合のような方位依存性を持つ相互作用の発現部位を決定する役割を持つ内部自由度を持った原子モデルを用いて、各原子の配位状態を識別し、複数の配位状態の重ね合わせ状態を表現することによって、配位状態及びそれに伴うポテンシャル定数の変化を表現する原子間相互作用モデルを用いる分子動力学シミュレーション手法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は半導体集積回路製作等に用いる材料設計支援シミュレーション技術、加工プロセスのシミュレーション技術、あるいは、蛋白質等の生体分子やその他有機化合物質の分子シミュレーション技術に係り、特にナノメートルスケールのサイズを有する数千原子規模以上の系を対象とする大規模分子シミュレーション方法、シミュレーション装置、シミュレーションプログラムを格納した記録媒体に関する。
【背景技術】
【０００２】
半導体集積回路の高集積化に伴い、ナノメートルオーダーの極めて高い精度の微細加工技術が要求されている。例えば、金属−絶縁膜−半導体型電界効果トランジスタの絶縁膜部分は僅か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.
【非特許文献６】D. Brenner, Physical Review B Vol.42, (1990) p.9458.
【非特許文献７】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ポテンシャルやBrennerポテンシャルなど、クラスタ汎関数型に分類されるポテンシャルで一部実現されている。しかし、いずれも、特定の２〜３種類の元素が混在した系に適用範囲が限定され、また、想定しうるあらゆる配位状態を正確に再現できるわけではない。しかも、パラメータの決定が非常に困難で、パラメータフィットの際に想定した状態と異なる系へ適用した場合の予測能力も未知数であった。
【００１１】
上記問題点に鑑み、本発明は、多元素混在系においても結合組み換えを伴う構造変化を正確に再現する、大規模分子動力学シミュレーション装置を提供することを目的とする。
【課題を解決するための手段】
【００１２】
上記課題を解決するため、本発明によれば、原子核座標及び内部座標の情報を格納する座標情報格納部と、各原子の配位状態に依存して異なる値をとるポテンシャル定数情報の格納部と、前記原子核座標及び内部座標に基いて原子間の結合次数を計算する動的結合次数評価部と、前記結合次数に基づいて原子の配位状態を判定する配位状態判定部と、前記原子核座標、内部座標、前記結合次数及び配位状態に基づいて全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部と、前記ポテンシャルエネルギーを１階微分処理することによって、各原子核座標および内部座標に作用する力を計算する１階微分計算部と、原子核座標および各原子の共有結合手の方向を決定する内部自由度に加わる力に基づいて、指定した微小時間間における移動量を計算して前記座標情報を更新し、これに基づいて前記結合次数計算、ポテンシャルエネルギー計算、及び微分計算を指定された回数繰り返しながら、原子核座標と内部自由度変数の時間発展に関する分子動力学計算を実行するシュミレーション部と、前記シュミレーション結果を出力する出力部とを有することを特徴とする分子動力学シミュレーション装置が提供される。
【００１３】
まず、本発明に含まれる原子間相互作用モデルの特徴は、結合手を露に記述しその運動まで表現する点である。そしてその結合手の重なり度に応じて原子間の結合次数が動的に変化し、これによって共有結合の滑らかな繋ぎ変えを実現できる。この動的に変化する結合次数の情報は、各原子の配位状態を判別することに使える。
【００１４】
本発明に含まれる原子間相互作用モデルによれば、各原子の配位状態、即ち、各原子がもつ全ての結合手の結合相手の組み合わせに応じて、異なるポテンシャル定数を使い分けられるようになる。これにより、例えば同じSi-O結合でも、そのSi原子やO原子の他の結合手の結合相手の組み合わせに応じて、異なる結合エネルギーを割り当てることができるようになる。更に、１重結合、２重結合、３重結合の違いまで表現できる。また、結合手の数を予め標準的な結合手の数より多く用意しておけば、超原子価状態も表現できるようになるため、新しい結合を徐々に形成しつつ、古い結合手が徐々に消滅していく、というリアルな化学反応経路を再現できるようになる。
本発明における原子間相互作用ポテンシャルの定義：
【００１５】
【数１】
2008052308.tif 000002
【００１６】
第１項は１体項で、各原子の内部自由度の調整や、イオン化エネルギー等を表す役目をする。第２項は２体項で、等方的に働く反発項U repと、共有結合が生じた部位にのみ働く引力項U atrからなる。残りは３体項、４体項である。４体項は、有機高分子への応用を視野に入れ、捩れ角歪みを表現するための項である。
【００１７】
１体項U1の関数形は次のように定義する。
【００１８】
【数２】
2008052308.tif 000003
【００１９】
f1が１体ポテンシャルの実体で、これは、原子iが、ni本の結合手を持っていたとして、その結合相手が、j1,j2,…である時に適用されるポテンシャルである。当然、この結合相手が変わるとf1も変わるが、その結合相手の組み合わせ全てについての重ね合わせでU1を定義している。この重ね合わせの係数を決めているのが動的結合次数項σα,jαである。これは、α番目の結合手が、その相手であるjαとの間の結合次数を表している。分母は、この重ね合わせの規格化因子である。
【００２０】
このように、あらゆる配位状態に応じた１体ポテンシャルf1を用意しておき、その重ね合わせで定義することで、ある配位状態から別の配位状態への遷移を滑らかに起こせる仕掛けになっている。この多重シグマを見ると、高速化が難しい印象を持つが、現在の配位状態に対応しない場合のΠσの項はゼロになりポテンシャルエネルギーに寄与しないため、有効な配位状態を効率よく列挙するプログラムを作成すれば、力の計算は原子数Ｎのオーダーで実現できる。
【００２１】
２体項引力項、３体項、基本的にU1と同じである。
【００２２】
【数３】
2008052308.tif 000004
【００２３】
【数４】
2008052308.tif 000005
【００２４】
３体項は１体項と同じく、中心の原子の配位状態の重ね合わせでポテンシャルエネルギーを表現する。２体項は、中央の結合を挟む両側の原子の配位状態に依存する。
【００２５】
動的結合次数項σは、注目しているダイナミックボンドが、相手の原子とどの程度軌道のオーバーラップを作っているかに応じて変化する関数である。
【発明を実施するための最良の形態】
【００２６】
（シミュレーション装置）
図１は、本発明のシミュレーション装置の機能的な構成を示すブロック図である。このシミュレーション装置は、図１に示すように、データ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０を少なくとも具備し、一連の多体分子シミュレーションの機能を有した処理制御部４１と、操作者からのデータや命令などの入力を受け付ける入力部４２と、シミュレーション結果を出力する出力部４５と、原子間相互作用ポテンシャルの定数、および原子核座標と各原子の内部座標などの入力データなどを格納したデータ記憶部４３と、シミュレーションプログラムなどを格納したプログラム記憶部４４とから少なくとも構成されている。
【００２７】
データ解析部１０は、構造評価部１１、統計情報計算部１２、構造可視化部１３等のような、原子座標およびその時系列データの解析機能を少なくとも備えている。例えば、構造評価部１１は、原子間相互作用計算部の原子座標情報管理部の情報から、動径分布関数等を計算する機能を有する。また、原子間相互作用計算部の動的結合次数評価部で得られる各原子間の結合次数情報から、結合欠損等の結晶欠陥解析を行う機能を有する。統計情報計算部１２は、原子間相互作用計算部の動的結合次数評価部で得られる各原子間の結合次数情報から、リング状結合ネットワークの統計解析を行う機能を有する。更に、シミュレーション実行部の古典分子動力学計算部で得られる各原子核座標の時系列情報から、拡散係数などの輸送係数を計算することができる。構造可視化部１３は、原子間相互作用計算部の原子核座標情報管理部に格納されている原子核座標および動的結合次数評価部で得られる各原子間の結合次数情報から、結合手を含む３次元コンピュータ・グラフィックスを生成し、画像ファイルやディスプレイに表示させる機能を有する。結合次数に応じて原子間結合手の太さを変えて表示させることができ、各原子間の結合の強さを直観的に把握できる。
【００２８】
シミュレーション実行部２０は、古典分子動力学計算部２１と、構造最適化計算部２２を少なくとも具備している。古典分子動力学計算部２１は、原子間相互作用計算部３０で計算される原子核座標および各原子の内部自由度に加わる力に基づいて、指定した微小時間における移動量を計算して座標情報を更新し、これを指定された回数繰り返して、原子核座標と内部自由度変数の時間発展を計算する。構造最適化計算部２２は、原子間相互作用計算部３０で計算される原子核座標および各原子の内部自由度に加わる力に基づいて、エネルギーを最小化する方向に移動させることにより、入力した初期座標配置に近い局所的な安定構造を探索する。
【００２９】
原子間相互作用計算部３０は、原子核座標情報管理部３１、内部座標情報管理部３２、配位数依存ポテンシャル定数情報管理部３３、動的結合次数評価部３４、配位状態判定部３５、ポテンシャルエネルギー計算部３６、１階微分計算部３７を少なくとも有する。原子核座標情報管理部３１は、原子核座標およびその時間微分の数値を保持する。内部座標情報管理部３２は、各原子の共有結合手の方向を決定する内部自由度変数およびその時間微分の数値を保持する。配位数依存ポテンシャル定数情報管理部３３は、配位状態に依存して使い分けられるポテンシャル定数の数値を分類し保持する。動的結合次数評価部３４は、原子核座標と内部座標の値から各原子間の結合次数を評価する機能を有する。配位状態判定部３５は、動的結合次数評価部３４で評価された結合次数から、各原子の配位状態を判定する。ポテンシャルエネルギー計算部３６は、原子核座標と内部座標の値、動的結合次数評価部３４で評価された結合次数、及び配位状態判定部３５で定められた配位状態から、系の全ポテンシャルエネルギーを計算する。１階微分計算部３７は、原子核座標と内部座標の値、さらに動的結合次数評価部３４で評価された結合次数から、原子核および各原子の内部自由度に加わる力を計算する。
【００３０】
図１において、入力部４２はキーボード、マウス、ライトペン、フレキシブルディスク装置、またはＬＡＮ端子などで構成される。データ記憶部４３およびプログラム記憶部４４は、周知の磁気テープ、磁気ドラム、磁気ディスク、光ディスク、光磁気ディスク、あるいはＲＯＭ、ＲＡＭなどの半導体メモリー等を用いた記録部である。また出力部４５は、ディスプレイ装置やプリンタ装置などにより構成されている。図１に示した処理制御部４１、データ記憶部４３、およびプログラム記憶部４４等はＣＰＵ，およびこのＣＰＵに接続されたＲＯＭ、ＲＡＭ、磁気ディスクなどの記憶装置を含む通常のコンピュータシステムで構成すればよい。図１に示した処理制御部４１で実行される各処理の入力データはデータ記憶部４３に格納され、プログラム命令はプログラム記憶部４４に記憶される。そしてこれらの入力データやプログラム命令は必要に応じてＣＰＵに読み込まれ、演算処理が実行されるとともに、データ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０で発生した数値情報などのデータはデータ記憶部４３に格納される。
【００３１】
次に、原子間相互作用計算部３０の動作について説明する。原子間相互作用計算部３０においては、まず動的結合次数評価部３４により各原子間の結合次数が評価される。その結果を用いて配位状態判定部３５において、各原子の配位状態が判定され、ポテンシャルエネルギー計算部３６において系の全ポテンシャルエネルギーが計算され、１階微分計算分３７において各原子核および内部自由度に加わる力が計算される。
【００３２】
各原子間の結合次数を判定する動的結合次数評価部３４の機能的な構成を示すブロック図を図２に示す。図２に示すように、動的結合次数評価部３４は、結合手の重なり判定部３４１、及びスイッチング関数部３４２を少なくとも有している。
【００３３】
図２において結合手の重なり判定部３4１は、原子核座標情報管理部３１及び内部座標情報管理部３２より座標情報を読取り、原子iのα番目の結合手と原子jとの間の結合手の重なり度を表すスカラ量σiα,jを、図２の式（＊１）乃至式（＊４）で計算する。式（＊１）乃至式（＊４）でｉ，ｊは各原子を示す指数（インデックス）で、αは原子ｉに属する結合手、βは原子ｊもしくはkに属する結合手を示す指数である。また式（＊４）で定義される関数ψの添え字Ｉは、原子iの種類を示す。ＣＩ，ｐＩ，ｎＩは、Iで指定される原子の種類に固有の定数である。ｓ，θ，φは、原子の内部自由度変数を示す。式（＊１）に見られるそれらの添え字は、各内部自由度変数が属する原子の指数と結合手の指数を示している。ｒi jは原子iと原子ｊの間の距離を示し、ｅi jは原子iの原子核位置を始点とし原子ｊの原子核位置を終点とするベクトルの単位ベクトルを示す。結合手の重なり度は、式（＊１）に示すように、原子iのα番目の結合手と原子ｊの間の各結合手の組み合わせについて式（＊３）で示される関数Ｓを計算し、その総和として定義される。もし相手の原子jとして、何れの原子とも結合していない未結合手の状態を表す記号Ｄが指定された場合、式（＊２）で示される関数f Dで未結合手状態の割合を計算する。式（＊２）及び（＊３）において、関数ｆｃはスイッチング関数部３４２で評価され、それは式（＊５）で定義される。近距離では１、遠距離では０となる連続関数である。式（＊５）Ｒ，Ｄは、スイッチングを起こす距離を規定するパラメータで、それらは式（＊２）及び（＊３）において、２原子の組み合わせ毎に指定できる。
【００３４】
以上の動的結合次数評価部３４における結合次数の決定方法を直観的に把握するための模式図を図３に示す。原子ｉのα番目の結合手と原子ｊ間のβ番目の結合手の間の関数Ｓは、原子ｉのα番目の結合手と原子ｊ間のβ番目の結合手の重なり度に応じて０から１まで連続的に増加する。この結合手の重なり度は、図３のＡ及びＢに示すように、原子ｉと原子ｊ間の距離と、それらの結合手の相対的な向きによって定まる。図３のＡ及びＢに描かれている結合手の形状は、式（＊４）で定義される関数ψを正確に反映している。即ち、式（＊４）で定義された関数ψを、引数である単位ベクトルｅの方位を変え、対応する関数値を単位ベクトルｅの向きをもつ動径の長さで表すと、動径の先端をプロットして描かれる曲面は、Ａ及びＢに示したような偏った亜鈴形を成す。もし、２原子が相互作用できる距離に近づいても、Ａに示したように、互いの結合手に重なりが生じていない場合には、式（＊１）乃至式（＊４）より、重なりは負の値となり、関数Ｓは０と判定される。一方、Ｂに示したように、互いの結合手に重なりが生じている場合には、式（＊１）乃至式（＊４）より、重なりは正の値を持つと評価され、関数Ｓは１と判定される。尚、図３には描かれていないが、関数Ｓは重なり度の単調増加関数ではなく、重なり度が一定値を超えると再び０に減少する式（＊３）で定義されている。これは、相手の原子を完全に覆い尽くすほど異常に結合手が伸びた場合、相手の原子の更に先に存在する原子とも同時に結合して不自然な安定化を起こす可能性があるため、これを防ぐための処置として、重なり度が一定値を超えると再びＳが０に減少するようになっている。
【００３５】
系の全ポテンシャルエネルギーを計算するポテンシャルエネルギー計算部３６の機能的な構成を示すブロック図を図４に示す。図４に示すように、ポテンシャルエネルギー計算部３６は、１体項計算部３６１、２体項計算部３６２、３体項計算部３６３及び全エネルギー合算部３６４を少なくとも有している。
【００３６】
図５に、１体項計算部３６１の機能を示す。１体項計算部３６１は、内部座標情報管理部３２より座標情報を読取り、式（＊７）に従って、内部座標変数のみに依存する項u1と、配位状態判定部３４で列挙された原子iがとる全ての配位状態について１体ポテンシャルエネルギーf1の平均した項を加えて１体項U1の値を定める。u1は、各結合手の第一自由度変数sの値から、式（＊８）で計算する。式（＊８）でｉは原子を示す指数で、αは原子ｉに属する結合手を示す指数である。nｉは、原子ｉに属する結合手の本数を示す。ａ，ｍは、原子iの種類に依存する定数である。f1は、指定された配位状態によってきまる１体ポテンシャルエネルギーで、配位数依存ポテンシャル情報管理部３３から値を読み取る。平均操作をする際、各配位状態におけるf1には重み付け係数として、指定された配位状態における各結合手の動的結合次数σの積を掛ける。
【００３７】
２体項計算部３６２の機能を図６に示す。２体項計算部３６２では、原子核座標情報管理部３１より座標情報を読取り、２原子間距離のみで決まる反発項U2repを計算する。また、配位状態判定部３５で列挙された原子iおよび原子ｊがとる全ての配位状態の組み合わせについて２体引力ポテンシャルエネルギーf2atrを求め、これら平均を式（＊１０）に従って計算し、これを引力項U2atrの値とする。平均操作をする際、各配位状態におけるf2atrには重み付け係数として、指定された配位状態における原子iおよびｊの全ての結合手の動的結合次数σの積を掛ける。f2repおよびf2atrは、式（＊１１）に示される２体ポテンシャルエネルギーｆ２typeで評価される。ｆ２の上付き添え字typeは、２体ポテンシャルエネルギーの種類を指定する標識で、本シミュレーション装置では、等方的に働く反発相互作用を表す添え字repと、結合が組まれた原子間のみに作用する引力相互作用を表す添え字atrの２種類が用意されている。２体ポテンシャルエネルギーｆ２typeは、式（＊１１）に示すように、式（＊１２）で定義される関数fsw2の２項の和で定義される。式（＊１２）で、ｒは２原子間距離を示し、A，ｐ，ａはパラメータを示す。これらのパラメータは２つの原子の種類Ｉ，Ｊと２体ポテンシャルの種類typeによって異なるの値を取る。f2atrとして使用される場合、A，ｐ，ａは更に配位状態にも依存して異なる値をとる。特にａはカットオフ距離と呼び、距離の次元を有し、引力相互作用が完全に打ち切られて０になる距離を指定する。
【００３８】
図７において３体項計算部３６３は、３原子間の位置関係で決まる３体ポテンシャルエネルギーU3を計算する。原子核座標情報管理部３１より座標情報を読取り、配位状態判定部３５で列挙された原子iがとる全ての配位状態の組み合わせについて３体ポテンシャルエネルギーf3を求め、これらの平均を式（＊１３）に従って計算し、これを３体項U3の値とする。平均操作をする際、各配位状態におけるf3には重み付け係数として、指定された配位状態における各結合手の動的結合次数σの積を掛ける。３体ポテンシャルエネルギーｆ3は、式（＊１４）に示すように、式（＊１５）で定義される３体関数Λと式（＊１６）で定義される3体関数Θの積で定義される。３体関数Λは、原子i−ｊ間の距離ｒｉｊと原子i−ｋ間の距離ｒｉｋの関数で、λ，γ，ａはパラメータである。各パラメータの下付き添え字は、３体項を構成する原子の種類を指定する。例えば、下付き添え字JIKは、Ｉで指定される原子の種類を頂点とし、それがJおよびＫで指定される原子と結合を組んで結合角を構成している場合を表す。また、γ，ａの上付き添え字は、上記結合角を張る２つの結合の一方を指定するための標識である。ａはカットオフ距離を表し、３体ポテンシャルエネルギーが作用する範囲を指定する。また、3体関数Θは、結合角の角度θjikの関数で、θ0，αはパラメータを示す。各パラメータの下付き添え字は、３体項を構成する原子の種類を指定する。特にパラメータθ0JIKは、Ｉで指定される原子の種類を頂点とし、それがJおよびＫで指定される原子と結合を組むことで構成される結合角の再安定角度を決定する。以上のパラメータは全て、f3の引数l1,l2,…l3で指定される配位状態に応じて異なる値をとる。
【００３９】
図４において、全エネルギー合算部３６４は、１体項計算部３６１、２体項計算部３６２、及び３体項計算部３６３で計算されたU１，U２rep，U２atr，U3を式（＊６）のように合計し、系の全ポテンシャルエネルギーΦを計算する。式（＊６）において、ｉ，ｊ，ｋは各原子を示す指数で、Ｎは系に含まれる全粒子数を表す。
【００４０】
各原子核および内部自由度に加わる力を計算する、１階微分計算部３７の機能的な構成を示すブロック図を図８に示す。図８に示すように、１階微分計算部３７は、原子核座標に関する微分実行部３７１、及び内部座標に関する微分実行部３７２を少なくとも有している。
【００４１】
図８において原子核座標に関する微分実行部３７１は、動的結合次数評価部３４及びポテンシャルエネルギー計算部３６で得られる計算結果を利用し、式（＊１７）より全ポテンシャルエネルギーの各原子核位置に関する１階偏微分係数を計算する。式（＊１７）において、ｉ，ｊ，ｋは、各原子を示す指数である。
【００４２】
図８において内部座標に関する微分実行部３７２は、動的結合次数評価部３４及びポテンシャルエネルギー計算部３６で得られる計算結果を利用し、式（＊１８）より全ポテンシャルエネルギーの各内部座標変数に関する１階偏微分係数を計算する。式（＊１８）において、ξiαは、式（＊１）で示した、原子ｉに属する結合手αの、３つの内部自由度変数siα，θiα，φiαをまとめたベクトル変数である。siα，θiα，φiαはこの順番で、ξiαの第一要素、第二要素、第三要素に対応する。
【００４３】
１階微分計算部３７で原子核座標及び内部座標に関するポテンシャルの１階微分が計算されると、それらを質量で除算することで、直ちに原子核及び内部座標変数に加わる力を求めることができる。ただし、ここで言う質量とは、原子核に加わる力を求める際はkgの次元を有する原子核の質量そのものであるが、内部座標変数に加わる力を求める際は、内部座標変数用に別途指定される、エネルギーの次元を有する仮想的な質量を用いる。質量は一般に、運動する自由度の変化の容易さを決めるパラメータであり、その次元はkgに限らない。本シミュレーション装置において内部座標変数を原子核に素早く追従させたい場合には、原子核の質量に比べて十分小さな質量を指定する。ただし、内部座標変数に対してあまり小さな質量を指定すると、計算の精度低下を招くため、扱う系に応じて適切な質量の値を使用者が指定することになる。
【００４４】
以上のように、原子間相互作用計算部３０によれば、原子核同様、各原子の共有結合手の方向を決定する内部自由度変数にも力が働くため、シミュレーション実行部２０によって、その力に基づいて結合手の方向を変化させることができる。これにより、各原子間の結合次数をシミュレーションの最中に連続的に変化させることができ、共有結合の組み換えを伴う分子シミュレーションを行うことが可能となる。このとき、結合次数に応じて配位状態判定部３５で判定される各原子の配位状態も変化するため、配位状態に依存してポテンシャルエネルギーを連続的に変化させることが可能となる。この拡張に伴って増える計算は、内部座標情報管理部３２及び配位数依存ポテンシャル定数情報管理部３３のメモリ操作、動的結合次数評価部３４、配位上体判定部３４に止まるため、従来の古典分子動力学計算と同様、数千個以上の原子から構成される大規模系のシミュレーションが、短時間で効率的に実行できる。
【００４５】
また、原子間相互作用計算部３０によれば、２体引力相互作用および３体相互作用の発現を結合が組まれている原子間に限定することができる。このため、企図しない原子間に余計な引力やポテンシャル障壁が形成される事を防ぐことができる。更に、結合を組んでいない原子間の２体引力相互作用を打ち消すように３体相互作用項を調節する必要がなくなるため、２体ポテンシャル定数と３体ポテンシャル定数を独立に決定することができ、多元素混在系用のポテンシャル定数を用意することが可能となる。
【００４６】
更に、原子間相互作用計算部３０によれば、各原子の配位状態、即ち、全ての結合手の結合相手の組み合わせに応じて異なるポテンシャル定数を割り当てることができる。これにより、例えば同じSi-O結合でも、そのSi原子やO原子がもつ他の結合手の結合相手の組み合わせに応じて、そのSi-O結合に異なる結合エネルギーを割り当てることができるようになる。この方法であれば、１重結合、２重結合、３重結合の違いまで表現できる。また、結合手の数を予め標準的な結合手の数より多く用意しておけば、超原子価状態も表現できるようになるため、新しい結合を徐々に形成しつつ、古い結合手が徐々に消滅していく、というリアルな共有結合組み換え反応経路を古典的分子動力学計算の枠組みで再現できるようになる。
【００４７】
（シミュレーション方法）
図１に示した原子間相互作用計算部３０を用いて分子動力学シミュレーションを実行するには、所望の統計アンサンブルを生成するように設計されたオイラー−ラグランジュ方程式の数値積分を行えばよい。また、構造最適化計算を実行するには、原子間相互作用計算部３０で計算される力の情報を用いて、共役勾配法等、多次元空間上のポテンシャルの極小点探索アルゴリズムを実行すればよい。以下において、分子動力学計算及び構造最適化計算を行うための具体例、即ち、シミュレーション方法１及びシミュレーション方法２を示す。
【００４８】
シミュレーン方法１：図９のフローチャートを参照して、本発明のシミュレーション方法１を説明する。
（イ）まず、ステップＳ１１０で、分子動力学計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標、配位状態依存ポテンシャル定数を入力部４２から読み込む。ここで分子動力学計算を実行する上で必要となる制御データとは、数値積分を実行する際の1ステップ当たりの時間刻み、全ステップ数など、分子動力学計算を実行する上で必要となるデータである。また温度制御を使用するかどうかの選択や、その際の設定温度の数値なども含む。分子動力学計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標の情報は基本的に全て入力部４２より読み込むが、入力が省略された情報については、必要に応じて自動的に補完する。例えば、内部座標情報が省略された場合、原子核座標情報に基づいて、各原子の最も近い原子から順に結合を組むように、内部自由度の座標を自動的に決定する。
（ロ）次にステップＳ１２０で、原子核座標及び内部座標の情報から、各原子間の動的結合次数を評価する。
（ハ）次にステップＳ１３０で、各原子の配位状態を全て列挙する。
（ニ）次にステップＳ１４０で、上記配位状態に属する２体ポテンシャルおよび３体ポテンシャルの定数が完備されているか判定する。完備されていない時はステップＳ１５０に進み、不足しているポテンシャル定数を要求するエラーメッセージを出力し、計算の途中データを記録して終了する。必要なポテンシャル定数が揃っている時はステップＳ１６０に進む。
（ホ）ステップＳ１６０ではポテンシャルエネルギー、１階微分を計算する。
（ヘ）ステップＳ１７０では原子座標、内部座標に関する１階微分の値から力を計算し、これをオイラー−ラグランジュ方程式に代入して各座標の加速度を計算し、分子動力学計算を１ステップ進める。
（ト）ステップＳ１８０では入力データで指定された全ステップ数に達したか判定する。達していない場合はＳ１２０に戻り、達している場合には計算結果を出力して終了する。
【００４９】
シミュレーション方法２：図１０のフローチャートを参照して、本発明のシミュレーション方法２を説明する。
（イ）まず、ステップＳ２１０で、構造最適化計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標、配位状態依存ポテンシャル定数を入力部４２から読み込む。ここで構造最適化計算を実行する上で必要となる制御データとは、構造最適化計算の収束判定条件など、構造最適化計算を実行する上で必要となるデータである。構造最適化計算を実行する上で必要となる制御データ、初期原子核座標及び初期内部座標の情報は基本的に全て入力部４２より読み込むが、入力が省略された情報については、必要に応じて自動的に補完する。例えば、内部座標情報が省略された場合、原子核座標情報に基づいて、各原子の最も近い原子から順に結合を組むように、内部自由度の座標を自動的に決定する。
（ロ）次にステップＳ２２０で、原子核座標及び内部座標の情報から、各原子間の動的結合次数を評価し、原子間結合のトポロジーを定める。
（ハ）次にステップＳ２３０で、各原子の配位状態を全て列挙する。
（ニ）次にステップＳ１４０で、上記配位状態に属する２体ポテンシャルおよび３体ポテンシャルの定数が完備されているか判定する。完備されていない時はステップＳ２５０に進み、不足しているポテンシャル定数を要求するエラーメッセージを出力し、計算の途中データを記録して終了する。必要なポテンシャル定数が揃っている時はステップＳ２６０に進む。
（ホ）ステップＳ２６０ではポテンシャルエネルギー、１階微分を計算する。
（ヘ）ステップＳ２７０では、次の１次元極小点探索の方向を決定する。ここで言う方向とは、原子核座標および各原子の内部自由度座標を併せた全自由度から構成される多次元空間内における方向ベクトルを指す。最初にこのステップに来た場合には１階微分から得られる力の方向を指定し、２回目以降にこのステップに来た場合は、既に探索した方向に共役な方向を、共役勾配法の手続きで計算する。
（ト）ステップＳ２８０では、前のステップで指定された、原子核座標および各原子の内部自由度座標を併せた全自由度から構成される多次元空間内における方向ベクトル上で、１次元の極小点探索を実行する。
（チ）ステップＳ２９０では、構造最適化計算の収束判定条件が満たされるか判定する。収束判定条件とは、具体的には、ステップＳ２８０の前後で、原子核座標および内部自由度座標の移動幅を計算し、その合計が入力データで指定した範囲に収まっていることを指す。収束判定条件が満たされない時は、ステップＳ２２０に戻る。収束判定条件が満たされた時は、計算結果を出力して終了する。
【００５０】
以上のシミュレーション方法１及びシミュレーション方法２によれば、共有結合の所在を明示しながら、その滑らかな繋ぎ変えを可能にする分子動力学計算及び構造最適化計算を行うことができる。共有結合の発現を明示された原子間に限定できるため、多元素が混在する複雑な系の、配位状態に応じて変化できる高度なポテンシャルを用意することが可能となる。つまり、シミュレーション方法１及びシミュレーション方法２によれば、多元素混在系において、共有結合の組み換えを伴う化学変化の分子シミュレーションを、各原子の配位状態に応じてポテンシャルエネルギーを調節しながら、数千原子以上の大きな規模で実行できる。
【００５１】
（シミュレーションプログラムを格納した記録媒体）
図９及び図１０に示すシミュレーション方法を含め、図１に示したデータ解析部１０、シミュレーション実行部２０、原子間相互作用計算部３０におけるシミュレーション方法を実現するためのプログラムはコンピュータ読取り可能な記録媒体に保存しておいても良い。この記録媒体をコンピュータシステムによって読み込ませ、図１に示したプログラム記憶部４４に格納し、このプログラムを処理制御部４１で実行して本シミュレーション方法を実現することもできる。ここで、記録媒体とは、例えば磁気ディスク、光ディスク、光磁気ディスク、磁気テープなどのプログラムを記録することができるような媒体などが含まれる。コンピュータの外部メモリも、ここで言う記録媒体に含まれる。
【００５２】
本発明のように様々な材料に応用できる本分子動力学法は、今後の分子シミュレーションに新たな潮流を引き起こす程の波及効果が期待される。パラメータ決定が容易になるため，酸化物材料，生体分子など，今後のナノテクノロジー，バイオテクノロジーに必要な様々な系の大規模分子シミュレーションが次々に可能になると思われる。Car-Parrinelloの第一原理分子動力学計算は，構成元素を選ばぬ汎用性と信頼性の高さから大変広く利用されているが，計算量の多さから扱える原子数は通常数十原子程度に限られ，スーパーコンピュータを独占使用しても千原子規模の分子動力学計算すら難しい。これに対し本発明を持ちいれれば，通常普及しているPCで数万原子規模の大規模計算が可能であり，64bitマルチプロセッサ機を導入すれば数十万〜百万原子規模も扱うことができる。自然に結合組替えを起こすことが可能となるため，これまでコンフォーメーション変化の再現が中心だった蛋白質分子のシミュレーションで，加水分解・脱水縮合などの化学反応まで表現することが可能となる。DNAやペプチド結合の組替え，ADPからATPが合成される反応など，重要な生体化学反応の多くが加水分解・脱水縮合反応であり，これらを古典的に扱えるようになれば，利用者は爆発的に増えると予想される。このように本発明は新たな科学技術計算ソフトフェアの市場を開拓する潜在力を秘めている。
【図面の簡単な説明】
【００５３】
【図１】本発明の実施の形態に係るシミュレーション装置の機能的な構成を示すブロック図である。
【図２】図１に示した動的結合次数部３４の機能的な構成を示すブロック図である。
【図３】動的結合次数評価部３４における結合次数の決定方法を模式的に示した図である。
【図４】図１に示したポテンシャルエネルギー計算部３６の機能的な構成を示すブロック図である。
【図５】図４に示した１体項計算部３６１の機能的な構成を示す図である。
【図６】図４に示した２体項計算部３６２の機能的な構成を示す図である。
【図７】図４に示した３体項計算部３６３の機能的な構成を示す図である。
【図８】図１に示した１階微分計算部３７の機能的な構成を示すブロック図である。
【図９】本発明のシミュレーション方法１を示すフローチャートである。
【図１０】本発明のシミュレーション方法２を示すフローチャートである。
【符号の説明】
【００５４】
１０  データ解析部
１１  構造評価部
１２  統計情報計算部
１３  構造可視化部
２０  シミュレーション実行部
２１  古典分子動力学計算部
２２  構造最適化計算部
３０  原子間相互作用計算部
３１  原子核座標情報管理部
３２  内部座標情報管理部
３３  配位数依存ポテンシャル定数情報管理部
３４  動的結合次数評価部
３５  配位状態判定部
３６  ポテンシャルエネルギー計算部
３７  １階微分計算部
４１  処理制御部
４２  入力部
４３  データ記憶部
４４  プログラム記憶部
４５  出力部
３４１  結合手の重なり判定部
３４２  スイッチング関数部
３６１  １体項計算部
３６２  ２体項計算部
３６３  ３体項計算部
３６４  全エネルギー合算部
３７１  原子核座標に関する微分実行部
３７２  内部座標に関する微分実行部
【図１】
2008052308.tif 000006
【図２】
2008052308.tif 000007
【図３】
2008052308.tif 000008
【図４】
2008052308.tif 000009
【図５】
2008052308.tif 000010
【図６】
2008052308.tif 000011
【図７】
2008052308.tif 000012
【図８】
2008052308.tif 000013
【図９】
2008052308.tif 000014
【図１０】
2008052308.tif 000015