# | Line 50 | Line 50 | |
---|---|---|
50 | #include "brains/ForceManager.hpp" | |
51 | #include "primitives/Molecule.hpp" | |
52 | #include "UseTheForce/doForces_interface.h" | |
53 | + | #define __C |
54 | + | #include "UseTheForce/DarkSide/fInteractionMap.h" |
55 | #include "utils/simError.h" | |
56 | + | #include "primitives/Bond.hpp" |
57 | + | #include "primitives/Bend.hpp" |
58 | namespace oopse { | |
59 | ||
60 | void ForceManager::calcForces(bool needPotential, bool needStress) { | |
61 | < | |
61 | > | |
62 | if (!info_->isFortranInitialized()) { | |
63 | info_->update(); | |
64 | } | |
65 | < | |
65 | > | |
66 | preCalculation(); | |
67 | ||
68 | calcShortRangeInteraction(); | |
69 | ||
70 | calcLongRangeInteraction(needPotential, needStress); | |
71 | ||
72 | < | postCalculation(); |
73 | < | |
72 | > | postCalculation(needStress); |
73 | > | |
74 | } | |
75 | < | |
75 | > | |
76 | void ForceManager::preCalculation() { | |
77 | SimInfo::MoleculeIterator mi; | |
78 | Molecule* mol; | |
# | Line 79 | Line 83 | namespace oopse { | |
83 | ||
84 | // forces are zeroed here, before any are accumulated. | |
85 | // NOTE: do not rezero the forces in Fortran. | |
86 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
86 | > | |
87 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
88 | > | mol = info_->nextMolecule(mi)) { |
89 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
90 | atom->zeroForcesAndTorques(); | |
91 | } | |
92 | < | |
92 | > | |
93 | //change the positions of atoms which belong to the rigidbodies | |
94 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
94 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
95 | > | rb = mol->nextRigidBody(rbIter)) { |
96 | rb->zeroForcesAndTorques(); | |
97 | } | |
98 | + | |
99 | } | |
100 | ||
101 | + | // Zero out the stress tensor |
102 | + | tau *= 0.0; |
103 | + | |
104 | } | |
105 | < | |
105 | > | |
106 | void ForceManager::calcShortRangeInteraction() { | |
107 | Molecule* mol; | |
108 | RigidBody* rb; | |
# | Line 103 | Line 114 | namespace oopse { | |
114 | Molecule::BondIterator bondIter;; | |
115 | Molecule::BendIterator bendIter; | |
116 | Molecule::TorsionIterator torsionIter; | |
117 | + | RealType bondPotential = 0.0; |
118 | + | RealType bendPotential = 0.0; |
119 | + | RealType torsionPotential = 0.0; |
120 | ||
121 | //calculate short range interactions | |
122 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
122 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
123 | > | mol = info_->nextMolecule(mi)) { |
124 | ||
125 | //change the positions of atoms which belong to the rigidbodies | |
126 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
127 | < | rb->updateAtoms(); |
126 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
127 | > | rb = mol->nextRigidBody(rbIter)) { |
128 | > | rb->updateAtoms(); |
129 | } | |
130 | ||
131 | < | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
132 | < | bond->calcForce(); |
131 | > | for (bond = mol->beginBond(bondIter); bond != NULL; |
132 | > | bond = mol->nextBond(bondIter)) { |
133 | > | bond->calcForce(); |
134 | > | bondPotential += bond->getPotential(); |
135 | } | |
136 | ||
137 | < | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
138 | < | bend->calcForce(); |
137 | > | for (bend = mol->beginBend(bendIter); bend != NULL; |
138 | > | bend = mol->nextBend(bendIter)) { |
139 | > | |
140 | > | RealType angle; |
141 | > | bend->calcForce(angle); |
142 | > | RealType currBendPot = bend->getPotential(); |
143 | > | bendPotential += bend->getPotential(); |
144 | > | std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend); |
145 | > | if (i == bendDataSets.end()) { |
146 | > | BendDataSet dataSet; |
147 | > | dataSet.prev.angle = dataSet.curr.angle = angle; |
148 | > | dataSet.prev.potential = dataSet.curr.potential = currBendPot; |
149 | > | dataSet.deltaV = 0.0; |
150 | > | bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet)); |
151 | > | }else { |
152 | > | i->second.prev.angle = i->second.curr.angle; |
153 | > | i->second.prev.potential = i->second.curr.potential; |
154 | > | i->second.curr.angle = angle; |
155 | > | i->second.curr.potential = currBendPot; |
156 | > | i->second.deltaV = fabs(i->second.curr.potential - |
157 | > | i->second.prev.potential); |
158 | > | } |
159 | } | |
160 | < | |
161 | < | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
162 | < | torsion->calcForce(); |
163 | < | } |
164 | < | |
160 | > | |
161 | > | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; |
162 | > | torsion = mol->nextTorsion(torsionIter)) { |
163 | > | RealType angle; |
164 | > | torsion->calcForce(angle); |
165 | > | RealType currTorsionPot = torsion->getPotential(); |
166 | > | torsionPotential += torsion->getPotential(); |
167 | > | std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion); |
168 | > | if (i == torsionDataSets.end()) { |
169 | > | TorsionDataSet dataSet; |
170 | > | dataSet.prev.angle = dataSet.curr.angle = angle; |
171 | > | dataSet.prev.potential = dataSet.curr.potential = currTorsionPot; |
172 | > | dataSet.deltaV = 0.0; |
173 | > | torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet)); |
174 | > | }else { |
175 | > | i->second.prev.angle = i->second.curr.angle; |
176 | > | i->second.prev.potential = i->second.curr.potential; |
177 | > | i->second.curr.angle = angle; |
178 | > | i->second.curr.potential = currTorsionPot; |
179 | > | i->second.deltaV = fabs(i->second.curr.potential - |
180 | > | i->second.prev.potential); |
181 | > | } |
182 | > | } |
183 | } | |
184 | ||
185 | < | double shortRangePotential = 0.0; |
186 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
131 | < | shortRangePotential += mol->getPotential(); |
132 | < | } |
133 | < | |
185 | > | RealType shortRangePotential = bondPotential + bendPotential + |
186 | > | torsionPotential; |
187 | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
188 | curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; | |
189 | + | curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential; |
190 | + | curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential; |
191 | + | curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential; |
192 | + | |
193 | } | |
194 | < | |
195 | < | void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) { |
194 | > | |
195 | > | void ForceManager::calcLongRangeInteraction(bool needPotential, |
196 | > | bool needStress) { |
197 | Snapshot* curSnapshot; | |
198 | DataStorage* config; | |
199 | < | double* frc; |
200 | < | double* pos; |
201 | < | double* trq; |
202 | < | double* A; |
203 | < | double* electroFrame; |
204 | < | double* rc; |
199 | > | RealType* frc; |
200 | > | RealType* pos; |
201 | > | RealType* trq; |
202 | > | RealType* A; |
203 | > | RealType* electroFrame; |
204 | > | RealType* rc; |
205 | > | RealType* particlePot; |
206 | ||
207 | //get current snapshot from SimInfo | |
208 | curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
209 | < | |
209 | > | |
210 | //get array pointers | |
211 | config = &(curSnapshot->atomData); | |
212 | frc = config->getArrayPointer(DataStorage::dslForce); | |
# | Line 155 | Line 214 | namespace oopse { | |
214 | trq = config->getArrayPointer(DataStorage::dslTorque); | |
215 | A = config->getArrayPointer(DataStorage::dslAmat); | |
216 | electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame); | |
217 | + | particlePot = config->getArrayPointer(DataStorage::dslParticlePot); |
218 | ||
219 | //calculate the center of mass of cutoff group | |
220 | SimInfo::MoleculeIterator mi; | |
# | Line 165 | Line 225 | namespace oopse { | |
225 | std::vector<Vector3d> rcGroup; | |
226 | ||
227 | if(info_->getNCutoffGroups() > 0){ | |
228 | < | |
229 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
230 | < | for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
228 | > | |
229 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
230 | > | mol = info_->nextMolecule(mi)) { |
231 | > | for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
232 | > | cg = mol->nextCutoffGroup(ci)) { |
233 | cg->getCOM(com); | |
234 | rcGroup.push_back(com); | |
235 | } | |
# | Line 175 | Line 237 | namespace oopse { | |
237 | ||
238 | rc = rcGroup[0].getArrayPointer(); | |
239 | } else { | |
240 | < | // center of mass of the group is the same as position of the atom if cutoff group does not exist |
240 | > | // center of mass of the group is the same as position of the atom |
241 | > | // if cutoff group does not exist |
242 | rc = pos; | |
243 | } | |
244 | < | |
244 | > | |
245 | //initialize data before passing to fortran | |
246 | < | double longRangePotential = 0.0; |
247 | < | Mat3x3d tau; |
246 | > | RealType longRangePotential[LR_POT_TYPES]; |
247 | > | RealType lrPot = 0.0; |
248 | > | Vector3d totalDipole; |
249 | short int passedCalcPot = needPotential; | |
250 | short int passedCalcStress = needStress; | |
251 | int isError = 0; | |
252 | ||
253 | < | doForceLoop( pos, |
254 | < | rc, |
255 | < | A, |
256 | < | electroFrame, |
257 | < | frc, |
258 | < | trq, |
259 | < | tau.getArrayPointer(), |
260 | < | &longRangePotential, |
261 | < | &passedCalcPot, |
262 | < | &passedCalcStress, |
263 | < | &isError ); |
264 | < | |
253 | > | for (int i=0; i<LR_POT_TYPES;i++){ |
254 | > | longRangePotential[i]=0.0; //Initialize array |
255 | > | } |
256 | > | |
257 | > | doForceLoop(pos, |
258 | > | rc, |
259 | > | A, |
260 | > | electroFrame, |
261 | > | frc, |
262 | > | trq, |
263 | > | tau.getArrayPointer(), |
264 | > | longRangePotential, |
265 | > | particlePot, |
266 | > | &passedCalcPot, |
267 | > | &passedCalcStress, |
268 | > | &isError ); |
269 | > | |
270 | if( isError ){ | |
271 | sprintf( painCave.errMsg, | |
272 | "Error returned from the fortran force calculation.\n" ); | |
273 | painCave.isFatal = 1; | |
274 | simError(); | |
275 | } | |
276 | < | |
276 | > | for (int i=0; i<LR_POT_TYPES;i++){ |
277 | > | lrPot += longRangePotential[i]; //Quick hack |
278 | > | } |
279 | > | |
280 | > | // grab the simulation box dipole moment if specified |
281 | > | if (info_->getCalcBoxDipole()){ |
282 | > | getAccumulatedBoxDipole(totalDipole.getArrayPointer()); |
283 | > | |
284 | > | curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0); |
285 | > | curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1); |
286 | > | curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2); |
287 | > | } |
288 | > | |
289 | //store the tau and long range potential | |
290 | < | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; |
291 | < | curSnapshot->statData.setTau(tau); |
290 | > | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; |
291 | > | curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; |
292 | > | curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; |
293 | } | |
294 | ||
295 | < | |
296 | < | void ForceManager::postCalculation() { |
295 | > | |
296 | > | void ForceManager::postCalculation(bool needStress) { |
297 | SimInfo::MoleculeIterator mi; | |
298 | Molecule* mol; | |
299 | Molecule::RigidBodyIterator rbIter; | |
300 | RigidBody* rb; | |
301 | + | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
302 | ||
303 | // collect the atomic forces onto rigid bodies | |
304 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
305 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
306 | < | rb->calcForcesAndTorques(); |
304 | > | |
305 | > | for (mol = info_->beginMolecule(mi); mol != NULL; |
306 | > | mol = info_->nextMolecule(mi)) { |
307 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
308 | > | rb = mol->nextRigidBody(rbIter)) { |
309 | > | if (needStress) { |
310 | > | Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial(); |
311 | > | tau += rbTau; |
312 | > | } else{ |
313 | > | rb->calcForcesAndTorques(); |
314 | > | } |
315 | } | |
316 | } | |
317 | ||
318 | + | if (needStress) { |
319 | + | #ifdef IS_MPI |
320 | + | Mat3x3d tmpTau(tau); |
321 | + | MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(), |
322 | + | 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
323 | + | #endif |
324 | + | curSnapshot->statData.setTau(tau); |
325 | + | } |
326 | } | |
327 | ||
328 | } //end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |