JP5332875B2

JP5332875B2 - 分子安定構造探索装置、分子安定構造探索プログラム、および分子安定構造探索方法

Info

Publication number: JP5332875B2
Application number: JP2009107672A
Authority: JP
Inventors: 義明谷田
Original assignee: Fujitsu Ltd
Current assignee: Fujitsu Ltd
Priority date: 2009-04-27
Filing date: 2009-04-27
Publication date: 2013-11-06
Anticipated expiration: 2029-04-27
Also published as: JP2010257297A

Description

本件は、分子安定構造探索装置、分子安定構造探索プログラム、および分子安定構造探索方法に関する。分子安定構造探索装置は、ポテンシャルエネルギーが低く安定した分子構造を探索する装置である。また、分子安定構造探索プログラムは、コンピュータをそのような分子安定構造探索装置として動作させるプログラムである。また、分子安定構造探索方法は、上記のような安定した分子構造を探索する方法である。

複数の原子が配列された分子構造を有する分子は、その分子を構成する各原子の位置に応じたポテンシャルエネルギーを有している。そして、分子構造は、分子が有するポテンシャルエネルギーが低い程安定した構造となる。

従来、安定した分子構造を探索する方法として、シミュレーティドアニリーリング（ＳＡ：ＳｉｍｕｌａｔｅｄＡｎｎｅａｌｉｎｇ）法と呼ばれる方法が知られている（例えば、特許文献１参照。）。ＳＡ法は、分子をデータ上で加熱した上で原子の位置を変えてポテンシャルエネルギーが低い構造を探りつつ、徐々に冷却して、最終的にポテンシャルエネルギーが極小となる安定構造を見つける方法である。

ここで、原子を移動させようとしたときに、その原子の移動が、ポテンシャルエネルギーのエネルギー障壁によって妨げられることがしばしばある。そして、このようなエネルギー障壁を超えた原子の移動によって安定構造を見つけることができることが多い。上記のＳＡ法では、このようなエネルギー障壁を越えるためには、分子を、そのエネルギー障壁を上回る高温で加熱する必要がある。ＳＡ法はシミュレーションなので、加熱温度をいくらでも高く設定することは可能である。しかしながら、そのような高温での加熱は、シミュレーション上での分子構造を根本的に変えてしまい、処理対象の分子についての安定構造の探索を不可能にしてしまうおそれがある。

そこで、上記のようなエネルギー障壁に関わる問題を回避して安定構造を探索する量子アニーリング（ＱＡ：ＱｕａｎｔｕｍＡｎｎｅａｌｉｎｇ）法と呼ばれる方法が提案されている。このＱＡ法は、ＳＡ法とは異なり、トンネル効果によってエネルギー障壁を見かけ上減じることで、このエネルギー障壁を越えた安定構造の探索が可能となっている。

ＱＡ法では、まず、このトンネル効果の度合いを示すパラメータであるトンネリング項が大きめに設定された上で、安定構造の探索が行われる。そして、ＱＡ法では、ポテンシャルエネルギーの低下に応じてトンネリング項を減じる等といった操作が行われて、最終的に、ポテンシャルエネルギーが極小となる安定構造が探索される。

特開平９−１０６３９２号公報

しかしながら、上記のＱＡ法では、シミュレーションにおいて分子構造を変える度に、上記のトンネリング項を用いた複雑な演算を繰り返す必要があり、分子の安定構造を効率的に探索することが難しいという問題がある。

本件は上記事情に鑑み、次のような分子安定構造探索装置、分子安定構造探索プログラム、および分子安定構造探索方法を提供することを目的とするものである。分子安定構造探索装置は、分子の安定構造を効率的に探索することができる装置である。また、分子安定構造探索プログラムは、コンピュータをそのような分子安定構造探索装置として動作させるプログラムである。また、分子安定構造探索方法は、分子の安定構造を効率的に探索することができる方法である。

上記目的を達成する分子安定構造探索装置の基本形態は、原子移動部と、移動判定部と、エネルギー算出部と、探索部とを備えている。

原子移動部は、複数の原子それぞれが各座標点に配置され、原子の座標移動に伴ってポテンシャルエネルギーが変化する分子について、その分子を構成する原子のうちの移動対象原子の座標をデータ上で移動させるものである。

移動判定部は、上記原子移動部での今回の移動対象原子の移動が、ポテンシャルエネルギーのエネルギー障壁によりその移動対象原子の移動が妨げられる第１種の移動であるかその第１種の移動を除く第２種の移動であるかを判定するものである。

エネルギー算出部は、上記移動対象原子の移動が、上記第１種の移動であるか上記第２種の移動であるかに応じて、それぞれ、次のような第１のポテンシャル関数、および第２のポテンシャル関数に基づいて、上記分子のポテンシャルエネルギーを算出するものである。第１のポテンシャル関数は、上記分子を構成する複数の原子の座標点に依存するとともに上記エネルギー障壁を減じるトンネリング項を含む関数である。また、第２のポテンシャル関数は、その第１のポテンシャル関数からそのトンネリング項を除いた関数である。

探索部は、次のような処理により、上記分子のポテンシャルエネルギーが極小となる、その分子を構成する複数の原子の座標点を探索するものである。即ち、この探索部は、上記原子移動部による移動対象原子の移動と、上記エネルギー算出部による上記分子のポテンシャルエネルギーの算出を交互に繰り返す。そして、この探索部は、その繰返しを、上記第１種の移動の場合は算出されたポテンシャルエネルギーの低下に応じて上記第１のポテンシャル関数のトンネリング項の影響をその影響の消滅を含んで減じながら実行する。

また、上記目的を達成する分子安定構造探索プログラムの基本形態は、コンピュータ内で実行され、そのコンピュータ上に、上記原子移動部と、上記移動判定部と、上記エネルギー算出部と、上記探索部とを構築するものである。

また、上記目的を達成する分子安定構造探索方法の基本形態は、原子移動過程と、移動判定過程と、エネルギー算出過程と、探索過程とを備えている。

原子移動過程は、複数の原子それぞれが各座標点に配置され、原子の座標移動に伴ってポテンシャルエネルギーが変化する分子について、その分子を構成する原子のうちの移動対象原子の座標をデータ上で移動させる過程である。

移動判定過程は、上記原子移動過程での今回の移動対象原子の移動が、ポテンシャルエネルギーのエネルギー障壁によりその移動対象原子の移動が妨げられる第１種の移動であるかその第１種の移動を除く第２種の移動であるかを判定する過程である。

エネルギー算出過程は、上記移動対象原子の移動が、上記第１種の移動であるか上記第２種の移動であるかに応じて、上記分子のポテンシャルエネルギーを算出する過程である。このエネルギー算出過程では、上記第１種の移動であるか上記第２種の移動であるかに応じて、このポテンシャルエネルギーが、それぞれ、次のような第１のポテンシャル関数、および第２のポテンシャル関数に基づいて算出される。第１のポテンシャル関数は、上記分子を構成する複数の原子の座標点に依存するとともに上記エネルギー障壁を減じるトンネリング項を含む関数である。また、第２のポテンシャル関数は、その第１のポテンシャル関数からそのトンネリング項を除いた関数である。

探索過程は、次のような処理により、上記分子のポテンシャルエネルギーが極小となる、その分子を構成する複数の原子の座標点を探索する過程である。即ち、この探索過程では、上記原子移動過程による移動対象原子の移動と、上記エネルギー算出過程による上記分子のポテンシャルエネルギーの算出が交互に繰り返される。そして、この探索過程では、その繰返しが、上記第１種の移動の場合は算出されたポテンシャルエネルギーの低下に応じて上記第１のポテンシャル関数のトンネリング項の影響をその影響の消滅を含んで減じながら実行される。

本件によれば、分子の安定構造を効率的に探索することができる。

基本形態について上述した分子安定構造探索装置の具体的な一実施形態として動作するコンピュータを示す図である。図１に示すコンピュータのハードウェア構成図である。基本形態について説明した分子安定構造探索プログラムの具体的な実施形態が記憶されたＣＤ−ＲＯＭを示す概念図である。基本形態について説明した分子安定構造探索装置の具体的な実施形態を示すブロック図である。図４の分子安定構造探索装置で実行される分子安定構造探索処理の流れを示すフローチャートにおけるステップＳ１０１からステップＳ１１３までの処理を示す図である。図４の分子安定構造探索装置で実行される分子安定構造探索処理の流れを示すフローチャートにおけるステップＳ１１４以降の処理を示す図である。二面角の変更を伴った原子の移動を示す模式図である。図５のステップＳ２００の処理を詳細に示す図である。単純な分子構造の一例を示す図である。図９の分子Ｍ１に対する分子安定構造探索処理により安定構造が探索される様子を模式的に示す図である。第２実施形態において、二面角の変更を伴う移動であった場合に移動後エネルギーＵ’を算出する処理の流れを示すフローチャートである。２種類のトンネリング項の使い分けについて説明する図である。

以下、上記に基本形態について説明した分子安定構造探索装置、分子安定構造探索プログラム、および分子安定構造探索方法の具体的な実施形態について、図面を参照して説明する。

まず、第１実施形態について説明する。

図１は、基本形態について上述した分子安定構造探索装置の具体的な一実施形態として動作するコンピュータを示す図である。また、図２は、図１に示すコンピュータのハードウェア構成図である。

このコンピュータ１００は、本体部１１０、画像表示装置１２０、キーボード１３０、およびマウス１４０を備えている。

本体部１１０は、後述するＣＰＵ１１１やハードディスク装置（ＨＤＤ）１１３等を内蔵し各種処理を実行する装置である。また、この本体部１１０は、ＣＤ−ＲＯＭ２００が装填されるＣＤ−ＲＯＭ装填口１１０Ａを有している。画像表示装置１２０は、本体部１１０からの指示によりモニタ１２１に画像表示を行う装置である。キーボード１３０は、このコンピュータ１００にユーザの指示や文字情報を入力するためのものである。マウス１４０は、モニタ１２１上の任意の位置を指定することによりその位置に応じた指示をこのコンピュータ１００に入力するものである。

本体部１１０の内部には、図２に示すように、ＣＰＵ１１１、主メモリ１１２、ＨＤＤ１１３、ＣＤ−ＲＯＭドライブ１１４、およびＩ／Ｏインタフェース１１５が内蔵されている。

ＣＰＵ１１１は、各種プログラムを実行するものである。主メモリ１１２は、ＨＤＤ１１３に格納されたプログラムが読み出されてＣＰＵ１１１での実行のために展開されるメモリである。ＨＤＤ１１３は、各種プログラムやデータが保存されるメモリである。ＣＤ−ＲＯＭドライブ１１４は、ＣＤ−ＲＯＭ装填口１１０Ａから装填されたＣＤ−ＲＯＭ２００にアクセスするものである。Ｉ／Ｏインタフェース１１５は、このコンピュータ１００が外部機器との間でデータのやりとりを行うためのインタフェースである。これらの各種要素、図１にも示す画像表示装置１２０、キーボード１３０、およびマウス１４０は、バス１１６を介して相互に接続されている。

以上に説明したコンピュータ１００に、後述の分子安定構造探索プログラムがインストールされると、このコンピュータ１００が分子安定構造探索装置として動作する。

図３は、基本形態について説明した分子安定構造探索プログラムの具体的な実施形態が記憶されたＣＤ−ＲＯＭを示す概念図である。

このＣＤ−ＲＯＭ２００に記憶され、コンピュータ１００を分子安定構造探索装置として動作させる分子安定構造探索プログラム３００が、基本形態について上述した分子安定構造探索プログラムの具体的な一実施形態に相当する。この分子安定構造探索プログラム３００は、初期構造データ入力部３１０と、原子移動部３２０と、移動判定部３３０と、エネルギー算出部３４０と、探索部３５０と、探索結果保存／表示部３６０とを有する。この分子安定構造探索プログラム３００の各要素の詳細については後述する。

図４は、基本形態について説明した分子安定構造探索装置の具体的な実施形態を示すブロック図である。

この図４に示す分子安定構造探索装置４００は、図３の分子安定構造探索プログラム３００が図１のコンピュータ１００にインストールされて実行されることによって構成されるものである。この分子安定構造探索装置４００は、初期構造データ入力部４１０と、原子移動部４２０と、移動判定部４３０と、エネルギー算出部４４０と、探索部４５０と、探索結果保存／表示部４６０とを有する。

図３の分子安定構造探索プログラム３００が、図１のコンピュータ１００にインストールされると、その分子安定構造探索プログラム３００の各要素が、分子安定構造探索装置４００の各要素を構築する。即ち、分子安定構造探索プログラム３００の初期構造データ入力部３１０が、分子安定構造探索装置４００の初期構造データ入力部４１０を構築する。また、分子安定構造探索プログラム３００の原子移動部３２０が、分子安定構造探索装置４００の原子移動部４２０を構築する。また、分子安定構造探索プログラム３００の移動判定部３３０が、分子安定構造探索装置４００の移動判定部４３０を構築する。また、分子安定構造探索プログラム３００のエネルギー算出部３４０が、分子安定構造探索装置４００のエネルギー算出部４４０を構築する。また、分子安定構造探索プログラム３００の探索部３５０が、分子安定構造探索装置４００の探索部４５０を構築する。そして、分子安定構造探索プログラム３００の探索結果保存／表示部３６０が、分子安定構造探索装置４００の探索結果保存／表示部４６０を構築する。

ここで、これら分子安定構造探索装置４００の各要素は、コンピュータのハードウェアとそのコンピュータで実行されるＯＳやアプリケーションプログラムとの組合せで構成されている。これに対し、図３の分子安定構造探索プログラム３００の各要素は、アプリケーションプログラムのみによって構成されている。

本実施形態の原子移動部４２０が、上述の基本形態における原子移動部の一例に相当する。また、移動判定部４３０が、上述の基本形態における移動判定部の一例に相当する。また、エネルギー算出部４４０が、上述の基本形態におけるエネルギー算出部の一例に相当する。また、探索部４５０が、上述の基本形態における探索部の一例に相当する。

以下、図４の分子安定構造探索装置４００の各要素の詳細について、この分子安定構造探索装置４００で実行される分子安定構造探索処理の流れに沿って説明する。また、この分子安定構造探索装置４００の各要素を説明することによって、図３に示す分子安定構造探索プログラム３００の各要素も併せて説明する。

図５は、図４の分子安定構造探索装置で実行される分子安定構造探索処理の流れを示すフローチャートにおけるステップＳ１０１からステップＳ１１３までの処理を示す図である。また、図６は、図４の分子安定構造探索装置で実行される分子安定構造探索処理の流れを示すフローチャートにおけるステップＳ１１４以降の処理を示す図である。これら図５および図６のフローチャートが示す分子安定構造探索処理が、基本形態について上述した分子安定構造探索方法の具体的な一実施形態に相当する。

図５および図６のフローチャートが示す分子安定構造探索処理は、処理対象の分子のポテンシャルエネルギーが最も低くなる分子構造（安定構造）を、モンテカルロシミュレーションによって各原子を移動しながら探索する処理である。

この分子安定構造探索処理は、ユーザによる所定の操作画面に対する操作によって、分子安定構造探索処理の開始が指示されるとスタートする。処理がスタートすると、初期構造データ入力部４１０が、処理対象の分子について、各原子が所定の初期位置に配置された初期構造を表わす初期構造データの、ユーザ操作による分子安定構造探索装置４００への入力を受け付ける（ステップＳ１０１）。本実施形態では、この初期構造データは、処理対象の分子を構成する各原子の初期位置を、所定の仮想三次元空間上での座標で定義したものである。尚、この初期構造データは、仮想三次元空間上で分子構造を構築する公知の設計支援ツール等によって予め構築されたものである。

次に、ステップＳ１０１で入力された初期構造データに基づいて、エネルギー算出部４４０が、その初期構造データに対応した分子が有するポテンシャルエネルギーを算出する（ステップＳ１０２）。

分子のポテンシャルエネルギーは、その分子を構成する複数の原子それぞれが有するポテンシャルエネルギーの総和となる。また、各原子のポテンシャルエネルギーは、隣接する原子との相互作用を含んだものとなる。このような各原子のポテンシャルエネルギーをＥｉ（Ｘｉ，Ｙｉ，Ｚｉ）とすると、分子のポテンシャルエネルギーＵは、以下の（１）式で表される。ここで、沿え字の「ｉ」は、分子を構成する複数の原始それぞれに対応する識別番号である。そして、「（Ｘｉ，Ｙｉ，Ｚｉ）」は、識別番号「ｉ」の原子の、仮想三次元空間上での位置座標である。

尚、各原子のポテンシャルエネルギーＥｉ（Ｘｉ，Ｙｉ，Ｚｉ）の算出については、公知であるのでここでは説明を割愛する。

ステップＳ１０２で分子のポテンシャルエネルギーを算出すると、エネルギー算出部４４０は、算出結果を探索部４５０に渡す。

次に、原子移動部４２０が、モンテカルロシミュレーションによる移動対象原子を、分子を構成する複数の原始の中から所定の選択確率で選択する（ステップＳ１０３）。また、分子を構成する複数の原始それぞれは、分子構造に起因する各原子の自由度に応じて動くことができる。そして、各原子が動くことができる可動方向は、その原子の自由度に応じた数だけ存在する。ステップＳ１０３では、原子移動部４２０は、移動対象原子を選択するとともに、その原子について存在する可動方向の中から、その原子の移動方向を所定の選択確率で選択する。

移動対象原子と移動方向が決まると、原子移動部４２０は、その原子を、その移動方向に、所定距離ｄＲだけ移動する（ステップＳ１０４）。

上記のステップS１０３の処理とステップＳ１０４の処理とを合わせた処理が、上述の基本形態における原子移動過程の一例に相当する。

次に、移動判定部４３０が、ステップＳ１０４での原子の移動が、以下に説明する二面角の変更を伴った移動である否かを判定する（ステップＳ１０５）。このステップＳ１０５の処理が、上述の基本形態における移動判定過程の一例に相当する。

図７は、二面角の変更を伴った原子の移動を示す模式図である。

この図７には、一例として、第１から第４の原子Ａ１，…，Ａ４が、第１から第３の結合軸Ｂ１，Ｂ２，Ｂ３で互いに結ばれてなる分子構造において、第２の原子Ａ２が二面角の変更を伴って動かされる様子が示されている。そして、図７のパート（Ａ）には、移動前の分子構造が示されている。また、図７のパート（Ｂ）には、移動後の分子構造が示されている。

この図７の例では、第１の原子Ａ１と第２の原子Ａ２と第３の原子Ａ３とが形作る面、および、第１の原子Ａ１と第３の原子Ａ３と第４の原子Ａ４とが形作る面がなす二面角が、第２の原子Ａ２の移動によって「θ１」から「θ２」に変更されている。

たんぱく質の分子や薬剤候補化合物の分子等では、図７に一例を示すような二面角の変更を伴う原子の移動は、他の原子との相互作用が大きく、分子のポテンシャルエネルギーの複雑で振幅が大きい増減を招くことが多い。一方で、原子間の距離を伸ばしたり縮めたりする移動は、他の原子との相互作用が小さく、この移動に伴う分子のポテンシャルエネルギーの増減は、概ね振幅が小さいことが多い。

図５のステップＳ１０５では、ステップＳ１０４での原子の移動が、処理対象の分子のポテンシャルエネルギーの複雑で振幅が大きい増減を招くことがある、二面角の変更を伴った移動であるか否かが判定される。

まず、移動判定部４４０が二面角の変更を伴った移動ではないと判定した場合（ステップＳ１０５におけるＮＯ判定）、エネルギー算出部４４０が次のような算出処理を実行する（ステップＳ１０６）。このステップＳ１０６では、エネルギー算出部４４０は、移動後の原子位置が反映された構造データに対応した分子のポテンシャルエネルギー（移動後エネルギー）Ｕ’を、上記の（１）式によって算出する。

一方、移動判定部４４０が二面角の変更を伴った移動であると判定した場合（ステップＳ１０５におけるＹＥＳ判定）、エネルギー算出部４４０は、移動後エネルギーＵ’を、以下に説明する一連の処理によって算出する（ステップＳ２００）。

図８は、図５のステップＳ２００の処理を詳細に示す図である。

この図８に示すように、エネルギー算出部４４０は、まず、移動対象原子について、上記のステップＳ１０４での移動前のポテンシャルエネルギーＥＡと、移動後のポテンシャルエネルギーＥＢとを算出する（ステップＳ２０１）。この算出は、ステップＳ１０２やステップＳ１０６の処理における原子のポテンシャルエネルギーＥｉ（Ｘｉ，Ｙｉ，Ｚｉ）の算出と同様に、公知の手法によって行われる。

次に、エネルギー算出部４４０は、以下に説明するトンネリング項Ｑを用いたハミルトニアン（トンネリング項導入ハミルトニアン）Ｈ’を作成し、そのトンネリング項導入ハミルトニアンＨ’の対角化を行う（ステップＳ２０２）。

このステップＳ２０２の処理では、エネルギー算出部４４０は、まず、上記の移動前後のポテンシャルエネルギーＥＡ，ＥＢを対角項とした二行二列の対角行列を作成する。さらに、その対角行列における「０」の項を、所定の値を有するトンネリング項Ｑに置き換えることで、図８に示す二行二列のトンネリング項導入ハミルトニアンＨ’を得る。そして、エネルギー算出部４４０は、この二行二列のトンネリング項導入ハミルトニアンＨ’の対角化によって、トンネリング項導入ハミルトニアンＨ’を、次の（２）式に示す形に変形する。

この（２）式が示すトンネリング項導入ハミルトニアンＨ’における対角項のうち左上の項が、移動対象原子についての、トンネリング項Ｑが導入された、移動前のポテンシャルエネルギーＥＡ’である。また、対角項のうち右下の項が、トンネリング項Ｑが導入された、移動後のポテンシャルエネルギーＥＢ’である。

エネルギー算出部４４０は、上記のようなトンネリング項導入ハミルトニアンＨ’の対角化を実行すると、その対角化済みのトンネリング項導入ハミルトニアンＨ’から右下の項を抽出する行列演算を行う（ステップＳ２０３）。このステップＳ２０３の処理により、エネルギー算出部４４０は、移動対象原子についての、トンネリング項Ｑが導入された、移動後のポテンシャルエネルギーＥＢ’を得る。

次に、エネルギー算出部４４０は、処理対象の分子について、移動対象原子以外の原子のポテンシャルエネルギーの総和Ｕ１を算出する（ステップＳ２０４）。このステップＳ２０４の算出は、上記のステップＳ１０２と同等な処理によって行われる。ただし、移動対象原子に隣接する原子に対する、この移動対象原子からの相互作用については、上記のトンネリング項Ｑが導入された、移動後のポテンシャルエネルギーＥＢ’に基づく相互作用が適用される。

ここで、移動対象原子のポテンシャルエネルギーが増加した場合には、その移動対象原子に隣接する原子に対する、この移動対象原子からの相互作用も増加する。このとき、上記の（２）式から、トンネリング項Ｑが導入された、移動後のポテンシャルエネルギーＥＢ’は、トンネリング項Ｑが導入されていない移動後のポテンシャルエネルギーＥＢよりも低くなることが分かる。つまり、移動対象原子の移動によって、その移動対象原子のポテンシャルエネルギーが増加したとしても、上記のようなトンネリング項Ｑの導入により、その移動対象原子のポテンシャルエネルギーの増し分が減じられる。これにより、その移動対象原子に隣接する原子に対する、この移動対象原子からの相互作用も、トンネリング項Ｑの導入により減じられることとなる。その結果、上記の移動対象原子以外の原子のポテンシャルエネルギーの総和Ｕ１の増し分が減じられる。

次に、エネルギー算出部４４０は、ステップＳ２０４で算出された移動対象原子以外の原子のポテンシャルエネルギーの総和Ｕ１に、移動対象原子について上記のステップＳ２０３で得られた移動後のポテンシャルエネルギーＥＢ’を加える（ステップＳ２０５）。この加算により、移動後の構造データに対応した分子について、トンネリング項Ｑが導入された移動後エネルギーＵ’が得られる。

ここで、上述したように、移動対象原子について上記のステップＳ２０３で得られた移動後のポテンシャルエネルギーＥＢ’は、トンネリング項Ｑの導入により低くなる。そして、移動対象原子以外の原子のポテンシャルエネルギーの総和Ｕ１については、トンネリング項Ｑの導入により増し分が減じられる。その結果、このステップＳ２０５で得られる移動後エネルギーＵ’についても増し分が減じられることとなる。このことは、移動対象原子の、二面角の変更を伴った移動がもたらすポテンシャルエネルギーの増減の振幅が緩和され、延いては、そのような移動によって生じるエネルギー障壁が減じられることを意味する。

以上に説明したステップＳ２００の一連の処理で使われる数式を総合したものが、上述の基本形態にいう第１のポテンシャル関数の一例に相当する。また、図５のステップＳ１０６の処理で使われる（１）式が、上述の基本形態にいう第２のポテンシャル関数の一例に相当する。

図５のステップＳ１０６の処理、あるいは、以上に説明したステップＳ２００の一連の処理によって得られた移動後エネルギーＵ’は、エネルギー算出部４４０から探索部４５０に渡される。

上記のステップＳ１０６の処理とステップＳ２００の処理とを合わせた処理が、上述の基本形態におけるエネルギー算出過程の一例に相当する。

移動後エネルギーＵ’を受け取ると、探索部４５０は、移動後エネルギーＵ’から移動前の分子のポテンシャルエネルギーＵを差し引いて差分エネルギーｄＵを求める（ステップＳ１０７）。そして、探索部４５０は、その差分エネルギーｄＵが負の値となっているか否か、即ち、移動後エネルギーＵ’が移動前のポテンシャルエネルギーＵよりも小さくなっているか否かを判定する（ステップＳ１０８）。

差分エネルギーｄＵが負の値となっていることは、移動対象原子が移動された後の分子が、移動前よりもポテンシャルエネルギーが低い安定した状態に移行していることを意味する。このため、差分エネルギーｄＵが負の値となっている場合には（ステップＳ１０８におけるＹＥＳ判定）、探索部４５０は、移動対象原子の位置として、分子全体を安定させる、移動後の位置を採用する（ステップＳ１０９）。そして、探索部４５０は、処理対象の分子のポテンシャルエネルギーＵとして、ステップＳ１０６あるいはステップＳ２００で得られた移動後エネルギーＵ’を採用する（ステップＳ１１０）。

一方、差分エネルギーｄＵが「０」以上の値となっていることは、原子が移動された後の分子が、移動前よりもポテンシャルエネルギーが高く不安定な状態に移行していることを意味する。ただし、ステップＳ１０４での所定距離ｄＲだけの移動で不安定な状態に移行したとしても、移動対象原子をさらに移動させたときにポテンシャルエネルギーが低い安定した状態に移行する可能性は未だ残っている。

そのため、差分エネルギーｄＵが「０」以上の値となっている場合には（ステップＳ１０８におけるＮＯ判定）、探索部４５０は、まず、差分エネルギーｄＵをパラメータとした確率関数ＥＸＰ（−ｄＵ／Ｂ）に、ステップＳ１０７で得られた差分エネルギーｄＵを代入する。そして、探索部４５０は、その代入結果を、「０」から「１．０」までの間での一様乱数と比較する（ステップＳ１１１）。

そして、上記の代入結果が一様乱数よりも大きい場合には（ステップＳ１１１におけるＹＥＳ判定）、ステップＳ１０９に進んで、探索部４５０は、移動後の原子の位置の採用と、移動後エネルギーＵ’の採用とを行う。

一方、上記の代入結果が一様乱数以下の場合には（ステップＳ１１１におけるＮＯ判定）、探索部４５０は、移動前の原子の位置の採用（ステップＳ１１２）と、移動前のポテンシャルエネルギーＵの採用（ステップＳ１１３）とを行う。

以上、説明したように、移動対象原子の移動について、移動後の原子の位置、即ち、移動後の分子における分子構造と、移動前のポテンシャルエネルギーＵの決定が終了すると、処理は図６のステップＳ１１４に進む。

このステップＳ１１４では、探索部４５０は、処理対象の分子についてのモンテカルロシミュレーションが終了したか否かを判定する（ステップＳ１１４）。本実施形態では、この判定は、処理対象の分子中で移動可能な全ての方向についての所定ステップ分の原子の移動が、その処理対象の分子を構成する全ての原子について行われたか否かを判定することで行われる。

未終了であると判定された場合には（ステップＳ１１４におけるＮＯ判定）、探索部４５０は、この時点で求められている分子のポテンシャルエネルギーＵが、所定の閾値よりも小さいか否かを判定する（ステップＳ１１５）。

そして、分子のポテンシャルエネルギーＵが閾値以上の場合には（ステップＳ１１４におけるＮＯ判定）、図５のステップＳ１０３以降の処理が再度実行される。本実施形態では、このような処理の再実行は、原子移動部４２０、およびエネルギー算出部４４０に対する探索部４５０の指示によって行われる。

一方、分子のポテンシャルエネルギーＵが閾値よりも小さい場合には（ステップＳ１１４におけるＹＥＳ判定）、まず、探索部４５０が、上記のトンネリング項Ｑを「０」に設定する（ステップＳ１１６）。そして、探索部４５０の指示により、図５のステップＳ１０３以降の処理が再度実行される。このステップＳ１１６の処理を経た場合には、上記のステップＳ２００では、ステップＳ１０６の処理と同等な処理が実行されることとなる。

以上に説明したステップＳ１０３からステップＳ１１６までの処理が、探索部４５０によってモンテカルロシミュレーションが終了したと判定される（ステップＳ１１４におけるＹＥＳ判定）まで繰り返される。そして、モンテカルロシミュレーションが終了すると、探索部４５０は、処理対象の分子について最終的に得られた分子構造を表わす構造データを、探索結果保存／表示部４６０に渡す。

上述のステップＳ１０７からステップＳ１１６までの処理が、上述の基本形態における探索過程の一例に相当する。

探索結果保存／表示部４６０は、最終的な構造データを受け取ると、その構造データを所定のメモリに保存するとともに、図１のモニタ１２１に、その構造データに対応した分子の画像を表示する。この構造データの表示をもって、図５および図６のフローチャートが示す分子安定構造探索処理が終了する。

以下、ここまでに説明した分子安定構造探索処理によって、分子の安定構造を効率的に探索することができることについて、単純な分子構造を例に挙げて説明する。

図９は、単純な分子構造の一例を示す図である。

この図９には、一連の識別番号ｉが割り当てられた９つの原子Ａ１，…，Ａ９が、ジグザグに結合されてなる分子Ｍ１が示されている。さらに、説明を単純化するために、この図９の分子Ｍ１では、識別番号ｉが「３」の原子Ａ３のみが移動可能であり、しかも、自由度は「１」で、図中の矢印Ｄ１が示す方向のみに移動可能であるとする。さらに、この矢印Ｄ１が示す移動方向が、二面角の変更を伴う移動であるとする。

また、この識別番号ｉが「３」の原子Ａ３では、図９中で実線で描かれた位置が初期位置Ｐ１である。そして、矢印Ｄ１が示す移動方向上の、図９中で点線で描かれた位置に、この分子Ｍ１のポテンシャルエネルギーが最低となる安定位置Ｐ２が存在すると仮定する。

この図９の分子Ｍ１に対して、図５および図６の分子安定構造探索処理が施されると、識別番号ｉが「３」の原子Ａ３が安定位置Ｐ２に配置された安定構造が、以下に説明するように探索されることとなる。

図１０は、図９の分子Ｍ１に対する分子安定構造探索処理により安定構造が探索される様子を模式的に示す図である。

この図１０のパート（Ａ）には、識別番号ｉが「３」の原子Ａ３を、図９の矢印Ｄ１が示す移動方向に連続的に移動したときの、分子Ｍ１のポテンシャルエネルギーＵの変化を示す第１曲線Ｌ１が、上記の初期位置Ｐ１に配置された原子Ａ３とともに示されている。また、図１０のパート（Ｂ）には、この第１曲線Ｌ１が、上記の安定位置Ｐ２に配置された原子Ａ３とともに示されている。

上述したように二面角の変更を伴う原子の移動は、分子のポテンシャルエネルギーの複雑で振幅の大きな増減を招くことが多い。図１０のパート（Ａ）およびパート（Ｂ）には、このようなポテンシャルエネルギーの複雑で振幅の大きな増減が第１曲線Ｌ１で示されている。

ここで、第１曲線Ｌ１が示すポテンシャルエネルギーの変化には、ポテンシャルエネルギーが最も低くなる基底状態の他に、この基底状態に比べれば高いが、ポテンシャルエネルギーが低くなった準安定状態が複数存在している。そして、準安定状態どうしの間がエネルギー障壁で隔てられている。このとき、移動対象原子Ａ３の初期位置Ｐ１が、図１０のパート（Ａ）のように準安定状態に対応していたとする。一方で、上記の安定位置Ｐ２は、図１０のパート（Ｂ）のように初期位置Ｐ１とは別の基底状態に対応した位置となっている。

上述したように、図５および図６の分子安定構造探索処理では、基本的には、所定距離ｄＲだけ原子を移動した後の移動後エネルギーＵ’が、移動前のポテンシャルエネルギーＵよりも小さくなっていれば、その移動後の原子の位置が採用される。また、移動後エネルギーＵ’が、移動前のポテンシャルエネルギーＵ以上の場合には、確率関数ＥＸＰ（−ｄＵ／Ｂ）に差分エネルギーｄＵを代入した代入結果が、上記の一様乱数よりも大きいときに移動後の原子の位置が採用される。

ここで、図１０のパート（Ａ）のように移動前の原子Ａ３が、いずれかの準安定状態に対応した位置（例えば、初期位置Ｐ１）に配置されていたとする。

また、本実施形態とは異なり、ポテンシャルエネルギーの複雑で振幅の大きな増減を招くことが多い二面角の変更を伴った原子の移動についても、移動後エネルギーＵ’を上述のトンネリング項を用いずに算出する場合について考える。

この場合、原子Ａ３を移動した後の移動後エネルギーＵ’は、エネルギー障壁に阻まれて必ず移動前のポテンシャルエネルギーＵよりも大きくなる。また、上述したように二面角の変更を伴って行われる原子の移動に伴うポテンシャルエネルギーの変化では増減の振幅が大きい。従って、原子の移動前後についての差分エネルギーｄＵが大きくなり、確率関数ＥＸＰ（−ｄＵ／Ｂ）に差分エネルギーｄＵを代入した代入結果が小さくなる。その結果、移動後の原子の位置が採用されにくくなり、安定位置Ｐ２に原子が配置された安定構造が見つかり難くなってしまう。

これに対し、本実施形態では、このような二面角の変更を伴った原子の移動については、上述したようにトンネリング項を用いることで、ポテンシャルエネルギーの増減の振幅を見かけ上緩和させた上で移動後エネルギーＵ’が算出される。

図１０のパート（Ｃ）には、分子Ｍ１のポテンシャルエネルギーＵの、増減の振幅が緩和された変化を示す第２曲線Ｌ２が示されている。

本実施形態におけるこのような算出では、エネルギー障壁が見かけ上低くなっている。そのため、ある準安定状態に対応する位置から原子Ａ３を移動したときのポテンシャルエネルギーＵの移動後の増分が小さくなる。即ち、原子Ａ３の移動前後についての差分エネルギーｄＵが小さくなり、確率関数ＥＸＰ（−ｄＵ／Ｂ）に差分エネルギーｄＵを代入した代入結果が大きくなる。その結果、移動後の原子Ａ３の位置が採用され易くなり、安定位置Ｐ２までの移動を阻むエネルギー障壁を超え易くなる。そして、本実施形態では、ある程度まで分子のポテンシャルエネルギーＵが下がった段階で、安定位置Ｐ２の周辺位置に対応したエネルギー障壁を超えたものとみなしトンネリング項が「０」に設定される。これにより、安定位置Ｐ２の周辺に至った原子Ａ３が、その安定位置Ｐ２を通り過ぎてしまうことが防がれる。本実施形態によれば、このようなトンネリング項を用いた演算処理により、原子Ａ３が安定位置Ｐ２に配置された安定構造を確実に見つけることができる。

また、本実施形態では、トンネリング項を用いた上述したような複雑な演算処理は、二面角の変更を伴った原子の移動についてのみ行われる。その結果、このような複雑な演算処理の実行が必要最小限に抑えられるので効率的である。即ち、本実施形態によれば、分子の安定構造を効率的に探索することができる。

そして、上述したように、たんぱく質の分子や薬剤候補化合物の分子等では、二面角の変更を伴う原子の移動は、他の原子との相互作用が大きく、分子のポテンシャルエネルギーの複雑で振幅が大きい増減を招くことが多い。このため、トンネリング項を用いた演算処理の処理対象を、二面角の変更を伴う原子の移動としたことは、上記のような分子の安定構造の探索に非常に有効となる。

このことは、上述の基本形態に対し、以下に説明する応用形態が好適であることを意味している。この応用形態では、上記移動判定部は、上記移動対象原子の移動が上記分子の分子構造における二面角の変更を伴う移動であるときに、その移動対象原子の移動を上記第１種の移動であると判定するものとなっている。

本実施形態の移動判定部４３０は、この応用形態における移動判定部の一例にも相当している。

また、本実施形態では、上述したように、ある程度まで分子のポテンシャルエネルギーＵが下がった段階でトンネリング項が「０」に設定される。一般的なＱＡ法では、ポテンシャルエネルギーの低下に伴ってトンネリング項を徐々に下げていく手法がとられることが多い。これに対し、本実施形態では、分子構造が、安定構造に近付くぎりぎりまで、ある大きさのトンネリング項を維持することで、安定構造探索のスピードが速められている。

このことは、上述の基本形態に対し、以下に説明する応用形態が好適であることを意味している。この応用形態では、上記探索部は、上記分子のポテンシャルエネルギーとして、所定の閾値よりも小さいポテンシャルエネルギーが算出された場合に、上記第１のポテンシャル関数のトンネリング項の影響を消滅させるものとなっている。

本実施形態の探索部４５０は、この応用形態における探索部の一例にも相当している。

次に、第２実施形態について説明する。

この第２実施形態は、移動対象原子の移動が、二面角の変更を伴う移動であった場合に移動後エネルギーＵ’を算出する処理が、上述の第１実施形態とは異なっている。

一方で、この第２実施形態は、分子安定構造探索プログラムや分子安定構造探索装置の基本的な構成については上述の第１実施形態の分子安定構造探索プログラム３００や分子安定構造探索装置４００と同等である。また、この第２実施形態は、分子安定構造探索処理の基本的な流れについても、上述の第１実施形態における分子安定構造探索処理と同等である。

そこで、以下では、第２実施形態の分子安定構造探索プログラム、分子安定構造探索装置、および分子安定構造探索処理の基本的な流れについての図示と重複説明とを割愛し、第１実施形態との上記の相違点に注目した説明を行う。尚、以下では、第１実施形態の分子安定構造探索装置４００を示す図４、および第１実施形態の分子安定構造探索処理の基本的な流れ示す図５とを適宜に参照する。

図１１は、第２実施形態において、二面角の変更を伴う移動であった場合に移動後エネルギーＵ’を算出する処理の流れを示すフローチャートである。

尚、この図１１では、第１実施形態においてこのような処理の流れを示す図８のフローチャートと同等なステップについては、図８と同じ符号が付されている。以下では、この同等なステップについては重複説明を省略する。

本実施形態では、エネルギー算出部４４０が、二面角の変更を伴う移動について移動後エネルギーＵ’を算出する場合に、互いに大きさが異なる第１および第２の２種類のトンネリング項を適宜に使い分ける。即ち、エネルギー算出部４４０は、処理対象の分子を二分したときの各原子群それぞれに対応する互いに値が異なる２種類のトンネリング項を有し、移動対象原子が属する原子群に対応したトンネリング項を使って算出する。

この図１１のフローチャートでは、上述のステップＳ２０１の処理の後、エネルギー算出部４４０は、まず、移動対象原子が、上記の２つの原子群のうちの一方の原子群である所定原子群に属しているか否かを判定する（ステップＳ３０１）。そして、移動対象原子が、その所定原子群に属さずもう一方の原子群に属している場合には（ステップＳ３０１におけるＮＯ判定）、第１のトンネリング項Ｑ１が採用される（ステップＳ３０２）。一方、移動対象原子が、その所定原子群に属している場合には（ステップＳ３０１におけるＹＥＳ判定）、第１のトンネリング項Ｑ１よりも値が大きな第２のトンネリング項Ｑ２が採用される（ステップＳ３０３）。

そして、ステップＳ３０２あるいはステップＳ３０３で採用されたトンネリング項を使って、以降のステップＳ２０２からステップＳ２０５までの処理が実行されて、移動後エネルギーＵ’が求められる。

本実施形態において、上記のように、移動対象原子が、分子中のどの原子群に属しているかによって２種類のトンネリング項を使い分けることには、以下に説明する意味がある。

図１２は、２種類のトンネリング項の使い分けについて説明する図である。

この図１２には、第１から第１３までの１３個の原子Ａ１，…，Ａ１３が、ジグザグに結合されてなる分子Ｍ２が示されている。さらに、この図１２には、この分子Ｍ２を二分する第１および第２の原子群Ｍ２＿１，Ｍ２＿２が示されている。そして、例えば第３の原子Ａ３のように、第１の原子群Ｍ２＿１内の原子については第１のトンネリング項Ｑ１が使われる。また、第９の原子Ａ９のように、第２の原子群Ｍ２＿２内の原子については、第１のトンネリング項Ｑ１よりも値が大きい第２のトンネリング項Ｑ２が使われる。

ここで、たんぱく質の分子や薬剤候補化合物の分子等では、原子の移動に伴うポテンシャルエネルギーの増減の振幅が、分子中の他の原子に比べて著しく大きくなる箇所が存在することがある。そこで、本実施形態では、ポテンシャルエネルギーの増減の振幅が特に大きい原子群（ここでの例では第２の原子群Ｍ２＿２）内の原子について、相対的に値が大きいトンネリング項Ｑ２が使われることとなっている。これにより、原子の移動に伴うポテンシャルエネルギーの増減の振幅が大きくなる箇所の原子についても、そのような増減の振幅が十分に緩和されることとなる。

一方で、ポテンシャルエネルギーの増減の振幅を緩和し過ぎると、移動の前後でポテンシャルエネルギーの明確な差が小さくなり過ぎてしまい、却って安定構造の探索が難しくなってしまう。しかしながら、本実施形態では、上記のような相対的に値が大きいトンネリング項Ｑ２の使用が、ポテンシャルエネルギーの増減の振幅が特に大きい原子群に限定されていることから、上記のような問題が回避されることとなっている。

尚、分子中のどの原子群が、ポテンシャルエネルギーの増減の振幅が特に大きくなるかについては、実際の分子を使った実験等により経験的に把握される。そして、上述の初期構造データの入力時に、分子構造を二分したときの各原子群を示す情報と、各原子群とトンネリング項とを互いに対応付ける対応情報とが、初期構造データと共に入力される。そして、その入力された各情報に基づいて、上記のようなトンネリング項の使い分けが実行される。

本実施形態では、このトンネリング項の使い分けにより、分子を構成する各原始に対してポテンシャルエネルギーの増減の振幅が適切に緩和された、良好な安定構造の探索処理が実行されることとなる。

このことは、上述の基本形態に対し、以下に説明する応用形態が好適であることを意味する。この応用形態では、上記エネルギー算出部は、上記分子を複数の原子群に分けたときの各原子群それぞれに対応付けられた互いに値が異なる複数のトンネリング項を有している。そして、このエネルギー算出部は、上記移動対象原子の移動が上記第１種の移動である場合において、該分子のポテンシャルエネルギーを、上記移動対象原子が属する原子群に対応したトンネリング項を使って算出する。

本実施形態のエネルギー算出部４４０は、この応用形態におけるエネルギー算出部の一例に相当する。

尚、上記では、基本形態について上述した分子安定構造探索装置の具体的な実施形態として、モンテカルロシミュレーションによって原子を１つずつ移動させて分子構造を変えながら安定構造を探索する形態を例示した。しかしながら、基本形態について上述した分子安定構造探索装置は、この実施形態に限るものではなく、例えば、分子軌道計算を利用して複数の原子を一度に移動させて分子構造を変えながら安定構造を探索する形態であっても良い。

また、上記では、上述の基本形態における探索部の一例として、分子のポテンシャルエネルギーが所定の閾値よりも小さくなるとトンネリング項を「０」に設定する探索部４５０を例示した。しかしながら、上述の基本形態における探索部はこれに限るものではなく、例えば、分子のポテンシャルエネルギーの低下に伴って、トンネリング項を徐々に「０」近付けるものであっても良い。

また、上記では、複数のトンネリング項を使い分けるタイプの応用形態におけるエネルギー算出部の一例として、２種類のトンネリング項を使い分けるエネルギー算出部４４０を例示した。しかしながら、上述のタイプの応用形態におけるエネルギー算出部はこれに限るものではなく、例えば、３種類以上のトンネリング項を使い分けるものであっても良い。

１００コンピュータ
１１０本体部
１１１ＣＰＵ
１１２主メモリ
１１３ＨＤＤ
１１４ＣＤ−ＲＯＭドライブ
１１５Ｉ／Ｏインタフェース
１１０ＡＣＤ−ＲＯＭ装填口
１２０画像表示装置
１２１モニタ
１３０キーボード
１４０マウス
２００ＣＤ−ＲＯＭ
３００分子安定構造探索プログラム
３１０初期構造データ入力部
３２０原子移動部
３３０移動判定部
３４０エネルギー算出部
３５０探索部
３６０探索結果保存／表示部
４００分子安定構造探索装置
４１０初期構造データ入力部
４２０原子移動部
４３０移動判定部
４４０エネルギー算出部
４５０探索部
４６０探索結果保存／表示部

Claims

複数の原子それぞれが各座標点に配置され、原子の座標移動に伴ってポテンシャルエネルギーが変化する分子について、該分子を構成する原子のうちの移動対象原子の座標をデータ上で移動させる原子移動部と、
前記原子移動部での今回の移動対象原子の移動が、ポテンシャルエネルギーのエネルギー障壁により該移動対象原子の移動が妨げられる第１種の移動であるか該第１種の移動を除く第２種の移動であるかを判定する移動判定部と、
前記移動対象原子の移動が、前記第１種の移動であるか前記第２種の移動であるかに応じて、それぞれ、前記分子を構成する複数の原子の座標点に依存するとともに前記エネルギー障壁を減じるトンネリング項を含む第１のポテンシャル関数、および該第１のポテンシャル関数から該トンネリング項を除いた第２のポテンシャル関数に基づいて、該分子のポテンシャルエネルギーを算出するエネルギー算出部と、
前記原子移動部による移動対象原子の移動と、前記エネルギー算出部による前記分子のポテンシャルエネルギーの算出を、前記第１種の移動の場合は算出されたポテンシャルエネルギーの低下に応じて前記第１のポテンシャル関数のトンネリング項の影響を該影響の消滅を含んで減じながら、交互に繰り返すことにより、前記分子のポテンシャルエネルギーが極小となる、該分子を構成する複数の原子の座標点を探索する探索部とを備えたことを特徴とする分子安定構造探索装置。
前記移動判定部は、前記移動対象原子の移動が前記分子の分子構造における二面角の変更を伴う移動であるときに、該移動対象原子の移動を前記第１種の移動であると判定するものであることを特徴とする請求項１記載の分子安定構造探索装置。
前記探索部は、前記分子のポテンシャルエネルギーとして、所定の閾値よりも小さいポテンシャルエネルギーが算出された場合に、前記第１のポテンシャル関数のトンネリング項の影響を消滅させるものであることを特徴とする請求項１又は２記載の分子安定構造探索装置。
前記エネルギー算出部は、前記分子を複数の原子群に分けたときの各原子群それぞれに対応付けられた互いに値が異なる複数のトンネリング項を有しており、前記移動対象原子の移動が前記第１種の移動である場合において、該分子のポテンシャルエネルギーを、前記移動対象原子が属する原子群に対応したトンネリング項を使って算出するものであることを特徴とする請求項１から３のうちいずれか１項記載の分子安定構造探索装置。
コンピュータ内で実行され、該コンピュータ上に、
複数の原子それぞれが各座標点に配置され、原子の座標移動に伴ってポテンシャルエネルギーが変化する分子について、該分子を構成する原子のうちの移動対象原子の座標をデータ上で移動させる原子移動部と、
前記原子移動部での今回の移動対象原子の移動が、ポテンシャルエネルギーのエネルギー障壁により該移動対象原子の移動が妨げられる第１種の移動であるか該第１種の移動を除く第２種の移動であるかを判定する移動判定部と、
前記移動対象原子の移動が、前記第１種の移動であるか前記第２種の移動であるかに応じて、それぞれ、前記分子を構成する複数の原子の座標点に依存するとともに前記エネルギー障壁を減じるトンネリング項を含む第１のポテンシャル関数、および該第１のポテンシャル関数から該トンネリング項を除いた第２のポテンシャル関数に基づいて、該分子のポテンシャルエネルギーを算出するエネルギー算出部と、
前記原子移動部による移動対象原子の移動と、前記エネルギー算出部による前記分子のポテンシャルエネルギーの算出を、前記第１種の移動の場合は算出されたポテンシャルエネルギーの低下に応じて前記第１のポテンシャル関数のトンネリング項の影響を該影響の消滅を含んで減じながら、交互に繰り返すことにより、前記分子のポテンシャルエネルギーが極小となる、該分子を構成する複数の原子の座標点を探索する探索部とを構築することを特徴とする分子安定構造探索プログラム。
複数の原子それぞれが各座標点に配置され、原子の座標移動に伴ってポテンシャルエネルギーが変化する分子について、該分子を構成する原子のうちの移動対象原子の座標をデータ上で移動させる原子移動過程と、
前記原子移動過程での今回の移動対象原子の移動が、ポテンシャルエネルギーのエネルギー障壁により該移動対象原子の移動が妨げられる第１種の移動であるか該第１種の移動を除く第２種の移動であるかを判定する移動判定過程と、
前記移動対象原子の移動が、前記第１種の移動であるか前記第２種の移動であるかに応じて、それぞれ、前記分子を構成する複数の原子の座標点に依存するとともに前記エネルギー障壁を減じるトンネリング項を含む第１のポテンシャル関数、および該第１のポテンシャル関数から該トンネリング項を除いた第２のポテンシャル関数に基づいて、該分子のポテンシャルエネルギーを算出するエネルギー算出過程と、
前記原子移動過程による移動対象原子の移動と、前記エネルギー算出過程による前記分子のポテンシャルエネルギーの算出を、前記第１種の移動の場合は算出されたポテンシャルエネルギーの低下に応じて前記第１のポテンシャル関数のトンネリング項の影響を該影響の消滅を含んで減じながら、交互に繰り返すことにより、前記分子のポテンシャルエネルギーが極小となる、該分子を構成する複数の原子の座標点を探索する探索過程とを備えたことを特徴とする分子安定構造探索方法。