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 1581 by gezelter, Mon Jun 13 22:13:12 2011 UTC vs.
Revision 1587 by gezelter, Fri Jul 8 20:25:32 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    
# Line 71 | Line 75 | namespace OpenMD {
75    /**
76     * setupCutoffs
77     *
78 <   * Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
78 >   * Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
79 >   * and cutoffPolicy
80     *
81     * cutoffRadius : realType
82     *  If the cutoffRadius was explicitly set, use that value.
# Line 88 | Line 93 | namespace OpenMD {
93     * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
94     *      If cutoffPolicy was explicitly set, use that choice.
95     *      If cutoffPolicy was not explicitly set, use TRADITIONAL
96 +   *
97 +   * switchingRadius : realType
98 +   *  If the cutoffMethod was set to SWITCHED:
99 +   *      If the switchingRadius was explicitly set, use that value
100 +   *          (but do a sanity check first).
101 +   *      If the switchingRadius was not explicitly set: use 0.85 *
102 +   *      cutoffRadius_
103 +   *  If the cutoffMethod was not set to SWITCHED:
104 +   *      Set switchingRadius equal to cutoffRadius for safety.
105     */
106    void ForceManager::setupCutoffs() {
107      
# Line 123 | Line 137 | namespace OpenMD {
137          painCave.severity = OPENMD_INFO;
138          simError();
139        }
126 –      fDecomp_->setUserCutoff(rCut_);
140      }
141  
142 +    fDecomp_->setUserCutoff(rCut_);
143 +    interactionMan_->setCutoffRadius(rCut_);
144 +
145      map<string, CutoffMethod> stringToCutoffMethod;
146      stringToCutoffMethod["HARD"] = HARD;
147      stringToCutoffMethod["SWITCHED"] = SWITCHED;
# Line 196 | Line 212 | namespace OpenMD {
212        simError();
213        cutoffPolicy_ = TRADITIONAL;        
214      }
199 –    fDecomp_->setCutoffPolicy(cutoffPolicy_);
200 –  }
215  
216 <  /**
217 <   * setupSwitching
204 <   *
205 <   * Sets the values of switchingRadius and
206 <   *  If the switchingRadius was explicitly set, use that value (but check it)
207 <   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
208 <   */
209 <  void ForceManager::setupSwitching() {
210 <    Globals* simParams_ = info_->getSimParams();
211 <
216 >    fDecomp_->setCutoffPolicy(cutoffPolicy_);
217 >        
218      // create the switching function object:
219 +
220      switcher_ = new SwitchingFunction();
221 <    
222 <    if (simParams_->haveSwitchingRadius()) {
223 <      rSwitch_ = simParams_->getSwitchingRadius();
224 <      if (rSwitch_ > rCut_) {        
221 >  
222 >    if (cutoffMethod_ == SWITCHED) {
223 >      if (simParams_->haveSwitchingRadius()) {
224 >        rSwitch_ = simParams_->getSwitchingRadius();
225 >        if (rSwitch_ > rCut_) {        
226 >          sprintf(painCave.errMsg,
227 >                  "ForceManager::setupCutoffs: switchingRadius (%f) is larger "
228 >                  "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
229 >          painCave.isFatal = 1;
230 >          painCave.severity = OPENMD_ERROR;
231 >          simError();
232 >        }
233 >      } else {      
234 >        rSwitch_ = 0.85 * rCut_;
235          sprintf(painCave.errMsg,
236 <                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
237 <                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
238 <        painCave.isFatal = 1;
239 <        painCave.severity = OPENMD_ERROR;
236 >                "ForceManager::setupCutoffs: No value was set for the switchingRadius.\n"
237 >                "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
238 >                "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
239 >        painCave.isFatal = 0;
240 >        painCave.severity = OPENMD_WARNING;
241          simError();
242        }
243 <    } else {      
244 <      rSwitch_ = 0.85 * rCut_;
245 <      sprintf(painCave.errMsg,
246 <              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
247 <              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
248 <              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
249 <      painCave.isFatal = 0;
250 <      painCave.severity = OPENMD_WARNING;
251 <      simError();
252 <    }          
243 >    } else {
244 >      if (simParams_->haveSwitchingRadius()) {
245 >        map<string, CutoffMethod>::const_iterator it;
246 >        string theMeth;
247 >        for (it = stringToCutoffMethod.begin();
248 >             it != stringToCutoffMethod.end(); ++it) {
249 >          if (it->second == cutoffMethod_) {
250 >            theMeth = it->first;
251 >            break;
252 >          }
253 >        }
254 >        sprintf(painCave.errMsg,
255 >                "ForceManager::setupCutoffs: the cutoffMethod (%s)\n"
256 >                "\tis not set to SWITCHED, so switchingRadius value\n"
257 >                "\twill be ignored for this simulation\n", theMeth.c_str());
258 >        painCave.isFatal = 0;
259 >        painCave.severity = OPENMD_WARNING;
260 >        simError();
261 >      }
262 >
263 >      rSwitch_ = rCut_;
264 >    }
265      
266      // Default to cubic switching function.
267      sft_ = cubic;
# Line 258 | Line 288 | namespace OpenMD {
288      }
289      switcher_->setSwitchType(sft_);
290      switcher_->setSwitch(rSwitch_, rCut_);
291 +    interactionMan_->setSwitchingRadius(rSwitch_);
292    }
293    
294    void ForceManager::initialize() {
# Line 272 | Line 303 | namespace OpenMD {
303        // query them for suggested cutoff values
304  
305        setupCutoffs();
275 –      setupSwitching();
306  
307        info_->prepareTopology();      
308      }
# Line 518 | Line 548 | namespace OpenMD {
548      InteractionData idat;
549      SelfData sdat;
550      RealType mf;
521 –    potVec pot(0.0);
522 –    potVec longRangePotential(0.0);
551      RealType lrPot;
552      RealType vpair;
553 +    potVec longRangePotential(0.0);
554 +    potVec workPot(0.0);
555  
556      int loopStart, loopEnd;
557  
558      idat.vdwMult = &vdwMult;
559      idat.electroMult = &electroMult;
560 <    idat.pot = &pot;
560 >    idat.pot = &workPot;
561 >    sdat.pot = fDecomp_->getEmbeddingPotential();
562      idat.vpair = &vpair;
563      idat.f1 = &f1;
564      idat.sw = &sw;
565 <
565 >    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
566 >    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
567 >    
568      loopEnd = PAIR_LOOP;
569      if (info_->requiresPrepair() ) {
570        loopStart = PREPAIR_LOOP;
571      } else {
572        loopStart = PAIR_LOOP;
573      }
574 <    
574 >  
575      for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
576      
577        if (iLoop == loopStart) {
# Line 564 | Line 597 | namespace OpenMD {
597          if (rgrpsq < rCutSq) {
598            idat.rcut = &cuts.first;
599            if (iLoop == PAIR_LOOP) {
600 <            vij *= 0.0;
600 >            vij = 0.0;
601              fij = V3Zero;
602            }
603            
# Line 581 | Line 614 | namespace OpenMD {
614              for (vector<int>::iterator jb = atomListColumn.begin();
615                   jb != atomListColumn.end(); ++jb) {              
616                atom2 = (*jb);
617 <              
617 >            
618                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
586 –                
619                  vpair = 0.0;
620 +                workPot = 0.0;
621                  f1 = V3Zero;
622  
623                  fDecomp_->fillInteractionData(idat, atom1, atom2);
# Line 624 | Line 657 | namespace OpenMD {
657              if (in_switching_region) {
658                swderiv = vij * dswdr / rgrp;
659                fg = swderiv * d_grp;
627 –
660                fij += fg;
661  
662                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
# Line 683 | Line 715 | namespace OpenMD {
715              fDecomp_->fillSelfData(sdat, atom1);
716              interactionMan_->doPreForce(sdat);
717            }
718 <
718 >          
719 >          
720            fDecomp_->distributeIntermediateData();        
721          }
722        }
# Line 691 | Line 724 | namespace OpenMD {
724      }
725      
726      fDecomp_->collectData();
694 –    
695 –    if ( info_->requiresSkipCorrection() ) {
696 –      
697 –      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
698 –
699 –        vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 );
727          
701 –        for (vector<int>::iterator jb = skipList.begin();
702 –             jb != skipList.end(); ++jb) {        
703 –    
704 –          atom2 = (*jb);
705 –          fDecomp_->fillSkipData(idat, atom1, atom2);
706 –          interactionMan_->doSkipCorrection(idat);
707 –
708 –        }
709 –      }
710 –    }
711 –    
728      if (info_->requiresSelfCorrection()) {
729  
730        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
# Line 718 | Line 734 | namespace OpenMD {
734  
735      }
736  
737 <    longRangePotential = fDecomp_->getLongRangePotential();
737 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
738 >      *(fDecomp_->getPairwisePotential());
739 >
740      lrPot = longRangePotential.sum();
741  
742      //store the tau and long range potential    

Diff Legend

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