--- trunk/src/nonbonded/Electrostatic.cpp 2013/07/12 17:38:06 1900 +++ trunk/src/nonbonded/Electrostatic.cpp 2014/04/14 18:32:51 1981 @@ -40,10 +40,15 @@ * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ +#ifdef IS_MPI +#include +#endif + #include #include #include +#include #include "nonbonded/Electrostatic.hpp" #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" @@ -55,6 +60,8 @@ #include "utils/PhysicalConstants.hpp" #include "math/erfc.hpp" #include "math/SquareMatrix.hpp" +#include "primitives/Molecule.hpp" +#include "flucq/FluctuatingChargeForces.hpp" namespace OpenMD { @@ -64,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(); @@ -191,13 +210,12 @@ namespace OpenMD { simError(); } - if (screeningMethod_ == DAMPED) { + if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) { if (!simParams_->haveDampingAlpha()) { // first set a cutoff dependent alpha value // we assume alpha depends linearly with rcut from 0 to 20.5 ang dampingAlpha_ = 0.425 - cutoffRadius_* 0.02; - if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0; - + if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0; // throw warning sprintf( painCave.errMsg, "Electrostatic::initialize: dampingAlpha was not specified in the\n" @@ -213,6 +231,7 @@ namespace OpenMD { haveDampingAlpha_ = true; } + Etypes.clear(); Etids.clear(); FQtypes.clear(); @@ -262,9 +281,6 @@ namespace OpenMD { b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; - //selfMult1_ = - 2.0 * a2 * invArootPi; - //selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0; - //selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0; // Half the Smith self piece: selfMult1_ = - a2 * invArootPi; selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0; @@ -288,11 +304,13 @@ namespace OpenMD { db0c_3 = 3.0*r*b2c - r2*r*b3c; db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; - - selfMult1_ -= b0c; - selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; - selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; + if (summationMethod_ != esm_EWALD_FULL) { + selfMult1_ -= b0c; + selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; + selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; + } + // working variables for the splines: RealType ri, ri2; RealType b0, b1, b2, b3, b4, b5; @@ -328,14 +346,6 @@ namespace OpenMD { vector v31v, v32v; vector v41v, v42v, v43v; - /* - vector dv01v; - vector dv11v; - vector dv21v, dv22v; - vector dv31v, dv32v; - vector dv41v, dv42v, dv43v; - */ - for (int i = 1; i < np_ + 1; i++) { r = RealType(i) * dx; rv.push_back(r); @@ -499,6 +509,7 @@ namespace OpenMD { case esm_SWITCHING_FUNCTION: case esm_HARD: + case esm_EWALD_FULL: v01 = f; v11 = g; @@ -557,7 +568,6 @@ namespace OpenMD { break; - case esm_EWALD_FULL: case esm_EWALD_PME: case esm_EWALD_SPME: default : @@ -586,17 +596,6 @@ namespace OpenMD { v41v.push_back(v41); v42v.push_back(v42); v43v.push_back(v43); - /* - dv01v.push_back(dv01); - dv11v.push_back(dv11); - dv21v.push_back(dv21); - dv22v.push_back(dv22); - dv31v.push_back(dv31); - dv32v.push_back(dv32); - dv41v.push_back(dv41); - dv42v.push_back(dv42); - dv43v.push_back(dv43); - */ } // construct the spline structures and fill them with the values we've @@ -621,27 +620,6 @@ namespace OpenMD { v43s = new CubicSpline(); v43s->addPoints(rv, v43v); - /* - dv01s = new CubicSpline(); - dv01s->addPoints(rv, dv01v); - dv11s = new CubicSpline(); - dv11s->addPoints(rv, dv11v); - dv21s = new CubicSpline(); - dv21s->addPoints(rv, dv21v); - dv22s = new CubicSpline(); - dv22s->addPoints(rv, dv22v); - dv31s = new CubicSpline(); - dv31s->addPoints(rv, dv31v); - dv32s = new CubicSpline(); - dv32s->addPoints(rv, dv32v); - dv41s = new CubicSpline(); - dv41s->addPoints(rv, dv41v); - dv42s = new CubicSpline(); - dv42s->addPoints(rv, dv42v); - dv43s = new CubicSpline(); - dv43s->addPoints(rv, dv43v); - */ - haveElectroSplines_ = true; initialized_ = true; @@ -715,7 +693,8 @@ namespace OpenMD { FQtids[atid] = fqtid; Jij[fqtid].resize(nFlucq_); - // Now, iterate over all known fluctuating and add to the coulomb integral map: + // Now, iterate over all known fluctuating and add to the + // coulomb integral map: std::set::iterator it; for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { @@ -801,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); @@ -922,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: @@ -942,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) { @@ -1023,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: @@ -1083,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); @@ -1114,7 +1091,6 @@ namespace OpenMD { // + 4.0 * cross(rhat, QbQar) Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - } } @@ -1177,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_) { @@ -1202,6 +1179,8 @@ namespace OpenMD { case esm_SHIFTED_FORCE: case esm_SHIFTED_POTENTIAL: + case esm_TAYLOR_SHIFTED: + case esm_EWALD_FULL: if (i_is_Charge) self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); if (i_is_Dipole) @@ -1226,5 +1205,336 @@ namespace OpenMD { // 12 angstroms seems to be a reasonably good guess for most // cases. return 12.0; + } + + + void Electrostatic::ReciprocalSpaceSum(RealType& pot) { + + RealType kPot = 0.0; + RealType kVir = 0.0; + + const RealType mPoleConverter = 0.20819434; // converts from the + // internal units of + // Debye (for dipoles) + // or Debye-angstroms + // (for quadrupoles) to + // electron angstroms or + // electron-angstroms^2 + + const RealType eConverter = 332.0637778; // convert the + // Charge-Charge + // electrostatic + // interactions into kcal / + // mol assuming distances + // are measured in + // angstroms. + + Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat(); + Vector3d box = hmat.diagonals(); + RealType boxMax = box.max(); + + //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); + int kMax = 7; + int kSqMax = kMax*kMax + 2; + + int kLimit = kMax+1; + int kLim2 = 2*kMax+1; + int kSqLim = kSqMax; + + vector AK(kSqLim+1, 0.0); + RealType xcl = 2.0 * M_PI / box.x(); + RealType ycl = 2.0 * M_PI / box.y(); + RealType zcl = 2.0 * M_PI / box.z(); + RealType rcl = 2.0 * M_PI / boxMax; + RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z()); + + if(dampingAlpha_ < 1.0e-12) return; + + RealType ralph = -0.25/pow(dampingAlpha_,2); + + // Calculate and store exponential factors + + vector > elc; + vector > emc; + vector > enc; + vector > els; + vector > ems; + vector > ens; + + int nMax = info_->getNAtoms(); + + 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++) { + 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; + + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + i = atom->getLocalIndex(); + r = atom->getPos(); + info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); + + tt.Vmul(t, r); + + 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; + + 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++) { + 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]; + } + } + } + + // Calculate and store AK coefficients: + + RealType eksq = 1.0; + RealType expf = 0.0; + if (ralph < 0.0) expf = exp(ralph*rcl*rcl); + for (i = 1; i <= kSqLim; i++) { + RealType rksq = float(i)*rcl*rcl; + eksq = expf*eksq; + AK[i] = eConverter * eksq/rksq; + } + + /* + * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz) + * the values of ll, mm and nn are selected so that the symmetry of + * reciprocal lattice is taken into account i.e. the following + * rules apply. + * + * ll ranges over the values 0 to kMax only. + * + * mm ranges over 0 to kMax when ll=0 and over + * -kMax to kMax otherwise. + * nn ranges over 1 to kMax when ll=mm=0 and over + * -kMax to kMax otherwise. + * + * Hence the result of the summation must be doubled at the end. + */ + + std::vector clm(nMax, 0.0); + std::vector slm(nMax, 0.0); + std::vector ckr(nMax, 0.0); + std::vector skr(nMax, 0.0); + std::vector ckc(nMax, 0.0); + std::vector cks(nMax, 0.0); + std::vector dkc(nMax, 0.0); + std::vector dks(nMax, 0.0); + std::vector qkc(nMax, 0.0); + 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; + rl = xcl * float(ll); + for (int mmm = mMin; mmm <= kLim2; mmm++) { + int mm = mmm - kLimit; + int m = abs(mm) + 1; + rm = ycl * float(mm); + // Set temporary products of exponential terms + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + i = atom->getLocalIndex(); + if(mm < 0) { + 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]=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; + rn = zcl * float(nn); + // Test on magnitude of k vector: + int kk=ll*ll + mm*mm + nn*nn; + if(kk <= kSqLim) { + 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)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + i = atom->getLocalIndex(); + + if (nn < 0) { + 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]*enc[n][i]-slm[i]*ens[n][i]; + skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i]; + } + } + } + + // Calculate scalar and vector products for each site: + + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + i = atom->getLocalIndex(); + int atid = atom->getAtomType()->getIdent(); + data = ElectrostaticMap[Etids[atid]]; + + if (data.is_Charge) { + C = data.fixedCharge; + if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); + ckc[i] = C * ckr[i]; + cks[i] = C * skr[i]; + } + + if (data.is_Dipole) { + 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) { + 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]; + } + } + } + + // calculate vector sums + + 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_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-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: + + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + i = atom->getLocalIndex(); + 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)); + + 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] ); + } + if (data.is_Quadrupole) { + atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] ); + } + } + } + } + } + nMin = 1; + } + mMin = 1; + } + pot += kPot; } }