--- trunk/src/brains/SimInfo.cpp	2006/08/30 18:42:29	1024
+++ trunk/src/brains/SimInfo.cpp	2006/12/29 20:21:53	1103
@@ -59,6 +59,7 @@
 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
 #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
 #include "UseTheForce/doForces_interface.h"
+#include "UseTheForce/DarkSide/neighborLists_interface.h"
 #include "UseTheForce/DarkSide/electrostatic_interface.h"
 #include "UseTheForce/DarkSide/switcheroo_interface.h"
 #include "utils/MemoryUtils.hpp"
@@ -67,6 +68,7 @@
 #include "io/ForceFieldOptions.hpp"
 #include "UseTheForce/ForceField.hpp"
 
+
 #ifdef IS_MPI
 #include "UseTheForce/mpiComponentPlan.h"
 #include "UseTheForce/DarkSide/simParallel_interface.h"
@@ -600,6 +602,8 @@ namespace oopse {
     //setup fortran force field
     /** @deprecate */    
     int isError = 0;
+    
+    setupCutoff();
     
     setupElectrostaticSummationMethod( isError );
     setupSwitchingFunction();
@@ -611,9 +615,6 @@ namespace oopse {
       painCave.isFatal = 1;
       simError();
     }
-  
-    
-    setupCutoff();
 
     calcNdf();
     calcNdfRaw();
@@ -670,13 +671,14 @@ namespace oopse {
     // set the useRF logical
     useRF = 0;
     useSF = 0;
+    useSP = 0;
 
 
     if (simParams_->haveElectrostaticSummationMethod()) {
       std::string myMethod = simParams_->getElectrostaticSummationMethod();
       toUpper(myMethod);
       if (myMethod == "REACTION_FIELD"){
-        useRF=1;
+        useRF = 1;
       } else if (myMethod == "SHIFTED_FORCE"){
 	useSF = 1;
       } else if (myMethod == "SHIFTED_POTENTIAL"){
@@ -784,21 +786,6 @@ namespace oopse {
     fInfo_.SIM_uses_SF = useSF;
     fInfo_.SIM_uses_SP = useSP;
     fInfo_.SIM_uses_BoxDipole = useBoxDipole;
-
-    if( myMethod == "REACTION_FIELD") {
-      
-      if (simParams_->haveDielectric()) {
-	fInfo_.dielect = simParams_->getDielectric();
-      } else {
-	sprintf(painCave.errMsg,
-		"SimSetup Error: No Dielectric constant was set.\n"
-		"\tYou are trying to use Reaction Field without"
-		"\tsetting a dielectric constant!\n");
-	painCave.isFatal = 1;
-	simError();
-      }      
-    }
-
   }
 
   void SimInfo::setupFortranSim() {
@@ -883,6 +870,13 @@ namespace oopse {
 	     "succesfully sent the simulation information to fortran.\n");
     MPIcheckPoint();
 #endif // is_mpi
+
+    // Setup number of neighbors in neighbor list if present
+    if (simParams_->haveNeighborListNeighbors()) {
+      setNeighbors(simParams_->getNeighborListNeighbors());
+    }
+   
+
   }
 
 
@@ -1074,10 +1068,8 @@ namespace oopse {
     int sm = UNDAMPED;
     RealType alphaVal;
     RealType dielectric;
-
+    
     errorOut = isError;
-    alphaVal = simParams_->getDampingAlpha();
-    dielectric = simParams_->getDielectric();
 
     if (simParams_->haveElectrostaticSummationMethod()) {
       std::string myMethod = simParams_->getElectrostaticSummationMethod();
@@ -1094,8 +1086,17 @@ namespace oopse {
 	    if (myMethod == "SHIFTED_FORCE") {            
 	      esm = SHIFTED_FORCE;
 	    } else {
-	      if (myMethod == "REACTION_FIELD") {	      
+	      if (myMethod == "REACTION_FIELD") {
 		esm = REACTION_FIELD;
+		dielectric = simParams_->getDielectric();
+		if (!simParams_->haveDielectric()) {
+		  // throw warning
+		  sprintf( painCave.errMsg,
+			   "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
+			   "\tA default value of %f will be used for the dielectric.\n", dielectric);
+		  painCave.isFatal = 0;
+		  simError();
+		}
 	      } else {
 		// throw error        
 		sprintf( painCave.errMsg,
@@ -1122,13 +1123,22 @@ namespace oopse {
 	if (myScreen == "DAMPED") {
 	  sm = DAMPED;
 	  if (!simParams_->haveDampingAlpha()) {
-	    //throw error
+	    // first set a cutoff dependent alpha value
+	    // we assume alpha depends linearly with rcut from 0 to 20.5 ang
+	    alphaVal = 0.5125 - rcut_* 0.025;
+	    // for values rcut > 20.5, alpha is zero
+	    if (alphaVal < 0) alphaVal = 0;
+
+	    // throw warning
 	    sprintf( painCave.errMsg,
 		     "SimInfo warning: dampingAlpha was not specified in the input file.\n"
-                     "\tA default value of %f (1/ang) will be used.\n", alphaVal);
+                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
 	    painCave.isFatal = 0;
 	    simError();
+	  } else {
+	    alphaVal = simParams_->getDampingAlpha();
 	  }
+	  
 	} else {
 	  // throw error        
 	  sprintf( painCave.errMsg,
@@ -1443,8 +1453,45 @@ namespace oopse {
   
   void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
     IOIndexToIntegrableObject= v;
+  }
+
+  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes 
+     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
+     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to 
+     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
+  */
+  void SimInfo::getGyrationalVolume(RealType &volume){
+    Mat3x3d intTensor;
+    RealType det;
+    Vector3d dummyAngMom; 
+    RealType sysconstants;
+    RealType geomCnst;
+
+    geomCnst = 3.0/2.0;
+    /* Get the inertial tensor and angular momentum for free*/
+    getInertiaTensor(intTensor,dummyAngMom);
+    
+    det = intTensor.determinant();
+    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
+    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
+    return;
   }
 
+  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
+    Mat3x3d intTensor;
+    Vector3d dummyAngMom; 
+    RealType sysconstants;
+    RealType geomCnst;
+
+    geomCnst = 3.0/2.0;
+    /* Get the inertial tensor and angular momentum for free*/
+    getInertiaTensor(intTensor,dummyAngMom);
+    
+    detI = intTensor.determinant();
+    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
+    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
+    return;
+  }
 /*
    void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
       assert( v.size() == nAtoms_ + nRigidBodies_);