added Makefile + changes to appendix and sim chapter
[lectures/latex.git] / posic / thesis / d_tersoff.tex
1 \chapter{Force evaluation for the three body Tersoff potential}
2 \label{app:d_tersoff}
3
4   \section{Form of the Tersoff potential and its derivative}
5
6 The Tersoff potential is of the form
7 \begin{eqnarray}
8 E & = & \sum_i E_i = \frac{1}{2} \sum_{i \ne j} V_{ij} \textrm{ ,} \\
9 V_{ij} & = & f_C(r_{ij}) [ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) ] \textrm{ .}
10 \end{eqnarray}
11 The repulsive $f_R$ and attractive $f_A$ part is given by
12 \begin{eqnarray}
13 f_R(r_{ij}) & = & A_{ij} \exp (- \lambda_{ij} r_{ij} ) \textrm{ ,} \\
14 f_A(r_{ij}) & = & -B_{ij} \exp (- \mu_{ij} r_{ij} ) \textrm{ .}
15 \end{eqnarray}
16 The bond order function $b_{ij}$ is
17 \begin{equation}
18 b_{ij} = \chi_{ij} (1 + \beta_i^{n_i} \zeta^{n_i}_{ij})^{-1/2n_i}
19 \end{equation}
20 with
21 \begin{eqnarray}
22 \zeta_{ij} & = & \sum_{k \ne i,j} f_C (r_{ik}) \omega_{ik} g(\theta_{ijk}) \textrm{ ,}\\
23 g(\theta_{ijk}) & = & 1 + c_i^2/d_i^2 - c_i^2/[d_i^2 + (h_i - \cos \theta_{ijk})^2] \textrm{ .}
24 \end{eqnarray}
25 The cutoff function $f_C$ is taken to be
26 \begin{equation}
27 f_C(r_{ij}) = \left\{
28   \begin{array}{ll}
29     1, & r_{ij} < R_{ij} \\
30     \frac{1}{2} + \frac{1}{2} \cos \Big[ \pi (r_{ij} - R_{ij})/(S_{ij} - R_{ij}) \Big], &
31  R_{ij} < r_{ij} < S_{ij} \\
32     0, & r_{ij} > S_{ij}
33   \end{array} \right.
34 \end{equation}
35 with $\theta_{ijk}$ being the bond angle between bonds $ij$ and $ik$ as shown in Figure \ref{img:tersoff_angle}.\\
36 \\
37 For a three body potential, if $V_{ij}$ is not equal to $V_{ji}$, the derivative is of the form
38 \begin{equation}
39 \nabla_{{\bf r}_i} E = \frac{1}{2} \big[ \sum_j ( \nabla_{{\bf r}_i} V_{ij} + \nabla_{{\bf r}_i} V_{ji} ) + \sum_k \sum_j \nabla_{{\bf r}_i} V_{jk} \big] \textrm{ .}
40 \end{equation}
41 In the following all the necessary derivatives to calculate $\nabla_{{\bf r}_i} E$ are done.
42
43   \section{Derivative of $V_{ij}$ with respect to ${\bf r}_i$}
44
45 \begin{eqnarray}
46 \nabla_{{\bf r}_i} V_{ij} & = & \nabla_{{\bf r}_i} f_C(r_{ij}) \big[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) \big] + \nonumber \\
47  & & + f_C(r_{ij}) \big[ \nabla_{{\bf r}_i} f_R(r_{ij}) + b_{ij} \nabla_{{\bf r}_i} f_A(r_{ij}) + f_A(r_{ij}) \nabla_{{\bf r}_i} b_{ij} \big]
48 \end{eqnarray}
49 \begin{eqnarray}
50 \nabla_{{\bf r}_i} f_R(r_{ij}) & = & A_{ij} \lambda_{ij} \frac{{\bf r}_{ij}}{r_{ij}} \exp(-\lambda_{ij} r_{ij}) \\
51 \nabla_{{\bf r}_i} f_A(r_{ij}) & = & - B_{ij} \mu_{ij} \frac{{\bf r}_{ij}}{r_{ij}} \exp(-\mu_{ij} r_{ij})
52 \end{eqnarray}
53 \begin{equation}
54 \nabla_{{\bf r}_i} f_C(r_{ij}) = \left\{
55   \begin{array}{ll}
56     \frac{1}{2} \sin \Big( \frac{\pi(r_{ij}-R_{ij})}{S_{ij}-R_{ij}} \Big) \frac{\pi}{S_{ij}-R_{ij}} \frac{{\bf r}_{ij}}{r_{ij}}, & R_{ij} < r_{ij} < S_{ij} \\
57     0, & \textrm{else.}
58   \end{array} \right.
59 \end{equation}
60 \begin{eqnarray}
61 \nabla_{{\bf r}_i} b_{ij} &= & - \frac{\chi_{ij}}{2} (1+\beta^{n_i} \zeta_{ij}^{n_i})^{-\frac{1}{2n_i}-1} \beta^{n_i} \zeta_{ij}^{n_i-1} \nabla_{{\bf r}_i} \zeta_{ij} \\
62 \nabla_{{\bf r}_i} \zeta_{ij} & = & \sum_{k \neq i,j} \big( g(\theta_{ijk}) \nabla_{{\bf r}_i} f_C(r_{ik}) + f_C(r_{ik}) \nabla_{{\bf r}_i} g(\theta_{ijk}) \big) \\
63 \nabla_{{\bf r}_i} g(\theta_{ijk}) & = & - \frac{2(h_i-\cos\theta_{ijk})c_i^2}{\big[d_i^2 + (h_i - \cos\theta_{ijk})^2\big]^2} \nabla_{{\bf r}_i} (\cos\theta_{ijk})
64 \end{eqnarray}
65 \begin{eqnarray}
66 \nabla_{{\bf r}_i} \cos \theta_{ijk} & = & \nabla_{{\bf r}_i} \Big( \frac{{\bf r}_{ij} {\bf r}_{ik}}{r_{ij} r_{ik}} \Big) \nonumber \\
67  & = & \Big[ \frac{\cos\theta_{ijk}}{r_{ij}^2} - \frac{1}{r_{ij} r_{ik}} \Big] {\bf r}_{ij} + \Big[ \frac{\cos\theta_{ijk}}{r_{ik}^2} - \frac{1}{r_{ij} r_{ik}} \Big] {\bf r}_{ik}
68 \end{eqnarray}
69
70   \section{Derivative of $V_{ji}$ with respect to ${\bf r}_i$}
71
72 \begin{eqnarray}
73 \nabla_{{\bf r}_i} V_{ji} & = & \nabla_{{\bf r}_i} f_C(r_{ji}) \big[ f_R(r_{ji}) + b_{ji} f_A(r_{ji}) \big] + \nonumber \\
74  & & + f_C(r_{ji}) \big[ \nabla_{{\bf r}_i} f_R(r_{ji}) + b_{ji} \nabla_{{\bf r}_i} f_A(r_{ji}) + f_A(r_{ji}) \nabla_{{\bf r}_i} b_{ji} \big]
75 \end{eqnarray}
76 \begin{eqnarray}
77 \nabla_{{\bf r}_i} f_R(r_{ji}) & = & - A_{ji} \lambda_{ji} \frac{{\bf r}_{ji}}{r_{ji}} \exp(-\lambda_{ji} r_{ji}) = \nabla_{{\bf r}_i} f_R(r_{ij}) \\
78 \nabla_{{\bf r}_i} f_A(r_{ji}) & = & + B_{ji} \mu_{ji} \frac{{\bf r}_{ji}}{r_{ji}} \exp(-\mu_{ji} r_{ji}) = \nabla_{{\bf r}_i} f_A(r_{ij})
79 \end{eqnarray}
80 \begin{equation}
81 \nabla_{{\bf r}_i} f_C(r_{ij}) = \nabla_{{\bf r}_i} f_C(r_{ij}) = \left\{
82   \begin{array}{ll}
83     - \frac{1}{2} \sin \Big( \frac{\pi(r_{ji}-R_{ji})}{S_{ji}-R_{ji}} \Big) \frac{\pi}{S_{ji}-R_{ji}} \frac{{\bf r}_{ji}}{r_{ji}}, & R_{ji} < r_{ji} < S_{ji} \\
84     0, & \textrm{else.}
85   \end{array} \right.
86 \end{equation}
87 \begin{eqnarray}
88 \nabla_{{\bf r}_i} b_{ji} &= & - \frac{\chi_{ji}}{2} (1+\beta^{n_j} \zeta_{ji}^{n_j})^{-\frac{1}{2n_j}-1} \beta^{n_j} \zeta_{ji}^{n_j-1} \nabla_{{\bf r}_i} \zeta_{ji} \\
89 \nabla_{{\bf r}_i} \zeta_{ji} & = & \sum_{k \neq j,i} \big( g(\theta_{jik}) \nabla_{{\bf r}_i} f_C(r_{jk}) + f_C(r_{jk}) \nabla_{{\bf r}_i} g(\theta_{jik}) \big) \nonumber \\
90  & = & \sum_{k \neq j,i} f_C(r_{jk}) \nabla_{{\bf r}_i} g(\theta_{jik}) \quad \Big(\textrm{Reason: }\nabla_{{\bf r}_i} f_C(r_{jk}) = 0\Big) \\
91 \nabla_{{\bf r}_i} g(\theta_{jik}) & = & - \frac{2(h_j-\cos\theta_{jik})c_j^2}{\big[d_j^2 + (h_j - \cos\theta_{jik})^2\big]^2} \nabla_{{\bf r}_i} (\cos\theta_{jik})
92 \end{eqnarray}
93 \begin{eqnarray}
94 \nabla_{{\bf r}_i} \cos \theta_{jik} & = & \nabla_{{\bf r}_i} \Big( \frac{{\bf r}_{ji} {\bf r}_{jk}}{r_{ji} r_{jk}} \Big) \nonumber \\
95  & = & \frac{1}{r_{ji} r_{jk}} {\bf r}_{jk} - \frac{\cos\theta_{jik}}{r_{ji}^2} {\bf r}_{ji}
96 \end{eqnarray}
97
98   \section{Derivative of $V_{jk}$ with respect to ${\bf r}_i$}
99
100 \begin{eqnarray}
101 \nabla_{{\bf r}_i} V_{jk} & = & f_C(r_{jk}) f_A(r_{jk}) \nabla_{{\bf r}_i} b_{jk} \\
102 \nabla_{{\bf r}_i} b_{jk} & = & - \frac{\chi_{jk}}{2} (1+\beta^{n_j} \zeta_{jk}^{n_j})^{-\frac{1}{2n_j}-1} \beta^{n_j} \zeta_{jk}^{n_j-1} \nabla_{{\bf r}_i} \zeta_{jk} \\
103 \nabla_{{\bf r}_i} \zeta_{jk} & = & \sum_{l \neq j,k} \big( g(\theta_{jkl}) \nabla_{{\bf r}_i} f_C(r_{jl}) + f_C(r_{jl}) \nabla_{{\bf r}_i} g(\theta_{jkl}) \big) \nonumber \\
104  & = & f_C(r_{ji}) \nabla_{{\bf r}_i} g(\theta_{jki}) + g(\theta_{jki}) \nabla_{{\bf r}_i} f_C(r_{ji}) \\
105 \nabla_{{\bf r}_i} g(\theta_{jki}) & = & - \frac{2(h_j-\cos\theta_{jki})c_j^2}{\big[d_j^2 + (h_j - \cos\theta_{jki})^2\big]^2} \nabla_{{\bf r}_i} (\cos\theta_{jki}) \\
106 \nabla_{{\bf r}_i} \cos \theta_{jki} & = & \nabla_{{\bf r}_i} \Big( \frac{{\bf r}_{jk} {\bf r}_{ji}}{r_{jk} r_{ji}} \Big) \nonumber \\
107  & = & \frac{1}{r_{jk} r_{ji}} {\bf r}_{jk} - \frac{\cos\theta_{jki}}{r_{ji}^2} {\bf r}_{ji}
108 \end{eqnarray}
109
110   \section{Implementation issues}
111
112 As seen in the last sections the derivatives of $V_{ij}$, $V_{ji}$ and $V_{jk}$
113 with respect to ${\bf r}_i$ are necessary to compute the forces for atom $i$.
114 According to this, for every triple $(ijk)$ the derivatives of the three
115 potential contributions, denoted by $V_{ijk}$, $V_{jik}$ and $V_{jki}$
116 have to be computed.
117 In simulation, however, it is not practical to evaluate all three potential
118 derivatives for each $(ijk)$ triple.
119 The $V_{jik}$ and $V_{jki}$ potential and its derivatives will be calculated
120 in subsequent loops anyways.
121 To avoid multiple computation of the same potential derivatives
122 the force contributions for atom $j$ and $k$ due to the $V_{ijk}$ contribution
123 have to be considered by calculating the derivatives of $V_{ijk}$
124 with respect to ${\bf r}_j$ and ${\bf r}_k$
125 inside the loop of the $(ijk)$ triple
126 in addition to the derivative with respect to ${\bf r}_i$.
127 This poses a more convenient method to obtain the forces
128 keeping in mind that all the necessary force contributions for atom $i$
129 are calculated and added in subsequent loops.
130
131 The following symmetry considerations help to obtain the 
132
133   \subsection{Derivative of $V_{ij}$ with respect to ${\bf r}_j$}
134
135 \begin{eqnarray}
136  \nabla_{{\bf r}_j} V_{ij} & = &
137  \nabla_{{\bf r}_j} f_C(r_{ij}) \big[ f_R(r_{ij}) +
138  b_{ij} f_A(r_{ij}) \big] + \nonumber \\
139  & & + f_C(r_{ij}) \big[ \nabla_{{\bf r}_j} f_R(r_{ij}) +
140  b_{ij} \nabla_{{\bf r}_j} f_A(r_{ij}) +
141  f_A(r_{ij}) \nabla_{{\bf r}_j} b_{ij} \big]
142 \end{eqnarray}
143 Using the equality $\nabla_{{\bf r}_i} r_{ij}=-\nabla_{{\bf r}_j} r_{ij}$
144 the following relations are valid:
145 \begin{eqnarray}
146  \nabla_{{\bf r}_j} f_R(r_{ij}) &=& - \nabla_{{\bf r}_i} f_R(r_{ij}) \\
147  \nabla_{{\bf r}_j} f_A(r_{ij}) &=& - \nabla_{{\bf r}_i} f_A(r_{ij}) \\
148  \nabla_{{\bf r}_j} f_C(r_{ij}) &=& - \nabla_{{\bf r}_i} f_C(r_{ij})
149 \end{eqnarray}
150 The pair contributions .... easy.
151 Now having a look at $b_{ij}$.
152 \begin{eqnarray}
153  \nabla_{{\bf r}_j}\cos\theta_{ijk} &=&
154  \nabla_{{\bf r}_j}\Big(\frac{{\bf r}_{ij}{\bf }r_{ik}}{r_{ij}r_{ik}}\Big)
155  \nonumber \\
156  &=& \frac{1}{r_{ij}r_{ik}}{\bf r}_{ik} -
157      \frac{\cos\theta_{ijk}}{r_{ij}^2}{\bf r}_{ij}
158 \end{eqnarray}
159
160   \subsection{Derivative of $V_{ij}$ with respect to ${\bf r}_k$}
161
162 The derivative of $V_{ij}$ with respect to ${\bf r}_k$ just consists of the
163 single term
164 \begin{eqnarray}
165  \nabla_{{\bf r}_k} V_{ij} & = &
166  f_C(r_{ij})f_A(r_{ij})\nabla_{{\bf r}_{k}}b_{ij}
167 \end{eqnarray}
168 since the derivatives of the functions only depending on atom $i$ and $j$
169 vanish.
170 \begin{eqnarray}
171  \nabla_{{\bf r}_k} f_R(r_{ij}) &=& 0 \\
172  \nabla_{{\bf r}_k} f_A(r_{ij}) &=& 0 \\
173  \nabla_{{\bf r}_k} f_C(r_{ij}) &=& 0
174 \end{eqnarray}
175 Now look at $b_{ij}$, not only angle important here!
176 \begin{eqnarray}
177  \nabla_{{\bf r}_k}\zeta_{ij} &=&
178  g(\theta_{ijk})\nabla_{{\bf r}_k}f_C(r_{ik}) +
179  f_C(r_{ik})\nabla_{{\bf r}_k}g(\theta_{ijk}) \\
180  \nabla_{{\bf r}_k}f_C(r_{ik}) &=& - \nabla_{{\bf r}_i}f_C(r_{ik}) \\
181  \nabla_{{\bf r}_k}\cos\theta_{ijk} &=&
182  \nabla_{{\bf r}_k}\Big(\frac{{\bf r}_{ij}{\bf }r_{ik}}{r_{ij}r_{ik}}\Big)
183  \nonumber \\
184  &=&\frac{1}{r_{ij}r_{ik}}{\bf r}_{ij} -
185  \frac{\cos\theta_{ijk}}{r_{ik}^2}{\bf r}_{ik}
186 \end{eqnarray}
187
188   \subsection{Code realization}
189
190 The implementation of the force evaluation shown in the following
191 is applied to the potential designed by Erhard and Albe.
192 There are slight differences comparted to the original potential by Tersoff:
193 \begin{itemize}
194  \item Difference in sign of the attractive part.
195  \item $c$, $d$ and $h$ values depend on atom $k$ in addition to atom $i$.
196  \item Difference in sign of the $\cos\theta_{ijk}$ term.
197  \item There are no parameters $\beta$ and $\chi$.
198  \item The exponent of the $b$ term is constantly $-\frac{1}{2}$.
199 \end{itemize}
200 These differences actually slightly ease code realization.
201
202 \begin{figure}
203 \renewcommand\labelitemi{}
204 \renewcommand\labelitemii{}
205 \renewcommand\labelitemiii{}
206 LOOP i \{
207 \begin{itemize}
208  \item // nop (only used in orig. Tersoff)
209  \item LOOP j \{
210        \begin{itemize}
211         \item // nop (only used in orig. Tersoff)
212        \end{itemize}
213  \item \}
214  \item LOOP j \{
215        \begin{itemize}
216         \item $\zeta_{ij}=0$
217         \item set $S_{ij}$ (cutoff)
218         \item calculate: $r_{ij}$, $r_{ij}^2$
219         \item IF $r_{ij} > S_{ij}$ THEN CONTINUE
220         \item
221         \item LOOP k \{
222               \begin{itemize}
223                \item set $ik$-depending values
224                \item calculate: $r_{ik}$, $r_{ik}^2$
225                \item IF $r_{ik} > S_{ik}$ THEN CONTINUE
226                \item calculate: $\theta_{ijk}$, $\cos(\theta_{ijk})$,
227                                 $dg(\theta_{ijk})$, $g(\theta)$,
228                                 $f_C(r_{ik})$, $df_C(r_{ik})$
229                \item $\zeta_{ij}=\zeta_{ij}+f_C(r_{ik})g(\theta)$
230               \end{itemize}
231         \item \}
232         \item
233         \item calculate: $f_C(r_{ij})$, $df_C(r_{ij})$, $f_A(r_{ij})$,
234                          $df_A(r_{ij})$, $f_R(r_{ij})$, $df_R(r_{ij})$,
235                          $b_{ij}$, $db_{ij}$
236         \item calculate:
237 $
238 F=-\frac{1}{2}\big(
239 \nabla_{{\bf r}_i} f_C(r_{ij}) \big[ f_R(r_{ij}) - b_{ij} f_A(r_{ij}) \big] +
240 f_C(r_{ij}) \big[ \nabla_{{\bf r}_i} f_R(r_{ij}) -
241                   b_{ij} \nabla_{{\bf r}_i} f_A(r_{ij}) \big]\big)
242 $
243         \item $F_{Atom\, i} = F_{Atom\, i} + F$
244         \item $F_{Atom\, j} = F_{Atom\, j} - F$
245         \item $E=E+\frac{1}{2}f_C(r_{ij})[f_R(r_{ij})-b_{ij}f_A(r_{ij})]$
246         \item $d\zeta_{ij}=\frac{1}{2}f_A(r_{ij})f_C(r_{ij})db_{ij}$
247         \item
248         \item LOOP k \{
249               \begin{itemize}
250                \item calculate: $\nabla_{{\bf r}_i}\cos\theta_{ijk}$,
251                                 $\nabla_{{\bf r}_j}\cos\theta_{ijk}$,
252                                 $\nabla_{{\bf r}_k}\cos\theta_{ijk}$
253                \item $
254 F_{Atom\, i}+= d\zeta_{ij}\big(
255 g(\theta_{ijk})\nabla_{{\bf r}_i}f_C(r_{ik}) +
256 f_C(r_{ik})dg(\theta_{ijk})\nabla_{{\bf r}_i}\cos\theta_{ijk}\big)$
257                \item $
258 F_{Atom\, j}+= d\zeta_{ij}
259 f_C(r_{ik})dg(\theta_{ijk})\nabla_{{\bf r}_j}\cos\theta_{ijk}$
260                \item $
261 F_{Atom\, k}+= d\zeta_{ij}\big(
262 g(\theta_{ijk})\nabla_{{\bf r}_k}f_C(r_{ik}) +
263 f_C(r_{ik})dg(\theta_{ijk})\nabla_{{\bf r}_k}\cos\theta_{ijk}\big)$
264               \end{itemize}
265         \item \}
266        \end{itemize}
267  \item \}
268 \end{itemize}
269 \}
270 \caption{Implementation of the force evaluation for Tersoff like bond-order
271          potentials using pseudocode.}
272 \end{figure}
273