| 107 |
|
debyeToCm_ = 3.33564095198e-30; |
| 108 |
|
|
| 109 |
|
// number of points for electrostatic splines |
| 110 |
< |
np_ = 100; |
| 110 |
> |
np_ = 1000; |
| 111 |
|
|
| 112 |
|
// variables to handle different summation methods for long-range |
| 113 |
|
// electrostatics: |
| 499 |
|
v33v.push_back(v33); |
| 500 |
|
v34v.push_back(v34); |
| 501 |
|
v35v.push_back(v35); |
| 502 |
< |
|
| 502 |
> |
|
| 503 |
|
v41v.push_back(v41); |
| 504 |
|
v42v.push_back(v42); |
| 505 |
|
v43v.push_back(v43); |
| 693 |
|
RealType pref; |
| 694 |
|
|
| 695 |
|
RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars |
| 696 |
+ |
RealType rQaQbr; |
| 697 |
|
Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors |
| 698 |
< |
Vector3d rQaQb, QaQbr, QaxQb; |
| 698 |
> |
Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; |
| 699 |
|
Mat3x3d QaQb; // Cross-interaction matrices |
| 700 |
|
|
| 701 |
|
RealType U(0.0); // Potential |
| 793 |
|
v46 = v46s->getValueAt( *(idat.rij) ); |
| 794 |
|
} |
| 795 |
|
|
| 795 |
– |
|
| 796 |
|
// calculate the single-site contributions (fields, etc). |
| 797 |
|
|
| 798 |
|
if (a_is_Charge) { |
| 1010 |
|
Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32); |
| 1011 |
|
} |
| 1012 |
|
if (b_is_Quadrupole) { |
| 1013 |
< |
pref = pre44_ * *(idat.electroMult); |
| 1013 |
> |
pref = pre44_ * *(idat.electroMult); // yes |
| 1014 |
|
QaQb = Q_a * Q_b; |
| 1015 |
|
trQaQb = QaQb.trace(); |
| 1016 |
|
rQaQb = rhat * QaQb; |
| 1017 |
< |
QaQbr = QaQb * rhat; |
| 1017 |
> |
QaQbr = QaQb * rhat; |
| 1018 |
|
QaxQb = cross(Q_a, Q_b); |
| 1019 |
< |
|
| 1020 |
< |
U += pref * (trQa * trQb + 2.0*trQaQb) * v41; |
| 1021 |
< |
U += pref * (trQa*rdQbr + trQb*rdQar + 4.0*dot(rQa, Qbr)) * v42; |
| 1019 |
> |
rQaQbr = dot(rQa, Qbr); |
| 1020 |
> |
rQaxQbr = cross(rQa, Qbr); |
| 1021 |
> |
|
| 1022 |
> |
U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; |
| 1023 |
> |
U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; |
| 1024 |
|
U += pref * (rdQar * rdQbr) * v43; |
| 1025 |
|
|
| 1026 |
< |
F += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44; |
| 1027 |
< |
F += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44; |
| 1028 |
< |
F += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44; |
| 1027 |
< |
F += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat |
| 1028 |
< |
+ 4.0*dot(rQa, Qbr)*rhat)*v45; |
| 1029 |
< |
F += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45; |
| 1030 |
< |
F += pref * (rdQar*rdQbr*rhat) * v46; |
| 1026 |
> |
F += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44; |
| 1027 |
> |
F += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45; |
| 1028 |
> |
F += pref * rhat * (rdQar * rdQbr)*v46; |
| 1029 |
|
|
| 1030 |
< |
Ta += pref * (-4.0 * QaxQb * v41); |
| 1031 |
< |
Ta += pref * (-2.0*trQb*cross(rQa, rhat) |
| 1032 |
< |
+ 4.0*cross(rhat, QaQbr) |
| 1033 |
< |
- 4.0*cross(rQa, Qbr)) * v42; |
| 1030 |
> |
F += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44; |
| 1031 |
> |
F += pref * 4.0 * (rQaQb + QaQbr)*v44; |
| 1032 |
> |
F += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45; |
| 1033 |
> |
|
| 1034 |
> |
Ta += pref * (- 4.0 * QaxQb * v41); |
| 1035 |
> |
Ta += pref * (- 2.0 * trQb * cross(rQa, rhat) |
| 1036 |
> |
+ 4.0 * cross(rhat, QaQbr) |
| 1037 |
> |
- 4.0 * rQaxQbr) * v42; |
| 1038 |
|
Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43; |
| 1039 |
|
|
| 1038 |
– |
Tb += pref * (4.0 * QaxQb * v41); |
| 1039 |
– |
Tb += pref * (-2.0*trQa*cross(rQb, rhat) |
| 1040 |
– |
- 4.0*cross(rQaQb, rhat) |
| 1041 |
– |
+ 4.0*cross(rQa, Qbr))*v42; |
| 1042 |
– |
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1040 |
|
|
| 1041 |
+ |
Tb += pref * (+ 4.0 * QaxQb * v41); |
| 1042 |
+ |
Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) |
| 1043 |
+ |
- 4.0 * cross(rQaQb, rhat) |
| 1044 |
+ |
+ 4.0 * rQaxQbr) * v42; |
| 1045 |
+ |
// Possible replacement for line 2 above: |
| 1046 |
+ |
// + 4.0 * cross(rhat, QbQar) |
| 1047 |
+ |
|
| 1048 |
+ |
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1049 |
+ |
|
| 1050 |
+ |
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1051 |
|
} |
| 1052 |
|
} |
| 1053 |
|
|
