# | 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/Bend.hpp" |
57 | + | #include "primitives/Bend.hpp" |
58 | namespace oopse { | |
59 | ||
60 | + | struct BendOrderStruct { |
61 | + | Bend* bend; |
62 | + | BendDataSet dataSet; |
63 | + | }; |
64 | + | struct TorsionOrderStruct { |
65 | + | Torsion* torsion; |
66 | + | TorsionDataSet dataSet; |
67 | + | }; |
68 | + | |
69 | + | bool BendSortFunctor(const BendOrderStruct& b1, const BendOrderStruct& b2) { |
70 | + | return b1.dataSet.deltaV < b2.dataSet.deltaV; |
71 | + | } |
72 | + | |
73 | + | bool TorsionSortFunctor(const TorsionOrderStruct& t1, const TorsionOrderStruct& t2) { |
74 | + | return t1.dataSet.deltaV < t2.dataSet.deltaV; |
75 | + | } |
76 | + | |
77 | void ForceManager::calcForces(bool needPotential, bool needStress) { | |
78 | ||
79 | if (!info_->isFortranInitialized()) { | |
# | Line 66 | Line 87 | namespace oopse { | |
87 | calcLongRangeInteraction(needPotential, needStress); | |
88 | ||
89 | postCalculation(); | |
90 | < | |
90 | > | |
91 | > | std::vector<BendOrderStruct> bendOrderStruct; |
92 | > | for(std::map<Bend*, BendDataSet>::iterator i = bendDataSets.begin(); i != bendDataSets.end(); ++i) { |
93 | > | BendOrderStruct tmp; |
94 | > | tmp.bend= const_cast<Bend*>(i->first); |
95 | > | tmp.dataSet = i->second; |
96 | > | bendOrderStruct.push_back(tmp); |
97 | > | } |
98 | > | |
99 | > | std::vector<TorsionOrderStruct> torsionOrderStruct; |
100 | > | for(std::map<Torsion*, TorsionDataSet>::iterator j = torsionDataSets.begin(); j != torsionDataSets.end(); ++j) { |
101 | > | TorsionOrderStruct tmp; |
102 | > | tmp.torsion = const_cast<Torsion*>(j->first); |
103 | > | tmp.dataSet = j->second; |
104 | > | torsionOrderStruct.push_back(tmp); |
105 | > | } |
106 | > | |
107 | > | std::sort(bendOrderStruct.begin(), bendOrderStruct.end(), std::ptr_fun(BendSortFunctor)); |
108 | > | std::sort(torsionOrderStruct.begin(), torsionOrderStruct.end(), std::ptr_fun(TorsionSortFunctor)); |
109 | > | std::cout << "bend" << std::endl; |
110 | > | for (std::vector<BendOrderStruct>::iterator k = bendOrderStruct.begin(); k != bendOrderStruct.end(); ++k) { |
111 | > | Bend* bend = k->bend; |
112 | > | std::cout << "atom1=" <<bend->getAtomA()->getGlobalIndex() << ",atom2 = "<< bend->getAtomB()->getGlobalIndex() << ",atom3="<<bend->getAtomC()->getGlobalIndex() << " "; |
113 | > | 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; |
114 | > | } |
115 | > | std::cout << "torsio" << std::endl; |
116 | > | for (std::vector<TorsionOrderStruct>::iterator l = torsionOrderStruct.begin(); l != torsionOrderStruct.end(); ++l) { |
117 | > | Torsion* torsion = l->torsion; |
118 | > | std::cout << "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 103 | Line 155 | namespace oopse { | |
155 | Molecule::BondIterator bondIter;; | |
156 | Molecule::BendIterator bendIter; | |
157 | Molecule::TorsionIterator torsionIter; | |
158 | + | double bondPotential = 0.0; |
159 | + | double bendPotential = 0.0; |
160 | + | double 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 | + | //int i =0; |
176 | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | |
177 | < | bend->calcForce(); |
177 | > | //std::cout << i++ << "\t"; |
178 | > | double angle; |
179 | > | bend->calcForce(angle); |
180 | > | double 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 | > | double angle; |
200 | > | torsion->calcForce(angle); |
201 | > | double 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; |
130 | < | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
131 | < | shortRangePotential += mol->getPotential(); |
132 | < | } |
133 | < | |
221 | > | double 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) { | |
# | Line 180 | Line 272 | namespace oopse { | |
272 | } | |
273 | ||
274 | //initialize data before passing to fortran | |
275 | < | double longRangePotential = 0.0; |
275 | > | double longRangePotential[LR_POT_TYPES]; |
276 | > | double lrPot = 0.0; |
277 | > | |
278 | Mat3x3d tau; | |
279 | short int passedCalcPot = needPotential; | |
280 | short int passedCalcStress = needStress; | |
281 | int isError = 0; | |
282 | ||
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, | |
# | Line 193 | Line 291 | namespace oopse { | |
291 | frc, | |
292 | trq, | |
293 | tau.getArrayPointer(), | |
294 | < | &longRangePotential, |
294 | > | longRangePotential, |
295 | &passedCalcPot, | |
296 | &passedCalcStress, | |
297 | &isError ); | |
# | Line 204 | Line 302 | namespace oopse { | |
302 | painCave.isFatal = 1; | |
303 | simError(); | |
304 | } | |
305 | + | for (int i=0; i<LR_POT_TYPES;i++){ |
306 | + | lrPot += longRangePotential[i]; //Quick hack |
307 | + | } |
308 | ||
309 | //store the tau and long range potential | |
310 | < | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; |
310 | > | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; |
311 | > | curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; |
312 | > | curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; |
313 | > | |
314 | curSnapshot->statData.setTau(tau); | |
315 | } | |
316 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |