1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 #include "meam.h"
15
16 #include "math_special.h"
17
18 #include <cmath>
19 #include <algorithm>
20
21 using namespace LAMMPS_NS;
22
23 void
meam_force(int i,int eflag_global,int eflag_atom,int vflag_global,int vflag_atom,double * eng_vdwl,double * eatom,int,int * type,int * fmap,double ** scale,double ** x,int numneigh,int * firstneigh,int numneigh_full,int * firstneigh_full,int fnoffset,double ** f,double ** vatom,double * virial)24 MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int vflag_atom,
25 double* eng_vdwl, double* eatom, int /*ntype*/, int* type, int* fmap,
26 double** scale, double** x, int numneigh, int* firstneigh, int numneigh_full,
27 int* firstneigh_full, int fnoffset, double** f, double** vatom, double *virial)
28 {
29 int j, jn, k, kn, kk, m, n, p, q;
30 int nv2, nv3, elti, eltj, eltk, ind;
31 int eflag_either = eflag_atom || eflag_global;
32 int vflag_either = vflag_atom || vflag_global;
33 double xitmp, yitmp, zitmp, delij[3], rij2, rij, rij3;
34 double v[6], fi[3], fj[3];
35 double third, sixth;
36 double pp, dUdrij, dUdsij, dUdrijm[3], force, forcem;
37 double recip, phi, phip;
38 double sij;
39 double a1, a1i, a1j, a2, a2i, a2j;
40 double a3i, a3j;
41 double shpi[3], shpj[3];
42 double ai, aj, ro0i, ro0j, invrei, invrej;
43 double rhoa0j, drhoa0j, rhoa0i, drhoa0i;
44 double rhoa1j, drhoa1j, rhoa1i, drhoa1i;
45 double rhoa2j, drhoa2j, rhoa2i, drhoa2i;
46 double a3, a3a, rhoa3j, drhoa3j, rhoa3i, drhoa3i;
47 double drho0dr1, drho0dr2, drho0ds1, drho0ds2;
48 double drho1dr1, drho1dr2, drho1ds1, drho1ds2;
49 double drho1drm1[3], drho1drm2[3];
50 double drho2dr1, drho2dr2, drho2ds1, drho2ds2;
51 double drho2drm1[3], drho2drm2[3];
52 double drho3dr1, drho3dr2, drho3ds1, drho3ds2;
53 double drho3drm1[3], drho3drm2[3];
54 double dt1dr1, dt1dr2, dt1ds1, dt1ds2;
55 double dt2dr1, dt2dr2, dt2ds1, dt2ds2;
56 double dt3dr1, dt3dr2, dt3ds1, dt3ds2;
57 double drhodr1, drhodr2, drhods1, drhods2, drhodrm1[3], drhodrm2[3];
58 double arg;
59 double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3;
60 double dsij1, dsij2, force1, force2;
61 double t1i, t2i, t3i, t1j, t2j, t3j;
62 double scaleij;
63
64 third = 1.0 / 3.0;
65 sixth = 1.0 / 6.0;
66
67 // Compute forces atom i
68
69 elti = fmap[type[i]];
70 if (elti < 0) return;
71
72 xitmp = x[i][0];
73 yitmp = x[i][1];
74 zitmp = x[i][2];
75
76 // Treat each pair
77 for (jn = 0; jn < numneigh; jn++) {
78 j = firstneigh[jn];
79 eltj = fmap[type[j]];
80 scaleij = scale[type[i]][type[j]];
81
82 if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) {
83
84 sij = scrfcn[fnoffset + jn] * fcpair[fnoffset + jn];
85 delij[0] = x[j][0] - xitmp;
86 delij[1] = x[j][1] - yitmp;
87 delij[2] = x[j][2] - zitmp;
88 rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
89 if (rij2 < this->cutforcesq) {
90 rij = sqrt(rij2);
91 recip = 1.0 / rij;
92
93 // Compute phi and phip
94 ind = this->eltind[elti][eltj];
95 pp = rij * this->rdrar;
96 kk = (int)pp;
97 kk = std::min(kk, this->nrar - 2);
98 pp = pp - kk;
99 pp = std::min(pp, 1.0);
100 phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + this->phirar[ind][kk];
101 phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk];
102
103 if (eflag_either != 0) {
104 double phi_sc = phi * scaleij;
105 if (eflag_global != 0)
106 *eng_vdwl = *eng_vdwl + phi_sc * sij;
107 if (eflag_atom != 0) {
108 eatom[i] = eatom[i] + 0.5 * phi_sc * sij;
109 eatom[j] = eatom[j] + 0.5 * phi_sc * sij;
110 }
111 }
112
113 // write(1,*) "force_meamf: phi: ",phi
114 // write(1,*) "force_meamf: phip: ",phip
115
116 // Compute pair densities and derivatives
117 invrei = 1.0 / this->re_meam[elti][elti];
118 ai = rij * invrei - 1.0;
119 ro0i = this->rho0_meam[elti];
120 rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai);
121 drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i;
122 rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai);
123 drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i;
124 rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai);
125 drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i;
126 rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai);
127 drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i;
128
129 if (elti != eltj) {
130 invrej = 1.0 / this->re_meam[eltj][eltj];
131 aj = rij * invrej - 1.0;
132 ro0j = this->rho0_meam[eltj];
133 rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj);
134 drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j;
135 rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj);
136 drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j;
137 rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj);
138 drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j;
139 rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj);
140 drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j;
141 } else {
142 rhoa0j = rhoa0i;
143 drhoa0j = drhoa0i;
144 rhoa1j = rhoa1i;
145 drhoa1j = drhoa1i;
146 rhoa2j = rhoa2i;
147 drhoa2j = drhoa2i;
148 rhoa3j = rhoa3i;
149 drhoa3j = drhoa3i;
150 }
151
152 const double t1mi = this->t1_meam[elti];
153 const double t2mi = this->t2_meam[elti];
154 const double t3mi = this->t3_meam[elti];
155 const double t1mj = this->t1_meam[eltj];
156 const double t2mj = this->t2_meam[eltj];
157 const double t3mj = this->t3_meam[eltj];
158
159 if (this->ialloy == 1) {
160 rhoa1j *= t1mj;
161 rhoa2j *= t2mj;
162 rhoa3j *= t3mj;
163 rhoa1i *= t1mi;
164 rhoa2i *= t2mi;
165 rhoa3i *= t3mi;
166 drhoa1j *= t1mj;
167 drhoa2j *= t2mj;
168 drhoa3j *= t3mj;
169 drhoa1i *= t1mi;
170 drhoa2i *= t2mi;
171 drhoa3i *= t3mi;
172 }
173
174 nv2 = 0;
175 nv3 = 0;
176 arg1i1 = 0.0;
177 arg1j1 = 0.0;
178 arg1i2 = 0.0;
179 arg1j2 = 0.0;
180 arg1i3 = 0.0;
181 arg1j3 = 0.0;
182 arg3i3 = 0.0;
183 arg3j3 = 0.0;
184 for (n = 0; n < 3; n++) {
185 for (p = n; p < 3; p++) {
186 for (q = p; q < 3; q++) {
187 arg = delij[n] * delij[p] * delij[q] * this->v3D[nv3];
188 arg1i3 = arg1i3 + arho3[i][nv3] * arg;
189 arg1j3 = arg1j3 - arho3[j][nv3] * arg;
190 nv3 = nv3 + 1;
191 }
192 arg = delij[n] * delij[p] * this->v2D[nv2];
193 arg1i2 = arg1i2 + arho2[i][nv2] * arg;
194 arg1j2 = arg1j2 + arho2[j][nv2] * arg;
195 nv2 = nv2 + 1;
196 }
197 arg1i1 = arg1i1 + arho1[i][n] * delij[n];
198 arg1j1 = arg1j1 - arho1[j][n] * delij[n];
199 arg3i3 = arg3i3 + arho3b[i][n] * delij[n];
200 arg3j3 = arg3j3 - arho3b[j][n] * delij[n];
201 }
202
203 // rho0 terms
204 drho0dr1 = drhoa0j * sij;
205 drho0dr2 = drhoa0i * sij;
206
207 // rho1 terms
208 a1 = 2 * sij / rij;
209 drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1;
210 drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1;
211 a1 = 2.0 * sij / rij;
212 for (m = 0; m < 3; m++) {
213 drho1drm1[m] = a1 * rhoa1j * arho1[i][m];
214 drho1drm2[m] = -a1 * rhoa1i * arho1[j][m];
215 }
216
217 // rho2 terms
218 a2 = 2 * sij / rij2;
219 drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - 2.0 / 3.0 * arho2b[i] * drhoa2j * sij;
220 drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - 2.0 / 3.0 * arho2b[j] * drhoa2i * sij;
221 a2 = 4 * sij / rij2;
222 for (m = 0; m < 3; m++) {
223 drho2drm1[m] = 0.0;
224 drho2drm2[m] = 0.0;
225 for (n = 0; n < 3; n++) {
226 drho2drm1[m] = drho2drm1[m] + arho2[i][this->vind2D[m][n]] * delij[n];
227 drho2drm2[m] = drho2drm2[m] - arho2[j][this->vind2D[m][n]] * delij[n];
228 }
229 drho2drm1[m] = a2 * rhoa2j * drho2drm1[m];
230 drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m];
231 }
232
233 // rho3 terms
234 rij3 = rij * rij2;
235 a3 = 2 * sij / rij3;
236 a3a = 6.0 / 5.0 * sij / rij;
237 drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - a3a * (drhoa3j - rhoa3j / rij) * arg3i3;
238 drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - a3a * (drhoa3i - rhoa3i / rij) * arg3j3;
239 a3 = 6 * sij / rij3;
240 a3a = 6 * sij / (5 * rij);
241 for (m = 0; m < 3; m++) {
242 drho3drm1[m] = 0.0;
243 drho3drm2[m] = 0.0;
244 nv2 = 0;
245 for (n = 0; n < 3; n++) {
246 for (p = n; p < 3; p++) {
247 arg = delij[n] * delij[p] * this->v2D[nv2];
248 drho3drm1[m] = drho3drm1[m] + arho3[i][this->vind3D[m][n][p]] * arg;
249 drho3drm2[m] = drho3drm2[m] + arho3[j][this->vind3D[m][n][p]] * arg;
250 nv2 = nv2 + 1;
251 }
252 }
253 drho3drm1[m] = (a3 * drho3drm1[m] - a3a * arho3b[i][m]) * rhoa3j;
254 drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i;
255 }
256
257 // Compute derivatives of weighting functions t wrt rij
258 t1i = t_ave[i][0];
259 t2i = t_ave[i][1];
260 t3i = t_ave[i][2];
261 t1j = t_ave[j][0];
262 t2j = t_ave[j][1];
263 t3j = t_ave[j][2];
264
265 if (this->ialloy == 1) {
266
267 a1i = fdiv_zero(drhoa0j * sij, tsq_ave[i][0]);
268 a1j = fdiv_zero(drhoa0i * sij, tsq_ave[j][0]);
269 a2i = fdiv_zero(drhoa0j * sij, tsq_ave[i][1]);
270 a2j = fdiv_zero(drhoa0i * sij, tsq_ave[j][1]);
271 a3i = fdiv_zero(drhoa0j * sij, tsq_ave[i][2]);
272 a3j = fdiv_zero(drhoa0i * sij, tsq_ave[j][2]);
273
274 dt1dr1 = a1i * (t1mj - t1i * MathSpecial::square(t1mj));
275 dt1dr2 = a1j * (t1mi - t1j * MathSpecial::square(t1mi));
276 dt2dr1 = a2i * (t2mj - t2i * MathSpecial::square(t2mj));
277 dt2dr2 = a2j * (t2mi - t2j * MathSpecial::square(t2mi));
278 dt3dr1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
279 dt3dr2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
280
281 } else if (this->ialloy == 2) {
282
283 dt1dr1 = 0.0;
284 dt1dr2 = 0.0;
285 dt2dr1 = 0.0;
286 dt2dr2 = 0.0;
287 dt3dr1 = 0.0;
288 dt3dr2 = 0.0;
289
290 } else {
291
292 ai = 0.0;
293 if (!iszero(rho0[i]))
294 ai = drhoa0j * sij / rho0[i];
295 aj = 0.0;
296 if (!iszero(rho0[j]))
297 aj = drhoa0i * sij / rho0[j];
298
299 dt1dr1 = ai * (t1mj - t1i);
300 dt1dr2 = aj * (t1mi - t1j);
301 dt2dr1 = ai * (t2mj - t2i);
302 dt2dr2 = aj * (t2mi - t2j);
303 dt3dr1 = ai * (t3mj - t3i);
304 dt3dr2 = aj * (t3mi - t3j);
305 }
306
307 // Compute derivatives of total density wrt rij, sij and rij(3)
308 get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpi);
309 get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj);
310
311 drhodr1 = dgamma1[i] * drho0dr1 +
312 dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +
313 dt3dr1 * rho3[i] + t3i * drho3dr1) -
314 dgamma3[i] * (shpi[0] * dt1dr1 + shpi[1] * dt2dr1 + shpi[2] * dt3dr1);
315 drhodr2 = dgamma1[j] * drho0dr2 +
316 dgamma2[j] * (dt1dr2 * rho1[j] + t1j * drho1dr2 + dt2dr2 * rho2[j] + t2j * drho2dr2 +
317 dt3dr2 * rho3[j] + t3j * drho3dr2) -
318 dgamma3[j] * (shpj[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2);
319 for (m = 0; m < 3; m++) {
320 drhodrm1[m] = 0.0;
321 drhodrm2[m] = 0.0;
322 drhodrm1[m] = dgamma2[i] * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]);
323 drhodrm2[m] = dgamma2[j] * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]);
324 }
325
326 // Compute derivatives wrt sij, but only if necessary
327 if (!iszero(dscrfcn[fnoffset + jn])) {
328 drho0ds1 = rhoa0j;
329 drho0ds2 = rhoa0i;
330 a1 = 2.0 / rij;
331 drho1ds1 = a1 * rhoa1j * arg1i1;
332 drho1ds2 = a1 * rhoa1i * arg1j1;
333 a2 = 2.0 / rij2;
334 drho2ds1 = a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j;
335 drho2ds2 = a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arho2b[j] * rhoa2i;
336 a3 = 2.0 / rij3;
337 a3a = 6.0 / (5.0 * rij);
338 drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3;
339 drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3;
340
341 if (this->ialloy == 1) {
342 a1i = fdiv_zero(rhoa0j, tsq_ave[i][0]);
343 a1j = fdiv_zero(rhoa0i, tsq_ave[j][0]);
344 a2i = fdiv_zero(rhoa0j, tsq_ave[i][1]);
345 a2j = fdiv_zero(rhoa0i, tsq_ave[j][1]);
346 a3i = fdiv_zero(rhoa0j, tsq_ave[i][2]);
347 a3j = fdiv_zero(rhoa0i, tsq_ave[j][2]);
348
349 dt1ds1 = a1i * (t1mj - t1i * MathSpecial::square(t1mj));
350 dt1ds2 = a1j * (t1mi - t1j * MathSpecial::square(t1mi));
351 dt2ds1 = a2i * (t2mj - t2i * MathSpecial::square(t2mj));
352 dt2ds2 = a2j * (t2mi - t2j * MathSpecial::square(t2mi));
353 dt3ds1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
354 dt3ds2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
355
356 } else if (this->ialloy == 2) {
357
358 dt1ds1 = 0.0;
359 dt1ds2 = 0.0;
360 dt2ds1 = 0.0;
361 dt2ds2 = 0.0;
362 dt3ds1 = 0.0;
363 dt3ds2 = 0.0;
364
365 } else {
366
367 ai = 0.0;
368 if (!iszero(rho0[i]))
369 ai = rhoa0j / rho0[i];
370 aj = 0.0;
371 if (!iszero(rho0[j]))
372 aj = rhoa0i / rho0[j];
373
374 dt1ds1 = ai * (t1mj - t1i);
375 dt1ds2 = aj * (t1mi - t1j);
376 dt2ds1 = ai * (t2mj - t2i);
377 dt2ds2 = aj * (t2mi - t2j);
378 dt3ds1 = ai * (t3mj - t3i);
379 dt3ds2 = aj * (t3mi - t3j);
380 }
381
382 drhods1 = dgamma1[i] * drho0ds1 +
383 dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 +
384 dt3ds1 * rho3[i] + t3i * drho3ds1) -
385 dgamma3[i] * (shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);
386 drhods2 = dgamma1[j] * drho0ds2 +
387 dgamma2[j] * (dt1ds2 * rho1[j] + t1j * drho1ds2 + dt2ds2 * rho2[j] + t2j * drho2ds2 +
388 dt3ds2 * rho3[j] + t3j * drho3ds2) -
389 dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
390 }
391
392 // Compute derivatives of energy wrt rij, sij and rij[3]
393 dUdrij = phip * sij + frhop[i] * drhodr1 + frhop[j] * drhodr2;
394 dUdsij = 0.0;
395 if (!iszero(dscrfcn[fnoffset + jn])) {
396 dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2;
397 }
398 for (m = 0; m < 3; m++) {
399 dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
400 }
401 if (!isone(scaleij)) {
402 dUdrij *= scaleij;
403 dUdsij *= scaleij;
404 dUdrijm[0] *= scaleij;
405 dUdrijm[1] *= scaleij;
406 dUdrijm[2] *= scaleij;
407 }
408
409 // Add the part of the force due to dUdrij and dUdsij
410
411 force = dUdrij * recip + dUdsij * dscrfcn[fnoffset + jn];
412 for (m = 0; m < 3; m++) {
413 forcem = delij[m] * force + dUdrijm[m];
414 f[i][m] = f[i][m] + forcem;
415 f[j][m] = f[j][m] - forcem;
416 }
417
418 // Tabulate per-atom virial as symmetrized stress tensor
419
420 if (vflag_either) {
421 fi[0] = delij[0] * force + dUdrijm[0];
422 fi[1] = delij[1] * force + dUdrijm[1];
423 fi[2] = delij[2] * force + dUdrijm[2];
424 v[0] = -0.5 * (delij[0] * fi[0]);
425 v[1] = -0.5 * (delij[1] * fi[1]);
426 v[2] = -0.5 * (delij[2] * fi[2]);
427 v[3] = -0.25 * (delij[0] * fi[1] + delij[1] * fi[0]);
428 v[4] = -0.25 * (delij[0] * fi[2] + delij[2] * fi[0]);
429 v[5] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]);
430
431 if (vflag_global) {
432 for (m = 0; m < 6; m++) {
433 virial[m] += 2.0*v[m];
434 }
435 }
436 if (vflag_atom) {
437 for (m = 0; m < 6; m++) {
438 vatom[i][m] += v[m];
439 vatom[j][m] += v[m];
440 }
441 }
442 }
443
444 // Now compute forces on other atoms k due to change in sij
445
446 if (iszero(sij) || isone(sij)) continue; //: cont jn loop
447
448 double dxik(0), dyik(0), dzik(0);
449 double dxjk(0), dyjk(0), dzjk(0);
450
451 for (kn = 0; kn < numneigh_full; kn++) {
452 k = firstneigh_full[kn];
453 eltk = fmap[type[k]];
454 if (k != j && eltk >= 0) {
455 double xik, xjk, cikj, sikj, dfc, a;
456 double dCikj1, dCikj2;
457 double delc, rik2, rjk2;
458
459 sij = scrfcn[jn+fnoffset] * fcpair[jn+fnoffset];
460 const double Cmax = this->Cmax_meam[elti][eltj][eltk];
461 const double Cmin = this->Cmin_meam[elti][eltj][eltk];
462
463 dsij1 = 0.0;
464 dsij2 = 0.0;
465 if (!iszero(sij) && !isone(sij)) {
466 const double rbound = rij2 * this->ebound_meam[elti][eltj];
467 delc = Cmax - Cmin;
468 dxjk = x[k][0] - x[j][0];
469 dyjk = x[k][1] - x[j][1];
470 dzjk = x[k][2] - x[j][2];
471 rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
472 if (rjk2 <= rbound) {
473 dxik = x[k][0] - x[i][0];
474 dyik = x[k][1] - x[i][1];
475 dzik = x[k][2] - x[i][2];
476 rik2 = dxik * dxik + dyik * dyik + dzik * dzik;
477 if (rik2 <= rbound) {
478 xik = rik2 / rij2;
479 xjk = rjk2 / rij2;
480 a = 1 - (xik - xjk) * (xik - xjk);
481 if (!iszero(a)) {
482 cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
483 if (cikj >= Cmin && cikj <= Cmax) {
484 cikj = (cikj - Cmin) / delc;
485 sikj = dfcut(cikj, dfc);
486 dCfunc2(rij2, rik2, rjk2, dCikj1, dCikj2);
487 a = sij / delc * dfc / sikj;
488 dsij1 = a * dCikj1;
489 dsij2 = a * dCikj2;
490 }
491 }
492 }
493 }
494 }
495
496 if (!iszero(dsij1) || !iszero(dsij2)) {
497 force1 = dUdsij * dsij1;
498 force2 = dUdsij * dsij2;
499
500 f[i][0] += force1 * dxik;
501 f[i][1] += force1 * dyik;
502 f[i][2] += force1 * dzik;
503 f[j][0] += force2 * dxjk;
504 f[j][1] += force2 * dyjk;
505 f[j][2] += force2 * dzjk;
506 f[k][0] -= force1 * dxik + force2 * dxjk;
507 f[k][1] -= force1 * dyik + force2 * dyjk;
508 f[k][2] -= force1 * dzik + force2 * dzjk;
509
510 // Tabulate per-atom virial as symmetrized stress tensor
511
512 if (vflag_either) {
513 fi[0] = force1 * dxik;
514 fi[1] = force1 * dyik;
515 fi[2] = force1 * dzik;
516 fj[0] = force2 * dxjk;
517 fj[1] = force2 * dyjk;
518 fj[2] = force2 * dzjk;
519 v[0] = -third * (dxik * fi[0] + dxjk * fj[0]);
520 v[1] = -third * (dyik * fi[1] + dyjk * fj[1]);
521 v[2] = -third * (dzik * fi[2] + dzjk * fj[2]);
522 v[3] = -sixth * (dxik * fi[1] + dxjk * fj[1] + dyik * fi[0] + dyjk * fj[0]);
523 v[4] = -sixth * (dxik * fi[2] + dxjk * fj[2] + dzik * fi[0] + dzjk * fj[0]);
524 v[5] = -sixth * (dyik * fi[2] + dyjk * fj[2] + dzik * fi[1] + dzjk * fj[1]);
525
526 if (vflag_global) {
527 for (m = 0; m < 6; m++) {
528 virial[m] += 3.0*v[m];
529 }
530 }
531
532 if (vflag_atom) {
533 for (m = 0; m < 6; m++) {
534 vatom[i][m] += v[m];
535 vatom[j][m] += v[m];
536 vatom[k][m] += v[m];
537 }
538 }
539 }
540 }
541 }
542 // end of k loop
543 }
544 }
545 }
546 // end of j loop
547 }
548 }
549