--- trunk/src/nonbonded/Electrostatic.cpp 2013/08/07 15:24:16 1925 +++ trunk/src/nonbonded/Electrostatic.cpp 2014/04/14 18:32:51 1981 @@ -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(); @@ -766,7 +780,6 @@ namespace OpenMD { // Excluded potential that is still computed for fluctuating charges excluded_Pot= 0.0; - // some variables we'll need independent of electrostatic type: ri = 1.0 / *(idat.rij); @@ -887,18 +900,19 @@ namespace OpenMD { Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or + rdQbr * rhat * (dv22 - 2.0*v22or)); } - + + if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; } - + if (a_is_Charge) { if (b_is_Charge) { pref = pre11_ * *(idat.electroMult); U += C_a * C_b * pref * v01; F += C_a * C_b * pref * dv01 * rhat; - + // If this is an excluded pair, there are still indirect // interactions via the reaction field we must worry about: @@ -907,22 +921,21 @@ namespace OpenMD { indirect_Pot += rfContrib; indirect_F += rfContrib * 2.0 * ri * rhat; } - + // Fluctuating charge forces are handled via Coulomb integrals // for excluded pairs (i.e. those connected via bonds) and // with the standard charge-charge interaction otherwise. - if (idat.excluded) { + if (idat.excluded) { if (a_is_Fluctuating || b_is_Fluctuating) { coulInt = J->getValueAt( *(idat.rij) ); - if (a_is_Fluctuating) dUdCa += coulInt * C_b; - if (b_is_Fluctuating) dUdCb += coulInt * C_a; - excluded_Pot += C_a * C_b * coulInt; - } + if (a_is_Fluctuating) dUdCa += C_b * coulInt; + if (b_is_Fluctuating) dUdCb += C_a * coulInt; + } } else { if (a_is_Fluctuating) dUdCa += C_b * pref * v01; - if (a_is_Fluctuating) dUdCb += C_a * pref * v01; - } + if (b_is_Fluctuating) dUdCb += C_a * pref * v01; + } } if (b_is_Dipole) { @@ -988,7 +1001,6 @@ namespace OpenMD { F -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat; Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb); - // Even if we excluded this pair from direct interactions, we // still have the reaction-field-mediated dipole-dipole // interaction: @@ -1048,7 +1060,7 @@ namespace OpenMD { trQaQb = QaQb.trace(); rQaQb = rhat * QaQb; QaQbr = QaQb * rhat; - QaxQb = cross(Q_a, Q_b); + QaxQb = mCross(Q_a, Q_b); rQaQbr = dot(rQa, Qbr); rQaxQbr = cross(rQa, Qbr); @@ -1079,7 +1091,6 @@ namespace OpenMD { // + 4.0 * cross(rhat, QbQar) Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - } } @@ -1142,10 +1153,11 @@ namespace OpenMD { if (i_is_Fluctuating) { C_a += *(sdat.flucQ); - // dVdFQ is really a force, so this is negative the derivative - *(sdat.dVdFQ) -= *(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_) { @@ -1248,7 +1260,6 @@ namespace OpenMD { vector > els; vector > ems; vector > ens; - int nMax = info_->getNAtoms(); @@ -1271,7 +1282,6 @@ namespace OpenMD { Vector3d t( 2.0 * M_PI ); t.Vdiv(t, box); - SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; int i; @@ -1441,11 +1451,10 @@ namespace OpenMD { dks[i] = dk * skr[i]; } if (data.is_Quadrupole) { - Q = atom->getQuadrupole(); - Q *= mPoleConverter; - Qk = Q * kVec; + 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]; } @@ -1462,18 +1471,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: @@ -1504,10 +1513,14 @@ namespace OpenMD { RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) -ckr[i]*(ckss+dkcs-qkss)); RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) - +skr[i]*(ckss+dkcs-qkss)); + +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] ); }