--- 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 03:31:50 1432 @@ -7,17 +7,20 @@ #include #include -#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[]){ @@ -34,8 +37,9 @@ int main( int argc, char* argv[]){ BaseLattice* simpleLat; int numMol; double latticeConstant; - vector lc; + vector lc; double mass; + const double rhoConvertConst = 1.661; double density; int nx, ny, nz; double Hmat[3][3]; @@ -45,46 +49,58 @@ 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()->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(); @@ -109,7 +125,7 @@ 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); @@ -117,14 +133,17 @@ int main( int argc, char* argv[]){ 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,30 +151,17 @@ 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 + //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,8 +172,8 @@ int main( int argc, char* argv[]){ cerr << "error in creating SimInfo" << endl; exit(1); } - - newSimSetup = new SimSetup(newInfo); + + newSimSetup = new SimSetup(); if(newSimSetup == NULL){ cerr << "error in creating SimSetup" << endl; exit(1); @@ -177,17 +183,30 @@ 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++) - newInfo.atoms[i]->setCoords(); + 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); + 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; @@ -204,7 +223,7 @@ 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], &(newInfo->molecules[curMolIndex++])); } } } @@ -215,8 +234,7 @@ int main( int argc, char* argv[]){ cerr << "error in creating DumpWriter" << endl; exit(1); } - writer->writeFinal(); - + writer->writeFinal(0); //delete objects if(!oldInfo) @@ -241,17 +259,20 @@ 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()) - 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); + if(line.find("nMol") < line.size()){ + + sprintf(buffer, "\t\tnMol = %d;", numMol); + newMdFile << buffer << endl; + } else newMdFile << buffer << endl; @@ -262,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; igetAtom(i); + totMass += myFF->getAtomTypeMass(currAtomStamp->getType()); + } + + return totMass; +}