--- trunk/src/nonbonded/Electrostatic.cpp 2013/08/23 15:59:23 1933 +++ trunk/src/nonbonded/Electrostatic.cpp 2014/04/29 17:32:31 1993 @@ -40,6 +40,10 @@ * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ +#ifdef IS_MPI +#include +#endif + #include #include @@ -57,9 +61,7 @@ #include "math/erfc.hpp" #include "math/SquareMatrix.hpp" #include "primitives/Molecule.hpp" -#ifdef IS_MPI -#include -#endif +#include "flucq/FluctuatingChargeForces.hpp" namespace OpenMD { @@ -69,8 +71,20 @@ namespace OpenMD { haveDampingAlpha_(false), haveDielectric_(false), haveElectroSplines_(false) - {} + { + flucQ_ = new FluctuatingChargeForces(info_); + } + void Electrostatic::setForceField(ForceField *ff) { + forceField_ = ff; + flucQ_->setForceField(forceField_); + } + + void Electrostatic::setSimulatedAtomTypes(set &simtypes) { + simTypes_ = simtypes; + flucQ_->setSimulatedAtomTypes(simTypes_); + } + void Electrostatic::initialize() { Globals* simParams_ = info_->getSimParams(); @@ -754,6 +768,8 @@ namespace OpenMD { Tb.zero(); // Torque on site b Ea.zero(); // Electric field at site a Eb.zero(); // Electric field at site b + Pa = 0.0; // Site potential at site a + Pb = 0.0; // Site potential at site b dUdCa = 0.0; // fluctuating charge force at site a dUdCb = 0.0; // fluctuating charge force at site a @@ -828,8 +844,9 @@ namespace OpenMD { if (idat.excluded) { *(idat.skippedCharge2) += C_a; } else { - // only do the field if we're not excluded: + // only do the field and site potentials if we're not excluded: Eb -= C_a * pre11_ * dv01 * rhat; + Pb += C_a * pre11_ * v01; } } @@ -837,8 +854,11 @@ namespace OpenMD { D_a = *(idat.dipole1); rdDa = dot(rhat, D_a); rxDa = cross(rhat, D_a); - if (!idat.excluded) + if (!idat.excluded) { Eb -= pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a); + Pb += pre12_ * v11 * rdDa; + } + } if (a_is_Quadrupole) { @@ -848,9 +868,11 @@ namespace OpenMD { rQa = rhat * Q_a; rdQar = dot(rhat, Qar); rxQar = cross(rhat, Qar); - if (!idat.excluded) + if (!idat.excluded) { Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or + rdQar * rhat * (dv22 - 2.0*v22or)); + Pb += pre14_ * (v21 * trQa + v22 * rdQar); + } } if (b_is_Charge) { @@ -864,6 +886,8 @@ namespace OpenMD { } else { // only do the field if we're not excluded: Ea += C_b * pre11_ * dv01 * rhat; + Pa += C_b * pre11_ * v01; + } } @@ -871,8 +895,10 @@ namespace OpenMD { D_b = *(idat.dipole2); rdDb = dot(rhat, D_b); rxDb = cross(rhat, D_b); - if (!idat.excluded) + if (!idat.excluded) { Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b); + Pa += pre12_ * v11 * rdDb; + } } if (b_is_Quadrupole) { @@ -882,9 +908,11 @@ namespace OpenMD { rQb = rhat * Q_b; rdQbr = dot(rhat, Qbr); rxQbr = cross(rhat, Qbr); - if (!idat.excluded) + if (!idat.excluded) { Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or + rdQbr * rhat * (dv22 - 2.0*v22or)); + Pa += pre14_ * (v21 * trQb + v22 * rdQbr); + } } @@ -1085,6 +1113,11 @@ namespace OpenMD { *(idat.eField2) += Eb * *(idat.electroMult); } + if (idat.doSitePotential) { + *(idat.sPot1) += Pa * *(idat.electroMult); + *(idat.sPot2) += Pb * *(idat.electroMult); + } + if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw); if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw); @@ -1139,9 +1172,11 @@ namespace OpenMD { if (i_is_Fluctuating) { C_a += *(sdat.flucQ); - *(sdat.flucQfrc) -= *(sdat.flucQ) * data.hardness + data.electronegativity; - (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * - (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); + + flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ), + (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY], + *(sdat.flucQfrc) ); + } switch (summationMethod_) { @@ -1244,7 +1279,6 @@ namespace OpenMD { vector > els; vector > ems; vector > ens; - int nMax = info_->getNAtoms(); @@ -1267,7 +1301,6 @@ namespace OpenMD { Vector3d t( 2.0 * M_PI ); t.Vdiv(t, box); - SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; int i; @@ -1440,7 +1473,7 @@ namespace OpenMD { Q = atom->getQuadrupole() * mPoleConverter; Qk = Q * kVec; qk = dot(kVec, Qk); - qxk[i] = cross(kVec, Qk); + qxk[i] = -cross(kVec, Qk); qkc[i] = qk * ckr[i]; qks[i] = qk * skr[i]; } @@ -1457,18 +1490,18 @@ namespace OpenMD { qkss = std::accumulate(qks.begin(),qks.end(),0.0); #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE, - MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE, - MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); #endif // Accumulate potential energy and virial contribution: @@ -1502,7 +1535,11 @@ namespace OpenMD { +skr[i]*(ckss+dkcs-qkss)); atom->addFrc( 4.0 * rvol * qfrc * kVec ); - + + if (atom->isFluctuatingCharge()) { + atom->addFlucQFrc( - 2.0 * rvol * qtrq2 ); + } + if (data.is_Dipole) { atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] ); }