--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/07/15 17:57:04	614
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/07/16 21:30:56	626
@@ -78,8 +78,6 @@ void SimSetup::createSim( void ){
 
 void SimSetup::createSim( void ){
 
-  MakeStamps *the_stamps;
-  Globals* the_globals;
   int i, j, k, globalAtomIndex;
   
   // gather all of the information from the Bass file
@@ -90,26 +88,10 @@ void SimSetup::createSim( void ){
 
   sysObjectsCreation();
 
-
-
-  // initialize the arrays
-
- 
-
-  makeMolecules();
-  info->identArray = new int[info->n_atoms];
-  for(i=0; i<info->n_atoms; i++){
-    info->identArray[i] = the_atoms[i]->getIdent();
-  }
-  
-
   // check on the post processing info
   
   finalInfoCheck();
 
-
-
-
   // initialize the system coordinates
 
   initSystemCoords();
@@ -119,229 +101,27 @@ void SimSetup::createSim( void ){
 
   makeOutNames();
    
- 
-  
-
-  
-
-  
-
-
   // make the integrator
   
+  makeIntegrator();
   
-  NVT*  myNVT = NULL;
-  NPTi* myNPTi = NULL;
-  NPTf* myNPTf = NULL;
-  NPTim* myNPTim = NULL;
-  NPTfm* myNPTfm = NULL;
-
-  switch( ensembleCase ){
-
-  case NVE_ENS:
-    new NVE( info, the_ff );
-    break;
-
-  case NVT_ENS:
-    myNVT = new NVT( info, the_ff );
-    myNVT->setTargetTemp(the_globals->getTargetTemp());
-
-    if (the_globals->haveTauThermostat()) 
-      myNVT->setTauThermostat(the_globals->getTauThermostat());
-
-    else {
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use the NVT\n"
-               "    ensemble, you must set tauThermostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    break;
-
-  case NPTi_ENS:
-    myNPTi = new NPTi( info, the_ff );
-    myNPTi->setTargetTemp( the_globals->getTargetTemp() );
-
-    if (the_globals->haveTargetPressure())
-      myNPTi->setTargetPressure(the_globals->getTargetPressure());
-    else {
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use a constant pressure\n"
-               "    ensemble, you must set targetPressure in the BASS file.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    
-    if( the_globals->haveTauThermostat() )
-      myNPTi->setTauThermostat( the_globals->getTauThermostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauThermostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-
-    if( the_globals->haveTauBarostat() )
-      myNPTi->setTauBarostat( the_globals->getTauBarostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauBarostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    break;
-
-  case NPTf_ENS:
-    myNPTf = new NPTf( info, the_ff );
-    myNPTf->setTargetTemp( the_globals->getTargetTemp());
-
-    if (the_globals->haveTargetPressure())
-      myNPTf->setTargetPressure(the_globals->getTargetPressure());
-    else {
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use a constant pressure\n"
-               "    ensemble, you must set targetPressure in the BASS file.\n");
-      painCave.isFatal = 1;
-      simError();
-    }    
-
-    if( the_globals->haveTauThermostat() )
-      myNPTf->setTauThermostat( the_globals->getTauThermostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauThermostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-
-    if( the_globals->haveTauBarostat() )
-      myNPTf->setTauBarostat( the_globals->getTauBarostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauBarostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    break;
-    
-  case NPTim_ENS:
-    myNPTim = new NPTim( info, the_ff );
-    myNPTim->setTargetTemp( the_globals->getTargetTemp());
-
-    if (the_globals->haveTargetPressure())
-      myNPTim->setTargetPressure(the_globals->getTargetPressure());
-    else {
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use a constant pressure\n"
-               "    ensemble, you must set targetPressure in the BASS file.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    
-    if( the_globals->haveTauThermostat() )
-      myNPTim->setTauThermostat( the_globals->getTauThermostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauThermostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-
-    if( the_globals->haveTauBarostat() )
-      myNPTim->setTauBarostat( the_globals->getTauBarostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauBarostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    break;
-
-  case NPTfm_ENS:
-    myNPTfm = new NPTfm( info, the_ff );
-    myNPTfm->setTargetTemp( the_globals->getTargetTemp());
-
-    if (the_globals->haveTargetPressure())
-      myNPTfm->setTargetPressure(the_globals->getTargetPressure());
-    else {
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use a constant pressure\n"
-               "    ensemble, you must set targetPressure in the BASS file.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    
-    if( the_globals->haveTauThermostat() )
-      myNPTfm->setTauThermostat( the_globals->getTauThermostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauThermostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-
-    if( the_globals->haveTauBarostat() )
-      myNPTfm->setTauBarostat( the_globals->getTauBarostat() );
-    else{
-      sprintf( painCave.errMsg,
-               "SimSetup error: If you use an NPT\n"
-               "    ensemble, you must set tauBarostat.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    break;
-
-  default:
-    sprintf( painCave.errMsg,
-	     "SimSetup Error. Unrecognized ensemble in case statement.\n");
-    painCave.isFatal = 1;
-    simError();
-  }
-
-
 #ifdef IS_MPI
   mpiSim->mpiRefresh();
 #endif
 
   // initialize the Fortran
 
+  initFortran();
 
-  info->refreshSim();
-  
-  if( !strcmp( info->mixingRule, "standard") ){
-    the_ff->initForceField( LB_MIXING_RULE );
-  }
-  else if( !strcmp( info->mixingRule, "explicit") ){
-    the_ff->initForceField( EXPLICIT_MIXING_RULE );
-  }
-  else{
-    sprintf( painCave.errMsg,
-	     "SimSetup Error: unknown mixing rule -> \"%s\"\n",
-	     info->mixingRule );
-    painCave.isFatal = 1;
-    simError();
-  }
 
 
-#ifdef IS_MPI
-  strcpy( checkPointMsg, 
-	  "Successfully intialized the mixingRule for Fortran." );
-  MPIcheckPoint();
-#endif // is_mpi
 }
 
 
 void SimSetup::makeMolecules( void ){
 
   int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
-  molInit info;
+  molInit molInfo;
   DirectionalAtom* dAtom;
   LinkedAssign* extras;
   LinkedAssign* current_extra;
@@ -370,32 +150,32 @@ void SimSetup::makeMolecules( void ){
     
     stampID = the_molecules[i].getStampID();
 
-    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;
+    molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
+    molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
+    molInfo.nBends    = comp_stamps[stampID]->getNBends();
+    molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
+    molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
 
-    info.myAtoms = &the_atoms[atomOffset];
-    info.myExcludes = &the_excludes[excludeOffset];
-    info.myBonds = new Bond*[info.nBonds];
-    info.myBends = new Bend*[info.nBends];
-    info.myTorsions = new Torsion*[info.nTorsions];
+    molInfo.myAtoms = &the_atoms[atomOffset];
+    molInfo.myExcludes = &the_excludes[excludeOffset];
+    molInfo.myBonds = new Bond*[molInfo.nBonds];
+    molInfo.myBends = new Bend*[molInfo.nBends];
+    molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
 
-    theBonds = new bond_pair[info.nBonds];
-    theBends = new bend_set[info.nBends];
-    theTorsions = new torsion_set[info.nTorsions];
+    theBonds = new bond_pair[molInfo.nBonds];
+    theBends = new bend_set[molInfo.nBends];
+    theTorsions = new torsion_set[molInfo.nTorsions];
     
     // make the Atoms
     
-    for(j=0; j<info.nAtoms; j++){
+    for(j=0; j<molInfo.nAtoms; j++){
       
       currentAtom = comp_stamps[stampID]->getAtom( j );
       if( currentAtom->haveOrientation() ){
 	
 	dAtom = new DirectionalAtom(j + atomOffset);
 	info->n_oriented++;
-	info.myAtoms[j] = dAtom;
+	molInfo.myAtoms[j] = dAtom;
 	
 	ux = currentAtom->getOrntX();
 	uy = currentAtom->getOrntY();
@@ -413,19 +193,19 @@ void SimSetup::makeMolecules( void ){
 	dAtom->setSUz( uz );
       }
       else{
-	info.myAtoms[j] = new GeneralAtom(j + atomOffset);
+	molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
       }
-      info.myAtoms[j]->setType( currentAtom->getType() );
+      molInfo.myAtoms[j]->setType( currentAtom->getType() );
     
 #ifdef IS_MPI
       
-      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
+      molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
       
 #endif // is_mpi
     } 
     
     // make the bonds
-    for(j=0; j<info.nBonds; j++){
+    for(j=0; j<molInfo.nBonds; j++){
       
       currentBond = comp_stamps[stampID]->getBond( j );
       theBonds[j].a = currentBond->getA() + atomOffset;
@@ -452,10 +232,10 @@ void SimSetup::makeMolecules( void ){
       the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
     }
-    excludeOffset += info.nBonds;
+    excludeOffset += molInfo.nBonds;
 
     //make the bends
-    for(j=0; j<info.nBends; j++){
+    for(j=0; j<molInfo.nBends; j++){
       
       currentBend = comp_stamps[stampID]->getBend( j );
       theBends[j].a = currentBend->getA() + atomOffset;
@@ -538,9 +318,9 @@ void SimSetup::makeMolecules( void ){
       the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
     }
-    excludeOffset += info.nBends;
+    excludeOffset += molInfo.nBends;
 
-    for(j=0; j<info.nTorsions; j++){
+    for(j=0; j<molInfo.nTorsions; j++){
       
       currentTorsion = comp_stamps[stampID]->getTorsion( j );
       theTorsions[j].a = currentTorsion->getA() + atomOffset;
@@ -568,21 +348,21 @@ void SimSetup::makeMolecules( void ){
       the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
     }
-    excludeOffset += info.nTorsions;
+    excludeOffset += molInfo.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 );
+    the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
+    the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
+    the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
+    the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
 
 
-    the_molecules[i].initialize( info );
+    the_molecules[i].initialize( molInfo );
 
 
-    atomOffset += info.nAtoms;
+    atomOffset += molInfo.nAtoms;
     delete[] theBonds;
     delete[] theBends;
     delete[] theTorsions;
@@ -618,9 +398,9 @@ void SimSetup::initFromBass( void ){
     have_extra =1;
 
     n_cells = (int)temp3 - 1;
-    cellx = info->boxLx / temp3;
-    celly = info->boxLy / temp3;
-    cellz = info->boxLz / temp3;
+    cellx = info->boxL[0] / temp3;
+    celly = info->boxL[1] / temp3;
+    cellz = info->boxL[2] / 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 );
@@ -635,9 +415,9 @@ void SimSetup::initFromBass( void ){
   }
   else{
     n_cells = (int)temp3;
-    cellx = info->boxLx / temp3;
-    celly = info->boxLy / temp3;
-    cellz = info->boxLz / temp3;
+    cellx = info->boxL[0] / temp3;
+    celly = info->boxL[1] / temp3;
+    cellz = info->boxL[2] / temp3;
   }
 
   current_mol = 0;
@@ -784,24 +564,25 @@ void SimSetup::gatherInfo( void ){
 
 
 void SimSetup::gatherInfo( void ){
+  int i,j,k;
 
   ensembleCase = -1;
   ffCase = -1;
 
   // get the stamps and globals;
-  the_stamps = stamps;
-  the_globals = globals;
+  stamps = stamps;
+  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();
+  info->target_temp = globals->getTargetTemp();
+  info->dt = globals->getDt();
+  info->run_time = globals->getRunTime();
+  n_components = globals->getNComponents();
 
 
   // get the forceField
 
-  strcpy( force_field, the_globals->getForceField() );
+  strcpy( force_field, globals->getForceField() );
 
   if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
   else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
@@ -815,7 +596,7 @@ void SimSetup::gatherInfo( void ){
 
   // get the ensemble
 
-  strcpy( ensemble, the_globals->getEnsemble() );
+  strcpy( ensemble, globals->getEnsemble() );
 
   if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
   else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
@@ -838,17 +619,17 @@ void SimSetup::gatherInfo( void ){
 
   // get the mixing rule
 
-  strcpy( info->mixingRule, the_globals->getMixingRule() );
-  info->usePBC = the_globals->getPBC();
+  strcpy( info->mixingRule, globals->getMixingRule() );
+  info->usePBC = globals->getPBC();
 	
   
   // get the components and calculate the tot_nMol and indvidual n_mol
  
-  the_components = the_globals->getComponents();
+  the_components = globals->getComponents();
   components_nmol = new int[n_components];
 
 
-  if( !the_globals->haveNMol() ){
+  if( !globals->haveNMol() ){
     // we don't have the total number of molecules, so we assume it is
     // given in each component
 
@@ -881,43 +662,43 @@ void SimSetup::gatherInfo( void ){
 
   // set the status, sample, and thermal kick times
   
-  if( the_globals->haveSampleTime() ){
-    info->sampleTime = the_globals->getSampleTime();
+  if( globals->haveSampleTime() ){
+    info->sampleTime = globals->getSampleTime();
     info->statusTime = info->sampleTime;
     info->thermalTime = info->sampleTime;
   }
   else{
-    info->sampleTime = the_globals->getRunTime();
+    info->sampleTime = globals->getRunTime();
     info->statusTime = info->sampleTime;
     info->thermalTime = info->sampleTime;
   }
 
-  if( the_globals->haveStatusTime() ){
-    info->statusTime = the_globals->getStatusTime();
+  if( globals->haveStatusTime() ){
+    info->statusTime = globals->getStatusTime();
   }
 
-  if( the_globals->haveThermalTime() ){
-    info->thermalTime = the_globals->getThermalTime();
+  if( globals->haveThermalTime() ){
+    info->thermalTime = globals->getThermalTime();
   }
 
   // check for the temperature set flag
 
-  if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
+  if( globals->haveTempSet() ) info->setTemp = 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();
+  if( globals->haveBox() ){
+    boxVector[0] = globals->getBox();
+    boxVector[1] = globals->getBox();
+    boxVector[2] = globals->getBox();
     
     info->setBox( boxVector );
   }
-  else if( the_globals->haveDensity() ){
+  else if( globals->haveDensity() ){
 
     double vol;
-    vol = (double)tot_nmol / the_globals->getDensity();
+    vol = (double)tot_nmol / globals->getDensity();
      boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
      boxVector[1] = boxVector[0];
      boxVector[2] = boxVector[0];
@@ -925,29 +706,29 @@ void SimSetup::gatherInfo( void ){
     info->setBox( boxVector );
   }
   else{
-    if( !the_globals->haveBoxX() ){
+    if( !globals->haveBoxX() ){
       sprintf( painCave.errMsg,
 	       "SimSetup error, no periodic BoxX size given.\n" );
       painCave.isFatal = 1;
       simError();
     }
-    boxVector[0] = the_globals->getBoxX();
+    boxVector[0] = globals->getBoxX();
 
-    if( !the_globals->haveBoxY() ){
+    if( !globals->haveBoxY() ){
       sprintf( painCave.errMsg,
 	       "SimSetup error, no periodic BoxY size given.\n" );
       painCave.isFatal = 1;
       simError();
     }
-    boxVector[1] = the_globals->getBoxY();
+    boxVector[1] = globals->getBoxY();
 
-    if( !the_globals->haveBoxZ() ){
+    if( !globals->haveBoxZ() ){
       sprintf( painCave.errMsg,
 	       "SimSetup error, no periodic BoxZ size given.\n" );
       painCave.isFatal = 1;
       simError();
     }
-    boxVector[2] = the_globals->getBoxZ();
+    boxVector[2] = globals->getBoxZ();
 
     info->setBox( boxVector );
   }
@@ -977,15 +758,16 @@ void SimSetup::finalInfoCheck( void ){
   }
   
 #ifdef IS_MPI
-  int myUse = usesDipoles
-  MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
+  int myUse = usesDipoles;
+  MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
 #endif //is_mpi
 
+  double theEcr, theEst;
 
-  if (the_globals->getUseRF() ) {
+  if (globals->getUseRF() ) {
     info->useReactionField = 1;
     
-    if( !the_globals->haveECR() ){
+    if( !globals->haveECR() ){
       sprintf( painCave.errMsg,
      	       "SimSetup Warning: using default value of 1/2 the smallest "
                "box length for the electrostaticCutoffRadius.\n"
@@ -993,27 +775,29 @@ void SimSetup::finalInfoCheck( void ){
       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;
+      smallest = info->boxL[0];
+      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
+      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
+      theEcr = 0.5 * smallest;
     } else {
-      info->ecr        = the_globals->getECR();
+      theEcr = globals->getECR();
     }
 
-    if( !the_globals->haveEST() ){
+    if( !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;
+      theEst = 0.05 * theEcr;
     } else {
-      info->est        = the_globals->getEST();
+      theEst= globals->getEST();
     }
+
+    info->setEcr( theEcr, theEst );
     
-    if(!the_globals->haveDielectric() ){
+    if(!globals->haveDielectric() ){
       sprintf( painCave.errMsg,
      	       "SimSetup Error: You are trying to use Reaction Field without"
                "setting a dielectric constant!\n"
@@ -1021,39 +805,41 @@ void SimSetup::finalInfoCheck( void ){
       painCave.isFatal = 1;
       simError();
     }
-    info->dielectric = the_globals->getDielectric();  
+    info->dielectric = 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;
+      if( !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->boxL[0];
+	if (info->boxL[1] <= smallest) smallest = info->boxL[1];
+	if (info->boxL[2] <= smallest) smallest = info->boxL[2];
+	theEcr = 0.5 * smallest;
       } else {
-        info->ecr        = the_globals->getECR();
+	theEcr = 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;
+      if( !globals->haveEST() ){
+	sprintf( painCave.errMsg,
+		 "SimSetup Warning: using default value of 0.05 * the "
+		 "electrostaticCutoffRadius for the "
+		 "electrostaticSkinThickness\n"
+		 );
+	painCave.isFatal = 0;
+	simError();
+	theEst = 0.05 * theEcr;
       } else {
-        info->est        = the_globals->getEST();
+	theEst= globals->getEST();
       }
+
+      info->setEcr( theEcr, theEst );
     }
   }  
 
@@ -1066,13 +852,13 @@ void SimSetup::initSystemCoords( void ){
 
 void SimSetup::initSystemCoords( void ){
 
- if( the_globals->haveInitialConfig() ){
+ if( globals->haveInitialConfig() ){
  
      InitializeFromFile* fileInit;
 #ifdef IS_MPI // is_mpi
      if( worldRank == 0 ){
 #endif //is_mpi
-   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
+   fileInit = new InitializeFromFile( globals->getInitialConfig() );
 #ifdef IS_MPI
      }else fileInit = new InitializeFromFile( NULL );
 #endif
@@ -1113,8 +899,8 @@ void SimSetup::makeOutNames( void ){
   if( worldRank == 0 ){
 #endif // is_mpi
     
-    if( the_globals->haveFinalConfig() ){
-      strcpy( info->finalName, the_globals->getFinalConfig() );
+    if( globals->haveFinalConfig() ){
+      strcpy( info->finalName, globals->getFinalConfig() );
     }
     else{
       strcpy( info->finalName, inFileName );
@@ -1197,6 +983,8 @@ void SimSetup::sysObjectsCreation( void ){
 
 void SimSetup::sysObjectsCreation( void ){
 
+  int i;
+
   // create the forceField
 
   createFF();
@@ -1219,7 +1007,14 @@ void SimSetup::sysObjectsCreation( void ){
   
   makeSysArrays();
 
+  // make and initialize the molecules (all but atomic coordinates)
   
+  makeMolecules();
+  info->identArray = new int[info->n_atoms];
+  for(i=0; i<info->n_atoms; i++){
+    info->identArray[i] = the_atoms[i]->getIdent();
+  }
+  
 
 
 }
@@ -1253,6 +1048,8 @@ void SimSetup::compList( void ){
 
 
 void SimSetup::compList( void ){
+
+  int i;
 
   comp_stamps = new MoleculeStamp*[n_components];
 
@@ -1279,7 +1076,7 @@ void SimSetup::compList( void ){
       
       // extract the component from the list;
       
-      currentStamp = the_stamps->extractMolStamp( id );
+      currentStamp = stamps->extractMolStamp( id );
       if( currentStamp == NULL ){
 	sprintf( painCave.errMsg,
 		 "SimSetup error: Component \"%s\" was not found in the "
@@ -1303,7 +1100,9 @@ void SimSetup::calcSysValues( void ){
 }
 
 void SimSetup::calcSysValues( void ){
+  int i, j, k;
 
+
   tot_atoms = 0;
   tot_bonds = 0;
   tot_bends = 0;
@@ -1333,6 +1132,7 @@ void SimSetup::mpiMolDivide( void ){
 
 void SimSetup::mpiMolDivide( void ){
   
+  int i, j, k;
   int localMol, allMol;
   int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
 
@@ -1402,7 +1202,9 @@ void SimSetup::makeSysArrays( void ){
 
 
 void SimSetup::makeSysArrays( void ){
+  int i, j, k;
 
+
   // create the atom and short range interaction arrays
 
   Atom::createArrays(info->n_atoms);
@@ -1476,3 +1278,209 @@ void SimSetup::makeSysArrays( void ){
   the_ff->setSimInfo( info );
 
 }
+
+void SimSetup::makeIntegrator( void ){
+
+  NVT*  myNVT = NULL;
+  NPTi* myNPTi = NULL;
+  NPTf* myNPTf = NULL;
+  NPTim* myNPTim = NULL;
+  NPTfm* myNPTfm = NULL;
+
+  switch( ensembleCase ){
+
+  case NVE_ENS:
+    new NVE( info, the_ff );
+    break;
+
+  case NVT_ENS:
+    myNVT = new NVT( info, the_ff );
+    myNVT->setTargetTemp(globals->getTargetTemp());
+
+    if (globals->haveTauThermostat()) 
+      myNVT->setTauThermostat(globals->getTauThermostat());
+
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use the NVT\n"
+               "    ensemble, you must set tauThermostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    break;
+
+  case NPTi_ENS:
+    myNPTi = new NPTi( info, the_ff );
+    myNPTi->setTargetTemp( globals->getTargetTemp() );
+
+    if (globals->haveTargetPressure())
+      myNPTi->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+               "    ensemble, you must set targetPressure in the BASS file.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    
+    if( globals->haveTauThermostat() )
+      myNPTi->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauThermostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+
+    if( globals->haveTauBarostat() )
+      myNPTi->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauBarostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    break;
+
+  case NPTf_ENS:
+    myNPTf = new NPTf( info, the_ff );
+    myNPTf->setTargetTemp( globals->getTargetTemp());
+
+    if (globals->haveTargetPressure())
+      myNPTf->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+               "    ensemble, you must set targetPressure in the BASS file.\n");
+      painCave.isFatal = 1;
+      simError();
+    }    
+
+    if( globals->haveTauThermostat() )
+      myNPTf->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauThermostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+
+    if( globals->haveTauBarostat() )
+      myNPTf->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauBarostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    break;
+    
+  case NPTim_ENS:
+    myNPTim = new NPTim( info, the_ff );
+    myNPTim->setTargetTemp( globals->getTargetTemp());
+
+    if (globals->haveTargetPressure())
+      myNPTim->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+               "    ensemble, you must set targetPressure in the BASS file.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    
+    if( globals->haveTauThermostat() )
+      myNPTim->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauThermostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+
+    if( globals->haveTauBarostat() )
+      myNPTim->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauBarostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    break;
+
+  case NPTfm_ENS:
+    myNPTfm = new NPTfm( info, the_ff );
+    myNPTfm->setTargetTemp( globals->getTargetTemp());
+
+    if (globals->haveTargetPressure())
+      myNPTfm->setTargetPressure(globals->getTargetPressure());
+    else {
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use a constant pressure\n"
+               "    ensemble, you must set targetPressure in the BASS file.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    
+    if( globals->haveTauThermostat() )
+      myNPTfm->setTauThermostat( globals->getTauThermostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauThermostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+
+    if( globals->haveTauBarostat() )
+      myNPTfm->setTauBarostat( globals->getTauBarostat() );
+    else{
+      sprintf( painCave.errMsg,
+               "SimSetup error: If you use an NPT\n"
+               "    ensemble, you must set tauBarostat.\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    break;
+
+  default:
+    sprintf( painCave.errMsg,
+	     "SimSetup Error. Unrecognized ensemble in case statement.\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+
+}
+
+void SimSetup::initFortran( void ){
+
+  info->refreshSim();
+  
+  if( !strcmp( info->mixingRule, "standard") ){
+    the_ff->initForceField( LB_MIXING_RULE );
+  }
+  else if( !strcmp( info->mixingRule, "explicit") ){
+    the_ff->initForceField( EXPLICIT_MIXING_RULE );
+  }
+  else{
+    sprintf( painCave.errMsg,
+	     "SimSetup Error: unknown mixing rule -> \"%s\"\n",
+	     info->mixingRule );
+    painCave.isFatal = 1;
+    simError();
+  }
+
+
+#ifdef IS_MPI
+  strcpy( checkPointMsg, 
+	  "Successfully intialized the mixingRule for Fortran." );
+  MPIcheckPoint();
+#endif // is_mpi
+
+}