--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/21 17:42:12	378
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/04/10 20:08:56	489
@@ -66,6 +66,7 @@ void SimSetup::createSim( void ){
 
   MakeStamps *the_stamps;
   Globals* the_globals;
+  ExtendedSystem* the_extendedsystem;
   int i, j;
 
   // get the stamps and globals;
@@ -80,18 +81,85 @@ 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() );
+
+  // get the ensemble and set up an extended system if we need it:
   strcpy( ensemble, the_globals->getEnsemble() );
+  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();
+    }
+
+  } else if ( !strcasecmp( ensemble, "NVE") ) {
+  } else {
+    sprintf( painCave.errMsg,
+	     "SimSetup Warning. Unrecognized Ensemble -> %s, "
+             "reverting to NVE for this simulation.\n",
+	     ensemble );
+    painCave.isFatal = 0;
+    simError();
+    strcpy( ensemble, "NVE" );
+  }  
   strcpy( simnfo->ensemble, ensemble );
 
   strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
   simnfo->usePBC = the_globals->getPBC();
 	  
-
-
-  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();
+  int usesDipoles = 0;
+  if( !strcmp( force_field, "TraPPE_Ex" ) ){
+    the_ff = new TraPPE_ExFF();
+    usesDipoles = 1;
+  }
+  else if( !strcasecmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
   else{
     sprintf( painCave.errMsg,
 	     "SimSetup Error. Unrecognized force field -> %s\n",
@@ -242,12 +310,13 @@ void SimSetup::createSim( void ){
 
   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;
@@ -259,8 +328,7 @@ void SimSetup::createSim( void ){
 
     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();
@@ -279,7 +347,7 @@ void SimSetup::createSim( void ){
   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,59 +372,69 @@ 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;
+  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 );
+      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;
 
 
   // get some of the tricky things that may still be in the globals
 
-  if( simnfo->n_dipoles ){
-
-    if( !the_globals->haveRRF() ){
-      sprintf( painCave.errMsg,
-	       "SimSetup Error, system has dipoles, but no rRF was set.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    if( !the_globals->haveDielectric() ){
-      sprintf( painCave.errMsg,
-	       "SimSetup Error, system has dipoles, but no"
-	       " dielectric was set.\n" );
-      painCave.isFatal = 1;
-      simError();
-    }
-
-    simnfo->rRF        = the_globals->getRRF();
-    simnfo->dielectric = the_globals->getDielectric();
-  }
-
-#ifdef IS_MPI
-  strcpy( checkPointMsg, "rRf and dielectric check out" );
-  MPIcheckPoint();
-#endif // is_mpi
   
   if( the_globals->haveBox() ){
     simnfo->box_x = the_globals->getBox();
@@ -407,28 +485,94 @@ void SimSetup::createSim( void ){
 
   the_ff->setSimInfo( simnfo );
 
-  makeAtoms();
+  makeMolecules();
   simnfo->identArray = new int[simnfo->n_atoms];
+  simnfo->molMembershipArray = new int[simnfo->n_atoms];
   for(i=0; i<simnfo->n_atoms; i++){
     simnfo->identArray[i] = the_atoms[i]->getIdent();
   }
-  
-  if( tot_bonds ){
-    makeBonds();
+  for(i=0; i< simnfo->n_mol; i++) {
+    the_molecules[i].atomicRollCall(simnfo->molMembershipArray);
   }
+  
+  if (the_globals->getUseRF() ) {
+    simnfo->useReactionField = 1;
+  
+    if( !the_globals->haveECR() ){
+      sprintf( painCave.errMsg,
+     	       "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;
+      simError();
+      double smallest;
+      smallest = simnfo->box_x;
+      if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
+      if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
+      simnfo->ecr = 0.5 * smallest;
+    } else {
+      simnfo->ecr        = the_globals->getECR();
+    }
 
-  if( tot_bends ){
-    makeBends();
-  }
-
-  if( tot_torsions ){
-    makeTorsions();
-  }
-
-
-
-
+    if( !the_globals->haveEST() ){
+      sprintf( painCave.errMsg,
+     	       "SimSetup Warning: using default value of 0.05 * the "
+               "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
+               );
+      painCave.isFatal = 0;
+      simError();
+      simnfo->est = 0.05 * simnfo->ecr;
+    } else {
+      simnfo->est        = the_globals->getEST();
+    }
+    
+    if(!the_globals->haveDielectric() ){
+      sprintf( painCave.errMsg,
+     	       "SimSetup Error: You are trying to use Reaction Field without"
+               "setting a dielectric constant!\n"
+               );
+      painCave.isFatal = 1;
+      simError();
+    }
+    simnfo->dielectric = the_globals->getDielectric();  
+  } else {
+    if (usesDipoles) {
+      
+      if( !the_globals->haveECR() ){
+        sprintf( painCave.errMsg,
+                 "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;
+        simError();
+        double smallest;
+        smallest = simnfo->box_x;
+        if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
+        if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
+        simnfo->ecr = 0.5 * smallest;
+      } else {
+        simnfo->ecr        = the_globals->getECR();
+      }
+      
+      if( !the_globals->haveEST() ){
+        sprintf( painCave.errMsg,
+                 "SimSetup Warning: using default value of 5%% of the "
+                 "electrostaticCutoffRadius for the "
+                 "electrostaticSkinThickness\n"
+                 );
+        painCave.isFatal = 0;
+        simError();
+        simnfo->est = 0.05 * simnfo->ecr;
+      } else {
+        simnfo->est        = the_globals->getEST();
+      }
+    }
+  }  
 
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "electrostatic parameters check out" );
+  MPIcheckPoint();
+#endif // is_mpi
 
  if( the_globals->haveInitialConfig() ){
  
@@ -587,15 +731,24 @@ void SimSetup::createSim( void ){
 
 //   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( force_field, "TraPPE_Ex" ) ){
+    new Symplectic(simnfo, the_ff, the_extendedsystem);
+  }
+  else if( !strcmp( force_field, "LJ" ) ){
+    new Verlet( *simnfo, the_ff, the_extendedsystem );
+  }
+  else {
+    std::cerr << "I'm a bug.\n";
+    fprintf( stderr, "Ima bug. stderr %s\n", force_field);
+  }
+#ifdef IS_MPI
+  mpiSim->mpiRefresh();
+#endif
 
-
   // initialize the Fortran
-  
+
+
   simnfo->refreshSim();
   
   if( !strcmp( simnfo->mixingRule, "standard") ){
@@ -620,358 +773,266 @@ void SimSetup::makeAtoms( 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;
+  LinkedAssign* extras;
+  LinkedAssign* current_extra;
+  AtomStamp* currentAtom;
+  BondStamp* currentBond;
+  BendStamp* currentBend;
+  TorsionStamp* currentTorsion;
 
-  lMolIndex = 0;
-  molIndex = 0;
-  index = 0;
-  for( i=0; i<n_components; i++ ){
+  bond_pair* theBonds;
+  bend_set* theBends;
+  torsion_set* theTorsions;
 
-    for( j=0; j<components_nmol[i]; j++ ){
+  
+  //init the forceField paramters
 
-#ifdef IS_MPI
-      if( mpiSim->getMyMolStart() <= molIndex &&
-	  molIndex <= mpiSim->getMyMolEnd() ){
-#endif // is_mpi	
+  the_ff->readParams();
 
-	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 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();
-}
-
-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
+    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;
 
-	  the_excludes[index*2] =     
-	    the_atoms[exI]->getGlobalIndex() + 1;
-	  the_excludes[index*2 + 1] = 
-	    the_atoms[exJ]->getGlobalIndex() + 1;
+    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];
 
-#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();
+    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() ){
 	
-#ifdef IS_MPI
+	dAtom = new DirectionalAtom(j + atomOffset);
+	simnfo->n_oriented++;
+	info.myAtoms[j] = dAtom;
+	
+	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;
+	    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;
-	  }
+	  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
-
-	  the_excludes[(index + tot_bonds)*2] =     
-	    the_atoms[exI]->getGlobalIndex() + 1;
-	  the_excludes[(index + tot_bonds)*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)*2] =     exI + 1; 
-	  the_excludes[(index + tot_bonds)*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
-	  
-	  
-	  // increment the index and repeat;
-	  index++;
-	}
-	offset += comp_stamps[i]->getNAtoms();
-	
-#ifdef IS_MPI
-      }
-#endif //is_mpi
-
-      molIndex++;
     }
-  }
+    excludeOffset += info.nTorsions;
+
+    
+    // send the arrays off to the forceField for init.
 
-#ifdef IS_MPI
-  sprintf( checkPointMsg,
-	   "Successfully created the bends list.\n" );
-  MPIcheckPoint();
-#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 );
 
-  the_ff->initializeBends( the_bends );
-}
 
-void SimSetup::makeTorsions( void ){
+    the_molecules[i].initialize( info );
 
-  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++ ){
+    atomOffset += info.nAtoms;
+    delete[] theBonds;
+    delete[] theBends;
+    delete[] theTorsions;
+  }
 
-    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 ){