--- trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/21 17:42:12	378
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2003/03/26 20:22:02	407
@@ -333,30 +333,6 @@ void SimSetup::createSim( void ){
 
   // get some of the tricky things that may still be in the globals
 
-  if( simnfo->n_dipoles ){
-
-    if( !the_globals->haveRRF() ){
-      sprintf( painCave.errMsg,
-	       "SimSetup Error, system has dipoles, but no rRF was set.\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-    if( !the_globals->haveDielectric() ){
-      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();
@@ -426,9 +402,84 @@ void SimSetup::createSim( void ){
   }
 
 
-
+  if (the_globals->getUseRF() ) {
+    simnfo->useReactionField = 1;
+  
+    if( !the_globals->haveECR() ){
+      sprintf( painCave.errMsg,
+     	       "SimSetup Warning: using default value of 1/2 the smallest "
+               "box length for the electrostaticCutoffRadius.\n"
+               "I hope you have a very fast processor!\n");
+      painCave.isFatal = 0;
+      simError();
+      double smallest;
+      smallest = simnfo->box_x;
+      if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
+      if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
+      simnfo->ecr = 0.5 * smallest;
+    } else {
+      simnfo->ecr        = the_globals->getECR();
+    }
 
+    if( !the_globals->haveEST() ){
+      sprintf( painCave.errMsg,
+     	       "SimSetup Warning: using default value of 0.05 * the "
+               "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
+               );
+      painCave.isFatal = 0;
+      simError();
+      simnfo->est = 0.05 * simnfo->ecr;
+    } else {
+      simnfo->est        = the_globals->getEST();
+    }
+    
+    if(!the_globals->haveDielectric() ){
+      sprintf( painCave.errMsg,
+     	       "SimSetup Error: You are trying to use Reaction Field without"
+               "setting a dielectric constant!\n"
+               );
+      painCave.isFatal = 1;
+      simError();
+    }
+    simnfo->dielectric = the_globals->getDielectric();  
+  } else {
+    if (simnfo->n_dipoles) {
+      
+      if( !the_globals->haveECR() ){
+        sprintf( painCave.errMsg,
+                 "SimSetup Warning: using default value of 1/2 the smallest"
+                 "box length for the electrostaticCutoffRadius.\n"
+                 "I hope you have a very fast processor!\n");
+        painCave.isFatal = 0;
+        simError();
+        double smallest;
+        smallest = simnfo->box_x;
+        if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
+        if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
+        simnfo->ecr = 0.5 * smallest;
+      } else {
+        simnfo->ecr        = the_globals->getECR();
+      }
+      
+      if( !the_globals->haveEST() ){
+        sprintf( painCave.errMsg,
+                 "SimSetup Warning: using default value of 5% of the"
+                 "electrostaticCutoffRadius for the "
+                 "electrostaticSkinThickness\n"
+                 );
+        painCave.isFatal = 0;
+        simError();
+        simnfo->est = 0.05 * simnfo->ecr;
+      } else {
+        simnfo->est        = the_globals->getEST();
+      }
+    }
+  }  
 
+#ifdef IS_MPI
+  strcpy( checkPointMsg, "electrostatic parameters check out" );
+  MPIcheckPoint();
+#endif // is_mpi
 
  if( the_globals->haveInitialConfig() ){
  
@@ -620,6 +671,128 @@ void SimSetup::createSim( void ){
 #endif // is_mpi
 }
 
+
+void SimSetup::makeMolecules( void ){
+
+  int i, j, exI, exJ, tempEx, stampID, atomOffset;
+  molInit info;
+  DirectionalAtom* dAtom;
+  AtomStamp* currentAtom;
+  BondStamp* currentBond;
+  BendStamp* currentBend;
+  TorsionStamp* currentTorsion;
+  
+  //init the forceField paramters
+
+  the_ff->readParams();
+
+  
+  // init the molecules
+
+  atomOffset = 0;
+  for(i=0; i<simnfo->n_mol; i++){
+    
+    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.myAtoms = &the_atoms[atomOffset];
+    info.myBonds = new Bond*[info.nBonds];
+    info.myBends = new Bend*[info.nBends];
+    info.myTorsions = new Torsions*[info.nTorsions];
+
+    theBonds = new bond_pair[info.nBonds];
+    theBends = new bend_set[info.nBends];
+    theTorsions = new torsion_set[info.nTorsions];
+    
+    // make the Atoms
+    
+    for(j=0; j<info.nAtoms; j++){
+      
+      currentAtom = theComponents[stampID]->getAtom( j );
+      if( currentAtom->haveOrientation() ){
+	
+	dAtom = new DirectionalAtom(j + atomOffset);
+	simnfo->n_oriented++;
+	info.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{
+	info.myAtoms[j] = new GeneralAtom(j + atomOffset);
+      }
+      info.myAtoms[j]->setType( currentAtom->getType() );
+    
+#ifdef IS_MPI
+      
+      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
+      
+#endif // is_mpi
+    } 
+    
+    // make the bonds
+    for(j=0; j<nBonds; j++){
+      
+      currentBond = comp_stamps[stampID]->getBond( j );
+      theBonds[j].a = currentBond->getA() + atomOffset;
+      theBonds[j].b = currentBond->getB() + atomOffset;
+
+      exI = theBonds[i].a;
+      exJ = theBonds[i].b;
+
+      // exclude_I must always be the smaller of the pair
+      if( exI > exJ ){
+	tempEx = exI;
+	exI = exJ;
+	exJ = tempEx;
+      }
+#ifdef IS_MPI
+      
+      the_excludes[index*2] =     
+	the_atoms[exI]->getGlobalIndex() + 1;
+      the_excludes[index*2 + 1] = 
+	the_atoms[exJ]->getGlobalIndex() + 1;
+      
+#else  // isn't MPI
+      
+      the_excludes[index*2] =     exI + 1; 
+      the_excludes[index*2 + 1] = exJ + 1;
+      // fortran index from 1 (hence the +1 in the indexing)
+
+#endif  //is_mpi
+    
+  }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 void SimSetup::makeAtoms( void ){
 
   int i, j, k, index;