OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
MoleculeCreator.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 MoleculeCreator.cpp
47 * @author tlin
48 * @date 11/04/2004
49 * @version 1.0
50 */
51
53
54#include <cassert>
55#include <set>
56#include <typeinfo>
57
59#include "primitives/GhostTorsion.hpp"
60#include "types/AtomType.hpp"
61#include "types/BendTypeParser.hpp"
62#include "types/BondTypeParser.hpp"
64#include "types/InversionTypeParser.hpp"
65#include "types/TorsionTypeParser.hpp"
66#include "utils/StringUtils.hpp"
67#include "utils/simError.h"
68
69namespace OpenMD {
70
71 Molecule* MoleculeCreator::createMolecule(ForceField* ff,
72 MoleculeStamp* molStamp,
73 int stampId, int globalIndex,
74 LocalIndexManager* localIndexMan) {
75 Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName(),
76 molStamp->getRegion());
77
78 // create atoms
79 Atom* atom;
80 AtomStamp* currentAtomStamp;
81 size_t nAtom = molStamp->getNAtoms();
82 for (size_t i = 0; i < nAtom; ++i) {
83 currentAtomStamp = molStamp->getAtomStamp(i);
84 atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
85 mol->addAtom(atom);
86 }
87
88 // create rigidbodies
89 RigidBody* rb;
90 RigidBodyStamp* currentRigidBodyStamp;
91 size_t nRigidbodies = molStamp->getNRigidBodies();
92
93 for (size_t i = 0; i < nRigidbodies; ++i) {
94 currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
95 rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
96 mol->addRigidBody(rb);
97 }
98
99 // create bonds
100 Bond* bond;
101 BondStamp* currentBondStamp;
102 size_t nBonds = molStamp->getNBonds();
103
104 for (size_t i = 0; i < nBonds; ++i) {
105 currentBondStamp = molStamp->getBondStamp(i);
106 bond = createBond(ff, mol, currentBondStamp, localIndexMan);
107 mol->addBond(bond);
108 }
109
110 // create bends
111 Bend* bend;
112 BendStamp* currentBendStamp;
113 size_t nBends = molStamp->getNBends();
114 for (size_t i = 0; i < nBends; ++i) {
115 currentBendStamp = molStamp->getBendStamp(i);
116 bend = createBend(ff, mol, currentBendStamp, localIndexMan);
117 mol->addBend(bend);
118 }
119
120 // create torsions
121 Torsion* torsion;
122 TorsionStamp* currentTorsionStamp;
123 size_t nTorsions = molStamp->getNTorsions();
124 for (size_t i = 0; i < nTorsions; ++i) {
125 currentTorsionStamp = molStamp->getTorsionStamp(i);
126 torsion = createTorsion(ff, mol, currentTorsionStamp, localIndexMan);
127 mol->addTorsion(torsion);
128 }
129
130 // create inversions
131 Inversion* inversion;
132 InversionStamp* currentInversionStamp;
133 size_t nInversions = molStamp->getNInversions();
134 for (size_t i = 0; i < nInversions; ++i) {
135 currentInversionStamp = molStamp->getInversionStamp(i);
136 inversion =
137 createInversion(ff, mol, currentInversionStamp, localIndexMan);
138 if (inversion != NULL) { mol->addInversion(inversion); }
139 }
140
141 // create cutoffGroups
142 CutoffGroup* cutoffGroup;
143 CutoffGroupStamp* currentCutoffGroupStamp;
144 size_t nCutoffGroups = molStamp->getNCutoffGroups();
145 for (size_t i = 0; i < nCutoffGroups; ++i) {
146 currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
147 cutoffGroup =
148 createCutoffGroup(mol, currentCutoffGroupStamp, localIndexMan);
149 mol->addCutoffGroup(cutoffGroup);
150 }
151
152 // every free atom is a cutoff group
153 std::vector<Atom*> freeAtoms;
154 std::vector<Atom*>::iterator ai;
155 std::vector<Atom*>::iterator fai;
156
157 // add all atoms into allAtoms set
158 for (atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
159 freeAtoms.push_back(atom);
160 }
161
162 Molecule::CutoffGroupIterator ci;
163 CutoffGroup* cg;
164
165 for (cg = mol->beginCutoffGroup(ci); cg != NULL;
166 cg = mol->nextCutoffGroup(ci)) {
167 for (atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
168 // erase the atoms belong to cutoff groups from freeAtoms vector
169 freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom),
170 freeAtoms.end());
171 }
172 }
173
174 // loop over the free atoms and then create one cutoff group for
175 // every single free atom
176
177 for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
178 cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
179 mol->addCutoffGroup(cutoffGroup);
180 }
181
182 // create bonded constraintPairs:
183 createConstraintPair(mol);
184
185 // create non-bonded constraintPairs
186 for (std::size_t i = 0; i < molStamp->getNConstraints(); ++i) {
187 ConstraintStamp* cStamp = molStamp->getConstraintStamp(i);
188 Atom* atomA;
189 Atom* atomB;
190
191 atomA = mol->getAtomAt(cStamp->getA());
192 atomB = mol->getAtomAt(cStamp->getB());
193 assert(atomA && atomB);
194
195 bool printConstraintForce = false;
196
197 if (cStamp->havePrintConstraintForce()) {
198 printConstraintForce = cStamp->getPrintConstraintForce();
199 }
200
201 if (!cStamp->haveConstrainedDistance()) {
202 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
203 "Constraint Error: A non-bond constraint was specified\n"
204 "\twithout providing a value for the constrainedDistance.\n");
205 painCave.isFatal = 1;
206 simError();
207 } else {
208 RealType distance = cStamp->getConstrainedDistance();
209 ConstraintElem* consElemA = new ConstraintElem(atomA);
210 ConstraintElem* consElemB = new ConstraintElem(atomB);
211 ConstraintPair* cPair = new ConstraintPair(
212 consElemA, consElemB, distance, printConstraintForce);
213 mol->addConstraintPair(cPair);
214 }
215 }
216
217 // now create the constraint elements:
218
219 createConstraintElem(mol);
220
221 // Does this molecule stamp define a total constrained charge value?
222 // If so, let the created molecule know about it.
223 if (molStamp->haveConstrainTotalCharge()) {
224 mol->setConstrainTotalCharge(molStamp->getConstrainTotalCharge());
225 }
226
227 // The construction of this molecule is finished:
228 mol->complete();
229
230 return mol;
231 }
232
234 LocalIndexManager* localIndexMan) {
235 AtomType* atomType;
236 Atom* atom;
237
238 atomType = ff->getAtomType(stamp->getType());
239 if (atomType == NULL) {
240 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
241 "Can not find Matching Atom Type for[%s]",
242 stamp->getType().c_str());
243
244 painCave.isFatal = 1;
245 simError();
246 }
247
248 // below code still have some kind of hard-coding smell
249 if (atomType->isDirectional()) {
250 DirectionalAtom* dAtom;
251 dAtom = new DirectionalAtom(atomType);
252 atom = dAtom;
253 } else {
254 atom = new Atom(atomType);
255 }
256
257 atom->setLocalIndex(localIndexMan->getNextAtomIndex());
258
259 return atom;
260 }
261
262 RigidBody* MoleculeCreator::createRigidBody(
263 MoleculeStamp* molStamp, Molecule* mol, RigidBodyStamp* rbStamp,
264 LocalIndexManager* localIndexMan) {
265 Atom* atom;
266 size_t nAtoms;
267 Vector3d refCoor;
268 AtomStamp* atomStamp;
269
270 RigidBody* rb = new RigidBody();
271 nAtoms = rbStamp->getNMembers();
272 for (std::size_t i = 0; i < nAtoms; ++i) {
273 // rbStamp->getMember(i) return the local index of current atom
274 // inside the molecule. It is not the same as local index of
275 // atom which is the index of atom at DataStorage class
276 atom = mol->getAtomAt(rbStamp->getMemberAt(i));
277 atomStamp = molStamp->getAtomStamp(rbStamp->getMemberAt(i));
278 rb->addAtom(atom, atomStamp);
279 }
280
281 // after all of the atoms are added, we need to calculate the
282 // reference coordinates
283 rb->calcRefCoords();
284
285 // set the local index of this rigid body, global index will be set later
286 rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
287
288 // The rule for naming a rigidbody is: MoleculeName_RB_Integer
289 // The first part is the name of the molecule
290 // The second part is always fixed as "RB"
291 // The third part is the index of the rigidbody defined in meta-data file
292 // For example, Butane_RB_0 is a valid rigid body name of butane molecule
293
294 std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
295 rb->setType(mol->getType() + "_RB_" + s.c_str());
296 return rb;
297 }
298
299 Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
300 BondStamp* stamp,
301 LocalIndexManager* localIndexMan) {
302 BondTypeParser btParser;
303 BondType* bondType = NULL;
304 Atom* atomA;
305 Atom* atomB;
306
307 atomA = mol->getAtomAt(stamp->getA());
308 atomB = mol->getAtomAt(stamp->getB());
309
310 assert(atomA && atomB);
311
312 if (stamp->hasOverride()) {
313 try {
314 bondType = btParser.parseTypeAndPars(stamp->getOverrideType(),
315 stamp->getOverridePars());
316 } catch (OpenMDException& e) {
317 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
318 "MoleculeCreator Error: %s "
319 "for molecule %s\n",
320 e.what(), mol->getType().c_str());
321 painCave.isFatal = 1;
322 simError();
323 }
324
325 } else {
326 bondType = ff->getBondType(atomA->getType(), atomB->getType());
327
328 if (bondType == NULL) {
329 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
330 "Can not find Matching Bond Type for[%s, %s]",
331 atomA->getType().c_str(), atomB->getType().c_str());
332
333 painCave.isFatal = 1;
334 simError();
335 }
336 }
337
338 Bond* bond = new Bond(atomA, atomB, bondType);
339
340 // set the local index of this bond, the global index will be set later
341 bond->setLocalIndex(localIndexMan->getNextBondIndex());
342
343 // The rule for naming a bond is: MoleculeName_Bond_Integer
344 // The first part is the name of the molecule
345 // The second part is always fixed as "Bond"
346 // The third part is the index of the bond defined in meta-data file
347 // For example, Butane_bond_0 is a valid Bond name in a butane molecule
348
349 std::string s = OpenMD_itoa(mol->getNBonds(), 10);
350 bond->setName(mol->getType() + "_Bond_" + s.c_str());
351 return bond;
352 }
353
354 Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
355 BendStamp* stamp,
356 LocalIndexManager* localIndexMan) {
357 BendTypeParser btParser;
358 BendType* bendType = NULL;
359 Bend* bend = NULL;
360
361 std::vector<int> bendAtoms = stamp->getMembers();
362 if (bendAtoms.size() == 3) {
363 Atom* atomA = mol->getAtomAt(bendAtoms[0]);
364 Atom* atomB = mol->getAtomAt(bendAtoms[1]);
365 Atom* atomC = mol->getAtomAt(bendAtoms[2]);
366
367 assert(atomA && atomB && atomC);
368
369 if (stamp->hasOverride()) {
370 try {
371 bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
372 stamp->getOverridePars());
373 } catch (OpenMDException& e) {
374 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
375 "MoleculeCreator Error: %s "
376 "for molecule %s\n",
377 e.what(), mol->getType().c_str());
378 painCave.isFatal = 1;
379 simError();
380 }
381 } else {
382 bendType =
383 ff->getBendType(atomA->getType().c_str(), atomB->getType().c_str(),
384 atomC->getType().c_str());
385
386 if (bendType == NULL) {
387 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
388 "Can not find Matching Bend Type for[%s, %s, %s]",
389 atomA->getType().c_str(), atomB->getType().c_str(),
390 atomC->getType().c_str());
391
392 painCave.isFatal = 1;
393 simError();
394 }
395 }
396
397 bend = new Bend(atomA, atomB, atomC, bendType);
398
399 } else if (bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
400 int ghostIndex = stamp->getGhostVectorSource();
401 int normalIndex =
402 ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
403 Atom* normalAtom = mol->getAtomAt(normalIndex);
404 DirectionalAtom* ghostAtom =
405 dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
406 if (ghostAtom == NULL) {
407 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
408 "Can not cast Atom to DirectionalAtom");
409 painCave.isFatal = 1;
410 simError();
411 }
412
413 if (stamp->hasOverride()) {
414 try {
415 bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
416 stamp->getOverridePars());
417 } catch (OpenMDException& e) {
418 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
419 "MoleculeCreator Error: %s "
420 "for molecule %s\n",
421 e.what(), mol->getType().c_str());
422 painCave.isFatal = 1;
423 simError();
424 }
425 } else {
426 bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(),
427 "GHOST");
428
429 if (bendType == NULL) {
430 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
431 "Can not find Matching Bend Type for[%s, %s, %s]",
432 normalAtom->getType().c_str(), ghostAtom->getType().c_str(),
433 "GHOST");
434
435 painCave.isFatal = 1;
436 simError();
437 }
438 }
439
440 bend = new GhostBend(normalAtom, ghostAtom, bendType);
441 }
442
443 // set the local index of this bend, the global index will be set later
444 bend->setLocalIndex(localIndexMan->getNextBendIndex());
445
446 // The rule for naming a bend is: MoleculeName_Bend_Integer
447 // The first part is the name of the molecule
448 // The second part is always fixed as "Bend"
449 // The third part is the index of the bend defined in meta-data file
450 // For example, Butane_Bend_0 is a valid Bend name in a butane molecule
451
452 std::string s = OpenMD_itoa(mol->getNBends(), 10);
453 bend->setName(mol->getType() + "_Bend_" + s.c_str());
454 return bend;
455 }
456
457 Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
458 TorsionStamp* stamp,
459 LocalIndexManager* localIndexMan) {
460 TorsionTypeParser ttParser;
461 TorsionType* torsionType = NULL;
462 Torsion* torsion = NULL;
463
464 std::vector<int> torsionAtoms = stamp->getMembers();
465 if (torsionAtoms.size() < 3) { return torsion; }
466
467 Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
468 Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
469 Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
470
471 if (torsionAtoms.size() == 4) {
472 Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
473
474 assert(atomA && atomB && atomC && atomD);
475
476 if (stamp->hasOverride()) {
477 try {
478 torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
479 stamp->getOverridePars());
480 } catch (OpenMDException& e) {
481 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
482 "MoleculeCreator Error: %s "
483 "for molecule %s\n",
484 e.what(), mol->getType().c_str());
485 painCave.isFatal = 1;
486 simError();
487 }
488 } else {
489 torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
490 atomC->getType(), atomD->getType());
491 if (torsionType == NULL) {
492 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
493 "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
494 atomA->getType().c_str(), atomB->getType().c_str(),
495 atomC->getType().c_str(), atomD->getType().c_str());
496
497 painCave.isFatal = 1;
498 simError();
499 }
500 }
501
502 torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);
503 } else {
504 DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(
505 mol->getAtomAt(stamp->getGhostVectorSource()));
506 if (dAtom == NULL) {
507 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
508 "Can not cast Atom to DirectionalAtom");
509 painCave.isFatal = 1;
510 simError();
511 }
512
513 if (stamp->hasOverride()) {
514 try {
515 torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
516 stamp->getOverridePars());
517 } catch (OpenMDException& e) {
518 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
519 "MoleculeCreator Error: %s "
520 "for molecule %s\n",
521 e.what(), mol->getType().c_str());
522 painCave.isFatal = 1;
523 simError();
524 }
525 } else {
526 torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
527 atomC->getType(), "GHOST");
528
529 if (torsionType == NULL) {
530 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
531 "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
532 atomA->getType().c_str(), atomB->getType().c_str(),
533 atomC->getType().c_str(), "GHOST");
534
535 painCave.isFatal = 1;
536 simError();
537 }
538 }
539
540 torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);
541 }
542
543 // set the local index of this torsion, the global index will be set later
544 torsion->setLocalIndex(localIndexMan->getNextTorsionIndex());
545
546 // The rule for naming a torsion is: MoleculeName_Torsion_Integer
547 // The first part is the name of the molecule
548 // The second part is always fixed as "Torsion"
549 // The third part is the index of the torsion defined in meta-data file
550 // For example, Butane_Torsion_0 is a valid Torsion name in a
551 // butane molecule
552
553 std::string s = OpenMD_itoa(mol->getNTorsions(), 10);
554 torsion->setName(mol->getType() + "_Torsion_" + s.c_str());
555 return torsion;
556 }
557
558 Inversion* MoleculeCreator::createInversion(
559 ForceField* ff, Molecule* mol, InversionStamp* stamp,
560 LocalIndexManager* localIndexMan) {
561 InversionTypeParser itParser;
562 InversionType* inversionType = NULL;
563 Inversion* inversion = NULL;
564
565 int center = stamp->getCenter();
566 std::vector<int> satellites = stamp->getSatellites();
567 if (satellites.size() != 3) { return inversion; }
568
569 Atom* atomA = mol->getAtomAt(center);
570 Atom* atomB = mol->getAtomAt(satellites[0]);
571 Atom* atomC = mol->getAtomAt(satellites[1]);
572 Atom* atomD = mol->getAtomAt(satellites[2]);
573
574 assert(atomA && atomB && atomC && atomD);
575
576 if (stamp->hasOverride()) {
577 try {
578 inversionType = itParser.parseTypeAndPars(stamp->getOverrideType(),
579 stamp->getOverridePars());
580 } catch (OpenMDException& e) {
581 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
582 "MoleculeCreator Error: %s "
583 "for molecule %s\n",
584 e.what(), mol->getType().c_str());
585 painCave.isFatal = 1;
586 simError();
587 }
588 } else {
589 inversionType = ff->getInversionType(atomA->getType(), atomB->getType(),
590 atomC->getType(), atomD->getType());
591
592 if (inversionType == NULL) {
593 snprintf(
594 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
595 "No Matching Inversion Type for[%s, %s, %s, %s]\n"
596 "\t(May not be a problem: not all inversions are parametrized)\n",
597 atomA->getType().c_str(), atomB->getType().c_str(),
598 atomC->getType().c_str(), atomD->getType().c_str());
599
600 painCave.isFatal = 0;
601 painCave.severity = OPENMD_INFO;
602 simError();
603 }
604 }
605 if (inversionType != NULL) {
606 inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
607
608 // set the local index of this inversion, the global index will
609 // be set later
610 inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
611
612 // The rule for naming an inversion is: MoleculeName_Inversion_Integer
613 // The first part is the name of the molecule
614 // The second part is always fixed as "Inversion"
615 // The third part is the index of the inversion defined in meta-data file
616 // For example, Benzene_Inversion_0 is a valid Inversion name in a
617 // Benzene molecule
618
619 std::string s = OpenMD_itoa(mol->getNInversions(), 10);
620 inversion->setName(mol->getType() + "_Inversion_" + s.c_str());
621 return inversion;
622 } else {
623 return NULL;
624 }
625 }
626
627 CutoffGroup* MoleculeCreator::createCutoffGroup(
628 Molecule* mol, CutoffGroupStamp* stamp,
629 LocalIndexManager* localIndexMan) {
630 size_t nAtoms;
631 CutoffGroup* cg;
632 Atom* atom;
633 cg = new CutoffGroup();
634
635 nAtoms = stamp->getNMembers();
636 for (size_t i = 0; i < nAtoms; ++i) {
637 atom = mol->getAtomAt(stamp->getMemberAt(i));
638 assert(atom);
639 cg->addAtom(atom);
640 }
641
642 // set the local index of this cutoffGroup, global index will be set later
643 cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
644
645 return cg;
646 }
647
648 CutoffGroup* MoleculeCreator::createCutoffGroup(
649 Molecule*, Atom* atom, LocalIndexManager* localIndexMan) {
650 CutoffGroup* cg;
651 cg = new CutoffGroup();
652 cg->addAtom(atom);
653
654 // set the local index of this cutoffGroup, global index will be set later
655 cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
656
657 return cg;
658 }
659
660 void MoleculeCreator::createConstraintPair(Molecule* mol) {
661 // add bond constraints
662 Molecule::BondIterator bi;
663 Bond* bond;
664 ConstraintPair* cPair;
665
666 for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
667 BondType* bt = bond->getBondType();
668
669 if (typeid(FixedBondType) == typeid(*bt)) {
670 FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
671
672 ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
673 ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());
674 cPair = new ConstraintPair(consElemA, consElemB,
675 fbt->getEquilibriumBondLength(), false);
676 mol->addConstraintPair(cPair);
677 }
678 }
679
680 // rigidbody -- rigidbody constraint is not support yet
681 }
682
683 void MoleculeCreator::createConstraintElem(Molecule* mol) {
684 ConstraintPair* consPair;
685 Molecule::ConstraintPairIterator cpi;
686 std::set<StuntDouble*> sdSet;
687 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
688 consPair = mol->nextConstraintPair(cpi)) {
689 StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();
690 if (sdSet.find(sdA) == sdSet.end()) {
691 sdSet.insert(sdA);
692 mol->addConstraintElem(new ConstraintElem(sdA));
693 }
694
695 StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();
696 if (sdSet.find(sdB) == sdSet.end()) {
697 sdSet.insert(sdB);
698 mol->addConstraintElem(new ConstraintElem(sdB));
699 }
700 }
701 }
702} // namespace OpenMD
AtomType * getAtomType()
Returns the AtomType of this Atom.
Definition Atom.hpp:83
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
"utils/LocalIndexManager.hpp"
virtual Atom * createAtom(ForceField *ff, Molecule *, AtomStamp *stamp, LocalIndexManager *localIndexMan)
Create an atom by its stamp.
size_t getNRigidBodies()
Returns the total number of rigid bodies in this molecule.
Definition Molecule.hpp:176
std::string getType()
Returns the name of the molecule.
Definition Molecule.hpp:120
void setLocalIndex(int index)
Sets the local index of this stuntDouble.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Real distance(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the distance between two DynamicVectors.