(57)【特許請求の範囲】
【請求項１】
多数の粒子を含む系について各粒子の運動方程式を数値的に解くことによりその系を解析する解析装置であって、前記解析装置は、入力部と演算部と記憶部と出力部とを備え、前記演算部は、モデル生成部とエネルギ調整部と分子動力学計算部とを備え、
前記モデル生成部は、前記入力部から入力されるユーザの入力情報に基づいて又は前記記憶部から読み込むことにより、粒子ごとに単一の又は複数のポテンシャルエネルギ関数が設定された前記多数の粒子を含む系の解析モデルを生成し、
前記複数のポテンシャルエネルギ関数は粒子間の複数の結合状態にそれぞれ対応し、
前記エネルギ調整部は、分子動力学計算中に粒子間の結合状態に応じて前記複数のポテンシャルエネルギ関数のいずれかを選択することにより、その結合に関与する粒子のポテンシャルエネルギ関数を変更し、
前記エネルギ調整部は、ポテンシャルエネルギ関数の変更によるポテンシャルエネルギの変化量をその結合に関与する粒子の運動エネルギに振り替えることにより、変更前後でエネルギ保存則が満たされるよう前記系のもつエネルギを調整し、
前記分子動力学計算部は、調整されたポテンシャルエネルギ及び運動エネルギから導かれる運動方程式を数値的に解くことにより各粒子の位置及び速度の少なくとも一方を計算し、
前記出力部は、計算された各粒子の位置及び速度の少なくとも一方を、又はそれに基づき前記演算部により演算される前記多数の粒子を含む系の物理量を、出力することを特徴とする解析装置。
【請求項２】
前記エネルギ調整部は、前記分子動力学計算部の演算結果に基づき更新された粒子間距離に基づいて結合状態を判別することを特徴とする請求項１に記載の解析装置。
【請求項３】
前記多数の粒子を含む系は、粒子間の複数の結合状態に各々が対応する複数のポテンシャルエネルギ関数が定義されている第１の粒子群と、単一のポテンシャルエネルギ関数が定義されている第２の粒子群と、を含み、
前記エネルギ調整部は、前記第１の粒子群について粒子間の結合状態に応じてポテンシャルエネルギ関数を選択し、
前記分子動力学計算部は、前記第１及び第２の粒子群について粒子の位置及び速度の少なくとも一方を求めることを特徴とする請求項１または２に記載の解析装置。
【請求項４】
多数の粒子を含む系について各粒子の運動方程式を数値的に解くことによりその系を解析する解析方法であって、本方法は、入力部と演算部と記憶部と出力部とを備える解析装置において実行され、前記演算部は、モデル生成部とエネルギ調整部と分子動力学計算部とを備え、本方法は、
前記モデル生成部が、前記入力部から入力されるユーザの入力情報に基づいて又は前記記憶部から読み込むことにより、粒子ごとに単一の又は複数のポテンシャルエネルギ関数が設定された前記多数の粒子を含む系の解析モデルを生成することを備え、
前記複数のポテンシャルエネルギ関数は粒子間の複数の結合状態にそれぞれ対応し、
本方法はさらに、
前記エネルギ調整部が、分子動力学計算中に粒子間の結合状態に応じて前記複数のポテンシャルエネルギ関数のいずれかを選択することによりポテンシャルエネルギ関数を変更することと、
前記エネルギ調整部が、ポテンシャルエネルギ関数の変更によるポテンシャルエネルギの変化量をその結合に関与する粒子の運動エネルギに振り替えることにより、変更前後でエネルギ保存則が満たされるよう前記系のもつエネルギを調整することと、
前記分子動力学計算部が、調整されたポテンシャルエネルギ及び運動エネルギから導かれる運動方程式を数値的に解くことにより各粒子の位置及び速度の少なくとも一方を計算することと、
前記出力部が、計算された各粒子の位置及び速度の少なくとも一方を、又はそれに基づき前記演算部により演算される前記多数の粒子を含む系の物理量を、出力することと、を備えることを特徴とする解析方法。
【請求項５】
請求項４に記載の方法をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、多数の粒子を含む系について各粒子の運動方程式を数値的に解くことによりその系を解析する解析装置及び解析方法に関する。
【背景技術】
【０００２】
特許文献１には、繰り込み群の手法を応用した分子動力学を用いたシミュレーション方法が記載されている。この方法によれば、計算を実行する系よりも粒子数の多い系での物理量を求めることができるので、従来シミュレーションが困難であったスケールの問題に対して精度の良いシミュレーションを行うことが可能となる。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００６−２８５８６６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
古典分子動力学法においては一般に、粒子間にはたらく相互作用が所与のポテンシャル関数として設定または仮定される。単純には例えば粒子間の結合がバネ接続であるなどというように仮定して、その結合力を表すためのポテンシャル関数が定義される。このため、原子間結合が新たに生成または解離する化学反応を分子動力学計算に組み入れるのは容易でない。量子力学を用いる第一原理分子動力学法によりこれに対処する手法もありうるが、演算負荷が非常に大きいので実際上は粒子数のごく小さい系への適用に制限されている。
【０００５】
そこで、本発明は、粒子間の結合状態の変化を分子動力学計算に比較的容易に組み入れることを可能とする解析方法及び解析装置を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明のある態様の解析装置は、多数の粒子を含む系について各粒子の運動方程式を数値的に解くことによりその系を解析する。この装置は、粒子間の相互作用を表すポテンシャルエネルギ関数を該粒子間の結合状態に応じて変更し、変更前後のポテンシャルエネルギの変化を運動エネルギの調整により相殺するエネルギ調整部と、調整されたポテンシャルエネルギ及び運動エネルギから導かれる運動方程式を数値的に解くことにより各粒子の位置及び速度の少なくとも一方を求める分子動力学計算部と、を備える。
【０００７】
本発明の別の態様は、解析方法である。この方法は、多数の粒子を含む系について各粒子の運動方程式を数値的に解くことによりその系を解析する解析方法であって、粒子間の相互作用を表すポテンシャルエネルギ関数を該粒子間の結合状態に応じて変更し、変更前後のポテンシャルエネルギの変化を運動エネルギの調整により相殺することを含む。
【０００８】
なお、以上の構成要素の任意の組み合わせ、本発明の表現を方法、装置、システム、コンピュータプログラム、データ構造、記録媒体等の間で変換したものもまた、本発明の態様として有効である。
【発明の効果】
【０００９】
本発明によれば、粒子間の結合状態の変化を分子動力学計算に反映させることができる。
【図面の簡単な説明】
【００１０】
【図１】本発明の一実施形態に係る解析処理を説明するためのフローチャートである。
【図２】本発明の一実施形態に係る解析装置の構成を模式的に示す図である。
【図３】本発明の一実施形態に係る分子動力学計算処理を説明するためのフローチャートである。
【図４】本発明の一実施形態に係る結合状態判別処理を説明するためのフローチャートである。
【図５】本発明の一実施形態に係る結合状態判別処理を概念的に説明するための図である。
【図６】本発明の一実施形態に係るポテンシャル切替処理を概念的に説明するための図である。
【図７】本発明の一実施形態に係るポテンシャル切替処理を概念的に説明するための図である。
【図８】本発明の一実施形態に係る計算例を示す図である。
【図９】本発明の一実施形態に係る計算例を示す図である。
【図１０】本発明の一実施形態に係る他の計算例を示す図である。
【発明を実施するための形態】
【００１１】
本発明の一実施形態に係る解析装置は、粒子間の相互作用を表すポテンシャルエネルギ関数を当該粒子間の結合状態に応じて変更し、各粒子について変更前後でエネルギー収支を整合させる。すなわち、変更前後でエネルギ保存則を満たすよう系のもつエネルギを調整する。例えば、ポテンシャルエネルギの変化を運動エネルギの調整により相殺する。このようにして、粒子間の結合状態の変化を分子動力学計算に直接反映させることができる。例えば、原子間の化学結合の生成及び解離を解析モデルに直接導入することができる。よって、粒子間に化学反応が生じる系について分子動力学法によるシミュレーションを行うことができる。
【００１２】
図１は、本発明の一実施形態に係る解析処理を説明するためのフローチャートである。図１に示される処理は、一実施例においては解析装置１０（図２参照）により実行される。解析装置１０はまず、解析モデルを生成する（Ｓ１０）。解析装置１０は、ユーザからの入力情報に基づいて解析モデルを作成してもよいし、既に作成されている解析モデルを読み込むことによりモデルを取得してもよい。この解析モデルは、分子動力学法に適する粒子モデルであり、Ｎ個の粒子からなる粒子系Ｓとして構成される。一実施例においては粒子は原子または分子である。
【００１３】
解析装置１０は、粒子系Ｓにおいて粒子間の相互作用を表すポテンシャルエネルギ関数（以下では適宜、ポテンシャル関数または単にポテンシャルとも称する）を設定する（Ｓ１２）。一実施例においてはポテンシャル関数は、２粒子間の距離を変数とする二体ポテンシャルの形式を有する。粒子系Ｓは、粒子間の結合状態が分子動力学計算中に変化しうる粒子群を含んでもよい。この粒子群については、複数の結合状態の各々に対応して複数のポテンシャル関数が設定される。後述するように、分子動力学計算中に結合状態に応じて複数のポテンシャル関数のうちいずれかが選択され、そのポテンシャル関数へと変更される。解析装置１０は例えば、粒子間の距離に基づいて結合状態を判別する。
【００１４】
一実施例においては、粒子群は互いに化学結合をする粒子を含んでおり、その結合前の状態における粒子間の相互作用を表すポテンシャルと、結合後の状態における粒子間の相互作用を表すポテンシャルとが設定される。例えば、粒子系Ｓは、互いに結合して分子を構成する原子を粒子として含んでいてもよい。この場合、分子を形成する前後それぞれの原子間相互作用が設定される。例えば、分子を構成していない状態の原子（すなわち未反応の原子）と他の未反応の原子との相互作用を表すポテンシャルと、分子を構成する原子（すなわち反応済の原子）と他の反応済の原子との相互作用を表すポテンシャルとを含む複数のポテンシャルが設定される。
【００１５】
解析装置１０は、粒子系Ｓに対し分子動力学計算を実行する（Ｓ１４）。粒子系Ｓにおいて各粒子は、与えられたポテンシャルのもとで刻々と位置及び速度を変化させながら自由に運動する。解析装置１０は、各粒子の運動方程式を数値的に解くことにより各粒子の位置及び速度の少なくとも一方を求める。解析装置１０は、計算結果に基づき粒子系Ｓのもつ物理量を求める。解析装置１０は、得られた結果を出力して処理を終了する。
【００１６】
また、一実施例においては、解析装置１０は、繰り込み群の手法を応用した分子動力学（以下では繰り込み群分子動力学とも称する）を粒子系Ｓに適用してもよい。この場合、分子動力学計算の実行に先立って、解析装置１０は、繰り込みにより粒子系Ｓを、それよりも粒子数の少ない系Ｓ’に対応づけてもよい。分子動力学計算は、粒子系Ｓに代えて、繰り込まれた系Ｓ’に対して行われる。解析装置１０は、繰り込まれた系Ｓ’の物理量を変換することにより粒子系Ｓのもつ物理量を求める。
【００１７】
図２は、本発明の一実施形態に係る解析装置１０の構成を模式的に示す図である。解析装置１０は演算部１２を備える。図２には、演算部１２が処理する主な機能ブロックが示されている。演算部１２は、モデル生成部１４と、エネルギ調整部１６と、分子動力学計算部１８と、を含んで構成される。
【００１８】
モデル生成部１４は、解析モデルを生成する処理を実行する。モデル生成部１４は例えば、ユーザからの入力情報に基づいて解析モデルを作成する。あるいは、モデル生成部１４は既に作成されている解析モデルを例えば記憶部２０から読み込むことによりモデルを取得してもよい。エネルギ調整部１６は、後述するように、ポテンシャル関数の変更に伴うエネルギ収支の調整を行う。分子動力学計算部１８は、粒子系Ｓまたは繰り込まれた系Ｓ’に対し分子動力学計算を実行する。
【００１９】
以上の各機能ブロックは、各種演算処理を実行するＣＰＵ、データの格納やプログラム実行のためのワークエリアとして利用されるＲＡＭなどのハードウェア、およびソフトウェアの連携によって実現される。したがって、これらの機能ブロックはハードウェアおよびソフトウェアの組み合わせによって様々な形で実現することができる。
【００２０】
解析装置１０はさらに、記憶部２０、入力部２２、及び出力部２４を含む。記憶部２０には、本発明の一実施例を実施するための分子動力学計算プログラム、エネルギ調整プログラム、及び粒子モデル生成プログラム等が格納されている。また、記憶部２０には、これらプログラムの実行に要する各種のデータが格納されている。演算部１２は記憶部２０から必要なプログラム及びデータを読み込んで処理を実行する。
【００２１】
入力部２２は、演算部１２で実行される処理に関係するユーザの入力を受けるためのキーボード、マウス等の入力装置を含む。入力部２２は、フレキシブルディスク、ＣＤ、ＤＶＤ等の記録媒体から入力を受けるよう構成されていてもよい。出力部２４は、ディスプレイ等の表示機器やプリンタ等を含む出力装置である。出力部２４は、適宜の記録媒体に結果を出力するよう構成されていてもよい。解析装置１０は例えば公知のパソコンであってもよく、演算部１２はＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等を含んで構成されてもよい。
【００２２】
図３は、本発明の一実施形態に係る分子動力学計算処理を説明するためのフローチャートである。図示されるように、この処理は、粒子の結合状態判別処理（Ｓ２０）と、結合状態に応じたエネルギ操作処理（Ｓ２６）と、を含む。結合状態判別処理（Ｓ２０）は、時刻ｔの粒子の挙動に基づき粒子間の結合状態を判別することを含む。エネルギ操作処理（Ｓ２６）は、その結合状態に応じてポテンシャル関数を変更し、ポテンシャルエネルギの変動分を運動エネルギに振り向けることを含む。
【００２３】
典型的な分子動力学計算処理における各粒子の追跡処理に結合状態判別処理（Ｓ２０）及びエネルギ操作処理（Ｓ２６）が組み込まれている。ここで追跡処理とは、各粒子の運動方程式を用いて時刻ｔの位置及び速度から微小時間Δｔ経過後の時刻ｔ＋Δｔの位置及び速度へと更新する処理である。
【００２４】
エネルギ調整部１６は、結合状態判別処理を実行する（Ｓ２０）。反復処理の第１回目の処理である場合には各粒子の初期状態に基づき粒子間の結合状態が判別される。第２回目以降の処理である場合には、前回得られた各粒子の挙動に基づき粒子間の結合状態が判別される。エネルギ調整部１６は、分子動力学計算部１８の演算結果に基づき更新された粒子間距離に基づいて結合状態を判別する。結合状態判別処理（Ｓ２０）は、結合解離判定処理（Ｓ２２）と、再結合判定処理（Ｓ２４）と、を含む。ここで粒子の挙動は、各粒子の位置、及びそこから導かれる各粒子間の距離を含む。
【００２５】
結合解離判定処理（Ｓ２２）は、粒子間距離を含む粒子の挙動に基づいて既存の粒子間結合を維持すべきか否かを判定することを含む。上述のように、前回演算された粒子間距離に基づいて結合を維持すべきか否かが判定される。一実施例においては、粒子間距離が所定のしきい値以下である場合にその粒子間の結合は維持すべきと判定され、粒子間距離が当該しきい値を超える場合にその結合状態を維持すべきでないと判定される。一実施例においては、結合状態は結合の有無の二段階に設定されており、結合状態を維持すべきでない場合にはその結合は解消される。
【００２６】
再結合判定処理（Ｓ２４）は、反応可能な状態にある粒子が他の粒子と結合するか否かを判定することを含む。例えば、結合解離判定処理（Ｓ２２）において結合が解消された粒子について他の粒子と結合するか否かを判定する。また、前回の追跡処理から引き続き反応可能な状態にある粒子について他の粒子と結合するか否かを判定してもよい。
【００２７】
再結合判定処理（Ｓ２４）は、反応可能な状態にある粒子の結合相手となる他の粒子を探索することを含んでもよい。結合可能な他の粒子が存在する場合にはその粒子との間に新たな結合が生成される。結合可能な他の粒子が存在しない場合には未反応の粒子であると判定され、次回の再結合判定処理で反応相手が再度探索される。一実施例においては、互いに反応可能な状態にある粒子間の距離が所定のしきい値以下である場合にその粒子間に新たな結合が生成され、粒子間距離が当該しきい値を超える場合に未反応の状態に維持される。ここで、既に結合状態にある粒子についてはしきい値より接近していても再結合の相手とすべきではない。よって、結合済の粒子は探索の対象外とする反応可能対象の選別処理を実行してもよい。この選別処理をしたうえで再結合判定処理を実行することにより、再結合判定処理の計算負荷を軽減することができる。
【００２８】
なお、一実施例においては、粒子間の結合状態は結合の有無の二段階のいずれかに限られず、粒子間の結合状態は相対的に強い結合と弱い結合とを含む多段階に設定されていてもよい。この場合、結合状態判別処理（Ｓ２０）においては、例えば粒子間距離に基づいて複数の結合状態のうちいずれをとるべきかが選択されてもよい。
【００２９】
また、結合状態判別処理（Ｓ２０）は、２粒子間の結合状態の判別には限られず、３粒子間または４粒子間等の多粒子間での結合状態を判別してもよい。この場合例えば、粒子間距離に基づく判定を拡張して、結合に関与する粒子のうち１つを含む微小領域に他の粒子が含まれるか否かによって結合の有無を判別してもよい。多段階の結合状態が設定されている場合には、それに対応させて多段階に微小領域を設定してもよい。
【００３０】
エネルギ調整部１６はさらに、エネルギ操作処理を実行する（Ｓ２６）。エネルギ操作処理（Ｓ２６）は、ポテンシャル変更処理（Ｓ２８）と、エネルギ調整処理（Ｓ３０）と、を含む。ポテンシャル変更処理（Ｓ２８）は、結合状態判別処理（Ｓ２０）により更新された粒子間の結合状態に対応するポテンシャル関数を各粒子について選択することを含む。設定されている複数の結合状態の各々に対応して複数のポテンシャル関数が予め定義され例えば記憶部２０に記憶されている。
【００３１】
エネルギ調整処理（Ｓ３０）は、ポテンシャル関数の変更によって生じたポテンシャルエネルギの変化量をその結合に関与する粒子の運動エネルギに振り替えることを含む。一実施例においては、エネルギ調整処理（Ｓ３０）は、ポテンシャル関数の変更前後のポテンシャルエネルギ差を、その結合に関与する各粒子の運動エネルギに均等に分配することを含む。このようにして、ポテンシャル関数の変更によるエネルギ収支を運動エネルギの調整により相殺することができる。後述のように、粒子間の結合状態の変化がポテンシャルエネルギ及び運動エネルギを通じて各粒子の運動方程式へと反映される。
【００３２】
分子動力学計算部１８は、例えば以下に述べる手法により粒子系Ｓの各粒子の位置及び速度を追跡する（Ｓ３２）。分子動力学計算部１８は、調整されたポテンシャルエネルギ及び運動エネルギから導かれる運動方程式を数値的に解くことにより各粒子の位置及び速度を求める。
【００３３】
Ｎ個の粒子からなる粒子系ＳのハミルトニアンＨは次の式（１）で表される。ここで、各粒子の質量がｍであり、粒子ｊの運動量ベクトルがｐｊであり、粒子ｊの位置ベクトル（位置座標）がｑｊであり、粒子間の相互作用ポテンシャルエネルギがφである。
【００３４】
【数１】
0005495910.tif 000002
【００３５】
ハミルトニアンＨを、ハミルトンの正準方程式に代入することにより、各粒子に対する運動方程式（２）及び（３）が得られる。
【００３６】
【数２】
0005495910.tif 000003
【００３７】
【数３】
0005495910.tif 000004
【００３８】
ここで式（３）において、ｖｊは粒子ｊの速度ベクトルである。分子動力学では、粒子系を構成する各粒子に対し、式（２）及び（３）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトル（または速度ベクトル）及び位置ベクトルが求められる。
【００３９】
一実施例においては、次式に従って数値積分を行い、各粒子の速度ベクトルｖ及び位置ベクトルｒが更新される。なお、速度ベクトルｖ及び位置ベクトルｒの両方を更新することは必須ではなく、最終的に求めるべき物理量によりいずれか一方のみを更新するようにしてもよい。
【００４０】
【数４】
0005495910.tif 000005
【００４１】
【数５】
0005495910.tif 000006
【００４２】
なお、数値積分には、Ｖｅｒｌｅｔ法が用いられることが多い。Ｖｅｒｌｅｔ法については、例えば、Ｊ．Ｍ．Ｔｈｉｊｓｓｅｎ，“Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ”，Ｃａｍｂｒｉｄｇｅ  Ｕｎｉｖｅｒｓｉｔｙ  Ｐｒｅｓｓ（１９９９）のｐ．１７５に説明されている。分子動力学計算で得られた各粒子の位置、速度に基づき、系の種々の物理量を算出することができる。
【００４３】
また、分子動力学計算部１８または演算部１２の繰り込み処理部は、繰り込み群分子動力学を粒子系Ｓに適用してもよい。一実施例においては、元の粒子系Ｓにおける粒子の質量、速度、及び位置を次式の変換則により、繰り込まれた粒子系Ｓ’における粒子の質量、速度、及び位置に変換してもよい。ここで、定数αは繰り込み因子である。なお、繰り込み群分子動力学については、特許文献１に更なる詳細が説明されている。
【００４４】
【数６】
0005495910.tif 000007
【００４５】
繰り込まれた粒子系Ｓ’においては、各粒子の速度ベクトルｖ’及び位置ベクトルｒ’は次式を用いて更新される。
【００４６】
【数７】
0005495910.tif 000008
【００４７】
図３に示されるように、演算部１２は所定の計算終了条件を満たしているかを判断する（Ｓ３４）。計算終了条件は例えば、所望の解析時間にわたる演算が終了したことを含む。計算終了条件が満たされていない場合には（Ｓ３４のＮ）、演算部１２は上述の結合状態判別処理（Ｓ２０）、エネルギ操作処理（Ｓ２６）、及び追跡処理（Ｓ３２）を再度実行する。計算終了条件が満たされている場合には（Ｓ３４のＹ）、演算部１２はこの反復処理を終了する。続いて、演算部１２は、各粒子の位置及び速度に基づき粒子系Ｓのもつ物理量（例えば温度など）を求める（Ｓ３６）。解析装置１０は、得られた結果を出力して処理を終了する（Ｓ３８）。
【００４８】
本発明の一実施例においては、粒子系Ｓのすべての粒子間の相互作用について複数のポテンシャルが設定され結合状態に応じて選択されてもよい。あるいは、粒子系Ｓには、粒子間の結合状態が分子動力学計算中に変化しない粒子群が含まれてもよい。この粒子群については、単一のポテンシャル関数が設定される。
【００４９】
よって、一実施例においては、粒子系Ｓは、粒子間の複数の結合状態に各々が対応する複数のポテンシャルエネルギ関数が定義されている第１の粒子群と、単一のポテンシャルエネルギ関数が定義されている第２の粒子群と、を含んでもよい。エネルギ調整部は、第１の粒子群について粒子間の結合状態に応じてポテンシャルエネルギ関数を選択してもよい。分子動力学計算部は、第１及び第２の粒子群について粒子の位置及び速度の少なくとも一方を求めてもよい。
【００５０】
一実施例においては、第１の粒子群は、解析モデルにおいて粒子間の化学反応が許容されている反応粒子群であってもよい。第２の粒子群は、解析モデルにおいて粒子間の化学反応が禁止されている非反応粒子群であってもよい。非反応粒子群は、例えば弾性体または剛体を表してもよい。反応粒子群の各粒子は、非反応粒子群が表す物体により画定される領域または空間において、ポテンシャルエネルギ関数の表す相互作用のもとで自由に運動可能であってもよい。反応粒子群は例えば気体または液体を表してもよい。反応粒子群は同種の原子を含んでもよいし、異種の原子を含んでもよい。解析モデルは、それぞれ異なるポテンシャルエネルギ関数が定義されている複数種類の反応粒子群を含んでもよい。解析モデルは、それぞれ異なるポテンシャルエネルギ関数が定義されている複数種類の非反応粒子群を含んでもよい。
【００５１】
反応粒子群により解析モデルに化学反応を直接導入することができる。非反応粒子群により物体を模擬することができる。このようにして、反応性を有する物質と、その反応の作用を受ける物体との相互作用を模擬することができる。これら反応粒子群及び非反応粒子群をもつ解析モデルを用いることにより、化学反応と機構との連成解析を容易に行うことができる。つまり、化学反応により生じたエネルギーを受けて動作する機構の運動解析を、統一された解析モデルにより実行することができる。従来は化学反応を機構解析に組み込むことが困難であったのに対して、上述のポテンシャル選択及びエネルギ調整を通じて反応と機構との連成解析を１つの解析モデルにより実行することができる。
【００５２】
図４は、本発明の一実施形態に係る結合状態判別処理を説明するためのフローチャートである。この結合状態判別処理は、結合解消処理（Ｓ４０）と、再結合候補選定処理（Ｓ４２）と、再結合原子特定処理（Ｓ４４）と、を含む。結合状態判別処理に続いて、ポテンシャル切替処理（Ｓ４６）が実行される。これらの処理はエネルギ調整部１６により実行される。
【００５３】
図５は、本発明の一実施形態に係る結合状態判別処理を概念的に説明するための図である。図５の左側及び右側はそれぞれ結合状態判別処理の実行前及び実行後を示す。ここでは、代表例として二原子分子を例に挙げる。結合状態として、２つの原子が結合した状態と、結合していない単原子の状態との２状態を考える。図５は、解析モデルのある局所領域の時刻ｔにおける各原子ａ１〜ａ１０の位置を示している。単原子状態の原子を白丸で示し、分子を構成する原子を黒丸で示す。図５の左側に示されるように、結合状態判別処理の実行前において、原子ａ１、ａ２、ａ９、ａ１０は単原子である。原子ａ３は原子ａ８と、原子ａ４は原子ａ６と、原子ａ５は原子ａ７と、二原子分子を形成している。一実施例においては各原子は共通の元素であり、例えば水素原子である。この場合、原子間の結合は共有結合である。
【００５４】
図４に示される結合解消処理（Ｓ４０）は、前回の演算で更新された各原子の位置に基づき分子の結合の維持または解消を判定することを含む。具体的には、エネルギ調整部１６は、各分子の原子間距離がしきい値を超えるか否かを判定する。このしきい値は、原子間の結合が許容される最大距離であり、解析的方法、実験的方法または経験的観点から適宜定めることができる。一実施例においては、エネルギ調整部１６は、原子間距離がしきい値以下である場合にその分子は結合を維持すると判定し、原子間距離が当該しきい値を超える場合にその結合を解消する。
【００５５】
再結合候補選定処理（Ｓ４２）は、ある粒子に近接する粒子を再結合の候補としてリスト化することを含む。遠方に位置する粒子との再結合は非現実的である。このように再結合の候補を近接粒子に制限することにより、再結合の相手となる粒子の探索処理を短時間で行うことができる。一実施例においては、エネルギ調整部１６は、単原子状態の各原子について所定の探索範囲にある原子をリスト化する。この探索範囲は、上述の結合しきい値よりは十分に大きな値を有する。なお、系の粒子数が少ない場合などのように探索処理の計算負荷が小さい場合には、この再結合候補選定処理を省略してもよい。
【００５６】
再結合対象特定処理（Ｓ４４）は、対象原子の再結合の相手となる原子を特定することを含む。一実施例においては、エネルギ調整部１６は、所定の特定条件が満たされる原子が存在する場合にはその原子を再結合の相手として特定し、特定条件が満たされる原子が存在しない場合には対象原子を未反応の単原子として取り扱う。特定条件は、対象原子との原子間距離が上述の結合判定しきい値以内であること、互いに最も近接している原子であること、及び、他の原子と結合して分子を形成してはいないこと、を含んでもよい。
【００５７】
図５の右側に示されるように、結合状態判別処理が実行されることにより、原子ａ３と原子ａ８との結合は解消される。これは、原子ａ３と原子ａ８との距離が結合判定しきい値を超えたためである。また、原子ａ９と原子ａ１０との結合が新たに生成される。原子間距離が結合判定しきい値以内であり、互いに最も近接しており、かつ、原子ａ９及び原子ａ１０はともに他の原子と未結合であるからである。一方、原子ａ５と原子ａ６との間には結合は生成されない。原子ａ５と原子ａ６とは原子間距離が結合判定しきい値以内であり、かつ互いに最も近接しているが、既に他の原子ａ４、ａ７と分子を形成しているからである。
【００５８】
このように、一実施例においては、一度結合した分子は原子間距離がしきい値を超えて離れるまでは結合が解消されずに維持されるようにしてもよい。周囲の原子または分子に結合状態が影響されないとの仮定をおくことにより、結合状態の判別処理を上述のように比較的簡単なアルゴリズムで実現することができる。
【００５９】
図４に示されるように、続いてポテンシャル切替処理が実行される（Ｓ４６）。エネルギ調整部１６は、分子が解消された原子間の相互作用を表すポテンシャル関数を、結合状態を表すポテンシャル関数から単原子間の相互作用を表すポテンシャル関数へと切り替える。一方、エネルギ調整部１６は、新たに結合が生成された原子間のポテンシャル関数にを、単原子間の相互作用を表すポテンシャル関数から結合状態を表すポテンシャル関数に切り替える。例えば、図５に示される原子ａ５と原子ａ６との相互作用については、単原子間の相互作用を表すポテンシャル関数から共有結合を表すポテンシャル関数へと切り替えられる。一実施例においては、単原子間の相互作用としてファンデルワールス力を仮定してもよい。
【００６０】
図６及び図７は、本発明の一実施形態に係るポテンシャル切替処理を概念的に説明するための図である。図６及び図７においては簡単な代表例として、１次元に配列された３つの原子ａｊ、ａｉ、ａｋと、原子間のポテンシャルとを示している。原子ａｊ、ａｉ、ａｋは例えばともに水素原子である。３つの原子ａｊ、ａｉ、ａｋの結合状態を図中の右上に示す。図６は時刻ｔにおける結合状態を示し、図７はその微小時間経過後の時刻ｔ＋Δｔにおける結合状態を示す。原子ａｊ、ａｉ間の距離をｒｉｊと表記し、原子ａｉ、ａｋ間の距離をｒｉｋと表記する。また、原子間距離Ｒは、原子間結合が生成または解消されるときの値を表す。結合している原子間の距離がしきい値Ｒを超えた場合に、その原子間結合が解消されると判定される。また、未結合の原子間距離がしきい値Ｒを下回ったことを１つの条件としてその原子間に新たに結合が生成される。
【００６１】
図６及び図７に示されるグラフの縦軸はポテンシャルφであり、横軸は原子間距離ｒである。原子ａｉ及び原子ａｊ間のポテンシャルφｉ−ｊを図の左側に示し、原子ａｉ及び原子ａｋ間のポテンシャルφｉ−ｋを図の右側に示す。実線で示されるポテンシャル関数φ１は原子間の結合、例えば水素分子の共有結合を表し、破線で示されるポテンシャル関数φ２は非結合原子間のファンデルワールス力を表す。この例では原子ａｊ、ａｉ、ａｋが同種の原子であるので、図示されるようにφｉ−ｊとφｉ−ｋとでφ１及びφ２はそれぞれ共通の関数である。ポテンシャル関数φｉ−ｊ、φｉ−ｋは、原子間結合がある場合には結合ポテンシャルφ１が選択され、結合がない場合には非結合ポテンシャルφ２が選択される。
【００６２】
一実施例においては、原子間結合を表すポテンシャル関数は第一原理計算に基づき求めてもよい。複数または多数の原子間距離ｒの値について第一原理計算によりポテンシャルの値を求め、この計算結果を近似する関数をポテンシャル関数としてもよい。近似関数は、古典論で使用される形式のポテンシャル関数であってもよい。例えば、共有結合を表す場合にはモース（Ｍｏｒｓｅ）ポテンシャルを用いてもよい。また、ファンデルワールス力を表す場合にはモース（Ｍｏｒｓｅ）ポテンシャルを用いてもよいし、レナード−ジョーンズ（Ｌｅｎｎａｒｄ−Ｊｏｎｅｓ）ポテンシャルを用いてもよい。すなわち、第一原理計算から得られる結果に近似関数をフィッティングさせるよう近似関数に含まれるパラメタを設定することにより、ポテンシャル関数を定めてもよい。第一原理計算については公知の手法を適宜用いることができるので、詳細の説明は省略する。
【００６３】
図６に示されるように、時刻ｔにおいては原子ａｉ、ａｋが結合しており、原子ａｊ、ａｉは結合していない。よって、原子ａｉ、ａｋ間のポテンシャルφｉ−ｋは結合ポテンシャルφ１が選択され、原子ａｊ、ａｉ間のポテンシャルφｉ−ｊは非結合ポテンシャルφ２が選択されている。この状態において各原子の位置及び速度の更新演算を１回行った状態が図７に示される。
【００６４】
図７に示される時刻ｔ＋Δｔにおいては、原子ａｋと結合していた原子ａｉが原子ａｊの近傍へと移動している。原子ａｉ、ａｋ間の距離ｒｉｋは結合判定しきい値Ｒを超えているので、原子ａｉ、ａｋ間の結合は解消される。原子ａｊと原子ａｉとが新たに結合する。この場合、原子ａｊ、ａｉ間のポテンシャルφｉ−ｊは非結合ポテンシャルφ２から結合ポテンシャルφ１へと切り替えられる（図中の矢印Ａ）。このとき、切替前のポテンシャルエネルギφ２（ｒｉｊ）と切替後のポテンシャルエネルギφ１（ｒｉｊ）とは異なっておりポテンシャルエネルギ差ΔＥ（＝φ２（ｒｉｊ）−φ１（ｒｉｊ））が生じる。
【００６５】
実際の系では、このポテンシャルエネルギ差ΔＥを得るまでに原子は加速され、運動エネルギをもっているはずである。よって、本解析においてもこの余剰のエネルギ差ΔＥを各原子に与えエネルギ収支の整合をとることとする。具体的には例えば、結合に関与した各原子に均等にエネルギ差ΔＥを分配する。図７に示される例においては、原子ａｊ及び原子ａｉのそれぞれの運動エネルギにΔＥ／２を付加する。すなわち、切替後の原子の運動エネルギｐ’２／２ｍは、切替前の原子の運動エネルギｐ２／２ｍと付加エネルギΔＥ／２との和で与えられる。ここで、ｐ’及びｐはそれぞれ切替前後の原子の運動量を表し、ｍは原子の質量を表す。
【００６６】
このようにして、エネルギ収支を合わせつつ結合状態の変化を反映させることができる。切替後のポテンシャルエネルギ及び運動エネルギを用いて粒子の位置及び速度を追跡することができる。
【００６７】
計算例を示す。図８は、本発明の一実施形態に用いた解析モデルを示す。図８の左側は初期状態を示し、図８の右側は計算終了時の状態を示す。図９は、図８に示す系の温度変化の計算結果を示す。図８及び図９は、１．０ｎｍ×１．０ｎｍ×１．０ｎｍの立方体空間に４つの水素原子を均等に分布させた状態を初期状態として、水素原子から水素分子への反応を計算した例を示す。ポテンシャル関数はモースポテンシャルを用いている。初期状態の温度は１５Ｋとし、１００ピコ秒（ｐｓ）の実時間に相当するシミュレーションを行っている。図８の右側に示されるように、初期状態から１００ピコ秒経過後には２つの水素分子が生成されている。
【００６８】
図９に示されるように、系の温度は概ね２段階に増加している。まず、初期状態から１０〜２０ピコ秒経過したときに約１５０００Ｋへと温度が増加し、初期状態から３０〜４０ピコ秒経過したときに約３００００Ｋへと温度が増加している。この段階的な温度変化はそれぞれ１つ目の水素分子の生成と２つ目の水素分子の生成に対応している。なお温度が振動的に大きく変化しているのは、系の粒子数が４つと少ないためである。粒子数を多くすることにより平均化され、振動的な温度変動は小さくなる。
【００６９】
５０ピコ秒以降計算終了時である１００ピコ秒までは２つ目の水素分子が既に生成されており、定常状態にあると言える。そこで、５０〜１００ピコ秒の平均温度を反応後の系の温度とする。算出された反応後の温度は約２８５００Ｋであり、これに基づき得られる反応熱Ｑは約７１０ｋＪ／ｍｏｌであった。反応熱Ｑは反応前後の温度差ΔＴと定積モル比熱との積である。ただし定積モル比熱Ｃｖは（３／２）Ｒと仮定した。Ｒは気体定数である。一方、公知文献によれば、水素原子から水素分子の生成反応（Ｈ＋Ｈ→Ｈ２）の反応熱は約４３５．８８ｋＪ／ｍｏｌである（例えば、化学ハンドブック編集委員会編、「化学ハンドブック」、オーム社）。
【００７０】
本実施例による解析によって得られた反応熱の値と文献に示された値とは数値のオーダーは合っている。よって、定性的に見て、本実施例は良好な結果を出していると考えられる。定量的な精度についても、例えば解析に用いる各種パラメタ（例えばポテンシャル関数のパラメタ）を調整することにより、さらに向上させることが可能であると考えられる。
【００７１】
第２の計算例を図１０に示す。この例は、水素により駆動されるロータリーエンジンをモデルとして反応と機構との連成解析を計算したものである。ロータリーエンジンモデル３０は、回転部３２と、回転部３２の運動を内部空間に拘束する壁部３４と、回転部３２に回転力を作用させる水素３６と、を含む。ロータリーエンジンモデル３０は、回転部３２、壁部３４、及び水素３６のそれぞれを表す粒子群により構成される粒子モデルである。当初生成された粒子モデルに対し繰り込み処理をして粒子数を低減させた粒子モデルが図１０に示されている。
【００７２】
回転部３２及び壁部３４の粒子については弾性体であるとして、それらの粒子間の相互作用は図１０及び次式（８）に示す２次のポテンシャル関数φａを設定している。ｒｉｊは粒子ｉと粒子ｊとの距離を示し、初期状態における粒子ｉ、ｊ間の距離がｒ０である。弾性体粒子と水素原子を表す粒子との相互作用は、図１０及び次式（９）に示すポテンシャル関数φｂを設定している。ポテンシャル関数φｂは、ポテンシャル関数φａで与えられる斥力のみを残し、引力をゼロとしたものである。水素原子を表す粒子間の相互作用は、ポテンシャル関数φｃである。ポテンシャル関数φｃは上述のように、結合状態に応じてφ１とφ２とで切り替えられる。
【００７３】
【数８】
0005495910.tif 000009
【００７４】
【数９】
0005495910.tif 000010
【００７５】
図１０には、初期状態４０と、計算の結果得られた所定時間経過後の第２の状態４２、さらに所定時間経過後の第３乃至第５の状態４４、４６、４８を示す。反応部分についてはポテンシャル関数の切替を許容し、機構部分についてはポテンシャル関数を固定して繰り込み分子動力学法を適用することにより、図示されるように、反応と機構との連成解析を実現することができる。
【符号の説明】
【００７６】
１０  解析装置、  １２  演算部、  １４  モデル生成部、  １６  エネルギ調整部、  １８  分子動力学計算部、  ２０  記憶部、  ２２  入力部、  ２４  出力部。
【図１】
0005495910.tif 000011
【図２】
0005495910.tif 000012
【図３】
0005495910.tif 000013
【図４】
0005495910.tif 000014
【図５】
0005495910.tif 000015
【図６】
0005495910.tif 000016
【図７】
0005495910.tif 000017
【図８】
0005495910.tif 000018
【図９】
0005495910.tif 000019
【図１０】
0005495910.tif 000020