--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/05/27 18:59:17	1203
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/06/01 15:57:30	1211
@@ -10,8 +10,9 @@
 #include "Integrator.hpp"
 #include "simError.h"
 #include "RigidBody.hpp"
-//#include "ConjugateMinimizer.hpp"
 #include "OOPSEMinimizer.hpp"
+//#include "ConstraintElement.hpp"
+//#include "ConstraintPair.hpp"
 
 #ifdef IS_MPI
 #include "mpiBASS.h"
@@ -198,6 +199,12 @@ void SimSetup::makeMolecules(void){
   char* molName;
   char rbName[100];
 
+  //ConstraintPair* consPair; //constraint pair
+  //ConstraintElement* consElement1;  //first element of constraint pair
+  //ConstraintElement* consElement2;  //second element of constraint pair
+  //int whichRigidBody; 
+  //int consAtomIndex;  //index of constraint atom in rigid body's atom array
+  //vector<pair<int, int> > jointAtoms;
   //init the forceField paramters
 
   the_ff->readParams();
@@ -210,6 +217,7 @@ void SimSetup::makeMolecules(void){
     the_ff->setSimInfo(&(info[k]));
 
     atomOffset = 0;
+    groupOffset = 0;
 
     for (i = 0; i < info[k].n_mol; i++){
       stampID = info[k].molecules[i].getStampID();
@@ -498,6 +506,7 @@ void SimSetup::makeMolecules(void){
         nMembers = currentCutoffGroup->getNMembers(); 
 
         myCutoffGroup = new CutoffGroup();
+        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
         
         for (int cg = 0; cg < nMembers; cg++) {
 
@@ -506,13 +515,23 @@ void SimSetup::makeMolecules(void){
 
           // tempI is atom numbering on local processor
           tempI = molI + atomOffset;
-          
+
+#ifdef IS_MPI
+          globalID = info[k].atoms[tempI]->getGlobalIndex()
+#else 
+          globalID = info[k].atoms[tempI]->getIndex();
+#endif
+
+          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; 
+
           myCutoffGroup->addAtom(info[k].atoms[tempI]);          
 
           cutoffAtomSet.insert(tempI);
         }
-
+       
         molInfo.myCutoffGroups.push_back(myCutoffGroup);
+        groupOffset++;
+
       }//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
@@ -522,14 +541,19 @@ void SimSetup::makeMolecules(void){
         if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
           myCutoffGroup = new CutoffGroup();
           myCutoffGroup->addAtom(molInfo.myAtoms[j]);
+          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
+#ifdef IS_MPI
+          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex()
+#else 
+          globalID = info[k].atoms[atomOffset + j]->getIndex();
+#endif
+          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; 
           molInfo.myCutoffGroups.push_back(myCutoffGroup);
+          groupOffset++;
         }
           
       }
 
-              
-
-
       // After this is all set up, scan through the atoms to 
       // see if they can be added to the integrableObjects:
 
@@ -560,8 +584,54 @@ void SimSetup::makeMolecules(void){
         info[k].integrableObjects.push_back(mySD);      
         molInfo.myIntegrableObjects.push_back(mySD);
       }
-    
+
+
+    /*
+
+      //creat ConstraintPair.
+      molInfo.myConstraintPair.clear();
       
+      for (j = 0; j < molInfo.nBonds; j++){
+
+        //if both atoms are in the same rigid body, just skip it
+        currentBond = comp_stamps[stampID]->getBond(j);
+        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
+
+          tempI = currentBond->getA() + atomOffset;
+          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
+            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
+          else
+             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);       
+
+          tempJ =  currentBond->getB() + atomOffset;
+          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
+            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
+          else
+             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);     
+
+          consPair = new DistanceConstraintPair(consElement1, consElement2);
+          molInfo.myConstraintPairs.push_back(consPair);
+        }
+      }  
+      
+      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair
+      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
+        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
+          
+          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
+
+          for(size_t m = 0; m < jointAtoms.size(); m++){          
+            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
+            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
+
+            consPair = new JointConstraintPair(consElement1, consElement2);  
+            molInfo.myConstraintPairs.push_back(consPair);            
+          }
+
+        }
+      }
+      
+*/      
       // send the arrays off to the forceField for init.
       
       the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
@@ -1500,8 +1570,10 @@ void SimSetup::mpiMolDivide(void){
   int i, j, k;
   int localMol, allMol;
   int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
-  int local_rigid;
+  int local_rigid, local_groups;
   vector<int> globalMolIndex;
+  int ncutgroups, atomsingroups, ngroupsinstamp;
+  CutoffGroupStamp* cg;
 
   mpiSim = new mpiSimulation(info);
 
@@ -1521,6 +1593,7 @@ void SimSetup::mpiMolDivide(void){
   local_bends = 0;
   local_torsions = 0;
   local_rigid = 0;
+  local_groups = 0;
   globalAtomCounter = 0;
 
   for (i = 0; i < n_components; i++){
@@ -1531,6 +1604,17 @@ void SimSetup::mpiMolDivide(void){
         local_bends += comp_stamps[i]->getNBends();
         local_torsions += comp_stamps[i]->getNTorsions();
         local_rigid += comp_stamps[i]->getNRigidBodies();
+
+        ncutgroups = comp_stamps[i]->getNCutoffGroups();
+        atomsingroups = 0;
+        for (k=0; k < ncutgroups; k++) {
+          cg = comp_stamps[i]->getCutoffGroup(k);
+          atomsingroups += cg->getNMembers();
+        }
+        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + 
+          ncutgroups;
+        local_groups += ngroupsinstamp;    
+
         localMol++;
       }      
       for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
@@ -1545,7 +1629,6 @@ void SimSetup::mpiMolDivide(void){
 
   info[0].n_atoms = mpiSim->getNAtomsLocal();  
   
-
   if (local_atoms != info[0].n_atoms){
     sprintf(painCave.errMsg,
             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
@@ -1555,6 +1638,16 @@ void SimSetup::mpiMolDivide(void){
     simError();
   }
 
+  info[0].ngroup = mpiSim->getNGroupsLocal();   
+  if (local_groups != info[0].ngroup){
+    sprintf(painCave.errMsg,
+            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
+            "\tlocalGroups (%d) are not equal.\n",
+            info[0].ngroup, local_groups);
+    painCave.isFatal = 1;
+    simError();
+  }
+  
   info[0].n_bonds = local_bonds;
   info[0].n_bends = local_bends;
   info[0].n_torsions = local_torsions;