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 1584 by gezelter, Fri Jun 17 20:16:35 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 58 | 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    
69    ForceManager::ForceManager(SimInfo * info) : info_(info) {
70 <
71 < #ifdef IS_MPI
72 <    fDecomp_ = new ForceMatrixDecomposition(info_);
68 < #else
69 <    // fDecomp_ = new ForceSerialDecomposition(info);
70 < #endif
70 >    forceField_ = info_->getForceField();
71 >    interactionMan_ = new InteractionManager();
72 >    fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
73    }
74 +
75 +  /**
76 +   * setupCutoffs
77 +   *
78 +   * Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
79 +   *
80 +   * cutoffRadius : realType
81 +   *  If the cutoffRadius was explicitly set, use that value.
82 +   *  If the cutoffRadius was not explicitly set:
83 +   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
84 +   *      No electrostatic atoms?  Poll the atom types present in the
85 +   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
86 +   *      Use the maximum suggested value that was found.
87 +   *
88 +   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
89 +   *      If cutoffMethod was explicitly set, use that choice.
90 +   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
91 +   *
92 +   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
93 +   *      If cutoffPolicy was explicitly set, use that choice.
94 +   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
95 +   */
96 +  void ForceManager::setupCutoffs() {
97 +    
98 +    Globals* simParams_ = info_->getSimParams();
99 +    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
100 +    
101 +    if (simParams_->haveCutoffRadius()) {
102 +      rCut_ = simParams_->getCutoffRadius();
103 +    } else {      
104 +      if (info_->usesElectrostaticAtoms()) {
105 +        sprintf(painCave.errMsg,
106 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
107 +                "\tOpenMD will use a default value of 12.0 angstroms"
108 +                "\tfor the cutoffRadius.\n");
109 +        painCave.isFatal = 0;
110 +        painCave.severity = OPENMD_INFO;
111 +        simError();
112 +        rCut_ = 12.0;
113 +      } else {
114 +        RealType thisCut;
115 +        set<AtomType*>::iterator i;
116 +        set<AtomType*> atomTypes;
117 +        atomTypes = info_->getSimulatedAtomTypes();        
118 +        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
119 +          thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
120 +          rCut_ = max(thisCut, rCut_);
121 +        }
122 +        sprintf(painCave.errMsg,
123 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
124 +                "\tOpenMD will use %lf angstroms.\n",
125 +                rCut_);
126 +        painCave.isFatal = 0;
127 +        painCave.severity = OPENMD_INFO;
128 +        simError();
129 +      }
130 +    }
131 +
132 +    fDecomp_->setUserCutoff(rCut_);
133 +    interactionMan_->setCutoffRadius(rCut_);
134 +
135 +    map<string, CutoffMethod> stringToCutoffMethod;
136 +    stringToCutoffMethod["HARD"] = HARD;
137 +    stringToCutoffMethod["SWITCHED"] = SWITCHED;
138 +    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
139 +    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
140    
141 <  void ForceManager::calcForces() {
141 >    if (simParams_->haveCutoffMethod()) {
142 >      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
143 >      map<string, CutoffMethod>::iterator i;
144 >      i = stringToCutoffMethod.find(cutMeth);
145 >      if (i == stringToCutoffMethod.end()) {
146 >        sprintf(painCave.errMsg,
147 >                "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
148 >                "\tShould be one of: "
149 >                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
150 >                cutMeth.c_str());
151 >        painCave.isFatal = 1;
152 >        painCave.severity = OPENMD_ERROR;
153 >        simError();
154 >      } else {
155 >        cutoffMethod_ = i->second;
156 >      }
157 >    } else {
158 >      sprintf(painCave.errMsg,
159 >              "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
160 >              "\tOpenMD will use SHIFTED_FORCE.\n");
161 >      painCave.isFatal = 0;
162 >      painCave.severity = OPENMD_INFO;
163 >      simError();
164 >      cutoffMethod_ = SHIFTED_FORCE;        
165 >    }
166 >
167 >    map<string, CutoffPolicy> stringToCutoffPolicy;
168 >    stringToCutoffPolicy["MIX"] = MIX;
169 >    stringToCutoffPolicy["MAX"] = MAX;
170 >    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
171 >
172 >    std::string cutPolicy;
173 >    if (forceFieldOptions_.haveCutoffPolicy()){
174 >      cutPolicy = forceFieldOptions_.getCutoffPolicy();
175 >    }else if (simParams_->haveCutoffPolicy()) {
176 >      cutPolicy = simParams_->getCutoffPolicy();
177 >    }
178 >
179 >    if (!cutPolicy.empty()){
180 >      toUpper(cutPolicy);
181 >      map<string, CutoffPolicy>::iterator i;
182 >      i = stringToCutoffPolicy.find(cutPolicy);
183 >
184 >      if (i == stringToCutoffPolicy.end()) {
185 >        sprintf(painCave.errMsg,
186 >                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
187 >                "\tShould be one of: "
188 >                "MIX, MAX, or TRADITIONAL\n",
189 >                cutPolicy.c_str());
190 >        painCave.isFatal = 1;
191 >        painCave.severity = OPENMD_ERROR;
192 >        simError();
193 >      } else {
194 >        cutoffPolicy_ = i->second;
195 >      }
196 >    } else {
197 >      sprintf(painCave.errMsg,
198 >              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
199 >              "\tOpenMD will use TRADITIONAL.\n");
200 >      painCave.isFatal = 0;
201 >      painCave.severity = OPENMD_INFO;
202 >      simError();
203 >      cutoffPolicy_ = TRADITIONAL;        
204 >    }
205 >    fDecomp_->setCutoffPolicy(cutoffPolicy_);
206 >  }
207 >
208 >  /**
209 >   * setupSwitching
210 >   *
211 >   * Sets the values of switchingRadius and
212 >   *  If the switchingRadius was explicitly set, use that value (but check it)
213 >   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
214 >   */
215 >  void ForceManager::setupSwitching() {
216 >    Globals* simParams_ = info_->getSimParams();
217 >
218 >    // create the switching function object:
219 >    switcher_ = new SwitchingFunction();
220      
221 +    if (simParams_->haveSwitchingRadius()) {
222 +      rSwitch_ = simParams_->getSwitchingRadius();
223 +      if (rSwitch_ > rCut_) {        
224 +        sprintf(painCave.errMsg,
225 +                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
226 +                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
227 +        painCave.isFatal = 1;
228 +        painCave.severity = OPENMD_ERROR;
229 +        simError();
230 +      }
231 +    } else {      
232 +      rSwitch_ = 0.85 * rCut_;
233 +      sprintf(painCave.errMsg,
234 +              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
235 +              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
236 +              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
237 +      painCave.isFatal = 0;
238 +      painCave.severity = OPENMD_WARNING;
239 +      simError();
240 +    }          
241 +    
242 +    // Default to cubic switching function.
243 +    sft_ = cubic;
244 +    if (simParams_->haveSwitchingFunctionType()) {
245 +      string funcType = simParams_->getSwitchingFunctionType();
246 +      toUpper(funcType);
247 +      if (funcType == "CUBIC") {
248 +        sft_ = cubic;
249 +      } else {
250 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
251 +          sft_ = fifth_order_poly;
252 +        } else {
253 +          // throw error        
254 +          sprintf( painCave.errMsg,
255 +                   "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
256 +                   "\tswitchingFunctionType must be one of: "
257 +                   "\"cubic\" or \"fifth_order_polynomial\".",
258 +                   funcType.c_str() );
259 +          painCave.isFatal = 1;
260 +          painCave.severity = OPENMD_ERROR;
261 +          simError();
262 +        }          
263 +      }
264 +    }
265 +    switcher_->setSwitchType(sft_);
266 +    switcher_->setSwitch(rSwitch_, rCut_);
267 +    interactionMan_->setSwitchingRadius(rSwitch_);
268 +  }
269 +  
270 +  void ForceManager::initialize() {
271 +
272      if (!info_->isTopologyDone()) {
273        info_->update();
274        interactionMan_->setSimInfo(info_);
275        interactionMan_->initialize();
276 <      swfun_ = interactionMan_->getSwitchingFunction();
277 <      fDecomp_->distributeInitialData();
278 <      info_->prepareTopology();
276 >
277 >      // We want to delay the cutoffs until after the interaction
278 >      // manager has set up the atom-atom interactions so that we can
279 >      // query them for suggested cutoff values
280 >
281 >      setupCutoffs();
282 >      setupSwitching();
283 >
284 >      info_->prepareTopology();      
285      }
286 +
287 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
288      
289 +    // Force fields can set options on how to scale van der Waals and electrostatic
290 +    // interactions for atoms connected via bonds, bends and torsions
291 +    // in this case the topological distance between atoms is:
292 +    // 0 = topologically unconnected
293 +    // 1 = bonded together
294 +    // 2 = connected via a bend
295 +    // 3 = connected via a torsion
296 +    
297 +    vdwScale_.reserve(4);
298 +    fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
299 +
300 +    electrostaticScale_.reserve(4);
301 +    fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
302 +
303 +    vdwScale_[0] = 1.0;
304 +    vdwScale_[1] = fopts.getvdw12scale();
305 +    vdwScale_[2] = fopts.getvdw13scale();
306 +    vdwScale_[3] = fopts.getvdw14scale();
307 +    
308 +    electrostaticScale_[0] = 1.0;
309 +    electrostaticScale_[1] = fopts.getelectrostatic12scale();
310 +    electrostaticScale_[2] = fopts.getelectrostatic13scale();
311 +    electrostaticScale_[3] = fopts.getelectrostatic14scale();    
312 +    
313 +    fDecomp_->distributeInitialData();
314 +
315 +    initialized_ = true;
316 +
317 +  }
318 +
319 +  void ForceManager::calcForces() {
320 +    
321 +    if (!initialized_) initialize();
322 +
323      preCalculation();  
324      shortRangeInteractions();
325      longRangeInteractions();
326 <    postCalculation();
88 <    
326 >    postCalculation();    
327    }
328    
329    void ForceManager::preCalculation() {
# Line 245 | Line 483 | namespace OpenMD {
483    
484    void ForceManager::longRangeInteractions() {
485  
248 –    // some of this initial stuff will go away:
486      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
487      DataStorage* config = &(curSnapshot->atomData);
488      DataStorage* cgConfig = &(curSnapshot->cgData);
252 –    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
253 –    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
254 –    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
255 –    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
256 –    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
257 –    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
258 –    RealType* rc;    
489  
490 <    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
491 <      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
490 >    //calculate the center of mass of cutoff group
491 >
492 >    SimInfo::MoleculeIterator mi;
493 >    Molecule* mol;
494 >    Molecule::CutoffGroupIterator ci;
495 >    CutoffGroup* cg;
496 >
497 >    if(info_->getNCutoffGroups() > 0){      
498 >      for (mol = info_->beginMolecule(mi); mol != NULL;
499 >           mol = info_->nextMolecule(mi)) {
500 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
501 >            cg = mol->nextCutoffGroup(ci)) {
502 >          cg->updateCOM();
503 >        }
504 >      }      
505      } else {
506        // center of mass of the group is the same as position of the atom  
507        // if cutoff group does not exist
508 <      rc = pos;
508 >      cgConfig->position = config->position;
509      }
267 –    
268 –    //initialize data before passing to fortran
269 –    RealType longRangePotential[N_INTERACTION_FAMILIES];
270 –    RealType lrPot = 0.0;
271 –    int isError = 0;
510  
511 <    // 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 <
278 <    // new stuff starts here:
279 <
511 >    fDecomp_->zeroWorkArrays();
512      fDecomp_->distributeData();
513 <
514 <    int cg1, cg2, atom1, atom2;
515 <    Vector3d d_grp, dag;
516 <    RealType rgrpsq, rgrp;
513 >    
514 >    int cg1, cg2, atom1, atom2, topoDist;
515 >    Vector3d d_grp, dag, d;
516 >    RealType rgrpsq, rgrp, r2, r;
517 >    RealType electroMult, vdwMult;
518      RealType vij;
519 <    Vector3d fij, fg;
520 <    pair<int, int> gtypes;
519 >    Vector3d fij, fg, f1;
520 >    tuple3<RealType, RealType, RealType> cuts;
521      RealType rCutSq;
522      bool in_switching_region;
523      RealType sw, dswdr, swderiv;
# Line 292 | Line 525 | namespace OpenMD {
525      InteractionData idat;
526      SelfData sdat;
527      RealType mf;
528 +    RealType lrPot;
529 +    RealType vpair;
530 +    potVec longRangePotential(0.0);
531 +    potVec workPot(0.0);
532  
533      int loopStart, loopEnd;
534  
535 +    idat.vdwMult = &vdwMult;
536 +    idat.electroMult = &electroMult;
537 +    idat.pot = &workPot;
538 +    sdat.pot = fDecomp_->getEmbeddingPotential();
539 +    idat.vpair = &vpair;
540 +    idat.f1 = &f1;
541 +    idat.sw = &sw;
542 +    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
543 +    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
544 +    
545      loopEnd = PAIR_LOOP;
546      if (info_->requiresPrepair() ) {
547        loopStart = PREPAIR_LOOP;
548      } else {
549        loopStart = PAIR_LOOP;
550      }
551 <
552 <    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
553 <      
551 >  
552 >    for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
553 >    
554        if (iLoop == loopStart) {
555          bool update_nlist = fDecomp_->checkNeighborList();
556          if (update_nlist)
557            neighborList = fDecomp_->buildNeighborList();
558 <      }
559 <
558 >      }      
559 >        
560        for (vector<pair<int, int> >::iterator it = neighborList.begin();
561               it != neighborList.end(); ++it) {
562 <        
562 >                
563          cg1 = (*it).first;
564          cg2 = (*it).second;
565 +        
566 +        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
567  
319 –        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
568          d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
569          curSnapshot->wrapVector(d_grp);        
570          rgrpsq = d_grp.lengthSquare();
323 –        rCutSq = groupCutoffMap[gtypes].first;
571  
572 +        rCutSq = cuts.second;
573 +
574          if (rgrpsq < rCutSq) {
575 <          *(idat.rcut) = groupCutoffMap[gtypes].second;
575 >          idat.rcut = &cuts.first;
576            if (iLoop == PAIR_LOOP) {
577              vij *= 0.0;
578              fij = V3Zero;
579            }
580            
581 <          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
582 <                                                  rgrp);              
581 >          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
582 >                                                     rgrp);
583 >              
584            atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
585            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
586  
# Line 341 | Line 591 | namespace OpenMD {
591              for (vector<int>::iterator jb = atomListColumn.begin();
592                   jb != atomListColumn.end(); ++jb) {              
593                atom2 = (*jb);
594 <              
594 >            
595                if (!fDecomp_->skipAtomPair(atom1, atom2)) {
596 +                vpair = 0.0;
597 +                workPot = 0.0;
598 +                f1 = V3Zero;
599 +
600 +                fDecomp_->fillInteractionData(idat, atom1, atom2);
601                  
602 <                idat = fDecomp_->fillInteractionData(atom1, atom2);
602 >                topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
603 >                vdwMult = vdwScale_[topoDist];
604 >                electroMult = electrostaticScale_[topoDist];
605  
606                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
607 <                  *(idat.d) = d_grp;
608 <                  *(idat.r2) = rgrpsq;
607 >                  idat.d = &d_grp;
608 >                  idat.r2 = &rgrpsq;
609                  } else {
610 <                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
611 <                  curSnapshot->wrapVector( *(idat.d) );
612 <                  *(idat.r2) = idat.d->lengthSquare();
610 >                  d = fDecomp_->getInteratomicVector(atom1, atom2);
611 >                  curSnapshot->wrapVector( d );
612 >                  r2 = d.lengthSquare();
613 >                  idat.d = &d;
614 >                  idat.r2 = &r2;
615                  }
616                  
617 <                *(idat.rij) = sqrt( *(idat.r2) );
617 >                r = sqrt( *(idat.r2) );
618 >                idat.rij = &r;
619                
620                  if (iLoop == PREPAIR_LOOP) {
621                    interactionMan_->doPrePair(idat);
622                  } else {
623                    interactionMan_->doPair(idat);
624 <                  vij += *(idat.vpair);
625 <                  fij += *(idat.f1);
626 <                  tau -= outProduct( *(idat.d), *(idat.f1));
624 >                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
625 >                  vij += vpair;
626 >                  fij += f1;
627 >                  tau -= outProduct( *(idat.d), f1);
628                  }
629                }
630              }
# Line 429 | Line 690 | namespace OpenMD {
690            fDecomp_->collectIntermediateData();
691  
692            for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
693 <            sdat = fDecomp_->fillSelfData(atom1);
693 >            fDecomp_->fillSelfData(sdat, atom1);
694              interactionMan_->doPreForce(sdat);
695            }
696 <
696 >          
697 >          
698            fDecomp_->distributeIntermediateData();        
699          }
700        }
# Line 445 | Line 707 | namespace OpenMD {
707        
708        for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
709  
710 <        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
710 >        vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 );
711          
712          for (vector<int>::iterator jb = skipList.begin();
713               jb != skipList.end(); ++jb) {        
714      
715            atom2 = (*jb);
716 <          idat = fDecomp_->fillSkipData(atom1, atom2);
716 >          fDecomp_->fillSkipData(idat, atom1, atom2);
717            interactionMan_->doSkipCorrection(idat);
718 +          fDecomp_->unpackSkipData(idat, atom1, atom2);
719  
720          }
721        }
# Line 461 | Line 724 | namespace OpenMD {
724      if (info_->requiresSelfCorrection()) {
725  
726        for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
727 <        sdat = fDecomp_->fillSelfData(atom1);
727 >        fDecomp_->fillSelfData(sdat, atom1);
728          interactionMan_->doSelfCorrection(sdat);
729        }
730  
731      }
732  
733 <    // dangerous to iterate over enums, but we'll live on the edge:
734 <    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
735 <      lrPot += longRangePotential[i]; //Quick hack
736 <    }
737 <        
733 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
734 >      *(fDecomp_->getPairwisePotential());
735 >
736 >    lrPot = longRangePotential.sum();
737 >
738      //store the tau and long range potential    
739      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
740      curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];

Diff Legend

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