いくつかの実装の詳細¶

この章では、いくつかの実装の詳細について説明します。これはまだ完全ではありませんが、そうでない場合、理解するのが難しいと思われるいくつかの点を明確にするために、この章を作成しました。

単一のビリアル・サマリ (Gromacs)¶

virial \(\Xi\) は、完全なテンソル形式で次のように表すことができます:

(13)¶\[\Xi~=~-\frac{1}{2}~\sum_{i < j}^N~\mathbf{r}_{ij}\otimes\mathbf{F}_{ij}\]

\(\otimes\) は、2 つのベクトルの直積を表します。 [1] MD プログラムの内部ループでこれを計算する場合、9 回の乗算と 9 回の加算が必要です。 [2]

ここに、内ループ 177 から、ビオールの計算を抽出する方法を示す。

バイラル¶

周期境界条件を持つシステムでは、バイラルを計算する際には、周期性を考慮する必要があります。

(14)¶\[\Xi~=~-\frac{1}{2}~\sum_{i < j}^{N}~\mathbf{r}_{ij}^n\otimes\mathbf{F}_{ij}\]

ここで、\(\mathbf{r}_{ij}^n\) は、原子 \(i\) から原子 \(j\) への 最も近いイメージ の距離ベクトルを表します。この定義では、原子 \(i\) の位置ベクトル \(\mathbf{r}_i\) に シフトベクトル \(\delta_i\) を加えます。したがって、差ベクトル \(\mathbf{r}_{ij}^n\) は次のようになります。

(15)¶\[\mathbf{r}_{ij}^n~=~\mathbf{r}_i+\delta_i-\mathbf{r}_j\]

または、短縮形として:

(16)¶\[\mathbf{r}_{ij}^n~=~\mathbf{r}_i^n-\mathbf{r}_j\]

トリキリン系のシステムでは、\(i\) の 27 種類の可能性があります。truncated octahedron を使用する場合、15 種類の可能性があります。

非相互作用力からのバイラル法¶

ここでは、*非相互作用力*ルーチンにおける単一の合計バイラルを導出する方法を示します。|Gromacs|固有のいくつかの考慮事項を以下に示します。

短距離相互作用を計算する際には、*最小画像規約*を適用し、各近傍の最も近い画像のみを考慮します。特に、粒子とそれ自身の周期的な画像との間の相互作用は決して考慮されません。以下に示すすべての式において、これは i ≡ j を意味します。
一般的に、\(i\) または \(j\) の粒子を隣接するセルに移動させることで、最も近い相互作用を得ることができます（移動：math:delta_{ij}）。ただし、最小画像規約では、中心セルの粒子に対して最大27種類の異なるシフトが存在し、典型的な（非常に短距離の）生体分子相互作用の場合、各粒子に対して通常、ほんのわずかな異なるシフトのみが関与します。さらに、各相互作用は、1つのシフトでのみ存在することが可能です。
GROMACS の非相互作用の場合、この方法を用いて、各粒子の近傍リストを複数の独立したリストに分割します。各リストには、粒子の \(i\) に対して一定のオフセット \(\delta_i\) が設定されます。これをシフトの和として表現します（ここでインデックス \(s\) を使用します）。この場合、各粒子相互作用は、この和の項のいずれか1つにのみ寄与し、オフセットは粒子 \(j\) に依存しなくなります。もし和に複雑な依存関係 \(s\) が含まれていない場合、この和は単純に \(i\) と/または \(j\) の和に還元されます。
いくつかの計算を簡略化するために、\(j<i\) のような和を、すべての粒子に対する二重の和に置き換え、2で割ります（ただし、\(i \ eq j\)）。

上記のバイラル定義から、次に得られる

(17)¶\[\begin{split}\begin{aligned} \Xi &~=~&-{\frac{1}{2}}~\sum_{i < j}^{N}~{\mathbf r}^n_{ij} \otimes {\mathbf F}_{ij} \nonumber \\ &~=~&-{\frac{1}{2}}~\sum_{i < j}^{N}~\left( {\mathbf r}_i + \delta_{ij} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{j=1}^{N}~\left( {\mathbf r}_i + \delta_{ij} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} - {\mathbf r}_j \right) \otimes {\mathbf F}_{ij,s} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{i=}^{N}~\sum_{s}~\sum_{j=1}^{N}~\left( \left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} -{\mathbf r}_j \otimes {\mathbf F}_{ij,s} \right) \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^N ~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} + {\frac{1}{4}}\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^{N} {\mathbf r}_j \otimes {\mathbf F}_{ij,s} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{i=1}^{N}~\sum_{s}~\sum_{j=1}^N ~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{ij,s} + {\frac{1}{4}}\sum_{i=1}^{N}~\sum_{j=1}^{N} {\mathbf r}_j \otimes {\mathbf F}_{ij} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes ~\sum_{j=1}^N {\mathbf F}_{ij,s} + {\frac{1}{4}}\sum_{j=1}^N {\mathbf r}_j \otimes \sum_{i=1}^{N} {\mathbf F}_{ij} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes ~\sum_{j=1}^N {\mathbf F}_{ij,s} - {\frac{1}{4}}\sum_{j=1}^N {\mathbf r}_j \otimes \sum_{i=1}^{N} {\mathbf F}_{ji} \nonumber \\ &~=~&-{\frac{1}{4}}~\sum_{s}~\sum_{i=1}^{N}~\left( {\mathbf r}_i + \delta_{i,s} \right) \otimes {\mathbf F}_{i,s} - {\frac{1}{4}}\sum_{j=1}^N~{\mathbf r}_j \otimes {\mathbf F}_{j} \nonumber \\ &~=~&-{\frac{1}{4}}~\left(\sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} + \sum_{j=1}^N~{\mathbf r}_j \otimes {\mathbf F}_{j} \right) - {\frac{1}{4}}\sum_{s}~\sum_{i=1}^{N} \delta_{i,s} \otimes {\mathbf F}_{i,s} \nonumber \\ &~=~&-{\frac{1}{2}}\sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} -{\frac{1}{4}}\sum_{s}~\sum_{i=1}^{N}~\delta_{i,s} \otimes {\mathbf F}_{i,s} \nonumber \\ &~=~&-{\frac{1}{2}}\sum_{i=1}^{N}~{\mathbf r}_i \otimes {\mathbf F}_{i} -{\frac{1}{4}}\sum_{s}~\delta_{s} \otimes {\mathbf F}_{s} \nonumber \\ &~=~&\Xi_0 + \Xi_1\end{aligned}\end{split}\]

第2の最終段階では、各シフトベクトルが粒子の座標 \(i\) に依存しないという性質を利用しました。したがって、各シフトベクトルに対応するすべての力を合計し、その後、カーネルの外で異なるシフトベクトルの合計を使用することが可能です。また、以下のものを利用しました。

(18)¶\[\begin{split}\begin{aligned} \mathbf{F}_i&~=~&\sum_{j=1}^N~\mathbf{F}_{ij} \\ \mathbf{F}_j&~=~&\sum_{i=1}^N~\mathbf{F}_{ji} \end{aligned}\end{split}\]

これは、\(i\) と \(j\) の間の総力を表します。ここで、ニュートンの第3法則を使用しています。

(19)¶\[\mathbf{F}_{ij}~=~-\mathbf{F}_{ji}\]

実装においては、シフトを含む項を2倍にする必要があります：delta_i。同様に、いくつかの箇所では、シフト依存の力をすべてのシフトに対して合計して、相互作用または粒子ごとの総力を算出しています。

これにより、総バイラル：math:Xi は、粒子に関する単一の和である成分：math:Xi_0 と、粒子シフトの影響を表す成分：math:Xi_1 に分割されます。：math:Xi_1 は、異なるシフトベクトルの和のみです。

分子内シフト (mol-shift)¶

「結合力とSHAKEを使用する場合、周期性を保存するための*mol-shift*リストを作成できます。これは、配列`mshift`に、各原子に対応する`shiftvec`配列内のインデックスを格納することで実現します。」

そのようなリストを生成するためのアルゴリズムは、分子内の各粒子をグラフ内のビーズとして、結合をエッジとして扱うことで、グラフ理論から導き出すことができます。

双方向グラフとして、結合と原子を表現する
すべての原子を白にする
白い原子を黒色に変え（原子：math:i）、中央の箱に入れます。
i の周囲の、現在白またはグレーになっているすべての隣接要素を、
選択肢として、グレーの原子（原子：math:j）を選び、その黒い隣接原子のいずれかとの適切な周期性を設定し、それを黒くする。
j の周囲のすべての、現在白またはグレーになっている隣接要素をすべて変更する。
もし灰色のアトムが残っている場合は、[5]に進んでください。
もし白い原子が残っている場合は、[3]に進んでください。

このアルゴリズムを使用することで、

結合力の計算とSHAKEの最適化を同時に行う
単一和法で、結合力を利用してバイラルを再計算する

双方向グラフとして、結合の表現を見つけ出す。

コバレントな結合からのバイラル¶

コバレント結合の力がビラールに寄与するため、以下の式が成り立ちます。

(20)¶\[\begin{split}\begin{aligned} b &~=~& \|\mathbf{r}_{ij}^n\| \\ V_b &~=~& \frac{1}{2} k_b(b-b_0)^2 \\ \mathbf{F}_i &~=~& -\nabla V_b \\ &~=~& k_b(b-b_0)\frac{\mathbf{r}_{ij}^n}{b} \\ \mathbf{F}_j &~=~& -\mathbf{F}_i\end{aligned}\end{split}\]

したがって、結合からのバイラル寄与は次のようになります:

(21)¶\[\begin{split}\begin{aligned} \Xi_b &~=~& -\frac{1}{2}(\mathbf{r}_i^n\otimes\mathbf{F}_i~+~\mathbf{r}_j\otimes\mathbf{F}_j) \\ &~=~& -\frac{1}{2}\mathbf{r}_{ij}^n\otimes\mathbf{F}_i\end{aligned}\end{split}\]

Virial (SHAKE から)¶

重要な貢献の一つは、SHAKEアルゴリズムにおける「シャック」です。粒子に作用する力を「シャック」することで、制約を満たすことができます。この力がアルゴリズムから得られない場合（標準のSHAKEの場合など）、後で（*leap-frog*を使用する場合）次のように計算できます。

(22)¶\[\begin{split}\begin{aligned} \Delta\mathbf{r}_i&~=~&{\mathbf{r}_i}(t+{\Delta t})- [\mathbf{r}_i(t)+{\bf v}_i(t-\frac{{\Delta t}}{2}){\Delta t}+\frac{\mathbf{F}_i}{m_i}{\Delta t}^2] \\ {\bf G}_i&~=~&\frac{m_i{\Delta}{\mathbf{r}_i}}{{\Delta t}^2i}\end{aligned}\end{split}\]

これは一般的なケースでは役に立ちません。周期性が不要な場合にのみ（例えば、硬水の場合）使用できます。それ以外の場合は、SHAKEの内部ループにバイラル計算を追加する必要があります。

「適用可能な場合は、バイラルは以下の式で計算できます:」

(23)¶\[\Xi~=~-\frac{1}{2}\sum_i^{N_c}~\mathbf{r}_i\otimes\mathbf{F}_i\]

ここで、N_c は制約された原子の数を表します。

最適化¶

ここでは、並列処理以外の|Gromacs|で使用されているいくつかのアルゴリズム最適化について説明します。

内部ループ（水）¶

GROMACS は、水分子と他の原子との間の非相互作用的な相互作用を計算するために、特別な内部ループを使用し、さらに水分子同士の相互作用を計算するための別のループセットを使用します。このループは、2種類の水モデルに対して高度に最適化されています。3サイトモデル（SPC 80 に類似したもの）の場合、たとえば、次のように記述できます。

この分子には3つの原子が含まれています。
分子全体が単一の電荷グループとして機能します。
最初の原子には、レンナー・ジョーンズ（第 lj 節を参照）とクーロン（第 coul 節を参照）の相互作用があります。
原子2と3は、クーロン相互作用のみを持ち、電荷は等しい。

これらのループも、SPC/E：<refBerendsen87> と TIP3P：<refJorgensen83> のような水モデルにも使用できます。また、TIP4P：<refJorgensen83> に類似した4サイトの水モデルにも使用できます。

この分子には4つの原子が含まれています。
分子全体が単一の電荷グループとして機能します。
最初の原子には、Lennard-Jones相互作用（第 lj 節を参照）のみが存在します。
原子2と3は、クーロン力（第 coul 参照）のみの影響を受け、電荷は等しい。
アトム4は、コクーロン相互作用のみを持ちます。

これらの実装の利点は、単一のループ内でより多くの浮動小数点演算を実行できることです。これにより、一部のコンパイラはコードの実行順序をより効率的に設定できます。しかし、最先端のコンパイラでも、実行順序の最適化には問題があることが判明しました。そのため、最適なパフォーマンスを得るには、手動での調整が必要となります。これには、共通部分式の削除や、コードの配置変更などが含まれます。