力の計算

今回の問題では、以下のレナード・ジョーンズ型ポテンシャルが粒子間に作用すると仮定しています。 粒子i, j間に作用する力 $(F_x^{ij}, F_y^{ij})$は$\nabla U(r)$なので

$$F_x^{ij} = \frac{\partial}{\partial x}U(r) = 12(x_i - x_j)(\frac{1}{r_{ij}^{14}}-\frac{1}{r_{ij}^{8}})$$ (4)

$$F_y^{ij} = \frac{\partial}{\partial y}U(r) = 12(y_i - y_j)(\frac{1}{r_{ij}^{14}}-\frac{1}{r_{ij}^{8}})$$ (5)

となります。なお、このときポテンシャル中の$\epsilon$と$\sigma$は1に規格化されています。この式を用いて以下の手順で各粒子が受ける力を計算します。

1. 配列force[200][2]を用意する。型に注意。また、配列の長さは粒子の数と同じで次元は2（x, y方向の力を格納）である。
2. 配列を0で初期化する。
3. 以下のプロセスをペアリストの数だけ繰り返す。
4. ペアリストから粒子番号(i, j)を抽出する。
5. 力の作用する範囲 $R_{cutoff} < R_{max}$にある場合だけ力を計算する。例えば、$R_{cutoff}$ = 2.5。$R_{max}$内の粒子全てを考慮しない理由を考えてみよ。
6. 一方の粒子iに作用する力$F^{ij}$は座標値の差$\Delta x$と距離$r$から計算できる。
7. 他方の粒子jには反作用の力$-F^{ij}$が作用する。
8. それぞれの力を該当する粒子の配列forceに足していく。例えばx方向であればforce[i][0]とforce[j][0]に、y方向であればforce[i][1]とforce[j][1]に加えていく。

ペアリストの数だけ力を足していくプロセスが完了すれば、各粒子に働く力が配列forceに格納されている。以上の内容も関数として実装します。