--- trunk/src/nonbonded/Electrostatic.cpp 2013/08/07 15:24:16 1925 +++ trunk/src/nonbonded/Electrostatic.cpp 2015/03/07 21:41:51 2071 @@ -40,6 +40,10 @@ * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ +#ifdef IS_MPI +#include +#endif + #include #include @@ -57,20 +61,31 @@ #include "math/erfc.hpp" #include "math/SquareMatrix.hpp" #include "primitives/Molecule.hpp" -#ifdef IS_MPI -#include -#endif +#include "flucq/FluctuatingChargeForces.hpp" namespace OpenMD { Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), - forceField_(NULL), info_(NULL), haveCutoffRadius_(false), haveDampingAlpha_(false), haveDielectric_(false), - haveElectroSplines_(false) - {} + haveElectroSplines_(false), + info_(NULL), forceField_(NULL) + + { + 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(); @@ -249,7 +264,7 @@ namespace OpenMD { RealType b0c, b1c, b2c, b3c, b4c, b5c; RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5; - RealType a2, expTerm, invArootPi; + RealType a2, expTerm, invArootPi(0.0); RealType r = cutoffRadius_; RealType r2 = r * r; @@ -754,6 +769,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 @@ -766,7 +783,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); @@ -829,8 +845,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; } } @@ -838,8 +855,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) { @@ -849,9 +869,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) { @@ -865,6 +887,8 @@ namespace OpenMD { } else { // only do the field if we're not excluded: Ea += C_b * pre11_ * dv01 * rhat; + Pa += C_b * pre11_ * v01; + } } @@ -872,8 +896,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) { @@ -883,22 +909,25 @@ 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); + } } - + + 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 +936,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 +1016,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 +1075,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,13 +1106,17 @@ namespace OpenMD { // + 4.0 * cross(rhat, QbQar) Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - } } if (idat.doElectricField) { *(idat.eField1) += Ea * *(idat.electroMult); *(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); @@ -1134,7 +1165,7 @@ namespace OpenMD { bool i_is_Quadrupole = data.is_Quadrupole; bool i_is_Fluctuating = data.is_Fluctuating; RealType C_a = data.fixedCharge; - RealType self(0.0), preVal, DdD, trQ, trQQ; + RealType self(0.0), preVal, DdD(0.0), trQ, trQQ; if (i_is_Dipole) { DdD = data.dipole.lengthSquare(); @@ -1142,10 +1173,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 +1280,6 @@ namespace OpenMD { vector > els; vector > ems; vector > ens; - int nMax = info_->getNAtoms(); @@ -1271,7 +1302,6 @@ namespace OpenMD { Vector3d t( 2.0 * M_PI ); t.Vdiv(t, box); - SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; int i; @@ -1441,11 +1471,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 +1491,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 +1533,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] ); } @@ -1524,4 +1557,120 @@ namespace OpenMD { } pot += kPot; } + + void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded, + RealType &spot1, RealType &spot2) { + + if (!initialized_) { + cerr << "initializing\n"; + initialize(); + cerr << "done\n"; + } + + const RealType mPoleConverter = 0.20819434; + + AtomType* atype1 = a1->getAtomType(); + AtomType* atype2 = a2->getAtomType(); + int atid1 = atype1->getIdent(); + int atid2 = atype2->getIdent(); + data1 = ElectrostaticMap[Etids[atid1]]; + data2 = ElectrostaticMap[Etids[atid2]]; + + Pa = 0.0; // Site potential at site a + Pb = 0.0; // Site potential at site b + + Vector3d d = a2->getPos() - a1->getPos(); + info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d); + RealType rij = d.length(); + // some variables we'll need independent of electrostatic type: + + RealType ri = 1.0 / rij; + rhat = d * ri; + + + if ((rij >= cutoffRadius_) || excluded) { + spot1 = 0.0; + spot2 = 0.0; + return; + } + + // logicals + + a_is_Charge = data1.is_Charge; + a_is_Dipole = data1.is_Dipole; + a_is_Quadrupole = data1.is_Quadrupole; + a_is_Fluctuating = data1.is_Fluctuating; + + b_is_Charge = data2.is_Charge; + b_is_Dipole = data2.is_Dipole; + b_is_Quadrupole = data2.is_Quadrupole; + b_is_Fluctuating = data2.is_Fluctuating; + + // Obtain all of the required radial function values from the + // spline structures: + + + if (a_is_Charge || b_is_Charge) { + v01 = v01s->getValueAt(rij); + } + if (a_is_Dipole || b_is_Dipole) { + v11 = v11s->getValueAt(rij); + v11or = ri * v11; + } + if (a_is_Quadrupole || b_is_Quadrupole) { + v21 = v21s->getValueAt(rij); + v22 = v22s->getValueAt(rij); + v22or = ri * v22; + } + + if (a_is_Charge) { + C_a = data1.fixedCharge; + + if (a_is_Fluctuating) { + C_a += a1->getFlucQPos(); + } + + Pb += C_a * pre11_ * v01; + } + + if (a_is_Dipole) { + D_a = a1->getDipole() * mPoleConverter; + rdDa = dot(rhat, D_a); + Pb += pre12_ * v11 * rdDa; + } + + if (a_is_Quadrupole) { + Q_a = a1->getQuadrupole() * mPoleConverter; + trQa = Q_a.trace(); + Qar = Q_a * rhat; + rdQar = dot(rhat, Qar); + Pb += pre14_ * (v21 * trQa + v22 * rdQar); + } + + if (b_is_Charge) { + C_b = data2.fixedCharge; + + if (b_is_Fluctuating) + C_b += a2->getFlucQPos(); + + Pa += C_b * pre11_ * v01; + } + + if (b_is_Dipole) { + D_b = a2->getDipole() * mPoleConverter; + rdDb = dot(rhat, D_b); + Pa += pre12_ * v11 * rdDb; + } + + if (b_is_Quadrupole) { + Q_a = a2->getQuadrupole() * mPoleConverter; + trQb = Q_b.trace(); + Qbr = Q_b * rhat; + rdQbr = dot(rhat, Qbr); + Pa += pre14_ * (v21 * trQb + v22 * rdQbr); + } + + spot1 = Pa; + spot2 = Pb; + } }