| 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, QbxQa, rQaxQbr, QbQar, rQbxQar; |
| 699 |
< |
Mat3x3d QaQb, QbQa; // Cross-interaction matrices |
| 698 |
> |
Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; |
| 699 |
> |
Mat3x3d QaQb; // Cross-interaction matrices |
| 700 |
|
|
| 701 |
|
RealType U(0.0); // Potential |
| 702 |
|
Vector3d F(0.0); // Force |
| 793 |
|
v46 = v46s->getValueAt( *(idat.rij) ); |
| 794 |
|
} |
| 795 |
|
|
| 796 |
– |
|
| 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; |
| 1016 |
– |
QbQa = Q_b * Q_a; |
| 1015 |
|
trQaQb = QaQb.trace(); |
| 1016 |
|
rQaQb = rhat * QaQb; |
| 1017 |
< |
QbQar = QbQa * rhat; |
| 1020 |
< |
QaQbr = QaQb * rhat; |
| 1017 |
> |
QaQbr = QaQb * rhat; |
| 1018 |
|
QaxQb = cross(Q_a, Q_b); |
| 1022 |
– |
QbxQa = cross(Q_b, Q_a); |
| 1019 |
|
rQaQbr = dot(rQa, Qbr); |
| 1020 |
|
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; |
| 1021 |
|
|
| 1022 |
|
U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; |
| 1023 |
|
U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; |
| 1040 |
|
|
| 1041 |
|
Tb += pref * (+ 4.0 * QaxQb * v41); |
| 1042 |
|
Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) |
| 1043 |
< |
+ 4.0 * cross(rhat, QbQar) |
| 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 |
|
|
| 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 |
– |
|
| 1050 |
|
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1051 |
|
} |
| 1052 |
|
} |
