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

Comparing:
trunk/src/brains/ForceManager.cpp (file contents), Revision 1448 by gezelter, Thu Jun 17 14:58:49 2010 UTC vs.
branches/development/src/brains/ForceManager.cpp (file contents), Revision 1540 by gezelter, Mon Jan 17 21:34:36 2011 UTC

# Line 57 | Line 57
57   #include "primitives/Bend.hpp"
58   #include "primitives/Torsion.hpp"
59   #include "primitives/Inversion.hpp"
60 – namespace OpenMD {
60  
61 <  void ForceManager::calcForces(bool needPotential, bool needStress) {
61 > namespace OpenMD {
62 >  
63 >  ForceManager::ForceManager(SimInfo * info) : info_(info),
64 >                                               NBforcesInitialized_(false) {
65 >  }
66 >
67 >  void ForceManager::calcForces() {
68      
69 +
70      if (!info_->isFortranInitialized()) {
71        info_->update();
72 +      nbiMan_->setSimInfo(info_);
73 +      nbiMan_->initialize();    
74 +      info_->setupFortran();
75      }
76      
77      preCalculation();
78      
79      calcShortRangeInteraction();
80  
81 <    calcLongRangeInteraction(needPotential, needStress);
81 >    calcLongRangeInteraction();
82  
83 <    postCalculation(needStress);
83 >    postCalculation();
84      
85    }
86    
# Line 82 | Line 91 | namespace OpenMD {
91      Atom* atom;
92      Molecule::RigidBodyIterator rbIter;
93      RigidBody* rb;
94 +    Molecule::CutoffGroupIterator ci;
95 +    CutoffGroup* cg;
96      
97      // forces are zeroed here, before any are accumulated.
98      // NOTE: do not rezero the forces in Fortran.
# Line 97 | Line 108 | namespace OpenMD {
108             rb = mol->nextRigidBody(rbIter)) {
109          rb->zeroForcesAndTorques();
110        }        
111 <          
111 >
112 >      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
113 >        std::cerr << "should not see me \n";
114 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
115 >            cg = mol->nextCutoffGroup(ci)) {
116 >          //calculate the center of mass of cutoff group
117 >          cg->updateCOM();
118 >        }
119 >      }      
120      }
121 <    
121 >  
122      // Zero out the stress tensor
123      tau *= 0.0;
124      
# Line 223 | Line 242 | namespace OpenMD {
242      
243    }
244    
245 <  void ForceManager::calcLongRangeInteraction(bool needPotential,
227 <                                              bool needStress) {
245 >  void ForceManager::calcLongRangeInteraction() {
246      Snapshot* curSnapshot;
247      DataStorage* config;
248 +    DataStorage* cgConfig;
249      RealType* frc;
250      RealType* pos;
251      RealType* trq;
# Line 240 | Line 259 | namespace OpenMD {
259      
260      //get array pointers
261      config = &(curSnapshot->atomData);
262 +    cgConfig = &(curSnapshot->cgData);
263      frc = config->getArrayPointer(DataStorage::dslForce);
264      pos = config->getArrayPointer(DataStorage::dslPosition);
265      trq = config->getArrayPointer(DataStorage::dslTorque);
# Line 247 | Line 267 | namespace OpenMD {
267      electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
268      particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
269  
270 <    //calculate the center of mass of cutoff group
271 <    SimInfo::MoleculeIterator mi;
272 <    Molecule* mol;
253 <    Molecule::CutoffGroupIterator ci;
254 <    CutoffGroup* cg;
255 <    Vector3d com;
256 <    std::vector<Vector3d> rcGroup;
257 <    
258 <    if(info_->getNCutoffGroups() > 0){
259 <      
260 <      for (mol = info_->beginMolecule(mi); mol != NULL;
261 <           mol = info_->nextMolecule(mi)) {
262 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
263 <            cg = mol->nextCutoffGroup(ci)) {
264 <          cg->getCOM(com);
265 <          rcGroup.push_back(com);
266 <        }
267 <      }// end for (mol)
268 <      
269 <      rc = rcGroup[0].getArrayPointer();
270 >    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
271 >      std::cerr << "should not see me \n";
272 >      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
273      } else {
274        // center of mass of the group is the same as position of the atom  
275        // if cutoff group does not exist
# Line 276 | Line 279 | namespace OpenMD {
279      //initialize data before passing to fortran
280      RealType longRangePotential[LR_POT_TYPES];
281      RealType lrPot = 0.0;
279 –    Vector3d totalDipole;
280 –    short int passedCalcPot = needPotential;
281 –    short int passedCalcStress = needStress;
282      int isError = 0;
283  
284      for (int i=0; i<LR_POT_TYPES;i++){
# Line 294 | Line 294 | namespace OpenMD {
294                  tau.getArrayPointer(),
295                  longRangePotential,
296                  particlePot,
297 –                &passedCalcPot,
298 –                &passedCalcStress,
297                  &isError );
298      
299      if( isError ){
# Line 307 | Line 305 | namespace OpenMD {
305      for (int i=0; i<LR_POT_TYPES;i++){
306        lrPot += longRangePotential[i]; //Quick hack
307      }
308 <    
311 <    // grab the simulation box dipole moment if specified
312 <    if (info_->getCalcBoxDipole()){
313 <      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
314 <      
315 <      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
316 <      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
317 <      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
318 <    }
319 <    
308 >        
309      //store the tau and long range potential    
310      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
311      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
# Line 324 | Line 313 | namespace OpenMD {
313    }
314  
315    
316 <  void ForceManager::postCalculation(bool needStress) {
316 >  void ForceManager::postCalculation() {
317      SimInfo::MoleculeIterator mi;
318      Molecule* mol;
319      Molecule::RigidBodyIterator rbIter;
# Line 337 | Line 326 | namespace OpenMD {
326           mol = info_->nextMolecule(mi)) {
327        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
328             rb = mol->nextRigidBody(rbIter)) {
329 <        if (needStress) {          
330 <          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
342 <          tau += rbTau;
343 <        } else{
344 <          rb->calcForcesAndTorques();
345 <        }
329 >        Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
330 >        tau += rbTau;
331        }
332      }
333 <
349 <    if (needStress) {
333 >    
334   #ifdef IS_MPI
335 <      Mat3x3d tmpTau(tau);
336 <      MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
337 <                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
335 >    Mat3x3d tmpTau(tau);
336 >    MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
337 >                  9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
338   #endif
339 <      curSnapshot->statData.setTau(tau);
356 <    }
339 >    curSnapshot->statData.setTau(tau);
340    }
341  
342   } //end namespace OpenMD

Diff Legend

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