--- trunk/src/brains/SimInfo.cpp	2005/05/30 14:01:52	555
+++ trunk/src/brains/SimInfo.cpp	2005/09/18 20:45:38	610
@@ -52,7 +52,10 @@
 #include "brains/SimInfo.hpp"
 #include "math/Vector3.hpp"
 #include "primitives/Molecule.hpp"
+#include "UseTheForce/fCutoffPolicy.h"
+#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
 #include "UseTheForce/doForces_interface.h"
+#include "UseTheForce/DarkSide/electrostatic_interface.h"
 #include "UseTheForce/notifyCutoffs_interface.h"
 #include "utils/MemoryUtils.hpp"
 #include "utils/simError.h"
@@ -462,7 +465,9 @@ namespace oopse {
     //setup fortran force field
     /** @deprecate */    
     int isError = 0;
-    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
+    
+    setupElectrostaticSummationMethod( isError );
+
     if(isError){
       sprintf( painCave.errMsg,
 	       "ForceField error: There was an error initializing the forceField in fortran.\n" );
@@ -518,7 +523,6 @@ namespace oopse {
     int useElectrostatics = 0;
     //usePBC and useRF are from simParams
     int usePBC = simParams_->getPBC();
-    int useRF = simParams_->getUseRF();
 
     //loop over all of the atom types
     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
@@ -581,9 +585,6 @@ namespace oopse {
     temp = useFLARB;
     MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
 
-    temp = useRF;
-    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
-    
 #endif
 
     fInfo_.SIM_uses_PBC = usePBC;    
@@ -598,7 +599,6 @@ namespace oopse {
     fInfo_.SIM_uses_EAM = useEAM;
     fInfo_.SIM_uses_Shapes = useShape;
     fInfo_.SIM_uses_FLARB = useFLARB;
-    fInfo_.SIM_uses_RF = useRF;
 
     if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
 
@@ -830,14 +830,87 @@ namespace oopse {
     }
   }
 
-  void SimInfo::setupCutoff() {
+  void SimInfo::setupCutoff() {    
     getCutoff(rcut_, rsw_);    
     double rnblist = rcut_ + 1; // skin of neighbor list
 
     //Pass these cutoff radius etc. to fortran. This function should be called once and only once
-    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
+    
+    int cp =  TRADITIONAL_CUTOFF_POLICY;
+    if (simParams_->haveCutoffPolicy()) {
+      std::string myPolicy = simParams_->getCutoffPolicy();
+      if (myPolicy == "MIX") {
+        cp = MIX_CUTOFF_POLICY;
+      } else {
+        if (myPolicy == "MAX") {
+          cp = MAX_CUTOFF_POLICY;
+        } else {
+          if (myPolicy == "TRADITIONAL") {            
+            cp = TRADITIONAL_CUTOFF_POLICY;
+          } else {
+            // throw error        
+            sprintf( painCave.errMsg,
+                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
+            painCave.isFatal = 1;
+            simError();
+          }     
+        }           
+      }
+    }
+    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
+    // also send cutoff notification to electrostatics
+    setElectrostaticCutoffRadius(&rcut_);
   }
 
+  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
+     
+    int errorOut;
+    int esm =  NONE;
+    double alphaVal;
+    double dielectric;
+
+    errorOut = isError;
+    alphaVal = simParams_->getDampingAlpha();
+    dielectric = simParams_->getDielectric();
+
+    if (simParams_->haveElectrostaticSummationMethod()) {
+      std::string myMethod = simParams_->getElectrostaticSummationMethod();
+      if (myMethod == "NONE") {
+        esm = NONE;
+      } else {
+        if (myMethod == "UNDAMPED_WOLF") {
+          esm = UNDAMPED_WOLF;
+        } else {
+          if (myMethod == "DAMPED_WOLF") {            
+            esm = DAMPED_WOLF;
+	    if (!simParams_->haveDampingAlpha()) {
+	      //throw error
+	      sprintf( painCave.errMsg,
+		       "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal);
+	      painCave.isFatal = 0;
+	      simError();
+	    }
+          } else {
+	    if (myMethod == "REACTION_FIELD") {
+	      esm = REACTION_FIELD;
+	    } else {
+	      // throw error        
+	      sprintf( painCave.errMsg,
+		       "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() );
+	      painCave.isFatal = 1;
+	      simError();
+	    }     
+	  }           
+	}
+      }
+    }
+    // let's pass some summation method variables to fortran
+    setElectrostaticSummationMethod( &esm );
+    setDampedWolfAlpha( &alphaVal );
+    setReactionFieldDielectric( &dielectric );
+    initFortranFF( &esm, &errorOut );
+  }
+
   void SimInfo::addProperty(GenericData* genData) {
     properties_.addProperty(genData);  
   }
@@ -981,6 +1054,11 @@ namespace oopse {
    
    /* 
    Return intertia tensor for entire system and angular momentum Vector.
+
+
+       [  Ixx -Ixy  -Ixz ]
+  J =| -Iyx  Iyy  -Iyz |
+       [ -Izx -Iyz   Izz ]
     */
 
    void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
@@ -1032,7 +1110,7 @@ namespace oopse {
       inertiaTensor(0,1) = -xy;
       inertiaTensor(0,2) = -xz;
       inertiaTensor(1,0) = -xy;
-      inertiaTensor(2,0) = xx + zz;
+      inertiaTensor(1,1) = xx + zz;
       inertiaTensor(1,2) = -yz;
       inertiaTensor(2,0) = -xz;
       inertiaTensor(2,1) = -yz;
@@ -1060,17 +1138,18 @@ namespace oopse {
       SimInfo::MoleculeIterator i;
       Molecule* mol;
       
-      Vector3d thisq(0.0);
-      Vector3d thisv(0.0);
+      Vector3d thisr(0.0);
+      Vector3d thisp(0.0);
       
-      double thisMass = 0.0;
+      double thisMass;
       
       for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {         
-         thisq = mol->getCom()-com;
-         thisv = mol->getComVel()-comVel;
-         thisMass = mol->getMass();
-         angularMomentum += cross( thisq, thisv ) * thisMass;
+        thisMass = mol->getMass(); 
+	thisr = mol->getCom()-com;
+	thisp = (mol->getComVel()-comVel)*thisMass;
          
+	angularMomentum += cross( thisr, thisp );
+         
       }  
        
 #ifdef IS_MPI