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 1590 by gezelter, Mon Jul 11 01:39:49 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 81 | 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 123 | Line 138 | namespace OpenMD {
138          painCave.severity = OPENMD_INFO;
139          simError();
140        }
126 –      fDecomp_->setUserCutoff(rCut_);
141      }
142  
143 +    fDecomp_->setUserCutoff(rCut_);
144 +    interactionMan_->setCutoffRadius(rCut_);
145 +
146      map<string, CutoffMethod> stringToCutoffMethod;
147      stringToCutoffMethod["HARD"] = HARD;
148      stringToCutoffMethod["SWITCHED"] = SWITCHED;
# Line 196 | Line 213 | namespace OpenMD {
213        simError();
214        cutoffPolicy_ = TRADITIONAL;        
215      }
199 –    fDecomp_->setCutoffPolicy(cutoffPolicy_);
200 –  }
216  
217 <  /**
218 <   * 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 <
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 258 | 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 270 | 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
273 –
306        setupCutoffs();
275 –      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 333 | 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 351 | Line 384 | namespace OpenMD {
384          }
385        }      
386      }
387 <  
387 >    
388      // Zero out the stress tensor
389 <    tau *= 0.0;
389 >    tau = Mat3x3d(0.0);
390      
391    }
392    
# Line 405 | 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 518 | Line 552 | namespace OpenMD {
552      InteractionData idat;
553      SelfData sdat;
554      RealType mf;
521 –    potVec pot(0.0);
522 –    potVec longRangePotential(0.0);
555      RealType lrPot;
556      RealType vpair;
557 +    potVec longRangePotential(0.0);
558 +    potVec workPot(0.0);
559  
560      int loopStart, loopEnd;
561  
562      idat.vdwMult = &vdwMult;
563      idat.electroMult = &electroMult;
564 <    idat.pot = &pot;
564 >    idat.pot = &workPot;
565 >    sdat.pot = fDecomp_->getEmbeddingPotential();
566      idat.vpair = &vpair;
567      idat.f1 = &f1;
568      idat.sw = &sw;
569 <
569 >    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
570 >    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
571 >    
572      loopEnd = PAIR_LOOP;
573      if (info_->requiresPrepair() ) {
574        loopStart = PREPAIR_LOOP;
575      } else {
576        loopStart = PAIR_LOOP;
577      }
578 <    
578 >  
579      for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
580      
581        if (iLoop == loopStart) {
# Line 564 | 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 581 | Line 618 | namespace OpenMD {
618              for (vector<int>::iterator jb = atomListColumn.begin();
619                   jb != atomListColumn.end(); ++jb) {              
620                atom2 = (*jb);
621 <              
621 >            
622                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
586 –                
623                  vpair = 0.0;
624 +                workPot = 0.0;
625                  f1 = V3Zero;
626  
627                  fDecomp_->fillInteractionData(idat, atom1, atom2);
# Line 624 | Line 661 | namespace OpenMD {
661              if (in_switching_region) {
662                swderiv = vij * dswdr / rgrp;
663                fg = swderiv * d_grp;
627 –
664                fij += fg;
665  
666                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
# Line 676 | Line 712 | namespace OpenMD {
712        }
713  
714        if (iLoop == PREPAIR_LOOP) {
715 <        if (info_->requiresPrepair()) {            
715 >        if (info_->requiresPrepair()) {
716 >
717            fDecomp_->collectIntermediateData();
718  
719            for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
# Line 684 | Line 721 | namespace OpenMD {
721              interactionMan_->doPreForce(sdat);
722            }
723  
724 <          fDecomp_->distributeIntermediateData();        
724 >          fDecomp_->distributeIntermediateData();
725 >
726          }
727        }
728  
729      }
730      
731      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 );
732          
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 –    
733      if (info_->requiresSelfCorrection()) {
734  
735        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
# Line 718 | Line 739 | namespace OpenMD {
739  
740      }
741  
742 <    longRangePotential = fDecomp_->getLongRangePotential();
742 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
743 >      *(fDecomp_->getPairwisePotential());
744 >
745      lrPot = longRangePotential.sum();
746  
747      //store the tau and long range potential    

Diff Legend

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