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 1569 by gezelter, Thu May 26 13:55:04 2011 UTC vs.
Revision 1577 by gezelter, Wed Jun 8 20:26:56 2011 UTC

# Line 47 | Line 47
47   * @version 1.0
48   */
49  
50 +
51   #include "brains/ForceManager.hpp"
52   #include "primitives/Molecule.hpp"
53   #define __OPENMD_C
# Line 62 | Line 63 | namespace OpenMD {
63   namespace OpenMD {
64    
65    ForceManager::ForceManager(SimInfo * info) : info_(info) {
66 <
66 < #ifdef IS_MPI
66 >    forceField_ = info_->getForceField();
67      fDecomp_ = new ForceMatrixDecomposition(info_);
68 < #else
69 <    // fDecomp_ = new ForceSerialDecomposition(info);
70 < #endif
68 >    interactionMan_ = new InteractionManager();
69    }
70 +
71 +  /**
72 +   * setupCutoffs
73 +   *
74 +   * Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
75 +   *
76 +   * cutoffRadius : realType
77 +   *  If the cutoffRadius was explicitly set, use that value.
78 +   *  If the cutoffRadius was not explicitly set:
79 +   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
80 +   *      No electrostatic atoms?  Poll the atom types present in the
81 +   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
82 +   *      Use the maximum suggested value that was found.
83 +   *
84 +   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
85 +   *      If cutoffMethod was explicitly set, use that choice.
86 +   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
87 +   *
88 +   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
89 +   *      If cutoffPolicy was explicitly set, use that choice.
90 +   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
91 +   */
92 +  void ForceManager::setupCutoffs() {
93 +    
94 +    Globals* simParams_ = info_->getSimParams();
95 +    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
96 +    
97 +    if (simParams_->haveCutoffRadius()) {
98 +      rCut_ = simParams_->getCutoffRadius();
99 +    } else {      
100 +      if (info_->usesElectrostaticAtoms()) {
101 +        sprintf(painCave.errMsg,
102 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
103 +                "\tOpenMD will use a default value of 12.0 angstroms"
104 +                "\tfor the cutoffRadius.\n");
105 +        painCave.isFatal = 0;
106 +        painCave.severity = OPENMD_INFO;
107 +        simError();
108 +        rCut_ = 12.0;
109 +      } else {
110 +        RealType thisCut;
111 +        set<AtomType*>::iterator i;
112 +        set<AtomType*> atomTypes;
113 +        atomTypes = info_->getSimulatedAtomTypes();        
114 +        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
115 +          thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
116 +          rCut_ = max(thisCut, rCut_);
117 +        }
118 +        sprintf(painCave.errMsg,
119 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
120 +                "\tOpenMD will use %lf angstroms.\n",
121 +                rCut_);
122 +        painCave.isFatal = 0;
123 +        painCave.severity = OPENMD_INFO;
124 +        simError();
125 +      }            
126 +    }
127 +
128 +    map<string, CutoffMethod> stringToCutoffMethod;
129 +    stringToCutoffMethod["HARD"] = HARD;
130 +    stringToCutoffMethod["SWITCHED"] = SWITCHED;
131 +    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
132 +    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
133    
134 <  void ForceManager::calcForces() {
134 >    if (simParams_->haveCutoffMethod()) {
135 >      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
136 >      map<string, CutoffMethod>::iterator i;
137 >      i = stringToCutoffMethod.find(cutMeth);
138 >      if (i == stringToCutoffMethod.end()) {
139 >        sprintf(painCave.errMsg,
140 >                "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
141 >                "\tShould be one of: "
142 >                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
143 >                cutMeth.c_str());
144 >        painCave.isFatal = 1;
145 >        painCave.severity = OPENMD_ERROR;
146 >        simError();
147 >      } else {
148 >        cutoffMethod_ = i->second;
149 >      }
150 >    } else {
151 >      sprintf(painCave.errMsg,
152 >              "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
153 >              "\tOpenMD will use SHIFTED_FORCE.\n");
154 >      painCave.isFatal = 0;
155 >      painCave.severity = OPENMD_INFO;
156 >      simError();
157 >      cutoffMethod_ = SHIFTED_FORCE;        
158 >    }
159 >
160 >    map<string, CutoffPolicy> stringToCutoffPolicy;
161 >    stringToCutoffPolicy["MIX"] = MIX;
162 >    stringToCutoffPolicy["MAX"] = MAX;
163 >    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
164 >
165 >    std::string cutPolicy;
166 >    if (forceFieldOptions_.haveCutoffPolicy()){
167 >      cutPolicy = forceFieldOptions_.getCutoffPolicy();
168 >    }else if (simParams_->haveCutoffPolicy()) {
169 >      cutPolicy = simParams_->getCutoffPolicy();
170 >    }
171 >
172 >    if (!cutPolicy.empty()){
173 >      toUpper(cutPolicy);
174 >      map<string, CutoffPolicy>::iterator i;
175 >      i = stringToCutoffPolicy.find(cutPolicy);
176 >
177 >      if (i == stringToCutoffPolicy.end()) {
178 >        sprintf(painCave.errMsg,
179 >                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
180 >                "\tShould be one of: "
181 >                "MIX, MAX, or TRADITIONAL\n",
182 >                cutPolicy.c_str());
183 >        painCave.isFatal = 1;
184 >        painCave.severity = OPENMD_ERROR;
185 >        simError();
186 >      } else {
187 >        cutoffPolicy_ = i->second;
188 >      }
189 >    } else {
190 >      sprintf(painCave.errMsg,
191 >              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
192 >              "\tOpenMD will use TRADITIONAL.\n");
193 >      painCave.isFatal = 0;
194 >      painCave.severity = OPENMD_INFO;
195 >      simError();
196 >      cutoffPolicy_ = TRADITIONAL;        
197 >    }
198 >  }
199 >
200 >  /**
201 >   * setupSwitching
202 >   *
203 >   * Sets the values of switchingRadius and
204 >   *  If the switchingRadius was explicitly set, use that value (but check it)
205 >   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
206 >   */
207 >  void ForceManager::setupSwitching() {
208 >    Globals* simParams_ = info_->getSimParams();
209 >
210 >    // create the switching function object:
211 >    switcher_ = new SwitchingFunction();
212      
213 +    if (simParams_->haveSwitchingRadius()) {
214 +      rSwitch_ = simParams_->getSwitchingRadius();
215 +      if (rSwitch_ > rCut_) {        
216 +        sprintf(painCave.errMsg,
217 +                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
218 +                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
219 +        painCave.isFatal = 1;
220 +        painCave.severity = OPENMD_ERROR;
221 +        simError();
222 +      }
223 +    } else {      
224 +      rSwitch_ = 0.85 * rCut_;
225 +      sprintf(painCave.errMsg,
226 +              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
227 +              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
228 +              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
229 +      painCave.isFatal = 0;
230 +      painCave.severity = OPENMD_WARNING;
231 +      simError();
232 +    }          
233 +    
234 +    // Default to cubic switching function.
235 +    sft_ = cubic;
236 +    if (simParams_->haveSwitchingFunctionType()) {
237 +      string funcType = simParams_->getSwitchingFunctionType();
238 +      toUpper(funcType);
239 +      if (funcType == "CUBIC") {
240 +        sft_ = cubic;
241 +      } else {
242 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
243 +          sft_ = fifth_order_poly;
244 +        } else {
245 +          // throw error        
246 +          sprintf( painCave.errMsg,
247 +                   "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
248 +                   "\tswitchingFunctionType must be one of: "
249 +                   "\"cubic\" or \"fifth_order_polynomial\".",
250 +                   funcType.c_str() );
251 +          painCave.isFatal = 1;
252 +          painCave.severity = OPENMD_ERROR;
253 +          simError();
254 +        }          
255 +      }
256 +    }
257 +    switcher_->setSwitchType(sft_);
258 +    switcher_->setSwitch(rSwitch_, rCut_);
259 +  }
260 +  
261 +  void ForceManager::initialize() {
262 +
263      if (!info_->isTopologyDone()) {
264        info_->update();
265        interactionMan_->setSimInfo(info_);
266        interactionMan_->initialize();
267 <      swfun_ = interactionMan_->getSwitchingFunction();
268 <      fDecomp_->distributeInitialData();
269 <      info_->prepareTopology();
267 >
268 >      // We want to delay the cutoffs until after the interaction
269 >      // manager has set up the atom-atom interactions so that we can
270 >      // query them for suggested cutoff values
271 >
272 >      setupCutoffs();
273 >      setupSwitching();
274 >
275 >      info_->prepareTopology();      
276      }
277 +
278 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
279      
280 +    // Force fields can set options on how to scale van der Waals and electrostatic
281 +    // interactions for atoms connected via bonds, bends and torsions
282 +    // in this case the topological distance between atoms is:
283 +    // 0 = topologically unconnected
284 +    // 1 = bonded together
285 +    // 2 = connected via a bend
286 +    // 3 = connected via a torsion
287 +    
288 +    vdwScale_.reserve(4);
289 +    fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
290 +
291 +    electrostaticScale_.reserve(4);
292 +    fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
293 +
294 +    vdwScale_[0] = 1.0;
295 +    vdwScale_[1] = fopts.getvdw12scale();
296 +    vdwScale_[2] = fopts.getvdw13scale();
297 +    vdwScale_[3] = fopts.getvdw14scale();
298 +    
299 +    electrostaticScale_[0] = 1.0;
300 +    electrostaticScale_[1] = fopts.getelectrostatic12scale();
301 +    electrostaticScale_[2] = fopts.getelectrostatic13scale();
302 +    electrostaticScale_[3] = fopts.getelectrostatic14scale();    
303 +    
304 +    fDecomp_->distributeInitialData();
305 +
306 +    initialized_ = true;
307 +
308 +  }
309 +
310 +  void ForceManager::calcForces() {
311 +    
312 +    if (!initialized_) initialize();
313 +
314      preCalculation();  
315      shortRangeInteractions();
316      longRangeInteractions();
317 <    postCalculation();
88 <    
317 >    postCalculation();    
318    }
319    
320    void ForceManager::preCalculation() {
# Line 265 | Line 494 | namespace OpenMD {
494        rc = pos;
495      }
496      
268 –    //initialize data before passing to fortran
269 –    RealType longRangePotential[N_INTERACTION_FAMILIES];
270 –    RealType lrPot = 0.0;
271 –    int isError = 0;
272 –
273 –    // dangerous to iterate over enums, but we'll live on the edge:
274 –    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
275 –      longRangePotential[i]=0.0; //Initialize array
276 –    }
277 –
497      // new stuff starts here:
498 <
498 >    fDecomp_->zeroWorkArrays();
499      fDecomp_->distributeData();
500  
501      int cg1, cg2, atom1, atom2;
# Line 284 | Line 503 | namespace OpenMD {
503      RealType rgrpsq, rgrp;
504      RealType vij;
505      Vector3d fij, fg;
506 <    pair<int, int> gtypes;
506 >    tuple3<RealType, RealType, RealType> cuts;
507      RealType rCutSq;
508      bool in_switching_region;
509      RealType sw, dswdr, swderiv;
# Line 292 | Line 511 | namespace OpenMD {
511      InteractionData idat;
512      SelfData sdat;
513      RealType mf;
514 +    potVec pot(0.0);
515 +    potVec longRangePotential(0.0);
516 +    RealType lrPot;
517  
518      int loopStart, loopEnd;
519  
# Line 315 | Line 537 | namespace OpenMD {
537          
538          cg1 = (*it).first;
539          cg2 = (*it).second;
540 +        
541 +        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
542  
319 –        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
543          d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
544          curSnapshot->wrapVector(d_grp);        
545          rgrpsq = d_grp.lengthSquare();
323 –        rCutSq = groupCutoffMap[gtypes].first;
546  
547 +        rCutSq = cuts.second;
548 +
549          if (rgrpsq < rCutSq) {
550 <          *(idat.rcut) = groupCutoffMap[gtypes].second;
550 >          *(idat.rcut) = cuts.first;
551            if (iLoop == PAIR_LOOP) {
552              vij *= 0.0;
553              fij = V3Zero;
554            }
555            
556 <          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
557 <                                                  rgrp);              
556 >          in_switching_region = switcher_->getSwitch(rgrpsq, *(idat.sw), dswdr,
557 >                                                     rgrp);
558 >              
559            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
560            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
561  
# Line 344 | Line 569 | namespace OpenMD {
569                
570                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
571                  
572 +                pot *= 0.0;
573 +
574                  idat = fDecomp_->fillInteractionData(atom1, atom2);
575 +                *(idat.pot) = pot;
576  
577                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
578                    *(idat.d) = d_grp;
# Line 361 | Line 589 | namespace OpenMD {
589                    interactionMan_->doPrePair(idat);
590                  } else {
591                    interactionMan_->doPair(idat);
592 +                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
593                    vij += *(idat.vpair);
594                    fij += *(idat.f1);
595                    tau -= outProduct( *(idat.d), *(idat.f1));
# Line 427 | Line 656 | namespace OpenMD {
656        if (iLoop == PREPAIR_LOOP) {
657          if (info_->requiresPrepair()) {            
658            fDecomp_->collectIntermediateData();
659 <          atomListLocal = fDecomp_->getAtomList();
660 <          for (vector<int>::iterator ia = atomListLocal.begin();
432 <               ia != atomListLocal.end(); ++ia) {              
433 <            atom1 = (*ia);            
659 >
660 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
661              sdat = fDecomp_->fillSelfData(atom1);
662              interactionMan_->doPreForce(sdat);
663            }
664 +
665            fDecomp_->distributeIntermediateData();        
666          }
667        }
# Line 442 | Line 670 | namespace OpenMD {
670      
671      fDecomp_->collectData();
672      
673 <    if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) {
674 <      atomListLocal = fDecomp_->getAtomList();
675 <      for (vector<int>::iterator ia = atomListLocal.begin();
448 <           ia != atomListLocal.end(); ++ia) {              
449 <        atom1 = (*ia);    
673 >    if ( info_->requiresSkipCorrection() ) {
674 >      
675 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
676  
677 <        if (info_->requiresSkipCorrection()) {
678 <          vector<int> skipList = fDecomp_->getSkipsForAtom(atom1);
679 <          for (vector<int>::iterator jb = skipList.begin();
680 <               jb != skipList.end(); ++jb) {              
681 <            atom2 = (*jb);
682 <            idat = fDecomp_->fillSkipData(atom1, atom2);
683 <            interactionMan_->doSkipCorrection(idat);
684 <          }
677 >        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
678 >        
679 >        for (vector<int>::iterator jb = skipList.begin();
680 >             jb != skipList.end(); ++jb) {        
681 >    
682 >          atom2 = (*jb);
683 >          idat = fDecomp_->fillSkipData(atom1, atom2);
684 >          interactionMan_->doSkipCorrection(idat);
685 >
686          }
460 –          
461 –        if (info_->requiresSelfCorrection()) {
462 –          sdat = fDecomp_->fillSelfData(atom1);
463 –          interactionMan_->doSelfCorrection(sdat);
464 –        }
687        }
688      }
689 +    
690 +    if (info_->requiresSelfCorrection()) {
691  
692 <    // dangerous to iterate over enums, but we'll live on the edge:
693 <    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
694 <      lrPot += longRangePotential[i]; //Quick hack
692 >      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
693 >        sdat = fDecomp_->fillSelfData(atom1);
694 >        interactionMan_->doSelfCorrection(sdat);
695 >      }
696 >
697      }
698 <        
698 >
699 >    longRangePotential = fDecomp_->getLongRangePotential();
700 >    lrPot = longRangePotential.sum();
701 >
702      //store the tau and long range potential    
703      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
704      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];

Diff Legend

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