--- trunk/src/applications/randomBuilder/randomBuilder.cpp	2006/04/25 22:59:27	950
+++ trunk/src/applications/randomBuilder/randomBuilder.cpp	2006/10/10 18:34:12	1067
@@ -42,7 +42,7 @@
  *
  *  Created by Charles F. Vardeman II on 10 Apr 2006.
  *  @author  Charles F. Vardeman II
- *  @version $Id: randomBuilder.cpp,v 1.1 2006-04-25 22:59:27 chuckv Exp $
+ *  @version $Id: randomBuilder.cpp,v 1.6 2006-10-10 18:34:12 gezelter Exp $
  *
  */
 
@@ -73,36 +73,29 @@ void createMdFile(const std::string&oldMdFileName, 
 
 void createMdFile(const std::string&oldMdFileName, 
                   const std::string&newMdFileName,
-                  int components,int* numMol);
+                  std::vector<int> nMol);
 
 int main(int argc, char *argv []) {
 
-  //register force fields
+  // register force fields
   registerForceFields();
   registerLattice();
     
   gengetopt_args_info args_info;
   std::string latticeType;
   std::string inputFileName;
-  std::string outPrefix;
-  std::string outMdFileName;
-  std::string outInitFileName;
+  std::string outputFileName;
   Lattice *simpleLat;
-  int* numMol;
-  double latticeConstant;
-  std::vector<double> lc;
-  double mass;
-  const double rhoConvertConst = 1.661;
-  double density;
-  int nx,
-    ny,
-    nz;
+  RealType latticeConstant;
+  std::vector<RealType> lc;
+  const RealType rhoConvertConst = 1.661;
+  RealType density;
+  int nx, ny, nz;
   Mat3x3d hmat;
   MoLocator *locator;
   std::vector<Vector3d> latticePos;
   std::vector<Vector3d> latticeOrt;
-  int numMolPerCell;
-  int curMolIndex;
+  int nMolPerCell;
   DumpWriter *writer;
 
   // parse command line arguments
@@ -112,7 +105,7 @@ int main(int argc, char *argv []) {
   density = args_info.density_arg;
 
   //get lattice type
-  latticeType = UpperCase(args_info.latticetype_arg);
+  latticeType = "FCC";
 
   simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
     
@@ -122,125 +115,188 @@ int main(int argc, char *argv []) {
     painCave.isFatal = 1;
     simError();
   }
+  nMolPerCell = simpleLat->getNumSitesPerCell();
 
-  //get the number of unit cell
+  //get the number of unit cells in each direction:
+
   nx = args_info.nx_arg;
 
   if (nx <= 0) {
-    std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl;
-    exit(1);
+    sprintf(painCave.errMsg, "The number of unit cells in the x direction "
+            "must be greater than 0.");
+    painCave.isFatal = 1;
+    simError();
   }
 
   ny = args_info.ny_arg;
 
   if (ny <= 0) {
-    std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl;
-    exit(1);
+    sprintf(painCave.errMsg, "The number of unit cells in the y direction "
+            "must be greater than 0.");
+    painCave.isFatal = 1;
+    simError();
   }
 
   nz = args_info.nz_arg;
 
   if (nz <= 0) {
-    std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl;
-    exit(1);
+    sprintf(painCave.errMsg, "The number of unit cells in the z direction "
+            "must be greater than 0.");
+    painCave.isFatal = 1;
+    simError();
   }
 
+  int nSites = nMolPerCell * nx * ny * nz;
+
   //get input file name
   if (args_info.inputs_num)
     inputFileName = args_info.inputs[0];
   else {
-    std::cerr << "You must specify a input file name.\n" << std::endl;
-    cmdline_parser_print_help();
-    exit(1);
+    sprintf(painCave.errMsg, "No input .md file name was specified "
+            "on the command line");
+    painCave.isFatal = 1;
+    simError();
   }
 
   //parse md file and set up the system
+
   SimCreator oldCreator;
   SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
   Globals* simParams = oldInfo->getSimParams();
 
-  int nComponents =simParams->getNComponents();
-  if (oldInfo->getNMoleculeStamp()>= 2) {
-    std::cerr << "can not build the system with more than two components"
-	      << std::endl;
-    exit(1);
-  }
+  // Calculate lattice constant (in Angstroms)
 
-  //get mass of molecule. 
+  std::vector<Component*> components = simParams->getComponents();
+  std::vector<RealType> molFractions;
+  std::vector<RealType> molecularMasses;
+  std::vector<int> nMol;
+  int nComponents = components.size();
 
-  mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
+  if (nComponents == 1) {
+    molFractions.push_back(1.0);    
+  } else {
+    if (args_info.molFraction_given == nComponents) {
+      for (int i = 0; i < nComponents; i++) {
+        molFractions.push_back(args_info.molFraction_arg[i]);
+      }
+    } else if (args_info.molFraction_given == nComponents-1) {
+      RealType remainingFraction = 1.0;
+      for (int i = 0; i < nComponents-1; i++) {
+        molFractions.push_back(args_info.molFraction_arg[i]);
+        remainingFraction -= molFractions[i];
+      }
+      molFractions.push_back(remainingFraction);
+    } else {    
+      sprintf(painCave.errMsg, "randomBuilder can't figure out molFractions "
+              "for all of the components in the <MetaData> block.");
+      painCave.isFatal = 1;
+      simError();
+    }
+  }
 
-  //creat lattice
-  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
+  // do some sanity checking:
+  
+  RealType totalFraction = 0.0;
 
-  if (simpleLat == NULL) {
-    std::cerr << "Error in creating lattice" << std::endl;
-    exit(1);
+  for (int i = 0; i < nComponents; i++) {
+    if (molFractions.at(i) < 0.0) {
+      sprintf(painCave.errMsg, "One of the requested molFractions was"
+              " less than zero!");
+      painCave.isFatal = 1;
+      simError();
+    }
+    if (molFractions.at(i) > 1.0) {
+      sprintf(painCave.errMsg, "One of the requested molFractions was"
+              " greater than one!");
+      painCave.isFatal = 1;
+      simError();
+    }
+    totalFraction += molFractions.at(i);
   }
+  if (abs(totalFraction - 1.0) > 1e-6) {
+    sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
+    painCave.isFatal = 1;
+    simError();
+  }
 
-  numMolPerCell = simpleLat->getNumSitesPerCell();
+  int remaining = nSites;
+  for (int i=0; i < nComponents-1; i++) {    
+    nMol.push_back(int((RealType)nSites * molFractions.at(i)));
+    remaining -= nMol.at(i);
+  }
+  nMol.push_back(remaining);
 
-  //calculate lattice constant (in Angstrom)
-  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
-			1.0 / 3.0);
+  // recompute actual mol fractions and perform final sanity check:
 
-  //set lattice constant
+  int totalMolecules = 0;
+  RealType totalMass = 0.0;
+  for (int i=0; i < nComponents; i++) {
+    molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
+    totalMolecules += nMol.at(i);
+    molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
+                                         oldInfo->getForceField()));
+    totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
+  }
+  RealType avgMass = totalMass / (RealType) totalMolecules;
+
+  if (totalMolecules != nSites) {
+    sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
+            "to the number of lattice sites!");
+    painCave.isFatal = 1;
+    simError();
+  }
+     
+  latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
+			(RealType)(1.0 / 3.0));
+  
+  // Set the lattice constant
+  
   lc.push_back(latticeConstant);
   simpleLat->setLatticeConstant(lc);
+  
+  // Calculate the lattice sites and fill the lattice vector.
 
-  //calculate the total number of molecules
-  int totMol = nx * ny * nz * numMolPerCell;
+  // Get the standard orientations of the cell sites
 
-  // Calculate the lattice sites and fill lattice vector.
-  vector<Vector3d> latticeSites;
+  latticeOrt = simpleLat->getLatticePointsOrt();
+
+  vector<Vector3d> sites;
+  vector<Vector3d> orientations;
   
   for(int i = 0; i < nx; i++) {
     for(int j = 0; j < ny; j++) {
       for(int k = 0; k < nz; k++) {
 
-	//get the position of the cell sites
+	// Get the position of the cell sites
+
 	simpleLat->getLatticePointsPos(latticePos, i, j, k);
 
-	for(int l = 0; l < numMolPerCell; l++) {
-	  latticeSites.push_back(latticePos[l]);
+	for(int l = 0; l < nMolPerCell; l++) {
+	  sites.push_back(latticePos[l]);
+          orientations.push_back(latticeOrt[l]);
 	}
       }
     }
   }
-
-  int numSites = latticeSites.size();
-
-  numMol = new int[nComponents];
-  if (nComponents != args_info.molFraction_given && nComponents != 1){
-    std::cerr << "Number of components does not equal molFraction occurances." << std::endl;
-    exit(1);
-  }
-  int totComponents = 0;
-  for (int i = 0;i<nComponents-1;i++){ /* Figure out Percent for each component */
-    numMol[i] = int((double)numSites * args_info.molFraction_arg[i]);
-    std::cout<<numMol[i]<<std::endl;
-    totComponents += numMol[i];
-  }
-  numMol[nComponents-1] = numSites - totComponents;
   
+  outputFileName = args_info.output_arg;
 
-  outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
-  outMdFileName = outPrefix + ".md";
+  // create a new .md file on the fly which corrects the number of molecules
 
-  //creat new .md file on fly which corrects the number of molecule     
-  createMdFile(inputFileName, outMdFileName, nComponents,numMol);
+  createMdFile(inputFileName, outputFileName, nMol);
 
+  if (oldInfo != NULL)
+    delete oldInfo;
 
+  // We need to read in the new SimInfo object, then Parse the 
+  // md file and set up the system
 
-  //determine the output file names  
-  if (args_info.output_given){
-    outInitFileName = args_info.output_arg;
-  }else{
-    outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
-  }
+  SimCreator newCreator;
+  SimInfo* newInfo = newCreator.createSim(outputFileName, false);
 
-  //fill Hmat
-  hmat(0, 0)= nx * latticeConstant;
+  // fill Hmat
+
+  hmat(0, 0) = nx * latticeConstant;
   hmat(0, 1) = 0.0;
   hmat(0, 2) = 0.0;
 
@@ -252,83 +308,63 @@ int main(int argc, char *argv []) {
   hmat(2, 1) = 0.0;
   hmat(2, 2) = nz * latticeConstant;
 
-  //set Hmat
-  oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
+  // Set Hmat
 
-  //place the molecules
+  newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
 
-  curMolIndex = 0;
+  // place the molecules
 
-  //get the orientation of the cell sites
-  //for the same type of molecule in same lattice, it will not change
-  latticeOrt = simpleLat->getLatticePointsOrt();
+  // Randomize a vector of ints:
 
- 
-  /* Randomize position vector */
-  std::random_shuffle(latticeSites.begin(), latticeSites.end());
+  vector<int> ids;
+  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
+  std::random_shuffle(ids.begin(), ids.end());
 
-
-  if (oldInfo != NULL)
-    delete oldInfo;
-  
-  
-  // We need to read in new siminfo object.	
-  //parse md file and set up the system
-  //SimCreator NewCreator;
-  SimCreator newCreator;
-  SimInfo* NewInfo = newCreator.createSim(outMdFileName, false);
-  
-  /* create Molocators */
-  locator = new MoLocator(NewInfo->getMoleculeStamp(0), NewInfo->getForceField());
-  
-  
-  
   Molecule* mol;
-  SimInfo::MoleculeIterator mi;
-  mol = NewInfo->beginMolecule(mi);
   int l = 0;
-  for (mol = NewInfo->beginMolecule(mi); mol != NULL; mol = NewInfo->nextMolecule(mi)) {
-    locator->placeMol(latticeSites[l], latticeOrt[l], mol);
-    l++;
+  for (int i = 0; i < nComponents; i++){
+    locator = new MoLocator(newInfo->getMoleculeStamp(i), 
+                            newInfo->getForceField());
+    for (int n = 0; n < nMol.at(i); n++) {
+      mol = newInfo->getMoleculeByGlobalIndex(l);
+      locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
+      l++;
+    }
   }
+  
+  // Create DumpWriter and write out the coordinates
 
+  writer = new DumpWriter(newInfo, outputFileName);
 
-
-  //create dumpwriter and write out the coordinates
-  oldInfo->setFinalConfigFileName(outInitFileName);
-  writer = new DumpWriter(oldInfo);
-
   if (writer == NULL) {
-    std::cerr << "error in creating DumpWriter" << std::endl;
-    exit(1);
+    sprintf(painCave.errMsg, "error in creating DumpWriter");
+    painCave.isFatal = 1;
+    simError();
   }
 
   writer->writeDump();
-  std::cout << "new initial configuration file: " << outInitFileName
-	    << " is generated." << std::endl;
 
-  //delete objects
+  // deleting the writer will put the closing at the end of the dump file.
 
-  //delete oldInfo and oldSimSetup
-  if (oldInfo != NULL)
-    delete oldInfo;
+  delete writer;
 
-  if (writer != NULL)
-    delete writer;
-    
-  delete simpleLat;
-
+  sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been "
+          "generated.\n", outputFileName.c_str());
+  painCave.isFatal = 0;
+  simError();
   return 0;
 }
 
-void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
-                  int components,int* numMol) {
+void createMdFile(const std::string&oldMdFileName, 
+                  const std::string&newMdFileName, 
+                  std::vector<int> nMol) {
   ifstream oldMdFile;
   ofstream newMdFile;
   const int MAXLEN = 65535;
   char buffer[MAXLEN];
   
   //create new .md file based on old .md file
+
   oldMdFile.open(oldMdFileName.c_str());
   newMdFile.open(newMdFileName.c_str());
   
@@ -339,8 +375,8 @@ void createMdFile(const std::string&oldMdFileName, con
     
     //correct molecule number
     if (strstr(buffer, "nMol") != NULL) {
-      if(i<components){
-	sprintf(buffer, "\tnMol = %i;", numMol[i]);				
+      if(i<nMol.size()){
+	sprintf(buffer, "\tnMol = %i;", nMol.at(i));
 	newMdFile << buffer << std::endl;
 	i++;
       }