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 #include "memory.h"
18
19 #include <cmath>
20
21 using namespace LAMMPS_NS;
22
23 void
meam_dens_setup(int atom_nmax,int nall,int n_neigh)24 MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
25 {
26 int i, j;
27
28 // grow local arrays if necessary
29
30 if (atom_nmax > nmax) {
31 memory->destroy(rho);
32 memory->destroy(rho0);
33 memory->destroy(rho1);
34 memory->destroy(rho2);
35 memory->destroy(rho3);
36 memory->destroy(frhop);
37 memory->destroy(gamma);
38 memory->destroy(dgamma1);
39 memory->destroy(dgamma2);
40 memory->destroy(dgamma3);
41 memory->destroy(arho2b);
42 memory->destroy(arho1);
43 memory->destroy(arho2);
44 memory->destroy(arho3);
45 memory->destroy(arho3b);
46 memory->destroy(t_ave);
47 memory->destroy(tsq_ave);
48
49 nmax = atom_nmax;
50
51 memory->create(rho, nmax, "pair:rho");
52 memory->create(rho0, nmax, "pair:rho0");
53 memory->create(rho1, nmax, "pair:rho1");
54 memory->create(rho2, nmax, "pair:rho2");
55 memory->create(rho3, nmax, "pair:rho3");
56 memory->create(frhop, nmax, "pair:frhop");
57 memory->create(gamma, nmax, "pair:gamma");
58 memory->create(dgamma1, nmax, "pair:dgamma1");
59 memory->create(dgamma2, nmax, "pair:dgamma2");
60 memory->create(dgamma3, nmax, "pair:dgamma3");
61 memory->create(arho2b, nmax, "pair:arho2b");
62 memory->create(arho1, nmax, 3, "pair:arho1");
63 memory->create(arho2, nmax, 6, "pair:arho2");
64 memory->create(arho3, nmax, 10, "pair:arho3");
65 memory->create(arho3b, nmax, 3, "pair:arho3b");
66 memory->create(t_ave, nmax, 3, "pair:t_ave");
67 memory->create(tsq_ave, nmax, 3, "pair:tsq_ave");
68 }
69
70 if (n_neigh > maxneigh) {
71 memory->destroy(scrfcn);
72 memory->destroy(dscrfcn);
73 memory->destroy(fcpair);
74 maxneigh = n_neigh;
75 memory->create(scrfcn, maxneigh, "pair:scrfcn");
76 memory->create(dscrfcn, maxneigh, "pair:dscrfcn");
77 memory->create(fcpair, maxneigh, "pair:fcpair");
78 }
79
80 // zero out local arrays
81
82 for (i = 0; i < nall; i++) {
83 rho0[i] = 0.0;
84 arho2b[i] = 0.0;
85 arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
86 for (j = 0; j < 6; j++)
87 arho2[i][j] = 0.0;
88 for (j = 0; j < 10; j++)
89 arho3[i][j] = 0.0;
90 arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
91 t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
92 tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
93 }
94 }
95
96 void
meam_dens_init(int i,int ntype,int * type,int * fmap,double ** x,int numneigh,int * firstneigh,int numneigh_full,int * firstneigh_full,int fnoffset)97 MEAM::meam_dens_init(int i, int ntype, int* type, int* fmap, double** x,
98 int numneigh, int* firstneigh,
99 int numneigh_full, int* firstneigh_full, int fnoffset)
100 {
101 // Compute screening function and derivatives
102 getscreen(i, &scrfcn[fnoffset], &dscrfcn[fnoffset], &fcpair[fnoffset], x, numneigh, firstneigh,
103 numneigh_full, firstneigh_full, ntype, type, fmap);
104
105 // Calculate intermediate density terms to be communicated
106 calc_rho1(i, ntype, type, fmap, x, numneigh, firstneigh, &scrfcn[fnoffset], &fcpair[fnoffset]);
107 }
108
109 // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
110
111 void
getscreen(int i,double * scrfcn,double * dscrfcn,double * fcpair,double ** x,int numneigh,int * firstneigh,int numneigh_full,int * firstneigh_full,int,int * type,int * fmap)112 MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
113 int* firstneigh, int numneigh_full, int* firstneigh_full, int /*ntype*/, int* type, int* fmap)
114 {
115 int jn, j, kn, k;
116 int elti, eltj, eltk;
117 double xitmp, yitmp, zitmp, delxij, delyij, delzij, rij2, rij;
118 double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/;
119 double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/;
120 double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj;
121 double Cmin, Cmax, delc, /*ebound,*/ a, coef1, coef2;
122 double dCikj;
123 double rnorm, fc, dfc, drinv;
124
125 drinv = 1.0 / this->delr_meam;
126 elti = fmap[type[i]];
127 if (elti < 0) return;
128
129 xitmp = x[i][0];
130 yitmp = x[i][1];
131 zitmp = x[i][2];
132
133 for (jn = 0; jn < numneigh; jn++) {
134 j = firstneigh[jn];
135
136 eltj = fmap[type[j]];
137 if (eltj < 0) continue;
138
139 // First compute screening function itself, sij
140 xjtmp = x[j][0];
141 yjtmp = x[j][1];
142 zjtmp = x[j][2];
143 delxij = xjtmp - xitmp;
144 delyij = yjtmp - yitmp;
145 delzij = zjtmp - zitmp;
146 rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
147
148 if (rij2 > this->cutforcesq) {
149 dscrfcn[jn] = 0.0;
150 scrfcn[jn] = 0.0;
151 fcpair[jn] = 0.0;
152 continue;
153 }
154
155 const double rbound = this->ebound_meam[elti][eltj] * rij2;
156 rij = sqrt(rij2);
157 rnorm = (this->cutforce - rij) * drinv;
158 sij = 1.0;
159
160 // if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
161 for (kn = 0; kn < numneigh_full; kn++) {
162 k = firstneigh_full[kn];
163 if (k == j) continue;
164 eltk = fmap[type[k]];
165 if (eltk < 0) continue;
166
167 xktmp = x[k][0];
168 yktmp = x[k][1];
169 zktmp = x[k][2];
170
171 delxjk = xktmp - xjtmp;
172 delyjk = yktmp - yjtmp;
173 delzjk = zktmp - zjtmp;
174 rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
175 if (rjk2 > rbound) continue;
176
177 delxik = xktmp - xitmp;
178 delyik = yktmp - yitmp;
179 delzik = zktmp - zitmp;
180 rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
181 if (rik2 > rbound) continue;
182
183 xik = rik2 / rij2;
184 xjk = rjk2 / rij2;
185 a = 1 - (xik - xjk) * (xik - xjk);
186 // if a < 0, then ellipse equation doesn't describe this case and
187 // atom k can't possibly screen i-j
188 if (a <= 0.0) continue;
189
190 cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
191 Cmax = this->Cmax_meam[elti][eltj][eltk];
192 Cmin = this->Cmin_meam[elti][eltj][eltk];
193 if (cikj >= Cmax) continue;
194 // note that cikj may be slightly negative (within numerical
195 // tolerance) if atoms are colinear, so don't reject that case here
196 // (other negative cikj cases were handled by the test on "a" above)
197 else if (cikj <= Cmin) {
198 sij = 0.0;
199 break;
200 } else {
201 delc = Cmax - Cmin;
202 cikj = (cikj - Cmin) / delc;
203 sikj = fcut(cikj);
204 }
205 sij *= sikj;
206 }
207
208 fc = dfcut(rnorm, dfc);
209 fcij = fc;
210 dfcij = dfc * drinv;
211
212 // Now compute derivatives
213 dscrfcn[jn] = 0.0;
214 sfcij = sij * fcij;
215 if (!iszero(sfcij) && !isone(sfcij)) {
216 for (kn = 0; kn < numneigh_full; kn++) {
217 k = firstneigh_full[kn];
218 if (k == j) continue;
219 eltk = fmap[type[k]];
220 if (eltk < 0) continue;
221
222 delxjk = x[k][0] - xjtmp;
223 delyjk = x[k][1] - yjtmp;
224 delzjk = x[k][2] - zjtmp;
225 rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
226 if (rjk2 > rbound) continue;
227
228 delxik = x[k][0] - xitmp;
229 delyik = x[k][1] - yitmp;
230 delzik = x[k][2] - zitmp;
231 rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
232 if (rik2 > rbound) continue;
233
234 xik = rik2 / rij2;
235 xjk = rjk2 / rij2;
236 a = 1 - (xik - xjk) * (xik - xjk);
237 // if a < 0, then ellipse equation doesn't describe this case and
238 // atom k can't possibly screen i-j
239 if (a <= 0.0) continue;
240
241 cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
242 Cmax = this->Cmax_meam[elti][eltj][eltk];
243 Cmin = this->Cmin_meam[elti][eltj][eltk];
244 if (cikj >= Cmax) {
245 continue;
246 // Note that cikj may be slightly negative (within numerical
247 // tolerance) if atoms are colinear, so don't reject that case
248 // here
249 // (other negative cikj cases were handled by the test on "a"
250 // above)
251 // Note that we never have 0<cikj<Cmin here, else sij=0
252 // (rejected above)
253 } else {
254 delc = Cmax - Cmin;
255 cikj = (cikj - Cmin) / delc;
256 sikj = dfcut(cikj, dfikj);
257 coef1 = dfikj / (delc * sikj);
258 dCikj = dCfunc(rij2, rik2, rjk2);
259 dscrfcn[jn] = dscrfcn[jn] + coef1 * dCikj;
260 }
261 }
262 coef1 = sfcij;
263 coef2 = sij * dfcij / rij;
264 dscrfcn[jn] = dscrfcn[jn] * coef1 - coef2;
265 }
266
267 scrfcn[jn] = sij;
268 fcpair[jn] = fcij;
269 }
270 }
271
272 // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
273
274 void
calc_rho1(int i,int,int * type,int * fmap,double ** x,int numneigh,int * firstneigh,double * scrfcn,double * fcpair)275 MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
276 double* scrfcn, double* fcpair)
277 {
278 int jn, j, m, n, p, elti, eltj;
279 int nv2, nv3;
280 double xtmp, ytmp, ztmp, delij[3], rij2, rij, sij;
281 double ai, aj, rhoa0j, rhoa1j, rhoa2j, rhoa3j, A1j, A2j, A3j;
282 // double G,Gbar,gam,shp[3+1];
283 double ro0i, ro0j;
284 double rhoa0i, rhoa1i, rhoa2i, rhoa3i, A1i, A2i, A3i;
285
286 elti = fmap[type[i]];
287 xtmp = x[i][0];
288 ytmp = x[i][1];
289 ztmp = x[i][2];
290 for (jn = 0; jn < numneigh; jn++) {
291 if (!iszero(scrfcn[jn])) {
292 j = firstneigh[jn];
293 sij = scrfcn[jn] * fcpair[jn];
294 delij[0] = x[j][0] - xtmp;
295 delij[1] = x[j][1] - ytmp;
296 delij[2] = x[j][2] - ztmp;
297 rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
298 if (rij2 < this->cutforcesq) {
299 eltj = fmap[type[j]];
300 rij = sqrt(rij2);
301 ai = rij / this->re_meam[elti][elti] - 1.0;
302 aj = rij / this->re_meam[eltj][eltj] - 1.0;
303 ro0i = this->rho0_meam[elti];
304 ro0j = this->rho0_meam[eltj];
305 rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj) * sij;
306 rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj) * sij;
307 rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj) * sij;
308 rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj) * sij;
309 rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai) * sij;
310 rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai) * sij;
311 rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai) * sij;
312 rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai) * sij;
313 if (this->ialloy == 1) {
314 rhoa1j = rhoa1j * this->t1_meam[eltj];
315 rhoa2j = rhoa2j * this->t2_meam[eltj];
316 rhoa3j = rhoa3j * this->t3_meam[eltj];
317 rhoa1i = rhoa1i * this->t1_meam[elti];
318 rhoa2i = rhoa2i * this->t2_meam[elti];
319 rhoa3i = rhoa3i * this->t3_meam[elti];
320 }
321 rho0[i] = rho0[i] + rhoa0j;
322 rho0[j] = rho0[j] + rhoa0i;
323 // For ialloy = 2, use single-element value (not average)
324 if (this->ialloy != 2) {
325 t_ave[i][0] = t_ave[i][0] + this->t1_meam[eltj] * rhoa0j;
326 t_ave[i][1] = t_ave[i][1] + this->t2_meam[eltj] * rhoa0j;
327 t_ave[i][2] = t_ave[i][2] + this->t3_meam[eltj] * rhoa0j;
328 t_ave[j][0] = t_ave[j][0] + this->t1_meam[elti] * rhoa0i;
329 t_ave[j][1] = t_ave[j][1] + this->t2_meam[elti] * rhoa0i;
330 t_ave[j][2] = t_ave[j][2] + this->t3_meam[elti] * rhoa0i;
331 }
332 if (this->ialloy == 1) {
333 tsq_ave[i][0] = tsq_ave[i][0] + this->t1_meam[eltj] * this->t1_meam[eltj] * rhoa0j;
334 tsq_ave[i][1] = tsq_ave[i][1] + this->t2_meam[eltj] * this->t2_meam[eltj] * rhoa0j;
335 tsq_ave[i][2] = tsq_ave[i][2] + this->t3_meam[eltj] * this->t3_meam[eltj] * rhoa0j;
336 tsq_ave[j][0] = tsq_ave[j][0] + this->t1_meam[elti] * this->t1_meam[elti] * rhoa0i;
337 tsq_ave[j][1] = tsq_ave[j][1] + this->t2_meam[elti] * this->t2_meam[elti] * rhoa0i;
338 tsq_ave[j][2] = tsq_ave[j][2] + this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i;
339 }
340 arho2b[i] = arho2b[i] + rhoa2j;
341 arho2b[j] = arho2b[j] + rhoa2i;
342
343 A1j = rhoa1j / rij;
344 A2j = rhoa2j / rij2;
345 A3j = rhoa3j / (rij2 * rij);
346 A1i = rhoa1i / rij;
347 A2i = rhoa2i / rij2;
348 A3i = rhoa3i / (rij2 * rij);
349 nv2 = 0;
350 nv3 = 0;
351 for (m = 0; m < 3; m++) {
352 arho1[i][m] = arho1[i][m] + A1j * delij[m];
353 arho1[j][m] = arho1[j][m] - A1i * delij[m];
354 arho3b[i][m] = arho3b[i][m] + rhoa3j * delij[m] / rij;
355 arho3b[j][m] = arho3b[j][m] - rhoa3i * delij[m] / rij;
356 for (n = m; n < 3; n++) {
357 arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n];
358 arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n];
359 nv2 = nv2 + 1;
360 for (p = n; p < 3; p++) {
361 arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p];
362 arho3[j][nv3] = arho3[j][nv3] - A3i * delij[m] * delij[n] * delij[p];
363 nv3 = nv3 + 1;
364 }
365 }
366 }
367 }
368 }
369 }
370 }
371
372