--- trunk/mdtools/interface_implementation/SimSetup.cpp	2002/10/09 22:29:40	131
+++ trunk/mdtools/interface_implementation/SimSetup.cpp	2003/01/30 15:20:21	253
@@ -6,14 +6,22 @@
 #include "parse_me.h"
 #include "LRI.hpp"
 #include "Integrator.hpp"
+#include "simError.h"
+
 #ifdef IS_MPI
-#include "mpiInterface.h"
+#include "mpiBASS.h"
+#include "mpiSimulation.hpp"
 #include "bassDiag.hpp"
 #endif
 
 SimSetup::SimSetup(){
   stamps = new MakeStamps();
   globals = new Globals();
+  
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "SimSetup creation successful" );
+  MPIcheckPoint();
+#endif // IS_MPI
 }
 
 SimSetup::~SimSetup(){
@@ -23,14 +31,23 @@ void SimSetup::parseFile( char* fileName ){
 
 void SimSetup::parseFile( char* fileName ){
 
-  inFileName = fileName;
-  set_interface_stamps( stamps, globals );
 #ifdef IS_MPI
-  mpiEventInit();
+  if( worldRank == 0 ){
+#endif // is_mpi
+    
+    inFileName = fileName;
+    set_interface_stamps( stamps, globals );
+    
+#ifdef IS_MPI
+    mpiEventInit();
 #endif
-  yacc_BASS( fileName );
+
+    yacc_BASS( fileName );
+
 #ifdef IS_MPI
-  throwMPIEvent(NULL);
+    throwMPIEvent(NULL);
+  }
+  else receiveParse();
 #endif
 
 }
@@ -40,6 +57,7 @@ void SimSetup::receiveParse(void){
 
     set_interface_stamps( stamps, globals );
     mpiEventInit();
+    MPIcheckPoint();
     mpiEventLoop();
 
 }
@@ -51,11 +69,12 @@ void SimSetup::createSim( void ){
   delete dumpMe;
 }
 #endif
+
 void SimSetup::createSim( void ){
 
   MakeStamps *the_stamps;
   Globals* the_globals;
-  int i;
+  int i, j;
 
   // get the stamps and globals;
   the_stamps = stamps;
@@ -74,12 +93,20 @@ void SimSetup::createSim( void ){
   if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
   else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
   else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
+  else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
   else{
-    std::cerr<< "SimSetup Error. Unrecognized force field -> "
-	     << force_field << "\n";
-    exit(8);
+    sprintf( painCave.errMsg,
+	     "SimSetup Error. Unrecognized force field -> %s\n",
+	     force_field );
+    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];
@@ -94,9 +121,11 @@ void SimSetup::createSim( void ){
 
       if( !the_components[i]->haveNMol() ){
 	// we have a problem
-	std::cerr << "SimSetup Error. No global NMol or component NMol"
-		  << " given. Cannot calculate the number of atoms.\n";
-	exit( 8 );
+	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();
@@ -104,32 +133,81 @@ void SimSetup::createSim( void ){
     }
   }
   else{
-    std::cerr << "NOT A SUPPORTED FEATURE\n";
-
-//     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
-
+    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
+    
   }
 
+#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++ ){
 
-    comp_stamps[i] =
-      the_stamps->getMolecule( the_components[i]->getType() );
+    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
+  
 
 
+
   // caclulate the number of atoms, bonds, bends and torsions
 
   tot_atoms = 0;
@@ -137,10 +215,10 @@ void SimSetup::createSim( void ){
   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_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();
   }
 
@@ -151,16 +229,82 @@ void SimSetup::createSim( void ){
   simnfo->n_bends = tot_bends;
   simnfo->n_torsions = tot_torsions;
   simnfo->n_SRI = tot_SRI;
+  simnfo->n_mol = tot_nmol;
+
+  
+#ifdef IS_MPI
+
+  // divide the molecules among processors here.
+  
+  mpiSim = new mpiSimulation( simnfo );
+  
+  globalIndex = mpiSim->divideLabor();
+
+  // set up the local variables 
+  
+  int localMol, allMol;
+  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
+  
+  allMol = 0;
+  localMol = 0;
+  local_atoms = 0;
+  local_bonds = 0;
+  local_bends = 0;
+  local_torsions = 0;
+  for( i=0; i<n_components; i++ ){
+
+    for( j=0; j<components_nmol[i]; j++ ){
+      
+      if( mpiSim->getMyMolStart() <= allMol &&
+	  allMol <= mpiSim->getMyMolEnd() ){
+	
+	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++;
+      }      
+      allMol++;
+    }
+  }
+  local_SRI = local_bonds + local_bends + local_torsions;
+  
+
+  simnfo->n_atoms = mpiSim->getMyNlocal();  
+  
+  if( local_atoms != simnfo->n_atoms ){
+    sprintf( painCave.errMsg,
+	     "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
+	     " localAtom (%d) are note equal.\n",
+	     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
 
-  the_atoms = new Atom*[tot_atoms];
-  the_molecules = new Molecule[tot_nmol];
+  Atom::createArrays(simnfo->n_atoms);
+  the_atoms = new Atom*[simnfo->n_atoms];
+  the_molecules = new Molecule[simnfo->n_mol];
 
 
-  if( tot_SRI ){
-    the_sris = new SRI*[tot_SRI];
-    the_excludes = new ex_pair[tot_SRI];
+  if( simnfo->n_SRI ){
+    the_sris = new SRI*[simnfo->n_SRI];
+    the_excludes = new ex_pair[simnfo->n_SRI];
   }
 
   // set the arrays into the SimInfo object
@@ -171,44 +315,33 @@ void SimSetup::createSim( void ){
   simnfo->excludes = the_excludes;
 
 
-  // initialize the arrays
-
-  the_ff->setSimInfo( simnfo );
-
-  makeAtoms();
-
-  if( tot_bonds ){
-    makeBonds();
-  }
-
-  if( tot_bends ){
-    makeBends();
-  }
-
-  if( tot_torsions ){
-    makeTorsions();
-  }
-
-  //  makeMolecules();
-
   // get some of the tricky things that may still be in the globals
 
   if( simnfo->n_dipoles ){
 
     if( !the_globals->haveRRF() ){
-      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
-      exit(8);
+      sprintf( painCave.errMsg,
+	       "SimSetup Error, system has dipoles, but no rRF was set.\n");
+      painCave.isFatal = 1;
+      simError();
     }
     if( !the_globals->haveDielectric() ){
-      std::cerr << "SimSetup Error, system has dipoles, but no"
-		<< " dielectric was set.\n";
-      exit(8);
+      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();
     simnfo->box_y = the_globals->getBox();
@@ -224,118 +357,188 @@ void SimSetup::createSim( void ){
   }
   else{
     if( !the_globals->haveBoxX() ){
-      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
-      exit(8);
+      sprintf( painCave.errMsg,
+	       "SimSetup error, no periodic BoxX size given.\n" );
+      painCave.isFatal = 1;
+      simError();
     }
     simnfo->box_x = the_globals->getBoxX();
 
     if( !the_globals->haveBoxY() ){
-      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
-      exit(8);
+      sprintf( painCave.errMsg,
+	       "SimSetup error, no periodic BoxY size given.\n" );
+      painCave.isFatal = 1;
+      simError();
     }
     simnfo->box_y = the_globals->getBoxY();
 
     if( !the_globals->haveBoxZ() ){
-      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
-      exit(8);
+      sprintf( painCave.errMsg,
+	       "SimSetup error, no periodic BoxZ size given.\n" );
+      painCave.isFatal = 1;
+      simError();
     }
     simnfo->box_z = the_globals->getBoxZ();
   }
 
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "Box size set up" );
+  MPIcheckPoint();
+#endif // is_mpi
 
-//   if( the_globals->haveInitialConfig() ){
-//        InitializeFromFile* fileInit;
-//     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
 
-//     fileInit->read_xyz( simnfo ); // default velocities on
+  // initialize the arrays
 
-//     delete fileInit;
-//   }
-//   else{
+  the_ff->setSimInfo( simnfo );
 
-    initFromBass();
+  makeAtoms();
 
+  if( tot_bonds ){
+    makeBonds();
+  }
 
-//   }
+  if( tot_bends ){
+    makeBends();
+  }
 
-//   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" );
-//       }
-//     }
-//   }
+  if( tot_torsions ){
+    makeTorsions();
+  }
 
-//   // 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" );
-//     }
-//   }
 
 
-  // set the status, sample, and themal kick times
 
+ 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
+
+
+   
+ 
+  
+
+  
+#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;
@@ -360,13 +563,15 @@ void SimSetup::createSim( void ){
   if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
 
 
-  // make the longe range forces and the integrator
+//   // make the longe range forces and the integrator
 
-  new AllLong( simnfo );
+//   new AllLong( simnfo );
 
-  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
+  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
   if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
   if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
+  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
+
 }
 
 void SimSetup::makeAtoms( void ){
@@ -375,90 +580,120 @@ void SimSetup::makeAtoms( void ){
   double ux, uy, uz, uSqr, u;
   AtomStamp* current_atom;
   DirectionalAtom* dAtom;
-  int molIndex, molStart, molEnd, nMemb;
+  int molIndex, molStart, molEnd, nMemb, lMolIndex;
 
-
+  lMolIndex = 0;
   molIndex = 0;
   index = 0;
   for( i=0; i<n_components; i++ ){
 
     for( j=0; j<components_nmol[i]; j++ ){
 
-      molStart = index;
-      nMemb = comp_stamps[i]->getNAtoms();
-      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
+#ifdef IS_MPI
+      if( mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= mpiSim->getMyMolEnd() ){
+#endif // is_mpi	
 
-	current_atom = comp_stamps[i]->getAtom( k );
-	if( current_atom->haveOrientation() ){
-
-	  dAtom = new DirectionalAtom;
-	  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;
-	}
-	the_atoms[index]->setType( current_atom->getType() );
-	the_atoms[index]->setIndex( index );
-
-	// increment the index and repeat;
-	index++;
-      }
+	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++;
 
-      molEnd = index -1;
-      the_molecules[molIndex].setNMembers( nMemb );
-      the_molecules[molIndex].setStartAtom( molStart );
-      the_molecules[molIndex].setEndAtom( molEnd );
+#ifdef IS_MPI
+      }
+#endif //is_mpi
+      
       molIndex++;
-
     }
   }
 
+#ifdef IS_MPI
+    for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
+    
+    delete[] globalIndex;
+#endif IS_MPI
+	  
   the_ff->initializeAtoms();
 }
 
 void SimSetup::makeBonds( void ){
 
-  int i, j, k, index, offset;
+  int i, j, k, index, offset, molIndex;
   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++ ){
 
-      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;
-
-	the_excludes[index].i = the_bonds[index].a;
-	the_excludes[index].j = the_bonds[index].b;
-
-	// increment the index and repeat;
-	index++;
+#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;
+	  
+	  the_excludes[index].i = the_bonds[index].a;
+	  the_excludes[index].j = the_bonds[index].b;
+	  
+	  // increment the index and repeat;
+	  index++;
+	}
+	offset += comp_stamps[i]->getNAtoms();
+	
+#ifdef IS_MPI
       }
-      offset += comp_stamps[i]->getNAtoms();
-    }
+#endif is_mpi
+      
+      molIndex++;
+    }      
   }
 
   the_ff->initializeBonds( the_bonds );
@@ -466,31 +701,43 @@ void SimSetup::makeBends( void ){
 
 void SimSetup::makeBends( void ){
 
-  int i, j, k, index, offset;
+  int i, j, k, index, offset, molIndex;
   bend_set* the_bends;
   BendStamp* current_bend;
 
   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++ ){
 
-      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
+#ifdef IS_MPI
+      if( mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= mpiSim->getMyMolEnd() ){
+#endif // is_mpi	
 
-	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;
-
-	the_excludes[index + tot_bonds].i = the_bends[index].a;
-	the_excludes[index + tot_bonds].j = the_bends[index].c;
-
-	// increment the index and repeat;
-	index++;
+	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
+	  
+	  current_bend = comp_stamps[i]->getBend( k );
+	  the_bends[index].a = current_bend->getA() + offset;
+	  the_bends[index].b = current_bend->getB() + offset;
+	  the_bends[index].c = current_bend->getC() + offset;
+	  
+	  the_excludes[index + tot_bonds].i = the_bends[index].a;
+	  the_excludes[index + tot_bonds].j = the_bends[index].c;
+	  
+	  // increment the index and repeat;
+	  index++;
+	}
+	offset += comp_stamps[i]->getNAtoms();
+	
+#ifdef IS_MPI
       }
-      offset += comp_stamps[i]->getNAtoms();
+#endif //is_mpi
+
+      molIndex++;
     }
   }
 
@@ -499,17 +746,23 @@ void SimSetup::makeTorsions( void ){
 
 void SimSetup::makeTorsions( void ){
 
-  int i, j, k, index, offset;
+  int i, j, k, index, offset, molIndex;
   torsion_set* the_torsions;
   TorsionStamp* current_torsion;
 
   the_torsions = new torsion_set[tot_torsions];
   index = 0;
   offset = 0;
+  molIndex = 0;
   for( i=0; i<n_components; i++ ){
 
     for( j=0; j<components_nmol[i]; j++ ){
 
+#ifdef IS_MPI
+      if( mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= mpiSim->getMyMolEnd() ){
+#endif // is_mpi	
+
       for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
 
 	current_torsion = comp_stamps[i]->getTorsion( k );
@@ -525,6 +778,12 @@ void SimSetup::makeTorsions( void ){
 	index++;
       }
       offset += comp_stamps[i]->getNAtoms();
+
+#ifdef IS_MPI
+      }
+#endif //is_mpi      
+
+      molIndex++;
     }
   }
 
@@ -558,8 +817,11 @@ void SimSetup::initFromBass( void ){
     n_per_extra = (int)ceil( temp1 );
 
     if( n_per_extra > 4){
-      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
-      exit(8);
+      sprintf( painCave.errMsg,
+	       "SimSetup error. There has been an error in constructing"
+	       " the non-complete lattice.\n" );
+      painCave.isFatal = 1;
+      simError();
     }
   }
   else{
@@ -664,12 +926,15 @@ void SimSetup::makeElement( double x, double y, double
 
     current_atom = comp_stamps[current_comp]->getAtom( k );
     if( !current_atom->havePosition() ){
-      std::cerr << "Component " << comp_stamps[current_comp]->getID()
-		<< ", atom " << current_atom->getType()
-		<< " does not have a position specified.\n"
-		<< "The initialization routine is unable to give a start"
-		<< " position.\n";
-      exit(8);
+      sprintf( painCave.errMsg,
+	       "SimSetup:initFromBass error.\n"
+	       "\tComponent %s, atom %s does not have a position specified.\n"
+	       "\tThe initialization routine is unable to give a start"
+	       " position.\n",
+	       comp_stamps[current_comp]->getID(),
+	       current_atom->getType() );
+      painCave.isFatal = 1;
+      simError();
     }
 
     the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );