| 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 | < | } | 
| 61 | < |  | 
| 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_); | 
| 72 | < | //integrableObjects_ don't own the objects | 
| 59 | > | } | 
| 60 | > |  | 
| 61 | > | Molecule::~Molecule() { | 
| 62 | > |  | 
| 63 | > | MemoryUtils::deletePointers(atoms_); | 
| 64 | > | MemoryUtils::deletePointers(bonds_); | 
| 65 | > | MemoryUtils::deletePointers(bends_); | 
| 66 | > | MemoryUtils::deletePointers(torsions_); | 
| 67 | > | MemoryUtils::deletePointers(rigidBodies_); | 
| 68 | > | MemoryUtils::deletePointers(cutoffGroups_); | 
| 69 | > | MemoryUtils::deletePointers(constraintPairs_); | 
| 70 | > | MemoryUtils::deletePointers(constraintElems_); | 
| 71 | > | // integrableObjects_ don't own the objects | 
| 72 |  | integrableObjects_.clear(); | 
| 73 |  |  | 
| 74 | < | } | 
| 75 | < |  | 
| 76 | < | void Molecule::addAtom(Atom* atom) { | 
| 74 | > | } | 
| 75 | > |  | 
| 76 | > | void Molecule::addAtom(Atom* atom) { | 
| 77 |  | if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) { | 
| 78 | < | atoms_.push_back(atom); | 
| 78 | > | atoms_.push_back(atom); | 
| 79 |  | } | 
| 80 | < | } | 
| 81 | < |  | 
| 82 | < | void Molecule::addBond(Bond* bond) { | 
| 80 | > | } | 
| 81 | > |  | 
| 82 | > | void Molecule::addBond(Bond* bond) { | 
| 83 |  | if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) { | 
| 84 | < | bonds_.push_back(bond); | 
| 84 | > | bonds_.push_back(bond); | 
| 85 |  | } | 
| 86 | < | } | 
| 87 | < |  | 
| 88 | < | void Molecule::addBend(Bend* bend) { | 
| 86 | > | } | 
| 87 | > |  | 
| 88 | > | void Molecule::addBend(Bend* bend) { | 
| 89 |  | if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) { | 
| 90 | < | bends_.push_back(bend); | 
| 90 | > | bends_.push_back(bend); | 
| 91 |  | } | 
| 92 | < | } | 
| 93 | < |  | 
| 94 | < | void Molecule::addTorsion(Torsion* torsion) { | 
| 95 | < | if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) { | 
| 96 | < | torsions_.push_back(torsion); | 
| 92 | > | } | 
| 93 | > |  | 
| 94 | > | void Molecule::addTorsion(Torsion* torsion) { | 
| 95 | > | if (std::find(torsions_.begin(), torsions_.end(), torsion) == | 
| 96 | > | torsions_.end()) { | 
| 97 | > | torsions_.push_back(torsion); | 
| 98 |  | } | 
| 99 | < | } | 
| 100 | < |  | 
| 101 | < | void Molecule::addRigidBody(RigidBody *rb) { | 
| 102 | < | if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) { | 
| 103 | < | rigidBodies_.push_back(rb); | 
| 99 | > | } | 
| 100 | > |  | 
| 101 | > | void Molecule::addRigidBody(RigidBody *rb) { | 
| 102 | > | if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == | 
| 103 | > | rigidBodies_.end()) { | 
| 104 | > | rigidBodies_.push_back(rb); | 
| 105 |  | } | 
| 106 | < | } | 
| 107 | < |  | 
| 108 | < | void Molecule::addCutoffGroup(CutoffGroup* cp) { | 
| 109 | < | if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) { | 
| 110 | < | cutoffGroups_.push_back(cp); | 
| 106 | > | } | 
| 107 | > |  | 
| 108 | > | void Molecule::addCutoffGroup(CutoffGroup* cp) { | 
| 109 | > | if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == | 
| 110 | > | cutoffGroups_.end()) { | 
| 111 | > | cutoffGroups_.push_back(cp); | 
| 112 | > | } | 
| 113 | > | } | 
| 114 | > |  | 
| 115 | > | void Molecule::addConstraintPair(ConstraintPair* cp) { | 
| 116 | > | if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == | 
| 117 | > | constraintPairs_.end()) { | 
| 118 | > | constraintPairs_.push_back(cp); | 
| 119 | > | } | 
| 120 | > | } | 
| 121 | > |  | 
| 122 | > | void Molecule::addConstraintElem(ConstraintElem* cp) { | 
| 123 | > | if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == | 
| 124 | > | constraintElems_.end()) { | 
| 125 | > | constraintElems_.push_back(cp); | 
| 126 |  | } | 
| 127 | < |  | 
| 128 | < | } | 
| 129 | < |  | 
| 114 | < | void Molecule::addConstraintPair(ConstraintPair* cp) { | 
| 115 | < | if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) { | 
| 116 | < | constraintPairs_.push_back(cp); | 
| 117 | < | } | 
| 118 | < |  | 
| 119 | < | } | 
| 120 | < |  | 
| 121 | < | void Molecule::addConstraintElem(ConstraintElem* cp) { | 
| 122 | < | if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) { | 
| 123 | < | constraintElems_.push_back(cp); | 
| 124 | < | } | 
| 125 | < |  | 
| 126 | < | } | 
| 127 | < |  | 
| 128 | < | void Molecule::complete() { | 
| 127 | > | } | 
| 128 | > |  | 
| 129 | > | void Molecule::complete() { | 
| 130 |  |  | 
| 131 |  | std::set<Atom*> rigidAtoms; | 
| 132 |  | RigidBody* rb; | 
| 134 |  |  | 
| 135 |  |  | 
| 136 |  | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { | 
| 137 | < | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); | 
| 137 | > | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); | 
| 138 |  | } | 
| 139 | < |  | 
| 139 | > |  | 
| 140 |  | Atom* atom; | 
| 141 |  | AtomIterator ai; | 
| 142 |  | for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { | 
| 143 | < |  | 
| 144 | < | if (rigidAtoms.find(*ai) == rigidAtoms.end()) { | 
| 145 | < | //if an atom does not belong to a rigid body, it is an integrable object | 
| 146 | < | integrableObjects_.push_back(*ai); | 
| 147 | < | } | 
| 147 | < | } | 
| 143 | > |  | 
| 144 | > | if (rigidAtoms.find(*ai) == rigidAtoms.end()) { | 
| 145 | > |  | 
| 146 | > | // If an atom does not belong to a rigid body, it is an | 
| 147 | > | // integrable object | 
| 148 |  |  | 
| 149 | + | integrableObjects_.push_back(*ai); | 
| 150 | + | } | 
| 151 | + | } | 
| 152 | + |  | 
| 153 |  | //find all free atoms (which do not belong to rigid bodies) | 
| 154 | < | //performs the "difference" operation from set theory,  the output range contains a copy of every | 
| 155 | < | //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in | 
| 156 | < | //[rigidAtoms.begin(), rigidAtoms.end()). | 
| 154 | > | // performs the "difference" operation from set theory, the output | 
| 155 | > | // range contains a copy of every element that is contained in | 
| 156 | > | // [allAtoms.begin(), allAtoms.end()) and not contained in | 
| 157 | > | // [rigidAtoms.begin(), rigidAtoms.end()). | 
| 158 |  | //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(), | 
| 159 |  | //                        std::back_inserter(integrableObjects_)); | 
| 160 |  |  | 
| 169 |  | integrableObjects_.push_back(rb); | 
| 170 |  | } | 
| 171 |  | //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); | 
| 172 | < | } | 
| 172 | > | } | 
| 173 |  |  | 
| 174 | < | double Molecule::getMass() { | 
| 174 | > | RealType Molecule::getMass() { | 
| 175 |  | StuntDouble* sd; | 
| 176 |  | std::vector<StuntDouble*>::iterator i; | 
| 177 | < | double mass = 0.0; | 
| 178 | < |  | 
| 179 | < | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 180 | < | mass += sd->getMass(); | 
| 177 | > | RealType mass = 0.0; | 
| 178 | > |  | 
| 179 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = | 
| 180 | > | nextIntegrableObject(i)){ | 
| 181 | > | mass += sd->getMass(); | 
| 182 |  | } | 
| 183 | + |  | 
| 184 | + | return mass; | 
| 185 | + | } | 
| 186 |  |  | 
| 187 | < | return mass; | 
| 179 | < |  | 
| 180 | < | } | 
| 181 | < |  | 
| 182 | < | Vector3d Molecule::getCom() { | 
| 187 | > | Vector3d Molecule::getCom() { | 
| 188 |  | StuntDouble* sd; | 
| 189 |  | std::vector<StuntDouble*>::iterator i; | 
| 190 |  | Vector3d com; | 
| 191 | < | double totalMass = 0; | 
| 192 | < | double mass; | 
| 191 | > | RealType totalMass = 0; | 
| 192 | > | RealType mass; | 
| 193 |  |  | 
| 194 | < | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 195 | < | mass = sd->getMass(); | 
| 196 | < | totalMass += mass; | 
| 197 | < | com += sd->getPos() * mass; | 
| 194 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = | 
| 195 | > | nextIntegrableObject(i)){ | 
| 196 | > | mass = sd->getMass(); | 
| 197 | > | totalMass += mass; | 
| 198 | > | com += sd->getPos() * mass; | 
| 199 |  | } | 
| 200 | < |  | 
| 200 | > |  | 
| 201 |  | com /= totalMass; | 
| 202 |  |  | 
| 203 |  | return com; | 
| 204 | < | } | 
| 204 | > | } | 
| 205 |  |  | 
| 206 | < | void Molecule::moveCom(const Vector3d& delta) { | 
| 206 | > | void Molecule::moveCom(const Vector3d& delta) { | 
| 207 |  | StuntDouble* sd; | 
| 208 |  | std::vector<StuntDouble*>::iterator i; | 
| 209 |  |  | 
| 210 | < | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 211 | < | sd->setPos(sd->getPos() + delta); | 
| 212 | < | } | 
| 210 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = | 
| 211 | > | nextIntegrableObject(i)){ | 
| 212 | > | sd->setPos(sd->getPos() + delta); | 
| 213 | > | } | 
| 214 | > | } | 
| 215 |  |  | 
| 216 | < | } | 
| 209 | < |  | 
| 210 | < | Vector3d Molecule::getComVel() { | 
| 216 | > | Vector3d Molecule::getComVel() { | 
| 217 |  | StuntDouble* sd; | 
| 218 |  | std::vector<StuntDouble*>::iterator i; | 
| 219 |  | Vector3d velCom; | 
| 220 | < | double totalMass = 0; | 
| 221 | < | double mass; | 
| 220 | > | RealType totalMass = 0; | 
| 221 | > | RealType mass; | 
| 222 |  |  | 
| 223 | < | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ | 
| 224 | < | mass = sd->getMass(); | 
| 225 | < | totalMass += mass; | 
| 226 | < | velCom += sd->getVel() * mass; | 
| 223 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = | 
| 224 | > | nextIntegrableObject(i)){ | 
| 225 | > | mass = sd->getMass(); | 
| 226 | > | totalMass += mass; | 
| 227 | > | velCom += sd->getVel() * mass; | 
| 228 |  | } | 
| 229 | < |  | 
| 229 | > |  | 
| 230 |  | velCom /= totalMass; | 
| 231 | < |  | 
| 231 | > |  | 
| 232 |  | return velCom; | 
| 233 | < | } | 
| 233 | > | } | 
| 234 |  |  | 
| 235 | < | double Molecule::getPotential() { | 
| 235 | > | RealType Molecule::getPotential() { | 
| 236 |  |  | 
| 237 |  | Bond* bond; | 
| 238 |  | Bend* bend; | 
| 241 |  | Molecule::BendIterator  bendIter; | 
| 242 |  | Molecule::TorsionIterator  torsionIter; | 
| 243 |  |  | 
| 244 | < | double potential = 0.0; | 
| 244 | > | RealType potential = 0.0; | 
| 245 |  |  | 
| 246 |  | for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { | 
| 247 | < | potential += bond->getPotential(); | 
| 247 | > | potential += bond->getPotential(); | 
| 248 |  | } | 
| 249 |  |  | 
| 250 |  | for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { | 
| 251 | < | potential += bend->getPotential(); | 
| 251 | > | potential += bend->getPotential(); | 
| 252 |  | } | 
| 253 |  |  | 
| 254 | < | for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) { | 
| 255 | < | potential += torsion->getPotential(); | 
| 254 | > | for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = | 
| 255 | > | nextTorsion(torsionIter)) { | 
| 256 | > | potential += torsion->getPotential(); | 
| 257 |  | } | 
| 258 | < |  | 
| 258 | > |  | 
| 259 |  | return potential; | 
| 260 | < |  | 
| 261 | < | } | 
| 262 | < |  | 
| 263 | < | std::ostream& operator <<(std::ostream& o, Molecule& mol) { | 
| 260 | > |  | 
| 261 | > | } | 
| 262 | > |  | 
| 263 | > | std::ostream& operator <<(std::ostream& o, Molecule& mol) { | 
| 264 |  | o << std::endl; | 
| 265 |  | o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; | 
| 266 |  | o << mol.getNAtoms() << " atoms" << std::endl; | 
| 272 |  | o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; | 
| 273 |  | o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; | 
| 274 |  | return o; | 
| 275 | < | } | 
| 276 | < |  | 
| 275 | > | } | 
| 276 | > |  | 
| 277 |  | }//end namespace oopse |