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 1576 by gezelter, Wed Jun 8 16:05:07 2011 UTC vs.
Revision 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC

# Line 59 | Line 59
59   #include "nonbonded/NonBondedInteraction.hpp"
60   #include "parallel/ForceMatrixDecomposition.hpp"
61  
62 + #include <cstdio>
63 + #include <iostream>
64 + #include <iomanip>
65 +
66   using namespace std;
67   namespace OpenMD {
68    
69    ForceManager::ForceManager(SimInfo * info) : info_(info) {
70      forceField_ = info_->getForceField();
71 <    fDecomp_ = new ForceMatrixDecomposition(info_);
71 >    interactionMan_ = new InteractionManager();
72 >    fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
73    }
74  
75    /**
# Line 121 | Line 126 | namespace OpenMD {
126          painCave.isFatal = 0;
127          painCave.severity = OPENMD_INFO;
128          simError();
129 <      }            
129 >      }
130      }
131  
132 +    fDecomp_->setUserCutoff(rCut_);
133 +    interactionMan_->setCutoffRadius(rCut_);
134 +
135      map<string, CutoffMethod> stringToCutoffMethod;
136      stringToCutoffMethod["HARD"] = HARD;
137      stringToCutoffMethod["SWITCHED"] = SWITCHED;
# Line 194 | Line 202 | namespace OpenMD {
202        simError();
203        cutoffPolicy_ = TRADITIONAL;        
204      }
205 +    fDecomp_->setCutoffPolicy(cutoffPolicy_);
206    }
207  
208    /**
# Line 205 | Line 214 | namespace OpenMD {
214     */
215    void ForceManager::setupSwitching() {
216      Globals* simParams_ = info_->getSimParams();
217 +
218 +    // create the switching function object:
219 +    switcher_ = new SwitchingFunction();
220      
221      if (simParams_->haveSwitchingRadius()) {
222        rSwitch_ = simParams_->getSwitchingRadius();
223        if (rSwitch_ > rCut_) {        
224          sprintf(painCave.errMsg,
225 <                "ForceManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
226 <                rSwitch_, rCut_);
225 >                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
226 >                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
227          painCave.isFatal = 1;
228          painCave.severity = OPENMD_ERROR;
229          simError();
# Line 227 | Line 239 | namespace OpenMD {
239        simError();
240      }          
241      
242 +    // Default to cubic switching function.
243 +    sft_ = cubic;
244      if (simParams_->haveSwitchingFunctionType()) {
245        string funcType = simParams_->getSwitchingFunctionType();
246        toUpper(funcType);
# Line 250 | Line 264 | namespace OpenMD {
264      }
265      switcher_->setSwitchType(sft_);
266      switcher_->setSwitch(rSwitch_, rCut_);
267 +    interactionMan_->setSwitchingRadius(rSwitch_);
268    }
269    
270    void ForceManager::initialize() {
# Line 468 | Line 483 | namespace OpenMD {
483    
484    void ForceManager::longRangeInteractions() {
485  
471 –    // some of this initial stuff will go away:
486      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
487      DataStorage* config = &(curSnapshot->atomData);
488      DataStorage* cgConfig = &(curSnapshot->cgData);
475 –    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
476 –    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
477 –    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
478 –    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
479 –    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
480 –    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
481 –    RealType* rc;    
489  
490 <    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
491 <      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
490 >    //calculate the center of mass of cutoff group
491 >
492 >    SimInfo::MoleculeIterator mi;
493 >    Molecule* mol;
494 >    Molecule::CutoffGroupIterator ci;
495 >    CutoffGroup* cg;
496 >
497 >    if(info_->getNCutoffGroups() > 0){      
498 >      for (mol = info_->beginMolecule(mi); mol != NULL;
499 >           mol = info_->nextMolecule(mi)) {
500 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
501 >            cg = mol->nextCutoffGroup(ci)) {
502 >          cg->updateCOM();
503 >        }
504 >      }      
505      } else {
506        // center of mass of the group is the same as position of the atom  
507        // if cutoff group does not exist
508 <      rc = pos;
508 >      cgConfig->position = config->position;
509      }
510 <    
491 <    // new stuff starts here:
510 >
511      fDecomp_->zeroWorkArrays();
512      fDecomp_->distributeData();
513 <
514 <    int cg1, cg2, atom1, atom2;
515 <    Vector3d d_grp, dag;
516 <    RealType rgrpsq, rgrp;
513 >    
514 >    int cg1, cg2, atom1, atom2, topoDist;
515 >    Vector3d d_grp, dag, d;
516 >    RealType rgrpsq, rgrp, r2, r;
517 >    RealType electroMult, vdwMult;
518      RealType vij;
519 <    Vector3d fij, fg;
519 >    Vector3d fij, fg, f1;
520      tuple3<RealType, RealType, RealType> cuts;
521      RealType rCutSq;
522      bool in_switching_region;
# Line 505 | Line 525 | namespace OpenMD {
525      InteractionData idat;
526      SelfData sdat;
527      RealType mf;
508 –    potVec pot(0.0);
509 –    potVec longRangePotential(0.0);
528      RealType lrPot;
529 +    RealType vpair;
530 +    potVec longRangePotential(0.0);
531 +    potVec workPot(0.0);
532  
533      int loopStart, loopEnd;
534  
535 +    idat.vdwMult = &vdwMult;
536 +    idat.electroMult = &electroMult;
537 +    idat.pot = &workPot;
538 +    sdat.pot = fDecomp_->getEmbeddingPotential();
539 +    idat.vpair = &vpair;
540 +    idat.f1 = &f1;
541 +    idat.sw = &sw;
542 +    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
543 +    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
544 +    
545      loopEnd = PAIR_LOOP;
546      if (info_->requiresPrepair() ) {
547        loopStart = PREPAIR_LOOP;
548      } else {
549        loopStart = PAIR_LOOP;
550      }
551 <
552 <    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
553 <      
551 >  
552 >    for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
553 >    
554        if (iLoop == loopStart) {
555          bool update_nlist = fDecomp_->checkNeighborList();
556          if (update_nlist)
557            neighborList = fDecomp_->buildNeighborList();
558 <      }
559 <
558 >      }      
559 >        
560        for (vector<pair<int, int> >::iterator it = neighborList.begin();
561               it != neighborList.end(); ++it) {
562 <        
562 >                
563          cg1 = (*it).first;
564          cg2 = (*it).second;
565          
# Line 541 | Line 572 | namespace OpenMD {
572          rCutSq = cuts.second;
573  
574          if (rgrpsq < rCutSq) {
575 <          *(idat.rcut) = cuts.first;
575 >          idat.rcut = &cuts.first;
576            if (iLoop == PAIR_LOOP) {
577              vij *= 0.0;
578              fij = V3Zero;
579            }
580            
581 <          in_switching_region = switcher_->getSwitch(rgrpsq, *(idat.sw), dswdr,
581 >          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
582                                                       rgrp);
583                
584            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
# Line 560 | Line 591 | namespace OpenMD {
591              for (vector<int>::iterator jb = atomListColumn.begin();
592                   jb != atomListColumn.end(); ++jb) {              
593                atom2 = (*jb);
594 <              
594 >            
595                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
596 +                vpair = 0.0;
597 +                workPot = 0.0;
598 +                f1 = V3Zero;
599 +
600 +                fDecomp_->fillInteractionData(idat, atom1, atom2);
601                  
602 <                pot *= 0.0;
602 >                topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
603 >                vdwMult = vdwScale_[topoDist];
604 >                electroMult = electrostaticScale_[topoDist];
605  
568 –                idat = fDecomp_->fillInteractionData(atom1, atom2);
569 –                *(idat.pot) = pot;
570 –
606                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
607 <                  *(idat.d) = d_grp;
608 <                  *(idat.r2) = rgrpsq;
607 >                  idat.d = &d_grp;
608 >                  idat.r2 = &rgrpsq;
609                  } else {
610 <                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
611 <                  curSnapshot->wrapVector( *(idat.d) );
612 <                  *(idat.r2) = idat.d->lengthSquare();
610 >                  d = fDecomp_->getInteratomicVector(atom1, atom2);
611 >                  curSnapshot->wrapVector( d );
612 >                  r2 = d.lengthSquare();
613 >                  idat.d = &d;
614 >                  idat.r2 = &r2;
615                  }
616                  
617 <                *(idat.rij) = sqrt( *(idat.r2) );
617 >                r = sqrt( *(idat.r2) );
618 >                idat.rij = &r;
619                
620                  if (iLoop == PREPAIR_LOOP) {
621                    interactionMan_->doPrePair(idat);
622                  } else {
623                    interactionMan_->doPair(idat);
624                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
625 <                  vij += *(idat.vpair);
626 <                  fij += *(idat.f1);
627 <                  tau -= outProduct( *(idat.d), *(idat.f1));
625 >                  vij += vpair;
626 >                  fij += f1;
627 >                  tau -= outProduct( *(idat.d), f1);
628                  }
629                }
630              }
# Line 652 | Line 690 | namespace OpenMD {
690            fDecomp_->collectIntermediateData();
691  
692            for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
693 <            sdat = fDecomp_->fillSelfData(atom1);
693 >            fDecomp_->fillSelfData(sdat, atom1);
694              interactionMan_->doPreForce(sdat);
695            }
696 <
696 >          
697 >          
698            fDecomp_->distributeIntermediateData();        
699          }
700        }
# Line 668 | Line 707 | namespace OpenMD {
707        
708        for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
709  
710 <        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
710 >        vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 );
711          
712          for (vector<int>::iterator jb = skipList.begin();
713               jb != skipList.end(); ++jb) {        
714      
715            atom2 = (*jb);
716 <          idat = fDecomp_->fillSkipData(atom1, atom2);
716 >          fDecomp_->fillSkipData(idat, atom1, atom2);
717            interactionMan_->doSkipCorrection(idat);
718 +          fDecomp_->unpackSkipData(idat, atom1, atom2);
719  
720          }
721        }
# Line 684 | Line 724 | namespace OpenMD {
724      if (info_->requiresSelfCorrection()) {
725  
726        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
727 <        sdat = fDecomp_->fillSelfData(atom1);
727 >        fDecomp_->fillSelfData(sdat, atom1);
728          interactionMan_->doSelfCorrection(sdat);
729        }
730  
731      }
732  
733 <    longRangePotential = fDecomp_->getLongRangePotential();
733 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
734 >      *(fDecomp_->getPairwisePotential());
735 >
736      lrPot = longRangePotential.sum();
737  
738      //store the tau and long range potential    

Diff Legend

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