--- trunk/src/brains/ForceManager.cpp 2013/07/29 17:55:17 1915 +++ trunk/src/brains/ForceManager.cpp 2014/11/01 14:12:16 2033 @@ -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 @@ -393,10 +393,7 @@ namespace OpenMD { switcher_->setSwitchType(sft_); switcher_->setSwitch(rSwitch_, rCut_); } - - - void ForceManager::initialize() { if (!info_->isTopologyDone()) { @@ -405,9 +402,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(); @@ -417,19 +414,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); @@ -447,8 +445,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); } @@ -637,14 +635,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(); @@ -666,6 +665,7 @@ namespace OpenMD { DataStorage* config = &(curSnapshot->atomData); DataStorage* cgConfig = &(curSnapshot->cgData); + //calculate the center of mass of cutoff group SimInfo::MoleculeIterator mi; @@ -709,10 +709,14 @@ 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; @@ -728,13 +732,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 || cutoffMethod_ == TAYLOR_SHIFTED) ? 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; @@ -755,7 +763,7 @@ namespace OpenMD { } for (vector >::iterator it = neighborList_.begin(); - it != neighborList_.end(); ++it) { + it != neighborList_.end(); ++it) { cg1 = (*it).first; cg2 = (*it).second; @@ -769,12 +777,13 @@ namespace OpenMD { 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, @@ -826,7 +835,7 @@ namespace OpenMD { r = sqrt( *(idat.r2) ); idat.rij = &r; - + if (iLoop == PREPAIR_LOOP) { interactionMan_->doPrePair(idat); } else { @@ -926,7 +935,9 @@ namespace OpenMD { // collects pairwise information fDecomp_->collectData(); if (cutoffMethod_ == EWALD_FULL) { - interactionMan_->doReciprocalSpaceSum(); + interactionMan_->doReciprocalSpaceSum(reciprocalPotential); + + curSnapshot->setReciprocalPotential(reciprocalPotential); } if (info_->requiresSelfCorrection()) { @@ -945,11 +956,10 @@ namespace OpenMD { curSnapshot->setLongRangePotential(longRangePotential); curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) + - *(fDecomp_->getExcludedPotential())); + *(fDecomp_->getExcludedPotential())); } - void ForceManager::postCalculation() { vector::iterator pi; @@ -975,44 +985,44 @@ 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); + RealType vol = curSnapshot->getVolume(); + RealType Elrc(0.0); + RealType Wlrc(0.0); - set::iterator i; - set::iterator j; + set::iterator i; + set::iterator j; - RealType n_i, n_j; - RealType rho_i, rho_j; - pair LRI; + RealType n_i, n_j; + RealType rho_i, rho_j; + pair LRI; - for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { + 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; + n_j = RealType(info_->getGlobalCountOfType(*j)); + rho_j = n_j / vol; - LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) ); + LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) ); - Elrc += n_i * rho_j * LRI.first; - Wlrc -= rho_i * rho_j * LRI.second; + Elrc += n_i * rho_j * LRI.first; + Wlrc -= rho_i * rho_j * LRI.second; } - } - Elrc *= 2.0 * NumericConstant::PI; - Wlrc *= 2.0 * NumericConstant::PI; + } + 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); + RealType lrp = curSnapshot->getLongRangePotential(); + curSnapshot->setLongRangePotential(lrp + Elrc); + stressTensor += Wlrc * SquareMatrix3::identity(); + curSnapshot->setStressTensor(stressTensor); */ }