--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/26 20:22:02	407
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/27 17:32:03	420
@@ -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,388 +768,150 @@ void SimSetup::makeMolecules( void ){
 	exJ = tempEx;
       }
 #ifdef IS_MPI
-      
-      the_excludes[index*2] =     
-	the_atoms[exI]->getGlobalIndex() + 1;
-      the_excludes[index*2 + 1] = 
-	the_atoms[exJ]->getGlobalIndex() + 1;
+      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[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
-    
-  }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-void SimSetup::makeAtoms( void ){
-
-  int i, j, k, index;
-  double ux, uy, uz, uSqr, u;
-  AtomStamp* current_atom;
-
-  DirectionalAtom* dAtom;
-  int molIndex, molStart, molEnd, nMemb, lMolIndex;
+    }
+    excludeOffset += info.nBonds;
 
-  lMolIndex = 0;
-  molIndex = 0;
-  index = 0;
-  for( i=0; i<n_components; i++ ){
-
-    for( j=0; j<components_nmol[i]; j++ ){
-
-#ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
-
-	molStart = index;
-	nMemb = comp_stamps[i]->getNAtoms();
-	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
+    //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;
 	  
-	  current_atom = comp_stamps[i]->getAtom( k );
-	  if( current_atom->haveOrientation() ){
+      if( currentBend->haveExtras() ){
 	    
-	    dAtom = new DirectionalAtom(index);
-	    simnfo->n_oriented++;
-	    the_atoms[index] = dAtom;
+	extras = current_bend->getExtras();
+	current_extra = extras;
 	    
-	    ux = current_atom->getOrntX();
-	    uy = current_atom->getOrntY();
-	    uz = current_atom->getOrntZ();
-	    
-	    uSqr = (ux * ux) + (uy * uy) + (uz * uz);
-	    
-	    u = sqrt( uSqr );
-	    ux = ux / u;
-	    uy = uy / u;
-	    uz = uz / u;
-	    
-	    dAtom->setSUx( ux );
-	    dAtom->setSUy( uy );
-	    dAtom->setSUz( uz );
-	  }
-	  else{
-	    the_atoms[index] = new GeneralAtom(index);
-	  }
-	  the_atoms[index]->setType( current_atom->getType() );
-	  the_atoms[index]->setIndex( index );
-	  
-	  // increment the index and repeat;
-	  index++;
-	}
-	
-	molEnd = index -1;
-	the_molecules[lMolIndex].setNMembers( nMemb );
-	the_molecules[lMolIndex].setStartAtom( molStart );
-	the_molecules[lMolIndex].setEndAtom( molEnd );
-	the_molecules[lMolIndex].setStampID( i );
-	lMolIndex++;
-
-#ifdef IS_MPI
-      }
-#endif //is_mpi
-      
-      molIndex++;
-    }
-  }
-
-#ifdef IS_MPI
-    for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
-    
-    delete[] globalIndex;
-
-    mpiSim->mpiRefresh();
-#endif //IS_MPI
-	  
-  the_ff->initializeAtoms();
-}
-
-void SimSetup::makeBonds( void ){
-
-  int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
-  bond_pair* the_bonds;
-  BondStamp* current_bond;
-
-  the_bonds = new bond_pair[tot_bonds];
-  index = 0;
-  offset = 0;
-  molIndex = 0;
-
-  for( i=0; i<n_components; i++ ){
-
-    for( j=0; j<components_nmol[i]; j++ ){
-
-#ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
-	
-	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
-	  
-	  current_bond = comp_stamps[i]->getBond( k );
-	  the_bonds[index].a = current_bond->getA() + offset;
-	  the_bonds[index].b = current_bond->getB() + offset;
-
-	  exI = the_bonds[index].a;
-	  exJ = the_bonds[index].b;
-
-	  // exclude_I must always be the smaller of the pair
-	  if( exI > exJ ){
-	    tempEx = exI;
-	    exI = exJ;
-	    exJ = tempEx;
-	  }
-
-	  
-#ifdef IS_MPI
-
-	  the_excludes[index*2] =     
-	    the_atoms[exI]->getGlobalIndex() + 1;
-	  the_excludes[index*2 + 1] = 
-	    the_atoms[exJ]->getGlobalIndex() + 1;
-
-#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)
-#endif  //is_mpi
-	  
-	  // increment the index and repeat;
-	  index++;
-	}
-	offset += comp_stamps[i]->getNAtoms();
-	
-#ifdef IS_MPI
-      }
-#endif //is_mpi
-      
-      molIndex++;
-    }      
-  }
-
-  the_ff->initializeBonds( the_bonds );
-}
-
-void SimSetup::makeBends( void ){
-
-  int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
-  bend_set* the_bends;
-  BendStamp* current_bend;
-  LinkedAssign* extras;
-  LinkedAssign* current_extra;
-  
-
-  the_bends = new bend_set[tot_bends];
-  index = 0;
-  offset = 0;
-  molIndex = 0;
-  for( i=0; i<n_components; i++ ){
-
-    for( j=0; j<components_nmol[i]; j++ ){
-
-#ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
-
-	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
-	  
-	  current_bend = comp_stamps[i]->getBend( k );
-	  the_bends[index].a = current_bend->getA() + offset;
-	  the_bends[index].b = current_bend->getB() + offset;
-	  the_bends[index].c = current_bend->getC() + offset;
-	  
-	  if( current_bend->haveExtras() ){
-	    
-	    extras = current_bend->getExtras();
-	    current_extra = extras;
-	    
-	    while( current_extra != NULL ){
-	      if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
+	while( current_extra != NULL ){
+	  if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
 		
-		switch( current_extra->getType() ){
+	    switch( current_extra->getType() ){
+	      
+	    case 0:
+	      theBends[j].ghost =
+		current_extra->getInt() + atomOffset;
+	      theBends[j].isGhost = 1;
+	      break;
 		  
-		case 0:
-		  the_bends[index].ghost =
-		    current_extra->getInt() + offset;
-		  the_bends[index].isGhost = 1;
-		  break;
-		  
-		case 1:
-		  the_bends[index].ghost = 
-		    (int)current_extra->getDouble() + offset;
-		  the_bends[index].isGhost = 1;
-		  break;
-		  
-		default:
-		  sprintf( painCave.errMsg,
-			   "SimSetup Error: ghostVectorSource was neiter a "
-			   "double nor an int.\n"
-			   "-->Bend[%d] in %s\n",
-			   k, comp_stamps[i]->getID() );
-		  painCave.isFatal = 1;
-		  simError();
-		}
-	      }
+	    case 1:
+	      theBends[j].ghost = 
+		(int)current_extra->getDouble() + atomOffset;
+	      theBends[j].isGhost = 1;
+	      break;
 	      
-	      else{
-		
-		sprintf( painCave.errMsg,
-			 "SimSetup Error: unhandled bend assignment:\n"
-			 "    -->%s in Bend[%d] in %s\n",
-			 current_extra->getlhs(),
-			 k, comp_stamps[i]->getID() );
-		painCave.isFatal = 1;
-		simError();
-	      }
-	      
-	      current_extra = current_extra->getNext();
+	    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();
 	    }
 	  }
 	  
-	  if( !the_bends[index].isGhost ){
-	    
-	    exI = the_bends[index].a;
-	    exJ = the_bends[index].c;
-	  }
 	  else{
 	    
-	    exI = the_bends[index].a;
-	    exJ = the_bends[index].b;
+	    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();
 	  }
 	  
-	  // exclude_I must always be the smaller of the pair
-	  if( exI > exJ ){
-	    tempEx = exI;
-	    exI = exJ;
-	    exJ = tempEx;
-	  }
-
-
-#ifdef IS_MPI
-
-	  the_excludes[(index + tot_bonds)*2] =     
-	    the_atoms[exI]->getGlobalIndex() + 1;
-	  the_excludes[(index + tot_bonds)*2 + 1] = 
-	    the_atoms[exJ]->getGlobalIndex() + 1;
-	  
-#else  // isn't MPI
-	  
-	  the_excludes[(index + tot_bonds)*2] =     exI + 1; 
-	  the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1;
-	  // fortran index from 1 (hence the +1 in the indexing)
-#endif  //is_mpi
-	  
-	  
-	  // increment the index and repeat;
-	  index++;
+	  current_extra = current_extra->getNext();
 	}
-	offset += comp_stamps[i]->getNAtoms();
+      }
+	  
+      if( !theBends[j].isGhost ){
+	    
+	exI = theBends[j].a;
+	exJ = theBends[j].c;
+      }
+      else{
 	
-#ifdef IS_MPI
+	exI = theBends[j].a;
+	exJ = theBends[j].b;
       }
-#endif //is_mpi
-
-      molIndex++;
-    }
-  }
+      
+      // 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;
 
-#ifdef IS_MPI
-  sprintf( checkPointMsg,
-	   "Successfully created the bends list.\n" );
-  MPIcheckPoint();
-#endif // is_mpi
-  
+    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;
 
-  the_ff->initializeBends( the_bends );
-}
-
-void SimSetup::makeTorsions( void ){
-
-  int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
-  torsion_set* the_torsions;
-  TorsionStamp* current_torsion;
-
-  the_torsions = new torsion_set[tot_torsions];
-  index = 0;
-  offset = 0;
-  molIndex = 0;
-  for( i=0; i<n_components; i++ ){
-
-    for( j=0; j<components_nmol[i]; j++ ){
-
+      // exclude_I must always be the smaller of the pair
+      if( exI > exJ ){
+	tempEx = exI;
+	exI = exJ;
+	exJ = tempEx;
+      }
 #ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
-
-      for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
-
-	current_torsion = comp_stamps[i]->getTorsion( k );
-	the_torsions[index].a = current_torsion->getA() + offset;
-	the_torsions[index].b = current_torsion->getB() + offset;
-	the_torsions[index].c = current_torsion->getC() + offset;
-	the_torsions[index].d = current_torsion->getD() + offset;
-
-	exI = the_torsions[index].a;
-	exJ = the_torsions[index].d;
-
-	
-	// exclude_I must always be the smaller of the pair
-	if( exI > exJ ){
-	  tempEx = exI;
-	  exI = exJ;
-	  exJ = tempEx;
-	}
-
-
-#ifdef IS_MPI
-	
-	the_excludes[(index + tot_bonds + tot_bends)*2] =     
-	  the_atoms[exI]->getGlobalIndex() + 1;
-	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = 
-	  the_atoms[exJ]->getGlobalIndex() + 1;
-	
+      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[(index + tot_bonds + tot_bends)*2] =     exI + 1; 
-	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1;
-	// fortran indexes from 1 (hence the +1 in the indexing)
+      the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
-	
+    }
+    excludeOffset += info.nTorsions;
 
-	// increment the index and repeat;
-	index++;
-      }
-      offset += comp_stamps[i]->getNAtoms();
+    
+    // send the arrays off to the forceField for init.
 
-#ifdef IS_MPI
-      }
-#endif //is_mpi      
+    the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
+    the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
+    the_ff->initializeBends( info.nBends, info.myBends, theBends );
+    the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
 
-      molIndex++;
-    }
+
+    the_molecules[i].initialize( info );
+    atomOffset += info.nAtoms;
   }
 
-  the_ff->initializeTorsions( the_torsions );
+  // clean up the forcefield
+  the_ff->calcRcut();
+  the_ff->cleanMe();
 }
 
 void SimSetup::initFromBass( void ){