OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
EAMAtomTypesSectionParser.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "io/EAMAtomTypesSectionParser.hpp"
46
47#include <functional>
48#include <string>
49
50#include "brains/ForceField.hpp"
51#include "types/AtomType.hpp"
52#include "types/EAMAdapter.hpp"
53#include "utils/simError.h"
54
55#ifndef FIX_UNUSED
56#define FIX_UNUSED(X) (void)(X) /* avoid warnings for unused params */
57#endif
58
59using namespace std;
60namespace OpenMD {
61
62 EAMAtomTypesSectionParser::EAMAtomTypesSectionParser(
63 ForceFieldOptions& options) :
64 options_(options) {
65 setSectionName("EAMAtomTypes");
66 }
67
68 void EAMAtomTypesSectionParser::parseLine(ForceField& ff, const string& line,
69 int lineNo) {
70 eus_ = options_.getMetallicEnergyUnitScaling();
71 dus_ = options_.getDistanceUnitScaling();
72
73 StringTokenizer tokenizer(line);
74 int nTokens = tokenizer.countTokens();
75
76 if (tokenizer.countTokens() >= 2) {
77 std::string atomTypeName = tokenizer.nextToken();
78 std::string eamParameterType = tokenizer.nextToken();
79 nTokens -= 2;
80 AtomType* atomType = ff.getAtomType(atomTypeName);
81 if (atomType != NULL) {
82 EAMAdapter ea = EAMAdapter(atomType);
83 toUpper(eamParameterType);
84
85 if (eamParameterType == "FUNCFL") {
86 std::string funcflFile = tokenizer.nextToken();
87 parseFuncflFile(ff, ea, funcflFile, atomType->getIdent());
88 } else if (eamParameterType == "ZHOU" ||
89 eamParameterType == "ZHOU2001") {
90 if (nTokens < 20) {
91 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
92 "EAMAtomTypesSectionParser Error: "
93 "Not enough tokens at line %d\n",
94 lineNo);
95 painCave.isFatal = 1;
96 simError();
97 } else {
98 std::string lattice = tokenizer.nextToken();
99 toUpper(lattice);
100 RealType re = dus_ * tokenizer.nextTokenAsDouble();
101 RealType fe = tokenizer.nextTokenAsDouble();
102 RealType rhoe = tokenizer.nextTokenAsDouble();
103 RealType alpha = tokenizer.nextTokenAsDouble();
104 RealType beta = tokenizer.nextTokenAsDouble();
105 RealType A = eus_ * tokenizer.nextTokenAsDouble();
106 RealType B = eus_ * tokenizer.nextTokenAsDouble();
107 RealType kappa = tokenizer.nextTokenAsDouble();
108 RealType lambda = tokenizer.nextTokenAsDouble();
109 std::vector<RealType> Fn;
110 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
111 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
112 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
113 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
114 std::vector<RealType> F;
115 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
116 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
117 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
118 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
119 RealType eta = tokenizer.nextTokenAsDouble();
120 RealType Fe = eus_ * tokenizer.nextTokenAsDouble();
121
122 ea.makeZhou2001(lattice, re, fe, rhoe, alpha, beta, A, B, kappa,
123 lambda, Fn, F, eta, Fe);
124 }
125 } else if (eamParameterType == "ZHOU2004") {
126 if (nTokens < 23) {
127 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
128 "EAMAtomTypesSectionParser Error: "
129 "Not enough tokens at line %d\n",
130 lineNo);
131 painCave.isFatal = 1;
132 simError();
133 } else {
134 std::string lattice = tokenizer.nextToken();
135 toUpper(lattice);
136 RealType re = dus_ * tokenizer.nextTokenAsDouble();
137 RealType fe = tokenizer.nextTokenAsDouble();
138 RealType rhoe = tokenizer.nextTokenAsDouble();
139 RealType rhos = tokenizer.nextTokenAsDouble();
140 RealType alpha = tokenizer.nextTokenAsDouble();
141 RealType beta = tokenizer.nextTokenAsDouble();
142 RealType A = eus_ * tokenizer.nextTokenAsDouble();
143 RealType B = eus_ * tokenizer.nextTokenAsDouble();
144 RealType kappa = tokenizer.nextTokenAsDouble();
145 RealType lambda = tokenizer.nextTokenAsDouble();
146 std::vector<RealType> Fn;
147 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
148 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
149 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
150 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
151 std::vector<RealType> F;
152 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
153 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
154 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
155 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
156 RealType eta = tokenizer.nextTokenAsDouble();
157 RealType Fe = eus_ * tokenizer.nextTokenAsDouble();
158 RealType rhol = tokenizer.nextTokenAsDouble();
159 RealType rhoh = tokenizer.nextTokenAsDouble();
160
161 ea.makeZhou2004(lattice, re, fe, rhoe, rhos, alpha, beta, A, B,
162 kappa, lambda, Fn, F, eta, Fe, rhol, rhoh);
163 }
164
165 } else if (eamParameterType == "ZHOU2005") {
166 if (nTokens < 22) {
167 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
168 "EAMAtomTypesSectionParser Error: "
169 "Not enough tokens at line %d\n",
170 lineNo);
171 painCave.isFatal = 1;
172 simError();
173 } else {
174 std::string lattice = tokenizer.nextToken();
175 toUpper(lattice);
176 RealType re = dus_ * tokenizer.nextTokenAsDouble();
177 RealType fe = tokenizer.nextTokenAsDouble();
178 RealType rhoe = tokenizer.nextTokenAsDouble();
179 RealType rhos = tokenizer.nextTokenAsDouble();
180 RealType alpha = tokenizer.nextTokenAsDouble();
181 RealType beta = tokenizer.nextTokenAsDouble();
182 RealType A = eus_ * tokenizer.nextTokenAsDouble();
183 RealType B = eus_ * tokenizer.nextTokenAsDouble();
184 RealType kappa = tokenizer.nextTokenAsDouble();
185 RealType lambda = tokenizer.nextTokenAsDouble();
186 std::vector<RealType> Fn;
187 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
188 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
189 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
190 Fn.push_back(eus_ * tokenizer.nextTokenAsDouble());
191 std::vector<RealType> F;
192 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
193 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
194 F.push_back(eus_ * tokenizer.nextTokenAsDouble());
195 RealType F3minus = eus_ * tokenizer.nextTokenAsDouble();
196 RealType F3plus = eus_ * tokenizer.nextTokenAsDouble();
197 RealType eta = tokenizer.nextTokenAsDouble();
198 RealType Fe = eus_ * tokenizer.nextTokenAsDouble();
199
200 ea.makeZhou2005(lattice, re, fe, rhoe, rhos, alpha, beta, A, B,
201 kappa, lambda, Fn, F, F3minus, F3plus, eta, Fe);
202 }
203 } else if (eamParameterType == "ZHOU2005OXYGEN") {
204 if (nTokens < 36) {
205 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
206 "EAMAtomTypesSectionParser Error: "
207 "Not enough tokens at line %d\n",
208 lineNo);
209 painCave.isFatal = 1;
210 simError();
211 } else {
212 RealType re = dus_ * tokenizer.nextTokenAsDouble();
213 RealType fe = tokenizer.nextTokenAsDouble();
214 RealType alpha = tokenizer.nextTokenAsDouble();
215 RealType beta = tokenizer.nextTokenAsDouble();
216 RealType A = eus_ * tokenizer.nextTokenAsDouble();
217 RealType B = eus_ * tokenizer.nextTokenAsDouble();
218 RealType kappa = tokenizer.nextTokenAsDouble();
219 RealType lambda = tokenizer.nextTokenAsDouble();
220 RealType gamma = tokenizer.nextTokenAsDouble();
221 RealType nu = tokenizer.nextTokenAsDouble();
222
223 std::vector<RealType> OrhoLimits;
224 OrhoLimits.push_back(tokenizer.nextTokenAsDouble());
225 OrhoLimits.push_back(tokenizer.nextTokenAsDouble());
226 OrhoLimits.push_back(tokenizer.nextTokenAsDouble());
227 OrhoLimits.push_back(tokenizer.nextTokenAsDouble());
228 OrhoLimits.push_back(tokenizer.nextTokenAsDouble());
229
230 std::vector<RealType> OrhoE;
231 OrhoE.push_back(tokenizer.nextTokenAsDouble());
232 OrhoE.push_back(tokenizer.nextTokenAsDouble());
233 OrhoE.push_back(tokenizer.nextTokenAsDouble());
234 OrhoE.push_back(tokenizer.nextTokenAsDouble());
235 OrhoE.push_back(tokenizer.nextTokenAsDouble());
236
237 std::vector<std::vector<RealType>> OF;
238 OF.resize(5);
239 OF[0].push_back(eus_ * tokenizer.nextTokenAsDouble());
240 OF[0].push_back(eus_ * tokenizer.nextTokenAsDouble());
241 OF[0].push_back(eus_ * tokenizer.nextTokenAsDouble());
242 OF[0].push_back(eus_ * tokenizer.nextTokenAsDouble());
243
244 OF[1].push_back(eus_ * tokenizer.nextTokenAsDouble());
245 OF[1].push_back(eus_ * tokenizer.nextTokenAsDouble());
246 OF[1].push_back(eus_ * tokenizer.nextTokenAsDouble());
247
248 OF[2].push_back(eus_ * tokenizer.nextTokenAsDouble());
249 OF[2].push_back(eus_ * tokenizer.nextTokenAsDouble());
250 OF[2].push_back(eus_ * tokenizer.nextTokenAsDouble());
251
252 OF[3].push_back(eus_ * tokenizer.nextTokenAsDouble());
253 OF[3].push_back(eus_ * tokenizer.nextTokenAsDouble());
254 OF[3].push_back(eus_ * tokenizer.nextTokenAsDouble());
255
256 OF[4].push_back(eus_ * tokenizer.nextTokenAsDouble());
257 OF[4].push_back(eus_ * tokenizer.nextTokenAsDouble());
258 OF[4].push_back(eus_ * tokenizer.nextTokenAsDouble());
259
260 ea.makeZhou2005Oxygen(re, fe, alpha, beta, A, B, kappa, lambda,
261 gamma, nu, OrhoLimits, OrhoE, OF);
262 }
263 } else if (eamParameterType == "ZHOUROSE") {
264 if (nTokens < 10) {
265 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
266 "EAMAtomTypesSectionParser Error: "
267 "Not enough tokens at line %d\n",
268 lineNo);
269 painCave.isFatal = 1;
270 simError();
271 } else {
272 RealType re = dus_ * tokenizer.nextTokenAsDouble();
273 RealType fe = tokenizer.nextTokenAsDouble();
274 RealType rhoe = tokenizer.nextTokenAsDouble();
275 RealType alpha = tokenizer.nextTokenAsDouble();
276 RealType beta = tokenizer.nextTokenAsDouble();
277 RealType A = eus_ * tokenizer.nextTokenAsDouble();
278 RealType B = eus_ * tokenizer.nextTokenAsDouble();
279 RealType kappa = tokenizer.nextTokenAsDouble();
280 RealType lambda = tokenizer.nextTokenAsDouble();
281 RealType F0 = eus_ * tokenizer.nextTokenAsDouble();
282
283 ea.makeZhouRose(re, fe, rhoe, alpha, beta, A, B, kappa, lambda, F0);
284 }
285 } else if (eamParameterType == "OXYFUNCFL") {
286 if (nTokens < 9) {
287 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
288 "EAMAtomTypesSectionParser Error: "
289 "Not enough tokens at line %d\n",
290 lineNo);
291 painCave.isFatal = 1;
292 simError();
293 } else {
294 std::string funcflFile = tokenizer.nextToken();
295
296 RealType re = dus_ * tokenizer.nextTokenAsDouble();
297 RealType fe = tokenizer.nextTokenAsDouble();
298 RealType alpha = tokenizer.nextTokenAsDouble();
299 RealType beta = tokenizer.nextTokenAsDouble();
300 RealType A = eus_ * tokenizer.nextTokenAsDouble();
301 RealType B = eus_ * tokenizer.nextTokenAsDouble();
302 RealType kappa = tokenizer.nextTokenAsDouble();
303 RealType lambda = tokenizer.nextTokenAsDouble();
304
305 ifstrstream* ppfStream = ff.openForceFieldFile(funcflFile);
306 const int bufferSize = 65535;
307 char buffer[bufferSize];
308
309 // skip first line
310 ppfStream->getline(buffer, bufferSize);
311
312 std::vector<RealType> F;
313
314 int nrho = 0;
315 RealType drho = 0.0;
316 if (ppfStream->getline(buffer, bufferSize)) {
317 StringTokenizer tokenizer1(buffer);
318
319 if (tokenizer1.countTokens() == 2) {
320 nrho = tokenizer1.nextTokenAsInt();
321 drho = tokenizer1.nextTokenAsDouble();
322 } else {
323 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
324 "EAMAtomTypesSectionParser Error: "
325 "Not enough tokens\n");
326 painCave.isFatal = 1;
327 simError();
328 }
329 }
330
331 parseEAMArray(*ppfStream, F, nrho);
332 delete ppfStream;
333
334 // Convert to eV using energy unit scaling in force field:
335 std::transform(F.begin(), F.end(), F.begin(),
336 std::bind(std::multiplies<RealType>(), eus_,
337 std::placeholders::_1));
338
339 ea.makeOxygenFuncfl(re, fe, alpha, beta, A, B, kappa, lambda, drho,
340 nrho, F);
341 }
342 } else {
343 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
344 "EAMAtomTypesSectionParser Error: %s "
345 "is not a recognized EAM type\n",
346 eamParameterType.c_str());
347 painCave.isFatal = 1;
348 simError();
349 }
350
351 } else {
352 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
353 "EAMAtomTypesSectionParser Error: "
354 "Can not find AtomType [%s]\n",
355 atomTypeName.c_str());
356 painCave.isFatal = 1;
357 simError();
358 }
359
360 } else {
361 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
362 "EAMAtomTypesSectionParser Error: "
363 "Not enough tokens at line %d\n",
364 lineNo);
365 painCave.isFatal = 1;
366 simError();
367 }
368 }
369
370 void EAMAtomTypesSectionParser::parseFuncflFile(ForceField& ff, EAMAdapter ea,
371 const string& funcflFile,
372 int) {
373 ifstrstream* ppfStream = ff.openForceFieldFile(funcflFile);
374 const int bufferSize = 65535;
375 char buffer[bufferSize];
376
377 // skip first line
378 ppfStream->getline(buffer, bufferSize);
379
380 // The second line contains atomic number, atomic mass, a lattice
381 // constant and lattice type
382 int atomicNumber;
383 RealType atomicMass;
384 RealType latticeConstant(0.0);
385 std::string lattice;
386
387 // The third line is nrho, drho, nr, dr and rcut
388 int nrho(0);
389 RealType drho(0.0);
390 int nr(0);
391 RealType dr(0.0);
392 RealType rcut(0.0);
393 std::vector<RealType> F;
394 std::vector<RealType> Z;
395 std::vector<RealType> rho;
396
397 if (ppfStream->getline(buffer, bufferSize)) {
398 StringTokenizer tokenizer1(buffer);
399
400 if (tokenizer1.countTokens() >= 4) {
401 atomicNumber = tokenizer1.nextTokenAsInt();
402 atomicMass = tokenizer1.nextTokenAsDouble();
403 latticeConstant = tokenizer1.nextTokenAsDouble() * dus_;
404 lattice = tokenizer1.nextToken();
405 } else {
406 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
407 "EAMAtomTypesSectionParser Error: "
408 "Not enough tokens\n");
409 painCave.isFatal = 1;
410 simError();
411 }
412 }
413
414 FIX_UNUSED(atomicNumber);
415 FIX_UNUSED(atomicMass);
416
417 if (ppfStream->getline(buffer, bufferSize)) {
418 StringTokenizer tokenizer2(buffer);
419
420 if (tokenizer2.countTokens() >= 5) {
421 nrho = tokenizer2.nextTokenAsInt();
422 drho = tokenizer2.nextTokenAsDouble();
423 nr = tokenizer2.nextTokenAsInt();
424 dr = tokenizer2.nextTokenAsDouble() * dus_;
425 rcut = tokenizer2.nextTokenAsDouble() * dus_;
426 } else {
427 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
428 "EAMAtomTypesSectionParser Error: "
429 "Not enough tokens\n");
430 painCave.isFatal = 1;
431 simError();
432 }
433 }
434
435 parseEAMArray(*ppfStream, F, nrho);
436 parseEAMArray(*ppfStream, Z, nr);
437 parseEAMArray(*ppfStream, rho, nr);
438
439 // Convert to kcal/mol using energy unit scaling in force field:
440 std::transform(
441 F.begin(), F.end(), F.begin(),
442 std::bind(std::multiplies<RealType>(), eus_, std::placeholders::_1));
443
444 ea.makeFuncfl(latticeConstant, lattice, nrho, drho, nr, dr, rcut, Z, rho,
445 F);
446
447 delete ppfStream;
448 }
449
450 void EAMAtomTypesSectionParser::parseEAMArray(istream& input,
451 std::vector<RealType>& array,
452 int num) {
453 const int dataPerLine = 5;
454 if (num % dataPerLine != 0) {}
455
456 int nlinesToRead = num / dataPerLine;
457
458 const int bufferSize = 65535;
459 char buffer[bufferSize];
460 int lineCount = 0;
461
462 while (lineCount < nlinesToRead && input.getline(buffer, bufferSize)) {
463 StringTokenizer tokenizer(buffer);
464 if (tokenizer.countTokens() >= dataPerLine) {
465 for (int i = 0; i < dataPerLine; ++i) {
466 array.push_back(tokenizer.nextTokenAsDouble());
467 }
468 } else {
469 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
470 "EAMAtomTypesSectionParser Error: "
471 "Not enough tokens\n");
472 painCave.isFatal = 1;
473 simError();
474 }
475 ++lineCount;
476 }
477
478 if (lineCount < nlinesToRead) {
479 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
480 "EAMAtomTypesSectionParser Error: "
481 "Not enough lines to read\n");
482 painCave.isFatal = 1;
483 simError();
484 }
485 }
486
487} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.