# | Line 50 | Line 50 | |
---|---|---|
50 | #include "brains/ForceManager.hpp" | |
51 | #include "primitives/Molecule.hpp" | |
52 | #include "UseTheForce/doForces_interface.h" | |
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/Torsion.hpp" |
59 | + | #include "primitives/Inversion.hpp" |
60 | namespace oopse { | |
61 | ||
62 | void ForceManager::calcForces(bool needPotential, bool needStress) { | |
63 | < | |
63 | > | |
64 | if (!info_->isFortranInitialized()) { | |
65 | info_->update(); | |
66 | } | |
67 | < | |
67 | > | |
68 | preCalculation(); | |
69 | ||
70 | calcShortRangeInteraction(); | |
71 | ||
72 | calcLongRangeInteraction(needPotential, needStress); | |
73 | ||
74 | < | postCalculation(); |
75 | < | |
74 | > | postCalculation(needStress); |
75 | > | |
76 | } | |
77 | < | |
77 | > | |
78 | void ForceManager::preCalculation() { | |
79 | SimInfo::MoleculeIterator mi; | |
80 | Molecule* mol; | |
# | Line 79 | Line 85 | namespace oopse { | |
85 | ||
86 | // forces are zeroed here, before any are accumulated. | |
87 | // NOTE: do not rezero the forces in Fortran. | |
88 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
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; rb = mol->nextRigidBody(rbIter)) { |
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 |
104 | + | tau *= 0.0; |
105 | + | |
106 | } | |
107 | < | |
107 | > | |
108 | void ForceManager::calcShortRangeInteraction() { | |
109 | Molecule* mol; | |
110 | RigidBody* rb; | |
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; mol = info_->nextMolecule(mi)) { |
127 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
128 | > | mol = info_->nextMolecule(mi)) { |
129 | ||
130 | //change the positions of atoms which belong to the rigidbodies | |
131 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
132 | < | rb->updateAtoms(); |
131 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
132 | > | rb = mol->nextRigidBody(rbIter)) { |
133 | > | rb->updateAtoms(); |
134 | } | |
135 | ||
136 | < | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
137 | < | bond->calcForce(); |
136 | > | for (bond = mol->beginBond(bondIter); bond != NULL; |
137 | > | bond = mol->nextBond(bondIter)) { |
138 | > | bond->calcForce(); |
139 | > | bondPotential += bond->getPotential(); |
140 | } | |
141 | ||
142 | < | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
143 | < | bend->calcForce(); |
142 | > | for (bend = mol->beginBend(bendIter); bend != NULL; |
143 | > | bend = mol->nextBend(bendIter)) { |
144 | > | |
145 | > | RealType angle; |
146 | > | bend->calcForce(angle); |
147 | > | RealType currBendPot = bend->getPotential(); |
148 | > | bendPotential += bend->getPotential(); |
149 | > | std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend); |
150 | > | if (i == bendDataSets.end()) { |
151 | > | BendDataSet dataSet; |
152 | > | dataSet.prev.angle = dataSet.curr.angle = angle; |
153 | > | dataSet.prev.potential = dataSet.curr.potential = currBendPot; |
154 | > | dataSet.deltaV = 0.0; |
155 | > | bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet)); |
156 | > | }else { |
157 | > | i->second.prev.angle = i->second.curr.angle; |
158 | > | i->second.prev.potential = i->second.curr.potential; |
159 | > | i->second.curr.angle = angle; |
160 | > | i->second.curr.potential = currBendPot; |
161 | > | i->second.deltaV = fabs(i->second.curr.potential - |
162 | > | i->second.prev.potential); |
163 | > | } |
164 | } | |
165 | + | |
166 | + | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; |
167 | + | torsion = mol->nextTorsion(torsionIter)) { |
168 | + | RealType angle; |
169 | + | torsion->calcForce(angle); |
170 | + | RealType currTorsionPot = torsion->getPotential(); |
171 | + | torsionPotential += torsion->getPotential(); |
172 | + | std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion); |
173 | + | if (i == torsionDataSets.end()) { |
174 | + | TorsionDataSet dataSet; |
175 | + | dataSet.prev.angle = dataSet.curr.angle = angle; |
176 | + | dataSet.prev.potential = dataSet.curr.potential = currTorsionPot; |
177 | + | dataSet.deltaV = 0.0; |
178 | + | torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet)); |
179 | + | }else { |
180 | + | i->second.prev.angle = i->second.curr.angle; |
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 (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
190 | < | torsion->calcForce(); |
191 | < | } |
192 | < | |
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 | > | } |
211 | > | } |
212 | } | |
213 | ||
214 | < | double shortRangePotential = 0.0; |
215 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
131 | < | shortRangePotential += mol->getPotential(); |
132 | < | } |
133 | < | |
214 | > | RealType shortRangePotential = bondPotential + bendPotential + |
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 | < | |
225 | < | void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) { |
224 | > | |
225 | > | void ForceManager::calcLongRangeInteraction(bool needPotential, |
226 | > | bool needStress) { |
227 | Snapshot* curSnapshot; | |
228 | DataStorage* config; | |
229 | < | double* frc; |
230 | < | double* pos; |
231 | < | double* trq; |
232 | < | double* A; |
233 | < | double* electroFrame; |
234 | < | double* rc; |
229 | > | RealType* frc; |
230 | > | RealType* pos; |
231 | > | RealType* trq; |
232 | > | RealType* A; |
233 | > | RealType* electroFrame; |
234 | > | RealType* rc; |
235 | > | RealType* particlePot; |
236 | ||
237 | //get current snapshot from SimInfo | |
238 | curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
239 | < | |
239 | > | |
240 | //get array pointers | |
241 | config = &(curSnapshot->atomData); | |
242 | frc = config->getArrayPointer(DataStorage::dslForce); | |
# | Line 155 | 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 163 | Line 253 | namespace oopse { | |
253 | CutoffGroup* cg; | |
254 | Vector3d com; | |
255 | std::vector<Vector3d> rcGroup; | |
256 | < | |
256 | > | |
257 | if(info_->getNCutoffGroups() > 0){ | |
258 | < | |
259 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
260 | < | for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
258 | > | |
259 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
260 | > | mol = info_->nextMolecule(mi)) { |
261 | > | for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
262 | > | cg = mol->nextCutoffGroup(ci)) { |
263 | cg->getCOM(com); | |
264 | rcGroup.push_back(com); | |
265 | } | |
# | Line 175 | Line 267 | namespace oopse { | |
267 | ||
268 | rc = rcGroup[0].getArrayPointer(); | |
269 | } else { | |
270 | < | // center of mass of the group is the same as position of the atom if cutoff group does not exist |
270 | > | // center of mass of the group is the same as position of the atom |
271 | > | // if cutoff group does not exist |
272 | rc = pos; | |
273 | } | |
274 | < | |
274 | > | |
275 | //initialize data before passing to fortran | |
276 | < | double longRangePotential = 0.0; |
277 | < | Mat3x3d tau; |
276 | > | RealType longRangePotential[LR_POT_TYPES]; |
277 | > | RealType lrPot = 0.0; |
278 | > | Vector3d totalDipole; |
279 | short int passedCalcPot = needPotential; | |
280 | short int passedCalcStress = needStress; | |
281 | int isError = 0; | |
282 | ||
283 | < | doForceLoop( pos, |
284 | < | rc, |
285 | < | A, |
286 | < | electroFrame, |
287 | < | frc, |
288 | < | trq, |
289 | < | tau.getArrayPointer(), |
290 | < | &longRangePotential, |
291 | < | &passedCalcPot, |
292 | < | &passedCalcStress, |
293 | < | &isError ); |
294 | < | |
283 | > | for (int i=0; i<LR_POT_TYPES;i++){ |
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 | > | particlePot, |
296 | > | &passedCalcPot, |
297 | > | &passedCalcStress, |
298 | > | &isError ); |
299 | > | |
300 | if( isError ){ | |
301 | sprintf( painCave.errMsg, | |
302 | "Error returned from the fortran force calculation.\n" ); | |
303 | painCave.isFatal = 1; | |
304 | simError(); | |
305 | } | |
306 | < | |
306 | > | for (int i=0; i<LR_POT_TYPES;i++){ |
307 | > | lrPot += longRangePotential[i]; //Quick hack |
308 | > | } |
309 | > | |
310 | > | // grab the simulation box dipole moment if specified |
311 | > | if (info_->getCalcBoxDipole()){ |
312 | > | getAccumulatedBoxDipole(totalDipole.getArrayPointer()); |
313 | > | |
314 | > | curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0); |
315 | > | curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1); |
316 | > | curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2); |
317 | > | } |
318 | > | |
319 | //store the tau and long range potential | |
320 | < | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; |
321 | < | curSnapshot->statData.setTau(tau); |
320 | > | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; |
321 | > | curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; |
322 | > | curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; |
323 | } | |
324 | ||
325 | < | |
326 | < | void ForceManager::postCalculation() { |
325 | > | |
326 | > | void ForceManager::postCalculation(bool needStress) { |
327 | SimInfo::MoleculeIterator mi; | |
328 | Molecule* mol; | |
329 | Molecule::RigidBodyIterator rbIter; | |
330 | RigidBody* rb; | |
331 | + | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
332 | ||
333 | // collect the atomic forces onto rigid bodies | |
334 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
335 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
336 | < | rb->calcForcesAndTorques(); |
334 | > | |
335 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
336 | > | mol = info_->nextMolecule(mi)) { |
337 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
338 | > | rb = mol->nextRigidBody(rbIter)) { |
339 | > | if (needStress) { |
340 | > | Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial(); |
341 | > | tau += rbTau; |
342 | > | } else{ |
343 | > | rb->calcForcesAndTorques(); |
344 | > | } |
345 | } | |
346 | } | |
347 | ||
348 | + | if (needStress) { |
349 | + | #ifdef IS_MPI |
350 | + | Mat3x3d tmpTau(tau); |
351 | + | MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(), |
352 | + | 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
353 | + | #endif |
354 | + | curSnapshot->statData.setTau(tau); |
355 | + | } |
356 | } | |
357 | ||
358 | } //end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |