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

Comparing trunk/src/brains/ForceManager.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 749 by tim, Wed Nov 16 23:10:02 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 >  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()) {
80 <        info_->update();
80 >      info_->update();
81      }
82  
83      preCalculation();
# Line 66 | Line 87 | void ForceManager::calcForces(bool needPotential, bool
87      calcLongRangeInteraction(needPotential, needStress);
88  
89      postCalculation();
69 –        
70 – }
90  
91 < void ForceManager::preCalculation() {
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() {
125      SimInfo::MoleculeIterator mi;
126      Molecule* mol;
127      Molecule::AtomIterator ai;
# Line 80 | Line 132 | void ForceManager::preCalculation() {
132      // forces are zeroed here, before any are accumulated.
133      // NOTE: do not rezero the forces in Fortran.
134      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
135 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
136 <            atom->zeroForcesAndTorques();
137 <        }
135 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
136 >        atom->zeroForcesAndTorques();
137 >      }
138          
139 <        //change the positions of atoms which belong to the rigidbodies
140 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
141 <            rb->zeroForcesAndTorques();
142 <        }        
139 >      //change the positions of atoms which belong to the rigidbodies
140 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
141 >        rb->zeroForcesAndTorques();
142 >      }        
143      }
144      
145 < }
145 >  }
146  
147 < void ForceManager::calcShortRangeInteraction() {
147 >  void ForceManager::calcShortRangeInteraction() {
148      Molecule* mol;
149      RigidBody* rb;
150      Bond* bond;
# Line 103 | Line 155 | void ForceManager::calcShortRangeInteraction() {
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();
168 <        }
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();
168 >      }
169  
170 <        for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
171 <            bond->calcForce();
172 <        }
170 >      for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
171 >        bond->calcForce();
172 >        bondPotential += bond->getPotential();
173 >      }
174  
175 <        for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
176 <            bend->calcForce();
177 <        }
175 >      //int i =0;
176 >      for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
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();
200 <        }
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      
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 < }
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) {
230 >  void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
231      Snapshot* curSnapshot;
232      DataStorage* config;
233      double* frc;
# Line 163 | Line 255 | void ForceManager::calcLongRangeInteraction(bool needP
255      CutoffGroup* cg;
256      Vector3d com;
257      std::vector<Vector3d> rcGroup;
166 –    
167 –    if(info_->getNCutoffGroups() > 0){
258  
259 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
259 >    if(info_->getNCutoffGroups() > 0){
260 >
261 >      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
262          for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
263 <            cg->getCOM(com);
264 <            rcGroup.push_back(com);
263 >          cg->getCOM(com);
264 >          rcGroup.push_back(com);
265          }
266 <    }// end for (mol)
266 >      }// end for (mol)
267        
268 <        rc = rcGroup[0].getArrayPointer();
268 >      rc = rcGroup[0].getArrayPointer();
269      } else {
270 <        // center of mass of the group is the same as position of the atom  if cutoff group does not exist
271 <        rc = pos;
270 >      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
271 >      rc = pos;
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,
290 <            electroFrame,
291 <            frc,
292 <            trq,
293 <            tau.getArrayPointer(),
294 <            &longRangePotential,
295 <            &passedCalcPot,
296 <            &passedCalcStress,
297 <            &isError );
288 >                 rc,
289 >                 A,
290 >                 electroFrame,
291 >                 frc,
292 >                 trq,
293 >                 tau.getArrayPointer(),
294 >                 longRangePotential,
295 >                 &passedCalcPot,
296 >                 &passedCalcStress,
297 >                 &isError );
298  
299      if( isError ){
300 <        sprintf( painCave.errMsg,
301 <             "Error returned from the fortran force calculation.\n" );
302 <        painCave.isFatal = 1;
303 <        simError();
300 >      sprintf( painCave.errMsg,
301 >               "Error returned from the fortran force calculation.\n" );
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 < }
315 >  }
316  
317  
318 < void ForceManager::postCalculation() {
318 >  void ForceManager::postCalculation() {
319      SimInfo::MoleculeIterator mi;
320      Molecule* mol;
321      Molecule::RigidBodyIterator rbIter;
# Line 219 | Line 323 | void ForceManager::postCalculation() {
323      
324      // collect the atomic forces onto rigid bodies
325      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
326 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
327 <            rb->calcForcesAndTorques();
328 <        }
326 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
327 >        rb->calcForcesAndTorques();
328 >      }
329      }
330  
331 < }
331 >  }
332  
333   } //end namespace oopse

Diff Legend

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