| 1 | < | /* | 
| 1 | > | /* | 
| 2 |  | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 |  | * | 
| 4 |  | * The University of Notre Dame grants you ("Licensee") a | 
| 54 |  | #include "utils/simError.h" | 
| 55 |  |  | 
| 56 |  | namespace oopse { | 
| 57 | < | Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) | 
| 57 | > | Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) | 
| 58 |  | : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) { | 
| 59 |  |  | 
| 60 | < | } | 
| 60 | > | } | 
| 61 |  |  | 
| 62 | < | Molecule::~Molecule() { | 
| 62 | > | Molecule::~Molecule() { | 
| 63 |  |  | 
| 64 | < | MemoryUtils::deleteVectorOfPointer(atoms_); | 
| 65 | < | MemoryUtils::deleteVectorOfPointer(bonds_); | 
| 66 | < | MemoryUtils::deleteVectorOfPointer(bends_); | 
| 67 | < | MemoryUtils::deleteVectorOfPointer(torsions_); | 
| 68 | < | MemoryUtils::deleteVectorOfPointer(rigidBodies_); | 
| 69 | < | MemoryUtils::deleteVectorOfPointer(cutoffGroups_); | 
| 70 | < | MemoryUtils::deleteVectorOfPointer(constraintPairs_); | 
| 71 | < | MemoryUtils::deleteVectorOfPointer(constraintElems_); | 
| 64 | > | MemoryUtils::deletePointers(atoms_); | 
| 65 | > | MemoryUtils::deletePointers(bonds_); | 
| 66 | > | MemoryUtils::deletePointers(bends_); | 
| 67 | > | MemoryUtils::deletePointers(torsions_); | 
| 68 | > | MemoryUtils::deletePointers(rigidBodies_); | 
| 69 | > | MemoryUtils::deletePointers(cutoffGroups_); | 
| 70 | > | MemoryUtils::deletePointers(constraintPairs_); | 
| 71 | > | MemoryUtils::deletePointers(constraintElems_); | 
| 72 |  | //integrableObjects_ don't own the objects | 
| 73 |  | integrableObjects_.clear(); | 
| 74 |  |  | 
| 75 | < | } | 
| 75 | > | } | 
| 76 |  |  | 
| 77 | < | void Molecule::addAtom(Atom* atom) { | 
| 77 | > | void Molecule::addAtom(Atom* atom) { | 
| 78 |  | if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) { | 
| 79 | < | atoms_.push_back(atom); | 
| 79 | > | atoms_.push_back(atom); | 
| 80 |  | } | 
| 81 | < | } | 
| 81 | > | } | 
| 82 |  |  | 
| 83 | < | void Molecule::addBond(Bond* bond) { | 
| 83 | > | void Molecule::addBond(Bond* bond) { | 
| 84 |  | if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) { | 
| 85 | < | bonds_.push_back(bond); | 
| 85 | > | bonds_.push_back(bond); | 
| 86 |  | } | 
| 87 | < | } | 
| 87 | > | } | 
| 88 |  |  | 
| 89 | < | void Molecule::addBend(Bend* bend) { | 
| 89 | > | void Molecule::addBend(Bend* bend) { | 
| 90 |  | if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) { | 
| 91 | < | bends_.push_back(bend); | 
| 91 | > | bends_.push_back(bend); | 
| 92 |  | } | 
| 93 | < | } | 
| 93 | > | } | 
| 94 |  |  | 
| 95 | < | void Molecule::addTorsion(Torsion* torsion) { | 
| 95 | > | void Molecule::addTorsion(Torsion* torsion) { | 
| 96 |  | if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) { | 
| 97 | < | torsions_.push_back(torsion); | 
| 97 | > | torsions_.push_back(torsion); | 
| 98 |  | } | 
| 99 | < | } | 
| 99 | > | } | 
| 100 |  |  | 
| 101 | < | void Molecule::addRigidBody(RigidBody *rb) { | 
| 101 | > | void Molecule::addRigidBody(RigidBody *rb) { | 
| 102 |  | if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) { | 
| 103 | < | rigidBodies_.push_back(rb); | 
| 103 | > | rigidBodies_.push_back(rb); | 
| 104 |  | } | 
| 105 | < | } | 
| 105 | > | } | 
| 106 |  |  | 
| 107 | < | void Molecule::addCutoffGroup(CutoffGroup* cp) { | 
| 107 | > | void Molecule::addCutoffGroup(CutoffGroup* cp) { | 
| 108 |  | if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) { | 
| 109 | < | cutoffGroups_.push_back(cp); | 
| 109 | > | cutoffGroups_.push_back(cp); | 
| 110 |  | } | 
| 111 |  |  | 
| 112 | < | } | 
| 112 | > | } | 
| 113 |  |  | 
| 114 | < | void Molecule::addConstraintPair(ConstraintPair* cp) { | 
| 114 | > | void Molecule::addConstraintPair(ConstraintPair* cp) { | 
| 115 |  | if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) { | 
| 116 | < | constraintPairs_.push_back(cp); | 
| 116 | > | constraintPairs_.push_back(cp); | 
| 117 |  | } | 
| 118 |  |  | 
| 119 | < | } | 
| 119 | > | } | 
| 120 |  |  | 
| 121 | < | void Molecule::addConstraintElem(ConstraintElem* cp) { | 
| 121 | > | void Molecule::addConstraintElem(ConstraintElem* cp) { | 
| 122 |  | if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) { | 
| 123 | < | constraintElems_.push_back(cp); | 
| 123 | > | constraintElems_.push_back(cp); | 
| 124 |  | } | 
| 125 |  |  | 
| 126 | < | } | 
| 126 | > | } | 
| 127 |  |  | 
| 128 | < | void Molecule::complete() { | 
| 128 | > | void Molecule::complete() { | 
| 129 |  |  | 
| 130 |  | std::set<Atom*> rigidAtoms; | 
| 131 |  | RigidBody* rb; | 
| 133 |  |  | 
| 134 |  |  | 
| 135 |  | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { | 
| 136 | < | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); | 
| 136 | > | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); | 
| 137 |  | } | 
| 138 |  |  | 
| 139 |  | Atom* atom; | 
| 140 |  | AtomIterator ai; | 
| 141 |  | for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { | 
| 142 |  |  | 
| 143 | < | if (rigidAtoms.find(*ai) == rigidAtoms.end()) { | 
| 144 | < | //if an atom does not belong to a rigid body, it is an integrable object | 
| 145 | < | integrableObjects_.push_back(*ai); | 
| 146 | < | } | 
| 143 | > | if (rigidAtoms.find(*ai) == rigidAtoms.end()) { | 
| 144 | > | //if an atom does not belong to a rigid body, it is an integrable object | 
| 145 | > | integrableObjects_.push_back(*ai); | 
| 146 | > | } | 
| 147 |  | } | 
| 148 |  |  | 
| 149 |  | //find all free atoms (which do not belong to rigid bodies) | 
| 164 |  | integrableObjects_.push_back(rb); | 
| 165 |  | } | 
| 166 |  | //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); | 
| 167 | < | } | 
| 167 | > | } | 
| 168 |  |  | 
| 169 | < | double Molecule::getMass() { | 
| 169 | > | RealType Molecule::getMass() { | 
| 170 |  | StuntDouble* sd; | 
| 171 |  | std::vector<StuntDouble*>::iterator i; | 
| 172 | < | double mass = 0.0; | 
| 172 | > | RealType mass = 0.0; | 
| 173 |  |  | 
| 174 |  | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 175 | < | mass += sd->getMass(); | 
| 175 | > | mass += sd->getMass(); | 
| 176 |  | } | 
| 177 |  |  | 
| 178 |  | return mass; | 
| 179 |  |  | 
| 180 | < | } | 
| 180 | > | } | 
| 181 |  |  | 
| 182 | < | Vector3d Molecule::getCom() { | 
| 182 | > | Vector3d Molecule::getCom() { | 
| 183 |  | StuntDouble* sd; | 
| 184 |  | std::vector<StuntDouble*>::iterator i; | 
| 185 |  | Vector3d com; | 
| 186 | < | double totalMass = 0; | 
| 187 | < | double mass; | 
| 186 | > | RealType totalMass = 0; | 
| 187 | > | RealType mass; | 
| 188 |  |  | 
| 189 |  | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 190 | < | mass = sd->getMass(); | 
| 191 | < | totalMass += mass; | 
| 192 | < | com += sd->getPos() * mass; | 
| 190 | > | mass = sd->getMass(); | 
| 191 | > | totalMass += mass; | 
| 192 | > | com += sd->getPos() * mass; | 
| 193 |  | } | 
| 194 |  |  | 
| 195 |  | com /= totalMass; | 
| 196 |  |  | 
| 197 |  | return com; | 
| 198 | < | } | 
| 198 | > | } | 
| 199 |  |  | 
| 200 | < | void Molecule::moveCom(const Vector3d& delta) { | 
| 200 | > | void Molecule::moveCom(const Vector3d& delta) { | 
| 201 |  | StuntDouble* sd; | 
| 202 |  | std::vector<StuntDouble*>::iterator i; | 
| 203 |  |  | 
| 204 |  | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 205 | < | sd->setPos(sd->getPos() + delta); | 
| 205 | > | sd->setPos(sd->getPos() + delta); | 
| 206 |  | } | 
| 207 |  |  | 
| 208 | < | } | 
| 208 | > | } | 
| 209 |  |  | 
| 210 | < | Vector3d Molecule::getComVel() { | 
| 210 | > | Vector3d Molecule::getComVel() { | 
| 211 |  | StuntDouble* sd; | 
| 212 |  | std::vector<StuntDouble*>::iterator i; | 
| 213 |  | Vector3d velCom; | 
| 214 | < | double totalMass = 0; | 
| 215 | < | double mass; | 
| 214 | > | RealType totalMass = 0; | 
| 215 | > | RealType mass; | 
| 216 |  |  | 
| 217 |  | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 218 | < | mass = sd->getMass(); | 
| 219 | < | totalMass += mass; | 
| 220 | < | velCom += sd->getVel() * mass; | 
| 218 | > | mass = sd->getMass(); | 
| 219 | > | totalMass += mass; | 
| 220 | > | velCom += sd->getVel() * mass; | 
| 221 |  | } | 
| 222 |  |  | 
| 223 |  | velCom /= totalMass; | 
| 224 |  |  | 
| 225 |  | return velCom; | 
| 226 | < | } | 
| 226 | > | } | 
| 227 |  |  | 
| 228 | < | double Molecule::getPotential() { | 
| 228 | > | RealType Molecule::getPotential() { | 
| 229 |  |  | 
| 230 |  | Bond* bond; | 
| 231 |  | Bend* bend; | 
| 234 |  | Molecule::BendIterator  bendIter; | 
| 235 |  | Molecule::TorsionIterator  torsionIter; | 
| 236 |  |  | 
| 237 | < | double potential = 0.0; | 
| 237 | > | RealType potential = 0.0; | 
| 238 |  |  | 
| 239 |  | for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { | 
| 240 | < | potential += bond->getPotential(); | 
| 240 | > | potential += bond->getPotential(); | 
| 241 |  | } | 
| 242 |  |  | 
| 243 |  | for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { | 
| 244 | < | potential += bend->getPotential(); | 
| 244 | > | potential += bend->getPotential(); | 
| 245 |  | } | 
| 246 |  |  | 
| 247 |  | for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) { | 
| 248 | < | potential += torsion->getPotential(); | 
| 248 | > | potential += torsion->getPotential(); | 
| 249 |  | } | 
| 250 |  |  | 
| 251 |  | return potential; | 
| 252 |  |  | 
| 253 | < | } | 
| 253 | > | } | 
| 254 |  |  | 
| 255 | < | std::ostream& operator <<(std::ostream& o, Molecule& mol) { | 
| 255 | > | std::ostream& operator <<(std::ostream& o, Molecule& mol) { | 
| 256 |  | o << std::endl; | 
| 257 |  | o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; | 
| 258 |  | o << mol.getNAtoms() << " atoms" << std::endl; | 
| 264 |  | o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; | 
| 265 |  | o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; | 
| 266 |  | return o; | 
| 267 | < | } | 
| 267 | > | } | 
| 268 |  |  | 
| 269 |  | }//end namespace oopse |