--- trunk/src/brains/SimInfo.cpp	2006/05/17 21:51:42	963
+++ trunk/src/brains/SimInfo.cpp	2006/09/22 22:19:59	1050
@@ -53,6 +53,7 @@
 #include "brains/SimInfo.hpp"
 #include "math/Vector3.hpp"
 #include "primitives/Molecule.hpp"
+#include "primitives/StuntDouble.hpp"
 #include "UseTheForce/fCutoffPolicy.h"
 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
@@ -89,7 +90,7 @@ namespace oopse {
     nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
     nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
     nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
-    sman_(NULL), fortranInitialized_(false) {
+    sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) {
 
       MoleculeStamp* molStamp;
       int nMolWithSameStamp;
@@ -600,8 +601,11 @@ namespace oopse {
     /** @deprecate */    
     int isError = 0;
     
+    setupCutoff();
+    
     setupElectrostaticSummationMethod( isError );
     setupSwitchingFunction();
+    setupAccumulateBoxDipole();
 
     if(isError){
       sprintf( painCave.errMsg,
@@ -609,9 +613,6 @@ namespace oopse {
       painCave.isFatal = 1;
       simError();
     }
-  
-    
-    setupCutoff();
 
     calcNdf();
     calcNdfRaw();
@@ -661,6 +662,8 @@ namespace oopse {
     int usePBC = simParams_->getUsePeriodicBoundaryConditions();
     int useRF;
     int useSF;
+    int useSP;
+    int useBoxDipole;
     std::string myMethod;
 
     // set the useRF logical
@@ -671,14 +674,18 @@ namespace oopse {
     if (simParams_->haveElectrostaticSummationMethod()) {
       std::string myMethod = simParams_->getElectrostaticSummationMethod();
       toUpper(myMethod);
-      if (myMethod == "REACTION_FIELD") {
+      if (myMethod == "REACTION_FIELD"){
         useRF=1;
-      } else {
-	if (myMethod == "SHIFTED_FORCE") {
-	  useSF = 1;
-	}
+      } else if (myMethod == "SHIFTED_FORCE"){
+	useSF = 1;
+      } else if (myMethod == "SHIFTED_POTENTIAL"){
+	useSP = 1;
       }
     }
+    
+    if (simParams_->haveAccumulateBoxDipole()) 
+      if (simParams_->getAccumulateBoxDipole())
+	useBoxDipole = 1;
 
     //loop over all of the atom types
     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
@@ -749,8 +756,14 @@ namespace oopse {
     MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
 
     temp = useSF;
-    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
+    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);   
+
+    temp = useSP;
+    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
 
+    temp = useBoxDipole;
+    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); 
+
 #endif
 
     fInfo_.SIM_uses_PBC = usePBC;    
@@ -768,21 +781,8 @@ namespace oopse {
     fInfo_.SIM_uses_FLARB = useFLARB;
     fInfo_.SIM_uses_RF = useRF;
     fInfo_.SIM_uses_SF = useSF;
-
-    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();
-      }      
-    }
-
+    fInfo_.SIM_uses_SP = useSP;
+    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
   }
 
   void SimInfo::setupFortranSim() {
@@ -1058,10 +1058,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();
@@ -1078,8 +1076,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,
@@ -1106,10 +1113,16 @@ 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();
 	  }
@@ -1158,6 +1171,17 @@ namespace oopse {
 
     // send switching function notification to switcheroo
     setFunctionType(&ft);
+
+  }
+
+  void SimInfo::setupAccumulateBoxDipole() {    
+
+    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
+    if ( simParams_->haveAccumulateBoxDipole() ) 
+      if ( simParams_->getAccumulateBoxDipole() ) {
+	setAccumulateBoxDipole();
+	calcBoxDipole_ = true;
+      }
 
   }
 
@@ -1409,7 +1433,25 @@ namespace oopse {
       
       return angularMomentum;
    }
-   
    
+  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
+    return IOIndexToIntegrableObject.at(index);
+  }
+  
+  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
+    IOIndexToIntegrableObject= v;
+  }
+
+/*
+   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
+      assert( v.size() == nAtoms_ + nRigidBodies_);
+      sdByGlobalIndex_ = v;
+    }
+
+    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
+      //assert(index < nAtoms_ + nRigidBodies_);
+      return sdByGlobalIndex_.at(index);
+    }   
+*/   
 }//end namespace oopse