--- trunk/src/nonbonded/Electrostatic.cpp 2014/03/07 13:13:03 1974 +++ trunk/src/nonbonded/Electrostatic.cpp 2015/03/07 21:41:51 2071 @@ -66,11 +66,12 @@ namespace OpenMD { 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_); } @@ -263,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; @@ -768,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 @@ -842,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; } } @@ -851,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) { @@ -862,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) { @@ -878,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; + } } @@ -885,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) { @@ -896,9 +909,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); + } } @@ -1097,6 +1112,11 @@ namespace OpenMD { 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); @@ -1145,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(); @@ -1518,7 +1538,7 @@ namespace OpenMD { atom->addFrc( 4.0 * rvol * qfrc * kVec ); if (atom->isFluctuatingCharge()) { - atom->addFlucQFrc( 2.0 * rvol * qtrq2 ); + atom->addFlucQFrc( - 2.0 * rvol * qtrq2 ); } if (data.is_Dipole) { @@ -1536,5 +1556,121 @@ namespace OpenMD { mMin = 1; } 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; } }