--- trunk/OOPSE-1.0/utils/sysbuilder/simpleBuilder.cpp	2004/07/28 20:29:49	1429
+++ trunk/OOPSE-1.0/utils/sysbuilder/simpleBuilder.cpp	2004/07/29 03:31:50	1432
@@ -20,6 +20,7 @@ void createMdFile(const string& oldMdFileName, const s
 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[]){
 
@@ -38,6 +39,7 @@ int main( int argc, char* argv[]){
   double latticeConstant;
   vector<double> lc;
   double mass;
+  const double rhoConvertConst = 1.661;
   double density;
   int nx, ny, nz;
   double Hmat[3][3];
@@ -47,24 +49,37 @@ 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) ;
 
+  //process density
   if(!args_info.density_given && !args_info.ndensity_given){
-    cerr << "density or number density must be given" << endl;
-    return false;
+    cerr << "SimpleBuilder Error: density or number density must be given" << endl;
+    cmdline_parser_print_help();
+    exit(1);
   }
+  else if(args_info.density_given){
+    if( args_info.ndensity_given)
+      cerr << "SimpleBuilder Warning: both of density and number density are given, use density" << endl;
+    
+    density = args_info.density_arg;
+  }
+  else if(args_info.ndensity_given){
+    //convert number density to density 
+  }
 
+
+  //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()->toString() << endl;
     exit(1);
   }
-	
+
+  //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;
@@ -118,13 +133,17 @@ int main( int argc, char* argv[]){
     
   oldSimSetup->setSimInfo(oldInfo );
   oldSimSetup->parseFile(&inputFileName[0] );
-
-
+  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,30 +151,17 @@ int main( int argc, char* argv[]){
  		exit(1);
 	}
 
-
   numMolPerCell = simpleLat->getNumSitesPerCell();
-  //calculate lattice constant 
-  latticeConstant = pow(numMolPerCell * mass /density, 1.0/3.0);
   
+  //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
+  //calculate the total number of molecules
   numMol = args_info.nx_arg * args_info.ny_arg * args_info.nz_arg * numMolPerCell;
-
-  //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][2] = 0.0;
-
-  Hmat[2][0] = nz * latticeConstant;
-  Hmat[2][1] = 0.0;
-  Hmat[2][2] = 0.0;
   
   //creat new .md file on fly 
   createMdFile(inputFileName, outMdFileName, numMol);
@@ -166,7 +172,7 @@ int main( int argc, char* argv[]){
      cerr << "error in creating SimInfo" << endl;
      exit(1);
   }
-
+  
   newSimSetup = new SimSetup();  
   if(newSimSetup == NULL){
      cerr << "error in creating SimSetup" << endl;
@@ -177,9 +183,6 @@ int main( int argc, char* argv[]){
   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++)
@@ -188,6 +191,22 @@ int main( int argc, char* argv[]){
   //creat Molocator
   locator = new MoLocator(newInfo->compStamps[0], newSimSetup->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][2] = 0.0;
+
+  Hmat[2][0] = nz * latticeConstant;
+  Hmat[2][1] = 0.0;
+  Hmat[2][2] = 0.0;
+
+  //set Hmat
+  newInfo->setBoxM(Hmat);
+  
   //place the molecules
 
   curMolIndex = 0;
@@ -216,7 +235,6 @@ int main( int argc, char* argv[]){
     exit(1);    
   }
   writer->writeFinal(0);
-
   
   //delete objects
   if(!oldInfo)
@@ -265,3 +283,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;
+}