(57)【要約】物理系の動力学についての計算における陰的積分法で使用される解析的ヤコビアンを得る方法。この方法を用いて、数値的ヤコビアンの少なくとも２倍の有効桁数を有するヤコビアンが計算され得る。これはまた、運動方程式の非線形ステージ方程式を解くためにより少ない反復回数が必要とされるため、より効率的な陰的積分法を生じる。計算におけるスピードの向上は、分子モデリングにおいて非常に有用である。本発明はまた、ねじれ角、剛体多体系によってモデリングされ得る、分子の挙動、または任意の物理系をシミュレートするための、コンピュータコードを提供する。陰的積分器を使用する、コンピュータコードにおけるモジュールは、上記のような解析的ヤコビアンを含む。【特許請求の範囲】
【請求項１】
分子の挙動をモデリングする方法であって、該方法は、以下：
該分子について、ねじれ角、剛体多体系モデルを選択する工程であって、該モデルが、運動方程式を有する、工程；
陰的積分器を選択する工程；および
解析的ヤコビアンを生成して、該陰的積分器が該運動方程式を積分して、該分子の挙動の計算を得る、工程、
を包含する、方法。
【請求項２】
前記解析的ヤコビアンが、前記運動方程式の残差形式の解析的ヤコビアンから導出される、請求項１に記載の方法。
【請求項３】
請求項２に記載の方法であって、前記解析的ヤコビアンＪが、以下：
【表１】
2004527027.tif 000198
を含み、ここで、ｑは、一般化座標であり、ｕは一般化速度であり、Ｗはジョイントマップ行列であり、そしてＭは質量行列であり、そしてρｕは、前記運動方程式の動力学的残差であり、そしてｚは、−Ｍ−１ρｕ（ｑ、ｕ、０）である、方法。
【請求項４】
前記陰的積分器を選択する工程がＬ安定性積分器を含む、請求項３に記載の方法。
【請求項５】
物理系の挙動をシミュレートする方法であって、該方法は、ねじれ角、剛体モデルを用いて該物理系をモデリングする工程であって、該モデルは、運動方程式を有する工程；および
陰的積分器を用いて該運動方程式を積分する工程；
を包含し、該陰的積分器が解析的ヤコビアンを有し、該物理系の該挙動の計算を得る、方法。
【請求項６】
前記解析的ヤコビアンが前記運動方程式の残差形式の解析的ヤコビアンから導出される、請求項５に記載の方法。
【請求項７】
請求項６に記載の方法であって、前記解析的ヤコビアンＪは、以下：
【表２】
2004527027.tif 000199
を含み、ここで、ｑは、一般化座標であり、ｕは一般化速度であり、Ｗはジョイントマップ行列であり、そしてＭは質量行列であり、そしてρｕは、前記運動方程式の動力学的残差であり、そしてｚは、−Ｍ−１ρｕ（ｑ、ｕ、０）である、方法。
【請求項８】
前記陰的積分器がＬ安定性積分器を含む、請求項７に記載の方法。
【請求項９】
分子の挙動をシミュレーションするためのコンピュータコードであって、該コードが、以下：
該分子のねじれ角、剛体多体系モデルのための第１のモジュールであって、該モデルが、運動方程式を有する、第１のモジュール；および
陰的積分器が解析的ヤコビアンを使用して該運動方程式を積分し、該分子の挙動の計算を得る、第２のモジュール；
を備える、コンピュータコード。
【請求項１０】
前記解析的ヤコビアンが前記運動方程式の残差形式の解析的ヤコビアンから導出される、請求項９に記載のコンピュータコード。
【請求項１１】
請求項１０に記載のコンピュータコードであって、前記解析的ヤコビアンＪが、以下：
【表３】
2004527027.tif 000200
を含み、ここで、ｑは、一般化座標であり、ｕは一般化速度であり、Ｗはジョイントマップ行列であり、そしてＭは、質量行列であり、そしてρｕは、前記運動方程式の動力学的残差であり、そしてｚは、−Ｍ−１ρｕ（ｑ、ｕ、０）である、コンピュータコード。
【請求項１２】
前記陰的積分器がＬ安定性積分器である、請求項１１に記載のコンピュータコード。
【請求項１３】
物理系の挙動をシミュレーションするためのコンピュータコードであって、該コードは、以下：
該系のねじれ角、剛体多体系モデルのための第１のモジュールであって、該モデルが、運動方程式を有する、第１のモジュール；および
陰的積分器が解析的ヤコビアンを使用して該運動方程式を積分し、該分子の挙動の計算を得る、第２のモジュール；
を備える、コンピュータコード。
【請求項１４】
前記解析的ヤコビアンが前記運動方程式の残差形式の解析的ヤコビアンから導出される、請求項１３に記載のコンピュータコード。
【請求項１５】
請求項１４に記載のコンピュータコードであって、前記解析的ヤコビアンＪが、以下：
【表４】
2004527027.tif 000201
を含み、ここで、ｑは、一般化座標であり、ｕは一般化速度であり、Ｗはジョイントマップ行列であり、そしてＭは、質量行列であり、そしてρｕは、前記運動方程式の動力学的残差であり、そしてｚは、−Ｍ−１ρｕ（ｑ、ｕ、０）である、コンピュータコード。
【請求項１６】
前記陰的積分器がＬ安定性積分器を含む、請求項１５に記載のコンピュータコード。
【発明の詳細な説明】
【０００１】
（関連出願の引用参照）
本出願は、仮特許出願第６０／２４５，７３０号（２０００年１１月２日に出願）；および、さらに、第６０／２４５，６８８号（２０００年１１月２日に出願）；第６０／２４５，７３１号（２０００年１１月２日に出願）；ならびに第６０／２４５，７３４号（２０００年１１月２日に出願）の優先出願日の利益を享受し、これらの出願の全てを本明細書中で参考として援用する。
【０００２】
（発明の背景）
本発明は分子モデリングの分野に関し、より詳細には、高分子の動力学的モデリングおよび静的解析のためのコンピュータにインプリメントされた方法に関する。
【０００３】
分子モデリングの分野は、コンピュータによって多くの複雑な分子系の、その運動を首尾よくシミュレート（分子動力学法、すなわち（ＭＤ））し、エネルギー極小状態または静止状態（静的解析）を決定してきた。代表的な分子モデリング用途としては、酵素−リガンドドッキング、分子拡散、反応経路、相転移、およびタンパク質のフォールディングの研究が挙げられる。生命科学、ならびに、製薬、ポリマー、および化学産業の研究者が、複雑な分子の化学的過程の性質を理解し、それに従って新規の薬物および材料を設計するために、それらの技術を使用し始めている。本来、これらのツールの受容は、現実のものを表現するときの結果の精度、モデルリングされ得るべき分子系のサイズおよび複雑さ、および解が得られる速度を含む、いくつかの要因に基づく。多くの計算の精度は実験と比較され、一般的に特定の境界内では適切であると理解される。しかし、従来技術におけるこれらのツールの使用は、適度な大きさの分子または分子系でさえも、それらをモデリングして有用である十分な長さの分子の時間履歴を得るために膨大な計算パワーを、必要としてきた。
【０００４】
時間積分を含む分子モデリング課題についての計算の複雑さの２つの原因が存在する。
【０００５】
１．特定の分子モデルであって、構成原子の位置、速度、および質量特性、構成原子間の原子間力、ならびにそれらの原子と構成原子の周囲環境との間の相互作用を記述するために使用される、分子モデル；ならびに
２．経時的にモデルを前進させるために使用される、特定の数値的方法。
時刻は、最終時刻に達するまで、非常に短い間隔（これは、時間刻みと呼称される）分だけ、繰り返し進められる。
【０００６】
実質的な研究は、例えば、以下のように分子モデルについての計算負荷を低減することにおいて完結した：剛体仮定条件で高次のモードを拘束することによって、モデルの複雑性の低減すること、剛体性または伸縮性の下部構造分解（ｓｕｂｓｔｒｕｃｔｕｒｉｎｇ）によってモデルを単純化すること、Ｎ次動力学、有効性のある非明示的な溶媒モデル、および力場モデルについての多重極法（ｍｕｌｔｉｐｏｌｅ  ｍｅｔｈｏｄ）（例えば、市販のＭＢＯ（Ｎ）Ｄソフトウエアパッケージに関する米国特許第５，４２４，９６３号を参照のこと））。同時係属出願である、「ＭＥＴＨＯＤ  ＦＯＲ  ＬＡＲＧＥ  ＴＩＭＥＳＴＥＰＳ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題する、米国特許出願番号＿＿＿＿、および「ＭＥＴＨＯＤ  ＦＯＲ  ＲＥＳＩＤＵＡＬ  ＦＯＲＭ  ＩＮ  ＭＯＬＥＣＵＬＥＲ  ＭＯＤＥＬＩＮＧ」と題する、米国特許出願番号＿＿＿＿（これらの両方は、同日付けで出願され、そして、上で引用した仮特許出願の優先権を主張し、本譲渡受人に譲渡され、そしてこれらは、本明細書中において参考として援用される）によって、分子モデリングおよび数値的方法におけるさらなる改善が記載される。
【０００７】
分子シミュレーションが非常に遅い主要な原因は、現在の数値的方法が非常に小さな時間刻み（代表的に、１〜１０フェムト秒（１０−１５〜１０−１４秒））を必要することである。採用された各時間刻みは、特定の分子モデルの新しい状態（各原子についての位置および運動）の計算を必要とし、次いで、力の新しいセットの計算は、この新しい状態から生じる。例えば、高分子の複雑な挙動（例えば、タンパク質のフォールディング）の分子動力学シミュレーションは、代表的には、少なくとも１マイクロ秒〜数秒（数分でさえも）の時間範囲に及ぶ必要がある。現在、通常用いられている技術によって、これは、コンピュータシミュレーションにおいて１０９〜１０１６の時間刻みをとる必要性を生じる。この状態、および特に力の、刻み当たりの計算は、問題の規模が大きくなるにつれて非常に高くなる。今日利用可能な最高速のコンピュータを用いたとしても、数ヶ月、数年、または数百年のコンピュータ処理時間が、適度な大きさの系であってもこのような問題を解決するために必要とされる。
【０００８】
これまで、これらの小さな時間刻みは、分子結合の振動において見出される非常に高い振動数の存在下で精度を維持するための要求の、不可避な必要条件であることが分子動力学者によって広く信じられてきた。例えば、Ｌｅａｃｈによる、Ｍｏｌｅｃｕｌａｒ  Ｍｏｄｅｌｌｉｎｇ  Ｐｒｉｎｃｉｐｌｅｓ  ａｎｄ  Ａｐｐｌｉｃａｔｉｏｎｓ，１９９６，ｐ．３２８；Ｂｅｒｅｎｄｓｅｎによる、Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ：Ｃｈａｌｌｅｎｇｅｓ，Ｍｅｔｈｏｄｓ，Ｉｄｅａｓ；Ｄｅｕｆｌｈａｒｄら（編）による、Ｓｐｒｉｎｇｅｒ，１９９９，ｐｇ．１８；Ｒａｐａｐｏｒｔによる、Ｔｈｅ  Ａｒｔ  ｏｆ  Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ  Ｓｉｍｕｌａｔｉｏｎ，Ｃａｍｂｒｉｄｇｅ，１９９５（訂正と共に再版された１９９８年、ｐ．５７）；および米国特許第５，４２４，９６３号を参照のこと。
【０００９】
しかし、この常識的な確信は誤りである。数値解析のうちのコンピュータサイエンス下位区分によって、高振動数が存在するが、それがほとんど関心事とならない問題に対する、数値積分法の広範な理論が生まれてきた。これらの問題は、「スティッフ（ｓｔｉｆｆ）」問題と呼称される（例えば、ＨａｉｒｅｒおよびＷａｎｎｅｒ、Ｓｏｌｖｉｎｇ  Ｏｒｄｉｎａｒｙ  Ｄｉｆｆｅｒｅｎｔｉａｌ  Ｅｑｕａｔｉｏｎｓ  ＩＩ：Ｓｔｉｆｆ  ａｎｄ  Ｄｉｆｆｅｒｅｎｔｉａｌ−Ａｌｇｅｂｒａｉｃ  Ｐｒｏｂｌｅｍｓ，第２版、Ｓｐｒｉｎｇｅｒ，１９９６を参照のこと）。これらの場合においては、時間刻みを制限するものは、これらの積分法の安定性であって、その必要とされる解の精度ではない。無条件な安定性を有する積分法が存在し、このことは、理論上、これらの積分法が任意の大きな時間刻みを採用することができることを意味する。このような方法は、「Ｌ安定性」と呼称される数学的な特性を有する。所謂「陰的」な積分法のみが、Ｌ安定性であり得るが、ほとんどの陰的積分法が、現実的に「Ｌ安定性」ではない。上で引用した同時係属出願である、「ＭＥＴＨＯＤ  ＦＯＲ  ＬＡＲＧＥ  ＴＩＭＥＳＴＥＰＳ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題する、米国特許出願番号＿＿＿＿は、陰的であり、特にＬ安定性である積分法の使用を取り扱う。
【００１０】
本発明は、陰的積分器（特にＬ安定性積分器）が大きな時間刻みをとることを可能にするための別の重要な局面、すなわち、「ヤコビアン」と呼称される陰的積分法に必要とされる成分を計算するためのより正確な方法を取り扱う。
【００１１】
しかし、ＭＤシミュレーションに対する高振動数の分子振動の同じ問題が、ヤコビアンの計算についても問題を生じる。例えば、２つの原子の電子軌道関数が重なるために生成される斥力性のファンデルワールス力は、分子内で説明されなければならない。この量子力学的効果は、しばしば、所謂、Ｌｅｎｎａｒｄ−Ｊｏｎｅｓポテンシャル（Ｒａｐａｐｏｒｔ，ｏｐ．ｃｉｔ．）によって、近似され、このポテンシャルは、この斥力を、１／ｒ１３に比例するものとしてモデリングする。ここで、ｒは、原子の中心間距離である。これは、非常にスティッフ（ｓｔｉｆｆ）である原子間力であり、これは、分子動力学（ＭＤ）シミュレーションの特徴であり、そして、これは、シミュレーションの間に、時間を前進させるために使用される任意の数値積分スキームについて特別な困難を引き起こす。結果として、かつ、上述したように、殆どの従来技術のＭＤ積分スキームは、これらの非常にスティッフな斥力により生成される高振動数によって限定される時間刻みを有する。そして、特に、この原子間力のスティッフ性は、特定の必要とされるヤコビ行列を形成する困難を顕著に増強する。このようなヤコビ行列は、直上で引用した同時係属中の出願に記載されるように、任意の安定性陰的積分スキームの必要不可欠な構成要素である。
【００１２】
全ての汎用目的のシミュレーションコードは、以下のようにヤコビアンを数値的に計算するルーチンを提供する。
【００１３】
【数１】
2004527027.tif 000002
を積分するため所与の式について、所望のヤコビアンは、
【００１４】
【数２】
2004527027.tif 000003
であり、以下のように数値的計算される：
【００１５】
【数３】
2004527027.tif 000004
。
【００１６】
ここで、摂動Δｙは、十分な結果を与えるように選択される必要があり、そして特にスティッフ（ｓｔｉｆｆ）な系について選択することがで有り得ることに留意のこと。対照的に、解析的ヤコビアンは、直接的、すなわち、この場合アルゴリズム的に、所望の導関数の方程式について解くことによって、計算される。
【００１７】
線形化されたモデルは、通常は、単純な系に対して解析的に生成される。このような線形化は、通常、平衡解の周囲で実行される。ＡＣＳＬ（Ａｄｖａｎｃｅｄ  Ｃｏｎｔｉｎｕｏｕｓ  Ｓｉｍｕｌａｔｉｏｎ  Ｌａｎｇｕａｇｅ）（ＡＣＳＬ  Ｒｅｆｅｒｅｎｃｅ  Ｇｕｉｄｅ，  Ｍｉｔｃｈｅｌｌ  Ｇａｕｔｈｉｅｒ  ａｎｄ  Ａｓｓｏｃｉａｔｅｓ  ，１９９６）およびＳＰＩＣＥ（サーキットシミュレーションパッケージ）（Ｒ．Ｋｉｅｌｋｏｗｓｋｉ，Ｉｎｓｉｄｅ  ＳＰＩＣＥ，ＭｃＧｒａｗ−Ｈｉｌｌ，１９８８）のようなパッケージにおいて、線形化の前に平衡解析を行うことは、一般的である。このような線形化は、数値的に実行される。
【００１８】
対照的に、本発明のヤコビアンは、微分方程式の瞬間の解（非平衡）をについての線形化を表現し、そして、解析的に生成される。解析的ヤコビアンを生成するための別の従来技術のアプローチは、任意の方程式を記号的に微分し得るＡＤＩＦＯＲ（Ａｕｔｏｍａｔｉｃ  Ｄｉｆｆｅｒｅｎｔｉａｔｉｏｎ  ｏｆ  Ｆｏｒｔｒａｎ）（Ｃ．Ｂｉｓｃｈｏｆら、ＡＤＩＦＯＲ  ２．０  Ｕｓｅｒｓ’Ｇｕｉｄｅ  Ａｒｇｏｎｎｅ  Ｎａｔｉｏｎａｌ  Ｌａｂｏｒａｔｏｒｙ，１９９８）のような自動化微分ツールを使用することであることに留意するべきである。これらのツールは、実際問題として、本発明をインプリメントするために使用され得る。しかし、この方程式の構造は、その計算を最小化し、丸めおよび項の相殺に起因する誤差を避け、そして、解かれる問題のサイズを限定し得る「方程式スウェル（ｅｑｕａｔｉｏｎ  ｓｗｅｌｌ）」を回避するように適切に活用されなければならない。
【００１９】
むしろ、本発明は、分子モデルの運動方程式の有効な陰的積分法（Ｌ安定性積分器を含む）についての解析的ヤコビアンの計算を可能にする。
【００２０】
（発明の要旨）
本発明は、分子の挙動をモデリングする方法を提供する。この方法は、その分子について、ねじれ角、剛体多体系モデルを選択する工程であって、このモデルが、運動方程式を有する工程；陰的積分器を選択する工程；および解析的ヤコビアンを生成して、陰的積分器が運動方程式を積分して分子の挙動の計算を得る、工程を、有する。
解析的ヤコビアンＪは、運動方程式の残差形式の解析的ヤコビアンから導出され、そして、これは、以下：
【００２１】
【数４】
2004527027.tif 000005
のように記述される。ここで、ｑは一般化座標、ｕは一般化速度、Ｗはジョイントマップ行列であり、そしてＭは質量行列であり、そしてρｕは、上記運動方程式の動力学的残差であり、そしてｚは、−Ｍ−１ρｕ（ｑ、ｕ、０）である。この方法はまた、ねじれ角、剛体多体系によってモデリングされ得る任意の物理系について使用され得る。
【００２２】
本発明はまた、ねじれ角、剛体多体系によってモデリングされ得る、分子の挙動、または任意の物理系をシミュレートするための、コンピュータコードを提供する。陰的積分器を使用する、コンピュータコードにおけるモジュールは、上記のような解析的ヤコビアンを含む。
【００２３】
（特定の実施形態の説明）
（本発明の一般的な説明）
モデリングした分子系のシミュレーションにおいて時間を前進させるために使用される数値的方法は、積分器、または積分法と呼ばれる。この積分は、この分子系の支配的な運動方程式を差分化することによって、前進する。陰的積分器の場合、この工程は、結合型非線形代数方程式（「ステージ（ｓｔａｇｅ）」方程式）の集合を生じ、そして、これらの方程式の解は、次の時間刻みにおけるこの分子系の状態ベクトルを生成する。
【００２４】
本発明は、このステージ方程式の解法を支援する。原子間力が、短距離にわたって非常に迅速に変動するので、原子間のトラジェクトリーを正確に伝達することが難しい。原子の位置のわずかな誤差が、このステージ方程式を満たす際に大きな誤差を生じる。このヤコビアンがこのステージ方程式を反復的に解くために使用されるので、不正確なヤコビアンが正確な解からかけ離れた試験的な解（ｔｒｉａｌ  ｓｏｌｕｔｉｏｎ）を生じる。このことによって、この試験的な解の撤回に至り、そしてこのシミュレーションを妨げる。
【００２５】
数値的ヤコビアンは、その有効数字の半分のみで正確であり得る。解析的ヤコビアンは、しばしば、その最後の桁以外の全てにおいて正確であり得る。この結果の利点は、この積分器が、特定の間隔のシミュレーションを実施するためにずっと少ない時間刻みを取ることが可能であり、この積分法の理論的な安定性の完全な享受が可能となる。
【００２６】
このヤコビアンを生成することの容易性または困難性は、この支配的な方程式を生成するために使用される定式化に決定的に依存する。例えば、グローバル（ｇｌｏｂａｌ）デカルト定式化は、非常に限定された陽的座標依存性を有する方程式を生じる。このような定式化に対する解析的ヤコビアンの生成は周知である。他方で、内部座標（これにおいて、各分子サブユニットの位置は、それより前のサブユニットの位置に対して表現される）を独立変数として使用することは、大きな利点を有する。この方法は、内部座標の任意の選択でモデリングされた任意の分子について最も価値があり、そして、特に、残基間の「ねじれ角」を内部座標として用いる、タンパク質のモデルまたは他の高分子とともに使用される場合に、価値がある。ねじれ角で定式化された系についての解析的ヤコビアン生成するためのアルゴリズムは、考案するのが極めて困難である。しかし、本発明は、この課題を達成する。
【００２７】
本発明は、以下のような、一見、多項式アルゴリズムが与えられない問題に取り組む：内部座標、および特にねじれ角を使用する定式化のための解析的ヤコビアンの生成（これは、一般的に非実用的であると考えられている）。数値的ヤコビアンよりも正確であることに加え、この解析的ヤコビアンはまた、生成されることが（コンピュータパワーにおいて）安価である。本発明はまた、状態微分のヤコビアンが、計算トルク法（ｃｏｍｐｕｔｅｄ  ｔｏｒｑｕｅ  ｍｅｔｈｏｄ）のヤコビアンに逆行列を適用することによって計算されるという重要な結果を認識する。この結果は、コンピュータ処理時間およびアルゴリズムを構築する労力の有意な節約となる。
【００２８】
さらに、ねじれ角、内部座標を使用する多体系定式化についての解析的ヤコビアンを生成するこの方法は、一般的なＭＢＳの文献においてこれまで見出されていない。本発明は、任意のねじれ角ＭＢＳ定式化のために使用され得、このねじれ角ＭＢＳは、以下を含むが必ずしもこれらに限定されない、分子シミュレーション以外の他の多数の学問分野にも適用され得る：機械システム、ロボットシステム、乗物システム、またはヒンジ接続剛体のセットとして記述され得る他の任意のシステム。
【００２９】
（説明の概観）
好ましい実施形態は、いくつかの節に分割される。節の第１のセットは、ＭＤシミュレーションアーキテクチャ、多体系（ＭＢＳ）の定義、およびＭＢＳ方程式の残差形式を、引き続き説明のために記載する。節の次のセットは、ヤコビアンの定義、陰的積分法におけるヤコビアンの役割、および残差形式を用いた解析的ヤコビアンの計算法を、議論する。数値的ヤコビアンに対する解析的ヤコビアンの優れた精度および能力もまた示される。解析的ヤコビアンの計算法におけるさらなる有効性が議論され、特に、剛体ＭＢＳをこのヤコビ行列のサイズを「減縮」するために利用すること、およびこのＭＢＳのトポロジー構造を利用して不必要な計算を取り除くことを議論する。
【００３０】
常微分方程式（ＯＤＥ）を解くために、従来技術のほとんどは、直接形式、すなわち、以下：
【００３１】
【数５】
2004527027.tif 000006
（ここで、
【００３２】
【数６】
2004527027.tif 000007
は、
【００３３】
【数７】
2004527027.tif 000008
を意味する）
で表現された方程式を使用してきた。生体分子系についての運動方程式は、この形式（そして、これは直接形式と呼ばれる）に投入され得る。分子モデリングにおいて、本発明に対する全ての先行技術は、この直接形式を使用してきた。すなわち、
【００３４】
【数８】
2004527027.tif 000009
であり、ここで、ｑおよびｕは、それぞれ、一般化座標および一般化速度であり、その結果、従来のＯＤＥの解法が適用され得る。しかし、これは、次数（Ｎ）〜次数（Ｎ３）の次数（オーダー）の浮動小数点演算（これは、使用されるアルゴリズムに依存し、ここで、Ｎはこの系の自由度である）を費やすＭ（この系の質量を表す）の逆行列を必要とする。なぜならば、この方程式の固有の形式が、この一般化速度の微分の間の慣性カップリング（ｉｎｅｒｔｉａｌ  ｃｏｕｐｌｉｎｇ）を生じるためである。すなわち、この方程式は、以下：
【００３５】
【数９】
2004527027.tif 000010
の形式で、固有に生成され、ここで、質量行列Ｍは、一般化座標ｑに陽的に依存しており、すなわち、Ｍ＝Ｍ（ｑ）である。この事実は、この状態微分が、経時的に積分器によって必要とされる各時刻での質量行列を形成し、そして効率的にこの質量行列を因数分解することを要求している。
【００３６】
一般化ジョイントマップ行列（ｇｅｎｅｒａｌｉｚｅｄ  ｊｏｉｎｔ  ｍａｐ  ｍａｔｒｉｘ）Ｗは、ブロック対角であり、そしてこれはまた、座標に依存する（つまり、Ｗ＝Ｗ（ｑ）である）が、これは重大な計算的な損失を有さない。
【００３７】
本発明に従って、分子系の方程式を解くための解法は、残差形式で表現され、直接的にこの状態微分を生成する慣例的な工程を回避する。この残差形式は、以下の工程を有する：
１）この解法の変数の差分化（ｄｉｓｃｒｅｔｉｚａｔｉｏｎ）。この差分化の特別な形式が、時間において分子モデルを前進させるために使用される特定の陰的積分法によって、決定づけられる。陰的積分法は、この残差形式に続く。陰的積分法、特に、Ｌ安定性積分器、および他の高度に安定性である積分器（例えば、陰的Ｅｕｌｅｒ法、Ｒａｄａｕ５、ＳＤＩＲＫ３、ＳＤＩＲＫ４、他の陰的Ｒｕｎｇｅ−Ｋｕｔｔａ法およびＤＡＳＳＬ）または他の陰的多段法がまた、分子モデリングについての他の利点を提供する。例えば、上記に引用された、同一日付けで出願された、「ＭＥＴＨＯＤ  ＦＯＲ  ＬＡＲＧＥ  ＴＩＭＥＳＴＥＰＳ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題する、米国特許出願番号＿＿＿＿＿を参照のこと。具体的な単純な例として、陰的Ｅｕｌｅｒ積分法を使用する場合、この差分化は、以下に従う：
【００３８】
【数１０】
2004527027.tif 000011
ここで、ｈはその時間刻みである。
【００３９】
２）この残差方程式への置換：
【００４０】
【数１１】
2004527027.tif 000012
３）ｑｎおよびｕｎについてのこの得られた非線形代数方程式
【００４１】
【数１２】
2004527027.tif 000013
の解法。
【００４２】
運動学的残差ρｑは、陰的積分器によって生成された推定
【００４３】
【数１３】
2004527027.tif 000014
を、以下で詳説されるようなこの分子モデルのジョイントを決定するためのルーチンによって計算された微分と比較する。この残差の第２行は、動力学的残差ρｕであり、これは、推定された
【００４４】
【数１４】
2004527027.tif 000015
が運動方程式を満たす度合いを決定する。
【００４５】
系の質量行列Ｍおよびいわゆる「無バイアスヒンジトルク（ｂｉａｓ−ｆｒｅｅ  ｈｉｎｇｅ  ｔｏｒｑｕｅ）」ｆは、両方とも状態依存性である。無バイアスヒンジトルクは、残差ルーチンへと運ばれる計算された
【００４６】
【数１５】
2004527027.tif 000016
ベクトルが零（０）である場合に、動力学的残差によって生成される。一般的に、このヒンジ加速度は、適用された力、ジョイントトルク、および運動誘導性効果（例えば、Ｃｏｒｉｏｌｉｓ力、および遠心力）に応答するものである。この系が静止しており、かつジョイントトルクにのみ供されている場合、無バイアス状態にあると考えられる。現実的な入力を伴う実際の系は、バイアスのある入力と等価なジョイントトルクの集合を計算することによって無バイアス状態へと換算され得る。両方の集合は、同じヒンジ加速度を生じる。
【００４７】
残差形式の好ましい実施形態は、その分子モデルのための次数（Ｎ）ねじれ角、剛体について示される。以下の節によって、基礎的な定義から分子モデルが展開され、そして、このモデルの運動を計算するためにこのモデルがどのように使用しされるかを示す。第１に、分子モデリングシミュレーションのための、全コンピュータコードアーキテクチャが、記載される。次いで、次数（Ｎ）ねじれ角、剛体多体系が、使用された記数法、基準配置、物体間のジョイントの定義、一般化座標、および一般化速度に従って、導出される。動力学状態についてのこのアプローチは、Ｔ．Ｒ．Ｋａｎｅ（Ｄｙｎａｍｉｃｓ、第３版、１９７８）によって使用されたアプローチと類似している。
【００４８】
最後に、この残差形式の結果を使用した解析的ヤコビアンの好ましい実施形態が開発された。
【００４９】
（分子動力学シミュレーションアーキテクチャ）
本発明に従うモデリングモジュールのための、ソフトウェアの汎用システムアーキテクチャ４８およびそのいくつかのプロセスを、図１に示す。大きい矩形ブロックの各々は、ソフトウエアモジュールを表し、矢印はソフトウエアモジュール間を通過する情報を表す。このソフトウエアシステムアーキテクチャは、モデラー（ｍｏｄｅｌｅｒ）モジュール５０、生化学的コンポーネントモジュール５２、物理学的モデルモジュール５４、解析モジュール５６、および可視化モジュール５８を有する。これらのモジュールのいくつかの詳細は、以下に説明される。他のモジュールが、公で利用可能である。
【００５０】
モデラーモジュール５０は、ユーザが特定の分子系を規定する物理学的パラメータを入力するインターフェースを提供する。インターフェースはグラフィックまたはデータファイル入力（またはその両方）を有し得る。生物学的コンポーネントモジュール５２は、分子系の特定の数学モデルのためのモデラー入力を変換し、モデリングされる系の、分子、力場、および溶媒のそれぞれを、数学的モデリングするための、変換サブモジュール６０、６２、および６４に、分割される。例えば、Ｔｉｎｋｅｒ（Ｊａｙ  Ｐｏｎｄｅｒ，「ＴＩＮＫＥＲ  Ｕｓｅｒ’ｓ  Ｇｕｉｄｅ」，Ｖｅｒｓｉｏｎ３．８，２０００年１０月，Ｗａｓｈｉｎｇｔｏｎ  Ｕｎｉｖｅｒｓｉｔｙ，Ｓｔ．Ｌｏｕｉｓ，ＭＯ）を含む、利用可能ないくつかのモデラーおよび生物学的コンポーネントモジュールが存在する。
【００５１】
生物学的コンポーネントモジュール５２から変換された物理パラメータを用いて、物理学的モデルモジュール５４は、分子系を数学的に定義する。多体系サブモジュール６６が、モジュール５４の中心に存在する。物理モデルモジュール５４および可視化モジュール５８と通信する解析モジュール５６は、物理学的モデルモジュール５４によって規定された分子系の計算モデルに解を提供する。解析モジュール５６は、物理学的モデルモジュール５４の異なる方程式を積分する、積分器サブモジュール６８のセットからなる。積分器サブモジュール６８は、経時的に分子系を前進させ、さらに分子系の極小エネルギー配置を決定する際に使用される静的解析を提供する。解析モジュール５６および積分器サブモジュールジュール６８は、本発明の内容のほとんどを含み、以下で詳細に説明される。
【００５２】
可視化モジュール５８は、生物学的コンポーネントモジュール５２および解析モジュール５６からの入力情報を受信し、ユーザに分子系の三次元画像表示および分子系に対して得られた解を提供する。多くの可視化モジュールが現在利用可能であり、例えば、ＶＭＤ（Ａ．Ｄａｌｋｅら、「ＶＭＤ  Ｕｓｅｒ’ｓ  Ｇｕｉｄｅ」，Ｖｅｒｓｉｏｎ１．５，２０００年６月、Ｔｈｅｏｒｅｔｉｃａｌ  Ｂｉｏｐｈｙｓｉｃｓ  Ｇｒｏｕｐ，Ｕｎｉｖｅｒｓｉｔｙ  ｏｆ  Ｉｌｌｉｎｏｉｓ，Ｕｒｂａｎａ，Ｉｌｌｉｎｏｉｓ）である。
【００５３】
記載されたソフトウェアコードは、Ｐｅｎｔｉｕｍ（登録商標）ＩＩＩマイクロプロセッサおよびＰｅｎｔｉｕｍ（登録商標）ＩＶマイクロプロセッサ（これえらは、Ｉｎｔｅｌ  Ｃｏｒｐｏｒａｔｉｏｎ（Ｓｎｔａ  Ｃｌａｒａ、Ｃａｌｉｆｏｎｉａ）によって製造される）を搭載したＰＣのような、従来のパーソナルコンピュータにて実行される。このことは、計算を実行させるためにスーパーコンピュータを使用する、分子モデリングにおける多くの現在の努力と対照をなす。糖残、さらなるスピードの改善が、より速いコンピュータ上で記載したソフトウェアを実行させることによって、得られる。
【００５４】
（分子モデルおよび多体系の説明）
サブモジュール６８における、以下で説明される積分器は、多体系（ｍｕｌｔｉｂｏｄｙ  ｓｙｓｔｅｍ；ＭＢＳ）に関する分子モデルの運動を記述する方程式のセットについて演算する。以下に詳細に説明された積分法の計算を支援するために、本発明に従って、ねじれ角、剛体モデルが、対象とする分子系を記述するために使用される。内部座標（選択された、一般化座標および一般化速度）が分子の状態を記述するために使用され得る。
【００５５】
このＭＢＳは、モデリングされるされる分子系を構成する、原子および有効な剛体性の結合の抽象化であり、かつ、シミュレーションによって扱われる問題についての重要な特徴を失うことなく、実際の物理学系、その環境下にある分子を単純化するように選択される。図１に示された一般的なシステムアーキテクチャに関して、このＭＢＳは、原子間の静電荷または他のエネルギー的相互作用を含まないだけでなく、この分子を囲む溶媒のモデルも含まない。力場は、生化学的コンポーネントのサブモジュール６２でモデリングされ、そして、溶媒は、生化学的コンポーネントのサブモジュール６４でモデリングされる。
【００５６】
図２は、対象分子のＭＢＳのツリー構造を示す。このＭＢＳの基本的な抽象化は、ヒンジ接続された剛体１７０の１つ以上の集合の抽象化である。剛体は、その物体を構成する全ての粒子が互いに対して固定された位置を有する物理的な物体の数学的抽象化である。伸縮または他の相対運動は許容されない。ヒンジ接続は、２つの剛体間の許容される相対運動を規定する数学的抽象化である。これらの剛体およびヒンジ接続の例は、以下に説明される。
【００５７】
基本物体（ｂａｓｅ  ｂｏｄｙ）１７２と呼ばれる１以上の物体は、運動力学的状態が、グラウンド１７４上の基準点に対して直接参照される点において固有の状態を有する。系のグラフは、１つ以上の「ツリー」である。ツリーの重要な特性は、任意の物体から任意の他の物体への経路が固有であり、すなわち、このグラフはループを含まないということである。ツリー内の物体は、ｎ個存在する（この出発点（ｂａｓｅ）は、標識１を有する）。ツリー内の物体は、規則的な標識が割り当てられ、このことによって、これらの標識が、基本物体から任意のリーフ物体（ｌｅａｆ  ｂｏｄｙ）１７６までの任意の経路上で決して減少しないことが意味される。リーフ物体は、単一の他の物体のみに接続される物体である。規則的な標識は、標識ｎをリーフ物体１７８の１つ（少なくとも１つでなければならない）に割り当てることによって達成され得る。この物体がグラフから除去される場合、その場合のツリーは、ｎ−１物体を有する。次いで、標識ｎ−１は、そのリーフ物体１８０のうちの１つに割り当てられ、全ての物体が標識化されるまで、このプロセスが繰り返される。これはまた、系内の任意の残っているツリーに対して行われる。
【００５８】
物体間の関係を維持することを支援するために、整数型の関数が使用され、系の各物体について内側物体（ｉｎｂｏａｒｄ  ｂｏｄｙ）を記録する。各出発点（ｂａｓｅ）対しては、内側物体が基準（ｇｒｏｕｎｄ）であり、物体ｋ（１８４）についてのペアレント物体すなわち内側物体１８２（ｉ）は、ｉ＝ｉｎｂ（ｋ）として参照される。さらに、記号Ｎは、慣性系、すなわち基準系１７４をいう。上付き文字Ｏは、座標原点（０，０，０）をいう。
【００５９】
ある点から別の点までのベクトルについての記号は２点の名前を含む。従って、ｒＰＱは、点Ｐから点Ｑまでのベクトルである。基準系内の点の速度を表すベクトルは、以下のように点および基準系の名称を含む：ＮｖＰ。以後導入される特定の記号は、２つの基準系を関連付ける。この場合では、記号は２つの系の名前を含む。従って、ｉＣｋは、系ｉにおける系ｋの配向のための方向余弦行列である。この記号はそのペアレント系における代表的物体のための方向余弦行列を指す。従って、ｉＣｋ（ｊ）は、問題にしている現実の物体ｊを示す。左および右の上付き文字は、物体のインデックスを変更しない。これはまた、他の記号についても真である。アスタリスクは、転置を示す（例えばＨ＊（ｋ））。ベクトルの上のチルダ記号は、３×３歪対称クロス乗積、すなわち
【００６０】
【数１６】
2004527027.tif 000017
を表す。
【００６１】
【数１７】
2004527027.tif 000018
は、ｉ×ｉの恒等行列であり、そして
【００６２】
【数１８】
2004527027.tif 000019
は、長さｉの零ベクトルであり、
【００６３】
【数１９】
2004527027.tif 000020
は、ｉ×ｉの零行列である。
【００６４】
（モデルの剛体）
図３は、ＭＢＳのサンプル「ツリー」の基準配置１９０を示す。１つを超えるツリーが許容される。各物体の点は、そのヒンジポイントＱとして表される。例えば、点Ｑｋ  １８６は、物体ｋ  １８４のヒンジポイントである。座標軸の固定された集合は、慣性系１９８に確立される。ＭＢＳの任意の配置は、その基準配置１９０として選択される。この配置における慣性座標軸の写像は、各物体内で物体に固定された軸の集合を確立するために使用される。基準配置では、各ヒンジポイントＱはＰ（そのペアレント物体（または拡張された物体）の点）と一致する。各物体に対して、点Ｐは、物体の内側ヒンジポイントと呼ばれる。そのため、物体ｋ  １８４についての内側ヒンジポイント  Ｐｋ  １８８は、ペアレント物体ｉ  １８２に固定された点である。各基本物体の内側ヒンジポイントは、グラウンドに固定された点Ｏ  １９２である。図２に示された拡大図は、点Ｑｋ  １８６が、本体ｋ  １８４に固定され、点Ｐｋ１８８はペアレント物体ｉ  １８２に固定されることをより明確に示す。
【００６５】
ヒンジポイントの位置は、各物体についての定ベクトルである
【００６６】
【数２０】
2004527027.tif 000021
１９４を規定し、これはまた
【００６７】
【数２１】
2004527027.tif 000022
と記載され得る。物体ｋについてのベクトルは、そのペアレント物体ｉに固定される。このベクトルは、物体ｉのヒンジポイントから物体ｋの内側ヒンジポイントの範囲にある。ベクトル
【００６８】
【数２２】
2004527027.tif 000023
１９６は、慣性系原点から第１の基本体の内側ヒンジポイント（またグラウンド内に固定された点）の範囲にあり、
【００６９】
【数２３】
2004527027.tif 000024
と記載され得る。
【００７０】
物体に対して、ｍ（ｋ）、
【００７１】
【数２４】
2004527027.tif 000025
は、そのヒンジ点Ｑｋに対する物体ｋの質量の特性を規定する。これらはそれぞれ、質量、第１質量モーメント、その物体の座標系にある、そのヒンジポイントについての物体の慣性行列である。粒子の分布を構成する剛体について、その質量特性は、予備処理モジュールによって計算される定数である。それらの計算の詳細は、ＫａｎｅによるＴ．Ｒ．Ｄｙｎａｍｉｃｓ，第３版、１９７８年１月、Ｓｔａｎｆｏｒｄ  Ｕｎｉｖｅｒｓｉｔｙ，Ｓｔａｎｆｏｒｄ，ＣＡ等の標準的な参考文献において見出され得る。
【００７２】
ヒンジポイントＱｋに対する物体ｋの空間的な慣性Ｍ（ｋ）は、対称６×６行列
【００７３】
【数２５】
2004527027.tif 000026
によって与えられる。
【００７４】
この系における各ジョイントは幾何学的データによって説明される。例えば、ピンジョイントは、そのジョイントによって接続された２つの物体に固定された軸によって特徴付けられる。ジョイントついての具体的なデータは、その型に依存する。番号ｎ、ｉｎｂ関数、系の質量特性、ベクトル
【００７５】
【数２６】
2004527027.tif 000027
、およびジョイント幾何学データ（ジョイント型を含む）は系パラメータを構成する。
【００７６】
（モデルのジョイントおよび一般化座標）
図４Ａ〜図４Ｃは、ＭＢＳの好適な実施形態のジョイントの定義（スライダージョイント１００、ピンジョイント１０２、およびボールジョイント１０４）を示す。各ジョイントは、内側ヒンジポイントＰｋ  １０８に関してヒンジポイントＱｋ  １０６の並進変位または回転変位を可能にする。これらの変位は、物体ｋについての一般化座標ｑ（ｋ）１１０によりパラメータ表示される。ちなみに、一般化座標は一般化量の一例であり、回転特性および並進特性の両方を有する量を指すことに留意されるべきである。例えば、点において作用する一般化力は、力ベクトルおよびトルクベクトルの両方からなる。スライダージョイント１００についての一般化座標ｑ（ｋ）は滑り変位ｘ  １１２である。ピンジョイント１０２についての一般化座標ｑ（ｋ）は、角変位θ  １１４である。ボールジョイント１０４についての一般化座標ｑ（ｋ）は、Ｅｕｌｅｒパラメータ（ε１，ε２，ε３，ε４）１１６である。
【００７７】
ジョイントの各々は、ピン、スライダー、またはボールジョイント；あるいはこれらのジョイントの組み合わせであり得る。フリージョイント、Ｕジョイント、シリンダージョイント、およびベアリングジョイントを含むがこれらに限定されない、多くの他のジョイント型が可能である。例えば、基本物体内側ヒンジポイントから基本物体ヒンジポイントまでのベクトルの慣性測定数、ｑ（ｋ）＝（ｘ，ｙ，ｚ）は、３つの直交スライダージョイントとして、グラウンドにおける基本物体の変位を表現する。フリージョイントは、ボールジョイントと結合された３つの直交スライダージョイントからなり、全部で６つの自由度を有する。
【００７８】
全ての物体に対する一般化座標の集合は、系の一般化座標であるベクトル
【００７９】
【数２７】
2004527027.tif 000028
を含む。
【００８０】
特定のジョイントについての一般化座標が与えられると、以下の２つの量：
ジョイント並進ベクトル
【００８１】
【数２８】
2004527027.tif 000029
およびペアレントにおける本体ｋに対する方向余弦行列
ｉＣｋ（ｋ）
が形成される。
この並進ベクトル
【００８２】
【数２９】
2004527027.tif 000030
は、ペアレント物体の座標系における物体ｋの内側ヒンジポイントＰから物体ｋのヒンジポイントＱへのベクトルを表す。これらの計算の詳細は、ジョイント型に依存し、容易に導出される。この記述の目的のために、系の一般化座標で与えられた
【００８３】
【数３０】
2004527027.tif 000031
およびｉＣｋ（ｋ）を生成する関数の呼び出しが、前提となる。
【００８４】
導入されたように、各物体についてのヒンジポイントの選択は任意である。しかし、賢明な選択によって、問題が大幅に簡略化される。例えば、ピンジョイントに対して、このヒンジポイントは、ジョイントの軸上にある点として選択されるべきである。この選択に対して、点Ｐおよび点Ｑは、ジョイント角度の全ての値に対して一致したままであるため、このジョイントの並進は０である。点Ｑがこの軸からの離れて選択される場合、点ＰおよびＱは、以下：
【００８５】
【数３１】
2004527027.tif 000032
のように互いに対して移動し、ここで、
【００８６】
【数３２】
2004527027.tif 000033
はジョイント軸の単位ベクトルであり、θはそのジョイント角度であり、
【００８７】
【数３３】
2004527027.tif 000034
は、軸上の任意の点から点Ｑへのベクトルである。
【００８８】
ピンジョイントおよびボールジョイントについて、軸上の点がヒンジポイントとして常に選択される。これらのジョイントに対して、並進ベクトル
【００８９】
【数３４】
2004527027.tif 000035
は０である。
【００９０】
スライダージョイントについて、並進ベクトル
【００９１】
【数３５】
2004527027.tif 000036
は
【００９２】
【数３６】
2004527027.tif 000037
である。
【００９３】
ピンに対する方向余弦行列は、
【００９４】
【数３７】
2004527027.tif 000038
である。スライダーに対する方向余弦行列は、
【００９５】
【数３８】
2004527027.tif 000039
である。
【００９６】
（モデルの一般化速度）
ペアレントｉにおいて測定された物体ｋのヒンジポイントの一般化速度を、一般化速度の集合ｕ（ｋ）によってパラメータ表示されたｉＶｋ（ｋ）とすると
【００９７】
【数３９】
2004527027.tif 000040
であり、
ここで、行列Ｈ（ｋ）は、このジョイントに対してジョイントマップと呼ばれる。これは６×ｎｕ（ｋ）の行列であり、ここで、ｎｕ（ｋ）は、そのジョイントに対する自由度の数（ピンまたはスライダに対して１、ボールに対して３、フリージョイントに対して６）である。Ｈ（ｋ）は、一般的には座標ｑに対する依存性を有し得る。ジョイントについての一般化速度が与えられると、ジョイントマップは、チャイルド（子）物体の系において表現されたジョイントの線形速度および角速度を生成する。これらのジョイントとして、本発明者らは、
【００９８】
【数４０】
2004527027.tif 000041
を使用する。
この全ての物体についての一般化速度の集合は、ベクトル
【００９９】
【数４１】
2004527027.tif 000042
、すなわち、この系についての一般化座標を含む。前述のように、（ｑ，ｕ）および特定のジョイント型が与えられた場合、ベクトルｉＶｋ（ｋ）を生成し得る関数の呼び出しが前提となる。以下の微分：
【０１００】
【数４２】
2004527027.tif 000043
を計算し得る関数の呼び出しがまた、前提となる。このルーチンは、以下の一般化位置座標の時間微分：
【０１０１】
【数４３】
2004527027.tif 000044
を生成し、ここで、Ｗ（ｑ）は、
【０１０２】
【数４４】
2004527027.tif 000045
およびｕに関連するブロック対角行列であり、以下のジョイント型：
【０１０３】
【数４５】
2004527027.tif 000046
およびフリーのジョイントに依存するブロックの各々は、３つのスライダージョイントおよび１つのボールジョントの組み合わせである。ボールジョイントについて、３つのｕと関連する４つの
【０１０４】
【数４６】
2004527027.tif 000047
（Ｅｕｌｅｒパラメーターの微分）が存在することに留意のこと。
【０１０５】
同様に、ｉＡｋ（ｋ）、すなわち、そのペアレント（ｐａｒｅｎｔ）の物体ｋのヒンジポイントの一般化加速度は、以下：
【０１０６】
【数４７】
2004527027.tif 000048
によって与えられる。
【０１０７】
ここで計算されるのは、この記述目的のための、分子系の、これらの一般化座標
【０１０８】
【数４８】
2004527027.tif 000049
、一般化速度
【０１０９】
【数４９】
2004527027.tif 000050
、内部座標である。この代表的な慣性座標（ｘ，ｙ，ｚ）およびこれらの慣性座標系における速度を用いて取り組むよりも、対象分子系についての計算が軽減される。
【０１１０】
（運動方程式の計算）
記載のように例示的な剛体多体系、ねじれ角モデルが与えられると、運動方程式が計算され得る。本発明に従って、このＭＢＳ分子モデルの運動は、残差形式によって、決定される。この残差形式法は、「第１」運動学的計算と呼称される計算を必要とし、「第２」運動学的計算とこの「第１」運動学的計算を区別し、ここで、第２の運動学的計算は、直接形式（これは、比較目的のためにこの説明に含まれる）によってさらに必要とされる。
【０１１１】
（分子モデルのための第１運動学的計算）
第１の運動学的計算において、この分子系の内部座標、すなわち、
【０１１２】
【数５０】
2004527027.tif 000051
およびこの系のパラメータが与えられると、分子モデルの各剛体ｋについての、以下の位置、速度および加速度の運動学的状態が計算される（ちなみに、第１運動学的計算は、残差形式法について実施された場合、ここで、
【０１１３】
【数５１】
2004527027.tif 000052
が、推定された解として、投入され、次いで、この積分法がこの解を、正確な解へと精緻化することに留意すべきである。対照的に、
【０１１４】
【数５２】
2004527027.tif 000053
が、直接形式について使用される場合には、零（０）に設定される。このことは、後の２つ方法の説明において明確に示される）。
【０１１５】
物体ｋの各々について、
【０１１６】
【数５３】
2004527027.tif 000054
を計算する。これらの計算は、反復的に実行され、基準の物体の各々から開始され、その細部へと進行する。
【０１１７】
ＮＣｋ（ｋ）、すなわち、基準（ｇｒｏｕｎｄ）におけるその物体ｋについての方向余弦行例が、以下のように定義される：
ＮＣｋ（１）＝ｉＣｋ（１）
ＮＣｋ（ｋ）＝ＮＣｋ（ｋ）ｉＣｋ（ｋ），ｋ＝２，．．．ｎ，ｉ＝ｉｎｂ（ｋ）
。
ｉＣｋ（ｋ）は、上記のジョイントルーチンから得られる。
【０１１８】
そのペアレント系において表現される場合、物体ｋのペアレントのヒンジポイントＱｉからの、物体ｋのペアレントのヒンジポイントＱｋへの位置ベクトルである、
【０１１９】
【数５４】
2004527027.tif 000055
は、以下：
【０１２０】
【数５５】
2004527027.tif 000056
のように定義される。
【０１２１】
【数５６】
2004527027.tif 000057
は、ジョイントルーチンから得られる。
【０１２２】
グローバル系として表現される場合、慣性系の原点Ｏから物体ｋのヒンジポイントＱｋへの位置ベクトルである
【０１２３】
【数５７】
2004527027.tif 000058
は、以下：
【０１２４】
【数５８】
2004527027.tif 000059
のように定義される。
【０１２５】
物体ｋについての剛体の変換演算子ｉφｋは、以下：
【０１２６】
【数５９】
2004527027.tif 000060
のように定義される。
【０１２７】
物体ｋの系で表現される場合、そのヒンジポイントでの物体ｋについて空間ベクトルＶ（ｋ）は、以下：
【０１２８】
【数６０】
2004527027.tif 000061
のように定義される。
【０１２９】
物体ｋの系で表現される場合の、そのヒンジポイントにおける、物体ｋについての空間的加速度Ａ（ｋ）は、以下：
【０１３０】
【数６１】
2004527027.tif 000062
で定義され、ここで、
【０１３１】
【数６２】
2004527027.tif 000063
である。
当然、これらの計算は、所望される場合、全て、単一の経路で計算され得る。
【０１３２】
時間刻み１増分についてのこれらの工程を完結した後に、ＭＢＳは、任意の物体の任意のポイントについて、（一般化）位置、（一般化）速度または（一般化）加速度の情報を計算するための運動論的要求に助力を提供する。これは、標準的な剛体式を使用して、その物体についてのヒンジ量に関する任意のポイントについて必要とされる情報を計算することによってなされる。
【０１３３】
（残差計算）
上記の第１運動学的計算を使用して、残差形式法のための残差計算が決定され得る。残差形式およびその分子モデリングへの適用の詳細な説明は、「ＭＥＴＨＨＯＤ  ＦＯＲ  ＲＥＳＩＤＵＡＬ  ＦＯＲＭ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題される、同一日に出願された上述の同時係属中の米国特許出願番号において見出される。この計算は、上で与えられたベクトル
【０１３４】
【数６３】
2004527027.tif 000064
の２つの部分からなる。第１の部分は、ρｑと呼称される、運動学的残差であり、そして、第２の部分は、ρｕと呼称される、動力学的残差である。この運動学的残差は、
【０１３５】
【数６４】
2004527027.tif 000065
（これは、（陰的）積分サブモジュール６６から投入される）と、各ジョイントによって計算された導関数との間の差から計算される：
【０１３６】
【数６５】
2004527027.tif 000066
。
【０１３７】
動力学的残差がまた、計算される。分子モデルの所定状態、すなわち、所定の
【０１３８】
【数６６】
2004527027.tif 000067
および系のパラメータを用いて開始して、プログラムルーチンは、ＭＢＳの「環境」をモデリングする。このようなルーチンは、コンピュータモデリング分野の当業者にとって容易に利用可能であり、作製可能であり得る。このルーチンは、積分サブモジュール６６の形態によって決定され、そして通過する、値（ｑ，ｕ）をとり、（状態依存的な）そのヒンジポイントＱｋにおける物体ｋについて適応された空間的力
【０１３９】
【数６７】
2004527027.tif 000068
、物体ｋについてのヒンジトルクσ（ｋ）を返す。その後、物体ｋについての一般化速度ｕ（ｋ）に関するこの動力学的残差ρｕ（ｋ）は、以下の工程１〜３：  １．投入された（ｐａｓｓｅｄ−ｉｎ）状態の値
【０１４０】
【数６８】
2004527027.tif 000069
を用いて、上記のように残差形式によって分子モデルのための計算を実施する
２．ｎ個の物体を有するモデルにおける各物体ｋについての空間的負荷収支
【０１４１】
【数６９】
2004527027.tif 000070
の生成：
【０１４２】
【数７０】
2004527027.tif 000071
３．ρｕ（ｋ）の計算
【０１４３】
【数７１】
2004527027.tif 000072
によって、計算される。
【０１４４】
残差形式法は、微分方程式が満たす程度を評価する。零（０）の残差は、適応された力が、慣性力とバランスをとることを示す。しかし、このことは、この系が静的平衡にあることを意味せず、それよりもむしろ、適用された力が、状態（ｑ，ｕ）にある系に対して適用された場合に、所定の
【０１４５】
【数７２】
2004527027.tif 000073
を再形成することを示す。この残差は、適用された力および慣性力とバランスをとるために必要とされる、さらなるヒンジトルクとして解釈され得る。文献によると、この方法は、逆方向動力学（ｉｎｖｅｒｓｅ  ｄｙｎａｍｉｃｓ）、または計算されたトルクの方法のいずれかとして公知である。これは、
【０１４６】
【数７３】
2004527027.tif 000074
が全て規定された場合を支配する。この点において、残差形式について必要とされる計算は、完全である。この残差ρｑおよびρｕは、積分器サブモジュール６８の中の陰的積分器によって直接的に使用される。
【０１４７】
（分子モデルのための第２運動学的計算）
直接形式法を実行するために、第１の運動学的計算が必要とされる。これらのさらなる計算は、第２運動学的計算と称する。Ｐ（ｋ）、Ｄ（ｋ）、ｉΨｋ（ｋ）、ｉＫｋ（ｋ）の値が、以下のように計算する：
１．上述したような残差形式によって分子モデルについての計算を実施（すなわち、第１運動学的計算）
２．各物体ｋの有節（ａｒｔｉｃｕｌａｔｅｄ）物体の慣性Ｐ（ｋ）の初期化。
【０１４８】
Ｐ（ｋ）＝Ｍ（ｋ）、ｋ＝１，．．．，ｎ
３．以下の対象の生成
【０１４９】
【数７４】
2004527027.tif 000075
これらの量の関数依存性は、一般化座標であるｑついてのみである。従って、この第１運動学的計算は、第２運動学的計算を見越してプログラムされる。
【０１５０】
（前進動力学計算）
最後に、
【０１５１】
【数７５】
2004527027.tif 000076
がこの分子の、以下のように、内側、次いで、外側に、及ぶことよって、計算される：
【０１５２】
【数７６】
2004527027.tif 000077
。
第１運動学的計算および第２運動学的計算、ならびに前進的動力学計算を使用して、この直接形式法が利用可能である。
【０１５３】
（運動方程式のための直接形式法）
直接形式法は、現在の状態（ｑ，ｕ）を取り、そしてその微分
【０１５４】
【数７７】
2004527027.tif 000078
を上記アルゴリズムを使用して計算し、このアルゴリズムは、時間を前進する積分法によって使用される。
【０１５５】
状態（ｑ，ｕ）が与えられると、
以下：
【０１５６】
【数７８】
2004527027.tif 000079
を計算する。
１．上記のようにジョイント特異的なルーチンを使用して
【０１５７】
【数７９】
2004527027.tif 000080
を計算する
２．上記の第１運動学的計算を、
【０１５８】
【数８０】
2004527027.tif 000081
を用いて実施する
３．上記のように、残差ρｕを生成する
４．残差ρｕ＝−ρｕを否定する
５．第２運動学的計算を実施する
６．上記前進動力学計算を使用して、
【０１５９】
【数８１】
2004527027.tif 000082
を計算する。
【０１６０】
直接形式法は、この系について作用する適用された力に応答してヒンジ加速度
【０１６１】
【数８２】
2004527027.tif 000083
を生成する。図５は、残差形式法および直接形式法の計算工程を概説する。
（陰的積分法におけるヤコビアン）
分子（例えば、タンパク質）をモデリングするＭＤ方程式は、多体系（ＭＢＳ）としてインプリメントされる。これらの方程式は、Ｎｅｗｔｏｎの法則を表し、そして、以下の微分方程式：
【０１６２】
【数８３】
2004527027.tif 000084
の集合によって、表現される。この微分方程式は、Ｎ次多体系動力学法のスイートを使用して、インプリメントされる。経時的にこの方程式を前進させるために、本発明に従って、数値積分の陰的方法、特にＬ安定性陰的積分法（例えば、陰的Ｅｕｌｅｒ法、Ｒａｄａｕ５法およびＳＤＩＲＫ３法）が使用され得る。
【０１６３】
この積分過程の重要な成分は、この微分方程式のヤコビアンの生成である。これは、以下：
【０１６４】
【数８４】
2004527027.tif 000085
である。この関数ｆは、それ自体、陽的な定式よりはむしろ、アルゴリズムによって計算さるので、このヤコビアン計算は、実質的な量の仕事を表現する。最も単純なアプローチにおいて、このヤコビアンは、この微分ルーチンを区別することによって数値的に形成され得る。これは、扱いにくい操作である。なぜならば、このヤコビアンの質は、丸め誤差と打ち切り誤差との間の諸条件の兼ね合いであるからである。結果における、代表的に半分の実際的な精度が、異なるスキームの中で十分な摂動サイズを選択するこのよって保持される。実際には、しかし、これは、実行することは困難である。
【０１６５】
しかし、この支配的な方程式の構造を活用して、ヤコビアンの計算を改善し得る。この例示的な多体系動力学法は、これを例示する。関連するアルゴリズムは数値的な方法がこの式を実行するために使用されたとしても、正確な微分を計算する。得られる微分は、丸め誤差および検討中の多体系の条件付けに依存する量による誤差の範囲にある。しかし、何一つ近似は、この方程式のレベルには含まれていない。
【０１６６】
一般に、Ｇ、つまり、この陰的積分器のＮｅｗｔｏｎループにおいて使用される反復行列は、Ｇ＝Ｅ−αＪの形式を有し、ここで、Ｅは、単位行列であり、αは、時間刻みのスカラー関数である。陰的積分器の記載について、同一日付けで出願された、「ＭＥＴＨＯＤ  ＦＯＲ  ＬＡＲＧＥ  ＴＩＭＥＳＴＥＰＳ  ＩＮ  ＭＯＬＥＣＵＬＡＲ  ＭＯＤＥＬＩＮＧ」と題する、上述の参照した米国特許出願番号＿＿＿＿＿を参照のこと。刻み幅の変化は、Ｇを因数分解（ｒｅｆａｃｔｏｒ）することを要求するが、しかし、Ｊを再形成（ｒｅｆｏｒｍ）することは必要としない。Ｊを再度形成することは、ヤコビアンが新状態で必要とされる場合のみ、必要とされる。Ｇは、Ｎｅｗｔｏｎループ内での線形性の副次的問題において使用される。以下のように解かれる：
【０１６７】
【数８５】
2004527027.tif 000086
ここで、ｒ（ｙ）は、特定の陰的積分法についての残差関数である。
【０１６８】
後ほど示されるように、Ｊは特別な構造を有しており、この構造は、Ｇに受け継がれている。これは、Ｇを用いた方程式の解法が、この構造が活用された場合に、少ない損失で実施され得ることを意味する。
【０１６９】
ヤコビアンの質は、積分器における打ち切りから得られる非線形性方程式を解く能力に影響を与える。Ｎｅｗｔｏｎループが解けないことによって、追跡工程の撤回、および積分時間刻み工程の減少を要求する。この時間刻みは、Ｎｅｗｔｏｎループにおける失敗よりも、むしろ精度によって制御されるべきある。
【０１７０】
（解析的ヤコビアンの計算）
ヤコビアンＪは、運動方程式の線形化を表す行列である。通常は、動力学的系を支配する方程式は、平衡状態、あるいは、定常運動状態の周囲で線形化される。この場合、方程式が、任意の状態の周りで線形化され、そして、全ての可能性のある寄与する項が、展開されるべきである。Ｊを以下のその分割に関してＪを記述することは、慣例的である；
【０１７１】
【数８６】
2004527027.tif 000087
。
【０１７２】
（ＪｑｑおよびＪｑｕの構造）
この
【０１７３】
【数８７】
2004527027.tif 000088
の方程式は、
【０１７４】
【数８８】
2004527027.tif 000089
であり、ここで、行列Ｗは、ブロック対角の構造を有している。各ブロックは、そのジョイント型に依存する。ピンジョイントおよびスライダージョイントは、１×１恒等ブロックを生じる。ボールジョイントは、Ｅｕｌｅｒパラメータおよび物体の角速度測定数（一般化速度）に関してＥｕｌｅｒパラメータ微分を表現する４×３ブロックを生じる。
【０１７５】
上記の
【０１７６】
【数８９】
2004527027.tif 000090
の方程式から、ヤコビ行列の２つの分割が見出される：
【０１７７】
【数９０】
2004527027.tif 000091
。
【０１７８】
これらの方程式は、象徴的な目的のみであると解釈される。実際において、行列Ｗを陽的に生成する必要はない。非ゼロのブロック対角要素は、運動学的残差についての前述の節で議論されたように満たされる。
【０１７９】
（ＪｕｑおよびＪｕｕの構造）
微分
【０１８０】
【数９１】
2004527027.tif 000092
は、さらに複雑である。なぜならば、
【０１８１】
【数９２】
2004527027.tif 000093
であるので、
【０１８２】
【数９３】
2004527027.tif 000094
が計算されなければならない（ρｕが以前に展開された動力学的ルーチンの残差であることに留意のこと）。ここで、数値解析の分野からの重要な結果は、この逆行列の微分を回避するのに使用される。
【０１８３】
Ａ（ｘ）ｙ（ｘ）＝ｂ（ｘ）であるとすると、ｙ（ｘ）＝Ａ−１（ｘ）ｂ（ｘ）あると記載し得る。ｙ（ｘ０）＝ｚが既知であり、かつ、
【０１８４】
【数９４】
2004527027.tif 000095
の値が得られなければならない場合、本発明者らは、
【０１８５】
【数９５】
2004527027.tif 000096
を有し、ここで、ｚ＝Ａ−１ｂは、この方程式の右辺を形成するときに固定される。上述の多体系ルーチンに適用された場合、
【０１８６】
【数９６】
2004527027.tif 000097
であり、ここで、
【０１８７】
【数９７】
2004527027.tif 000098
である。この結果は、Ｍ−１（これは、超行列（ｈｙｐｅｒｍａｔｒｉｘ）である）の微分を計算することを回避する。この逆行列は、「引き出される」。上記の方程式において、「ｘ」は、一般化座標ｑか、または一般化速度ｕであり、計算されるヤコビアン分割に依存する。
【０１８８】
要約すると、ブロックＪｕｑおよびＪｕｕを計算するために、３つの工程が続く：
１．所定の（ｑ，ｕ）は、直接法を使用して
【０１８９】
【数９８】
2004527027.tif 000099
を計算する。これは、単純に
【０１９０】
【数９９】
2004527027.tif 000100
を現在の状態から計算することをいい、そして、これを変数ｚとして保存する。
【０１９１】
２．動力学的残差ルーチンの解析的ヤコビアンを計算する。この工程において、行列
【０１９２】
【数１００】
2004527027.tif 000101
が形成される。この量は、明らかに工程１で計算されたベクトルｚに依存する。この工程における残差ρｕの数値的値は、発明者らは、この運動方程式の共通解のまわりでヤコビアンを計算しているので、各要素に対して零（０）であることを留意のこと。この残差ヤコビアンの分割ＪｕｑおよびＪｕｕは、上記「ｘ」について「ｑ」および「ｕ」で置換することによって得られる。
３．工程２の結果を、質量行列を用いて後退代入法することによって、所望のブロックを得る。この後退代入法操作は、残差ベクトルを、
【０１９３】
【数１０１】
2004527027.tif 000102
ベクトルへと処理することによって、直接形式ルーチンにおいて解決される。第２運動学的工程は、１度のみ実施される必要がある。なぜならば、この後退代入法は、この状態の名目上の値で実施されるからである。実際に、この第２運動学的ルーチンは、ｚを計算している間に、工程１で呼び出され、そしてこの変数はキャッシュに入れられる。
【０１９４】
言葉で表すと、本発明の微分ルーチンのヤコビアンは、本発明者らの残差ルーチンのヤコビアンの後退代入法によって実施され得る。この残差ヤコビアンは、この微分ヤコビアンよりもずっと容易に計算される。上記の工程２および工程３は、以下の節で導出される。
【０１９５】
（残差ヤコビアン）
残差ヤコビアンの計算は、力学についての残差形式に密接に関連し、ここで以下のように要約する：
１．位置および速度に依存する運動データを計算する、外側への（ｏｕｔｂｏａｒｄ）パスを実行する。
２．原子間力を発生させ、そしてこれらを物体に作用する空間的負荷に強化する、力ルーチンを呼び出す。
３．加速度レベル値を（入力した
【０１９６】
【数１０２】
2004527027.tif 000103
を使用して）計算し、そして慣性力を工程２からの空間的負荷と組み合わせる、別の運動パスを実行する。
４．各結合における残差を発生させる内向パスを実行する。このパスは、問題の結合の外側物体の「ネスト」について、（空間的）慣性力と適用された力との合力を再帰的に計算する。残差は、結合マップデータによって与えられる、結合の自由度の合力の投影である。
【０１９７】
高レベルにおいて、残差の計算は、２種類の力（「運動力」および外力）に依存すると考えられ得る。運動力は、多体系によって直接計算される。外力は、種々の原子間力（例えば、静電力および溶媒）を計算する力モデリングルーチンから、多体系に対して利用可能である。ヤコビアンを計算する場合に、類似の手順に従う。多体系は、運動力のヤコビアンを構築し、そしてこれを外力のヤコビアンと組み合わせる。
【０１９８】
（力場の作用への分配）
分子に作用し得る多くの力が存在し、そしてこれらの力は、種々の固有座標系（これらは、この特定の力の「効果」に対して最も好都合である）において計算され得る。例えば、静電力の項は、多極方法および球座標を使用して計算され得、共有結合の項は、ねじれ角および結合角に関して計算され得、そして溶媒力の項は、グローバルデカルト座標において計算され得る。残差の計算の間に、これらの力は、それらの固有座標系からＭＢＳ座標系へと変換される。
【０１９９】
同じ交換が、ヤコビアンを計算するために起こる。ネイティブのヤコビアンは、その固有座標において、ＭＢＳ座標にされる。このことは、固有座標とＭＢＳ一般化座標との間での変換のために、連鎖法則の使用を必要とする。各効果が、その関数値およびヤコビアンを同時に計算することが重要である。なぜなら、同じ項の多くが、それぞれ計算される必要があるからである。各効果は、空間負荷Ｔｅｆｆｅｃｔ（ｋ）のセットに変換され、ここで、ｋは、系における一般的な物体の指数である。これらの効果の合計は、記号Ｔ（ｋ）を与えられる。
【０２００】
（ヤコビアンを残差ヤコビアンにする効果）
高レベルにおいて、残差ルーチンは、以前に、以下の等式：
ρｕ  ＝  Ｈφ（ＭＡ−Ｔ）
から実行された。この実行は、次数（Ｎ）の方法（これは、上記等式からすぐに明らかである）を使用する。この等式において、Ｔは、多体系の旋回に作用する空間的負荷のベクトルであり、ここで、各要素は、空間的負荷である（１つの力および１つのトルクからなる６ベクトル）。これは実際に、慣性負荷または純粋なヒンジ負荷以外の全ての効果を表す。括弧内の項は、各物体に対する負荷のバランスを表す。第一項は、慣性力であり、そして次の項は、空間的負荷である。Ｍ（ｋ）Ａ（ｋ）は、代表的な物体に対する空間的慣性力である。これは、物体の質量特性および物体の旋回の空間的な加速度から構築される。この空間的な加速度は、残差ルーチンが先行する力学ルーチンによって実行される前に、計算される。演算子Ｈφは、次数（Ｎ）の内向パスを実行するルーチンにおいて、実行される。
【０２０１】
上記等式によって実行される計算の詳細についての全てを知ることなしにさえ、
【０２０２】
【数１０３】
2004527027.tif 000104
（空間的負荷ヤコビアン（効果ヤコビアン））の、
【０２０３】
【数１０４】
2004527027.tif 000105
（残差ヤコビアン）に対する寄与が、すぐに推定され得る：
【０２０４】
【数１０５】
2004527027.tif 000106
ここで、．．．は、項が効果ヤコビアンを含まないことを表す。再度、「ｘ」についてのｑまたはｕは、ヤコビアンのどの分配が計算されるかに依存して、置換される。
【０２０５】
残差ヤコビアンにおける効果ヤコビアンの役割は、マルチ物体等式における効果の役割と同じである。このことは、Ｔが、残差ρｕに対して、
【０２０６】
【数１０６】
2004527027.tif 000107
の列が
【０２０７】
【数１０７】
2004527027.tif 000108
の列に寄与する様式と同じ様式で寄与することを意味する。両方が、同じ演算子Ｈφによって処理される。これは、非常に重要な点である。なぜなら、このことは、ヤコビアン計算のこの部分に対して新たな方法が必要とされないことを意味するからである。（異なる方法は、
【０２０８】
【数１０８】
2004527027.tif 000109
を得ることを必要とする）。
【０２０９】
従って、効果ヤコビアンを考慮して、その寄与は、この効果ヤコビアンに対してそのもとの残差ルーチンを実行すること、この効果ヤコビアンをマルチ列セットの空間的負荷ベクトルとして処理することによって、残差ヤコビアンに組み立てられる。このことは、この等式の直線性の線形性である結果である。
【０２１０】
残差ヤコビアンの列は、残差ベクトルが順方向力学ルーチンにおいて果たすものと同じ役割を、微分ヤコビアンルーチンにおいて果たす。力学ルーチンは、このルーチンが受信するデータベクトルに対して後退代入法（ｂａｃｋ−ｓｏｌｖｅ）を実行し、そしてどのデータであるかを知る必要がなく、そのデータに対してどの演算が実行されているかのみを知る必要がある。これは、全てのルーチンに当てはまる。
【０２１１】
．．．の項を追加することによって、連鎖法則を使用し、等式全体を以下のように示す：
【０２１２】
【数１０９】
2004527027.tif 000110
このレベルにおいて、ヤコビアンに対して４つの寄与（慣性力、空間的な力、ならびに演算子φおよびＨにおける変化に起因する寄与）が存在する。量ｚおよびｙは、それぞれ（ＭＡ−Ｔ）およびφｚを表し、これらは、最後の２つの項を詳説する間、一定に維持される。これらの項の数値は、残差計算から既に利用可能である。上記等式についての別の観察は、演算子φがｑのみに依存し、そしてｕには依存しないことである。従って、この項は、分配Ｊｕｕが計算されている間、一定のままである。同様に、空間負荷は、その連続した効果に分離しえ、これらのいくつかは、ｕに依存しない。一般に、このことは、マルチ物体等式の知識が、計算を最適化するために利用され得る。
【０２１３】
ここまでで、一旦項が計算されたら
【０２１４】
【数１１０】
2004527027.tif 000111
をどう取り扱うかが記載されてきたが、この項をどのように形成するかの説明はなされていない。これらの詳細を、以下の節に示す。
【０２１５】
（効果ヤコビアンの計算および残差ヤコビアンとの組み合わせ）
ここまでは、ヤコビアン計算の高レベルの説明がなされた。これらの計算は、それを行うためのアルゴリズムの特色が非常に強い。先に記載された順方向力学ルーチンについてもまた存在したような、作業のための非常に異なる段階が存在する。そこで、原子間力の計算は、明らかに、ボトルネックの工程である。しかし、力を取り扱うためのマルチ物体の等式におけるオーバーヘッドは、かなり小さい。この場合、呼び出しが、力ルーチンに対してなされ、そしてこのルーチン内において起こることは、無視された。ヤコビアンのこととなると、この局面はさほど正確ではない。
【０２１６】
呼び出しが、効果ヤコビアンを得るために依然としてなされるが、効果ヤコビアンが残差ヤコビアンに組み立てられ得る前に必要とされる、多数の処理が存在する。力場を扱ってヤコビアンを作成する詳細は、次の節に網羅されており、そして静電力を組み込む例が開発される。他の全ての付加は、類似の発達に従う。
【０２１７】
（例の効果としての静電力）
静電力の基本的な前提は、２つの荷電粒子間の力が、以下：
【０２１８】
【数１１１】
2004527027.tif 000112
であることである。これは、古典的な逆二乗則である。Ｆｉｊ（粒子ｊに起因して粒子ｉに作用する力）は、これらの粒子の電荷に依存し、反対に荷電した粒子については引力であり、同じ荷電については反発力である。記号
【０２１９】
【数１１２】
2004527027.tif 000113
は、粒子ｉから粒子ｊへの方向の単位ベクトルである；ｒｉｊは、これらの粒子間の距離である；ｋは、力の強度に関連する、単位に依存する定数である。もちろん、粒子ｊに作用する力は、粒子ｉに作用力に等しく、そしてその逆である。従って、上で与えたクーロンの法則によって相互作用する粒子の収集を考慮すると、各粒子に作用する正味の力は、対ごとの力を合計することによって計算される。この例に対して、これらの力は、グローバルデカルト座標において計算される。
【０２２０】
原子間力が使用できる状態で、マルチ物体力が発生し得る。各物体の粒子に作用する力の系は、各物体の旋回において作用する空間的負荷によって置換される。原子間力は、物体を固定した基礎において最初に表現され、次いで、その力が結合する特定の原子のステーション座標を使用して、その旋回にシフトされる。
【０２２１】
Ｆｉｊは、ベクトルである。ｄｒｉｊに対するＦｉｊ（粒子の相対位置の小さな変化）の微分は、
【０２２２】
【数１１３】
2004527027.tif 000114
（テンソル）である。これは：
【０２２３】
【数１１４】
2004527027.tif 000115
となるようなものである。クーロン力に対して、テンソルは：
【０２２４】
【数１１５】
2004527027.tif 000116
である。
【０２２５】
（いくつかの観測）
１．力ヤコビアンは、大きさｎａｔｏｍｓ×ｎａｔｏｍｓの大きさの行列である。各要素は、３×３のテンソルである。（ｉ，ｊ）ブロックは、原子ｊの位置の小さな変化に対する原子ｉへの力の微分を与える。一般に、全ての力のモデルが、分析処理のための固有ヤコビアン方法を支持するために必要とされる。
２．力ヤコビアンについての格納要件は、急速に非実用的になる。このことは、原子間の対ごとに作用する全ての力のヤコビアンが、物体間の対ごとに作用するように、減少するかまたは「縮小」される場合の、界面の「縮小」の概念を導く。
３．クーロン力は対ごとの相互作用であるので、各力が、ヤコビアン全体における２つのブロックに寄与する。従って、各力は対照的な費用で処理され、そしてヤコビアン全体は、原子の数の二乗（すなわち、次数（Ｎ２））に比例する費用で計算される。このことは、力自体の計算費用と同じである。これは、解析的ヤコビアンを計算するために非常に良い結果である。数値のヤコビアンは、その状態の要素が摂動するごとに、新たな力の計算を必要とする。これは、数値のヤコビアンの費用に、三乗の成長（すなわち、次数（Ｎ３））を導く。従って、解析的ヤコビアンは、数値のヤコビアンより計算が安価であり、そしてより正確である。
４．ヤコビアンの計算は、力を計算するルーチンと同じルーチンにおいて、好都合になされる（同時計算）。しかし、これは代表的に、力の計算よりはるかに少ない頻度でなされる必要がある。従って、フラッグが使用されて、必要な場合にのみヤコビアンの計算を誘発する。
【０２２６】
（変位勾配へのカップリング）
（固有の）力ヤコビアンが得られたら、そのヤコビアンをさらに処理することが必要である。このことは、多体系が相対座標を使用して形成されるという事実に起因する。連鎖法則が各原子間力Ｆ（ｋ，ｉ）に適用されて、「変位勾配の計算」と称される。これは、物体ｋに属する原子ｉに対するグローバルデカルト力を表す。
【０２２７】
【数１１６】
2004527027.tif 000117
ここで、ｒ（ｐ，ｓ）は、物体「ｐ」上の原子「ｓ」の位置である。
【０２２８】
この合計の第一項は、計算されたばかりの力ヤコビアンの要素を選択する。量
【０２２９】
【数１１７】
2004527027.tif 000118
は、変位勾配の要素である。代表的な項は、一般化座標における小さな変化に起因する、原子の位置の変化を与える。この項は、厳密には、力の計算とは無関係の運動量であることに注目のこと。従って、力ヤコビアンは一回計算され得、次いで、多体系における各座標に対する連鎖法則によって、連続的に再処理され得る。この工程は、行列ベクトルの多重を必要とする。なぜなら、
【０２３０】
【数１１８】
2004527027.tif 000119
は、ｎａｔｏｍｓの実体（各３ベクトル）を有する列ベクトルであり、そしてヤコビアンは、ｎａｔｏｍｓ×ｎａｔｏｍｓの正方行列であり、ここで、各要素は３×３のテンソルであるからである。
【０２３１】
この計算を改善することが可能である。なぜなら、変位ヤコビアンにおける実体の多くは、０であることが既知であるからである。このことは、特定のヒンジを漸増させることによって、その系の全ての原子が変位するのではなく、変位したヒンジの外側でありかつ問題のヒンジから外れていない原子のみが変位するという事実に起因する。例えば、基本の物体を回転させることによって、その系の全ての原子の変化が誘導される。しかし、任意の末端物体においてねじれ角を摂動させることにより、その末端物体に属する原子のみの変化が誘導される。従って、大まかに言えば、およそ半分の仕事が、計算を最適化させることによって節約され得る。この減少は、厳密に知識に基づくアプローチから生じる。
【０２３２】
（界面の縮小）
Ｔ（ｋ）（物体ｋに対する空間的負荷）を形成するプロセスにおいて、この負荷は、物体ｋに属する原子に対して作用する原子間力から生じる。各力は、グローバル座標から局所座標に変換され、次いで物体旋回にシフトする。この手順の簡潔な言及は、以下である：
【０２３３】
【数１１９】
2004527027.tif 000120
演算子φ（ｋ，ｉ）は：
【０２３４】
【数１２０】
2004527027.tif 000121
である。ここで、ρ（ｋ，ｉ）は、物体ｋ上の原子ｉの固定されたステーション座標である。新たな値Ｔ（ｋ，ｉ）が現れることに注目のこと。これは、原子間力が、原子ｉの原子位置における空間的負荷に変わっただけのものである。
【０２３５】
【数１２１】
2004527027.tif 000122
原子の空間的負荷の第一の要素は、０である。なぜなら、個々の原子に対する力場によって、トルクが付与されないからである。
【０２３６】
ここで、Ｔ（ｋ）は、原子間力を物体の空間的負荷に関連付ける。従って、この等式の微分は、差示的原子間力と差示的空間的負荷とを関連付ける：
【０２３７】
【数１２２】
2004527027.tif 000123
。この等式の第二項Ｔ２（ｋ）は、この節の最後に議論される。なぜなら、この項は、空間的負荷を含むが、付加微分を含まないからである。このことは、この項が一般に、空間的負荷がどのように計算されるかを心配せずに処理され得ることを意味する。
【０２３８】
Ｔ（ｋ）の定義をＴ１（ｋ）に置換すると：
【０２３９】
【数１２３】
2004527027.tif 000124
であり、ここで、記号
【０２４０】
【数１２４】
2004527027.tif 000125
は、現在、行が減少された力ヤコビアンの要素である。このヤコビアンは、差示的（物体）空間的負荷を、差示的原子変位に直接関連付ける。再度、ｒ（ｐ，ｓ）は、物体ｐ上の原子ｓのグローバルな位置を表す。項
【０２４１】
【数１２５】
2004527027.tif 000126
は、各原子間力ヤコビアン要素を、原子のφ（ｋ，ｉ）行列によって重み付けされた、減少したヤコビアンにおける行き先要素に合計することによって、形成される。行が減少したヤコビアンの各要素は、６×３行列である。従って、力ヤコビアンの行は、縮小されている。この縮小は、以下の記録において明らかである：分子は物体指数のみを有し、一方で分母は物体指数と原子指数との両方を有する。１つの物体あたりの原子の数に依存して、行の減少は、差示的な空間的負荷が形成されなければならない場合に、格納時間と実行時間との両方を節約する。
【０２４２】
行を減少させる手順は、残差ヤコビアンを計算する前に１回のみなされる必要があることに、注目のこと。残差を実行するオーバーヘッドは、形成されなければならない、より小さな行列ベクトル積の減少した費用によって、オフセットを超える。
【０２４３】
【数１２６】
2004527027.tif 000127
を形成する際に、原子間力ヤコビアンの要素を節約する必要がないことに注目のこと。すなわち、各要素
【０２４４】
【数１２７】
2004527027.tif 000128
は、その
【０２４５】
【数１２８】
2004527027.tif 000129
に対する寄与が計算される間に利用可能であることのみが必要である。従って、大きな力ヤコビアンの１より多くの要素は、同時には必要とされない。
【０２４６】
代表的な原子ｒ（ｐ，ｓ）のグローバル座標は、ｒ（ｐ）（物体ｐの旋回のグローバル座標）、およびρ（ｐ，ｓ）（原子のステーション座標）の観点で計算される：
ｒ（ｑ，ｓ）  ＝  ｒ（ｐ）＋ＮＣｋ（ｐ）ρ（ｐ，ｓ）
微分することによって、本発明者らは、以下：
【０２４７】
【数１２９】
2004527027.tif 000130
を見出した（いくらかは、有限回転の運動から生じる）。
【０２４８】
この等式を、さらなる等式λ（ｐ，ｓ）＝λ（ｐ）で拡大し、そしてｗ（空間的微分）
【０２４９】
【数１３０】
2004527027.tif 000131
を定義すると、その結果は、以下である：
【０２５０】
【数１３１】
2004527027.tif 000132
ベクトルλは、各物体に対する配向の変化の割合を発生させる場合に、中断され得る。これは、空間における各点においてポテンシャルが変化し得るという意味で、場の量である。純粋な回転（変形なし）を受ける剛体については、この回転によって影響を受ける各物体に対して、これは一定である。λを計算するための次数（Ｎ）アルゴリズムは、厳密には、後の節に記載される。
【０２５１】
上記等式をＴ１（ｋ）についての等式に適用することによって、最後の換算公式：
【０２５２】
【数１３２】
2004527027.tif 000133
が得られる。
換算ヤコビアン
【０２５３】
【数１３３】
2004527027.tif 000134
の各要素は、ここで６×６行列である。還元ヤコビアンの要素は、物体の旋回における空間的負荷を、その系における別の旋回において起こる空間的微分（旋回の変位勾配へのカップリングの後）に関連付ける。
【０２５４】
行の減少は、各物体における全ての原子間力を強化し、その系における全ての原子の変位微分にカップリングされた、空間的負荷の微分を残した。列の減少は、全ての原子変位微分を強化し、空間的旋回微分にカップリングされた、空間的負荷微分を残した。従って、力ヤコビアンは、「ネイティブ」な形態に変わる。還元ヤコビアンを用いる作業は、空間的負荷微分の計算を、１つの物体あたりの原子の数のおよそ２乗で加速させる。このスピードアップは、１００倍以上に容易に近付き得る。
【０２５５】
この節は、例として静電力を使用して、原子レベルでどのように力ヤコビアンが構築されるかを示す。このヤコビアンは、差示的原子間力を差示的原子変位に関連付ける。界面収縮の概念を導入して、差示的空間負荷は、行を減少させた力ヤコビアンを介して、差示的原子変位に関連付けられることが示される。計算に対する別の改善は、最後に、完全に収縮したヤコビアンを生じ、これは、差示的食うか敵変位に対する差して器空間的負荷が繰るかも。
【０２５６】
ここで、力ヤコビアン
【０２５７】
【数１３４】
2004527027.tif 000135
が、上記の空間的変位勾配にカップリングされなければならない。行列ベクトル多重度は、この工程を実行し、そしてその系における物体の数（原子の数ではない）と規模を合わせる。
第２項Ｔ２（ｋ）は、以下：
【０２５８】
【数１３５】
2004527027.tif 000136
によって与えられ、そしてもとの空間的負荷Ｔ（ｋ）および演算子φ（ｋ，ｉ）を含む：
【０２５９】
【数１３６】
2004527027.tif 000137
は、基部物体から外向きに、繰り返し計算され、そして
【０２６０】
【数１３７】
2004527027.tif 000138
ここで、ｄｑ（ｋ）は、次の節で定義され、そして
【０２６１】
【数１３８】
2004527027.tif 000139
（物体間方向余弦行列の偏微分）は、物体を接続するジョイントの型の関数である：
【０２６２】
【数１３９】
2004527027.tif 000140
次の節において、力ヤコビアンを、慣性力ヤコビアンと組み合わせて、残差ルーチンのヤコビアンを最終的に形成する。
【０２６３】
（残差ヤコビアン）
先行する節は、力ヤコビアン
【０２６４】
【数１４０】
2004527027.tif 000141
を形成し、これは、空間的力の微分を形成するために、空間的変位勾配をカップリングしなければならない。この節は、空間的変位勾配の形成および残差ルーチンのヤコビアンの形成を記載する。
【０２６５】
ヤコビアン全体をカップリングするための残差アルゴリズムが、記載される。ヤコビアンアルゴリズムは、正確には、ヤコビアンを計算するように設定されない。自動差示ルーチンにおいて代表的であるように、これは、行列ベクトル積Ｊｕｑｄｑ＋Ｊｕｕｄｕを、ベクトルｄｑおよびｄｕについての任意の入力値について計算する。実際に、ヤコビアンを計算するために、「ヤコビアンルーチン」が効果的に、一連のブールベクトル（１つの入力が１に設定され、そして他の全ての入力場０に設定されたベクトル）を用いて繰り返し呼び出される。各呼び出しが、ヤコビアンの対応する列を発生させる。これらの工程のうちのいくつかは、残差形式方法または直接形式方法（順方向力学計算）について既に計算されたかまたは計算されているが、本明細書中では明瞭にするために再形成されていることに、注目のこと。
１．所定の（ｑ，ｕ）が、直接形式方法を使用して、
【０２６６】
【数１４１】
2004527027.tif 000142
を計算する。また、
【０２６７】
【数１４２】
2004527027.tif 000143
を設定し、そしてＡ（ｋ）、次いでρｕを再計算し、これは、
【０２６８】
【数１４３】
2004527027.tif 000144
を再計算する。
２．収縮を実行して、完全に行および列が減少した力ヤコビアン
【０２６９】
【数１４４】
2004527027.tif 000145
を、「界面構築」の節に記載されるように計算する：
【０２７０】
【数１４５】
2004527027.tif 000146
。
以下の工程３〜１０は、Ｊｕｑの列を満たすために使用される：
３．
【０２７１】
【数１４６】
2004527027.tif 000147
の項の解析的ヤコビアン分配
【０２７２】
【数１４７】
2004527027.tif 000148
を、第１運動学的計算のために必要なルーチンと類似の結合ルーチンを使用して計算する。
４．ｑ（位置量の微分）および空間的勾配についての項を計算する：
先に記載した方法を使用して、特定の結合特異的場を満たした。これらの量は、ｉＣｋ（ｋ）（物体内方向余弦行列）、
【０２７３】
【数１４７Ａ】
2004527027.tif 000149
（各物体についてのスパンベクトル）、およびＨ（ｋ）（各物体の内側への（ｉｎｂｏａｒｄ）結合についてのジョイントマップ）からなった。これらの量の各々の微分についていうために、添え字ｄがその記号の名称に付加されて、この参照を一般的にする。従って、ｉｄＣｋ（ｋ）とは、方向余弦行列ｉＣｋ（ｋ）の微分を意味する。
【０２７４】
各物体間方向余弦行列（および全ての結合特異的）量は、個々の結合の一般化座標のみに依存する。従って、ｉｄＣｋ（ｋ）は、この微分が物体ｋについての座標のいずれかに対して微分がとられる場合にのみ、０ではない。試験中の反復を適切に「シード」するために、ベクトルｄｑが、このルーチンに入力される。ヤコビアン計算のために、本発明者らは、１つの入力を１に、そして他の全ての入力を０に設定した。次いで、必要とされた予備的な量を、代表的なループによって作製した：
【０２７５】
【数１４８】
2004527027.tif 000150
方向余弦行列の偏微分を、分析的に生成し、そして上記「界面収縮」の節で表すように表示する。これらの偏微分は、計算されているヤコビアンの特定の列に依存しない。ｄｑの特定の入力を１に設定し、そして残り全てを０に設定することは、シード量の正しいサブセットの消滅の効果を有する。
【０２７６】
【数１４９】
2004527027.tif 000151
（物体間スパンベクトルの偏微分）は、
【０２７７】
【数１５０】
2004527027.tif 000152
によって与えられる。
λ（ｋ）は、本明細書中で、親の物体に接続する、物体のスライド軸をいう。
【０２７８】
【数１５１】
2004527027.tif 000153
（結合マップの偏微分）は、
【０２７９】
【数１５２】
2004527027.tif 000154
である。
上記偏微分の定義を用いて、反復は、以下のループでシードされる：
【０２８０】
【数１５３】
2004527027.tif 000155
これらのループの実行の後に、全ての物体は、
【０２８１】
【数１５４】
2004527027.tif 000156
を有し、これらの物体間微分量は利用可能である。
【０２８２】
空間的変位勾配計算において必要とされる１つの新たな量もまた、計算される。これは、界面収縮の節からのλ（ｋ）であり、摂動した結合の外側の、各物体についての配向の変化の割合を発生させる、回転軸である。ここで、この変数は、記号ｄθ（ｋ）を与えられ、各物体についての差示的な回転軸は、以下である：
【０２８３】
【数１５５】
2004527027.tif 000157
一連のオイラーパラメータに対する任意の摂動は、純粋な回転を生じないので、ヤコビアンの対応する列を計算する場合に、列の収縮は使用され得ない。行が減少した力ヤコビアンは、依然として使用される（そして使用されなければならない）。
【０２８４】
反復をシードした後に、
【０２８５】
【数１５６】
2004527027.tif 000158
が計算される：
【０２８６】
【数１５７】
2004527027.tif 000159
５．ｑ（速度の微分）を計算する：
結合角の変化に起因する、結合速度の変化の割合を計算するループは、以下のプロセスで開始する：
【０２８７】
【数１５８】
2004527027.tif 000160
この量は、結合角の変化に起因する、結合速度の変化の割合である。明らかに、マップが座標依存性を含む結合についてのみ、この量は０ではない。自由結合について、一般化された速度は、結合の配向に依存する相対線形速度を生じる。
【０２８８】
ｉｄｖｋ（ｋ）、ｄＶ（ｋ）を計算した後に、各物体の空間的速度の微分が計算される。これは、以下のループによってなされる：
【０２８９】
【数１５９】
2004527027.tif 000161
６．力ヤコビアンを空間的変位勾配にカップリングして、Ｔ１（ｋ）を計算する
【０２９０】
【数１６０】
2004527027.tif 000162
７．力ヤコビアンＴ２（ｋ）の第二項を計算して、Ｔ１（ｋ）に添える：
【０２９１】
【数１６１】
2004527027.tif 000163
８．ｑ（加速度関連項の微分）を計算する：
再度、このプロセスは、
【０２９２】
【数１６２】
2004527027.tif 000164
を計算するループで開始する：
【０２９３】
【数１６３】
2004527027.tif 000165
これは、座標の変化に起因する、結合加速度の変化である。次いで、ｄＡ（ｋ）（各物体の空間的加速度の微分）が、以下のように計算される：
【０２９４】
【数１６４】
2004527027.tif 000166
ここで、→は、６ベクトルが２つの３ベクトルに分解されることを表し、
Ｋ＝２〜ｎｂｏｄについて、
【０２９５】
【数１６５】
2004527027.tif 000167
ここで
【０２９６】
【数１６６】
2004527027.tif 000168
によって導入される記号は、ｄＡ（ｋ）の計算の後には必要とされない一時的な変数を意味する
空間的加速度の微分を計算した後に、
【０２９７】
【数１６７】
2004527027.tif 000169
（空間的慣性力の微分）の計算が実行される：
【０２９８】
【数１６８】
2004527027.tif 000170
９．ｄρ（ｋ）（物体ｋに対する結合残差微分）を計算する：
【０２９９】
【数１６９】
2004527027.tif 000171
このルーチンを実行した後に、ｄρ（ｋ）に格納された値は、残差ヤコビアン
【０３００】
【数１７０】
2004527027.tif 000172
の新たな列である。
１０．質量行列を用いて、先行する工程の結果
【０３０１】
【数１７１】
2004527027.tif 000173
を後退代入法に供して、所望の
【０３０２】
【数１７２】
2004527027.tif 000174
を得る。
後退代入法の演算は、残差ベクトルを
【０３０３】
【数１７３】
2004527027.tif 000175
ベクトルに処理することによって、直接形式方法において達成される。第２運動学的計算は、ヤコビアン全体に対して１回のみ実施される必要がある。なぜなら、この後退代入法は、その状態の正常値においてなされるからである。実際に、第２運動ルーチンは、ｚを計算している間に、変数が依然として隠されるように、工程１において呼び出されなければならない。
【０３０４】
以下の工程１１〜１３は、Ｊｕｕの列を満たすために使用される：
１１．ｕ（速度の微分）を計算する：
このルーチンは、入力ベクトルｄｕを採用し、そしてｉｄｖｋ（ｋ）＝Ｈ＊（ｋ）ｄｕ（ｋ）を計算する。次いで、ｄＶ（ｋ）（各物体の空間的速度の微分）が計算される：
【０３０５】
【数１７４】
2004527027.tif 000176
１２．速度により誘導される微分
【０３０６】
【数１７５】
2004527027.tif 000177
を計算する。本明細書中に提示されるように、このルーチンは、速度に依存する外部負荷がない場合に対して特殊化される。残っている項は、内力のみの変化に起因するものである。外部負荷の変化が存在する場合でさえも、以前のように、ｄＴ（ｋ）の寄与を含むために必要とされるのみである。
【０３０７】
【数１７６】
2004527027.tif 000178
空間的加速度の微分
【０３０８】
【数１７７】
2004527027.tif 000179
を計算した後に、空間的慣性力の微分が計算される：
【０３０９】
【数１７８】
2004527027.tif 000180
１３．ｄρ（ｋ）（物体ｋについての結合残差微分）を計算する：
【０３１０】
【数１７９】
2004527027.tif 000181
このルーチンを実行した後に、ｄρ（ｋ）に格納される値は、残差ヤコビアン
【０３１１】
【数１８０】
2004527027.tif 000182
の新たな列である。
１４．質量行列を用いて、先行する工程の結果
【０３１２】
【数１８１】
2004527027.tif 000183
を後退代入法し、所望の
【０３１３】
【数１８２】
2004527027.tif 000184
を得る。
【０３１４】
後退代入法の演算は、残差ベクトルを
【０３１５】
【数１８３】
2004527027.tif 000185
ベクトルに処理することによって、直接方法において達成される。第２運動学的計算は、１回のみ実施される必要がある。なぜなら、この後退代入法は、その状態の正常値においてなされるからである。実際に、第二運動ルーチンは、ｚを計算している間に、変数が依然として隠されるように、工程１において呼び出されなければならない。
【０３１６】
上記工程は、力がｑへの依存のみを有する限り、解析的ヤコビアンの計算を完了する。これは、全ての原子間力がポテンシャル関数から誘導される、古典的な状況に当てはまる。速度依存性の力（例えば、単純な粘度ダンピング）に当てはまるためには、上記工程のいくつかが、以下のように改変される必要がある：
上記工程２において、本発明者らはまた、収縮した速度ヤコビアン
【０３１７】
【数１８４】
2004527027.tif 000186
を計算する必要がある。これは、ブロックの対角であり、これもまた計算されなければならない。
【０３１８】
上記工程６において、Ｔ１（ｋ）の計算は、収縮された速度ヤコビアンで増大されなければならない：
【０３１９】
【数１８５】
2004527027.tif 000187
次いで、工程１１の後に工程が追加される。この工程は、工程１１ａと呼ばれる。この新たな工程は、ｄＴ（ｋ）を計算する：
【０３２０】
【数１８６】
2004527027.tif 000188
上記工程１２を実行する間、
【０３２１】
【数１８７】
2004527027.tif 000189
についての最後のループは、速度依存性の力の微分ｄＴ（ｋ）を除くことによって、改変される：
【０３２２】
【数１８８】
2004527027.tif 000190
これらの工程の残りは、同じままである。
【０３２３】
図６は、上で詳細に記載された、解析的ヤコビアン方法の実行工程を要約する。
【０３２４】
図７は、例示的なＭＤ系についての解析的ヤコビアンの正確さに対する、数値ヤコビアンの正確さのプロットを示す。摂動が完全に選択された最良の場合において、数字ヤコビアンからの一般化座標（ｑ）および一般化速度（ｕ）についての正確さの数字（線１５２によって示される）は、依然として、線１５０によって示される解析的ヤコビアンのものの半分のみであった。
【０３２５】
（さらなる実施形態）
本発明は、上記に記載された実施例に加えて、多くの実施形態を含む。以下の列挙は、他の実施形態および適用を有する。
【０３２６】
・ヤコビアンに含まれる力の次数
ヤコビアンに含まれる力の任意の次数としては、次数（Ｎ）、次数（Ｎ２）、次数（Ｎ３）、および次数（Ｎ４）が挙げられるが、これらに限定されない。次数（Ｎ）の力場の例は、次数（Ｎ２）である直接方法ではなく、第一の多極展開方法を使用する、静電力場である（例えば、Ｇｒｅｅｎｇａａｒｄ，Ｔｈｅ  Ｒａｐｉｄ  Ｅｖａｌｕａｔｉｏｎ  ｏｆ  Ｐｏｔｅｎｔｉａｌ  Ｆｉｅｌｄｓ  ｉｎ  Ｐａｒｔｉｃｌｅ  Ｓｙｓｔｅｍｓ，Ｍａｓｓａｃｈｕｓｅｔｔｓ  Ｉｎｓｔｉｔｕｔｅ  ｏｆ  Ｔｅｃｈｎｏｌｏｇｙ  Ｄｉｓｓｅｒｔａｔｉｏｎ，１９８８を参照のこと）。
【０３２７】
・直接形式のための解析的ヤコビアン
支配的な等式が直接形式である場合、いわゆる「順方向力学」形態の等式が得られる。この形式において、これらの等式は、状態ベクトルおよび適用される効果を処理し、そしてその系においてモデル化された結合の各々において加速度を発生させる。
【０３２８】
【数１８９】
2004527027.tif 000191
次いで、ヤコビアンは、状態ベクトルの要素に関して、加速度の偏微分を表す。好ましい実施形態は、これらの偏微分を計算するための、いくつかのアルゴリズム方法を示す。これらの方法は正確であり、そして微分を形成するために数値的近似を利用しない。
【０３２９】
直接形式は、ヤコビアンの
【０３３０】
【数１９０】
2004527027.tif 000192
分配
【０３３１】
【数１９１】
2004527027.tif 000193
を、
【０３３２】
【数１９２】
2004527027.tif 000194
関数を計算する関数の対アルゴリズム（ａｌｇｏｒｉｔｈｍｉｃ  ｃｏｕｎｔｅｒｐａｒｔ）を使用することによって、生じる。
【０３３３】
Ｍ−１ｆの計算が、次数（Ｎ）の浮動小数点演算（ｆｌｏｐｓ）において、質量行列の演算子の逆を利用することによって、達成される：
【０３３４】
【数１９３】
2004527027.tif 000195
ここで、全てのブロック行列は、第２運動学的計算において先に規定され、そして各因子は、次数（Ｎ）ｆｌｏｐｓにおいてｎベクトルに適用され得る演算子を表す。
【０３３５】
【数１９４】
2004527027.tif 000196
であるので、連鎖法則から、
【０３３６】
【数１９５】
2004527027.tif 000197
であり、Ｊｕｕについても同様である。
【０３３７】
従って、本発明は、上記項の各々を閉じた形態で計算するアルゴリズムを作成する際の、直接形式の等式のために使用され得る。
【０３３８】
従って、本発明は、多くの利点を提供する。解析的ヤコビアンは、よりずっと正確である（有効数字の桁の２倍）。直接形式ではなく残差形式からのカップリングは、よりずっと効率的である。「原子の数」から「物体の数」への、行および列の「収縮」は、力ヤコビアン行列の大きさを減少させる。ヤコビアン計算は、各列が摂動した場合には過剰に高い次数であるよりむしろ、力の計算と同じ次数である。従って、力が次数（Ｎ３）の演算子で計算される場合、例えば、数値ヤコビアンは、次数（Ｎ４）を必要とし、一方で解析的ヤコビアンは、次数（Ｎ３）の演算子のみを必要とする。ヤコビアン計算におけるループ構造の範囲を制御することによって、計算はなおさらに減少され得る（外側物体についてちょうど計算する）。
【０３３９】
従って、上述のものは、本発明の実施形態の完全な説明であるが、種々の改変、変更および均等物が、作製および使用され得ることが、明らかであるはずである。従って、上記説明は、添付の特許請求の範囲の組み合わせおよび境界によって規定される、本発明の範囲を限定するとはみなされるべきではない。
【図面の簡単な説明】
【図１】
図１は、本発明に従うソフトウエアシステムアーキテクチャのブロックモジュール図である。
【図２】
図２は、本発明に従う分子モデルの多体系のツリー構造を示す。
【図３】
図３は、図２の多体系の基準配置を示す。
【図４Ａ】
図４Ａは、図２の多体系における２つの物体間のスライダージョイントを示す。
【図４Ｂ】
図４Ｂは、図２の多体系における２つの物体間のピンジョイントを示す。
【図４Ｃ】
図４Ｃは、図２の多体系における２つの物体間のボールジョイントを示す。
【図５】
図５は、解析的ヤコビアン計算について使用される残差形式法についての一般的な計算工程を示す。
【図６】
図６は、解析的ヤコビアン計算のための一般的計算工程を要約するチャートである。
【図７】
図７は、数値的ヤコビアンを上回る、解析的ヤコビアンの有効数字のプロットである。
【図１】
2004527027.tif 000202
【図５】
2004527027.tif 000203
【図６】
2004527027.tif 000204
【図７】
2004527027.tif 000205
【国際公開パンフレット】
2004527027.tif 000206
2004527027.tif 000207
2004527027.tif 000208
2004527027.tif 000209
2004527027.tif 000210
2004527027.tif 000211
2004527027.tif 000212
2004527027.tif 000213
2004527027.tif 000214
2004527027.tif 000215
2004527027.tif 000216
2004527027.tif 000217
2004527027.tif 000218
2004527027.tif 000219
2004527027.tif 000220
2004527027.tif 000221
2004527027.tif 000222
2004527027.tif 000223
2004527027.tif 000224
2004527027.tif 000225
2004527027.tif 000226
2004527027.tif 000227
2004527027.tif 000228
2004527027.tif 000229
2004527027.tif 000230
2004527027.tif 000231
2004527027.tif 000232
2004527027.tif 000233
2004527027.tif 000234
2004527027.tif 000235
2004527027.tif 000236
2004527027.tif 000237
2004527027.tif 000238
2004527027.tif 000239
2004527027.tif 000240
2004527027.tif 000241
2004527027.tif 000242
2004527027.tif 000243
2004527027.tif 000244
2004527027.tif 000245
2004527027.tif 000246
2004527027.tif 000247
2004527027.tif 000248
2004527027.tif 000249
2004527027.tif 000250
2004527027.tif 000251
2004527027.tif 000252
2004527027.tif 000253
2004527027.tif 000254
2004527027.tif 000255
2004527027.tif 000256
2004527027.tif 000257
2004527027.tif 000258
2004527027.tif 000259
2004527027.tif 000260
【国際調査報告】
2004527027.tif 000261