逆誤差伝播法による微分のコーディング・ノート#3

Tersoffポテンシャル実装のテスト

#2に続いて，3体間ポテンシャルについて逆誤差伝播法（back propagation）を応用した微分のコーディングを紹介する．

3体間ポテンシャルの中でも複雑なTersoffポテンシャルに対して，Pythonでの実装を行う．

Tersoffポテンシャル関数

系の全ポテンシャルエネルギー

E=\frac{1}{2}\sum_{j\ne i}\sum_i f_c(r_{ij})\left\{ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) \right\}

カットオフ関数

f_c(r) = \left\{ 
    \begin{array}{ll}
    1 & (r<R-D>) \\
    \frac{1}{2}-\frac{1}{2}\sin\left\{\frac{\pi}{2}\frac{r-R}{D} \right\} & (R-D\le r \le R+D)\\
    0 & (r>R+D)
    \end{array}
\right.

反発項（二体間）

f_R(r_{ij}) = A_{ij}\exp(-\lambda_1 r_{ij})

引力項（二体間）

f_A(r_{ij}) = -B_{ij}\exp(-\lambda_2 r_{ij})

ボンドオーダー項（三体間）

b_{ij} = (1+\beta_i^n \zeta_{ij}^n)^{-1/(2n)}

\zeta_{ij} = \sum_{k\ne i,j} \zeta_{jik}

\zeta_{jik} = f_c(r_{ik})g(\theta_{jik})\exp\left\{\lambda_3^3(r_{ij}-r_{ik})^3\right\}

（ここで$\theta_{jik}$は角$jik$，すなわちi-jボンドとi-kボンドがなす角とする）

ボンド角依存項

g(\theta) = 1 + \frac{c^2}{d^2} - \frac{c^2}{d^2+(h-\cos\theta)^2}

なお，原子$i$あたりのエネルギー$E_i$を

E_i = \frac{1}{2}\sum_{j\ne i}f_c(r_{ij})\{f_R(r_{ij})+b_{ij}f_A(r_{ij})\}

と定義しておく($E=\sum_i E_i$)．

計算の流れ

エネルギーを計算するためには

1. $r_{ij} = |\mathbf{r}_i-\mathbf{r}_j|$
2. $f_c^{ij} := f_c(r_{ij})$, $f_A^{ij} := f_A(r_{ij})$, $f_R^{ij} := f_R(r_{ij})$
3. $r_{ik} = |\mathbf{r}_i-\mathbf{r}_k|$
4. $r_{jk} = |\mathbf{r}_k-\mathbf{r}_j|$
5. $\cos\theta_{jik} = (r_{ij}^2+r_{ik}^2-r_{jk}^2)/(2r_{ij}r_{ik})$
6. $g(\theta_{jik})$
7. $\zeta_{jik}$
8. $\zeta_{ji} = \sum_{k}\zeta_{jik}$
9. $b_{ij}$
10. $E_i = \frac{1}{2}\sum_{j\ne i}f_c(r_{ij}){f_R(r_{ij})+b_{ij}f_A(r_{ij})}$
11. $E=\sum_i E_i$

の順に計算する．1~8が$i,j$の二重ループで，その中に$k$のループ（3~6）が入る形になる．

各変数での微分

#1, #2でやったように，各変数に対する$E$の微分を$^*$をつけて表すことにする．

E_i^* = \frac{\partial E}{\partial E_i} = 1

e_{ij}^* = \frac{\partial E}{\partial e_{ij}} = \sum_k E_k^* \frac{\partial E_k}{\partial e_{ij}} = \frac{1}{2} E_i^* = \frac{1}{2}

b_{ij}^* = \frac{\partial E}{\partial b_{ij}} = e_{ij}^*\frac{\partial e_{ij}}{\partial b_{ij}} = e_{ij}^* f_c(r_{ij})f_A(r_{ij})

\zeta_{ij}^* = b_{ij}^*\frac{\partial b_{ij}}{\partial\zeta_{ij}}=b_{ij}^*\left(-\frac{1}{2n}(1+\beta_i^n\zeta_{ij}^n)^{-\frac{1}{2n}-1}\right)\cdot n\beta_i^n \zeta_{ij}^{n-1}

\zeta_{jik}^* = \zeta_{ij}^* \frac{\partial\zeta_{ij}}{\partial\zeta_{jik}} = \zeta_{ij}^*

g_{jik}^* = \zeta_{jik}^* \frac{\partial\zeta_{jik}}{\partial g_{jik}}=\zeta_{jik}^* f_c(r_{ik})\exp\{\lambda_3^3(r_{ij}-r_{ik})^3\}

$r_{ij}^*$ の項
ここからがやや難しい（注意しないとハマる）かもしれない．

r_{ij}^* = {f_c^{ij}}^*\frac{\partial f_c^{ij}}{\partial r_{ij}}
 + {f_R^{ij}}^*\frac{\partial f_R^{ij}}{\partial r_{ij}}
 + {f_A^{ij}}^*\frac{\partial f_A^{ij}}{\partial r_{ij}}
 + b_{ij}^*\frac{\partial b_{ij}}{\partial r_{ij}}
 + \sum_{k\ne i,j}b_{ik}^*\frac{\partial b_{ik}}{\partial r_{ij}}

特に，右辺第4項・第5項は見落としがちかもしれない．これらの項は，$\zeta_{jik}$の関数形が$r_{ij}, r_{ik}$を変数としていることにより現れる．なお，$\zeta_{jik}$は$\theta_{jik}$の関数でもあるが，その寄与は後で出てくる．

蛇足だが，ここで$\frac{\partial\theta_{jik}}{\partial r_{ij}}$を考えてはいけない．（え，そんなの当り前だろうって？失礼しました．）

ここで，右辺第4項・第5項をもう少し進めておく．第4項は

b_{ij}^*\frac{\partial b_{ij}}{\partial r_{ij}} 
= \sum_{k\ne i,j}\zeta_{jik}^*\frac{\partial\zeta_{jik}}{\partial r_{ij}}
= \zeta_{ij}^*\sum_{k\ne i,j}\frac{\partial\zeta_{jik}}{\partial r_{ij}}

となる．ここで

\frac{\partial\zeta_{jik}}{\partial r_{ij}}=3\lambda_3^3(r_{ij}-r_{ik})^2\zeta_{jik}

である．

第5項は

\sum_{k\ne i,j}b_{ik}^*\frac{\partial b_{ik}}{\partial r_{ij}} 
=\sum_{k\ne i,j}b_{ik}^*\frac{\partial b_{ik}}{\partial\zeta_{ik}}\frac{\partial\zeta_{ik}}{\partial r_{ij}} 
= \sum_{k\ne i,j}\zeta_{ik}^*\frac{\partial\zeta_{ik}}{\partial r_{ij}}

とできるが，$\zeta_{ik}=\sum_{l\ne i,k}\zeta_{kil}$なので$l=j$のときのみ偏微分がnonzeroとなるため，結局

\sum_{k\ne i,j}\zeta_{ik}^*\frac{\partial\zeta_{kij}}{\partial r_{ij}}

になるが，これをこのまま用いて $r_{ij}^*$ へ $\sum_k b_{ik}^*\frac{\partial b_{ik}}{\partial r_{ij}}=\sum_k \zeta_{ik}^*\frac{\partial\zeta_{kij}}{\partial r_{ij}}$ の項を足しこむ，とするのは非効率である．なぜなら，kについてのループを書き，その中で $\zeta_{ik}^*\frac{\partial\zeta_{kij}}{\partial r_{ij}}$ を求める操作は，$\zeta_{ik}^*$ が k ごとに変わるため，いちいち $\zeta_{ik}^*$ を計算しなおさなくてはならないことになる．ではどうするか？

原子 i の neighbor について j のループをとって $r_{ij}^*$ を全て求めるのであるから，k と j のループを入れ替えても問題はない．これを，k が j の内側ループになるように書き換えると，$r_{ik}^*$ へ $\zeta_{ij}^*\frac{\partial\zeta_{jik}}{\partial r_{ik}}$ を足しこむ，という操作

r_{ik}^* \leftarrow \zeta_{ij}^*\frac{\partial\zeta_{jik}}{\partial r_{ik}}

をまず k についてループし，kループの外側で j についてループさせればよい．このようにすれば，kループの中では $\zeta_{ij}^*$ は不変なので，効率的に計算できる．

ここで，

\frac{\partial\zeta_{kij}}{\partial r_{ij}}=-3\lambda_3^3(r_{ij}-r_{ik})^2\zeta_{jik}+\frac{f_c'(r_{ik})}{f_c(r_{ik})}\zeta_{jik}

である．

$(\cos\theta_{jik})^*$ の項

(\cos\theta_{jik})^* 
= g_{jik}^* \frac{\partial g_{jik}}{\partial(\cos\theta_{jik})}
= g_{jik}^* \frac{-2c^2(h-\cos\theta_{jik})}{\{d^2+(h-\cos\theta_{jik})^2\}^2}

$\mathbf{R}_i^*$ の項

\mathbf{R}_l^* = \frac{\partial E}{\partial\mathbf{R}_l}
=\sum_{i,j}r_{ij}^*\frac{\partial r_{ij}}{\partial\mathbf{R}_l}
+\sum_{i,k}r_{ik}^*\frac{\partial r_{ik}}{\partial\mathbf{R}_l}
+\sum_{i,j,k}\cos\theta_{jik}^*\frac{\partial\cos\theta_{jik}}{\partial\mathbf{R}_l}\\
+\sum_{i,j,k}\cos\theta_{ijk}^*\frac{\partial\cos\theta_{ijk}}{\partial\mathbf{R}_l}
+\sum_{i,j,k}\cos\theta_{jki}^*\frac{\partial\cos\theta_{jki}}{\partial\mathbf{R}_l}

であるから

\mathbf{R}_i^* 
=\sum_{j\ne i}r_{ij}^*\frac{\partial r_{ij}}{\partial\mathbf{R}_i}
+\sum_{j\ne i}r_{ji}^*\frac{\partial r_{ji}}{\partial\mathbf{R}_i}\\
+\sum_{j\ne i}\sum_{k\ne i,j}\left(
\cos\theta_{jik}^*\frac{\partial\cos\theta_{jik}}{\partial\mathbf{R}_i}
+\cos\theta_{ijk}^*\frac{\partial\cos\theta_{ijk}}{\partial\mathbf{R}_i}
+\cos\theta_{jki}^*\frac{\partial\cos\theta_{jki}}{\partial\mathbf{R}_i}
\right)

その他関数の微分

\frac{\partial f_c(r)}{\partial r} = \left\{
   \begin{array}{ll}
   -\frac{\pi}{4D}\cos\{\frac{\pi}{2}\frac{r-R}{D}\} & (R-D\le r\le R+D)\\
   0 & (r<R-D, R+D<r)
   \end{array} 
\right.

\frac{\partial f_R(r)}{\partial r} = -\lambda_1 f_R(r)

\frac{\partial f_A(r)}{\partial r} = -\lambda_2 f_A(r)

Pythonでの実装

基本事項

#1, #2を踏まえ，関数作成のポイントを整理しておく．

1. $\mathbf{x}$より$\mathbf{X}$を求める（$\mathbf{X}$が$\mathbf{x}$の関数である）とき，これをBeadXのforward処理として実装する．
2. $\mathbf{x}$，$\mathbf{X}$の要素数がそれぞれ$m$個，$n$個ならば，a=BeadX(m,n)としてインスタンスを発行する．
3. この関数のbackward処理は，$\mathbf{X}^*$（$\mathbf{X}^*=\frac{\partial E}{\partial\mathbf{X}}$なので$n$次）を入力として$\mathbf{x}^*$（$m$次）を出力する．

用意する関数

• BeadR(6,1): (forward処理で) $\mathbf{r}_i, \mathbf{r}_j \to r_{ij}$
• を計算する（同様に，$r_{ik}, r_{jk}$の計算にも使う）．rij=BeadR(6,1)で発行，ri = [$r_i^x, r_i^y, r_i^z$], rj = [$r_j^x, r_j^y, r_j^z$] などとしておき，rij.forward(ri,rj)のように使う．
• BeadCT(12,1): $\mathbf{r}_i, \mathbf{r}_j, \mathbf{r}_k, r_{ij}, r_{ik}, r_{jk} \to \cos\theta_{jik}$．forwardの入力としては$\mathbf{r}_i, \mathbf{r}_j, \mathbf{r}_k$で十分なのであるが，この関数を使う前に$r_{ij}$等はrij関数などにより求められているため，同じ計算をBeadCTの中で行うのではなく，入力値として要求する設計とする．backward処理は$\mathbf{R}_i^*, \mathbf{R}_j^*, \mathbf{R}_k^*$を出力（$r_{ij}^*$等は計算しない）．
• BeadFC(1,1): $r\to f_c(r)$
• BeadFR(1,1): $r\to f_R(r)$
• BeadFA(1,1): $r\to f_A(r)$
• BeadG(1,1): $\cos\theta_{jik}\to g_{jik}$
• BeadZ(3,1): $g_{jik},r_{ij},r_{ik}\to\zeta_{jik}$
• BeadZIJ(1,1): $\zeta_{jik}\to\zeta_{ij}$ ($\sum_k$をとる)
• BeadBIJ(1,1): $\zeta_{ij}\to b_{ij}$
• BeadEIJ(4,1): $f_c^{ij},f_R^{ij},b_{ij},f_A^{ij}\to e_{ij}$
• BeadEI(1,1): $e_{ij}\to E_i$（$\sum_j$をとる）
• BeadE(1,1): $E_i\to E$（$\sum_i$をとる）

Tersoffパラメータとカットオフ関数

Phys. Rev. B 38 (1988) 9902でSi(B)として紹介されているパラメータを採用する．また，カットオフ関数$f_c(r)$とその微分を定義しておく．

import numpy as np

Ter_R = 3.0
Ter_D = 0.2
Ter_A = 3.2647e3
Ter_B = 9.5373e2
Ter_lam1 = 3.2394
Ter_lam2 = 1.3258
Ter_lam3 = 1.3258
Ter_bet = 3.3675e-1
Ter_n = 2.2956e1
Ter_c = 4.8381
Ter_d = 2.0417
Ter_h = 0.0000

def fcut(r): # cutoff
    if (r < Ter_R-Ter_D):
        f = 1
    elif (r > Ter_R+Ter_D):
        f = 0
    else:
        f = 0.5 - 0.5 * np.sin(np.pi/2.0 * (r-Ter_R)/Ter_D)
    return f

def dfcut(r): # derivative of fc
    if (r < Ter_R-Ter_D):
        df = 0
    elif (r > Ter_R+Ter_D):
        df = 0
    else:
        df = -np.pi/4.0/Ter_D * np.cos(np.pi/2.0 *(r-Ter_R)/Ter_D)
    return df

親クラスの準備

コンストラクタとリセット機能を書いておく．ninsがforwad処理の入力変数の次元，noutが出力変数のそれを表す．従って，backward処理についてはこれらが逆転する．

class Bead(object):
    def __init__(self, nins, nout):
        self.nins = nins
        self.nout = nout
        if (nins == 1):
            self.ar_in = 0
            self.ar_gout = 0
        else:
            self.ar_gout = np.zeros(nins)
            self.ar_in = np.zeros(nins)
        if (nout == 1):
            self.ar_out = 0
            self.ar_gin = 0
        else:
            self.ar_out = np.zeros(nout)
            self.ar_gin = np.zeros(nout)
    def reset(self):
        if (self.nout == 1):
            self.ar_out = 0
        else:
            self.ar_out[:] = 0
    def resetg(self):
        if (self.nins == 1):
            self.ar_gout = 0
        else:
            self.ar_gout[:] = 0

BeadR関数の実装

BeadRのforward, backward処理を実装する．forwardのインプットとして${\mathbf{r}_i, \mathbf{r}_j}={r_i^x,r_i^y,r_i^z,r_j^x,r_j^y,r_j^z}$の６次元ベクトルをとる．

アウトプットは$r_{ij}$とスカラーなので１次元となる．

r_{ij} = |\mathbf{r}_i-\mathbf{r}_j|

\mathbf{R}_l^* = \sum_i \sum_{j\ne i} r_{ij}^*\frac{\partial r_{ij}}{\partial \mathbf{R}_i}

であったから，

\mathbf{R}_i^* = \sum_k \sum_{l\ne k} r_{kl}^*\frac{\partial r_{kl}}{\partial \mathbf{R}_i}
=\sum_{j\ne i}r_{ij}^*\frac{\partial r_{ij}}{\partial \mathbf{R}_i}
+\sum_{j\ne i}r_{ji}^*\frac{\partial r_{ji}}{\partial \mathbf{R}_i}

各原子の$\mathbf{R}^*$を求めるには$i$と$j$の二重ループを回すことになるが，その際に

r_{ij}^*\frac{\partial r_{ij}}{\partial \mathbf{R}_i} \to \mathbf{R}_i^*\quad に

r_{ij}^*\frac{\partial r_{ij}}{\partial \mathbf{R}_j} \to \mathbf{R}_j^*\quad に

加算していけばよい．したがって

\frac{\partial r_{ij}}{\partial \mathbf{R}_l} = 
\begin{cases}
(\mathbf{r}_i-\mathbf{r}_j)/r_{ij} \quad (l = i)\\
(\mathbf{r}_j-\mathbf{r}_i)/r_{ij} \quad (l = j)\\
0 \quad (l\ne i, l\ne j)
\end{cases}

を用いると，以下のように書ける．

class BeadR(Bead):
    def forward(self, arri_in, arrj_in): # takes [r_i(x,y,z),r_j(x,y,z)] (6dim.) as input
        for k in range(3):
            self.ar_in[k  ] = arri_in[k]
            self.ar_in[k+3] = arrj_in[k]
        dx = arrj_in[0] - arri_in[0]
        dy = arrj_in[1] - arri_in[1]
        dz = arrj_in[2] - arri_in[2]        
        self.ar_out = np.sqrt(dx*dx + dy*dy + dz*dz)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        xx = ar_gin / self.ar_out 
        self.ar_gout[0] = xx * (self.ar_in[0] - self.ar_in[3]) 
        self.ar_gout[1] = xx * (self.ar_in[1] - self.ar_in[4]) 
        self.ar_gout[2] = xx * (self.ar_in[2] - self.ar_in[5]) 
        self.ar_gout[3] = -self.ar_gout[0]
        self.ar_gout[4] = -self.ar_gout[1]
        self.ar_gout[5] = -self.ar_gout[2]
        return self.ar_gout

なお，ここで得られる$\mathbf{R}^*$はボンド長変化に因るものであり，BeadCTのbackwardで得られるボンド角依存項と足し合わせる必要がある．

BeadCT関数の実装

forward処理では，$\mathbf{r}_i, \mathbf{r}_j, \mathbf{r}_k, r_{ij}, r_{ik}, r_{jk}$をinputとし，$\cos\theta_{jik}$をoutputとして生成する．

$r_{ij},r_{ik},r_{jk}$は$\mathbf{r}_i, \mathbf{r}_j, \mathbf{r}_k$から計算できるため，これらを全てinputとして要求するのは冗長ではある．だが，この関数を使う前に$r_{ij}$等はrij関数などにより求められているため，同じ計算をBeadCTの中で行うのではなく，入力値として要求する設計とする．

なお，backward処理は$\mathbf{R}_i^*, \mathbf{R}_j^*, \mathbf{R}_k^*$を出力（$r_{ij}^*$等は計算しない）する．くどいようだが，ここで得られるのは$\mathbf{R}^*$のうち角度依存項によるものであるので，BeadR関数のbackward処理で得られる$\mathbf{R}^*$と足し合わせる必要がある．

forward処理は上述の $\cos\theta_{jik} = \frac{r_{ij}^2+r_{ik}^2-r_{jk}}{2r_{ij}r_{ik}}$ で得られる．backward処理は$\cos\theta_{jik}^*，\cos\theta_{ijk}^*，\cos\theta_{jki}^*$と

\frac{\partial\cos\theta_{jik}}{\partial\mathbf{R}_i},\quad
\frac{\partial\cos\theta_{ijk}}{\partial\mathbf{R}_i},\quad
\frac{\partial\cos\theta_{jki}}{\partial\mathbf{R}_i}

より計算可能であるが，そのようにまともに（？）計算すると，$\mathbf{R}_i^*$を求めるために$\cos\theta_{jik}，\cos\theta_{ijk}，\cos\theta_{jki}$やその微分を全て求めねばならず，効率が悪い（無闇にループを書かなくてはいけなくなる）．
そうではなく，$\cos\theta_{jik}$（やその微分）が求められたときに，その寄与を$\mathbf{R}_i^*,\mathbf{R}_j^*,\mathbf{R}_k^*$に足しこむ，とした方がループ構造がシンプルにでき，効率的である．すなわち，

\frac{\partial\cos\theta_{jik}}{\partial\mathbf{R}_i}
=\frac{\partial r_{ij}}{\partial\mathbf{R}_i}\frac{\partial\cos\theta_{jik}}{\partial r_{ij}}
+\frac{\partial r_{ik}}{\partial\mathbf{R}_i}\frac{\partial\cos\theta_{jik}}{\partial r_{ik}}

\frac{\partial\cos\theta_{jik}}{\partial\mathbf{R}_j}
=\frac{\partial r_{ij}}{\partial\mathbf{R}_j}\frac{\partial\cos\theta_{jik}}{\partial r_{ij}}
+\frac{\partial r_{jk}}{\partial\mathbf{R}_j}\frac{\partial\cos\theta_{jik}}{\partial r_{jk}}

\frac{\partial\cos\theta_{jik}}{\partial\mathbf{R}_k}
=\frac{\partial r_{ik}}{\partial\mathbf{R}_k}\frac{\partial\cos\theta_{jik}}{\partial r_{ik}}
+\frac{\partial r_{jk}}{\partial\mathbf{R}_k}\frac{\partial\cos\theta_{jik}}{\partial r_{jk}}

$a:=r_{jk},b:=r_{ij},c:=r_{ik}$とすると

\frac{\partial\cos\theta_{jik}}{\partial r_{ij}}
=\frac{1}{c}-\frac{1}{b}\cos\theta_{jik}

\frac{\partial\cos\theta_{jik}}{\partial r_{ik}}
=\frac{1}{b}-\frac{1}{c}\cos\theta_{jik}

\frac{\partial\cos\theta_{jik}}{\partial r_{jk}}
=-\frac{a}{bc}

以上より，BeadCT関数を次のように実装する．

class BeadCT(Bead):
    def forward(self, ar_in): # takes r_i,r_j,r_k,r_ij,r_ik,r_jk as input
        self.ar_in = ar_in
        b = ar_in[9]; c = ar_in[10]; a = ar_in[11]
        self.ar_out = (b*b+c*c-a*a)/(2.0*b*c)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        b = self.ar_in[9]; c = self.ar_in[10]; a = self.ar_in[11]
        dcosjik_drij = 1.0/c - 1.0/b*self.ar_out
        dcosjik_drik = 1.0/b - 1.0/c*self.ar_out
        dcosjik_drjk = -a/b/c
        drij_dRi = [(x-y)/b for (x,y) in zip(self.ar_in[0:3],self.ar_in[3:6])]
        drik_dRi = [(x-y)/c for (x,y) in zip(self.ar_in[0:3],self.ar_in[6:9])]
        ri_g = [x*dcosjik_drij+y*dcosjik_drik for (x,y) in zip(drij_dRi,drik_dRi)]
        drij_dRj = [(y-x)/b for (x,y) in zip(self.ar_in[0:3],self.ar_in[3:6])]
        drjk_dRj = [(x-y)/a for (x,y) in zip(self.ar_in[3:6],self.ar_in[6:9])]
        rj_g = [x*dcosjik_drij+y*dcosjik_drjk for (x,y) in zip(drij_dRj,drjk_dRj)]
        drik_dRk = [(y-x)/c for (x,y) in zip(self.ar_in[0:3],self.ar_in[6:9])]
        drjk_dRk = [(y-x)/a for (x,y) in zip(self.ar_in[3:6],self.ar_in[6:9])]
        rk_g = [x*dcosjik_drik+y*dcosjik_drjk for (x,y) in zip(drik_dRk,drjk_dRk)]
        xx = ri_g + rj_g + rk_g
        xx.append(0); xx.append(0); xx.append(0)
        self.ar_gout = [n*ar_gin for n in xx]
        return self.ar_gout

BeadFC, BeadFR, BeadFA各関数の実装

これは簡単．

class BeadFC(Bead):
    def forward(self, ar_in): # takes r_ij as input
        self.ar_in = ar_in
        self.ar_out = fcut(ar_in)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gout = ar_gin * dfcut(self.ar_in)
        return self.ar_gout

class BeadFR(Bead):
    def forward(self, ar_in): # takes r_ij as input
        self.ar_in = ar_in
        self.ar_out = Ter_A * np.exp(-Ter_lam1*ar_in)
        return self.ar_out
    def backward(self, ar_gin):
        dfr_drij = -Ter_lam1 * self.ar_out
        self.ar_gout = ar_gin * dfr_drij
        return self.ar_gout

class BeadFA(Bead):
    def forward(self, ar_in): # takes r_ij as input
        self.ar_in = ar_in
        self.ar_out = -Ter_B * np.exp(-Ter_lam2*ar_in)
        return self.ar_out
    def backward(self, ar_gin):
        dfa_drij = -Ter_lam2 * self.ar_out
        self.ar_gout = ar_gin * dfa_drij
        return self.ar_gout

BeadG関数の実装

これも$\cos\theta\to g(\theta)$への変換なので比較的簡単．

class BeadG(Bead):
    def forward(self, ar_in): # takes cosjik as input
        self.ar_in = ar_in
        self.ar_out = 1.0 + (Ter_c/Ter_d)**2 - Ter_c**2 / (Ter_d**2 + (Ter_h - ar_in)**2)
        return self.ar_out
    def backward(self, ar_gin):
        dg_dcos = -2.0 * Ter_c**2 * (Ter_h - self.ar_in) / (Ter_d**2 + (Ter_h - self.ar_in)**2)**2
        self.ar_gout = ar_gin * dg_dcos
        return self.ar_gout

BeadZ関数の実装

forward処理は $g_{jik},r_{ij},r_{ik}\to\zeta_{jik}$ である．backward処理で $g_{jik}^*,r_{ij}^*,r_{ik}^*$ を出力する．

class BeadZ(Bead):
    def forward(self, ar_in): # takes g_jik,r_ij,r_ik as input
        self.ar_in = ar_in
        g = ar_in[0]; rij = ar_in[1]; rik = ar_in[2]
        if (fcut(rik) != 0.0):
            self.ar_out = fcut(rik) * g * np.exp((Ter_lam3 * (rij - rik))**3)
        else:
            self.ar_out = 0.0
        return self.ar_out
    def backward(self, ar_gin):
        rij = self.ar_in[1]; rik = self.ar_in[2]
        if (fcut(rik) != 0.0):
            self.ar_gout[0] = ar_gin * fcut(rik) * np.exp(Ter_lam3**3 * (rij - rik)**3)
        else:
            self.ar_gout[0] = 0.0
        self.ar_gout[1] =  ar_gin * 3.0*Ter_lam3**3*(rij-rik)**2 * self.ar_out
        self.ar_gout[2] = -self.ar_gout[1]
        if (fcut(rik) != 0.0 and dfcut(rik) != 0.0):
            self.ar_gout[2] += ar_gin * (self.ar_out/fcut(rik)*dfcut(rik))
        return self.ar_gout

BeadZIJ関数の実装

k ループを作りsummationを取るだけのものなので簡単．

class BeadZIJ(Bead):
    def forward(self, ar_in): # takes zet_jik as input
        self.ar_in = ar_in
        self.ar_out += self.ar_in
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        self.ar_gout = self.ar_gin
        return self.ar_gout

BeadBIJ関数の実装

$\zeta_{ij}\to b_{ij}$ の関数であり，これも難しくない．

class BeadBIJ(Bead):
    def forward(self, ar_in): # takes z_ij as input
        self.ar_in = ar_in
        self.ar_out = np.power(1.0 + (Ter_bet*ar_in)**Ter_n, -0.5/Ter_n)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        x = (-0.5/Ter_n) * self.ar_out / (1.0 + (Ter_bet*self.ar_in)**Ter_n)
        self.ar_gout = ar_gin * x * (Ter_n * Ter_bet**Ter_n * self.ar_in**(Ter_n-1))
        return self.ar_gout

BeadEIJ関数の実装

forward処理は $f_c,f_R,b_{ij},f_A\to e_{ij}$ であり，backward処理は当然ながら $f_c^*,f_R^*,b_{ij}^*,f_A^*\leftarrow e_{ij}^*$ であることに注意する．

class BeadEIJ(Bead):
    def forward(self, ar_in): # takes fc,fr,bij,fa as input
        self.ar_in = ar_in
        self.ar_out = ar_in[0] * (ar_in[1] + ar_in[2] * ar_in[3])
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        fc_g  = ar_gin * (self.ar_in[1] + self.ar_in[2]*self.ar_in[3])
        fr_g  = ar_gin * self.ar_in[0]
        bij_g = ar_gin * self.ar_in[0] * self.ar_in[3]
        fa_g  = ar_gin * self.ar_in[0] * self.ar_in[2]
        self.ar_gout = [fc_g,fr_g,bij_g,fa_g]
        return self.ar_gout

BeadEI関数の実装

class BeadEI(Bead):
    def forward(self, ar_in): # takes eij as input
        self.ar_in = ar_in
        self.ar_out += 0.5 * self.ar_in
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        dedp = 0.5
        self.ar_gout = ar_gin * dedp
        return self.ar_gout

BeadE関数の実装

class BeadE(Bead): # takes E_i as input
    def forward(self, ar_in):
        self.ar_in = ar_in
        self.ar_out += self.ar_in
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        self.ar_gout = 1.0
        return self.ar_gout

実装した関数を用いてエネルギーと力を計算

所々，backward処理を行う前にforward計算をやり直している部分があることに注意．（もう少し賢いやり方があるかもしれないが・・・）

rij = BeadR(6,1)
rik = BeadR(6,1)
rjk = BeadR(6,1)
cosjik  = BeadCT(12,1)
fc  = BeadFC(1,1)
fr  = BeadFR(1,1)
fa  = BeadFA(1,1)
g   = BeadG(1,1)
zet = BeadZ(3,1)
zij = BeadZIJ(1,1)
bij = BeadBIJ(1,1)
eij = BeadEIJ(4,1)
ei  = BeadEI(1,1)
e   = BeadE(1,1)

e.reset()
for i in range(n):
    e.backward(1.0) #E_i*
    ei.reset()
    for j in range(n):
        if (j != i):
            ri = [rx[i],ry[i],rz[i]]
            rj = [rx[j],ry[j],rz[j]]
            rij.forward(ri,rj)
            fc.forward(rij.ar_out)
            fr.forward(rij.ar_out)
            fa.forward(rij.ar_out)
            zij.reset()
            rij_g = 0.0
            for k in range(n):
                if (k != i and k != j):
                    rk = [rx[k],ry[k],rz[k]]
                    rik.forward(ri,rk)
                    rjk.forward(rj,rk)
                    xx = ri+rj+rk
                    xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                    cosjik.forward(xx)
                    g.forward(cosjik.ar_out)
                    xx = [g.ar_out,rij.ar_out,rik.ar_out]
                    zet.forward(xx)
                    zij.forward(zet.ar_out)
            bij.forward(zij.ar_out)
            xx = [fc.ar_out,fr.ar_out,bij.ar_out,fa.ar_out]
            eij.forward(xx)
            ei.forward(eij.ar_out)
            #
            ei.backward(e.ar_gout) #eij*
            xx = eij.backward(ei.ar_gout) #bij*
            fc_g  = xx[0]; fr_g  = xx[1]; bij_g = xx[2]; fa_g  = xx[3] 
            bij.backward(eij.ar_gout[2]) #zij*
            zjik_g = zij.backward(bij.ar_gout) #zjik*
            rij_g += fc.backward(fc_g)
            rij_g += fr.backward(fr_g)
            rij_g += fa.backward(fa_g)
            for k in range(n):
                if (k != i and k != j):
                    # cosjik to be recalculated
                    rk = [rx[k],ry[k],rz[k]]
                    rik.forward(ri,rk)
                    rjk.forward(rj,rk)
                    xx = ri+rj+rk
                    xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                    cosjik.forward(xx)
                    g.forward(cosjik.ar_out)
                    xx = [g.ar_out,rij.ar_out,rik.ar_out]
                    zet.forward(xx)
                    #
                    xx = zet.backward(zjik_g) #gjik*,rij*,rik*
                    gjik_g = xx[0]; rij_g += xx[1]; rik_g = xx[2] # treatment of rik
                    #redo cosjik.forward..
                    rk = [rx[k],ry[k],rz[k]]
                    rik.forward(ri,rk)
                    rjk.forward(rj,rk)
                    # treatment of rik
                    xxx = rik.backward(rik_g)
                    fx[i] -= xxx[0]; fy[i] -= xxx[1]; fz[i] -= xxx[2]
                    fx[k] -= xxx[3]; fy[k] -= xxx[4]; fz[k] -= xxx[5]
                    #
                    xx = ri+rj+rk
                    xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                    cosjik.forward(xx)
                    #then backward..
                    cosjik_g = g.backward(zet.ar_gout[0]) #cosjik*
                    xx = cosjik.backward(g.ar_gout) #Ri*,Rj*,Rk*,rij*,rik*,rjk*
                    fx[i] -= xx[0]; fy[i] -= xx[1]; fz[i] -= xx[2]
                    fx[j] -= xx[3]; fy[j] -= xx[4]; fz[j] -= xx[5]
                    fx[k] -= xx[6]; fy[k] -= xx[7]; fz[k] -= xx[8]
            xxx = rij.backward(rij_g)
            fx[i] -= xxx[0]; fy[i] -= xxx[1]; fz[i] -= xxx[2]
            fx[j] -= xxx[3]; fy[j] -= xxx[4]; fz[j] -= xxx[5]
    ee = e.forward(ei.ar_out)
print('Etot = ',ee)

作成したコードとテスト

コードはこちら．

bptest04_tersoff.py

'''
Testing back-propagation for derivative calculation
04: Tersoff potential
r_ij := |r_j - r_i|
E_i := 1/2 * sum_{j\ne i} * f_c(r_ij) * { f_R(r_ij) + b_ij f_A(r_ij) }
E := \sum_i E_i
         1 (r<R-D)
f_c(r) = 1/2 - 1/2 * sin{pi/2 * (r-R)/D} (R-D <= r <= R+D)
         0 (r>R+D)
f_R(r) =  A * exp(-\lambda_1 r)
f_A(r) = -B * exp(-\lambda_2 r)
b_ij = (1 + \beta_i^n * \zeta_ij^n)^(-1/2n)
\zeta_ij = \sum_{k\ne i,j} f_c(r_ik) * g(\theta_ijk) * exp{\lambda_3^3(r_ij-r_ik)^3}
g(\theta) = 1 + c^2/d^2 - c^2 / {d^2 + (h-costheta)^2}
'''

import numpy as np

ev = 1.6021892e-19
Ter_R = 3.0
Ter_D = 0.2
Ter_A = 3.2647e3
Ter_B = 9.5373e2
Ter_lam1 = 3.2394
Ter_lam2 = 1.3258
Ter_lam3 = 1.3258
Ter_bet = 3.3675e-1
Ter_n = 2.2956e1
Ter_c = 4.8381
Ter_d = 2.0417
Ter_h = 0.0000

def fcut(r): # cutoff
    if (r < Ter_R-Ter_D):
        f = 1
    elif (r > Ter_R+Ter_D):
        f = 0
    else:
        f = 0.5 - 0.5 * np.sin(np.pi/2.0 * (r-Ter_R)/Ter_D)
    return f

def dfcut(r): # derivative of fc
    if (r < Ter_R-Ter_D):
        df = 0
    elif (r > Ter_R+Ter_D):
        df = 0
    else:
        df = -np.pi/4.0/Ter_D * np.cos(np.pi/2.0 *(r-Ter_R)/Ter_D)
    return df

def fex(a,lam,r) : # A * exp(-\lambda * r)
    f = a * np.exp(-lam*r)
    return f

def dfex(a,lam,r) : # derivative of fex
    df = -lam*a*np.exp(-lam*r)
    return f

class Bead(object):
    def __init__(self, nins, nout):
        self.nins = nins
        self.nout = nout
        if (nins == 1):
            self.ar_in = 0
            self.ar_gout = 0
        else:
            self.ar_gout = np.zeros(nins)
            self.ar_in = np.zeros(nins)
        if (nout == 1):
            self.ar_out = 0
            self.ar_gin = 0
        else:
            self.ar_out = np.zeros(nout)
            self.ar_gin = np.zeros(nout)
    def reset(self):
        if (self.nout == 1):
            self.ar_out = 0
        else:
            self.ar_out[:] = 0
    def resetg(self):
        if (self.nins == 1):
            self.ar_gout = 0
        else:
            self.ar_gout[:] = 0

class BeadR(Bead):
    def forward(self, arri_in, arrj_in): # takes [r_i(x,y,z), r_j(x,y,z)] (6dim.) as input
        for k in range(3):
            self.ar_in[k  ] = arri_in[k]
            self.ar_in[k+3] = arrj_in[k]
        dx = arrj_in[0] - arri_in[0]
        dy = arrj_in[1] - arri_in[1]
        dz = arrj_in[2] - arri_in[2]        
        self.ar_out = np.sqrt(dx*dx + dy*dy + dz*dz)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        xx = ar_gin / self.ar_out 
        self.ar_gout[0] = xx * (self.ar_in[0] - self.ar_in[3]) 
        self.ar_gout[1] = xx * (self.ar_in[1] - self.ar_in[4]) 
        self.ar_gout[2] = xx * (self.ar_in[2] - self.ar_in[5]) 
        self.ar_gout[3] = -self.ar_gout[0]
        self.ar_gout[4] = -self.ar_gout[1]
        self.ar_gout[5] = -self.ar_gout[2]
        return self.ar_gout

class BeadCT(Bead):
    def forward(self, ar_in): # takes r_i,r_j,r_k,r_ij,r_ik,r_jk as input
        self.ar_in = ar_in
        b = ar_in[9]; c = ar_in[10]; a = ar_in[11]
        self.ar_out = (b*b+c*c-a*a)/(2.0*b*c)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        b = self.ar_in[9]; c = self.ar_in[10]; a = self.ar_in[11]
        dcosjik_drij = 1.0/c - 1.0/b*self.ar_out
        dcosjik_drik = 1.0/b - 1.0/c*self.ar_out
        dcosjik_drjk = -a/b/c
        drij_dRi = [(x-y)/b for (x,y) in zip(self.ar_in[0:3],self.ar_in[3:6])]
        drik_dRi = [(x-y)/c for (x,y) in zip(self.ar_in[0:3],self.ar_in[6:9])]
        ri_g = [x*dcosjik_drij+y*dcosjik_drik for (x,y) in zip(drij_dRi,drik_dRi)]
        drij_dRj = [(y-x)/b for (x,y) in zip(self.ar_in[0:3],self.ar_in[3:6])]
        drjk_dRj = [(x-y)/a for (x,y) in zip(self.ar_in[3:6],self.ar_in[6:9])]
        rj_g = [x*dcosjik_drij+y*dcosjik_drjk for (x,y) in zip(drij_dRj,drjk_dRj)]
        drik_dRk = [(y-x)/c for (x,y) in zip(self.ar_in[0:3],self.ar_in[6:9])]
        drjk_dRk = [(y-x)/a for (x,y) in zip(self.ar_in[3:6],self.ar_in[6:9])]
        rk_g = [x*dcosjik_drik+y*dcosjik_drjk for (x,y) in zip(drik_dRk,drjk_dRk)]
        xx = ri_g + rj_g + rk_g
        xx.append(0); xx.append(0); xx.append(0)
        self.ar_gout = [n*ar_gin for n in xx]
        return self.ar_gout

class BeadFC(Bead):
    def forward(self, ar_in): # takes r_ij as input
        self.ar_in = ar_in
        self.ar_out = fcut(ar_in)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gout = ar_gin * dfcut(self.ar_in)
        return self.ar_gout

class BeadFR(Bead):
    def forward(self, ar_in): # takes r_ij as input
        self.ar_in = ar_in
        self.ar_out = Ter_A * np.exp(-Ter_lam1*ar_in)
        return self.ar_out
    def backward(self, ar_gin):
        dfr_drij = -Ter_lam1 * self.ar_out
        self.ar_gout = ar_gin * dfr_drij
        return self.ar_gout

class BeadFA(Bead):
    def forward(self, ar_in): # takes r_ij as input
        self.ar_in = ar_in
        self.ar_out = -Ter_B * np.exp(-Ter_lam2*ar_in)
        return self.ar_out
    def backward(self, ar_gin):
        dfa_drij = -Ter_lam2 * self.ar_out
        self.ar_gout = ar_gin * dfa_drij
        return self.ar_gout

class BeadG(Bead):
    def forward(self, ar_in): # takes cosjik as input
        self.ar_in = ar_in
        self.ar_out = 1.0 + (Ter_c/Ter_d)**2 - Ter_c**2 / (Ter_d**2 + (Ter_h - ar_in)**2)
        return self.ar_out
    def backward(self, ar_gin):
        dg_dcos = -2.0 * Ter_c**2 * (Ter_h - self.ar_in) / (Ter_d**2 + (Ter_h - self.ar_in)**2)**2
        self.ar_gout = ar_gin * dg_dcos
        return self.ar_gout

class BeadZ(Bead):
    def forward(self, ar_in): # takes g_jik,r_ij,r_ik as input
        self.ar_in = ar_in
        g = ar_in[0]; rij = ar_in[1]; rik = ar_in[2]
        if (fcut(rik) != 0.0):
            self.ar_out = fcut(rik) * g * np.exp((Ter_lam3 * (rij - rik))**3)
        else:
            self.ar_out = 0.0
        return self.ar_out
    def backward(self, ar_gin):
        rij = self.ar_in[1]; rik = self.ar_in[2]
        if (fcut(rik) != 0.0):
            self.ar_gout[0] = ar_gin * fcut(rik) * np.exp(Ter_lam3**3 * (rij - rik)**3)
        else:
            self.ar_gout[0] = 0.0
        self.ar_gout[1] =  ar_gin * 3.0*Ter_lam3**3*(rij-rik)**2 * self.ar_out
        self.ar_gout[2] = -self.ar_gout[1]
        if (fcut(rik) != 0.0 and dfcut(rik) != 0.0):
            self.ar_gout[2] += ar_gin * (self.ar_out/fcut(rik)*dfcut(rik))
        return self.ar_gout

class BeadZIJ(Bead):
    def forward(self, ar_in): # takes zet_jik as input
        self.ar_in = ar_in
        self.ar_out += self.ar_in
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        self.ar_gout = self.ar_gin
        return self.ar_gout

class BeadBIJ(Bead):
    def forward(self, ar_in): # takes z_ij as input
        self.ar_in = ar_in
        self.ar_out = np.power(1.0 + (Ter_bet*ar_in)**Ter_n, -0.5/Ter_n)
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        x = (-0.5/Ter_n) * self.ar_out / (1.0 + (Ter_bet*self.ar_in)**Ter_n)
        self.ar_gout = ar_gin * x * (Ter_n * Ter_bet**Ter_n * self.ar_in**(Ter_n-1))
        return self.ar_gout

class BeadEIJ(Bead):
    def forward(self, ar_in): # takes fc,fr,bij,fa as input
        self.ar_in = ar_in
        self.ar_out = ar_in[0] * (ar_in[1] + ar_in[2] * ar_in[3])
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        fc_g  = ar_gin * (self.ar_in[1] + self.ar_in[2]*self.ar_in[3])
        fr_g  = ar_gin * self.ar_in[0]
        bij_g = ar_gin * self.ar_in[0] * self.ar_in[3]
        fa_g  = ar_gin * self.ar_in[0] * self.ar_in[2]
        self.ar_gout = [fc_g,fr_g,bij_g,fa_g]
        return self.ar_gout

class BeadEI(Bead):
    def forward(self, ar_in): # takes eij as input
        self.ar_in = ar_in
        self.ar_out += 0.5 * self.ar_in
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        dedp = 0.5
        self.ar_gout = ar_gin * dedp
        return self.ar_gout

class BeadE(Bead): # takes E_i as input
    def forward(self, ar_in):
        self.ar_in = ar_in
        self.ar_out += self.ar_in
        return self.ar_out
    def backward(self, ar_gin):
        self.ar_gin = ar_gin
        self.ar_gout = 1.0
        return self.ar_gout

rx = [] # m
ry = []
rz = []
rx_org = []
ry_org = []
rz_org = []
vx = [] # m/s
vy = []
vz = []
fx = [] # N
fy = []
fz = []
fx_n = [] # for calculation by numerical derivative
fy_n = []
fz_n = []
f = open('initial3d_Ter.d', 'r')
n = 0 # to count the number of total atoms
for line in f:
    xy = line.split()
    rx += [float(xy[0])]
    ry += [float(xy[1])]
    rz += [float(xy[2])]
    rx_org += [float(xy[0])]
    ry_org += [float(xy[1])]
    rz_org += [float(xy[2])]
    vx += [float(xy[3])]
    vy += [float(xy[4])]
    vz += [float(xy[5])]
    fx += [0]
    fy += [0]
    fz += [0]
    fx_n += [0]
    fy_n += [0]
    fz_n += [0]
    n += 1
print("number of atoms = ",n)

rij = BeadR(6,1)
rik = BeadR(6,1)
rjk = BeadR(6,1)
cosjik  = BeadCT(12,1)
fc  = BeadFC(1,1)
fr  = BeadFR(1,1)
fa  = BeadFA(1,1)
g   = BeadG(1,1)
zet = BeadZ(3,1)
zij = BeadZIJ(1,1)
bij = BeadBIJ(1,1)
eij = BeadEIJ(4,1)
ei  = BeadEI(1,1)
e   = BeadE(1,1)

e.reset()
for i in range(n):
    e.backward(1.0) #E_i*
    ei.reset()
    for j in range(n):
        if (j != i):
            ri = [rx[i],ry[i],rz[i]]
            rj = [rx[j],ry[j],rz[j]]
            rij.forward(ri,rj)
            fc.forward(rij.ar_out)
            fr.forward(rij.ar_out)
            fa.forward(rij.ar_out)
            zij.reset()
            rij_g = 0.0
            for k in range(n):
                if (k != i and k != j):
                    rk = [rx[k],ry[k],rz[k]]
                    rik.forward(ri,rk)
                    rjk.forward(rj,rk)
                    xx = ri+rj+rk
                    xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                    cosjik.forward(xx)
                    g.forward(cosjik.ar_out)
                    xx = [g.ar_out,rij.ar_out,rik.ar_out]
                    zet.forward(xx)
                    zij.forward(zet.ar_out)
            bij.forward(zij.ar_out)
            xx = [fc.ar_out,fr.ar_out,bij.ar_out,fa.ar_out]
            eij.forward(xx)
            ei.forward(eij.ar_out)
            #
            ei.backward(e.ar_gout) #eij*
            xx = eij.backward(ei.ar_gout) #bij* (using backward of BeadEIJ)
            fc_g  = xx[0]; fr_g  = xx[1]; bij_g = xx[2]; fa_g  = xx[3] 
            bij.backward(eij.ar_gout[2]) #zij*
            zjik_g = zij.backward(bij.ar_gout) #zjik*
            rij_g += fc.backward(fc_g)
            rij_g += fr.backward(fr_g)
            rij_g += fa.backward(fa_g)
            for k in range(n):
                if (k != i and k != j):
                    # POINT!!
                    rk = [rx[k],ry[k],rz[k]]
                    rik.forward(ri,rk)
                    rjk.forward(rj,rk)
                    xx = ri+rj+rk
                    xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                    cosjik.forward(xx)
                    g.forward(cosjik.ar_out)
                    xx = [g.ar_out,rij.ar_out,rik.ar_out]
                    zet.forward(xx)
                    #
                    xx = zet.backward(zjik_g) #gjik*,rij*,rik*
                    gjik_g = xx[0]; rij_g += xx[1]; rik_g = xx[2] # POINT!! (treatment of rik)
                    #redo cosjik.forward..
                    rk = [rx[k],ry[k],rz[k]]
                    rik.forward(ri,rk)
                    rjk.forward(rj,rk)
                    # POINT!! (treatment of rik)
                    xxx = rik.backward(rik_g)
                    fx[i] -= xxx[0]; fy[i] -= xxx[1]; fz[i] -= xxx[2]
                    fx[k] -= xxx[3]; fy[k] -= xxx[4]; fz[k] -= xxx[5]
                    #
                    xx = ri+rj+rk
                    xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                    cosjik.forward(xx)
                    #then backward..
                    cosjik_g = g.backward(zet.ar_gout[0]) #cosjik*
                    xx = cosjik.backward(g.ar_gout) #Ri*,Rj*,Rk*,rij*,rik*,rjk*
                    fx[i] -= xx[0]; fy[i] -= xx[1]; fz[i] -= xx[2]
                    fx[j] -= xx[3]; fy[j] -= xx[4]; fz[j] -= xx[5]
                    fx[k] -= xx[6]; fy[k] -= xx[7]; fz[k] -= xx[8]
            xxx = rij.backward(rij_g)
            fx[i] -= xxx[0]; fy[i] -= xxx[1]; fz[i] -= xxx[2]
            fx[j] -= xxx[3]; fy[j] -= xxx[4]; fz[j] -= xxx[5]
    ee = e.forward(ei.ar_out)
print('Etot = ',ee)

# numerical derivative
e0 = ee
ddx = 0.001
for k in range(n):
    for l in range(3):
        # plus
        for i in range(n):
            rx[i] = rx_org[i]
            ry[i] = ry_org[i]
            rz[i] = rz_org[i]
        if (l==0):
            rx[k] = rx_org[k] + ddx
        if (l==1):
            ry[k] = ry_org[k] + ddx
        if (l==2):
            rz[k] = rz_org[k] + ddx
        e.reset()
        for i in range(n):
            ei.reset()
            for j in range(n):
                if (j != i):
                    ri = [rx[i],ry[i],rz[i]]
                    rj = [rx[j],ry[j],rz[j]]
                    rij.forward(ri,rj)
                    fc.forward(rij.ar_out)
                    fr.forward(rij.ar_out)
                    fa.forward(rij.ar_out)
                    zij.reset()
                    rij_g = 0.0
                    for kk in range(n):
                        if (kk != i and kk != j):
                            rk = [rx[kk],ry[kk],rz[kk]]
                            rik.forward(ri,rk)
                            rjk.forward(rj,rk)
                            xx = ri+rj+rk
                            xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                            cosjik.forward(xx)
                            g.forward(cosjik.ar_out)
                            xx = [g.ar_out,rij.ar_out,rik.ar_out]
                            zet.forward(xx)
                            zij.forward(zet.ar_out)
                    bij.forward(zij.ar_out)
                    xx = [fc.ar_out,fr.ar_out,bij.ar_out,fa.ar_out]
                    eij.forward(xx)
                    ei.forward(eij.ar_out)
            eep = e.forward(ei.ar_out)
        # minus
        for i in range(n):
            rx[i] = rx_org[i]
            ry[i] = ry_org[i]
            rz[i] = rz_org[i]
        if (l==0):
            rx[k] = rx_org[k] - ddx
        if (l==1):
            ry[k] = ry_org[k] - ddx
        if (l==2):
            rz[k] = rz_org[k] - ddx
        e.reset()
        for i in range(n):
            ei.reset()
            for j in range(n):
                if (j != i):
                    ri = [rx[i],ry[i],rz[i]]
                    rj = [rx[j],ry[j],rz[j]]
                    rij.forward(ri,rj)
                    fc.forward(rij.ar_out)
                    fr.forward(rij.ar_out)
                    fa.forward(rij.ar_out)
                    zij.reset()
                    rij_g = 0.0
                    for kk in range(n):
                        if (kk != i and kk != j):
                            rk = [rx[kk],ry[kk],rz[kk]]
                            rik.forward(ri,rk)
                            rjk.forward(rj,rk)
                            xx = ri+rj+rk
                            xx.append(rij.ar_out); xx.append(rik.ar_out); xx.append(rjk.ar_out)
                            cosjik.forward(xx)
                            g.forward(cosjik.ar_out)
                            xx = [g.ar_out,rij.ar_out,rik.ar_out]
                            zet.forward(xx)
                            zij.forward(zet.ar_out)
                    bij.forward(zij.ar_out)
                    xx = [fc.ar_out,fr.ar_out,bij.ar_out,fa.ar_out]
                    eij.forward(xx)
                    ei.forward(eij.ar_out)
            eem = e.forward(ei.ar_out)
        numderiv = -(eep-eem)/ddx/2.0
        if (l==0):
            fx_n[k] = numderiv
        if (l==1):
            fy_n[k] = numderiv
        if (l==2):
            fz_n[k] = numderiv
for i in range(n):
    print('Analytical:',i,fx[i],fy[i],fz[i])
    print('Numerical :',i,fx_n[i],fy_n[i],fz_n[i])

このコードは，系のエネルギーと原子に働く力を求めるとともに，原子に微小変位を加えることでエネルギーの数値微分を行い，それによって求めた力の値と比べることで，微分がうまく行っていることを確認している．
原子の配置は以下のファイルで与えられる．

initial3d_Ter.d

3.0 2.0 0.0 0 0 0
5.0 2.0 0.0 0 0 0
4.0 3.0 2.0 0 0 0
2.4 1.6 1.0 0 0 0

左からx,y,z座標（単位はÅ）を与えている．4〜6番コラムは初期速度を与えるために用意しているがここでは用いられない．上記の例では原子数は4である．

実行結果

number of atoms =  4
Etot =  -269.3394974652807
Analytical: 0 126.71182641298195 32.495350662534534 -8.55921044903075
Numerical : 0 126.71161555047661 32.49532862065507 -8.55927491184616
Analytical: 1 -137.24788412101856 14.26757417134369 50.6967763138305
Numerical : 1 -137.24789247572744 14.267580700561666 50.69680315415326
Analytical: 2 -38.54319669778745 -71.72987606445713 -96.67701946018232
Numerical : 2 -38.54316748359565 -71.72986382244062 -96.67708448523626
Analytical: 3 49.07925440582404 24.96695123057889 54.53945359538258
Numerical : 3 49.07913418702492 24.96699702840033 54.53948254029228

back propagationによる微分で求めた力（Analytical）と，数値微分により求めたもの（Numerical）がほぼ一致していることがわかる．