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 1503 by gezelter, Sat Oct 2 19:54:41 2010 UTC vs.
Revision 1575 by gezelter, Fri Jun 3 21:39:49 2011 UTC

# Line 49 | Line 49
49  
50   #include "brains/ForceManager.hpp"
51   #include "primitives/Molecule.hpp"
52 – #include "UseTheForce/doForces_interface.h"
52   #define __OPENMD_C
54 – #include "UseTheForce/DarkSide/fInteractionMap.h"
53   #include "utils/simError.h"
54   #include "primitives/Bond.hpp"
55   #include "primitives/Bend.hpp"
56   #include "primitives/Torsion.hpp"
57   #include "primitives/Inversion.hpp"
58 + #include "nonbonded/NonBondedInteraction.hpp"
59 + #include "parallel/ForceMatrixDecomposition.hpp"
60  
61 + using namespace std;
62   namespace OpenMD {
63    
64 <  ForceManager::ForceManager(SimInfo * info) : info_(info),
65 <                                               NBforcesInitialized_(false) {
64 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
65 >
66 >    fDecomp_ = new ForceMatrixDecomposition(info_);
67    }
68 <
68 >  
69    void ForceManager::calcForces() {
70      
71 <    if (!info_->isFortranInitialized()) {
71 >    if (!info_->isTopologyDone()) {
72        info_->update();
73 +      interactionMan_->setSimInfo(info_);
74 +      interactionMan_->initialize();
75 +      swfun_ = interactionMan_->getSwitchingFunction();
76 +      info_->prepareTopology();
77 +      fDecomp_->distributeInitialData();
78      }
79      
80 <    preCalculation();
81 <    
82 <    calcShortRangeInteraction();
76 <
77 <    calcLongRangeInteraction();
78 <
80 >    preCalculation();  
81 >    shortRangeInteractions();
82 >    longRangeInteractions();
83      postCalculation();
84      
85    }
# Line 87 | Line 91 | namespace OpenMD {
91      Atom* atom;
92      Molecule::RigidBodyIterator rbIter;
93      RigidBody* rb;
94 +    Molecule::CutoffGroupIterator ci;
95 +    CutoffGroup* cg;
96      
97      // forces are zeroed here, before any are accumulated.
92 –    // NOTE: do not rezero the forces in Fortran.
98      
99      for (mol = info_->beginMolecule(mi); mol != NULL;
100           mol = info_->nextMolecule(mi)) {
# Line 102 | Line 107 | namespace OpenMD {
107             rb = mol->nextRigidBody(rbIter)) {
108          rb->zeroForcesAndTorques();
109        }        
110 <          
110 >
111 >      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
112 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
113 >            cg = mol->nextCutoffGroup(ci)) {
114 >          //calculate the center of mass of cutoff group
115 >          cg->updateCOM();
116 >        }
117 >      }      
118      }
119 <    
119 >  
120      // Zero out the stress tensor
121      tau *= 0.0;
122      
123    }
124    
125 <  void ForceManager::calcShortRangeInteraction() {
125 >  void ForceManager::shortRangeInteractions() {
126      Molecule* mol;
127      RigidBody* rb;
128      Bond* bond;
# Line 152 | Line 164 | namespace OpenMD {
164          RealType currBendPot = bend->getPotential();          
165          
166          bendPotential += bend->getPotential();
167 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
167 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
168          if (i == bendDataSets.end()) {
169            BendDataSet dataSet;
170            dataSet.prev.angle = dataSet.curr.angle = angle;
171            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
172            dataSet.deltaV = 0.0;
173 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
173 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
174          }else {
175            i->second.prev.angle = i->second.curr.angle;
176            i->second.prev.potential = i->second.curr.potential;
# Line 175 | Line 187 | namespace OpenMD {
187          torsion->calcForce(angle);
188          RealType currTorsionPot = torsion->getPotential();
189          torsionPotential += torsion->getPotential();
190 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
190 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
191          if (i == torsionDataSets.end()) {
192            TorsionDataSet dataSet;
193            dataSet.prev.angle = dataSet.curr.angle = angle;
194            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
195            dataSet.deltaV = 0.0;
196 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
196 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
197          }else {
198            i->second.prev.angle = i->second.curr.angle;
199            i->second.prev.potential = i->second.curr.potential;
# Line 191 | Line 203 | namespace OpenMD {
203                                     i->second.prev.potential);
204          }      
205        }      
206 <
206 >      
207        for (inversion = mol->beginInversion(inversionIter);
208             inversion != NULL;
209             inversion = mol->nextInversion(inversionIter)) {
# Line 199 | Line 211 | namespace OpenMD {
211          inversion->calcForce(angle);
212          RealType currInversionPot = inversion->getPotential();
213          inversionPotential += inversion->getPotential();
214 <        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
214 >        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
215          if (i == inversionDataSets.end()) {
216            InversionDataSet dataSet;
217            dataSet.prev.angle = dataSet.curr.angle = angle;
218            dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
219            dataSet.deltaV = 0.0;
220 <          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
220 >          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
221          }else {
222            i->second.prev.angle = i->second.curr.angle;
223            i->second.prev.potential = i->second.curr.potential;
# Line 224 | Line 236 | namespace OpenMD {
236      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
237      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
238      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
239 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
228 <    
239 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
240    }
241    
242 <  void ForceManager::calcLongRangeInteraction() {
232 <    Snapshot* curSnapshot;
233 <    DataStorage* config;
234 <    RealType* frc;
235 <    RealType* pos;
236 <    RealType* trq;
237 <    RealType* A;
238 <    RealType* electroFrame;
239 <    RealType* rc;
240 <    RealType* particlePot;
241 <    
242 <    //get current snapshot from SimInfo
243 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
244 <    
245 <    //get array pointers
246 <    config = &(curSnapshot->atomData);
247 <    frc = config->getArrayPointer(DataStorage::dslForce);
248 <    pos = config->getArrayPointer(DataStorage::dslPosition);
249 <    trq = config->getArrayPointer(DataStorage::dslTorque);
250 <    A   = config->getArrayPointer(DataStorage::dslAmat);
251 <    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
252 <    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
242 >  void ForceManager::longRangeInteractions() {
243  
244 <    //calculate the center of mass of cutoff group
245 <    SimInfo::MoleculeIterator mi;
246 <    Molecule* mol;
247 <    Molecule::CutoffGroupIterator ci;
248 <    CutoffGroup* cg;
249 <    Vector3d com;
250 <    std::vector<Vector3d> rcGroup;
251 <    
252 <    if(info_->getNCutoffGroups() > 0){
253 <      
254 <      for (mol = info_->beginMolecule(mi); mol != NULL;
255 <           mol = info_->nextMolecule(mi)) {
256 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
257 <            cg = mol->nextCutoffGroup(ci)) {
268 <          cg->getCOM(com);
269 <          rcGroup.push_back(com);
270 <        }
271 <      }// end for (mol)
272 <      
273 <      rc = rcGroup[0].getArrayPointer();
244 >    // some of this initial stuff will go away:
245 >    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
246 >    DataStorage* config = &(curSnapshot->atomData);
247 >    DataStorage* cgConfig = &(curSnapshot->cgData);
248 >    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
249 >    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
250 >    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
251 >    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
252 >    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
253 >    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
254 >    RealType* rc;    
255 >
256 >    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
257 >      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
258      } else {
259        // center of mass of the group is the same as position of the atom  
260        // if cutoff group does not exist
261        rc = pos;
262      }
263      
264 <    //initialize data before passing to fortran
265 <    RealType longRangePotential[LR_POT_TYPES];
266 <    RealType lrPot = 0.0;
267 <    int isError = 0;
264 >    // new stuff starts here:
265 >    fDecomp_->zeroWorkArrays();
266 >    fDecomp_->distributeData();
267 >
268 >    int cg1, cg2, atom1, atom2;
269 >    Vector3d d_grp, dag;
270 >    RealType rgrpsq, rgrp;
271 >    RealType vij;
272 >    Vector3d fij, fg;
273 >    pair<int, int> gtypes;
274 >    RealType rCutSq;
275 >    bool in_switching_region;
276 >    RealType sw, dswdr, swderiv;
277 >    vector<int> atomListColumn, atomListRow, atomListLocal;
278 >    InteractionData idat;
279 >    SelfData sdat;
280 >    RealType mf;
281 >    potVec pot(0.0);
282 >    potVec longRangePotential(0.0);
283 >    RealType lrPot;
284  
285 <    for (int i=0; i<LR_POT_TYPES;i++){
286 <      longRangePotential[i]=0.0; //Initialize array
285 >    int loopStart, loopEnd;
286 >
287 >    loopEnd = PAIR_LOOP;
288 >    if (info_->requiresPrepair() ) {
289 >      loopStart = PREPAIR_LOOP;
290 >    } else {
291 >      loopStart = PAIR_LOOP;
292      }
293 <    
294 <    doForceLoop(pos,
295 <                rc,
296 <                A,
297 <                electroFrame,
298 <                frc,
299 <                trq,
300 <                tau.getArrayPointer(),
301 <                longRangePotential,
302 <                particlePot,
303 <                &isError );
304 <    
305 <    if( isError ){
306 <      sprintf( painCave.errMsg,
307 <               "Error returned from the fortran force calculation.\n" );
308 <      painCave.isFatal = 1;
309 <      simError();
310 <    }
311 <    for (int i=0; i<LR_POT_TYPES;i++){
312 <      lrPot += longRangePotential[i]; //Quick hack
293 >
294 >    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
295 >      
296 >      if (iLoop == loopStart) {
297 >        bool update_nlist = fDecomp_->checkNeighborList();
298 >        if (update_nlist)
299 >          neighborList = fDecomp_->buildNeighborList();
300 >      }
301 >
302 >      for (vector<pair<int, int> >::iterator it = neighborList.begin();
303 >             it != neighborList.end(); ++it) {
304 >        
305 >        cg1 = (*it).first;
306 >        cg2 = (*it).second;
307 >
308 >        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
309 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
310 >        curSnapshot->wrapVector(d_grp);        
311 >        rgrpsq = d_grp.lengthSquare();
312 >        rCutSq = groupCutoffMap[gtypes].first;
313 >
314 >        if (rgrpsq < rCutSq) {
315 >          *(idat.rcut) = groupCutoffMap[gtypes].second;
316 >          if (iLoop == PAIR_LOOP) {
317 >            vij *= 0.0;
318 >            fij = V3Zero;
319 >          }
320 >          
321 >          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
322 >                                                  rgrp);              
323 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
324 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
325 >
326 >          for (vector<int>::iterator ia = atomListRow.begin();
327 >               ia != atomListRow.end(); ++ia) {            
328 >            atom1 = (*ia);
329 >            
330 >            for (vector<int>::iterator jb = atomListColumn.begin();
331 >                 jb != atomListColumn.end(); ++jb) {              
332 >              atom2 = (*jb);
333 >              
334 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
335 >                
336 >                pot *= 0.0;
337 >
338 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
339 >                *(idat.pot) = pot;
340 >
341 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
342 >                  *(idat.d) = d_grp;
343 >                  *(idat.r2) = rgrpsq;
344 >                } else {
345 >                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
346 >                  curSnapshot->wrapVector( *(idat.d) );
347 >                  *(idat.r2) = idat.d->lengthSquare();
348 >                }
349 >                
350 >                *(idat.rij) = sqrt( *(idat.r2) );
351 >              
352 >                if (iLoop == PREPAIR_LOOP) {
353 >                  interactionMan_->doPrePair(idat);
354 >                } else {
355 >                  interactionMan_->doPair(idat);
356 >                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
357 >                  vij += *(idat.vpair);
358 >                  fij += *(idat.f1);
359 >                  tau -= outProduct( *(idat.d), *(idat.f1));
360 >                }
361 >              }
362 >            }
363 >          }
364 >
365 >          if (iLoop == PAIR_LOOP) {
366 >            if (in_switching_region) {
367 >              swderiv = vij * dswdr / rgrp;
368 >              fg = swderiv * d_grp;
369 >
370 >              fij += fg;
371 >
372 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
373 >                tau -= outProduct( *(idat.d), fg);
374 >              }
375 >          
376 >              for (vector<int>::iterator ia = atomListRow.begin();
377 >                   ia != atomListRow.end(); ++ia) {            
378 >                atom1 = (*ia);                
379 >                mf = fDecomp_->getMassFactorRow(atom1);
380 >                // fg is the force on atom ia due to cutoff group's
381 >                // presence in switching region
382 >                fg = swderiv * d_grp * mf;
383 >                fDecomp_->addForceToAtomRow(atom1, fg);
384 >
385 >                if (atomListRow.size() > 1) {
386 >                  if (info_->usesAtomicVirial()) {
387 >                    // find the distance between the atom
388 >                    // and the center of the cutoff group:
389 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
390 >                    tau -= outProduct(dag, fg);
391 >                  }
392 >                }
393 >              }
394 >              for (vector<int>::iterator jb = atomListColumn.begin();
395 >                   jb != atomListColumn.end(); ++jb) {              
396 >                atom2 = (*jb);
397 >                mf = fDecomp_->getMassFactorColumn(atom2);
398 >                // fg is the force on atom jb due to cutoff group's
399 >                // presence in switching region
400 >                fg = -swderiv * d_grp * mf;
401 >                fDecomp_->addForceToAtomColumn(atom2, fg);
402 >
403 >                if (atomListColumn.size() > 1) {
404 >                  if (info_->usesAtomicVirial()) {
405 >                    // find the distance between the atom
406 >                    // and the center of the cutoff group:
407 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
408 >                    tau -= outProduct(dag, fg);
409 >                  }
410 >                }
411 >              }
412 >            }
413 >            //if (!SIM_uses_AtomicVirial) {
414 >            //  tau -= outProduct(d_grp, fij);
415 >            //}
416 >          }
417 >        }
418 >      }
419 >
420 >      if (iLoop == PREPAIR_LOOP) {
421 >        if (info_->requiresPrepair()) {            
422 >          fDecomp_->collectIntermediateData();
423 >
424 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
425 >            sdat = fDecomp_->fillSelfData(atom1);
426 >            interactionMan_->doPreForce(sdat);
427 >          }
428 >
429 >          fDecomp_->distributeIntermediateData();        
430 >        }
431 >      }
432 >
433      }
434 +    
435 +    fDecomp_->collectData();
436 +    
437 +    if ( info_->requiresSkipCorrection() ) {
438 +      
439 +      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
440 +
441 +        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
442          
443 +        for (vector<int>::iterator jb = skipList.begin();
444 +             jb != skipList.end(); ++jb) {        
445 +    
446 +          atom2 = (*jb);
447 +          idat = fDecomp_->fillSkipData(atom1, atom2);
448 +          interactionMan_->doSkipCorrection(idat);
449 +
450 +        }
451 +      }
452 +    }
453 +    
454 +    if (info_->requiresSelfCorrection()) {
455 +
456 +      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
457 +        sdat = fDecomp_->fillSelfData(atom1);
458 +        interactionMan_->doSelfCorrection(sdat);
459 +      }
460 +
461 +    }
462 +
463 +    longRangePotential = fDecomp_->getLongRangePotential();
464 +    lrPot = longRangePotential.sum();
465 +
466      //store the tau and long range potential    
467      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
468 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
469 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
468 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
469 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
470    }
471  
472    

Diff Legend

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