| 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; |
| 699 |
< |
Mat3x3d QaQb; // Cross-interaction matrices |
| 698 |
> |
Vector3d rQaQb, QaQbr, QaxQb, QbxQa, rQaxQbr, QbQar, rQbxQar; |
| 699 |
> |
Mat3x3d QaQb, QbQa; // Cross-interaction matrices |
| 700 |
|
|
| 701 |
|
RealType U(0.0); // Potential |
| 702 |
|
Vector3d F(0.0); // Force |
| 1013 |
|
if (b_is_Quadrupole) { |
| 1014 |
|
pref = pre44_ * *(idat.electroMult); |
| 1015 |
|
QaQb = Q_a * Q_b; |
| 1016 |
+ |
QbQa = Q_b * Q_a; |
| 1017 |
|
trQaQb = QaQb.trace(); |
| 1018 |
|
rQaQb = rhat * QaQb; |
| 1019 |
+ |
QbQar = QbQa * rhat; |
| 1020 |
|
QaQbr = QaQb * rhat; |
| 1021 |
|
QaxQb = cross(Q_a, Q_b); |
| 1022 |
< |
|
| 1023 |
< |
U += pref * (trQa * trQb + 2.0*trQaQb) * v41; |
| 1024 |
< |
U += pref * (trQa*rdQbr + trQb*rdQar + 4.0*dot(rQa, Qbr)) * v42; |
| 1025 |
< |
U += pref * (rdQar * rdQbr) * v43; |
| 1022 |
> |
QbxQa = cross(Q_b, Q_a); |
| 1023 |
> |
rQaQbr = dot(rQa, Qbr); |
| 1024 |
> |
rQaxQbr = cross(rQa, Qbr); |
| 1025 |
> |
//rQbxQar = cross(rQb, Qar); |
| 1026 |
|
|
| 1027 |
< |
F += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44; |
| 1028 |
< |
F += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44; |
| 1029 |
< |
F += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44; |
| 1030 |
< |
F += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat |
| 1031 |
< |
+ 4.0*dot(rQa, Qbr)*rhat)*v45; |
| 1032 |
< |
F += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45; |
| 1030 |
< |
F += pref * (rdQar*rdQbr*rhat) * v46; |
| 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 |
< |
Ta += pref * (-4.0 * QaxQb * v41); |
| 1035 |
< |
Ta += pref * (-2.0*trQb*cross(rQa, rhat) |
| 1034 |
< |
+ 4.0*cross(rhat, QaQbr) |
| 1035 |
< |
- 4.0*cross(rQa, Qbr)) * v42; |
| 1036 |
< |
Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43; |
| 1034 |
> |
// cerr << "Aa:\n"; |
| 1035 |
> |
// cerr << *(idat.A1) << "\n\n"; |
| 1036 |
|
|
| 1037 |
< |
Tb += pref * (4.0 * QaxQb * v41); |
| 1038 |
< |
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; |
| 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; |
| 1064 |
+ |
|
| 1065 |
+ |
U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; |
| 1066 |
+ |
U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; |
| 1067 |
+ |
U += pref * (rdQar * rdQbr) * v43; |
| 1068 |
+ |
|
| 1069 |
+ |
F += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44; |
| 1070 |
+ |
F += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45; |
| 1071 |
+ |
F += pref * rhat * (rdQar * rdQbr)*v46; |
| 1072 |
+ |
|
| 1073 |
+ |
F += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44; |
| 1074 |
+ |
F += pref * 4.0 * (rQaQb + QaQbr)*v44; |
| 1075 |
+ |
F += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45; |
| 1076 |
+ |
|
| 1077 |
+ |
Ta += pref * (- 4.0 * QaxQb * v41); |
| 1078 |
+ |
Ta += pref * (- 2.0 * trQb * cross(rQa, rhat) |
| 1079 |
+ |
+ 4.0 * cross(rhat, QaQbr) |
| 1080 |
+ |
- 4.0 * rQaxQbr) * v42; |
| 1081 |
+ |
Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43; |
| 1082 |
+ |
|
| 1083 |
+ |
|
| 1084 |
+ |
Tb += pref * (+ 4.0 * QaxQb * v41); |
| 1085 |
+ |
Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) |
| 1086 |
+ |
+ 4.0 * cross(rhat, QbQar) |
| 1087 |
+ |
+ 4.0 * rQaxQbr) * v42; |
| 1088 |
+ |
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1089 |
+ |
|
| 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 |
+ |
|
| 1096 |
+ |
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1097 |
|
} |
| 1098 |
|
} |
| 1099 |
|
|
