1 /// **************************************************************************
2 // tersoff_zbl_extra.h
3 // -------------------
4 // Trung Dac Nguyen
5 //
6 // Device code for Tersoff math routines
7 //
8 // __________________________________________________________________________
9 // This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
10 // __________________________________________________________________________
11 //
12 // begin :
13 // email : ndactrung@gmail.com
14 // ***************************************************************************/*
15
16 #ifndef LAL_TERSOFF_ZBL_EXTRA_H
17 #define LAL_TERSOFF_ZBL_EXTRA_H
18
19 #if defined(NV_KERNEL) || defined(USE_HIP)
20 #include "lal_aux_fun1.h"
21 #else
22 #endif
23
24 #define MY_PI (numtyp)3.14159265358979323846
25 #define MY_PI2 (numtyp)1.57079632679489661923
26 #define MY_PI4 (numtyp)0.78539816339744830962
27
28 /* ---------------------------------------------------------------------- */
29
vec3_dot(const numtyp x[3],const numtyp y[3])30 ucl_inline numtyp vec3_dot(const numtyp x[3], const numtyp y[3])
31 {
32 return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]);
33 }
34
vec3_add(const numtyp x[3],const numtyp y[3],numtyp z[3])35 ucl_inline void vec3_add(const numtyp x[3], const numtyp y[3], numtyp z[3])
36 {
37 z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2];
38 }
39
vec3_scale(const numtyp k,const numtyp x[3],numtyp y[3])40 ucl_inline void vec3_scale(const numtyp k, const numtyp x[3], numtyp y[3])
41 {
42 y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2];
43 }
44
vec3_scaleadd(const numtyp k,const numtyp x[3],const numtyp y[3],numtyp z[3])45 ucl_inline void vec3_scaleadd(const numtyp k, const numtyp x[3],
46 const numtyp y[3], numtyp z[3])
47 {
48 z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2];
49 }
50
51 /* ---------------------------------------------------------------------- */
52
ters_gijk(const numtyp costheta,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma)53 ucl_inline numtyp ters_gijk(const numtyp costheta,
54 const numtyp param_c,
55 const numtyp param_d,
56 const numtyp param_h,
57 const numtyp param_gamma)
58 {
59 const numtyp ters_c = param_c * param_c;
60 const numtyp ters_d = param_d * param_d;
61 const numtyp hcth = param_h - costheta;
62 return param_gamma*((numtyp)1.0 + ters_c*ucl_recip(ters_d) -
63 ters_c *ucl_recip(ters_d + hcth*hcth));
64 }
65
66 /* ---------------------------------------------------------------------- */
67
ters_gijk_d(const numtyp costheta,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma)68 ucl_inline numtyp ters_gijk_d(const numtyp costheta,
69 const numtyp param_c,
70 const numtyp param_d,
71 const numtyp param_h,
72 const numtyp param_gamma)
73 {
74 const numtyp ters_c = param_c * param_c;
75 const numtyp ters_d = param_d * param_d;
76 const numtyp hcth = param_h - costheta;
77 const numtyp numerator = (numtyp)-2.0 * ters_c * hcth;
78 const numtyp denominator = ucl_recip(ters_d + hcth*hcth);
79 return param_gamma*numerator*denominator*denominator;
80 }
81
82 /* ---------------------------------------------------------------------- */
83
costheta_d(const numtyp rij_hat[3],const numtyp rij,const numtyp rik_hat[3],const numtyp rik,numtyp * dri,numtyp * drj,numtyp * drk)84 ucl_inline void costheta_d(const numtyp rij_hat[3],
85 const numtyp rij,
86 const numtyp rik_hat[3],
87 const numtyp rik,
88 numtyp *dri,
89 numtyp *drj,
90 numtyp *drk)
91 {
92 // first element is derivative wrt Ri, second wrt Rj, third wrt Rk
93
94 numtyp cos_theta = vec3_dot(rij_hat,rik_hat);
95
96 vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj);
97 vec3_scale(ucl_recip(rij),drj,drj);
98 vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk);
99 vec3_scale(ucl_recip(rik),drk,drk);
100 vec3_add(drj,drk,dri);
101 vec3_scale((numtyp)-1.0,dri,dri);
102 }
103
104 /* ---------------------------------------------------------------------- */
105
ters_fc(const numtyp r,const numtyp param_bigr,const numtyp param_bigd)106 ucl_inline numtyp ters_fc(const numtyp r,
107 const numtyp param_bigr,
108 const numtyp param_bigd)
109 {
110 if (r < param_bigr-param_bigd) return (numtyp)1.0;
111 if (r > param_bigr+param_bigd) return (numtyp)0.0;
112 return (numtyp)0.5*((numtyp)1.0 - sin(MY_PI2*(r - param_bigr)/param_bigd));
113 }
114
115 /* ---------------------------------------------------------------------- */
116
ters_fc_d(const numtyp r,const numtyp param_bigr,const numtyp param_bigd)117 ucl_inline numtyp ters_fc_d(const numtyp r,
118 const numtyp param_bigr,
119 const numtyp param_bigd)
120 {
121 if (r < param_bigr-param_bigd) return (numtyp)0.0;
122 if (r > param_bigr+param_bigd) return (numtyp)0.0;
123 return -(MY_PI4/param_bigd) * cos(MY_PI2*(r - param_bigr)/param_bigd);
124 }
125
126 /* ---------------------------------------------------------------------- */
127
F_fermi(const numtyp r,const numtyp param_ZBLcut,const numtyp param_ZBLexpscale)128 ucl_inline numtyp F_fermi(const numtyp r, const numtyp param_ZBLcut,
129 const numtyp param_ZBLexpscale)
130 {
131 return ucl_recip((numtyp)1.0+ucl_exp(-param_ZBLexpscale*(r-param_ZBLcut)));
132 }
133
134 /* ---------------------------------------------------------------------- */
135
F_fermi_d(const numtyp r,const numtyp param_ZBLcut,const numtyp param_ZBLexpscale)136 ucl_inline numtyp F_fermi_d(const numtyp r, const numtyp param_ZBLcut,
137 const numtyp param_ZBLexpscale)
138 {
139 numtyp a = ucl_exp(-param_ZBLexpscale*(r-param_ZBLcut));
140 numtyp b = (numtyp)1.0 + a;
141 return param_ZBLexpscale*a*ucl_recip(b*b);
142 }
143
144 /* ---------------------------------------------------------------------- */
145
ters_fa(const numtyp r,const numtyp param_bigb,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_lam2,const numtyp param_ZBLcut,const numtyp param_ZBLexpscale)146 ucl_inline numtyp ters_fa(const numtyp r,
147 const numtyp param_bigb,
148 const numtyp param_bigr,
149 const numtyp param_bigd,
150 const numtyp param_lam2,
151 const numtyp param_ZBLcut,
152 const numtyp param_ZBLexpscale)
153 {
154 if (r > param_bigr + param_bigd) return (numtyp)0.0;
155 return -param_bigb * ucl_exp(-param_lam2 * r) *
156 ters_fc(r,param_bigr,param_bigd)*F_fermi(r,param_ZBLcut,param_ZBLexpscale);
157 }
158
159 /* ---------------------------------------------------------------------- */
160
ters_fa_d(const numtyp r,const numtyp param_bigb,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_lam2,const numtyp param_ZBLcut,const numtyp param_ZBLexpscale)161 ucl_inline numtyp ters_fa_d(const numtyp r,
162 const numtyp param_bigb,
163 const numtyp param_bigr,
164 const numtyp param_bigd,
165 const numtyp param_lam2,
166 const numtyp param_ZBLcut,
167 const numtyp param_ZBLexpscale)
168 {
169 if (r > param_bigr + param_bigd) return (numtyp)0.0;
170 numtyp f = F_fermi(r,param_ZBLcut,param_ZBLexpscale);
171 return param_bigb * ucl_exp(-param_lam2 * r) *
172 (param_lam2 * ters_fc(r,param_bigr,param_bigd) * f -
173 ters_fc_d(r,param_bigr,param_bigd) * f -
174 ters_fc(r,param_bigr,param_bigd) * F_fermi_d(r,param_ZBLcut,param_ZBLexpscale));
175 }
176
177 /* ---------------------------------------------------------------------- */
178
ters_bij(const numtyp zeta,const numtyp param_beta,const numtyp param_powern,const numtyp param_c1,const numtyp param_c2,const numtyp param_c3,const numtyp param_c4)179 ucl_inline numtyp ters_bij(const numtyp zeta,
180 const numtyp param_beta,
181 const numtyp param_powern,
182 const numtyp param_c1,
183 const numtyp param_c2,
184 const numtyp param_c3,
185 const numtyp param_c4)
186 {
187 numtyp tmp = param_beta * zeta;
188 if (tmp > param_c1) return ucl_rsqrt(tmp);
189 if (tmp > param_c2)
190 return ((numtyp)1.0 - ucl_powr(tmp,-param_powern) /
191 ((numtyp)2.0*param_powern))*ucl_rsqrt(tmp);
192 if (tmp < param_c4) return (numtyp)1.0;
193 if (tmp < param_c3)
194 return (numtyp)1.0 - ucl_powr(tmp,param_powern)/((numtyp)2.0*param_powern);
195 return ucl_powr((numtyp)1.0 + ucl_powr(tmp,param_powern),
196 (numtyp)-1.0/((numtyp)2.0*param_powern));
197 }
198
199 /* ---------------------------------------------------------------------- */
200
ters_bij_d(const numtyp zeta,const numtyp param_beta,const numtyp param_powern,const numtyp param_c1,const numtyp param_c2,const numtyp param_c3,const numtyp param_c4)201 ucl_inline numtyp ters_bij_d(const numtyp zeta,
202 const numtyp param_beta,
203 const numtyp param_powern,
204 const numtyp param_c1,
205 const numtyp param_c2,
206 const numtyp param_c3,
207 const numtyp param_c4)
208 {
209 numtyp tmp = param_beta * zeta;
210 if (tmp > param_c1)
211 return param_beta * (numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5);
212 if (tmp > param_c2)
213 return param_beta * ((numtyp)-0.5*ucl_powr(tmp,(numtyp)-1.5) *
214 // error in negligible 2nd term fixed 9/30/2015
215 // (1.0 - 0.5*(1.0 + 1.0/(2.0*param->powern)) *
216 ((numtyp)1.0 - ((numtyp)1.0 + (numtyp)1.0 /((numtyp)2.0 * param_powern)) *
217 ucl_powr(tmp,-param_powern)));
218 if (tmp < param_c4) return (numtyp)0.0;
219 if (tmp < param_c3)
220 return (numtyp)-0.5*param_beta * ucl_powr(tmp,param_powern-(numtyp)1.0);
221
222 numtyp tmp_n = ucl_powr(tmp,param_powern);
223 return (numtyp)-0.5 * ucl_powr((numtyp)1.0+tmp_n, (numtyp) -
224 (numtyp)1.0-((numtyp)1.0 / ((numtyp)2.0 * param_powern)))*tmp_n / zeta;
225 }
226
227 /* ---------------------------------------------------------------------- */
228
ters_zetaterm_d(const numtyp prefactor,const numtyp rij_hat[3],const numtyp rij,const numtyp rik_hat[3],const numtyp rik,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,numtyp dri[3],numtyp drj[3],numtyp drk[3])229 ucl_inline void ters_zetaterm_d(const numtyp prefactor,
230 const numtyp rij_hat[3],
231 const numtyp rij,
232 const numtyp rik_hat[3],
233 const numtyp rik,
234 const numtyp param_bigr,
235 const numtyp param_bigd,
236 const numtyp param_powermint,
237 const numtyp param_lam3,
238 const numtyp param_c,
239 const numtyp param_d,
240 const numtyp param_h,
241 const numtyp param_gamma,
242 numtyp dri[3],
243 numtyp drj[3],
244 numtyp drk[3])
245 {
246 numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
247 numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
248
249 fc = ters_fc(rik,param_bigr,param_bigd);
250 dfc = ters_fc_d(rik,param_bigr,param_bigd);
251
252 numtyp t = param_lam3*(rij-rik);
253 if ((int)param_powermint == 3) tmp = t*t*t;
254 else tmp = t;
255
256 if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
257 else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
258 else ex_delr = ucl_exp(tmp);
259
260 if ((int)param_powermint == 3)
261 ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
262 else ex_delr_d = param_lam3 * ex_delr;
263
264 cos_theta = vec3_dot(rij_hat,rik_hat);
265 gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
266 gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
267 costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
268
269 // compute the derivative wrt Ri
270 // dri = -dfc*gijk*ex_delr*rik_hat;
271 // dri += fc*gijk_d*ex_delr*dcosdri;
272 // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
273
274 vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
275 vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
276 vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
277 vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
278 vec3_scale(prefactor,dri,dri);
279
280 // compute the derivative wrt Rj
281 // drj = fc*gijk_d*ex_delr*dcosdrj;
282 // drj += fc*gijk*ex_delr_d*rij_hat;
283
284 vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
285 vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
286 vec3_scale(prefactor,drj,drj);
287
288 // compute the derivative wrt Rk
289 // drk = dfc*gijk*ex_delr*rik_hat;
290 // drk += fc*gijk_d*ex_delr*dcosdrk;
291 // drk += -fc*gijk*ex_delr_d*rik_hat;
292
293 vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
294 vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
295 vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
296 vec3_scale(prefactor,drk,drk);
297 }
298
ters_zetaterm_d_fi(const numtyp prefactor,const numtyp rij_hat[3],const numtyp rij,const numtyp rik_hat[3],const numtyp rik,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,numtyp dri[3])299 ucl_inline void ters_zetaterm_d_fi(const numtyp prefactor,
300 const numtyp rij_hat[3],
301 const numtyp rij,
302 const numtyp rik_hat[3],
303 const numtyp rik,
304 const numtyp param_bigr,
305 const numtyp param_bigd,
306 const numtyp param_powermint,
307 const numtyp param_lam3,
308 const numtyp param_c,
309 const numtyp param_d,
310 const numtyp param_h,
311 const numtyp param_gamma,
312 numtyp dri[3])
313 {
314 numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
315 numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
316
317 fc = ters_fc(rik,param_bigr,param_bigd);
318 dfc = ters_fc_d(rik,param_bigr,param_bigd);
319
320 numtyp t = param_lam3*(rij-rik);
321 if ((int)param_powermint == 3) tmp = t*t*t;
322 else tmp = t;
323
324 if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
325 else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
326 else ex_delr = ucl_exp(tmp);
327
328 if ((int)param_powermint == 3)
329 ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
330 else ex_delr_d = param_lam3 * ex_delr;
331
332 cos_theta = vec3_dot(rij_hat,rik_hat);
333 gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
334 gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
335 costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
336
337 // compute the derivative wrt Ri
338 // dri = -dfc*gijk*ex_delr*rik_hat;
339 // dri += fc*gijk_d*ex_delr*dcosdri;
340 // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
341
342 vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
343 vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
344 vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
345 vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
346 vec3_scale(prefactor,dri,dri);
347 }
348
ters_zetaterm_d_fj(const numtyp prefactor,const numtyp rij_hat[3],const numtyp rij,const numtyp rik_hat[3],const numtyp rik,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,numtyp drj[3])349 ucl_inline void ters_zetaterm_d_fj(const numtyp prefactor,
350 const numtyp rij_hat[3],
351 const numtyp rij,
352 const numtyp rik_hat[3],
353 const numtyp rik,
354 const numtyp param_bigr,
355 const numtyp param_bigd,
356 const numtyp param_powermint,
357 const numtyp param_lam3,
358 const numtyp param_c,
359 const numtyp param_d,
360 const numtyp param_h,
361 const numtyp param_gamma,
362 numtyp drj[3])
363 {
364 numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,cos_theta,tmp;
365 numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
366
367 fc = ters_fc(rik,param_bigr,param_bigd);
368
369 numtyp t = param_lam3*(rij-rik);
370 if ((int)param_powermint == 3) tmp = t*t*t;
371 else tmp = t;
372
373 if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
374 else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
375 else ex_delr = ucl_exp(tmp);
376
377 if ((int)param_powermint == 3)
378 ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
379 else ex_delr_d = param_lam3 * ex_delr;
380
381 cos_theta = vec3_dot(rij_hat,rik_hat);
382 gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
383 gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
384 costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
385
386 // compute the derivative wrt Rj
387 // drj = fc*gijk_d*ex_delr*dcosdrj;
388 // drj += fc*gijk*ex_delr_d*rij_hat;
389
390 vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
391 vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
392 vec3_scale(prefactor,drj,drj);
393 }
394
ters_zetaterm_d_fk(const numtyp prefactor,const numtyp rij_hat[3],const numtyp rij,const numtyp rik_hat[3],const numtyp rik,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,numtyp drk[3])395 ucl_inline void ters_zetaterm_d_fk(const numtyp prefactor,
396 const numtyp rij_hat[3],
397 const numtyp rij,
398 const numtyp rik_hat[3],
399 const numtyp rik,
400 const numtyp param_bigr,
401 const numtyp param_bigd,
402 const numtyp param_powermint,
403 const numtyp param_lam3,
404 const numtyp param_c,
405 const numtyp param_d,
406 const numtyp param_h,
407 const numtyp param_gamma,
408 numtyp drk[3])
409 {
410 numtyp gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
411 numtyp dcosdri[3],dcosdrj[3],dcosdrk[3];
412
413 fc = ters_fc(rik,param_bigr,param_bigd);
414 dfc = ters_fc_d(rik,param_bigr,param_bigd);
415
416 numtyp t = param_lam3*(rij-rik);
417 if ((int)param_powermint == 3) tmp = t*t*t;
418 else tmp = t;
419
420 if (tmp > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
421 else if (tmp < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
422 else ex_delr = ucl_exp(tmp);
423
424 if ((int)param_powermint == 3)
425 ex_delr_d = (numtyp)3.0*param_lam3*t*t*ex_delr;
426 else ex_delr_d = param_lam3 * ex_delr;
427
428 cos_theta = vec3_dot(rij_hat,rik_hat);
429 gijk = ters_gijk(cos_theta,param_c,param_d,param_h,param_gamma);
430 gijk_d = ters_gijk_d(cos_theta,param_c,param_d,param_h,param_gamma);
431 costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
432
433 // compute the derivative wrt Rk
434 // drk = dfc*gijk*ex_delr*rik_hat;
435 // drk += fc*gijk_d*ex_delr*dcosdrk;
436 // drk += -fc*gijk*ex_delr_d*rik_hat;
437
438 vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
439 vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
440 vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
441 vec3_scale(prefactor,drk,drk);
442 }
443
444 /* ---------------------------------------------------------------------- */
445
repulsive(const numtyp param_bigr,const numtyp param_bigd,const numtyp param_lam1,const numtyp param_biga,const numtyp param_Z_i,const numtyp param_Z_j,const numtyp param_ZBLcut,const numtyp param_ZBLexpscale,const numtyp global_e,const numtyp global_a_0,const numtyp global_epsilon_0,const numtyp rsq,const int eflag,numtyp * ans)446 ucl_inline void repulsive(const numtyp param_bigr,
447 const numtyp param_bigd,
448 const numtyp param_lam1,
449 const numtyp param_biga,
450 const numtyp param_Z_i,
451 const numtyp param_Z_j,
452 const numtyp param_ZBLcut,
453 const numtyp param_ZBLexpscale,
454 const numtyp global_e,
455 const numtyp global_a_0,
456 const numtyp global_epsilon_0,
457 const numtyp rsq,
458 const int eflag,
459 numtyp *ans)
460 {
461 numtyp r,tmp_fc,tmp_fc_d,tmp_exp;
462
463 // Tersoff repulsive portion
464
465 r = ucl_sqrt(rsq);
466 tmp_fc = ters_fc(r,param_bigr,param_bigd);
467 tmp_fc_d = ters_fc_d(r,param_bigr,param_bigd);
468 tmp_exp = ucl_exp(-param_lam1 * r);
469
470 numtyp fforce_ters = param_biga * tmp_exp * (tmp_fc_d - tmp_fc*param_lam1);
471 numtyp eng_ters = tmp_fc * param_biga * tmp_exp;
472
473 // ZBL repulsive portion
474
475 numtyp esq = global_e*global_e;
476 numtyp a_ij = ((numtyp)0.8854*global_a_0) /
477 (ucl_powr(param_Z_i,(numtyp)0.23) + ucl_powr(param_Z_j,(numtyp)0.23));
478 numtyp premult = (param_Z_i * param_Z_j * esq)/((numtyp)4.0*MY_PI*global_epsilon_0);
479 numtyp r_ov_a = r/a_ij;
480 numtyp t1 = (numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a);
481 numtyp t2 = (numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a);
482 numtyp t3 = (numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a);
483 numtyp t4 = (numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a);
484 numtyp phi = t1 + t2 + t3 + t4;
485 numtyp dphi = (numtyp)-3.2*t1 - (numtyp)0.9423*t2 - (numtyp)0.4029*t3 -
486 (numtyp)0.2016*t4;
487 dphi *= ucl_recip(a_ij);
488 /*
489 numtyp phi = (numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a) +
490 (numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a) +
491 (numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a) +
492 (numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a);
493 numtyp dphi = ucl_recip(a_ij) * ((numtyp)-3.2*(numtyp)0.1818*ucl_exp((numtyp)-3.2*r_ov_a) -
494 (numtyp)0.9423*(numtyp)0.5099*ucl_exp((numtyp)-0.9423*r_ov_a) -
495 (numtyp)0.4029*(numtyp)0.2802*ucl_exp((numtyp)-0.4029*r_ov_a) -
496 (numtyp)0.2016*(numtyp)0.02817*ucl_exp((numtyp)-0.2016*r_ov_a));
497 */
498 numtyp rinv = ucl_recip(r);
499 numtyp fforce_ZBL = premult*(-phi)/rsq + premult*dphi*rinv;
500 numtyp eng_ZBL = premult*rinv*phi;
501
502 // combine two parts with smoothing by Fermi-like function
503 // ans[0] = fforce
504 numtyp f = F_fermi(r,param_ZBLcut,param_ZBLexpscale);
505 numtyp f_d = F_fermi_d(r,param_ZBLcut,param_ZBLexpscale);
506 ans[0] = -(-f_d * eng_ZBL + ((numtyp)1.0 - f)*fforce_ZBL + f_d*eng_ters +
507 f*fforce_ters) * rinv;
508
509 // ans[1] = eng
510 if (eflag) ans[1] = ((numtyp)1.0 - f)*eng_ZBL + f*eng_ters;
511 }
512
513 /* ---------------------------------------------------------------------- */
514
zeta(const numtyp param_powermint,const numtyp param_lam3,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,const numtyp rsqij,const numtyp rsqik,const numtyp4 delrij,const numtyp4 delrik)515 ucl_inline numtyp zeta(const numtyp param_powermint,
516 const numtyp param_lam3,
517 const numtyp param_bigr,
518 const numtyp param_bigd,
519 const numtyp param_c,
520 const numtyp param_d,
521 const numtyp param_h,
522 const numtyp param_gamma,
523 const numtyp rsqij,
524 const numtyp rsqik,
525 const numtyp4 delrij,
526 const numtyp4 delrik)
527 {
528 numtyp rij,rik,costheta,arg,ex_delr;
529
530 rij = ucl_sqrt(rsqij);
531 rik = ucl_sqrt(rsqik);
532 costheta = (delrij.x*delrik.x + delrij.y*delrik.y +
533 delrij.z*delrik.z) / (rij*rik);
534
535 numtyp t = param_lam3*(rij-rik);
536 if ((int)param_powermint == 3) arg = t*t*t;
537 else arg = t;
538
539 if (arg > (numtyp)69.0776) ex_delr = (numtyp)1.e30;
540 else if (arg < (numtyp)-69.0776) ex_delr = (numtyp)0.0;
541 else ex_delr = ucl_exp(arg);
542
543 return ters_fc(rik,param_bigr,param_bigd) *
544 ters_gijk(costheta,param_c, param_d, param_h, param_gamma) * ex_delr;
545 }
546
547 /* ---------------------------------------------------------------------- */
548
force_zeta(const numtyp param_bigb,const numtyp param_bigr,const numtyp param_bigd,const numtyp param_lam2,const numtyp param_beta,const numtyp param_powern,const numtyp param_c1,const numtyp param_c2,const numtyp param_c3,const numtyp param_c4,const numtyp param_ZBLcut,const numtyp param_ZBLexpscale,const numtyp rsq,const numtyp zeta_ij,const int eflag,numtyp fpfeng[4])549 ucl_inline void force_zeta(const numtyp param_bigb,
550 const numtyp param_bigr,
551 const numtyp param_bigd,
552 const numtyp param_lam2,
553 const numtyp param_beta,
554 const numtyp param_powern,
555 const numtyp param_c1,
556 const numtyp param_c2,
557 const numtyp param_c3,
558 const numtyp param_c4,
559 const numtyp param_ZBLcut,
560 const numtyp param_ZBLexpscale,
561 const numtyp rsq,
562 const numtyp zeta_ij,
563 const int eflag,
564 numtyp fpfeng[4])
565 {
566 numtyp r,fa,fa_d,bij;
567
568 r = ucl_sqrt(rsq);
569 fa = ters_fa(r,param_bigb,param_bigr,param_bigd,param_lam2,param_ZBLcut,param_ZBLexpscale);
570 fa_d = ters_fa_d(r,param_bigb,param_bigr,param_bigd,param_lam2,param_ZBLcut,param_ZBLexpscale);
571 bij = ters_bij(zeta_ij,param_beta,param_powern,
572 param_c1,param_c2, param_c3, param_c4);
573 fpfeng[0] = (numtyp)0.5*bij*fa_d * ucl_recip(r); // fforce
574 fpfeng[1] = (numtyp)-0.5*fa * ters_bij_d(zeta_ij,param_beta, param_powern,
575 param_c1,param_c2, param_c3, param_c4); // prefactor
576 if (eflag) fpfeng[2] = (numtyp)0.5*bij*fa; // eng
577 }
578
579 /* ----------------------------------------------------------------------
580 attractive term
581 use param_ij cutoff for rij test
582 use param_ijk cutoff for rik test
583 ------------------------------------------------------------------------- */
584
attractive(const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,const numtyp prefactor,const numtyp rij,const numtyp rijinv,const numtyp rik,const numtyp rikinv,const numtyp delrij[3],const numtyp delrik[3],numtyp fi[3],numtyp fj[3],numtyp fk[3])585 ucl_inline void attractive(const numtyp param_bigr,
586 const numtyp param_bigd,
587 const numtyp param_powermint,
588 const numtyp param_lam3,
589 const numtyp param_c,
590 const numtyp param_d,
591 const numtyp param_h,
592 const numtyp param_gamma,
593 const numtyp prefactor,
594 const numtyp rij,
595 const numtyp rijinv,
596 const numtyp rik,
597 const numtyp rikinv,
598 const numtyp delrij[3],
599 const numtyp delrik[3],
600 numtyp fi[3],
601 numtyp fj[3],
602 numtyp fk[3])
603 {
604 numtyp rij_hat[3],rik_hat[3];
605 vec3_scale(rijinv,delrij,rij_hat);
606 vec3_scale(rikinv,delrik,rik_hat);
607 ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,
608 param_bigr, param_bigd, param_powermint, param_lam3,
609 param_c, param_d, param_h, param_gamma, fi, fj, fk);
610 }
611
attractive_fi(const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,const numtyp prefactor,const numtyp rij,const numtyp rijinv,const numtyp rik,const numtyp rikinv,const numtyp delrij[3],const numtyp delrik[3],numtyp fi[3])612 ucl_inline void attractive_fi(const numtyp param_bigr,
613 const numtyp param_bigd,
614 const numtyp param_powermint,
615 const numtyp param_lam3,
616 const numtyp param_c,
617 const numtyp param_d,
618 const numtyp param_h,
619 const numtyp param_gamma,
620 const numtyp prefactor,
621 const numtyp rij,
622 const numtyp rijinv,
623 const numtyp rik,
624 const numtyp rikinv,
625 const numtyp delrij[3],
626 const numtyp delrik[3],
627 numtyp fi[3])
628 {
629 numtyp rij_hat[3],rik_hat[3];
630 vec3_scale(rijinv,delrij,rij_hat);
631 vec3_scale(rikinv,delrik,rik_hat);
632 ters_zetaterm_d_fi(prefactor,rij_hat,rij,rik_hat,rik,
633 param_bigr, param_bigd, param_powermint, param_lam3,
634 param_c, param_d, param_h, param_gamma, fi);
635 }
636
attractive_fj(const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,const numtyp prefactor,const numtyp rij,const numtyp rijinv,const numtyp rik,const numtyp rikinv,const numtyp delrij[3],const numtyp delrik[3],numtyp fj[3])637 ucl_inline void attractive_fj(const numtyp param_bigr,
638 const numtyp param_bigd,
639 const numtyp param_powermint,
640 const numtyp param_lam3,
641 const numtyp param_c,
642 const numtyp param_d,
643 const numtyp param_h,
644 const numtyp param_gamma,
645 const numtyp prefactor,
646 const numtyp rij,
647 const numtyp rijinv,
648 const numtyp rik,
649 const numtyp rikinv,
650 const numtyp delrij[3],
651 const numtyp delrik[3],
652 numtyp fj[3])
653 {
654 numtyp rij_hat[3],rik_hat[3];
655 vec3_scale(rijinv,delrij,rij_hat);
656 vec3_scale(rikinv,delrik,rik_hat);
657 ters_zetaterm_d_fj(prefactor,rij_hat,rij,rik_hat,rik,
658 param_bigr, param_bigd, param_powermint, param_lam3,
659 param_c, param_d, param_h, param_gamma, fj);
660 }
661
attractive_fk(const numtyp param_bigr,const numtyp param_bigd,const numtyp param_powermint,const numtyp param_lam3,const numtyp param_c,const numtyp param_d,const numtyp param_h,const numtyp param_gamma,const numtyp prefactor,const numtyp rij,const numtyp rijinv,const numtyp rik,const numtyp rikinv,const numtyp delrij[3],const numtyp delrik[3],numtyp fk[3])662 ucl_inline void attractive_fk(const numtyp param_bigr,
663 const numtyp param_bigd,
664 const numtyp param_powermint,
665 const numtyp param_lam3,
666 const numtyp param_c,
667 const numtyp param_d,
668 const numtyp param_h,
669 const numtyp param_gamma,
670 const numtyp prefactor,
671 const numtyp rij,
672 const numtyp rijinv,
673 const numtyp rik,
674 const numtyp rikinv,
675 const numtyp delrij[3],
676 const numtyp delrik[3],
677 numtyp fk[3])
678 {
679 numtyp rij_hat[3],rik_hat[3];
680 vec3_scale(rijinv,delrij,rij_hat);
681 vec3_scale(rikinv,delrik,rik_hat);
682 ters_zetaterm_d_fk(prefactor,rij_hat,rij,rik_hat,rik,
683 param_bigr, param_bigd, param_powermint, param_lam3,
684 param_c, param_d, param_h, param_gamma, fk);
685 }
686
687
688 #endif
689
690
691