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