| 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) { |
| 293 |
< |
selfMult1_ *= 2.0; |
| 294 |
< |
selfMult2_ *= 2.0; |
| 295 |
< |
selfMult4_ *= 2.0; |
| 296 |
< |
} else { |
| 292 |
> |
if (summationMethod_ != esm_EWALD_FULL) { |
| 293 |
|
selfMult1_ -= b0c; |
| 294 |
|
selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; |
| 295 |
|
selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; |
| 1219 |
|
Vector3d box = hmat.diagonals(); |
| 1220 |
|
RealType boxMax = box.max(); |
| 1221 |
|
|
| 1226 |
– |
cerr << "da = " << dampingAlpha_ << " rc = " << cutoffRadius_ << "\n"; |
| 1227 |
– |
cerr << "boxMax = " << boxMax << "\n"; |
| 1222 |
|
//int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); |
| 1223 |
< |
int kMax = 5; |
| 1230 |
< |
cerr << "kMax = " << kMax << "\n"; |
| 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; |
| 1275 |
|
r = atom->getPos(); |
| 1276 |
|
info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); |
| 1277 |
|
|
| 1284 |
– |
// Shift so that all coordinates are in the range [0,L]: |
| 1285 |
– |
|
| 1286 |
– |
r += box/2.0; |
| 1287 |
– |
|
| 1278 |
|
tt.Vmul(t, r); |
| 1279 |
|
|
| 1290 |
– |
//cerr << "tt = " << tt << "\n"; |
| 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 |
+ |
|
| 1286 |
|
u = eCos[2][i]; |
| 1287 |
|
w = eSin[2][i]; |
| 1288 |
|
|
| 1289 |
|
for(int l = 3; l <= kLimit; l++) { |
| 1290 |
< |
a.Vmul(eCos[l-1][i], u); |
| 1291 |
< |
b.Vmul(eSin[l-1][i], w); |
| 1292 |
< |
eCos[l][i] = a - b; |
| 1293 |
< |
a.Vmul(eSin[l-1][i], u); |
| 1294 |
< |
b.Vmul(eCos[l-1][i], w); |
| 1295 |
< |
eSin[l][i] = a + b; |
| 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; |
| 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 |
|
} |
| 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]; |
| 1463 |
|
MPI::SUM); |
| 1464 |
|
#endif |
| 1465 |
|
|
| 1466 |
< |
// Accumulate potential energy and virial contribution: |
| 1466 |
> |
// Accumulate potential energy and virial contribution: |
| 1467 |
|
|
| 1468 |
|
kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) |
| 1469 |
|
+ (ckcs-dkss-qkcs)*(ckcs-dkss-qkss)); |
| 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 |
< |
|
| 1495 |
< |
|
| 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; |
