| 40 |
|
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 42 |
|
|
| 43 |
+ |
#ifdef IS_MPI |
| 44 |
+ |
#include <mpi.h> |
| 45 |
+ |
#endif |
| 46 |
+ |
|
| 47 |
|
#include <stdio.h> |
| 48 |
|
#include <string.h> |
| 49 |
|
|
| 50 |
|
#include <cmath> |
| 51 |
+ |
#include <numeric> |
| 52 |
|
#include "nonbonded/Electrostatic.hpp" |
| 53 |
|
#include "utils/simError.h" |
| 54 |
|
#include "types/NonBondedInteractionType.hpp" |
| 60 |
|
#include "utils/PhysicalConstants.hpp" |
| 61 |
|
#include "math/erfc.hpp" |
| 62 |
|
#include "math/SquareMatrix.hpp" |
| 63 |
+ |
#include "primitives/Molecule.hpp" |
| 64 |
+ |
#include "flucq/FluctuatingChargeForces.hpp" |
| 65 |
|
|
| 66 |
|
namespace OpenMD { |
| 67 |
|
|
| 71 |
|
haveDampingAlpha_(false), |
| 72 |
|
haveDielectric_(false), |
| 73 |
|
haveElectroSplines_(false) |
| 74 |
< |
{} |
| 74 |
> |
{ |
| 75 |
> |
flucQ_ = new FluctuatingChargeForces(info_); |
| 76 |
> |
} |
| 77 |
|
|
| 78 |
+ |
void Electrostatic::setForceField(ForceField *ff) { |
| 79 |
+ |
forceField_ = ff; |
| 80 |
+ |
flucQ_->setForceField(forceField_); |
| 81 |
+ |
} |
| 82 |
+ |
|
| 83 |
+ |
void Electrostatic::setSimulatedAtomTypes(set<AtomType*> &simtypes) { |
| 84 |
+ |
simTypes_ = simtypes; |
| 85 |
+ |
flucQ_->setSimulatedAtomTypes(simTypes_); |
| 86 |
+ |
} |
| 87 |
+ |
|
| 88 |
|
void Electrostatic::initialize() { |
| 89 |
|
|
| 90 |
|
Globals* simParams_ = info_->getSimParams(); |
| 210 |
|
simError(); |
| 211 |
|
} |
| 212 |
|
|
| 213 |
< |
if (screeningMethod_ == DAMPED) { |
| 213 |
> |
if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) { |
| 214 |
|
if (!simParams_->haveDampingAlpha()) { |
| 215 |
|
// first set a cutoff dependent alpha value |
| 216 |
|
// we assume alpha depends linearly with rcut from 0 to 20.5 ang |
| 217 |
|
dampingAlpha_ = 0.425 - cutoffRadius_* 0.02; |
| 218 |
< |
if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0; |
| 200 |
< |
|
| 218 |
> |
if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0; |
| 219 |
|
// throw warning |
| 220 |
|
sprintf( painCave.errMsg, |
| 221 |
|
"Electrostatic::initialize: dampingAlpha was not specified in the\n" |
| 231 |
|
haveDampingAlpha_ = true; |
| 232 |
|
} |
| 233 |
|
|
| 234 |
+ |
|
| 235 |
|
Etypes.clear(); |
| 236 |
|
Etids.clear(); |
| 237 |
|
FQtypes.clear(); |
| 281 |
|
b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; |
| 282 |
|
b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; |
| 283 |
|
b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; |
| 265 |
– |
//selfMult1_ = - 2.0 * a2 * invArootPi; |
| 266 |
– |
//selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0; |
| 267 |
– |
//selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0; |
| 284 |
|
// Half the Smith self piece: |
| 285 |
|
selfMult1_ = - a2 * invArootPi; |
| 286 |
|
selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0; |
| 304 |
|
db0c_3 = 3.0*r*b2c - r2*r*b3c; |
| 305 |
|
db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; |
| 306 |
|
db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; |
| 291 |
– |
|
| 292 |
– |
selfMult1_ -= b0c; |
| 293 |
– |
selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; |
| 294 |
– |
selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; |
| 307 |
|
|
| 308 |
+ |
if (summationMethod_ != esm_EWALD_FULL) { |
| 309 |
+ |
selfMult1_ -= b0c; |
| 310 |
+ |
selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; |
| 311 |
+ |
selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; |
| 312 |
+ |
} |
| 313 |
+ |
|
| 314 |
|
// working variables for the splines: |
| 315 |
|
RealType ri, ri2; |
| 316 |
|
RealType b0, b1, b2, b3, b4, b5; |
| 346 |
|
vector<RealType> v31v, v32v; |
| 347 |
|
vector<RealType> v41v, v42v, v43v; |
| 348 |
|
|
| 331 |
– |
/* |
| 332 |
– |
vector<RealType> dv01v; |
| 333 |
– |
vector<RealType> dv11v; |
| 334 |
– |
vector<RealType> dv21v, dv22v; |
| 335 |
– |
vector<RealType> dv31v, dv32v; |
| 336 |
– |
vector<RealType> dv41v, dv42v, dv43v; |
| 337 |
– |
*/ |
| 338 |
– |
|
| 349 |
|
for (int i = 1; i < np_ + 1; i++) { |
| 350 |
|
r = RealType(i) * dx; |
| 351 |
|
rv.push_back(r); |
| 509 |
|
|
| 510 |
|
case esm_SWITCHING_FUNCTION: |
| 511 |
|
case esm_HARD: |
| 512 |
+ |
case esm_EWALD_FULL: |
| 513 |
|
|
| 514 |
|
v01 = f; |
| 515 |
|
v11 = g; |
| 568 |
|
|
| 569 |
|
break; |
| 570 |
|
|
| 560 |
– |
case esm_EWALD_FULL: |
| 571 |
|
case esm_EWALD_PME: |
| 572 |
|
case esm_EWALD_SPME: |
| 573 |
|
default : |
| 596 |
|
v41v.push_back(v41); |
| 597 |
|
v42v.push_back(v42); |
| 598 |
|
v43v.push_back(v43); |
| 589 |
– |
/* |
| 590 |
– |
dv01v.push_back(dv01); |
| 591 |
– |
dv11v.push_back(dv11); |
| 592 |
– |
dv21v.push_back(dv21); |
| 593 |
– |
dv22v.push_back(dv22); |
| 594 |
– |
dv31v.push_back(dv31); |
| 595 |
– |
dv32v.push_back(dv32); |
| 596 |
– |
dv41v.push_back(dv41); |
| 597 |
– |
dv42v.push_back(dv42); |
| 598 |
– |
dv43v.push_back(dv43); |
| 599 |
– |
*/ |
| 599 |
|
} |
| 600 |
|
|
| 601 |
|
// construct the spline structures and fill them with the values we've |
| 620 |
|
v43s = new CubicSpline(); |
| 621 |
|
v43s->addPoints(rv, v43v); |
| 622 |
|
|
| 624 |
– |
/* |
| 625 |
– |
dv01s = new CubicSpline(); |
| 626 |
– |
dv01s->addPoints(rv, dv01v); |
| 627 |
– |
dv11s = new CubicSpline(); |
| 628 |
– |
dv11s->addPoints(rv, dv11v); |
| 629 |
– |
dv21s = new CubicSpline(); |
| 630 |
– |
dv21s->addPoints(rv, dv21v); |
| 631 |
– |
dv22s = new CubicSpline(); |
| 632 |
– |
dv22s->addPoints(rv, dv22v); |
| 633 |
– |
dv31s = new CubicSpline(); |
| 634 |
– |
dv31s->addPoints(rv, dv31v); |
| 635 |
– |
dv32s = new CubicSpline(); |
| 636 |
– |
dv32s->addPoints(rv, dv32v); |
| 637 |
– |
dv41s = new CubicSpline(); |
| 638 |
– |
dv41s->addPoints(rv, dv41v); |
| 639 |
– |
dv42s = new CubicSpline(); |
| 640 |
– |
dv42s->addPoints(rv, dv42v); |
| 641 |
– |
dv43s = new CubicSpline(); |
| 642 |
– |
dv43s->addPoints(rv, dv43v); |
| 643 |
– |
*/ |
| 644 |
– |
|
| 623 |
|
haveElectroSplines_ = true; |
| 624 |
|
|
| 625 |
|
initialized_ = true; |
| 693 |
|
FQtids[atid] = fqtid; |
| 694 |
|
Jij[fqtid].resize(nFlucq_); |
| 695 |
|
|
| 696 |
< |
// Now, iterate over all known fluctuating and add to the coulomb integral map: |
| 696 |
> |
// Now, iterate over all known fluctuating and add to the |
| 697 |
> |
// coulomb integral map: |
| 698 |
|
|
| 699 |
|
std::set<int>::iterator it; |
| 700 |
|
for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { |
| 780 |
|
// Excluded potential that is still computed for fluctuating charges |
| 781 |
|
excluded_Pot= 0.0; |
| 782 |
|
|
| 804 |
– |
|
| 783 |
|
// some variables we'll need independent of electrostatic type: |
| 784 |
|
|
| 785 |
|
ri = 1.0 / *(idat.rij); |
| 900 |
|
Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or |
| 901 |
|
+ rdQbr * rhat * (dv22 - 2.0*v22or)); |
| 902 |
|
} |
| 903 |
< |
|
| 903 |
> |
|
| 904 |
> |
|
| 905 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 906 |
|
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 907 |
|
} |
| 908 |
< |
|
| 908 |
> |
|
| 909 |
|
if (a_is_Charge) { |
| 910 |
|
|
| 911 |
|
if (b_is_Charge) { |
| 912 |
|
pref = pre11_ * *(idat.electroMult); |
| 913 |
|
U += C_a * C_b * pref * v01; |
| 914 |
|
F += C_a * C_b * pref * dv01 * rhat; |
| 915 |
< |
|
| 915 |
> |
|
| 916 |
|
// If this is an excluded pair, there are still indirect |
| 917 |
|
// interactions via the reaction field we must worry about: |
| 918 |
|
|
| 921 |
|
indirect_Pot += rfContrib; |
| 922 |
|
indirect_F += rfContrib * 2.0 * ri * rhat; |
| 923 |
|
} |
| 924 |
< |
|
| 924 |
> |
|
| 925 |
|
// Fluctuating charge forces are handled via Coulomb integrals |
| 926 |
|
// for excluded pairs (i.e. those connected via bonds) and |
| 927 |
|
// with the standard charge-charge interaction otherwise. |
| 928 |
|
|
| 929 |
< |
if (idat.excluded) { |
| 929 |
> |
if (idat.excluded) { |
| 930 |
|
if (a_is_Fluctuating || b_is_Fluctuating) { |
| 931 |
|
coulInt = J->getValueAt( *(idat.rij) ); |
| 932 |
< |
if (a_is_Fluctuating) dUdCa += coulInt * C_b; |
| 933 |
< |
if (b_is_Fluctuating) dUdCb += coulInt * C_a; |
| 934 |
< |
excluded_Pot += C_a * C_b * coulInt; |
| 956 |
< |
} |
| 932 |
> |
if (a_is_Fluctuating) dUdCa += C_b * coulInt; |
| 933 |
> |
if (b_is_Fluctuating) dUdCb += C_a * coulInt; |
| 934 |
> |
} |
| 935 |
|
} else { |
| 936 |
|
if (a_is_Fluctuating) dUdCa += C_b * pref * v01; |
| 937 |
< |
if (a_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 938 |
< |
} |
| 937 |
> |
if (b_is_Fluctuating) dUdCb += C_a * pref * v01; |
| 938 |
> |
} |
| 939 |
|
} |
| 940 |
|
|
| 941 |
|
if (b_is_Dipole) { |
| 1001 |
|
F -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat; |
| 1002 |
|
Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); |
| 1003 |
|
Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb); |
| 1026 |
– |
|
| 1004 |
|
// Even if we excluded this pair from direct interactions, we |
| 1005 |
|
// still have the reaction-field-mediated dipole-dipole |
| 1006 |
|
// interaction: |
| 1060 |
|
trQaQb = QaQb.trace(); |
| 1061 |
|
rQaQb = rhat * QaQb; |
| 1062 |
|
QaQbr = QaQb * rhat; |
| 1063 |
< |
QaxQb = cross(Q_a, Q_b); |
| 1063 |
> |
QaxQb = mCross(Q_a, Q_b); |
| 1064 |
|
rQaQbr = dot(rQa, Qbr); |
| 1065 |
|
rQaxQbr = cross(rQa, Qbr); |
| 1066 |
|
|
| 1091 |
|
// + 4.0 * cross(rhat, QbQar) |
| 1092 |
|
|
| 1093 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1117 |
– |
|
| 1094 |
|
} |
| 1095 |
|
} |
| 1096 |
|
|
| 1153 |
|
|
| 1154 |
|
if (i_is_Fluctuating) { |
| 1155 |
|
C_a += *(sdat.flucQ); |
| 1156 |
< |
// dVdFQ is really a force, so this is negative the derivative |
| 1157 |
< |
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1158 |
< |
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * |
| 1159 |
< |
(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); |
| 1156 |
> |
|
| 1157 |
> |
flucQ_->getSelfInteraction(sdat.atid, *(sdat.flucQ), |
| 1158 |
> |
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY], |
| 1159 |
> |
*(sdat.flucQfrc) ); |
| 1160 |
> |
|
| 1161 |
|
} |
| 1162 |
|
|
| 1163 |
|
switch (summationMethod_) { |
| 1179 |
|
|
| 1180 |
|
case esm_SHIFTED_FORCE: |
| 1181 |
|
case esm_SHIFTED_POTENTIAL: |
| 1182 |
+ |
case esm_TAYLOR_SHIFTED: |
| 1183 |
+ |
case esm_EWALD_FULL: |
| 1184 |
|
if (i_is_Charge) |
| 1185 |
|
self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); |
| 1186 |
|
if (i_is_Dipole) |
| 1205 |
|
// 12 angstroms seems to be a reasonably good guess for most |
| 1206 |
|
// cases. |
| 1207 |
|
return 12.0; |
| 1208 |
+ |
} |
| 1209 |
+ |
|
| 1210 |
+ |
|
| 1211 |
+ |
void Electrostatic::ReciprocalSpaceSum(RealType& pot) { |
| 1212 |
+ |
|
| 1213 |
+ |
RealType kPot = 0.0; |
| 1214 |
+ |
RealType kVir = 0.0; |
| 1215 |
+ |
|
| 1216 |
+ |
const RealType mPoleConverter = 0.20819434; // converts from the |
| 1217 |
+ |
// internal units of |
| 1218 |
+ |
// Debye (for dipoles) |
| 1219 |
+ |
// or Debye-angstroms |
| 1220 |
+ |
// (for quadrupoles) to |
| 1221 |
+ |
// electron angstroms or |
| 1222 |
+ |
// electron-angstroms^2 |
| 1223 |
+ |
|
| 1224 |
+ |
const RealType eConverter = 332.0637778; // convert the |
| 1225 |
+ |
// Charge-Charge |
| 1226 |
+ |
// electrostatic |
| 1227 |
+ |
// interactions into kcal / |
| 1228 |
+ |
// mol assuming distances |
| 1229 |
+ |
// are measured in |
| 1230 |
+ |
// angstroms. |
| 1231 |
+ |
|
| 1232 |
+ |
Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat(); |
| 1233 |
+ |
Vector3d box = hmat.diagonals(); |
| 1234 |
+ |
RealType boxMax = box.max(); |
| 1235 |
+ |
|
| 1236 |
+ |
//int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); |
| 1237 |
+ |
int kMax = 7; |
| 1238 |
+ |
int kSqMax = kMax*kMax + 2; |
| 1239 |
+ |
|
| 1240 |
+ |
int kLimit = kMax+1; |
| 1241 |
+ |
int kLim2 = 2*kMax+1; |
| 1242 |
+ |
int kSqLim = kSqMax; |
| 1243 |
+ |
|
| 1244 |
+ |
vector<RealType> AK(kSqLim+1, 0.0); |
| 1245 |
+ |
RealType xcl = 2.0 * M_PI / box.x(); |
| 1246 |
+ |
RealType ycl = 2.0 * M_PI / box.y(); |
| 1247 |
+ |
RealType zcl = 2.0 * M_PI / box.z(); |
| 1248 |
+ |
RealType rcl = 2.0 * M_PI / boxMax; |
| 1249 |
+ |
RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z()); |
| 1250 |
+ |
|
| 1251 |
+ |
if(dampingAlpha_ < 1.0e-12) return; |
| 1252 |
+ |
|
| 1253 |
+ |
RealType ralph = -0.25/pow(dampingAlpha_,2); |
| 1254 |
+ |
|
| 1255 |
+ |
// Calculate and store exponential factors |
| 1256 |
+ |
|
| 1257 |
+ |
vector<vector<RealType> > elc; |
| 1258 |
+ |
vector<vector<RealType> > emc; |
| 1259 |
+ |
vector<vector<RealType> > enc; |
| 1260 |
+ |
vector<vector<RealType> > els; |
| 1261 |
+ |
vector<vector<RealType> > ems; |
| 1262 |
+ |
vector<vector<RealType> > ens; |
| 1263 |
+ |
|
| 1264 |
+ |
int nMax = info_->getNAtoms(); |
| 1265 |
+ |
|
| 1266 |
+ |
elc.resize(kLimit+1); |
| 1267 |
+ |
emc.resize(kLimit+1); |
| 1268 |
+ |
enc.resize(kLimit+1); |
| 1269 |
+ |
els.resize(kLimit+1); |
| 1270 |
+ |
ems.resize(kLimit+1); |
| 1271 |
+ |
ens.resize(kLimit+1); |
| 1272 |
+ |
|
| 1273 |
+ |
for (int j = 0; j < kLimit+1; j++) { |
| 1274 |
+ |
elc[j].resize(nMax); |
| 1275 |
+ |
emc[j].resize(nMax); |
| 1276 |
+ |
enc[j].resize(nMax); |
| 1277 |
+ |
els[j].resize(nMax); |
| 1278 |
+ |
ems[j].resize(nMax); |
| 1279 |
+ |
ens[j].resize(nMax); |
| 1280 |
+ |
} |
| 1281 |
+ |
|
| 1282 |
+ |
Vector3d t( 2.0 * M_PI ); |
| 1283 |
+ |
t.Vdiv(t, box); |
| 1284 |
+ |
|
| 1285 |
+ |
SimInfo::MoleculeIterator mi; |
| 1286 |
+ |
Molecule::AtomIterator ai; |
| 1287 |
+ |
int i; |
| 1288 |
+ |
Vector3d r; |
| 1289 |
+ |
Vector3d tt; |
| 1290 |
+ |
|
| 1291 |
+ |
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1292 |
+ |
mol = info_->nextMolecule(mi)) { |
| 1293 |
+ |
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1294 |
+ |
atom = mol->nextAtom(ai)) { |
| 1295 |
+ |
|
| 1296 |
+ |
i = atom->getLocalIndex(); |
| 1297 |
+ |
r = atom->getPos(); |
| 1298 |
+ |
info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); |
| 1299 |
+ |
|
| 1300 |
+ |
tt.Vmul(t, r); |
| 1301 |
+ |
|
| 1302 |
+ |
elc[1][i] = 1.0; |
| 1303 |
+ |
emc[1][i] = 1.0; |
| 1304 |
+ |
enc[1][i] = 1.0; |
| 1305 |
+ |
els[1][i] = 0.0; |
| 1306 |
+ |
ems[1][i] = 0.0; |
| 1307 |
+ |
ens[1][i] = 0.0; |
| 1308 |
+ |
|
| 1309 |
+ |
elc[2][i] = cos(tt.x()); |
| 1310 |
+ |
emc[2][i] = cos(tt.y()); |
| 1311 |
+ |
enc[2][i] = cos(tt.z()); |
| 1312 |
+ |
els[2][i] = sin(tt.x()); |
| 1313 |
+ |
ems[2][i] = sin(tt.y()); |
| 1314 |
+ |
ens[2][i] = sin(tt.z()); |
| 1315 |
+ |
|
| 1316 |
+ |
for(int l = 3; l <= kLimit; l++) { |
| 1317 |
+ |
elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i]; |
| 1318 |
+ |
emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i]; |
| 1319 |
+ |
enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i]; |
| 1320 |
+ |
els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i]; |
| 1321 |
+ |
ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i]; |
| 1322 |
+ |
ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i]; |
| 1323 |
+ |
} |
| 1324 |
+ |
} |
| 1325 |
+ |
} |
| 1326 |
+ |
|
| 1327 |
+ |
// Calculate and store AK coefficients: |
| 1328 |
+ |
|
| 1329 |
+ |
RealType eksq = 1.0; |
| 1330 |
+ |
RealType expf = 0.0; |
| 1331 |
+ |
if (ralph < 0.0) expf = exp(ralph*rcl*rcl); |
| 1332 |
+ |
for (i = 1; i <= kSqLim; i++) { |
| 1333 |
+ |
RealType rksq = float(i)*rcl*rcl; |
| 1334 |
+ |
eksq = expf*eksq; |
| 1335 |
+ |
AK[i] = eConverter * eksq/rksq; |
| 1336 |
+ |
} |
| 1337 |
+ |
|
| 1338 |
+ |
/* |
| 1339 |
+ |
* Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz) |
| 1340 |
+ |
* the values of ll, mm and nn are selected so that the symmetry of |
| 1341 |
+ |
* reciprocal lattice is taken into account i.e. the following |
| 1342 |
+ |
* rules apply. |
| 1343 |
+ |
* |
| 1344 |
+ |
* ll ranges over the values 0 to kMax only. |
| 1345 |
+ |
* |
| 1346 |
+ |
* mm ranges over 0 to kMax when ll=0 and over |
| 1347 |
+ |
* -kMax to kMax otherwise. |
| 1348 |
+ |
* nn ranges over 1 to kMax when ll=mm=0 and over |
| 1349 |
+ |
* -kMax to kMax otherwise. |
| 1350 |
+ |
* |
| 1351 |
+ |
* Hence the result of the summation must be doubled at the end. |
| 1352 |
+ |
*/ |
| 1353 |
+ |
|
| 1354 |
+ |
std::vector<RealType> clm(nMax, 0.0); |
| 1355 |
+ |
std::vector<RealType> slm(nMax, 0.0); |
| 1356 |
+ |
std::vector<RealType> ckr(nMax, 0.0); |
| 1357 |
+ |
std::vector<RealType> skr(nMax, 0.0); |
| 1358 |
+ |
std::vector<RealType> ckc(nMax, 0.0); |
| 1359 |
+ |
std::vector<RealType> cks(nMax, 0.0); |
| 1360 |
+ |
std::vector<RealType> dkc(nMax, 0.0); |
| 1361 |
+ |
std::vector<RealType> dks(nMax, 0.0); |
| 1362 |
+ |
std::vector<RealType> qkc(nMax, 0.0); |
| 1363 |
+ |
std::vector<RealType> qks(nMax, 0.0); |
| 1364 |
+ |
std::vector<Vector3d> dxk(nMax, V3Zero); |
| 1365 |
+ |
std::vector<Vector3d> qxk(nMax, V3Zero); |
| 1366 |
+ |
RealType rl, rm, rn; |
| 1367 |
+ |
Vector3d kVec; |
| 1368 |
+ |
Vector3d Qk; |
| 1369 |
+ |
Mat3x3d k2; |
| 1370 |
+ |
RealType ckcs, ckss, dkcs, dkss, qkcs, qkss; |
| 1371 |
+ |
int atid; |
| 1372 |
+ |
ElectrostaticAtomData data; |
| 1373 |
+ |
RealType C, dk, qk; |
| 1374 |
+ |
Vector3d D; |
| 1375 |
+ |
Mat3x3d Q; |
| 1376 |
+ |
|
| 1377 |
+ |
int mMin = kLimit; |
| 1378 |
+ |
int nMin = kLimit + 1; |
| 1379 |
+ |
for (int l = 1; l <= kLimit; l++) { |
| 1380 |
+ |
int ll = l - 1; |
| 1381 |
+ |
rl = xcl * float(ll); |
| 1382 |
+ |
for (int mmm = mMin; mmm <= kLim2; mmm++) { |
| 1383 |
+ |
int mm = mmm - kLimit; |
| 1384 |
+ |
int m = abs(mm) + 1; |
| 1385 |
+ |
rm = ycl * float(mm); |
| 1386 |
+ |
// Set temporary products of exponential terms |
| 1387 |
+ |
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1388 |
+ |
mol = info_->nextMolecule(mi)) { |
| 1389 |
+ |
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1390 |
+ |
atom = mol->nextAtom(ai)) { |
| 1391 |
+ |
|
| 1392 |
+ |
i = atom->getLocalIndex(); |
| 1393 |
+ |
if(mm < 0) { |
| 1394 |
+ |
clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i]; |
| 1395 |
+ |
slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i]; |
| 1396 |
+ |
} else { |
| 1397 |
+ |
clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i]; |
| 1398 |
+ |
slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i]; |
| 1399 |
+ |
} |
| 1400 |
+ |
} |
| 1401 |
+ |
} |
| 1402 |
+ |
for (int nnn = nMin; nnn <= kLim2; nnn++) { |
| 1403 |
+ |
int nn = nnn - kLimit; |
| 1404 |
+ |
int n = abs(nn) + 1; |
| 1405 |
+ |
rn = zcl * float(nn); |
| 1406 |
+ |
// Test on magnitude of k vector: |
| 1407 |
+ |
int kk=ll*ll + mm*mm + nn*nn; |
| 1408 |
+ |
if(kk <= kSqLim) { |
| 1409 |
+ |
kVec = Vector3d(rl, rm, rn); |
| 1410 |
+ |
k2 = outProduct(kVec, kVec); |
| 1411 |
+ |
// Calculate exp(ikr) terms |
| 1412 |
+ |
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1413 |
+ |
mol = info_->nextMolecule(mi)) { |
| 1414 |
+ |
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1415 |
+ |
atom = mol->nextAtom(ai)) { |
| 1416 |
+ |
i = atom->getLocalIndex(); |
| 1417 |
+ |
|
| 1418 |
+ |
if (nn < 0) { |
| 1419 |
+ |
ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i]; |
| 1420 |
+ |
skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i]; |
| 1421 |
+ |
|
| 1422 |
+ |
} else { |
| 1423 |
+ |
ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i]; |
| 1424 |
+ |
skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i]; |
| 1425 |
+ |
} |
| 1426 |
+ |
} |
| 1427 |
+ |
} |
| 1428 |
+ |
|
| 1429 |
+ |
// Calculate scalar and vector products for each site: |
| 1430 |
+ |
|
| 1431 |
+ |
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1432 |
+ |
mol = info_->nextMolecule(mi)) { |
| 1433 |
+ |
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1434 |
+ |
atom = mol->nextAtom(ai)) { |
| 1435 |
+ |
i = atom->getLocalIndex(); |
| 1436 |
+ |
int atid = atom->getAtomType()->getIdent(); |
| 1437 |
+ |
data = ElectrostaticMap[Etids[atid]]; |
| 1438 |
+ |
|
| 1439 |
+ |
if (data.is_Charge) { |
| 1440 |
+ |
C = data.fixedCharge; |
| 1441 |
+ |
if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); |
| 1442 |
+ |
ckc[i] = C * ckr[i]; |
| 1443 |
+ |
cks[i] = C * skr[i]; |
| 1444 |
+ |
} |
| 1445 |
+ |
|
| 1446 |
+ |
if (data.is_Dipole) { |
| 1447 |
+ |
D = atom->getDipole() * mPoleConverter; |
| 1448 |
+ |
dk = dot(D, kVec); |
| 1449 |
+ |
dxk[i] = cross(D, kVec); |
| 1450 |
+ |
dkc[i] = dk * ckr[i]; |
| 1451 |
+ |
dks[i] = dk * skr[i]; |
| 1452 |
+ |
} |
| 1453 |
+ |
if (data.is_Quadrupole) { |
| 1454 |
+ |
Q = atom->getQuadrupole() * mPoleConverter; |
| 1455 |
+ |
Qk = Q * kVec; |
| 1456 |
+ |
qk = dot(kVec, Qk); |
| 1457 |
+ |
qxk[i] = -cross(kVec, Qk); |
| 1458 |
+ |
qkc[i] = qk * ckr[i]; |
| 1459 |
+ |
qks[i] = qk * skr[i]; |
| 1460 |
+ |
} |
| 1461 |
+ |
} |
| 1462 |
+ |
} |
| 1463 |
+ |
|
| 1464 |
+ |
// calculate vector sums |
| 1465 |
+ |
|
| 1466 |
+ |
ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); |
| 1467 |
+ |
ckss = std::accumulate(cks.begin(),cks.end(),0.0); |
| 1468 |
+ |
dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0); |
| 1469 |
+ |
dkss = std::accumulate(dks.begin(),dks.end(),0.0); |
| 1470 |
+ |
qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); |
| 1471 |
+ |
qkss = std::accumulate(qks.begin(),qks.end(),0.0); |
| 1472 |
+ |
|
| 1473 |
+ |
#ifdef IS_MPI |
| 1474 |
+ |
MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE, |
| 1475 |
+ |
MPI_SUM, MPI_COMM_WORLD); |
| 1476 |
+ |
MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE, |
| 1477 |
+ |
MPI_SUM, MPI_COMM_WORLD); |
| 1478 |
+ |
MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE, |
| 1479 |
+ |
MPI_SUM, MPI_COMM_WORLD); |
| 1480 |
+ |
MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE, |
| 1481 |
+ |
MPI_SUM, MPI_COMM_WORLD); |
| 1482 |
+ |
MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE, |
| 1483 |
+ |
MPI_SUM, MPI_COMM_WORLD); |
| 1484 |
+ |
MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE, |
| 1485 |
+ |
MPI_SUM, MPI_COMM_WORLD); |
| 1486 |
+ |
#endif |
| 1487 |
+ |
|
| 1488 |
+ |
// Accumulate potential energy and virial contribution: |
| 1489 |
+ |
|
| 1490 |
+ |
kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) |
| 1491 |
+ |
+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs)); |
| 1492 |
+ |
|
| 1493 |
+ |
kVir += 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss |
| 1494 |
+ |
+4.0*(ckss*dkcs-ckcs*dkss) |
| 1495 |
+ |
+3.0*(dkcs*dkcs+dkss*dkss) |
| 1496 |
+ |
-6.0*(ckss*qkss+ckcs*qkcs) |
| 1497 |
+ |
+8.0*(dkss*qkcs-dkcs*qkss) |
| 1498 |
+ |
+5.0*(qkss*qkss+qkcs*qkcs)); |
| 1499 |
+ |
|
| 1500 |
+ |
// Calculate force and torque for each site: |
| 1501 |
+ |
|
| 1502 |
+ |
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1503 |
+ |
mol = info_->nextMolecule(mi)) { |
| 1504 |
+ |
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1505 |
+ |
atom = mol->nextAtom(ai)) { |
| 1506 |
+ |
|
| 1507 |
+ |
i = atom->getLocalIndex(); |
| 1508 |
+ |
atid = atom->getAtomType()->getIdent(); |
| 1509 |
+ |
data = ElectrostaticMap[Etids[atid]]; |
| 1510 |
+ |
|
| 1511 |
+ |
RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs) |
| 1512 |
+ |
- (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); |
| 1513 |
+ |
RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) |
| 1514 |
+ |
-ckr[i]*(ckss+dkcs-qkss)); |
| 1515 |
+ |
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) |
| 1516 |
+ |
+skr[i]*(ckss+dkcs-qkss)); |
| 1517 |
+ |
|
| 1518 |
+ |
atom->addFrc( 4.0 * rvol * qfrc * kVec ); |
| 1519 |
+ |
|
| 1520 |
+ |
if (atom->isFluctuatingCharge()) { |
| 1521 |
+ |
atom->addFlucQFrc( - 2.0 * rvol * qtrq2 ); |
| 1522 |
+ |
} |
| 1523 |
+ |
|
| 1524 |
+ |
if (data.is_Dipole) { |
| 1525 |
+ |
atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] ); |
| 1526 |
+ |
} |
| 1527 |
+ |
if (data.is_Quadrupole) { |
| 1528 |
+ |
atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] ); |
| 1529 |
+ |
} |
| 1530 |
+ |
} |
| 1531 |
+ |
} |
| 1532 |
+ |
} |
| 1533 |
+ |
} |
| 1534 |
+ |
nMin = 1; |
| 1535 |
+ |
} |
| 1536 |
+ |
mMin = 1; |
| 1537 |
+ |
} |
| 1538 |
+ |
pot += kPot; |
| 1539 |
|
} |
| 1540 |
|
} |
