--- trunk/src/brains/ForceManager.cpp 2013/06/17 18:28:30 1880 +++ trunk/src/brains/ForceManager.cpp 2014/10/31 18:57:19 2032 @@ -57,7 +57,7 @@ #include "primitives/Torsion.hpp" #include "primitives/Inversion.hpp" #include "nonbonded/NonBondedInteraction.hpp" -#include "perturbations/ElectricField.hpp" +#include "perturbations/UniformField.hpp" #include "parallel/ForceMatrixDecomposition.hpp" #include @@ -99,7 +99,7 @@ namespace OpenMD { * Use the maximum suggested value that was found. * * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED, - * or SHIFTED_POTENTIAL) + * SHIFTED_POTENTIAL, or EWALD_FULL) * If cutoffMethod was explicitly set, use that choice. * If cutoffMethod was not explicitly set, use SHIFTED_FORCE * @@ -172,6 +172,7 @@ namespace OpenMD { 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()); @@ -182,7 +183,7 @@ namespace OpenMD { "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" "\tShould be one of: " "HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n" - "\tor SHIFTED_FORCE\n", + "\tSHIFTED_FORCE, or EWALD_FULL\n", cutMeth.c_str()); painCave.isFatal = 1; painCave.severity = OPENMD_ERROR; @@ -228,13 +229,15 @@ namespace OpenMD { 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" || - myMethod == "TAYLOR_SHIFTED") { + myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") { if (simParams_->haveSwitchingRadius()){ sprintf(painCave.errMsg, "ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n" @@ -402,9 +405,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(); @@ -414,19 +417,20 @@ namespace OpenMD { 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); @@ -444,8 +448,8 @@ 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); } @@ -634,14 +638,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(); @@ -663,6 +668,7 @@ namespace OpenMD { DataStorage* config = &(curSnapshot->atomData); DataStorage* cgConfig = &(curSnapshot->cgData); + //calculate the center of mass of cutoff group SimInfo::MoleculeIterator mi; @@ -695,7 +701,7 @@ namespace OpenMD { RealType vij; Vector3d fij, fg, f1; tuple3 cuts; - RealType rCutSq; + RealType rCut, rCutSq, rListSq; bool in_switching_region; RealType sw, dswdr, swderiv; vector atomListColumn, atomListRow; @@ -706,14 +712,19 @@ namespace OpenMD { RealType dVdFQ1(0.0); RealType dVdFQ2(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); + vector::iterator ia, jb; int loopStart, loopEnd; - + + idat.rcut = &rCut; idat.vdwMult = &vdwMult; idat.electroMult = &electroMult; idat.pot = &workPot; @@ -724,13 +735,16 @@ 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 || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false; idat.doParticlePot = doParticlePot_; idat.doElectricField = doElectricField_; + idat.doSitePotential = doSitePotential_; sdat.doParticlePot = doParticlePot_; loopEnd = PAIR_LOOP; @@ -746,31 +760,32 @@ namespace OpenMD { if (update_nlist) { if (!usePeriodicBoundaryConditions_) Mat3x3d bbox = thermo->getBoundingBox(); - neighborList = fDecomp_->buildNeighborList(); + fDecomp_->buildNeighborList(neighborList_); } } - for (vector >::iterator it = neighborList.begin(); - it != neighborList.end(); ++it) { + for (vector >::iterator it = neighborList_.begin(); + it != neighborList_.end(); ++it) { cg1 = (*it).first; cg2 = (*it).second; - cuts = fDecomp_->getGroupCutoffs(cg1, cg2); + fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq); d_grp = fDecomp_->getIntergroupVector(cg1, cg2); - curSnapshot->wrapVector(d_grp); + // already wrapped in the getIntergroupVector call: + // 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.zero(); eField1.zero(); eField2.zero(); + sPot1 = 0.0; + sPot2 = 0.0; } in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, @@ -822,7 +837,7 @@ namespace OpenMD { r = sqrt( *(idat.r2) ); idat.rij = &r; - + if (iLoop == PREPAIR_LOOP) { interactionMan_->doPrePair(idat); } else { @@ -921,6 +936,11 @@ namespace OpenMD { // 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++) { @@ -942,7 +962,6 @@ namespace OpenMD { } - void ForceManager::postCalculation() { vector::iterator pi; @@ -968,8 +987,8 @@ 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);