OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
ForceField.cpp
Go to the documentation of this file.
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/**
49 * @file ForceField.cpp
50 * @author tlin
51 * @date 11/04/2004
52 * @version 1.0
53 */
54
55#include "brains/ForceField.hpp"
56
57#include <algorithm>
58#include <tuple>
59
60#include "io/AtomTypesSectionParser.hpp"
61#include "io/BaseAtomTypesSectionParser.hpp"
62#include "io/BendTypesSectionParser.hpp"
63#include "io/BondTypesSectionParser.hpp"
64#include "io/ChargeAtomTypesSectionParser.hpp"
65#include "io/DirectionalAtomTypesSectionParser.hpp"
66#include "io/EAMAtomTypesSectionParser.hpp"
67#include "io/FluctuatingChargeAtomTypesSectionParser.hpp"
68#include "io/GayBerneAtomTypesSectionParser.hpp"
69#include "io/InversionTypesSectionParser.hpp"
70#include "io/LennardJonesAtomTypesSectionParser.hpp"
71#include "io/MultipoleAtomTypesSectionParser.hpp"
72#include "io/NonBondedInteractionsSectionParser.hpp"
73#include "io/OptionSectionParser.hpp"
74#include "io/PolarizableAtomTypesSectionParser.hpp"
75#include "io/SCAtomTypesSectionParser.hpp"
76#include "io/ShapeAtomTypesSectionParser.hpp"
77#include "io/StickyAtomTypesSectionParser.hpp"
78#include "io/StickyPowerAtomTypesSectionParser.hpp"
79#include "io/TorsionTypesSectionParser.hpp"
80#include "io/UFFAtomTypesSectionParser.hpp"
81#include "types/EAMAdapter.hpp"
82#include "types/GayBerneAdapter.hpp"
83#include "types/LennardJonesAdapter.hpp"
84#include "types/StickyAdapter.hpp"
85#include "types/SuttonChenAdapter.hpp"
86#include "utils/simError.h"
87
88namespace OpenMD {
89
90 ForceField::ForceField(std::string ffName) : wildCardAtomTypeName_("X") {
91 char* tempPath;
92 tempPath = getenv("FORCE_PARAM_PATH");
93
94 if (tempPath == NULL) {
95 // convert a macro from compiler to a string in c++
96 STR_DEFINE(ffPath_, FRC_PATH);
97 } else {
98 ffPath_ = tempPath;
99 }
100
101 setForceFieldFileName(ffName + ".frc");
102
103 /**
104 * The order of adding section parsers is important.
105 *
106 * OptionSectionParser must come first to set options for other
107 * parsers
108 *
109 * DirectionalAtomTypesSectionParser should be added before
110 * AtomTypesSectionParser, and these two section parsers will
111 * actually create "real" AtomTypes (AtomTypesSectionParser will
112 * create AtomType and DirectionalAtomTypesSectionParser will
113 * create DirectionalAtomType, which is a subclass of AtomType and
114 * should come first).
115 *
116 * Other AtomTypes Section Parsers will not create the "real"
117 * AtomType, they only add and set some attributes of the AtomType
118 * (via the Adapters). Thus ordering of these is not important.
119 * AtomTypesSectionParser should be added before other atom type
120 *
121 * The order of BondTypesSectionParser, BendTypesSectionParser and
122 * TorsionTypesSectionParser, etc. are not important.
123 */
124
125 spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
126 spMan_.push_back(new BaseAtomTypesSectionParser());
127 spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
128 spMan_.push_back(new AtomTypesSectionParser());
129
130 spMan_.push_back(
131 new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
132 spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
133 spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
134 spMan_.push_back(
135 new FluctuatingChargeAtomTypesSectionParser(forceFieldOptions_));
136 spMan_.push_back(new PolarizableAtomTypesSectionParser(forceFieldOptions_));
137 spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
138 spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
139 spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
140 spMan_.push_back(new UFFAtomTypesSectionParser(forceFieldOptions_));
141 spMan_.push_back(new ShapeAtomTypesSectionParser(forceFieldOptions_));
142 spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
143 spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
144
145 spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
146 spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
147 spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
148 spMan_.push_back(new InversionTypesSectionParser(forceFieldOptions_));
149
150 spMan_.push_back(
151 new NonBondedInteractionsSectionParser(forceFieldOptions_));
152 }
153
154 void ForceField::parse(const std::string& filename) {
155 ifstrstream* ffStream;
156
157 ffStream = openForceFieldFile(filename);
158
159 spMan_.parse(*ffStream, *this);
160
161 ForceField::AtomTypeContainer::MapTypeIterator i;
162 AtomType* at;
163
164 for (at = atomTypeCont_.beginType(i); at != NULL;
165 at = atomTypeCont_.nextType(i)) {
166 // useBase sets the responsibilities, and these have to be done
167 // after the atomTypes and Base types have all been scanned:
168
169 std::vector<AtomType*> ayb = at->allYourBase();
170 if (ayb.size() > 1) {
171 for (int j = ayb.size() - 1; j > 0; j--) {
172 ayb[j - 1]->useBase(ayb[j]);
173 }
174 }
175 }
176
177 delete ffStream;
178 }
179
180 /**
181 * getAtomType by string
182 *
183 * finds the requested atom type in this force field using the string
184 * name of the atom type.
185 */
186 AtomType* ForceField::getAtomType(const std::string& at) {
187 std::vector<std::string> keys;
188 keys.push_back(at);
189 return atomTypeCont_.find(keys);
190 }
191
192 /**
193 * getAtomType by ident
194 *
195 * finds the requested atom type in this force field using the
196 * integer ident instead of the string name of the atom type.
197 */
198 AtomType* ForceField::getAtomType(int ident) {
199 std::string at = atypeIdentToName.find(ident)->second;
200 return getAtomType(at);
201 }
202
203 BondType* ForceField::getBondType(const std::string& at1,
204 const std::string& at2) {
205 std::vector<std::string> keys;
206 keys.push_back(at1);
207 keys.push_back(at2);
208
209 // try exact match first
210 BondType* bondType = bondTypeCont_.find(keys);
211 if (bondType) {
212 return bondType;
213 } else {
214 AtomType* atype1;
215 AtomType* atype2;
216 std::vector<std::string> at1key;
217 at1key.push_back(at1);
218 atype1 = atomTypeCont_.find(at1key);
219
220 std::vector<std::string> at2key;
221 at2key.push_back(at2);
222 atype2 = atomTypeCont_.find(at2key);
223
224 // query atom types for their chains of responsibility
225 std::vector<AtomType*> at1Chain = atype1->allYourBase();
226 std::vector<AtomType*> at2Chain = atype2->allYourBase();
227
228 std::vector<AtomType*>::iterator i;
229 std::vector<AtomType*>::iterator j;
230
231 int ii = 0;
232 int jj = 0;
233 int bondTypeScore;
234
235 std::vector<std::pair<int, std::vector<std::string>>> foundBonds;
236
237 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
238 jj = 0;
239 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
240 bondTypeScore = ii + jj;
241
242 std::vector<std::string> myKeys;
243 myKeys.push_back((*i)->getName());
244 myKeys.push_back((*j)->getName());
245
246 BondType* bondType = bondTypeCont_.find(myKeys);
247 if (bondType) {
248 foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
249 }
250 jj++;
251 }
252 ii++;
253 }
254
255 if (!foundBonds.empty()) {
256 // sort the foundBonds by the score:
257 std::sort(foundBonds.begin(), foundBonds.end());
258
259 std::vector<std::string> theKeys = foundBonds[0].second;
260
261 BondType* bestType = bondTypeCont_.find(theKeys);
262
263 return bestType;
264 } else {
265 // if no exact match found, try wild card match
266 return bondTypeCont_.find(keys, wildCardAtomTypeName_);
267 }
268 }
269 }
270
271 BendType* ForceField::getBendType(const std::string& at1,
272 const std::string& at2,
273 const std::string& at3) {
274 std::vector<std::string> keys;
275 keys.push_back(at1);
276 keys.push_back(at2);
277 keys.push_back(at3);
278
279 // try exact match first
280 BendType* bendType = bendTypeCont_.find(keys);
281 if (bendType) {
282 return bendType;
283 } else {
284 AtomType* atype1;
285 AtomType* atype2;
286 AtomType* atype3;
287 std::vector<std::string> at1key;
288 at1key.push_back(at1);
289 atype1 = atomTypeCont_.find(at1key);
290
291 std::vector<std::string> at2key;
292 at2key.push_back(at2);
293 atype2 = atomTypeCont_.find(at2key);
294
295 std::vector<std::string> at3key;
296 at3key.push_back(at3);
297 atype3 = atomTypeCont_.find(at3key);
298
299 // query atom types for their chains of responsibility
300 std::vector<AtomType*> at1Chain = atype1->allYourBase();
301 std::vector<AtomType*> at2Chain = atype2->allYourBase();
302 std::vector<AtomType*> at3Chain = atype3->allYourBase();
303
304 std::vector<AtomType*>::iterator i;
305 std::vector<AtomType*>::iterator j;
306 std::vector<AtomType*>::iterator k;
307
308 int ii = 0;
309 int jj = 0;
310 int kk = 0;
311 int IKscore;
312
313 std::vector<std::tuple<int, int, std::vector<std::string>>> foundBends;
314
315 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
316 ii = 0;
317 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
318 kk = 0;
319 for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
320 IKscore = ii + kk;
321
322 std::vector<std::string> myKeys;
323 myKeys.push_back((*i)->getName());
324 myKeys.push_back((*j)->getName());
325 myKeys.push_back((*k)->getName());
326
327 BendType* bendType = bendTypeCont_.find(myKeys);
328 if (bendType) {
329 foundBends.push_back(std::make_tuple(jj, IKscore, myKeys));
330 }
331 kk++;
332 }
333 ii++;
334 }
335 jj++;
336 }
337
338 if (!foundBends.empty()) {
339 std::sort(foundBends.begin(), foundBends.end());
340 std::vector<std::string> theKeys = std::get<2>(foundBends[0]);
341
342 BendType* bestType = bendTypeCont_.find(theKeys);
343 return bestType;
344 } else {
345 // if no exact match found, try wild card match
346 return bendTypeCont_.find(keys, wildCardAtomTypeName_);
347 }
348 }
349 }
350
351 TorsionType* ForceField::getTorsionType(const std::string& at1,
352 const std::string& at2,
353 const std::string& at3,
354 const std::string& at4) {
355 std::vector<std::string> keys;
356 keys.push_back(at1);
357 keys.push_back(at2);
358 keys.push_back(at3);
359 keys.push_back(at4);
360
361 // try exact match first
362 TorsionType* torsionType = torsionTypeCont_.find(keys);
363 if (torsionType) {
364 return torsionType;
365 } else {
366 AtomType* atype1;
367 AtomType* atype2;
368 AtomType* atype3;
369 AtomType* atype4;
370 std::vector<std::string> at1key;
371 at1key.push_back(at1);
372 atype1 = atomTypeCont_.find(at1key);
373
374 std::vector<std::string> at2key;
375 at2key.push_back(at2);
376 atype2 = atomTypeCont_.find(at2key);
377
378 std::vector<std::string> at3key;
379 at3key.push_back(at3);
380 atype3 = atomTypeCont_.find(at3key);
381
382 std::vector<std::string> at4key;
383 at4key.push_back(at4);
384 atype4 = atomTypeCont_.find(at4key);
385
386 // query atom types for their chains of responsibility
387 std::vector<AtomType*> at1Chain = atype1->allYourBase();
388 std::vector<AtomType*> at2Chain = atype2->allYourBase();
389 std::vector<AtomType*> at3Chain = atype3->allYourBase();
390 std::vector<AtomType*> at4Chain = atype4->allYourBase();
391
392 std::vector<AtomType*>::iterator i;
393 std::vector<AtomType*>::iterator j;
394 std::vector<AtomType*>::iterator k;
395 std::vector<AtomType*>::iterator l;
396
397 int ii = 0;
398 int jj = 0;
399 int kk = 0;
400 int ll = 0;
401 int ILscore;
402 int JKscore;
403
404 std::vector<std::tuple<int, int, std::vector<std::string>>> foundTorsions;
405
406 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
407 kk = 0;
408 for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
409 ii = 0;
410 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
411 ll = 0;
412 for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
413 ILscore = ii + ll;
414 JKscore = jj + kk;
415
416 std::vector<std::string> myKeys;
417 myKeys.push_back((*i)->getName());
418 myKeys.push_back((*j)->getName());
419 myKeys.push_back((*k)->getName());
420 myKeys.push_back((*l)->getName());
421
422 TorsionType* torsionType = torsionTypeCont_.find(myKeys);
423 if (torsionType) {
424 foundTorsions.push_back(
425 std::make_tuple(JKscore, ILscore, myKeys));
426 }
427 ll++;
428 }
429 ii++;
430 }
431 kk++;
432 }
433 jj++;
434 }
435
436 if (!foundTorsions.empty()) {
437 std::sort(foundTorsions.begin(), foundTorsions.end());
438 std::vector<std::string> theKeys = std::get<2>(foundTorsions[0]);
439
440 TorsionType* bestType = torsionTypeCont_.find(theKeys);
441 return bestType;
442 } else {
443 // if no exact match found, try wild card match
444
445 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
446 kk = 0;
447 for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
448 ii = 0;
449 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
450 ll = 0;
451 for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
452 ILscore = ii + ll;
453 JKscore = jj + kk;
454
455 std::vector<std::string> myKeys;
456 myKeys.push_back((*i)->getName());
457 myKeys.push_back((*j)->getName());
458 myKeys.push_back((*k)->getName());
459 myKeys.push_back((*l)->getName());
460
461 TorsionType* torsionType =
462 torsionTypeCont_.find(myKeys, wildCardAtomTypeName_);
463 if (torsionType) {
464 foundTorsions.push_back(
465 std::make_tuple(JKscore, ILscore, myKeys));
466 }
467 ll++;
468 }
469 ii++;
470 }
471 kk++;
472 }
473 jj++;
474 }
475
476 if (!foundTorsions.empty()) {
477 std::sort(foundTorsions.begin(), foundTorsions.end());
478 std::vector<std::string> theKeys = std::get<2>(foundTorsions[0]);
479 TorsionType* bestType =
480 torsionTypeCont_.find(theKeys, wildCardAtomTypeName_);
481
482 return bestType;
483 } else {
484 return NULL;
485 }
486 }
487 }
488 }
489
490 InversionType* ForceField::getInversionType(const std::string& at1,
491 const std::string& at2,
492 const std::string& at3,
493 const std::string& at4) {
494 std::vector<std::string> keys;
495 keys.push_back(at1);
496 keys.push_back(at2);
497 keys.push_back(at3);
498 keys.push_back(at4);
499
500 // try exact match first
501 InversionType* inversionType =
502 inversionTypeCont_.permutedFindSkippingFirstElement(keys);
503 if (inversionType) {
504 return inversionType;
505 } else {
506 AtomType* atype1;
507 AtomType* atype2;
508 AtomType* atype3;
509 AtomType* atype4;
510 std::vector<std::string> at1key;
511 at1key.push_back(at1);
512 atype1 = atomTypeCont_.find(at1key);
513
514 std::vector<std::string> at2key;
515 at2key.push_back(at2);
516 atype2 = atomTypeCont_.find(at2key);
517
518 std::vector<std::string> at3key;
519 at3key.push_back(at3);
520 atype3 = atomTypeCont_.find(at3key);
521
522 std::vector<std::string> at4key;
523 at4key.push_back(at4);
524 atype4 = atomTypeCont_.find(at4key);
525
526 // query atom types for their chains of responsibility
527 std::vector<AtomType*> at1Chain = atype1->allYourBase();
528 std::vector<AtomType*> at2Chain = atype2->allYourBase();
529 std::vector<AtomType*> at3Chain = atype3->allYourBase();
530 std::vector<AtomType*> at4Chain = atype4->allYourBase();
531
532 std::vector<AtomType*>::iterator i;
533 std::vector<AtomType*>::iterator j;
534 std::vector<AtomType*>::iterator k;
535 std::vector<AtomType*>::iterator l;
536
537 int ii = 0;
538 int jj = 0;
539 int kk = 0;
540 int ll = 0;
541 int Iscore;
542 int JKLscore;
543
544 std::vector<std::tuple<int, int, std::vector<std::string>>>
545 foundInversions;
546
547 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
548 kk = 0;
549 for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
550 ii = 0;
551 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
552 ll = 0;
553 for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
554 Iscore = ii;
555 JKLscore = jj + kk + ll;
556
557 std::vector<std::string> myKeys;
558 myKeys.push_back((*i)->getName());
559 myKeys.push_back((*j)->getName());
560 myKeys.push_back((*k)->getName());
561 myKeys.push_back((*l)->getName());
562
563 InversionType* inversionType =
564 inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
565 if (inversionType) {
566 foundInversions.push_back(
567 std::make_tuple(Iscore, JKLscore, myKeys));
568 }
569 ll++;
570 }
571 ii++;
572 }
573 kk++;
574 }
575 jj++;
576 }
577
578 if (!foundInversions.empty()) {
579 std::sort(foundInversions.begin(), foundInversions.end());
580 std::vector<std::string> theKeys = std::get<2>(foundInversions[0]);
581
582 InversionType* bestType =
583 inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
584
585 return bestType;
586 } else {
587 // if no exact match found, try wild card match
588 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
589 kk = 0;
590 for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
591 ii = 0;
592 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
593 ll = 0;
594 for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
595 Iscore = ii;
596 JKLscore = jj + kk + ll;
597
598 std::vector<std::string> myKeys;
599 myKeys.push_back((*i)->getName());
600 myKeys.push_back((*j)->getName());
601 myKeys.push_back((*k)->getName());
602 myKeys.push_back((*l)->getName());
603 InversionType* inversionType =
604 inversionTypeCont_.permutedFindSkippingFirstElement(
605 myKeys, wildCardAtomTypeName_);
606 if (inversionType) {
607 foundInversions.push_back(
608 std::make_tuple(Iscore, JKLscore, myKeys));
609 }
610 ll++;
611 }
612 ii++;
613 }
614 kk++;
615 }
616 jj++;
617 }
618 if (!foundInversions.empty()) {
619 std::sort(foundInversions.begin(), foundInversions.end());
620 std::vector<std::string> theKeys = std::get<2>(foundInversions[0]);
621 InversionType* bestType =
622 inversionTypeCont_.permutedFindSkippingFirstElement(
623 theKeys, wildCardAtomTypeName_);
624 return bestType;
625 } else {
626 return NULL;
627 }
628 }
629 }
630 }
631
632 NonBondedInteractionType* ForceField::getNonBondedInteractionType(
633 const std::string& at1, const std::string& at2) {
634 std::vector<std::string> keys;
635 keys.push_back(at1);
636 keys.push_back(at2);
637
638 // try exact match first
639 NonBondedInteractionType* nbiType =
640 nonBondedInteractionTypeCont_.find(keys);
641 if (nbiType) {
642 return nbiType;
643 } else {
644 AtomType* atype1;
645 AtomType* atype2;
646 std::vector<std::string> at1key;
647 at1key.push_back(at1);
648 atype1 = atomTypeCont_.find(at1key);
649
650 std::vector<std::string> at2key;
651 at2key.push_back(at2);
652 atype2 = atomTypeCont_.find(at2key);
653
654 // query atom types for their chains of responsibility
655 std::vector<AtomType*> at1Chain = atype1->allYourBase();
656 std::vector<AtomType*> at2Chain = atype2->allYourBase();
657
658 std::vector<AtomType*>::iterator i;
659 std::vector<AtomType*>::iterator j;
660
661 int ii = 0;
662 int jj = 0;
663 int nbiTypeScore;
664
665 std::vector<std::pair<int, std::vector<std::string>>> foundNBI;
666
667 for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
668 jj = 0;
669 for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
670 nbiTypeScore = ii + jj;
671
672 std::vector<std::string> myKeys;
673 myKeys.push_back((*i)->getName());
674 myKeys.push_back((*j)->getName());
675
676 NonBondedInteractionType* nbiType =
677 nonBondedInteractionTypeCont_.find(myKeys);
678 if (nbiType) {
679 foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
680 }
681 jj++;
682 }
683 ii++;
684 }
685
686 if (!foundNBI.empty()) {
687 // sort the foundNBI by the score:
688 std::sort(foundNBI.begin(), foundNBI.end());
689 std::vector<std::string> theKeys = foundNBI[0].second;
690
691 NonBondedInteractionType* bestType =
692 nonBondedInteractionTypeCont_.find(theKeys);
693 return bestType;
694 } else {
695 // if no exact match found, try wild card match
696 return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
697 }
698 }
699 }
700
701 BondType* ForceField::getExactBondType(const std::string& at1,
702 const std::string& at2) {
703 std::vector<std::string> keys;
704 keys.push_back(at1);
705 keys.push_back(at2);
706 return bondTypeCont_.find(keys);
707 }
708
709 BendType* ForceField::getExactBendType(const std::string& at1,
710 const std::string& at2,
711 const std::string& at3) {
712 std::vector<std::string> keys;
713 keys.push_back(at1);
714 keys.push_back(at2);
715 keys.push_back(at3);
716 return bendTypeCont_.find(keys);
717 }
718
719 TorsionType* ForceField::getExactTorsionType(const std::string& at1,
720 const std::string& at2,
721 const std::string& at3,
722 const std::string& at4) {
723 std::vector<std::string> keys;
724 keys.push_back(at1);
725 keys.push_back(at2);
726 keys.push_back(at3);
727 keys.push_back(at4);
728 return torsionTypeCont_.find(keys);
729 }
730
731 InversionType* ForceField::getExactInversionType(const std::string& at1,
732 const std::string& at2,
733 const std::string& at3,
734 const std::string& at4) {
735 std::vector<std::string> keys;
736 keys.push_back(at1);
737 keys.push_back(at2);
738 keys.push_back(at3);
739 keys.push_back(at4);
740 return inversionTypeCont_.find(keys);
741 }
742
743 NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(
744 const std::string& at1, const std::string& at2) {
745 std::vector<std::string> keys;
746 keys.push_back(at1);
747 keys.push_back(at2);
748 return nonBondedInteractionTypeCont_.find(keys);
749 }
750
751 bool ForceField::addAtomType(const std::string& at, AtomType* atomType) {
752 std::vector<std::string> keys;
753 keys.push_back(at);
754 atypeIdentToName[atomType->getIdent()] = at;
755 return atomTypeCont_.add(keys, atomType);
756 }
757
758 bool ForceField::replaceAtomType(const std::string& at, AtomType* atomType) {
759 std::vector<std::string> keys;
760 keys.push_back(at);
761 atypeIdentToName[atomType->getIdent()] = at;
762 return atomTypeCont_.replace(keys, atomType);
763 }
764
765 bool ForceField::addBondType(const std::string& at1, const std::string& at2,
766 BondType* bondType) {
767 std::vector<std::string> keys;
768 keys.push_back(at1);
769 keys.push_back(at2);
770 return bondTypeCont_.add(keys, bondType);
771 }
772
773 bool ForceField::addBendType(const std::string& at1, const std::string& at2,
774 const std::string& at3, BendType* bendType) {
775 std::vector<std::string> keys;
776 keys.push_back(at1);
777 keys.push_back(at2);
778 keys.push_back(at3);
779 return bendTypeCont_.add(keys, bendType);
780 }
781
782 bool ForceField::addTorsionType(const std::string& at1,
783 const std::string& at2,
784 const std::string& at3,
785 const std::string& at4,
786 TorsionType* torsionType) {
787 std::vector<std::string> keys;
788 keys.push_back(at1);
789 keys.push_back(at2);
790 keys.push_back(at3);
791 keys.push_back(at4);
792 return torsionTypeCont_.add(keys, torsionType);
793 }
794
795 bool ForceField::addInversionType(const std::string& at1,
796 const std::string& at2,
797 const std::string& at3,
798 const std::string& at4,
799 InversionType* inversionType) {
800 std::vector<std::string> keys;
801 keys.push_back(at1);
802 keys.push_back(at2);
803 keys.push_back(at3);
804 keys.push_back(at4);
805 return inversionTypeCont_.add(keys, inversionType);
806 }
807
808 bool ForceField::addNonBondedInteractionType(
809 const std::string& at1, const std::string& at2,
810 NonBondedInteractionType* nbiType) {
811 std::vector<std::string> keys;
812 keys.push_back(at1);
813 keys.push_back(at2);
814 return nonBondedInteractionTypeCont_.add(keys, nbiType);
815 }
816
817 RealType ForceField::getRcutFromAtomType(AtomType* at) {
818 RealType rcut(0.0);
819
820 LennardJonesAdapter lja = LennardJonesAdapter(at);
821 if (lja.isLennardJones()) { rcut = 2.5 * lja.getSigma(); }
822 EAMAdapter ea = EAMAdapter(at);
823 if (ea.isEAM()) {
824 switch (ea.getEAMType()) {
825 case eamFuncfl:
826 rcut = max(rcut, ea.getRcut());
827 break;
828 case eamZhou2001:
829 case eamZhou2004:
830 rcut = max(rcut, ea.getLatticeConstant() * sqrt(10.0) / 2.0);
831 break;
832 default:
833 break;
834 }
835 }
836 SuttonChenAdapter sca = SuttonChenAdapter(at);
837 if (sca.isSuttonChen()) { rcut = max(rcut, 2.0 * sca.getAlpha()); }
838 GayBerneAdapter gba = GayBerneAdapter(at);
839 if (gba.isGayBerne()) {
840 rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL()));
841 }
842 StickyAdapter sa = StickyAdapter(at);
843 if (sa.isSticky()) { rcut = max(rcut, max(sa.getRu(), sa.getRup())); }
844
845 return rcut;
846 }
847
848 ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
849 std::string forceFieldFilename(filename);
850
851 ifstrstream* ffStream = new ifstrstream();
852
853 // Try to open the force field file in current directory first
854 ffStream->open(forceFieldFilename.c_str());
855
856 if (!ffStream->is_open()) {
857 // If current directory does not contain the force field file,
858 // try to open it in ffPath_:
859 forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
860 ffStream->open(forceFieldFilename.c_str());
861
862 if (!ffStream->is_open()) {
863 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
864 "Error opening the force field parameter file:\n"
865 "\t%s\n"
866 "\tHave you tried setting the FORCE_PARAM_PATH environment "
867 "variable?\n",
868 forceFieldFilename.c_str());
869 painCave.severity = OPENMD_ERROR;
870 painCave.isFatal = 1;
871 simError();
872 }
873 }
874 return ffStream;
875 }
876} // namespace OpenMD
OpenMD::AtomType
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:69
OpenMD::AtomTypesSectionParser
"io/AtomTypesSectionParser.hpp"
Definition AtomTypesSectionParser.hpp:59
OpenMD::BaseAtomTypesSectionParser
"io/BaseAtomTypesSectionParser.hpp"
Definition BaseAtomTypesSectionParser.hpp:59
OpenMD::BendTypesSectionParser
"io/BendTypesSectionParser.hpp"
Definition BendTypesSectionParser.hpp:60
OpenMD::BondType
BondType class is responsible for calculating the force and energy of the bond.
Definition BondType.hpp:67
OpenMD::BondTypesSectionParser
"io/BondTypesSectionParser.hpp"
Definition BondTypesSectionParser.hpp:60
OpenMD::EAMAtomTypesSectionParser
"io/EAMAtomTypesSectionParser.hpp"
Definition EAMAtomTypesSectionParser.hpp:61
OpenMD::ForceField::ForceField
ForceField(std::string ffName)
Definition ForceField.cpp:90
OpenMD::ForceField::getAtomType
AtomType * getAtomType(const std::string &at)
getAtomType by string
Definition ForceField.cpp:186
OpenMD::InversionTypesSectionParser
"io/InversionTypesSectionParser.hpp"
Definition InversionTypesSectionParser.hpp:60
OpenMD::NonBondedInteractionType
NonBondedInteractionType class is responsible for keeping track of static (unchanging) parameters for...
Definition NonBondedInteractionType.hpp:85
OpenMD::OptionSectionParser
"io/OptionSectionParser.hpp"
Definition OptionSectionParser.hpp:61
OpenMD::SCAtomTypesSectionParser
"io/SCAtomTypesSectionParser.hpp"
Definition SCAtomTypesSectionParser.hpp:60
OpenMD::StickyAtomTypesSectionParser
"io/StickyAtomTypesSectionParser.hpp"
Definition StickyAtomTypesSectionParser.hpp:60
OpenMD::StickyPowerAtomTypesSectionParser
StickyPowerAtomTypesSectionParser.hpp "io/StickyAtomTypesSectionParser.hpp".
Definition StickyPowerAtomTypesSectionParser.hpp:60
OpenMD::TorsionTypesSectionParser
"io/TorsionTypesSectionParser.hpp"
Definition TorsionTypesSectionParser.hpp:60
OpenMD::TypeContainer::find
ElemPtr find(KeyType &keys)
Exact Match.
Definition TypeContainer.hpp:110
OpenMD::ifstrstream
ifstrstream class provides a stream interface to read data from files.
Definition ifstrstream.hpp:86
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:63