1 /* -*- c++ -*- ----------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 #ifndef LMP_MEAM_H
15 #define LMP_MEAM_H
16
17 #include "math_const.h" // IWYU pragma: export
18
19 #include <cmath>
20 #include <string>
21
22 #define maxelt 5
23
24 namespace LAMMPS_NS {
25 class Memory;
26
27 typedef enum { FCC, BCC, HCP, DIM, DIA, DIA3, B1, C11, L12, B2, CH4, LIN, ZIG, TRI } lattice_t;
28
29 class MEAM {
30 public:
31 MEAM(Memory *mem);
32 ~MEAM();
33
34 private:
35 Memory *memory;
36
37 // cutforce = force cutoff
38 // cutforcesq = force cutoff squared
39
40 double cutforce, cutforcesq;
41
42 // Ec_meam = cohesive energy
43 // re_meam = nearest-neighbor distance
44 // B_meam = bulk modulus
45 // ielt_meam = atomic number of element
46 // A_meam = adjustable parameter
47 // alpha_meam = sqrt(9*Omega*B/Ec)
48 // rho0_meam = density scaling parameter
49 // delta_meam = heat of formation for alloys
50 // beta[0-3]_meam = electron density constants
51 // t[0-3]_meam = coefficients on densities in Gamma computation
52 // rho_ref_meam = background density for reference structure
53 // ibar_meam(i) = selection parameter for Gamma function for elt i,
54 // lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j)
55 // neltypes = maximum number of element type defined
56 // eltind = index number of pair (similar to Voigt notation; ij = ji)
57 // phir = pair potential function array
58 // phirar[1-6] = spline coeffs
59 // attrac_meam = attraction parameter in Rose energy
60 // repuls_meam = repulsion parameter in Rose energy
61 // nn2_meam = 1 if second nearest neighbors are to be computed, else 0
62 // zbl_meam = 1 if zbl potential for small r to be use, else 0
63 // emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0
64 // bkgd_dyn = 1 if reference densities follows Dynamo, else 0
65 // Cmin_meam, Cmax_meam = min and max values in screening cutoff
66 // rc_meam = cutoff distance for meam
67 // delr_meam = cutoff region for meam
68 // ebound_meam = factor giving maximum boundary of sceen fcn ellipse
69 // augt1 = flag for whether t1 coefficient should be augmented
70 // ialloy = flag for newer alloy formulation (as in dynamo code)
71 // mix_ref_t = flag to recover "old" way of computing t in reference config
72 // erose_form = selection parameter for form of E_rose function
73 // gsmooth_factor = factor determining length of G smoothing region
74 // vind[23]D = Voight notation index maps for 2 and 3D
75 // v2D,v3D = array of factors to apply for Voight notation
76
77 // nr,dr = pair function discretization parameters
78 // nrar,rdrar = spline coeff array parameters
79
80 // theta = angle between three atoms in line, zigzag, and trimer reference structures
81 // stheta_meam = sin(theta/2) in radian used in line, zigzag, and trimer reference structures
82 // ctheta_meam = cos(theta/2) in radian used in line, zigzag, and trimer reference structures
83
84 double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
85 double A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt];
86 double delta_meam[maxelt][maxelt];
87 double beta0_meam[maxelt], beta1_meam[maxelt];
88 double beta2_meam[maxelt], beta3_meam[maxelt];
89 double t0_meam[maxelt], t1_meam[maxelt];
90 double t2_meam[maxelt], t3_meam[maxelt];
91 double rho_ref_meam[maxelt];
92 int ibar_meam[maxelt], ielt_meam[maxelt];
93 lattice_t lattce_meam[maxelt][maxelt];
94 int nn2_meam[maxelt][maxelt];
95 int zbl_meam[maxelt][maxelt];
96 int eltind[maxelt][maxelt];
97 int neltypes;
98
99 double **phir;
100
101 double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5, **phirar6;
102
103 double attrac_meam[maxelt][maxelt], repuls_meam[maxelt][maxelt];
104
105 double Cmin_meam[maxelt][maxelt][maxelt];
106 double Cmax_meam[maxelt][maxelt][maxelt];
107 double rc_meam, delr_meam, ebound_meam[maxelt][maxelt];
108 int augt1, ialloy, mix_ref_t, erose_form;
109 int emb_lin_neg, bkgd_dyn;
110 double gsmooth_factor;
111
112 int vind2D[3][3], vind3D[3][3][3]; // x-y-z to Voigt-like index
113 int v2D[6], v3D[10]; // multiplicity of Voigt index (i.e. [1] -> xy+yx = 2
114
115 int nr, nrar;
116 double dr, rdrar;
117
118 public:
119 int nmax;
120 double *rho, *rho0, *rho1, *rho2, *rho3, *frhop;
121 double *gamma, *dgamma1, *dgamma2, *dgamma3, *arho2b;
122 double **arho1, **arho2, **arho3, **arho3b, **t_ave, **tsq_ave;
123
124 int maxneigh;
125 double *scrfcn, *dscrfcn, *fcpair;
126
127 //angle for trimer, zigzag, line reference structures
128 double stheta_meam[maxelt][maxelt];
129 double ctheta_meam[maxelt][maxelt];
130
131 protected:
132 // meam_funcs.cpp
133
134 //-----------------------------------------------------------------------------
135 // Cutoff function
136 //
fcut(const double xi)137 static double fcut(const double xi)
138 {
139 double a;
140 if (xi >= 1.0)
141 return 1.0;
142 else if (xi <= 0.0)
143 return 0.0;
144 else {
145 // ( 1.d0 - (1.d0 - xi)**4 )**2, but with better codegen
146 a = 1.0 - xi;
147 a *= a;
148 a *= a;
149 a = 1.0 - a;
150 return a * a;
151 }
152 }
153
154 //-----------------------------------------------------------------------------
155 // Cutoff function and derivative
156 //
dfcut(const double xi,double & dfc)157 static double dfcut(const double xi, double &dfc)
158 {
159 double a, a3, a4, a1m4;
160 if (xi >= 1.0) {
161 dfc = 0.0;
162 return 1.0;
163 } else if (xi <= 0.0) {
164 dfc = 0.0;
165 return 0.0;
166 } else {
167 a = 1.0 - xi;
168 a3 = a * a * a;
169 a4 = a * a3;
170 a1m4 = 1.0 - a4;
171
172 dfc = 8 * a1m4 * a3;
173 return a1m4 * a1m4;
174 }
175 }
176
177 //-----------------------------------------------------------------------------
178 // Derivative of Cikj w.r.t. rij
179 // Inputs: rij,rij2,rik2,rjk2
180 //
dCfunc(const double rij2,const double rik2,const double rjk2)181 static double dCfunc(const double rij2, const double rik2, const double rjk2)
182 {
183 double rij4, a, asq, b, denom;
184
185 rij4 = rij2 * rij2;
186 a = rik2 - rjk2;
187 b = rik2 + rjk2;
188 asq = a * a;
189 denom = rij4 - asq;
190 denom = denom * denom;
191 return -4 * (-2 * rij2 * asq + rij4 * b + asq * b) / denom;
192 }
193
194 //-----------------------------------------------------------------------------
195 // Derivative of Cikj w.r.t. rik and rjk
196 // Inputs: rij,rij2,rik2,rjk2
197 //
dCfunc2(const double rij2,const double rik2,const double rjk2,double & dCikj1,double & dCikj2)198 static void dCfunc2(const double rij2, const double rik2, const double rjk2, double &dCikj1,
199 double &dCikj2)
200 {
201 double rij4, rik4, rjk4, a, denom;
202
203 rij4 = rij2 * rij2;
204 rik4 = rik2 * rik2;
205 rjk4 = rjk2 * rjk2;
206 a = rik2 - rjk2;
207 denom = rij4 - a * a;
208 denom = denom * denom;
209 dCikj1 = 4 * rij2 * (rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom;
210 dCikj2 = 4 * rij2 * (rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom;
211 }
212
213 double G_gam(const double gamma, const int ibar, int &errorflag) const;
214 double dG_gam(const double gamma, const int ibar, double &dG) const;
215 static double zbl(const double r, const int z1, const int z2);
216 double embedding(const double A, const double Ec, const double rhobar, double &dF) const;
217 static double erose(const double r, const double re, const double alpha, const double Ec,
218 const double repuls, const double attrac, const int form);
219
220 static void get_shpfcn(const lattice_t latt, const double sthe, const double cthe,
221 double (&s)[3]);
222
223 static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, const double sthe,
224 double &a, double &S);
225 static int get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S);
226
227 protected:
228 void meam_checkindex(int, int, int, int *, int *);
229 void getscreen(int i, double *scrfcn, double *dscrfcn, double *fcpair, double **x, int numneigh,
230 int *firstneigh, int numneigh_full, int *firstneigh_full, int ntype, int *type,
231 int *fmap);
232 void calc_rho1(int i, int ntype, int *type, int *fmap, double **x, int numneigh, int *firstneigh,
233 double *scrfcn, double *fcpair);
234
235 void alloyparams();
236 void compute_pair_meam();
237 double phi_meam(double, int, int);
238 double phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b,
239 const double r, const double arat);
240 void compute_reference_density();
241 void get_tavref(double *, double *, double *, double *, double *, double *, double, double,
242 double, double, double, double, double, int, int, lattice_t);
243 void get_sijk(double, int, int, int, double *);
244 void get_densref(double, int, int, double *, double *, double *, double *, double *, double *,
245 double *, double *);
246 void interpolate_meam(int);
247
248 public:
249 // clang-format off
250 //-----------------------------------------------------------------------------
251 // convert lattice spec to lattice_t
252 // only use single-element lattices if single=true
253 // return false on failure
254 // return true and set lat on success
str_to_lat(const std::string & str,bool single,lattice_t & lat)255 static bool str_to_lat(const std::string & str, bool single, lattice_t& lat)
256 {
257 if (str == "fcc") lat = FCC;
258 else if (str == "bcc") lat = BCC;
259 else if (str == "hcp") lat = HCP;
260 else if (str == "dim") lat = DIM;
261 else if (str == "dia") lat = DIA;
262 else if (str == "dia3") lat = DIA3;
263 else if (str == "lin") lat = LIN;
264 else if (str == "zig") lat = ZIG;
265 else if (str == "tri") lat = TRI;
266 else {
267 if (single)
268 return false;
269
270 if (str == "b1") lat = B1;
271 else if (str == "c11") lat = C11;
272 else if (str == "l12") lat = L12;
273 else if (str == "b2") lat = B2;
274 else if (str == "ch4") lat = CH4;
275 else if (str == "lin") lat =LIN;
276 else if (str == "zig") lat = ZIG;
277 else if (str == "tri") lat = TRI;
278 else return false;
279 }
280 return true;
281 }
282 // clang-format on
283 static int get_Zij(const lattice_t latt);
284 void meam_setup_global(int nelt, lattice_t *lat, int *ielement, double *atwt, double *alpha,
285 double *b0, double *b1, double *b2, double *b3, double *alat, double *esub,
286 double *asub, double *t0, double *t1, double *t2, double *t3,
287 double *rozero, int *ibar);
288 void meam_setup_param(int which, double value, int nindex, int *index /*index(3)*/,
289 int *errorflag);
290 void meam_setup_done(double *cutmax);
291 void meam_dens_setup(int atom_nmax, int nall, int n_neigh);
292 void meam_dens_init(int i, int ntype, int *type, int *fmap, double **x, int numneigh,
293 int *firstneigh, int numneigh_full, int *firstneigh_full, int fnoffset);
294 void meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom,
295 double *eng_vdwl, double *eatom, int ntype, int *type, int *fmap,
296 double **scale, int &errorflag);
297 void meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int vflag_atom,
298 double *eng_vdwl, double *eatom, int ntype, int *type, int *fmap, double **scale,
299 double **x, int numneigh, int *firstneigh, int numneigh_full,
300 int *firstneigh_full, int fnoffset, double **f, double **vatom, double *virial);
301 };
302
303 // Functions we need for compat
304
iszero(const double f)305 static inline bool iszero(const double f)
306 {
307 return fabs(f) < 1e-20;
308 }
309
isone(const double f)310 static inline bool isone(const double f)
311 {
312 return fabs(f - 1.0) < 1e-20;
313 }
314
315 // Helper functions
316
fdiv_zero(const double n,const double d)317 static inline double fdiv_zero(const double n, const double d)
318 {
319 if (iszero(d)) return 0.0;
320 return n / d;
321 }
322
323 } // namespace LAMMPS_NS
324 #endif
325