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 1593 by gezelter, Fri Jul 15 21:35:14 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;
390      
# 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 492 | Line 526 | namespace OpenMD {
526             mol = info_->nextMolecule(mi)) {
527          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
528              cg = mol->nextCutoffGroup(ci)) {
529 +          cerr << "branch1\n";
530 +          cerr << "globind = " << cg->getGlobalIndex() << "\n";
531            cg->updateCOM();
532          }
533        }      
534      } else {
535        // center of mass of the group is the same as position of the atom  
536        // if cutoff group does not exist
537 +      cerr << "branch2\n";
538        cgConfig->position = config->position;
539      }
540  
# Line 518 | Line 555 | namespace OpenMD {
555      InteractionData idat;
556      SelfData sdat;
557      RealType mf;
521 –    potVec pot(0.0);
522 –    potVec longRangePotential(0.0);
558      RealType lrPot;
559      RealType vpair;
560 +    potVec longRangePotential(0.0);
561 +    potVec workPot(0.0);
562  
563      int loopStart, loopEnd;
564  
565      idat.vdwMult = &vdwMult;
566      idat.electroMult = &electroMult;
567 <    idat.pot = &pot;
567 >    idat.pot = &workPot;
568 >    sdat.pot = fDecomp_->getEmbeddingPotential();
569      idat.vpair = &vpair;
570      idat.f1 = &f1;
571      idat.sw = &sw;
572 <
572 >    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
573 >    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
574 >    
575      loopEnd = PAIR_LOOP;
576      if (info_->requiresPrepair() ) {
577        loopStart = PREPAIR_LOOP;
578      } else {
579        loopStart = PAIR_LOOP;
580      }
581 <    
581 >  
582      for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
583      
584        if (iLoop == loopStart) {
# Line 564 | Line 604 | namespace OpenMD {
604          if (rgrpsq < rCutSq) {
605            idat.rcut = &cuts.first;
606            if (iLoop == PAIR_LOOP) {
607 <            vij *= 0.0;
607 >            vij = 0.0;
608              fij = V3Zero;
609            }
610            
# Line 581 | Line 621 | namespace OpenMD {
621              for (vector<int>::iterator jb = atomListColumn.begin();
622                   jb != atomListColumn.end(); ++jb) {              
623                atom2 = (*jb);
624 <              
624 >
625                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
586 –                
626                  vpair = 0.0;
627 +                workPot = 0.0;
628                  f1 = V3Zero;
629  
630                  fDecomp_->fillInteractionData(idat, atom1, atom2);
# Line 596 | Line 636 | namespace OpenMD {
636                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
637                    idat.d = &d_grp;
638                    idat.r2 = &rgrpsq;
639 +                  cerr << "dgrp = " << d_grp << "\n";
640                  } else {
641                    d = fDecomp_->getInteratomicVector(atom1, atom2);
642                    curSnapshot->wrapVector( d );
643                    r2 = d.lengthSquare();
644 +                  cerr << "datm = " << d<< "\n";
645                    idat.d = &d;
646                    idat.r2 = &r2;
647                  }
648                  
649 +                cerr << "idat.d = " << *(idat.d) << "\n";
650                  r = sqrt( *(idat.r2) );
651                  idat.rij = &r;
652                
# Line 612 | Line 655 | namespace OpenMD {
655                  } else {
656                    interactionMan_->doPair(idat);
657                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
658 +
659 +                  cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << "\n";
660                    vij += vpair;
661                    fij += f1;
662                    tau -= outProduct( *(idat.d), f1);
# Line 624 | Line 669 | namespace OpenMD {
669              if (in_switching_region) {
670                swderiv = vij * dswdr / rgrp;
671                fg = swderiv * d_grp;
627 –
672                fij += fg;
673  
674                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
# Line 676 | Line 720 | namespace OpenMD {
720        }
721  
722        if (iLoop == PREPAIR_LOOP) {
723 <        if (info_->requiresPrepair()) {            
723 >        if (info_->requiresPrepair()) {
724 >
725            fDecomp_->collectIntermediateData();
726  
727            for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
# Line 684 | Line 729 | namespace OpenMD {
729              interactionMan_->doPreForce(sdat);
730            }
731  
732 <          fDecomp_->distributeIntermediateData();        
732 >          fDecomp_->distributeIntermediateData();
733 >
734          }
735        }
736  
737      }
738      
739      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 );
740          
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 –    
741      if (info_->requiresSelfCorrection()) {
742  
743        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
# Line 718 | Line 747 | namespace OpenMD {
747  
748      }
749  
750 <    longRangePotential = fDecomp_->getLongRangePotential();
750 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
751 >      *(fDecomp_->getPairwisePotential());
752 >
753      lrPot = longRangePotential.sum();
754  
755      //store the tau and long range potential    

Diff Legend

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