# | 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 | > | std::cout << "bend" << std::endl; |
112 | > | for (std::vector<BendOrderStruct>::iterator k = bendOrderStruct.begin(); k != bendOrderStruct.end(); ++k) { |
113 | > | Bend* bend = k->bend; |
114 | > | std::cout << "atom1=" <<bend->getAtomA()->getGlobalIndex() << ",atom2 = "<< bend->getAtomB()->getGlobalIndex() << ",atom3="<<bend->getAtomC()->getGlobalIndex() << " "; |
115 | > | 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; |
116 | > | } |
117 | > | std::cout << "torsio" << std::endl; |
118 | > | for (std::vector<TorsionOrderStruct>::iterator l = torsionOrderStruct.begin(); l != torsionOrderStruct.end(); ++l) { |
119 | > | Torsion* torsion = l->torsion; |
120 | > | std::cout << "atom1=" <<torsion->getAtomA()->getGlobalIndex() << ",atom2 = "<< torsion->getAtomB()->getGlobalIndex() << ",atom3="<<torsion->getAtomC()->getGlobalIndex() << ",atom4="<<torsion->getAtomD()->getGlobalIndex()<< " "; |
121 | > | 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; |
122 | > | } |
123 | > | */ |
124 | } | |
125 | ||
126 | void ForceManager::preCalculation() { | |
# | Line 105 | Line 157 | namespace oopse { | |
157 | Molecule::BondIterator bondIter;; | |
158 | Molecule::BendIterator bendIter; | |
159 | Molecule::TorsionIterator torsionIter; | |
160 | + | double bondPotential = 0.0; |
161 | + | double bendPotential = 0.0; |
162 | + | double torsionPotential = 0.0; |
163 | ||
164 | //calculate short range interactions | |
165 | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { | |
166 | ||
167 | //change the positions of atoms which belong to the rigidbodies | |
168 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
169 | < | rb->updateAtoms(); |
169 | > | rb->updateAtoms(); |
170 | } | |
171 | ||
172 | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
173 | < | bond->calcForce(); |
173 | > | bond->calcForce(); |
174 | > | bondPotential += bond->getPotential(); |
175 | } | |
176 | + | |
177 | ||
178 | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { | |
122 | – | bend->calcForce(); |
123 | – | } |
179 | ||
180 | < | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
181 | < | torsion->calcForce(); |
180 | > | double angle; |
181 | > | bend->calcForce(angle); |
182 | > | double currBendPot = bend->getPotential(); |
183 | > | bendPotential += bend->getPotential(); |
184 | > | std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend); |
185 | > | if (i == bendDataSets.end()) { |
186 | > | BendDataSet dataSet; |
187 | > | dataSet.prev.angle = dataSet.curr.angle = angle; |
188 | > | dataSet.prev.potential = dataSet.curr.potential = currBendPot; |
189 | > | dataSet.deltaV = 0.0; |
190 | > | bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet)); |
191 | > | }else { |
192 | > | i->second.prev.angle = i->second.curr.angle; |
193 | > | i->second.prev.potential = i->second.curr.potential; |
194 | > | i->second.curr.angle = angle; |
195 | > | i->second.curr.potential = currBendPot; |
196 | > | i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential); |
197 | > | } |
198 | } | |
199 | ||
129 | – | } |
130 | – | |
131 | – | |
132 | – | double bondPotential = 0.0; |
133 | – | double bendPotential = 0.0; |
134 | – | double torsionPotential = 0.0; |
135 | – | |
136 | – | for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
137 | – | |
138 | – | for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
139 | – | bondPotential += bond->getPotential(); |
140 | – | } |
141 | – | |
142 | – | for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
143 | – | bendPotential += bend->getPotential(); |
144 | – | } |
145 | – | |
200 | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
201 | + | double angle; |
202 | + | torsion->calcForce(angle); |
203 | + | double currTorsionPot = torsion->getPotential(); |
204 | torsionPotential += torsion->getPotential(); | |
205 | + | std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion); |
206 | + | if (i == torsionDataSets.end()) { |
207 | + | TorsionDataSet dataSet; |
208 | + | dataSet.prev.angle = dataSet.curr.angle = angle; |
209 | + | dataSet.prev.potential = dataSet.curr.potential = currTorsionPot; |
210 | + | dataSet.deltaV = 0.0; |
211 | + | torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet)); |
212 | + | }else { |
213 | + | i->second.prev.angle = i->second.curr.angle; |
214 | + | i->second.prev.potential = i->second.curr.potential; |
215 | + | i->second.curr.angle = angle; |
216 | + | i->second.curr.potential = currTorsionPot; |
217 | + | i->second.deltaV = fabs(i->second.curr.potential - i->second.prev.potential); |
218 | + | } |
219 | } | |
220 | ||
221 | < | } |
222 | < | |
221 | > | } |
222 | > | |
223 | double shortRangePotential = bondPotential + bendPotential + torsionPotential; | |
224 | Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); | |
225 | curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; | |
# | Line 215 | Line 286 | namespace oopse { | |
286 | longRangePotential[i]=0.0; //Initialize array | |
287 | } | |
288 | ||
218 | – | |
219 | – | |
289 | doForceLoop( pos, | |
290 | rc, | |
291 | A, | |
# | Line 241 | Line 310 | namespace oopse { | |
310 | ||
311 | //store the tau and long range potential | |
312 | curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; | |
244 | – | // curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; |
313 | curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; | |
314 | curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; | |
315 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |