力の計算プログラムの記述

Stilliger-Weberポテンシャル

$$ V = \varepsilon \sum_{i}\sum_{j>i} f_2(r_{ij}) + \varepsilon \sum_i\sum_{j \neq i}\sum_{\stackrel{k>j}{k \neq i}} h(r_{ij}, r_{ik}, \theta_{jik}) $$

で、原子\(i\)が受ける力は

$$\begin{eqnarray}
\boldsymbol{F}_i & = & – \varepsilon \sum_{j\neq i} \frac{\partial f_2(r_{ij})}{\partial r_{ij}}\frac{ \boldsymbol{r}_i-\boldsymbol{r}_j }{r_{ij}}\\
&&- \varepsilon \sum_{j \neq i} \sum_{k > j} \left ( A_{jik} \frac{\boldsymbol{r}_i – \boldsymbol{r}_j}{r_{ij}} + B_{jik} \frac{\boldsymbol{r}_i – \boldsymbol{r}_k}{r_{ik}} \right ) \\
& & – \varepsilon \sum_{j \neq i}\sum_{k \neq i} \left ( C_{ijk} \frac{\boldsymbol{r}_i – \boldsymbol{r}_j}{r_{ij}} + D_{ijk} \frac{\boldsymbol{r}_i – \boldsymbol{r}_k}{r_{ik}}\right ) \\
\end{eqnarray} $$

と書けることがわかった。問題はこれをどうプログラムコードにするかである。

力の計算ルーチンのターム指向記述で解説したように、ターム指向の記述に従って書いていけばそれほど難しくない。

1. ポテンシャルエネルギーを計算するコードを記述する。

$$ \begin{eqnarray}
&& \bigl \{ V \;+\!\!=\; \varepsilon \sum_{i}\sum_{j>i} f_2(r_{ij}) + \varepsilon \sum_i\sum_{j \neq i}\sum_{\stackrel{k>j}{k \neq i}} h(r_{ij}, r_{ik}, \theta_{jik}) \bigr \}\\
&=& \mathop{\large \Lambda}_{i} \mathop{\large \Lambda}_{j>i} \bigl\{ V \;+\!\!=\; \varepsilon f_2(r_{ij})\bigr \} + \mathop{\large \Lambda}_{i} \mathop{\large \Lambda}_{j \neq i} \mathop{\large \Lambda}_{\stackrel{k>j}{k \neq i}} \bigl\{ V \;+\!\!=\; \varepsilon h(r_{ij}, r_{ik}, \theta_{jik})\bigr \}
\end{eqnarray}$$

\(\mathop{\large \Lambda}\)はC言語のforループを表す本稿独自の記法である。詳しくは力の計算ルーチンのターム指向記述を参照のこと。