--- branches/development/src/nonbonded/Electrostatic.cpp	2012/11/27 21:13:48	1814
+++ branches/development/src/nonbonded/Electrostatic.cpp	2013/01/07 20:05:43	1821
@@ -107,7 +107,7 @@ namespace OpenMD {
     debyeToCm_ = 3.33564095198e-30;
     
     // number of points for electrostatic splines
-    np_ = 100;
+    np_ = 1000;
     
     // variables to handle different summation methods for long-range 
     // electrostatics:
@@ -499,7 +499,7 @@ namespace OpenMD {
       v33v.push_back(v33);
       v34v.push_back(v34);
       v35v.push_back(v35);
-
+      
       v41v.push_back(v41);
       v42v.push_back(v42);
       v43v.push_back(v43);
@@ -693,9 +693,10 @@ namespace OpenMD {
     RealType pref;
 
     RealType DadDb, trQaQb, DadQbr, DbdQar;       // Cross-interaction scalars
+    RealType rQaQbr;
     Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors
-    Vector3d rQaQb, QaQbr, QaxQb;
-    Mat3x3d  QaQb;                                // Cross-interaction matrices
+    Vector3d rQaQb, QaQbr, QaxQb, QbxQa, rQaxQbr, QbQar, rQbxQar;
+    Mat3x3d  QaQb, QbQa;                          // Cross-interaction matrices
 
     RealType U(0.0);  // Potential
     Vector3d F(0.0);  // Force
@@ -1012,35 +1013,87 @@ namespace OpenMD {
       if (b_is_Quadrupole) {
         pref = pre44_ * *(idat.electroMult);
         QaQb = Q_a * Q_b;
+        QbQa = Q_b * Q_a;
         trQaQb = QaQb.trace();
         rQaQb = rhat * QaQb;
+        QbQar = QbQa * rhat;
         QaQbr = QaQb * rhat;        
         QaxQb = cross(Q_a, Q_b);
-
-        U  += pref * (trQa * trQb + 2.0*trQaQb) * v41;
-        U  += pref * (trQa*rdQbr + trQb*rdQar  + 4.0*dot(rQa, Qbr)) * v42;
-        U  += pref * (rdQar * rdQbr) * v43;
+        QbxQa = cross(Q_b, Q_a);
+        rQaQbr = dot(rQa, Qbr);
+        rQaxQbr = cross(rQa, Qbr);
+        //rQbxQar = cross(rQb, Qar);
 
-        F  += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44;
-        F  += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44;
-        F  += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44;
-        F  += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat 
-                      + 4.0*dot(rQa, Qbr)*rhat)*v45;
-        F  += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45;
-        F  += pref * (rdQar*rdQbr*rhat) * v46;
+
+        // First part of potential and associated forces:
+        // U += pref * (trQa * trQb) * v41;
+        // F += pref * rhat * (trQa * trQb) * v44;
+        // Ta += 0.0;
+        // Tb += 0.0;
 
-        Ta += pref * (-4.0 * QaxQb  * v41);
-        Ta += pref * (-2.0*trQb*cross(rQa, rhat) 
-                      + 4.0*cross(rhat, QaQbr) 
-                      - 4.0*cross(rQa, Qbr)) * v42;
-        Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
+        // cerr << "Aa:\n";
+        // cerr << *(idat.A1)  << "\n\n";
 
-        Tb += pref * (4.0 * QaxQb * v41);
-        Tb += pref * (-2.0*trQa*cross(rQb, rhat) 
-                      - 4.0*cross(rQaQb, rhat) 
-                      + 4.0*cross(rQa, Qbr))*v42;
-        Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;        
+        // cerr << "Ab:\n";
+        // cerr << *(idat.A2)  << "\n\n";
 
+        // cerr << "Qa:\n";
+        // cerr << Q_a << "\n\n";
+        // cerr << "Qb:\n";
+        // cerr << Q_b << "\n\n";
+
+         
+        // Second part of potential and associated forces:
+        // Probably suspect (particularly torques):
+        // cerr << "trQaQb = " << trQaQb << "\n";
+
+        // U += pref * 2.0 * trQaQb * v41;
+        // F += pref * rhat * (2.0 * trQaQb) * v44;
+        // Ta += pref * (-4.0*QaxQb) * v41;
+        // Tb += pref * (-4.0*QbxQa) * v41;
+
+        // cerr << "QaxQb = " << QaxQb << "\n";
+        // cerr << "QbxQa = " << QbxQa << "\n";
+
+        // Third part of potential:
+        // U += pref * (trQa * rdQbr) * v42;
+        // F += pref * rhat * (trQa * rdQbr) * v45;
+        // F += pref * 2.0 * (trQa * rQb) * v44;
+        // Ta += 0.0;
+        // Tb += pref * 2.0 * trQa * cross(rhat, Qbr) * v42;
+        
+        U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
+        U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
+        U  += pref * (rdQar * rdQbr) * v43;
+
+        F  += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*v44;
+        F  += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr)*v45;
+        F  += pref * rhat * (rdQar * rdQbr)*v46;
+
+        F  += pref * 2.0 * (trQb*rQa + trQa*rQb)*v44;
+        F  += pref * 4.0 * (rQaQb + QaQbr)*v44;
+        F  += pref * 2.0 * (rQa*rdQbr + rdQar*rQb)*v45;
+
+        Ta += pref * (- 4.0 * QaxQb * v41);
+        Ta += pref * (- 2.0 * trQb * cross(rQa, rhat) 
+                      + 4.0 * cross(rhat, QaQbr) 
+                      - 4.0 * rQaxQbr) * v42;
+        Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
+
+
+        Tb += pref * (+ 4.0 * QaxQb * v41);
+        Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) 
+                      + 4.0 * cross(rhat, QbQar) 
+                      + 4.0 * rQaxQbr) * v42;
+        Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
+
+        // Tb += pref * (+ 4.0 * QaxQb * v41);
+        // Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) 
+        //               - 4.0 * cross(rQaQb, rhat) 
+        //               + 4.0 * rQaxQbr) * v42;
+        // Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
+
+        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
       }
     }