--- trunk/mdtools/interface_implementation/SimSetup.cpp	2002/12/06 21:20:42	198
+++ trunk/mdtools/interface_implementation/SimSetup.cpp	2002/12/12 21:21:59	206
@@ -73,7 +73,7 @@ void SimSetup::createSim( void ){
 
   MakeStamps *the_stamps;
   Globals* the_globals;
-  int i;
+  int i, j;
 
   // get the stamps and globals;
   the_stamps = stamps;
@@ -213,10 +213,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();
   }
 
@@ -234,10 +234,62 @@ void SimSetup::createSim( void ){
 
   // divide the molecules among processors here.
   
-  new mpiSimulation( simnfo );
+  mpiSimulation* mpiSim = new mpiSimulation( simnfo );
   
-  simnfo->mpiSim->divideLabor( n_components, comp_stamps, components_nmol );
+  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
   
 
@@ -353,15 +405,21 @@ void SimSetup::createSim( void ){
 
 
 
-//   if( the_globals->haveInitialConfig() ){
-//        InitializeFromFile* fileInit;
-//     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
-
-//     fileInit->read_xyz( simnfo ); // default velocities on
+ 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{
+   delete fileInit;
+ }
+ else{
 
 #ifdef IS_MPI 
 
@@ -376,8 +434,10 @@ void SimSetup::createSim( void ){
 
   initFromBass();
 
-#endif // is_mpi
- 
+
+#endif
+ }
+
 #ifdef IS_MPI
   strcpy( checkPointMsg, "Successfully read in the initial configuration" );
   MPIcheckPoint();
@@ -387,7 +447,7 @@ void SimSetup::createSim( void ){
    
  
   
-  //   }
+
   
 #ifdef IS_MPI
   if( worldRank == 0 ){
@@ -513,57 +573,68 @@ 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( simnfo->mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= simnfo->mpiSim->getMyMolEnd() ){
+#endif // is_mpi	
 
-	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++;
-      }
+	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++;
-
     }
   }
 
@@ -572,31 +643,43 @@ void SimSetup::makeBonds( void ){
 
 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( simnfo->mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= simnfo->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 );
@@ -604,31 +687,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( simnfo->mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= simnfo->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++;
     }
   }
 
@@ -637,17 +732,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( simnfo->mpiSim->getMyMolStart() <= molIndex &&
+	  molIndex <= simnfo->mpiSim->getMyMolEnd() ){
+#endif // is_mpi	
+
       for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
 
 	current_torsion = comp_stamps[i]->getTorsion( k );
@@ -663,6 +764,12 @@ void SimSetup::makeTorsions( void ){
 	index++;
       }
       offset += comp_stamps[i]->getNAtoms();
+
+#ifdef IS_MPI
+      }
+#endif //is_mpi      
+
+      molIndex++;
     }
   }