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 1570 by gezelter, Thu May 26 21:56:04 2011 UTC vs.
Revision 1576 by gezelter, Wed Jun 8 16:05:07 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    }
69 +
70 +  /**
71 +   * setupCutoffs
72 +   *
73 +   * Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
74 +   *
75 +   * cutoffRadius : realType
76 +   *  If the cutoffRadius was explicitly set, use that value.
77 +   *  If the cutoffRadius was not explicitly set:
78 +   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
79 +   *      No electrostatic atoms?  Poll the atom types present in the
80 +   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
81 +   *      Use the maximum suggested value that was found.
82 +   *
83 +   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
84 +   *      If cutoffMethod was explicitly set, use that choice.
85 +   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
86 +   *
87 +   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
88 +   *      If cutoffPolicy was explicitly set, use that choice.
89 +   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
90 +   */
91 +  void ForceManager::setupCutoffs() {
92 +    
93 +    Globals* simParams_ = info_->getSimParams();
94 +    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
95 +    
96 +    if (simParams_->haveCutoffRadius()) {
97 +      rCut_ = simParams_->getCutoffRadius();
98 +    } else {      
99 +      if (info_->usesElectrostaticAtoms()) {
100 +        sprintf(painCave.errMsg,
101 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
102 +                "\tOpenMD will use a default value of 12.0 angstroms"
103 +                "\tfor the cutoffRadius.\n");
104 +        painCave.isFatal = 0;
105 +        painCave.severity = OPENMD_INFO;
106 +        simError();
107 +        rCut_ = 12.0;
108 +      } else {
109 +        RealType thisCut;
110 +        set<AtomType*>::iterator i;
111 +        set<AtomType*> atomTypes;
112 +        atomTypes = info_->getSimulatedAtomTypes();        
113 +        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
114 +          thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
115 +          rCut_ = max(thisCut, rCut_);
116 +        }
117 +        sprintf(painCave.errMsg,
118 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
119 +                "\tOpenMD will use %lf angstroms.\n",
120 +                rCut_);
121 +        painCave.isFatal = 0;
122 +        painCave.severity = OPENMD_INFO;
123 +        simError();
124 +      }            
125 +    }
126 +
127 +    map<string, CutoffMethod> stringToCutoffMethod;
128 +    stringToCutoffMethod["HARD"] = HARD;
129 +    stringToCutoffMethod["SWITCHED"] = SWITCHED;
130 +    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
131 +    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
132    
133 <  void ForceManager::calcForces() {
133 >    if (simParams_->haveCutoffMethod()) {
134 >      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
135 >      map<string, CutoffMethod>::iterator i;
136 >      i = stringToCutoffMethod.find(cutMeth);
137 >      if (i == stringToCutoffMethod.end()) {
138 >        sprintf(painCave.errMsg,
139 >                "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
140 >                "\tShould be one of: "
141 >                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
142 >                cutMeth.c_str());
143 >        painCave.isFatal = 1;
144 >        painCave.severity = OPENMD_ERROR;
145 >        simError();
146 >      } else {
147 >        cutoffMethod_ = i->second;
148 >      }
149 >    } else {
150 >      sprintf(painCave.errMsg,
151 >              "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
152 >              "\tOpenMD will use SHIFTED_FORCE.\n");
153 >      painCave.isFatal = 0;
154 >      painCave.severity = OPENMD_INFO;
155 >      simError();
156 >      cutoffMethod_ = SHIFTED_FORCE;        
157 >    }
158 >
159 >    map<string, CutoffPolicy> stringToCutoffPolicy;
160 >    stringToCutoffPolicy["MIX"] = MIX;
161 >    stringToCutoffPolicy["MAX"] = MAX;
162 >    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
163 >
164 >    std::string cutPolicy;
165 >    if (forceFieldOptions_.haveCutoffPolicy()){
166 >      cutPolicy = forceFieldOptions_.getCutoffPolicy();
167 >    }else if (simParams_->haveCutoffPolicy()) {
168 >      cutPolicy = simParams_->getCutoffPolicy();
169 >    }
170 >
171 >    if (!cutPolicy.empty()){
172 >      toUpper(cutPolicy);
173 >      map<string, CutoffPolicy>::iterator i;
174 >      i = stringToCutoffPolicy.find(cutPolicy);
175 >
176 >      if (i == stringToCutoffPolicy.end()) {
177 >        sprintf(painCave.errMsg,
178 >                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
179 >                "\tShould be one of: "
180 >                "MIX, MAX, or TRADITIONAL\n",
181 >                cutPolicy.c_str());
182 >        painCave.isFatal = 1;
183 >        painCave.severity = OPENMD_ERROR;
184 >        simError();
185 >      } else {
186 >        cutoffPolicy_ = i->second;
187 >      }
188 >    } else {
189 >      sprintf(painCave.errMsg,
190 >              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
191 >              "\tOpenMD will use TRADITIONAL.\n");
192 >      painCave.isFatal = 0;
193 >      painCave.severity = OPENMD_INFO;
194 >      simError();
195 >      cutoffPolicy_ = TRADITIONAL;        
196 >    }
197 >  }
198 >
199 >  /**
200 >   * setupSwitching
201 >   *
202 >   * Sets the values of switchingRadius and
203 >   *  If the switchingRadius was explicitly set, use that value (but check it)
204 >   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
205 >   */
206 >  void ForceManager::setupSwitching() {
207 >    Globals* simParams_ = info_->getSimParams();
208      
209 +    if (simParams_->haveSwitchingRadius()) {
210 +      rSwitch_ = simParams_->getSwitchingRadius();
211 +      if (rSwitch_ > rCut_) {        
212 +        sprintf(painCave.errMsg,
213 +                "ForceManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
214 +                rSwitch_, rCut_);
215 +        painCave.isFatal = 1;
216 +        painCave.severity = OPENMD_ERROR;
217 +        simError();
218 +      }
219 +    } else {      
220 +      rSwitch_ = 0.85 * rCut_;
221 +      sprintf(painCave.errMsg,
222 +              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
223 +              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
224 +              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
225 +      painCave.isFatal = 0;
226 +      painCave.severity = OPENMD_WARNING;
227 +      simError();
228 +    }          
229 +    
230 +    if (simParams_->haveSwitchingFunctionType()) {
231 +      string funcType = simParams_->getSwitchingFunctionType();
232 +      toUpper(funcType);
233 +      if (funcType == "CUBIC") {
234 +        sft_ = cubic;
235 +      } else {
236 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
237 +          sft_ = fifth_order_poly;
238 +        } else {
239 +          // throw error        
240 +          sprintf( painCave.errMsg,
241 +                   "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
242 +                   "\tswitchingFunctionType must be one of: "
243 +                   "\"cubic\" or \"fifth_order_polynomial\".",
244 +                   funcType.c_str() );
245 +          painCave.isFatal = 1;
246 +          painCave.severity = OPENMD_ERROR;
247 +          simError();
248 +        }          
249 +      }
250 +    }
251 +    switcher_->setSwitchType(sft_);
252 +    switcher_->setSwitch(rSwitch_, rCut_);
253 +  }
254 +  
255 +  void ForceManager::initialize() {
256 +
257      if (!info_->isTopologyDone()) {
258        info_->update();
259        interactionMan_->setSimInfo(info_);
260        interactionMan_->initialize();
261 <      swfun_ = interactionMan_->getSwitchingFunction();
262 <      fDecomp_->distributeInitialData();
263 <      info_->prepareTopology();
261 >
262 >      // We want to delay the cutoffs until after the interaction
263 >      // manager has set up the atom-atom interactions so that we can
264 >      // query them for suggested cutoff values
265 >
266 >      setupCutoffs();
267 >      setupSwitching();
268 >
269 >      info_->prepareTopology();      
270      }
271 +
272 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
273 +    
274 +    // Force fields can set options on how to scale van der Waals and electrostatic
275 +    // interactions for atoms connected via bonds, bends and torsions
276 +    // in this case the topological distance between atoms is:
277 +    // 0 = topologically unconnected
278 +    // 1 = bonded together
279 +    // 2 = connected via a bend
280 +    // 3 = connected via a torsion
281 +    
282 +    vdwScale_.reserve(4);
283 +    fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
284 +
285 +    electrostaticScale_.reserve(4);
286 +    fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
287 +
288 +    vdwScale_[0] = 1.0;
289 +    vdwScale_[1] = fopts.getvdw12scale();
290 +    vdwScale_[2] = fopts.getvdw13scale();
291 +    vdwScale_[3] = fopts.getvdw14scale();
292      
293 +    electrostaticScale_[0] = 1.0;
294 +    electrostaticScale_[1] = fopts.getelectrostatic12scale();
295 +    electrostaticScale_[2] = fopts.getelectrostatic13scale();
296 +    electrostaticScale_[3] = fopts.getelectrostatic14scale();    
297 +    
298 +    fDecomp_->distributeInitialData();
299 +
300 +    initialized_ = true;
301 +
302 +  }
303 +
304 +  void ForceManager::calcForces() {
305 +    
306 +    if (!initialized_) initialize();
307 +
308      preCalculation();  
309      shortRangeInteractions();
310      longRangeInteractions();
311 <    postCalculation();
88 <    
311 >    postCalculation();    
312    }
313    
314    void ForceManager::preCalculation() {
# Line 265 | Line 488 | namespace OpenMD {
488        rc = pos;
489      }
490      
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 –
491      // new stuff starts here:
492 <
492 >    fDecomp_->zeroWorkArrays();
493      fDecomp_->distributeData();
494  
495      int cg1, cg2, atom1, atom2;
# Line 284 | Line 497 | namespace OpenMD {
497      RealType rgrpsq, rgrp;
498      RealType vij;
499      Vector3d fij, fg;
500 <    pair<int, int> gtypes;
500 >    tuple3<RealType, RealType, RealType> cuts;
501      RealType rCutSq;
502      bool in_switching_region;
503      RealType sw, dswdr, swderiv;
# Line 292 | Line 505 | namespace OpenMD {
505      InteractionData idat;
506      SelfData sdat;
507      RealType mf;
508 +    potVec pot(0.0);
509 +    potVec longRangePotential(0.0);
510 +    RealType lrPot;
511  
512      int loopStart, loopEnd;
513  
# Line 315 | Line 531 | namespace OpenMD {
531          
532          cg1 = (*it).first;
533          cg2 = (*it).second;
534 +        
535 +        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
536  
319 –        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
537          d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
538          curSnapshot->wrapVector(d_grp);        
539          rgrpsq = d_grp.lengthSquare();
323 –        rCutSq = groupCutoffMap[gtypes].first;
540  
541 +        rCutSq = cuts.second;
542 +
543          if (rgrpsq < rCutSq) {
544 <          *(idat.rcut) = groupCutoffMap[gtypes].second;
544 >          *(idat.rcut) = cuts.first;
545            if (iLoop == PAIR_LOOP) {
546              vij *= 0.0;
547              fij = V3Zero;
548            }
549            
550 <          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
551 <                                                  rgrp);              
550 >          in_switching_region = switcher_->getSwitch(rgrpsq, *(idat.sw), dswdr,
551 >                                                     rgrp);
552 >              
553            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
554            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
555  
# Line 344 | Line 563 | namespace OpenMD {
563                
564                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
565                  
566 +                pot *= 0.0;
567 +
568                  idat = fDecomp_->fillInteractionData(atom1, atom2);
569 +                *(idat.pot) = pot;
570  
571                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
572                    *(idat.d) = d_grp;
# Line 361 | Line 583 | namespace OpenMD {
583                    interactionMan_->doPrePair(idat);
584                  } else {
585                    interactionMan_->doPair(idat);
586 +                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
587                    vij += *(idat.vpair);
588                    fij += *(idat.f1);
589                    tau -= outProduct( *(idat.d), *(idat.f1));
# Line 467 | Line 690 | namespace OpenMD {
690  
691      }
692  
693 <    // dangerous to iterate over enums, but we'll live on the edge:
694 <    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
695 <      lrPot += longRangePotential[i]; //Quick hack
473 <    }
474 <        
693 >    longRangePotential = fDecomp_->getLongRangePotential();
694 >    lrPot = longRangePotential.sum();
695 >
696      //store the tau and long range potential    
697      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
698      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];

Diff Legend

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