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 using namespace LAMMPS_NS;
17
18 void
meam_dens_final(int nlocal,int eflag_either,int eflag_global,int eflag_atom,double * eng_vdwl,double * eatom,int,int * type,int * fmap,double ** scale,int & errorflag)19 MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,
20 double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, int& errorflag)
21 {
22 int i, elti;
23 int m;
24 double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z;
25 double denom, rho_bkgd, Fl;
26 double scaleii;
27
28 // Complete the calculation of density
29
30 for (i = 0; i < nlocal; i++) {
31 elti = fmap[type[i]];
32 if (elti >= 0) {
33 scaleii = scale[type[i]][type[i]];
34 rho1[i] = 0.0;
35 rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i];
36 rho3[i] = 0.0;
37 for (m = 0; m < 3; m++) {
38 rho1[i] = rho1[i] + arho1[i][m] * arho1[i][m];
39 rho3[i] = rho3[i] - 3.0 / 5.0 * arho3b[i][m] * arho3b[i][m];
40 }
41 for (m = 0; m < 6; m++) {
42 rho2[i] = rho2[i] + this->v2D[m] * arho2[i][m] * arho2[i][m];
43 }
44 for (m = 0; m < 10; m++) {
45 rho3[i] = rho3[i] + this->v3D[m] * arho3[i][m] * arho3[i][m];
46 }
47
48 if (rho0[i] > 0.0) {
49 if (this->ialloy == 1) {
50 t_ave[i][0] = fdiv_zero(t_ave[i][0], tsq_ave[i][0]);
51 t_ave[i][1] = fdiv_zero(t_ave[i][1], tsq_ave[i][1]);
52 t_ave[i][2] = fdiv_zero(t_ave[i][2], tsq_ave[i][2]);
53 } else if (this->ialloy == 2) {
54 t_ave[i][0] = this->t1_meam[elti];
55 t_ave[i][1] = this->t2_meam[elti];
56 t_ave[i][2] = this->t3_meam[elti];
57 } else {
58 t_ave[i][0] = t_ave[i][0] / rho0[i];
59 t_ave[i][1] = t_ave[i][1] / rho0[i];
60 t_ave[i][2] = t_ave[i][2] / rho0[i];
61 }
62 }
63
64 gamma[i] = t_ave[i][0] * rho1[i] + t_ave[i][1] * rho2[i] + t_ave[i][2] * rho3[i];
65
66 if (rho0[i] > 0.0) {
67 gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
68 }
69
70 Z = get_Zij(this->lattce_meam[elti][elti]);
71
72 G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
73 if (errorflag != 0)
74 return;
75
76 get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shp);
77
78 if (this->ibar_meam[elti] <= 0) {
79 Gbar = 1.0;
80 dGbar = 0.0;
81 } else {
82 if (this->mix_ref_t == 1) {
83 gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
84 } else {
85 gam = (this->t1_meam[elti] * shp[0] + this->t2_meam[elti] * shp[1] + this->t3_meam[elti] * shp[2]) /
86 (Z * Z);
87 }
88 Gbar = G_gam(gam, this->ibar_meam[elti], errorflag);
89 }
90 rho[i] = rho0[i] * G;
91
92 if (this->mix_ref_t == 1) {
93 if (this->ibar_meam[elti] <= 0) {
94 Gbar = 1.0;
95 dGbar = 0.0;
96 } else {
97 gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
98 Gbar = dG_gam(gam, this->ibar_meam[elti], dGbar);
99 }
100 rho_bkgd = this->rho0_meam[elti] * Z * Gbar;
101 } else {
102 if (this->bkgd_dyn == 1) {
103 rho_bkgd = this->rho0_meam[elti] * Z;
104 } else {
105 rho_bkgd = this->rho_ref_meam[elti];
106 }
107 }
108 rhob = rho[i] / rho_bkgd;
109 denom = 1.0 / rho_bkgd;
110
111 G = dG_gam(gamma[i], this->ibar_meam[elti], dG);
112
113 dgamma1[i] = (G - 2 * dG * gamma[i]) * denom;
114
115 if (!iszero(rho0[i])) {
116 dgamma2[i] = (dG / rho0[i]) * denom;
117 } else {
118 dgamma2[i] = 0.0;
119 }
120
121 // dgamma3 is nonzero only if we are using the "mixed" rule for
122 // computing t in the reference system (which is not correct, but
123 // included for backward compatibility
124 if (this->mix_ref_t == 1) {
125 dgamma3[i] = rho0[i] * G * dGbar / (Gbar * Z * Z) * denom;
126 } else {
127 dgamma3[i] = 0.0;
128 }
129
130 Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);
131
132 if (eflag_either != 0) {
133 Fl *= scaleii;
134 if (eflag_global != 0) {
135 *eng_vdwl = *eng_vdwl + Fl;
136 }
137 if (eflag_atom != 0) {
138 eatom[i] = eatom[i] + Fl;
139 }
140 }
141 }
142 }
143 }
144
145