--- trunk/src/nonbonded/Electrostatic.cpp 2013/08/05 21:46:11 1924 +++ 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,8 +61,8 @@ #include "math/erfc.hpp" #include "math/SquareMatrix.hpp" #include "primitives/Molecule.hpp" +#include "flucq/FluctuatingChargeForces.hpp" - namespace OpenMD { Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), @@ -67,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(); @@ -763,7 +779,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: @@ -885,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: @@ -905,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) { @@ -986,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: @@ -1046,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); @@ -1077,7 +1091,6 @@ namespace OpenMD { // + 4.0 * cross(rhat, QbQar) Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - } } @@ -1140,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_) { @@ -1194,7 +1208,7 @@ namespace OpenMD { } - void Electrostatic::ReciprocalSpaceSum(potVec& pot) { + void Electrostatic::ReciprocalSpaceSum(RealType& pot) { RealType kPot = 0.0; RealType kVir = 0.0; @@ -1240,31 +1254,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 +1299,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 +1363,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 +1391,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 +1416,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 +1434,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 +1463,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 +1505,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 +1535,6 @@ namespace OpenMD { } mMin = 1; } - cerr << "kPot = " << kPot << "\n"; - pot[ELECTROSTATIC_FAMILY] += kPot; + pot += kPot; } }