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 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2469 by tim, Fri Dec 2 15:38:03 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 39 | Line 39
39   * such damages.
40   */
41  
42 < /**
43 <  * @file ForceManager.cpp
44 <  * @author tlin
45 <  * @date 11/09/2004
46 <  * @time 10:39am
47 <  * @version 1.0
48 <  */
42 > /**
43 > * @file ForceManager.cpp
44 > * @author tlin
45 > * @date 11/09/2004
46 > * @time 10:39am
47 > * @version 1.0
48 > */
49  
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 < void ForceManager::calcForces(bool needPotential, bool needStress) {
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()) {
81 <        info_->update();
81 >      info_->update();
82      }
83  
84      preCalculation();
# Line 66 | Line 88 | void ForceManager::calcForces(bool needPotential, bool
88      calcLongRangeInteraction(needPotential, needStress);
89  
90      postCalculation();
69 –        
70 – }
91  
92 < void ForceManager::preCalculation() {
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() {
127      SimInfo::MoleculeIterator mi;
128      Molecule* mol;
129      Molecule::AtomIterator ai;
# Line 80 | Line 134 | void ForceManager::preCalculation() {
134      // forces are zeroed here, before any are accumulated.
135      // NOTE: do not rezero the forces in Fortran.
136      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
137 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
138 <            atom->zeroForcesAndTorques();
139 <        }
137 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
138 >        atom->zeroForcesAndTorques();
139 >      }
140          
141 <        //change the positions of atoms which belong to the rigidbodies
142 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
143 <            rb->zeroForcesAndTorques();
144 <        }        
141 >      //change the positions of atoms which belong to the rigidbodies
142 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
143 >        rb->zeroForcesAndTorques();
144 >      }        
145      }
146      
147 < }
147 >  }
148  
149 < void ForceManager::calcShortRangeInteraction() {
149 >  void ForceManager::calcShortRangeInteraction() {
150      Molecule* mol;
151      RigidBody* rb;
152      Bond* bond;
# Line 103 | Line 157 | void ForceManager::calcShortRangeInteraction() {
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();
170 <        }
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();
170 >      }
171  
172 <        for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
173 <            bond->calcForce();
174 <        }
172 >      for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
173 >        bond->calcForce();
174 >        bondPotential += bond->getPotential();
175 >      }
176  
119 –        for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
120 –            bend->calcForce();
121 –        }
177  
178 <        for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
124 <            torsion->calcForce();
125 <        }
178 >      for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
179  
180 <    }
181 <    
182 <    double  shortRangePotential = 0.0;
183 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
184 <        shortRangePotential += mol->getPotential();
185 <    }
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  
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 +    
223 +    double  shortRangePotential = bondPotential + bendPotential + torsionPotential;    
224      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
225      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
226 < }
226 >    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
227 >    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
228 >    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
229 >    
230 >  }
231  
232 < void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
232 >  void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
233      Snapshot* curSnapshot;
234      DataStorage* config;
235      double* frc;
# Line 163 | Line 257 | void ForceManager::calcLongRangeInteraction(bool needP
257      CutoffGroup* cg;
258      Vector3d com;
259      std::vector<Vector3d> rcGroup;
166 –    
167 –    if(info_->getNCutoffGroups() > 0){
260  
261 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
261 >    if(info_->getNCutoffGroups() > 0){
262 >
263 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
264          for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
265 <            cg->getCOM(com);
266 <            rcGroup.push_back(com);
265 >          cg->getCOM(com);
266 >          rcGroup.push_back(com);
267          }
268 <    }// end for (mol)
268 >      }// end for (mol)
269        
270 <        rc = rcGroup[0].getArrayPointer();
270 >      rc = rcGroup[0].getArrayPointer();
271      } else {
272 <        // center of mass of the group is the same as position of the atom  if cutoff group does not exist
273 <        rc = pos;
272 >      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
273 >      rc = pos;
274      }
275    
276      //initialize data before passing to fortran
277 <    double longRangePotential = 0.0;
277 >    double longRangePotential[LR_POT_TYPES];
278 >    double lrPot = 0.0;
279 >    
280      Mat3x3d tau;
281      short int passedCalcPot = needPotential;
282      short int passedCalcStress = needStress;
283      int isError = 0;
284  
285 +    for (int i=0; i<LR_POT_TYPES;i++){
286 +      longRangePotential[i]=0.0; //Initialize array
287 +    }
288 +
289      doForceLoop( pos,
290 <            rc,
291 <            A,
292 <            electroFrame,
293 <            frc,
294 <            trq,
295 <            tau.getArrayPointer(),
296 <            &longRangePotential,
297 <            &passedCalcPot,
298 <            &passedCalcStress,
299 <            &isError );
290 >                 rc,
291 >                 A,
292 >                 electroFrame,
293 >                 frc,
294 >                 trq,
295 >                 tau.getArrayPointer(),
296 >                 longRangePotential,
297 >                 &passedCalcPot,
298 >                 &passedCalcStress,
299 >                 &isError );
300  
301      if( isError ){
302 <        sprintf( painCave.errMsg,
303 <             "Error returned from the fortran force calculation.\n" );
304 <        painCave.isFatal = 1;
305 <        simError();
302 >      sprintf( painCave.errMsg,
303 >               "Error returned from the fortran force calculation.\n" );
304 >      painCave.isFatal = 1;
305 >      simError();
306      }
307 +    for (int i=0; i<LR_POT_TYPES;i++){
308 +      lrPot += longRangePotential[i]; //Quick hack
309 +    }
310  
311      //store the tau and long range potential    
312 <    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
312 >    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
313 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
314 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
315 >
316      curSnapshot->statData.setTau(tau);
317 < }
317 >  }
318  
319  
320 < void ForceManager::postCalculation() {
320 >  void ForceManager::postCalculation() {
321      SimInfo::MoleculeIterator mi;
322      Molecule* mol;
323      Molecule::RigidBodyIterator rbIter;
# Line 219 | Line 325 | void ForceManager::postCalculation() {
325      
326      // collect the atomic forces onto rigid bodies
327      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
328 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
329 <            rb->calcForcesAndTorques();
330 <        }
328 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
329 >        rb->calcForcesAndTorques();
330 >      }
331      }
332  
333 < }
333 >  }
334  
335   } //end namespace oopse

Diff Legend

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