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 1544 by gezelter, Fri Mar 18 19:31:52 2011 UTC vs.
Revision 1545 by gezelter, Fri Apr 8 21:25:19 2011 UTC

# Line 60 | Line 60
60   #include "parallel/ForceDecomposition.hpp"
61   //#include "parallel/SerialDecomposition.hpp"
62  
63 + using namespace std;
64   namespace OpenMD {
65    
66 <  ForceManager::ForceManager(SimInfo * info) : info_(info),
67 <                                               NBforcesInitialized_(false) {
66 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
67 >
68   #ifdef IS_MPI
69      decomp_ = new ForceDecomposition(info_);
70   #else
71 <    //  decomp_ = new SerialDecomposition(info);
71 >    // decomp_ = new SerialDecomposition(info);
72   #endif
73    }
74 <
74 >  
75    void ForceManager::calcForces() {
76      
76 –
77      if (!info_->isFortranInitialized()) {
78        info_->update();
79        nbiMan_->setSimInfo(info_);
80        nbiMan_->initialize();
81 +      swfun_ = nbiMan_->getSwitchingFunction();
82        decomp_->distributeInitialData();
83        info_->setupFortran();
84      }
# Line 100 | Line 101 | namespace OpenMD {
101      CutoffGroup* cg;
102      
103      // forces are zeroed here, before any are accumulated.
103 –    // NOTE: do not rezero the forces in Fortran.
104      
105      for (mol = info_->beginMolecule(mi); mol != NULL;
106           mol = info_->nextMolecule(mi)) {
# Line 115 | Line 115 | namespace OpenMD {
115        }        
116  
117        if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
118 –        std::cerr << "should not see me \n";
118          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
119              cg = mol->nextCutoffGroup(ci)) {
120            //calculate the center of mass of cutoff group
# Line 171 | Line 170 | namespace OpenMD {
170          RealType currBendPot = bend->getPotential();          
171          
172          bendPotential += bend->getPotential();
173 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
173 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
174          if (i == bendDataSets.end()) {
175            BendDataSet dataSet;
176            dataSet.prev.angle = dataSet.curr.angle = angle;
177            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
178            dataSet.deltaV = 0.0;
179 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
179 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
180          }else {
181            i->second.prev.angle = i->second.curr.angle;
182            i->second.prev.potential = i->second.curr.potential;
# Line 194 | Line 193 | namespace OpenMD {
193          torsion->calcForce(angle);
194          RealType currTorsionPot = torsion->getPotential();
195          torsionPotential += torsion->getPotential();
196 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
196 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
197          if (i == torsionDataSets.end()) {
198            TorsionDataSet dataSet;
199            dataSet.prev.angle = dataSet.curr.angle = angle;
200            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
201            dataSet.deltaV = 0.0;
202 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
202 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
203          }else {
204            i->second.prev.angle = i->second.curr.angle;
205            i->second.prev.potential = i->second.curr.potential;
# Line 210 | Line 209 | namespace OpenMD {
209                                     i->second.prev.potential);
210          }      
211        }      
212 <
212 >      
213        for (inversion = mol->beginInversion(inversionIter);
214             inversion != NULL;
215             inversion = mol->nextInversion(inversionIter)) {
# Line 218 | Line 217 | namespace OpenMD {
217          inversion->calcForce(angle);
218          RealType currInversionPot = inversion->getPotential();
219          inversionPotential += inversion->getPotential();
220 <        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
220 >        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
221          if (i == inversionDataSets.end()) {
222            InversionDataSet dataSet;
223            dataSet.prev.angle = dataSet.curr.angle = angle;
224            dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
225            dataSet.deltaV = 0.0;
226 <          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
226 >          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
227          }else {
228            i->second.prev.angle = i->second.curr.angle;
229            i->second.prev.potential = i->second.curr.potential;
# Line 243 | Line 242 | namespace OpenMD {
242      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
243      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
244      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
245 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
247 <    
245 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
246    }
247    
248    void ForceManager::calcLongRangeInteraction() {
251 –    Snapshot* curSnapshot;
252 –    DataStorage* config;
253 –    DataStorage* cgConfig;
254 –    RealType* frc;
255 –    RealType* pos;
256 –    RealType* trq;
257 –    RealType* A;
258 –    RealType* electroFrame;
259 –    RealType* rc;
260 –    RealType* particlePot;
261 –    
262 –    //get current snapshot from SimInfo
263 –    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
264 –    
265 –    //get array pointers
266 –    config = &(curSnapshot->atomData);
267 –    cgConfig = &(curSnapshot->cgData);
268 –    frc = config->getArrayPointer(DataStorage::dslForce);
269 –    pos = config->getArrayPointer(DataStorage::dslPosition);
270 –    trq = config->getArrayPointer(DataStorage::dslTorque);
271 –    A   = config->getArrayPointer(DataStorage::dslAmat);
272 –    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
273 –    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
249  
250 +    // some of this initial stuff will go away:
251 +    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
252 +    DataStorage* config = &(curSnapshot->atomData);
253 +    DataStorage* cgConfig = &(curSnapshot->cgData);
254 +    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
255 +    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
256 +    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
257 +    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
258 +    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
259 +    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
260 +    RealType* rc;    
261 +
262      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
276 –      std::cerr << "should not see me \n";
263        rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
264      } else {
265        // center of mass of the group is the same as position of the atom  
# Line 289 | Line 275 | namespace OpenMD {
275      for (int i=0; i<LR_POT_TYPES;i++){
276        longRangePotential[i]=0.0; //Initialize array
277      }
278 <    
278 >
279 >    // new stuff starts here:
280 >
281      decomp_->distributeData();
282 +
283 +    int cg1, cg2;
284 +    Vector3d d_grp;
285 +    RealType rgrpsq, rgrp;
286 +    RealType vij;
287 +    Vector3d fij, fg;
288 +    pair<int, int> gtypes;
289 +    RealType rCutSq;
290 +    bool in_switching_region;
291 +    RealType sw, dswdr, swderiv;
292 +    vector<int> atomListI;
293 +    vector<int> atomListJ;
294 +    InteractionData idat;
295  
296 <    int nLoops = 1;
296 <    for (int iLoop = 0; iLoop < nLoops; iLoop++) {
297 <      doForceLoop(pos,
298 <                  rc,
299 <                  A,
300 <                  electroFrame,
301 <                  frc,
302 <                  trq,
303 <                  tau.getArrayPointer(),
304 <                  longRangePotential,
305 <                  particlePot,
306 <                  &isError );  
296 >    int loopStart, loopEnd;
297  
298 <      if (nLoops > 1) {
299 <        decomp_->collectIntermediateData();
300 <        decomp_->distributeIntermediateData();
301 <      }
298 >    loopEnd = PAIR_LOOP;
299 >    if (info_->requiresPrepair_) {
300 >      loopStart = PREPAIR_LOOP;
301 >    } else {
302 >      loopStart = PAIR_LOOP;
303      }
313 –  
314 –    decomp_->collectData();
304  
305 <    if( isError ){
306 <      sprintf( painCave.errMsg,
307 <               "Error returned from the fortran force calculation.\n" );
308 <      painCave.isFatal = 1;
309 <      simError();
305 >    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
306 >      
307 >      if (iLoop == loopStart) {
308 >        bool update_nlist = decomp_->checkNeighborList();
309 >        if (update_nlist)
310 >          neighborList = decomp_->buildNeighborList();
311 >      }
312 >
313 >      for (vector<pair<int, int> >::iterator it = neighborList.begin();
314 >             it != neighborList.end(); ++it) {
315 >        
316 >        cg1 = (*it).first;
317 >        cg2 = (*it).second;
318 >
319 >        gtypes = decomp_->getGroupTypes(cg1, cg2);
320 >        d_grp  = decomp_->getIntergroupVector(cg1, cg2);
321 >        curSnapshot->wrapVector(d_grp);        
322 >        rgrpsq = d_grp.lengthSquare();
323 >        rCutSq = groupCutoffMap(gtypes).first;
324 >
325 >        if (rgrpsq < rCutSq) {
326 >          idat.rcut = groupCutoffMap(gtypes).second;
327 >          if (iLoop == PAIR_LOOP) {
328 >            vij = 0.0;
329 >            fij = V3Zero;
330 >          }
331 >          
332 >          in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, idat.dswdr, rgrp);    
333 >          
334 >          atomListI = decomp_->getAtomsInGroupI(cg1);
335 >          atomListJ = decomp_->getAtomsInGroupJ(cg2);
336 >
337 >          for (vector<int>::iterator ia = atomListI.begin();
338 >               ia != atomListI.end(); ++ia) {            
339 >            atom1 = (*ia);
340 >            
341 >            for (vector<int>::iterator jb = atomListJ.begin();
342 >                 jb != atomListJ.end(); ++jb) {              
343 >              atom2 = (*jb);
344 >              
345 >              if (!decomp_->skipAtomPair(atom1, atom2)) {
346 >                
347 >                if (atomListI.size() == 1 && atomListJ.size() == 1) {
348 >                  idat.d = d_grp;
349 >                  idat.r2 = rgrpsq;
350 >                } else {
351 >                  idat.d = decomp_->getInteratomicVector(atom1, atom2);
352 >                  curSnapshot->wrapVector(idat.d);
353 >                  idat.r2 = idat.d.lengthSquare();
354 >                }
355 >                
356 >                idat.r = sqrt(idat.r2);
357 >                decomp_->fillInteractionData(atom1, atom2, idat);
358 >                
359 >                if (iLoop == PREPAIR_LOOP) {
360 >                  interactionMan_->doPrePair(idat);
361 >                } else {
362 >                  interactionMan_->doPair(idat);
363 >                  vij += idat.vpair;
364 >                  fij += idat.f1;
365 >                  tau -= outProduct(idat.d, idat.f);
366 >                }
367 >              }
368 >            }
369 >          }
370 >
371 >          if (iLoop == PAIR_LOOP) {
372 >            if (in_switching_region) {
373 >              swderiv = vij * dswdr / rgrp;
374 >              fg = swderiv * d_grp;
375 >
376 >              fij += fg;
377 >
378 >              if (atomListI.size() == 1 && atomListJ.size() == 1) {
379 >                tau -= outProduct(idat.d, fg);
380 >              }
381 >          
382 >              for (vector<int>::iterator ia = atomListI.begin();
383 >                   ia != atomListI.end(); ++ia) {            
384 >                atom1 = (*ia);                
385 >                mf = decomp_->getMfactI(atom1);
386 >                // fg is the force on atom ia due to cutoff group's
387 >                // presence in switching region
388 >                fg = swderiv * d_grp * mf;
389 >                decomp_->addForceToAtomI(atom1, fg);
390 >
391 >                if (atomListI.size() > 1) {
392 >                  if (info_->usesAtomicVirial_) {
393 >                    // find the distance between the atom
394 >                    // and the center of the cutoff group:
395 >                    dag = decomp_->getAtomToGroupVectorI(atom1, cg1);
396 >                    tau -= outProduct(dag, fg);
397 >                  }
398 >                }
399 >              }
400 >              for (vector<int>::iterator jb = atomListJ.begin();
401 >                   jb != atomListJ.end(); ++jb) {              
402 >                atom2 = (*jb);
403 >                mf = decomp_->getMfactJ(atom2);
404 >                // fg is the force on atom jb due to cutoff group's
405 >                // presence in switching region
406 >                fg = -swderiv * d_grp * mf;
407 >                decomp_->addForceToAtomJ(atom2, fg);
408 >
409 >                if (atomListJ.size() > 1) {
410 >                  if (info_->usesAtomicVirial_) {
411 >                    // find the distance between the atom
412 >                    // and the center of the cutoff group:
413 >                    dag = decomp_->getAtomToGroupVectorJ(atom2, cg2);
414 >                    tau -= outProduct(dag, fg);
415 >                  }
416 >                }
417 >              }
418 >            }
419 >            //if (!SIM_uses_AtomicVirial) {
420 >            //  tau -= outProduct(d_grp, fij);
421 >            //}
422 >          }
423 >        }
424 >      }
425 >
426 >      if (iLoop == PREPAIR_LOOP) {
427 >        if (info_->requiresPrepair_) {            
428 >          decomp_->collectIntermediateData();
429 >          atomList = decomp_->getAtomList();
430 >          for (vector<int>::iterator ia = atomList.begin();
431 >               ia != atomList.end(); ++ia) {              
432 >            atom1 = (*ia);            
433 >            decomp_->populateSelfData(atom1, SelfData sdat);
434 >            interactionMan_->doPreForce(sdat);
435 >          }
436 >          decomp_->distributeIntermediateData();        
437 >        }
438 >      }
439 >
440      }
441 +    
442 +    decomp_->collectData();
443 +    
444 +    if (info_->requiresSkipCorrection_ || info_->requiresSelfCorrection_) {
445 +      atomList = decomp_->getAtomList();
446 +      for (vector<int>::iterator ia = atomList.begin();
447 +           ia != atomList.end(); ++ia) {              
448 +        atom1 = (*ia);    
449 +
450 +        if (info_->requiresSkipCorrection_) {
451 +          vector<int> skipList = decomp_->getSkipsForAtom(atom1);
452 +          for (vector<int>::iterator jb = skipList.begin();
453 +               jb != skipList.end(); ++jb) {              
454 +            atom2 = (*jb);
455 +            decomp_->populateSkipData(atom1, atom2, InteractionData idat);
456 +            interactionMan_->doSkipCorrection(idat);
457 +          }
458 +        }
459 +          
460 +        if (info_->requiresSelfCorrection_) {
461 +          decomp_->populateSelfData(atom1, SelfData sdat);
462 +          interactionMan_->doSelfCorrection(sdat);
463 +      }
464 +      
465 +      
466 +    }
467 +
468      for (int i=0; i<LR_POT_TYPES;i++){
469        lrPot += longRangePotential[i]; //Quick hack
470      }

Diff Legend

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