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 1551 by gezelter, Thu Apr 28 18:38:21 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 {
# Line 66 | Line 64 | namespace OpenMD {
64    ForceManager::ForceManager(SimInfo * info) : info_(info) {
65  
66   #ifdef IS_MPI
67 <    decomp_ = new ForceDecomposition(info_);
67 >    fDecomp_ = new ForceMatrixDecomposition(info_);
68   #else
69 <    // decomp_ = new SerialDecomposition(info);
69 >    // fDecomp_ = new ForceSerialDecomposition(info);
70   #endif
71    }
72    
# Line 79 | Line 77 | namespace OpenMD {
77        interactionMan_->setSimInfo(info_);
78        interactionMan_->initialize();
79        swfun_ = interactionMan_->getSwitchingFunction();
80 <      decomp_->distributeInitialData();
80 >      fDecomp_->distributeInitialData();
81        info_->setupFortran();
82      }
83      
# Line 268 | Line 266 | namespace OpenMD {
266      }
267      
268      //initialize data before passing to fortran
269 <    RealType longRangePotential[LR_POT_TYPES];
269 >    RealType longRangePotential[N_INTERACTION_FAMILIES];
270      RealType lrPot = 0.0;
271      int isError = 0;
272  
273 <    for (int i=0; i<LR_POT_TYPES;i++){
273 >    // dangerous to iterate over enums, but we'll live on the edge:
274 >    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
275        longRangePotential[i]=0.0; //Initialize array
276      }
277  
278      // new stuff starts here:
279  
280 <    decomp_->distributeData();
280 >    fDecomp_->distributeData();
281  
282      int cg1, cg2, atom1, atom2;
283      Vector3d d_grp, dag;
284      RealType rgrpsq, rgrp;
285 <    Vector<RealType, 4> vij;
285 >    RealType vij;
286      Vector3d fij, fg;
287      pair<int, int> gtypes;
288      RealType rCutSq;
289      bool in_switching_region;
290      RealType sw, dswdr, swderiv;
291 <    vector<int> atomListI, atomListJ, atomList;
291 >    vector<int> atomListColumn, atomListRow, atomListLocal;
292      InteractionData idat;
293      SelfData sdat;
294      RealType mf;
# Line 306 | Line 305 | namespace OpenMD {
305      for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
306        
307        if (iLoop == loopStart) {
308 <        bool update_nlist = decomp_->checkNeighborList();
308 >        bool update_nlist = fDecomp_->checkNeighborList();
309          if (update_nlist)
310 <          neighborList = decomp_->buildNeighborList();
310 >          neighborList = fDecomp_->buildNeighborList();
311        }
312  
313        for (vector<pair<int, int> >::iterator it = neighborList.begin();
# Line 317 | Line 316 | namespace OpenMD {
316          cg1 = (*it).first;
317          cg2 = (*it).second;
318  
319 <        gtypes = decomp_->getGroupTypes(cg1, cg2);
320 <        d_grp  = decomp_->getIntergroupVector(cg1, cg2);
319 >        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
320 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
321          curSnapshot->wrapVector(d_grp);        
322          rgrpsq = d_grp.lengthSquare();
323          rCutSq = groupCutoffMap[gtypes].first;
# Line 331 | Line 330 | namespace OpenMD {
330            }
331            
332            in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp);              
333 <          atomListI = decomp_->getAtomsInGroupI(cg1);
334 <          atomListJ = decomp_->getAtomsInGroupJ(cg2);
333 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
334 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
335  
336 <          for (vector<int>::iterator ia = atomListI.begin();
337 <               ia != atomListI.end(); ++ia) {            
336 >          for (vector<int>::iterator ia = atomListRow.begin();
337 >               ia != atomListRow.end(); ++ia) {            
338              atom1 = (*ia);
339              
340 <            for (vector<int>::iterator jb = atomListJ.begin();
341 <                 jb != atomListJ.end(); ++jb) {              
340 >            for (vector<int>::iterator jb = atomListColumn.begin();
341 >                 jb != atomListColumn.end(); ++jb) {              
342                atom2 = (*jb);
343                
344 <              if (!decomp_->skipAtomPair(atom1, atom2)) {
344 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
345                  
346 <                idat = decomp_->fillInteractionData(atom1, atom2);
346 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
347  
348 <                if (atomListI.size() == 1 && atomListJ.size() == 1) {
348 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
349                    idat.d = d_grp;
350                    idat.r2 = rgrpsq;
351                  } else {
352 <                  idat.d = decomp_->getInteratomicVector(atom1, atom2);
352 >                  idat.d = fDecomp_->getInteratomicVector(atom1, atom2);
353                    curSnapshot->wrapVector(idat.d);
354                    idat.r2 = idat.d.lengthSquare();
355                  }
# Line 376 | Line 375 | namespace OpenMD {
375  
376                fij += fg;
377  
378 <              if (atomListI.size() == 1 && atomListJ.size() == 1) {
378 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
379                  tau -= outProduct(idat.d, fg);
380                }
381            
382 <              for (vector<int>::iterator ia = atomListI.begin();
383 <                   ia != atomListI.end(); ++ia) {            
382 >              for (vector<int>::iterator ia = atomListRow.begin();
383 >                   ia != atomListRow.end(); ++ia) {            
384                  atom1 = (*ia);                
385 <                mf = decomp_->getMfactI(atom1);
385 >                mf = fDecomp_->getMfactRow(atom1);
386                  // fg is the force on atom ia due to cutoff group's
387                  // presence in switching region
388                  fg = swderiv * d_grp * mf;
389 <                decomp_->addForceToAtomI(atom1, fg);
389 >                fDecomp_->addForceToAtomRow(atom1, fg);
390  
391 <                if (atomListI.size() > 1) {
391 >                if (atomListRow.size() > 1) {
392                    if (info_->usesAtomicVirial()) {
393                      // find the distance between the atom
394                      // and the center of the cutoff group:
395 <                    dag = decomp_->getAtomToGroupVectorI(atom1, cg1);
395 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
396                      tau -= outProduct(dag, fg);
397                    }
398                  }
399                }
400 <              for (vector<int>::iterator jb = atomListJ.begin();
401 <                   jb != atomListJ.end(); ++jb) {              
400 >              for (vector<int>::iterator jb = atomListColumn.begin();
401 >                   jb != atomListColumn.end(); ++jb) {              
402                  atom2 = (*jb);
403 <                mf = decomp_->getMfactJ(atom2);
403 >                mf = fDecomp_->getMfactColumn(atom2);
404                  // fg is the force on atom jb due to cutoff group's
405                  // presence in switching region
406                  fg = -swderiv * d_grp * mf;
407 <                decomp_->addForceToAtomJ(atom2, fg);
407 >                fDecomp_->addForceToAtomColumn(atom2, fg);
408  
409 <                if (atomListJ.size() > 1) {
409 >                if (atomListColumn.size() > 1) {
410                    if (info_->usesAtomicVirial()) {
411                      // find the distance between the atom
412                      // and the center of the cutoff group:
413 <                    dag = decomp_->getAtomToGroupVectorJ(atom2, cg2);
413 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
414                      tau -= outProduct(dag, fg);
415                    }
416                  }
# Line 426 | Line 425 | namespace OpenMD {
425  
426        if (iLoop == PREPAIR_LOOP) {
427          if (info_->requiresPrepair()) {            
428 <          decomp_->collectIntermediateData();
429 <          atomList = decomp_->getAtomList();
430 <          for (vector<int>::iterator ia = atomList.begin();
431 <               ia != atomList.end(); ++ia) {              
428 >          fDecomp_->collectIntermediateData();
429 >          atomListLocal = fDecomp_->getAtomList();
430 >          for (vector<int>::iterator ia = atomListLocal.begin();
431 >               ia != atomListLocal.end(); ++ia) {              
432              atom1 = (*ia);            
433 <            sdat = decomp_->fillSelfData(atom1);
433 >            sdat = fDecomp_->fillSelfData(atom1);
434              interactionMan_->doPreForce(sdat);
435            }
436 <          decomp_->distributeIntermediateData();        
436 >          fDecomp_->distributeIntermediateData();        
437          }
438        }
439  
440      }
441      
442 <    decomp_->collectData();
442 >    fDecomp_->collectData();
443      
444      if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) {
445 <      atomList = decomp_->getAtomList();
446 <      for (vector<int>::iterator ia = atomList.begin();
447 <           ia != atomList.end(); ++ia) {              
445 >      atomListLocal = fDecomp_->getAtomList();
446 >      for (vector<int>::iterator ia = atomListLocal.begin();
447 >           ia != atomListLocal.end(); ++ia) {              
448          atom1 = (*ia);    
449  
450          if (info_->requiresSkipCorrection()) {
451 <          vector<int> skipList = decomp_->getSkipsForAtom(atom1);
451 >          vector<int> skipList = fDecomp_->getSkipsForAtom(atom1);
452            for (vector<int>::iterator jb = skipList.begin();
453                 jb != skipList.end(); ++jb) {              
454              atom2 = (*jb);
455 <            idat = decomp_->fillSkipData(atom1, atom2);
455 >            idat = fDecomp_->fillSkipData(atom1, atom2);
456              interactionMan_->doSkipCorrection(idat);
457            }
458          }
459            
460          if (info_->requiresSelfCorrection()) {
461 <          sdat = decomp_->fillSelfData(atom1);
461 >          sdat = fDecomp_->fillSelfData(atom1);
462            interactionMan_->doSelfCorrection(sdat);
463 +        }
464        }
465 –      
466 –      
465      }
466  
467 <    for (int i=0; i<LR_POT_TYPES;i++){
467 >    // dangerous to iterate over enums, but we'll live on the edge:
468 >    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
469        lrPot += longRangePotential[i]; //Quick hack
470      }
471          
472      //store the tau and long range potential    
473      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
474 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
475 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
474 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
475 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
476    }
477  
478    

Diff Legend

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