--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/26 21:24:08	411
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/26 21:50:33	412
@@ -307,20 +307,20 @@ void SimSetup::createSim( void ){
 
 
   if( simnfo->n_SRI ){
-    the_sris = new SRI*[simnfo->n_SRI];
-    the_excludes = new int[2 * simnfo->n_SRI];
+    Exclude::createArray(simnfo->n_SRI);
+    the_excludes = new Exclude*[simnfo->n_SRI];
     simnfo->globalExcludes = new int;
     simnfo->n_exclude = tot_SRI;
   }
   else{
     
-    the_excludes = new int[2];
-    the_excludes[0] = 0;
-    the_excludes[1] = 0;
+    Exclude::createArray( 1 );
+    the_excludes = new Exclude*;
+    the_excludes[0] = new Exclude(0);
+    the_excludes[0]->setPair( 0,0 );
     simnfo->globalExcludes = new int;
     simnfo->globalExcludes[0] = 0;
-
-    simnfo->n_exclude = 1;
+    simnfo->n_exclude = 0;
   }
 
   // set the arrays into the SimInfo object
@@ -674,9 +674,11 @@ void SimSetup::makeMolecules( void ){
 
 void SimSetup::makeMolecules( void ){
 
-  int i, j, exI, exJ, tempEx, stampID, atomOffset;
+  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
   molInit info;
   DirectionalAtom* dAtom;
+  LinkedAssign* extras;
+  LinkedAssign* current_extra;
   AtomStamp* currentAtom;
   BondStamp* currentBond;
   BendStamp* currentBend;
@@ -690,6 +692,7 @@ void SimSetup::makeMolecules( void ){
   // init the molecules
 
   atomOffset = 0;
+  excludeOffset = 0;
   for(i=0; i<simnfo->n_mol; i++){
     
     stampID = the_molecules[i].getStampID();
@@ -698,8 +701,10 @@ void SimSetup::makeMolecules( void ){
     info.nBonds    = comp_stamps[stampID]->getNBonds();
     info.nBends    = comp_stamps[stampID]->getNBends();
     info.nTorsions = comp_stamps[stampID]->getNTorsions();
-    
+    info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
+
     info.myAtoms = &the_atoms[atomOffset];
+    info.myExcludes = &the_excludes[excludeOffset];
     info.myBonds = new Bond*[info.nBonds];
     info.myBends = new Bend*[info.nBends];
     info.myTorsions = new Torsions*[info.nTorsions];
@@ -747,7 +752,7 @@ void SimSetup::makeMolecules( void ){
     } 
     
     // make the bonds
-    for(j=0; j<nBonds; j++){
+    for(j=0; j<info.nBonds; j++){
       
       currentBond = comp_stamps[stampID]->getBond( j );
       theBonds[j].a = currentBond->getA() + atomOffset;
@@ -763,30 +768,135 @@ void SimSetup::makeMolecules( void ){
 	exJ = tempEx;
       }
 #ifdef IS_MPI
+      tempEx = exI;
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
+      tempEx = exJ;
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
       
-      the_excludes[index*2] =     
-	the_atoms[exI]->getGlobalIndex() + 1;
-      the_excludes[index*2 + 1] = 
-	the_atoms[exJ]->getGlobalIndex() + 1;
-      
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
 #else  // isn't MPI
-      
-      the_excludes[index*2] =     exI + 1; 
-      the_excludes[index*2 + 1] = exJ + 1;
-      // fortran index from 1 (hence the +1 in the indexing)
-
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
-    
-  }
+    }
+    excludeOffset += info.nBonds;
 
-
-
-
-
-
+    //make the bends
+    for(j=0; j<info.nBends; j++){
+      
+      currentBend = comp_stamps[stampID]->getBend( j );
+      theBends[j].a = currentBend->getA() + atomOffset;
+      theBends[j].b = currentBend->getB() + atomOffset;
+      theBends[j].c = currentBend->getC() + atomOffset;
+	  
+      if( currentBend->haveExtras() ){
+	    
+	extras = current_bend->getExtras();
+	current_extra = extras;
+	    
+	while( current_extra != NULL ){
+	  if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
+		
+	    switch( current_extra->getType() ){
+	      
+	    case 0:
+	      theBends[j].ghost =
+		current_extra->getInt() + atomOffset;
+	      theBends[j].isGhost = 1;
+	      break;
+		  
+	    case 1:
+	      theBends[j].ghost = 
+		(int)current_extra->getDouble() + atomOffset;
+	      theBends[j].isGhost = 1;
+	      break;
+	      
+	    default:
+	      sprintf( painCave.errMsg,
+		       "SimSetup Error: ghostVectorSource was neiter a "
+		       "double nor an int.\n"
+		       "-->Bend[%d] in %s\n",
+		       j, comp_stamps[stampID]->getID() );
+	      painCave.isFatal = 1;
+	      simError();
+	    }
+	  }
+	  
+	  else{
+	    
+	    sprintf( painCave.errMsg,
+		     "SimSetup Error: unhandled bend assignment:\n"
+		     "    -->%s in Bend[%d] in %s\n",
+		     current_extra->getlhs(),
+		     j, comp_stamps[stampID]->getID() );
+	    painCave.isFatal = 1;
+	    simError();
+	  }
+	  
+	  current_extra = current_extra->getNext();
+	}
+      }
+	  
+      if( !theBends[j].isGhost ){
+	    
+	exI = theBends[j].a;
+	exJ = theBends[j].c;
+      }
+      else{
+	
+	exI = theBends[j].a;
+	exJ = theBends[j].b;
+      }
+      
+      // exclude_I must always be the smaller of the pair
+      if( exI > exJ ){
+	tempEx = exI;
+	exI = exJ;
+	exJ = tempEx;
+      }
+#ifdef IS_MPI
+      tempEx = exI;
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
+      tempEx = exJ;
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
+      
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
+#else  // isn't MPI
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
+#endif  //is_mpi
+    }
+    excludeOffset += info.nBends;
 
+    for(j=0; j<info.nTorsions; j++){
+      
+      currentTorsion = comp_stamps[stampID]->getTorsion( j );
+      theTorsions[j].a = currentTorsion->getA() + atomOffset;
+      theTorsions[j].b = currentTorsion->getB() + atomOffset;
+      theTorsions[j].c = currentTorsion->getC() + atomOffset;
+      theTorsions[j].d = currentTorsion->getD() + atomOffset;
+      
+      exI = theTorsions[j].a;
+      exJ = theTorsions[j].d;
 
+      // exclude_I must always be the smaller of the pair
+      if( exI > exJ ){
+	tempEx = exI;
+	exI = exJ;
+	exJ = tempEx;
+      }
+#ifdef IS_MPI
+      tempEx = exI;
+      exI = the_atoms[tempEx]->getGlobalIndex() + 1;
+      tempEx = exJ;
+      exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
+      
+      the_excludes[j+excludeOffset]->setPair( exI, exJ );
+#else  // isn't MPI
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
+#endif  //is_mpi
+    }
+    excludeOffset += info.nTorsions;
 
+