OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
NonBondedInteractionsSectionParser.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 following paper when you publish your work:
33 *
34 * [1] Drisko et al., J. Open Source Softw. 9, 7004 (2024).
35 *
36 * Good starting points for code and simulation methodology are:
37 *
38 * [2] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
39 * [3] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
40 * [4] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
41 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
42 * [6] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
43 * [7] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
44 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
45 * [9] Drisko & Gezelter, J. Chem. Theory Comput. 20, 4986-4997 (2024).
46 */
47
48#include "io/NonBondedInteractionsSectionParser.hpp"
49
50#include "brains/ForceField.hpp"
51#include "types/AtomType.hpp"
52#include "types/BuckinghamInteractionType.hpp"
53#include "types/EAMInteractionType.hpp"
54#include "types/InversePowerSeriesInteractionType.hpp"
55#include "types/LennardJonesInteractionType.hpp"
56#include "types/MAWInteractionType.hpp"
57#include "types/MieInteractionType.hpp"
58#include "types/MorseInteractionType.hpp"
59#include "types/RepulsivePowerInteractionType.hpp"
60#include "utils/simError.h"
61
62namespace OpenMD {
63
64 NonBondedInteractionsSectionParser::NonBondedInteractionsSectionParser(
65 ForceFieldOptions& options) :
66 options_(options) {
67 setSectionName("NonBondedInteractions");
68
69 stringToEnumMap_["MAW"] = MAW;
70 stringToEnumMap_["ShiftedMorse"] = ShiftedMorse;
71 stringToEnumMap_["LennardJones"] = LennardJones;
72 stringToEnumMap_["RepulsiveMorse"] = RepulsiveMorse;
73 stringToEnumMap_["RepulsivePower"] = RepulsivePower;
74 stringToEnumMap_["Mie"] = Mie;
75 stringToEnumMap_["Buckingham"] = Buckingham;
76 stringToEnumMap_["EAMTable"] = EAMTable;
77 stringToEnumMap_["EAMZhou"] = EAMZhou;
78 stringToEnumMap_["EAMOxides"] = EAMOxides;
79 stringToEnumMap_["InversePowerSeries"] = InversePowerSeries;
80 }
81
82 void NonBondedInteractionsSectionParser::parseLine(ForceField& ff,
83 const std::string& line,
84 int lineNo) {
85 StringTokenizer tokenizer(line);
86 int nTokens = tokenizer.countTokens();
87 if (nTokens < 3) {
88 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
89 "NonBondedInteractionsSectionParser Error: Not enough tokens at "
90 "line %d\n",
91 lineNo);
92 painCave.isFatal = 1;
93 simError();
94 }
95
96 meus_ = options_.getMetallicEnergyUnitScaling();
97 eus_ = options_.getEnergyUnitScaling();
98 dus_ = options_.getDistanceUnitScaling();
99
100 std::string at1 = tokenizer.nextToken();
101 std::string at2 = tokenizer.nextToken();
102 std::string itype = tokenizer.nextToken();
103
104 NonBondedInteractionTypeEnum nbit = getNonBondedInteractionTypeEnum(itype);
105 nTokens -= 3;
106 NonBondedInteractionType* interactionType = NULL;
107
108 // switch is a nightmare to maintain
109 switch (nbit) {
110 case MAW:
111 if (nTokens != 5) {
112 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
113 "NonBondedInteractionsSectionParser Error: Token number "
114 "mismatch at line "
115 "%d. 8 tokens expected. \n",
116 lineNo);
117 painCave.isFatal = 1;
118 simError();
119 } else {
120 RealType r_e = dus_ * tokenizer.nextTokenAsDouble();
121 RealType D_e = eus_ * tokenizer.nextTokenAsDouble();
122 RealType beta = tokenizer.nextTokenAsDouble() / dus_;
123 RealType ca1 = tokenizer.nextTokenAsDouble();
124 RealType cb1 = tokenizer.nextTokenAsDouble();
125 interactionType = new MAWInteractionType(D_e, beta, r_e, ca1, cb1);
126 }
127 break;
128
129 case ShiftedMorse:
130 if (nTokens != 3) {
131 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
132 "NonBondedInteractionsSectionParser Error: Token number "
133 "mismatch at line "
134 "%d. 6 tokens expected. \n",
135 lineNo);
136 painCave.isFatal = 1;
137 simError();
138 } else {
139 RealType r0 = dus_ * tokenizer.nextTokenAsDouble();
140 RealType D0 = eus_ * tokenizer.nextTokenAsDouble();
141 RealType beta0 = tokenizer.nextTokenAsDouble() / dus_;
142 interactionType = new MorseInteractionType(D0, beta0, r0, mtShifted);
143 }
144 break;
145
146 case RepulsiveMorse:
147 if (nTokens != 3) {
148 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
149 "NonBondedInteractionsSectionParser Error: Token number "
150 "mismatch at line "
151 "%d. 6 tokens expected. \n",
152 lineNo);
153 painCave.isFatal = 1;
154 simError();
155 } else {
156 RealType r0 = dus_ * tokenizer.nextTokenAsDouble();
157 RealType D0 = eus_ * tokenizer.nextTokenAsDouble();
158 RealType beta0 = tokenizer.nextTokenAsDouble() / dus_;
159 interactionType = new MorseInteractionType(D0, beta0, r0, mtRepulsive);
160 }
161 break;
162
163 case LennardJones:
164 if (nTokens != 2) {
165 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
166 "NonBondedInteractionsSectionParser Error: Token number "
167 "mismatch at line "
168 "%d. 5 tokens expected. \n",
169 lineNo);
170 painCave.isFatal = 1;
171 simError();
172 } else {
173 RealType sigma = dus_ * tokenizer.nextTokenAsDouble();
174 RealType epsilon = eus_ * tokenizer.nextTokenAsDouble();
175 interactionType = new LennardJonesInteractionType(sigma, epsilon);
176 }
177 break;
178
179 case RepulsivePower:
180 if (nTokens < 3) {
181 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
182 "NonBondedInteractionsSectionParser Error: Token number "
183 "mismatch at line "
184 "%d. 6 tokens expected. \n",
185 lineNo);
186 painCave.isFatal = 1;
187 simError();
188 } else {
189 RealType sigma = dus_ * tokenizer.nextTokenAsDouble();
190 RealType epsilon = eus_ * tokenizer.nextTokenAsDouble();
191 int nRep = tokenizer.nextTokenAsInt();
192 interactionType =
193 new RepulsivePowerInteractionType(sigma, epsilon, nRep);
194 }
195 break;
196
197 case Mie:
198 if (nTokens != 4) {
199 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
200 "NonBondedInteractionsSectionParser Error: Token number "
201 "mismatch at line "
202 "%d. 7 tokens expected. \n",
203 lineNo);
204 painCave.isFatal = 1;
205 simError();
206 } else {
207 RealType sigma = dus_ * tokenizer.nextTokenAsDouble();
208 RealType epsilon = eus_ * tokenizer.nextTokenAsDouble();
209 int nRep = tokenizer.nextTokenAsInt();
210 int mAtt = tokenizer.nextTokenAsInt();
211 interactionType = new MieInteractionType(sigma, epsilon, nRep, mAtt);
212 }
213 break;
214
215 case Buckingham:
216 if (nTokens < 4) {
217 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
218 "NonBondedInteractionsSectionParser Error: Not enough tokens "
219 "at line %d\n",
220 lineNo);
221 painCave.isFatal = 1;
222 simError();
223 } else {
224 std::string btype = tokenizer.nextToken();
225 toUpper(btype);
226
227 RealType A = eus_ * tokenizer.nextTokenAsDouble();
228 RealType B = tokenizer.nextTokenAsDouble() / dus_;
229 RealType C =
230 tokenizer.nextTokenAsDouble(); // should also have a scaling
231 RealType sigma = 0.0;
232 RealType epsilon = 0.0;
233
234 if (btype.compare("MODIFIED")) {
235 sigma = dus_ * tokenizer.nextTokenAsDouble();
236 epsilon = eus_ * tokenizer.nextTokenAsDouble();
237 interactionType = new BuckinghamInteractionType(A, B, C, sigma,
238 epsilon, btModified);
239
240 } else if (btype.compare("TRADITIONAL")) {
241 interactionType =
242 new BuckinghamInteractionType(A, B, C, btTraditional);
243 } else {
244 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
245 "NonBondedInteractionsSectionParser Error: Unknown "
246 "Buckingham Type at "
247 "line %d\n",
248 lineNo);
249 painCave.isFatal = 1;
250 simError();
251 }
252 }
253 break;
254
255 case EAMZhou:
256 if (nTokens != 7) {
257 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
258 "NonBondedInteractionsSectionParser Error: Token number "
259 "mismatch at line "
260 "%d. 10 tokens expected. \n",
261 lineNo);
262 painCave.isFatal = 1;
263 simError();
264 } else {
265 RealType re = dus_ * tokenizer.nextTokenAsDouble();
266 RealType alpha = tokenizer.nextTokenAsDouble();
267 RealType beta = tokenizer.nextTokenAsDouble();
268 // Because EAM is a metallic potential, we'll use the metallic
269 // unit scaling for these two parameters
270 RealType A = meus_ * tokenizer.nextTokenAsDouble();
271 RealType B = meus_ * tokenizer.nextTokenAsDouble();
272 RealType kappa = tokenizer.nextTokenAsDouble();
273 RealType lambda = tokenizer.nextTokenAsDouble();
274
275 interactionType =
276 new EAMInteractionType(re, alpha, beta, A, B, kappa, lambda);
277 }
278 break;
279
280 case EAMOxides:
281 if (nTokens != 5) {
282 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
283 "NonBondedInteractionsSectionParser Error: Token number "
284 "mismatch at line "
285 "%d. 8 tokens expected. \n",
286 lineNo);
287 painCave.isFatal = 1;
288 simError();
289 } else {
290 RealType re = dus_ * tokenizer.nextTokenAsDouble();
291 RealType alpha = tokenizer.nextTokenAsDouble();
292 // Because EAM is a metallic potential, we'll use the metallic
293 // unit scaling for these two parameters
294 RealType A = meus_ * tokenizer.nextTokenAsDouble();
295 RealType Ci = tokenizer.nextTokenAsDouble();
296 RealType Cj = tokenizer.nextTokenAsDouble();
297
298 interactionType = new EAMInteractionType(re, alpha, A, Ci, Cj);
299 }
300 break;
301
302 case InversePowerSeries:
303 if (nTokens < 2 || nTokens % 2 != 0) {
304 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
305 "NonBondedInteractionsSectionParser Error: Not enough tokens "
306 "at line %d\n",
307 lineNo);
308 painCave.isFatal = 1;
309 simError();
310 } else {
311 std::vector<std::pair<int, RealType>> series;
312 int nPairs = nTokens / 2;
313 int power;
314 RealType coefficient;
315
316 for (int i = 0; i < nPairs; ++i) {
317 power = tokenizer.nextTokenAsInt();
318 coefficient = tokenizer.nextTokenAsDouble() * eus_ * pow(dus_, power);
319 series.push_back(std::make_pair(power, coefficient));
320 }
321 interactionType = new InversePowerSeriesInteractionType(series);
322 }
323 break;
324
325 case Unknown:
326 default:
327 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
328 "NonBondedInteractionsSectionParser Error: Unknown Interaction "
329 "Type at "
330 "line %d\n",
331 lineNo);
332 painCave.isFatal = 1;
333 simError();
334
335 break;
336 }
337
338 if (interactionType != NULL) {
339 ff.addNonBondedInteractionType(at1, at2, interactionType);
340 }
341 }
342
343 NonBondedInteractionsSectionParser::NonBondedInteractionTypeEnum
344 NonBondedInteractionsSectionParser::getNonBondedInteractionTypeEnum(
345 const std::string& str) {
346 std::map<std::string, NonBondedInteractionTypeEnum>::iterator i;
347 i = stringToEnumMap_.find(str);
348
349 return i == stringToEnumMap_.end() ? Unknown : i->second;
350 }
351
352} // namespace OpenMD
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:63