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 1550 by gezelter, Wed Apr 27 21:49:59 2011 UTC vs.
Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC

# Line 55 | Line 55
55   #include "primitives/Bend.hpp"
56   #include "primitives/Torsion.hpp"
57   #include "primitives/Inversion.hpp"
58 – #include "parallel/ForceMatrixDecomposition.hpp"
58   #include "nonbonded/NonBondedInteraction.hpp"
59 + #include "parallel/ForceMatrixDecomposition.hpp"
60  
61   using namespace std;
62   namespace OpenMD {
63    
64    ForceManager::ForceManager(SimInfo * info) : info_(info) {
65  
66 – #ifdef IS_MPI
66      fDecomp_ = new ForceMatrixDecomposition(info_);
68 – #else
69 –    // fDecomp_ = new ForceSerialDecomposition(info);
70 – #endif
67    }
68    
69    void ForceManager::calcForces() {
70      
71 <    if (!info_->isFortranInitialized()) {
71 >    if (!info_->isTopologyDone()) {
72        info_->update();
73        interactionMan_->setSimInfo(info_);
74        interactionMan_->initialize();
75        swfun_ = interactionMan_->getSwitchingFunction();
76 +      info_->prepareTopology();
77        fDecomp_->distributeInitialData();
81 –      info_->setupFortran();
78      }
79      
80      preCalculation();  
# Line 323 | Line 319 | namespace OpenMD {
319          rCutSq = groupCutoffMap[gtypes].first;
320  
321          if (rgrpsq < rCutSq) {
322 <          idat.rcut = groupCutoffMap[gtypes].second;
322 >          *(idat.rcut) = groupCutoffMap[gtypes].second;
323            if (iLoop == PAIR_LOOP) {
324              vij *= 0.0;
325              fij = V3Zero;
326            }
327            
328 <          in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp);              
328 >          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
329 >                                                  rgrp);              
330            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
331            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
332  
# Line 346 | Line 343 | namespace OpenMD {
343                  idat = fDecomp_->fillInteractionData(atom1, atom2);
344  
345                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
346 <                  idat.d = d_grp;
347 <                  idat.r2 = rgrpsq;
346 >                  *(idat.d) = d_grp;
347 >                  *(idat.r2) = rgrpsq;
348                  } else {
349 <                  idat.d = fDecomp_->getInteratomicVector(atom1, atom2);
350 <                  curSnapshot->wrapVector(idat.d);
351 <                  idat.r2 = idat.d.lengthSquare();
349 >                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
350 >                  curSnapshot->wrapVector( *(idat.d) );
351 >                  *(idat.r2) = idat.d->lengthSquare();
352                  }
353                  
354 <                idat.rij = sqrt(idat.r2);
354 >                *(idat.rij) = sqrt( *(idat.r2) );
355                
356                  if (iLoop == PREPAIR_LOOP) {
357                    interactionMan_->doPrePair(idat);
358                  } else {
359                    interactionMan_->doPair(idat);
360 <                  vij += idat.vpair;
361 <                  fij += idat.f1;
362 <                  tau -= outProduct(idat.d, idat.f1);
360 >                  vij += *(idat.vpair);
361 >                  fij += *(idat.f1);
362 >                  tau -= outProduct( *(idat.d), *(idat.f1));
363                  }
364                }
365              }
# Line 376 | Line 373 | namespace OpenMD {
373                fij += fg;
374  
375                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
376 <                tau -= outProduct(idat.d, fg);
376 >                tau -= outProduct( *(idat.d), fg);
377                }
378            
379                for (vector<int>::iterator ia = atomListRow.begin();
380                     ia != atomListRow.end(); ++ia) {            
381                  atom1 = (*ia);                
382 <                mf = fDecomp_->getMfactRow(atom1);
382 >                mf = fDecomp_->getMassFactorRow(atom1);
383                  // fg is the force on atom ia due to cutoff group's
384                  // presence in switching region
385                  fg = swderiv * d_grp * mf;
# Line 400 | Line 397 | namespace OpenMD {
397                for (vector<int>::iterator jb = atomListColumn.begin();
398                     jb != atomListColumn.end(); ++jb) {              
399                  atom2 = (*jb);
400 <                mf = fDecomp_->getMfactColumn(atom2);
400 >                mf = fDecomp_->getMassFactorColumn(atom2);
401                  // fg is the force on atom jb due to cutoff group's
402                  // presence in switching region
403                  fg = -swderiv * d_grp * mf;
# Line 426 | Line 423 | namespace OpenMD {
423        if (iLoop == PREPAIR_LOOP) {
424          if (info_->requiresPrepair()) {            
425            fDecomp_->collectIntermediateData();
426 <          atomListLocal = fDecomp_->getAtomList();
427 <          for (vector<int>::iterator ia = atomListLocal.begin();
431 <               ia != atomListLocal.end(); ++ia) {              
432 <            atom1 = (*ia);            
426 >
427 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
428              sdat = fDecomp_->fillSelfData(atom1);
429              interactionMan_->doPreForce(sdat);
430            }
431 +
432            fDecomp_->distributeIntermediateData();        
433          }
434        }
# Line 441 | Line 437 | namespace OpenMD {
437      
438      fDecomp_->collectData();
439      
440 <    if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) {
441 <      atomListLocal = fDecomp_->getAtomList();
442 <      for (vector<int>::iterator ia = atomListLocal.begin();
447 <           ia != atomListLocal.end(); ++ia) {              
448 <        atom1 = (*ia);    
440 >    if ( info_->requiresSkipCorrection() ) {
441 >      
442 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
443  
444 <        if (info_->requiresSkipCorrection()) {
445 <          vector<int> skipList = fDecomp_->getSkipsForAtom(atom1);
446 <          for (vector<int>::iterator jb = skipList.begin();
447 <               jb != skipList.end(); ++jb) {              
448 <            atom2 = (*jb);
449 <            idat = fDecomp_->fillSkipData(atom1, atom2);
450 <            interactionMan_->doSkipCorrection(idat);
451 <          }
444 >        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
445 >        
446 >        for (vector<int>::iterator jb = skipList.begin();
447 >             jb != skipList.end(); ++jb) {        
448 >    
449 >          atom2 = (*jb);
450 >          idat = fDecomp_->fillSkipData(atom1, atom2);
451 >          interactionMan_->doSkipCorrection(idat);
452 >
453          }
459 –          
460 –        if (info_->requiresSelfCorrection()) {
461 –          sdat = fDecomp_->fillSelfData(atom1);
462 –          interactionMan_->doSelfCorrection(sdat);
463 –        }
454        }
455      }
456 +    
457 +    if (info_->requiresSelfCorrection()) {
458  
459 +      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
460 +        sdat = fDecomp_->fillSelfData(atom1);
461 +        interactionMan_->doSelfCorrection(sdat);
462 +      }
463 +
464 +    }
465 +
466      // dangerous to iterate over enums, but we'll live on the edge:
467      for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
468        lrPot += longRangePotential[i]; //Quick hack

Diff Legend

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