--- trunk/src/brains/MoleculeCreator.cpp 2005/03/08 21:06:49 403 +++ trunk/src/brains/MoleculeCreator.cpp 2014/12/02 22:11:04 2046 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,42 +28,58 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ - /** - * @file MoleculeCreator.cpp - * @author tlin - * @date 11/04/2004 - * @time 13:44am - * @version 1.0 - */ +/** + * @file MoleculeCreator.cpp + * @author tlin + * @date 11/04/2004 + * @version 1.0 + */ #include +#include #include #include "brains/MoleculeCreator.hpp" #include "primitives/GhostBend.hpp" #include "primitives/GhostTorsion.hpp" -#include "types/DirectionalAtomType.hpp" +#include "types/AtomType.hpp" #include "types/FixedBondType.hpp" +#include "types/BondTypeParser.hpp" +#include "types/BendTypeParser.hpp" +#include "types/TorsionTypeParser.hpp" +#include "types/InversionTypeParser.hpp" #include "utils/simError.h" #include "utils/StringUtils.hpp" -namespace oopse { +namespace OpenMD { + + Molecule* MoleculeCreator::createMolecule(ForceField* ff, + MoleculeStamp *molStamp, + int stampId, int globalIndex, + LocalIndexManager* localIndexMan) { + Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName(), + molStamp->getRegion() ); -Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp, - int stampId, int globalIndex, LocalIndexManager* localIndexMan) { - - Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID()); - //create atoms Atom* atom; AtomStamp* currentAtomStamp; int nAtom = molStamp->getNAtoms(); for (int i = 0; i < nAtom; ++i) { - currentAtomStamp = molStamp->getAtom(i); - atom = createAtom(ff, mol, currentAtomStamp, localIndexMan); - mol->addAtom(atom); + currentAtomStamp = molStamp->getAtomStamp(i); + atom = createAtom(ff, mol, currentAtomStamp, localIndexMan); + mol->addAtom(atom); } //create rigidbodies @@ -81,20 +88,21 @@ Molecule* MoleculeCreator::createMolecule(ForceField* int nRigidbodies = molStamp->getNRigidBodies(); for (int i = 0; i < nRigidbodies; ++i) { - currentRigidBodyStamp = molStamp->getRigidBody(i); - rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan); - mol->addRigidBody(rb); + currentRigidBodyStamp = molStamp->getRigidBodyStamp(i); + rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, + localIndexMan); + mol->addRigidBody(rb); } - + //create bonds Bond* bond; BondStamp* currentBondStamp; int nBonds = molStamp->getNBonds(); for (int i = 0; i < nBonds; ++i) { - currentBondStamp = molStamp->getBond(i); - bond = createBond(ff, mol, currentBondStamp); - mol->addBond(bond); + currentBondStamp = molStamp->getBondStamp(i); + bond = createBond(ff, mol, currentBondStamp, localIndexMan); + mol->addBond(bond); } //create bends @@ -102,9 +110,9 @@ Molecule* MoleculeCreator::createMolecule(ForceField* BendStamp* currentBendStamp; int nBends = molStamp->getNBends(); for (int i = 0; i < nBends; ++i) { - currentBendStamp = molStamp->getBend(i); - bend = createBend(ff, mol, currentBendStamp); - mol->addBend(bend); + currentBendStamp = molStamp->getBendStamp(i); + bend = createBend(ff, mol, currentBendStamp, localIndexMan); + mol->addBend(bend); } //create torsions @@ -112,120 +120,157 @@ Molecule* MoleculeCreator::createMolecule(ForceField* TorsionStamp* currentTorsionStamp; int nTorsions = molStamp->getNTorsions(); for (int i = 0; i < nTorsions; ++i) { - currentTorsionStamp = molStamp->getTorsion(i); - torsion = createTorsion(ff, mol, currentTorsionStamp); - mol->addTorsion(torsion); + currentTorsionStamp = molStamp->getTorsionStamp(i); + torsion = createTorsion(ff, mol, currentTorsionStamp, localIndexMan); + mol->addTorsion(torsion); } + //create inversions + Inversion* inversion; + InversionStamp* currentInversionStamp; + int nInversions = molStamp->getNInversions(); + for (int i = 0; i < nInversions; ++i) { + currentInversionStamp = molStamp->getInversionStamp(i); + inversion = createInversion(ff, mol, currentInversionStamp, + localIndexMan); + if (inversion != NULL ) { + mol->addInversion(inversion); + } + } + //create cutoffGroups CutoffGroup* cutoffGroup; CutoffGroupStamp* currentCutoffGroupStamp; int nCutoffGroups = molStamp->getNCutoffGroups(); for (int i = 0; i < nCutoffGroups; ++i) { - currentCutoffGroupStamp = molStamp->getCutoffGroup(i); - cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp); - mol->addCutoffGroup(cutoffGroup); + currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i); + cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp, + localIndexMan); + mol->addCutoffGroup(cutoffGroup); } //every free atom is a cutoff group - std::set allAtoms; - Molecule::AtomIterator ai; + std::vector freeAtoms; + std::vector::iterator ai; + std::vector::iterator fai; //add all atoms into allAtoms set - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - allAtoms.insert(atom); + for(atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) { + freeAtoms.push_back(atom); } Molecule::CutoffGroupIterator ci; CutoffGroup* cg; - std::set cutoffAtoms; - //add all of the atoms belong to cutoff groups into cutoffAtoms set - for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { - - for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { - cutoffAtoms.insert(atom); - } - + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { + + for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { + //erase the atoms belong to cutoff groups from freeAtoms vector + freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom), + freeAtoms.end()); + } } - //find all free atoms (which do not belong to cutoff groups) - //performs the "difference" operation from set theory, the output range contains a copy of every - //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in - //[cutoffAtoms.begin(), cutoffAtoms.end()). - std::vector freeAtoms; - std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(), - std::back_inserter(freeAtoms)); + // loop over the free atoms and then create one cutoff group for + // every single free atom + + for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) { + cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan); + mol->addCutoffGroup(cutoffGroup); + } - if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) { - //Some atoms in rigidAtoms are not in allAtoms, something must be wrong - sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule"); + //create bonded constraintPairs: + createConstraintPair(mol); + //create non-bonded constraintPairs + for (int i = 0; i < molStamp->getNConstraints(); ++i) { + ConstraintStamp* cStamp = molStamp->getConstraintStamp(i); + Atom* atomA; + Atom* atomB; + + atomA = mol->getAtomAt(cStamp->getA()); + atomB = mol->getAtomAt(cStamp->getB()); + assert( atomA && atomB ); + + RealType distance; + bool printConstraintForce = false; + + if (!cStamp->haveConstrainedDistance()) { + sprintf(painCave.errMsg, + "Constraint Error: A non-bond constraint was specified\n" + "\twithout providing a value for the constrainedDistance.\n"); painCave.isFatal = 1; - simError(); + simError(); + } else { + distance = cStamp->getConstrainedDistance(); + } + + if (cStamp->havePrintConstraintForce()) { + printConstraintForce = cStamp->getPrintConstraintForce(); + } + + ConstraintElem* consElemA = new ConstraintElem(atomA); + ConstraintElem* consElemB = new ConstraintElem(atomB); + ConstraintPair* cPair = new ConstraintPair(consElemA, consElemB, distance, + printConstraintForce); + mol->addConstraintPair(cPair); } - //loop over the free atoms and then create one cutoff group for every single free atom - std::vector::iterator fai; + // now create the constraint elements: - for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) { - cutoffGroup = createCutoffGroup(mol, *fai); - mol->addCutoffGroup(cutoffGroup); - } - //create constraints - createConstraintPair(mol); createConstraintElem(mol); + // Does this molecule stamp define a total constrained charge value? + // If so, let the created molecule know about it. + + if (molStamp->haveConstrainTotalCharge() ) { + mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() ); + } + //the construction of this molecule is finished mol->complete(); - + return mol; -} + } -Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp, - LocalIndexManager* localIndexMan) { + Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, + AtomStamp* stamp, + LocalIndexManager* localIndexMan) { AtomType * atomType; Atom* atom; atomType = ff->getAtomType(stamp->getType()); - + if (atomType == NULL) { - sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]", - stamp->getType()); + sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]", + stamp->getType().c_str()); - painCave.isFatal = 1; - simError(); + painCave.isFatal = 1; + simError(); } - + //below code still have some kind of hard-coding smell if (atomType->isDirectional()){ - - DirectionalAtomType* dAtomType = dynamic_cast(atomType); - - if (dAtomType == NULL) { - sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType"); - painCave.isFatal = 1; - simError(); - } - - DirectionalAtom* dAtom; - dAtom = new DirectionalAtom(dAtomType); - atom = dAtom; + DirectionalAtom* dAtom; + dAtom = new DirectionalAtom(atomType); + atom = dAtom; } else{ - atom = new Atom(atomType); + atom = new Atom(atomType); } atom->setLocalIndex(localIndexMan->getNextAtomIndex()); return atom; -} - -RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol, - RigidBodyStamp* rbStamp, - LocalIndexManager* localIndexMan) { + } + + RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, + Molecule* mol, + RigidBodyStamp* rbStamp, + LocalIndexManager* localIndexMan){ Atom* atom; int nAtoms; Vector3d refCoor; @@ -234,207 +279,391 @@ RigidBody* MoleculeCreator::createRigidBody(MoleculeSt RigidBody* rb = new RigidBody(); nAtoms = rbStamp->getNMembers(); for (int i = 0; i < nAtoms; ++i) { - //rbStamp->getMember(i) return the local index of current atom inside the molecule. - //It is not the same as local index of atom which is the index of atom at DataStorage class - atom = mol->getAtomAt(rbStamp->getMember(i)); - atomStamp= molStamp->getAtom(rbStamp->getMember(i)); - rb->addAtom(atom, atomStamp); + //rbStamp->getMember(i) return the local index of current atom + //inside the molecule. It is not the same as local index of + //atom which is the index of atom at DataStorage class + atom = mol->getAtomAt(rbStamp->getMemberAt(i)); + atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i)); + rb->addAtom(atom, atomStamp); } - //after all of the atoms are added, we need to calculate the reference coordinates + //after all of the atoms are added, we need to calculate the + //reference coordinates rb->calcRefCoords(); //set the local index of this rigid body, global index will be set later rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex()); - //the rule for naming rigidbody MoleculeName_RB_Integer - //The first part is the name of the molecule - //The second part is alway fixed as "RB" - //The third part is the index of the rigidbody defined in meta-data file - //For example, Butane_RB_0 is a valid rigid body name of butane molecule - /**@todo replace itoa by lexi_cast */ - std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10); - rb->setType(mol->getType() + "_RB_" + s.c_str()); + // The rule for naming a rigidbody is: MoleculeName_RB_Integer + // The first part is the name of the molecule + // The second part is always fixed as "RB" + // The third part is the index of the rigidbody defined in meta-data file + // For example, Butane_RB_0 is a valid rigid body name of butane molecule + std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10); + rb->setType(mol->getType() + "_RB_" + s.c_str()); return rb; -} + } -Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) { - BondType* bondType; + Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, + BondStamp* stamp, + LocalIndexManager* localIndexMan) { + BondTypeParser btParser; + BondType* bondType = NULL; Atom* atomA; Atom* atomB; - + atomA = mol->getAtomAt(stamp->getA()); atomB = mol->getAtomAt(stamp->getB()); - - assert( atomA && atomB); - bondType = ff->getBondType(atomA->getType(), atomB->getType()); + assert( atomA && atomB); - if (bondType == NULL) { - sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]", - atomA->getType().c_str(), - atomB->getType().c_str()); + if (stamp->hasOverride()) { + try { + bondType = btParser.parseTypeAndPars(stamp->getOverrideType(), + stamp->getOverridePars() ); + } + catch( OpenMDException e) { + sprintf(painCave.errMsg, "MoleculeCreator Error: %s " + "for molecule %s\n", + e.what(), mol->getType().c_str() ); painCave.isFatal = 1; simError(); - } - return new Bond(atomA, atomB, bondType); -} + } -Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) { - bool isGhostBend = false; - int ghostIndex; + } else { + bondType = ff->getBondType(atomA->getType(), atomB->getType()); - - // - if (stamp->haveExtras()){ - LinkedAssign* extras = stamp->getExtras(); - LinkedAssign* currentExtra = extras; - - while (currentExtra != NULL){ - if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){ - switch (currentExtra->getType()){ - case 0: - ghostIndex = currentExtra->getInt(); - isGhostBend = true; - break; - - default: - sprintf(painCave.errMsg, - "SimSetup Error: ghostVectorSource must be an int.\n"); - painCave.isFatal = 1; - simError(); - } - } else{ - sprintf(painCave.errMsg, - "SimSetup Error: unhandled bend assignment:\n"); - painCave.isFatal = 1; - simError(); - } - currentExtra = currentExtra->getNext(); - } + if (bondType == NULL) { + sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]", + atomA->getType().c_str(), + atomB->getType().c_str()); + painCave.isFatal = 1; + simError(); + } } + + Bond* bond = new Bond(atomA, atomB, bondType); - if (isGhostBend) { + //set the local index of this bond, the global index will be set later + bond->setLocalIndex(localIndexMan->getNextBondIndex()); - int indexA = stamp->getA(); - int indexB= stamp->getB(); + // The rule for naming a bond is: MoleculeName_Bond_Integer + // The first part is the name of the molecule + // The second part is always fixed as "Bond" + // The third part is the index of the bond defined in meta-data file + // For example, Butane_bond_0 is a valid Bond name in a butane molecule - assert(indexA != indexB); + std::string s = OpenMD_itoa(mol->getNBonds(), 10); + bond->setName(mol->getType() + "_Bond_" + s.c_str()); + return bond; + } + + Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, + BendStamp* stamp, + LocalIndexManager* localIndexMan) { + BendTypeParser btParser; + BendType* bendType = NULL; + Bend* bend = NULL; + + std::vector bendAtoms = stamp->getMembers(); + if (bendAtoms.size() == 3) { + Atom* atomA = mol->getAtomAt(bendAtoms[0]); + Atom* atomB = mol->getAtomAt(bendAtoms[1]); + Atom* atomC = mol->getAtomAt(bendAtoms[2]); + + assert( atomA && atomB && atomC ); - int normalIndex; - if (indexA == ghostIndex) { - normalIndex = indexB; - } else if (indexB == ghostIndex) { - normalIndex = indexA; - } + if (stamp->hasOverride()) { - Atom* normalAtom = mol->getAtomAt(normalIndex) ; - DirectionalAtom* ghostAtom = dynamic_cast(mol->getAtomAt(ghostIndex)); - if (ghostAtom == NULL) { - sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom"); - painCave.isFatal = 1; - simError(); + try { + bendType = btParser.parseTypeAndPars(stamp->getOverrideType(), + stamp->getOverridePars() ); } - - BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST"); - - if (bendType == NULL) { - sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]", - normalAtom->getType().c_str(), - ghostAtom->getType().c_str(), - "GHOST"); - - painCave.isFatal = 1; - simError(); + catch( OpenMDException e) { + sprintf(painCave.errMsg, "MoleculeCreator Error: %s " + "for molecule %s\n", + e.what(), mol->getType().c_str() ); + painCave.isFatal = 1; + simError(); } + } else { - return new GhostBend(normalAtom, ghostAtom, bendType); - - } else { - - Atom* atomA = mol->getAtomAt(stamp->getA()); - Atom* atomB = mol->getAtomAt(stamp->getB()); - Atom* atomC = mol->getAtomAt(stamp->getC()); - - assert( atomA && atomB && atomC); - - BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType()); - - if (bendType == NULL) { - sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]", - atomA->getType().c_str(), - atomB->getType().c_str(), - atomC->getType().c_str()); + bendType = ff->getBendType(atomA->getType().c_str(), + atomB->getType().c_str(), + atomC->getType().c_str()); + + if (bendType == NULL) { + sprintf(painCave.errMsg, + "Can not find Matching Bend Type for[%s, %s, %s]", + atomA->getType().c_str(), + atomB->getType().c_str(), + atomC->getType().c_str()); + + painCave.isFatal = 1; + simError(); + } + } + + bend = new Bend(atomA, atomB, atomC, bendType); + + } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) { + int ghostIndex = stamp->getGhostVectorSource(); + int normalIndex = ghostIndex != bendAtoms[0] ? + bendAtoms[0] : bendAtoms[1]; + Atom* normalAtom = mol->getAtomAt(normalIndex) ; + DirectionalAtom* ghostAtom = dynamic_cast(mol->getAtomAt(ghostIndex)); + if (ghostAtom == NULL) { + sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom"); + painCave.isFatal = 1; + simError(); + } - painCave.isFatal = 1; - simError(); + if (stamp->hasOverride()) { + + try { + bendType = btParser.parseTypeAndPars(stamp->getOverrideType(), + stamp->getOverridePars() ); } + catch( OpenMDException e) { + sprintf(painCave.errMsg, "MoleculeCreator Error: %s " + "for molecule %s\n", + e.what(), mol->getType().c_str() ); + painCave.isFatal = 1; + simError(); + } + } else { + + bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), + "GHOST"); + + if (bendType == NULL) { + sprintf(painCave.errMsg, + "Can not find Matching Bend Type for[%s, %s, %s]", + normalAtom->getType().c_str(), + ghostAtom->getType().c_str(), + "GHOST"); + + painCave.isFatal = 1; + simError(); + } + } + + bend = new GhostBend(normalAtom, ghostAtom, bendType); + + } + + //set the local index of this bend, the global index will be set later + bend->setLocalIndex(localIndexMan->getNextBendIndex()); + + // The rule for naming a bend is: MoleculeName_Bend_Integer + // The first part is the name of the molecule + // The second part is always fixed as "Bend" + // The third part is the index of the bend defined in meta-data file + // For example, Butane_Bend_0 is a valid Bend name in a butane molecule + + std::string s = OpenMD_itoa(mol->getNBends(), 10); + bend->setName(mol->getType() + "_Bend_" + s.c_str()); + return bend; + } - return new Bend(atomA, atomB, atomC, bendType); + Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, + TorsionStamp* stamp, + LocalIndexManager* localIndexMan) { + + TorsionTypeParser ttParser; + TorsionType* torsionType = NULL; + Torsion* torsion = NULL; + + std::vector torsionAtoms = stamp->getMembers(); + if (torsionAtoms.size() < 3) { + return torsion; } -} -Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) { + Atom* atomA = mol->getAtomAt(torsionAtoms[0]); + Atom* atomB = mol->getAtomAt(torsionAtoms[1]); + Atom* atomC = mol->getAtomAt(torsionAtoms[2]); - Atom* atomA = mol->getAtomAt(stamp->getA()); - Atom* atomB = mol->getAtomAt(stamp->getB()); - Atom* atomC = mol->getAtomAt(stamp->getC()); - Torsion* torsion; + if (torsionAtoms.size() == 4) { + Atom* atomD = mol->getAtomAt(torsionAtoms[3]); - if (stamp->getD() != -1) { - Atom* atomD = mol->getAtomAt(stamp->getD()); + assert(atomA && atomB && atomC && atomD ); - assert(atomA && atomB && atomC && atomD); + if (stamp->hasOverride()) { - TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), - atomC->getType(), atomD->getType()); + try { + torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(), + stamp->getOverridePars() ); + } + catch( OpenMDException e) { + sprintf(painCave.errMsg, "MoleculeCreator Error: %s " + "for molecule %s\n", + e.what(), mol->getType().c_str() ); + painCave.isFatal = 1; + simError(); + } + } else { + + torsionType = ff->getTorsionType(atomA->getType(), + atomB->getType(), + atomC->getType(), + atomD->getType()); if (torsionType == NULL) { - sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]", - atomA->getType().c_str(), - atomB->getType().c_str(), - atomC->getType().c_str(), - atomD->getType().c_str()); - - painCave.isFatal = 1; - simError(); + sprintf(painCave.errMsg, + "Can not find Matching Torsion Type for[%s, %s, %s, %s]", + atomA->getType().c_str(), + atomB->getType().c_str(), + atomC->getType().c_str(), + atomD->getType().c_str()); + + painCave.isFatal = 1; + simError(); } + } + + torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType); + } else { + + DirectionalAtom* dAtom = dynamic_cast(mol->getAtomAt(stamp->getGhostVectorSource())); + if (dAtom == NULL) { + sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom"); + painCave.isFatal = 1; + simError(); + } + + if (stamp->hasOverride()) { - torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType); + try { + torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(), + stamp->getOverridePars() ); + } + catch( OpenMDException e) { + sprintf(painCave.errMsg, "MoleculeCreator Error: %s " + "for molecule %s\n", + e.what(), mol->getType().c_str() ); + painCave.isFatal = 1; + simError(); + } + } else { + torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), + atomC->getType(), "GHOST"); + + if (torsionType == NULL) { + sprintf(painCave.errMsg, + "Can not find Matching Torsion Type for[%s, %s, %s, %s]", + atomA->getType().c_str(), + atomB->getType().c_str(), + atomC->getType().c_str(), + "GHOST"); + + painCave.isFatal = 1; + simError(); + } + } + + torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType); } - else { - DirectionalAtom* dAtom = dynamic_cast(atomC); - if (dAtom == NULL) { - sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom"); - painCave.isFatal = 1; - simError(); - } + //set the local index of this torsion, the global index will be set later + torsion->setLocalIndex(localIndexMan->getNextTorsionIndex()); + + // The rule for naming a torsion is: MoleculeName_Torsion_Integer + // The first part is the name of the molecule + // The second part is always fixed as "Torsion" + // The third part is the index of the torsion defined in meta-data file + // For example, Butane_Torsion_0 is a valid Torsion name in a + // butane molecule - TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), - atomC->getType(), "GHOST"); + std::string s = OpenMD_itoa(mol->getNTorsions(), 10); + torsion->setName(mol->getType() + "_Torsion_" + s.c_str()); + return torsion; + } - if (torsionType == NULL) { - sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]", - atomA->getType().c_str(), - atomB->getType().c_str(), - atomC->getType().c_str(), - "GHOST"); + Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol, + InversionStamp* stamp, + LocalIndexManager* localIndexMan) { - painCave.isFatal = 1; - simError(); - } + InversionTypeParser itParser; + InversionType* inversionType = NULL; + Inversion* inversion = NULL; + + int center = stamp->getCenter(); + std::vector satellites = stamp->getSatellites(); + if (satellites.size() != 3) { + return inversion; + } + + Atom* atomA = mol->getAtomAt(center); + Atom* atomB = mol->getAtomAt(satellites[0]); + Atom* atomC = mol->getAtomAt(satellites[1]); + Atom* atomD = mol->getAtomAt(satellites[2]); + + assert(atomA && atomB && atomC && atomD); + + if (stamp->hasOverride()) { + + try { + inversionType = itParser.parseTypeAndPars(stamp->getOverrideType(), + stamp->getOverridePars() ); + } + catch( OpenMDException e) { + sprintf(painCave.errMsg, "MoleculeCreator Error: %s " + "for molecule %s\n", + e.what(), mol->getType().c_str() ); + painCave.isFatal = 1; + simError(); + } + } else { + + inversionType = ff->getInversionType(atomA->getType(), + atomB->getType(), + atomC->getType(), + atomD->getType()); + + if (inversionType == NULL) { + sprintf(painCave.errMsg, + "No Matching Inversion Type for[%s, %s, %s, %s]\n" + "\t(May not be a problem: not all inversions are parametrized)\n", + atomA->getType().c_str(), + atomB->getType().c_str(), + atomC->getType().c_str(), + atomD->getType().c_str()); - torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + } } + if (inversionType != NULL) { + + inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType); + + // set the local index of this inversion, the global index will + // be set later + inversion->setLocalIndex(localIndexMan->getNextInversionIndex()); + + // The rule for naming an inversion is: MoleculeName_Inversion_Integer + // The first part is the name of the molecule + // The second part is always fixed as "Inversion" + // The third part is the index of the inversion defined in meta-data file + // For example, Benzene_Inversion_0 is a valid Inversion name in a + // Benzene molecule - return torsion; -} + std::string s = OpenMD_itoa(mol->getNInversions(), 10); + inversion->setName(mol->getType() + "_Inversion_" + s.c_str()); + return inversion; + } else { + return NULL; + } + } + -CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) { + CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, + CutoffGroupStamp* stamp, + LocalIndexManager* localIndexMan) { int nAtoms; CutoffGroup* cg; Atom* atom; @@ -442,71 +671,73 @@ CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu nAtoms = stamp->getNMembers(); for (int i =0; i < nAtoms; ++i) { - atom = mol->getAtomAt(stamp->getMember(i)); - assert(atom); - cg->addAtom(atom); + atom = mol->getAtomAt(stamp->getMemberAt(i)); + assert(atom); + cg->addAtom(atom); } - + + //set the local index of this cutoffGroup, global index will be set later + cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex()); + return cg; -} - -CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) { + } + + CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom, + LocalIndexManager* localIndexMan) { CutoffGroup* cg; cg = new CutoffGroup(); cg->addAtom(atom); + + //set the local index of this cutoffGroup, global index will be set later + cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex()); + return cg; -} + } -void MoleculeCreator::createConstraintPair(Molecule* mol) { + void MoleculeCreator::createConstraintPair(Molecule* mol) { //add bond constraints Molecule::BondIterator bi; Bond* bond; + ConstraintPair* cPair; + for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) { - BondType* bt = bond->getBondType(); + BondType* bt = bond->getBondType(); - //class Parent1 {}; - //class Child1 : public Parent {}; - //class Child2 : public Parent {}; - //Child1* ch1 = new Child1(); - //Child2* ch2 = dynamic_cast(ch1); - //the dynamic_cast is succeed in above line. A compiler bug? + if (typeid(FixedBondType) == typeid(*bt)) { + FixedBondType* fbt = dynamic_cast(bt); - if (typeid(FixedBondType) == typeid(*bt)) { - FixedBondType* fbt = dynamic_cast(bt); - - ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA()); - ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB()); - ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength()); - mol->addConstraintPair(consPair); - } + ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA()); + ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB()); + cPair = new ConstraintPair(consElemA, consElemB, + fbt->getEquilibriumBondLength(), false); + mol->addConstraintPair(cPair); + } } //rigidbody -- rigidbody constraint is not support yet -} + } -void MoleculeCreator::createConstraintElem(Molecule* mol) { + void MoleculeCreator::createConstraintElem(Molecule* mol) { ConstraintPair* consPair; Molecule::ConstraintPairIterator cpi; std::set sdSet; - for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) { + for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; + consPair = mol->nextConstraintPair(cpi)) { - StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble(); - if (sdSet.find(sdA) == sdSet.end()){ - sdSet.insert(sdA); - mol->addConstraintElem(new ConstraintElem(sdA)); - } - - StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble(); - if (sdSet.find(sdB) == sdSet.end()){ - sdSet.insert(sdB); - mol->addConstraintElem(new ConstraintElem(sdB)); - } - + StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble(); + if (sdSet.find(sdA) == sdSet.end()){ + sdSet.insert(sdA); + mol->addConstraintElem(new ConstraintElem(sdA)); + } + + StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble(); + if (sdSet.find(sdB) == sdSet.end()){ + sdSet.insert(sdB); + mol->addConstraintElem(new ConstraintElem(sdB)); + } } - + } } - -}