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