| 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 | 
| 39 |  | * such damages. | 
| 40 |  | */ | 
| 41 |  |  | 
| 42 | < | /** | 
| 43 | < | * @file ForceManager.cpp | 
| 44 | < | * @author tlin | 
| 45 | < | * @date 11/09/2004 | 
| 46 | < | * @time 10:39am | 
| 47 | < | * @version 1.0 | 
| 48 | < | */ | 
| 42 | > | /** | 
| 43 | > | * @file ForceManager.cpp | 
| 44 | > | * @author tlin | 
| 45 | > | * @date 11/09/2004 | 
| 46 | > | * @time 10:39am | 
| 47 | > | * @version 1.0 | 
| 48 | > | */ | 
| 49 |  |  | 
| 50 |  | #include "brains/ForceManager.hpp" | 
| 51 |  | #include "primitives/Molecule.hpp" | 
| 53 |  | #include "utils/simError.h" | 
| 54 |  | namespace oopse { | 
| 55 |  |  | 
| 56 | < | void ForceManager::calcForces(bool needPotential, bool needStress) { | 
| 56 | > | void ForceManager::calcForces(bool needPotential, bool needStress) { | 
| 57 |  |  | 
| 58 |  | if (!info_->isFortranInitialized()) { | 
| 59 | < | info_->update(); | 
| 59 | > | info_->update(); | 
| 60 |  | } | 
| 61 |  |  | 
| 62 |  | preCalculation(); | 
| 67 |  |  | 
| 68 |  | postCalculation(); | 
| 69 |  |  | 
| 70 | < | } | 
| 70 | > | } | 
| 71 |  |  | 
| 72 | < | void ForceManager::preCalculation() { | 
| 72 | > | void ForceManager::preCalculation() { | 
| 73 |  | SimInfo::MoleculeIterator mi; | 
| 74 |  | Molecule* mol; | 
| 75 |  | Molecule::AtomIterator ai; | 
| 80 |  | // forces are zeroed here, before any are accumulated. | 
| 81 |  | // NOTE: do not rezero the forces in Fortran. | 
| 82 |  | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 83 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 84 | < | atom->zeroForcesAndTorques(); | 
| 85 | < | } | 
| 83 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 84 | > | atom->zeroForcesAndTorques(); | 
| 85 | > | } | 
| 86 |  |  | 
| 87 | < | //change the positions of atoms which belong to the rigidbodies | 
| 88 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 89 | < | rb->zeroForcesAndTorques(); | 
| 90 | < | } | 
| 87 | > | //change the positions of atoms which belong to the rigidbodies | 
| 88 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 89 | > | rb->zeroForcesAndTorques(); | 
| 90 | > | } | 
| 91 |  | } | 
| 92 |  |  | 
| 93 | < | } | 
| 93 | > | } | 
| 94 |  |  | 
| 95 | < | void ForceManager::calcShortRangeInteraction() { | 
| 95 | > | void ForceManager::calcShortRangeInteraction() { | 
| 96 |  | Molecule* mol; | 
| 97 |  | RigidBody* rb; | 
| 98 |  | Bond* bond; | 
| 107 |  | //calculate short range interactions | 
| 108 |  | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 109 |  |  | 
| 110 | < | //change the positions of atoms which belong to the rigidbodies | 
| 111 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 112 | < | rb->updateAtoms(); | 
| 113 | < | } | 
| 110 | > | //change the positions of atoms which belong to the rigidbodies | 
| 111 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 112 | > | rb->updateAtoms(); | 
| 113 | > | } | 
| 114 |  |  | 
| 115 | < | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 116 | < | bond->calcForce(); | 
| 117 | < | } | 
| 115 | > | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | 
| 116 | > | bond->calcForce(); | 
| 117 | > | } | 
| 118 |  |  | 
| 119 | < | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | 
| 120 | < | bend->calcForce(); | 
| 121 | < | } | 
| 119 | > | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | 
| 120 | > | bend->calcForce(); | 
| 121 | > | } | 
| 122 |  |  | 
| 123 | < | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 124 | < | torsion->calcForce(); | 
| 125 | < | } | 
| 123 | > | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | 
| 124 | > | torsion->calcForce(); | 
| 125 | > | } | 
| 126 |  |  | 
| 127 |  | } | 
| 128 |  |  | 
| 129 |  | double  shortRangePotential = 0.0; | 
| 130 |  | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 131 | < | shortRangePotential += mol->getPotential(); | 
| 131 | > | shortRangePotential += mol->getPotential(); | 
| 132 |  | } | 
| 133 |  |  | 
| 134 |  | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | 
| 135 |  | curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; | 
| 136 | < | } | 
| 136 | > | } | 
| 137 |  |  | 
| 138 | < | void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) { | 
| 138 | > | void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) { | 
| 139 |  | Snapshot* curSnapshot; | 
| 140 |  | DataStorage* config; | 
| 141 |  | double* frc; | 
| 166 |  |  | 
| 167 |  | if(info_->getNCutoffGroups() > 0){ | 
| 168 |  |  | 
| 169 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 169 | > | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 170 |  | for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 171 | < | cg->getCOM(com); | 
| 172 | < | rcGroup.push_back(com); | 
| 171 | > | cg->getCOM(com); | 
| 172 | > | rcGroup.push_back(com); | 
| 173 |  | } | 
| 174 | < | }// end for (mol) | 
| 174 | > | }// end for (mol) | 
| 175 |  |  | 
| 176 | < | rc = rcGroup[0].getArrayPointer(); | 
| 176 | > | rc = rcGroup[0].getArrayPointer(); | 
| 177 |  | } else { | 
| 178 | < | // center of mass of the group is the same as position of the atom  if cutoff group does not exist | 
| 179 | < | rc = pos; | 
| 178 | > | // center of mass of the group is the same as position of the atom  if cutoff group does not exist | 
| 179 | > | rc = pos; | 
| 180 |  | } | 
| 181 |  |  | 
| 182 |  | //initialize data before passing to fortran | 
| 187 |  | int isError = 0; | 
| 188 |  |  | 
| 189 |  | doForceLoop( pos, | 
| 190 | < | rc, | 
| 191 | < | A, | 
| 192 | < | electroFrame, | 
| 193 | < | frc, | 
| 194 | < | trq, | 
| 195 | < | tau.getArrayPointer(), | 
| 196 | < | &longRangePotential, | 
| 197 | < | &passedCalcPot, | 
| 198 | < | &passedCalcStress, | 
| 199 | < | &isError ); | 
| 190 | > | rc, | 
| 191 | > | A, | 
| 192 | > | electroFrame, | 
| 193 | > | frc, | 
| 194 | > | trq, | 
| 195 | > | tau.getArrayPointer(), | 
| 196 | > | &longRangePotential, | 
| 197 | > | &passedCalcPot, | 
| 198 | > | &passedCalcStress, | 
| 199 | > | &isError ); | 
| 200 |  |  | 
| 201 |  | if( isError ){ | 
| 202 | < | sprintf( painCave.errMsg, | 
| 203 | < | "Error returned from the fortran force calculation.\n" ); | 
| 204 | < | painCave.isFatal = 1; | 
| 205 | < | simError(); | 
| 202 | > | sprintf( painCave.errMsg, | 
| 203 | > | "Error returned from the fortran force calculation.\n" ); | 
| 204 | > | painCave.isFatal = 1; | 
| 205 | > | simError(); | 
| 206 |  | } | 
| 207 |  |  | 
| 208 |  | //store the tau and long range potential | 
| 209 |  | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; | 
| 210 |  | curSnapshot->statData.setTau(tau); | 
| 211 | < | } | 
| 211 | > | } | 
| 212 |  |  | 
| 213 |  |  | 
| 214 | < | void ForceManager::postCalculation() { | 
| 214 | > | void ForceManager::postCalculation() { | 
| 215 |  | SimInfo::MoleculeIterator mi; | 
| 216 |  | Molecule* mol; | 
| 217 |  | Molecule::RigidBodyIterator rbIter; | 
| 219 |  |  | 
| 220 |  | // collect the atomic forces onto rigid bodies | 
| 221 |  | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | 
| 222 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 223 | < | rb->calcForcesAndTorques(); | 
| 224 | < | } | 
| 222 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 223 | > | rb->calcForcesAndTorques(); | 
| 224 | > | } | 
| 225 |  | } | 
| 226 |  |  | 
| 227 | < | } | 
| 227 | > | } | 
| 228 |  |  | 
| 229 |  | } //end namespace oopse |