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
15 #include "meam.h"
16
17 #include "math_const.h"
18
19 #include <algorithm>
20 #include <cmath>
21
22 using namespace LAMMPS_NS;
23 using namespace MathConst;
24
25 // do a sanity check on index parameters
26 void
meam_checkindex(int num,int lim,int nidx,int * idx,int * ierr)27 MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr)
28 {
29 //: idx[0..2]
30 *ierr = 0;
31 if (nidx < num) {
32 *ierr = 2;
33 return;
34 }
35
36 for (int i = 0; i < num; i++) {
37 if ((idx[i] < 0) || (idx[i] >= lim)) {
38 *ierr = 3;
39 return;
40 }
41 }
42 }
43
44 // The "which" argument corresponds to the index of the "keyword" array
45 // in pair_meam.cpp:
46 //
47 // 0 = Ec_meam
48 // 1 = alpha_meam
49 // 2 = rho0_meam
50 // 3 = delta_meam
51 // 4 = lattce_meam
52 // 5 = attrac_meam
53 // 6 = repuls_meam
54 // 7 = nn2_meam
55 // 8 = Cmin_meam
56 // 9 = Cmax_meam
57 // 10 = rc_meam
58 // 11 = delr_meam
59 // 12 = augt1
60 // 13 = gsmooth_factor
61 // 14 = re_meam
62 // 15 = ialloy
63 // 16 = mixture_ref_t
64 // 17 = erose_form
65 // 18 = zbl_meam
66 // 19 = emb_lin_neg
67 // 20 = bkgd_dyn
68 // 21 = theta
69
70 // The returned errorflag has the following meanings:
71
72 // 0 = no error
73 // 1 = "which" out of range / invalid keyword
74 // 2 = not enough indices given
75 // 3 = an element index is out of range
76
77 void
meam_setup_param(int which,double value,int nindex,int * index,int * errorflag)78 MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag)
79 {
80 //: index[0..2]
81 int i1, i2;
82 lattice_t vlat;
83 *errorflag = 0;
84
85 switch (which) {
86 // 0 = Ec_meam
87 case 0:
88 meam_checkindex(2, neltypes, nindex, index, errorflag);
89 if (*errorflag != 0)
90 return;
91 this->Ec_meam[index[0]][index[1]] = value;
92 break;
93
94 // 1 = alpha_meam
95 case 1:
96 meam_checkindex(2, neltypes, nindex, index, errorflag);
97 if (*errorflag != 0)
98 return;
99 this->alpha_meam[index[0]][index[1]] = value;
100 break;
101
102 // 2 = rho0_meam
103 case 2:
104 meam_checkindex(1, neltypes, nindex, index, errorflag);
105 if (*errorflag != 0)
106 return;
107 this->rho0_meam[index[0]] = value;
108 break;
109
110 // 3 = delta_meam
111 case 3:
112 meam_checkindex(2, neltypes, nindex, index, errorflag);
113 if (*errorflag != 0)
114 return;
115 this->delta_meam[index[0]][index[1]] = value;
116 break;
117
118 // 4 = lattce_meam
119 case 4:
120 meam_checkindex(2, neltypes, nindex, index, errorflag);
121 if (*errorflag != 0)
122 return;
123 vlat = (lattice_t)value;
124
125 this->lattce_meam[index[0]][index[1]] = vlat;
126 break;
127
128 // 5 = attrac_meam
129 case 5:
130 meam_checkindex(2, neltypes, nindex, index, errorflag);
131 if (*errorflag != 0)
132 return;
133 this->attrac_meam[index[0]][index[1]] = value;
134 break;
135
136 // 6 = repuls_meam
137 case 6:
138 meam_checkindex(2, neltypes, nindex, index, errorflag);
139 if (*errorflag != 0)
140 return;
141 this->repuls_meam[index[0]][index[1]] = value;
142 break;
143
144 // 7 = nn2_meam
145 case 7:
146 meam_checkindex(2, neltypes, nindex, index, errorflag);
147 if (*errorflag != 0)
148 return;
149 i1 = std::min(index[0], index[1]);
150 i2 = std::max(index[0], index[1]);
151 this->nn2_meam[i1][i2] = (int)value;
152 break;
153
154 // 8 = Cmin_meam
155 case 8:
156 meam_checkindex(3, neltypes, nindex, index, errorflag);
157 if (*errorflag != 0)
158 return;
159 this->Cmin_meam[index[0]][index[1]][index[2]] = value;
160 break;
161
162 // 9 = Cmax_meam
163 case 9:
164 meam_checkindex(3, neltypes, nindex, index, errorflag);
165 if (*errorflag != 0)
166 return;
167 this->Cmax_meam[index[0]][index[1]][index[2]] = value;
168 break;
169
170 // 10 = rc_meam
171 case 10:
172 this->rc_meam = value;
173 break;
174
175 // 11 = delr_meam
176 case 11:
177 this->delr_meam = value;
178 break;
179
180 // 12 = augt1
181 case 12:
182 this->augt1 = (int)value;
183 break;
184
185 // 13 = gsmooth
186 case 13:
187 this->gsmooth_factor = value;
188 break;
189
190 // 14 = re_meam
191 case 14:
192 meam_checkindex(2, neltypes, nindex, index, errorflag);
193 if (*errorflag != 0)
194 return;
195 this->re_meam[index[0]][index[1]] = value;
196 break;
197
198 // 15 = ialloy
199 case 15:
200 this->ialloy = (int)value;
201 break;
202
203 // 16 = mixture_ref_t
204 case 16:
205 this->mix_ref_t = (int)value;
206 break;
207
208 // 17 = erose_form
209 case 17:
210 this->erose_form = (int)value;
211 break;
212
213 // 18 = zbl_meam
214 case 18:
215 meam_checkindex(2, neltypes, nindex, index, errorflag);
216 if (*errorflag != 0)
217 return;
218 i1 = std::min(index[0], index[1]);
219 i2 = std::max(index[0], index[1]);
220 this->zbl_meam[i1][i2] = (int)value;
221 break;
222
223 // 19 = emb_lin_neg
224 case 19:
225 this->emb_lin_neg = (int)value;
226 break;
227
228 // 20 = bkgd_dyn
229 case 20:
230 this->bkgd_dyn = (int)value;
231 break;
232
233 // 21 = theta
234 // see alloyparams(void) in meam_setup_done.cpp
235 case 21:
236 // const double PI = 3.141592653589793238463;
237 meam_checkindex(2, neltypes, nindex, index, errorflag);
238 if (*errorflag != 0)
239 return;
240 i1 = std::min(index[0], index[1]);
241 i2 = std::max(index[0], index[1]);
242 // we don't use theta, instead stheta and ctheta
243 this->stheta_meam[i1][i2] = sin(value/2*MY_PI/180.0);
244 this->ctheta_meam[i1][i2] = cos(value/2*MY_PI/180.0);
245 break;
246
247 default:
248 *errorflag = 1;
249 }
250 }
251