(57)【要約】【課題】触媒反応などの進行に伴って電子状態が変化する場合でも、古典力学の範囲で高速に触媒反応をシミュレーションする。【解決手段】本シミュレーション方法及びプログラムは、例えば触媒などの結晶表面上のガス反応を分子動力学法で解析する際に好適に利用される。本方法では、分子動力学モデルを構成する各原子の局所状態が判定され、局所状態に基づいて、前記各原子に関するポテンシャル関数が変更される。そして、変更されたポテンシャル関数を用いて、分子動力学法を用いたアルゴリズムにより各原子の状態が判定され、各原子の状態に基づいて分子動力学モデルの状態が判定される。よって、吸着・脱離や化学反応などの前後で電子状態が変化する場合でも、その局所状態に応じてポテンシャル関数が変更され、分子動力学法が適用される。よって、状態変化を伴う反応を高速にシミュレーションすることができる。【選択図】図１
2007316976.tif 000002
【特許請求の範囲】
【請求項１】
分子動力学法を用いたシミュレーション方法であって、
分子動力学モデルの初期状態を作成するステップと、
前記分子動力学モデルを構成する各原子の局所状態を判定するステップと、
前記局所状態に基づいて、前記各原子に関するポテンシャル関数を変更するステップと、
変更されたポテンシャル関数を用いて、分子動力学法を用いたアルゴリズムにより、前記各原子の状態を判定するステップと、
前記各原子の状態に基づいて、前記分子動力学モデルの状態を判定するステップと、を有することを特徴とするシミュレーション方法。
【請求項２】
分子動力学法を用いたシミュレーションプログラムであって、コンピュータ上で実行されることにより、
分子動力学モデルの初期状態を作成するステップと、
前記分子動力学モデルを構成する各原子の局所状態を判定するステップと、
前記局所状態に基づいて、前記各原子に関するポテンシャル関数を変更するステップと、
変更されたポテンシャル関数を用いて、分子動力学法を用いたアルゴリズムにより、前記各原子の状態を判定するステップと、
前記各原子の状態に基づいて、前記分子動力学モデルの状態を判定するステップと、を実行することを特徴とするシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学法を用いて、系を構成する原子の状態を解析するシミュレーション手法に関する。
【背景技術】
【０００２】
触媒などの結晶表面上のガス反応を分子動力学法（ＭＤ法：Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ  ｍｅｔｈｏｄ）により解析することが行われている。触媒表面上のガス反応では、化学吸着・反応・脱離など電子状態変化を伴う素過程が最重要である。触媒反応では吸着・脱離・化学反応などの前後で電子状態が変化するため、原子に働く力を導出するポテンシャル関数を、反応の全行程を通じて同一の関数で記述することができない。即ち、反応全過程を通して原子に働く力を導出できるポテンシャル関数は存在しないため、原子に働く力を古典的分子動力学でシミュレーションすることは不可能であった。
【０００３】
一方、第１原理分子動力学法では強力な計算機能力と膨大な計算時間を要するため、実験との対応可能なモデルスケールと時間スケールで色々な熱力学条件下の触媒表面反応シミュレーションを行うことは不可能である。
【０００４】
なお、物質に対する分子の吸着・吸収過程を、モンテカルロ法の分子シミュレーションにより解析する手法が特許文献１に記載されている。
【０００５】
【特許文献１】特開平１０−２７１８７号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
本発明は、上記のような課題を解決するためになされたものであり、反応などの進行に伴って電子状態が変化する場合でも、古典力学の範囲で高速に反応をシミュレーションすることを目的とする。
【課題を解決するための手段】
【０００７】
本発明の１つの観点では、分子動力学法を用いたシミュレーション方法は、分子動力学モデルの初期状態を作成するステップと、前記分子動力学モデルを構成する各原子の局所状態を判定するステップと、前記局所状態に基づいて、前記各原子に関するポテンシャル関数を変更するステップと、変更されたポテンシャル関数を用いて、分子動力学法を用いたアルゴリズムにより、前記各原子の状態を判定するステップと、前記各原子の状態に基づいて、前記分子動力学モデルの状態を判定するステップと、を有する。
【０００８】
本発明の他の観点では、分子動力学法を用いたシミュレーションプログラムは、コンピュータ上で実行されることにより、分子動力学モデルの初期状態を作成するステップと、前記分子動力学モデルを構成する各原子の局所状態を判定するステップと、前記局所状態に基づいて、前記各原子に関するポテンシャル関数を変更するステップと、変更されたポテンシャル関数を用いて、分子動力学法を用いたアルゴリズムにより、前記各原子の状態を判定するステップと、前記各原子の状態に基づいて、前記分子動力学モデルの状態を判定するステップと、を実行する。
【０００９】
上記のシミュレーション方法及びプログラムは、例えば触媒などの結晶表面上のガス反応を分子動力学法で解析する際に好適に利用される。本方法では、分子動力学モデルを構成する各原子の局所状態が判定され、局所状態に基づいて、前記各原子に関するポテンシャル関数が変更される。そして、変更されたポテンシャル関数を用いて、分子動力学法を用いたアルゴリズムにより各原子の状態が判定され、各原子の状態に基づいて分子動力学モデルの状態が判定される。よって、吸着・脱離や化学反応などの前後で電子状態が変化する場合でも、その局所状態に応じてポテンシャル関数が変更され、分子動力学法が適用される。よって、状態変化を伴う反応を高速にシミュレーションすることができる。
【発明を実施するための最良の形態】
【００１０】
以下、図面を参照して本発明の好適な実施の形態について説明する。
【００１１】
図１は、本発明による表面反応シミュレーションプログラムの概略フローチャートである。図１に示すように、この表面反応シミュレーションプログラムでは、まず、解析の対象となる体系のモデル構成、熱力学的状態など（以下、「ＭＤモデル」と呼ぶ。）が入力、設定される。また、ＭＤシミュレーションのプログラムパラメータなどを含むシミュレーション条件が入力、設定される（ステップＳ１）。
【００１２】
次に、設定された条件に基づき、解析の対象となる体系の初期状態が設定される。ここで、体系の初期状態とは、体系を構成する複数の原子各々の初期位置、初期速度などを含む。具体的には、今回のシミュレーションが新規計算であるか否かが判定され（ステップＳ２）、新規計算であれば、ステップＳ１で設定されたＭＤモデルの初期状態が作成される（ステップＳ３）。一方、新規計算でなければ、前回のシミュレーションにより得られた最終結果を読み込み、これを今回のシミュレーションにおける初期状態として使用する（ステップＳ４）。
【００１３】
次に、境界条件を考慮して原子位置を修正する（ステップＳ５）。次に、各原子に働くエネルギー・力を計算する（ステップＳ６）。詳しくは、まず、ＭＤモデルを構成する各原子について、近傍の局所状況（原子配置、触媒表面からの距離）やその保持時間などを用いて状態判定を行う。次に、状態判定結果に基づいて各原子の状態に応じた状態フラグを設定し、この状態フラグに基づいて各原子についてポテンシャル関数を選定する。そして、そのポテンシャル関数を用いて、当該原子に働くエネルギー・力を計算する。
【００１４】
次に、各原子について得られたエネルギー・力に基づいて、ＭＤアルゴリズムを実行して各原子の位置、速度の変化を計算し、各原子の位置、速度を更新する（ステップＳ７）。次に、ＭＤアルゴリズムによる計算結果に基づいて、体系の熱力学的状態を計算し、ステップＳ１で設定したＭＤモデルの状態と比較し、状態の修正が必要か否かを判定する（ステップＳ８）。修正が必要であれば、構成原子数、温度、圧力などを修正する（ステップＳ９）。こうして得られた計算結果は、必要に応じてファイルに記録され、解析され、表示装置などに可視化された状態で出力される（ステップＳ１０）。
【００１５】
以上のうち、ステップＳ５〜Ｓ７がメインのＭＤステップを構成し、ステップＳ５〜Ｓ１０が１つのＭＤサイクルを構成する。
【実施例】
【００１６】
次に、本発明の表面反応シミュレーションプログラムの好適な実施例について説明する。この実施例は、本プログラムにより、白金（Ｐｔ）表面上での一酸化炭素（ＣＯ）の酸化反応を例として、分子動力学（ＭＤ）シミュレーションを行う例である。
【００１７】
まず、本実施例のプログラムの概略について説明する。図２は、本実施例による計算の流れの概要を示すフローチャートである。
【００１８】
本実施例のプログラムでは、まず、解析の対象となる体系に対応するＭＤモデルの構成及びシミュレーション条件が設定される（ステップＳ１１）。なお、ＭＤモデルの構成は、具体的には、結晶部分の形状、ガス部分の形状・組成、温度・圧力条件などを含む。
【００１９】
次に、ＭＤモデルの初期状態が作成される（ステップＳ１２）。詳しくは、今回のシミュレーションが新規計算の場合、結晶部及びガス部に対し、シミュレーション条件に適合した位置及び速度の原子が生成される（ステップＳ１３）。一方、今回のシミュレーションが新規計算でない場合、即ち、以前のシミュレーション結果が存在する場合には、前回のシミュレーションの最終状態を読込み、初期状態とする（ステップＳ１４）。
【００２０】
次に、境界条件を考慮し（ステップＳ１５）、原子の局所状態判定を行う（ステップＳ１５）。この原子の局所状態を判定するステップは本実施例のプログラムにおける特徴部分である。具体的には、反応、吸着など、電子状態の変化を伴う素過程を原子の局所状態から判定して状態フラグを変更し、各原子に対するポテンシャル関数を選定する。なお、ポテンシャル関数は、各種の状態フラグを含む基本の関数式として用意されている。よって、各原子の局所状態に応じて状態フラグを設定、変更することにより、その局所状態に適合したポテンシャル関数を得ることができる。即ち、反応、吸着などの素過程において電子状態が変化した場合でも、変化後の原子の局所状態に適合したポテンシャル関数を得ることができる。
【００２１】
こうして、各原子の局所状態に適合したポテンシャル関数が得られると、当該ポテンシャル関数を用いて、各原子のエネルギー、各原子が受ける力を計算する（ステップＳ１７）。次に、ＭＤアルゴリズムにより、原子の速度、位置を更新する（ステップＳ１８）。
【００２２】
次に、ＭＤモデルが示す体系の状態判定を行う（ステップＳ１９）。詳しくは、吸着、反応によるガス相の組成変化が生じるため、ガス相中の原子数をモニタし、原子数がステップＳ１１で設定された設定条件から外れた場合に、原料ガス原子の生成、反応ガス分子の消滅を行う。また、ＭＤモデル中の温度・圧力をモニタし、設定条件からのずれを補正する（ステップＳ１９）。最後に、こうして得られたモデル状態データをファイルに記録するとともに、必要に応じて解析、可視化する（ステップＳ２０）。
【００２３】
次に、本実施例のプログラムの主要部分について詳細に説明する。
【００２４】
図３は、新規計算である場合のＭＤモデルの初期状態作成処理、即ちステップＳ１２及びＳ１３に相当する処理のフローチャートである。具体的には、まず、ＭＤセル形状が設定される（ステップＳ１２−１）。なお、ＭＤセルとは、解析の対象となる体系に対応するＭＤモデルを構成する単位を言う。この実施例では、全粒子数は、白金粒子数、酸素原子数及び一酸化炭素原子数の和として設定される。
【００２５】
次に、結晶部の初期状態が設定される（ステップＳ１２−２）。具体的には、１〜Ｎptの粒子番号が付与されたＮpt個のＰｔ粒子（原子）について、元素フラグ、状態フラグが付与され、粒子座標が設定される。
【００２６】
次に、ガス部の初期状態が設定される（ステップＳ１２−３）。具体的には、全粒子数のうち、Ｐｔ粒子の数を除いた数（即ち、酸素及び一酸化炭素原子の数）に相当する粒子に対して、粒子番号、元素フラグ及び状態フラグが付与される。なお、ガス部は、結晶部との間に厚みｄの空間をあけて設定される。
【００２７】
次に、メインＭＤステップ部分、即ち図２におけるステップＳ１５〜Ｓ１８について詳しく説明する。図４は、メインＭＤステップ部分のフローチャートである。図示のように、メインＭＤステップ部分は、境界条件を考慮して粒子位置を修正するステップＳ１５−１と、２体力の場合に粒子のエネルギー・力を計算するステップＳ１７−１と、３体力の場合に粒子のエネルギー・力を計算するステップＳ１７−２と、ＭＤアルゴリズムを実行するステップＳ１８−１とを有する。境界条件を考慮するステップＳ１５−１は全粒子について実行される。
【００２８】
図５は、２体力の場合のエネルギー・力の計算ステップの詳細を示す。なお、２体力の場合とは、２つの粒子間の関係に基づいてエネルギー・力を計算する場合を言う。まず、１つの粒子（「ｉ−粒子」とする。）を特定し、その粒子に付与されている元素フラグを取得する（ステップＳ１７−１−１）。次に、他の１つの粒子（「ｊ−粒子」とする。）
を特定し、両粒子間の状態判定を行う（ステップＳ１７−１−２）。具体的には、両粒子間の距離を所定の吸着距離及び反応距離と比較し、その結果に基づいて両粒子の状態フラグを変更し、その状態フラグに基づいてポテンシャル関数を選定する。そして、選定されたポテンシャル関数を用いて、ｉ−粒子に働く力を計算する。この処理を、図５に示すように、全粒子について、即ち、ｉ−粒子及びｊ−粒子が全粒子数に対応する粒子番号に至るまで繰り返す。こうして、各粒子について２体力のエネルギー・力の計算が終了する。
【００２９】
図６は、３体力の場合のエネルギー・力の計算ステップの詳細を示す。なお、２体力の場合とは、３つの粒子間の関係に基づいてエネルギー・力を計算する場合を言う。まず、１つの粒子（「ｉ−粒子」とする。）と、他の１つの粒子（「ｊ−粒子」とする。）を選択し、それらの距離に基づいて状態判定を行う（ステップＳ１７−２−１）。そして、状態判定の結果（例えば、吸着距離にある、反応距離にあるなど）に応じて両粒子の状態フラグを変更し、適当なポテンシャル関数を選定する（ステップＳ１７−２−２、Ｓ１７−２−４）。そして、選定されたポテンシャル関数を用いて、ポテンシャルエネルギー、即ちｉ−粒子に働く力が計算される（ステップＳ１７−２−３、Ｓ１７−２−５）。
【００３０】
次に、ｉ−粒子とｊ−粒子の状態に基づいて３体力があるか否かが判定され、３体力の場合の処理が実行される（ステップＳ１７−２−６）。具体的には、第３の粒子（「ｋ−粒子」とする。）が選択され、３つの粒子間の状態、即ち、吸着状態にあるか、反応状態にあるかが判定される（ステップＳ１７−２−７）。そして、状態判定の結果（例えば、吸着状態にある、反応状態にあるなど）に応じて３つの粒子の状態フラグを変更し、適当なポテンシャル関数を選定する（ステップＳ１７−２−８、Ｓ１７−２−１０）。そして、選定されたポテンシャル関数を用いて、ポテンシャルエネルギー、即ちｉ−粒子に働く力が計算される（ステップＳ１７−２−９、Ｓ１７−２−１１）。
【００３１】
図６に示すように、以上の処理を全粒子について、即ち、ｉ−粒子、ｊ−粒子及びｋ−粒子が全粒子数に対応する粒子番号に至るまで繰り返すことにより、各粒子についてのエネルギー・力が計算される。
【００３２】
図７は、図２のステップＳ１８に示すＭＤアルゴリズムの実行ステップの詳細を示す。具体的には、ベルレアルゴリズムなどのＭＤアルゴリズムを利用し、ステップＳ１７で計算された各粒子のエネルギー・力に基づいて、各粒子の速度、位置を更新する（ステップＳ１８−１）。そして、更新後の各粒子の速度、位置をステップＳ１１で設定したシミュレーション条件と比較し、粒子の生成・消滅、温度・圧力の変更などを行う（ステップＳ１８−２）。
【００３３】
図８は、本実施例のプログラムの特徴部分である状態判定の方法の一例を示す。触媒表面上でのガス反応の中で電子状態変化を伴う素過程としては図示のようにガス吸着、ガス脱離、化学反応（本例では一酸化炭素ＣＯの酸化反応）がある。
【００３４】
状態判定においては、運動を考える原子の局所的状況により場合分けし、例えば図９に示すように状態フラグや、ポテンシャル関数を切り替えることにより、電子状態変化を伴う素過程を古典力学の範囲で取り扱えるようにする。なお、ポテンシャル関数の選定やポテンシャルパラメータの設定は、分子軌道計算（ＭＯ）やバンド計算などにより行われる。
【００３５】
図９の例では、酸素原子（Ｏ）及び一酸化炭素原子（ＣＯ）のガス成分については、ガス成分の位置がガス部と結晶部との間の空間外にある場合は気相状態として、状態フラグＫ(i)＝１とする。ガス成分の位置がガス部と結晶部の空間内にある場合及びガス粒子と白金粒子の距離が空間距離ｄad以下である場合は吸着状態であるとして、状態フラグＫ(i)＝２とする。ガス粒子と白金粒子の距離が反応距離ｒacより小さい場合は反応状態として、状態フラグＫ(i)＝３とする。
【００３６】
また、白金の触媒表面については、白金粒子とガス粒子との距離が空間距離ｄadより大きい場合は未吸着状態として、状態フラグＫ(k)＝１とする。白金粒子とガス粒子との距離が空間距離ｄad以下である場合は吸着状態として、状態フラグＫ(k)＝２とする。また、白金粒子の周囲の状態フラグＫ(l)＝２のときは、状態フラグＫ(k)＝３とする。
【００３７】
このように、本発明では、触媒反応の進行に伴う状態変化に応じて各粒子の状態フラグを変更し、ポテンシャル関数を選定するので、電子状態の変化を伴う反応を高速にシミュレーションすることが可能となる。
【図面の簡単な説明】
【００３８】
【図１】本発明の表面反応シミュレーションプログラムのフローチャートである。
【図２】本発明の実施例による計算の流れの概要を示すフローチャートである。
【図３】ＭＤモデルの初期状態作成ステップのフローチャートである。
【図４】メインＭＤステップのフローチャートである。
【図５】２体力の場合のエネルギー・力の計算ステップのフローチャートである。
【図６】３体力の場合のエネルギー・力の計算ステップのフローチャートである。
【図７】ＭＤアルゴリズム実行ステップのフローチャートである。
【図８】本プログラムによる状態判定方法の例を示す。
【図９】状態フラグの変更方法の例を示す。
【図１】
2007316976.tif 000003
【図２】
2007316976.tif 000004
【図３】
2007316976.tif 000005
【図４】
2007316976.tif 000006
【図５】
2007316976.tif 000007
【図６】
2007316976.tif 000008
【図７】
2007316976.tif 000009
【図８】
2007316976.tif 000010
【図９】
2007316976.tif 000011