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"
67#include "utils/simError.h"
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());
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);
90 RigidBodyStamp* currentRigidBodyStamp;
91 size_t nRigidbodies = molStamp->getNRigidBodies();
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);
101 BondStamp* currentBondStamp;
102 size_t nBonds = molStamp->getNBonds();
104 for (
size_t i = 0; i < nBonds; ++i) {
105 currentBondStamp = molStamp->getBondStamp(i);
106 bond = createBond(ff, mol, currentBondStamp, localIndexMan);
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);
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);
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);
137 createInversion(ff, mol, currentInversionStamp, localIndexMan);
138 if (inversion != NULL) { mol->addInversion(inversion); }
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);
148 createCutoffGroup(mol, currentCutoffGroupStamp, localIndexMan);
149 mol->addCutoffGroup(cutoffGroup);
153 std::vector<Atom*> freeAtoms;
154 std::vector<Atom*>::iterator ai;
155 std::vector<Atom*>::iterator fai;
158 for (atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
159 freeAtoms.push_back(atom);
162 Molecule::CutoffGroupIterator ci;
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)) {
169 freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom),
177 for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
178 cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
179 mol->addCutoffGroup(cutoffGroup);
183 createConstraintPair(mol);
186 for (std::size_t i = 0; i < molStamp->getNConstraints(); ++i) {
187 ConstraintStamp* cStamp = molStamp->getConstraintStamp(i);
191 atomA = mol->getAtomAt(cStamp->getA());
192 atomB = mol->getAtomAt(cStamp->getB());
193 assert(atomA && atomB);
195 bool printConstraintForce =
false;
197 if (cStamp->havePrintConstraintForce()) {
198 printConstraintForce = cStamp->getPrintConstraintForce();
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;
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);
219 createConstraintElem(mol);
223 if (molStamp->haveConstrainTotalCharge()) {
224 mol->setConstrainTotalCharge(molStamp->getConstrainTotalCharge());
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());
244 painCave.isFatal = 1;
249 if (atomType->isDirectional()) {
254 atom =
new Atom(atomType);
257 atom->setLocalIndex(localIndexMan->getNextAtomIndex());
262 RigidBody* MoleculeCreator::createRigidBody(
271 nAtoms = rbStamp->getNMembers();
272 for (std::size_t i = 0; i < nAtoms; ++i) {
276 atom = mol->getAtomAt(rbStamp->getMemberAt(i));
277 atomStamp = molStamp->getAtomStamp(rbStamp->getMemberAt(i));
278 rb->addAtom(atom, atomStamp);
295 rb->setType(mol->
getType() +
"_RB_" + s.c_str());
299 Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
301 LocalIndexManager* localIndexMan) {
302 BondTypeParser btParser;
303 BondType* bondType = NULL;
307 atomA = mol->getAtomAt(stamp->getA());
308 atomB = mol->getAtomAt(stamp->getB());
310 assert(atomA && atomB);
312 if (stamp->hasOverride()) {
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 "
320 e.what(), mol->getType().c_str());
321 painCave.isFatal = 1;
326 bondType = ff->getBondType(atomA->getType(), atomB->getType());
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());
333 painCave.isFatal = 1;
338 Bond* bond =
new Bond(atomA, atomB, bondType);
341 bond->setLocalIndex(localIndexMan->getNextBondIndex());
349 std::string s = OpenMD_itoa(mol->getNBonds(), 10);
350 bond->setName(mol->getType() +
"_Bond_" + s.c_str());
354 Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
356 LocalIndexManager* localIndexMan) {
357 BendTypeParser btParser;
358 BendType* bendType = NULL;
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]);
367 assert(atomA && atomB && atomC);
369 if (stamp->hasOverride()) {
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 "
377 e.what(), mol->getType().c_str());
378 painCave.isFatal = 1;
383 ff->getBendType(atomA->getType().c_str(), atomB->getType().c_str(),
384 atomC->getType().c_str());
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());
392 painCave.isFatal = 1;
397 bend =
new Bend(atomA, atomB, atomC, bendType);
399 }
else if (bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
400 int ghostIndex = stamp->getGhostVectorSource();
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;
413 if (stamp->hasOverride()) {
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 "
421 e.what(), mol->getType().c_str());
422 painCave.isFatal = 1;
426 bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(),
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(),
435 painCave.isFatal = 1;
440 bend =
new GhostBend(normalAtom, ghostAtom, bendType);
444 bend->setLocalIndex(localIndexMan->getNextBendIndex());
452 std::string s = OpenMD_itoa(mol->getNBends(), 10);
453 bend->setName(mol->getType() +
"_Bend_" + s.c_str());
457 Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
459 LocalIndexManager* localIndexMan) {
460 TorsionTypeParser ttParser;
461 TorsionType* torsionType = NULL;
462 Torsion* torsion = NULL;
464 std::vector<int> torsionAtoms = stamp->getMembers();
465 if (torsionAtoms.size() < 3) {
return torsion; }
467 Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
468 Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
469 Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
471 if (torsionAtoms.size() == 4) {
472 Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
474 assert(atomA && atomB && atomC && atomD);
476 if (stamp->hasOverride()) {
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 "
484 e.what(), mol->getType().c_str());
485 painCave.isFatal = 1;
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());
497 painCave.isFatal = 1;
502 torsion =
new Torsion(atomA, atomB, atomC, atomD, torsionType);
504 DirectionalAtom* dAtom =
dynamic_cast<DirectionalAtom*
>(
505 mol->getAtomAt(stamp->getGhostVectorSource()));
507 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
508 "Can not cast Atom to DirectionalAtom");
509 painCave.isFatal = 1;
513 if (stamp->hasOverride()) {
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 "
521 e.what(), mol->getType().c_str());
522 painCave.isFatal = 1;
526 torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
527 atomC->getType(),
"GHOST");
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");
535 painCave.isFatal = 1;
540 torsion =
new GhostTorsion(atomA, atomB, dAtom, torsionType);
544 torsion->setLocalIndex(localIndexMan->getNextTorsionIndex());
553 std::string s = OpenMD_itoa(mol->getNTorsions(), 10);
554 torsion->setName(mol->getType() +
"_Torsion_" + s.c_str());
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;
565 int center = stamp->getCenter();
566 std::vector<int> satellites = stamp->getSatellites();
567 if (satellites.size() != 3) {
return inversion; }
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]);
574 assert(atomA && atomB && atomC && atomD);
576 if (stamp->hasOverride()) {
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 "
584 e.what(), mol->getType().c_str());
585 painCave.isFatal = 1;
589 inversionType = ff->getInversionType(atomA->getType(), atomB->getType(),
590 atomC->getType(), atomD->getType());
592 if (inversionType == NULL) {
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());
600 painCave.isFatal = 0;
601 painCave.severity = OPENMD_INFO;
605 if (inversionType != NULL) {
606 inversion =
new Inversion(atomA, atomB, atomC, atomD, inversionType);
610 inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
619 std::string s = OpenMD_itoa(mol->getNInversions(), 10);
620 inversion->setName(mol->getType() +
"_Inversion_" + s.c_str());
627 CutoffGroup* MoleculeCreator::createCutoffGroup(
628 Molecule* mol, CutoffGroupStamp* stamp,
629 LocalIndexManager* localIndexMan) {
633 cg =
new CutoffGroup();
635 nAtoms = stamp->getNMembers();
636 for (
size_t i = 0; i < nAtoms; ++i) {
637 atom = mol->getAtomAt(stamp->getMemberAt(i));
643 cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
648 CutoffGroup* MoleculeCreator::createCutoffGroup(
649 Molecule*, Atom* atom, LocalIndexManager* localIndexMan) {
651 cg =
new CutoffGroup();
655 cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
660 void MoleculeCreator::createConstraintPair(Molecule* mol) {
662 Molecule::BondIterator bi;
664 ConstraintPair* cPair;
666 for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
667 BondType* bt = bond->getBondType();
669 if (
typeid(FixedBondType) ==
typeid(*bt)) {
670 FixedBondType* fbt =
dynamic_cast<FixedBondType*
>(bt);
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);
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()) {
692 mol->addConstraintElem(
new ConstraintElem(sdA));
695 StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();
696 if (sdSet.find(sdB) == sdSet.end()) {
698 mol->addConstraintElem(
new ConstraintElem(sdB));
AtomType * getAtomType()
Returns the AtomType of this Atom.
AtomType is what OpenMD looks to for unchanging data about an atom.
"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.
std::string getType()
Returns the name of the molecule.
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.