1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 http://lammps.sandia.gov, 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 Contributing author: Aidan Thompson (SNL) - original Tersoff implementation
15 David Farrell (NWU) - ZBL addition
16
17 Modified for use with KIM by Tobias Brink (2020).
18 ------------------------------------------------------------------------- */
19
20 #include "pair_tersoff_zbl.hpp"
21
22 #include <cmath>
23 #include <fstream>
24 #include <sstream>
25 #include <iomanip>
26
27 using namespace model_driver_Tersoff;
28 using namespace std;
29
30 /* ---------------------------------------------------------------------- */
31
PairTersoffZBL(const string & parameter_file,int n_spec,map<string,int> type_map,double energy_conv,double inv_energy_conv,double length_conv,double inv_length_conv,double charge_conv)32 PairTersoffZBL::PairTersoffZBL(const string& parameter_file,
33 int n_spec,
34 map<string,int> type_map,
35 // Conversion factors.
36 double energy_conv,
37 double inv_energy_conv,
38 double length_conv,
39 double inv_length_conv,
40 double charge_conv)
41 : PairTersoff(n_spec, type_map),
42 kim_params_zbl(n_spec),
43 params_zbl_2(n_spec, n_spec),
44 global_a_0(0.529 * length_conv),
45 global_epsilon_0(0.00552635
46 * charge_conv * charge_conv
47 * inv_energy_conv * inv_length_conv),
48 global_e(1.0 * charge_conv),
49 global_e_sq(global_e * global_e)
50 {
51 // Read parameter file.
52 std::fstream f(parameter_file.c_str(), std::ios_base::in);
53 read_params(f, type_map, energy_conv, length_conv, inv_length_conv);
54 }
55
~PairTersoffZBL()56 PairTersoffZBL::~PairTersoffZBL() {}
57
58 /* ---------------------------------------------------------------------- */
59
60 /* Read parameters.
61
62 Reads the parameters from an input file and stores
63 them. Automatically calls prepare_params() to check validity of
64 parameters and pre-compute some values.
65
66 This is a bit of a copy-paste from pair_tersoff.cpp, but this way
67 is still easiest.
68 */
read_params(istream & infile,std::map<string,int> type_map,double energy_conv,double length_conv,double inv_length_conv)69 void PairTersoffZBL::read_params(istream& infile, std::map<string,int> type_map,
70 double energy_conv,
71 double length_conv,
72 double inv_length_conv) {
73 // Strip comment lines.
74 stringstream buffer;
75 string line;
76 while(getline(infile, line))
77 buffer << line.substr(0, line.find('#'));
78 // Read in parameters.
79 Array3D<bool> got_interaction(n_spec,n_spec,n_spec);
80 got_interaction = false;
81 Params2 temp_params2;
82 Params3 temp_params3;
83 ParamsZBL2 temp_params_zbl_2;
84 double m; // m is an integer but some input files use floating point
85 // notation, so we need to read into a double variable,
86 // otherwise the C++ standard library chokes on the input.
87 string type_i, type_j, type_k;
88 double c, d; // Those are only stored in the published KIM parameters
89 // since they are not needed for computation, we use
90 // precomputed c², d² and c²/d² instead.
91 double Z_i, Z_j; // Equivalent for these.
92 while (buffer >> type_i
93 >> type_j
94 >> type_k
95 >> m
96 >> temp_params3.gamma
97 >> temp_params3.lam3
98 >> c
99 >> d
100 >> temp_params3.h // costheta0
101 >> temp_params2.n
102 >> temp_params2.beta
103 >> temp_params2.lam2
104 >> temp_params2.B
105 >> temp_params3.R
106 >> temp_params3.D
107 >> temp_params2.lam1
108 >> temp_params2.A
109 >> Z_i
110 >> Z_j
111 >> temp_params_zbl_2.ZBLcut
112 >> temp_params_zbl_2.ZBLexpscale) {
113 // Convert m to integer.
114 temp_params3.m = m;
115 if (abs(m - temp_params3.m) > 1e-8)
116 throw runtime_error("m must be an integer");
117 // Unit conversion.
118 temp_params2.A *= energy_conv;
119 temp_params2.B *= energy_conv;
120 temp_params2.lam1 *= inv_length_conv;
121 temp_params2.lam2 *= inv_length_conv;
122 temp_params3.lam3 *= inv_length_conv;
123 temp_params3.R *= length_conv;
124 temp_params3.D *= length_conv;
125 temp_params_zbl_2.ZBLcut *= length_conv;
126 temp_params_zbl_2.ZBLexpscale *= inv_length_conv;
127 // Get the atom type indices.
128 std::map<std::string,int>::const_iterator it;
129 int i,j,k;
130 it = type_map.find(type_i);
131 if (it != type_map.end())
132 i = it->second;
133 else
134 throw runtime_error("Unknown species: " + type_i);
135 it = type_map.find(type_j);
136 if (it != type_map.end())
137 j = it->second;
138 else
139 throw runtime_error("Unknown species: " + type_j);
140 it = type_map.find(type_k);
141 if (it != type_map.end())
142 k = it->second;
143 else
144 throw runtime_error("Unknown species: " + type_k);
145 // Check if parameters were given twice.
146 if (got_interaction(i,j,k))
147 throw runtime_error("Interaction " + type_i +
148 "-" + type_j +
149 "-" + type_k +
150 " is defined twice!");
151 // All OK, store.
152 got_interaction(i,j,k) = true;
153 if (j == k) {
154 temp_params2.R = temp_params3.R;
155 temp_params2.D = temp_params3.D;
156 params2(i,j) = temp_params2;
157 params_zbl_2(i,j) = temp_params_zbl_2;
158 kim_params_zbl.Z_i(i,j) = Z_i;
159 kim_params_zbl.Z_j(i,j) = Z_j;
160 }
161 params3(i,j,k) = temp_params3;
162 kim_params.c(i,j,k) = c;
163 kim_params.d(i,j,k) = d;
164 }
165 if (!got_interaction.all())
166 throw runtime_error("Not all interactions were set!");
167 // Check parameters and pre-compute values.
168 prepare_params();
169 // Copy to KIM-published data.
170 kim_params.from_params(params2, params3);
171 kim_params_zbl.from_params(params_zbl_2);
172 }
173
update_params()174 void PairTersoffZBL::update_params() {
175 kim_params.to_params(params2, params3);
176 kim_params_zbl.to_params(params_zbl_2);
177 prepare_params();
178 }
179
prepare_params()180 void PairTersoffZBL::prepare_params() {
181 PairTersoff::prepare_params();
182
183 for (int i = 0; i != n_spec; ++i) {
184 const string type_i = to_spec.at(i);
185 for (int j = 0; j != n_spec; ++j) {
186 const string type_j = to_spec.at(j);
187 ParamsZBL2& temp_params_zbl_2 = params_zbl_2(i,j);
188 if (kim_params_zbl.Z_i(i,j) < 1)
189 throw runtime_error("Parameter Z_i ("
190 + type_i +
191 "-" + type_j +
192 ") may not be smaller than one.");
193 if (kim_params_zbl.Z_j(i,j) < 1)
194 throw runtime_error("Parameter Z_j ("
195 + type_i +
196 "-" + type_j +
197 ") may not be smaller than one.");
198 if (temp_params_zbl_2.ZBLcut < 0)
199 throw runtime_error("Parameter ZBLcut ("
200 + type_i +
201 "-" + type_j +
202 ") may not be smaller than one.");
203 if (temp_params_zbl_2.ZBLexpscale < 0)
204 throw runtime_error("Parameter ZBLexpscale ("
205 + type_i +
206 "-" + type_j +
207 ") may not be smaller than one.");
208 // Pre-compute values.
209 temp_params_zbl_2.a =
210 0.8854 * global_a_0 / (pow(kim_params_zbl.Z_i(i,j), 0.23)
211 + pow(kim_params_zbl.Z_j(i,j), 0.23));
212 temp_params_zbl_2.premult =
213 (kim_params_zbl.Z_i(i,j) * kim_params_zbl.Z_j(i,j) * global_e_sq)
214 /
215 (4.0 * pi * global_epsilon_0);
216 }
217 }
218 }
219
write_params(ofstream & outfile)220 void PairTersoffZBL::write_params(ofstream& outfile) {
221 // Set maximum precision.
222 outfile << setprecision(16);
223 // Write.
224 for (int i = 0; i < n_spec; ++i) {
225 const string type_i = to_spec.at(i);
226 for (int j = 0; j < n_spec; ++j) {
227 const string type_j = to_spec.at(j);
228 for (int k = 0; k < n_spec; ++k) {
229 const string type_k = to_spec.at(k);
230 outfile << type_i << " " << type_j << " " << type_k << " ";
231 outfile << kim_params.m(i,j,k) << " ";
232 outfile << kim_params.gamma(i,j,k) << " ";
233 outfile << kim_params.lam3(i,j,k) << " ";
234 outfile << kim_params.c(i,j,k) << " ";
235 outfile << kim_params.d(i,j,k) << " ";
236 outfile << kim_params.h(i,j,k) << " ";
237 if (j == k) {
238 // Two-body parameters are only taken from j == k, so in
239 // order to make that clear, we write zeros otherwise.
240 outfile << kim_params.n(i,j) << " ";
241 outfile << kim_params.beta(i,j) << " ";
242 outfile << kim_params.lam2(i,j) << " ";
243 outfile << kim_params.B(i,j) << " ";
244 } else {
245 outfile << "0 0 0 0 ";
246 }
247 outfile << kim_params.R(i,j,k) << " ";
248 outfile << kim_params.D(i,j,k) << " ";
249 if (j == k) {
250 // Two-body parameters are only taken from j == k, so in
251 // order to make that clear, we write zeros otherwise.
252 outfile << kim_params.lam1(i,j) << " ";
253 outfile << kim_params.A(i,j) << " ";
254 outfile << kim_params_zbl.Z_i(i,j) << " ";
255 outfile << kim_params_zbl.Z_j(i,j) << " ";
256 outfile << kim_params_zbl.ZBLcut(i,j) << " ";
257 outfile << kim_params_zbl.ZBLexpscale(i,j) << endl;
258 } else {
259 outfile << "0 0 0 0 0 0" << endl;
260 }
261 }
262 }
263 }
264 }
265
266
267 /* ---------------------------------------------------------------------- */
268
repulsive(double r,double fc,double fc_d,int itype,int jtype,bool eflag,double & eng) const269 double PairTersoffZBL::repulsive(double r, double fc, double fc_d,
270 int itype, int jtype,
271 bool eflag, double &eng) const
272 {
273 const double lam1 = params2(itype,jtype).lam1;
274 const double A = params2(itype,jtype).A;
275
276 // Tersoff repulsive portion
277
278 const double tmp_exp = exp(-lam1 * r);
279 const double eng_ters = fc * A * tmp_exp;
280 const double fforce_ters = A * tmp_exp * (fc_d - fc*lam1); // note, will be multiplied with -1/r at the end.
281
282 // ZBL repulsive portion
283
284 const double a_ij = params_zbl_2(itype,jtype).a;
285 const double premult = params_zbl_2(itype,jtype).premult;
286 const double ZBLexpscale = params_zbl_2(itype,jtype).ZBLexpscale;
287 const double ZBLcut = params_zbl_2(itype,jtype).ZBLcut;
288
289 const double r_over_a = r / a_ij;
290 const double exp1 = 0.1818 * exp(-3.2 * r_over_a);
291 const double exp2 = 0.5099 * exp(-0.9423 * r_over_a);
292 const double exp3 = 0.2802 * exp(-0.4029 * r_over_a);
293 const double exp4 = 0.02817 * exp(-0.2016 * r_over_a);
294 const double phi = exp1 + exp2 + exp3 + exp4;
295 const double dphi = (1.0/a_ij) * (-3.2 * exp1 -
296 0.9423 * exp2 -
297 0.4029 * exp3 -
298 0.2016 * exp4);
299 const double eng_ZBL = premult * (1.0/r) * phi;
300 const double fforce_ZBL = premult * -phi / pow2(r) + premult * dphi / r;
301
302 // combine two parts with smoothing by Fermi-like function
303
304 const double f_F = F_fermi(r, ZBLexpscale, ZBLcut);
305 const double f_F_d = F_fermi_d(r, ZBLexpscale, ZBLcut);
306
307 if (eflag) {
308 eng = (1.0 - f_F) * eng_ZBL + f_F * eng_ters;
309 }
310
311 return -(-f_F_d * eng_ZBL +
312 (1.0 - f_F) * fforce_ZBL +
313 f_F_d * eng_ters + f_F * fforce_ters) / r;
314 }
315
316 /* ---------------------------------------------------------------------- */
317
ters_fa(double r,double fc,int itype,int jtype) const318 double PairTersoffZBL::ters_fa(double r, double fc,
319 int itype, int jtype) const
320 {
321 if (fc == 0.0) return 0.0;
322
323 const double B = params2(itype,jtype).B;
324 const double lam2 = params2(itype,jtype).lam2;
325 const double ZBLexpscale = params_zbl_2(itype,jtype).ZBLexpscale;
326 const double ZBLcut = params_zbl_2(itype,jtype).ZBLcut;
327
328 return -B * exp(-lam2 * r) * fc * F_fermi(r, ZBLexpscale, ZBLcut);
329 }
330
331 /* ---------------------------------------------------------------------- */
332
ters_fa_d(double r,double fc,double fc_d,int itype,int jtype) const333 double PairTersoffZBL::ters_fa_d(double r, double fc, double fc_d,
334 int itype, int jtype) const
335 {
336 if (fc == 0.0) return 0.0;
337
338 const double B = params2(itype,jtype).B;
339 const double lam2 = params2(itype,jtype).lam2;
340 const double ZBLexpscale = params_zbl_2(itype,jtype).ZBLexpscale;
341 const double ZBLcut = params_zbl_2(itype,jtype).ZBLcut;
342
343 const double f_F = F_fermi(r, ZBLexpscale, ZBLcut);
344 const double f_F_d = F_fermi_d(r, ZBLexpscale, ZBLcut);
345
346 return B * exp(-lam2 * r) * (lam2 * fc * f_F - fc_d * f_F - fc * f_F_d);
347 }
348
349 /* ----------------------------------------------------------------------
350 Fermi-like smoothing function
351 ------------------------------------------------------------------------- */
352
F_fermi(double r,double ZBLexpscale,double ZBLcut) const353 double PairTersoffZBL::F_fermi(double r, double ZBLexpscale, double ZBLcut) const
354 {
355 return 1.0 / (1.0 + exp(-ZBLexpscale * (r - ZBLcut)));
356 }
357
358 /* ----------------------------------------------------------------------
359 Fermi-like smoothing function derivative with respect to r
360 ------------------------------------------------------------------------- */
361
F_fermi_d(double r,double ZBLexpscale,double ZBLcut) const362 double PairTersoffZBL::F_fermi_d(double r, double ZBLexpscale, double ZBLcut) const
363 {
364 return ZBLexpscale * exp(-ZBLexpscale * (r - ZBLcut)) /
365 pow2(1.0 + exp(-ZBLexpscale * (r - ZBLcut)));
366 }
367