--- trunk/src/brains/ForceManager.cpp 2013/07/01 21:09:37 1895 +++ trunk/src/brains/ForceManager.cpp 2013/08/07 15:24:16 1925 @@ -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" @@ -662,6 +665,7 @@ namespace OpenMD { Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); DataStorage* config = &(curSnapshot->atomData); DataStorage* cgConfig = &(curSnapshot->cgData); + //calculate the center of mass of cutoff group @@ -695,7 +699,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,6 +710,7 @@ 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); @@ -714,6 +719,7 @@ namespace OpenMD { int loopStart, loopEnd; + idat.rcut = &rCut; idat.vdwMult = &vdwMult; idat.electroMult = &electroMult; idat.pot = &workPot; @@ -756,17 +762,15 @@ namespace OpenMD { cg1 = (*it).first; cg2 = (*it).second; - cuts = fDecomp_->getGroupCutoffs(cg1, cg2); + fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq); d_grp = fDecomp_->getIntergroupVector(cg1, cg2); // 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(); @@ -823,7 +827,7 @@ namespace OpenMD { r = sqrt( *(idat.r2) ); idat.rij = &r; - + if (iLoop == PREPAIR_LOOP) { interactionMan_->doPrePair(idat); } else { @@ -922,6 +926,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++) { @@ -943,7 +952,6 @@ namespace OpenMD { } - void ForceManager::postCalculation() { vector::iterator pi;