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 749 by tim, Wed Nov 16 23:10:02 2005 UTC vs.
Revision 1245 by chuckv, Tue May 27 16:39:06 2008 UTC

# Line 53 | Line 53
53   #define __C
54   #include "UseTheForce/DarkSide/fInteractionMap.h"
55   #include "utils/simError.h"
56 + #include "primitives/Bond.hpp"
57   #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 –  
60    void ForceManager::calcForces(bool needPotential, bool needStress) {
61 <
61 >    
62      if (!info_->isFortranInitialized()) {
63        info_->update();
64      }
65 <
65 >    
66      preCalculation();
67      
68      calcShortRangeInteraction();
69  
70      calcLongRangeInteraction(needPotential, needStress);
71  
72 <    postCalculation();
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 <    }
72 >    postCalculation(needStress);
73      
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 –    
74    }
75 <
75 >  
76    void ForceManager::preCalculation() {
77      SimInfo::MoleculeIterator mi;
78      Molecule* mol;
# Line 131 | Line 83 | namespace oopse {
83      
84      // forces are zeroed here, before any are accumulated.
85      // NOTE: do not rezero the forces in Fortran.
86 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
86 >    
87 >    for (mol = info_->beginMolecule(mi); mol != NULL;
88 >         mol = info_->nextMolecule(mi)) {
89        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
90          atom->zeroForcesAndTorques();
91        }
92 <        
92 >          
93        //change the positions of atoms which belong to the rigidbodies
94 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
94 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
95 >           rb = mol->nextRigidBody(rbIter)) {
96          rb->zeroForcesAndTorques();
97        }        
98 +          
99      }
100      
101 +    // Zero out the stress tensor
102 +    tau *= 0.0;
103 +    
104    }
105 <
105 >  
106    void ForceManager::calcShortRangeInteraction() {
107      Molecule* mol;
108      RigidBody* rb;
# Line 155 | Line 114 | namespace oopse {
114      Molecule::BondIterator bondIter;;
115      Molecule::BendIterator  bendIter;
116      Molecule::TorsionIterator  torsionIter;
117 <    double bondPotential = 0.0;
118 <    double bendPotential = 0.0;
119 <    double torsionPotential = 0.0;
117 >    RealType bondPotential = 0.0;
118 >    RealType bendPotential = 0.0;
119 >    RealType torsionPotential = 0.0;
120  
121      //calculate short range interactions    
122 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
122 >    for (mol = info_->beginMolecule(mi); mol != NULL;
123 >         mol = info_->nextMolecule(mi)) {
124  
125        //change the positions of atoms which belong to the rigidbodies
126 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
127 <          rb->updateAtoms();
126 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
127 >           rb = mol->nextRigidBody(rbIter)) {
128 >        rb->updateAtoms();
129        }
130  
131 <      for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
131 >      for (bond = mol->beginBond(bondIter); bond != NULL;
132 >           bond = mol->nextBond(bondIter)) {
133          bond->calcForce();
134          bondPotential += bond->getPotential();
135        }
136  
137 <      //int i =0;
138 <      for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
139 <          //std::cout << i++ << "\t";
140 <          double angle;
141 <            bend->calcForce(angle);
142 <          double currBendPot = bend->getPotential();          
143 <            bendPotential += bend->getPotential();
144 <          std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
145 <          if (i == bendDataSets.end()) {
146 <            BendDataSet dataSet;
147 <            dataSet.prev.angle = dataSet.curr.angle = angle;
148 <            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
149 <            dataSet.deltaV = 0.0;
150 <            bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
151 <          }else {
152 <            i->second.prev.angle = i->second.curr.angle;
153 <            i->second.prev.potential = i->second.curr.potential;
154 <            i->second.curr.angle = angle;
155 <            i->second.curr.potential = currBendPot;
156 <            i->second.deltaV =  fabs(i->second.curr.potential -  i->second.prev.potential);
157 <          }
158 <      }
197 <
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 <          }      
137 >      for (bend = mol->beginBend(bendIter); bend != NULL;
138 >           bend = mol->nextBend(bendIter)) {
139 >        
140 >        RealType angle;
141 >        bend->calcForce(angle);
142 >        RealType currBendPot = bend->getPotential();          
143 >        bendPotential += bend->getPotential();
144 >        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
145 >        if (i == bendDataSets.end()) {
146 >          BendDataSet dataSet;
147 >          dataSet.prev.angle = dataSet.curr.angle = angle;
148 >          dataSet.prev.potential = dataSet.curr.potential = currBendPot;
149 >          dataSet.deltaV = 0.0;
150 >          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
151 >        }else {
152 >          i->second.prev.angle = i->second.curr.angle;
153 >          i->second.prev.potential = i->second.curr.potential;
154 >          i->second.curr.angle = angle;
155 >          i->second.curr.potential = currBendPot;
156 >          i->second.deltaV =  fabs(i->second.curr.potential -  
157 >                                   i->second.prev.potential);
158 >        }
159        }
160 <
160 >      
161 >      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
162 >           torsion = mol->nextTorsion(torsionIter)) {
163 >        RealType angle;
164 >        torsion->calcForce(angle);
165 >        RealType currTorsionPot = torsion->getPotential();
166 >        torsionPotential += torsion->getPotential();
167 >        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
168 >        if (i == torsionDataSets.end()) {
169 >          TorsionDataSet dataSet;
170 >          dataSet.prev.angle = dataSet.curr.angle = angle;
171 >          dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
172 >          dataSet.deltaV = 0.0;
173 >          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
174 >        }else {
175 >          i->second.prev.angle = i->second.curr.angle;
176 >          i->second.prev.potential = i->second.curr.potential;
177 >          i->second.curr.angle = angle;
178 >          i->second.curr.potential = currTorsionPot;
179 >          i->second.deltaV =  fabs(i->second.curr.potential -  
180 >                                   i->second.prev.potential);
181 >        }      
182 >      }      
183      }
184      
185 <    double  shortRangePotential = bondPotential + bendPotential + torsionPotential;    
185 >    RealType  shortRangePotential = bondPotential + bendPotential +
186 >      torsionPotential;    
187      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
188      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
189      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
# Line 226 | Line 191 | namespace oopse {
191      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
192      
193    }
194 <
195 <  void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
194 >  
195 >  void ForceManager::calcLongRangeInteraction(bool needPotential,
196 >                                              bool needStress) {
197      Snapshot* curSnapshot;
198      DataStorage* config;
199 <    double* frc;
200 <    double* pos;
201 <    double* trq;
202 <    double* A;
203 <    double* electroFrame;
204 <    double* rc;
199 >    RealType* frc;
200 >    RealType* pos;
201 >    RealType* trq;
202 >    RealType* A;
203 >    RealType* electroFrame;
204 >    RealType* rc;
205 >    RealType* particlePot;
206      
207      //get current snapshot from SimInfo
208      curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
209 <
209 >    
210      //get array pointers
211      config = &(curSnapshot->atomData);
212      frc = config->getArrayPointer(DataStorage::dslForce);
# Line 247 | Line 214 | namespace oopse {
214      trq = config->getArrayPointer(DataStorage::dslTorque);
215      A   = config->getArrayPointer(DataStorage::dslAmat);
216      electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
217 +    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
218  
219      //calculate the center of mass of cutoff group
220      SimInfo::MoleculeIterator mi;
# Line 255 | Line 223 | namespace oopse {
223      CutoffGroup* cg;
224      Vector3d com;
225      std::vector<Vector3d> rcGroup;
226 <
226 >    
227      if(info_->getNCutoffGroups() > 0){
228 <
229 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
230 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
228 >      
229 >      for (mol = info_->beginMolecule(mi); mol != NULL;
230 >           mol = info_->nextMolecule(mi)) {
231 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
232 >            cg = mol->nextCutoffGroup(ci)) {
233            cg->getCOM(com);
234            rcGroup.push_back(com);
235          }
# Line 267 | Line 237 | namespace oopse {
237        
238        rc = rcGroup[0].getArrayPointer();
239      } else {
240 <      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
240 >      // center of mass of the group is the same as position of the atom  
241 >      // if cutoff group does not exist
242        rc = pos;
243      }
273 –  
274 –    //initialize data before passing to fortran
275 –    double longRangePotential[LR_POT_TYPES];
276 –    double lrPot = 0.0;
244      
245 <    Mat3x3d tau;
245 >    //initialize data before passing to fortran
246 >    RealType longRangePotential[LR_POT_TYPES];
247 >    RealType lrPot = 0.0;
248 >    Vector3d totalDipole;
249      short int passedCalcPot = needPotential;
250      short int passedCalcStress = needStress;
251      int isError = 0;
# Line 283 | Line 253 | namespace oopse {
253      for (int i=0; i<LR_POT_TYPES;i++){
254        longRangePotential[i]=0.0; //Initialize array
255      }
256 <
257 <    doForceLoop( pos,
258 <                 rc,
259 <                 A,
260 <                 electroFrame,
261 <                 frc,
262 <                 trq,
263 <                 tau.getArrayPointer(),
264 <                 longRangePotential,
265 <                 &passedCalcPot,
266 <                 &passedCalcStress,
267 <                 &isError );
268 <
256 >    
257 >    doForceLoop(pos,
258 >                rc,
259 >                A,
260 >                electroFrame,
261 >                frc,
262 >                trq,
263 >                tau.getArrayPointer(),
264 >                longRangePotential,
265 >                particlePot,
266 >                &passedCalcPot,
267 >                &passedCalcStress,
268 >                &isError );
269 >    
270      if( isError ){
271        sprintf( painCave.errMsg,
272                 "Error returned from the fortran force calculation.\n" );
# Line 305 | Line 276 | namespace oopse {
276      for (int i=0; i<LR_POT_TYPES;i++){
277        lrPot += longRangePotential[i]; //Quick hack
278      }
279 <
279 >    
280 >    // grab the simulation box dipole moment if specified
281 >    if (info_->getCalcBoxDipole()){
282 >      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
283 >      
284 >      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
285 >      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
286 >      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
287 >    }
288 >    
289      //store the tau and long range potential    
290      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
291      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
292      curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
313 –
314 –    curSnapshot->statData.setTau(tau);
293    }
294  
295 <
296 <  void ForceManager::postCalculation() {
295 >  
296 >  void ForceManager::postCalculation(bool needStress) {
297      SimInfo::MoleculeIterator mi;
298      Molecule* mol;
299      Molecule::RigidBodyIterator rbIter;
300      RigidBody* rb;
301 +    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
302      
303      // collect the atomic forces onto rigid bodies
304 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
305 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
306 <        rb->calcForcesAndTorques();
304 >    
305 >    for (mol = info_->beginMolecule(mi); mol != NULL;
306 >         mol = info_->nextMolecule(mi)) {
307 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
308 >           rb = mol->nextRigidBody(rbIter)) {
309 >        if (needStress) {          
310 >          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
311 >          tau += rbTau;
312 >        } else{
313 >          rb->calcForcesAndTorques();
314 >        }
315        }
316      }
317  
318 +    if (needStress) {
319 + #ifdef IS_MPI
320 +      Mat3x3d tmpTau(tau);
321 +      MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
322 +                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
323 + #endif
324 +      curSnapshot->statData.setTau(tau);
325 +    }
326    }
327  
328   } //end namespace oopse

Diff Legend

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