--- trunk/OOPSE-1.0/utils/sysbuilder/simpleBuilder.cpp	2004/07/28 18:42:59	1427
+++ trunk/OOPSE-1.0/utils/sysbuilder/simpleBuilder.cpp	2004/07/29 19:12:00	1436
@@ -7,17 +7,20 @@
 #include <map>
 #include <fstream>
 
-#include "simError.h"
 #include "Globals.hpp"
 #include "SimInfo.hpp"
 #include "SimSetup.hpp"
-#include "simpleBuilderCmd.hpp"
+#include "simpleBuilderCmd.h"
 #include "StringUtils.hpp"
 #include "LatticeFactory.hpp"
 #include "Vector3d.hpp"
+#include "MoLocator.hpp"
+#include "Lattice.hpp"
 
 using namespace std;
 
+void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol);
+double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF);
 
 int main( int argc, char* argv[]){
 
@@ -28,14 +31,13 @@ int main( int argc, char* argv[]){
   string outMdFileName;
   string outInitFileName;
   SimInfo* oldInfo;
-  SimInfo* newInfo;
   SimSetup* oldSimSetup;
-  SimSetup* newSimSetup;
   BaseLattice* simpleLat;
   int numMol;
   double latticeConstant;
-  vector lc;
+  vector<double> lc;
   double mass;
+  const double rhoConvertConst = 1.661;
   double density;
   int nx, ny, nz;
   double Hmat[3][3];
@@ -45,63 +47,50 @@ int main( int argc, char* argv[]){
   int numMolPerCell;
   int curMolIndex;
   DumpWriter* writer;
-
+  
   // parse command line arguments
-
   if (cmdline_parser (argc, argv, &args_info) != 0)
     exit(1) ;
+  
+  density = args_info.density_arg;
 
-  if(!args_info.density_given && !args_info.ndensity_given){
-    cerr << "density or number density must be given" << endl;
-    return false;
-  }
-
+  //get lattice type
   latticeType = UpperCase(args_info.latticetype_arg);
   if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){
     cerr << latticeType << " is an invalid lattice type" << endl;
-    cerr << LatticeFactory::getInstance()->oString << endl;
+    cerr << LatticeFactory::getInstance()->toString() << endl;
     exit(1);
   }
-	
-  nx = args_info.nx_args;
+
+  //get the number of unit cell
+  nx = args_info.nx_arg;
   if(nx <= 0){
     cerr << "The number of unit cell in h direction must be greater than 0" << endl;
     exit(1);
   }
 
-  nx = args_info.nx_args;
-  if(nx <= 0){
+  ny = args_info.ny_arg;
+  if(ny <= 0){
     cerr << "The number of unit cell in l direction must be greater than 0" << endl;
     exit(1);
   }
 
-  nx = args_info.nx_args;
-  if(nx <= 0){
+  nz = args_info.nz_arg;
+  if(nz <= 0){
     cerr << "The number of unit cell in k direction must be greater than 0" << endl;
     exit(1);
   }
 	
   //get input file name
-  if (args_info.inputs_num) {
+  if (args_info.inputs_num) 
     inputFileName = args_info.inputs[0];
-    cerr << in_name << "\n";      }
   else {		 
     cerr <<"You must specify a input file name.\n" << endl;
     cmdline_parser_print_help();
     exit(1);
   }
 
-  //determine the output file names
   
-  if (args_info.output_given){
-    outPrefix = args_info.output_arg;
-  }
-  else
-    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
-  
-  outMdFileName = outPrefix + ".md";
-  outInitFileName = outPrefix + ".in";
-
   //parse md file and set up the system
   oldInfo = new SimInfo;
   if(oldInfo == NULL){
@@ -109,22 +98,26 @@ int main( int argc, char* argv[]){
      exit(1);
   }
 
-  oldSimSetup = new SimSetup(oldInfo);  
+  oldSimSetup = new SimSetup();  
   if(oldSimSetup == NULL){
      cerr << "error in creating SimSetup" << endl;
      exit(1);
   }
-    
+
+  oldSimSetup->suspendInit();
   oldSimSetup->setSimInfo(oldInfo );
   oldSimSetup->parseFile(&inputFileName[0] );
-  oldSimSetup->createSim();
-
-
-  if(oldInfo->nComponets >=2){
+  oldSimSetup->createSim(); 
+  
+  if(oldInfo->nComponents >=2){
       cerr << "can not build the system with more than two components" << endl;
       exit(1);
   }
-
+  
+  //get mass of molecule. 
+  //Due to the design of OOPSE, given atom type, we have to query forcefiled to get the mass
+  mass = getMolMass(oldInfo->compStamps[0], oldSimSetup->getForceField());
+  
   //creat lattice
 	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
 	if(simpleLat == NULL){
@@ -132,62 +125,61 @@ int main( int argc, char* argv[]){
  		exit(1);
 	}
 
-
-  numMolPerCell = simpleLat->getNpoints();
-  //calculate lattice constant 
-  latticeConstant = pow(numMolPerCell * mass /density, 1.0/3.0);
+  numMolPerCell = simpleLat->getNumSitesPerCell();
   
+  //calculate lattice constant (in Angstrom)
+  latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0);
+  
   //set lattice constant
   lc.push_back(latticeConstant);
   simpleLat->setLatticeConstant(lc);
   
-  //calculate the total
-  numMol = args_info.nx_arg * args_info.ny_arg * args_info.nz_arg * numMolPerCell;
+  //calculate the total number of molecules
+  numMol = nx * ny * nz * numMolPerCell;
+
+  if (oldInfo->n_mol != numMol){
+
+    outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
+    outMdFileName = outPrefix + ".md";
+
+    //creat new .md file on fly which corrects the number of molecule     
+    createMdFile(inputFileName, outMdFileName, numMol);
+    cerr << "SimpleBuilder Error: the number of molecule and the density are not matched" <<endl;
+    cerr << "A new .md file: " << outMdFileName << " is generated, use it to rerun the simpleBuilder" << endl;
+    exit(1);
+  }
+  
+  //determine the output file names  
+  if (args_info.output_given)
+    outInitFileName = args_info.output_arg;
+  else
+    outInitFileName = getPrefix(inputFileName.c_str())  + ".in";
+  
+  
+  //allocat memory for storing pos, vel and etc
+  oldInfo->getConfiguration()->createArrays(oldInfo->n_atoms);
+  for (int i = 0; i < oldInfo->n_atoms; i++)
+    oldInfo->atoms[i]->setCoords();  
 
+  //creat Molocator
+  locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->getForceField());
+
   //fill Hmat
   Hmat[0][0] = nx * latticeConstant;
   Hmat[0][1] = 0.0;
   Hmat[0][2] = 0.0;
 
-  Hmat[1][0] = ny * latticeConstant;
-  Hmat[1][1] = 0.0;
+  Hmat[1][0] = 0.0;
+  Hmat[1][1] = ny * latticeConstant;
   Hmat[1][2] = 0.0;
 
-  Hmat[2][0] = nz * latticeConstant;
+  Hmat[2][0] = 0.0;
   Hmat[2][1] = 0.0;
-  Hmat[2][2] = 0.0;
-  
-  //creat new .md file on fly 
-  createMdFile(inputFileName, outMdFileName, numMol);
+  Hmat[2][2] = nz * latticeConstant ;
 
-  //parse new .md file and set up the system
-  newInfo = new SimInfo;
-  if(newInfo == NULL){
-     cerr << "error in creating SimInfo" << endl;
-     exit(1);
-  }
-
-  newSimSetup = new SimSetup(newInfo);  
-  if(newSimSetup == NULL){
-     cerr << "error in creating SimSetup" << endl;
-     exit(1);
-  }
-    
-  newSimSetup->setSimInfo(newInfo );
-  newSimSetup->parseFile(&outMdFileName[0] );
-  newSimSetup->createSim(); 
-
   //set Hmat
-  newInfo->setBoxM(Hmat);
-
-  //allocat memory for storing pos, vel and etc
-  newInfo.getConfiguration()->createArrays(newInfo.n_atoms);
-  for (int i = 0; i < newInfo.n_atoms; i++)
-    newInfo.atoms[i]->setCoords();  
-
-  //creat Molocator
-  locator = new MoLocator(newInfo->compStamps[0], newSimSetup->the_ff);
-
+  oldInfo->setBoxM(Hmat);
+  
   //place the molecules
 
   curMolIndex = 0;
@@ -204,33 +196,29 @@ int main( int argc, char* argv[]){
           simpleLat->getLatticePointsPos(latticePos, i, j, k);
 
           for(int l = 0; l < numMolPerCell; l++)
-            locator->placeMol(latticePos[l], latticeOrt[l], newInfo->molecules[curMolIndex++]);
+            locator->placeMol(latticePos[l], latticeOrt[l], &(oldInfo->molecules[curMolIndex++]));
        }
     }
   }
 
   //create dumpwriter and write out the coordinates
-  writer = new DumpWriter( newInfo );
+  oldInfo->finalName = outInitFileName;
+  writer = new DumpWriter( oldInfo );
   if(writer == NULL){
     cerr << "error in creating DumpWriter" << endl;
     exit(1);    
   }
-  writer->writeFinal();
-
-  
+  writer->writeFinal(0);
+  cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl;
   //delete objects
-  if(!oldInfo)
+
+  //delete oldInfo and oldSimSetup
+  if(oldInfo != NULL)
      delete oldInfo;
   
-  if(!oldSimSetup)
-     delete oldSimSetup;
-
-  if(!newInfo)
-     delete newInfo;
+  if(oldSimSetup != NULL)
+     delete oldSimSetup; 
   
-  if(!newSimSetup)
-     delete newSimSetup;
-
   if (writer != NULL)
     delete writer;
   return 0;
@@ -242,16 +230,18 @@ void createMdFile(const string& oldMdFileName, const s
   const int MAXLEN = 65535;
   char buffer[MAXLEN];
 
-
   //create new .md file based on old .md file
-  oldMdFile.open(oldMdFileName.c_str())
-  newMdFile.open(newMdFileName);
+  oldMdFile.open(oldMdFileName.c_str());
+  newMdFile.open(newMdFileName.c_str());
 
   oldMdFile.getline(buffer, MAXLEN);
   while(!oldMdFile.eof()){
 
-    if(line.find("nMol") < line.size())
-      sprintf(buffer, "nMol = %s;", numMol);
+    //correct molecule number
+    if(strstr(buffer, "nMol") !=NULL){      
+      sprintf(buffer, "\t\tnMol = %d;", numMol);
+      newMdFile << buffer << endl;
+    }
     else
       newMdFile << buffer << endl;
 
@@ -262,3 +252,19 @@ void createMdFile(const string& oldMdFileName, const s
   newMdFile.close();
 
 }
+
+double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF){
+  int nAtoms;
+  AtomStamp* currAtomStamp;
+  double totMass;
+  
+  totMass = 0;
+  nAtoms = molStamp->getNAtoms();
+
+  for(size_t i=0; i<nAtoms; i++){
+    currAtomStamp = molStamp->getAtom(i);
+    totMass += myFF->getAtomTypeMass(currAtomStamp->getType());
+  }
+
+  return totMass;
+}