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