| 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: |
| 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, rQaxQbr, QbQar, rQbxQar; |
| 698 |
> |
Vector3d rQaQb, QaQbr, QaxQb, QbxQa, rQaxQbr, QbQar, rQbxQar; |
| 699 |
|
Mat3x3d QaQb, QbQa; // Cross-interaction matrices |
| 700 |
|
|
| 701 |
|
RealType U(0.0); // Potential |
| 1019 |
|
QbQar = QbQa * rhat; |
| 1020 |
|
QaQbr = QaQb * rhat; |
| 1021 |
|
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); |
| 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; |
| 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; |
| 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 * rQbxQar) * v42; |
| 1087 |
> |
+ 4.0 * rQaxQbr) * v42; |
| 1088 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1089 |
|
|
| 1090 |
|
// Tb += pref * (+ 4.0 * QaxQb * v41); |
| 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"; |
| 1096 |
> |
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1097 |
|
} |
| 1098 |
|
} |
| 1099 |
|
|
