--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/05/12 14:30:12	1163
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/05/27 18:59:17	1203
@@ -185,7 +185,9 @@ void SimSetup::makeMolecules(void){
   RigidBodyStamp* currentRigidBody;
   CutoffGroupStamp* currentCutoffGroup;
   CutoffGroup* myCutoffGroup;
-  
+  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
+  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
+
   bond_pair* theBonds;
   bend_set* theBends;
   torsion_set* theTorsions;
@@ -218,22 +220,23 @@ void SimSetup::makeMolecules(void){
       molInfo.nBends = comp_stamps[stampID]->getNBends();
       molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
       molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
-      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
+
+      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
       
       molInfo.myAtoms = &(info[k].atoms[atomOffset]);
 
       if (molInfo.nBonds > 0) 
-        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
+        molInfo.myBonds = new Bond*[molInfo.nBonds];
       else 
         molInfo.myBonds = NULL;
 
       if (molInfo.nBends > 0) 
-        molInfo.myBends = new (Bend *) [molInfo.nBends];
+        molInfo.myBends = new Bend*[molInfo.nBends];
       else 
         molInfo.myBends = NULL;
 
       if (molInfo.nTorsions > 0) 
-        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
+        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
       else 
         molInfo.myTorsions = NULL;
 
@@ -485,8 +488,11 @@ void SimSetup::makeMolecules(void){
       
 
       //create cutoff group for molecule
+
+      cutoffAtomSet.clear();
       molInfo.myCutoffGroups.clear();
-      for (j = 0; j < molInfo.nCutoffGroups; j++){
+      
+      for (j = 0; j < nCutoffGroups; j++){
 
         currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
         nMembers = currentCutoffGroup->getNMembers(); 
@@ -500,15 +506,30 @@ void SimSetup::makeMolecules(void){
 
           // tempI is atom numbering on local processor
           tempI = molI + atomOffset;
-
+          
           myCutoffGroup->addAtom(info[k].atoms[tempI]);          
+
+          cutoffAtomSet.insert(tempI);
         }
 
         molInfo.myCutoffGroups.push_back(myCutoffGroup);
       }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
-      
 
+      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
 
+      for(j = 0; j < molInfo.nAtoms; j++){
+
+        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
+          myCutoffGroup = new CutoffGroup();
+          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
+          molInfo.myCutoffGroups.push_back(myCutoffGroup);
+        }
+          
+      }
+
+              
+
+
       // After this is all set up, scan through the atoms to 
       // see if they can be added to the integrableObjects:
 
@@ -927,7 +948,37 @@ void SimSetup::gatherInfo(void){
 
     info[i].useInitXSstate = globals->getUseInitXSstate();
     info[i].orthoTolerance = globals->getOrthoBoxTolerance();
-    
+
+    // check for thermodynamic integration
+    if (globals->getUseThermInt()) {
+      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
+	info[i].useThermInt = globals->getUseThermInt();
+	info[i].thermIntLambda = globals->getThermIntLambda();
+	info[i].thermIntK = globals->getThermIntK();
+	
+	Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
+	info[i].restraint = myRestraint;
+      }
+      else {
+	sprintf(painCave.errMsg,
+		"SimSetup Error:\n"
+		"\tKeyword useThermInt was set to 'true' but\n"
+		"\tthermodynamicIntegrationLambda (and/or\n"
+		"\tthermodynamicIntegrationK) was not specified.\n"
+		"\tPlease provide a lambda value and k value in your .bass file.\n");
+	painCave.isFatal = 1;
+	simError();    
+      }
+    }
+    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
+	sprintf(painCave.errMsg,
+		"SimSetup Warning: If you want to use Thermodynamic\n"
+		"\tIntegration, set useThermInt to 'true' in your .bass file.\n"
+		"\tThe useThermInt keyword is 'false' by default, so your\n"
+		"\tlambda and/or k values are being ignored.\n");
+	painCave.isFatal = 0;
+	simError();   
+    }
   }
   
   //setup seed for random number generator
@@ -1240,6 +1291,28 @@ void SimSetup::makeOutNames(void){
         }
       }
 
+      strcpy(info[k].rawPotName, inFileName);
+      nameLength = strlen(info[k].rawPotName);
+      endTest = &(info[k].rawPotName[nameLength - 5]);
+      if (!strcmp(endTest, ".bass")){
+        strcpy(endTest, ".raw");
+      }
+      else if (!strcmp(endTest, ".BASS")){
+        strcpy(endTest, ".raw");
+      }
+      else{
+        endTest = &(info[k].rawPotName[nameLength - 4]);
+        if (!strcmp(endTest, ".bss")){
+          strcpy(endTest, ".raw");
+        }
+        else if (!strcmp(endTest, ".mdl")){
+          strcpy(endTest, ".raw");
+        }
+        else{
+          strcat(info[k].rawPotName, ".raw");
+        }
+      }
+
 #ifdef IS_MPI
 
     }
@@ -1377,21 +1450,33 @@ void SimSetup::calcSysValues(void){
 }
 
 void SimSetup::calcSysValues(void){
-  int i;
+  int i, j;
+  int ncutgroups, atomsingroups, ngroupsinstamp;
 
   int* molMembershipArray;
+  CutoffGroupStamp* cg;
 
   tot_atoms = 0;
   tot_bonds = 0;
   tot_bends = 0;
   tot_torsions = 0;
   tot_rigid = 0;
+  tot_groups = 0;
   for (i = 0; i < n_components; i++){
     tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
     tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
     tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
     tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
     tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
+
+    ncutgroups = comp_stamps[i]->getNCutoffGroups();
+    atomsingroups = 0;
+    for (j=0; j < ncutgroups; j++) {
+      cg = comp_stamps[i]->getCutoffGroup(j);
+      atomsingroups += cg->getNMembers();
+    }
+    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
+    tot_groups += components_nmol[i] * ngroupsinstamp;    
   }
   
   tot_SRI = tot_bonds + tot_bends + tot_torsions;
@@ -1404,7 +1489,7 @@ void SimSetup::calcSysValues(void){
     info[i].n_torsions = tot_torsions;
     info[i].n_SRI = tot_SRI;
     info[i].n_mol = tot_nmol;
-
+    info[i].ngroup = tot_groups;
     info[i].molMembershipArray = molMembershipArray;
   }
 }
@@ -1458,7 +1543,7 @@ void SimSetup::mpiMolDivide(void){
   }
   local_SRI = local_bonds + local_bends + local_torsions;
 
-  info[0].n_atoms = mpiSim->getMyNlocal();  
+  info[0].n_atoms = mpiSim->getNAtomsLocal();  
   
 
   if (local_atoms != info[0].n_atoms){
@@ -1506,7 +1591,7 @@ void SimSetup::makeSysArrays(void){
 
 
     molIndex = 0;
-    for (i = 0; i < mpiSim->getTotNmol(); i++){
+    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
       if (mol2proc[i] == worldRank){
         the_molecules[molIndex].setStampID(molCompType[i]);
         the_molecules[molIndex].setMyIndex(molIndex);