| 677 |
|
FQtids[atid] = fqtid; |
| 678 |
|
Jij[fqtid].resize(nFlucq_); |
| 679 |
|
|
| 680 |
< |
// Now, iterate over all known fluctuating and add to the coulomb integral map: |
| 680 |
> |
// Now, iterate over all known fluctuating and add to the |
| 681 |
> |
// coulomb integral map: |
| 682 |
|
|
| 683 |
|
std::set<int>::iterator it; |
| 684 |
|
for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { |
| 1166 |
|
case esm_SHIFTED_FORCE: |
| 1167 |
|
case esm_SHIFTED_POTENTIAL: |
| 1168 |
|
case esm_TAYLOR_SHIFTED: |
| 1169 |
+ |
case esm_EWALD_FULL: |
| 1170 |
|
if (i_is_Charge) |
| 1171 |
|
self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); |
| 1172 |
|
if (i_is_Dipole) |
| 1194 |
|
} |
| 1195 |
|
|
| 1196 |
|
|
| 1197 |
< |
void Electrostatic::ReciprocalSpaceSum () { |
| 1197 |
> |
void Electrostatic::ReciprocalSpaceSum(potVec& pot) { |
| 1198 |
|
|
| 1199 |
|
RealType kPot = 0.0; |
| 1200 |
|
RealType kVir = 0.0; |
| 1219 |
|
Vector3d box = hmat.diagonals(); |
| 1220 |
|
RealType boxMax = box.max(); |
| 1221 |
|
|
| 1222 |
< |
//int kMax = int(pow(dampingAlpha_,2)*cutoffRadius_ * boxMax / M_PI); |
| 1223 |
< |
const int kMax = 5; |
| 1222 |
> |
//int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); |
| 1223 |
> |
int kMax = 7; |
| 1224 |
|
int kSqMax = kMax*kMax + 2; |
| 1225 |
|
|
| 1226 |
|
int kLimit = kMax+1; |
| 1254 |
|
|
| 1255 |
|
Vector3d t( 2.0 * M_PI ); |
| 1256 |
|
t.Vdiv(t, box); |
| 1257 |
+ |
|
| 1258 |
|
|
| 1259 |
|
SimInfo::MoleculeIterator mi; |
| 1260 |
|
Molecule::AtomIterator ai; |
| 1263 |
|
Vector3d tt; |
| 1264 |
|
Vector3d w; |
| 1265 |
|
Vector3d u; |
| 1266 |
+ |
Vector3d a; |
| 1267 |
+ |
Vector3d b; |
| 1268 |
|
|
| 1269 |
|
for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; |
| 1270 |
|
mol = info_->nextMolecule(mi)) { |
| 1276 |
|
info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); |
| 1277 |
|
|
| 1278 |
|
tt.Vmul(t, r); |
| 1279 |
+ |
|
| 1280 |
|
|
| 1281 |
|
eCos[1][i] = Vector3d(1.0, 1.0, 1.0); |
| 1282 |
|
eSin[1][i] = Vector3d(0.0, 0.0, 0.0); |
| 1283 |
|
eCos[2][i] = Vector3d(cos(tt.x()), cos(tt.y()), cos(tt.z())); |
| 1284 |
|
eSin[2][i] = Vector3d(sin(tt.x()), sin(tt.y()), sin(tt.z())); |
| 1285 |
< |
u = 2.0 * eCos[1][i]; |
| 1286 |
< |
eCos[3][i].Vmul(u, eCos[2][i]); |
| 1287 |
< |
eSin[3][i].Vmul(u, eSin[2][i]); |
| 1285 |
> |
|
| 1286 |
> |
u = eCos[2][i]; |
| 1287 |
> |
w = eSin[2][i]; |
| 1288 |
|
|
| 1289 |
|
for(int l = 3; l <= kLimit; l++) { |
| 1290 |
< |
w.Vmul(u, eCos[l-1][i]); |
| 1291 |
< |
eCos[l][i] = w - eCos[l-2][i]; |
| 1292 |
< |
w.Vmul(u, eSin[l-1][i]); |
| 1293 |
< |
eSin[l][i] = w - eSin[l-2][i]; |
| 1290 |
> |
eCos[l][i].x() = eCos[l-1][i].x()*eCos[2][i].x() - eSin[l-1][i].x()*eSin[2][i].x(); |
| 1291 |
> |
eCos[l][i].y() = eCos[l-1][i].y()*eCos[2][i].y() - eSin[l-1][i].y()*eSin[2][i].y(); |
| 1292 |
> |
eCos[l][i].z() = eCos[l-1][i].z()*eCos[2][i].z() - eSin[l-1][i].z()*eSin[2][i].z(); |
| 1293 |
> |
|
| 1294 |
> |
eSin[l][i].x() = eSin[l-1][i].x()*eCos[2][i].x() + eCos[l-1][i].x()*eSin[2][i].x(); |
| 1295 |
> |
eSin[l][i].y() = eSin[l-1][i].y()*eCos[2][i].y() + eCos[l-1][i].y()*eSin[2][i].y(); |
| 1296 |
> |
eSin[l][i].z() = eSin[l-1][i].z()*eCos[2][i].z() + eCos[l-1][i].z()*eSin[2][i].z(); |
| 1297 |
> |
|
| 1298 |
> |
|
| 1299 |
> |
// a.Vmul(eCos[l-1][i], u); |
| 1300 |
> |
// b.Vmul(eSin[l-1][i], w); |
| 1301 |
> |
// eCos[l][i] = a - b; |
| 1302 |
> |
// a.Vmul(eSin[l-1][i], u); |
| 1303 |
> |
// b.Vmul(eCos[l-1][i], w); |
| 1304 |
> |
// eSin[l][i] = a + b; |
| 1305 |
> |
|
| 1306 |
|
} |
| 1307 |
|
} |
| 1308 |
|
} |
| 1350 |
|
int mMin = kLimit; |
| 1351 |
|
int nMin = kLimit + 1; |
| 1352 |
|
for (int l = 1; l <= kLimit; l++) { |
| 1353 |
< |
int ll =l - 1; |
| 1353 |
> |
int ll = l - 1; |
| 1354 |
|
RealType rl = xcl * float(ll); |
| 1355 |
|
for (int mmm = mMin; mmm <= kLim2; mmm++) { |
| 1356 |
|
int mm = mmm - kLimit; |
| 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() |
| 1373 |
|
- eSin[l][i].x()*eSin[m][i].y(); |
| 1374 |
|
slm[i] = eSin[l][i].x()*eCos[m][i].y() |
| 1375 |
< |
+ eSin[m][i].y()*eCos[l][i].x(); |
| 1375 |
> |
+ eSin[m][i].y()*eCos[l][i].x(); |
| 1376 |
|
} |
| 1377 |
|
} |
| 1378 |
|
} |
| 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) { |
| 1423 |
|
Vector3d D = atom->getDipole() * mPoleConverter; |
| 1424 |
< |
RealType dk = dot(kVec, D); |
| 1425 |
< |
dxk[i] = cross(kVec, D); |
| 1424 |
> |
RealType dk = dot(D, kVec); |
| 1425 |
> |
dxk[i] = cross(D, kVec); |
| 1426 |
|
dkc[i] = dk * ckr[i]; |
| 1427 |
|
dks[i] = dk * skr[i]; |
| 1428 |
|
} |
| 1429 |
|
if (data.is_Quadrupole) { |
| 1430 |
|
Mat3x3d Q = atom->getQuadrupole(); |
| 1431 |
|
Q *= mPoleConverter; |
| 1432 |
< |
RealType qk = -( Q * k2 ).trace(); |
| 1432 |
> |
RealType qk = - doubleDot(Q, k2); |
| 1433 |
> |
// RealType qk = -( Q * k2 ).trace(); |
| 1434 |
|
qxk[i] = -2.0 * cross(k2, Q); |
| 1435 |
|
qkc[i] = qk * ckr[i]; |
| 1436 |
|
qks[i] = qk * skr[i]; |
| 1437 |
|
} |
| 1438 |
|
} |
| 1439 |
|
} |
| 1440 |
< |
|
| 1440 |
> |
|
| 1441 |
|
// calculate vector sums |
| 1442 |
|
|
| 1443 |
|
RealType ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); |
| 1446 |
|
RealType dkss = std::accumulate(dks.begin(),dks.end(),0.0); |
| 1447 |
|
RealType qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); |
| 1448 |
|
RealType qkss = std::accumulate(qks.begin(),qks.end(),0.0); |
| 1449 |
+ |
|
| 1450 |
|
|
| 1451 |
|
#ifdef IS_MPI |
| 1452 |
|
MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, |
| 1464 |
|
#endif |
| 1465 |
|
|
| 1466 |
|
// Accumulate potential energy and virial contribution: |
| 1467 |
< |
|
| 1448 |
< |
//cerr << "l, m, n = " << l << " " << m << " " << n << "\n"; |
| 1449 |
< |
cerr << "kVec = " << kVec << "\n"; |
| 1450 |
< |
cerr << "ckss = " << ckss << " ckcs = " << ckcs << "\n"; |
| 1467 |
> |
|
| 1468 |
|
kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) |
| 1469 |
|
+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkss)); |
| 1470 |
< |
//cerr << "kspace pot = " << kPot << "\n"; |
| 1470 |
> |
|
| 1471 |
|
kVir -= 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss |
| 1472 |
|
+4.0*(ckss*dkcs-ckcs*dkss) |
| 1473 |
|
+3.0*(dkcs*dkcs+dkss*dkss) |
| 1487 |
|
ElectrostaticAtomData data = ElectrostaticMap[Etids[atid]]; |
| 1488 |
|
|
| 1489 |
|
RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs) |
| 1490 |
< |
- (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); |
| 1490 |
> |
- (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); |
| 1491 |
|
RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) |
| 1492 |
|
-ckr[i]*(ckss+dkcs-qkss)); |
| 1493 |
|
RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs)+ |
| 1494 |
|
skr[i]*(ckss+dkcs-qkss)); |
| 1495 |
< |
|
| 1479 |
< |
|
| 1495 |
> |
|
| 1496 |
|
atom->addFrc( 4.0 * rvol * qfrc * kVec ); |
| 1497 |
|
|
| 1498 |
|
if (data.is_Dipole) { |
| 1505 |
|
} |
| 1506 |
|
} |
| 1507 |
|
} |
| 1508 |
+ |
nMin = 1; |
| 1509 |
|
} |
| 1510 |
+ |
mMin = 1; |
| 1511 |
|
} |
| 1512 |
+ |
cerr << "kPot = " << kPot << "\n"; |
| 1513 |
+ |
pot[ELECTROSTATIC_FAMILY] += kPot; |
| 1514 |
|
} |
| 1515 |
|
} |
