first sent away beta version
[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 \cite{tersoff_m} 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 cut-off 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 \subsection{Derivative of $V_{ij}$ with respect to ${\bf r}_j$}
132
133 \begin{eqnarray}
134  \nabla_{{\bf r}_j} V_{ij} & = &
135  \nabla_{{\bf r}_j} f_C(r_{ij}) \big[ f_R(r_{ij}) +
136  b_{ij} f_A(r_{ij}) \big] + \nonumber \\
137  & & + f_C(r_{ij}) \big[ \nabla_{{\bf r}_j} f_R(r_{ij}) +
138  b_{ij} \nabla_{{\bf r}_j} f_A(r_{ij}) +
139  f_A(r_{ij}) \nabla_{{\bf r}_j} b_{ij} \big]
140 \end{eqnarray}
141 Using the equality $\nabla_{{\bf r}_i} r_{ij}=-\nabla_{{\bf r}_j} r_{ij}$
142 the following relations are valid:
143 \begin{eqnarray}
144  \nabla_{{\bf r}_j} f_R(r_{ij}) &=& - \nabla_{{\bf r}_i} f_R(r_{ij}) \\
145  \nabla_{{\bf r}_j} f_A(r_{ij}) &=& - \nabla_{{\bf r}_i} f_A(r_{ij}) \\
146  \nabla_{{\bf r}_j} f_C(r_{ij}) &=& - \nabla_{{\bf r}_i} f_C(r_{ij})
147 \end{eqnarray}
148 The pair contributions are, thus, easily obtained.
149 The contribution of the bond order term is given by:
150 \begin{eqnarray}
151  \nabla_{{\bf r}_j}\cos\theta_{ijk} &=&
152  \nabla_{{\bf r}_j}\Big(\frac{{\bf r}_{ij}{\bf }r_{ik}}{r_{ij}r_{ik}}\Big)
153  \nonumber \\
154  &=& \frac{1}{r_{ij}r_{ik}}{\bf r}_{ik} -
155      \frac{\cos\theta_{ijk}}{r_{ij}^2}{\bf r}_{ij}
156 \end{eqnarray}
157
158 \subsection{Derivative of $V_{ij}$ with respect to ${\bf r}_k$}
159
160 The derivative of $V_{ij}$ with respect to ${\bf r}_k$ just consists of the
161 single term
162 \begin{eqnarray}
163  \nabla_{{\bf r}_k} V_{ij} & = &
164  f_C(r_{ij})f_A(r_{ij})\nabla_{{\bf r}_{k}}b_{ij}
165 \end{eqnarray}
166 since the derivatives of the functions only depending on atom $i$ and $j$
167 vanish.
168 \begin{eqnarray}
169  \nabla_{{\bf r}_k} f_R(r_{ij}) &=& 0 \\
170  \nabla_{{\bf r}_k} f_A(r_{ij}) &=& 0 \\
171  \nabla_{{\bf r}_k} f_C(r_{ij}) &=& 0
172 \end{eqnarray}
173 Concerning $b_{ij}$, in addition to the angular term, the derivative of the cut-off function has to be considered.
174 \begin{eqnarray}
175  \nabla_{{\bf r}_k}\zeta_{ij} &=&
176  g(\theta_{ijk})\nabla_{{\bf r}_k}f_C(r_{ik}) +
177  f_C(r_{ik})\nabla_{{\bf r}_k}g(\theta_{ijk}) \\
178  \nabla_{{\bf r}_k}f_C(r_{ik}) &=& - \nabla_{{\bf r}_i}f_C(r_{ik}) \\
179  \nabla_{{\bf r}_k}\cos\theta_{ijk} &=&
180  \nabla_{{\bf r}_k}\Big(\frac{{\bf r}_{ij}{\bf r}_{ik}}{r_{ij}r_{ik}}\Big)
181  \nonumber \\
182  &=&\frac{1}{r_{ij}r_{ik}}{\bf r}_{ij} -
183  \frac{\cos\theta_{ijk}}{r_{ik}^2}{\bf r}_{ik}
184 \end{eqnarray}
185
186   \subsection{Code realization}
187
188 The implementation of the force evaluation shown in the following is applied to the potential designed by Erhard and Albe \cite{albe_sic_pot}.
189 There are slight differences compared to the original potential by Tersoff:
190 \begin{itemize}
191  \item Difference in sign of the attractive part.
192  \item $c$, $d$ and $h$ values depend on atom $k$ in addition to atom $i$.
193  \item Difference in sign of the $\cos\theta_{ijk}$ term.
194  \item There are no parameters $\beta$ and $\chi$.
195  \item The exponent of the $b$ term is constantly $-\frac{1}{2}$.
196 \end{itemize}
197 These differences actually slightly ease code realization.
198 The respective flow chart is displayed in Fig.~\ref{fig:flowchart}.
199
200 \begin{figure}
201 \renewcommand\labelitemi{}
202 \renewcommand\labelitemii{}
203 \renewcommand\labelitemiii{}
204 {\small
205 \fbox{\begin{minipage}{\textwidth}
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}$ (cut-off)
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 \end{minipage}}
271 }
272 \caption{Flow chart of the force evaluation for Tersoff-like bond order potentials using pseudocode.}
273 \label{fig:flowchart}
274 \end{figure}
275