JP6840668B2 - 結合された粒子間の相互作用の抑制 - Google Patents

Description

本開示は、分子動力学のシミュレーションに関し、特に、粒子対間の相互作用の不必要な計算を回避し、粒子を表すデータのキューへの流れおよびキューからの流れを管理することに関する。

分子動力学シミュレーション装置は、互いに相互作用する粒子からなる系の進化をシミュレートするものである。粒子間の相互作用は離れた位置で作用する力によって生じる。そのような力の一例は重力である。別の例は電気力である。

粒子動力学のシミュレーションでは、通常、粒子に作用する全ての力を合計し、ニュートンの法則を使用して、それらの力の合計の結果として粒子の動きを決定する。原則、各粒子は、系内の他の全ての粒子からの力を受ける。しかしながら、これらの力の大きさは、一般には、粒子間の距離の２乗で減少する。したがって、実際には、粒子が差異を生じさせるのに十分近い場合にのみ計算を実行することが好ましい。

公知のシミュレーション装置の構造および機能の詳細については、その内容が参照により本明細書に組み込まれる特許文献１および特許文献２において見ることができる。

米国特許第８，１２６，９５６号明細書 国際公開第２００６／１１３８２５号パンフレット

本発明は、バウアーズ（Bowers）らによる「DETERMINING COMPUTATION UNITS FOR COMPUTING MULTIPLE BODY INTERACTIONS」を名称とし、２０１２年２月２８日に発行された特許文献１に記載の装置および方法の改良に関する。これらの改良は、粒子間の不要な計算の抑制と、シミュレートされた粒子のキューへの流れの管理の改良とを含む。

本発明は、距離以外にも粒子間の相互作用の計算を抑制する他の理由が存在するという認識に部分的に基づいている。すなわち、粒子が互いに非常に接近している場合であっても、それらの粒子間の相互作用の計算を抑制する理由が存在し得る。したがって、本発明は、このような粒子対を識別し、そのような粒子対間の相互作用の計算を抑制し、それによって分子動力学シミュレーション装置の計算性能を向上させる方法を提供する。

第１の態様において、本発明は、コンピュータによって２つの粒子間のトポロジ距離が閾値未満であると判定することを含む方法を特徴とする。
いくつかの実施形態は、この決定に基づいて、コンピュータによる２つの粒子間の相互作用の計算を抑制する方法を含む。

他の実施形態は、粒子間の関係を表すコンピュータ読み取り可能なグラフを提供することを含む。トポロジ距離は、グラフ内の２つの粒子間のエッジ数に少なくとも部分的に依存する。ある実施形態において、粒子は第１の原子と第２の原子を含む。また、いくつかの実施形態において、２つの粒子間の相互作用の計算を抑制することは、計算が抑制されない場合に、分子動力学シミュレーションシステムによって実行される計算を抑制することを含む。

他の実施形態は、２つの粒子間のトポロジ距離が閾値未満であると判定することが、２つの粒子のうちの第１の粒子と、その一部が第１の粒子である粒子セットとの間のトポロジの関係を示す第１のトポロジ識別子により第１の粒子をタグ付けすることを含む。ある実施形態では、第１の粒子が原子であり、粒子セットが粒子対間の共有結合よりなる分子である。

追加の実施形態は、２つの粒子のうちの第２の粒子を第２のトポロジ識別子によりタグ付けすること、第１および第２のトポロジ識別子に基づいて第１および第２の粒子間のトポロジ距離を決定することを含む。これらの実施形態の中には、第１および第２のトポロジ識別子の各々が主鎖（backbone）識別子を含むもの、第１および第２のトポロジ識別子の各々が少なくとも１つの側鎖（side-chain）識別子を含むもの、第１および第２のトポロジ識別子の各々が末端（termination）フラグを含むもの、第１および第２のトポロジ識別子の各々が、主鎖識別子、一次側鎖識別子、二次側鎖識別子、および末端フラグを含むものが含まれる。

トポロジ距離を計算する種々の方法が特許請求の範囲に含まれる。例えば、ある粒子について、第１および第２のトポロジ識別子の主鎖識別子間の差の大きさと、第１および第２のトポロジ識別子の一次側鎖識別子間の差の大きさと、第１および第２のトポロジ識別子の二次側鎖識別子間の差の大きさと、第１および第２のトポロジ識別子の末端識別子間の差の大きさとを合計することによって、トポロジ距離を計算することを含む。他の実施形態は、第１および第２のトポロジ識別子の一次側鎖識別子間の差の大きさと、第１および第２のトポロジ識別子の二次側鎖識別子間の差の大きさと、第１および第２のトポロジ識別子の末端識別子間の差の大きさとを合計することによってトポロジ距離を計算することを含む。さらに他の実施形態は、第１および第２のトポロジ識別子の二次側鎖識別子間の差の大きさと、第１および第２のトポロジ識別子の末端識別子間の差の大きさとを合計することによってトポロジ距離を計算することを含む。

いくつかの実施形態は、グラフが環状部分を有すると判定し、グラフ内に表された粒子にトポロジ識別子を割り当てる前に環状部分を除去するようにグラフを修正する追加ステップを含む。

別の態様において、本発明は、コンピュータシステムにより実行されることで上述した方法のいずれかを実行するソフトウェアを符号化した非一時的かつ有形のコンピュータ可読媒体を特徴とする。

別の態様において、本発明は、上述した方法のうちのいずれかを実行するように構成されたコンピュータシステムを備える装置を特徴とする。
別の態様において、本発明は、２つの粒子間のトポロジ距離が閾値未満であると判定することを含む方法を特徴とする。

さらに別の態様において、本発明は、上述した方法のいずれかを実行するように構成された分子動力学シミュレーション装置を特徴とする。
第２の態様において、本発明は、ペアワイズ点相互作用（pairwise-point-interaction）モジュールアレイへの粒子の流れを管理する方法を特徴とする。このような方法は、第１の粒子セットを第１のキュー内に受け取ることを含む。第１の粒子セットは、現在の期間中にペアワイズ点相互作用モジュールアレイに受け渡される全ての粒子を含む第２の粒子セットのうちの適切なサブセットである。方法は、第２の粒子セットから全ての粒子を受け取る前に、第１の粒子セットの粒子を第１のキューからアレイへと受け渡すことを含む。

一実施形態において、方法は、第２の粒子セットからの追加の粒子が第１のキュー内に受け取られるときに第１の粒子セットの粒子をアレイにロードし続けることを含む。
別の実施形態において、方法は、第３の粒子セットを第２のキュー内に受け取ることを含む。この第３の粒子セットは、現在の期間中にアレイにのみロードされる全ての粒子を含む。そのような実施形態において、第１の粒子セットの粒子を第１のキューからアレイへと受け渡すことは、第３の粒子セットの全ての粒子がアレイにロードされた後にのみ生じる。

さらに別の実施形態において、方法は、第３の粒子セットを第２および第３のキュー内に受け取ることを含む。第３の粒子セットは、現在の期間中にアレイにロードされる全ての粒子を含む。この実施形態において、第１の粒子セットの粒子を第１のキューからアレイへと受け渡すことは、第３の粒子セットの全ての粒子がアレイにロードされた後にのみ生じる。

さらなる実施形態において、方法は、第４の粒子セットを第３のキュー内に受け取ることを含む。第４の粒子セットは、アレイにロードされるとともにストリーミングされる。この実施形態において、第４の粒子セットの粒子のストリーミングは、第４の粒子セットの粒子のロードが完了した後にのみ開始される。

追加の実施形態は、第１のキューが論理キューであることを含む。このような実施形態は、複数の選択された物理キューから第１および第２の粒子セットを選択することを含む。

別の態様において、本発明は、コンピュータシステムにより実行されることで本発明の第２の態様に関連する方法を実行するソフトウェアを符号化した非一時的かつ有形のコンピュータ可読媒体を含む。

別の態様において、本発明は、本発明の第２の態様に関連する方法のうちのいずれか一つを実行するように構成されたコンピュータシステムを備える装置を含む。
別の態様において、本発明は、本発明の第２の態様に関連する方法のうちのいずれか一つを実行するように構成された分子動力学シミュレーション装置を含む。

本発明の上記および他の態様は、以下の詳細な説明および添付図面から明らかとなる。

例示的な分子動力学シミュレーション装置におけるいくつかのノードと、典型的なノードの構造を示す図。 フレックスタイルの構造を示す図。 図１の相互作用タイルの構造を示す図。 分子の構造式を示す図。 分子のトポロジ識別子を示す図。 分子の環状グラフをどのように修正してトポロジ識別子によるその原子のタグ付けに適応させるかを示す図。 図３の高スループット相互作用サブシステムタイルにおけるキューを示す図。 関数の区分近似に使用するためのテーブルを示す図。 図８のテーブルと実質的に同じ区分近似を実行するための論理を伴うより小さなテーブルを示す図。

分子動力学シミュレーションは、力に応答して粒子の運動をシミュレートすることに関する。これらの力の多くは短距離の力であるため、粒子に関する殆どの計算は近くの粒子との相互作用に限られている。したがって、近接粒子に関する計算は、多くの場合、他の近接粒子に関する計算とは独立して実行することができる。この特性は、並列処理に適している。

この固有の並列性の利点を利用するため、図１に示す分子動力学シミュレーション用のシミュレーション装置１０は、ネットワークにより接続された複数のノード１２を考慮する。これらのノードは集合的にシミュレーション空間の体積を表し、各ノードはそのシミュレーション空間の特定の部分に対応する。ネットワーク内におけるノード間のトポロジの関係は、シミュレーション空間におけるそれらの空間的関係に対応する。このシミュレーション体積は、シミュレーション装置１０によって計算される互いに相互作用する粒子によって占有される。特定の一実施形態において、シミュレーション体積はトロイダルである。しかしながら、プリズムのような他の体積も可能である。

固有の並列性のために、シミュレーション体積をノードボックスに分割することは有益であり、各ノードボックスはノード１２のちの一つによって処理される。異なるノード間で計算を割り当てる方法の説明は、ショー（Shaw）による、「ORTHOGONAL METHOD」を名称とし、２０１０年４月７日に発行された米国特許第７，７０７，０１６号に記載されており、その内容が参照により本明細書に組み込まれる。

各ノード１２は、ペアワイズ相互作用を実行するように特に設計された第１のハードウェア要素と、潜在的に相互作用する粒子を第１のハードウェア要素に提供する機能を有する第２のハードウェア要素との組み合わせとして実装される。一実施形態では、各ノード１２は、１６個のフレックスタイル１４と、２つの相互作用タイル１６と、ホストインターフェース１８と、論理アナライザ２０とを有する特定用途向け集積回路として実装される。

ホストインターフェース１８は、ＰＣＩリンクを介して外部ホストとの通信を提供する。論理アナライザ２０は、デバッグのためのノードアクティビティを捕捉して記憶するために主に使用される。また、各ノード１２は、３つの局所座標方向の各々において隣接ノード間のデータ伝送のための通信インタフェース２２を含む。ノード１２内において、そのノード１２の構成要素間のデータ伝送は、オンチップメッシュネットワークによって実行される。

図２に示されるように、各フレックスタイル１４は複数のジオメトリコア２４を有し、各ジオメトリコア２４は関連するメモリ２６を有する。図示する特定の実施形態では、このようなジオメトリコア２４が４つ存在する。しかしながら、他の数のジオメトリコアも使用可能である。また、各フレックスタイル１４は、ノード１２の他の構成要素との通信を可能にするネットワークインタフェース２８と、細粒度イベント駆動型（fine-grained event-driven）計算のためのハードウェアサポートを提供するディスパッチユニット３０とを含む。ディスパッチユニット３０については、グロースマン（Grossman）による、「EVENT-DRIVEN COMPUTATION」を名称とし、２０１４年９月１８日に出願された米国特許出願公開第２０１４／０２８２５７６号に記載されており、その内容が参照により本明細書に組み込まれる。また、フレックスタイルには、４つのジオメトリコア２４全てに利用可能な共通メモリ３２も含まれている。

図３に示すように、各相互作用タイル１６は、相互作用コントローラ３４と、簡易フレックスタイル４４と、ペアワイズ点相互作用モジュール（pairwise-point-interaction module：ＰＰＩＭ）アレイ５２とを特徴とする。ＰＰＩＭアレイ５２は、複数のＰＰＩＭ５４から構成される。図示する実施形態では、ＰＰＩＭアレイ５２は３８個のＰＰＩＭ５４を有する。しかしながら、他の数のＰＰＩＭも使用可能である。

相互作用タイル１６の特に相互作用コントローラ３４は、オンチップメッシュネットワークを介してフレックスタイル１４から粒子および格子点を受け取る。そして、これらの粒子および格子点をローカルメモリ３８に格納されたキュー３６にエンキューする。相互作用コントローラ３４の動作は、簡易フレックスタイル４４内のジオメトリコア４２から受信した命令４０によって制御される。

ジオメトリコア４２に加えて、簡易フレックスタイル４４は、ローカルメモリ４６と、ディスパッチユニット４８と、ネットワークインタフェース５０とを有し、これらの全ては、図２に関連して説明したものと同様の機能を果たす。

シミュレーション装置１０は、小さな時間ステップで全ての粒子間力を繰り返し計算し積分することによって、粒子の集合の進化をシミュレートする。各時間ステップの開始時に、各フレックスタイル１４は、オンチップネットワークを使用して、同じノードおよび他のノード上の相互作用タイル１６にパケットを送信する。これらのパケットは相互作用タイル１６が最終的に計算によって明らかにする相互作用粒子に関する情報を含む。

上述したように、各ノード１２は、そのノードボックス内の粒子に関する計算を担う。リソース制限のために、各ノードボックスを複数のサブボックスにさらに分割する必要があり得る。各サブボックスはいくつかの粒子集団を有する。この粒子集団は、他の粒子によって加えられる力に応じて、シミュレーション体積内を移動する粒子の結果として経時的に変動する。各相互作用タイル１６は、フレックスタイル１４から複数のサブボックスのいくつかの可変数の粒子を受け取る。また、各相互作用タイル１６は、フレックスタイル１４から複数のカウントパケットを受信する。各カウントパケットは、各サブボックスからの予想される粒子数を通知する。これらのカウントパケットは、相互作用タイル１６によって用いられ、全てのサブボックスから全ての粒子を受信したときを決定する。

粒子の動きを正確にシミュレートするためには、非常に近い時間間隔で粒子の相互作用を評価する必要があることが分かっている。典型的なシミュレーションでは、数フェムト秒ごとのシミュレーション時間ごとに粒子相互作用が評価される。

多くの興味深い事象がミリ秒以上の時間スケールで生じるか、数十万個の粒子が関わって生じるか、あるいはその両方が起こる。この種の事象の例には、生物学的高分子が相互作用する生化学系において生じるものが含まれる。

このような長い時間スケールで事象が展開される系のシミュレーションは大量のリアルタイム処理を必要とする。これは、シミュレーション時間と実際の時間との変換率が現在１０億分の１のオーダーであるためである。したがって、シミュレート時間をフェムト秒だけ進めるのに必要な相互作用を計算するのにマイクロ秒費やす必要がある。これは速く見えるが、この問題に注視すると、わずか１ミリ秒のリアルタイムシミュレーションをこの変換率で行うには、シミュレーション装置１０は１００万秒動作する必要があり、これは１１日間と少しに亘って連続計算するものとなる。

以下の説明では、粒子は、特定のハードウェアに格納されているかストリーミングされているものとする。なお、「粒子」は、実際の粒子ではなく、粒子を表す情報またはシミュレートされた粒子を意味する。

粒子は、特定の順序で相互作用タイル１６に到着する。粒子が到着すると、粒子ディレクタ７４は、粒子を図７に示すような異なる先入れ先出し（ＦＩＦＯ）キュー７６に入れる。一実施形態では、このようなキュー７６が２^７個ある。一実施形態では、粒子ディレクタ７４は、プログラム可能なルックアップテーブルとして実装される。

粒子は不特定の順序で到着するので、不特定の順序でキュー７６に入れられる。しかしながら、キュー７６が使用される順序は予め分かっている。また、各キューは、予期される粒子分のサブボックスの数を把握するように予めプログラムされている。また、カウントパケットがフレックスタイル１４から受信されているので、各キューは、サブボックス毎に到着する粒子数も決定することができる。したがって、各キュー７６は、予期される粒子の数およびそれら粒子が全て到着したか否かを把握する。

キュー７６は、第１のキュー７８、第２のキュー８０、および第３のキュー８２に分割される。ＰＰＩＭアレイ５２にロードされるべき粒子は第１のキュー７８に入れられる。ロードされる粒子を通過してストリーミングされるべき粒子は、第２のキュー８０に入れられる。ロードされストリーミングされるべき粒子は、第３のキュー８２に入れられる。

例えば、第１の空間体積内にのみ存在する粒子は第１のキュー７８に入れられ、第１のキュー７８からＰＰＩＭアレイ５２にロードされる。また、第２の空間体積内にのみ存在する粒子は第２のキュー８０に入れられ、第２のキュー８０からＰＰＩＭアレイ５２にストリーミングされる。そして、第１および第２の空間体積の交点にある粒子は第３のキュー８２に入れられ、第３のキュー８２からＰＰＩＭアレイ５２にロードされストリーミングされる。

一実施形態では、相互作用コントローラ３４は、第１および第３のキュー７８，８２が満たされるまで待機する。上述したように、フレックスタイル１４は粒子とともにカウントパケットを送信しているので、この情報が利用可能である。次いで、相互作用コントローラ３４は、第１および第３のキュー７８，８２の全ての粒子をラウンドロビン方式でＰＰＩＭ５４にロードする。粒子がロードされると、ＰＰＩＭアレイ５２は、第２のキュー８０、さらには第３のキュー８２からの粒子のストリーミングを受け入れる準備が整う。キューがその粒子をＰＰＩＭアレイ５２に入れることを開始すべきか否か、または待機すべきか否かの制御は、ソフトウェア命令によって実行される。

ロードが完了したことを検出すると、相互作用コントローラ３４は、第２のキュー８０が満たされているか否かを判定する。ＰＰＩＭアレイ５２へのロードは完了していないので、この時点によって第３のキュー８２が満たされていることは把握されている。第２のキュー８０が満たされると、相互作用コントローラ３４は、粒子を第２のキュー８０および第３のキュー８２から引き出す。その後、これらの粒子は、ＰＰＩＭアレイ５２の全てのＰＰＩＭ５４を通ってストリーミングされ始める。そして、ＰＰＩＭ５４は、第１および第３のキュー７８，８２におけるロードされた粒子と、第２および第３のキュー８０，８２におけるストリーミングされた粒子との間の相互作用を計算する処理を開始する。

ＰＰＩＭアレイ５２への粒子の流れを管理する別の方法において、相互作用コントローラ３４は、第１および第３のキュー７８，８２が満たされるまで敢えて待機しない。代わりに、相互作用コントローラ３４は、可能な限り早く、第１および第３のキュー７８，８２からの粒子をＰＰＩＭ５４にロードする。これにより、ＰＰＩＭアレイ５２は、粒子のロードを開始する前に、全ての粒子が到着するまで待機する必要がなくなるため、時間が節約される。

粒子が全てＰＰＩＭアレイ５２にロードされると、相互作用コントローラ３４は、第２のキュー８０および第３のキュー８２からＰＰＩＭアレイ５２への粒子のストリーミングを開始する。これは、第２のキュー８０がストリーミング用に指定された粒子をまだ全て受信していない場合にも行われる。

別の実施形態では、図７に示されたキューを、全ての物理キューのセットのうちのサブセットからキューを選択的に決定することによりアセンブルされる論理キューとして考慮することが有用である。これらの論理キューに入れられる粒子は、物理キューのセットのうちのサブセットを選択することによって取得される。選択すべきサブセットの決定は、対象となっているサブボックスの近傍に依存する。これにより、必要に応じて物理キューを柔軟に再利用することが可能となる。

ＰＰＩＭ５４にストリーミングされる間に、第２および第３のキュー８０，８２からの各ストリーミング粒子は、第１および第３のキュー７８，８２からの各ロード粒子と接近する。このような各接近時に、相互作用を計算すべきか否かを判定する必要がある。この判定が「肯定」である場合、ロード粒子とストリーミング粒子との間の相互作用から生じる力が計算される。これらの力は、計算されるごとに累積される。ストリーミング粒子がＰＰＩＭアレイ５２を通過すると、粒子に作用する全ての力が出力５６に送られる。

従って、「これら２つの粒子間の相互作用を計算すべきか」との論題に答えるための試験を考慮しなければならない。
相互作用を計算すべきか否かを判定するための１つの試験は、粒子間の距離が計算すべき価値があるほど十分に近いか否かを調べることである。２つの粒子が離れすぎていると、相互作用は計算されない。

しかしながら、この試験は良い近似ではあるが、粒子同士が近接している場合であっても粒子間の相互作用を計算すべきではない場合があることによって複雑になる。このような複雑さは、原子が共有結合して分子を形成するという理由により、分子動力学において生じる。その場合、これらの粒子を結合で保持する力は、シミュレートされる粒子間力を小さなものとし易い。既知のシミュレーション装置では、これらの相互作用は計算されるが、後に訂正パイプラインで削除される。

この欠陥を是正して訂正パイプラインの必要性を排除するために、各粒子は、その粒子と他の近くの粒子との間のトポロジの関係の性質を伝達するトポロジ識別子に関連付けられる。一般性を失うことなく、このトポロジ識別子は、分子を形成するために他の原子と結合する原子、特にタンパク質や脂質のような生物学的高分子に関連して議論される。

図４に示すように、分子５８は、ノードが原子６０を表し、エッジが原子６０間の結合６２を表すグラフとして表すことができる。分子５８内の２つの原子６０を分離するエッジの数を用いて原子６０間のトポロジ距離を規定することができる。そして、２つの原子間の相互作用の計算を抑制する基準としてこのトポロジ距離を用いることができる。その結果、分子内で隣り合っている２つの原子について相互作用は計算されない。しかしながら、２つの原子が高分子の両側の端にある場合、相互作用を計算することができる。

図５に示すように、トポロジ識別子６４は、原子６０に固有の短い記述子である。この一意性は異なる分子にまたがっている。
トポロジ識別子６４は、分子５８内の原子６０間のトポロジの関係を符号化する。したがって、２つの原子６０のトポロジ識別子６４を比較することによって、２つの原子６０間のトポロジ距離を規定することができる。これにより、通常であれば相互作用が計算される程度に十分に原子同士が近づいている場合であっても、原子６０間の相互作用が排除されるべきか否かを高精度に判定することが可能となる。

このため、トポロジ識別子６４を比較することにより、通常であれば訂正パイプラインで行われた殆どの訂正を回避し、無駄な計算も回避することができる。以下に説明するように、計算が実行されるべきではない場合でも計算が実行される特定の場合がある。しかしながら、個別のハードウェア訂正パイプラインを有する代わりに、ソフトウェアによって計算の訂正が行われる場合はごくわずかである。

本明細書において説明するトポロジ識別子６４の使用は、多くの分子５８において原子の主鎖（backbone）がその主鎖から分岐した一次側鎖（primary side-chain）を有するという特徴に依存する。これらの一次側鎖は二次側鎖（secondary side-chain）を有し得る。また、これらの二次側鎖は三次側鎖を有し得る。このように側鎖は無限に存在し得る。しかしながら、対象となる殆どの分子５８は、主鎖から分岐した一次側鎖と一次側鎖から分岐した二次側鎖とを有する主鎖を有することが見出されている。したがって、実際的な実施形態では、一次側鎖および二次側鎖のみを考慮することが必要である。

本明細書において説明する実施形態では、各原子６０に整数四重項を割り当てる。四重項の要素は、主鎖、一次側鎖、二次側鎖、および末端（termination）フラグを識別する。一般的には、トポロジ識別子は、Ｎ＋２個の要素を有する整数タプルであり、Ｎは、考慮すべき側鎖のレベルの数である。本実施形態では、一次側鎖および二次側鎖のみが考慮されるのでＮ＝２である。

任意の原子６０について、四重項の第１の要素は、その原子の関連する主鎖原子を識別する主鎖識別子である。第１の要素を割り当てるために、原子６０自体が主鎖原子である場合、その原子６０はそれ自身に関連するとみなされる。最終的に同じ主鎖原子に接続する側鎖の全ての原子は、同じ主鎖識別子を有する。

四重項の第２の要素は、原子６０に関連する一次側鎖を識別する一次側鎖識別子である。一次側鎖は、主鎖に直接結合される原子６０を有する。同じ側鎖に存在するか、またはその同じ側鎖に接続された側鎖に存在する全ての原子は、同じ一次側鎖識別子を有する。

四重項の第３の要素は、原子に関連する第２のレベルの側鎖を規定する二次側鎖識別子である。
最後に、第４の要素は、原子が末端原子であるか否かを識別する末端フラグである。末端原子は、一つの共有結合のみによって分子の残りに結合しているものである。本明細書で使用する「共有結合」という用語は、結合に関与する電子軌道の数とは無関係であり、したがって、二重結合や三重結合を含む。

上述した方法は本質的に再帰的であり、末端フラグと主鎖識別子との間に適切な要素を単に追加することによって任意の数の側鎖に拡張することができることは明らかである。
図示された実施形態では、トポロジ識別子は整数四重項（ｎ，ｍ，ｋ，ｔ）である。主鎖原子は（ｎ，０，０，０）として識別される。ｎ番目の主鎖原子から離れた一次側鎖中の原子は（ｎ，ｍ，０，０）として識別され、ここで、ｍは、その原子と主鎖原子との間の鎖に沿った距離を表す１以上の整数である。二次側鎖の原子は（ｎ，ｍ，ｋ，０）として識別され、ここで、ｋは、その原子と二次側鎖が一次側と交差する原子（ｎ，ｍ，０，０）との間の距離を表す１以上の整数である。隣接するものが１つのみである末端原子の末端フラグは１に設定される。したがって、隣接する原子が原子（ｎ，ｍ，ｋ，０）である末端原子は、そのトポロジ識別子６４として整数四重項（ｎ，ｍ，ｋ，１）を有する。

図５は、上述した規則によって割り当てられた関連するトポロジ識別子を有する原子を備えた高分子の一部に対応するグラフを示す。図示のグラフは、４〜８まで付番された５つの主鎖原子を示している。すなわち、上述した規則に従って、これらの主鎖原子（ｎ，０，０，０）にはｎ＝４からｎ＝８まで番号が付けられている。

主鎖原子４に結合された原子は１つの隣接原子しか有さないため末端原子である。規則に従って、その四重項は、末端フラグが１に設定されていることを除いて、その隣接原子の四重項、つまり主鎖原子８の四重項と同じである。また、末端原子は、同じ規則に従って構築された対応する四重項を備える主鎖原子６，７に結合されていることが分かる。

主鎖原子５の外側の側鎖にある原子の識別子を検査することによって、側鎖に対するトポロジ識別子の割り当てが分かる。図示のように、最終的に主鎖原子５に接続する全ての原子の各四重項はｎ＝５を有する。同じ一次側鎖の全ての原子は同じｍの値を有するが、同じ二次側鎖の全ての原子は同じｋの値を有する。

整数四重極（ｎ_１，ｍ_１，ｋ_１，ｔ_１）および（ｎ_２，ｍ_２，ｋ_２，ｔ_２）を有する２つの原子間のトポロジ距離は３つの場合に分けられる。
第１の場合、原子は異なる主鎖原子に結合されており、このため、異なる主鎖識別子を有する。これは、ｎ_１≠ｎ_２であることを意味する。このような場合、距離は、主鎖識別子間の差の大きさを取って、それを２つの整数四重子の残りの要素全ての合計に加えることによって、
｜ｎ_１−ｎ_２｜＋ｍ_１＋ｍ_２＋ｋ_１＋ｋ_２＋ｔ_１＋ｔ_２
として得られる。

第２の場合、２つの主鎖識別子は同じであるが、２つの原子は異なる一次側鎖上にある。したがって、ｎ_１＝ｎ_２であり、ｍ_１≠ｍ_２である。その場合、距離は第１の場合と同じであるが、一次側鎖識別子を加算する代わりに、その差の大きさを評価して、
｜ｎ_１−ｎ_２｜＋｜ｍ_１−ｍ_２｜＋ｋ_１＋ｋ_２＋ｔ_１＋ｔ_２
として得られる。

第３の場合、２つの原子は一次側鎖上にあるが、それらは異なる二次側鎖上にある。これは、ｎ_１＝ｎ_２、ｍ_１＝ｍ_２であり、ｋ_１≠ｋ_２であることを意味する。その場合、距離は第２の場合と同じ方法で計算されるが、二次側鎖識別子を加算する代わりに、その差の大きさを評価して、
｜ｎ_１−ｎ_２｜＋｜ｍ_１−ｍ_２｜＋｜ｋ_１−ｋ_２｜＋ｔ_１＋ｔ_２
として得られる。

トポロジ識別子６４を割り当てる上述の方法は、分子５８のグラフが、最大一つの側鎖が任意の原子から発する非環状グラフであると仮定したものである。これは、対象とする殆どの分子に当てはまるが、例外もある。そのような場合、この条件が満たされるまで、分子のグラフから少数のエッジが除去される。

図６は、環状グラフを有する分子６８の例を示す。分子のグラフが環状であるかまたは環状部分を含む場合、トポロジ識別子を割り当てる上述の方法は機能しない。この場合に適応するために、グラフから２つのエッジ７０，７２を除去して環状部分を開放し、対応する非環状分子７４を生成する必要がある。この不規則性に対処するためにソフトウェアで後に修正が行われる。

原子のトポロジ識別子およびその位置は、ＰＰＩＭアレイ５２を通過する過程で、その位置に束ねられる。その結果、符号化はできるだけコンパクトでなければならない。
実際の問題として、対象とする分子において、殆どの側鎖は短い。殆どの脂質やタンパク質について、一次側鎖識別子を符号化するには３ビットで十分であり、二次側鎖識別子を符号化するには１ビットで十分である。末端フラグを符号化するには更に１ビット必要となる。したがって、残りのビットを使用して主鎖識別子を符号化することができる。

符号化を実行するのに十分なビットが割り当てられていないために、全ての主鎖識別子を符号化するには化学系が大きすぎる場合がある。あるいは、利用可能なビットで側鎖を符号化できない場合がある。分子が環を持つ場合などのように、これらのいずれの場合においてもソフトウェアで修正する必要がある。

粒子間力などの量を評価する過程で、引数の関数を評価することを含む計算を実行することがしばしば必要となる。関数、特に超越関数の評価は時間のかかる作業である。この作業を高速化するために、ルックアップテーブル内の特定の引数の関数の値を単純に検索することが知られている。しかしながら、近似誤差が小さいテーブルは非常に大きくなる。別のアプローチは、関数の定義域を複数の部分に分割し、三次形式などのパラメトリック形式を使用して所望の関数を近似することである。そのような場合、テーブルは、定義域領域から三次多項式の係数へのマッピングを提供する。定義域の不均一区画を使用することで、三次多項式が関数と良好に一致しない定義域の部分、例えば近似されるべき関数が速く変化している定義域の部分に、より細かい区画を提供することが可能となる。逆に、一致が良好な場合は、より粗い区画を使用することができる。

関数の値がその定義域の或る領域について急激に変化するために、三次多項式を使用した場合でも、その定義域の部分の関数を適切にモデル化するために、三次係数のルックアップテーブルが非常に多くのエントリを必要とする場合がある。

図８は、三次多項式によって区分的に近似された関数のルックアップテーブル８４を示す。特定のｘ_ｏの値について関数ｕ（ｘ_ｏ）の値を得るために、ｘ_ｏを含む定義域の部分に対応するエントリ８６を探す。このエントリ８６は三次多項式の係数を生成し、適切な組合せ論理８８を用いて評価することにより関数ｕ（ｘ_ｏ）の値を得ることができる。

当然ながら、近似は区分的なものであるため、多項式の係数は近似関数の定義域全体にわたって変化する。このため、テーブル８４には複数のエントリがなければならない。図８から分かるように、その定義域の特定の領域については、近似関数が急激に変化し、テーブル内のエントリは非常に小さい領域に対してのみ有効である。もちろん、この問題はテーブルのエントリを増やし、各エントリが関数の定義域のより小さな部分をカバーすることで容易に解決することができる。しかし、このようなテーブルは大量のダイスペースを消費するのでこれを避けることが望ましい。

図９に示される別の実施形態では、テーブル８４は、近似関数が急激に変化するような近似関数の定義域の部分についても、より少ないエントリを有する。特定のｘ_ｏの値について関数ｕ（ｘ_ｏ）の近似値を得るために、ｘ_ｏを含む定義域の部分に対応するエントリ８６を探す。このエントリ８６は三次多項式の係数を生成する。これらの係数が適切な組合せ論理８８によって使用されることにより、引数ｘ_ｏにおける関数の値の第１の近似、すなわちｕ〜（ｘ_ｏ）（〜はチルダを表す）を生成することができる。この第１の近似を向上させるために、引数ｘ_ｏは、その引数に基づいて１つまたは複数のパラメータを出力するルックアップテーブル８９にも供給される。次いで、この１つまたは複数のパラメータは増加／減衰論理９０に供給される。増加／減衰論理９０は、論理で容易に実行され第２の近似をもたらし得る機能を有する。増加／遅延論理９０の出力は関数ｆ（ｘ_ｏ）の値である。この値と第１の近似ｕ〜（ｘ_ｏ）は乗算器９２に供給され、その出力が第２の近似ｕ（ｘ_ｏ）＝ｕ〜（ｘ_ｏ）ｆ（ｘ_ｏ）となる。好ましい実施形態では、ｆ（ｘ_ｏ）はｘ_ｏ ^ｋ・ｅｘｐ（−ｘ_ｏ ^ｍ）であり、ここで、ｋおよびｍは、特定のｘ_ｏの値についてルックアップテーブル８９から得られるパラメータである。

シャン（Shan）らによる、「GRID BASED COMPUTATION FOR MULTIPLE BODY SIMULATION」を名称とし、２００９年４月２８日に発行された米国特許第７，５２６，４１５号に記載されているＧＳＥ法を用いた電荷拡散および力補間に関連してＰＰＩＭアレイ５２を使用することが知られており、その内容が参照により本明細書に組み込まれる。このシャンに適用される方法の数学的基礎は、シャン（Shan）らによる「GAUSSIAN SPLIT EWALD: A FAST EWALD MESH METHOD FOR MOLECULAR SIMULATION」、J. Chem. Phy. １２２０５４１０１（２００５）に記載されており、その内容が参照により本明細書に組み込まれる。粒子および格子位置がＰＰＩＭアレイ５２にどのように分布するかの詳細は、ショー（Shaw）らによる、「PARALLEL COMPUTER ARCHITECTURE FOR COMPUTATION OF PARTICLE INTERACTIONS」を名称とし、２０１３年４月１１日に公開された米国特許出願公開第２０１３／００９１３４１号に記載されており、その内容が参照により本明細書に組み込まれる。

全ての格子位置がＰＰＩＭアレイ５２に渡されるこれらに記載された方法の改良において、代替方法は、格子位置がランダムに配置されるのではなく、実際にはある空間規則性を有することを利用する。この規則性を利用することにより、選択された位置をアレイに渡し、それらの選択された位置から格子位置を導出することが可能となる。

例えば、Ｘをｍ個の格子位置（ｘ_１，ｘ_２，…ｘ_ｍ）の集合とし、各ｘ_ｉは、シミュレーション空間に適した次元を有する位置ベクトルであるとする。従来技術の方法によれば、点格子位置をＸに渡すために、全てのｍ点をＰＰＩＭアレイ５２に渡す。

改良された方法では、ｎ個の位置（ｙ_１，ｙ_２，…ｙ_ｎ）の集合Ｙが存在する。ここで、ｎ＜ｍである。また、Ｘ＝Ｒ（Ｙ）となる規則Ｒも存在する。したがって、ＸではなくＹをＰＰＩＭアレイ５２に渡し、ＰＰＩＭアレイ５２に規則Ｒを使用してＸを再構成する。これは圧縮技術に相当し、圧縮の程度はｍ：ｎの割合に依存する。

例えば、集合Ｘ内の２つのベクトルが（ｘ，ｙ，ｚ_１）および（ｘ，ｙ，ｚ_２）であった場合は、ベクトル（ｘ，ｙ，（ｚ_１＋ｚ_２）／２）を含む集合Ｙを単純に渡す。格子間隔が分かっている場合は、集合Ｘから元の２つのベクトルを導くことができる。あるいは、集合Ｙは、集合Ｘからの１つおきの点と等しくてもよく、この場合、Ｙのベクトルの適切な座標に適切な格子間隔を加えることによって元の集合Ｘを再構成することができる。

一つの方法では、位置の集合ＹがＰＰＩＭアレイ５２にロードされ、原子がＰＰＩＭアレイ５２にストリーミングされる。この手順は電荷拡散を実行するために使用される。別の実施形態では、規則が逆転され、原子がＰＰＩＭアレイ５２にロードされる一方、位置の集合ＹがＰＰＩＭアレイ５２にストリーミングされる。この手順は力補間を実行するために使用される。

シミュレーション装置１０の一実施形態では、処理のために相互作用タイル１６に到達する粒子は、常に現在の時間ステップと関連付けられる。しかしながら、別の実施形態は、粒子データを伝送するデータパケットに関連付けられた位相ビットを相互作用タイル１６に導入する。位相ビットの値は、特定の時間ステップに関連付けられる。これにより、２つの異なる時間ステップで粒子を区別する方法が提供される。その結果、相互作用タイル１６が２つの時間ステップに関連付けられたデータパケットを受信することが可能となる。

位相ビットを提供する実施形態では、相互作用コントローラ３４は、位相ビットの各値に関連付けられたキュー３６を管理する。粒子を受け取ると、粒子の位相ビットが検査され、粒子が位相ビットに適したキューに入れられる。したがって、位相ビットは、相互作用タイル１６が異なる時間ステップに関連付けられたデータを受信することを可能にし、それによって相互作用タイル１６とフレックスタイル１４とを同期させる必要性を排除する。動作中、現在の時間ステップに対応するキューのみがＰＰＩＭ５４にロードされる。

本発明およびその好ましい実施形態を説明したが、新規であり特許で保護される対象は特許請求の範囲に記載される。

Claims (20)

1. コンピュータにより２つの粒子間のトポロジ関係を示すトポロジ識別子に基づいて２つの粒子間のトポロジ距離を計算すること、
   前記コンピュータにより前記トポロジ距離が閾値未満であると判定すること、
   前記トポロジ距離が閾値未満である場合に、前記コンピュータによる前記２つの粒子間の相互作用の計算を抑制すること、
   前記トポロジ距離が閾値以上である場合に、前記コンピュータによる前記２つの粒子間の相互作用の計算を行うことを備える方法。
2. 前記２つの粒子間の関係を表すコンピュータ読み取り可能なグラフを提供することをさらに備え、前記トポロジ距離は、前記グラフ内の前記２つの粒子間のエッジ数に少なくとも部分的に依存する、請求項１に記載の方法。
3. 前記２つの粒子は第１の原子と第２の原子であり、前記２つの粒子間の相互作用の計算を抑制することは、計算が抑制されない場合に、分子動力学シミュレーションシステムによって実行される計算を抑制することを含む、請求項１に記載の方法。
4. 前記２つの粒子間のトポロジ距離が閾値未満であると判定することは、前記２つの粒子のうちの第１の粒子と、その一部が前記第１の粒子である粒子セットとの間のトポロジの関係を示す第１のトポロジ識別子により前記第１の粒子をタグ付けすることを含む、請求項１に記載の方法。
5. 前記第１の粒子が原子であり、前記粒子セットが粒子対間の共有結合よりなる分子である、請求項４に記載の方法。
6. 前記２つの粒子のうちの第２の粒子を第２のトポロジ識別子によりタグ付けすること、
   前記第１および第２のトポロジ識別子に基づいて前記第１および第２の粒子間のトポロジ距離を決定すること、をさらに備える請求項４に記載の方法。
7. 前記第１および第２のトポロジ識別子の各々が主鎖識別子を含む、請求項６に記載の方法。
8. 前記第１および第２のトポロジ識別子の各々が少なくとも１つの側鎖識別子を含む、請求項６に記載の方法。
9. 前記第１および第２のトポロジ識別子の各々が末端フラグを含む、請求項６に記載の方法。
10. 前記第１および第２のトポロジ識別子の各々が、主鎖識別子、一次側鎖識別子、二次側鎖識別子、および末端フラグのうちの少なくとも１つを含む、請求項６に記載の方法。
11. 前記第１および第２のトポロジ識別子に基づいてトポロジ距離を計算することをさらに備え、前記トポロジ距離を計算することは、前記第１および第２のトポロジ識別子の主鎖識別子間の差の大きさと、前記第１および第２のトポロジ識別子の一次側鎖識別子間の差の大きさと、前記第１および第２のトポロジ識別子の二次側鎖識別子間の差の大きさと、前記第１および第２のトポロジ識別子の末端識別子間の差の大きさとを合計することを含む、請求項６に記載の方法。
12. 前記第１および第２のトポロジ識別子に基づいてトポロジ距離を計算することをさらに備え、前記トポロジ距離を計算することは、前記第１および第２のトポロジ識別子の一次側鎖識別子間の差の大きさと、前記第１および第２のトポロジ識別子の二次側鎖識別子間の差の大きさと、前記第１および第２のトポロジ識別子の末端識別子間の差の大きさとを合計することを含む、請求項６に記載の方法。
13. 前記第１および第２のトポロジ識別子に基づいてトポロジ距離を計算することをさらに備え、前記トポロジ距離を計算することは、前記第１および第２のトポロジ識別子の二次側鎖識別子間の差の大きさと、前記第１および第２のトポロジ識別子の末端識別子間の差の大きさとを合計することを含む、請求項６に記載の方法。
14. 前記グラフが環状部分を有すると判定すること、
    前記グラフ内に表された粒子にトポロジ識別子を割り当てる前に前記環状部分を除去するように前記グラフを修正すること、をさらに備える請求項２に記載の方法。
15. 分子動力学シミュレーション装置において、相互作用コントローラからペアワイズ点相互作用モジュールのアレイへの粒子の流れを管理する方法であって、
    前記相互作用コントローラが、第１の粒子セットを第１のキュー内に受け取ることであって、前記第１の粒子セットは、現在の期間中にペアワイズ点相互作用モジュールのアレイに受け渡される全ての粒子を含む第２の粒子セットのうちの適切なサブセットである、前記第１の粒子セットを第１のキュー内に受け取ること、
    前記相互作用コントローラが、前記第２の粒子セットから全ての粒子を受け取る前に、前記第１の粒子セットの粒子を前記第１のキューから前記アレイへと受け渡すこと、を備える方法。
16. 前記第２の粒子セットからの追加の粒子が前記第１のキュー内に受け取られるときに前記第１の粒子セットの粒子を前記アレイにロードし続けることをさらに備える請求項１５に記載の方法。
17. 第３の粒子セットを第２のキュー内に受け取ることをさらに備え、前記第３の粒子セットは前記現在の期間中に前記アレイにのみロードされる全ての粒子を含み、
    前記第１の粒子セットの粒子を前記第１のキューから前記アレイへと受け渡すことは、前記第３の粒子セットの全ての粒子が前記アレイにロードされた後にのみ生じる、請求項１５に記載の方法。
18. 第３の粒子セットを第２および第３のキュー内に受け取ることをさらに備え、前記第３の粒子セットは前記現在の期間中に前記アレイにロードされる全ての粒子を含み、
    前記第１の粒子セットの粒子を前記第１のキューから前記アレイへと受け渡すことは、前記第３の粒子セットの全ての粒子が前記アレイにロードされた後にのみ生じる、請求項１５に記載の方法。
19. 第４の粒子セットを第３のキュー内に受け取ることをさらに備え、前記第４の粒子セットは、前記アレイにロードされるとともにストリーミングされ、前記第４の粒子セットの粒子のストリーミングは、前記第４の粒子セットの粒子のロードが完了した後にのみ開始される、請求項１５に記載の方法。
20. 前記第１のキューが論理キューであり、複数の選択された物理キューから前記第１および第２の粒子セットを選択することを含む、請求項１５に記載の方法。