--- trunk/src/nonbonded/Electrostatic.cpp 2013/08/05 21:46:11 1924 +++ 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,18 +61,31 @@ #include "math/erfc.hpp" #include "math/SquareMatrix.hpp" #include "primitives/Molecule.hpp" +#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(); @@ -247,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; @@ -752,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 @@ -764,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); @@ -827,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; } } @@ -836,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) { @@ -847,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) { @@ -863,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; + } } @@ -870,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) { @@ -881,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: @@ -905,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) { @@ -986,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: @@ -1046,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); @@ -1077,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); @@ -1132,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(); @@ -1140,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_) { @@ -1194,7 +1228,7 @@ namespace OpenMD { } - void Electrostatic::ReciprocalSpaceSum(potVec& pot) { + void Electrostatic::ReciprocalSpaceSum(RealType& pot) { RealType kPot = 0.0; RealType kVir = 0.0; @@ -1240,31 +1274,39 @@ namespace OpenMD { // Calculate and store exponential factors - vector > eCos; - vector > eSin; + vector > elc; + vector > emc; + vector > enc; + vector > els; + vector > ems; + vector > ens; int nMax = info_->getNAtoms(); - eCos.resize(kLimit+1); - eSin.resize(kLimit+1); + elc.resize(kLimit+1); + emc.resize(kLimit+1); + enc.resize(kLimit+1); + els.resize(kLimit+1); + ems.resize(kLimit+1); + ens.resize(kLimit+1); + for (int j = 0; j < kLimit+1; j++) { - eCos[j].resize(nMax); - eSin[j].resize(nMax); + elc[j].resize(nMax); + emc[j].resize(nMax); + enc[j].resize(nMax); + els[j].resize(nMax); + ems[j].resize(nMax); + ens[j].resize(nMax); } Vector3d t( 2.0 * M_PI ); t.Vdiv(t, box); - SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; int i; Vector3d r; Vector3d tt; - Vector3d w; - Vector3d u; - Vector3d a; - Vector3d b; for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { @@ -1277,32 +1319,27 @@ namespace OpenMD { tt.Vmul(t, r); - - eCos[1][i] = Vector3d(1.0, 1.0, 1.0); - eSin[1][i] = Vector3d(0.0, 0.0, 0.0); - eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z())); - eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z())); + elc[1][i] = 1.0; + emc[1][i] = 1.0; + enc[1][i] = 1.0; + els[1][i] = 0.0; + ems[1][i] = 0.0; + ens[1][i] = 0.0; - u = eCos[2][i]; - w = eSin[2][i]; + elc[2][i] = cos(tt.x()); + emc[2][i] = cos(tt.y()); + enc[2][i] = cos(tt.z()); + els[2][i] = sin(tt.x()); + ems[2][i] = sin(tt.y()); + ens[2][i] = sin(tt.z()); for(int l = 3; l <= kLimit; l++) { - eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x(); - eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y(); - eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z(); - - eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x(); - eSin[l][i].y() = eSin[l-1][i].y()*eCos[2][i].y() + eCos[l-1][i].y()*eSin[2][i].y(); - eSin[l][i].z() = eSin[l-1][i].z()*eCos[2][i].z() + eCos[l-1][i].z()*eSin[2][i].z(); - - - // a.Vmul(eCos[l-1][i], u); - // b.Vmul(eSin[l-1][i], w); - // eCos[l][i] = a - b; - // a.Vmul(eSin[l-1][i], u); - // b.Vmul(eCos[l-1][i], w); - // eSin[l][i] = a + b; - + elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i]; + emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i]; + enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i]; + els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i]; + ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i]; + ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i]; } } } @@ -1346,16 +1383,26 @@ namespace OpenMD { std::vector qks(nMax, 0.0); std::vector dxk(nMax, V3Zero); std::vector qxk(nMax, V3Zero); - + RealType rl, rm, rn; + Vector3d kVec; + Vector3d Qk; + Mat3x3d k2; + RealType ckcs, ckss, dkcs, dkss, qkcs, qkss; + int atid; + ElectrostaticAtomData data; + RealType C, dk, qk; + Vector3d D; + Mat3x3d Q; + int mMin = kLimit; int nMin = kLimit + 1; for (int l = 1; l <= kLimit; l++) { int ll = l - 1; - RealType rl = xcl * float(ll); + rl = xcl * float(ll); for (int mmm = mMin; mmm <= kLim2; mmm++) { int mm = mmm - kLimit; int m = abs(mm) + 1; - RealType rm = ycl * float(mm); + rm = ycl * float(mm); // Set temporary products of exponential terms for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { @@ -1364,27 +1411,23 @@ namespace OpenMD { i = atom->getLocalIndex(); if(mm < 0) { - clm[i] = eCos[l][i].x()*eCos[m][i].y() - + eSin[l][i].x()*eSin[m][i].y(); - slm[i] = eSin[l][i].x()*eCos[m][i].y() - - eSin[m][i].y()*eCos[l][i].x(); + clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i]; + slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i]; } else { - clm[i] = eCos[l][i].x()*eCos[m][i].y() - - eSin[l][i].x()*eSin[m][i].y(); - slm[i] = eSin[l][i].x()*eCos[m][i].y() - + eSin[m][i].y()*eCos[l][i].x(); + clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i]; + slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i]; } } } for (int nnn = nMin; nnn <= kLim2; nnn++) { int nn = nnn - kLimit; int n = abs(nn) + 1; - RealType rn = zcl * float(nn); + rn = zcl * float(nn); // Test on magnitude of k vector: int kk=ll*ll + mm*mm + nn*nn; if(kk <= kSqLim) { - Vector3d kVec = Vector3d(rl, rm, rn); - Mat3x3d k2 = outProduct(kVec, kVec); + kVec = Vector3d(rl, rm, rn); + k2 = outProduct(kVec, kVec); // Calculate exp(ikr) terms for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { @@ -1393,11 +1436,12 @@ namespace OpenMD { i = atom->getLocalIndex(); if (nn < 0) { - ckr[i]=clm[i]*eCos[n][i].z()+slm[i]*eSin[n][i].z(); - skr[i]=slm[i]*eCos[n][i].z()-clm[i]*eSin[n][i].z(); + ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i]; + skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i]; + } else { - ckr[i]=clm[i]*eCos[n][i].z()-slm[i]*eSin[n][i].z(); - skr[i]=slm[i]*eCos[n][i].z()+clm[i]*eSin[n][i].z(); + ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i]; + skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i]; } } } @@ -1410,28 +1454,27 @@ namespace OpenMD { atom = mol->nextAtom(ai)) { i = atom->getLocalIndex(); int atid = atom->getAtomType()->getIdent(); - ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; + data = ElectrostaticMap[Etids[atid]]; if (data.is_Charge) { - RealType C = data.fixedCharge; + C = data.fixedCharge; if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); ckc[i] = C * ckr[i]; cks[i] = C * skr[i]; } if (data.is_Dipole) { - Vector3d D = atom->getDipole() * mPoleConverter; - RealType dk = dot(D, kVec); + D = atom->getDipole() * mPoleConverter; + dk = dot(D, kVec); dxk[i] = cross(D, kVec); dkc[i] = dk * ckr[i]; dks[i] = dk * skr[i]; } if (data.is_Quadrupole) { - Mat3x3d Q = atom->getQuadrupole(); - Q *= mPoleConverter; - RealType qk = - doubleDot(Q, k2); - // RealType qk = -( Q * k2 ).trace(); - qxk[i] = -2.0 * cross(k2, Q); + Q = atom->getQuadrupole() * mPoleConverter; + Qk = Q * kVec; + qk = dot(kVec, Qk); + qxk[i] = -cross(kVec, Qk); qkc[i] = qk * ckr[i]; qks[i] = qk * skr[i]; } @@ -1440,40 +1483,39 @@ namespace OpenMD { // calculate vector sums - RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); - RealType ckss = std::accumulate(cks.begin(),cks.end(),0.0); - RealType dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0); - RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0); - RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); - RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0); - + ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); + ckss = std::accumulate(cks.begin(),cks.end(),0.0); + dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0); + dkss = std::accumulate(dks.begin(),dks.end(),0.0); + qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); + 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: - kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) - + (ckcs-dkss-qkcs)*(ckcs-dkss-qkss)); + kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) + + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs)); - kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss - +4.0*(ckss*dkcs-ckcs*dkss) - +3.0*(dkcs*dkcs+dkss*dkss) - -6.0*(ckss*qkss+ckcs*qkcs) - +8.0*(dkss*qkcs-dkcs*qkss) - +5.0*(qkss*qkss+qkcs*qkcs)); + kVir += 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss + +4.0*(ckss*dkcs-ckcs*dkss) + +3.0*(dkcs*dkcs+dkss*dkss) + -6.0*(ckss*qkss+ckcs*qkcs) + +8.0*(dkss*qkcs-dkcs*qkss) + +5.0*(qkss*qkss+qkcs*qkcs)); // Calculate force and torque for each site: @@ -1483,18 +1525,22 @@ namespace OpenMD { atom = mol->nextAtom(ai)) { i = atom->getLocalIndex(); - int atid = atom->getAtomType()->getIdent(); - ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; - + atid = atom->getAtomType()->getIdent(); + data = ElectrostaticMap[Etids[atid]]; + RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs) - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); 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)); + RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) + +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] ); } @@ -1509,7 +1555,122 @@ namespace OpenMD { } mMin = 1; } - cerr << "kPot = " << kPot << "\n"; - pot[ELECTROSTATIC_FAMILY] += kPot; + 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; + } }