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 1545 by gezelter, Fri Apr 8 21:25:19 2011 UTC vs.
Revision 1571 by gezelter, Fri May 27 16:45:44 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 <      nbiMan_->setSimInfo(info_);
74 <      nbiMan_->initialize();
75 <      swfun_ = nbiMan_->getSwitchingFunction();
76 <      decomp_->distributeInitialData();
77 <      info_->setupFortran();
73 >      interactionMan_->setSimInfo(info_);
74 >      interactionMan_->initialize();
75 >      swfun_ = interactionMan_->getSwitchingFunction();
76 >      info_->prepareTopology();
77 >      fDecomp_->distributeInitialData();
78      }
79      
80      preCalculation();  
81 <    calcShortRangeInteraction();
82 <    calcLongRangeInteraction();
81 >    shortRangeInteractions();
82 >    longRangeInteractions();
83      postCalculation();
84      
85    }
# Line 128 | Line 122 | namespace OpenMD {
122      
123    }
124    
125 <  void ForceManager::calcShortRangeInteraction() {
125 >  void ForceManager::shortRangeInteractions() {
126      Molecule* mol;
127      RigidBody* rb;
128      Bond* bond;
# Line 245 | Line 239 | namespace OpenMD {
239      curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
240    }
241    
242 <  void ForceManager::calcLongRangeInteraction() {
242 >  void ForceManager::longRangeInteractions() {
243  
244      // some of this initial stuff will go away:
245      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 268 | Line 262 | namespace OpenMD {
262      }
263      
264      //initialize data before passing to fortran
265 <    RealType longRangePotential[LR_POT_TYPES];
265 >    RealType longRangePotential[N_INTERACTION_FAMILIES];
266      RealType lrPot = 0.0;
267      int isError = 0;
268  
269 <    for (int i=0; i<LR_POT_TYPES;i++){
269 >    // dangerous to iterate over enums, but we'll live on the edge:
270 >    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
271        longRangePotential[i]=0.0; //Initialize array
272      }
273  
274      // new stuff starts here:
275  
276 <    decomp_->distributeData();
276 >    fDecomp_->distributeData();
277  
278 <    int cg1, cg2;
279 <    Vector3d d_grp;
278 >    int cg1, cg2, atom1, atom2;
279 >    Vector3d d_grp, dag;
280      RealType rgrpsq, rgrp;
281      RealType vij;
282      Vector3d fij, fg;
# Line 289 | Line 284 | namespace OpenMD {
284      RealType rCutSq;
285      bool in_switching_region;
286      RealType sw, dswdr, swderiv;
287 <    vector<int> atomListI;
293 <    vector<int> atomListJ;
287 >    vector<int> atomListColumn, atomListRow, atomListLocal;
288      InteractionData idat;
289 +    SelfData sdat;
290 +    RealType mf;
291  
292      int loopStart, loopEnd;
293  
294      loopEnd = PAIR_LOOP;
295 <    if (info_->requiresPrepair_) {
295 >    if (info_->requiresPrepair() ) {
296        loopStart = PREPAIR_LOOP;
297      } else {
298        loopStart = PAIR_LOOP;
# Line 305 | Line 301 | namespace OpenMD {
301      for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
302        
303        if (iLoop == loopStart) {
304 <        bool update_nlist = decomp_->checkNeighborList();
304 >        bool update_nlist = fDecomp_->checkNeighborList();
305          if (update_nlist)
306 <          neighborList = decomp_->buildNeighborList();
306 >          neighborList = fDecomp_->buildNeighborList();
307        }
308  
309        for (vector<pair<int, int> >::iterator it = neighborList.begin();
# Line 316 | Line 312 | namespace OpenMD {
312          cg1 = (*it).first;
313          cg2 = (*it).second;
314  
315 <        gtypes = decomp_->getGroupTypes(cg1, cg2);
316 <        d_grp  = decomp_->getIntergroupVector(cg1, cg2);
315 >        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
316 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
317          curSnapshot->wrapVector(d_grp);        
318          rgrpsq = d_grp.lengthSquare();
319 <        rCutSq = groupCutoffMap(gtypes).first;
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;
324 >            vij *= 0.0;
325              fij = V3Zero;
326            }
327            
328 <          in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, idat.dswdr, rgrp);    
329 <          
330 <          atomListI = decomp_->getAtomsInGroupI(cg1);
331 <          atomListJ = decomp_->getAtomsInGroupJ(cg2);
328 >          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
329 >                                                  rgrp);              
330 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
331 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
332  
333 <          for (vector<int>::iterator ia = atomListI.begin();
334 <               ia != atomListI.end(); ++ia) {            
333 >          for (vector<int>::iterator ia = atomListRow.begin();
334 >               ia != atomListRow.end(); ++ia) {            
335              atom1 = (*ia);
336              
337 <            for (vector<int>::iterator jb = atomListJ.begin();
338 <                 jb != atomListJ.end(); ++jb) {              
337 >            for (vector<int>::iterator jb = atomListColumn.begin();
338 >                 jb != atomListColumn.end(); ++jb) {              
339                atom2 = (*jb);
340                
341 <              if (!decomp_->skipAtomPair(atom1, atom2)) {
341 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
342                  
343 <                if (atomListI.size() == 1 && atomListJ.size() == 1) {
344 <                  idat.d = d_grp;
345 <                  idat.r2 = rgrpsq;
343 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
344 >
345 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
346 >                  *(idat.d) = d_grp;
347 >                  *(idat.r2) = rgrpsq;
348                  } else {
349 <                  idat.d = decomp_->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.r = sqrt(idat.r2);
355 <                decomp_->fillInteractionData(atom1, atom2, idat);
358 <                
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.f);
360 >                  vij += *(idat.vpair);
361 >                  fij += *(idat.f1);
362 >                  tau -= outProduct( *(idat.d), *(idat.f1));
363                  }
364                }
365              }
# Line 375 | Line 372 | namespace OpenMD {
372  
373                fij += fg;
374  
375 <              if (atomListI.size() == 1 && atomListJ.size() == 1) {
376 <                tau -= outProduct(idat.d, fg);
375 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
376 >                tau -= outProduct( *(idat.d), fg);
377                }
378            
379 <              for (vector<int>::iterator ia = atomListI.begin();
380 <                   ia != atomListI.end(); ++ia) {            
379 >              for (vector<int>::iterator ia = atomListRow.begin();
380 >                   ia != atomListRow.end(); ++ia) {            
381                  atom1 = (*ia);                
382 <                mf = decomp_->getMfactI(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;
386 <                decomp_->addForceToAtomI(atom1, fg);
386 >                fDecomp_->addForceToAtomRow(atom1, fg);
387  
388 <                if (atomListI.size() > 1) {
389 <                  if (info_->usesAtomicVirial_) {
388 >                if (atomListRow.size() > 1) {
389 >                  if (info_->usesAtomicVirial()) {
390                      // find the distance between the atom
391                      // and the center of the cutoff group:
392 <                    dag = decomp_->getAtomToGroupVectorI(atom1, cg1);
392 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
393                      tau -= outProduct(dag, fg);
394                    }
395                  }
396                }
397 <              for (vector<int>::iterator jb = atomListJ.begin();
398 <                   jb != atomListJ.end(); ++jb) {              
397 >              for (vector<int>::iterator jb = atomListColumn.begin();
398 >                   jb != atomListColumn.end(); ++jb) {              
399                  atom2 = (*jb);
400 <                mf = decomp_->getMfactJ(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;
404 <                decomp_->addForceToAtomJ(atom2, fg);
404 >                fDecomp_->addForceToAtomColumn(atom2, fg);
405  
406 <                if (atomListJ.size() > 1) {
407 <                  if (info_->usesAtomicVirial_) {
406 >                if (atomListColumn.size() > 1) {
407 >                  if (info_->usesAtomicVirial()) {
408                      // find the distance between the atom
409                      // and the center of the cutoff group:
410 <                    dag = decomp_->getAtomToGroupVectorJ(atom2, cg2);
410 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
411                      tau -= outProduct(dag, fg);
412                    }
413                  }
# Line 424 | Line 421 | namespace OpenMD {
421        }
422  
423        if (iLoop == PREPAIR_LOOP) {
424 <        if (info_->requiresPrepair_) {            
425 <          decomp_->collectIntermediateData();
426 <          atomList = decomp_->getAtomList();
427 <          for (vector<int>::iterator ia = atomList.begin();
428 <               ia != atomList.end(); ++ia) {              
432 <            atom1 = (*ia);            
433 <            decomp_->populateSelfData(atom1, SelfData sdat);
424 >        if (info_->requiresPrepair()) {            
425 >          fDecomp_->collectIntermediateData();
426 >
427 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
428 >            sdat = fDecomp_->fillSelfData(atom1);
429              interactionMan_->doPreForce(sdat);
430            }
431 <          decomp_->distributeIntermediateData();        
431 >
432 >          fDecomp_->distributeIntermediateData();        
433          }
434        }
435  
436      }
437      
438 <    decomp_->collectData();
438 >    fDecomp_->collectData();
439      
440 <    if (info_->requiresSkipCorrection_ || info_->requiresSelfCorrection_) {
441 <      atomList = decomp_->getAtomList();
442 <      for (vector<int>::iterator ia = atomList.begin();
447 <           ia != atomList.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 = decomp_->getSkipsForAtom(atom1);
446 <          for (vector<int>::iterator jb = skipList.begin();
447 <               jb != skipList.end(); ++jb) {              
448 <            atom2 = (*jb);
449 <            decomp_->populateSkipData(atom1, atom2, InteractionData idat);
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 –          decomp_->populateSelfData(atom1, SelfData sdat);
462 –          interactionMan_->doSelfCorrection(sdat);
454        }
464 –      
465 –      
455      }
456 +    
457 +    if (info_->requiresSelfCorrection()) {
458  
459 <    for (int i=0; i<LR_POT_TYPES;i++){
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
469      }
470          
471      //store the tau and long range potential    
472      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
473 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
474 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
473 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
474 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
475    }
476  
477    

Diff Legend

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