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 1546 by gezelter, Sun Apr 10 15:16:39 2011 UTC vs.
Revision 1575 by gezelter, Fri Jun 3 21:39:49 2011 UTC

# Line 49 | Line 49
49  
50   #include "brains/ForceManager.hpp"
51   #include "primitives/Molecule.hpp"
52 – #include "UseTheForce/doForces_interface.h"
52   #define __OPENMD_C
54 – #include "UseTheForce/DarkSide/fInteractionMap.h"
53   #include "utils/simError.h"
54   #include "primitives/Bond.hpp"
55   #include "primitives/Bend.hpp"
56   #include "primitives/Torsion.hpp"
57   #include "primitives/Inversion.hpp"
58 < #include "parallel/ForceDecomposition.hpp"
59 < //#include "parallel/SerialDecomposition.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
69 <    decomp_ = new ForceDecomposition(info_);
70 < #else
71 <    // decomp_ = new SerialDecomposition(info);
72 < #endif
66 >    fDecomp_ = new ForceMatrixDecomposition(info_);
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 <      decomp_->distributeInitialData();
77 <      info_->setupFortran();
76 >      info_->prepareTopology();
77 >      fDecomp_->distributeInitialData();
78      }
79      
80      preCalculation();  
# Line 267 | Line 261 | namespace OpenMD {
261        rc = pos;
262      }
263      
270 –    //initialize data before passing to fortran
271 –    RealType longRangePotential[LR_POT_TYPES];
272 –    RealType lrPot = 0.0;
273 –    int isError = 0;
274 –
275 –    for (int i=0; i<LR_POT_TYPES;i++){
276 –      longRangePotential[i]=0.0; //Initialize array
277 –    }
278 –
264      // new stuff starts here:
265 <
266 <    decomp_->distributeData();
265 >    fDecomp_->zeroWorkArrays();
266 >    fDecomp_->distributeData();
267  
268      int cg1, cg2, atom1, atom2;
269      Vector3d d_grp, dag;
270      RealType rgrpsq, rgrp;
271 <    Vector<RealType, 4> vij;
271 >    RealType vij;
272      Vector3d fij, fg;
273      pair<int, int> gtypes;
274      RealType rCutSq;
275      bool in_switching_region;
276      RealType sw, dswdr, swderiv;
277 <    vector<int> atomListI, atomListJ, atomList;
277 >    vector<int> atomListColumn, atomListRow, atomListLocal;
278      InteractionData idat;
279      SelfData sdat;
280      RealType mf;
281 +    potVec pot(0.0);
282 +    potVec longRangePotential(0.0);
283 +    RealType lrPot;
284  
285      int loopStart, loopEnd;
286  
# Line 306 | Line 294 | namespace OpenMD {
294      for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
295        
296        if (iLoop == loopStart) {
297 <        bool update_nlist = decomp_->checkNeighborList();
297 >        bool update_nlist = fDecomp_->checkNeighborList();
298          if (update_nlist)
299 <          neighborList = decomp_->buildNeighborList();
299 >          neighborList = fDecomp_->buildNeighborList();
300        }
301  
302        for (vector<pair<int, int> >::iterator it = neighborList.begin();
# Line 317 | Line 305 | namespace OpenMD {
305          cg1 = (*it).first;
306          cg2 = (*it).second;
307  
308 <        gtypes = decomp_->getGroupTypes(cg1, cg2);
309 <        d_grp  = decomp_->getIntergroupVector(cg1, cg2);
308 >        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
309 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
310          curSnapshot->wrapVector(d_grp);        
311          rgrpsq = d_grp.lengthSquare();
312          rCutSq = groupCutoffMap[gtypes].first;
313  
314          if (rgrpsq < rCutSq) {
315 <          idat.rcut = groupCutoffMap[gtypes].second;
315 >          *(idat.rcut) = groupCutoffMap[gtypes].second;
316            if (iLoop == PAIR_LOOP) {
317              vij *= 0.0;
318              fij = V3Zero;
319            }
320            
321 <          in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp);              
322 <          atomListI = decomp_->getAtomsInGroupI(cg1);
323 <          atomListJ = decomp_->getAtomsInGroupJ(cg2);
321 >          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
322 >                                                  rgrp);              
323 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
324 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
325  
326 <          for (vector<int>::iterator ia = atomListI.begin();
327 <               ia != atomListI.end(); ++ia) {            
326 >          for (vector<int>::iterator ia = atomListRow.begin();
327 >               ia != atomListRow.end(); ++ia) {            
328              atom1 = (*ia);
329              
330 <            for (vector<int>::iterator jb = atomListJ.begin();
331 <                 jb != atomListJ.end(); ++jb) {              
330 >            for (vector<int>::iterator jb = atomListColumn.begin();
331 >                 jb != atomListColumn.end(); ++jb) {              
332                atom2 = (*jb);
333                
334 <              if (!decomp_->skipAtomPair(atom1, atom2)) {
334 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
335                  
336 <                idat = decomp_->fillInteractionData(atom1, atom2);
336 >                pot *= 0.0;
337  
338 <                if (atomListI.size() == 1 && atomListJ.size() == 1) {
339 <                  idat.d = d_grp;
340 <                  idat.r2 = rgrpsq;
338 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
339 >                *(idat.pot) = pot;
340 >
341 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
342 >                  *(idat.d) = d_grp;
343 >                  *(idat.r2) = rgrpsq;
344                  } else {
345 <                  idat.d = decomp_->getInteratomicVector(atom1, atom2);
346 <                  curSnapshot->wrapVector(idat.d);
347 <                  idat.r2 = idat.d.lengthSquare();
345 >                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
346 >                  curSnapshot->wrapVector( *(idat.d) );
347 >                  *(idat.r2) = idat.d->lengthSquare();
348                  }
349                  
350 <                idat.rij = sqrt(idat.r2);
350 >                *(idat.rij) = sqrt( *(idat.r2) );
351                
352                  if (iLoop == PREPAIR_LOOP) {
353                    interactionMan_->doPrePair(idat);
354                  } else {
355                    interactionMan_->doPair(idat);
356 <                  vij += idat.vpair;
357 <                  fij += idat.f1;
358 <                  tau -= outProduct(idat.d, idat.f1);
356 >                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
357 >                  vij += *(idat.vpair);
358 >                  fij += *(idat.f1);
359 >                  tau -= outProduct( *(idat.d), *(idat.f1));
360                  }
361                }
362              }
# Line 376 | Line 369 | namespace OpenMD {
369  
370                fij += fg;
371  
372 <              if (atomListI.size() == 1 && atomListJ.size() == 1) {
373 <                tau -= outProduct(idat.d, fg);
372 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
373 >                tau -= outProduct( *(idat.d), fg);
374                }
375            
376 <              for (vector<int>::iterator ia = atomListI.begin();
377 <                   ia != atomListI.end(); ++ia) {            
376 >              for (vector<int>::iterator ia = atomListRow.begin();
377 >                   ia != atomListRow.end(); ++ia) {            
378                  atom1 = (*ia);                
379 <                mf = decomp_->getMfactI(atom1);
379 >                mf = fDecomp_->getMassFactorRow(atom1);
380                  // fg is the force on atom ia due to cutoff group's
381                  // presence in switching region
382                  fg = swderiv * d_grp * mf;
383 <                decomp_->addForceToAtomI(atom1, fg);
383 >                fDecomp_->addForceToAtomRow(atom1, fg);
384  
385 <                if (atomListI.size() > 1) {
385 >                if (atomListRow.size() > 1) {
386                    if (info_->usesAtomicVirial()) {
387                      // find the distance between the atom
388                      // and the center of the cutoff group:
389 <                    dag = decomp_->getAtomToGroupVectorI(atom1, cg1);
389 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
390                      tau -= outProduct(dag, fg);
391                    }
392                  }
393                }
394 <              for (vector<int>::iterator jb = atomListJ.begin();
395 <                   jb != atomListJ.end(); ++jb) {              
394 >              for (vector<int>::iterator jb = atomListColumn.begin();
395 >                   jb != atomListColumn.end(); ++jb) {              
396                  atom2 = (*jb);
397 <                mf = decomp_->getMfactJ(atom2);
397 >                mf = fDecomp_->getMassFactorColumn(atom2);
398                  // fg is the force on atom jb due to cutoff group's
399                  // presence in switching region
400                  fg = -swderiv * d_grp * mf;
401 <                decomp_->addForceToAtomJ(atom2, fg);
401 >                fDecomp_->addForceToAtomColumn(atom2, fg);
402  
403 <                if (atomListJ.size() > 1) {
403 >                if (atomListColumn.size() > 1) {
404                    if (info_->usesAtomicVirial()) {
405                      // find the distance between the atom
406                      // and the center of the cutoff group:
407 <                    dag = decomp_->getAtomToGroupVectorJ(atom2, cg2);
407 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
408                      tau -= outProduct(dag, fg);
409                    }
410                  }
# Line 426 | Line 419 | namespace OpenMD {
419  
420        if (iLoop == PREPAIR_LOOP) {
421          if (info_->requiresPrepair()) {            
422 <          decomp_->collectIntermediateData();
423 <          atomList = decomp_->getAtomList();
424 <          for (vector<int>::iterator ia = atomList.begin();
425 <               ia != atomList.end(); ++ia) {              
433 <            atom1 = (*ia);            
434 <            sdat = decomp_->fillSelfData(atom1);
422 >          fDecomp_->collectIntermediateData();
423 >
424 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
425 >            sdat = fDecomp_->fillSelfData(atom1);
426              interactionMan_->doPreForce(sdat);
427            }
428 <          decomp_->distributeIntermediateData();        
428 >
429 >          fDecomp_->distributeIntermediateData();        
430          }
431        }
432  
433      }
434      
435 <    decomp_->collectData();
435 >    fDecomp_->collectData();
436      
437 <    if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) {
438 <      atomList = decomp_->getAtomList();
439 <      for (vector<int>::iterator ia = atomList.begin();
448 <           ia != atomList.end(); ++ia) {              
449 <        atom1 = (*ia);    
437 >    if ( info_->requiresSkipCorrection() ) {
438 >      
439 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
440  
441 <        if (info_->requiresSkipCorrection()) {
442 <          vector<int> skipList = decomp_->getSkipsForAtom(atom1);
443 <          for (vector<int>::iterator jb = skipList.begin();
444 <               jb != skipList.end(); ++jb) {              
445 <            atom2 = (*jb);
446 <            idat = decomp_->fillSkipData(atom1, atom2);
447 <            interactionMan_->doSkipCorrection(idat);
448 <          }
441 >        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
442 >        
443 >        for (vector<int>::iterator jb = skipList.begin();
444 >             jb != skipList.end(); ++jb) {        
445 >    
446 >          atom2 = (*jb);
447 >          idat = fDecomp_->fillSkipData(atom1, atom2);
448 >          interactionMan_->doSkipCorrection(idat);
449 >
450          }
460 –          
461 –        if (info_->requiresSelfCorrection()) {
462 –          sdat = decomp_->fillSelfData(atom1);
463 –          interactionMan_->doSelfCorrection(sdat);
451        }
465 –      
466 –      
452      }
453 +    
454 +    if (info_->requiresSelfCorrection()) {
455  
456 <    for (int i=0; i<LR_POT_TYPES;i++){
457 <      lrPot += longRangePotential[i]; //Quick hack
456 >      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
457 >        sdat = fDecomp_->fillSelfData(atom1);
458 >        interactionMan_->doSelfCorrection(sdat);
459 >      }
460 >
461      }
462 <        
462 >
463 >    longRangePotential = fDecomp_->getLongRangePotential();
464 >    lrPot = longRangePotential.sum();
465 >
466      //store the tau and long range potential    
467      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
468 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
469 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
468 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
469 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
470    }
471  
472    

Diff Legend

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