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

Comparing branches/devel_omp/src/brains/ForceManager.cpp (file contents):
Revision 1607 by mciznick, Wed Jul 27 14:26:53 2011 UTC vs.
Revision 1608 by mciznick, Tue Aug 9 01:58:56 2011 UTC

+        preCalculation();
+        shortRangeInteractions();
+        //    longRangeInteractions();
-<
+        longRangeInteractionsRapaport();
->
+        //      longRangeInteractionsRapaport();
->
+        longRangeInteractionsParallel();
+        postCalculation();
+        curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
+        curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
+        curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
-+
+}
-+
-+
+void ForceManager::longRangeInteractionsParallel() {
-+
-+
+        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() << ":" << __LINE__ << "\n";
-+
+                                cg->updateCOM();
-+
-+
+                                //                              cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
-+
+                                //                                              << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
-+
+                                //                                              << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
-+
+                        }
-+
+                }
-+
+        } else
-+
+        {
-+
+                // center of mass of the group is the same as position of the atom
-+
+                // if cutoff group does not exist
-+
+                //              cerr << ":" << __LINE__ << "branch2\n";
-+
+                cgConfig->position = config->position;
-+
+        }
-+
-+
+        fDecomp_->zeroWorkArrays();
-+
+        fDecomp_->distributeData();
-+
-+
+        int 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;
-+
-+
+        InteractionDataPrv idatPrv;
-+
-+
+        SelfData sdat;
-+
+        RealType mf;
-+
+        RealType lrPot;
-+
+        RealType vpair;
-+
+        potVec longRangePotential(0.0);
-+
+        potVec workPot(0.0);
-+
-+
+        int loopStart, loopEnd;
-+
+        sdat.pot = fDecomp_->getEmbeddingPotential();
-+
-+
+        vector<CutoffGroup *> cgs;
-+
-+
+        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
-+
+        {
-+
+                for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
-+
+                {
-+
+                        cgs.push_back(cg);
-+
+                }
-+
+        }
-+
-+
+        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();
-+
+                }
-+
-+
+                vector<CutoffGroup *>::iterator cg1;
-+
+                vector<CutoffGroup *>::iterator cg2;
-+
-+
+//              int nThreads = 2;
-+
+                int chunkSize = cgs.size() / (omp_get_num_threads() * 20);
-+
-+
+                //              printf("before omp loop\n");
-+
+#pragma omp parallel /*num_threads(nThreads)*/ default(none) shared(curSnapshot, iLoop, cgs, chunkSize) \
-+
+        private(cg1, cg2, cuts, d_grp, rgrpsq, rCutSq, idatPrv, vij, fij, in_switching_region, dswdr, rgrp, \
-+
+                        atomListRow, atomListColumn, atom1, atom2, topoDist, d, r2, swderiv, fg, mf, dag)
-+
+                {
-+
+                        idatPrv.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
-+
+                        idatPrv.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
-+
-+
+//                      printf("Thread %d\n", omp_get_thread_num());
-+
+#pragma omp for schedule(dynamic, chunkSize)
-+
+                        for (cg1 = cgs.begin(); cg1 < cgs.end(); ++cg1)
-+
+                        {
-+
+                                for (cg2 = neighborMatW[(*cg1)->getGlobalIndex()].begin(); cg2 < neighborMatW[(*cg1)->getGlobalIndex()].end(); ++cg2)
-+
+                                {
-+
-+
+                                        cuts = fDecomp_->getGroupCutoffs((*cg1)->getGlobalIndex(), (*cg2)->getGlobalIndex());
-+
-+
+                                        d_grp = fDecomp_->getIntergroupVector((*cg1), (*cg2));
-+
+                                        curSnapshot->wrapVector(d_grp);
-+
+                                        rgrpsq = d_grp.lengthSquare();
-+
-+
+                                        rCutSq = cuts.second;
-+
-+
+//                                      printf("Thread %d\tcg1:%d\tcg2:%d d_grp\tx:%f\ty:%f\tz:%f\trgrpsq:%f\n", omp_get_thread_num(), (*cg1)->getGlobalIndex(), (*cg2)->getGlobalIndex(), d_grp.x(), d_grp.y(), d_grp.z(), rgrpsq);
-+
-+
+                                        if (rgrpsq < rCutSq)
-+
+                                        {
-+
+                                                idatPrv.rcut = cuts.first;
-+
+                                                if (iLoop == PAIR_LOOP)
-+
+                                                {
-+
+                                                        vij = 0.0;
-+
+                                                        fij = V3Zero;
-+
+                                                }
-+
-+
+                                                in_switching_region = switcher_->getSwitch(rgrpsq, /*sw*/idatPrv.sw, dswdr, rgrp);
-+
-+
+//                                              printf("in_switching_region:%d\trgrpsq:%f\t*idatPrv.sw:%f\tdswdr:%f\trgrp:%f\n", (in_switching_region == false ? 0 : 1), rgrpsq, idatPrv.sw, dswdr, rgrp);
-+
-+
+                                                atomListRow = fDecomp_->getAtomsInGroupRow((*cg1)->getGlobalIndex());
-+
+                                                atomListColumn = fDecomp_->getAtomsInGroupColumn((*cg2)->getGlobalIndex());
-+
-+
+                                                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);
-+
-+
+//                                                              printf("atom1:%d atom2:%d\n", atom1, atom2);
-+
+                                                                if (!fDecomp_->skipAtomPair(atom1, atom2))
-+
+                                                                {
-+
+                                                                        idatPrv.vpair = 0.0;
-+
+                                                                        idatPrv.pot = 0.0;
-+
+                                                                        idatPrv.f1 = V3Zero;
-+
-+
+                                                                        fDecomp_->fillInteractionDataOMP(idatPrv, atom1, atom2);
-+
-+
+                                                                        topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
-+
+                                                                        idatPrv.vdwMult = vdwScale_[topoDist];
-+
+                                                                        idatPrv.electroMult = electrostaticScale_[topoDist];
-+
-+
+//                                                                      printf("topoDist:%d\tidatPrv.vdwMult:%f\tidatPrv.electroMult:%f\n", topoDist, idatPrv.vdwMult, idatPrv.electroMult);
-+
-+
+                                                                        if (atomListRow.size() == 1 && atomListColumn.size() == 1)
-+
+                                                                        {
-+
+                                                                                idatPrv.d = d_grp;
-+
+                                                                                idatPrv.r2 = rgrpsq;
-+
+                                                                                //                                                                              cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
-+
+                                                                        } else
-+
+                                                                        {
-+
+                                                                                d = fDecomp_->getInteratomicVector(atom1, atom2);
-+
+                                                                                curSnapshot->wrapVector(d);
-+
+                                                                                r2 = d.lengthSquare();
-+
+                                                                                //                                                                              cerr << "datm = " << d << ":" << __LINE__ << "\n";
-+
+                                                                                idatPrv.d = d;
-+
+                                                                                idatPrv.r2 = r2;
-+
+                                                                        }
-+
-+
+//                                                                      printf("idatPrv.d x:%f\ty:%f\tz:%f\tidatPrv.r2:%f\n", (idatPrv.d).x(), (idatPrv.d).y(), (idatPrv.d).z(), idatPrv.r2);
-+
-+
+                                                                        //                                                                      cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
-+
+                                                                        idatPrv.rij = sqrt((idatPrv.r2));
-+
+                                                                        //                                                                      cerr << "idat.rij = " << *(idat.rij) << "\n";
-+
-+
+                                                                        #pragma omp critical
-+
+                                                                        {
-+
+                                                                                interactionMan_->initializeOMP();
-+
+                                                                        }
-+
-+
+                                                                        if (iLoop == PREPAIR_LOOP)
-+
+                                                                        {
-+
+                                                                                interactionMan_->doPrePairOMP(idatPrv);
-+
+                                                                        } else
-+
+                                                                        {
-+
+                                                                                interactionMan_->doPairOMP(idatPrv);
-+
+                                                                                fDecomp_->unpackInteractionDataOMP(idatPrv, atom1, atom2);
-+
-+
+                                                                                //                                                                              cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
-+
+//                                                                              printf("d x:%f y:%f z:%f vpair:%f f1 x:%f y:%f z:%f\n", idatPrv.d.x(), idatPrv.d.y(), idatPrv.d.z(), idatPrv.vpair, idatPrv.f1.x(), idatPrv.f1.y(), idatPrv.f1.z());
-+
+                                                                                #pragma omp critical
-+
+                                                                                {
-+
+                                                                                        vij += idatPrv.vpair;
-+
+                                                                                        fij += idatPrv.f1;
-+
+                                                                                        tau -= outProduct(idatPrv.d, idatPrv.f1);
-+
-+
+//                                                                                      printf("vij:%f fij x:%f y:%f z:%f\n", vij, fij.x(), fij.y(), fij.z());
-+
+                                                                                }
-+
+                                                                        }
-+
+                                                                }
-+
+                                                        }
-+
+                                                }
-+
-+
+                                                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)
-+
+                                                                {
-+
+                                                                        #pragma omp critical
-+
+                                                                        {
-+
+                                                                                tau -= outProduct(idatPrv.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_->addForceToAtomRowOMP(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)->getGlobalIndex());
-+
+                                                                                        #pragma omp critical
-+
+                                                                                        {
-+
+                                                                                                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)->getGlobalIndex());
-+
+                                                                                        #pragma omp critical
-+
+                                                                                        {
-+
+                                                                                                tau -= outProduct(dag, fg);
-+
+                                                                                        }
-+
+                                                                                }
-+
+                                                                        }
-+
+                                                                }
-+
+                                                        }
-+
+                                                        //if (!SIM_uses_AtomicVirial) {
-+
+                                                        //  tau -= outProduct(d_grp, fij);
-+
+                                                        //}
-+
+                                                }
-+
+                                        }
-+
+                                }
-+
+                        }// END: omp for loop
-+
+                        //              printf("after omp loop\n");
-+
+                }
-+
-+
+                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::longRangeInteractionsRapaport() {
+                        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci))
-<
+//                              cerr << "branch1\n";
-<
+//                              cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
->
+                                //                              cerr << "branch1\n";
->
+                                //                              cerr << "globind = " << cg->getGlobalIndex() << ":" << __LINE__ << "\n";
+                                cg->updateCOM();
-<
+//                              cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
-<
+//                                              << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
-<
+//                                              << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
->
+                                //                              cerr << "gbI: " << cg->getGlobalIndex() << " locI: " << cg->getLocalIndex() << " x: "
->
+                                //                                              << cgConfig->position[cg->getLocalIndex()].x() << " y: " << cgConfig->position[cg->getLocalIndex()].y()
->
+                                //                                              << " z: " << cgConfig->position[cg->getLocalIndex()].z() << "\n";
+        } else
+                // center of mass of the group is the same as position of the atom
+                // if cutoff group does not exist
-<
+//              cerr << ":" << __LINE__ << "branch2\n";
->
+                //              cerr << ":" << __LINE__ << "branch2\n";
+                cgConfig->position = config->position;
+                                neighborMatW = fDecomp_->buildLayerBasedNeighborList();
-–
+                //              printf("before omp loop\n");
-–
+                //#pragma omp parallel for num_threads(3) default(none) shared(curSnapshot, idat, iLoop, sw, cerr) \
-–
+        private(i, j, cg1, cg2, cuts, d_grp, rgrpsq, rCutSq, vij, fij, in_switching_region, rgrp, dswdr, atomListRow, atomListColumn, atom1, atom2, vpair, workPot, f1, topoDist, vdwMult, electroMult, d, r2, r, swderiv, fg, mf, dag)
+                for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi))
+                        for (cg1 = mol->beginCutoffGroup(ci); cg1 != NULL; cg1 = mol->nextCutoffGroup(ci))
+                                //                      printf("Thread %d executes loop iteration %d\n", omp_get_thread_num(), i);
-<
+                                for (vector<CutoffGroup *>::iterator cg2 = neighborMatW[cg1->getGlobalIndex()].begin(); cg2 != neighborMatW[cg1->getGlobalIndex()].end(); ++cg2)
->
+                                for (vector<CutoffGroup *>::iterator cg2 = neighborMatW[cg1->getGlobalIndex()].begin(); cg2
->
+                                                != neighborMatW[cg1->getGlobalIndex()].end(); ++cg2)
+                                        cuts = fDecomp_->getGroupCutoffs(cg1->getGlobalIndex(), (*cg2)->getGlobalIndex());
+                                                                                idat.d = &d_grp;
+                                                                                idat.r2 = &rgrpsq;
-<
+//                                                                              cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
->
+                                                                                //                                                                              cerr << "dgrp = " << d_grp << ":" << __LINE__ << "\n";
+                                                                        } else
+                                                                                d = fDecomp_->getInteratomicVector(atom1, atom2);
+                                                                                curSnapshot->wrapVector(d);
+                                                                                r2 = d.lengthSquare();
-<
+//                                                                              cerr << "datm = " << d << ":" << __LINE__ << "\n";
->
+                                                                                //                                                                              cerr << "datm = " << d << ":" << __LINE__ << "\n";
+                                                                                idat.d = &d;
+                                                                                idat.r2 = &r2;
-<
+//                                                                      cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
->
+                                                                        //                                                                      cerr << "idat.d = " << *(idat.d) << ":" << __LINE__ << "\n";
+                                                                        r = sqrt(*(idat.r2));
+                                                                        idat.rij = &r;
-<
+//                                                                      cerr << "idat.rij = " << *(idat.rij) << "\n";
->
+                                                                        //                                                                      cerr << "idat.rij = " << *(idat.rij) << "\n";
+                                                                        if (iLoop == PREPAIR_LOOP)
+                                                                                interactionMan_->doPair(idat);
+                                                                                fDecomp_->unpackInteractionData(idat, atom1, atom2);
-<
+//                                                                              cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
->
+                                                                                //                                                                              cerr << "d = " << *(idat.d) << "\tv=" << vpair << "\tf=" << f1 << ":" << __LINE__ << "\n";
+                                                                                vij += vpair;
+                                                                                fij += f1;
+                                                                                tau -= outProduct(*(idat.d), f1);
-<
+                }// END: omp for loop
-<
+                //              printf("after omp loop\n");
->
+                }
+                if (iLoop == PREPAIR_LOOP)

Diff Legend

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

Comparing branches/devel_omp/src/brains/ForceManager.cpp (file contents): Revision 1607 by mciznick, Wed Jul 27 14:26:53 2011 UTC vs. Revision 1608 by mciznick, Tue Aug 9 01:58:56 2011 UTC

Diff Legend

Comparing branches/devel_omp/src/brains/ForceManager.cpp (file contents):
Revision 1607 by mciznick, Wed Jul 27 14:26:53 2011 UTC vs.
Revision 1608 by mciznick, Tue Aug 9 01:58:56 2011 UTC