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

Comparing branches/development/src/brains/ForceManager.cpp (file contents):
Revision 1540 by gezelter, Mon Jan 17 21:34:36 2011 UTC vs.
Revision 1549 by gezelter, Wed Apr 27 18:38:15 2011 UTC

# Line 57 | Line 57
57   #include "primitives/Bend.hpp"
58   #include "primitives/Torsion.hpp"
59   #include "primitives/Inversion.hpp"
60 + #include "parallel/ForceMatrixDecomposition.hpp"
61 + //#include "parallel/ForceSerialDecomposition.hpp"
62  
63 + using namespace std;
64   namespace OpenMD {
65    
66 <  ForceManager::ForceManager(SimInfo * info) : info_(info),
67 <                                               NBforcesInitialized_(false) {
66 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
67 >
68 > #ifdef IS_MPI
69 >    fDecomp_ = new ForceMatrixDecomposition(info_);
70 > #else
71 >    // fDecomp_ = new ForceSerialDecomposition(info);
72 > #endif
73    }
74 <
74 >  
75    void ForceManager::calcForces() {
76      
69 –
77      if (!info_->isFortranInitialized()) {
78        info_->update();
79 <      nbiMan_->setSimInfo(info_);
80 <      nbiMan_->initialize();    
79 >      interactionMan_->setSimInfo(info_);
80 >      interactionMan_->initialize();
81 >      swfun_ = interactionMan_->getSwitchingFunction();
82 >      fDecomp_->distributeInitialData();
83        info_->setupFortran();
84      }
85      
86 <    preCalculation();
87 <    
88 <    calcShortRangeInteraction();
80 <
81 <    calcLongRangeInteraction();
82 <
86 >    preCalculation();  
87 >    shortRangeInteractions();
88 >    longRangeInteractions();
89      postCalculation();
90      
91    }
# Line 95 | Line 101 | namespace OpenMD {
101      CutoffGroup* cg;
102      
103      // forces are zeroed here, before any are accumulated.
98 –    // NOTE: do not rezero the forces in Fortran.
104      
105      for (mol = info_->beginMolecule(mi); mol != NULL;
106           mol = info_->nextMolecule(mi)) {
# Line 110 | Line 115 | namespace OpenMD {
115        }        
116  
117        if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
113 –        std::cerr << "should not see me \n";
118          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
119              cg = mol->nextCutoffGroup(ci)) {
120            //calculate the center of mass of cutoff group
# Line 124 | Line 128 | namespace OpenMD {
128      
129    }
130    
131 <  void ForceManager::calcShortRangeInteraction() {
131 >  void ForceManager::shortRangeInteractions() {
132      Molecule* mol;
133      RigidBody* rb;
134      Bond* bond;
# Line 166 | Line 170 | namespace OpenMD {
170          RealType currBendPot = bend->getPotential();          
171          
172          bendPotential += bend->getPotential();
173 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
173 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
174          if (i == bendDataSets.end()) {
175            BendDataSet dataSet;
176            dataSet.prev.angle = dataSet.curr.angle = angle;
177            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
178            dataSet.deltaV = 0.0;
179 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
179 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
180          }else {
181            i->second.prev.angle = i->second.curr.angle;
182            i->second.prev.potential = i->second.curr.potential;
# Line 189 | Line 193 | namespace OpenMD {
193          torsion->calcForce(angle);
194          RealType currTorsionPot = torsion->getPotential();
195          torsionPotential += torsion->getPotential();
196 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
196 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
197          if (i == torsionDataSets.end()) {
198            TorsionDataSet dataSet;
199            dataSet.prev.angle = dataSet.curr.angle = angle;
200            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
201            dataSet.deltaV = 0.0;
202 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
202 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
203          }else {
204            i->second.prev.angle = i->second.curr.angle;
205            i->second.prev.potential = i->second.curr.potential;
# Line 205 | Line 209 | namespace OpenMD {
209                                     i->second.prev.potential);
210          }      
211        }      
212 <
212 >      
213        for (inversion = mol->beginInversion(inversionIter);
214             inversion != NULL;
215             inversion = mol->nextInversion(inversionIter)) {
# Line 213 | Line 217 | namespace OpenMD {
217          inversion->calcForce(angle);
218          RealType currInversionPot = inversion->getPotential();
219          inversionPotential += inversion->getPotential();
220 <        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
220 >        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
221          if (i == inversionDataSets.end()) {
222            InversionDataSet dataSet;
223            dataSet.prev.angle = dataSet.curr.angle = angle;
224            dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
225            dataSet.deltaV = 0.0;
226 <          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
226 >          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
227          }else {
228            i->second.prev.angle = i->second.curr.angle;
229            i->second.prev.potential = i->second.curr.potential;
# Line 238 | Line 242 | namespace OpenMD {
242      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
243      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
244      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
245 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
242 <    
245 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
246    }
247    
248 <  void ForceManager::calcLongRangeInteraction() {
246 <    Snapshot* curSnapshot;
247 <    DataStorage* config;
248 <    DataStorage* cgConfig;
249 <    RealType* frc;
250 <    RealType* pos;
251 <    RealType* trq;
252 <    RealType* A;
253 <    RealType* electroFrame;
254 <    RealType* rc;
255 <    RealType* particlePot;
256 <    
257 <    //get current snapshot from SimInfo
258 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
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);
266 <    A   = config->getArrayPointer(DataStorage::dslAmat);
267 <    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
268 <    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
248 >  void ForceManager::longRangeInteractions() {
249  
250 +    // some of this initial stuff will go away:
251 +    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
252 +    DataStorage* config = &(curSnapshot->atomData);
253 +    DataStorage* cgConfig = &(curSnapshot->cgData);
254 +    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
255 +    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
256 +    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
257 +    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
258 +    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
259 +    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
260 +    RealType* rc;    
261 +
262      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
271 –      std::cerr << "should not see me \n";
263        rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
264      } else {
265        // center of mass of the group is the same as position of the atom  
# Line 284 | Line 275 | namespace OpenMD {
275      for (int i=0; i<LR_POT_TYPES;i++){
276        longRangePotential[i]=0.0; //Initialize array
277      }
278 <    
279 <    doForceLoop(pos,
280 <                rc,
281 <                A,
282 <                electroFrame,
283 <                frc,
284 <                trq,
285 <                tau.getArrayPointer(),
286 <                longRangePotential,
287 <                particlePot,
288 <                &isError );
289 <    
290 <    if( isError ){
291 <      sprintf( painCave.errMsg,
292 <               "Error returned from the fortran force calculation.\n" );
293 <      painCave.isFatal = 1;
294 <      simError();
278 >
279 >    // new stuff starts here:
280 >
281 >    fDecomp_->distributeData();
282 >
283 >    int cg1, cg2, atom1, atom2;
284 >    Vector3d d_grp, dag;
285 >    RealType rgrpsq, rgrp;
286 >    RealType vij;
287 >    Vector3d fij, fg;
288 >    pair<int, int> gtypes;
289 >    RealType rCutSq;
290 >    bool in_switching_region;
291 >    RealType sw, dswdr, swderiv;
292 >    vector<int> atomListColumn, atomListRow, atomListLocal;
293 >    InteractionData idat;
294 >    SelfData sdat;
295 >    RealType mf;
296 >
297 >    int loopStart, loopEnd;
298 >
299 >    loopEnd = PAIR_LOOP;
300 >    if (info_->requiresPrepair() ) {
301 >      loopStart = PREPAIR_LOOP;
302 >    } else {
303 >      loopStart = PAIR_LOOP;
304      }
305 +
306 +    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
307 +      
308 +      if (iLoop == loopStart) {
309 +        bool update_nlist = fDecomp_->checkNeighborList();
310 +        if (update_nlist)
311 +          neighborList = fDecomp_->buildNeighborList();
312 +      }
313 +
314 +      for (vector<pair<int, int> >::iterator it = neighborList.begin();
315 +             it != neighborList.end(); ++it) {
316 +        
317 +        cg1 = (*it).first;
318 +        cg2 = (*it).second;
319 +
320 +        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
321 +        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
322 +        curSnapshot->wrapVector(d_grp);        
323 +        rgrpsq = d_grp.lengthSquare();
324 +        rCutSq = groupCutoffMap[gtypes].first;
325 +
326 +        if (rgrpsq < rCutSq) {
327 +          idat.rcut = groupCutoffMap[gtypes].second;
328 +          if (iLoop == PAIR_LOOP) {
329 +            vij *= 0.0;
330 +            fij = V3Zero;
331 +          }
332 +          
333 +          in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp);              
334 +          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
335 +          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
336 +
337 +          for (vector<int>::iterator ia = atomListRow.begin();
338 +               ia != atomListRow.end(); ++ia) {            
339 +            atom1 = (*ia);
340 +            
341 +            for (vector<int>::iterator jb = atomListColumn.begin();
342 +                 jb != atomListColumn.end(); ++jb) {              
343 +              atom2 = (*jb);
344 +              
345 +              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
346 +                
347 +                idat = fDecomp_->fillInteractionData(atom1, atom2);
348 +
349 +                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
350 +                  idat.d = d_grp;
351 +                  idat.r2 = rgrpsq;
352 +                } else {
353 +                  idat.d = fDecomp_->getInteratomicVector(atom1, atom2);
354 +                  curSnapshot->wrapVector(idat.d);
355 +                  idat.r2 = idat.d.lengthSquare();
356 +                }
357 +                
358 +                idat.rij = sqrt(idat.r2);
359 +              
360 +                if (iLoop == PREPAIR_LOOP) {
361 +                  interactionMan_->doPrePair(idat);
362 +                } else {
363 +                  interactionMan_->doPair(idat);
364 +                  vij += idat.vpair;
365 +                  fij += idat.f1;
366 +                  tau -= outProduct(idat.d, idat.f1);
367 +                }
368 +              }
369 +            }
370 +          }
371 +
372 +          if (iLoop == PAIR_LOOP) {
373 +            if (in_switching_region) {
374 +              swderiv = vij * dswdr / rgrp;
375 +              fg = swderiv * d_grp;
376 +
377 +              fij += fg;
378 +
379 +              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
380 +                tau -= outProduct(idat.d, fg);
381 +              }
382 +          
383 +              for (vector<int>::iterator ia = atomListRow.begin();
384 +                   ia != atomListRow.end(); ++ia) {            
385 +                atom1 = (*ia);                
386 +                mf = fDecomp_->getMfactRow(atom1);
387 +                // fg is the force on atom ia due to cutoff group's
388 +                // presence in switching region
389 +                fg = swderiv * d_grp * mf;
390 +                fDecomp_->addForceToAtomRow(atom1, fg);
391 +
392 +                if (atomListRow.size() > 1) {
393 +                  if (info_->usesAtomicVirial()) {
394 +                    // find the distance between the atom
395 +                    // and the center of the cutoff group:
396 +                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
397 +                    tau -= outProduct(dag, fg);
398 +                  }
399 +                }
400 +              }
401 +              for (vector<int>::iterator jb = atomListColumn.begin();
402 +                   jb != atomListColumn.end(); ++jb) {              
403 +                atom2 = (*jb);
404 +                mf = fDecomp_->getMfactColumn(atom2);
405 +                // fg is the force on atom jb due to cutoff group's
406 +                // presence in switching region
407 +                fg = -swderiv * d_grp * mf;
408 +                fDecomp_->addForceToAtomColumn(atom2, fg);
409 +
410 +                if (atomListColumn.size() > 1) {
411 +                  if (info_->usesAtomicVirial()) {
412 +                    // find the distance between the atom
413 +                    // and the center of the cutoff group:
414 +                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
415 +                    tau -= outProduct(dag, fg);
416 +                  }
417 +                }
418 +              }
419 +            }
420 +            //if (!SIM_uses_AtomicVirial) {
421 +            //  tau -= outProduct(d_grp, fij);
422 +            //}
423 +          }
424 +        }
425 +      }
426 +
427 +      if (iLoop == PREPAIR_LOOP) {
428 +        if (info_->requiresPrepair()) {            
429 +          fDecomp_->collectIntermediateData();
430 +          atomListLocal = fDecomp_->getAtomList();
431 +          for (vector<int>::iterator ia = atomListLocal.begin();
432 +               ia != atomListLocal.end(); ++ia) {              
433 +            atom1 = (*ia);            
434 +            sdat = fDecomp_->fillSelfData(atom1);
435 +            interactionMan_->doPreForce(sdat);
436 +          }
437 +          fDecomp_->distributeIntermediateData();        
438 +        }
439 +      }
440 +
441 +    }
442 +    
443 +    fDecomp_->collectData();
444 +    
445 +    if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) {
446 +      atomListLocal = fDecomp_->getAtomList();
447 +      for (vector<int>::iterator ia = atomListLocal.begin();
448 +           ia != atomListLocal.end(); ++ia) {              
449 +        atom1 = (*ia);    
450 +
451 +        if (info_->requiresSkipCorrection()) {
452 +          vector<int> skipList = fDecomp_->getSkipsForAtom(atom1);
453 +          for (vector<int>::iterator jb = skipList.begin();
454 +               jb != skipList.end(); ++jb) {              
455 +            atom2 = (*jb);
456 +            idat = fDecomp_->fillSkipData(atom1, atom2);
457 +            interactionMan_->doSkipCorrection(idat);
458 +          }
459 +        }
460 +          
461 +        if (info_->requiresSelfCorrection()) {
462 +          sdat = fDecomp_->fillSelfData(atom1);
463 +          interactionMan_->doSelfCorrection(sdat);
464 +        }
465 +      }
466 +    }
467 +
468      for (int i=0; i<LR_POT_TYPES;i++){
469        lrPot += longRangePotential[i]; //Quick hack
470      }

Diff Legend

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