--- branches/development/src/nonbonded/Electrostatic.cpp	2013/01/07 20:05:43	1821
+++ branches/development/src/nonbonded/Electrostatic.cpp	2013/01/08 15:29:03	1822
@@ -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, QbxQa, 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,56 +1010,14 @@ 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);
-        QbxQa = cross(Q_b, Q_a);
         rQaQbr = dot(rQa, Qbr);
         rQaxQbr = cross(rQa, Qbr);
-        //rQbxQar = cross(rQb, Qar);
-
-
-        // First part of potential and associated forces:
-        // U += pref * (trQa * trQb) * v41;
-        // F += pref * rhat * (trQa * trQb) * v44;
-        // Ta += 0.0;
-        // Tb += 0.0;
-
-        // cerr << "Aa:\n";
-        // cerr << *(idat.A1)  << "\n\n";
-
-        // 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;
@@ -1083,16 +1040,13 @@ namespace OpenMD {
 
         Tb += pref * (+ 4.0 * QaxQb * v41);
         Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) 
-                      + 4.0 * cross(rhat, QbQar) 
+                      - 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";
       }
     }