| 106 |
|
angstromToM_ = 1.0e-10; |
| 107 |
|
debyeToCm_ = 3.33564095198e-30; |
| 108 |
|
|
| 109 |
< |
// number of points for electrostatic splines |
| 110 |
< |
np_ = 1000; |
| 109 |
> |
// Default number of points for electrostatic splines |
| 110 |
> |
np_ = 140; |
| 111 |
|
|
| 112 |
|
// variables to handle different summation methods for long-range |
| 113 |
|
// electrostatics: |
| 246 |
|
b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; |
| 247 |
|
b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; |
| 248 |
|
b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; |
| 249 |
< |
selfMult_ = b0c + 2.0 * a2 * invArootPi; |
| 249 |
> |
selfMult_ = b0c + a2 * invArootPi; |
| 250 |
|
} else { |
| 251 |
|
a2 = 0.0; |
| 252 |
|
b0c = 1.0 / r; |
| 279 |
|
// working variables for Taylor expansion: |
| 280 |
|
RealType rmRc, rmRc2, rmRc3, rmRc4; |
| 281 |
|
|
| 282 |
+ |
// Approximate using splines using a maximum of 0.1 Angstroms |
| 283 |
+ |
// between points. |
| 284 |
+ |
int nptest = int((cutoffRadius_ + 2.0) / 0.1); |
| 285 |
+ |
np_ = (np_ > nptest) ? np_ : nptest; |
| 286 |
+ |
|
| 287 |
|
// Add a 2 angstrom safety window to deal with cutoffGroups that |
| 288 |
|
// have charged atoms longer than the cutoffRadius away from each |
| 289 |
|
// other. Splining is almost certainly the best choice here. |
| 690 |
|
Vector3d D_a, D_b; // Dipoles (space-fixed) |
| 691 |
|
Mat3x3d Q_a, Q_b; // Quadrupoles (space-fixed) |
| 692 |
|
|
| 693 |
< |
RealType ri, ri2, ri3, ri4; // Distance utility scalars |
| 693 |
> |
RealType ri; // Distance utility scalar |
| 694 |
|
RealType rdDa, rdDb; // Dipole utility scalars |
| 695 |
|
Vector3d rxDa, rxDb; // Dipole utility vectors |
| 696 |
|
RealType rdQar, rdQbr, trQa, trQb; // Quadrupole utility scalars |
| 700 |
|
RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars |
| 701 |
|
RealType rQaQbr; |
| 702 |
|
Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors |
| 703 |
< |
Vector3d rQaQb, QaQbr, QaxQb, QbxQa, rQaxQbr, QbQar, rQbxQar; |
| 704 |
< |
Mat3x3d QaQb, QbQa; // Cross-interaction matrices |
| 703 |
> |
Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; |
| 704 |
> |
Mat3x3d QaQb; // Cross-interaction matrices |
| 705 |
|
|
| 706 |
|
RealType U(0.0); // Potential |
| 707 |
|
Vector3d F(0.0); // Force |
| 735 |
|
|
| 736 |
|
ri = 1.0 / *(idat.rij); |
| 737 |
|
Vector3d rhat = *(idat.d) * ri; |
| 733 |
– |
ri2 = ri * ri; |
| 738 |
|
|
| 739 |
|
// logicals |
| 740 |
|
|
| 796 |
|
v45 = v45s->getValueAt( *(idat.rij) ); |
| 797 |
|
v46 = v46s->getValueAt( *(idat.rij) ); |
| 798 |
|
} |
| 795 |
– |
|
| 799 |
|
|
| 800 |
|
// calculate the single-site contributions (fields, etc). |
| 801 |
|
|
| 1014 |
|
Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32); |
| 1015 |
|
} |
| 1016 |
|
if (b_is_Quadrupole) { |
| 1017 |
< |
pref = pre44_ * *(idat.electroMult); |
| 1017 |
> |
pref = pre44_ * *(idat.electroMult); // yes |
| 1018 |
|
QaQb = Q_a * Q_b; |
| 1016 |
– |
QbQa = Q_b * Q_a; |
| 1019 |
|
trQaQb = QaQb.trace(); |
| 1020 |
|
rQaQb = rhat * QaQb; |
| 1021 |
< |
QbQar = QbQa * rhat; |
| 1020 |
< |
QaQbr = QaQb * rhat; |
| 1021 |
> |
QaQbr = QaQb * rhat; |
| 1022 |
|
QaxQb = cross(Q_a, Q_b); |
| 1022 |
– |
QbxQa = cross(Q_b, Q_a); |
| 1023 |
|
rQaQbr = dot(rQa, Qbr); |
| 1024 |
|
rQaxQbr = cross(rQa, Qbr); |
| 1025 |
– |
//rQbxQar = cross(rQb, Qar); |
| 1026 |
– |
|
| 1027 |
– |
|
| 1028 |
– |
// First part of potential and associated forces: |
| 1029 |
– |
// U += pref * (trQa * trQb) * v41; |
| 1030 |
– |
// F += pref * rhat * (trQa * trQb) * v44; |
| 1031 |
– |
// Ta += 0.0; |
| 1032 |
– |
// Tb += 0.0; |
| 1033 |
– |
|
| 1034 |
– |
// cerr << "Aa:\n"; |
| 1035 |
– |
// cerr << *(idat.A1) << "\n\n"; |
| 1036 |
– |
|
| 1037 |
– |
// cerr << "Ab:\n"; |
| 1038 |
– |
// cerr << *(idat.A2) << "\n\n"; |
| 1039 |
– |
|
| 1040 |
– |
// cerr << "Qa:\n"; |
| 1041 |
– |
// cerr << Q_a << "\n\n"; |
| 1042 |
– |
// cerr << "Qb:\n"; |
| 1043 |
– |
// cerr << Q_b << "\n\n"; |
| 1044 |
– |
|
| 1045 |
– |
|
| 1046 |
– |
// Second part of potential and associated forces: |
| 1047 |
– |
// Probably suspect (particularly torques): |
| 1048 |
– |
// cerr << "trQaQb = " << trQaQb << "\n"; |
| 1049 |
– |
|
| 1050 |
– |
// U += pref * 2.0 * trQaQb * v41; |
| 1051 |
– |
// F += pref * rhat * (2.0 * trQaQb) * v44; |
| 1052 |
– |
// Ta += pref * (-4.0*QaxQb) * v41; |
| 1053 |
– |
// Tb += pref * (-4.0*QbxQa) * v41; |
| 1054 |
– |
|
| 1055 |
– |
// cerr << "QaxQb = " << QaxQb << "\n"; |
| 1056 |
– |
// cerr << "QbxQa = " << QbxQa << "\n"; |
| 1057 |
– |
|
| 1058 |
– |
// Third part of potential: |
| 1059 |
– |
// U += pref * (trQa * rdQbr) * v42; |
| 1060 |
– |
// F += pref * rhat * (trQa * rdQbr) * v45; |
| 1061 |
– |
// F += pref * 2.0 * (trQa * rQb) * v44; |
| 1062 |
– |
// Ta += 0.0; |
| 1063 |
– |
// Tb += pref * 2.0 * trQa * cross(rhat, Qbr) * v42; |
| 1025 |
|
|
| 1026 |
|
U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; |
| 1027 |
|
U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; |
| 1044 |
|
|
| 1045 |
|
Tb += pref * (+ 4.0 * QaxQb * v41); |
| 1046 |
|
Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) |
| 1047 |
< |
+ 4.0 * cross(rhat, QbQar) |
| 1047 |
> |
- 4.0 * cross(rQaQb, rhat) |
| 1048 |
|
+ 4.0 * rQaxQbr) * v42; |
| 1049 |
+ |
// Possible replacement for line 2 above: |
| 1050 |
+ |
// + 4.0 * cross(rhat, QbQar) |
| 1051 |
+ |
|
| 1052 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1053 |
|
|
| 1090 |
– |
// Tb += pref * (+ 4.0 * QaxQb * v41); |
| 1091 |
– |
// Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) |
| 1092 |
– |
// - 4.0 * cross(rQaQb, rhat) |
| 1093 |
– |
// + 4.0 * rQaxQbr) * v42; |
| 1094 |
– |
// Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1095 |
– |
|
| 1054 |
|
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1055 |
|
} |
| 1056 |
|
} |
| 1092 |
|
*(idat.t2) += *(idat.sw) * indirect_Tb; |
| 1093 |
|
} |
| 1094 |
|
return; |
| 1095 |
< |
} |
| 1095 |
> |
} |
| 1096 |
|
|
| 1097 |
|
void Electrostatic::calcSelfCorrection(SelfData &sdat) { |
| 1098 |
|
|
| 1135 |
|
|
| 1136 |
|
case esm_SHIFTED_FORCE: |
| 1137 |
|
case esm_SHIFTED_POTENTIAL: |
| 1138 |
< |
if (i_is_Charge) { |
| 1139 |
< |
self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; |
| 1138 |
> |
if (i_is_Charge) { |
| 1139 |
> |
self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; |
| 1140 |
|
(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; |
| 1141 |
|
} |
| 1142 |
|
break; |
