JPWO2005029385A1 - 分子シミュレーション方法及び装置 - Google Patents

Abstract

シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、ＱＭ空間に対して非経験的分子軌道法を適用し、ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション方法は、記憶部から、シミュレーション対象の分子または分子の一部を構造データを取り出してＱＭ空間及びＭＭ空間に分割する段階と、ＱＭ空間に関する非経験的分子軌道法における全エネルギー表式の一部を経験的ポテンシャルで置き換える段階と、を有する。

Description

本発明は、量子化学の手法によって分子シミュレーションを行うための方法及び装置に関し、特に、量子化学の理論手法のうち、非経験的分子軌道法と分子力場法を融合してひとつの理論体系として扱うＱＭ／ＭＭ(Quantum Mechanics/Molecular Mechanics)法による分子シミュレーション方法及び装置に関する。

量子化学理論の発展や計算機技術の進歩によって、計算により、化合物分子の構造、物性や、分子内の化学結合や分子軌道、電子状態などを精度よくシミュレーションできるようになってきた。そのような手法の中でも経験的パラメータに原則として依存しない非経験的分子軌道法は、過去数十年にわたる数多くの計算例から、その理論的枠組の正しさが実証されてきた。近年、非経験的分子軌道法の生体分子への適応が、創薬や機能性食品の基礎研究の観点から求められている。例えば、特定のタンパク質の活性部位に高度に選択的に結合するリガンドを発見することは、特定疾患に有効な薬剤の探索の上で極めて有効である。このようなリガンドは、そのタンパク質の活性を阻害するという意味では、インヒビターとして機能する。そこで、リガンドの候補となるべき極めて多数の化合物から実際の候補化合物を探索するスクリーニングが行われることになるが、その場合、実際に化学反応を起こさせることなく、計算機シミュレーションによって候補化合物をスクリーニングすることができれば、全体のスクリーニングに要する時間の大幅な短縮が望めるほか、大量合成が困難なタンパク質に対するリガンドの探索や、未だに化学合成例のないリガンドを対象にすることができるなど、大きな利点を得ることができる。そこで、非経験的分子軌道法は、計算機によるこのような分子シミュレーションを行うために利用されるようになってきている。

しかしながら、非経験的分子軌道法では、その全電子計算に必要とされる膨大な計算量から、計算できる分子の大きさは依然として中規模の分子に限定される状況にある。例えば、タンパク質とリガンドの複合体の非経験的分子軌道計算は、上述したように薬候補化合物の探索に重要な情報をもたらすが、複合体の全系の電子状態計算は、最新のスーパーコンピュータをにもってしても長時間を要するものであって、企業における研究開発の時間スケールでは一般に容認できるものではない。

ところで一般に、タンパク質では、それが非常にアミノ酸残基数が多いものであっても、リガンドなどと結合する部位は限られた場所であり、したがって、電子状態についての正確なシミュレーションはその部位及びその近傍についてだけ行えばよいと考えられている。そこで、多数の分子からなる化学系において、分子もしくは分子の一部を、注目している化学現象の起こるＱＭ(Quantum Mechanics)空間と、それ以外の２次的なＭＭ(Molecular Mechanics)空間とに分割し、前者を非経験的分子軌道法など量子力学的取り扱いにより、後者を分子力場法などの経験的ポテンシャルで記述する、ＱＭ／ＭＭ法が提唱されている[1]。この方法の利点は：
（１）注目しているＱＭ空間での化学現象は、量子力学に基づく方法により記述されるので、計算結果について高い信頼性が期待され；
（２）分子系の大部分を占めるＭＭ空間は、計算負荷の少ない経験的ポテンシャルにより記述され、古典力学的取扱いで処理されるので、計算時間が大幅に短縮される；
ということである。ＱＭ／ＭＭ法は、多数提案されている分子シミュレーションの理論の中でも、計算結果の信頼性維持と計算時間の短縮という、分子シミュレーションに対する相反する要請を満たす技術として、ここ数年、注目されている。

このように、従来の全電子計算に対する計算時間の短縮が、ＱＭ／ＭＭ法の特徴である。５００アミノ酸残基からなるタンパク質で計算時間を見積もると、全計算時間の９５％以上をＱＭ空間に対する計算に消費されている。５００アミノ酸残基の全体に対してＱＭ法を適用した場合には、その計算には、５万軌道の程度の原子軌道が必要である。計算時間は、用いる原子軌道の数の３乗に比例して増減する。ここでＱＭ／ＭＭ法を適用し、活性部位近傍のみをＱＭ空間として扱う場合、計算しなければならない原子軌道の数は５千軌道程度になり、１／１０００に計算時間が短縮される。特に、タンパク質＋基質複合体の安定構造や酵素の反応機構の計算には、このような計算を繰り返すことになるので、ＱＭ／ＭＭ法以外の手法によって意味のある分子シミュレーションを行おうとすることは、現実的でない。

図１は、タンパク質−リガンド複合体におけるＱＭ空間とＭＭ空間との分割を概念的に示す図である。この図で示すように、リガンド１がタンパク質２に結合している活性部位の近傍をＱＭ空間２１とし、それ以外をＭＭ空間２２としている。

しかしながら、ＱＭ／ＭＭ法では、ＱＭ空間に属する原子とＭＭ空間に属する原子との間の相互作用の記述方法に課題が残っている。ＱＭ空間に属する原子とＭＭ空間に属する原子が相互に共有結合を直接形成していない場合には、それらの相互作用をクーロン力やファンデルワールス力などで表現することにより、原理的な困難を生じることなく、ＱＭ空間とＭＭ空間との接続を取り扱うことが可能である。しかしながら、共有結合を形成している場合、その一方の原子がＱＭ空間に属し、他方がＭＭ空間に属することになり、深刻な問題が発生する。例えば，タンパク質の活性部位近傍をＱＭ空間にする場合、共有結合であるタンパク質の主鎖をＱＭ−ＭＭ接合面が横切ることになり、隣接する原子間のいわゆるσ（シグマ）結合が人為的に切断されることになる。

図２は、ＱＭ空間とＭＭ空間との分割を説明する図であって、タンパク質におけるペプチド結合の一例を示している。ここでは図示破線で示すように、Ｎ末端（図示左端）から数えて２番目のアミノ酸残基内において主鎖のＣ−Ｎ結合の位置で、ＱＭ空間とＭＭ空間とを分けている。このＱＭ／ＭＭ法に特徴的な空間分割では：
１）共有結合に関与している電子の一方をＱＭ空間の属する原子に帰属させるため、不対電子が人為的にその原子に発生し、その原子がラジカルになってしまう。その結果、系の化学的性質が大きく歪められ、ＱＭ空間の波動関数の正しい記述が得られないことになる；
２）ＭＭ空間に属する原子は、電荷をもつ質点として表現されるため、連続的な電子分布をもつＱＭ空間と質点からなる不連続なＭＭ空間が、ひとつの化学結合を境にして相互作用することになり、系の滑らかで連続的なポテンシャル関数が得られない。このため、ＱＭ空間とＭＭ空間との境界領域近傍の原子の挙動が正しく記述されない危険性がある；
という、理論上の問題を提起する。

前者については、リンクアトム(Link Atom)法と呼ばれる、不対電子に対して共有結合で結合する水素類似原子を人為的に導入する方法が、伝統的に用いられている[2]。しかしながら、リンクアトム法では、元々の系には存在しない水素類似原子を導入するため、余計なエネルギー項が出現し、その補正が必要となるという問題がある。また、導入された水素類似原子自体が、波動関数の性質に変化をもたらす可能性もある。最近では、ＭＭ空間に帰属された原子を、水素類似原子に置き直す工夫もされているが、後者の問題の解決策とはなっていない。

さらには、Bond Orbitalと呼ばれる局在した軌道を共有結合にあてがい、不対電子の発生を防ぐ方法も提案されている[3]。しかしながら、化学結合の軸方向に固定された軌道を用いるため、原子核の移動と共に軌道の回転などの補正が必要となり、理論が複雑になるという難点がある。

ＱＭ空間とＭＭ空間との接合部近傍のポテンシャルが連続的でないという問題は、構造最適化もしくは分子動力学による時間発展の計算を行った時、顕著になる。すなわち、近傍の原子核が、計算上、本来の位置から大きくずれて、例えば、タンパク質の活性部位の構造を正しく維持できなくなる可能性がある。

この課題を解決する方法として、ＯＮＩＯＭ(our own n-layered integrated molecular orbital + molecular mechanics method)が提唱されている[4]。この方法では、系全体をＭＭ領域として一旦扱うことで、境界領域でのポテンシャルの不連続性の問題は回避しているものの、ＱＭ空間の原子群を全く独立した系として別途計算するため、ＭＭ空間からの影響がＱＭ空間に自己無撞着的に反映されないという問題が指摘されている[3]。

以下、本明細書で引用した参考文献を列挙する。
[1] J. Gao, "Methods and Applications of Combined Quantum Mechanical and Molecular Mechanical Potentials" in "Reviews in Computational Chemistry", Vol. 7, K. B. Lipkowitz and D. B. Boyd, Editors, VCH Publishers, Inc. New York, 1996 [2] M. J. Field, P. A.Bash and M. Karplus, J. Comp. Chem., Vol. 11, 700 (1990) [3] D. M. Philipp, R. A. Friesner, J. Comp. Chem., Vo; 20, 1468 (1999) [4] T. Vreven, K. Morokuma, J. Comp. Chem., Vol. 21, 1419 (2000) [5] 米澤 康滋、高田 俊和、佐久間 俊広、中田 一人、中村 春木、生物物理、第４３巻、Supplement 1, B198（日本生物物理学会第４１会年会講演予稿集）(2003) [6] 電子構造の理論入門：新しい量子化学（上）、ｐ１５５、東京大学出版社、１９８７年 [7] P. Pulay, "Direct Use of Gradient for investigating Molecular Energy Surfaces", in "Modern Theoretical Chemistry 4", H. F. Scharfer, Editor, Plenum Press, New York and London, 1977 [8] J. A. Pople, R. Krishnan, H. B. Schlegel and J. S. Binkley, Int. J. Quant. Chem., Vol. 13, 225 (1979)

ＱＭ／ＭＭ法を、生体分子などに対する実用的な分子シミュレーション技術にするためには、前述したように、１）不対電子による人為的な波動関数の歪曲、２）不連続なポテンシャルによる分子構造の不安定性、の課題を解決しなければならない。

ＭＭ法で用いる経験ポテンシャルＥ(MM)は、(1)式で示すように、結合距離Ｅ_bonds、結合角Ｅ_angles、２面角Ｅ_torsions、非共有結合Ｅ_nonbondsで与えられる。

ここで、Ｅ^MM(０)は、平衡状態での全ＭＭエネルギーを表す。このため、全系を経験的ポテンシャルで記述すれば、異常な分子骨格が計算中に出現することはない。他方、ＱＭ空間自体は、量子力学的な手法により厳密に記述されているので、ＱＭ空間の分子構造が異常をきたすことは基本的にありえない。この議論から、ＱＭ空間の内側で、ＭＭ空間に接合する表層部分に、理論的に矛盾なく、ＭＭ空間から経験的ポテンシャルを染み込ませることができれば、分子構造の不安定性の問題を克服できることが分かる。その時、ＱＭ空間の末端の原子に発生する不対電子を、エネルギー表式に基づいて、論理的に隠蔽できることも、解決のための必須の条件である。

本発明者らは、ＱＭ領域とＭＭ領域の境界に局所化された分子軌道を配置し、この局所軌道をＱＭ領域に対してFrozen Orbitalとして計算に取り入れることによって、ＱＭ領域とＭＭ領域とのスムースかつ精密な接続を実現できることを既に予想している[5]。

そこで本発明の目的は、不対電子を隠蔽しつつ、ＭＭ空間に接合するＱＭ空間の表層に、(1)式で与えられる経験的ポテンシャルを染み込ませることのできる、論理的に整合性の取れた分子シミュレーション方法及び装置を提供することにある。

本発明の第１の様相に従えば、分子シミュレーション方法は、シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、ＱＭ空間に対して非経験的分子軌道法を適用し、ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション方法であって、記憶部から、シミュレーション対象の分子または分子の一部を構造データを取り出してＱＭ空間及びＭＭ空間に分割する段階と、ＱＭ空間に関する非経験的分子軌道法における全エネルギー表式の一部を経験的ポテンシャルで置き換える段階と、を有する。

本発明の第２の様相に従えば、分子シミュレーション方法は、シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、ＱＭ空間に対して非経験的分子軌道法を適用し、ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション方法であって、記憶部から、シミュレーション対象の分子または分子の一部を構造データを取り出してＱＭ空間及びＭＭ空間に分割し、さらに、ＱＭ空間を、ＭＭ空間に隣接する領域である表層ＱＭ領域と、表層ＱＭ領域以外の領域であるＱＭ領域とに分割する段階と、表層ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を得る段階と、ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を得る段階と、表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換える段階と、ＱＭ空間についての非経験的分子軌道法による全エネルギーを求める段階と、を有する。

本発明の第３の様相に従えば、分子シミュレーション装置は、シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、ＱＭ空間に対して非経験的分子軌道法を適用し、ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション装置であって、シミュレーション対象の分子または分子の一部の構造データを格納する記憶部と、記憶部から、シミュレーション対象の分子または分子の一部を構造データを取り出してＱＭ空間及びＭＭ空間に分割し、さらに、ＱＭ空間を、ＭＭ空間に隣接する領域である表層ＱＭ領域と、表層ＱＭ領域以外の領域であるＱＭ領域とに分割する領域分割部と、表層ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を求め、ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を求め、表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換え、ＱＭ空間についての非経験的分子軌道法による全エネルギーを求める第１の演算部と、を有する。

このように本発明では、ＱＭ／ＭＭ法による分子シミュレーションにおいて、ＱＭ法が適用される空間のうちＭＭ空間に接続する部分を表層ＱＭ領域とし、この表層ＱＭ領域には経験的ポテンシャルが染み込むようにしている。本発明の手法を用いることにより、シミュレーション計算の過程において、タンパク質の活性部位の近傍の構造、すなわちＱＭ空間とＭＭ空間とが接続する領域において、タンパク質の構造が正しく維持される。このように本発明によれば、ＱＭ空間とＭＭ空間とに分割したことに起因する不整合の問題を回避でき、精度の高い分子シミュレーションを行えるようになる。

タンパク質−リガンド複合体におけるＱＭ−ＭＭ空間分割を概念的に示す図である。 ＱＭ空間とＭＭ空間との分割の一例を説明する図である。 本発明の実施の一形態の分子シミュレーション装置の構成を示すブロック図である。 ＱＭ空間とＭＭ空間との接合部分の分子配列を拡大して示した図である。 分子シミュレーションの手順を示すフローチャートである。

本発明に基づく方法は、ＱＭ／ＭＭ法によって分子シミュレーションを行う際に、ＱＭ法が適用される空間の一部において、非経験的分子軌道法における全エネルギー表式の一部を経験的ポテンシャルに置き換えることにより、従来のＱＭ／ＭＭ法において発生するＱＭ−ＭＭ境界領域の分子構造不安定性を回避するものである。具体的には、ＱＭ／ＭＭ法の適用を前提として、シミュレーション対象の分子をＱＭ法が適用されるＱＭ空間と、ＭＭ法が適用されるＭＭ空間とに分割し、さらに、ＱＭ空間のうち、ＭＭ空間に近接する部分を表層ＱＭ領域とする。そして、表層ＱＭ領域においては、その全エネルギーの一部を経験的ポテンシャルに置き換えている。経験的ポテンシャルにどの程度置き換えるかは、外部からのパラメータで調整可能とすることが好ましい。さらに、表層ＱＭ領域では、その表層ＱＭ領域の一部を構成する分子もしくは分子の一部の波動関数を、局在化分子軌道で表現することが好ましい。このようにして本実施形態では、ＱＭ空間の一部に経験的ポテンシャルを持ち込んでいる。

まず、本実施形態に基づく方法の理論的側面を詳細に説明する。

いま、ＱＭ空間とＭＭ空間が接合している表層領域に分布する分子もしくは分子の部分構造に、局在している局在化分子軌道を考える。他方、表層以外のＱＭ空間の原子に属する電子を記述する分子軌道は、その空間全体に拡がる正準軌道(Canonical Orbital)で与えられるとする。これらの軌道を用いた、Hartree-Fock（ハートリーフォック）法の全エネルギーＥ(ＱＭ＋ＰＱＭ)は、(2)式で与えられる。

ただし、

である。ここで、Ｑは正準軌道の空間を、Ｐは局在化分子軌道の空間を意味している。後述する説明から明らかになるように、Ｑが表す空間は、経験的ポテンシャルが持ち込まれずに純粋にＱＭ法によって処理される領域（ＱＭ領域）であり、Ｐが表す空間は、ＱＭ法による計算が行われるものの経験的ポテンシャルが持ち込まれる領域（表層ＱＭ領域）である。また、(2)式は、分子軌道法の分野における通常の表記を採用しており、したがって、φ_iは分子軌道であって、

と表される。Ｒ_ABは原子核ＡとＢとの距離、Ｚ_A，Ｚ_Bはそれぞれ原子核Ａ，Ｂの電荷である。Hartree-Fock法の代わりに、多配置ＳＣＦ（自己無撞着場；Self-consistent Field）法、配置間相互作用法等のポストHartree-Fock法への拡張は、局在化分子軌道による部分を、いわゆる“Frozen Core Orbital”と考えれば、同様の取り扱いにより可能である。(2)式の最初の３項は、ＱＭ領域に属する原子核と電子のみに関係し、次の３項は、ＱＭ領域と表層ＱＭ領域との相互作用に関する項であり、(2)式の最後の３項は、表層ＱＭ領域に属する原子核と電子のみに関係していることが分かる。そこで、表層ＱＭ領域のみに関わる項、すなわち(2)式の最後の３項に着目し、これらの項に関してパラメータαを用いて、(3)式のように分割する。ただし０≦α≦１である。

ここではαと（１−α）を用いて形式的に分割しているだけであるので、この段階では、全エネルギーの値は、(2)式と全く同等である。いま、(3)式の最後の項すなわち係数αが乗じられている項について、平衡構造の近傍で、原子核の座標の微小変化Δｑに関して級数展開（テーラー展開）すると、(4)式のように表現される。

もし、ＱＭ空間に属する原子が全て安定構造の位置にあるとすると、言い換えれば平衡点にあるとすると、展開項数を無限に取れば、(3)式と(4)式は、再び等しい全エネルギーを与えることになる。

ここで、ＱＭ空間の表層領域に、経験的ポテンシャルによるエネルギー表式を組み込むために、(4)式の最終項を経験的ポテンシャルの項で置き換えて、(5)式のように、全エネルギー表式を変更する。

ＭＭ空間に属する原子との相互作用は，通常の手続きにしたがって、表層ＱＭ領域の原子との間でのみ、結合距離ｒ、結合角Θ、２面角Ψを用いて記述すれば、全く過不足なく、滑らかにＭＭ空間のポテンシャルと連続的につながる。(5)式の最後の項は、経験的ポテンシャルの項であるから、パラメータαは、表層ＱＭ領域においてＱＭによる項とＭＭによる項とをどのように混ぜるか、言い換えれば、表層ＱＭ領域において経験的ポテンシャルをどの程度“染み込ませるか”を示している。パラメータαの好適な値は、分子系や、分子シミュレーションの目的によって変化し得るものであるが、例えば、０．２程度に設定される。

このようにして、ＱＭ空間とＭＭ空間を合わせた全エネルギーＥ(Total)は、(6)式で表される。

もちろん、Ｅ(ＭＭ)は、ＭＭ空間に対して分子動力学的方法などにより経験的ポテンシャルを適用して得られたエネルギーである。ここで、ＱＭ空間原子とＭＭ空間原子の相互作用を表現する式は、注目している物理的化学的性質により異なるが、ここでは標準的な、ＱＭ空間原子に対するＭＭ空間原子のつくる静電的効果と、ＱＭ空間の電子分布の変化によりＭＭ空間原子に誘起される双極子能率で記述される効果とを示してある。

以上の説明においては、表層ＱＭ領域が１つのみ設定されているものとしているが、１つの分子において複数の表層ＱＭ領域が設定されている場合においても上述の議論を拡張することができる。その場合には、表層ＱＭ領域の相互間でのクーロン力の項及びファンデルワールス力の項も考慮しなければならない。そこで、(2)式を、複数の表層ＱＭ領域が設定されている場合に拡張する。その場合、(2)式の最初の３項は、ＱＭ領域に属する原子核と電子のみに関する項であるから、変更を要しない。(2)式における次の３項は、ＱＭ領域と表層ＱＭ領域との相互作用に関する項であるから、各表層ＱＭ領域ごとに独立して求めてそれらを加算すればよい。(2)式の最後の３項は、上述と同様の手順にしたがって、(7)式のように変形することができる。

ここでは、表層ＱＭ領域がＰＱＭ個存在するものとし、各表層ＱＭ領域ごとに、その表層ＱＭ領域においてＱＭによる項とＭＭによる項とをどのように混ぜるかを示すパラメータα_mが定められている。(7)式中のＥ₀ ^mは、(5)式中のＥ^MM(0)に対応するものであって、ｍ番目の表層ＱＭ領域についての平衡状態での全ＭＭエネルギーを表している。Ｅ^mは、(8)式のように表される。

(8)式は、クーロン相互作用とファンデルワールス相互作用を含む経験的ポテンシャルによる式であり、最初の項は結合距離ｒによる寄与を示し、次の項は結合角Θによる寄与を示し、３番目の項は二面角Ψによる寄与を示し、４番目の項は表層ＱＭ領域の相互間でのファンデルワールス力の項であり、最後の項は表層ＱＭ領域の相互間でのクーロン力の項である。このように、クーロン力及びファンデルワールス力による項を考慮しなければならないことを除けば、表層ＱＭ領域が複数ある場合も、表層ＱＭ領域が１つしかない場合と同様の取り扱いをすることができる。

次に、表層ＱＭ領域の末端原子について述べる。末端の原子については、切断したままでは不対電子の問題が発生するので、共有結合している相手側の、現存する原子を水素類似原子と見なしたリンクアトム法を採用する。この時、基底関数としては、原子価電子の原子軌道のみを用い、全電子計算によるマリケン(Mulliken)電荷を再現するように、Orbital Exponentを事前に調整しておく必要がある。

次に、前述の理論に基づく本実施形態の分子シミュレーション方法及び装置について、具体的に説明する。図３は、本発明の実施の一形態の分子シミュレーション装置の構成を示すブロック図である。

この分子シミュレーション装置は、シミュレーション対象となる分子の構造データとして、その分子を構成する各原子の座標データなどを使用する。分子シミュレーション装置は、対象分子の構造データなどを格納する初期データ格納部１１と、初期データ格納部１１に格納されている構造データをＱＭ空間とＭＭ空間とに分割し、さらにＱＭ空間については、ＭＭ空間との接続部である表層ＱＭ領域とそうでないＱＭ領域とに分割する領域分割部１２と、ＭＭ空間について、分子動力学法などのＭＭ法に基づいて分子シミュレーション演算を実行するＭＭ演算部１３と、ＱＭ空間について非経験的分子軌道法による分子シミュレーション演算を実行するＱＭ演算部１４と、上述したパラメータα（０≦α≦１）が入力するパラメータ入力部１５と、ＭＭ演算部１３及びＱＭ演算部１４での演算結果を総合して出力する出力部１６と、を備えている。ここでＱＭ演算部１４は、上述した(5)式に示すように、ＱＭ法による全エネルギーＥ(ＱＭ＋ＰＱＭ)を算出する。すなわちＱＭ演算部１４は、表層ＱＭ領域については、局在化分子軌道を用い、パラメータαで指定される割合で、ＱＭ法によるエネルギー成分と経験的ポテンシャルによる成分と重ね合わせて、計算を実行する。なおＱＭ領域に対しては、ＱＭ演算部１４は、通常のＱＭ法と同様に、正準軌道（拡がった分子軌道）を用いて計算を実行する。なお、ＭＭ演算部１３は、経験的ポテンシャルの計算を主として行い、ＱＭ演算部１４は、主として２電子積分の計算を行うものであるから、それぞれ、最適なハードウェア構成が異なっている。本実施形態の分子シミュレーション装置は、計算機クラスタとして構成するのに適したものであるが、その場合であっても、ＭＭ演算部１３及びＱＭ演算部１４は、それぞれ別のハードウェア構成の計算機から構成されるようにすることが好ましい。

図４は、ＱＭ空間（ＱＭ領域＋表層ＱＭ領域）とＭＭ空間との分割の一例を示している。ＱＭ領域と表層ＱＭ領域との境界は、例えば、図示左から２番目のアミノ酸残基におけるＣ−Ｃ結合の位置に設定されている。また、表層ＱＭ領域とＭＭ空間との境界は、図示左から４番目のアミノ酸残基のＣ−Ｃ結合に設定されている。これらの境界の位置は、分子シミュレーションの対象となる分子やシミュレーションの目的によって適宜に定められるものであるが、計算の精度を高めるためには、一重結合（σ結合）の位置に境界を定めることが好ましい。

次に、前述した理論的側面に基づき、本実施形態での実際の計算手順について、図５のフローチャートを利用して説明する。

まず、ステップ１０１において、初期データ格納部１１は、シミュレーション対象の分子の構造データ（座標データ）を初期データ格納部１１から読み出し、その分子をＱＭ空間とＭＭ空間とに分割し、さらに、ＱＭ空間をＱＭ領域と表層ＱＭ領域とに分割し、ＭＭ空間に接合する表層ＱＭ領域に属する分子もしくは分子の一部を仮想分子として切出す。ＭＭ空間についての構造データはＭＭ演算部１３に送られ、ステップ１０２において、ＭＭ演算部１３は、ＭＭ空間について、分子動力学法などの、経験的ポテンシャルに基づく手法によって分子シミュレーションを実行する。ＱＭ空間（ＱＭ領域及び表層ＱＭ領域）の構造データは、ＱＭ演算部１４に送られる。なお、ＭＭ計算を行う際に、表層ＱＭ領域内の、ＭＭ空間に隣接する原子の座標データが必要であるから、その座標データもＭＭ演算部１３に送られる。同様に、ＱＭ演算部１４には、ＭＭ空間内の原子のうち、表層ＱＭ領域に隣接する原子の座標データも送られる。

構造データが送られてくると、ステップ１０３において、ＱＭ演算部１４は、ＱＭ演算、具体的にはHartree-Fock計算を開始する。このとき、上述したように分子をＱＭ領域、表層ＱＭ領域及びＭＭ空間に分割したことに伴い、分子もしくは分子の一部において共有結合が切られ、その結果、分子もしくは分子の一部の両端または一端に不対電子が発生する。ＱＭ演算部１４は、上述したリンクアトム法により、水素類似原子の導入によりそのような不対電子を隠蔽し、まず、表層ＱＭ領域に関して、Hartree-Fock計算を行うようにする。

ステップ１０３の結果、正準軌道が得られるから、ステップ１０４において、その正準軌道が局在化分子軌道に変換され、ステップ１０５において、ＭＭ空間の反対側に位置する水素類似原子の分子軌道の係数を無視して、分子軌道の再規格化が行なわれる。その後、ＱＭ演算部１４は、ステップ１０６において、表層ＱＭ領域のこれらの局在化分子軌道の内、ＭＭ空間寄りに局在化した分子軌道を適宜複数選択してLowdinの直交化[6]を行い、規格直交化された局在化分子軌道基底を得る。ここで直交化を行うのは、水素類似原子（リンクアトム）の導入などの影響により、直交性が保証されていないからである。

ＱＭ空間内の電子のうち、ステップ１０３〜１０６での処理対象となったもの以外の電子、すなわちＱＭ領域の電子については、ＱＭ演算部１４は、ステップ１０７において、非経験的分子軌道法における通常の手順によって、初期分子軌道の係数を求める。ＱＭ演算部１４は、ステップ１０６とステップ１０７で得た分子軌道を用いて、反対称化された波動関数を作り、上述した(5)式の全エネルギーＥ(ＱＭ＋ＰＱＭ)を求める。このとき、(5)式に含まれるパラメータαが必要となるが、このパラメータαは、パラメータ入力部１５に外部から入力して、パラメータ入力部１５からＱＭ演算部１４に与えられる。その後、ＱＭ演算部１４は、ステップ１０９において、変分法に基づく自己無撞着（ＳＣＦ）法により、局在化分子軌道以外の分子軌道の係数を最適化し、ＱＭ空間の全エネルギーを求める。

以上の処理によって、ＭＭ空間の全エネルギーとＱＭ空間の全エネルギーとが求められたことになるから、出力部１６は、ステップ１１０において、両方の空間の全エネルギーを(6)式に基づいて組み合わせることにより、シミュレーション対象の分子全体のエネルギーＥ(Total)を出力する。

分子の安定構造や化学反応トラジェクトリーを求めるのに必要な原子核に働く力は、(6)式を原子核の座標で微分することにより求められる。この時、分子軌道の係数に対する偏微分項が現れるが、正準軌道に対する偏微分項はエネルギー勾配法[7]により、重なり積分の偏微分に置き換えられる。これに対し局在化分子軌道の係数に対する偏微分項の計算は必要となるが、ＣＰＨＦ(Coupled Perturbed Hartree-Fock)法[8]により、厳密に計算することができる。次いで、新しい原子核の座標を用いて、表層ＱＭ領域の分子もしくは分子の一部について、局在化分子軌道を新たに求める。この時、前回の局在化分子軌道との重なりを求め、最も重なりの大きい軌道を選択すれば、当初選択された局在化分子軌道に最も近い軌道が自動的に求められる。この局在化分子軌道の計算、Hartree-Fock法による全エネルギーと分子軌道係数の決定、原子核に働く力による新しい原子核の座標の算出、を繰り返せば、安定構造や化学反応トラジェクトリーの計算が行える。

以上説明した分子シミュレーション装置は、典型的には、計算機クラスタによって実現されるものであり、計算機クラスタには、初期データ格納部１１、領域分割部１２、パラメータ入力部１５及び出力部１６として機能する制御用の計算機と、ＭＭ演算部１３として機能する計算機または計算機群と、ＱＭ演算部１４として機能する計算機または計算機群が設けられる。これらの計算機は、それぞれその計算機が果たすべき機能を実行するためのプログラムを読み込むことによって、制御用の計算機、ＭＭ演算用の計算機あるいはＱＭ演算用の計算機として機能することになる。そのようなプログラムは、磁気テープ（ＭＴ）あるいはＣＤ−ＲＯＭなどの記録媒体を介して、あるいはネットワークを介して、計算機に読み込まれる。

比較的小規模な分子を対象とした場合には、単一の計算機を用いて、上述した分子シミュレーションを実行することも可能である。その場合には、上述した手順で分子シミュレーションを実行するための計算機プログラムを、スーパーコンピュータやパーソナルコンピュータなどの計算機に読み込ませ、そのプログラムを実行させればよい。分子シミュレーションを行うためのプログラムは、磁気テープやＣＤ−ＲＯＭなどの記録媒体によって、あるいはネットワークを介して計算機に読み込まれることになる。

上述したプログラム、そのようなプログラムを格納した記録媒体、そのようなプログラムからなるプログラムプロダクツも、本発明の範疇に含まれる。

タンパク質は、生体中で様々な機能を実現しているが、酵素という観点で捉えると、化学関連産業全般に関わる重要な化学物質である。したがって、本発明に基づいて、タンパク質と基質の相互作用を、量子力学的手法に則って、高い信頼性でシミュレーションできることは、医薬品や機能性食品の製造、化学工業分野全般、さらには、環境保全物質の開発など、広範囲の産業に有効な研究手法を提供することになる。

Claims (13)

1. シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、前記ＱＭ空間に対して非経験的分子軌道法を適用し、前記ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション方法であって、
   記憶部から、前記シミュレーション対象の分子または分子の一部を構造データを取り出して前記ＱＭ空間及び前記ＭＭ空間に分割する段階と、
   前記ＱＭ空間に関する非経験的分子軌道法における全エネルギー表式の一部を経験的ポテンシャルで置き換える段階と、
   を有する、分子シミュレーション方法。
2. シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、前記ＱＭ空間に対して非経験的分子軌道法を適用し、前記ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション方法であって、
   記憶部から、前記シミュレーション対象の分子または分子の一部の構造データを取り出して前記ＱＭ空間及び前記ＭＭ空間に分割し、さらに、前記ＱＭ空間を、前記ＭＭ空間に隣接する領域である表層ＱＭ領域と、前記表層ＱＭ領域以外の領域であるＱＭ領域とに分割する段階と、
   前記表層ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を得る段階と、
   前記ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を得る段階と、
   前記表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換える段階と、
   前記ＱＭ空間についての非経験的分子軌道法による全エネルギーを求める段階と、
   を有する、分子シミュレーション方法。
3. パラメータを入力する段階をさらに有し、前記表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換える割合が前記パラメータによって指定される、請求項２に記載の方法。
4. 複数の前記表層ＱＭ領域が設定され、前記各表層ＱＭ領域ごとに前記パラメータが入力され、表層ＱＭ領域相互間のクーロン相互作用及びファンデルワールス相互作用を含む経験的ポテンシャルが用いられる、請求項３に記載の方法。
5. 前記表層ＱＭ領域における波動関数を局在化分子軌道で表現する段階をさらに有する、請求項２乃至４のいずれか１項に記載の方法。
6. 分子構造最適化及び時間発展計算を行う段階をさらに有し、前記分子構造最適化及び時間発展計算に際し前記局在化分子軌道を更新する、請求項５に記載の方法。
7. シミュレーション対象の分子または分子の一部をＱＭ空間とＭＭ空間とに分割し、前記ＱＭ空間に対して非経験的分子軌道法を適用し、前記ＭＭ空間に対しては経験的ポテンシャルに基づく方法を適用して分子シミュレーションを行う分子シミュレーション装置であって、
   前記シミュレーション対象の分子または分子の一部の構造データを格納する記憶部と、
   前記記憶部から、前記シミュレーション対象の分子または分子の一部の構造データを取り出して前記ＱＭ空間及び前記ＭＭ空間に分割し、さらに、前記ＱＭ空間を、前記ＭＭ空間に隣接する領域である表層ＱＭ領域と、前記表層ＱＭ領域以外の領域であるＱＭ領域とに分割する領域分割部と、
   前記表層ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を求め、前記ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を求め、前記表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換え、前記ＱＭ空間についての非経験的分子軌道法による全エネルギーを求める第１の演算部と、
   を有する、分子シミュレーション装置。
8. パラメータを入力するパラメータ入力部をさらに有し、前記第１の演算部が前記表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換える割合が前記パラメータによって指定される、請求項７に記載の装置。
9. 複数の前記表層ＱＭ領域が設定され、前記各表層ＱＭ領域ごとに前記パラメータが入力され、表層ＱＭ領域相互間のクーロン相互作用及びファンデルワールス相互作用を含む経験的ポテンシャルが用いられる、請求項８に記載の装置。
10. 前記ＭＭ空間に関して経験的ポテンシャルによる方法に基づいて全エネルギーを求める第２の演算部と、
    前記第１の演算部で求めた全エネルギーと前記第２の演算部で求めた全エネルギーとから前記シミュレーション対象の分子または分子の一部の全エネルギーを算出する出力部と、
    をさらに有する、請求項７乃至９のいずれか１項に記載の装置。
11. コンピュータを、
    シミュレーション対象の分子または分子の一部の構造データを格納する記憶部、
    前記記憶部から、前記シミュレーション対象の分子または分子の一部を構造データを取り出して前記ＱＭ空間及び前記ＭＭ空間に分割し、さらに、前記ＱＭ空間を、前記ＭＭ空間に隣接する領域である表層ＱＭ領域と、前記表層ＱＭ領域以外の領域であるＱＭ領域とに分割する領域分割部、
    前記表層ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を求め、前記ＱＭ領域に関して非経験的分子軌道法による全エネルギー表式を求め、前記表層ＱＭ領域の全エネルギー表式の一部を経験的ポテンシャルに基づく項に置き換え、非経験的分子軌道法によって前記ＱＭ空間についての全エネルギーを求める第１の演算部、
    として機能させるプログラム。
12. 前記コンピュータを、さらに、
    前記ＭＭ空間に関して経験的ポテンシャルによる方法に基づいて全エネルギーを求める第２の演算部、
    前記第１の演算部で求めた全エネルギーと前記第２の演算部で求めた全エネルギーとから前記シミュレーション対象の分子または分子の一部の全エネルギーを算出する出力部、
    として機能させる、請求項１１に記載のプログラム。
13. コンピュータが読み取り可能な記録媒体であって、請求項１１または１２に記載のプログラムを格納した記録媒体。