# | Line 50 | Line 50 | |
---|---|---|
50 | #include "brains/ForceManager.hpp" | |
51 | #include "primitives/Molecule.hpp" | |
52 | #include "UseTheForce/doForces_interface.h" | |
53 | < | #define __C |
53 | > | #define __OOPSE_C |
54 | #include "UseTheForce/DarkSide/fInteractionMap.h" | |
55 | #include "utils/simError.h" | |
56 | + | #include "primitives/Bond.hpp" |
57 | #include "primitives/Bend.hpp" | |
58 | < | #include "primitives/Bend.hpp" |
58 | > | #include "primitives/Torsion.hpp" |
59 | > | #include "primitives/Inversion.hpp" |
60 | namespace oopse { | |
61 | ||
62 | void ForceManager::calcForces(bool needPotential, bool needStress) { | |
# | Line 83 | Line 85 | namespace oopse { | |
85 | ||
86 | // forces are zeroed here, before any are accumulated. | |
87 | // NOTE: do not rezero the forces in Fortran. | |
88 | < | |
88 | > | |
89 | for (mol = info_->beginMolecule(mi); mol != NULL; | |
90 | mol = info_->nextMolecule(mi)) { | |
91 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
92 | atom->zeroForcesAndTorques(); | |
93 | } | |
94 | < | |
94 | > | |
95 | //change the positions of atoms which belong to the rigidbodies | |
96 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
97 | rb = mol->nextRigidBody(rbIter)) { | |
98 | rb->zeroForcesAndTorques(); | |
99 | } | |
100 | + | |
101 | } | |
102 | ||
103 | // Zero out the stress tensor | |
# | Line 108 | Line 111 | namespace oopse { | |
111 | Bond* bond; | |
112 | Bend* bend; | |
113 | Torsion* torsion; | |
114 | + | Inversion* inversion; |
115 | SimInfo::MoleculeIterator mi; | |
116 | Molecule::RigidBodyIterator rbIter; | |
117 | Molecule::BondIterator bondIter;; | |
118 | Molecule::BendIterator bendIter; | |
119 | Molecule::TorsionIterator torsionIter; | |
120 | + | Molecule::InversionIterator inversionIter; |
121 | RealType bondPotential = 0.0; | |
122 | RealType bendPotential = 0.0; | |
123 | RealType torsionPotential = 0.0; | |
124 | + | RealType inversionPotential = 0.0; |
125 | ||
126 | //calculate short range interactions | |
127 | for (mol = info_->beginMolecule(mi); mol != NULL; | |
# | Line 175 | Line 181 | namespace oopse { | |
181 | i->second.prev.potential = i->second.curr.potential; | |
182 | i->second.curr.angle = angle; | |
183 | i->second.curr.potential = currTorsionPot; | |
184 | + | i->second.deltaV = fabs(i->second.curr.potential - |
185 | + | i->second.prev.potential); |
186 | + | } |
187 | + | } |
188 | + | |
189 | + | for (inversion = mol->beginInversion(inversionIter); |
190 | + | inversion != NULL; |
191 | + | inversion = mol->nextInversion(inversionIter)) { |
192 | + | RealType angle; |
193 | + | inversion->calcForce(angle); |
194 | + | RealType currInversionPot = inversion->getPotential(); |
195 | + | inversionPotential += inversion->getPotential(); |
196 | + | std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion); |
197 | + | if (i == inversionDataSets.end()) { |
198 | + | InversionDataSet dataSet; |
199 | + | dataSet.prev.angle = dataSet.curr.angle = angle; |
200 | + | dataSet.prev.potential = dataSet.curr.potential = currInversionPot; |
201 | + | dataSet.deltaV = 0.0; |
202 | + | inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet)); |
203 | + | }else { |
204 | + | i->second.prev.angle = i->second.curr.angle; |
205 | + | i->second.prev.potential = i->second.curr.potential; |
206 | + | i->second.curr.angle = angle; |
207 | + | i->second.curr.potential = currInversionPot; |
208 | i->second.deltaV = fabs(i->second.curr.potential - | |
209 | i->second.prev.potential); | |
210 | } | |
# | Line 182 | Line 212 | namespace oopse { | |
212 | } | |
213 | ||
214 | RealType shortRangePotential = bondPotential + bendPotential + | |
215 | < | torsionPotential; |
215 | > | torsionPotential + inversionPotential; |
216 | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
217 | curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; | |
218 | curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential; | |
219 | curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential; | |
220 | curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential; | |
221 | + | curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential; |
222 | ||
223 | } | |
224 | ||
# | Line 201 | Line 232 | namespace oopse { | |
232 | RealType* A; | |
233 | RealType* electroFrame; | |
234 | RealType* rc; | |
235 | + | RealType* particlePot; |
236 | ||
237 | //get current snapshot from SimInfo | |
238 | curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
# | Line 212 | Line 244 | namespace oopse { | |
244 | trq = config->getArrayPointer(DataStorage::dslTorque); | |
245 | A = config->getArrayPointer(DataStorage::dslAmat); | |
246 | electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame); | |
247 | + | particlePot = config->getArrayPointer(DataStorage::dslParticlePot); |
248 | ||
249 | //calculate the center of mass of cutoff group | |
250 | SimInfo::MoleculeIterator mi; | |
# | Line 251 | Line 284 | namespace oopse { | |
284 | longRangePotential[i]=0.0; //Initialize array | |
285 | } | |
286 | ||
287 | < | doForceLoop( pos, |
288 | < | rc, |
289 | < | A, |
290 | < | electroFrame, |
291 | < | frc, |
292 | < | trq, |
293 | < | tau.getArrayPointer(), |
294 | < | longRangePotential, |
295 | < | &passedCalcPot, |
296 | < | &passedCalcStress, |
297 | < | &isError ); |
298 | < | |
287 | > | doForceLoop(pos, |
288 | > | rc, |
289 | > | A, |
290 | > | electroFrame, |
291 | > | frc, |
292 | > | trq, |
293 | > | tau.getArrayPointer(), |
294 | > | longRangePotential, |
295 | > | particlePot, |
296 | > | &passedCalcPot, |
297 | > | &passedCalcStress, |
298 | > | &isError ); |
299 | > | |
300 | if( isError ){ | |
301 | sprintf( painCave.errMsg, | |
302 | "Error returned from the fortran force calculation.\n" ); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |