--- trunk/src/brains/ForceManager.cpp 2012/11/05 19:41:28 1809 +++ trunk/src/brains/ForceManager.cpp 2015/03/05 15:35:37 2067 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ @@ -44,7 +44,6 @@ * @file ForceManager.cpp * @author tlin * @date 11/09/2004 - * @time 10:39am * @version 1.0 */ @@ -58,7 +57,8 @@ #include "primitives/Torsion.hpp" #include "primitives/Inversion.hpp" #include "nonbonded/NonBondedInteraction.hpp" -#include "perturbations/ElectricField.hpp" +#include "perturbations/UniformField.hpp" +#include "perturbations/UniformGradient.hpp" #include "parallel/ForceMatrixDecomposition.hpp" #include @@ -68,17 +68,27 @@ namespace OpenMD { using namespace std; namespace OpenMD { - ForceManager::ForceManager(SimInfo * info) : info_(info) { + ForceManager::ForceManager(SimInfo * info) : initialized_(false), info_(info), + switcher_(NULL) { forceField_ = info_->getForceField(); interactionMan_ = new InteractionManager(); fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_); + thermo = new Thermo(info_); } + ForceManager::~ForceManager() { + perturbations_.clear(); + + delete switcher_; + delete interactionMan_; + delete fDecomp_; + delete thermo; + } + /** * setupCutoffs * - * Sets the values of cutoffRadius, switchingRadius, cutoffMethod, - * and cutoffPolicy + * Sets the values of cutoffRadius, switchingRadius, and cutoffMethod * * cutoffRadius : realType * If the cutoffRadius was explicitly set, use that value. @@ -88,15 +98,11 @@ namespace OpenMD { * simulation for suggested cutoff values (e.g. 2.5 * sigma). * Use the maximum suggested value that was found. * - * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, - * or SHIFTED_POTENTIAL) + * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED, + * SHIFTED_POTENTIAL, or EWALD_FULL) * If cutoffMethod was explicitly set, use that choice. * If cutoffMethod was not explicitly set, use SHIFTED_FORCE * - * cutoffPolicy : (one of MIX, MAX, TRADITIONAL) - * If cutoffPolicy was explicitly set, use that choice. - * If cutoffPolicy was not explicitly set, use TRADITIONAL - * * switchingRadius : realType * If the cutoffMethod was set to SWITCHED: * If the switchingRadius was explicitly set, use that value @@ -109,7 +115,6 @@ namespace OpenMD { void ForceManager::setupCutoffs() { Globals* simParams_ = info_->getSimParams(); - ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); int mdFileVersion; rCut_ = 0.0; //Needs a value for a later max() call; @@ -118,6 +123,13 @@ namespace OpenMD { else mdFileVersion = 0; + // We need the list of simulated atom types to figure out cutoffs + // as well as long range corrections. + + set::iterator i; + set atomTypes_; + atomTypes_ = info_->getSimulatedAtomTypes(); + if (simParams_->haveCutoffRadius()) { rCut_ = simParams_->getCutoffRadius(); } else { @@ -132,10 +144,7 @@ namespace OpenMD { rCut_ = 12.0; } else { RealType thisCut; - set::iterator i; - set atomTypes; - atomTypes = info_->getSimulatedAtomTypes(); - for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { thisCut = interactionMan_->getSuggestedCutoffRadius((*i)); rCut_ = max(thisCut, rCut_); } @@ -149,14 +158,17 @@ namespace OpenMD { } } - fDecomp_->setUserCutoff(rCut_); + fDecomp_->setCutoffRadius(rCut_); interactionMan_->setCutoffRadius(rCut_); + rCutSq_ = rCut_ * rCut_; map stringToCutoffMethod; stringToCutoffMethod["HARD"] = HARD; stringToCutoffMethod["SWITCHED"] = SWITCHED; stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; + stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED; + stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL; if (simParams_->haveCutoffMethod()) { string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); @@ -166,7 +178,8 @@ namespace OpenMD { sprintf(painCave.errMsg, "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" "\tShould be one of: " - "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", + "HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n" + "\tSHIFTED_FORCE, or EWALD_FULL\n", cutMeth.c_str()); painCave.isFatal = 1; painCave.severity = OPENMD_ERROR; @@ -210,12 +223,17 @@ namespace OpenMD { cutoffMethod_ = SHIFTED_POTENTIAL; } else if (myMethod == "SHIFTED_FORCE") { cutoffMethod_ = SHIFTED_FORCE; + } else if (myMethod == "TAYLOR_SHIFTED") { + cutoffMethod_ = TAYLOR_SHIFTED; + } else if (myMethod == "EWALD_FULL") { + cutoffMethod_ = EWALD_FULL; } if (simParams_->haveSwitchingRadius()) rSwitch_ = simParams_->getSwitchingRadius(); - if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE" || + myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") { if (simParams_->haveSwitchingRadius()){ sprintf(painCave.errMsg, "ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n" @@ -250,49 +268,8 @@ namespace OpenMD { } } } - } - } - - map stringToCutoffPolicy; - stringToCutoffPolicy["MIX"] = MIX; - stringToCutoffPolicy["MAX"] = MAX; - stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL; - - string cutPolicy; - if (forceFieldOptions_.haveCutoffPolicy()){ - cutPolicy = forceFieldOptions_.getCutoffPolicy(); - }else if (simParams_->haveCutoffPolicy()) { - cutPolicy = simParams_->getCutoffPolicy(); - } - - if (!cutPolicy.empty()){ - toUpper(cutPolicy); - map::iterator i; - i = stringToCutoffPolicy.find(cutPolicy); - - if (i == stringToCutoffPolicy.end()) { - sprintf(painCave.errMsg, - "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n" - "\tShould be one of: " - "MIX, MAX, or TRADITIONAL\n", - cutPolicy.c_str()); - painCave.isFatal = 1; - painCave.severity = OPENMD_ERROR; - simError(); - } else { - cutoffPolicy_ = i->second; } - } else { - sprintf(painCave.errMsg, - "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n" - "\tOpenMD will use TRADITIONAL.\n"); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); - cutoffPolicy_ = TRADITIONAL; } - - fDecomp_->setCutoffPolicy(cutoffPolicy_); // create the switching function object: @@ -370,12 +347,8 @@ namespace OpenMD { } switcher_->setSwitchType(sft_); switcher_->setSwitch(rSwitch_, rCut_); - interactionMan_->setSwitchingRadius(rSwitch_); } - - - void ForceManager::initialize() { if (!info_->isTopologyDone()) { @@ -384,9 +357,9 @@ namespace OpenMD { interactionMan_->setSimInfo(info_); interactionMan_->initialize(); - // We want to delay the cutoffs until after the interaction - // manager has set up the atom-atom interactions so that we can - // query them for suggested cutoff values + //! We want to delay the cutoffs until after the interaction + //! manager has set up the atom-atom interactions so that we can + //! query them for suggested cutoff values setupCutoffs(); info_->prepareTopology(); @@ -394,19 +367,22 @@ namespace OpenMD { doParticlePot_ = info_->getSimParams()->getOutputParticlePotential(); doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux(); if (doHeatFlux_) doParticlePot_ = true; + + doElectricField_ = info_->getSimParams()->getOutputElectricField(); + doSitePotential_ = info_->getSimParams()->getOutputSitePotential(); } ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); - // Force fields can set options on how to scale van der Waals and - // electrostatic interactions for atoms connected via bonds, bends - // and torsions in this case the topological distance between - // atoms is: - // 0 = topologically unconnected - // 1 = bonded together - // 2 = connected via a bend - // 3 = connected via a torsion + //! Force fields can set options on how to scale van der Waals and + //! electrostatic interactions for atoms connected via bonds, bends + //! and torsions in this case the topological distance between + //! atoms is: + //! 0 = topologically unconnected + //! 1 = bonded together + //! 2 = connected via a bend + //! 3 = connected via a torsion vdwScale_.reserve(4); fill(vdwScale_.begin(), vdwScale_.end(), 0.0); @@ -424,21 +400,29 @@ namespace OpenMD { electrostaticScale_[2] = fopts.getelectrostatic13scale(); electrostaticScale_[3] = fopts.getelectrostatic14scale(); - if (info_->getSimParams()->haveElectricField()) { - ElectricField* eField = new ElectricField(info_); + if (info_->getSimParams()->haveUniformField()) { + UniformField* eField = new UniformField(info_); perturbations_.push_back(eField); } - - fDecomp_->distributeInitialData(); - - initialized_ = true; - + if (info_->getSimParams()->haveUniformGradientStrength() || + info_->getSimParams()->haveUniformGradientDirection1() || + info_->getSimParams()->haveUniformGradientDirection2() ) { + UniformGradient* eGrad = new UniformGradient(info_); + perturbations_.push_back(eGrad); + } + + usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions(); + + fDecomp_->distributeInitialData(); + + initialized_ = true; + } - + void ForceManager::calcForces() { if (!initialized_) initialize(); - + preCalculation(); shortRangeInteractions(); longRangeInteractions(); @@ -612,14 +596,15 @@ namespace OpenMD { // Collect from all nodes. This should eventually be moved into a // SystemDecomposition, but this is a better place than in // Thermo to do the collection. - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1, - MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1, - MPI::REALTYPE, MPI::SUM); + + MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); @@ -637,7 +622,6 @@ namespace OpenMD { void ForceManager::longRangeInteractions() { - Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); DataStorage* config = &(curSnapshot->atomData); DataStorage* cgConfig = &(curSnapshot->cgData); @@ -648,8 +632,8 @@ namespace OpenMD { Molecule* mol; Molecule::CutoffGroupIterator ci; CutoffGroup* cg; - - if(info_->getNCutoffGroups() > 0){ + + if(info_->getNCutoffGroups() != info_->getNAtoms()){ for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { for(cg = mol->beginCutoffGroup(ci); cg != NULL; @@ -673,26 +657,30 @@ namespace OpenMD { RealType electroMult, vdwMult; RealType vij; Vector3d fij, fg, f1; - tuple3 cuts; - RealType rCutSq; bool in_switching_region; RealType sw, dswdr, swderiv; - vector atomListColumn, atomListRow, atomListLocal; + vector atomListColumn, atomListRow; InteractionData idat; SelfData sdat; RealType mf; RealType vpair; RealType dVdFQ1(0.0); RealType dVdFQ2(0.0); - Vector3d eField1(0.0); - Vector3d eField2(0.0); potVec longRangePotential(0.0); + RealType reciprocalPotential(0.0); potVec workPot(0.0); potVec exPot(0.0); + Vector3d eField1(0.0); + Vector3d eField2(0.0); + RealType sPot1(0.0); + RealType sPot2(0.0); + bool newAtom1; + vector::iterator ia, jb; int loopStart, loopEnd; - + + idat.rcut = &rCut_; idat.vdwMult = &vdwMult; idat.electroMult = &electroMult; idat.pot = &workPot; @@ -703,12 +691,17 @@ namespace OpenMD { idat.dVdFQ1 = &dVdFQ1; idat.dVdFQ2 = &dVdFQ2; idat.eField1 = &eField1; - idat.eField2 = &eField2; + idat.eField2 = &eField2; + idat.sPot1 = &sPot1; + idat.sPot2 = &sPot2; idat.f1 = &f1; idat.sw = &sw; idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false; - idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false; + idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || + cutoffMethod_ == TAYLOR_SHIFTED) ? true : false; idat.doParticlePot = doParticlePot_; + idat.doElectricField = doElectricField_; + idat.doSitePotential = doSitePotential_; sdat.doParticlePot = doParticlePot_; loopEnd = PAIR_LOOP; @@ -721,179 +714,192 @@ namespace OpenMD { if (iLoop == loopStart) { bool update_nlist = fDecomp_->checkNeighborList(); - if (update_nlist) - neighborList = fDecomp_->buildNeighborList(); - } - - for (vector >::iterator it = neighborList.begin(); - it != neighborList.end(); ++it) { - - cg1 = (*it).first; - cg2 = (*it).second; + if (update_nlist) { + if (!usePeriodicBoundaryConditions_) + Mat3x3d bbox = thermo->getBoundingBox(); + fDecomp_->buildNeighborList(neighborList_, point_); + } + } + + for (cg1 = 0; cg1 < int(point_.size()) - 1; cg1++) { - cuts = fDecomp_->getGroupCutoffs(cg1, cg2); + atomListRow = fDecomp_->getAtomsInGroupRow(cg1); + newAtom1 = true; + + for (int m2 = point_[cg1]; m2 < point_[cg1+1]; m2++) { - 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; - } + cg2 = neighborList_[m2]; - in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, - rgrp); - - atomListRow = fDecomp_->getAtomsInGroupRow(cg1); - atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2); - - if (doHeatFlux_) - gvel2 = fDecomp_->getGroupVelocityColumn(cg2); - - for (ia = atomListRow.begin(); - ia != atomListRow.end(); ++ia) { - atom1 = (*ia); - - for (jb = atomListColumn.begin(); - jb != atomListColumn.end(); ++jb) { - atom2 = (*jb); - - if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) { - - vpair = 0.0; - workPot = 0.0; - exPot = 0.0; - f1 = V3Zero; - dVdFQ1 = 0.0; - dVdFQ2 = 0.0; - - 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; - if (doHeatFlux_) - vel2 = gvel2; - } else { - d = fDecomp_->getInteratomicVector(atom1, atom2); - curSnapshot->wrapVector( d ); - r2 = d.lengthSquare(); - idat.d = &d; - idat.r2 = &r2; - if (doHeatFlux_) - vel2 = fDecomp_->getAtomVelocityColumn(atom2); + d_grp = fDecomp_->getIntergroupVector(cg1, cg2); + + // already wrapped in the getIntergroupVector call: + // curSnapshot->wrapVector(d_grp); + rgrpsq = d_grp.lengthSquare(); + + if (rgrpsq < rCutSq_) { + if (iLoop == PAIR_LOOP) { + vij = 0.0; + fij.zero(); + eField1.zero(); + eField2.zero(); + sPot1 = 0.0; + sPot2 = 0.0; + } + + in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, + rgrp); + + atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2); + + if (doHeatFlux_) + gvel2 = fDecomp_->getGroupVelocityColumn(cg2); + + for (ia = atomListRow.begin(); + ia != atomListRow.end(); ++ia) { + atom1 = (*ia); + + for (jb = atomListColumn.begin(); + jb != atomListColumn.end(); ++jb) { + atom2 = (*jb); + + if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) { + + vpair = 0.0; + workPot = 0.0; + exPot = 0.0; + f1.zero(); + dVdFQ1 = 0.0; + dVdFQ2 = 0.0; + + fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1); + + 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; + if (doHeatFlux_) + vel2 = gvel2; + } else { + d = fDecomp_->getInteratomicVector(atom1, atom2); + curSnapshot->wrapVector( d ); + r2 = d.lengthSquare(); + idat.d = &d; + idat.r2 = &r2; + if (doHeatFlux_) + vel2 = fDecomp_->getAtomVelocityColumn(atom2); + } + + r = sqrt( *(idat.r2) ); + idat.rij = &r; + + if (iLoop == PREPAIR_LOOP) { + interactionMan_->doPrePair(idat); + } else { + interactionMan_->doPair(idat); + fDecomp_->unpackInteractionData(idat, atom1, atom2); + vij += vpair; + fij += f1; + stressTensor -= outProduct( *(idat.d), f1); + if (doHeatFlux_) + fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2)); + } } - - r = sqrt( *(idat.r2) ); - idat.rij = &r; - - if (iLoop == PREPAIR_LOOP) { - interactionMan_->doPrePair(idat); - } else { - interactionMan_->doPair(idat); - fDecomp_->unpackInteractionData(idat, atom1, atom2); - vij += vpair; - fij += f1; - stressTensor -= outProduct( *(idat.d), f1); - if (doHeatFlux_) - fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2)); - } } } - } - - 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) { - if (!fDecomp_->skipAtomPair(atomListRow[0], - atomListColumn[0], - cg1, cg2)) { + + 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) { + if (!fDecomp_->skipAtomPair(atomListRow[0], + atomListColumn[0], + cg1, cg2)) { stressTensor -= outProduct( *(idat.d), fg); if (doHeatFlux_) fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2)); - } - } - - for (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); - stressTensor -= outProduct(dag, fg); - if (doHeatFlux_) - fDecomp_->addToHeatFlux( dag * dot(fg, vel2)); + } + } + + for (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); + stressTensor -= outProduct(dag, fg); + if (doHeatFlux_) + fDecomp_->addToHeatFlux( dag * dot(fg, vel2)); + } } } - } - for (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); - stressTensor -= outProduct(dag, fg); - if (doHeatFlux_) - fDecomp_->addToHeatFlux( dag * dot(fg, vel2)); + for (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); + stressTensor -= outProduct(dag, fg); + if (doHeatFlux_) + fDecomp_->addToHeatFlux( dag * dot(fg, vel2)); + } } } } + //if (!info_->usesAtomicVirial()) { + // stressTensor -= outProduct(d_grp, fij); + // if (doHeatFlux_) + // fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2)); + //} } - //if (!info_->usesAtomicVirial()) { - // stressTensor -= outProduct(d_grp, fij); - // if (doHeatFlux_) - // fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2)); - //} } } + newAtom1 = false; } - + if (iLoop == PREPAIR_LOOP) { if (info_->requiresPrepair()) { - + fDecomp_->collectIntermediateData(); - + for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { fDecomp_->fillSelfData(sdat, atom1); interactionMan_->doPreForce(sdat); } - + fDecomp_->distributeIntermediateData(); - + } } } - + // collects pairwise information fDecomp_->collectData(); + if (cutoffMethod_ == EWALD_FULL) { + interactionMan_->doReciprocalSpaceSum(reciprocalPotential); + + curSnapshot->setReciprocalPotential(reciprocalPotential); + } if (info_->requiresSelfCorrection()) { for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { @@ -911,11 +917,10 @@ namespace OpenMD { curSnapshot->setLongRangePotential(longRangePotential); curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) + - *(fDecomp_->getExcludedPotential())); + *(fDecomp_->getExcludedPotential())); } - void ForceManager::postCalculation() { vector::iterator pi; @@ -941,10 +946,46 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9, - MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif curSnapshot->setStressTensor(stressTensor); + if (info_->getSimParams()->getUseLongRangeCorrections()) { + /* + RealType vol = curSnapshot->getVolume(); + RealType Elrc(0.0); + RealType Wlrc(0.0); + + set::iterator i; + set::iterator j; + + RealType n_i, n_j; + RealType rho_i, rho_j; + pair LRI; + + for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { + n_i = RealType(info_->getGlobalCountOfType(*i)); + rho_i = n_i / vol; + for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) { + n_j = RealType(info_->getGlobalCountOfType(*j)); + rho_j = n_j / vol; + + LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) ); + + Elrc += n_i * rho_j * LRI.first; + Wlrc -= rho_i * rho_j * LRI.second; + } + } + Elrc *= 2.0 * NumericConstant::PI; + Wlrc *= 2.0 * NumericConstant::PI; + + RealType lrp = curSnapshot->getLongRangePotential(); + curSnapshot->setLongRangePotential(lrp + Elrc); + stressTensor += Wlrc * SquareMatrix3::identity(); + curSnapshot->setStressTensor(stressTensor); + */ + + } } -} //end namespace OpenMD +}