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 1583 by gezelter, Thu Jun 16 22:00:08 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 +
134      map<string, CutoffMethod> stringToCutoffMethod;
135      stringToCutoffMethod["HARD"] = HARD;
136      stringToCutoffMethod["SWITCHED"] = SWITCHED;
# Line 194 | Line 201 | namespace OpenMD {
201        simError();
202        cutoffPolicy_ = TRADITIONAL;        
203      }
204 +    fDecomp_->setCutoffPolicy(cutoffPolicy_);
205    }
206  
207    /**
# Line 205 | Line 213 | namespace OpenMD {
213     */
214    void ForceManager::setupSwitching() {
215      Globals* simParams_ = info_->getSimParams();
216 +
217 +    // create the switching function object:
218 +    switcher_ = new SwitchingFunction();
219      
220      if (simParams_->haveSwitchingRadius()) {
221        rSwitch_ = simParams_->getSwitchingRadius();
222        if (rSwitch_ > rCut_) {        
223          sprintf(painCave.errMsg,
224 <                "ForceManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
225 <                rSwitch_, rCut_);
224 >                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
225 >                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
226          painCave.isFatal = 1;
227          painCave.severity = OPENMD_ERROR;
228          simError();
# Line 227 | Line 238 | namespace OpenMD {
238        simError();
239      }          
240      
241 +    // Default to cubic switching function.
242 +    sft_ = cubic;
243      if (simParams_->haveSwitchingFunctionType()) {
244        string funcType = simParams_->getSwitchingFunctionType();
245        toUpper(funcType);
# Line 468 | Line 481 | namespace OpenMD {
481    
482    void ForceManager::longRangeInteractions() {
483  
471 –    // some of this initial stuff will go away:
484      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
485      DataStorage* config = &(curSnapshot->atomData);
486      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;    
487  
488 <    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
489 <      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
488 >    //calculate the center of mass of cutoff group
489 >
490 >    SimInfo::MoleculeIterator mi;
491 >    Molecule* mol;
492 >    Molecule::CutoffGroupIterator ci;
493 >    CutoffGroup* cg;
494 >
495 >    if(info_->getNCutoffGroups() > 0){      
496 >      for (mol = info_->beginMolecule(mi); mol != NULL;
497 >           mol = info_->nextMolecule(mi)) {
498 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
499 >            cg = mol->nextCutoffGroup(ci)) {
500 >          cg->updateCOM();
501 >        }
502 >      }      
503      } else {
504        // center of mass of the group is the same as position of the atom  
505        // if cutoff group does not exist
506 <      rc = pos;
506 >      cgConfig->position = config->position;
507      }
508 <    
491 <    // new stuff starts here:
508 >
509      fDecomp_->zeroWorkArrays();
510      fDecomp_->distributeData();
511 <
512 <    int cg1, cg2, atom1, atom2;
513 <    Vector3d d_grp, dag;
514 <    RealType rgrpsq, rgrp;
511 >    
512 >    int cg1, cg2, atom1, atom2, topoDist;
513 >    Vector3d d_grp, dag, d;
514 >    RealType rgrpsq, rgrp, r2, r;
515 >    RealType electroMult, vdwMult;
516      RealType vij;
517 <    Vector3d fij, fg;
517 >    Vector3d fij, fg, f1;
518      tuple3<RealType, RealType, RealType> cuts;
519      RealType rCutSq;
520      bool in_switching_region;
# Line 505 | Line 523 | namespace OpenMD {
523      InteractionData idat;
524      SelfData sdat;
525      RealType mf;
508 –    potVec pot(0.0);
509 –    potVec longRangePotential(0.0);
526      RealType lrPot;
527 +    RealType vpair;
528 +    potVec longRangePotential(0.0);
529 +    potVec workPot(0.0);
530  
531      int loopStart, loopEnd;
532  
533 +    idat.vdwMult = &vdwMult;
534 +    idat.electroMult = &electroMult;
535 +    idat.pot = &workPot;
536 +    sdat.pot = fDecomp_->getEmbeddingPotential();
537 +    idat.vpair = &vpair;
538 +    idat.f1 = &f1;
539 +    idat.sw = &sw;
540 +    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
541 +    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
542 +    
543      loopEnd = PAIR_LOOP;
544      if (info_->requiresPrepair() ) {
545        loopStart = PREPAIR_LOOP;
546      } else {
547        loopStart = PAIR_LOOP;
548      }
549 <
550 <    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
551 <      
549 >  
550 >    for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
551 >    
552        if (iLoop == loopStart) {
553          bool update_nlist = fDecomp_->checkNeighborList();
554          if (update_nlist)
555            neighborList = fDecomp_->buildNeighborList();
556 <      }
557 <
556 >      }      
557 >        
558        for (vector<pair<int, int> >::iterator it = neighborList.begin();
559               it != neighborList.end(); ++it) {
560 <        
560 >                
561          cg1 = (*it).first;
562          cg2 = (*it).second;
563          
# Line 541 | Line 570 | namespace OpenMD {
570          rCutSq = cuts.second;
571  
572          if (rgrpsq < rCutSq) {
573 <          *(idat.rcut) = cuts.first;
573 >          idat.rcut = &cuts.first;
574            if (iLoop == PAIR_LOOP) {
575              vij *= 0.0;
576              fij = V3Zero;
577            }
578            
579 <          in_switching_region = switcher_->getSwitch(rgrpsq, *(idat.sw), dswdr,
579 >          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
580                                                       rgrp);
581                
582            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
# Line 560 | Line 589 | namespace OpenMD {
589              for (vector<int>::iterator jb = atomListColumn.begin();
590                   jb != atomListColumn.end(); ++jb) {              
591                atom2 = (*jb);
592 <              
592 >            
593                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
594 +                vpair = 0.0;
595 +                workPot = 0.0;
596 +                f1 = V3Zero;
597 +
598 +                fDecomp_->fillInteractionData(idat, atom1, atom2);
599                  
600 <                pot *= 0.0;
600 >                topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
601 >                vdwMult = vdwScale_[topoDist];
602 >                electroMult = electrostaticScale_[topoDist];
603  
568 –                idat = fDecomp_->fillInteractionData(atom1, atom2);
569 –                *(idat.pot) = pot;
570 –
604                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
605 <                  *(idat.d) = d_grp;
606 <                  *(idat.r2) = rgrpsq;
605 >                  idat.d = &d_grp;
606 >                  idat.r2 = &rgrpsq;
607                  } else {
608 <                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
609 <                  curSnapshot->wrapVector( *(idat.d) );
610 <                  *(idat.r2) = idat.d->lengthSquare();
608 >                  d = fDecomp_->getInteratomicVector(atom1, atom2);
609 >                  curSnapshot->wrapVector( d );
610 >                  r2 = d.lengthSquare();
611 >                  idat.d = &d;
612 >                  idat.r2 = &r2;
613                  }
614                  
615 <                *(idat.rij) = sqrt( *(idat.r2) );
615 >                r = sqrt( *(idat.r2) );
616 >                idat.rij = &r;
617                
618                  if (iLoop == PREPAIR_LOOP) {
619                    interactionMan_->doPrePair(idat);
620                  } else {
621                    interactionMan_->doPair(idat);
622                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
623 <                  vij += *(idat.vpair);
624 <                  fij += *(idat.f1);
625 <                  tau -= outProduct( *(idat.d), *(idat.f1));
623 >                  vij += vpair;
624 >                  fij += f1;
625 >                  tau -= outProduct( *(idat.d), f1);
626                  }
627                }
628              }
# Line 652 | Line 688 | namespace OpenMD {
688            fDecomp_->collectIntermediateData();
689  
690            for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
691 <            sdat = fDecomp_->fillSelfData(atom1);
691 >            fDecomp_->fillSelfData(sdat, atom1);
692              interactionMan_->doPreForce(sdat);
693            }
694 <
694 >          
695 >          
696            fDecomp_->distributeIntermediateData();        
697          }
698        }
# Line 668 | Line 705 | namespace OpenMD {
705        
706        for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
707  
708 <        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
708 >        vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 );
709          
710          for (vector<int>::iterator jb = skipList.begin();
711               jb != skipList.end(); ++jb) {        
712      
713            atom2 = (*jb);
714 <          idat = fDecomp_->fillSkipData(atom1, atom2);
714 >          fDecomp_->fillSkipData(idat, atom1, atom2);
715            interactionMan_->doSkipCorrection(idat);
716 +          fDecomp_->unpackSkipData(idat, atom1, atom2);
717  
718          }
719        }
# Line 684 | Line 722 | namespace OpenMD {
722      if (info_->requiresSelfCorrection()) {
723  
724        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
725 <        sdat = fDecomp_->fillSelfData(atom1);
725 >        fDecomp_->fillSelfData(sdat, atom1);
726          interactionMan_->doSelfCorrection(sdat);
727        }
728  
729      }
730  
731 <    longRangePotential = fDecomp_->getLongRangePotential();
731 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
732 >      *(fDecomp_->getPairwisePotential());
733 >
734      lrPot = longRangePotential.sum();
735  
736      //store the tau and long range potential    

Diff Legend

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