[ViewVC] Diff of: OpenMD/branches/devel_omp/src/brains/ForceManager.cpp

Comparing branches/devel_omp/src/brains/ForceManager.cpp (file contents):
Revision 1594 by chuckv, Tue Jul 19 16:45:30 2011 UTC vs.
Revision 1595 by chuckv, Tue Jul 19 18:50:04 2011 UTC

+    preCalculation();
+    shortRangeInteractions();
-<
+    longRangeInteractions();
->
+//    longRangeInteractions();
->
+    longRangeInteractionsRapaport();
+    postCalculation();
+    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
-+
+  void ForceManager::longRangeInteractionsRapaport() {
-+
-+
+    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-+
+    DataStorage* config = &(curSnapshot->atomData);
-+
+    DataStorage* cgConfig = &(curSnapshot->cgData);
-+
-+
+    //calculate the center of mass of cutoff group
-+
-+
+    SimInfo::MoleculeIterator mi;
-+
+    Molecule* mol;
-+
+    Molecule::CutoffGroupIterator ci;
-+
+    CutoffGroup* cg;
-+
-+
+    if(info_->getNCutoffGroups() > 0){
-+
+      for (mol = info_->beginMolecule(mi); mol != NULL;
-+
+           mol = info_->nextMolecule(mi)) {
-+
+        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
-+
+            cg = mol->nextCutoffGroup(ci)) {
-+
+          cerr << "branch1\n";
-+
+          cerr << "globind = " << cg->getGlobalIndex() << "\n";
-+
+          cg->updateCOM();
-+
+        }
-+
+      }
-+
+    } else {
-+
+      // center of mass of the group is the same as position of the atom
-+
+      // if cutoff group does not exist
-+
+      cerr << "branch2\n";
-+
+      cgConfig->position = config->position;
-+
+    }
-+
-+
+    fDecomp_->zeroWorkArrays();
-+
+    fDecomp_->distributeData();
-+
-+
+    int cg1, cg2, atom1, atom2, topoDist;
-+
+    Vector3d d_grp, dag, d;
-+
+    RealType rgrpsq, rgrp, r2, r;
-+
+    RealType electroMult, vdwMult;
-+
+    RealType vij;
-+
+    Vector3d fij, fg, f1;
-+
+    tuple3<RealType, RealType, RealType> cuts;
-+
+    RealType rCutSq;
-+
+    bool in_switching_region;
-+
+    RealType sw, dswdr, swderiv;
-+
+    vector<int> atomListColumn, atomListRow, atomListLocal;
-+
+    InteractionData idat;
-+
+    SelfData sdat;
-+
+    RealType mf;
-+
+    RealType lrPot;
-+
+    RealType vpair;
-+
+    potVec longRangePotential(0.0);
-+
+    potVec workPot(0.0);
-+
-+
+    int loopStart, loopEnd;
-+
-+
+    idat.vdwMult = &vdwMult;
-+
+    idat.electroMult = &electroMult;
-+
+    idat.pot = &workPot;
-+
+    sdat.pot = fDecomp_->getEmbeddingPotential();
-+
+    idat.vpair = &vpair;
-+
+    idat.f1 = &f1;
-+
+    idat.sw = &sw;
-+
+    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
-+
+    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
-+
-+
+    loopEnd = PAIR_LOOP;
-+
+    if (info_->requiresPrepair() ) {
-+
+      loopStart = PREPAIR_LOOP;
-+
+    } else {
-+
+      loopStart = PAIR_LOOP;
-+
+    }
-+
-+
+    for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
-+
-+
+      if (iLoop == loopStart) {
-+
+        bool update_nlist = fDecomp_->checkNeighborList();
-+
+        if (update_nlist)
-+
+                neighborMatW = fDecomp_->buildLayerBasedNeighborList();
-+
+      }
-+
-+
+      int i;
-+
+#pragma omp parallel for num_threads(2) private(i)
-+
+      for(i = 0; i < neighborMatW.size(); ++i)
-+
+          for(vector<int>::iterator j = neighborMatW[i].begin(); j != neighborMatW[i].end(); ++j)
-+
+                  {
-+
+                         cg1 = i;
-+
+                         cg2 = *j;
-+
-+
+                        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
-+
-+
+                        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
-+
+                        curSnapshot->wrapVector(d_grp);
-+
+                        rgrpsq = d_grp.lengthSquare();
-+
-+
+                        rCutSq = cuts.second;
-+
-+
+                        if (rgrpsq < rCutSq) {
-+
+                          idat.rcut = &cuts.first;
-+
+                          if (iLoop == PAIR_LOOP) {
-+
+                                vij = 0.0;
-+
+                                fij = V3Zero;
-+
+                          }
-+
-+
+                          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
-+
+                                                                                                                 rgrp);
-+
-+
+                          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
-+
+                          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
-+
-+
+                          for (vector<int>::iterator ia = atomListRow.begin();
-+
+                                   ia != atomListRow.end(); ++ia) {
-+
+                                atom1 = (*ia);
-+
-+
+                                for (vector<int>::iterator jb = atomListColumn.begin();
-+
+                                         jb != atomListColumn.end(); ++jb) {
-+
+                                  atom2 = (*jb);
-+
-+
+                                  if (!fDecomp_->skipAtomPair(atom1, atom2)) {
-+
+                                        vpair = 0.0;
-+
+                                        workPot = 0.0;
-+
+                                        f1 = V3Zero;
-+
-+
+                                        fDecomp_->fillInteractionData(idat, atom1, atom2);
-+
-+
+                                        topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
-+
+                                        vdwMult = vdwScale_[topoDist];
-+
+                                        electroMult = electrostaticScale_[topoDist];
-+
-+
+                                        if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
-+
+                                          idat.d = &d_grp;
-+
+                                          idat.r2 = &rgrpsq;
-+
+                                          cerr << "dgrp = " << d_grp << "\n";
-+
+                                        } else {
-+
+                                          d = fDecomp_->getInteratomicVector(atom1, atom2);
-+
+                                          curSnapshot->wrapVector( d );
-+
+                                          r2 = d.lengthSquare();
-+
+                                          cerr << "datm = " << d<< "\n";
-+
+                                          idat.d = &d;
-+
+                                          idat.r2 = &r2;
-+
+                                        }
-+
-+
+                                        cerr << "idat.d = " << *(idat.d) << "\n";
-+
+                                        r = sqrt( *(idat.r2) );
-+
+                                        idat.rij = &r;
-+
-+
+                                        if (iLoop == PREPAIR_LOOP) {
-+
+                                          interactionMan_->doPrePair(idat);
-+
+                                        } else {
-+
+                                          interactionMan_->doPair(idat);
-+
+                                          fDecomp_->unpackInteractionData(idat, atom1, atom2);
-+
-+
+                                          cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << "\n";
-+
+                                          vij += vpair;
-+
+                                          fij += f1;
-+
+                                          tau -= outProduct( *(idat.d), f1);
-+
+                                        }
-+
+                                  }
-+
+                                }
-+
+                          }
-+
-+
+                          if (iLoop == PAIR_LOOP) {
-+
+                                if (in_switching_region) {
-+
+                                  swderiv = vij * dswdr / rgrp;
-+
+                                  fg = swderiv * d_grp;
-+
+                                  fij += fg;
-+
-+
+                                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
-+
+                                        tau -= outProduct( *(idat.d), fg);
-+
+                                  }
-+
-+
+                                  for (vector<int>::iterator ia = atomListRow.begin();
-+
+                                           ia != atomListRow.end(); ++ia) {
-+
+                                        atom1 = (*ia);
-+
+                                        mf = fDecomp_->getMassFactorRow(atom1);
-+
+                                        // fg is the force on atom ia due to cutoff group's
-+
+                                        // presence in switching region
-+
+                                        fg = swderiv * d_grp * mf;
-+
+                                        fDecomp_->addForceToAtomRow(atom1, fg);
-+
-+
+                                        if (atomListRow.size() > 1) {
-+
+                                          if (info_->usesAtomicVirial()) {
-+
+                                                // find the distance between the atom
-+
+                                                // and the center of the cutoff group:
-+
+                                                dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
-+
+                                                tau -= outProduct(dag, fg);
-+
+                                          }
-+
+                                        }
-+
+                                  }
-+
+                                  for (vector<int>::iterator jb = atomListColumn.begin();
-+
+                                           jb != atomListColumn.end(); ++jb) {
-+
+                                        atom2 = (*jb);
-+
+                                        mf = fDecomp_->getMassFactorColumn(atom2);
-+
+                                        // fg is the force on atom jb due to cutoff group's
-+
+                                        // presence in switching region
-+
+                                        fg = -swderiv * d_grp * mf;
-+
+                                        fDecomp_->addForceToAtomColumn(atom2, fg);
-+
-+
+                                        if (atomListColumn.size() > 1) {
-+
+                                          if (info_->usesAtomicVirial()) {
-+
+                                                // find the distance between the atom
-+
+                                                // and the center of the cutoff group:
-+
+                                                dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
-+
+                                                tau -= outProduct(dag, fg);
-+
+                                          }
-+
+                                        }
-+
+                                  }
-+
+                                }
-+
+                                //if (!SIM_uses_AtomicVirial) {
-+
+                                //  tau -= outProduct(d_grp, fij);
-+
+                                //}
-+
+                          }
-+
+                        }
-+
+                  }
-+
-+
+      if (iLoop == PREPAIR_LOOP) {
-+
+        if (info_->requiresPrepair()) {
-+
-+
+          fDecomp_->collectIntermediateData();
-+
-+
+          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
-+
+            fDecomp_->fillSelfData(sdat, atom1);
-+
+            interactionMan_->doPreForce(sdat);
-+
+          }
-+
-+
+          fDecomp_->distributeIntermediateData();
-+
-+
+        }
-+
+      }
-+
-+
+    }
-+
-+
+    fDecomp_->collectData();
-+
-+
+    if (info_->requiresSelfCorrection()) {
-+
-+
+      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
-+
+        fDecomp_->fillSelfData(sdat, atom1);
-+
+        interactionMan_->doSelfCorrection(sdat);
-+
+      }
-+
-+
+    }
-+
-+
+    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
-+
+      *(fDecomp_->getPairwisePotential());
-+
-+
+    lrPot = longRangePotential.sum();
-+
-+
+    //store the tau and long range potential
-+
+    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
-+
+    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
-+
+    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
-+
+  }
-+
+  void ForceManager::longRangeInteractions() {
+    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
+      if (iLoop == loopStart) {
+        bool update_nlist = fDecomp_->checkNeighborList();
-<
+        if (update_nlist)
->
+        if (update_nlist)
+          neighborList = fDecomp_->buildNeighborList();
-+
-<
+      for (vector<pair<int, int> >::iterator it = neighborList.begin();
-<
+             it != neighborList.end(); ++it) {
-<
->
+      for (vector<pair<int, int> >::iterator it = neighborList.begin();
->
+             it != neighborList.end(); ++it)
->
+      {
+        cg1 = (*it).first;
+        cg2 = (*it).second;

Diff Legend

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

Comparing branches/devel_omp/src/brains/ForceManager.cpp (file contents): Revision 1594 by chuckv, Tue Jul 19 16:45:30 2011 UTC vs. Revision 1595 by chuckv, Tue Jul 19 18:50:04 2011 UTC

Diff Legend

Comparing branches/devel_omp/src/brains/ForceManager.cpp (file contents):
Revision 1594 by chuckv, Tue Jul 19 16:45:30 2011 UTC vs.
Revision 1595 by chuckv, Tue Jul 19 18:50:04 2011 UTC