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 #include <algorithm>
21
22 using namespace LAMMPS_NS;
23
meam_setup_done(double * cutmax)24 void MEAM::meam_setup_done(double* cutmax)
25 {
26 int nv2, nv3, m, n, p;
27
28 // Force cutoff
29 this->cutforce = this->rc_meam;
30 this->cutforcesq = this->cutforce * this->cutforce;
31
32 // Pass cutoff back to calling program
33 *cutmax = this->cutforce;
34
35 // Augment t1 term
36 for (int i = 0; i < maxelt; i++)
37 this->t1_meam[i] = this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i];
38
39 // Compute off-diagonal alloy parameters
40 alloyparams();
41
42 // indices and factors for Voight notation
43 nv2 = 0;
44 nv3 = 0;
45 for (m = 0; m < 3; m++) {
46 for (n = m; n < 3; n++) {
47 this->vind2D[m][n] = nv2;
48 this->vind2D[n][m] = nv2;
49 nv2 = nv2 + 1;
50 for (p = n; p < 3; p++) {
51 this->vind3D[m][n][p] = nv3;
52 this->vind3D[m][p][n] = nv3;
53 this->vind3D[n][m][p] = nv3;
54 this->vind3D[n][p][m] = nv3;
55 this->vind3D[p][m][n] = nv3;
56 this->vind3D[p][n][m] = nv3;
57 nv3 = nv3 + 1;
58 }
59 }
60 }
61
62 this->v2D[0] = 1;
63 this->v2D[1] = 2;
64 this->v2D[2] = 2;
65 this->v2D[3] = 1;
66 this->v2D[4] = 2;
67 this->v2D[5] = 1;
68
69 this->v3D[0] = 1;
70 this->v3D[1] = 3;
71 this->v3D[2] = 3;
72 this->v3D[3] = 3;
73 this->v3D[4] = 6;
74 this->v3D[5] = 3;
75 this->v3D[6] = 1;
76 this->v3D[7] = 3;
77 this->v3D[8] = 3;
78 this->v3D[9] = 1;
79
80 nv2 = 0;
81 for (m = 0; m < this->neltypes; m++) {
82 for (n = m; n < this->neltypes; n++) {
83 this->eltind[m][n] = nv2;
84 this->eltind[n][m] = nv2;
85 nv2 = nv2 + 1;
86 }
87 }
88
89 // Compute background densities for reference structure
90 compute_reference_density();
91
92 // Compute pair potentials and setup arrays for interpolation
93 this->nr = 1000;
94 this->dr = 1.1 * this->rc_meam / this->nr;
95 compute_pair_meam();
96 }
97
98 // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
99 // Fill off-diagonal alloy parameters
alloyparams()100 void MEAM::alloyparams()
101 {
102
103 int i, j, k;
104 double eb;
105
106 // Loop over pairs
107 for (i = 0; i < this->neltypes; i++) {
108 for (j = 0; j < this->neltypes; j++) {
109 // Treat off-diagonal pairs
110 // If i>j, set all equal to i<j case (which has already been set,
111 // here or in the input file)
112 if (i > j) {
113 this->re_meam[i][j] = this->re_meam[j][i];
114 this->Ec_meam[i][j] = this->Ec_meam[j][i];
115 this->alpha_meam[i][j] = this->alpha_meam[j][i];
116 this->lattce_meam[i][j] = this->lattce_meam[j][i];
117 this->nn2_meam[i][j] = this->nn2_meam[j][i];
118 // theta for lin,tri,zig references
119 this->stheta_meam[i][j] = this->stheta_meam[j][i];
120 this->ctheta_meam[i][j] = this->ctheta_meam[j][i];
121 // If i<j and term is unset, use default values (e.g. mean of i-i and
122 // j-j)
123 } else if (j > i) {
124 if (iszero(this->Ec_meam[i][j])) {
125 if (this->lattce_meam[i][j] == L12)
126 this->Ec_meam[i][j] =
127 (3 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 4.0 - this->delta_meam[i][j];
128 else if (this->lattce_meam[i][j] == C11) {
129 if (this->lattce_meam[i][i] == DIA)
130 this->Ec_meam[i][j] =
131 (2 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
132 else
133 this->Ec_meam[i][j] =
134 (this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
135 } else
136 this->Ec_meam[i][j] = (this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 - this->delta_meam[i][j];
137 }
138 if (iszero(this->alpha_meam[i][j]))
139 this->alpha_meam[i][j] = (this->alpha_meam[i][i] + this->alpha_meam[j][j]) / 2.0;
140 if (iszero(this->re_meam[i][j]))
141 this->re_meam[i][j] = (this->re_meam[i][i] + this->re_meam[j][j]) / 2.0;
142 }
143 }
144 }
145
146 // Cmin[i][k][j] is symmetric in i-j, but not k. For all triplets
147 // where i>j, set equal to the i<j element. Likewise for Cmax.
148 for (i = 1; i < this->neltypes; i++) {
149 for (j = 0; j < i; j++) {
150 for (k = 0; k < this->neltypes; k++) {
151 this->Cmin_meam[i][j][k] = this->Cmin_meam[j][i][k];
152 this->Cmax_meam[i][j][k] = this->Cmax_meam[j][i][k];
153 }
154 }
155 }
156
157 // ebound gives the squared distance such that, for rik2 or rjk2>ebound,
158 // atom k definitely lies outside the screening function ellipse (so
159 // there is no need to calculate its effects). Here, compute it for all
160 // triplets [i][j][k] so that ebound[i][j] is the maximized over k
161 for (i = 0; i < this->neltypes; i++) {
162 for (j = 0; j < this->neltypes; j++) {
163 for (k = 0; k < this->neltypes; k++) {
164 eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / (4.0 * (this->Cmax_meam[i][j][k] - 1.0));
165 this->ebound_meam[i][j] = std::max(this->ebound_meam[i][j], eb);
166 }
167 }
168 }
169 }
170
171 //-----------------------------------------------------------------------
172 // compute MEAM pair potential for each pair of element types
173 //
174
compute_pair_meam()175 void MEAM::compute_pair_meam()
176 {
177 double r;
178 int j, a, b, nv2;
179 double astar, frac, phizbl;
180 int Z1, Z2;
181 double arat, rarat, scrn, scrn2;
182 double phiaa, phibb /*unused:,phitmp*/;
183 double C, s111, s112, s221, S11, S22;
184
185 // check for previously allocated arrays and free them
186 if (this->phir != nullptr)
187 memory->destroy(this->phir);
188 if (this->phirar != nullptr)
189 memory->destroy(this->phirar);
190 if (this->phirar1 != nullptr)
191 memory->destroy(this->phirar1);
192 if (this->phirar2 != nullptr)
193 memory->destroy(this->phirar2);
194 if (this->phirar3 != nullptr)
195 memory->destroy(this->phirar3);
196 if (this->phirar4 != nullptr)
197 memory->destroy(this->phirar4);
198 if (this->phirar5 != nullptr)
199 memory->destroy(this->phirar5);
200 if (this->phirar6 != nullptr)
201 memory->destroy(this->phirar6);
202
203 // allocate memory for array that defines the potential
204 memory->create(this->phir, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phir");
205
206 // allocate coeff memory
207
208 memory->create(this->phirar, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar");
209 memory->create(this->phirar1, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar1");
210 memory->create(this->phirar2, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar2");
211 memory->create(this->phirar3, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar3");
212 memory->create(this->phirar4, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar4");
213 memory->create(this->phirar5, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar5");
214 memory->create(this->phirar6, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar6");
215
216 // loop over pairs of element types
217 nv2 = 0;
218 for (a = 0; a < this->neltypes; a++) {
219 for (b = a; b < this->neltypes; b++) {
220 // loop over r values and compute
221 for (j = 0; j < this->nr; j++) {
222 r = j * this->dr;
223
224 this->phir[nv2][j] = phi_meam(r, a, b);
225
226 // if using second-nearest neighbor, solve recursive problem
227 // (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
228 if (this->nn2_meam[a][b] == 1) {
229 Z1 = get_Zij(this->lattce_meam[a][b]);
230 Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
231 this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn);
232
233 // The B1, B2, and L12 cases with NN2 have a trick to them; we need to
234 // compute the contributions from second nearest neighbors, like a-a
235 // pairs, but need to include NN2 contributions to those pairs as
236 // well.
237 if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
238 this->lattce_meam[a][b] == L12 || this->lattce_meam[a][b] == DIA) {
239 rarat = r * arat;
240
241 // phi_aa
242 phiaa = phi_meam(rarat, a, a);
243 Z1 = get_Zij(this->lattce_meam[a][a]);
244 Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
245 this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);
246 phiaa+= phi_meam_series(scrn, Z1, Z2, a, a, rarat, arat);
247
248 // phi_bb
249 phibb = phi_meam(rarat, b, b);
250 Z1 = get_Zij(this->lattce_meam[b][b]);
251 Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b],
252 this->Cmax_meam[b][b][b], this->stheta_meam[b][b], arat, scrn);
253 phibb+= phi_meam_series(scrn, Z1, Z2, b, b, rarat, arat);
254
255 if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
256 this->lattce_meam[a][b] == DIA) {
257 // Add contributions to the B1 or B2 potential
258 Z1 = get_Zij(this->lattce_meam[a][b]);
259 Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
260 this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn);
261 this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
262 Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[b][b][a],
263 this->Cmax_meam[b][b][a], this->stheta_meam[a][b], arat, scrn2);
264
265 this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
266
267 } else if (this->lattce_meam[a][b] == L12) {
268 // The L12 case has one last trick; we have to be careful to
269 // compute
270 // the correct screening between 2nd-neighbor pairs. 1-1
271 // second-neighbor pairs are screened by 2 type 1 atoms and
272 // two type
273 // 2 atoms. 2-2 second-neighbor pairs are screened by 4 type
274 // 1
275 // atoms.
276 C = 1.0;
277 get_sijk(C, a, a, a, &s111);
278 get_sijk(C, a, a, b, &s112);
279 get_sijk(C, b, b, a, &s221);
280 S11 = s111 * s111 * s112 * s112;
281 S22 = pow(s221, 4);
282 this->phir[nv2][j] = this->phir[nv2][j] - 0.75 * S11 * phiaa - 0.25 * S22 * phibb;
283 }
284
285 } else {
286 this->phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat);
287 }
288 }
289
290 // For Zbl potential:
291 // if astar <= -3
292 // potential is zbl potential
293 // else if -3 < astar < -1
294 // potential is linear combination with zbl potential
295 // endif
296 if (this->zbl_meam[a][b] == 1) {
297 astar = this->alpha_meam[a][b] * (r / this->re_meam[a][b] - 1.0);
298 if (astar <= -3.0)
299 this->phir[nv2][j] = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
300 else if (astar > -3.0 && astar < -1.0) {
301 frac = fcut(1 - (astar + 1.0) / (-3.0 + 1.0));
302 phizbl = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
303 this->phir[nv2][j] = frac * this->phir[nv2][j] + (1 - frac) * phizbl;
304 }
305 }
306 }
307
308 // call interpolation
309 interpolate_meam(nv2);
310
311 nv2 = nv2 + 1;
312 }
313 }
314 }
315
316 //----------------------------------------------------------------------c
317 // Compute MEAM pair potential for distance r, element types a and b
318 //
phi_meam(double r,int a,int b)319 double MEAM::phi_meam(double r, int a, int b)
320 {
321 /*unused:double a1,a2,a12;*/
322 double t11av, t21av, t31av, t12av, t22av, t32av;
323 double G1, G2, s1[3], s2[3], rho0_1, rho0_2;
324 double Gam1, Gam2, Z1, Z2;
325 double rhobar1, rhobar2, F1, F2, dF;
326 double rho01, rho11, rho21, rho31;
327 double rho02, rho12, rho22, rho32;
328 double scalfac, phiaa, phibb;
329 double Eu;
330 double arat, scrn, scrn2;
331 int Z12, errorflag;
332 int Z1nn, Z2nn;
333 lattice_t latta /*unused:,lattb*/;
334 double rho_bkgd1, rho_bkgd2;
335 double b11s, b22s;
336
337 double phi_m = 0.0;
338
339 // Equation numbers below refer to:
340 // I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615
341
342 // get number of neighbors in the reference structure
343 // Nref[i][j] = # of i's neighbors of type j
344 Z1 = get_Zij(this->lattce_meam[a][a]);
345 Z2 = get_Zij(this->lattce_meam[b][b]);
346 Z12 = get_Zij(this->lattce_meam[a][b]);
347
348 get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32);
349
350 // if densities are too small, numerical problems may result; just return zero
351 if (rho01 <= 1e-14 && rho02 <= 1e-14)
352 return 0.0;
353
354 // calculate average weighting factors for the reference structure
355 if (this->lattce_meam[a][b] == C11) {
356 if (this->ialloy == 2) {
357 t11av = this->t1_meam[a];
358 t12av = this->t1_meam[b];
359 t21av = this->t2_meam[a];
360 t22av = this->t2_meam[b];
361 t31av = this->t3_meam[a];
362 t32av = this->t3_meam[b];
363 } else {
364 scalfac = 1.0 / (rho01 + rho02);
365 t11av = scalfac * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02);
366 t12av = t11av;
367 t21av = scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02);
368 t22av = t21av;
369 t31av = scalfac * (this->t3_meam[a] * rho01 + this->t3_meam[b] * rho02);
370 t32av = t31av;
371 }
372 } else {
373 // average weighting factors for the reference structure, eqn. I.8
374 get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, this->t1_meam[a], this->t2_meam[a],
375 this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b,
376 this->lattce_meam[a][b]);
377 }
378
379 // for c11b structure, calculate background electron densities
380 if (this->lattce_meam[a][b] == C11) {
381 latta = this->lattce_meam[a][a];
382 if (latta == DIA) {
383 rhobar1 = MathSpecial::square((Z12 / 2) * (rho02 + rho01))
384 + t11av * MathSpecial::square(rho12 - rho11)
385 + t21av / 6.0 * MathSpecial::square(rho22 + rho21)
386 + 121.0 / 40.0 * t31av * MathSpecial::square(rho32 - rho31);
387 rhobar1 = sqrt(rhobar1);
388 rhobar2 = MathSpecial::square(Z12 * rho01) + 2.0 / 3.0 * t21av * MathSpecial::square(rho21);
389 rhobar2 = sqrt(rhobar2);
390 } else {
391 rhobar2 = MathSpecial::square((Z12 / 2) * (rho01 + rho02))
392 + t12av * MathSpecial::square(rho11 - rho12)
393 + t22av / 6.0 * MathSpecial::square(rho21 + rho22)
394 + 121.0 / 40.0 * t32av * MathSpecial::square(rho31 - rho32);
395 rhobar2 = sqrt(rhobar2);
396 rhobar1 = MathSpecial::square(Z12 * rho02) + 2.0 / 3.0 * t22av * MathSpecial::square(rho22);
397 rhobar1 = sqrt(rhobar1);
398 }
399 } else {
400 // for other structures, use formalism developed in Huang's paper
401 //
402 // composition-dependent scaling, equation I.7
403 // If using mixing rule for t, apply to reference structure; else
404 // use precomputed values
405 if (this->mix_ref_t == 1) {
406 if (this->ibar_meam[a] <= 0)
407 G1 = 1.0;
408 else {
409 get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], s1);
410 Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1);
411 G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
412 }
413 if (this->ibar_meam[b] <= 0)
414 G2 = 1.0;
415 else {
416 get_shpfcn(this->lattce_meam[b][b], this->stheta_meam[b][b], this->ctheta_meam[b][b], s2);
417 Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2);
418 G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
419 }
420 rho0_1 = this->rho0_meam[a] * Z1 * G1;
421 rho0_2 = this->rho0_meam[b] * Z2 * G2;
422 }
423 Gam1 = (t11av * rho11 + t21av * rho21 + t31av * rho31);
424 if (rho01 < 1.0e-14)
425 Gam1 = 0.0;
426 else
427 Gam1 = Gam1 / (rho01 * rho01);
428
429 Gam2 = (t12av * rho12 + t22av * rho22 + t32av * rho32);
430 if (rho02 < 1.0e-14)
431 Gam2 = 0.0;
432 else
433 Gam2 = Gam2 / (rho02 * rho02);
434
435 G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
436 G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
437 if (this->mix_ref_t == 1) {
438 rho_bkgd1 = rho0_1;
439 rho_bkgd2 = rho0_2;
440 } else {
441 if (this->bkgd_dyn == 1) {
442 rho_bkgd1 = this->rho0_meam[a] * Z1;
443 rho_bkgd2 = this->rho0_meam[b] * Z2;
444 } else {
445 rho_bkgd1 = this->rho_ref_meam[a];
446 rho_bkgd2 = this->rho_ref_meam[b];
447 }
448 }
449 rhobar1 = rho01 / rho_bkgd1 * G1;
450 rhobar2 = rho02 / rho_bkgd2 * G2;
451 }
452
453 // compute embedding functions, eqn I.5
454
455 F1 = embedding(this->A_meam[a], this->Ec_meam[a][a], rhobar1, dF);
456 F2 = embedding(this->A_meam[b], this->Ec_meam[b][b], rhobar2, dF);
457
458
459 // compute Rose function, I.16
460 Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], this->Ec_meam[a][b], this->repuls_meam[a][b],
461 this->attrac_meam[a][b], this->erose_form);
462
463 // calculate the pair energy
464 if (this->lattce_meam[a][b] == C11) {
465 latta = this->lattce_meam[a][a];
466 if (latta == DIA) {
467 phiaa = phi_meam(r, a, a);
468 phi_m = (3 * Eu - F2 - 2 * F1 - 5 * phiaa) / Z12;
469 } else {
470 phibb = phi_meam(r, b, b);
471 phi_m = (3 * Eu - F1 - 2 * F2 - 5 * phibb) / Z12;
472 }
473 } else if (this->lattce_meam[a][b] == L12) {
474 phiaa = phi_meam(r, a, a);
475 // account for second neighbor a-a potential here...
476 Z1nn = get_Zij(this->lattce_meam[a][a]);
477 Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
478 this->Cmax_meam[a][a][a], this->stheta_meam[a][b], arat, scrn);
479
480
481 phiaa += phi_meam_series(scrn, Z1nn, Z2nn, a, a, r, arat);
482 phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;
483
484 } else if (this->lattce_meam[a][b] == CH4) {
485 phi_m = (5 * Eu - F1 - 4*F2)/4;
486
487 } else if (this->lattce_meam[a][b] == ZIG) {
488 if (a==b) {
489 phi_m = (2 * Eu - F1 - F2) / Z12;
490 } else{
491 Z1 = get_Zij(this->lattce_meam[a][b]);
492 Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);
493 b11s = -Z2/Z1*scrn;
494 Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], scrn2);
495 b22s = -Z2/Z1*scrn2;
496
497 phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
498 phibb = phi_meam(2.0*this->stheta_meam[a][b]*r, b, b);
499 phi_m = (2.0*Eu - F1 - F2 + phiaa*b11s + phibb*b22s) / Z12;
500 }
501
502 } else if (this->lattce_meam[a][b] == TRI) {
503 if (a==b) {
504 phi_m = (3.0*Eu - 2.0*F1 - F2) / Z12;
505 } else {
506 Z1 = get_Zij(this->lattce_meam[a][b]);
507 Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);
508 b11s = -Z2/Z1*scrn;
509 phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
510 phi_m = (3.0*Eu - 2.0*F1 - F2 + phiaa*b11s) / Z12;
511 }
512
513 } else {
514 // potential is computed from Rose function and embedding energy
515 phi_m = (2 * Eu - F1 - F2) / Z12;
516 }
517
518 // if r = 0, just return 0
519 if (iszero(r)) {
520 phi_m = 0.0;
521 }
522
523 return phi_m;
524 }
525
526 //----------------------------------------------------------------------c
527 // Compute 2NN series terms for phi
528 // To avoid nan values of phir due to rapid decrease of b2nn^n or/and
529 // argument of phi_meam, i.e. r*arat^n, in some cases (3NN dia with low Cmin value)
530 //
phi_meam_series(const double scrn,const int Z1,const int Z2,const int a,const int b,const double r,const double arat)531 double MEAM::phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b,
532 const double r, const double arat)
533 {
534 double phi_sum = 0.0;
535 double b2nn, phi_val;
536 if (scrn > 0.0) {
537 b2nn = -Z2*scrn/Z1;
538 for (int n = 1; n <= 10; n++) {
539 phi_val = MathSpecial::powint(b2nn,n) * phi_meam(r * MathSpecial::powint(arat, n), a, b);
540 if (iszero(phi_val)) {
541 // once either term becomes zero at some point, all folliwng will also be zero
542 // necessary to avoid numerical error (nan or infty) due to exponential decay in phi_meam
543 break;
544 }
545 phi_sum += phi_val;
546 }
547 }
548 return phi_sum;
549 }
550
551 //----------------------------------------------------------------------c
552 // Compute background density for reference structure of each element
compute_reference_density()553 void MEAM::compute_reference_density()
554 {
555 int a, Z, Z2, errorflag;
556 double gam, Gbar, shp[3];
557 double rho0, rho0_2nn, arat, scrn;
558
559 // loop over element types
560 for (a = 0; a < this->neltypes; a++) {
561 Z = get_Zij(this->lattce_meam[a][a]);
562 if (this->ibar_meam[a] <= 0)
563 Gbar = 1.0;
564 else {
565 get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], shp);
566 gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z);
567 Gbar = G_gam(gam, this->ibar_meam[a], errorflag);
568 }
569
570 // The zeroth order density in the reference structure, with
571 // equilibrium spacing, is just the number of first neighbors times
572 // the rho0_meam coefficient...
573 rho0 = this->rho0_meam[a] * Z;
574
575 // ...unless we have unscreened second neighbors, in which case we
576 // add on the contribution from those (accounting for partial
577 // screening)
578 if (this->nn2_meam[a][a] == 1) {
579 Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
580 this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);
581 rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
582 rho0 = rho0 + Z2 * rho0_2nn * scrn;
583 }
584
585 this->rho_ref_meam[a] = rho0 * Gbar;
586 }
587 }
588
589 //------------------------------------------------------------------------------c
590 // Average weighting factors for the reference structure
get_tavref(double * t11av,double * t21av,double * t31av,double * t12av,double * t22av,double * t32av,double t11,double t21,double t31,double t12,double t22,double t32,double r,int a,int b,lattice_t latt)591 void MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, double* t22av, double* t32av,
592 double t11, double t21, double t31, double t12, double t22, double t32, double r, int a,
593 int b, lattice_t latt)
594 {
595 double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/;
596
597 // For ialloy = 2, no averaging is done
598 if (this->ialloy == 2) {
599 *t11av = t11;
600 *t21av = t21;
601 *t31av = t31;
602 *t12av = t12;
603 *t22av = t22;
604 *t32av = t32;
605 } else switch (latt) {
606 case FCC:
607 case BCC:
608 case DIA:
609 case DIA3:
610 case HCP:
611 case B1:
612 case DIM:
613 case B2:
614 case CH4:
615 case LIN:
616 case ZIG:
617 case TRI:
618 // all neighbors are of the opposite type
619 *t11av = t12;
620 *t21av = t22;
621 *t31av = t32;
622 *t12av = t11;
623 *t22av = t21;
624 *t32av = t31;
625 break;
626 default:
627 a1 = r / this->re_meam[a][a] - 1.0;
628 a2 = r / this->re_meam[b][b] - 1.0;
629 rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
630 rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
631 if (latt == L12) {
632 rho01 = 8 * rhoa01 + 4 * rhoa02;
633 *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
634 *t12av = t11;
635 *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01;
636 *t22av = t21;
637 *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01;
638 *t32av = t31;
639 } else {
640 // call error('Lattice not defined in get_tavref.')
641 }
642 }
643 }
644
645 //------------------------------------------------------------------------------c
get_sijk(double C,int i,int j,int k,double * sijk)646 void MEAM::get_sijk(double C, int i, int j, int k, double* sijk)
647 {
648 double x;
649 x = (C - this->Cmin_meam[i][j][k]) / (this->Cmax_meam[i][j][k] - this->Cmin_meam[i][j][k]);
650 *sijk = fcut(x);
651 }
652
653 //------------------------------------------------------------------------------c
654 // Calculate density functions, assuming reference configuration
get_densref(double r,int a,int b,double * rho01,double * rho11,double * rho21,double * rho31,double * rho02,double * rho12,double * rho22,double * rho32)655 void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, double* rho31,
656 double* rho02, double* rho12, double* rho22, double* rho32)
657 {
658 double a1, a2;
659 double s[3];
660 lattice_t lat;
661 int Zij,Zij2nn;
662 double rhoa01nn, rhoa02nn;
663 double rhoa01, rhoa11, rhoa21, rhoa31;
664 double rhoa02, rhoa12, rhoa22, rhoa32;
665 double arat, scrn, denom;
666 double C, s111, s112, s221, S11, S22;
667
668 a1 = r / this->re_meam[a][a] - 1.0;
669 a2 = r / this->re_meam[b][b] - 1.0;
670
671 rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
672 rhoa11 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1);
673 rhoa21 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1);
674 rhoa31 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1);
675 rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
676 rhoa12 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2);
677 rhoa22 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2);
678 rhoa32 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2);
679
680 lat = this->lattce_meam[a][b];
681
682 Zij = get_Zij(lat);
683
684 *rho11 = 0.0;
685 *rho21 = 0.0;
686 *rho31 = 0.0;
687 *rho12 = 0.0;
688 *rho22 = 0.0;
689 *rho32 = 0.0;
690
691 switch (lat) {
692 case FCC:
693 *rho01 = 12.0 * rhoa02;
694 *rho02 = 12.0 * rhoa01;
695 break;
696 case BCC:
697 *rho01 = 8.0 * rhoa02;
698 *rho02 = 8.0 * rhoa01;
699 break;
700 case B1:
701 *rho01 = 6.0 * rhoa02;
702 *rho02 = 6.0 * rhoa01;
703 break;
704 case DIA:
705 case DIA3:
706 *rho01 = 4.0 * rhoa02;
707 *rho02 = 4.0 * rhoa01;
708 *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
709 *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
710 break;
711 case HCP:
712 *rho01 = 12 * rhoa02;
713 *rho02 = 12 * rhoa01;
714 *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
715 *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
716 break;
717 case DIM:
718 get_shpfcn(DIM, 0, 0, s);
719 *rho01 = rhoa02;
720 *rho02 = rhoa01;
721 *rho11 = s[0] * rhoa12 * rhoa12;
722 *rho12 = s[0] * rhoa11 * rhoa11;
723 *rho21 = s[1] * rhoa22 * rhoa22;
724 *rho22 = s[1] * rhoa21 * rhoa21;
725 *rho31 = s[2] * rhoa32 * rhoa32;
726 *rho32 = s[2] * rhoa31 * rhoa31;
727 break;
728 case C11:
729 *rho01 = rhoa01;
730 *rho02 = rhoa02;
731 *rho11 = rhoa11;
732 *rho12 = rhoa12;
733 *rho21 = rhoa21;
734 *rho22 = rhoa22;
735 *rho31 = rhoa31;
736 *rho32 = rhoa32;
737 break;
738 case L12:
739 *rho01 = 8 * rhoa01 + 4 * rhoa02;
740 *rho02 = 12 * rhoa01;
741 if (this->ialloy == 1) {
742 *rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]);
743 denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]);
744 if (denom > 0.)
745 *rho21 = *rho21 / denom * *rho01;
746 } else
747 *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
748 break;
749 case B2:
750 *rho01 = 8.0 * rhoa02;
751 *rho02 = 8.0 * rhoa01;
752 break;
753 case CH4:
754 *rho01 = 4.0 * rhoa02; //in assumption that 'a' represent carbon
755 *rho02 = rhoa01; //in assumption that 'b' represent hydrogen
756
757 get_shpfcn(DIM, 0, 0, s); //H
758 *rho12 = s[0] * rhoa11 * rhoa11;
759 *rho22 = s[1] * rhoa21 * rhoa21;
760 *rho32 = s[2] * rhoa31 * rhoa31;
761
762 get_shpfcn(CH4, 0, 0, s); //C
763 *rho11 = s[0] * rhoa12 * rhoa12;
764 *rho21 = s[1] * rhoa22 * rhoa22;
765 *rho31 = s[2] * rhoa32 * rhoa32;
766 break;
767 case LIN:
768 *rho01 = rhoa02*Zij;
769 *rho02 = rhoa01*Zij;
770
771 get_shpfcn(LIN, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
772 *rho12 = s[0] * rhoa11 * rhoa11;
773 *rho22 = s[1] * rhoa21 * rhoa21;
774 *rho32 = s[2] * rhoa31 * rhoa31;
775 *rho11 = s[0] * rhoa12 * rhoa12;
776 *rho21 = s[1] * rhoa22 * rhoa22;
777 *rho31 = s[2] * rhoa32 * rhoa32;
778 break;
779 case ZIG:
780 *rho01 = rhoa02*Zij;
781 *rho02 = rhoa01*Zij;
782
783 get_shpfcn(ZIG, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
784 *rho12 = s[0] * rhoa11 * rhoa11;
785 *rho22 = s[1] * rhoa21 * rhoa21;
786 *rho32 = s[2] * rhoa31 * rhoa31;
787 *rho11 = s[0] * rhoa12 * rhoa12;
788 *rho21 = s[1] * rhoa22 * rhoa22;
789 *rho31 = s[2] * rhoa32 * rhoa32;
790 break;
791 case TRI:
792 *rho01 = rhoa02;
793 *rho02 = rhoa01*Zij;
794
795 get_shpfcn(TRI, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
796 *rho12 = s[0] * rhoa11 * rhoa11;
797 *rho22 = s[1] * rhoa21 * rhoa21;
798 *rho32 = s[2] * rhoa31 * rhoa31;
799 s[0] = 1.0;
800 s[1] = 2.0/3.0;
801 s[2] = 1.0 - 0.6*s[0];
802
803 *rho11 = s[0] * rhoa12 * rhoa12;
804 *rho21 = s[1] * rhoa22 * rhoa22;
805 *rho31 = s[2] * rhoa32 * rhoa32;
806 break;
807
808
809 // default:
810 // call error('Lattice not defined in get_densref.')
811 }
812
813 if (this->nn2_meam[a][b] == 1) {
814
815
816 Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b],
817 this->stheta_meam[a][b], arat, scrn);
818
819 a1 = arat * r / this->re_meam[a][a] - 1.0;
820 a2 = arat * r / this->re_meam[b][b] - 1.0;
821
822 rhoa01nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
823 rhoa02nn = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
824
825 if (lat == L12) {
826 // As usual, L12 thinks it's special; we need to be careful computing
827 // the screening functions
828 C = 1.0;
829 get_sijk(C, a, a, a, &s111);
830 get_sijk(C, a, a, b, &s112);
831 get_sijk(C, b, b, a, &s221);
832 S11 = s111 * s111 * s112 * s112;
833 S22 = s221 * s221 * s221 * s221;
834 *rho01 = *rho01 + 6 * S11 * rhoa01nn;
835 *rho02 = *rho02 + 6 * S22 * rhoa02nn;
836
837 } else {
838 // For other cases, assume that second neighbor is of same type,
839 // first neighbor may be of different type
840
841 *rho01 = *rho01 + Zij2nn * scrn * rhoa01nn;
842
843 // Assume Zij2nn and arat don't depend on order, but scrn might
844 Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a],
845 this->stheta_meam[a][b], arat, scrn);
846 *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;
847 }
848 }
849 }
850
851
interpolate_meam(int ind)852 void MEAM::interpolate_meam(int ind)
853 {
854 int j;
855 double drar;
856
857 // map to coefficient space
858
859 this->nrar = this->nr;
860 drar = this->dr;
861 this->rdrar = 1.0 / drar;
862
863 // phir interp
864 for (j = 0; j < this->nrar; j++) {
865 this->phirar[ind][j] = this->phir[ind][j];
866 }
867 this->phirar1[ind][0] = this->phirar[ind][1] - this->phirar[ind][0];
868 this->phirar1[ind][1] = 0.5 * (this->phirar[ind][2] - this->phirar[ind][0]);
869 this->phirar1[ind][this->nrar - 2] =
870 0.5 * (this->phirar[ind][this->nrar - 1] - this->phirar[ind][this->nrar - 3]);
871 this->phirar1[ind][this->nrar - 1] = 0.0;
872 for (j = 2; j < this->nrar - 2; j++) {
873 this->phirar1[ind][j] = ((this->phirar[ind][j - 2] - this->phirar[ind][j + 2]) +
874 8.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j - 1])) /
875 12.;
876 }
877
878 for (j = 0; j < this->nrar - 1; j++) {
879 this->phirar2[ind][j] = 3.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]) -
880 2.0 * this->phirar1[ind][j] - this->phirar1[ind][j + 1];
881 this->phirar3[ind][j] = this->phirar1[ind][j] + this->phirar1[ind][j + 1] -
882 2.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]);
883 }
884 this->phirar2[ind][this->nrar - 1] = 0.0;
885 this->phirar3[ind][this->nrar - 1] = 0.0;
886
887 for (j = 0; j < this->nrar; j++) {
888 this->phirar4[ind][j] = this->phirar1[ind][j] / drar;
889 this->phirar5[ind][j] = 2.0 * this->phirar2[ind][j] / drar;
890 this->phirar6[ind][j] = 3.0 * this->phirar3[ind][j] / drar;
891 }
892 }
893