--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/07/15 03:27:24	605
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/07/15 17:57:04	614
@@ -82,705 +82,51 @@ void SimSetup::createSim( void ){
   Globals* the_globals;
   int i, j, k, globalAtomIndex;
   
-  int ensembleCase;
-  int ffCase;
+  // gather all of the information from the Bass file
   
-  ensembleCase = -1;
-  ffCase = -1;
+  gatherInfo();
 
-  // get the stamps and globals;
-  the_stamps = stamps;
-  the_globals = globals;
+  // creation of complex system objects
 
-  // set the easy ones first
-  simnfo->target_temp = the_globals->getTargetTemp();
-  simnfo->dt = the_globals->getDt();
-  simnfo->run_time = the_globals->getRunTime();
+  sysObjectsCreation();
 
-  // 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() );
+  // initialize the arrays
 
-  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
-  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
-  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") ) 
-    ensembleCase = NPTi_ENS;
-  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
-  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
-  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_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 ){
-
-  case FF_DUFF:
-    the_ff = new DUFF();
-    usesDipoles = 1;
-    break;
-
-  case FF_LJ:
-    the_ff = new LJFF();
-    break;
-
-  default:
-    sprintf( painCave.errMsg,
-	     "SimSetup Error. Unrecognized force field in case statement.\n");
-    painCave.isFatal = 1;
-    simError();
-  }
-
-#ifdef IS_MPI
-  strcpy( checkPointMsg, "ForceField creation successful" );
-  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];
-  comp_stamps = new MoleculeStamp*[n_components];
-
-  if( !the_globals->haveNMol() ){
-    // we don't have the total number of molecules, so we assume it is
-    // given in each component
-
-    tot_nmol = 0;
-    for( i=0; i<n_components; i++ ){
-
-      if( !the_components[i]->haveNMol() ){
-	// we have a problem
-	sprintf( painCave.errMsg,
-		 "SimSetup Error. No global NMol or component NMol"
-		 " given. Cannot calculate the number of atoms.\n" );
-	painCave.isFatal = 1;
-	simError();
-      }
-
-      tot_nmol += the_components[i]->getNMol();
-      components_nmol[i] = the_components[i]->getNMol();
-    }
-  }
-  else{
-    sprintf( painCave.errMsg,
-	     "SimSetup error.\n"
-	     "\tSorry, the ability to specify total"
-	     " nMols and then give molfractions in the components\n"
-	     "\tis not currently supported."
-	     " Please give nMol in the components.\n" );
-    painCave.isFatal = 1;
-    simError();
-    
-    
-    //     tot_nmol = the_globals->getNMol();
-    
-    //   //we have the total number of molecules, now we check for molfractions
-    //     for( i=0; i<n_components; i++ ){
-    
-    //       if( !the_components[i]->haveMolFraction() ){
-    
-    // 	if( !the_components[i]->haveNMol() ){
-    // 	  //we have a problem
-    // 	  std::cerr << "SimSetup error. Neither molFraction nor "
-    // 		    << " nMol was given in component
-    
+  makeMolecules();
+  info->identArray = new int[info->n_atoms];
+  for(i=0; i<info->n_atoms; i++){
+    info->identArray[i] = the_atoms[i]->getIdent();
   }
-
-#ifdef IS_MPI
-  strcpy( checkPointMsg, "Have the number of components" );
-  MPIcheckPoint();
-#endif // is_mpi
-
-  // make an array of molecule stamps that match the components used.
-  // also extract the used stamps out into a separate linked list
-
-  simnfo->nComponents = n_components;
-  simnfo->componentsNmol = components_nmol;
-  simnfo->compStamps = comp_stamps;
-  simnfo->headStamp = new LinkedMolStamp();
   
-  char* id;
-  LinkedMolStamp* headStamp = simnfo->headStamp;
-  LinkedMolStamp* currentStamp = NULL;
-  for( i=0; i<n_components; i++ ){
 
-    id = the_components[i]->getType();
-    comp_stamps[i] = NULL;
-    
-    // check to make sure the component isn't already in the list
-
-    comp_stamps[i] = headStamp->match( id );
-    if( comp_stamps[i] == NULL ){
-      
-      // extract the component from the list;
-      
-      currentStamp = the_stamps->extractMolStamp( id );
-      if( currentStamp == NULL ){
-	sprintf( painCave.errMsg,
-		 "SimSetup error: Component \"%s\" was not found in the "
-		 "list of declared molecules\n",
-		 id );
-	painCave.isFatal = 1;
-	simError();
-      }
-      
-      headStamp->add( currentStamp );
-      comp_stamps[i] = headStamp->match( id );
-    }
-  }
-
-#ifdef IS_MPI
-  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
-  MPIcheckPoint();
-#endif // is_mpi
+  // check on the post processing info
   
+  finalInfoCheck();
 
 
 
-  // caclulate the number of atoms, bonds, bends and torsions
 
-  tot_atoms = 0;
-  tot_bonds = 0;
-  tot_bends = 0;
-  tot_torsions = 0;
-  for( i=0; i<n_components; i++ ){
-    
-    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
-    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
-    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
-    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
-  }
+  // initialize the system coordinates
 
-  tot_SRI = tot_bonds + tot_bends + tot_torsions;
-
-  simnfo->n_atoms = tot_atoms;
-  simnfo->n_bonds = tot_bonds;
-  simnfo->n_bends = tot_bends;
-  simnfo->n_torsions = tot_torsions;
-  simnfo->n_SRI = tot_SRI;
-  simnfo->n_mol = tot_nmol;
+  initSystemCoords();
   
-  simnfo->molMembershipArray = new int[tot_atoms];
 
-#ifdef IS_MPI
+  // make the output filenames
 
-  // 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;
-  local_atoms = 0;
-  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( mol2proc[allMol] == worldRank ){
-	
-	local_atoms +=    comp_stamps[i]->getNAtoms();
-	local_bonds +=    comp_stamps[i]->getNBonds();
-	local_bends +=    comp_stamps[i]->getNBends();
-	local_torsions += comp_stamps[i]->getNTorsions();
-	localMol++;
-      }      
-      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 not equal.\n",
-	     simnfo->n_atoms,
-	     local_atoms );
-    painCave.isFatal = 1;
-    simError();
-  }
-
-  simnfo->n_bonds = local_bonds;
-  simnfo->n_bends = local_bends;
-  simnfo->n_torsions = local_torsions;
-  simnfo->n_SRI = local_SRI;
-  simnfo->n_mol = localMol;
-
-  strcpy( checkPointMsg, "Passed nlocal consistency check." );
-  MPIcheckPoint();
-  
-  
-#endif // is_mpi
-  
-
-  // create the atom and short range interaction arrays
-
-  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
-
-#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 = simnfo->n_SRI;
-  }
-  else{
-    
-    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 = 0;
-  }
-
-  // set the arrays into the SimInfo object
-
-  simnfo->atoms = the_atoms;
-  simnfo->molecules = the_molecules;
-  simnfo->nGlobalExcludes = 0;
-  simnfo->excludes = the_excludes;
-
-
-  // get some of the tricky things that may still be in the globals
-
-  double boxVector[3];
-  if( the_globals->haveBox() ){
-    boxVector[0] = the_globals->getBox();
-    boxVector[1] = the_globals->getBox();
-    boxVector[2] = the_globals->getBox();
-    
-    simnfo->setBox( boxVector );
-  }
-  else if( the_globals->haveDensity() ){
-
-    double vol;
-    vol = (double)tot_nmol / the_globals->getDensity();
-     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
-     boxVector[1] = boxVector[0];
-     boxVector[2] = boxVector[0];
-
-    simnfo->setBox( boxVector );
-  }
-  else{
-    if( !the_globals->haveBoxX() ){
-      sprintf( painCave.errMsg,
-	       "SimSetup error, no periodic BoxX size given.\n" );
-      painCave.isFatal = 1;
-      simError();
-    }
-    boxVector[0] = the_globals->getBoxX();
-
-    if( !the_globals->haveBoxY() ){
-      sprintf( painCave.errMsg,
-	       "SimSetup error, no periodic BoxY size given.\n" );
-      painCave.isFatal = 1;
-      simError();
-    }
-    boxVector[1] = the_globals->getBoxY();
-
-    if( !the_globals->haveBoxZ() ){
-      sprintf( painCave.errMsg,
-	       "SimSetup error, no periodic BoxZ size given.\n" );
-      painCave.isFatal = 1;
-      simError();
-    }
-    boxVector[2] = the_globals->getBoxZ();
-
-    simnfo->setBox( boxVector );
-  }
-
-#ifdef IS_MPI
-  strcpy( checkPointMsg, "Box size set up" );
-  MPIcheckPoint();
-#endif // is_mpi
-
-
-  // initialize the arrays
-
-  the_ff->setSimInfo( simnfo );
-
-  makeMolecules();
-  simnfo->identArray = new int[simnfo->n_atoms];
-  for(i=0; i<simnfo->n_atoms; i++){
-    simnfo->identArray[i] = the_atoms[i]->getIdent();
-  }
-  
-  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->boxLx;
-      if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
-      if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
-      simnfo->ecr = 0.5 * smallest;
-    } else {
-      simnfo->ecr        = the_globals->getECR();
-    }
-
-    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->boxLx;
-        if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
-        if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
-        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() ){
- 
-     InitializeFromFile* fileInit;
-#ifdef IS_MPI // is_mpi
-     if( worldRank == 0 ){
-#endif //is_mpi
-   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
-#ifdef IS_MPI
-     }else fileInit = new InitializeFromFile( NULL );
-#endif
-   fileInit->read_xyz( simnfo ); // default velocities on
-
-   delete fileInit;
- }
- else{
-
-#ifdef IS_MPI 
-
-  // no init from bass
-  
-  sprintf( painCave.errMsg,
-	   "Cannot intialize a parallel simulation without an initial configuration file.\n" );
-  painCave.isFatal;
-  simError();
-  
-#else
-
-  initFromBass();
-
-
-#endif
- }
-
-#ifdef IS_MPI
-  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
-  MPIcheckPoint();
-#endif // is_mpi
-
-
+  makeOutNames();
    
  
   
 
   
-#ifdef IS_MPI
-  if( worldRank == 0 ){
-#endif // is_mpi
-    
-    if( the_globals->haveFinalConfig() ){
-      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
-    }
-    else{
-      strcpy( simnfo->finalName, inFileName );
-      char* endTest;
-      int nameLength = strlen( simnfo->finalName );
-      endTest = &(simnfo->finalName[nameLength - 5]);
-      if( !strcmp( endTest, ".bass" ) ){
-	strcpy( endTest, ".eor" );
-      }
-      else if( !strcmp( endTest, ".BASS" ) ){
-	strcpy( endTest, ".eor" );
-      }
-      else{
-	endTest = &(simnfo->finalName[nameLength - 4]);
-	if( !strcmp( endTest, ".bss" ) ){
-	  strcpy( endTest, ".eor" );
-	}
-	else if( !strcmp( endTest, ".mdl" ) ){
-	  strcpy( endTest, ".eor" );
-	}
-	else{
-	  strcat( simnfo->finalName, ".eor" );
-	}
-      }
-    }
-    
-    // make the sample and status out names
-    
-    strcpy( simnfo->sampleName, inFileName );
-    char* endTest;
-    int nameLength = strlen( simnfo->sampleName );
-    endTest = &(simnfo->sampleName[nameLength - 5]);
-    if( !strcmp( endTest, ".bass" ) ){
-      strcpy( endTest, ".dump" );
-    }
-    else if( !strcmp( endTest, ".BASS" ) ){
-      strcpy( endTest, ".dump" );
-    }
-    else{
-      endTest = &(simnfo->sampleName[nameLength - 4]);
-      if( !strcmp( endTest, ".bss" ) ){
-	strcpy( endTest, ".dump" );
-      }
-      else if( !strcmp( endTest, ".mdl" ) ){
-	strcpy( endTest, ".dump" );
-      }
-      else{
-	strcat( simnfo->sampleName, ".dump" );
-      }
-    }
-    
-    strcpy( simnfo->statusName, inFileName );
-    nameLength = strlen( simnfo->statusName );
-    endTest = &(simnfo->statusName[nameLength - 5]);
-    if( !strcmp( endTest, ".bass" ) ){
-      strcpy( endTest, ".stat" );
-    }
-    else if( !strcmp( endTest, ".BASS" ) ){
-      strcpy( endTest, ".stat" );
-    }
-    else{
-      endTest = &(simnfo->statusName[nameLength - 4]);
-      if( !strcmp( endTest, ".bss" ) ){
-	strcpy( endTest, ".stat" );
-      }
-      else if( !strcmp( endTest, ".mdl" ) ){
-	strcpy( endTest, ".stat" );
-      }
-      else{
-	strcat( simnfo->statusName, ".stat" );
-      }
-    }
-    
-#ifdef IS_MPI
-  }
-#endif // is_mpi
+
   
-  // set the status, sample, and themal kick times
-  
-  if( the_globals->haveSampleTime() ){
-    simnfo->sampleTime = the_globals->getSampleTime();
-    simnfo->statusTime = simnfo->sampleTime;
-    simnfo->thermalTime = simnfo->sampleTime;
-  }
-  else{
-    simnfo->sampleTime = the_globals->getRunTime();
-    simnfo->statusTime = simnfo->sampleTime;
-    simnfo->thermalTime = simnfo->sampleTime;
-  }
 
-  if( the_globals->haveStatusTime() ){
-    simnfo->statusTime = the_globals->getStatusTime();
-  }
 
-  if( the_globals->haveThermalTime() ){
-    simnfo->thermalTime = the_globals->getThermalTime();
-  }
-
-  // check for the temperature set flag
-
-  if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
-
-
   // make the integrator
   
   
@@ -793,11 +139,11 @@ void SimSetup::createSim( void ){
   switch( ensembleCase ){
 
   case NVE_ENS:
-    new NVE( simnfo, the_ff );
+    new NVE( info, the_ff );
     break;
 
   case NVT_ENS:
-    myNVT = new NVT( simnfo, the_ff );
+    myNVT = new NVT( info, the_ff );
     myNVT->setTargetTemp(the_globals->getTargetTemp());
 
     if (the_globals->haveTauThermostat()) 
@@ -813,8 +159,8 @@ void SimSetup::createSim( void ){
     break;
 
   case NPTi_ENS:
-    myNPTi = new NPTi( simnfo, the_ff );
-    myNPTi->setTargetTemp( the_globals->getTargetTemp());
+    myNPTi = new NPTi( info, the_ff );
+    myNPTi->setTargetTemp( the_globals->getTargetTemp() );
 
     if (the_globals->haveTargetPressure())
       myNPTi->setTargetPressure(the_globals->getTargetPressure());
@@ -848,7 +194,7 @@ void SimSetup::createSim( void ){
     break;
 
   case NPTf_ENS:
-    myNPTf = new NPTf( simnfo, the_ff );
+    myNPTf = new NPTf( info, the_ff );
     myNPTf->setTargetTemp( the_globals->getTargetTemp());
 
     if (the_globals->haveTargetPressure())
@@ -883,7 +229,7 @@ void SimSetup::createSim( void ){
     break;
     
   case NPTim_ENS:
-    myNPTim = new NPTim( simnfo, the_ff );
+    myNPTim = new NPTim( info, the_ff );
     myNPTim->setTargetTemp( the_globals->getTargetTemp());
 
     if (the_globals->haveTargetPressure())
@@ -918,7 +264,7 @@ void SimSetup::createSim( void ){
     break;
 
   case NPTfm_ENS:
-    myNPTfm = new NPTfm( simnfo, the_ff );
+    myNPTfm = new NPTfm( info, the_ff );
     myNPTfm->setTargetTemp( the_globals->getTargetTemp());
 
     if (the_globals->haveTargetPressure())
@@ -967,18 +313,18 @@ void SimSetup::createSim( void ){
   // initialize the Fortran
 
 
-  simnfo->refreshSim();
+  info->refreshSim();
   
-  if( !strcmp( simnfo->mixingRule, "standard") ){
+  if( !strcmp( info->mixingRule, "standard") ){
     the_ff->initForceField( LB_MIXING_RULE );
   }
-  else if( !strcmp( simnfo->mixingRule, "explicit") ){
+  else if( !strcmp( info->mixingRule, "explicit") ){
     the_ff->initForceField( EXPLICIT_MIXING_RULE );
   }
   else{
     sprintf( painCave.errMsg,
 	     "SimSetup Error: unknown mixing rule -> \"%s\"\n",
-	     simnfo->mixingRule );
+	     info->mixingRule );
     painCave.isFatal = 1;
     simError();
   }
@@ -1020,7 +366,7 @@ void SimSetup::makeMolecules( void ){
   
   atomOffset = 0;
   excludeOffset = 0;
-  for(i=0; i<simnfo->n_mol; i++){
+  for(i=0; i<info->n_mol; i++){
     
     stampID = the_molecules[i].getStampID();
 
@@ -1048,7 +394,7 @@ void SimSetup::makeMolecules( void ){
       if( currentAtom->haveOrientation() ){
 	
 	dAtom = new DirectionalAtom(j + atomOffset);
-	simnfo->n_oriented++;
+	info->n_oriented++;
 	info.myAtoms[j] = dAtom;
 	
 	ux = currentAtom->getOrntX();
@@ -1272,9 +618,9 @@ void SimSetup::initFromBass( void ){
     have_extra =1;
 
     n_cells = (int)temp3 - 1;
-    cellx = simnfo->boxLx / temp3;
-    celly = simnfo->boxLy / temp3;
-    cellz = simnfo->boxLz / temp3;
+    cellx = info->boxLx / temp3;
+    celly = info->boxLy / temp3;
+    cellz = info->boxLz / temp3;
     n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
     temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
     n_per_extra = (int)ceil( temp1 );
@@ -1289,9 +635,9 @@ void SimSetup::initFromBass( void ){
   }
   else{
     n_cells = (int)temp3;
-    cellx = simnfo->boxLx / temp3;
-    celly = simnfo->boxLy / temp3;
-    cellz = simnfo->boxLz / temp3;
+    cellx = info->boxLx / temp3;
+    celly = info->boxLy / temp3;
+    cellz = info->boxLz / temp3;
   }
 
   current_mol = 0;
@@ -1371,10 +717,10 @@ void SimSetup::initFromBass( void ){
   }
 
 
-  for( i=0; i<simnfo->n_atoms; i++ ){
-    simnfo->atoms[i]->set_vx( 0.0 );
-    simnfo->atoms[i]->set_vy( 0.0 );
-    simnfo->atoms[i]->set_vz( 0.0 );
+  for( i=0; i<info->n_atoms; i++ ){
+    info->atoms[i]->set_vx( 0.0 );
+    info->atoms[i]->set_vy( 0.0 );
+    info->atoms[i]->set_vz( 0.0 );
   }
 }
 
@@ -1433,5 +779,700 @@ void SimSetup::makeElement( double x, double y, double
 
     current_comp_mol = 0;
     current_comp++;
+  }
+}
+
+
+void SimSetup::gatherInfo( void ){
+
+  ensembleCase = -1;
+  ffCase = -1;
+
+  // get the stamps and globals;
+  the_stamps = stamps;
+  the_globals = globals;
+
+  // set the easy ones first
+  info->target_temp = the_globals->getTargetTemp();
+  info->dt = the_globals->getDt();
+  info->run_time = the_globals->getRunTime();
+  n_components = the_globals->getNComponents();
+
+
+  // get the forceField
+
+  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" ))      ensembleCase = NVE_ENS;
+  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
+  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") ) 
+    ensembleCase = NPTi_ENS;
+  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
+  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
+  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_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( info->ensemble, ensemble );
+
+  // get the mixing rule
+
+  strcpy( info->mixingRule, the_globals->getMixingRule() );
+  info->usePBC = the_globals->getPBC();
+	
+  
+  // get the components and calculate the tot_nMol and indvidual n_mol
+ 
+  the_components = the_globals->getComponents();
+  components_nmol = new int[n_components];
+
+
+  if( !the_globals->haveNMol() ){
+    // we don't have the total number of molecules, so we assume it is
+    // given in each component
+
+    tot_nmol = 0;
+    for( i=0; i<n_components; i++ ){
+
+      if( !the_components[i]->haveNMol() ){
+	// we have a problem
+	sprintf( painCave.errMsg,
+		 "SimSetup Error. No global NMol or component NMol"
+		 " given. Cannot calculate the number of atoms.\n" );
+	painCave.isFatal = 1;
+	simError();
+      }
+
+      tot_nmol += the_components[i]->getNMol();
+      components_nmol[i] = the_components[i]->getNMol();
+    }
+  }
+  else{
+    sprintf( painCave.errMsg,
+	     "SimSetup error.\n"
+	     "\tSorry, the ability to specify total"
+	     " nMols and then give molfractions in the components\n"
+	     "\tis not currently supported."
+	     " Please give nMol in the components.\n" );
+    painCave.isFatal = 1;
+    simError();
+  }
+
+  // set the status, sample, and thermal kick times
+  
+  if( the_globals->haveSampleTime() ){
+    info->sampleTime = the_globals->getSampleTime();
+    info->statusTime = info->sampleTime;
+    info->thermalTime = info->sampleTime;
+  }
+  else{
+    info->sampleTime = the_globals->getRunTime();
+    info->statusTime = info->sampleTime;
+    info->thermalTime = info->sampleTime;
+  }
+
+  if( the_globals->haveStatusTime() ){
+    info->statusTime = the_globals->getStatusTime();
+  }
+
+  if( the_globals->haveThermalTime() ){
+    info->thermalTime = the_globals->getThermalTime();
+  }
+
+  // check for the temperature set flag
+
+  if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
+
+  // get some of the tricky things that may still be in the globals
+
+  double boxVector[3];
+  if( the_globals->haveBox() ){
+    boxVector[0] = the_globals->getBox();
+    boxVector[1] = the_globals->getBox();
+    boxVector[2] = the_globals->getBox();
+    
+    info->setBox( boxVector );
+  }
+  else if( the_globals->haveDensity() ){
+
+    double vol;
+    vol = (double)tot_nmol / the_globals->getDensity();
+     boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
+     boxVector[1] = boxVector[0];
+     boxVector[2] = boxVector[0];
+
+    info->setBox( boxVector );
+  }
+  else{
+    if( !the_globals->haveBoxX() ){
+      sprintf( painCave.errMsg,
+	       "SimSetup error, no periodic BoxX size given.\n" );
+      painCave.isFatal = 1;
+      simError();
+    }
+    boxVector[0] = the_globals->getBoxX();
+
+    if( !the_globals->haveBoxY() ){
+      sprintf( painCave.errMsg,
+	       "SimSetup error, no periodic BoxY size given.\n" );
+      painCave.isFatal = 1;
+      simError();
+    }
+    boxVector[1] = the_globals->getBoxY();
+
+    if( !the_globals->haveBoxZ() ){
+      sprintf( painCave.errMsg,
+	       "SimSetup error, no periodic BoxZ size given.\n" );
+      painCave.isFatal = 1;
+      simError();
+    }
+    boxVector[2] = the_globals->getBoxZ();
+
+    info->setBox( boxVector );
   }
+
+
+    
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "Succesfully gathered all information from Bass\n" );
+  MPIcheckPoint();
+#endif // is_mpi
+
 }
+
+
+void SimSetup::finalInfoCheck( void ){
+  int index;
+  int usesDipoles;
+  
+
+  // check electrostatic parameters
+  
+  index = 0;
+  usesDipoles = 0;
+  while( (index < info->n_atoms) && !usesDipoles ){
+    usesDipoles = ((info->atoms)[index])->hasDipole();
+    index++;
+  }
+  
+#ifdef IS_MPI
+  int myUse = usesDipoles
+  MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
+#endif //is_mpi
+
+
+  if (the_globals->getUseRF() ) {
+    info->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 = info->boxLx;
+      if (info->boxLy <= smallest) smallest = info->boxLy;
+      if (info->boxLz <= smallest) smallest = info->boxLz;
+      info->ecr = 0.5 * smallest;
+    } else {
+      info->ecr        = the_globals->getECR();
+    }
+
+    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();
+      info->est = 0.05 * info->ecr;
+    } else {
+      info->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();
+    }
+    info->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 = info->boxLx;
+        if (info->boxLy <= smallest) smallest = info->boxLy;
+        if (info->boxLz <= smallest) smallest = info->boxLz;
+        info->ecr = 0.5 * smallest;
+      } else {
+        info->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();
+        info->est = 0.05 * info->ecr;
+      } else {
+        info->est        = the_globals->getEST();
+      }
+    }
+  }  
+
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "post processing checks out" );
+  MPIcheckPoint();
+#endif // is_mpi
+
+}
+
+void SimSetup::initSystemCoords( void ){
+
+ if( the_globals->haveInitialConfig() ){
+ 
+     InitializeFromFile* fileInit;
+#ifdef IS_MPI // is_mpi
+     if( worldRank == 0 ){
+#endif //is_mpi
+   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
+#ifdef IS_MPI
+     }else fileInit = new InitializeFromFile( NULL );
+#endif
+   fileInit->read_xyz( info ); // default velocities on
+
+   delete fileInit;
+ }
+ else{
+
+#ifdef IS_MPI 
+
+  // no init from bass
+  
+  sprintf( painCave.errMsg,
+	   "Cannot intialize a parallel simulation without an initial configuration file.\n" );
+  painCave.isFatal;
+  simError();
+  
+#else
+
+  initFromBass();
+
+
+#endif
+ }
+
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
+  MPIcheckPoint();
+#endif // is_mpi
+
+}
+
+
+void SimSetup::makeOutNames( void ){
+
+#ifdef IS_MPI
+  if( worldRank == 0 ){
+#endif // is_mpi
+    
+    if( the_globals->haveFinalConfig() ){
+      strcpy( info->finalName, the_globals->getFinalConfig() );
+    }
+    else{
+      strcpy( info->finalName, inFileName );
+      char* endTest;
+      int nameLength = strlen( info->finalName );
+      endTest = &(info->finalName[nameLength - 5]);
+      if( !strcmp( endTest, ".bass" ) ){
+	strcpy( endTest, ".eor" );
+      }
+      else if( !strcmp( endTest, ".BASS" ) ){
+	strcpy( endTest, ".eor" );
+      }
+      else{
+	endTest = &(info->finalName[nameLength - 4]);
+	if( !strcmp( endTest, ".bss" ) ){
+	  strcpy( endTest, ".eor" );
+	}
+	else if( !strcmp( endTest, ".mdl" ) ){
+	  strcpy( endTest, ".eor" );
+	}
+	else{
+	  strcat( info->finalName, ".eor" );
+	}
+      }
+    }
+    
+    // make the sample and status out names
+    
+    strcpy( info->sampleName, inFileName );
+    char* endTest;
+    int nameLength = strlen( info->sampleName );
+    endTest = &(info->sampleName[nameLength - 5]);
+    if( !strcmp( endTest, ".bass" ) ){
+      strcpy( endTest, ".dump" );
+    }
+    else if( !strcmp( endTest, ".BASS" ) ){
+      strcpy( endTest, ".dump" );
+    }
+    else{
+      endTest = &(info->sampleName[nameLength - 4]);
+      if( !strcmp( endTest, ".bss" ) ){
+	strcpy( endTest, ".dump" );
+      }
+      else if( !strcmp( endTest, ".mdl" ) ){
+	strcpy( endTest, ".dump" );
+      }
+      else{
+	strcat( info->sampleName, ".dump" );
+      }
+    }
+    
+    strcpy( info->statusName, inFileName );
+    nameLength = strlen( info->statusName );
+    endTest = &(info->statusName[nameLength - 5]);
+    if( !strcmp( endTest, ".bass" ) ){
+      strcpy( endTest, ".stat" );
+    }
+    else if( !strcmp( endTest, ".BASS" ) ){
+      strcpy( endTest, ".stat" );
+    }
+    else{
+      endTest = &(info->statusName[nameLength - 4]);
+      if( !strcmp( endTest, ".bss" ) ){
+	strcpy( endTest, ".stat" );
+      }
+      else if( !strcmp( endTest, ".mdl" ) ){
+	strcpy( endTest, ".stat" );
+      }
+      else{
+	strcat( info->statusName, ".stat" );
+      }
+    }
+    
+#ifdef IS_MPI
+  }
+#endif // is_mpi
+
+}
+
+
+void SimSetup::sysObjectsCreation( void ){
+
+  // create the forceField
+
+  createFF();
+
+  // extract componentList
+
+  compList();
+
+  // calc the number of atoms, bond, bends, and torsions
+
+  calcSysValues();
+
+#ifdef IS_MPI
+  // divide the molecules among the processors
+  
+  mpiMolDivide();
+#endif //is_mpi
+  
+  // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
+  
+  makeSysArrays();
+
+  
+
+
+}
+
+
+void SimSetup::createFF( void ){
+
+  switch( ffCase ){
+
+  case FF_DUFF:
+    the_ff = new DUFF();
+    break;
+
+  case FF_LJ:
+    the_ff = new LJFF();
+    break;
+
+  default:
+    sprintf( painCave.errMsg,
+	     "SimSetup Error. Unrecognized force field in case statement.\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "ForceField creation successful" );
+  MPIcheckPoint();
+#endif // is_mpi
+
+}
+
+
+void SimSetup::compList( void ){
+
+  comp_stamps = new MoleculeStamp*[n_components];
+
+  // make an array of molecule stamps that match the components used.
+  // also extract the used stamps out into a separate linked list
+
+  info->nComponents = n_components;
+  info->componentsNmol = components_nmol;
+  info->compStamps = comp_stamps;
+  info->headStamp = new LinkedMolStamp();
+  
+  char* id;
+  LinkedMolStamp* headStamp = info->headStamp;
+  LinkedMolStamp* currentStamp = NULL;
+  for( i=0; i<n_components; i++ ){
+
+    id = the_components[i]->getType();
+    comp_stamps[i] = NULL;
+    
+    // check to make sure the component isn't already in the list
+
+    comp_stamps[i] = headStamp->match( id );
+    if( comp_stamps[i] == NULL ){
+      
+      // extract the component from the list;
+      
+      currentStamp = the_stamps->extractMolStamp( id );
+      if( currentStamp == NULL ){
+	sprintf( painCave.errMsg,
+		 "SimSetup error: Component \"%s\" was not found in the "
+		 "list of declared molecules\n",
+		 id );
+	painCave.isFatal = 1;
+	simError();
+      }
+      
+      headStamp->add( currentStamp );
+      comp_stamps[i] = headStamp->match( id );
+    }
+  }
+
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
+  MPIcheckPoint();
+#endif // is_mpi
+
+
+}
+
+void SimSetup::calcSysValues( void ){
+
+  tot_atoms = 0;
+  tot_bonds = 0;
+  tot_bends = 0;
+  tot_torsions = 0;
+  for( i=0; i<n_components; i++ ){
+    
+    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
+    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
+    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
+    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
+  }
+
+  tot_SRI = tot_bonds + tot_bends + tot_torsions;
+
+  info->n_atoms = tot_atoms;
+  info->n_bonds = tot_bonds;
+  info->n_bends = tot_bends;
+  info->n_torsions = tot_torsions;
+  info->n_SRI = tot_SRI;
+  info->n_mol = tot_nmol;
+  
+  info->molMembershipArray = new int[tot_atoms];
+}
+
+
+#ifdef IS_MPI
+
+void SimSetup::mpiMolDivide( void ){
+  
+  int localMol, allMol;
+  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
+
+  mpiSim = new mpiSimulation( info );
+  
+  globalIndex = mpiSim->divideLabor();
+
+  // set up the local variables 
+  
+  mol2proc = mpiSim->getMolToProcMap();
+  molCompType = mpiSim->getMolComponentType();
+  
+  allMol = 0;
+  localMol = 0;
+  local_atoms = 0;
+  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( mol2proc[allMol] == worldRank ){
+	
+	local_atoms +=    comp_stamps[i]->getNAtoms();
+	local_bonds +=    comp_stamps[i]->getNBonds();
+	local_bends +=    comp_stamps[i]->getNBends();
+	local_torsions += comp_stamps[i]->getNTorsions();
+	localMol++;
+      }      
+      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
+        info->molMembershipArray[globalAtomIndex] = allMol;
+        globalAtomIndex++;
+      }
+
+      allMol++;      
+    }
+  }
+  local_SRI = local_bonds + local_bends + local_torsions;
+  
+  info->n_atoms = mpiSim->getMyNlocal();  
+  
+  if( local_atoms != info->n_atoms ){
+    sprintf( painCave.errMsg,
+	     "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
+	     " localAtom (%d) are not equal.\n",
+	     info->n_atoms,
+	     local_atoms );
+    painCave.isFatal = 1;
+    simError();
+  }
+
+  info->n_bonds = local_bonds;
+  info->n_bends = local_bends;
+  info->n_torsions = local_torsions;
+  info->n_SRI = local_SRI;
+  info->n_mol = localMol;
+
+  strcpy( checkPointMsg, "Passed nlocal consistency check." );
+  MPIcheckPoint();
+}
+  
+#endif // is_mpi
+
+
+void SimSetup::makeSysArrays( void ){
+
+  // create the atom and short range interaction arrays
+
+  Atom::createArrays(info->n_atoms);
+  the_atoms = new Atom*[info->n_atoms];
+  the_molecules = new Molecule[info->n_mol];
+  int molIndex;
+
+  // initialize the molecule's stampID's
+
+#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++) {
+        info->molMembershipArray[globalAtomIndex] = molIndex;
+        globalAtomIndex++;
+      }
+      molIndex++;
+    }
+  }
+    
+
+#endif // is_mpi
+
+
+  if( info->n_SRI ){
+    
+    Exclude::createArray(info->n_SRI);
+    the_excludes = new Exclude*[info->n_SRI];
+    for( int ex=0; ex<info->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
+    info->globalExcludes = new int;
+    info->n_exclude = info->n_SRI;
+  }
+  else{
+    
+    Exclude::createArray( 1 );
+    the_excludes = new Exclude*;
+    the_excludes[0] = new Exclude(0);
+    the_excludes[0]->setPair( 0,0 );
+    info->globalExcludes = new int;
+    info->globalExcludes[0] = 0;
+    info->n_exclude = 0;
+  }
+
+  // set the arrays into the SimInfo object
+
+  info->atoms = the_atoms;
+  info->molecules = the_molecules;
+  info->nGlobalExcludes = 0;
+  info->excludes = the_excludes;
+
+  the_ff->setSimInfo( info );
+
+}