| 289 |
|
db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; |
| 290 |
|
db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; |
| 291 |
|
|
| 292 |
< |
if (summationMethod_ != esm_EWALD_FULL) { |
| 292 |
> |
if (summationMethod_ == esm_EWALD_FULL) { |
| 293 |
> |
selfMult1_ *= 2.0; |
| 294 |
> |
selfMult2_ *= 2.0; |
| 295 |
> |
selfMult4_ *= 2.0; |
| 296 |
> |
} else { |
| 297 |
|
selfMult1_ -= b0c; |
| 298 |
|
selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; |
| 299 |
|
selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; |
| 681 |
|
FQtids[atid] = fqtid; |
| 682 |
|
Jij[fqtid].resize(nFlucq_); |
| 683 |
|
|
| 684 |
< |
// Now, iterate over all known fluctuating and add to the coulomb integral map: |
| 684 |
> |
// Now, iterate over all known fluctuating and add to the |
| 685 |
> |
// coulomb integral map: |
| 686 |
|
|
| 687 |
|
std::set<int>::iterator it; |
| 688 |
|
for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { |
| 1170 |
|
case esm_SHIFTED_FORCE: |
| 1171 |
|
case esm_SHIFTED_POTENTIAL: |
| 1172 |
|
case esm_TAYLOR_SHIFTED: |
| 1173 |
+ |
case esm_EWALD_FULL: |
| 1174 |
|
if (i_is_Charge) |
| 1175 |
|
self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); |
| 1176 |
|
if (i_is_Dipole) |
| 1198 |
|
} |
| 1199 |
|
|
| 1200 |
|
|
| 1201 |
< |
void Electrostatic::ReciprocalSpaceSum () { |
| 1201 |
> |
void Electrostatic::ReciprocalSpaceSum(potVec& pot) { |
| 1202 |
|
|
| 1203 |
|
RealType kPot = 0.0; |
| 1204 |
|
RealType kVir = 0.0; |
| 1223 |
|
Vector3d box = hmat.diagonals(); |
| 1224 |
|
RealType boxMax = box.max(); |
| 1225 |
|
|
| 1226 |
< |
//int kMax = int(pow(dampingAlpha_,2)*cutoffRadius_ * boxMax / M_PI); |
| 1227 |
< |
const int kMax = 5; |
| 1226 |
> |
cerr << "da = " << dampingAlpha_ << " rc = " << cutoffRadius_ << "\n"; |
| 1227 |
> |
cerr << "boxMax = " << boxMax << "\n"; |
| 1228 |
> |
//int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); |
| 1229 |
> |
int kMax = 5; |
| 1230 |
> |
cerr << "kMax = " << kMax << "\n"; |
| 1231 |
|
int kSqMax = kMax*kMax + 2; |
| 1232 |
|
|
| 1233 |
|
int kLimit = kMax+1; |
| 1269 |
|
Vector3d tt; |
| 1270 |
|
Vector3d w; |
| 1271 |
|
Vector3d u; |
| 1272 |
+ |
Vector3d a; |
| 1273 |
+ |
Vector3d b; |
| 1274 |
|
|
| 1275 |
|
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1276 |
|
mol = info_->nextMolecule(mi)) { |
| 1281 |
|
r = atom->getPos(); |
| 1282 |
|
info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); |
| 1283 |
|
|
| 1284 |
+ |
// Shift so that all coordinates are in the range [0,L]: |
| 1285 |
+ |
|
| 1286 |
+ |
r += box/2.0; |
| 1287 |
+ |
|
| 1288 |
|
tt.Vmul(t, r); |
| 1289 |
+ |
|
| 1290 |
+ |
//cerr << "tt = " << tt << "\n"; |
| 1291 |
|
|
| 1292 |
|
eCos[1][i] = Vector3d(1.0, 1.0, 1.0); |
| 1293 |
|
eSin[1][i] = Vector3d(0.0, 0.0, 0.0); |
| 1294 |
|
eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z())); |
| 1295 |
|
eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z())); |
| 1296 |
< |
u = 2.0 * eCos[1][i]; |
| 1297 |
< |
eCos[3][i].Vmul(u, eCos[2][i]); |
| 1281 |
< |
eSin[3][i].Vmul(u, eSin[2][i]); |
| 1296 |
> |
u = eCos[2][i]; |
| 1297 |
> |
w = eSin[2][i]; |
| 1298 |
|
|
| 1299 |
|
for(int l = 3; l <= kLimit; l++) { |
| 1300 |
< |
w.Vmul(u, eCos[l-1][i]); |
| 1301 |
< |
eCos[l][i] = w - eCos[l-2][i]; |
| 1302 |
< |
w.Vmul(u, eSin[l-1][i]); |
| 1303 |
< |
eSin[l][i] = w - eSin[l-2][i]; |
| 1300 |
> |
a.Vmul(eCos[l-1][i], u); |
| 1301 |
> |
b.Vmul(eSin[l-1][i], w); |
| 1302 |
> |
eCos[l][i] = a - b; |
| 1303 |
> |
a.Vmul(eSin[l-1][i], u); |
| 1304 |
> |
b.Vmul(eCos[l-1][i], w); |
| 1305 |
> |
eSin[l][i] = a + b; |
| 1306 |
|
} |
| 1307 |
|
} |
| 1308 |
|
} |
| 1366 |
|
if(mm < 0) { |
| 1367 |
|
clm[i] = eCos[l][i].x()*eCos[m][i].y() |
| 1368 |
|
+ eSin[l][i].x()*eSin[m][i].y(); |
| 1369 |
< |
slm[i] = eCos[l][i].x()*eCos[m][i].y() |
| 1369 |
> |
slm[i] = eSin[l][i].x()*eCos[m][i].y() |
| 1370 |
|
- eSin[m][i].y()*eCos[l][i].x(); |
| 1371 |
|
} else { |
| 1372 |
|
clm[i] = eCos[l][i].x()*eCos[m][i].y() |
| 1408 |
|
mol = info_->nextMolecule(mi)) { |
| 1409 |
|
for(Atom* atom = mol->beginAtom(ai); atom != NULL; |
| 1410 |
|
atom = mol->nextAtom(ai)) { |
| 1411 |
< |
i = atom->getGlobalIndex(); |
| 1411 |
> |
i = atom->getLocalIndex(); |
| 1412 |
|
int atid = atom->getAtomType()->getIdent(); |
| 1413 |
|
ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; |
| 1414 |
|
|
| 1416 |
|
RealType C = data.fixedCharge; |
| 1417 |
|
if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); |
| 1418 |
|
ckc[i] = C * ckr[i]; |
| 1419 |
< |
cks[i] = C * cks[i]; |
| 1419 |
> |
cks[i] = C * skr[i]; |
| 1420 |
|
} |
| 1421 |
|
|
| 1422 |
|
if (data.is_Dipole) { |
| 1436 |
|
} |
| 1437 |
|
} |
| 1438 |
|
} |
| 1439 |
< |
|
| 1439 |
> |
|
| 1440 |
|
// calculate vector sums |
| 1441 |
|
|
| 1442 |
|
RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); |
| 1445 |
|
RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0); |
| 1446 |
|
RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); |
| 1447 |
|
RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0); |
| 1448 |
+ |
|
| 1449 |
|
|
| 1450 |
|
#ifdef IS_MPI |
| 1451 |
|
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, |
| 1462 |
|
MPI::SUM); |
| 1463 |
|
#endif |
| 1464 |
|
|
| 1465 |
< |
// Accumulate potential energy and virial contribution: |
| 1466 |
< |
|
| 1448 |
< |
//cerr << "l, m, n = " << l << " " << m << " " << n << "\n"; |
| 1449 |
< |
cerr << "kVec = " << kVec << "\n"; |
| 1450 |
< |
cerr << "ckss = " << ckss << " ckcs = " << ckcs << "\n"; |
| 1465 |
> |
// Accumulate potential energy and virial contribution: |
| 1466 |
> |
|
| 1467 |
|
kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) |
| 1468 |
|
+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkss)); |
| 1469 |
< |
//cerr << "kspace pot = " << kPot << "\n"; |
| 1469 |
> |
|
| 1470 |
|
kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss |
| 1471 |
|
+4.0*(ckss*dkcs-ckcs*dkss) |
| 1472 |
|
+3.0*(dkcs*dkcs+dkss*dkss) |
| 1507 |
|
} |
| 1508 |
|
} |
| 1509 |
|
} |
| 1510 |
+ |
cerr << "kPot = " << kPot << "\n"; |
| 1511 |
+ |
pot[ELECTROSTATIC_FAMILY] += kPot; |
| 1512 |
|
} |
| 1513 |
|
} |
