--- trunk/src/brains/SimInfo.cpp	2006/09/22 22:19:59	1050
+++ trunk/src/brains/SimInfo.cpp	2007/04/20 18:15:48	1129
@@ -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"
@@ -90,7 +92,8 @@ 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), calcBoxDipole_(false) {
+    sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), 
+    useAtomicVirial_(true) {
 
       MoleculeStamp* molStamp;
       int nMolWithSameStamp;
@@ -664,18 +667,20 @@ namespace oopse {
     int useSF;
     int useSP;
     int useBoxDipole;
+
     std::string myMethod;
 
     // 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"){
@@ -687,6 +692,8 @@ namespace oopse {
       if (simParams_->getAccumulateBoxDipole())
 	useBoxDipole = 1;
 
+    useAtomicVirial_ = simParams_->getUseAtomicVirial();
+
     //loop over all of the atom types
     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
       useLennardJones |= (*i)->isLennardJones();
@@ -764,6 +771,9 @@ namespace oopse {
     temp = useBoxDipole;
     MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); 
 
+    temp = useAtomicVirial_;
+    MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); 
+
 #endif
 
     fInfo_.SIM_uses_PBC = usePBC;    
@@ -783,6 +793,7 @@ namespace oopse {
     fInfo_.SIM_uses_SF = useSF;
     fInfo_.SIM_uses_SP = useSP;
     fInfo_.SIM_uses_BoxDipole = useBoxDipole;
+    fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_;
   }
 
   void SimInfo::setupFortranSim() {
@@ -867,6 +878,14 @@ 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()) {
+      int nlistNeighbors = simParams_->getNeighborListNeighbors();
+      setNeighbors(&nlistNeighbors);
+    }
+   
+
   }
 
 
@@ -936,6 +955,10 @@ namespace oopse {
     // Check the cutoff policy
     int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
 
+    // Set LJ shifting bools to false
+    ljsp_ = false;
+    ljsf_ = false;
+
     std::string myPolicy;
     if (forceFieldOptions_.haveCutoffPolicy()){
       myPolicy = forceFieldOptions_.getCutoffPolicy();
@@ -999,9 +1022,19 @@ namespace oopse {
 	  simError();
 	}
       }
+
+      if (simParams_->haveElectrostaticSummationMethod()) {
+	std::string myMethod = simParams_->getElectrostaticSummationMethod();
+	toUpper(myMethod);
+        
+	if (myMethod == "SHIFTED_POTENTIAL") {
+	  ljsp_ = true;
+	} else if (myMethod == "SHIFTED_FORCE") {
+	  ljsf_ = true;
+	}
+      }
+      notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
       
-      notifyFortranCutoffs(&rcut_, &rsw_);
-      
     } else {
       
       // For electrostatic atoms, we'll assume a large safe value:
@@ -1017,7 +1050,15 @@ namespace oopse {
         if (simParams_->haveElectrostaticSummationMethod()) {
           std::string myMethod = simParams_->getElectrostaticSummationMethod();
           toUpper(myMethod);
-          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
+      
+      // For the time being, we're tethering the LJ shifted behavior to the
+      // electrostaticSummationMethod keyword options
+	  if (myMethod == "SHIFTED_POTENTIAL") {
+	    ljsp_ = true;
+	  } else if (myMethod == "SHIFTED_FORCE") {
+	    ljsf_ = true;
+	  }
+	  if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
             if (simParams_->haveSwitchingRadius()){
               sprintf(painCave.errMsg,
                       "SimInfo Warning: A value was set for the switchingRadius\n"
@@ -1040,7 +1081,9 @@ namespace oopse {
           simError();
           rsw_ = 0.85 * rcut_;
         }
-        notifyFortranCutoffs(&rcut_, &rsw_);
+
+        notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
+
       } else {
         // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
         // We'll punt and let fortran figure out the cutoffs later.
@@ -1125,7 +1168,10 @@ namespace oopse {
                      "\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,
@@ -1442,6 +1488,43 @@ namespace oopse {
     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_);