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 /* ----------------------------------------------------------------------
16 Contributing author: Aidan Thompson (SNL) - original Tersoff implementation
17 David Farrell (NWU) - ZBL addition
18 ------------------------------------------------------------------------- */
19
20 #include "pair_tersoff_zbl_omp.h"
21
22 #include "comm.h"
23 #include "error.h"
24 #include "math_const.h"
25 #include "math_special.h"
26 #include "memory.h"
27 #include "potential_file_reader.h"
28 #include "tokenizer.h"
29 #include "update.h"
30
31 #include <cmath>
32 #include <cstring>
33
34 using namespace LAMMPS_NS;
35 using namespace MathConst;
36 using namespace MathSpecial;
37
38 #define DELTA 4
39
40 /* ----------------------------------------------------------------------
41 Fermi-like smoothing function
42 ------------------------------------------------------------------------- */
43
F_fermi(const double r,const double expsc,const double cut)44 static double F_fermi(const double r, const double expsc, const double cut)
45 {
46 return 1.0 / (1.0 + exp(-expsc*(r-cut)));
47 }
48
49 /* ----------------------------------------------------------------------
50 Fermi-like smoothing function derivative with respect to r
51 ------------------------------------------------------------------------- */
52
F_fermi_d(const double r,const double expsc,const double cut)53 static double F_fermi_d(const double r, const double expsc, const double cut)
54 {
55 return expsc*exp(-expsc*(r-cut)) / square(1.0 + exp(-expsc*(r-cut)));
56 }
57
58 /* ---------------------------------------------------------------------- */
59
PairTersoffZBLOMP(LAMMPS * lmp)60 PairTersoffZBLOMP::PairTersoffZBLOMP(LAMMPS *lmp) : PairTersoffOMP(lmp)
61 {
62 // hard-wired constants in metal or real units
63 // a0 = Bohr radius
64 // epsilon0 = permittivity of vacuum = q / energy-distance units
65 // e = unit charge
66 // 1 Kcal/mole = 0.043365121 eV
67
68 if (strcmp(update->unit_style,"metal") == 0) {
69 global_a_0 = 0.529;
70 global_epsilon_0 = 0.00552635;
71 global_e = 1.0;
72 } else if (strcmp(update->unit_style,"real") == 0) {
73 global_a_0 = 0.529;
74 global_epsilon_0 = 0.00552635 * 0.043365121;
75 global_e = 1.0;
76 } else error->all(FLERR,"Pair tersoff/zbl requires metal or real units");
77 }
78
79 /* ---------------------------------------------------------------------- */
80
read_file(char * file)81 void PairTersoffZBLOMP::read_file(char *file)
82 {
83 memory->sfree(params);
84 params = nullptr;
85 nparams = maxparam = 0;
86
87 // open file on proc 0
88
89 if (comm->me == 0) {
90 PotentialFileReader reader(PairTersoff::lmp, file, "tersoff/zbl",
91 unit_convert_flag);
92 char * line;
93
94 // transparently convert units for supported conversions
95
96 int unit_convert = reader.get_unit_convert();
97 double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
98 unit_convert);
99
100 while ((line = reader.next_line(NPARAMS_PER_LINE))) {
101 try {
102 ValueTokenizer values(line);
103
104 std::string iname = values.next_string();
105 std::string jname = values.next_string();
106 std::string kname = values.next_string();
107
108 // ielement,jelement,kelement = 1st args
109 // if all 3 args are in element list, then parse this line
110 // else skip to next entry in file
111 int ielement, jelement, kelement;
112
113 for (ielement = 0; ielement < nelements; ielement++)
114 if (iname == elements[ielement]) break;
115 if (ielement == nelements) continue;
116 for (jelement = 0; jelement < nelements; jelement++)
117 if (jname == elements[jelement]) break;
118 if (jelement == nelements) continue;
119 for (kelement = 0; kelement < nelements; kelement++)
120 if (kname == elements[kelement]) break;
121 if (kelement == nelements) continue;
122
123 // load up parameter settings and error check their values
124
125 if (nparams == maxparam) {
126 maxparam += DELTA;
127 params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
128 "pair:params");
129 // make certain all addional allocated storage is initialized
130 // to avoid false positives when checking with valgrind
131
132 memset(params + nparams, 0, DELTA*sizeof(Param));
133 }
134
135 params[nparams].ielement = ielement;
136 params[nparams].jelement = jelement;
137 params[nparams].kelement = kelement;
138 params[nparams].powerm = values.next_double();
139 params[nparams].gamma = values.next_double();
140 params[nparams].lam3 = values.next_double();
141 params[nparams].c = values.next_double();
142 params[nparams].d = values.next_double();
143 params[nparams].h = values.next_double();
144 params[nparams].powern = values.next_double();
145 params[nparams].beta = values.next_double();
146 params[nparams].lam2 = values.next_double();
147 params[nparams].bigb = values.next_double();
148 params[nparams].bigr = values.next_double();
149 params[nparams].bigd = values.next_double();
150 params[nparams].lam1 = values.next_double();
151 params[nparams].biga = values.next_double();
152 params[nparams].Z_i = values.next_double();
153 params[nparams].Z_j = values.next_double();
154 params[nparams].ZBLcut = values.next_double();
155 params[nparams].ZBLexpscale = values.next_double();
156 params[nparams].powermint = int(params[nparams].powerm);
157
158 if (unit_convert) {
159 params[nparams].biga *= conversion_factor;
160 params[nparams].bigb *= conversion_factor;
161 }
162 } catch (TokenizerException &e) {
163 error->one(FLERR, e.what());
164 }
165
166 // currently only allow m exponent of 1 or 3
167 if (params[nparams].c < 0.0 ||
168 params[nparams].d < 0.0 ||
169 params[nparams].powern < 0.0 ||
170 params[nparams].beta < 0.0 ||
171 params[nparams].lam2 < 0.0 ||
172 params[nparams].bigb < 0.0 ||
173 params[nparams].bigr < 0.0 ||
174 params[nparams].bigd < 0.0 ||
175 params[nparams].bigd > params[nparams].bigr ||
176 params[nparams].lam1 < 0.0 ||
177 params[nparams].biga < 0.0 ||
178 params[nparams].powerm - params[nparams].powermint != 0.0 ||
179 (params[nparams].powermint != 3 &&
180 params[nparams].powermint != 1) ||
181 params[nparams].gamma < 0.0 ||
182 params[nparams].Z_i < 1.0 ||
183 params[nparams].Z_j < 1.0 ||
184 params[nparams].ZBLcut < 0.0 ||
185 params[nparams].ZBLexpscale < 0.0)
186 error->one(FLERR,"Illegal Tersoff parameter");
187
188 nparams++;
189 }
190 }
191
192 MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
193 MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
194
195 if (comm->me != 0) {
196 params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
197 }
198
199 MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
200 }
201
202 /* ---------------------------------------------------------------------- */
203
force_zeta(Param * param,double rsq,double zeta_ij,double & fforce,double & prefactor,int eflag,double & eng)204 void PairTersoffZBLOMP::force_zeta(Param *param, double rsq, double zeta_ij,
205 double &fforce, double &prefactor,
206 int eflag, double &eng)
207 {
208 double r,fa,fa_d,bij;
209
210 r = sqrt(rsq);
211
212 fa = (r > param->bigr + param->bigd) ? 0.0 :
213 -param->bigb * exp(-param->lam2 * r) * ters_fc(r,param) *
214 F_fermi(r,param->ZBLexpscale,param->ZBLcut);
215
216 fa_d = (r > param->bigr + param->bigd) ? 0.0 :
217 param->bigb * exp(-param->lam2 * r) *
218 (param->lam2 * ters_fc(r,param) *
219 F_fermi(r,param->ZBLexpscale,param->ZBLcut) -
220 ters_fc_d(r,param) * F_fermi(r,param->ZBLexpscale,param->ZBLcut)
221 - ters_fc(r,param) * F_fermi_d(r,param->ZBLexpscale,param->ZBLcut));
222
223 bij = ters_bij(zeta_ij,param);
224 fforce = 0.5*bij*fa_d;
225 prefactor = -0.5*fa * ters_bij_d(zeta_ij,param);
226 if (eflag) eng = 0.5*bij*fa;
227 }
228
229 /* ---------------------------------------------------------------------- */
230
repulsive(Param * param,double rsq,double & fforce,int eflag,double & eng)231 void PairTersoffZBLOMP::repulsive(Param *param, double rsq, double &fforce,
232 int eflag, double &eng)
233 {
234 double r,tmp_fc,tmp_fc_d,tmp_exp;
235
236 // Tersoff repulsive portion
237
238 r = sqrt(rsq);
239 tmp_fc = ters_fc(r,param);
240 tmp_fc_d = ters_fc_d(r,param);
241 tmp_exp = exp(-param->lam1 * r);
242 double fforce_ters = param->biga * tmp_exp * (tmp_fc_d - tmp_fc*param->lam1);
243 double eng_ters = tmp_fc * param->biga * tmp_exp;
244
245 // ZBL repulsive portion
246
247 double esq = square(global_e);
248 double a_ij = (0.8854*global_a_0) /
249 (pow(param->Z_i,0.23) + pow(param->Z_j,0.23));
250 double premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_PI*global_epsilon_0);
251 double r_ov_a = r/a_ij;
252 double phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) +
253 0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a);
254 double dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) -
255 0.9423*0.5099*exp(-0.9423*r_ov_a) -
256 0.4029*0.2802*exp(-0.4029*r_ov_a) -
257 0.2016*0.02817*exp(-0.2016*r_ov_a));
258 double fforce_ZBL = premult*-phi/rsq + premult*dphi/r;
259 double eng_ZBL = premult*(1.0/r)*phi;
260
261 // combine two parts with smoothing by Fermi-like function
262
263 fforce = -(-F_fermi_d(r,param->ZBLexpscale,param->ZBLcut) * eng_ZBL +
264 (1.0 - F_fermi(r,param->ZBLexpscale,param->ZBLcut))*fforce_ZBL +
265 F_fermi_d(r,param->ZBLexpscale,param->ZBLcut)*eng_ters +
266 F_fermi(r,param->ZBLexpscale,param->ZBLcut)*fforce_ters) / r;
267
268 if (eflag)
269 eng = (1.0 - F_fermi(r,param->ZBLexpscale,param->ZBLcut))*eng_ZBL +
270 F_fermi(r,param->ZBLexpscale,param->ZBLcut)*eng_ters;
271 }
272