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 1583 by gezelter, Thu Jun 16 22:00:08 2011 UTC vs.
Revision 1592 by gezelter, Tue Jul 12 20:33:14 2011 UTC

# Line 75 | 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 85 | Line 86 | namespace OpenMD {
86     *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
87     *      Use the maximum suggested value that was found.
88     *
89 <   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
89 >   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE,
90 >   *                        or SHIFTED_POTENTIAL)
91     *      If cutoffMethod was explicitly set, use that choice.
92     *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
93     *
94     * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
95     *      If cutoffPolicy was explicitly set, use that choice.
96     *      If cutoffPolicy was not explicitly set, use TRADITIONAL
97 +   *
98 +   * switchingRadius : realType
99 +   *  If the cutoffMethod was set to SWITCHED:
100 +   *      If the switchingRadius was explicitly set, use that value
101 +   *          (but do a sanity check first).
102 +   *      If the switchingRadius was not explicitly set: use 0.85 *
103 +   *      cutoffRadius_
104 +   *  If the cutoffMethod was not set to SWITCHED:
105 +   *      Set switchingRadius equal to cutoffRadius for safety.
106     */
107    void ForceManager::setupCutoffs() {
108      
# Line 130 | Line 141 | namespace OpenMD {
141      }
142  
143      fDecomp_->setUserCutoff(rCut_);
144 +    interactionMan_->setCutoffRadius(rCut_);
145  
146      map<string, CutoffMethod> stringToCutoffMethod;
147      stringToCutoffMethod["HARD"] = HARD;
# Line 201 | Line 213 | namespace OpenMD {
213        simError();
214        cutoffPolicy_ = TRADITIONAL;        
215      }
204 –    fDecomp_->setCutoffPolicy(cutoffPolicy_);
205 –  }
216  
217 <  /**
218 <   * setupSwitching
209 <   *
210 <   * Sets the values of switchingRadius and
211 <   *  If the switchingRadius was explicitly set, use that value (but check it)
212 <   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
213 <   */
214 <  void ForceManager::setupSwitching() {
215 <    Globals* simParams_ = info_->getSimParams();
216 <
217 >    fDecomp_->setCutoffPolicy(cutoffPolicy_);
218 >        
219      // create the switching function object:
220 +
221      switcher_ = new SwitchingFunction();
222 <    
223 <    if (simParams_->haveSwitchingRadius()) {
224 <      rSwitch_ = simParams_->getSwitchingRadius();
225 <      if (rSwitch_ > rCut_) {        
222 >  
223 >    if (cutoffMethod_ == SWITCHED) {
224 >      if (simParams_->haveSwitchingRadius()) {
225 >        rSwitch_ = simParams_->getSwitchingRadius();
226 >        if (rSwitch_ > rCut_) {        
227 >          sprintf(painCave.errMsg,
228 >                  "ForceManager::setupCutoffs: switchingRadius (%f) is larger "
229 >                  "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
230 >          painCave.isFatal = 1;
231 >          painCave.severity = OPENMD_ERROR;
232 >          simError();
233 >        }
234 >      } else {      
235 >        rSwitch_ = 0.85 * rCut_;
236          sprintf(painCave.errMsg,
237 <                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
238 <                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
239 <        painCave.isFatal = 1;
240 <        painCave.severity = OPENMD_ERROR;
237 >                "ForceManager::setupCutoffs: No value was set for the switchingRadius.\n"
238 >                "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
239 >                "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
240 >        painCave.isFatal = 0;
241 >        painCave.severity = OPENMD_WARNING;
242          simError();
243        }
244 <    } else {      
245 <      rSwitch_ = 0.85 * rCut_;
246 <      sprintf(painCave.errMsg,
247 <              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
248 <              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
249 <              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
250 <      painCave.isFatal = 0;
251 <      painCave.severity = OPENMD_WARNING;
252 <      simError();
253 <    }          
244 >    } else {
245 >      if (simParams_->haveSwitchingRadius()) {
246 >        map<string, CutoffMethod>::const_iterator it;
247 >        string theMeth;
248 >        for (it = stringToCutoffMethod.begin();
249 >             it != stringToCutoffMethod.end(); ++it) {
250 >          if (it->second == cutoffMethod_) {
251 >            theMeth = it->first;
252 >            break;
253 >          }
254 >        }
255 >        sprintf(painCave.errMsg,
256 >                "ForceManager::setupCutoffs: the cutoffMethod (%s)\n"
257 >                "\tis not set to SWITCHED, so switchingRadius value\n"
258 >                "\twill be ignored for this simulation\n", theMeth.c_str());
259 >        painCave.isFatal = 0;
260 >        painCave.severity = OPENMD_WARNING;
261 >        simError();
262 >      }
263 >
264 >      rSwitch_ = rCut_;
265 >    }
266      
267      // Default to cubic switching function.
268      sft_ = cubic;
# Line 263 | Line 289 | namespace OpenMD {
289      }
290      switcher_->setSwitchType(sft_);
291      switcher_->setSwitch(rSwitch_, rCut_);
292 +    interactionMan_->setSwitchingRadius(rSwitch_);
293    }
294    
295    void ForceManager::initialize() {
296  
297      if (!info_->isTopologyDone()) {
298 +
299        info_->update();
300        interactionMan_->setSimInfo(info_);
301        interactionMan_->initialize();
# Line 275 | Line 303 | namespace OpenMD {
303        // We want to delay the cutoffs until after the interaction
304        // manager has set up the atom-atom interactions so that we can
305        // query them for suggested cutoff values
278 –
306        setupCutoffs();
280 –      setupSwitching();
307  
308        info_->prepareTopology();      
309      }
310  
311      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
312      
313 <    // Force fields can set options on how to scale van der Waals and electrostatic
314 <    // interactions for atoms connected via bonds, bends and torsions
315 <    // in this case the topological distance between atoms is:
313 >    // Force fields can set options on how to scale van der Waals and
314 >    // electrostatic interactions for atoms connected via bonds, bends
315 >    // and torsions in this case the topological distance between
316 >    // atoms is:
317      // 0 = topologically unconnected
318      // 1 = bonded together
319      // 2 = connected via a bend
# Line 338 | Line 365 | namespace OpenMD {
365      
366      for (mol = info_->beginMolecule(mi); mol != NULL;
367           mol = info_->nextMolecule(mi)) {
368 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
368 >      for(atom = mol->beginAtom(ai); atom != NULL;
369 >          atom = mol->nextAtom(ai)) {
370          atom->zeroForcesAndTorques();
371        }
372 <          
372 >      
373        //change the positions of atoms which belong to the rigidbodies
374        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
375             rb = mol->nextRigidBody(rbIter)) {
376          rb->zeroForcesAndTorques();
377        }        
378 <
378 >      
379        if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
380          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
381              cg = mol->nextCutoffGroup(ci)) {
# Line 356 | Line 384 | namespace OpenMD {
384          }
385        }      
386      }
387 <  
387 >    
388      // Zero out the stress tensor
389      tau *= 0.0;
390      
# Line 410 | Line 438 | namespace OpenMD {
438            dataSet.prev.angle = dataSet.curr.angle = angle;
439            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
440            dataSet.deltaV = 0.0;
441 <          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
441 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend,
442 >                                                                  dataSet));
443          }else {
444            i->second.prev.angle = i->second.curr.angle;
445            i->second.prev.potential = i->second.curr.potential;
# Line 572 | Line 601 | namespace OpenMD {
601          if (rgrpsq < rCutSq) {
602            idat.rcut = &cuts.first;
603            if (iLoop == PAIR_LOOP) {
604 <            vij *= 0.0;
604 >            vij = 0.0;
605              fij = V3Zero;
606            }
607            
# Line 622 | Line 651 | namespace OpenMD {
651                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
652                    vij += vpair;
653                    fij += f1;
654 +                  cerr << "ats = " << atom1 << " " << atom2 << "d = " << *(idat.d) << " f1 = " << f1 << "\n";
655                    tau -= outProduct( *(idat.d), f1);
656                  }
657                }
# Line 632 | Line 662 | namespace OpenMD {
662              if (in_switching_region) {
663                swderiv = vij * dswdr / rgrp;
664                fg = swderiv * d_grp;
635 –
665                fij += fg;
666  
667                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
# Line 684 | Line 713 | namespace OpenMD {
713        }
714  
715        if (iLoop == PREPAIR_LOOP) {
716 <        if (info_->requiresPrepair()) {            
716 >        if (info_->requiresPrepair()) {
717 >
718            fDecomp_->collectIntermediateData();
719  
720            for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
721              fDecomp_->fillSelfData(sdat, atom1);
722              interactionMan_->doPreForce(sdat);
723            }
724 <          
725 <          
726 <          fDecomp_->distributeIntermediateData();        
724 >
725 >          fDecomp_->distributeIntermediateData();
726 >
727          }
728        }
729  
730      }
731      
732      fDecomp_->collectData();
703 –    
704 –    if ( info_->requiresSkipCorrection() ) {
705 –      
706 –      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
707 –
708 –        vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 );
733          
710 –        for (vector<int>::iterator jb = skipList.begin();
711 –             jb != skipList.end(); ++jb) {        
712 –    
713 –          atom2 = (*jb);
714 –          fDecomp_->fillSkipData(idat, atom1, atom2);
715 –          interactionMan_->doSkipCorrection(idat);
716 –          fDecomp_->unpackSkipData(idat, atom1, atom2);
717 –
718 –        }
719 –      }
720 –    }
721 –    
734      if (info_->requiresSelfCorrection()) {
735  
736        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          

Diff Legend

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