--- trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp	2005/10/11 21:57:22	653
+++ trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp	2006/10/14 20:21:26	1070
@@ -47,11 +47,11 @@
 #include <string>
 #include <map>
 #include <fstream>
+#include <algorithm>
 
 #include "config.h"
-
+#include "shapedLatticeSpherical.hpp"
 #include "nanoparticleBuilderCmd.h"
-#include "sphericalNanoparticle.hpp"
 #include "lattice/LatticeFactory.hpp"
 #include "utils/MoLocator.hpp"
 #include "lattice/Lattice.hpp"
@@ -67,7 +67,7 @@ void createMdFile(const std::string&oldMdFileName, 
 using namespace oopse;
 void createMdFile(const std::string&oldMdFileName, 
                   const std::string&newMdFileName,
-                  int numMol);
+                  std::vector<int> numMol);
 
 int main(int argc, char *argv []) {
   
@@ -78,156 +78,174 @@ int main(int argc, char *argv []) {
   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;
+  MoLocator* locator;
+  int* numMol;
+  int nComponents;
   double latticeConstant;
   std::vector<double> lc;
-  double mass;
+  double mass;										
   const double rhoConvertConst = 1.661;
   double density;
-  
-  
+  double particleRadius;
 
   Mat3x3d hmat;
-  MoLocator *locator;
-  sphericalNanoparticle *nanoparticle;
   std::vector<Vector3d> latticePos;
   std::vector<Vector3d> latticeOrt;
   int numMolPerCell;
   int nShells; /* Number of shells in nanoparticle*/
-  int numSites;
-  
+ 
   DumpWriter *writer;
   
-  // parse command line arguments
+  // Parse Command Line Arguments
   if (cmdline_parser(argc, argv, &args_info) != 0)
     exit(1);
-  
-	
-	
+         
   /* get lattice type */
-  latticeType = UpperCase(args_info.latticetype_arg);
-    
+  latticeType = "FCC";
+
   /* 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;
+      sprintf(painCave.errMsg, "No input .md file name was specified"
+	      "on the command line");
+    painCave.isFatal = 1;
     cmdline_parser_print_help();
-    exit(1);
+    simError();
   }
   
   /* parse md file and set up the system */
   SimCreator oldCreator;
   SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
   
-  nShells = 0;
-  if (args_info.coreShellRadius_given){
-    nShells = args_info.coreShellRadius_given;
-  }
+  latticeConstant = args_info.latticeCnst_arg;
+  particleRadius = args_info.radius_arg;
+  Globals* simParams = oldInfo->getSimParams();
   
-  nComponents = oldInfo->getNMoleculeStamp();
+  /* Create nanoparticle */
+  shapedLatticeSpherical nanoParticle(latticeConstant, latticeType,
+                                      particleRadius);
   
-  /* Check to see if we have enough components to build that many shells. */
-  if (nShells){
-    if (oldInfo->getNMoleculeStamp() != nShells) {
-      std::cerr << "Not enough components present in MD file to build specified number of shells"
-      << std::endl;
-      exit(1);
+  /* Build a lattice and get lattice points for this lattice constant */
+  vector<Vector3d> sites = nanoParticle.getSites();
+  vector<Vector3d> orientations = nanoParticle.getOrientations();
+
+  std::cout <<"nSites: " << sites.size() << std::endl;
+
+  /* Get number of lattice sites */
+  int nSites = sites.size();
+
+  std::vector<Component*> components = simParams->getComponents();
+  std::vector<RealType> molFractions;
+  std::vector<RealType> molecularMasses;
+  std::vector<int> nMol;
+  nComponents = components.size();
+
+ 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, "nanoparticleBuilder can't figure out molFractions "
+             "for all of the components in the <MetaData> block.");
+     painCave.isFatal = 1;
+     simError();
+   }
+ }
+ 
+ RealType totalFraction = 0.0;
+ 
+ /* Do some simple sanity checking*/
+
+  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);
   }
-  
-  
-  //creat lattice
-  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
-  
-  if (simpleLat == NULL) {
-    std::cerr << "Error in creating lattice" << std::endl;
-    exit(1);
+  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();
-  
-  /*calculate lattice constant (in Angstrom)
-  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
-                        1.0 / 3.0);*/
-  
-  latticeConstant = args_info.latticeCnst_arg;
-  particleRadius = args_info.radius_arg;
-  particleDiameter = 2.0 * particleRadius;
-  
-  /* set lattice constant */
-  lc.push_back(latticeConstant);
-  simpleLat->setLatticeConstant(lc);
-  
-  
-  /*determine the output file names*/
-  if (args_info.output_given)
-    outInitFileName = args_info.output_arg;
-  else
-    outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
-  
-  
-         
-  
-  
-  
-  /* create Molocators */
-  locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
-  
-  /* create a new spherical nanoparticle */
-  nanoparticle = new sphericalNanoparticle(particleRadius,latticeConstant);
-  /* Build a nanoparticle to see how many sites are there */
-  numSites = new int[nComponents]
-  nanoparticle.getNMol(numSites);
-  
-  numMol = new int[nComponents];
+
+  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);
+
+
+
+// recompute actual mol fractions and perform final sanity check:
+
+  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();
+  }
+     
+  vector<int> ids;
+  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
   /* Random particle is the default case*/
   if (!args_info.ShellRadius_given){
-    std::cout << "Creating a random nanoparticle" << std::endl;
-    /* Check to see if we have enough components */
-    if (nComponents != args_info.molFraction_given + 1){
-      std::cerr << "Number of components does not equal molFraction occurances." << std::endl;
-      exit 1;
-    }
-    int totComponents = 0;
-    for (int i = 0;i<nComponents-2;i++){ /* Figure out Percent for each component */
-      numMol[i] = int((double)numSites * args_info.molFraction_arg[i]);
-      totComponents += numMol[i];
-    }
-    numMol[nComponents-1] = numSites - totComponents;
-
+    /* do the iPod thing, Shuffle da vector */
+    std::random_shuffle(ids.begin(), ids.end());
   } else{ /*Handle core-shell with multiple components.*/
     std::cout << "Creating a core-shell nanoparticle." << std::endl;
     if (nComponents != args_info.ShellRadius_given + 1){
-      std::cerr << "Number of components does not equal ShellRadius occurances." << std::endl;
-      exit 1;
-    }
+      sprintf(painCave.errMsg, "Number of .md components "
+	      "does not match the number of shell radius specifications");
+      painCave.isFatal = 1;
+      simError();
+    }   
     
-    
-    
   }
 
-   //get the orientation of the cell sites
-  //for the same type of molecule in same lattice, it will not change
-   latticeOrt = simpleLat->getLatticePointsOrt();
   
   
   
-  // needed for writing out new md file.
+  outputFileName = args_info.output_arg;
+ 
   
-    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
-    outMdFileName = outPrefix + ".md";
+   //creat new .md file on fly which corrects the number of molecule     
+  createMdFile(inputFileName, outputFileName, nMol);
   
-    //creat new .md file on fly which corrects the number of molecule     
-    createMdFile(inputFileName, outMdFileName, numcomponents,numMol);
-  
   if (oldInfo != NULL)
     delete oldInfo;
   
@@ -235,89 +253,70 @@ int main(int argc, char *argv []) {
   // We need to read in new siminfo object.	
   //parse md file and set up the system
   //SimCreator NewCreator;
+  SimCreator newCreator;
+  SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
   
-  SimInfo* NewInfo = oldCreator.createSim(outMdFileName, false);
   
-  // This was so much fun the first time, lets do it again.
-  
+  // Place molecules
   Molecule* mol;
   SimInfo::MoleculeIterator mi;
   mol = NewInfo->beginMolecule(mi);
+  int l = 0;
 
-
-  for(int i = -nx; i < nx; i++) {
-     for(int j = -ny; j < ny; j++) {
-        for(int k = -nz; k < nz; k++) {
-           
-           //get the position of the cell sites
-           simpleLat->getLatticePointsPos(latticePos, i, j, k);
-           
-           for(int l = 0; l < numMolPerCell; l++) {
-#ifdef HAVE_CGAL              
-              if (myGeometry->isInsidePolyhedron(latticePos[l][0],latticePos[l][1],latticePos[l][2])){
-#endif                              
-                 if (mol != NULL) {
-                    locator->placeMol(latticePos[l], latticeOrt[l], mol);
-                 } else {
-                    std::cerr<<"Error in placing molecule " << std::endl;                    
-                 }
-                 mol = NewInfo->nextMolecule(mi);
-#ifdef HAVE_CGAL                
-              }
-#endif              
-           }
-        }
-     }
+  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++;
+    }
   }
   
 
   
   //fill Hmat
-  hmat(0, 0)= nx * latticeConstant;
+  hmat(0, 0)=  2.0*particleRadius;
   hmat(0, 1) = 0.0;
   hmat(0, 2) = 0.0;
   
   hmat(1, 0) = 0.0;
-  hmat(1, 1) = ny * latticeConstant;
+  hmat(1, 1) =  2.0*particleRadius;
   hmat(1, 2) = 0.0;
   
   hmat(2, 0) = 0.0;
   hmat(2, 1) = 0.0;
-  hmat(2, 2) = nz * latticeConstant;
+  hmat(2, 2) =  2.0*particleRadius;
   
   //set Hmat
   NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
   
   
   //create dumpwriter and write out the coordinates
-  NewInfo->setFinalConfigFileName(outInitFileName);
-  writer = new DumpWriter(NewInfo);
+  writer = new DumpWriter(NewInfo, outputFileName);
   
   if (writer == NULL) {
-    std::cerr << "error in creating DumpWriter" << std::endl;
-    exit(1);
+    sprintf(painCave.errMsg, "Error in creating dumpwrite object ");
+      painCave.isFatal = 1;
+      simError();
   }
   
   writer->writeDump();
-  std::cout << "new initial configuration file: " << outInitFileName
-            << " is generated." << std::endl;
-  
-  //delete objects
-  
-  //delete oldInfo and oldSimSetup
-  
-  if (NewInfo != NULL)
-    delete NewInfo;
-  
-  if (writer != NULL)
-    delete writer;
-  delete simpleLat;	
-  cmdline_parser_free(&args_info);
+
+  // deleting the writer will put the closing at the end of the dump file
+
+  delete writer;
+
+  // cleanup a by calling sim error.....
+  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) {
+                  std::vector<int> nMol) {
   ifstream oldMdFile;
   ofstream newMdFile;
   const int MAXLEN = 65535;
@@ -328,13 +327,17 @@ void createMdFile(const std::string&oldMdFileName, con
   newMdFile.open(newMdFileName.c_str());
   
   oldMdFile.getline(buffer, MAXLEN);
-  
+ 
+  int i = 0;
   while (!oldMdFile.eof()) {
     
     //correct molecule number
     if (strstr(buffer, "nMol") != NULL) {
-      sprintf(buffer, "\tnMol = %i;", numMol);				
-      newMdFile << buffer << std::endl;
+      if(i<nMol.size()){
+	sprintf(buffer, "\tnMol = %i;", nMol.at(i));				
+	newMdFile << buffer << std::endl;
+	i++;
+      }
     } else
       newMdFile << buffer << std::endl;