--- 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 19:12:00 1436 @@ -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[]){ @@ -30,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; double mass; + const double rhoConvertConst = 1.661; double density; int nx, ny, nz; double Hmat[3][3]; @@ -47,24 +47,22 @@ 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()->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; @@ -92,17 +90,7 @@ int main( int argc, char* argv[]){ 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){ @@ -115,16 +103,21 @@ int main( int argc, char* argv[]){ cerr << "error in creating SimSetup" << endl; exit(1); } - + + oldSimSetup->suspendInit(); 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,62 +125,61 @@ 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 - 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" <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(); - 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->getForceField()); - + 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(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; @@ -241,7 +229,6 @@ void createMdFile(const string& oldMdFileName, const s ofstream newMdFile; const int MAXLEN = 65535; char buffer[MAXLEN]; - string line; //create new .md file based on old .md file oldMdFile.open(oldMdFileName.c_str()); @@ -250,8 +237,8 @@ void createMdFile(const string& oldMdFileName, const s oldMdFile.getline(buffer, MAXLEN); while(!oldMdFile.eof()){ - if(line.find("nMol") < line.size()){ - + //correct molecule number + if(strstr(buffer, "nMol") !=NULL){ sprintf(buffer, "\t\tnMol = %d;", numMol); newMdFile << buffer << endl; } @@ -265,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; igetAtom(i); + totMass += myFF->getAtomTypeMass(currAtomStamp->getType()); + } + + return totMass; +}