ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/brains/ForceManager.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/brains/ForceManager.cpp (file contents):
Revision 2364 by tim, Thu Oct 13 22:26:47 2005 UTC vs.
Revision 2475 by gezelter, Fri Dec 2 20:10:49 2005 UTC

# Line 53 | Line 53 | namespace oopse {
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 +    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 +
176        for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
122 –        bend->calcForce();
123 –      }
177  
178 <      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
179 <        torsion->calcForce();
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  
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 –
198        for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
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 <
219 >    }
220 >    
221      double  shortRangePotential = bondPotential + bendPotential + torsionPotential;    
222      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
223      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
# Line 215 | Line 284 | namespace oopse {
284        longRangePotential[i]=0.0; //Initialize array
285      }
286  
218 –
219 –
287      doForceLoop( pos,
288                   rc,
289                   A,
# Line 241 | Line 308 | namespace oopse {
308  
309      //store the tau and long range potential    
310      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
311 <    //  curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
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  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines