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 1483 by gezelter, Tue Jul 27 21:17:31 2010 UTC vs.
Revision 1551 by gezelter, Thu Apr 28 18:38:21 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),
64 <                                               NBforcesInitialized_(false) {
65 <    lj_ = LJ::Instance();
66 <    lj_->setForceField(info_->getForceField());
64 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
65  
66 <    gb_ = GB::Instance();
67 <    gb_->setForceField(info_->getForceField());
68 <
69 <    eam_ = EAM::Instance();
70 <    eam_->setForceField(info_->getForceField());
66 > #ifdef IS_MPI
67 >    fDecomp_ = new ForceMatrixDecomposition(info_);
68 > #else
69 >    // fDecomp_ = new ForceSerialDecomposition(info);
70 > #endif
71    }
72 <
72 >  
73    void ForceManager::calcForces() {
74      
75      if (!info_->isFortranInitialized()) {
76        info_->update();
77 +      interactionMan_->setSimInfo(info_);
78 +      interactionMan_->initialize();
79 +      swfun_ = interactionMan_->getSwitchingFunction();
80 +      fDecomp_->distributeInitialData();
81 +      info_->setupFortran();
82      }
83      
84 <    preCalculation();
85 <    
86 <    calcShortRangeInteraction();
84 <
85 <    calcLongRangeInteraction();
86 <
84 >    preCalculation();  
85 >    shortRangeInteractions();
86 >    longRangeInteractions();
87      postCalculation();
88      
89    }
# Line 95 | Line 95 | namespace OpenMD {
95      Atom* atom;
96      Molecule::RigidBodyIterator rbIter;
97      RigidBody* rb;
98 +    Molecule::CutoffGroupIterator ci;
99 +    CutoffGroup* cg;
100      
101      // forces are zeroed here, before any are accumulated.
100 –    // NOTE: do not rezero the forces in Fortran.
102      
103      for (mol = info_->beginMolecule(mi); mol != NULL;
104           mol = info_->nextMolecule(mi)) {
# Line 110 | Line 111 | namespace OpenMD {
111             rb = mol->nextRigidBody(rbIter)) {
112          rb->zeroForcesAndTorques();
113        }        
114 <          
114 >
115 >      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
116 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
117 >            cg = mol->nextCutoffGroup(ci)) {
118 >          //calculate the center of mass of cutoff group
119 >          cg->updateCOM();
120 >        }
121 >      }      
122      }
123 <    
123 >  
124      // Zero out the stress tensor
125      tau *= 0.0;
126      
127    }
128    
129 <  void ForceManager::calcShortRangeInteraction() {
129 >  void ForceManager::shortRangeInteractions() {
130      Molecule* mol;
131      RigidBody* rb;
132      Bond* bond;
# Line 160 | Line 168 | namespace OpenMD {
168          RealType currBendPot = bend->getPotential();          
169          
170          bendPotential += bend->getPotential();
171 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
171 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
172          if (i == bendDataSets.end()) {
173            BendDataSet dataSet;
174            dataSet.prev.angle = dataSet.curr.angle = angle;
175            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
176            dataSet.deltaV = 0.0;
177 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
177 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
178          }else {
179            i->second.prev.angle = i->second.curr.angle;
180            i->second.prev.potential = i->second.curr.potential;
# Line 183 | Line 191 | namespace OpenMD {
191          torsion->calcForce(angle);
192          RealType currTorsionPot = torsion->getPotential();
193          torsionPotential += torsion->getPotential();
194 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
194 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
195          if (i == torsionDataSets.end()) {
196            TorsionDataSet dataSet;
197            dataSet.prev.angle = dataSet.curr.angle = angle;
198            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
199            dataSet.deltaV = 0.0;
200 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
200 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
201          }else {
202            i->second.prev.angle = i->second.curr.angle;
203            i->second.prev.potential = i->second.curr.potential;
# Line 199 | Line 207 | namespace OpenMD {
207                                     i->second.prev.potential);
208          }      
209        }      
210 <
210 >      
211        for (inversion = mol->beginInversion(inversionIter);
212             inversion != NULL;
213             inversion = mol->nextInversion(inversionIter)) {
# Line 207 | Line 215 | namespace OpenMD {
215          inversion->calcForce(angle);
216          RealType currInversionPot = inversion->getPotential();
217          inversionPotential += inversion->getPotential();
218 <        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
218 >        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
219          if (i == inversionDataSets.end()) {
220            InversionDataSet dataSet;
221            dataSet.prev.angle = dataSet.curr.angle = angle;
222            dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
223            dataSet.deltaV = 0.0;
224 <          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
224 >          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
225          }else {
226            i->second.prev.angle = i->second.curr.angle;
227            i->second.prev.potential = i->second.curr.potential;
# Line 232 | Line 240 | namespace OpenMD {
240      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
241      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
242      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
243 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
236 <    
243 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
244    }
245    
246 <  void ForceManager::calcLongRangeInteraction() {
240 <    Snapshot* curSnapshot;
241 <    DataStorage* config;
242 <    RealType* frc;
243 <    RealType* pos;
244 <    RealType* trq;
245 <    RealType* A;
246 <    RealType* electroFrame;
247 <    RealType* rc;
248 <    RealType* particlePot;
249 <    
250 <    //get current snapshot from SimInfo
251 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
252 <    
253 <    //get array pointers
254 <    config = &(curSnapshot->atomData);
255 <    frc = config->getArrayPointer(DataStorage::dslForce);
256 <    pos = config->getArrayPointer(DataStorage::dslPosition);
257 <    trq = config->getArrayPointer(DataStorage::dslTorque);
258 <    A   = config->getArrayPointer(DataStorage::dslAmat);
259 <    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
260 <    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
246 >  void ForceManager::longRangeInteractions() {
247  
248 <    //calculate the center of mass of cutoff group
249 <    SimInfo::MoleculeIterator mi;
250 <    Molecule* mol;
251 <    Molecule::CutoffGroupIterator ci;
252 <    CutoffGroup* cg;
253 <    Vector3d com;
254 <    std::vector<Vector3d> rcGroup;
255 <    
256 <    if(info_->getNCutoffGroups() > 0){
257 <      
258 <      for (mol = info_->beginMolecule(mi); mol != NULL;
259 <           mol = info_->nextMolecule(mi)) {
260 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
261 <            cg = mol->nextCutoffGroup(ci)) {
276 <          cg->getCOM(com);
277 <          rcGroup.push_back(com);
278 <        }
279 <      }// end for (mol)
280 <      
281 <      rc = rcGroup[0].getArrayPointer();
248 >    // some of this initial stuff will go away:
249 >    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
250 >    DataStorage* config = &(curSnapshot->atomData);
251 >    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;    
259 >
260 >    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
261 >      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
262      } else {
263        // center of mass of the group is the same as position of the atom  
264        // if cutoff group does not exist
# Line 286 | Line 266 | namespace OpenMD {
266      }
267      
268      //initialize data before passing to fortran
269 <    RealType longRangePotential[LR_POT_TYPES];
269 >    RealType longRangePotential[N_INTERACTION_FAMILIES];
270      RealType lrPot = 0.0;
291 –    Vector3d totalDipole;
271      int isError = 0;
272  
273 <    for (int i=0; i<LR_POT_TYPES;i++){
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 <    
278 <    doForceLoop(pos,
279 <                rc,
280 <                A,
281 <                electroFrame,
282 <                frc,
283 <                trq,
284 <                tau.getArrayPointer(),
285 <                longRangePotential,
286 <                particlePot,
287 <                &isError );
288 <    
289 <    if( isError ){
290 <      sprintf( painCave.errMsg,
291 <               "Error returned from the fortran force calculation.\n" );
292 <      painCave.isFatal = 1;
293 <      simError();
277 >
278 >    // new stuff starts here:
279 >
280 >    fDecomp_->distributeData();
281 >
282 >    int cg1, cg2, atom1, atom2;
283 >    Vector3d d_grp, dag;
284 >    RealType rgrpsq, rgrp;
285 >    RealType vij;
286 >    Vector3d fij, fg;
287 >    pair<int, int> gtypes;
288 >    RealType rCutSq;
289 >    bool in_switching_region;
290 >    RealType sw, dswdr, swderiv;
291 >    vector<int> atomListColumn, atomListRow, atomListLocal;
292 >    InteractionData idat;
293 >    SelfData sdat;
294 >    RealType mf;
295 >
296 >    int loopStart, loopEnd;
297 >
298 >    loopEnd = PAIR_LOOP;
299 >    if (info_->requiresPrepair() ) {
300 >      loopStart = PREPAIR_LOOP;
301 >    } else {
302 >      loopStart = PAIR_LOOP;
303      }
304 <    for (int i=0; i<LR_POT_TYPES;i++){
305 <      lrPot += longRangePotential[i]; //Quick hack
317 <    }
318 <    
319 <    // grab the simulation box dipole moment if specified
320 <    if (info_->getCalcBoxDipole()){
321 <      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
304 >
305 >    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
306        
307 <      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
308 <      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
309 <      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
307 >      if (iLoop == loopStart) {
308 >        bool update_nlist = fDecomp_->checkNeighborList();
309 >        if (update_nlist)
310 >          neighborList = fDecomp_->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 = fDecomp_->getGroupTypes(cg1, cg2);
320 >        d_grp  = fDecomp_->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, dswdr, rgrp);              
333 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
334 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
335 >
336 >          for (vector<int>::iterator ia = atomListRow.begin();
337 >               ia != atomListRow.end(); ++ia) {            
338 >            atom1 = (*ia);
339 >            
340 >            for (vector<int>::iterator jb = atomListColumn.begin();
341 >                 jb != atomListColumn.end(); ++jb) {              
342 >              atom2 = (*jb);
343 >              
344 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
345 >                
346 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
347 >
348 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
349 >                  idat.d = d_grp;
350 >                  idat.r2 = rgrpsq;
351 >                } else {
352 >                  idat.d = fDecomp_->getInteratomicVector(atom1, atom2);
353 >                  curSnapshot->wrapVector(idat.d);
354 >                  idat.r2 = idat.d.lengthSquare();
355 >                }
356 >                
357 >                idat.rij = sqrt(idat.r2);
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.f1);
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 (atomListRow.size() == 1 && atomListColumn.size() == 1) {
379 >                tau -= outProduct(idat.d, fg);
380 >              }
381 >          
382 >              for (vector<int>::iterator ia = atomListRow.begin();
383 >                   ia != atomListRow.end(); ++ia) {            
384 >                atom1 = (*ia);                
385 >                mf = fDecomp_->getMfactRow(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 >                fDecomp_->addForceToAtomRow(atom1, fg);
390 >
391 >                if (atomListRow.size() > 1) {
392 >                  if (info_->usesAtomicVirial()) {
393 >                    // find the distance between the atom
394 >                    // and the center of the cutoff group:
395 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
396 >                    tau -= outProduct(dag, fg);
397 >                  }
398 >                }
399 >              }
400 >              for (vector<int>::iterator jb = atomListColumn.begin();
401 >                   jb != atomListColumn.end(); ++jb) {              
402 >                atom2 = (*jb);
403 >                mf = fDecomp_->getMfactColumn(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 >                fDecomp_->addForceToAtomColumn(atom2, fg);
408 >
409 >                if (atomListColumn.size() > 1) {
410 >                  if (info_->usesAtomicVirial()) {
411 >                    // find the distance between the atom
412 >                    // and the center of the cutoff group:
413 >                    dag = fDecomp_->getAtomToGroupVectorColumn(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 >          fDecomp_->collectIntermediateData();
429 >          atomListLocal = fDecomp_->getAtomList();
430 >          for (vector<int>::iterator ia = atomListLocal.begin();
431 >               ia != atomListLocal.end(); ++ia) {              
432 >            atom1 = (*ia);            
433 >            sdat = fDecomp_->fillSelfData(atom1);
434 >            interactionMan_->doPreForce(sdat);
435 >          }
436 >          fDecomp_->distributeIntermediateData();        
437 >        }
438 >      }
439 >
440      }
441      
442 +    fDecomp_->collectData();
443 +    
444 +    if (info_->requiresSkipCorrection() || info_->requiresSelfCorrection()) {
445 +      atomListLocal = fDecomp_->getAtomList();
446 +      for (vector<int>::iterator ia = atomListLocal.begin();
447 +           ia != atomListLocal.end(); ++ia) {              
448 +        atom1 = (*ia);    
449 +
450 +        if (info_->requiresSkipCorrection()) {
451 +          vector<int> skipList = fDecomp_->getSkipsForAtom(atom1);
452 +          for (vector<int>::iterator jb = skipList.begin();
453 +               jb != skipList.end(); ++jb) {              
454 +            atom2 = (*jb);
455 +            idat = fDecomp_->fillSkipData(atom1, atom2);
456 +            interactionMan_->doSkipCorrection(idat);
457 +          }
458 +        }
459 +          
460 +        if (info_->requiresSelfCorrection()) {
461 +          sdat = fDecomp_->fillSelfData(atom1);
462 +          interactionMan_->doSelfCorrection(sdat);
463 +        }
464 +      }
465 +    }
466 +
467 +    // dangerous to iterate over enums, but we'll live on the edge:
468 +    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
469 +      lrPot += longRangePotential[i]; //Quick hack
470 +    }
471 +        
472      //store the tau and long range potential    
473      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
474 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
475 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
474 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
475 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
476    }
477  
478    

Diff Legend

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