--- trunk/src/nonbonded/Electrostatic.cpp 2013/07/12 17:38:06 1900 +++ trunk/src/nonbonded/Electrostatic.cpp 2013/08/07 15:24:16 1925 @@ -44,6 +44,7 @@ #include #include +#include #include "nonbonded/Electrostatic.hpp" #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" @@ -55,6 +56,10 @@ #include "utils/PhysicalConstants.hpp" #include "math/erfc.hpp" #include "math/SquareMatrix.hpp" +#include "primitives/Molecule.hpp" +#ifdef IS_MPI +#include +#endif namespace OpenMD { @@ -191,13 +196,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 +217,7 @@ namespace OpenMD { haveDampingAlpha_ = true; } + Etypes.clear(); Etids.clear(); FQtypes.clear(); @@ -262,9 +267,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 +290,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; @@ -327,14 +331,6 @@ namespace OpenMD { vector v21v, v22v; 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; @@ -499,6 +495,7 @@ namespace OpenMD { case esm_SWITCHING_FUNCTION: case esm_HARD: + case esm_EWALD_FULL: v01 = f; v11 = g; @@ -557,7 +554,6 @@ namespace OpenMD { break; - case esm_EWALD_FULL: case esm_EWALD_PME: case esm_EWALD_SPME: default : @@ -586,17 +582,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 @@ -620,27 +605,6 @@ namespace OpenMD { v42s->addPoints(rv, v42v); 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; @@ -715,7 +679,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) { @@ -1202,6 +1167,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) @@ -1227,4 +1194,334 @@ namespace OpenMD { // 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(); + Q *= 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::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); +#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 (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; + } }