--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/24 21:55:34	394
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/06/19 19:11:53	557
@@ -12,6 +12,16 @@
 #include "mpiSimulation.hpp"
 #endif
 
+// some defines for ensemble and Forcefield  cases
+
+#define NVE_ENS 0
+#define NVT_ENS 1
+#define NPT_ENS 2
+
+#define FF_DUFF 0
+#define FF_LJ   1
+
+
 SimSetup::SimSetup(){
   stamps = new MakeStamps();
   globals = new Globals();
@@ -66,7 +76,13 @@ void SimSetup::createSim( void ){
 
   MakeStamps *the_stamps;
   Globals* the_globals;
-  int i, j;
+  int i, j, k, globalAtomIndex;
+  
+  int ensembleCase;
+  int ffCase;
+  
+  ensembleCase = -1;
+  ffCase = -1;
 
   // get the stamps and globals;
   the_stamps = stamps;
@@ -80,22 +96,109 @@ void SimSetup::createSim( void ){
   // get the ones we know are there, yet still may need some work.
   n_components = the_globals->getNComponents();
   strcpy( force_field, the_globals->getForceField() );
+
+  if( !strcasecmp( force_field, "DUFF" ))    ffCase = FF_DUFF;
+  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
+  else{
+    sprintf( painCave.errMsg,
+	     "SimSetup Error. Unrecognized force field -> %s\n",
+	     force_field );
+    painCave.isFatal = 1;
+    simError();
+  }
+
+  // get the ensemble:
   strcpy( ensemble, the_globals->getEnsemble() );
+
+  if( !strcasecmp( ensemble, "NVE" ))      ffCase = NVE_ENS;
+  else if( !strcasecmp( ensemble, "NVT" )) ffCase = NVT_ENS;
+  else if( !strcasecmp( ensemble, "NPT" )) ffCase = NPT_ENS;
+  else{
+    sprintf( painCave.errMsg,
+	     "SimSetup Warning. Unrecognized Ensemble -> %s, "
+             "reverting to NVE for this simulation.\n",
+	     ensemble );
+    painCave.isFatal = 0;
+    simError();
+    strcpy( ensemble, "NVE" );
+    ensembleCase = NVE_ENS;
+  }  
   strcpy( simnfo->ensemble, ensemble );
 
+
+//   if( !strcasecmp( ensemble, "NPT" ) ) {
+//     the_extendedsystem = new ExtendedSystem( simnfo );
+//     the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
+//     if (the_globals->haveTargetPressure())
+//       the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
+//     else {
+//       sprintf( painCave.errMsg,
+//                "SimSetup error: If you use the constant pressure\n"
+//                "    ensemble, you must set targetPressure.\n"
+//                "    This was found in the BASS file.\n");
+//       painCave.isFatal = 1;
+//       simError();
+//     }
+
+//     if (the_globals->haveTauThermostat()) 
+//       the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
+//     else if (the_globals->haveQmass())
+//       the_extendedsystem->setQmass(the_globals->getQmass());
+//     else {
+//       sprintf( painCave.errMsg,
+//                "SimSetup error: If you use one of the constant temperature\n"
+//                "    ensembles, you must set either tauThermostat or qMass.\n"
+//                "    Neither of these was found in the BASS file.\n");
+//       painCave.isFatal = 1;
+//       simError();
+//     }
+
+//     if (the_globals->haveTauBarostat())
+//       the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
+//     else {
+//       sprintf( painCave.errMsg,
+//                "SimSetup error: If you use the constant pressure\n"
+//                "    ensemble, you must set tauBarostat.\n"
+//                "    This was found in the BASS file.\n");
+//       painCave.isFatal = 1;
+//       simError();
+//     }
+
+//   } else if ( !strcasecmp( ensemble, "NVT") ) {
+//     the_extendedsystem = new ExtendedSystem( simnfo );
+//     the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
+
+//     if (the_globals->haveTauThermostat()) 
+//       the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
+//     else if (the_globals->haveQmass())
+//       the_extendedsystem->setQmass(the_globals->getQmass());
+//     else {
+//       sprintf( painCave.errMsg,
+//                "SimSetup error: If you use one of the constant temperature\n"
+//                "    ensembles, you must set either tauThermostat or qMass.\n"
+//                "    Neither of these was found in the BASS file.\n");
+//       painCave.isFatal = 1;
+//       simError();
+//     }
+
   strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
   simnfo->usePBC = the_globals->getPBC();
 	  
+  int usesDipoles = 0;
+  switch( ffCase ){
 
-
-  if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
-  else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
-  else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
-  else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
-  else{
+  case FF_DUFF:
+    the_ff = new DUFF();
+    usesDipoles = 1;
+    break;
+
+  case FF_LJ:
+    the_ff = LJFF();
+    break;
+
+  default:
     sprintf( painCave.errMsg,
-	     "SimSetup Error. Unrecognized force field -> %s\n",
-	     force_field );
+	     "SimSetup Error. Unrecognized force field in case statement.\n");
     painCave.isFatal = 1;
     simError();
   }
@@ -105,8 +208,6 @@ void SimSetup::createSim( void ){
   MPIcheckPoint();
 #endif // is_mpi
 
-  
-
   // get the components and calculate the tot_nMol and indvidual n_mol
   the_components = the_globals->getComponents();
   components_nmol = new int[n_components];
@@ -230,24 +331,24 @@ void SimSetup::createSim( void ){
   simnfo->n_torsions = tot_torsions;
   simnfo->n_SRI = tot_SRI;
   simnfo->n_mol = tot_nmol;
-
   
+  simnfo->molMembershipArray = new int[tot_atoms];
+
 #ifdef IS_MPI
 
   // divide the molecules among processors here.
   
   mpiSim = new mpiSimulation( simnfo );
   
-  
-
   globalIndex = mpiSim->divideLabor();
 
-
-
   // set up the local variables 
   
   int localMol, allMol;
   int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
+
+  int* mol2proc = mpiSim->getMolToProcMap();
+  int* molCompType = mpiSim->getMolComponentType();
   
   allMol = 0;
   localMol = 0;
@@ -255,12 +356,14 @@ void SimSetup::createSim( void ){
   local_bonds = 0;
   local_bends = 0;
   local_torsions = 0;
+  globalAtomIndex = 0;
+
+
   for( i=0; i<n_components; i++ ){
 
     for( j=0; j<components_nmol[i]; j++ ){
       
-      if( mpiSim->getMyMolStart() <= allMol &&
-	  allMol <= mpiSim->getMyMolEnd() ){
+      if( mol2proc[allMol] == worldRank ){
 	
 	local_atoms +=    comp_stamps[i]->getNAtoms();
 	local_bonds +=    comp_stamps[i]->getNBonds();
@@ -268,18 +371,22 @@ void SimSetup::createSim( void ){
 	local_torsions += comp_stamps[i]->getNTorsions();
 	localMol++;
       }      
-      allMol++;
+      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
+        simnfo->molMembershipArray[globalAtomIndex] = allMol;
+        globalAtomIndex++;
+      }
+
+      allMol++;      
     }
   }
   local_SRI = local_bonds + local_bends + local_torsions;
   
-
   simnfo->n_atoms = mpiSim->getMyNlocal();  
   
   if( local_atoms != simnfo->n_atoms ){
     sprintf( painCave.errMsg,
 	     "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
-	     " localAtom (%d) are note equal.\n",
+	     " localAtom (%d) are not equal.\n",
 	     simnfo->n_atoms,
 	     local_atoms );
     painCave.isFatal = 1;
@@ -304,29 +411,68 @@ void SimSetup::createSim( void ){
   Atom::createArrays(simnfo->n_atoms);
   the_atoms = new Atom*[simnfo->n_atoms];
   the_molecules = new Molecule[simnfo->n_mol];
+  int molIndex;
 
+  // initialize the molecule's stampID's
 
-  if( simnfo->n_SRI ){
-    the_sris = new SRI*[simnfo->n_SRI];
-    the_excludes = new int[2 * simnfo->n_SRI];
+#ifdef IS_MPI
+  
+
+  molIndex = 0;
+  for(i=0; i<mpiSim->getTotNmol(); i++){
+    
+    if(mol2proc[i] == worldRank ){
+      the_molecules[molIndex].setStampID( molCompType[i] );
+      the_molecules[molIndex].setMyIndex( molIndex );
+      the_molecules[molIndex].setGlobalIndex( i );
+      molIndex++;
+    }
+  }
+
+#else // is_mpi
+  
+  molIndex = 0;
+  globalAtomIndex = 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++) {
+        simnfo->molMembershipArray[globalAtomIndex] = molIndex;
+        globalAtomIndex++;
+      }
+      molIndex++;
+    }
+  }
+    
+
+#endif // is_mpi
+
+
+  if( simnfo->n_SRI ){
+    
+    Exclude::createArray(simnfo->n_SRI);
+    the_excludes = new Exclude*[simnfo->n_SRI];
+    for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
     simnfo->globalExcludes = new int;
-    simnfo->n_exclude = tot_SRI;
+    simnfo->n_exclude = simnfo->n_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
 
   simnfo->atoms = the_atoms;
-  simnfo->sr_interactions = the_sris;
+  simnfo->molecules = the_molecules;
   simnfo->nGlobalExcludes = 0;
   simnfo->excludes = the_excludes;
 
@@ -383,25 +529,12 @@ void SimSetup::createSim( void ){
 
   the_ff->setSimInfo( simnfo );
 
-  makeAtoms();
+  makeMolecules();
   simnfo->identArray = new int[simnfo->n_atoms];
   for(i=0; i<simnfo->n_atoms; i++){
     simnfo->identArray[i] = the_atoms[i]->getIdent();
   }
   
-  if( tot_bonds ){
-    makeBonds();
-  }
-
-  if( tot_bends ){
-    makeBends();
-  }
-
-  if( tot_torsions ){
-    makeTorsions();
-  }
-
-
   if (the_globals->getUseRF() ) {
     simnfo->useReactionField = 1;
   
@@ -443,11 +576,11 @@ void SimSetup::createSim( void ){
     }
     simnfo->dielectric = the_globals->getDielectric();  
   } else {
-    if (simnfo->n_dipoles) {
+    if (usesDipoles) {
       
       if( !the_globals->haveECR() ){
         sprintf( painCave.errMsg,
-                 "SimSetup Warning: using default value of 1/2 the smallest"
+                 "SimSetup Warning: using default value of 1/2 the smallest "
                  "box length for the electrostaticCutoffRadius.\n"
                  "I hope you have a very fast processor!\n");
         painCave.isFatal = 0;
@@ -463,7 +596,7 @@ void SimSetup::createSim( void ){
       
       if( !the_globals->haveEST() ){
         sprintf( painCave.errMsg,
-                 "SimSetup Warning: using default value of 5% of the"
+                 "SimSetup Warning: using default value of 5%% of the "
                  "electrostaticCutoffRadius for the "
                  "electrostaticSkinThickness\n"
                  );
@@ -634,19 +767,23 @@ void SimSetup::createSim( void ){
   if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
 
 
-//   // make the longe range forces and the integrator
+  // make the integrator
 
-//   new AllLong( simnfo );
 
-  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
-  if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
-  if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
-  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
+  if( !strcmp( ensemble, "TraPPE_Ex" ) ){
+    new Symplectic(simnfo, the_ff, the_extendedsystem);
+  }
+  else if( !strcmp( force_field, "LJ" ) ){
+    new Verlet( *simnfo, the_ff, the_extendedsystem );
+  }
 
-
+#ifdef IS_MPI
+  mpiSim->mpiRefresh();
+#endif
 
   // initialize the Fortran
-  
+
+
   simnfo->refreshSim();
   
   if( !strcmp( simnfo->mixingRule, "standard") ){
@@ -671,358 +808,266 @@ void SimSetup::createSim( void ){
 #endif // is_mpi
 }
 
-void SimSetup::makeAtoms( void ){
 
-  int i, j, k, index;
-  double ux, uy, uz, uSqr, u;
-  AtomStamp* current_atom;
+void SimSetup::makeMolecules( void ){
 
+  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
+  molInit info;
   DirectionalAtom* dAtom;
-  int molIndex, molStart, molEnd, nMemb, lMolIndex;
-
-  lMolIndex = 0;
-  molIndex = 0;
-  index = 0;
-  for( i=0; i<n_components; i++ ){
-
-    for( j=0; j<components_nmol[i]; j++ ){
+  LinkedAssign* extras;
+  LinkedAssign* current_extra;
+  AtomStamp* currentAtom;
+  BondStamp* currentBond;
+  BendStamp* currentBend;
+  TorsionStamp* currentTorsion;
 
-#ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
+  bond_pair* theBonds;
+  bend_set* theBends;
+  torsion_set* theTorsions;
 
-	molStart = index;
-	nMemb = comp_stamps[i]->getNAtoms();
-	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
-	  
-	  current_atom = comp_stamps[i]->getAtom( k );
-	  if( current_atom->haveOrientation() ){
-	    
-	    dAtom = new DirectionalAtom(index);
-	    simnfo->n_oriented++;
-	    the_atoms[index] = dAtom;
-	    
-	    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++;
+  
+  //init the forceField paramters
+
+  the_ff->readParams();
+
+  
+  // init the atoms
 
-#ifdef IS_MPI
-      }
-#endif //is_mpi
-      
-      molIndex++;
-    }
-  }
-
-#ifdef IS_MPI
-    for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
+  double ux, uy, uz, u, uSqr;
+  
+  atomOffset = 0;
+  excludeOffset = 0;
+  for(i=0; i<simnfo->n_mol; i++){
     
-    delete[] globalIndex;
+    stampID = the_molecules[i].getStampID();
 
-    mpiSim->mpiRefresh();
-#endif //IS_MPI
-	  
-  the_ff->initializeAtoms();
-}
+    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
+    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;
 
-void SimSetup::makeBonds( void ){
+    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 Torsion*[info.nTorsions];
 
-  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	
+    theBonds = new bond_pair[info.nBonds];
+    theBends = new bend_set[info.nBends];
+    theTorsions = new torsion_set[info.nTorsions];
+    
+    // make the Atoms
+    
+    for(j=0; j<info.nAtoms; j++){
+      
+      currentAtom = comp_stamps[stampID]->getAtom( j );
+      if( currentAtom->haveOrientation() ){
 	
-	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();
+	dAtom = new DirectionalAtom(j + atomOffset);
+	simnfo->n_oriented++;
+	info.myAtoms[j] = dAtom;
 	
-#ifdef IS_MPI
+	ux = currentAtom->getOrntX();
+	uy = currentAtom->getOrntY();
+	uz = currentAtom->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 );
       }
-#endif //is_mpi
+      else{
+	info.myAtoms[j] = new GeneralAtom(j + atomOffset);
+      }
+      info.myAtoms[j]->setType( currentAtom->getType() );
+    
+#ifdef IS_MPI
       
-      molIndex++;
-    }      
-  }
+      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
+      
+#endif // is_mpi
+    } 
+    
+    // make the bonds
+    for(j=0; j<info.nBonds; j++){
+      
+      currentBond = comp_stamps[stampID]->getBond( j );
+      theBonds[j].a = currentBond->getA() + atomOffset;
+      theBonds[j].b = currentBond->getB() + atomOffset;
 
-  the_ff->initializeBonds( the_bonds );
-}
+      exI = theBonds[j].a;
+      exJ = theBonds[j].b;
 
-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++ ){
-
+      // 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	
+      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
 
-	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
+      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;
 	  
-	  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() ){
+      if( currentBend->haveExtras() ){
 	    
-	    extras = current_bend->getExtras();
-	    current_extra = extras;
+	extras = currentBend->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 neither 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;
-	  }
-	  
-	  // exclude_I must always be the smaller of the pair
-	  if( exI > exJ ){
-	    tempEx = exI;
-	    exI = exJ;
-	    exJ = tempEx;
+	    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();
 	  }
-
-
-#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
+      
+      // 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;
 
-      molIndex++;
+    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;
+
+    
+    // send the arrays off to the forceField for init.
+
+    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 );
 
-#ifdef IS_MPI
-  sprintf( checkPointMsg,
-	   "Successfully created the bends list.\n" );
-  MPIcheckPoint();
-#endif // is_mpi
-  
 
-  the_ff->initializeBends( the_bends );
-}
+    the_molecules[i].initialize( info );
 
-void SimSetup::makeTorsions( void ){
 
-  int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
-  torsion_set* the_torsions;
-  TorsionStamp* current_torsion;
+    atomOffset += info.nAtoms;
+    delete[] theBonds;
+    delete[] theBends;
+    delete[] theTorsions;
+  }
 
-  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++ ){
-
 #ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
+  sprintf( checkPointMsg, "all molecules initialized succesfully" );
+  MPIcheckPoint();
+#endif // is_mpi
 
-      for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
+  // clean up the forcefield
+  the_ff->calcRcut();
+  the_ff->cleanMe();
 
-	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;
-	
-#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)
-#endif  //is_mpi
-	
-
-	// increment the index and repeat;
-	index++;
-      }
-      offset += comp_stamps[i]->getNAtoms();
-
-#ifdef IS_MPI
-      }
-#endif //is_mpi      
-
-      molIndex++;
-    }
-  }
-
-  the_ff->initializeTorsions( the_torsions );
 }
 
 void SimSetup::initFromBass( void ){