効率的な多体シミュレーション方法
内部座標ヤコビアンを平方的な（内部座標の数Ｎに対してＮ２のオーダである）コストで計算するための高速かつ汎用的に適用可能な方法であって、分子モデル化方法の劇的な高速化に使用できる方法。一実施形態において本発明は、ペア‐ポテンシャルエネルギー項の制約なしにデカルト・ヘシアンのねじりヤコビアンへの変換に有用であり、コンピュータのメモリを高効率で使用しそのような計算を実行するに有用である方法とアルゴリズムを提供する。そのような方法とアルゴリズムは、変換を、例えば分子系などの任意の多体系をモデル化するために使用される内部座標の数に関して平方的な計算時間で行い得る。関連の実施の形態において、本発明は、ペア‐ポテンシャルエネルギー項の制約なしに剛性行列を計算するために有用であり、かつコンピュータのメモリを高効率で使用しつつそのような計算を実行するために有用である方法およびアルゴリズムを提供する。
【発明の詳細な説明】
【技術分野】
【０００１】
本願は、Ｄａｎ Ｒｏｓｅｎｔｈａｌによる仮特許出願第６０／４７７，２３７号（２００３年６月９日出願）の出願日とＤａｎ Ｒｏｓｅｎｔｈａｌによる仮特許出願第６０／５５２，２２２号（２００４年３月１０日出願）の出願日との優先権の利益を享受し、この２つの出願は、共に本明細書にその全体が援用される。
【０００２】
本発明は、数値的方法の分野に関し、さらに詳しくは、例えば多体機械系（ｍｕｌｔｉｂｏｄｙ ｍｅｃｈａｎｉｃａｌ ｓｙｓｔｅｍ）、とくには分子のモデル化の用途において使用される機械系など、機械系の運動方程式の解法との関連において使用される数値的方法に関する。本明細書にて言及する参考文献はすべて、あらゆる目的において、言及によって本明細書に組み込まれたものとする。
【背景技術】
【０００３】
分子のモデル化は、分子の構造および動的な挙動など、分子のさまざまな側面をシミュレートするために使用できるいくつかの技法を含んでいる。分子のモデル化の典型的な用途には、酵素‐リガンドの結合、分子の拡散、反応経路、相転移、およびたんぱく質の折り畳みの研究などが挙げられる。生物科学ならびに薬品、高分子化合物、および化学品産業における研究者らは、複雑な分子における化学プロセスの性質を理解し、それによって新規な薬剤または物質を設計するために、そのような技法を使用し始めつつある。当然ながら、これらのツールが受け入れられるか否かは、実体の表現における結果の正確さ、モデル化できる分子系のサイズおよび複雑さ、ならびに解を得る際の速度など、いくつかの要因にもとづいている。
【０００４】
分子または他の機械系のシミュレーションは、一般に、微分方程式の組を解くべく数値的方法を使用して実行される。したがって、そのようなシミュレーションの速度は、部分的には、使用される数値的方法の速度によって制約され、数値的方法を使用して解かれる計算およびアルゴリズムの速度によって制約される。機械的な（例えば、分子の）シミュレーションの計算速度は、系内の自由度（「ＤＯＦ（ｄｅｇｒｅｅｓ ｏｆ ｆｒｅｅｄｏｍ）」）の数への依存の程度の点で特徴付けられ、一般にＤＯＦの数は、系内の物体（例えば、原子）の数に比例する。例えば、「オーダ（Ｎ）」のアルゴリズムの計算時間が、系内のＤＯＦの数に線形に比例する一方で、「オーダ（Ｎ平方）」（オーダ（Ｎ２））のアルゴリズムの計算時間は平方の依存性を有し（すなわち、ＤＯＦの数の二乗に比例して増加する）、「オーダ（Ｎ３）」のアルゴリズムの計算時間は立方の依存性を有し、以下同様である。系のサイズ（したがって、ＤＯＦの数）が大きくなるにつれ、より高いオーダ（Ｎ）の依存性を有するアルゴリズムは、計算コストに関して急激に認容できなくなるほどになる可能性がある。したがって、可能な限りオーダ（Ｎ）への依存性が小さいアルゴリズムを開発することが、強く望まれる。
【０００５】
分子のシミュレーションは、一般に、デカルト座標（いずれも絶対基準座標系に関して表現される系内の各原子または物体についてのｘ、ｙ、ｚ座標）または内部座標（典型的には、二面角またはねじり角座標）のいずれかを使用して実行されている。内部座標を使用する利点は、使用されるアルゴリズムの少なくともいくつかにおいて複雑さが増すことと引き換えに、系内のＤＯＦの数を大きく減らす（例えば、デカルト座標の場合の約１／６〜１／８の間まで）ことができる点にある（例えば、Ｒｏｓｅｎｔｈａｌの２００２年８月８日付の国際公開パンフレットＷ００２／０６１６６２号、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」、またはＲｏｓｅｎｔｈａｌの２００３年１月２３日付の米国特許出願公開第ＵＳ２００３／００１８４５５Ａ１号明細書、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」を参照）。
【０００６】
典型的な２種類の分子モデル化技法は、時間にわたって分子系の運動をシミュレートする分子動力学（「ＭＤ」）アプローチ、および系を構成している原子または物体をランダムに変化させることによって分子または分子系について種々の状態を生成し、連続する各状態のエネルギーを評価するモンテカルロ法である（例えば、Ｌｅａｃｈ． Ａ． Ｒ．のＭｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ‐Ｐｒｉｎｃｉｐｌｅｓ ａｎｄ Ａｐｐｌｉｃａｔｉｏｎｓ ２ｎｄ Ｅｄ．、２００１年、Ｄｏｒｓｅｔ Ｐｒｅｓｓ、Ｄｏｒｃｈｅｓｔｅｒ、Ｄｏｒｓｅｔを参照）。
【０００７】
これらの、ならびに他の分子および機械モデル化技法は、構造的なエネルギー関数（位置および／または運動エネルギー）の一次（勾配、またはヤコビアン）および／または二次（へシアン）の導関数行列を利用するように定式化されることが多い。これらの導関数は、デカルト座標においては容易に計算可能であるが、内部座標において第１原理から計算することはきわめて困難である。したがって、内部座標についてエネルギー（Ｅ）または力（Ｆ）の関数のヤコビアンまたはヘシアンを見出したいと望む場合には、一般に、デカルトのヤコビアンまたはヘシアンを対応する内部座標（例えば、ねじり角座標）のヤコビアンまたはヘシアンに変換する必要がある。
【０００８】
デカルトのヤコビアンの計算コストは、オーダ（Ｎ２）と同じく低いが、デカルトのヤコビアンを内部座標のヤコビアンに変換するために一般に適用可能な従来技術の方法は、Ｎについて立方の依存性を有している（すなわち、オーダ（Ｎ３））。したがって、大規模な系について作業を行なうとき、変換のための計算がボトルネックとなる可能性がある。分子のモデル化の文脈においてヤコビアンの計算を高速に（オーダＮ２まで低く）するため、いくつかの方法が説明されている（例えば、Ｇｉｂｒａｔ、１９９７、「Ｆａｓｔ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｆｉｒｓｔ ａｎｄ ｓｅｃｏｎｄ ｄｅｒｉｖａｔｉｖｅｓ ｏｆ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｙ ｏｆ ｏｌｉｇｏｍｅｒｉｃ ｐｒｏｔｅｉｎｓ ｗｉｔｈ ｒｅｓｐｅｃｔ ｔｏ ａ ｓｅｔ ｏｆ ｖａｒｉａｂｌｅｓ ｍｉｘｉｎｇ ｄｉｈｅｄｒａｌ ａｎｇｌｅｓ ａｎｄ ｒｉｇｉｄ ｂｏｄｙ ｖａｒｉａｂｌｅｓ」、Ｊ． Ｃｈｅｍ Ｐｈｙｓ ９４：１２３４‐１２５６；Ｎｏｇｕｔｉ ＆ Ｇｏ、１９８３、Ｊ Ｐｈｙｓ． Ｓｏｃ． Ｊａｐ． ５２：３６８５）が、それらは以下に詳しく述べるとおり、或る限界を抱えている。
【０００９】
力またはエネルギーのヤコビアンを計算するための多くの技法は、連鎖法則を使用しており、例えば力のヤコビアンを内部座標について
【００１０】
【数６】
と表現することができ、ここでデカルト座標がｒとして表わされ、内部座標がｑとして表わされている。
【００１１】
【数７】
は、典型的には、Ｆそのものがポテンシャルエネルギーの勾配として得られるときは常にデカルト・ヘシアンであり、そうでない場合、
【００１２】
【数８】
はデカルト力ヤコビアンと称される。
【００１３】
【数９】
は、ねじり変位勾配と称される。陽解行列乗算法（「ｅｘｐｌｉｃｉｔ ｍａｔｒｉｘ ｍｕｌｔｉｐｌｉｃａｔｉｏｎｍｅｔｈｏｄ」）（例えば、両者ともここでの言及によってあらゆる目的について本明細書に組み込まれたものとするが、Ｒｏｓｅｎｔｈａｌの２００２年８月８日の国際公開パンフレットＷＯ０２／０６１６６２号、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」；またはＲｏｓｅｎｔｈａｌの２００３年１月２３日の米国特許出願公開第ＵＳ２００３／００１８４５５Ａ１号、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」を参照されたい）においては、デカルト・ヤコビアンおよび変位勾配の両者が形成され、陽解行列乗算が実行される。大規模な分子系においてはデカルト・ヤコビアンを保存するために通常でない量のメモリが必要とされるため、一般に、ヤコビアンをブロックに部分的に形成するため、或る特別な規定がなされる。行列の乗算のコストは、原子の数について立方的（オーダ（Ｎ３））である。使用される力の場に応じ、デカルト・ヤコビアンの形成のコストは、平方的（オーダ（Ｎ２））でありうる。変位勾配の形成のコストも平方的であり、したがってこの方法は、最終的には、行列乗算の立方的なコストによって支配される。
【００１４】
数値的摂動にもとづく方法（例えば、Ｌｙｎｅｓｓ ＆ Ｍｏｌｅｒ、１９６７、「Ｎｕｍｅｒｉｃａｌ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ ｏｆ ａｎａｌｙｔｉｃ ｆｕｎｃｔｉｏｎｓ」、ＳＩＡＭ Ｊ Ｎｕｍｅｒ． Ａｎａｌ． ４：２０２‐２１０を参照）も、すでに説明されている。使用されるシミュレーション・エンジンが原子間力を計算できる場合、その原子間力は、一般化座標の陰関数：Ｆ＝Ｆ（ｒ（ｑ））として眺めることができる。したがって、
【００１５】
【数１０】
を、例えば差分を使用して数値的に形成することができる。この仕組みの全体コストは、ヤコビアンの各列が少なくとも１つの力の計算を必要とするため、立方的である。この数値的な仕組みの精度も、丸めおよび切り捨てによる影響ゆえ、苦しいものになりうる。得られた結果の精度を回復させるため、複雑化技法（例えば、Ｍａｒｔｉｎｓら、２０００年、「Ａｎ ａｕｔｏｍａｔｅｄ ｍｅｔｈｏｄ ｆｏｒ ｓｅｎｓｉｔｉｖｉｔｙ ａｎａｌｙｓｉｓ ｕｓｉｎｇ ｃｏｍｐｌｅｘ ｖａｒｉａｂｌｅｓ」、ＡＩＡＡ ２０００−０６８９を参照）を使用することができるが、これは力の場の全体を複雑な算術を使用すべく書き換えることを必要とし、多くの用途において現実的でない要件である。
【００１６】
プログラム強化技法の適用から生じるフォワード・モードの微分（ｆｏｒｗａｒｄ ｍｏｄｅ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ）にもとづく方法も、先行技術において記載されている（例えば、Ｂｉｓｃｈｏｆら、１９９２年、「ＡＤＩＦＯＲ：Ｇｅｎｅｒａｔｉｎｇ Ｄｅｒｉｖａｔｉｖｅ Ｃｏｄｅｓ ｆｒｏｍ Ｆｏｒｔｒａｎ Ｐｒｏｇｒａｍｓ」、Ｓｃｉｅｎｔｉｆｉｃ Ｐｒｏｇｒａｍｍｉｎｇ １：１１−２９）。これらの方法の適用は、Ａｄｉｆｏｒ（Ａｒｇｏｎｎｅ Ｎａｔｉｏｎａｌ ＬａｂｏｒａｔｏｒｙおよびＲｉｃｅ Ｕｎｉｖｅｒｓｉｔｙ）などのコード生成プログラム、または手動のいずれかによって、行なうことができる。そのような方法の適用は、渡された（ｐａｓｓｅｄ‐ｉｎ）データ・ベクトルｚについての行列ベクトル積
【００１７】
【数１１】
を計算できるプログラムをもたらす。そのようなプログラムの良好な使用は、変位勾配を形成し、フォワード・モードの微分ルーチンに渡すことにある。これは、依然として立方的な全体コストをもたらすが、数値的差分法よりも約３倍高速であって、完全な精度が手に入る。
【００１８】
Ｎｏｇｕｔｉ ＆ Ｇｏ，１９８３年、（Ｊ． Ｐｈｙｓ． Ｓｏｃ． Ｊａｐ． ５２：３６８５）の方法を脚色してなるＧｉｂｒａｔの方法（Ｇｉｂｒａｔ、１９９７年、「Ｆａｓｔ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｆｉｒｓｔ ａｎｄ ｓｅｃｏｎｄ ｄｅｒｉｖａｔｉｖｅｓ ｏｆ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｙ ｏｆ ｏｌｉｇｏｍｅｒｉｃ ｐｒｏｔｅｉｎｓ ｗｉｔｈ ｒｅｓｐｅｃｔ ｔｏ ａ ｓｅｔ ｏｆ ｖａｒｉａｂｌｅｓ ｍｉｘｉｎｇ ｄｉｈｅｄｒａｌ ａｎｇｌｅｓ ａｎｄ ｒｉｇｉｄ ｂｏｄｙ ｖａｒｉａｂｌｅｓ」、Ｊ． Ｃｈｅｍ Ｐｈｙｓ ９４：１２３４−１２５６）においては、ヤコビアンおよびねじりヘシアンが、平方的なコストで生成されるが、ポテンシャルの組として表現できるエネルギー関数に限定されている。すなわち、１‐４ねじりポテンシャル（ねじりポテンシャルは、４原子項の加法和として表現されるため）、１‐３結合角度ペア、およびポテンシャルの組として専ら表現することが不可能である他のエネルギー関数など、単純なエネルギー関数には適用不可能である。さらに、真の多体ポテンシャルとして表現され、陰解法モデルを利用した分子シミュレーションに使用されることが多い一般化ボーン・ソルベント・アクセシビリティ（Ｇｅｎｅｒａｌｉｚｅｄ Ｂｏｒｎ Ｓｏｌｖｅｎｔ Ａｃｃｅｓｓｉｂｉｌｉｔｙ：「ＧＢＳＡ」）解法モデルなど、分子のモデル化に使用され、ポテンシャルの組として専ら表現することが不可能である他のいくつかのツールへの適用も不可能である。
【００１９】
本発明は、従来技術における前記およびその他の制限を克服するものであり、内部座標（例えば、ねじり角）ヤコビアンを平方的な（オーダ（Ｎ２）の）コストで計算するための高速かつ汎用的に適用可能な方法であって、幾多の用途の中でも、とりわけ分子のモデル化方法を劇的に高速化するために使用できる方法を提供する。
【発明の開示】
【課題を解決するための手段】
【００２０】
一実施の形態において、本発明は、ペア‐ポテンシャルエネルギー項の制約なしにデカルト・ヘシアンをねじりヤコビアンに変換するために有用であり、かつコンピュータのメモリを高効率で使用しつつそのような計算を実行するために有用である方法およびアルゴリズムを提供する。そのような方法およびアルゴリズムは、例えば分子系などの任意の多体系をモデル化するために使用されている内部座標の数について平方的な計算時間で変換を行なうことができる。関連の実施の形態において、本発明は、ペア‐ポテンシャルエネルギー項の制約なしに剛性行列を計算するために有用であり、かつコンピュータのメモリを高効率で使用しつつそのような計算を実行するために有用である方法およびアルゴリズムを提供する。
【発明を実施するための最良の形態】
【００２１】
（Ｉ．定義）
そのようではないと示されない限り、「約（ａｂｏｕｔ）」とは、記載の値±３０％を意味する。
【００２２】
そのようではないと示されない限り、「物体（ｂｏｄｙ）」とは、分子の表現の構成要素の文脈において、単一の質量または分子のモデル化の目的のための幾何構造として取り扱われる当該表現の単位として定義される。したがって、物体は、分子の個々の原子、原子の集まり、または質量からなる他の抽象系の表現であってよい。「剛体（ｒｉｇｉｄ ｂｏｄｙ）」は、剛であるとしてモデル化される物体である（すなわち、力が加わっても変形しない）。
【００２３】
そのようでないと示されない限り、「高速演算子の実装（ｆａｓｔ ｏｐｅｒａｔｏｒ ｉｍｐｌｅｍｅｎｔａｔｉｏｎ）」とは、積行列の要素が計算される行列の乗算の評価の方法であって、積行列の全体を計算および保存する必要なく「オンザフライ（ｏｎ‐ｔｈｅ‐ｆｌｙ）」で使用または解放される方法を意味する。したがって、高速オーダの実装は、標準的な行列の乗算（オーダ（Ｎ３））に比べ、実行速度および保存用メモリの両者の要件に関して大幅に効率的でありうる（オーダ（Ｎ２）と高い）。
【００２４】
そのようでないと示されない限り、「ジョイント（ｊｏｉｎｔ）」は、２つの物体間の接続である。一般的なジョイントの種類としては、ピン・ジョイント、スライダ・ジョイント、およびボール・ジョイントが挙げられるが、さらにこれらに限られるわけではないが、フリー・ジョイント、Ｕジョイント、円柱ジョイント、ベアリング・ジョイント、および前記いずれかの組み合わせなど、他にも多数の種類のジョイントが可能である。例えば、フリー・ジョイントは、直交する３つのスライダ・ジョイントをボール・ジョイントに組み合わせて構成され、完全な６自由度を有している。ジョイントの種類およびそれらの仕様についてのさらに詳細な検討を、例えば上記Ｒｏｓｅｎｔｈａｌに見つけることができる。
【００２５】
そのようでないと示されない限り、特定のジョイントまたはピボットにおける「凍結体（ｆｒｏｚｅｎ ｂｏｄｙ」」は、当該ピボットの外にあって剛体として一体に運動するすべての物体で構成される。
【００２６】
そのようでないと示されない限り、「計算実現可能なオーダ（ｃｏｍｐｕｔａｔｉｏｎａｌｌｙ‐ｆｅａｓｉｂｌｅ ｏｒｄｅｒ）」とは、タスクの特定の並びを最終的な結果を変化させることなく実行できる任意のオーダを指す。このコンセプトは、いくつかの方法においては特定のステップのオーダが、それらステップが実行され、かつ結果がそれらステップが当初提示されたオーダで実行された場合と実質的に同じである限りにおいて、重要ではないことから惹起されている。
【００２７】
そのようでないと示されない限り、用語「シリコン内（ｉｎ ｓｉｌｉｃｏ）」は、コンピュータを使用して実行される任意の方法またはプロセスを指す。
【００２８】
そのようでないと示されない限り、「分子（ｍｏｌｅｃｕｌｅ）」は、化学的結合によって接続された２つ以上の原子で形成される任意の巨視的構造である。分子の例としては、これらに限られるわけではないが、たんぱく質（例えば、抗体、受容体、など）、ペプチド、脂質、核酸（例えば、天然または合成のＤＮＡ、ＲＮＡ、ｇＤＮＡ、ｃＤＮＡ、ｍＲＮＡ、ｔＲＮＡ、など）、レクチン、糖（例えば、レクチン‐糖の複合体を形成している）、糖たんぱく、低分子、有機化合物、塩、金属、などの単原子または多原子の構造が挙げられる。
【００２９】
そのようでないと示されない限り、「表現（ｒｅｐｒｅｓｅｎｔａｔｉｏｎ）」とは、分子または分子構造に適用されるとき、例えばコンピュータ・シミュレーションにおいて機械によって使用するための分子または分子構造の抽象的説明を指す。例えば、分子の１つの表現は、分子を構成している原子または物体、あるいはそれらの抽象的代理の位置を集合的に定める座標の組である。分子の表現の例には、これらに限られるわけではないが、Ｘ線座標、「ｐｄｂ」（たんぱく質データバンク）ファイル、などが挙げられる。
【００３０】
そのようでないと示されない限り、「大きく（ｓｉｆｇｎｉｆｉｃａｎｔ）」という用語は、「大きく異なる」、「大きく抑制する」、または「大きく促す」に関して使用されるとき、比較される２つのグループ間の定量化可能なパラメータにおける差異であって、標準的な統計的検定を使用して統計的に有意である差異を指す。
【００３１】
そのようでないと示されない限り、「溶媒（ｓｏｌｖｅｎｔ）」は、溶質分子を含みうる任意の媒体を指す。溶媒の例としては、これらに限られるわけではないが、水および他の水性溶媒、ならびに有機溶媒（例えば、ＤＭＳＯ、脂質、アルコール、など）が挙げられる。溶媒は、一様であっても、非一様であってもよく、固体状、液体状、または気体状であってよい。
【００３２】
そのようでないと示されない限り、「ターゲット分子（ｔａｒｇｅｔ ｍｏｌｅｃｕｌｅ）」とは、分子シミュレーションの対象である主たる分子である。シミュレーションは、例えば互いに相互作用する２つ以上のたんぱく質のシミュレーションなど、２つ以上のターゲット分子を有することができる。
【００３３】
そのようでないと示されない限り、「汎用的に適用可能な方法（ｕｎｉｖｅｒｓａｌｌｙ‐ａｐｐｌｉｃａｂｌｅ ｍｅｔｈｏｄ）」とは、ポテンシャルの組における使用に限られず、分子のシミュレーションにおける使用に適した任意の方法である。
【００３４】
（ＩＩ．多体系（ＭＵＬＴＩ‐ＢＯＤＹ ＳＹＳＴＥＭＳ））
一実施の形態において、本発明の方法は、機械系の多体系（ＭＢＳ）の実装に関して使用される（例えば、ＭＢＳ分子シミュレーション；前記Ｒｏｓｅｎｔｈａｌを参照）。ＭＢＳの基本的な抽象化は、ジョイントで接続された剛体の集まりの抽象化である。物体のうちのベースと呼ばれる１つが、その運動が直接地面に参照付けられる特別な地位を有する。系のグラフは、ループを含みうる。ループ切断アルゴリズムが、循環するグラフをツリーおよび切断の組に還元する。切断は、ジョイントまたは物体において生じうる。これにより、すべてのジョイントがツリー系に置かれる。ループ切断アルゴリズムによって切断された物体は、溶接ジョイントを課すことによって回復される。このジョイントそのものは、６つの距離の制約の組へと分解される。溶接ジョイントが、元の物体の２つの部分を再組み立てする。
【００３５】
ツリーの重要な特性は、或る任意の物体から他の任意の物体への経路が一意である点にある。ツリー内の物体の数はｎであり、通常は、ベース物体にラベル１が割り当てられる。これは、ベース物体から任意の葉物体（ｌｅａｆ ｂｏｄｙ）へと到るあらゆる経路において、物体のラベルが減少することがないということを意味する。葉物体は、ただ１つの物体に接続された物体である。規則的なラベル付けを、葉物体（少なくとも１つ存在しなければならない）のうちの１つにラベルｎを割り当てることによって達成できる。この物体がグラフから除かれるならば、ｎ−１個の物体を有するツリーが残される。このツリーの葉物体のうちの１つにラベルｎ−１を割り当て、このプロセスを、すべての物体がラベル付けされるまで繰り返す。一般的な機械系のＭＢＳの実装の詳細は、すでに説明されている（例えば、Ｒｅｃｕｒｓｉｖｅ Ｍａｓｓ Ｍａｔｒｉｘ Ｆａｃｔｏｒｉｚａｔｉｏｎ ａｎｄ Ｉｎｖｅｒｓｉｏｎ、ＪＰＬ Ｐｕｂｌｉｃａｔｉｏｎ ８８‐１１、１９８８におけるＲｏｄｒｉｇｕｅｚおよびＫｒｅｕｔｚを参照）。分子系のＭＢＳの実装の詳細、ならびに関連の表記法も、やはり説明されている（例えば、前記Ｒｏｓｅｎｔｈａｌを参照）。
【００３６】
（ＩＩＩ．残差の計算）
典型的な実施の形態において、本発明の方法は、
【００３７】
【数１２】
の形式の運動方程式を有する多体系として表現された系（例えば、分子系）のシミュレーションに適用される。この方程式の右手側は残差ρであり、慣性力および力の場からのアクティブな力の寄与を含んでいる。残差のヤコビアンは、
【００３８】
【数１３】
である。本発明の方法は、アクティブな力の成分を表現している力の場からの残差への寄与を計算するために適用可能である。多体系は、原子間力の計算をトリガして、残差ベクトルの要素へと集めることができるデータ処理方法の集まりを具現化している。これらのデータ処理ステップは、特定の行列操作の点で高いレベルの表現を有することができるが、実際には、これらの操作の実際の実施は、典型的には、実行時間が問題の大きさと線形に比例するアルゴリズムのかたちである。それらは、いわゆるＯ（ｎ）法であって、先行技術に教示されている（例えば、前記Ｒｏｓｅｎｔｈａｌを参照）。しかしながら、本発明の一態様によれば、これらＯ（ｎ）法を、ヤコビアンの形成ステップに適用することができる。
【００３９】
例えば、残差は、演算子方程式
【００４０】
【数１４】
から形成される。
【００４１】
Ｆは、典型的には、３ｎａ×１のベクトルである。各エントリは、好ましくは、所与の原子に作用する全原子間力の大域的なデカルト成分で構成される。行列演算子Ｈ、Φ、およびＰは、Ｒｅｃｕｒｓｉｖｅ Ｍａｓｓ Ｍａｔｒｉｘ Ｆａｃｔｏｒｉｚａｔｉｏｎ ａｎｄ Ｉｎｖｅｒｓｉｏｎ、ＪＰＬ Ｐｕｂｌｉｃａｔｉｏｎ ８８‐１１、１９８８においてＲｏｄｒｉｇｕｅｚおよびＫｒｅｕｔｚによって記載されている。力分布行列Ｐは、多体系のピボットにおいて原子間力から作用している多体空間荷重を生成する。したがって、
【００４２】
【数１５】
と書くことができ、ここでＴは、典型的には６ｎｂ×１のベクトルである。特定の原子に作用する所与の力が、当該原子が位置する物体のピボット点に作用する力およびトルクへと写像される。このように、行列Ｐの所与の行が、特定の物体に作用する空間荷重に対応する。所与の行は、当該特定の物体上に位置する原子に対応するエントリについてのみ、ゼロでないブロックを有する。分子系の多体定式化の多くにおいて、各物体は、典型的には、約３つの原子を含み、あるいは表現している。したがって、行列Ｐは大きくかつスパースであり、明示的に形成される必要はない（例えば、前記Ｒｏｓｅｎｔｈａｌを参照）。むしろ、行列ベクトル積ＰＦを計算するアルゴリズムを使用することができる。そのような計算の実効時間は、原子の数に線形である（前記Ｒｏｓｅｎｔｈａｌ）。
【００４３】
多体の力学においては、さらに演算子の転置が有用な量であることが一般的である。この場合、演算子ＰＴが、空間速度（直線および角）を、各物体ピボットから、当該物体に位置する原子位置へと伝播させる。この行列の典型的な用途は、分子系内の原子の微分空間変位の計算の目的のためにある。
【００４４】
行列Φは、多体遷移演算子として知られている。この行列は、データ・ベクトルに作用し、ツリーのベース物体の下流の葉物体からのデータ・ベクトルの要素のシフトおよび累積の効果を有している。例えば、積ΦＴが、反応荷重の６ｎｂ×１のベクトルＲを生成する。所与の物体について、反応荷重の要素は、そのジョイントを中心とする当該物体の運動によって運動力学的に影響される物体に作用するそれらの力の物体内ピボットについての結果としての力およびトルクを表わしている。ΦＴの計算は、物体の数について線形である。
【００４５】
Φを理解するために、最初に行列εφを導入することが有用である。この（ブロック）行列は、スパースであってスーパー・ダイアゴナルである。ｉ番目の行において、物体ｉの子に対応する列のみが満たされている。図１に示した物体の図について、行列εφが図２に示されている。εφのゼロでない各要素は、６×６の行列である。要素ｉφｋが、空間荷重を子のピボットからその親のピボットへとシフトする。転置された要素は、空間速度を親ピボットから子ピボットへとシフトする。
【００４６】
今や、Φとεφとの間の関係を、
【００４７】
【数１６】
と記述することができる。
【００４８】
この方程式の意味は、図１に示した系に関して、各ピボットに加えられる空間荷重を考えて各ピボットにおける空間反応荷重の計算を説明する以下の説明を参照することによって、よりよく理解できるであろう。反応荷重は、以下のＯ（ｎ）スイープ（コードのシークエンス）によって計算することができる。
【００４９】
【数１７】
これらの方程式は、各ピボットにおける反応を累積させる。スイープは、葉物体から出発して、ベース物体へとツリーをさかのぼる。各物体が、自身の反応を自身の親へとシフトさせ、そこで兄弟からの寄与に加算される。これらの方程式が、先に計算された量への後方参照を行なうことに注意すべきである。例えば、Ｒ（６）はＲ（７）を参照し、Ｒ（７）がＲ（８）およびＲ（１２）を参照し、以下同様である。しかしながら、このアルゴリズムが提示のとおりに構成される場合、不満足な参照に遭遇することはない。演算子方程式Ｒ＝ΦＴと上記手順との間の対応は、以下のとおり説明できる。すなわち、この手順が、Ｒ＝εφＲ＋Ｔを示しており、ここから（Ｉ−εφ）Ｒ＝Ｔ、またはＲ＝（Ｉ−εφ）‐１Ｔ＝ΦＴである。他方の方向（方程式からコード・シークエンス）に進むためには、これらのステップを単純に反転できる。
【００５０】
【数１８】
Φの要素を露出したいと望む場合には、例えばＲの計算において繰り返しの置換を実行することによって後方参照を除き、次いでΦの行列要素を特定することができる。あるいは、アルゴリズム的な代替案を使用することも可能である。（Ｉ−εφ）Φ＝Ｉであるため、Φ＝εφΦ＋Ｉが得られ、この方程式は、Ｒ＝εφＲ＋Ｔと同じ構造を有している。Φの列に同じアルゴリズムを、最後の列から出発して最初へと後ろ向きに作業して適用できる。しかしながら、明らかであろうが、典型的には、演算子の作用のみが必要とされるため、Φの要素が必要になる状況が存在しない。したがって、Φ（密な上三角行列）による乗算を示すすべてのアルゴリズムが、スパースなスーパー・ダイアゴナル行列εφに関して進行できる。残差の計算を終えるため、行列ベクトル積ＨＲが、ベクトルＲの要素をジョイントの自由に投影する。行列Ｈは、ブロック対角である。各ブロックは、特定のジョイントについてのジョイント・マップで満たされている。ピン・ジョイントについてのマップは、１×６のベクトル［λ ０］であり、単位ベクトルλがピンの幾何学的性質を指定している。行列ベクトル積ＨＲは、非再帰的なドット積のシークエンスを代表するため、線形なコストで容易に計算される。演算子Ｈは、以下の形状を有しており、ここでｎは物体の数であり、ｎｕは一般化速度の数である。
【００５１】
【数１９】
以下の形状を有する対応する演算子
【００５２】
【数２０】
が存在し、ここでｎｑは一般化座標の数である。
【００５３】
【数２１】
前記に照らし、残差の計算を力の場によって計算された原子間力の要素へと適用されるＯ（ｎ）法のシークエンスとして見ることができることを、理解できるであろう。これらは、例えば、真空の項、溶媒の項（極性または非極性）、および速度関連の項を含む。これらの項目のいくつかは、平方的な計算時間を必要とする。例えば、多くの実装例において、静電およびＧＢＳＡ溶媒モデルは、一般に、計算時間の大部分の力の計算に消費する。したがって、残差コストは、原子の数について平方的である。これは、残差のヤコビアンが、最良でも平方的なコストであり、おそらくはより悪いことを意味する。
【００５４】
ρ＝ＨΦＰＦを考えると、
【００５５】
【数２２】
を
【００５６】
【数２３】
と表現することができる。
【００５７】
原子間力ベクトルＦは、第１および第２項の両者に現れる。しかしながら、第２項は、力ベクトルの解析的形態についての知識を必要とせず、力ベクトルの数値のみを必要とする。したがって、この項は、この項の計算の際に新しい力の評価は不要であるため、元のコストがＯ（ｎ）である項の微分を呈する。したがって、この項がヤコビアンにおいて平方的なコストのみを引き起こしうるということを実証できる。これは、例えばＢｉｃｈｏｆら、１９９４、（「Ａｕｔｏｍａｔｉｃ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ：ｏｂｔａｉｎｉｎｇ ｆａｓｔ ａｎｄ ｒｅｌｉａｂｌｅ ｄｅｒｉｖａｔｉｖｅｓ−ｆａｓｔ」Ｐｒｏｃ． ＳＩＡＭ Ｓｙｍｐｏｓｉｕｍ ｏｎ Ｃｏｎｔｒｏｌ Ｐｒｏｂｌｅｍｓ ｉｎ Ｉｎｄｕｓｔｒｙ Ａｒｇｏｎｎｅ、Ｐｒｅｐｒｉｎｔ ＭＣＳ−Ｐ４８４−１１９４）に記載のとおり、残差をフォワード・モードの微分を使用し、力ベクトルを既知の定数として取り扱って微分でき、これによりコストが（コール当たりの）元のプログラムのコストの定数倍によって拘束されるプログラムがもたらされることに注目することによって理解できる。定数は、プログラム中に現れる固有の関数の性質によって決まる。多体系は、算術的な操作に加え、三角法のただ１つの変関数のみを含んでいる。すなわち、一般的に、定数は、典型的には、元のプログラムのコストの１〜２倍のオーダ、またはそれ以下であり、第２項は、ヤコビアン全体について最大でも平方的なコストで容易に計算できる。
【００５８】
【数２４】
を含んでいるヤコビアンの第１項の計算は、以下で説明される。
【００５９】
【数２５】
本発明の一態様によれば、変位勾配を、Ｏ（ｎ）多体演算子Ｐ、Φ、Ｈに関して表現できる。
【００６０】
【数２６】
とくに多体系の実施に適用可能な一実施の形態においては、これらの演算子が、平方的なコストで変位勾配を形成するために使用される。これらの演算子を使用しない場合、コストは立方的であろう。このアプローチは、変位勾配が明示的には必要とされず、デカルト・ヤコビアンによる行列積のみが必要とされるという事実を利用している。積は、
【００６１】
【数２７】
であるため、Ｏ（ｎ）演算子に関して形成できる。
【００６２】
これは、残差の計算に関してすでに説明したものと同じ演算子である。したがって、この方法は、デカルト・ヤコビアンが立方的なコストで計算可能ならば、内部座標（例えば、ねじり）ヤコビアンを、列について線形なコストで、すなわち平方的な全体コストでもたらす。
【００６３】
上記関係を組み合わせて、
【００６４】
【数２８】
をもたらすことができ、ここで行列
【００６５】
【数２９】
である。行列Ｗが多体系において形成される場合、この行列の一般要素Ｗｉｊは、空間における多体系の各物体の配置を反映する。この要素が、物体ｊのピボットの空間的変位ゆえに物体ｉに作用する空間荷重の微分の計算に使用される。Ｗの要素の計算は、デカルト・ヘシアンへの接続がなされる所であり、方程式
【００６６】
【数３０】
を有している。
【００６７】
このアプローチのいくつかの特徴は、注目を有している。すなわち、（ｉ）この方程式がヘシアン要素を使用する。それは、物体ｊの原子リストに対する物体ｉの原子リストによって案内される。また、（ｉｉ）各物体が、典型的にはおおまかに３つの原子で構成されるため、これは単純な計算である。（ｉｉｉ）各原子が、ただ１つの物体に属するため、デカルト・ヘシアンの各要素は、ただ１つの要素Ｗｉｊにのみ寄与する。さらに、従来技術の方法（例えば、Ｇｉｂｒａｔ Ｊ． Ｆ．（１９９７）「Ｆａｓｔ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｔｈｅ ｆｉｒｓｔ ａｎｄ ｓｅｃｏｎｄ ｄｅｒｉｖａｔｉｖｅｓ ｏｆ ｔｈｅ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｙ ｏｆ ｏｌｉｇｏｍｅｒｉｃ ｐｒｏｔｅｉｎｓ ｗｉｔｈ ｒｅｓｐｅｃｔ ｔｏ ａ ｓｅｔ ｏｆ ｖａｒｉａｂｌｅｓ ｍｉｘｉｎｇ ｄｉｈｅｄｒａｌ ａｎｇｌｅｓ ａｎｄ ｒｉｇｉｄ‐ｂｏｄｙ ｖａｒｉａｂｌｅｓ」、Ｊ． Ｃｈｅｍ Ｐｈｙｓ ９４：１２３４‐１２５６）と異なり、上述の方法は、ヘシアン要素を直接使用する。したがって、ポテンシャルの組について制約はない。Ｇｉｂｒａｔ（前記）による論文中の対応する量Ｌｉｊは、物体ｊの原子に対する物体ｉの原子から直接的に構成されるが、一般には、Ｌｉｊ≠Ｗｉｊである。本明細書に記載するようなヘシアン要素の使用は、より単純な再帰シークエンスにつながる。
【００６８】
高速演算子の実装においては、全体の計算が、２つの構成されたダウンループ（ｄｏｗｎ‐ｌｏｏｐ）を呈する。項
【００６９】
【数３１】
が、以下の再帰的なアルゴリズムによって計算される。まず、Ｙ＝ΦＷ、Ｚ＝ＨＹとする。これは、
【００７０】
【数３２】
をもたらす。
【００７１】
すでに述べたとおり、行列εφはサブダイアゴナルであり、要素ｉφｋ（ｋ）で占められている。したがって、「積層された」方程式を見つけることができる。
【００７２】
【数３３】
これらの方程式において、Ｗ１は最初の（論理的な）行を表わしており、ＷｎはＷの最後の（論理的な）行を表わしている。これは、６×６ｎｂのブロックである。Ｗの各行は、先の行に関して計算される。Ｈがブロック・ダイアゴナルな構造を有しているため、Ｚ＝ＨＹの計算は瑣末である。
【００７３】
今や、第２の再帰を、以下のとおり表現できる。
【００７４】
【数３４】
演算子ＨΦが、Ｚの行に作用していると見ることができる。これは、Ｚを行のオーダで生成した第１の再帰と好都合に接合する。第２の再帰は、目的の行列が対称であるため、対角線までの所与の列の要素の計算のみを必要とする。
【００７５】
以上に照らし、原子間力ヤコビアンを収縮させてＷを形成できることを、理解できるであろう。全体アルゴリズムは、力のヤコビアンへの行（または、列アクセス）のみを必要とし、同時に行列全体にアクセスする必要はない。さらに、所望であれば第１および第２の再帰を散布できることも、理解できるであろう。このようにして、ＹまたはＺの要素を、それらが計算されるや否や解放することができる。したがって、１つのステップでＹまたはＺのすべてを形成する必要はない。
【００７６】
（Ｖ．メモリに関して効率的な計算）
前記の方法は、Ｗの対称性をさらに充分に利用することによって、計算を実行するために必要となるメモリの量を最小化するやり方で、実装することができる。メモリの必要性を少なくするための１つのアプローチは、デカルト・ヘシアンを大きなブロックで生成し、次いで各ブロックを処理する「スライス化」技法を使用することにある。このような方法は、メモリの使用を少なくするが実行時間が苦しくなりうるという意味で、実行可能であろう。以下で説明する方法は、実行速度を犠牲にすることなくメモリの必要性を大きく減らすことができる処理の仕組みをもたらすように設計されている。
【００７７】
一実施の形態において、この方法は、ねじりヤコビアンを
【００７８】
【数３５】
と表現することにもとづいている。
【００７９】
新しい行列Ωは、加えられた荷重Ｔに対して反応Ｒが有している関係と同じ関係を、Ｗに対して有している。行列Ｗは、多体系の物体についての空間へシアンである、空間に浮遊するが、ジョイントによって接続されていない物体について関連性があるのは、ヘシアンである。行列Ωの要素が、凍結体の「組」を結び付ける。各ピボットにおいて、凍結体は、ピボットの外のすべての物体であって剛体として一体に運動する物体で構成されている。
【００８０】
各凍結体のピボットにおける微分空間変位が、他のすべての凍結体のピボットにおける微分空間荷重を生じさせる。変位させられた物体を「ソース」物体として考え、荷重を受けている物体を「受け手」と考えると、各ソース物体が各受け手物体と相互作用することを理解できるであろう。しかしながら、ソース物体における空間的移動は、凍結体のピボットを起源とし、ΦＴの作用を通じてソース系内の他の物体に伝播する。各ソース物体が、結合行列Ｗを通じて受け手系の各物体上に微分空間荷重を生む。この空間荷重は、ＷΦＴによって与えられる。受け手系の空間荷重は、（ＷΦＴ）に作用する演算子Φの作用によって受け手ピボットへと総計される。
【００８１】
高速演算子の実装において、Ωは以下のとおり効率的に計算できる。
【００８２】
【数３６】
最後の方程式が、Ωの計算のための基礎をもたらしている。εφがスーパー・ダイアゴナルである（主たる対角線の上方のみがゼロでない）ため、項εφΩは、Ωの低い番号の行をより高い番号の行に結合させているが、その逆は真ではない。すなわち、所与の列において、要素は、Ωがεφで前もって乗算されるとき、列を上方へ移動する。同様に、項ΩεφＴは、低い番号の列をより高い番号の列に結合させるだけである。要素は、所与の行を横切って左方に移動する。第３の項は、初期の行および初期の列の両者を、より高い行およびより高い列に結合させているが、結合は、他の方向においては生じない。要素は、上、左、または上および左へと移動する。正味の結果は、ツリー系が規則的にラベル付けされている限り、Ωが常に、逆の列のオーダで計算されるというものである。
【００８３】
望ましいラベル付けは、εφの帯域幅を最小にするラベル付けである。εφの帯域幅は、行列の「リーチ」を制限し、Ωの計算において後方への参照がどの程度まで「後ろ」に延びることができるのかを決定する。
【００８４】
一例として、上記アルゴリズムを、以下の方程式を計算することによって図１に示した系の部分集合に適用できる。
【００８５】
【数３７】
最初の２つの方程式は、結合行列Ｗを通じての物体の直接的寄与のみを取り上げている。（１２、１０）の要素が、凍結体１０の空間変位による凍結体１２への空間荷重を計算する。凍結体１０は、実際の物体１０および１１で構成されている。１０のピボットにおける変位が、ｉφｋ＊（１１）を通じて１１に伝播する。この変位が、結合要素Ω（１２、１１）を通じて１２への空間荷重を生む。（１２、７）の要素は、さらにわずかに複雑である。凍結体７は、物体７ならびに凍結体８および１２で構成される。これらは、物体７の子に相当する凍結体である。物体７のピボットにおける変位が、その子へと伝播し、Ωのすでに計算された要素に結合し、物体１２のピボットへの荷重を生成する。
【００８６】
再帰なしで、子の凍結体を構成している物体の集まり全体にわたって計算が必要であろう。これは、計算が進むにつれて成長傾向のリストとなるであろう。再帰ありでは、子の凍結体のみを訪れるだけでよい。さらに、各物体が実際には少数の子しか有していない（通常は２つまたは３つであるが、より多い場合もある）ため、要素当たりのコストは比較的小さい。
【００８７】
以下で説明され、図３Ａ〜３Ｆに示されるアプローチは、２つの凍結体間の相互作用の計算を示している。説明されている操作は、実際には、前記方程式の計算の手順において「自動的」に実行されるが、これらの方法のよりよい理解および実施を可能にするために示されている。検討対象の物体およびそれらの子の追跡を容易にするため、下記の表１のような表が作成される。
【００８８】
【化１】
図３Ａに目を向けると、この図は、凍結体３１２と凍結体３１４との間の相互作用３１０を示している。この例ａ≧ｂにおいて、凍結体３１２（「受け手物体」）は、物体ａ、その子ｋ、およびａの１つの孫を含んでおり、凍結体３１４（「ソース物体」）は、物体ｂ、その子ｒ、およびｂの２つの孫を含んでいる。相互作用３１０は、任意の計算実現可能なオーダにて実行される図３Ｂ〜３Ｄに示すステップによる分解によって決定できる。例えば、最初に、２つの個々の物体ａおよびｂの間の直接項３１６が、図３Ｂに示すとおり決定される。次いで、図３Ｃに示すとおり、ｂ３２２の凍結した子と相互作用する凍結体ａ３２０によって生成された項３１８（ΩεφＴから得られる）が、加算される。必要となる最後の相互作用（３２４）は、図３Ｄに示すとおり、物体ｂ３２６とａの凍結した子であるｋ（３２８）との間の相互作用である。この項は、好ましくは、物体ｂが単独では凍結体ではないため、間接的に計算される。一実施の形態においては、ａの凍結した子に対する物体ｂの相互作用は、減算によって計算される。すなわち、最初に、図３Ｅに示されるとおり、凍結体ｂ３１４とａの凍結した子３２８との間の相互作用３３０が計算され、εφΩから得られる項を生じさせる。次いで、図３Ｆに示すとおり、ａの凍結した子とｂの凍結した子との間の相互作用３２２（項εφΩεφＴ）が、減算される。これによって、図３Ｄに示すとおりａの凍結した子に対するｂの相互作用がもたらされる。
【００８９】
一例として、前記操作を、図１に示した多体系の一部分の解析へと適用することができる。第１に、両方の物体が子を有している任意の要素が一般的な場合である。そのような要素の例は、Ω（８、６）である。
【００９０】
【化２】
【００９１】
【数３８】
対角要素は、非対角の要素と同じ方法によって処理される。例は、Ω（５、５）について与えられる。
【００９２】
【化３】
【００９３】
【数３９】
当然ながら、行列の対称性ゆえ、Ω（５、６）＝Ω（６、５）Ｔであることを理解すべきである。要約すると、各要素が、式
【００９４】
【数４０】
を使用して計算される。
【００９５】
表の構成によって案内され、行列を逆のオーダで計算できる。各ステップにおいて、Ｗの特定の要素がアクセスされる。これが、当該Ｗの要素へのアクセスが必要となる唯一の点である。したがって、好ましいことに、Ｗが、要素ごとのアクセスを提供する方法によって計算される。Ｗのいかなる要素をも保存する必要がない。Ｗはデカルト・ヘシアンから構成されるため、同様の検討がデカルト・ヘシアンにも当てはまる。
【００９６】
状況は、行列Ωについては少し異なっている。要素Ωｂ（ｉ）、ｂ（ｊ）が計算されるや否や、ねじりヤコビアンの（ｉ、ｊ）要素を生成できる。
【００９７】
【数４１】
Ωのこの要素のために使用される保存場所が、好ましくは、当該アルゴリズムにおいて後にそれへの後方参照が発生しないことが明らかになるまでは自動的には解放されないことに、注意すべきである。これは、可能な最大の参照を拘束しているεφの既知の帯域幅から判断できる。例えば、図１の多体系においては、物体７が物体１２についての親である。これは、Ωの５つの列を保存すれば充分であることを意味している。同様に、（７、２）要素は、要素（７、３）、（７、４）、（８、２）、（１２、２）、（８、３）、（８、４）、（１２、３）、および（１２、４）を参照する。これらの要素はすべて、（７、２）からの５列の帯の中に位置している。実際に、再計算を回避するためには、一般に、少数の列を保存すれば充分である。
【００９８】
（ＶＩＩ．コンピュータ・システム）
上記の計算を実行するため、コンピュータ・システムを、少なくとも１つのプロセッサと、プロセッサに前記操作を実行するよう指示するためのコンピュータ・コードを保持するための関連のメモリ・サブシステムとを備えて使用できる。図４が、プロセッサ４０１と、メモリ・サブシステム４０２と、入力／出力装置（キーボード、マウス、ディスプレイ、など）などの周辺機器４０３と、おそらくは計算を補助するためのコプロセッサ４０４と、ネットワークインターフェイス装置４０５とを、すべてバス４００で相互接続して有しているそのようなコンピュータ・システムの基本的なアーキテクチャを示している。メモリ・サブシステムは、アクセス待ちの時間が長くなる順番で、キャッシュメモリ、主メモリ、およびハードディスク・ドライブなどの恒久的保存メモリを、最適に備えている。計算の激しさの程度を考慮し、コンピュータ・システムが、計算を並列に実行するために複数の関連のメモリ・サブシステムとともに複数のプロセッサを備えてもよく、あるいは種々のコンピュータ要素を図４によって示されるような従来からのコンピュータ・アーキテクチャにてバスで接続して有するのではなく、コンピュータ・システムが、ネットワークによって互いに接続される複数のプロセッサおよび複数のメモリ・サブシステムによって形成されてもよいことを、理解すべきである。
【００９９】
以下の例は、本発明を説明するためのものであるが、いかなるかたちでも本発明を限定するものではない。
【０１００】
（実施例１）
（従来の方法および本発明の方法の計算速度の比較）
動的残差ρｕが、以下のとおり原子間力Ｆから計算される。
【０１０１】
【数４２】
ρｕの微分を、混合デカルト‐ねじり導関数行列
【０１０２】
【数４３】
を使用する
【０１０３】
【数４４】
を使用することによって、従来の方法を用いて得ることができる。
【０１０４】
しかしながら、本発明の方法を用いると、今や動的残差の微分を
【０１０５】
【数４５】
として計算できる。
【０１０６】
本明細書に詳細を説明した方法を、以下のように従来のからのアプローチに組み合わせた。便宜上、原子間力を２つの成分に、それらの導関数が２つの別個の項
【０１０７】
【数４６】
で計算されるような方法で分割した。
【０１０８】
【数４７】
したがって、ρｕの微分は、
【０１０９】
【数４８】
として記載され、第１が本発明の実施の形態によるデカルト・ヘシアンであり、第２が既に使用されている混合デカルト‐ねじりヘシアンである。
【０１１０】
（基本アルゴリズム）
以下の表現
【０１１１】
【数４９】
を、本明細書に記載の方法に従って計算し、既存の方法によって計算された
【０１１２】
【数５０】
の残りを加えた。
【０１１３】
これは、
【０１１４】
【数５１】
と書き換えられ、
１．「物体」ヘシアン
【０１１５】
【数５２】
の計算、および
２．以下の並びの演算子を使用するρｕ’の計算
【０１１６】
【数５３】
を可能にした。
【０１１７】
（省メモリのアルゴリズム）
省メモリのアルゴリズムは、Ｏ（Ｎ）中間メモリのみを保存する。動的残差ｑ‐導関数の表現を、以下のとおり書くことができる。
【０１１８】
【数５４】
行列ＷおよびΩを、効率的な要素ごとの方法で計算できる。すなわち、
子および親のインデックス間の差ｍを最小にするため物体を並べ直し、
行列Ωのｍ個の列のためにメモリを予約し、
【０１１９】
【数５５】
ここで、
【０１２０】
【数５６】
上記において、アインシュタインの総和規約が使用される。行列Ｗの要素が、必要に応じ
【０１２１】
【数５７】
から計算される。
【０１２２】
（実施例２）
（ヤコビ行列の計算）
ヤコビ行列が、４つのブロックで構成される。
【０１２３】
【数５８】
左下のブロック以外のすべては、計算が単純である。本明細書に記載の方法によれば、左下のブロックを、質量行列Ｍおよび動的残差ρｕから計算できる。
【０１２４】
【数５９】
最後の項は、本発明において説明したとおり、動的残差の一般化座標に関する導関数を含んでいる。
【０１２５】
（実施例３）
（剛性行列の計算）
系が、ｎｑ個の一般化座標およびｎｕ個の一般化速度を有するとする。制約または４元ジョイント（ｑｕａｔｅｒｎｉｏｎ ｊｏｉｎｔ）が存在しないならば、ｎｑ＝ｎｕである。この場合、一般化座標の独立組に関する最大階数（非特異）の剛性行列を、直接計算することができる。
【０１２６】
【数６０】
これは、ｎｕ×ｎｕの行列である。
【０１２７】
一方、４元および／または制約が存在する場合には、ｎｑ＞ｎｕである。この場合、剛性行列を２つのステップで計算することができ、最初のステップはヤコビアンの計算の部分集合である。
１．（従属）一般化座標の完全な組に関する特異剛性行列が、動的残差ρｕの導関数として計算される。
【０１２８】
【数６１】
これは、ｎｑ×ｎｑの行列である。
２．次いで、非特異剛性行列が、従属一般化座標を除くことによって、この行列から計算される。これは、いくつかある標準的なアプローチのいずれかを使用して行なうことができる。例えば、制約行列の特異値分解を使用することができる。Ｇが制約行列を表わしている場合、ＧＺ＝０を解くことによってＧのヌル空間Ｚを計算する。次いで、最大階数の剛性行列を、
【０１２９】
【数６２】
として計算できる。
剛性行列の用途の１つは、この技術分野においてよく知られているとおり、正規モードの計算である。
【０１３０】
（結果）
以上の実装例を、既存の方法および本発明の方法を使用するリゾチーム（＋／−オメガ）およびトリプシン（＋オメガ）の
【０１３１】
【数６３】
の計算に使用した。結果は、下記の表１にまとめられている。
【０１３２】
【化４】
速度（ＰＡ）は、前記従来の方法を使用した計算速度（単位は秒）であり、速度（ＭＩ）は、本発明の方法を使用した計算速度（単位は秒）である。本発明の方法が、これまで使用されている方法に比べて
【０１３３】
【数６４】
の計算速度を大きく向上させることができることを、理解できるであろう。
【０１３４】
表２は、基本アルゴリズム、および高速ヤコビアン・アルゴリズムの省メモリ（ＬＭ）実装からの結果を示している。
【０１３５】
【化５】
（拡大可能性）
デカルト・ヘシアン行列
【０１３６】
【数６５】
は大きくなる可能性があり、最終的に３２ビットのアドレス空間を超える可能性もある。この問題を解決するため、行列を部分に構成できる。このアルゴリズムは、ヘシアン行列の行について動作させ、Ｗを行ごとに構成し、ρｕ’の対応する列を構成することによって拡大可能である。
【０１３７】
以上、本発明の実施の形態を完全に説明したが、さまざまな変更、代案、および均等物が可能でありかつ利用可能であることは、明らかである。したがって、上記説明は、添付の特許請求の範囲の計量および境界によって定められる本発明の技術的範囲を限定するものとして理解すべきではない。
【図面の簡単な説明】
【０１３８】
【図１】多体系のツリー構造を示している。
【図２】図１に示した多体系についての行列εφである。
【図３Ａ】２つの凍結体間の相互作用を象徴的に示している。
【図３Ｂ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｃ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｄ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｅ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｆ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図４】本発明の方法の実行に有用であるコンピュータ・システムの図である。
【技術分野】
【０００１】
本願は、Ｄａｎ Ｒｏｓｅｎｔｈａｌによる仮特許出願第６０／４７７，２３７号（２００３年６月９日出願）の出願日とＤａｎ Ｒｏｓｅｎｔｈａｌによる仮特許出願第６０／５５２，２２２号（２００４年３月１０日出願）の出願日との優先権の利益を享受し、この２つの出願は、共に本明細書にその全体が援用される。
【０００２】
本発明は、数値的方法の分野に関し、さらに詳しくは、例えば多体機械系（ｍｕｌｔｉｂｏｄｙ ｍｅｃｈａｎｉｃａｌ ｓｙｓｔｅｍ）、とくには分子のモデル化の用途において使用される機械系など、機械系の運動方程式の解法との関連において使用される数値的方法に関する。本明細書にて言及する参考文献はすべて、あらゆる目的において、言及によって本明細書に組み込まれたものとする。
【背景技術】
【０００３】
分子のモデル化は、分子の構造および動的な挙動など、分子のさまざまな側面をシミュレートするために使用できるいくつかの技法を含んでいる。分子のモデル化の典型的な用途には、酵素‐リガンドの結合、分子の拡散、反応経路、相転移、およびたんぱく質の折り畳みの研究などが挙げられる。生物科学ならびに薬品、高分子化合物、および化学品産業における研究者らは、複雑な分子における化学プロセスの性質を理解し、それによって新規な薬剤または物質を設計するために、そのような技法を使用し始めつつある。当然ながら、これらのツールが受け入れられるか否かは、実体の表現における結果の正確さ、モデル化できる分子系のサイズおよび複雑さ、ならびに解を得る際の速度など、いくつかの要因にもとづいている。
【０００４】
分子または他の機械系のシミュレーションは、一般に、微分方程式の組を解くべく数値的方法を使用して実行される。したがって、そのようなシミュレーションの速度は、部分的には、使用される数値的方法の速度によって制約され、数値的方法を使用して解かれる計算およびアルゴリズムの速度によって制約される。機械的な（例えば、分子の）シミュレーションの計算速度は、系内の自由度（「ＤＯＦ（ｄｅｇｒｅｅｓ ｏｆ ｆｒｅｅｄｏｍ）」）の数への依存の程度の点で特徴付けられ、一般にＤＯＦの数は、系内の物体（例えば、原子）の数に比例する。例えば、「オーダ（Ｎ）」のアルゴリズムの計算時間が、系内のＤＯＦの数に線形に比例する一方で、「オーダ（Ｎ平方）」（オーダ（Ｎ２））のアルゴリズムの計算時間は平方の依存性を有し（すなわち、ＤＯＦの数の二乗に比例して増加する）、「オーダ（Ｎ３）」のアルゴリズムの計算時間は立方の依存性を有し、以下同様である。系のサイズ（したがって、ＤＯＦの数）が大きくなるにつれ、より高いオーダ（Ｎ）の依存性を有するアルゴリズムは、計算コストに関して急激に認容できなくなるほどになる可能性がある。したがって、可能な限りオーダ（Ｎ）への依存性が小さいアルゴリズムを開発することが、強く望まれる。
【０００５】
分子のシミュレーションは、一般に、デカルト座標（いずれも絶対基準座標系に関して表現される系内の各原子または物体についてのｘ、ｙ、ｚ座標）または内部座標（典型的には、二面角またはねじり角座標）のいずれかを使用して実行されている。内部座標を使用する利点は、使用されるアルゴリズムの少なくともいくつかにおいて複雑さが増すことと引き換えに、系内のＤＯＦの数を大きく減らす（例えば、デカルト座標の場合の約１／６〜１／８の間まで）ことができる点にある（例えば、Ｒｏｓｅｎｔｈａｌの２００２年８月８日付の国際公開パンフレットＷ００２／０６１６６２号、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」、またはＲｏｓｅｎｔｈａｌの２００３年１月２３日付の米国特許出願公開第ＵＳ２００３／００１８４５５Ａ１号明細書、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」を参照）。
【０００６】
典型的な２種類の分子モデル化技法は、時間にわたって分子系の運動をシミュレートする分子動力学（「ＭＤ」）アプローチ、および系を構成している原子または物体をランダムに変化させることによって分子または分子系について種々の状態を生成し、連続する各状態のエネルギーを評価するモンテカルロ法である（例えば、Ｌｅａｃｈ． Ａ． Ｒ．のＭｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ‐Ｐｒｉｎｃｉｐｌｅｓ ａｎｄ Ａｐｐｌｉｃａｔｉｏｎｓ ２ｎｄ Ｅｄ．、２００１年、Ｄｏｒｓｅｔ Ｐｒｅｓｓ、Ｄｏｒｃｈｅｓｔｅｒ、Ｄｏｒｓｅｔを参照）。
【０００７】
これらの、ならびに他の分子および機械モデル化技法は、構造的なエネルギー関数（位置および／または運動エネルギー）の一次（勾配、またはヤコビアン）および／または二次（へシアン）の導関数行列を利用するように定式化されることが多い。これらの導関数は、デカルト座標においては容易に計算可能であるが、内部座標において第１原理から計算することはきわめて困難である。したがって、内部座標についてエネルギー（Ｅ）または力（Ｆ）の関数のヤコビアンまたはヘシアンを見出したいと望む場合には、一般に、デカルトのヤコビアンまたはヘシアンを対応する内部座標（例えば、ねじり角座標）のヤコビアンまたはヘシアンに変換する必要がある。
【０００８】
デカルトのヤコビアンの計算コストは、オーダ（Ｎ２）と同じく低いが、デカルトのヤコビアンを内部座標のヤコビアンに変換するために一般に適用可能な従来技術の方法は、Ｎについて立方の依存性を有している（すなわち、オーダ（Ｎ３））。したがって、大規模な系について作業を行なうとき、変換のための計算がボトルネックとなる可能性がある。分子のモデル化の文脈においてヤコビアンの計算を高速に（オーダＮ２まで低く）するため、いくつかの方法が説明されている（例えば、Ｇｉｂｒａｔ、１９９７、「Ｆａｓｔ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｆｉｒｓｔ ａｎｄ ｓｅｃｏｎｄ ｄｅｒｉｖａｔｉｖｅｓ ｏｆ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｙ ｏｆ ｏｌｉｇｏｍｅｒｉｃ ｐｒｏｔｅｉｎｓ ｗｉｔｈ ｒｅｓｐｅｃｔ ｔｏ ａ ｓｅｔ ｏｆ ｖａｒｉａｂｌｅｓ ｍｉｘｉｎｇ ｄｉｈｅｄｒａｌ ａｎｇｌｅｓ ａｎｄ ｒｉｇｉｄ ｂｏｄｙ ｖａｒｉａｂｌｅｓ」、Ｊ． Ｃｈｅｍ Ｐｈｙｓ ９４：１２３４‐１２５６；Ｎｏｇｕｔｉ ＆ Ｇｏ、１９８３、Ｊ Ｐｈｙｓ． Ｓｏｃ． Ｊａｐ． ５２：３６８５）が、それらは以下に詳しく述べるとおり、或る限界を抱えている。
【０００９】
力またはエネルギーのヤコビアンを計算するための多くの技法は、連鎖法則を使用しており、例えば力のヤコビアンを内部座標について
【００１０】
【数６】
と表現することができ、ここでデカルト座標がｒとして表わされ、内部座標がｑとして表わされている。
【００１１】
【数７】
は、典型的には、Ｆそのものがポテンシャルエネルギーの勾配として得られるときは常にデカルト・ヘシアンであり、そうでない場合、
【００１２】
【数８】
はデカルト力ヤコビアンと称される。
【００１３】
【数９】
は、ねじり変位勾配と称される。陽解行列乗算法（「ｅｘｐｌｉｃｉｔ ｍａｔｒｉｘ ｍｕｌｔｉｐｌｉｃａｔｉｏｎｍｅｔｈｏｄ」）（例えば、両者ともここでの言及によってあらゆる目的について本明細書に組み込まれたものとするが、Ｒｏｓｅｎｔｈａｌの２００２年８月８日の国際公開パンフレットＷＯ０２／０６１６６２号、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」；またはＲｏｓｅｎｔｈａｌの２００３年１月２３日の米国特許出願公開第ＵＳ２００３／００１８４５５Ａ１号、「Ｍｅｔｈｏｄ ｆｏｒ Ａｎａｌｙｔｉｃａｌ Ｊａｃｏｂｉａｎ Ｃｏｍｐｕｔａｔｉｏｎ ｉｎ Ｍｏｌｅｃｕｌａｒ Ｍｏｄｅｌｉｎｇ」を参照されたい）においては、デカルト・ヤコビアンおよび変位勾配の両者が形成され、陽解行列乗算が実行される。大規模な分子系においてはデカルト・ヤコビアンを保存するために通常でない量のメモリが必要とされるため、一般に、ヤコビアンをブロックに部分的に形成するため、或る特別な規定がなされる。行列の乗算のコストは、原子の数について立方的（オーダ（Ｎ３））である。使用される力の場に応じ、デカルト・ヤコビアンの形成のコストは、平方的（オーダ（Ｎ２））でありうる。変位勾配の形成のコストも平方的であり、したがってこの方法は、最終的には、行列乗算の立方的なコストによって支配される。
【００１４】
数値的摂動にもとづく方法（例えば、Ｌｙｎｅｓｓ ＆ Ｍｏｌｅｒ、１９６７、「Ｎｕｍｅｒｉｃａｌ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ ｏｆ ａｎａｌｙｔｉｃ ｆｕｎｃｔｉｏｎｓ」、ＳＩＡＭ Ｊ Ｎｕｍｅｒ． Ａｎａｌ． ４：２０２‐２１０を参照）も、すでに説明されている。使用されるシミュレーション・エンジンが原子間力を計算できる場合、その原子間力は、一般化座標の陰関数：Ｆ＝Ｆ（ｒ（ｑ））として眺めることができる。したがって、
【００１５】
【数１０】
を、例えば差分を使用して数値的に形成することができる。この仕組みの全体コストは、ヤコビアンの各列が少なくとも１つの力の計算を必要とするため、立方的である。この数値的な仕組みの精度も、丸めおよび切り捨てによる影響ゆえ、苦しいものになりうる。得られた結果の精度を回復させるため、複雑化技法（例えば、Ｍａｒｔｉｎｓら、２０００年、「Ａｎ ａｕｔｏｍａｔｅｄ ｍｅｔｈｏｄ ｆｏｒ ｓｅｎｓｉｔｉｖｉｔｙ ａｎａｌｙｓｉｓ ｕｓｉｎｇ ｃｏｍｐｌｅｘ ｖａｒｉａｂｌｅｓ」、ＡＩＡＡ ２０００−０６８９を参照）を使用することができるが、これは力の場の全体を複雑な算術を使用すべく書き換えることを必要とし、多くの用途において現実的でない要件である。
【００１６】
プログラム強化技法の適用から生じるフォワード・モードの微分（ｆｏｒｗａｒｄ ｍｏｄｅ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ）にもとづく方法も、先行技術において記載されている（例えば、Ｂｉｓｃｈｏｆら、１９９２年、「ＡＤＩＦＯＲ：Ｇｅｎｅｒａｔｉｎｇ Ｄｅｒｉｖａｔｉｖｅ Ｃｏｄｅｓ ｆｒｏｍ Ｆｏｒｔｒａｎ Ｐｒｏｇｒａｍｓ」、Ｓｃｉｅｎｔｉｆｉｃ Ｐｒｏｇｒａｍｍｉｎｇ １：１１−２９）。これらの方法の適用は、Ａｄｉｆｏｒ（Ａｒｇｏｎｎｅ Ｎａｔｉｏｎａｌ ＬａｂｏｒａｔｏｒｙおよびＲｉｃｅ Ｕｎｉｖｅｒｓｉｔｙ）などのコード生成プログラム、または手動のいずれかによって、行なうことができる。そのような方法の適用は、渡された（ｐａｓｓｅｄ‐ｉｎ）データ・ベクトルｚについての行列ベクトル積
【００１７】
【数１１】
を計算できるプログラムをもたらす。そのようなプログラムの良好な使用は、変位勾配を形成し、フォワード・モードの微分ルーチンに渡すことにある。これは、依然として立方的な全体コストをもたらすが、数値的差分法よりも約３倍高速であって、完全な精度が手に入る。
【００１８】
Ｎｏｇｕｔｉ ＆ Ｇｏ，１９８３年、（Ｊ． Ｐｈｙｓ． Ｓｏｃ． Ｊａｐ． ５２：３６８５）の方法を脚色してなるＧｉｂｒａｔの方法（Ｇｉｂｒａｔ、１９９７年、「Ｆａｓｔ Ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｆｉｒｓｔ ａｎｄ ｓｅｃｏｎｄ ｄｅｒｉｖａｔｉｖｅｓ ｏｆ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｙ ｏｆ ｏｌｉｇｏｍｅｒｉｃ ｐｒｏｔｅｉｎｓ ｗｉｔｈ ｒｅｓｐｅｃｔ ｔｏ ａ ｓｅｔ ｏｆ ｖａｒｉａｂｌｅｓ ｍｉｘｉｎｇ ｄｉｈｅｄｒａｌ ａｎｇｌｅｓ ａｎｄ ｒｉｇｉｄ ｂｏｄｙ ｖａｒｉａｂｌｅｓ」、Ｊ． Ｃｈｅｍ Ｐｈｙｓ ９４：１２３４−１２５６）においては、ヤコビアンおよびねじりヘシアンが、平方的なコストで生成されるが、ポテンシャルの組として表現できるエネルギー関数に限定されている。すなわち、１‐４ねじりポテンシャル（ねじりポテンシャルは、４原子項の加法和として表現されるため）、１‐３結合角度ペア、およびポテンシャルの組として専ら表現することが不可能である他のエネルギー関数など、単純なエネルギー関数には適用不可能である。さらに、真の多体ポテンシャルとして表現され、陰解法モデルを利用した分子シミュレーションに使用されることが多い一般化ボーン・ソルベント・アクセシビリティ（Ｇｅｎｅｒａｌｉｚｅｄ Ｂｏｒｎ Ｓｏｌｖｅｎｔ Ａｃｃｅｓｓｉｂｉｌｉｔｙ：「ＧＢＳＡ」）解法モデルなど、分子のモデル化に使用され、ポテンシャルの組として専ら表現することが不可能である他のいくつかのツールへの適用も不可能である。
【００１９】
本発明は、従来技術における前記およびその他の制限を克服するものであり、内部座標（例えば、ねじり角）ヤコビアンを平方的な（オーダ（Ｎ２）の）コストで計算するための高速かつ汎用的に適用可能な方法であって、幾多の用途の中でも、とりわけ分子のモデル化方法を劇的に高速化するために使用できる方法を提供する。
【発明の開示】
【課題を解決するための手段】
【００２０】
一実施の形態において、本発明は、ペア‐ポテンシャルエネルギー項の制約なしにデカルト・ヘシアンをねじりヤコビアンに変換するために有用であり、かつコンピュータのメモリを高効率で使用しつつそのような計算を実行するために有用である方法およびアルゴリズムを提供する。そのような方法およびアルゴリズムは、例えば分子系などの任意の多体系をモデル化するために使用されている内部座標の数について平方的な計算時間で変換を行なうことができる。関連の実施の形態において、本発明は、ペア‐ポテンシャルエネルギー項の制約なしに剛性行列を計算するために有用であり、かつコンピュータのメモリを高効率で使用しつつそのような計算を実行するために有用である方法およびアルゴリズムを提供する。
【発明を実施するための最良の形態】
【００２１】
（Ｉ．定義）
そのようではないと示されない限り、「約（ａｂｏｕｔ）」とは、記載の値±３０％を意味する。
【００２２】
そのようではないと示されない限り、「物体（ｂｏｄｙ）」とは、分子の表現の構成要素の文脈において、単一の質量または分子のモデル化の目的のための幾何構造として取り扱われる当該表現の単位として定義される。したがって、物体は、分子の個々の原子、原子の集まり、または質量からなる他の抽象系の表現であってよい。「剛体（ｒｉｇｉｄ ｂｏｄｙ）」は、剛であるとしてモデル化される物体である（すなわち、力が加わっても変形しない）。
【００２３】
そのようでないと示されない限り、「高速演算子の実装（ｆａｓｔ ｏｐｅｒａｔｏｒ ｉｍｐｌｅｍｅｎｔａｔｉｏｎ）」とは、積行列の要素が計算される行列の乗算の評価の方法であって、積行列の全体を計算および保存する必要なく「オンザフライ（ｏｎ‐ｔｈｅ‐ｆｌｙ）」で使用または解放される方法を意味する。したがって、高速オーダの実装は、標準的な行列の乗算（オーダ（Ｎ３））に比べ、実行速度および保存用メモリの両者の要件に関して大幅に効率的でありうる（オーダ（Ｎ２）と高い）。
【００２４】
そのようでないと示されない限り、「ジョイント（ｊｏｉｎｔ）」は、２つの物体間の接続である。一般的なジョイントの種類としては、ピン・ジョイント、スライダ・ジョイント、およびボール・ジョイントが挙げられるが、さらにこれらに限られるわけではないが、フリー・ジョイント、Ｕジョイント、円柱ジョイント、ベアリング・ジョイント、および前記いずれかの組み合わせなど、他にも多数の種類のジョイントが可能である。例えば、フリー・ジョイントは、直交する３つのスライダ・ジョイントをボール・ジョイントに組み合わせて構成され、完全な６自由度を有している。ジョイントの種類およびそれらの仕様についてのさらに詳細な検討を、例えば上記Ｒｏｓｅｎｔｈａｌに見つけることができる。
【００２５】
そのようでないと示されない限り、特定のジョイントまたはピボットにおける「凍結体（ｆｒｏｚｅｎ ｂｏｄｙ」」は、当該ピボットの外にあって剛体として一体に運動するすべての物体で構成される。
【００２６】
そのようでないと示されない限り、「計算実現可能なオーダ（ｃｏｍｐｕｔａｔｉｏｎａｌｌｙ‐ｆｅａｓｉｂｌｅ ｏｒｄｅｒ）」とは、タスクの特定の並びを最終的な結果を変化させることなく実行できる任意のオーダを指す。このコンセプトは、いくつかの方法においては特定のステップのオーダが、それらステップが実行され、かつ結果がそれらステップが当初提示されたオーダで実行された場合と実質的に同じである限りにおいて、重要ではないことから惹起されている。
【００２７】
そのようでないと示されない限り、用語「シリコン内（ｉｎ ｓｉｌｉｃｏ）」は、コンピュータを使用して実行される任意の方法またはプロセスを指す。
【００２８】
そのようでないと示されない限り、「分子（ｍｏｌｅｃｕｌｅ）」は、化学的結合によって接続された２つ以上の原子で形成される任意の巨視的構造である。分子の例としては、これらに限られるわけではないが、たんぱく質（例えば、抗体、受容体、など）、ペプチド、脂質、核酸（例えば、天然または合成のＤＮＡ、ＲＮＡ、ｇＤＮＡ、ｃＤＮＡ、ｍＲＮＡ、ｔＲＮＡ、など）、レクチン、糖（例えば、レクチン‐糖の複合体を形成している）、糖たんぱく、低分子、有機化合物、塩、金属、などの単原子または多原子の構造が挙げられる。
【００２９】
そのようでないと示されない限り、「表現（ｒｅｐｒｅｓｅｎｔａｔｉｏｎ）」とは、分子または分子構造に適用されるとき、例えばコンピュータ・シミュレーションにおいて機械によって使用するための分子または分子構造の抽象的説明を指す。例えば、分子の１つの表現は、分子を構成している原子または物体、あるいはそれらの抽象的代理の位置を集合的に定める座標の組である。分子の表現の例には、これらに限られるわけではないが、Ｘ線座標、「ｐｄｂ」（たんぱく質データバンク）ファイル、などが挙げられる。
【００３０】
そのようでないと示されない限り、「大きく（ｓｉｆｇｎｉｆｉｃａｎｔ）」という用語は、「大きく異なる」、「大きく抑制する」、または「大きく促す」に関して使用されるとき、比較される２つのグループ間の定量化可能なパラメータにおける差異であって、標準的な統計的検定を使用して統計的に有意である差異を指す。
【００３１】
そのようでないと示されない限り、「溶媒（ｓｏｌｖｅｎｔ）」は、溶質分子を含みうる任意の媒体を指す。溶媒の例としては、これらに限られるわけではないが、水および他の水性溶媒、ならびに有機溶媒（例えば、ＤＭＳＯ、脂質、アルコール、など）が挙げられる。溶媒は、一様であっても、非一様であってもよく、固体状、液体状、または気体状であってよい。
【００３２】
そのようでないと示されない限り、「ターゲット分子（ｔａｒｇｅｔ ｍｏｌｅｃｕｌｅ）」とは、分子シミュレーションの対象である主たる分子である。シミュレーションは、例えば互いに相互作用する２つ以上のたんぱく質のシミュレーションなど、２つ以上のターゲット分子を有することができる。
【００３３】
そのようでないと示されない限り、「汎用的に適用可能な方法（ｕｎｉｖｅｒｓａｌｌｙ‐ａｐｐｌｉｃａｂｌｅ ｍｅｔｈｏｄ）」とは、ポテンシャルの組における使用に限られず、分子のシミュレーションにおける使用に適した任意の方法である。
【００３４】
（ＩＩ．多体系（ＭＵＬＴＩ‐ＢＯＤＹ ＳＹＳＴＥＭＳ））
一実施の形態において、本発明の方法は、機械系の多体系（ＭＢＳ）の実装に関して使用される（例えば、ＭＢＳ分子シミュレーション；前記Ｒｏｓｅｎｔｈａｌを参照）。ＭＢＳの基本的な抽象化は、ジョイントで接続された剛体の集まりの抽象化である。物体のうちのベースと呼ばれる１つが、その運動が直接地面に参照付けられる特別な地位を有する。系のグラフは、ループを含みうる。ループ切断アルゴリズムが、循環するグラフをツリーおよび切断の組に還元する。切断は、ジョイントまたは物体において生じうる。これにより、すべてのジョイントがツリー系に置かれる。ループ切断アルゴリズムによって切断された物体は、溶接ジョイントを課すことによって回復される。このジョイントそのものは、６つの距離の制約の組へと分解される。溶接ジョイントが、元の物体の２つの部分を再組み立てする。
【００３５】
ツリーの重要な特性は、或る任意の物体から他の任意の物体への経路が一意である点にある。ツリー内の物体の数はｎであり、通常は、ベース物体にラベル１が割り当てられる。これは、ベース物体から任意の葉物体（ｌｅａｆ ｂｏｄｙ）へと到るあらゆる経路において、物体のラベルが減少することがないということを意味する。葉物体は、ただ１つの物体に接続された物体である。規則的なラベル付けを、葉物体（少なくとも１つ存在しなければならない）のうちの１つにラベルｎを割り当てることによって達成できる。この物体がグラフから除かれるならば、ｎ−１個の物体を有するツリーが残される。このツリーの葉物体のうちの１つにラベルｎ−１を割り当て、このプロセスを、すべての物体がラベル付けされるまで繰り返す。一般的な機械系のＭＢＳの実装の詳細は、すでに説明されている（例えば、Ｒｅｃｕｒｓｉｖｅ Ｍａｓｓ Ｍａｔｒｉｘ Ｆａｃｔｏｒｉｚａｔｉｏｎ ａｎｄ Ｉｎｖｅｒｓｉｏｎ、ＪＰＬ Ｐｕｂｌｉｃａｔｉｏｎ ８８‐１１、１９８８におけるＲｏｄｒｉｇｕｅｚおよびＫｒｅｕｔｚを参照）。分子系のＭＢＳの実装の詳細、ならびに関連の表記法も、やはり説明されている（例えば、前記Ｒｏｓｅｎｔｈａｌを参照）。
【００３６】
（ＩＩＩ．残差の計算）
典型的な実施の形態において、本発明の方法は、
【００３７】
【数１２】
の形式の運動方程式を有する多体系として表現された系（例えば、分子系）のシミュレーションに適用される。この方程式の右手側は残差ρであり、慣性力および力の場からのアクティブな力の寄与を含んでいる。残差のヤコビアンは、
【００３８】
【数１３】
である。本発明の方法は、アクティブな力の成分を表現している力の場からの残差への寄与を計算するために適用可能である。多体系は、原子間力の計算をトリガして、残差ベクトルの要素へと集めることができるデータ処理方法の集まりを具現化している。これらのデータ処理ステップは、特定の行列操作の点で高いレベルの表現を有することができるが、実際には、これらの操作の実際の実施は、典型的には、実行時間が問題の大きさと線形に比例するアルゴリズムのかたちである。それらは、いわゆるＯ（ｎ）法であって、先行技術に教示されている（例えば、前記Ｒｏｓｅｎｔｈａｌを参照）。しかしながら、本発明の一態様によれば、これらＯ（ｎ）法を、ヤコビアンの形成ステップに適用することができる。
【００３９】
例えば、残差は、演算子方程式
【００４０】
【数１４】
から形成される。
【００４１】
Ｆは、典型的には、３ｎａ×１のベクトルである。各エントリは、好ましくは、所与の原子に作用する全原子間力の大域的なデカルト成分で構成される。行列演算子Ｈ、Φ、およびＰは、Ｒｅｃｕｒｓｉｖｅ Ｍａｓｓ Ｍａｔｒｉｘ Ｆａｃｔｏｒｉｚａｔｉｏｎ ａｎｄ Ｉｎｖｅｒｓｉｏｎ、ＪＰＬ Ｐｕｂｌｉｃａｔｉｏｎ ８８‐１１、１９８８においてＲｏｄｒｉｇｕｅｚおよびＫｒｅｕｔｚによって記載されている。力分布行列Ｐは、多体系のピボットにおいて原子間力から作用している多体空間荷重を生成する。したがって、
【００４２】
【数１５】
と書くことができ、ここでＴは、典型的には６ｎｂ×１のベクトルである。特定の原子に作用する所与の力が、当該原子が位置する物体のピボット点に作用する力およびトルクへと写像される。このように、行列Ｐの所与の行が、特定の物体に作用する空間荷重に対応する。所与の行は、当該特定の物体上に位置する原子に対応するエントリについてのみ、ゼロでないブロックを有する。分子系の多体定式化の多くにおいて、各物体は、典型的には、約３つの原子を含み、あるいは表現している。したがって、行列Ｐは大きくかつスパースであり、明示的に形成される必要はない（例えば、前記Ｒｏｓｅｎｔｈａｌを参照）。むしろ、行列ベクトル積ＰＦを計算するアルゴリズムを使用することができる。そのような計算の実効時間は、原子の数に線形である（前記Ｒｏｓｅｎｔｈａｌ）。
【００４３】
多体の力学においては、さらに演算子の転置が有用な量であることが一般的である。この場合、演算子ＰＴが、空間速度（直線および角）を、各物体ピボットから、当該物体に位置する原子位置へと伝播させる。この行列の典型的な用途は、分子系内の原子の微分空間変位の計算の目的のためにある。
【００４４】
行列Φは、多体遷移演算子として知られている。この行列は、データ・ベクトルに作用し、ツリーのベース物体の下流の葉物体からのデータ・ベクトルの要素のシフトおよび累積の効果を有している。例えば、積ΦＴが、反応荷重の６ｎｂ×１のベクトルＲを生成する。所与の物体について、反応荷重の要素は、そのジョイントを中心とする当該物体の運動によって運動力学的に影響される物体に作用するそれらの力の物体内ピボットについての結果としての力およびトルクを表わしている。ΦＴの計算は、物体の数について線形である。
【００４５】
Φを理解するために、最初に行列εφを導入することが有用である。この（ブロック）行列は、スパースであってスーパー・ダイアゴナルである。ｉ番目の行において、物体ｉの子に対応する列のみが満たされている。図１に示した物体の図について、行列εφが図２に示されている。εφのゼロでない各要素は、６×６の行列である。要素ｉφｋが、空間荷重を子のピボットからその親のピボットへとシフトする。転置された要素は、空間速度を親ピボットから子ピボットへとシフトする。
【００４６】
今や、Φとεφとの間の関係を、
【００４７】
【数１６】
と記述することができる。
【００４８】
この方程式の意味は、図１に示した系に関して、各ピボットに加えられる空間荷重を考えて各ピボットにおける空間反応荷重の計算を説明する以下の説明を参照することによって、よりよく理解できるであろう。反応荷重は、以下のＯ（ｎ）スイープ（コードのシークエンス）によって計算することができる。
【００４９】
【数１７】
これらの方程式は、各ピボットにおける反応を累積させる。スイープは、葉物体から出発して、ベース物体へとツリーをさかのぼる。各物体が、自身の反応を自身の親へとシフトさせ、そこで兄弟からの寄与に加算される。これらの方程式が、先に計算された量への後方参照を行なうことに注意すべきである。例えば、Ｒ（６）はＲ（７）を参照し、Ｒ（７）がＲ（８）およびＲ（１２）を参照し、以下同様である。しかしながら、このアルゴリズムが提示のとおりに構成される場合、不満足な参照に遭遇することはない。演算子方程式Ｒ＝ΦＴと上記手順との間の対応は、以下のとおり説明できる。すなわち、この手順が、Ｒ＝εφＲ＋Ｔを示しており、ここから（Ｉ−εφ）Ｒ＝Ｔ、またはＲ＝（Ｉ−εφ）‐１Ｔ＝ΦＴである。他方の方向（方程式からコード・シークエンス）に進むためには、これらのステップを単純に反転できる。
【００５０】
【数１８】
Φの要素を露出したいと望む場合には、例えばＲの計算において繰り返しの置換を実行することによって後方参照を除き、次いでΦの行列要素を特定することができる。あるいは、アルゴリズム的な代替案を使用することも可能である。（Ｉ−εφ）Φ＝Ｉであるため、Φ＝εφΦ＋Ｉが得られ、この方程式は、Ｒ＝εφＲ＋Ｔと同じ構造を有している。Φの列に同じアルゴリズムを、最後の列から出発して最初へと後ろ向きに作業して適用できる。しかしながら、明らかであろうが、典型的には、演算子の作用のみが必要とされるため、Φの要素が必要になる状況が存在しない。したがって、Φ（密な上三角行列）による乗算を示すすべてのアルゴリズムが、スパースなスーパー・ダイアゴナル行列εφに関して進行できる。残差の計算を終えるため、行列ベクトル積ＨＲが、ベクトルＲの要素をジョイントの自由に投影する。行列Ｈは、ブロック対角である。各ブロックは、特定のジョイントについてのジョイント・マップで満たされている。ピン・ジョイントについてのマップは、１×６のベクトル［λ ０］であり、単位ベクトルλがピンの幾何学的性質を指定している。行列ベクトル積ＨＲは、非再帰的なドット積のシークエンスを代表するため、線形なコストで容易に計算される。演算子Ｈは、以下の形状を有しており、ここでｎは物体の数であり、ｎｕは一般化速度の数である。
【００５１】
【数１９】
以下の形状を有する対応する演算子
【００５２】
【数２０】
が存在し、ここでｎｑは一般化座標の数である。
【００５３】
【数２１】
前記に照らし、残差の計算を力の場によって計算された原子間力の要素へと適用されるＯ（ｎ）法のシークエンスとして見ることができることを、理解できるであろう。これらは、例えば、真空の項、溶媒の項（極性または非極性）、および速度関連の項を含む。これらの項目のいくつかは、平方的な計算時間を必要とする。例えば、多くの実装例において、静電およびＧＢＳＡ溶媒モデルは、一般に、計算時間の大部分の力の計算に消費する。したがって、残差コストは、原子の数について平方的である。これは、残差のヤコビアンが、最良でも平方的なコストであり、おそらくはより悪いことを意味する。
【００５４】
ρ＝ＨΦＰＦを考えると、
【００５５】
【数２２】
を
【００５６】
【数２３】
と表現することができる。
【００５７】
原子間力ベクトルＦは、第１および第２項の両者に現れる。しかしながら、第２項は、力ベクトルの解析的形態についての知識を必要とせず、力ベクトルの数値のみを必要とする。したがって、この項は、この項の計算の際に新しい力の評価は不要であるため、元のコストがＯ（ｎ）である項の微分を呈する。したがって、この項がヤコビアンにおいて平方的なコストのみを引き起こしうるということを実証できる。これは、例えばＢｉｃｈｏｆら、１９９４、（「Ａｕｔｏｍａｔｉｃ ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ：ｏｂｔａｉｎｉｎｇ ｆａｓｔ ａｎｄ ｒｅｌｉａｂｌｅ ｄｅｒｉｖａｔｉｖｅｓ−ｆａｓｔ」Ｐｒｏｃ． ＳＩＡＭ Ｓｙｍｐｏｓｉｕｍ ｏｎ Ｃｏｎｔｒｏｌ Ｐｒｏｂｌｅｍｓ ｉｎ Ｉｎｄｕｓｔｒｙ Ａｒｇｏｎｎｅ、Ｐｒｅｐｒｉｎｔ ＭＣＳ−Ｐ４８４−１１９４）に記載のとおり、残差をフォワード・モードの微分を使用し、力ベクトルを既知の定数として取り扱って微分でき、これによりコストが（コール当たりの）元のプログラムのコストの定数倍によって拘束されるプログラムがもたらされることに注目することによって理解できる。定数は、プログラム中に現れる固有の関数の性質によって決まる。多体系は、算術的な操作に加え、三角法のただ１つの変関数のみを含んでいる。すなわち、一般的に、定数は、典型的には、元のプログラムのコストの１〜２倍のオーダ、またはそれ以下であり、第２項は、ヤコビアン全体について最大でも平方的なコストで容易に計算できる。
【００５８】
【数２４】
を含んでいるヤコビアンの第１項の計算は、以下で説明される。
【００５９】
【数２５】
本発明の一態様によれば、変位勾配を、Ｏ（ｎ）多体演算子Ｐ、Φ、Ｈに関して表現できる。
【００６０】
【数２６】
とくに多体系の実施に適用可能な一実施の形態においては、これらの演算子が、平方的なコストで変位勾配を形成するために使用される。これらの演算子を使用しない場合、コストは立方的であろう。このアプローチは、変位勾配が明示的には必要とされず、デカルト・ヤコビアンによる行列積のみが必要とされるという事実を利用している。積は、
【００６１】
【数２７】
であるため、Ｏ（ｎ）演算子に関して形成できる。
【００６２】
これは、残差の計算に関してすでに説明したものと同じ演算子である。したがって、この方法は、デカルト・ヤコビアンが立方的なコストで計算可能ならば、内部座標（例えば、ねじり）ヤコビアンを、列について線形なコストで、すなわち平方的な全体コストでもたらす。
【００６３】
上記関係を組み合わせて、
【００６４】
【数２８】
をもたらすことができ、ここで行列
【００６５】
【数２９】
である。行列Ｗが多体系において形成される場合、この行列の一般要素Ｗｉｊは、空間における多体系の各物体の配置を反映する。この要素が、物体ｊのピボットの空間的変位ゆえに物体ｉに作用する空間荷重の微分の計算に使用される。Ｗの要素の計算は、デカルト・ヘシアンへの接続がなされる所であり、方程式
【００６６】
【数３０】
を有している。
【００６７】
このアプローチのいくつかの特徴は、注目を有している。すなわち、（ｉ）この方程式がヘシアン要素を使用する。それは、物体ｊの原子リストに対する物体ｉの原子リストによって案内される。また、（ｉｉ）各物体が、典型的にはおおまかに３つの原子で構成されるため、これは単純な計算である。（ｉｉｉ）各原子が、ただ１つの物体に属するため、デカルト・ヘシアンの各要素は、ただ１つの要素Ｗｉｊにのみ寄与する。さらに、従来技術の方法（例えば、Ｇｉｂｒａｔ Ｊ． Ｆ．（１９９７）「Ｆａｓｔ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｔｈｅ ｆｉｒｓｔ ａｎｄ ｓｅｃｏｎｄ ｄｅｒｉｖａｔｉｖｅｓ ｏｆ ｔｈｅ ｃｏｎｆｏｒｍａｔｉｏｎａｌ ｅｎｅｒｇｙ ｏｆ ｏｌｉｇｏｍｅｒｉｃ ｐｒｏｔｅｉｎｓ ｗｉｔｈ ｒｅｓｐｅｃｔ ｔｏ ａ ｓｅｔ ｏｆ ｖａｒｉａｂｌｅｓ ｍｉｘｉｎｇ ｄｉｈｅｄｒａｌ ａｎｇｌｅｓ ａｎｄ ｒｉｇｉｄ‐ｂｏｄｙ ｖａｒｉａｂｌｅｓ」、Ｊ． Ｃｈｅｍ Ｐｈｙｓ ９４：１２３４‐１２５６）と異なり、上述の方法は、ヘシアン要素を直接使用する。したがって、ポテンシャルの組について制約はない。Ｇｉｂｒａｔ（前記）による論文中の対応する量Ｌｉｊは、物体ｊの原子に対する物体ｉの原子から直接的に構成されるが、一般には、Ｌｉｊ≠Ｗｉｊである。本明細書に記載するようなヘシアン要素の使用は、より単純な再帰シークエンスにつながる。
【００６８】
高速演算子の実装においては、全体の計算が、２つの構成されたダウンループ（ｄｏｗｎ‐ｌｏｏｐ）を呈する。項
【００６９】
【数３１】
が、以下の再帰的なアルゴリズムによって計算される。まず、Ｙ＝ΦＷ、Ｚ＝ＨＹとする。これは、
【００７０】
【数３２】
をもたらす。
【００７１】
すでに述べたとおり、行列εφはサブダイアゴナルであり、要素ｉφｋ（ｋ）で占められている。したがって、「積層された」方程式を見つけることができる。
【００７２】
【数３３】
これらの方程式において、Ｗ１は最初の（論理的な）行を表わしており、ＷｎはＷの最後の（論理的な）行を表わしている。これは、６×６ｎｂのブロックである。Ｗの各行は、先の行に関して計算される。Ｈがブロック・ダイアゴナルな構造を有しているため、Ｚ＝ＨＹの計算は瑣末である。
【００７３】
今や、第２の再帰を、以下のとおり表現できる。
【００７４】
【数３４】
演算子ＨΦが、Ｚの行に作用していると見ることができる。これは、Ｚを行のオーダで生成した第１の再帰と好都合に接合する。第２の再帰は、目的の行列が対称であるため、対角線までの所与の列の要素の計算のみを必要とする。
【００７５】
以上に照らし、原子間力ヤコビアンを収縮させてＷを形成できることを、理解できるであろう。全体アルゴリズムは、力のヤコビアンへの行（または、列アクセス）のみを必要とし、同時に行列全体にアクセスする必要はない。さらに、所望であれば第１および第２の再帰を散布できることも、理解できるであろう。このようにして、ＹまたはＺの要素を、それらが計算されるや否や解放することができる。したがって、１つのステップでＹまたはＺのすべてを形成する必要はない。
【００７６】
（Ｖ．メモリに関して効率的な計算）
前記の方法は、Ｗの対称性をさらに充分に利用することによって、計算を実行するために必要となるメモリの量を最小化するやり方で、実装することができる。メモリの必要性を少なくするための１つのアプローチは、デカルト・ヘシアンを大きなブロックで生成し、次いで各ブロックを処理する「スライス化」技法を使用することにある。このような方法は、メモリの使用を少なくするが実行時間が苦しくなりうるという意味で、実行可能であろう。以下で説明する方法は、実行速度を犠牲にすることなくメモリの必要性を大きく減らすことができる処理の仕組みをもたらすように設計されている。
【００７７】
一実施の形態において、この方法は、ねじりヤコビアンを
【００７８】
【数３５】
と表現することにもとづいている。
【００７９】
新しい行列Ωは、加えられた荷重Ｔに対して反応Ｒが有している関係と同じ関係を、Ｗに対して有している。行列Ｗは、多体系の物体についての空間へシアンである、空間に浮遊するが、ジョイントによって接続されていない物体について関連性があるのは、ヘシアンである。行列Ωの要素が、凍結体の「組」を結び付ける。各ピボットにおいて、凍結体は、ピボットの外のすべての物体であって剛体として一体に運動する物体で構成されている。
【００８０】
各凍結体のピボットにおける微分空間変位が、他のすべての凍結体のピボットにおける微分空間荷重を生じさせる。変位させられた物体を「ソース」物体として考え、荷重を受けている物体を「受け手」と考えると、各ソース物体が各受け手物体と相互作用することを理解できるであろう。しかしながら、ソース物体における空間的移動は、凍結体のピボットを起源とし、ΦＴの作用を通じてソース系内の他の物体に伝播する。各ソース物体が、結合行列Ｗを通じて受け手系の各物体上に微分空間荷重を生む。この空間荷重は、ＷΦＴによって与えられる。受け手系の空間荷重は、（ＷΦＴ）に作用する演算子Φの作用によって受け手ピボットへと総計される。
【００８１】
高速演算子の実装において、Ωは以下のとおり効率的に計算できる。
【００８２】
【数３６】
最後の方程式が、Ωの計算のための基礎をもたらしている。εφがスーパー・ダイアゴナルである（主たる対角線の上方のみがゼロでない）ため、項εφΩは、Ωの低い番号の行をより高い番号の行に結合させているが、その逆は真ではない。すなわち、所与の列において、要素は、Ωがεφで前もって乗算されるとき、列を上方へ移動する。同様に、項ΩεφＴは、低い番号の列をより高い番号の列に結合させるだけである。要素は、所与の行を横切って左方に移動する。第３の項は、初期の行および初期の列の両者を、より高い行およびより高い列に結合させているが、結合は、他の方向においては生じない。要素は、上、左、または上および左へと移動する。正味の結果は、ツリー系が規則的にラベル付けされている限り、Ωが常に、逆の列のオーダで計算されるというものである。
【００８３】
望ましいラベル付けは、εφの帯域幅を最小にするラベル付けである。εφの帯域幅は、行列の「リーチ」を制限し、Ωの計算において後方への参照がどの程度まで「後ろ」に延びることができるのかを決定する。
【００８４】
一例として、上記アルゴリズムを、以下の方程式を計算することによって図１に示した系の部分集合に適用できる。
【００８５】
【数３７】
最初の２つの方程式は、結合行列Ｗを通じての物体の直接的寄与のみを取り上げている。（１２、１０）の要素が、凍結体１０の空間変位による凍結体１２への空間荷重を計算する。凍結体１０は、実際の物体１０および１１で構成されている。１０のピボットにおける変位が、ｉφｋ＊（１１）を通じて１１に伝播する。この変位が、結合要素Ω（１２、１１）を通じて１２への空間荷重を生む。（１２、７）の要素は、さらにわずかに複雑である。凍結体７は、物体７ならびに凍結体８および１２で構成される。これらは、物体７の子に相当する凍結体である。物体７のピボットにおける変位が、その子へと伝播し、Ωのすでに計算された要素に結合し、物体１２のピボットへの荷重を生成する。
【００８６】
再帰なしで、子の凍結体を構成している物体の集まり全体にわたって計算が必要であろう。これは、計算が進むにつれて成長傾向のリストとなるであろう。再帰ありでは、子の凍結体のみを訪れるだけでよい。さらに、各物体が実際には少数の子しか有していない（通常は２つまたは３つであるが、より多い場合もある）ため、要素当たりのコストは比較的小さい。
【００８７】
以下で説明され、図３Ａ〜３Ｆに示されるアプローチは、２つの凍結体間の相互作用の計算を示している。説明されている操作は、実際には、前記方程式の計算の手順において「自動的」に実行されるが、これらの方法のよりよい理解および実施を可能にするために示されている。検討対象の物体およびそれらの子の追跡を容易にするため、下記の表１のような表が作成される。
【００８８】
【化１】
図３Ａに目を向けると、この図は、凍結体３１２と凍結体３１４との間の相互作用３１０を示している。この例ａ≧ｂにおいて、凍結体３１２（「受け手物体」）は、物体ａ、その子ｋ、およびａの１つの孫を含んでおり、凍結体３１４（「ソース物体」）は、物体ｂ、その子ｒ、およびｂの２つの孫を含んでいる。相互作用３１０は、任意の計算実現可能なオーダにて実行される図３Ｂ〜３Ｄに示すステップによる分解によって決定できる。例えば、最初に、２つの個々の物体ａおよびｂの間の直接項３１６が、図３Ｂに示すとおり決定される。次いで、図３Ｃに示すとおり、ｂ３２２の凍結した子と相互作用する凍結体ａ３２０によって生成された項３１８（ΩεφＴから得られる）が、加算される。必要となる最後の相互作用（３２４）は、図３Ｄに示すとおり、物体ｂ３２６とａの凍結した子であるｋ（３２８）との間の相互作用である。この項は、好ましくは、物体ｂが単独では凍結体ではないため、間接的に計算される。一実施の形態においては、ａの凍結した子に対する物体ｂの相互作用は、減算によって計算される。すなわち、最初に、図３Ｅに示されるとおり、凍結体ｂ３１４とａの凍結した子３２８との間の相互作用３３０が計算され、εφΩから得られる項を生じさせる。次いで、図３Ｆに示すとおり、ａの凍結した子とｂの凍結した子との間の相互作用３２２（項εφΩεφＴ）が、減算される。これによって、図３Ｄに示すとおりａの凍結した子に対するｂの相互作用がもたらされる。
【００８９】
一例として、前記操作を、図１に示した多体系の一部分の解析へと適用することができる。第１に、両方の物体が子を有している任意の要素が一般的な場合である。そのような要素の例は、Ω（８、６）である。
【００９０】
【化２】
【００９１】
【数３８】
対角要素は、非対角の要素と同じ方法によって処理される。例は、Ω（５、５）について与えられる。
【００９２】
【化３】
【００９３】
【数３９】
当然ながら、行列の対称性ゆえ、Ω（５、６）＝Ω（６、５）Ｔであることを理解すべきである。要約すると、各要素が、式
【００９４】
【数４０】
を使用して計算される。
【００９５】
表の構成によって案内され、行列を逆のオーダで計算できる。各ステップにおいて、Ｗの特定の要素がアクセスされる。これが、当該Ｗの要素へのアクセスが必要となる唯一の点である。したがって、好ましいことに、Ｗが、要素ごとのアクセスを提供する方法によって計算される。Ｗのいかなる要素をも保存する必要がない。Ｗはデカルト・ヘシアンから構成されるため、同様の検討がデカルト・ヘシアンにも当てはまる。
【００９６】
状況は、行列Ωについては少し異なっている。要素Ωｂ（ｉ）、ｂ（ｊ）が計算されるや否や、ねじりヤコビアンの（ｉ、ｊ）要素を生成できる。
【００９７】
【数４１】
Ωのこの要素のために使用される保存場所が、好ましくは、当該アルゴリズムにおいて後にそれへの後方参照が発生しないことが明らかになるまでは自動的には解放されないことに、注意すべきである。これは、可能な最大の参照を拘束しているεφの既知の帯域幅から判断できる。例えば、図１の多体系においては、物体７が物体１２についての親である。これは、Ωの５つの列を保存すれば充分であることを意味している。同様に、（７、２）要素は、要素（７、３）、（７、４）、（８、２）、（１２、２）、（８、３）、（８、４）、（１２、３）、および（１２、４）を参照する。これらの要素はすべて、（７、２）からの５列の帯の中に位置している。実際に、再計算を回避するためには、一般に、少数の列を保存すれば充分である。
【００９８】
（ＶＩＩ．コンピュータ・システム）
上記の計算を実行するため、コンピュータ・システムを、少なくとも１つのプロセッサと、プロセッサに前記操作を実行するよう指示するためのコンピュータ・コードを保持するための関連のメモリ・サブシステムとを備えて使用できる。図４が、プロセッサ４０１と、メモリ・サブシステム４０２と、入力／出力装置（キーボード、マウス、ディスプレイ、など）などの周辺機器４０３と、おそらくは計算を補助するためのコプロセッサ４０４と、ネットワークインターフェイス装置４０５とを、すべてバス４００で相互接続して有しているそのようなコンピュータ・システムの基本的なアーキテクチャを示している。メモリ・サブシステムは、アクセス待ちの時間が長くなる順番で、キャッシュメモリ、主メモリ、およびハードディスク・ドライブなどの恒久的保存メモリを、最適に備えている。計算の激しさの程度を考慮し、コンピュータ・システムが、計算を並列に実行するために複数の関連のメモリ・サブシステムとともに複数のプロセッサを備えてもよく、あるいは種々のコンピュータ要素を図４によって示されるような従来からのコンピュータ・アーキテクチャにてバスで接続して有するのではなく、コンピュータ・システムが、ネットワークによって互いに接続される複数のプロセッサおよび複数のメモリ・サブシステムによって形成されてもよいことを、理解すべきである。
【００９９】
以下の例は、本発明を説明するためのものであるが、いかなるかたちでも本発明を限定するものではない。
【０１００】
（実施例１）
（従来の方法および本発明の方法の計算速度の比較）
動的残差ρｕが、以下のとおり原子間力Ｆから計算される。
【０１０１】
【数４２】
ρｕの微分を、混合デカルト‐ねじり導関数行列
【０１０２】
【数４３】
を使用する
【０１０３】
【数４４】
を使用することによって、従来の方法を用いて得ることができる。
【０１０４】
しかしながら、本発明の方法を用いると、今や動的残差の微分を
【０１０５】
【数４５】
として計算できる。
【０１０６】
本明細書に詳細を説明した方法を、以下のように従来のからのアプローチに組み合わせた。便宜上、原子間力を２つの成分に、それらの導関数が２つの別個の項
【０１０７】
【数４６】
で計算されるような方法で分割した。
【０１０８】
【数４７】
したがって、ρｕの微分は、
【０１０９】
【数４８】
として記載され、第１が本発明の実施の形態によるデカルト・ヘシアンであり、第２が既に使用されている混合デカルト‐ねじりヘシアンである。
【０１１０】
（基本アルゴリズム）
以下の表現
【０１１１】
【数４９】
を、本明細書に記載の方法に従って計算し、既存の方法によって計算された
【０１１２】
【数５０】
の残りを加えた。
【０１１３】
これは、
【０１１４】
【数５１】
と書き換えられ、
１．「物体」ヘシアン
【０１１５】
【数５２】
の計算、および
２．以下の並びの演算子を使用するρｕ’の計算
【０１１６】
【数５３】
を可能にした。
【０１１７】
（省メモリのアルゴリズム）
省メモリのアルゴリズムは、Ｏ（Ｎ）中間メモリのみを保存する。動的残差ｑ‐導関数の表現を、以下のとおり書くことができる。
【０１１８】
【数５４】
行列ＷおよびΩを、効率的な要素ごとの方法で計算できる。すなわち、
子および親のインデックス間の差ｍを最小にするため物体を並べ直し、
行列Ωのｍ個の列のためにメモリを予約し、
【０１１９】
【数５５】
ここで、
【０１２０】
【数５６】
上記において、アインシュタインの総和規約が使用される。行列Ｗの要素が、必要に応じ
【０１２１】
【数５７】
から計算される。
【０１２２】
（実施例２）
（ヤコビ行列の計算）
ヤコビ行列が、４つのブロックで構成される。
【０１２３】
【数５８】
左下のブロック以外のすべては、計算が単純である。本明細書に記載の方法によれば、左下のブロックを、質量行列Ｍおよび動的残差ρｕから計算できる。
【０１２４】
【数５９】
最後の項は、本発明において説明したとおり、動的残差の一般化座標に関する導関数を含んでいる。
【０１２５】
（実施例３）
（剛性行列の計算）
系が、ｎｑ個の一般化座標およびｎｕ個の一般化速度を有するとする。制約または４元ジョイント（ｑｕａｔｅｒｎｉｏｎ ｊｏｉｎｔ）が存在しないならば、ｎｑ＝ｎｕである。この場合、一般化座標の独立組に関する最大階数（非特異）の剛性行列を、直接計算することができる。
【０１２６】
【数６０】
これは、ｎｕ×ｎｕの行列である。
【０１２７】
一方、４元および／または制約が存在する場合には、ｎｑ＞ｎｕである。この場合、剛性行列を２つのステップで計算することができ、最初のステップはヤコビアンの計算の部分集合である。
１．（従属）一般化座標の完全な組に関する特異剛性行列が、動的残差ρｕの導関数として計算される。
【０１２８】
【数６１】
これは、ｎｑ×ｎｑの行列である。
２．次いで、非特異剛性行列が、従属一般化座標を除くことによって、この行列から計算される。これは、いくつかある標準的なアプローチのいずれかを使用して行なうことができる。例えば、制約行列の特異値分解を使用することができる。Ｇが制約行列を表わしている場合、ＧＺ＝０を解くことによってＧのヌル空間Ｚを計算する。次いで、最大階数の剛性行列を、
【０１２９】
【数６２】
として計算できる。
剛性行列の用途の１つは、この技術分野においてよく知られているとおり、正規モードの計算である。
【０１３０】
（結果）
以上の実装例を、既存の方法および本発明の方法を使用するリゾチーム（＋／−オメガ）およびトリプシン（＋オメガ）の
【０１３１】
【数６３】
の計算に使用した。結果は、下記の表１にまとめられている。
【０１３２】
【化４】
速度（ＰＡ）は、前記従来の方法を使用した計算速度（単位は秒）であり、速度（ＭＩ）は、本発明の方法を使用した計算速度（単位は秒）である。本発明の方法が、これまで使用されている方法に比べて
【０１３３】
【数６４】
の計算速度を大きく向上させることができることを、理解できるであろう。
【０１３４】
表２は、基本アルゴリズム、および高速ヤコビアン・アルゴリズムの省メモリ（ＬＭ）実装からの結果を示している。
【０１３５】
【化５】
（拡大可能性）
デカルト・ヘシアン行列
【０１３６】
【数６５】
は大きくなる可能性があり、最終的に３２ビットのアドレス空間を超える可能性もある。この問題を解決するため、行列を部分に構成できる。このアルゴリズムは、ヘシアン行列の行について動作させ、Ｗを行ごとに構成し、ρｕ’の対応する列を構成することによって拡大可能である。
【０１３７】
以上、本発明の実施の形態を完全に説明したが、さまざまな変更、代案、および均等物が可能でありかつ利用可能であることは、明らかである。したがって、上記説明は、添付の特許請求の範囲の計量および境界によって定められる本発明の技術的範囲を限定するものとして理解すべきではない。
【図面の簡単な説明】
【０１３８】
【図１】多体系のツリー構造を示している。
【図２】図１に示した多体系についての行列εφである。
【図３Ａ】２つの凍結体間の相互作用を象徴的に示している。
【図３Ｂ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｃ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｄ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｅ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図３Ｆ】図３Ａの２つの凍結体間の相互作用の計算におけるステップを示している。
【図４】本発明の方法の実行に有用であるコンピュータ・システムの図である。
【特許請求の範囲】
【請求項１】
広く適用可能な多体シミュレーションを実行する方法であって、
動的残差の導関数（ＤＤＲ）をオーダ（Ｎ２）で一般化座標に関して計算することと、
該多体シミュレーションにおいて該導関数を使用することと
を包含する、方法。
【請求項２】
前記計算することが、
（ｉ）前記ＤＤＲを、多体演算子とその転置とを含む第１項と、１つ以上の追加の項とに分解することと、
（ｉｉ）高速演算子のインプリメンテーションを用いて該第１項を評価することと
を包含する、請求項１に記載の方法。
【請求項３】
前記ＤＤＲが、
【数１】
の形式をとる、請求項２に記載の方法。
【請求項４】
前記評価することが、物体ヘシアン（Ｗ）を特定することと、前記第１項をオーダ（Ｎ２）の表現に再定式化することとを包含する、請求項２に記載の方法。
【請求項５】
前記評価することが、前記物体ヘシアンの要素を局所的に決定することを包含する、請求項４に記載の方法。
【請求項６】
前記要素のそれぞれは１回だけ決定される、請求項５に記載の方法。
【請求項７】
前記評価することが、オメガ（Ω）についての表現を規定することを包含する、請求項２に記載の方法。
【請求項８】
前記Ωについての表現が、
【数２】
の形式をとる、請求項７に記載の方法。
【請求項９】
前記第１項が、剛性行列の計算において使用される、請求項２に記載の方法。
【請求項１０】
前記第１項が、ヤコビ行列の計算において使用される、請求項２に記載の方法。
【請求項１１】
前記多体シミュレーションが分子のシミュレーションである、請求項１に記載の方法。
【請求項１２】
前記多体シミュレーションがコンピュータを使用して実行される、請求項１に記載の方法。
【請求項１３】
前記一般化座標がねじり角座標である、請求項１に記載の方法。
【請求項１４】
多体系によるシミュレーションを実行する方法であって、
【数３】
を評価することによって、ねじりヤコビ行列を計算することと、
該ねじりヤコビ行列を使用することによって、空間における該多体系の物体の特徴を決定することと
を包含する、方法。
【請求項１５】
前記特徴が、前記物体の位置および速度を含む、請求項１４に記載の方法。
【請求項１６】
前記シミュレーションが分子のシミュレーションである、請求項１４に記載の方法。
【請求項１７】
前記使用することが、
【数４】
の形式を有するオメガ（Ω）についての表現を規定することを包含する、請求項１４に記載の方法。
【請求項１８】
多体系における物体の特徴を決定する方法であって、
【数５】
を評価することによって、剛性行列を計算することと、
該剛性行列を使用することによって、該特徴を決定することと
を包含する、方法。
【請求項１９】
前記多体系が、分子のシミュレーションにおける１つ以上の分子を表わす、請求項１８に記載の方法。
【請求項２０】
前記特徴が、前記系の正規モードを含む、請求項１８に記載の方法。
【請求項１】
広く適用可能な多体シミュレーションを実行する方法であって、
動的残差の導関数（ＤＤＲ）をオーダ（Ｎ２）で一般化座標に関して計算することと、
該多体シミュレーションにおいて該導関数を使用することと
を包含する、方法。
【請求項２】
前記計算することが、
（ｉ）前記ＤＤＲを、多体演算子とその転置とを含む第１項と、１つ以上の追加の項とに分解することと、
（ｉｉ）高速演算子のインプリメンテーションを用いて該第１項を評価することと
を包含する、請求項１に記載の方法。
【請求項３】
前記ＤＤＲが、
【数１】
の形式をとる、請求項２に記載の方法。
【請求項４】
前記評価することが、物体ヘシアン（Ｗ）を特定することと、前記第１項をオーダ（Ｎ２）の表現に再定式化することとを包含する、請求項２に記載の方法。
【請求項５】
前記評価することが、前記物体ヘシアンの要素を局所的に決定することを包含する、請求項４に記載の方法。
【請求項６】
前記要素のそれぞれは１回だけ決定される、請求項５に記載の方法。
【請求項７】
前記評価することが、オメガ（Ω）についての表現を規定することを包含する、請求項２に記載の方法。
【請求項８】
前記Ωについての表現が、
【数２】
の形式をとる、請求項７に記載の方法。
【請求項９】
前記第１項が、剛性行列の計算において使用される、請求項２に記載の方法。
【請求項１０】
前記第１項が、ヤコビ行列の計算において使用される、請求項２に記載の方法。
【請求項１１】
前記多体シミュレーションが分子のシミュレーションである、請求項１に記載の方法。
【請求項１２】
前記多体シミュレーションがコンピュータを使用して実行される、請求項１に記載の方法。
【請求項１３】
前記一般化座標がねじり角座標である、請求項１に記載の方法。
【請求項１４】
多体系によるシミュレーションを実行する方法であって、
【数３】
を評価することによって、ねじりヤコビ行列を計算することと、
該ねじりヤコビ行列を使用することによって、空間における該多体系の物体の特徴を決定することと
を包含する、方法。
【請求項１５】
前記特徴が、前記物体の位置および速度を含む、請求項１４に記載の方法。
【請求項１６】
前記シミュレーションが分子のシミュレーションである、請求項１４に記載の方法。
【請求項１７】
前記使用することが、
【数４】
の形式を有するオメガ（Ω）についての表現を規定することを包含する、請求項１４に記載の方法。
【請求項１８】
多体系における物体の特徴を決定する方法であって、
【数５】
を評価することによって、剛性行列を計算することと、
該剛性行列を使用することによって、該特徴を決定することと
を包含する、方法。
【請求項１９】
前記多体系が、分子のシミュレーションにおける１つ以上の分子を表わす、請求項１８に記載の方法。
【請求項２０】
前記特徴が、前記系の正規モードを含む、請求項１８に記載の方法。
【図１】
【図２】
【図３Ａ】
【図３Ｂ】
【図３Ｃ】
【図３Ｄ】
【図３Ｅ】
【図３Ｆ】
【図４】
【図２】
【図３Ａ】
【図３Ｂ】
【図３Ｃ】
【図３Ｄ】
【図３Ｅ】
【図３Ｆ】
【図４】
【公表番号】特表２００７−５１６５０８（Ｐ２００７−５１６５０８Ａ）
【公表日】平成１９年６月２１日（２００７．６．２１）
【国際特許分類】
【出願番号】特願２００６−５３３６６０（Ｐ２００６−５３３６６０）
【出願日】平成１６年６月９日（２００４．６．９）
【国際出願番号】ＰＣＴ／ＵＳ２００４／０１８３６８
【国際公開番号】ＷＯ２００４／１１１９１４
【国際公開日】平成１６年１２月２３日（２００４．１２．２３）
【出願人】（５０５４５４１４４）ルーカス ファーマシューティカルズ， インコーポレイテッド (3)
【Ｆターム（参考）】
【公表日】平成１９年６月２１日（２００７．６．２１）
【国際特許分類】
【出願日】平成１６年６月９日（２００４．６．９）
【国際出願番号】ＰＣＴ／ＵＳ２００４／０１８３６８
【国際公開番号】ＷＯ２００４／１１１９１４
【国際公開日】平成１６年１２月２３日（２００４．１２．２３）
【出願人】（５０５４５４１４４）ルーカス ファーマシューティカルズ， インコーポレイテッド (3)
【Ｆターム（参考）】
[ Back to top ]