--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/04/12 20:32:20	1097
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/05/11 20:33:41	1157
@@ -147,13 +147,6 @@ void SimSetup::createSim(void){
   // make the output filenames
 
   makeOutNames();
-
-  if (globals->haveMinimizer())
-    // make minimizer
-    makeMinimizer();
-  else
-    // make the integrator
-    makeIntegrator();
   
 #ifdef IS_MPI
   mpiSim->mpiRefresh();
@@ -162,12 +155,20 @@ void SimSetup::createSim(void){
   // initialize the Fortran
 
   initFortran();
+
+  if (globals->haveMinimizer())
+    // make minimizer
+    makeMinimizer();
+  else
+    // make the integrator
+    makeIntegrator();
+
 }
 
 
 void SimSetup::makeMolecules(void){
   int i, j, k;
-  int exI, exJ, exK, exL, slI;
+  int exI, exJ, exK, exL, slI, slJ;
   int tempI, tempJ, tempK, tempL;
   int molI;
   int stampID, atomOffset, rbOffset;
@@ -182,7 +183,9 @@ void SimSetup::makeMolecules(void){
   BendStamp* currentBend;
   TorsionStamp* currentTorsion;
   RigidBodyStamp* currentRigidBody;
-
+  CutoffGroupStamp* currentCutoffGroup;
+  CutoffGroup* myCutoffGroup;
+  
   bond_pair* theBonds;
   bend_set* theBends;
   torsion_set* theTorsions;
@@ -190,6 +193,8 @@ void SimSetup::makeMolecules(void){
   set<int> skipList;
 
   double phi, theta, psi;
+  char* molName;
+  char rbName[100];
 
   //init the forceField paramters
 
@@ -206,12 +211,14 @@ void SimSetup::makeMolecules(void){
 
     for (i = 0; i < info[k].n_mol; i++){
       stampID = info[k].molecules[i].getStampID();
+      molName = comp_stamps[stampID]->getID();
 
       molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
       molInfo.nBonds = comp_stamps[stampID]->getNBonds();
       molInfo.nBends = comp_stamps[stampID]->getNBends();
       molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
       molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
+      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
       
       molInfo.myAtoms = &(info[k].atoms[atomOffset]);
 
@@ -259,13 +266,13 @@ void SimSetup::makeMolecules(void){
         else{
 
           molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
+
         }
 
         molInfo.myAtoms[j]->setType(currentAtom->getType());
-
 #ifdef IS_MPI
 
-        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
+        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
 
 #endif // is_mpi
       } 
@@ -406,6 +413,9 @@ void SimSetup::makeMolecules(void){
         info[k].excludes->addPair(exK, exL);
       }
 
+      
+      molInfo.myRigidBodies.clear();
+      
       for (j = 0; j < molInfo.nRigidBodies; j++){
 
         currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
@@ -414,6 +424,9 @@ void SimSetup::makeMolecules(void){
         // Create the Rigid Body:
 
         myRB = new RigidBody();
+
+        sprintf(rbName,"%s_RB_%d", molName, j);
+        myRB->setType(rbName);
         
         for (rb1 = 0; rb1 < nMembers; rb1++) {
 
@@ -454,19 +467,80 @@ void SimSetup::makeMolecules(void){
             // used for the exclude list:
             
 #ifdef IS_MPI
-            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
-            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
+            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
+            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
 #else
-            exI = tempI + 1;
-            exJ = tempJ + 1;
+            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
+            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
 #endif
             
             info[k].excludes->addPair(exI, exJ);
             
           }
         }
+
+        molInfo.myRigidBodies.push_back(myRB);
+        info[k].rigidBodies.push_back(myRB);
       }
       
+
+      //creat cutoff group for molecule
+      molInfo.myCutoffGroups.clear();
+      for (j = 0; j < molInfo.nCutoffGroups; j++){
+
+        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
+        nMembers = currentCutoffGroup->getNMembers(); 
+
+        myCutoffGroup = new CutoffGroup();
+        
+        for (int cg = 0; cg < nMembers; cg++) {
+
+          // molI is atom numbering inside this molecule
+          molI = currentCutoffGroup->getMember(cg);     
+
+          // tempI is atom numbering on local processor
+          tempI = molI + atomOffset;
+
+          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
+        }
+
+        molInfo.myCutoffGroups.push_back(myCutoffGroup);
+      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
+      
+
+
+      // After this is all set up, scan through the atoms to 
+      // see if they can be added to the integrableObjects:
+
+      molInfo.myIntegrableObjects.clear();
+      
+
+      for (j = 0; j < molInfo.nAtoms; j++){
+
+#ifdef IS_MPI
+        slJ = molInfo.myAtoms[j]->getGlobalIndex();
+#else
+        slJ = j+atomOffset;
+#endif
+
+        // if they aren't on the skip list, then they can be integrated
+
+        if (skipList.find(slJ) == skipList.end()) {
+          mySD = (StuntDouble *) molInfo.myAtoms[j];
+          info[k].integrableObjects.push_back(mySD);
+          molInfo.myIntegrableObjects.push_back(mySD);
+        }
+      }
+
+      // all rigid bodies are integrated:
+
+      for (j = 0; j < molInfo.nRigidBodies; j++) {
+        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
+        info[k].integrableObjects.push_back(mySD);      
+        molInfo.myIntegrableObjects.push_back(mySD);
+      }
+    
+      
       // send the arrays off to the forceField for init.
       
       the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
@@ -482,28 +556,7 @@ void SimSetup::makeMolecules(void){
       delete[] theBonds;
       delete[] theBends;
       delete[] theTorsions;
-    }
-
-    // build up the integrableObjects vector:
-
-    for (i = 0; i < info[k].n_atoms; i++) {
-      
-#ifdef IS_MPI
-      slI = info[k].atoms[i]->getGlobalIndex();
-#else
-      slI = i;
-#endif
-
-      if (skipList.find(slI) == skipList.end()) {
-        mySD = (StuntDouble *) info[k].atoms[i];
-        info[k].integrableObjects.push_back(mySD);
-      }
-    }
-    for (i = 0; i < info[k].rigidBodies.size(); i++) {
-      mySD = (StuntDouble *) info[k].rigidBodies[i];
-      info[k].integrableObjects.push_back(mySD);      
-    }
-    
+    }    
   }
 
 #ifdef IS_MPI
@@ -513,13 +566,13 @@ void SimSetup::makeMolecules(void){
 
   // clean up the forcefield
 
-  if (!globals->haveLJrcut()){
+  if (!globals->haveRcut()){
 
     the_ff->calcRcut();
 
   } else {
     
-    the_ff->setRcut( globals->getLJrcut() );
+    the_ff->setRcut( globals->getRcut() );
   }
 
   the_ff->cleanMe();
@@ -809,7 +862,7 @@ void SimSetup::gatherInfo(void){
   }
 
   //check whether sample time, status time, thermal time and reset time are divisble by dt
-  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
+  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
     sprintf(painCave.errMsg,
             "Sample time is not divisible by dt.\n"
             "\tThis will result in samples that are not uniformly\n"
@@ -819,7 +872,7 @@ void SimSetup::gatherInfo(void){
     simError();    
   }
 
-  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
+  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
     sprintf(painCave.errMsg,
             "Status time is not divisible by dt.\n"
             "\tThis will result in status reports that are not uniformly\n"
@@ -855,12 +908,10 @@ void SimSetup::gatherInfo(void){
     if (globals->haveSampleTime()){
       info[i].sampleTime = globals->getSampleTime();
       info[i].statusTime = info[i].sampleTime;
-      info[i].thermalTime = info[i].sampleTime;
     }
     else{
       info[i].sampleTime = globals->getRunTime();
       info[i].statusTime = info[i].sampleTime;
-      info[i].thermalTime = info[i].sampleTime;
     }
 
     if (globals->haveStatusTime()){
@@ -869,6 +920,8 @@ void SimSetup::gatherInfo(void){
 
     if (globals->haveThermalTime()){
       info[i].thermalTime = globals->getThermalTime();
+    } else {
+      info[i].thermalTime = globals->getRunTime();
     }
 
     info[i].resetIntegrator = 0;
@@ -939,6 +992,7 @@ void SimSetup::finalInfoCheck(void){
 void SimSetup::finalInfoCheck(void){
   int index;
   int usesDipoles;
+  int usesCharges;
   int i;
 
   for (i = 0; i < nInfo; i++){
@@ -950,45 +1004,49 @@ void SimSetup::finalInfoCheck(void){
       usesDipoles = (info[i].atoms[index])->hasDipole();
       index++;
     }
-
+    index = 0;
+    usesCharges = 0;
+    while ((index < info[i].n_atoms) && !usesCharges){
+      usesCharges= (info[i].atoms[index])->hasCharge();
+      index++;
+    }
 #ifdef IS_MPI
     int myUse = usesDipoles;
     MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
 #endif //is_mpi
 
-    double theEcr, theEst;
+    double theRcut, theRsw;
 
     if (globals->getUseRF()){
       info[i].useReactionField = 1;
 
-      if (!globals->haveECR()){
+      if (!globals->haveRcut()){
         sprintf(painCave.errMsg,
-                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
+                "SimSetup Warning: No value was set for the cutoffRadius.\n"
                 "\tOOPSE will use a default value of 15.0 angstroms"
-                "\tfor the electrostaticCutoffRadius.\n");
+                "\tfor the cutoffRadius.\n");
         painCave.isFatal = 0;
         simError();
-        theEcr = 15.0;
+        theRcut = 15.0;
       }
       else{
-        theEcr = globals->getECR();
+        theRcut = globals->getRcut();
       }
 
-      if (!globals->haveEST()){
+      if (!globals->haveRsw()){
         sprintf(painCave.errMsg,
-                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
+                "SimSetup Warning: No value was set for switchingRadius.\n"
                 "\tOOPSE will use a default value of\n"
-                "\t0.05 * electrostaticCutoffRadius\n"
-                "\tfor the electrostaticSkinThickness\n");
+                "\t0.95 * cutoffRadius for the switchingRadius\n");
         painCave.isFatal = 0;
         simError();
-        theEst = 0.05 * theEcr;
+        theRsw = 0.95 * theRcut;
       }
       else{
-        theEst = globals->getEST();
+        theRsw = globals->getRsw();
       }
 
-      info[i].setDefaultEcr(theEcr, theEst);
+      info[i].setDefaultRcut(theRcut, theRsw);
 
       if (!globals->haveDielectric()){
         sprintf(painCave.errMsg,
@@ -1001,35 +1059,36 @@ void SimSetup::finalInfoCheck(void){
       info[i].dielectric = globals->getDielectric();
     }
     else{
-      if (usesDipoles){
-        if (!globals->haveECR()){
+      if (usesDipoles || usesCharges){
+
+        if (!globals->haveRcut()){
           sprintf(painCave.errMsg,
-                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
+                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
                   "\tOOPSE will use a default value of 15.0 angstroms"
-                  "\tfor the electrostaticCutoffRadius.\n");
-	  painCave.isFatal = 0;
-	  simError();
-	  theEcr = 15.0;
-        }
-        else{
-          theEcr = globals->getECR();
-        }
-	
-        if (!globals->haveEST()){
+                  "\tfor the cutoffRadius.\n");
+          painCave.isFatal = 0;
+          simError();
+          theRcut = 15.0;
+      }
+        else{
+          theRcut = globals->getRcut();
+        }
+        
+        if (!globals->haveRsw()){
           sprintf(painCave.errMsg,
-                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
+                  "SimSetup Warning: No value was set for switchingRadius.\n"
                   "\tOOPSE will use a default value of\n"
-                  "\t0.05 * electrostaticCutoffRadius\n"
-                  "\tfor the electrostaticSkinThickness\n");
+                  "\t0.95 * cutoffRadius for the switchingRadius\n");
           painCave.isFatal = 0;
           simError();
-          theEst = 0.05 * theEcr;
+          theRsw = 0.95 * theRcut;
         }
         else{
-          theEst = globals->getEST();
+          theRsw = globals->getRsw();
         }
+        
+        info[i].setDefaultRcut(theRcut, theRsw);
 	
-        info[i].setDefaultEcr(theEcr, theEst);
       }
     }
   }
@@ -1251,7 +1310,10 @@ void SimSetup::compList(void){
   LinkedMolStamp* headStamp = new LinkedMolStamp();
   LinkedMolStamp* currentStamp = NULL;
   comp_stamps = new MoleculeStamp * [n_components];
+  bool haveCutoffGroups;
 
+  haveCutoffGroups = false;
+  
   // make an array of molecule stamps that match the components used.
   // also extract the used stamps out into a separate linked list
 
@@ -1286,7 +1348,13 @@ void SimSetup::compList(void){
       headStamp->add(currentStamp);
       comp_stamps[i] = headStamp->match(id);
     }
+
+    if(comp_stamps[i]->getNCutoffGroups() > 0)
+      haveCutoffGroups = true;    
   }
+    
+  for (i = 0; i < nInfo; i++)
+    info[i].haveCutoffGroups = haveCutoffGroups;
 
 #ifdef IS_MPI
   strcpy(checkPointMsg, "Component stamps successfully extracted\n");
@@ -1334,10 +1402,13 @@ void SimSetup::mpiMolDivide(void){
   int localMol, allMol;
   int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
   int local_rigid;
+  vector<int> globalMolIndex;
 
   mpiSim = new mpiSimulation(info);
 
-  globalIndex = mpiSim->divideLabor();
+  mpiSim->divideLabor();
+  globalAtomIndex = mpiSim->getGlobalAtomIndex();
+  //globalMolIndex = mpiSim->getGlobalMolIndex();
 
   // set up the local variables 
 
@@ -1351,7 +1422,7 @@ void SimSetup::mpiMolDivide(void){
   local_bends = 0;
   local_torsions = 0;
   local_rigid = 0;
-  globalAtomIndex = 0;
+  globalAtomCounter = 0;
 
   for (i = 0; i < n_components; i++){
     for (j = 0; j < components_nmol[i]; j++){
@@ -1364,8 +1435,8 @@ void SimSetup::mpiMolDivide(void){
         localMol++;
       }      
       for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
-        info[0].molMembershipArray[globalAtomIndex] = allMol;
-        globalAtomIndex++;
+        info[0].molMembershipArray[globalAtomCounter] = allMol;
+        globalAtomCounter++;
       }
 
       allMol++;
@@ -1433,15 +1504,15 @@ void SimSetup::makeSysArrays(void){
 #else // is_mpi
 
     molIndex = 0;
-    globalAtomIndex = 0;
+    globalAtomCounter = 0;
     for (i = 0; i < n_components; i++){
       for (j = 0; j < components_nmol[i]; j++){
         the_molecules[molIndex].setStampID(i);
         the_molecules[molIndex].setMyIndex(molIndex);
         the_molecules[molIndex].setGlobalIndex(molIndex);
         for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
-          info[l].molMembershipArray[globalAtomIndex] = molIndex;
-          globalAtomIndex++;
+          info[l].molMembershipArray[globalAtomCounter] = molIndex;
+          globalAtomCounter++;
         }
         molIndex++;
       }
@@ -1458,7 +1529,7 @@ void SimSetup::makeSysArrays(void){
     info[l].atoms = the_atoms;
     info[l].molecules = the_molecules;
     info[l].nGlobalExcludes = 0;
-
+    
     the_ff->setSimInfo(info);
   }
 }