--- branches/development/src/nonbonded/Electrostatic.cpp	2012/12/18 16:23:13	1820
+++ branches/development/src/nonbonded/Electrostatic.cpp	2013/01/08 15:29:03	1822
@@ -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:
@@ -695,8 +695,8 @@ namespace OpenMD {
     RealType DadDb, trQaQb, DadQbr, DbdQar;       // Cross-interaction scalars
     RealType rQaQbr;
     Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors
-    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr, QbQar, rQbxQar;
-    Mat3x3d  QaQb, QbQa;                          // Cross-interaction matrices
+    Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr;
+    Mat3x3d  QaQb;                                // Cross-interaction matrices
 
     RealType U(0.0);  // Potential
     Vector3d F(0.0);  // Force
@@ -793,7 +793,6 @@ namespace OpenMD {
       v46 = v46s->getValueAt( *(idat.rij) );
     }
 
-
     // calculate the single-site contributions (fields, etc).
     
     if (a_is_Charge) {
@@ -1011,18 +1010,15 @@ namespace OpenMD {
         Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
       }
       if (b_is_Quadrupole) {
-        pref = pre44_ * *(idat.electroMult);
+        pref = pre44_ * *(idat.electroMult);  // yes
         QaQb = Q_a * Q_b;
-        QbQa = Q_b * Q_a;
         trQaQb = QaQb.trace();
         rQaQb = rhat * QaQb;
-        QbQar = QbQa * rhat;
-        QaQbr = QaQb * rhat;        
+        QaQbr = QaQb * rhat;
         QaxQb = cross(Q_a, Q_b);
         rQaQbr = dot(rQa, Qbr);
         rQaxQbr = cross(rQa, Qbr);
-        rQbxQar = cross(rQb, Qar);
-
+        
         U  += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
         U  += pref * (trQa * rdQbr + trQb * rdQar  + 4.0 * rQaQbr) * v42;
         U  += pref * (rdQar * rdQbr) * v43;
@@ -1044,17 +1040,14 @@ namespace OpenMD {
 
         Tb += pref * (+ 4.0 * QaxQb * v41);
         Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) 
-                      + 4.0 * cross(rhat, QbQar) 
-                      - 4.0 * rQbxQar) * v42;
+                      - 4.0 * cross(rQaQb, rhat) 
+                      + 4.0 * rQaxQbr) * v42;
+        // Possible replacement for line 2 above:
+        //             + 4.0 * cross(rhat, QbQar) 
+
         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";
+        //  cerr << " tsum = " << Ta + Tb - cross(  *(idat.d) , F ) << "\n";
       }
     }