--- trunk/src/primitives/Molecule.cpp 2005/03/07 22:39:33 398 +++ branches/development/src/primitives/Molecule.cpp 2012/05/22 21:55:31 1715 @@ -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,6 +28,16 @@ * 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, 24107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -53,217 +54,266 @@ #include "utils/MemoryUtils.hpp" #include "utils/simError.h" -namespace oopse { -Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) +namespace OpenMD { + Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) { - -} - -Molecule::~Molecule() { - + } + + Molecule::~Molecule() { + MemoryUtils::deletePointers(atoms_); MemoryUtils::deletePointers(bonds_); MemoryUtils::deletePointers(bends_); MemoryUtils::deletePointers(torsions_); + MemoryUtils::deletePointers(inversions_); MemoryUtils::deletePointers(rigidBodies_); MemoryUtils::deletePointers(cutoffGroups_); MemoryUtils::deletePointers(constraintPairs_); MemoryUtils::deletePointers(constraintElems_); - //integrableObjects_ don't own the objects + + // integrableObjects_ don't own the objects integrableObjects_.clear(); + fluctuatingCharges_.clear(); -} - -void Molecule::addAtom(Atom* atom) { + } + + void Molecule::addAtom(Atom* atom) { if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) { - atoms_.push_back(atom); + atoms_.push_back(atom); } -} - -void Molecule::addBond(Bond* bond) { + } + + void Molecule::addBond(Bond* bond) { if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) { - bonds_.push_back(bond); + bonds_.push_back(bond); } -} - -void Molecule::addBend(Bend* bend) { + } + + void Molecule::addBend(Bend* bend) { if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) { - bends_.push_back(bend); + bends_.push_back(bend); } -} - -void Molecule::addTorsion(Torsion* torsion) { - if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) { - torsions_.push_back(torsion); + } + + void Molecule::addTorsion(Torsion* torsion) { + if (std::find(torsions_.begin(), torsions_.end(), torsion) == + torsions_.end()) { + torsions_.push_back(torsion); } -} + } -void Molecule::addRigidBody(RigidBody *rb) { - if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) { - rigidBodies_.push_back(rb); + void Molecule::addInversion(Inversion* inversion) { + if (std::find(inversions_.begin(), inversions_.end(), inversion) == + inversions_.end()) { + inversions_.push_back(inversion); } -} - -void Molecule::addCutoffGroup(CutoffGroup* cp) { - if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) { - cutoffGroups_.push_back(cp); + } + + void Molecule::addRigidBody(RigidBody *rb) { + if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == + rigidBodies_.end()) { + rigidBodies_.push_back(rb); } - -} - -void Molecule::addConstraintPair(ConstraintPair* cp) { - if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) { - constraintPairs_.push_back(cp); + } + + void Molecule::addCutoffGroup(CutoffGroup* cp) { + if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == + cutoffGroups_.end()) { + cutoffGroups_.push_back(cp); + } + } + + void Molecule::addConstraintPair(ConstraintPair* cp) { + if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == + constraintPairs_.end()) { + constraintPairs_.push_back(cp); + } + } + + void Molecule::addConstraintElem(ConstraintElem* cp) { + if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == + constraintElems_.end()) { + constraintElems_.push_back(cp); } - -} - -void Molecule::addConstraintElem(ConstraintElem* cp) { - if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) { - constraintElems_.push_back(cp); - } - -} - -void Molecule::complete() { + } + + void Molecule::complete() { std::set rigidAtoms; + Atom* atom; + AtomIterator ai; RigidBody* rb; - std::vector::iterator rbIter; + RigidBodyIterator rbIter; - + // Get list of all the atoms that are part of rigid bodies + for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { - rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); + rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); } + + // add any atom that wasn't part of a rigid body to the list of integrableObjects - Atom* atom; - AtomIterator ai; for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { - - if (rigidAtoms.find(*ai) == rigidAtoms.end()) { - //if an atom does not belong to a rigid body, it is an integrable object - integrableObjects_.push_back(*ai); - } - } + + if (rigidAtoms.find(atom) == rigidAtoms.end()) { - //find all free atoms (which do not belong to rigid bodies) - //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 - //[rigidAtoms.begin(), rigidAtoms.end()). - //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(), - // std::back_inserter(integrableObjects_)); + // If an atom does not belong to a rigid body, it is an + // integrable object - //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) { - // //Some atoms in rigidAtoms are not in allAtoms, something must be wrong - // sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule"); - // - // painCave.isFatal = 1; - // simError(); - //} - for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { - integrableObjects_.push_back(rb); - } - //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); -} + integrableObjects_.push_back(atom); + } + } + + // then add the rigid bodies themselves to the integrableObjects -double Molecule::getMass() { - StuntDouble* sd; - std::vector::iterator i; - double mass = 0.0; + for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { + integrableObjects_.push_back(rb); + } - for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ - mass += sd->getMass(); + // find the atoms that are fluctuating charges and add them to the + // fluctuatingCharges_ vector + + for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { + if ( atom->isFluctuatingCharge() ) + fluctuatingCharges_.push_back( atom ); } - return mass; + } -} + RealType Molecule::getMass() { + StuntDouble* sd; + std::vector::iterator i; + RealType mass = 0.0; + + for (sd = beginIntegrableObject(i); sd != NULL; sd = + nextIntegrableObject(i)){ + mass += sd->getMass(); + } + + return mass; + } -Vector3d Molecule::getCom() { + Vector3d Molecule::getCom() { StuntDouble* sd; std::vector::iterator i; Vector3d com; - double totalMass = 0; - double mass; + RealType totalMass = 0; + RealType mass; - for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ - mass = sd->getMass(); - totalMass += mass; - com += sd->getPos() * mass; + for (sd = beginIntegrableObject(i); sd != NULL; sd = + nextIntegrableObject(i)){ + mass = sd->getMass(); + totalMass += mass; + com += sd->getPos() * mass; } - + com /= totalMass; return com; -} + } -void Molecule::moveCom(const Vector3d& delta) { + void Molecule::moveCom(const Vector3d& delta) { StuntDouble* sd; std::vector::iterator i; - for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ - sd->setPos(sd->getPos() + delta); - } + for (sd = beginIntegrableObject(i); sd != NULL; sd = + nextIntegrableObject(i)){ + sd->setPos(sd->getPos() + delta); + } + } -} - -Vector3d Molecule::getComVel() { + Vector3d Molecule::getComVel() { StuntDouble* sd; std::vector::iterator i; Vector3d velCom; - double totalMass = 0; - double mass; + RealType totalMass = 0; + RealType mass; - for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ - mass = sd->getMass(); - totalMass += mass; - velCom += sd->getVel() * mass; + for (sd = beginIntegrableObject(i); sd != NULL; sd = + nextIntegrableObject(i)){ + mass = sd->getMass(); + totalMass += mass; + velCom += sd->getVel() * mass; } - + velCom /= totalMass; - + return velCom; -} + } -double Molecule::getPotential() { + RealType Molecule::getPotential() { Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; Molecule::BondIterator bondIter;; Molecule::BendIterator bendIter; Molecule::TorsionIterator torsionIter; + Molecule::InversionIterator inversionIter; - double potential = 0.0; + RealType potential = 0.0; for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { - potential += bond->getPotential(); + potential += bond->getPotential(); } for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { - potential += bend->getPotential(); + potential += bend->getPotential(); } - for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) { - potential += torsion->getPotential(); + for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = + nextTorsion(torsionIter)) { + potential += torsion->getPotential(); } - + + for (inversion = beginInversion(inversionIter); torsion != NULL; + inversion = nextInversion(inversionIter)) { + potential += inversion->getPotential(); + } + return potential; + + } + + void Molecule::addProperty(GenericData* genData) { + properties_.addProperty(genData); + } -} + void Molecule::removeProperty(const std::string& propName) { + properties_.removeProperty(propName); + } -std::ostream& operator <<(std::ostream& o, Molecule& mol) { + void Molecule::clearProperties() { + properties_.clearProperties(); + } + + std::vector Molecule::getPropertyNames() { + return properties_.getPropertyNames(); + } + + std::vector Molecule::getProperties() { + return properties_.getProperties(); + } + + GenericData* Molecule::getPropertyByName(const std::string& propName) { + return properties_.getPropertyByName(propName); + } + + std::ostream& operator <<(std::ostream& o, Molecule& mol) { o << std::endl; o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; o << mol.getNAtoms() << " atoms" << std::endl; o << mol.getNBonds() << " bonds" << std::endl; o << mol.getNBends() << " bends" << std::endl; o << mol.getNTorsions() << " torsions" << std::endl; + o << mol.getNInversions() << " inversions" << std::endl; o << mol.getNRigidBodies() << " rigid bodies" << std::endl; - o << mol.getNIntegrableObjects() << "integrable objects" << std::endl; - o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; - o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; + o << mol.getNIntegrableObjects() << " integrable objects" << std::endl; + o << mol.getNCutoffGroups() << " cutoff groups" << std::endl; + o << mol.getNConstraintPairs() << " constraint pairs" << std::endl; + o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl; return o; -} - -}//end namespace oopse + } + +}//end namespace OpenMD