--- branches/development/src/nonbonded/Electrostatic.cpp	2012/12/18 16:23:13	1820
+++ branches/development/src/nonbonded/Electrostatic.cpp	2013/01/08 21:02:27	1823
@@ -106,8 +106,8 @@ namespace OpenMD {
     angstromToM_ = 1.0e-10;
     debyeToCm_ = 3.33564095198e-30;
     
-    // number of points for electrostatic splines
-    np_ = 100;
+    // Default number of points for electrostatic splines
+    np_ = 140;
     
     // variables to handle different summation methods for long-range 
     // electrostatics:
@@ -246,7 +246,7 @@ namespace OpenMD {
       b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
       b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
       b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
-      selfMult_ = b0c  +  2.0 * a2 * invArootPi;
+      selfMult_ = b0c + a2 * invArootPi;
     } else {
       a2 = 0.0;
       b0c = 1.0 / r;
@@ -279,6 +279,11 @@ namespace OpenMD {
     // working variables for Taylor expansion:
     RealType rmRc, rmRc2, rmRc3, rmRc4;
 
+    // Approximate using splines using a maximum of 0.1 Angstroms
+    // between points.
+    int nptest = int((cutoffRadius_ + 2.0) / 0.1);
+    np_ = (np_ > nptest) ? np_ : nptest;
+  
     // Add a 2 angstrom safety window to deal with cutoffGroups that
     // have charged atoms longer than the cutoffRadius away from each
     // other.  Splining is almost certainly the best choice here.
@@ -695,8 +700,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 +798,6 @@ namespace OpenMD {
       v46 = v46s->getValueAt( *(idat.rij) );
     }
 
-
     // calculate the single-site contributions (fields, etc).
     
     if (a_is_Charge) {
@@ -1011,18 +1015,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 +1045,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";
       }
     }
 
@@ -1095,7 +1093,7 @@ namespace OpenMD {
         *(idat.t2) += *(idat.sw) * indirect_Tb;
     }
     return;
-  }  
+  }
     
   void Electrostatic::calcSelfCorrection(SelfData &sdat) {
 
@@ -1138,8 +1136,8 @@ namespace OpenMD {
       
     case esm_SHIFTED_FORCE:
     case esm_SHIFTED_POTENTIAL:
-      if (i_is_Charge) {        
-        self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
+      if (i_is_Charge) {
+        self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
         (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
       }
       break;