--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/08/15 19:24:13	699
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/08/20 14:34:04	701
@@ -189,39 +189,39 @@ void SimSetup::makeMolecules( void ){
       // make the Atoms
     
       for(j=0; j<molInfo.nAtoms; j++){
-	
-	currentAtom = comp_stamps[stampID]->getAtom( j );
-	if( currentAtom->haveOrientation() ){
-	  
-	  dAtom = new DirectionalAtom( (j + atomOffset),
-				       info[k].getConfiguration() );
-	  info[k].n_oriented++;
-	  molInfo.myAtoms[j] = dAtom;
-	  
-	  ux = currentAtom->getOrntX();
-	  uy = currentAtom->getOrntY();
-	  uz = currentAtom->getOrntZ();
-	  
-	  uSqr = (ux * ux) + (uy * uy) + (uz * uz);
-	  
-	  u = sqrt( uSqr );
-	  ux = ux / u;
-	  uy = uy / u;
-	  uz = uz / u;
-	  
-	  dAtom->setSUx( ux );
-	  dAtom->setSUy( uy );
-	  dAtom->setSUz( uz );
-	}
-	else{
-	  molInfo.myAtoms[j] = new GeneralAtom( (j + atomOffset),
-						info[k].getConfiguration() );
-	}
-	molInfo.myAtoms[j]->setType( currentAtom->getType() );
+  
+  currentAtom = comp_stamps[stampID]->getAtom( j );
+  if( currentAtom->haveOrientation() ){
     
+    dAtom = new DirectionalAtom( (j + atomOffset),
+               info[k].getConfiguration() );
+    info[k].n_oriented++;
+    molInfo.myAtoms[j] = dAtom;
+    
+    ux = currentAtom->getOrntX();
+    uy = currentAtom->getOrntY();
+    uz = currentAtom->getOrntZ();
+    
+    uSqr = (ux * ux) + (uy * uy) + (uz * uz);
+    
+    u = sqrt( uSqr );
+    ux = ux / u;
+    uy = uy / u;
+    uz = uz / u;
+    
+    dAtom->setSUx( ux );
+    dAtom->setSUy( uy );
+    dAtom->setSUz( uz );
+  }
+  else{
+    molInfo.myAtoms[j] = new GeneralAtom( (j + atomOffset),
+            info[k].getConfiguration() );
+  }
+  molInfo.myAtoms[j]->setType( currentAtom->getType() );
+    
 #ifdef IS_MPI
       
-	molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
+  molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
       
 #endif // is_mpi
       } 
@@ -229,145 +229,145 @@ void SimSetup::makeMolecules( void ){
     // make the bonds
       for(j=0; j<molInfo.nBonds; j++){
       
-	currentBond = comp_stamps[stampID]->getBond( j );
-	theBonds[j].a = currentBond->getA() + atomOffset;
-	theBonds[j].b = currentBond->getB() + atomOffset;
-	
-	exI = theBonds[j].a;
-	exJ = theBonds[j].b;
-	
-	// exclude_I must always be the smaller of the pair
-	if( exI > exJ ){
-	  tempEx = exI;
-	  exI = exJ;
-	  exJ = tempEx;
-	}
-#ifdef IS_MPI
-	tempEx = exI;
-	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-	tempEx = exJ;
-	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-	
-	info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
+  currentBond = comp_stamps[stampID]->getBond( j );
+  theBonds[j].a = currentBond->getA() + atomOffset;
+  theBonds[j].b = currentBond->getB() + atomOffset;
+  
+  exI = theBonds[j].a;
+  exJ = theBonds[j].b;
+  
+  // exclude_I must always be the smaller of the pair
+  if( exI > exJ ){
+    tempEx = exI;
+    exI = exJ;
+    exJ = tempEx;
+  }
+#ifdef IS_MPI
+  tempEx = exI;
+  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+  tempEx = exJ;
+  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+  
+  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
 #else  // isn't MPI
-	
-	info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
+  
+  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
       }
       excludeOffset += molInfo.nBonds;
       
       //make the bends
       for(j=0; j<molInfo.nBends; j++){
-	
-	currentBend = comp_stamps[stampID]->getBend( j );
-	theBends[j].a = currentBend->getA() + atomOffset;
-	theBends[j].b = currentBend->getB() + atomOffset;
-	theBends[j].c = currentBend->getC() + atomOffset;
-	
-	if( currentBend->haveExtras() ){
-	  
-	  extras = currentBend->getExtras();
-	  current_extra = extras;
-	  
-	  while( current_extra != NULL ){
-	    if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
-	      
-	      switch( current_extra->getType() ){
-		
-	      case 0:
-		theBends[j].ghost =
-		  current_extra->getInt() + atomOffset;
-		theBends[j].isGhost = 1;
-		break;
-		
-	      case 1:
-		theBends[j].ghost = 
-		  (int)current_extra->getDouble() + atomOffset;
-		theBends[j].isGhost = 1;
-		break;
-		
-	      default:
-		sprintf( painCave.errMsg,
-			 "SimSetup Error: ghostVectorSource was neither a "
-			 "double nor an int.\n"
-			 "-->Bend[%d] in %s\n",
-			 j, comp_stamps[stampID]->getID() );
-		painCave.isFatal = 1;
-		simError();
-	      }
-	    }
-	    
-	    else{
-	      
-	      sprintf( painCave.errMsg,
-		       "SimSetup Error: unhandled bend assignment:\n"
-		       "    -->%s in Bend[%d] in %s\n",
-		       current_extra->getlhs(),
-		       j, comp_stamps[stampID]->getID() );
-	      painCave.isFatal = 1;
-	      simError();
-	    }
-	    
-	    current_extra = current_extra->getNext();
-	  }
-	}
-	
-	if( !theBends[j].isGhost ){
-	  
-	  exI = theBends[j].a;
-	  exJ = theBends[j].c;
-	}
-	else{
-	  
-	  exI = theBends[j].a;
-	  exJ = theBends[j].b;
-	}
-	
-	// exclude_I must always be the smaller of the pair
-	if( exI > exJ ){
-	  tempEx = exI;
-	  exI = exJ;
-	  exJ = tempEx;
-	}
+  
+  currentBend = comp_stamps[stampID]->getBend( j );
+  theBends[j].a = currentBend->getA() + atomOffset;
+  theBends[j].b = currentBend->getB() + atomOffset;
+  theBends[j].c = currentBend->getC() + atomOffset;
+  
+  if( currentBend->haveExtras() ){
+    
+    extras = currentBend->getExtras();
+    current_extra = extras;
+    
+    while( current_extra != NULL ){
+      if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
+        
+        switch( current_extra->getType() ){
+    
+        case 0:
+    theBends[j].ghost =
+      current_extra->getInt() + atomOffset;
+    theBends[j].isGhost = 1;
+    break;
+    
+        case 1:
+    theBends[j].ghost = 
+      (int)current_extra->getDouble() + atomOffset;
+    theBends[j].isGhost = 1;
+    break;
+    
+        default:
+    sprintf( painCave.errMsg,
+       "SimSetup Error: ghostVectorSource was neither a "
+       "double nor an int.\n"
+       "-->Bend[%d] in %s\n",
+       j, comp_stamps[stampID]->getID() );
+    painCave.isFatal = 1;
+    simError();
+        }
+      }
+      
+      else{
+        
+        sprintf( painCave.errMsg,
+           "SimSetup Error: unhandled bend assignment:\n"
+           "    -->%s in Bend[%d] in %s\n",
+           current_extra->getlhs(),
+           j, comp_stamps[stampID]->getID() );
+        painCave.isFatal = 1;
+        simError();
+      }
+      
+      current_extra = current_extra->getNext();
+    }
+  }
+  
+  if( !theBends[j].isGhost ){
+    
+    exI = theBends[j].a;
+    exJ = theBends[j].c;
+  }
+  else{
+    
+    exI = theBends[j].a;
+    exJ = theBends[j].b;
+  }
+  
+  // exclude_I must always be the smaller of the pair
+  if( exI > exJ ){
+    tempEx = exI;
+    exI = exJ;
+    exJ = tempEx;
+  }
 #ifdef IS_MPI
-	tempEx = exI;
-	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-	tempEx = exJ;
-	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+  tempEx = exI;
+  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+  tempEx = exJ;
+  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
       
-	info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
+  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
 #else  // isn't MPI
-	info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
+  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
       }
       excludeOffset += molInfo.nBends;
       
       for(j=0; j<molInfo.nTorsions; j++){
-	
-	currentTorsion = comp_stamps[stampID]->getTorsion( j );
-	theTorsions[j].a = currentTorsion->getA() + atomOffset;
-	theTorsions[j].b = currentTorsion->getB() + atomOffset;
-	theTorsions[j].c = currentTorsion->getC() + atomOffset;
-	theTorsions[j].d = currentTorsion->getD() + atomOffset;
-	
-	exI = theTorsions[j].a;
-	exJ = theTorsions[j].d;
-	
-	// exclude_I must always be the smaller of the pair
-	if( exI > exJ ){
-	  tempEx = exI;
-	  exI = exJ;
-	  exJ = tempEx;
-	}
+  
+  currentTorsion = comp_stamps[stampID]->getTorsion( j );
+  theTorsions[j].a = currentTorsion->getA() + atomOffset;
+  theTorsions[j].b = currentTorsion->getB() + atomOffset;
+  theTorsions[j].c = currentTorsion->getC() + atomOffset;
+  theTorsions[j].d = currentTorsion->getD() + atomOffset;
+  
+  exI = theTorsions[j].a;
+  exJ = theTorsions[j].d;
+  
+  // exclude_I must always be the smaller of the pair
+  if( exI > exJ ){
+    tempEx = exI;
+    exI = exJ;
+    exJ = tempEx;
+  }
 #ifdef IS_MPI
-	tempEx = exI;
-	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-	tempEx = exJ;
-	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-	
-	info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
+  tempEx = exI;
+  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+  tempEx = exJ;
+  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+  
+  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
 #else  // isn't MPI
-	info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
+  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
 #endif  //is_mpi
       }
       excludeOffset += molInfo.nTorsions;
@@ -436,8 +436,8 @@ void SimSetup::initFromBass( void ){
 
     if( n_per_extra > 4){
       sprintf( painCave.errMsg,
-	       "SimSetup error. There has been an error in constructing"
-	       " the non-complete lattice.\n" );
+         "SimSetup error. There has been an error in constructing"
+         " the non-complete lattice.\n" );
       painCave.isFatal = 1;
       simError();
     }
@@ -458,21 +458,21 @@ void SimSetup::initFromBass( void ){
     for( j=0; j < n_cells; j++ ){
       for( k=0; k < n_cells; k++ ){
 
-	makeElement( i * cellx,
-		     j * celly,
-		     k * cellz );
+  makeElement( i * cellx,
+         j * celly,
+         k * cellz );
 
-	makeElement( i * cellx + 0.5 * cellx,
-		     j * celly + 0.5 * celly,
-		     k * cellz );
+  makeElement( i * cellx + 0.5 * cellx,
+         j * celly + 0.5 * celly,
+         k * cellz );
 
-	makeElement( i * cellx,
-		     j * celly + 0.5 * celly,
-		     k * cellz + 0.5 * cellz );
+  makeElement( i * cellx,
+         j * celly + 0.5 * celly,
+         k * cellz + 0.5 * cellz );
 
-	makeElement( i * cellx + 0.5 * cellx,
-		     j * celly,
-		     k * cellz + 0.5 * cellz );
+  makeElement( i * cellx + 0.5 * cellx,
+         j * celly,
+         k * cellz + 0.5 * cellz );
       }
     }
   }
@@ -484,43 +484,43 @@ void SimSetup::initFromBass( void ){
     for( i=0; i < (n_cells+1) && !done; i++ ){
       for( j=0; j < (n_cells+1) && !done; j++ ){
 
-	if( i < n_cells ){
+  if( i < n_cells ){
 
-	  if( j < n_cells ){
-	    start_ndx = n_cells;
-	  }
-	  else start_ndx = 0;
-	}
-	else start_ndx = 0;
+    if( j < n_cells ){
+      start_ndx = n_cells;
+    }
+    else start_ndx = 0;
+  }
+  else start_ndx = 0;
 
-	for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
+  for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
 
-	  makeElement( i * cellx,
-		       j * celly,
-		       k * cellz );
-	  done = ( current_mol >= tot_nmol );
+    makeElement( i * cellx,
+           j * celly,
+           k * cellz );
+    done = ( current_mol >= tot_nmol );
 
-	  if( !done && n_per_extra > 1 ){
-	    makeElement( i * cellx + 0.5 * cellx,
-			 j * celly + 0.5 * celly,
-			 k * cellz );
-	    done = ( current_mol >= tot_nmol );
-	  }
+    if( !done && n_per_extra > 1 ){
+      makeElement( i * cellx + 0.5 * cellx,
+       j * celly + 0.5 * celly,
+       k * cellz );
+      done = ( current_mol >= tot_nmol );
+    }
 
-	  if( !done && n_per_extra > 2){
-	    makeElement( i * cellx,
-			 j * celly + 0.5 * celly,
-			 k * cellz + 0.5 * cellz );
-	    done = ( current_mol >= tot_nmol );
-	  }
+    if( !done && n_per_extra > 2){
+      makeElement( i * cellx,
+       j * celly + 0.5 * celly,
+       k * cellz + 0.5 * cellz );
+      done = ( current_mol >= tot_nmol );
+    }
 
-	  if( !done && n_per_extra > 3){
-	    makeElement( i * cellx + 0.5 * cellx,
-			 j * celly,
-			 k * cellz + 0.5 * cellz );
-	    done = ( current_mol >= tot_nmol );
-	  }
-	}
+    if( !done && n_per_extra > 3){
+      makeElement( i * cellx + 0.5 * cellx,
+       j * celly,
+       k * cellz + 0.5 * cellz );
+      done = ( current_mol >= tot_nmol );
+    }
+  }
       }
     }
   }
@@ -543,12 +543,12 @@ void SimSetup::makeElement( double x, double y, double
     current_atom = comp_stamps[current_comp]->getAtom( k );
     if( !current_atom->havePosition() ){
       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() );
+         "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();
     }
@@ -617,8 +617,8 @@ void SimSetup::gatherInfo( void ){
   else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
   else{
     sprintf( painCave.errMsg,
-	     "SimSetup Error. Unrecognized force field -> %s\n",
-	     force_field );
+       "SimSetup Error. Unrecognized force field -> %s\n",
+       force_field );
     painCave.isFatal = 1;
     simError();
   }
@@ -636,9 +636,9 @@ void SimSetup::gatherInfo( void ){
   else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
   else{
     sprintf( painCave.errMsg,
-	     "SimSetup Warning. Unrecognized Ensemble -> %s, "
+       "SimSetup Warning. Unrecognized Ensemble -> %s, "
              "reverting to NVE for this simulation.\n",
-	     ensemble );
+       ensemble );
     painCave.isFatal = 0;
     simError();
     strcpy( ensemble, "NVE" );
@@ -669,12 +669,12 @@ void SimSetup::gatherInfo( void ){
     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();
+  // 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();
@@ -683,11 +683,11 @@ void SimSetup::gatherInfo( void ){
   }
   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" );
+       "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();
   }
@@ -741,26 +741,26 @@ void SimSetup::gatherInfo( void ){
   }
     else{
       if( !globals->haveBoxX() ){
-	sprintf( painCave.errMsg,
-		 "SimSetup error, no periodic BoxX size given.\n" );
-	painCave.isFatal = 1;
-	simError();
+  sprintf( painCave.errMsg,
+     "SimSetup error, no periodic BoxX size given.\n" );
+  painCave.isFatal = 1;
+  simError();
       }
       boxVector[0] = globals->getBoxX();
       
       if( !globals->haveBoxY() ){
-	sprintf( painCave.errMsg,
-		 "SimSetup error, no periodic BoxY size given.\n" );
-	painCave.isFatal = 1;
-	simError();
+  sprintf( painCave.errMsg,
+     "SimSetup error, no periodic BoxY size given.\n" );
+  painCave.isFatal = 1;
+  simError();
       }
       boxVector[1] = globals->getBoxY();
       
       if( !globals->haveBoxZ() ){
-	sprintf( painCave.errMsg,
-		 "SimSetup error, no periodic BoxZ size given.\n" );
-	painCave.isFatal = 1;
-	simError();
+  sprintf( painCave.errMsg,
+     "SimSetup error, no periodic BoxZ size given.\n" );
+  painCave.isFatal = 1;
+  simError();
       }
       boxVector[2] = globals->getBoxZ();
       
@@ -803,78 +803,78 @@ void SimSetup::finalInfoCheck( void ){
       info[i].useReactionField = 1;
       
       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[i].boxL[0];
-	if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
-	if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
-	theEcr = 0.5 * smallest;
+  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[i].boxL[0];
+  if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
+  if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
+  theEcr = 0.5 * smallest;
       } else {
-	theEcr = globals->getECR();
+  theEcr = globals->getECR();
       }
       
       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;
+  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 {
-	theEst= globals->getEST();
+  theEst= globals->getEST();
       }
       
       info[i].setEcr( theEcr, theEst );
       
       if(!globals->haveDielectric() ){
-	sprintf( painCave.errMsg,
-		 "SimSetup Error: You are trying to use Reaction Field without"
-		 "setting a dielectric constant!\n"
-		 );
-	painCave.isFatal = 1;
-	simError();
+  sprintf( painCave.errMsg,
+     "SimSetup Error: You are trying to use Reaction Field without"
+     "setting a dielectric constant!\n"
+     );
+  painCave.isFatal = 1;
+  simError();
       }
       info[i].dielectric = globals->getDielectric();  
     } 
     else {
       if (usesDipoles) {
-	
-	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[i].boxL[0];
-	  if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
-	  if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
-	  theEcr = 0.5 * smallest;
-	} else {
-	  theEcr = globals->getECR();
-	}
-	
-	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 {
-	  theEst= globals->getEST();
-	}
-	
-	info[i].setEcr( theEcr, theEst );
+  
+  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[i].boxL[0];
+    if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
+    if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
+    theEcr = 0.5 * smallest;
+  } else {
+    theEcr = globals->getECR();
+  }
+  
+  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 {
+    theEst= globals->getEST();
+  }
+  
+  info[i].setEcr( theEcr, theEst );
       }
     }  
   }
@@ -919,7 +919,7 @@ void SimSetup::initSystemCoords( void ){
     // no init from bass
     
     sprintf( painCave.errMsg,
-	     "Cannot intialize a parallel simulation without an initial configuration file.\n" );
+       "Cannot intialize a parallel simulation without an initial configuration file.\n" );
     painCave.isFatal;
     simError();
     
@@ -951,31 +951,31 @@ void SimSetup::makeOutNames( void ){
 #endif // is_mpi
       
       if( globals->haveFinalConfig() ){
-	strcpy( info[k].finalName, globals->getFinalConfig() );
+  strcpy( info[k].finalName, globals->getFinalConfig() );
       }
       else{
-	strcpy( info[k].finalName, inFileName );
-	char* endTest;
-	int nameLength = strlen( info[k].finalName );
-	endTest = &(info[k].finalName[nameLength - 5]);
-	if( !strcmp( endTest, ".bass" ) ){
-	  strcpy( endTest, ".eor" );
-	}
-	else if( !strcmp( endTest, ".BASS" ) ){
-	  strcpy( endTest, ".eor" );
-	}
-	else{
-	  endTest = &(info[k].finalName[nameLength - 4]);
-	  if( !strcmp( endTest, ".bss" ) ){
-	    strcpy( endTest, ".eor" );
-	  }
-	  else if( !strcmp( endTest, ".mdl" ) ){
-	    strcpy( endTest, ".eor" );
-	  }
-	  else{
-	    strcat( info[k].finalName, ".eor" );
-	  }
-	}
+  strcpy( info[k].finalName, inFileName );
+  char* endTest;
+  int nameLength = strlen( info[k].finalName );
+  endTest = &(info[k].finalName[nameLength - 5]);
+  if( !strcmp( endTest, ".bass" ) ){
+    strcpy( endTest, ".eor" );
+  }
+  else if( !strcmp( endTest, ".BASS" ) ){
+    strcpy( endTest, ".eor" );
+  }
+  else{
+    endTest = &(info[k].finalName[nameLength - 4]);
+    if( !strcmp( endTest, ".bss" ) ){
+      strcpy( endTest, ".eor" );
+    }
+    else if( !strcmp( endTest, ".mdl" ) ){
+      strcpy( endTest, ".eor" );
+    }
+    else{
+      strcat( info[k].finalName, ".eor" );
+    }
+  }
       }
       
       // make the sample and status out names
@@ -985,44 +985,44 @@ void SimSetup::makeOutNames( void ){
       int nameLength = strlen( info[k].sampleName );
       endTest = &(info[k].sampleName[nameLength - 5]);
       if( !strcmp( endTest, ".bass" ) ){
-	strcpy( endTest, ".dump" );
+  strcpy( endTest, ".dump" );
       }
       else if( !strcmp( endTest, ".BASS" ) ){
-	strcpy( endTest, ".dump" );
+  strcpy( endTest, ".dump" );
       }
       else{
-	endTest = &(info[k].sampleName[nameLength - 4]);
-	if( !strcmp( endTest, ".bss" ) ){
-	  strcpy( endTest, ".dump" );
-	}
-	else if( !strcmp( endTest, ".mdl" ) ){
-	  strcpy( endTest, ".dump" );
-	}
-	else{
-	  strcat( info[k].sampleName, ".dump" );
-	}
+  endTest = &(info[k].sampleName[nameLength - 4]);
+  if( !strcmp( endTest, ".bss" ) ){
+    strcpy( endTest, ".dump" );
+  }
+  else if( !strcmp( endTest, ".mdl" ) ){
+    strcpy( endTest, ".dump" );
+  }
+  else{
+    strcat( info[k].sampleName, ".dump" );
+  }
       }
       
       strcpy( info[k].statusName, inFileName );
       nameLength = strlen( info[k].statusName );
       endTest = &(info[k].statusName[nameLength - 5]);
       if( !strcmp( endTest, ".bass" ) ){
-	strcpy( endTest, ".stat" );
+  strcpy( endTest, ".stat" );
       }
       else if( !strcmp( endTest, ".BASS" ) ){
-	strcpy( endTest, ".stat" );
+  strcpy( endTest, ".stat" );
       }
       else{
-	endTest = &(info[k].statusName[nameLength - 4]);
-	if( !strcmp( endTest, ".bss" ) ){
-	  strcpy( endTest, ".stat" );
-	}
-	else if( !strcmp( endTest, ".mdl" ) ){
-	  strcpy( endTest, ".stat" );
-	}
-	else{
-	  strcat( info[k].statusName, ".stat" );
-	}
+  endTest = &(info[k].statusName[nameLength - 4]);
+  if( !strcmp( endTest, ".bss" ) ){
+    strcpy( endTest, ".stat" );
+  }
+  else if( !strcmp( endTest, ".mdl" ) ){
+    strcpy( endTest, ".stat" );
+  }
+  else{
+    strcat( info[k].statusName, ".stat" );
+  }
       }
       
 #ifdef IS_MPI
@@ -1089,7 +1089,7 @@ void SimSetup::createFF( void ){
 
   default:
     sprintf( painCave.errMsg,
-	     "SimSetup Error. Unrecognized force field in case statement.\n");
+       "SimSetup Error. Unrecognized force field in case statement.\n");
     painCave.isFatal = 1;
     simError();
   }
@@ -1135,12 +1135,12 @@ void SimSetup::compList( void ){
       
       currentStamp = 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();
+  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 );
@@ -1219,12 +1219,12 @@ void SimSetup::mpiMolDivide( void ){
     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++;
+  
+  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[0].molMembershipArray[globalAtomIndex] = allMol;
@@ -1240,10 +1240,10 @@ void SimSetup::mpiMolDivide( void ){
   
   if( local_atoms != info[0].n_atoms ){
     sprintf( painCave.errMsg,
-	     "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
-	     " localAtom (%d) are not equal.\n",
-	     info[0].n_atoms,
-	     local_atoms );
+       "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
+       " localAtom (%d) are not equal.\n",
+       info[0].n_atoms,
+       local_atoms );
     painCave.isFatal = 1;
     simError();
   }
@@ -1286,10 +1286,10 @@ void SimSetup::makeSysArrays( void ){
     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++;
+  the_molecules[molIndex].setStampID( molCompType[i] );
+  the_molecules[molIndex].setMyIndex( molIndex );
+  the_molecules[molIndex].setGlobalIndex( i );
+  molIndex++;
       }
     }
     
@@ -1299,14 +1299,14 @@ void SimSetup::makeSysArrays( void ){
     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[l].molMembershipArray[globalAtomIndex] = molIndex;
-	  globalAtomIndex++;
-	}
-	molIndex++;
+  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[l].molMembershipArray[globalAtomIndex] = molIndex;
+    globalAtomIndex++;
+  }
+  molIndex++;
       }
     }
     
@@ -1319,7 +1319,7 @@ void SimSetup::makeSysArrays( void ){
       Exclude::createArray(info[l].n_SRI);
       the_excludes = new Exclude*[info[l].n_SRI];
       for( int ex=0; ex<info[l].n_SRI; ex++){
-	the_excludes[ex] = new Exclude(ex);
+  the_excludes[ex] = new Exclude(ex);
       }
       info[l].globalExcludes = new int;
       info[l].n_exclude = info[l].n_SRI;
@@ -1362,204 +1362,204 @@ void SimSetup::makeIntegrator( void ){
     switch( ensembleCase ){
       
     case NVE_ENS:
-	if (globals->haveZconstraints()){
-         setupZConstraint(info[k]);
-	   new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
-	}
+      if (globals->haveZconstraints()){
+        setupZConstraint(info[k]);
+        new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
+     }
 
-	else
+     else
         new NVE<RealIntegrator>( &(info[k]), the_ff );
       break;
       
     case NVT_ENS:
-	if (globals->haveZconstraints()){
-         setupZConstraint(info[k]);
-	   myNVT = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
-	}
-	else
+      if (globals->haveZconstraints()){
+        setupZConstraint(info[k]);
+        myNVT = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
+      }
+      else
         myNVT = new NVT<RealIntegrator>( &(info[k]), the_ff );
 
-      myNVT->setTargetTemp(globals->getTargetTemp());
+        myNVT->setTargetTemp(globals->getTargetTemp());
       
-      if (globals->haveTauThermostat()) 
-	myNVT->setTauThermostat(globals->getTauThermostat());
+        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;
+        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:
-	if (globals->haveZconstraints()){
-         setupZConstraint(info[k]);
-	   myNPTi = new ZConstraint<NPTi<RealIntegrator> >( &(info[k]), the_ff );
-	}
-	else
+      if (globals->haveZconstraints()){
+             setupZConstraint(info[k]);
+         myNPTi = new ZConstraint<NPTi<RealIntegrator> >( &(info[k]), the_ff );
+      }
+      else
         myNPTi = new NPTi<RealIntegrator>( &(info[k]), the_ff );
 
-	myNPTi->setTargetTemp( globals->getTargetTemp() );
-      
+      myNPTi->setTargetTemp( globals->getTargetTemp() );
+          
       if (globals->haveTargetPressure())
-	myNPTi->setTargetPressure(globals->getTargetPressure());
+        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();
+         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() );
+        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();
+         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() );
+        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;
+        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:
-	if (globals->haveZconstraints()){
-         setupZConstraint(info[k]);
-	   myNPTf = new ZConstraint<NPTf<RealIntegrator> >( &(info[k]), the_ff );
-	}
-	else
+      if (globals->haveZconstraints()){
+        setupZConstraint(info[k]);
+        myNPTf = new ZConstraint<NPTf<RealIntegrator> >( &(info[k]), the_ff );
+      }
+      else
         myNPTf = new NPTf<RealIntegrator>( &(info[k]), the_ff );
 
       myNPTf->setTargetTemp( globals->getTargetTemp());
-      
+          
       if (globals->haveTargetPressure())
-	myNPTf->setTargetPressure(globals->getTargetPressure());
+        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();
+        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() );
+        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();
+        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() );
+        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();
+        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:
-	if (globals->haveZconstraints()){
-         setupZConstraint(info[k]);
-	   myNPTim = new ZConstraint<NPTim<RealIntegrator> >( &(info[k]), the_ff );
-	}
-	else
+      if (globals->haveZconstraints()){
+        setupZConstraint(info[k]);
+        myNPTim = new ZConstraint<NPTim<RealIntegrator> >( &(info[k]), the_ff );
+      }
+      else
         myNPTim = new NPTim<RealIntegrator>( &(info[k]), the_ff );
 
-	myNPTim->setTargetTemp( globals->getTargetTemp());
-      
+        myNPTim->setTargetTemp( globals->getTargetTemp());
+          
       if (globals->haveTargetPressure())
-	myNPTim->setTargetPressure(globals->getTargetPressure());
+        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();
+        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() );
+        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();
+        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() );
+        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();
+        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:
-	if (globals->haveZconstraints()){
-         setupZConstraint(info[k]);
-	   myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >( &(info[k]), the_ff );
-	}
-	else
+      if (globals->haveZconstraints()){
+        setupZConstraint(info[k]);
+        myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >( &(info[k]), the_ff );
+      }
+      else
         myNPTfm = new NPTfm<RealIntegrator>( &(info[k]), the_ff );
 
-	myNPTfm->setTargetTemp( globals->getTargetTemp());
-      
+      myNPTfm->setTargetTemp( globals->getTargetTemp());
+
       if (globals->haveTargetPressure())
-	myNPTfm->setTargetPressure(globals->getTargetPressure());
+        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();
+        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() );
+        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();
+        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() );
+        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();
+        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");
+                 "SimSetup Error. Unrecognized ensemble in case statement.\n");
       painCave.isFatal = 1;
       simError();
     }
@@ -1578,8 +1578,8 @@ void SimSetup::initFortran( void ){
   }
   else{
     sprintf( painCave.errMsg,
-	     "SimSetup Error: unknown mixing rule -> \"%s\"\n",
-	     info[0].mixingRule );
+       "SimSetup Error: unknown mixing rule -> \"%s\"\n",
+       info[0].mixingRule );
     painCave.isFatal = 1;
     simError();
   }
@@ -1587,7 +1587,7 @@ void SimSetup::initFortran( void ){
 
 #ifdef IS_MPI
   strcpy( checkPointMsg, 
-	  "Successfully intialized the mixingRule for Fortran." );
+    "Successfully intialized the mixingRule for Fortran." );
   MPIcheckPoint();
 #endif // is_mpi
 
@@ -1595,99 +1595,99 @@ void SimSetup::setupZConstraint(SimInfo& theInfo)
 
 void SimSetup::setupZConstraint(SimInfo& theInfo)
 {
-    int nZConstraints;
-    ZconStamp** zconStamp;
-	
-    if(globals->haveZconstraintTime()){   
-      
-      //add sample time of z-constraint  into SimInfo's property list                    
-      DoubleData* zconsTimeProp = new DoubleData();
-      zconsTimeProp->setID(ZCONSTIME_ID);
-      zconsTimeProp->setData(globals->getZconsTime());
-      theInfo.addProperty(zconsTimeProp);
-    }
-    else{
-      sprintf( painCave.errMsg,
-	       "ZConstraint error: If you use an ZConstraint\n"
-	       " , you must set sample time.\n");
-      painCave.isFatal = 1;
-      simError();      
-    }
+  int nZConstraints;
+  ZconStamp** zconStamp;
 
-    //push zconsTol into siminfo, if user does not specify
-    //value for zconsTol, a default value will be used
-    DoubleData* zconsTol = new DoubleData();
-    zconsTol->setID(ZCONSTOL_ID);
-    if(globals->haveZconsTol()){
-      zconsTol->setData(globals->getZconsTol());
-    }
-	 else{
-	 	double defaultZConsTol = 0.01;
-      sprintf( painCave.errMsg,
-	       "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
-	       " , default value %f is used.\n", defaultZConsTol);
-      painCave.isFatal = 0;
-      simError();      
+  if(globals->haveZconstraintTime()){   
+    
+    //add sample time of z-constraint  into SimInfo's property list                    
+    DoubleData* zconsTimeProp = new DoubleData();
+    zconsTimeProp->setID(ZCONSTIME_ID);
+    zconsTimeProp->setData(globals->getZconsTime());
+    theInfo.addProperty(zconsTimeProp);
+  }
+  else{
+    sprintf( painCave.errMsg,
+       "ZConstraint error: If you use an ZConstraint\n"
+       " , you must set sample time.\n");
+    painCave.isFatal = 1;
+    simError();      
+  }
 
-      zconsTol->setData(defaultZConsTol);
-	 }
-    theInfo.addProperty(zconsTol);
+  //push zconsTol into siminfo, if user does not specify
+  //value for zconsTol, a default value will be used
+  DoubleData* zconsTol = new DoubleData();
+  zconsTol->setID(ZCONSTOL_ID);
+  if(globals->haveZconsTol()){
+    zconsTol->setData(globals->getZconsTol());
+  }
+  else{
+  double defaultZConsTol = 0.01;
+    sprintf( painCave.errMsg,
+       "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
+       " , default value %f is used.\n", defaultZConsTol);
+    painCave.isFatal = 0;
+    simError();      
 
-    //set Force Substraction Policy
-    StringData* zconsForcePolicy =  new StringData();
-    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
-		
-	 if(globals->haveZconsForcePolicy()){
-      zconsForcePolicy->setData(globals->getZconsForcePolicy());
-	 }	
-	 else{
-       sprintf( painCave.errMsg,
-               "ZConstraint Warning: User does not set force substraction policy, "
-               "average force substraction policy is used\n");
-       painCave.isFatal = 0;
-       simError(); 
-		 zconsForcePolicy->setData("BYNUMBER");
-	 }
-	 
-	 theInfo.addProperty(zconsForcePolicy);
-	 
-    //Determine the name of ouput file and add it into SimInfo's property list 
-    //Be careful, do not use inFileName, since it is a pointer which
-    //point to a string at master node, and slave nodes do not contain that string
-    
-    string zconsOutput(theInfo.finalName);
-    
-    zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
-    
-    StringData* zconsFilename = new StringData();
-    zconsFilename->setID(ZCONSFILENAME_ID);
-    zconsFilename->setData(zconsOutput);
-    
-    theInfo.addProperty(zconsFilename);
-		
-    //setup index, pos and other parameters of z-constraint molecules
-    nZConstraints = globals->getNzConstraints();
-    theInfo.nZconstraints = nZConstraints;
-	
-    zconStamp = globals->getZconStamp();
-    ZConsParaItem tempParaItem;
+    zconsTol->setData(defaultZConsTol);
+  }
+  theInfo.addProperty(zconsTol);
 
-    ZConsParaData* zconsParaData = new ZConsParaData();
-    zconsParaData->setID(ZCONSPARADATA_ID);
-   
-    for(int i = 0; i < nZConstraints; i++){
+  //set Force Substraction Policy
+  StringData* zconsForcePolicy =  new StringData();
+  zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
+  
+  if(globals->haveZconsForcePolicy()){
+    zconsForcePolicy->setData(globals->getZconsForcePolicy());
+  }  
+  else{
+     sprintf( painCave.errMsg,
+             "ZConstraint Warning: User does not set force substraction policy, "
+             "average force substraction policy is used\n");
+     painCave.isFatal = 0;
+     simError(); 
+     zconsForcePolicy->setData("BYNUMBER");
+  }
+ 
+ theInfo.addProperty(zconsForcePolicy);
+ 
+  //Determine the name of ouput file and add it into SimInfo's property list 
+  //Be careful, do not use inFileName, since it is a pointer which
+  //point to a string at master node, and slave nodes do not contain that string
+  
+  string zconsOutput(theInfo.finalName);
+  
+  zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
+  
+  StringData* zconsFilename = new StringData();
+  zconsFilename->setID(ZCONSFILENAME_ID);
+  zconsFilename->setData(zconsOutput);
+  
+  theInfo.addProperty(zconsFilename);
+  
+  //setup index, pos and other parameters of z-constraint molecules
+  nZConstraints = globals->getNzConstraints();
+  theInfo.nZconstraints = nZConstraints;
+
+  zconStamp = globals->getZconStamp();
+  ZConsParaItem tempParaItem;
+
+  ZConsParaData* zconsParaData = new ZConsParaData();
+  zconsParaData->setID(ZCONSPARADATA_ID);
+ 
+  for(int i = 0; i < nZConstraints; i++){
     tempParaItem.havingZPos = zconStamp[i]->haveZpos();
     tempParaItem.zPos = zconStamp[i]->getZpos();
     tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
     tempParaItem.kRatio = zconStamp[i]->getKratio();
 
     zconsParaData->addItem(tempParaItem);
-    }
+  }
 
-    //sort the parameters by index of molecules
-    zconsParaData->sortByIndex();
-	
-    //push data into siminfo, therefore, we can retrieve later
-    theInfo.addProperty(zconsParaData);
+  //sort the parameters by index of molecules
+  zconsParaData->sortByIndex();
+  
+  //push data into siminfo, therefore, we can retrieve later
+  theInfo.addProperty(zconsParaData);
       
 }