# | Line 53 | Line 53 | |
---|---|---|
53 | #define __C | |
54 | #include "UseTheForce/DarkSide/fInteractionMap.h" | |
55 | #include "utils/simError.h" | |
56 | + | #include "primitives/Bend.hpp" |
57 | + | #include "primitives/Bend.hpp" |
58 | namespace oopse { | |
59 | ||
60 | + | /* |
61 | + | struct BendOrderStruct { |
62 | + | Bend* bend; |
63 | + | BendDataSet dataSet; |
64 | + | }; |
65 | + | struct TorsionOrderStruct { |
66 | + | Torsion* torsion; |
67 | + | TorsionDataSet dataSet; |
68 | + | }; |
69 | + | |
70 | + | bool BendSortFunctor(const BendOrderStruct& b1, const BendOrderStruct& b2) { |
71 | + | return b1.dataSet.deltaV < b2.dataSet.deltaV; |
72 | + | } |
73 | + | |
74 | + | bool TorsionSortFunctor(const TorsionOrderStruct& t1, const TorsionOrderStruct& t2) { |
75 | + | return t1.dataSet.deltaV < t2.dataSet.deltaV; |
76 | + | } |
77 | + | */ |
78 | void ForceManager::calcForces(bool needPotential, bool needStress) { | |
79 | ||
80 | if (!info_->isFortranInitialized()) { | |
# | Line 68 | Line 88 | namespace oopse { | |
88 | calcLongRangeInteraction(needPotential, needStress); | |
89 | ||
90 | postCalculation(); | |
91 | < | |
91 | > | |
92 | > | /* |
93 | > | std::vector<BendOrderStruct> bendOrderStruct; |
94 | > | for(std::map<Bend*, BendDataSet>::iterator i = bendDataSets.begin(); i != bendDataSets.end(); ++i) { |
95 | > | BendOrderStruct tmp; |
96 | > | tmp.bend= const_cast<Bend*>(i->first); |
97 | > | tmp.dataSet = i->second; |
98 | > | bendOrderStruct.push_back(tmp); |
99 | > | } |
100 | > | |
101 | > | std::vector<TorsionOrderStruct> torsionOrderStruct; |
102 | > | for(std::map<Torsion*, TorsionDataSet>::iterator j = torsionDataSets.begin(); j != torsionDataSets.end(); ++j) { |
103 | > | TorsionOrderStruct tmp; |
104 | > | tmp.torsion = const_cast<Torsion*>(j->first); |
105 | > | tmp.dataSet = j->second; |
106 | > | torsionOrderStruct.push_back(tmp); |
107 | > | } |
108 | > | |
109 | > | std::sort(bendOrderStruct.begin(), bendOrderStruct.end(), std::ptr_fun(BendSortFunctor)); |
110 | > | std::sort(torsionOrderStruct.begin(), torsionOrderStruct.end(), std::ptr_fun(TorsionSortFunctor)); |
111 | > | for (std::vector<BendOrderStruct>::iterator k = bendOrderStruct.begin(); k != bendOrderStruct.end(); ++k) { |
112 | > | Bend* bend = k->bend; |
113 | > | std::cout << "Bend: atom1=" <<bend->getAtomA()->getGlobalIndex() << ",atom2 = "<< bend->getAtomB()->getGlobalIndex() << ",atom3="<<bend->getAtomC()->getGlobalIndex() << " "; |
114 | > | std::cout << "deltaV=" << k->dataSet.deltaV << ",p_theta=" << k->dataSet.prev.angle <<",p_pot=" << k->dataSet.prev.potential<< ",c_theta=" << k->dataSet.curr.angle << ", c_pot = " << k->dataSet.curr.potential <<std::endl; |
115 | > | } |
116 | > | for (std::vector<TorsionOrderStruct>::iterator l = torsionOrderStruct.begin(); l != torsionOrderStruct.end(); ++l) { |
117 | > | Torsion* torsion = l->torsion; |
118 | > | std::cout << "Torsion: atom1=" <<torsion->getAtomA()->getGlobalIndex() << ",atom2 = "<< torsion->getAtomB()->getGlobalIndex() << ",atom3="<<torsion->getAtomC()->getGlobalIndex() << ",atom4="<<torsion->getAtomD()->getGlobalIndex()<< " "; |
119 | > | std::cout << "deltaV=" << l->dataSet.deltaV << ",p_theta=" << l->dataSet.prev.angle <<",p_pot=" << l->dataSet.prev.potential<< ",c_theta=" << l->dataSet.curr.angle << ", c_pot = " << l->dataSet.curr.potential <<std::endl; |
120 | > | } |
121 | > | */ |
122 | } | |
123 | ||
124 | void ForceManager::preCalculation() { | |
# | Line 105 | Line 155 | namespace oopse { | |
155 | Molecule::BondIterator bondIter;; | |
156 | Molecule::BendIterator bendIter; | |
157 | Molecule::TorsionIterator torsionIter; | |
158 | + | RealType bondPotential = 0.0; |
159 | + | RealType bendPotential = 0.0; |
160 | + | RealType torsionPotential = 0.0; |
161 | ||
162 | //calculate short range interactions | |
163 | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | |
164 | ||
165 | //change the positions of atoms which belong to the rigidbodies | |
166 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
167 | < | rb->updateAtoms(); |
167 | > | rb->updateAtoms(); |
168 | } | |
169 | ||
170 | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
171 | < | bond->calcForce(); |
171 | > | bond->calcForce(); |
172 | > | bondPotential += bond->getPotential(); |
173 | } | |
174 | ||
175 | + | |
176 | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | |
177 | < | bend->calcForce(); |
177 | > | |
178 | > | RealType angle; |
179 | > | bend->calcForce(angle); |
180 | > | RealType currBendPot = bend->getPotential(); |
181 | > | bendPotential += bend->getPotential(); |
182 | > | std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend); |
183 | > | if (i == bendDataSets.end()) { |
184 | > | BendDataSet dataSet; |
185 | > | dataSet.prev.angle = dataSet.curr.angle = angle; |
186 | > | dataSet.prev.potential = dataSet.curr.potential = currBendPot; |
187 | > | dataSet.deltaV = 0.0; |
188 | > | bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet)); |
189 | > | }else { |
190 | > | i->second.prev.angle = i->second.curr.angle; |
191 | > | i->second.prev.potential = i->second.curr.potential; |
192 | > | i->second.curr.angle = angle; |
193 | > | i->second.curr.potential = currBendPot; |
194 | > | i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential); |
195 | > | } |
196 | } | |
197 | ||
198 | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
199 | < | torsion->calcForce(); |
199 | > | RealType angle; |
200 | > | torsion->calcForce(angle); |
201 | > | RealType currTorsionPot = torsion->getPotential(); |
202 | > | torsionPotential += torsion->getPotential(); |
203 | > | std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion); |
204 | > | if (i == torsionDataSets.end()) { |
205 | > | TorsionDataSet dataSet; |
206 | > | dataSet.prev.angle = dataSet.curr.angle = angle; |
207 | > | dataSet.prev.potential = dataSet.curr.potential = currTorsionPot; |
208 | > | dataSet.deltaV = 0.0; |
209 | > | torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet)); |
210 | > | }else { |
211 | > | i->second.prev.angle = i->second.curr.angle; |
212 | > | i->second.prev.potential = i->second.curr.potential; |
213 | > | i->second.curr.angle = angle; |
214 | > | i->second.curr.potential = currTorsionPot; |
215 | > | i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential); |
216 | > | } |
217 | } | |
218 | ||
219 | } | |
220 | ||
221 | < | double shortRangePotential = 0.0; |
132 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
133 | < | shortRangePotential += mol->getPotential(); |
134 | < | } |
135 | < | |
221 | > | RealType shortRangePotential = bondPotential + bendPotential + torsionPotential; |
222 | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
223 | curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; | |
224 | + | curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential; |
225 | + | curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential; |
226 | + | curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential; |
227 | + | |
228 | } | |
229 | ||
230 | void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) { | |
231 | Snapshot* curSnapshot; | |
232 | DataStorage* config; | |
233 | < | double* frc; |
234 | < | double* pos; |
235 | < | double* trq; |
236 | < | double* A; |
237 | < | double* electroFrame; |
238 | < | double* rc; |
233 | > | RealType* frc; |
234 | > | RealType* pos; |
235 | > | RealType* trq; |
236 | > | RealType* A; |
237 | > | RealType* electroFrame; |
238 | > | RealType* rc; |
239 | ||
240 | //get current snapshot from SimInfo | |
241 | curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
# | Line 182 | Line 272 | namespace oopse { | |
272 | } | |
273 | ||
274 | //initialize data before passing to fortran | |
275 | < | double longRangePotential[LR_POT_TYPES]; |
276 | < | double lrPot = 0.0; |
277 | < | |
275 | > | RealType longRangePotential[LR_POT_TYPES]; |
276 | > | RealType lrPot = 0.0; |
277 | > | Vector3d totalDipole; |
278 | Mat3x3d tau; | |
279 | short int passedCalcPot = needPotential; | |
280 | short int passedCalcStress = needStress; | |
# | Line 194 | Line 284 | namespace oopse { | |
284 | longRangePotential[i]=0.0; //Initialize array | |
285 | } | |
286 | ||
197 | – | |
198 | – | |
287 | doForceLoop( pos, | |
288 | rc, | |
289 | A, | |
# | Line 218 | Line 306 | namespace oopse { | |
306 | lrPot += longRangePotential[i]; //Quick hack | |
307 | } | |
308 | ||
309 | + | // grab the simulation box dipole moment if specified |
310 | + | if (info_->getCalcBoxDipole()){ |
311 | + | getAccumulatedBoxDipole(totalDipole.getArrayPointer()); |
312 | + | |
313 | + | curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0); |
314 | + | curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1); |
315 | + | curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2); |
316 | + | } |
317 | + | |
318 | //store the tau and long range potential | |
319 | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; | |
320 | < | // curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; |
320 | > | curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; |
321 | > | curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; |
322 | > | |
323 | curSnapshot->statData.setTau(tau); | |
324 | } | |
325 | ||
# | Line 236 | Line 335 | namespace oopse { | |
335 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
336 | rb->calcForcesAndTorques(); | |
337 | } | |
338 | < | } |
338 | > | } |
339 | ||
340 | } | |
341 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |