--- trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp 2006/01/19 22:54:23 875 +++ trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp 2006/10/13 20:16:59 1069 @@ -47,11 +47,11 @@ #include #include #include +#include #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 numMol); int main(int argc, char *argv []) { @@ -79,29 +79,31 @@ int main(int argc, char *argv []) { std::string latticeType; std::string inputFileName; std::string outPrefix; - std::string outMdFileName; + std::string outputFileName; std::string outInitFileName; Lattice *simpleLat; - int numMol; + MoLocator* locator; + int* numMol; + int nComponents; double latticeConstant; std::vector lc; - double mass; + double mass; + const double rhoConvertConst = 1.661; double density; + double particleRadius; Mat3x3d hmat; - MoLocator *locator; - sphericalNanoparticle *nanoparticle; std::vector latticePos; std::vector latticeOrt; int numMolPerCell; int nShells; /* Number of shells in nanoparticle*/ - int numSites; + DumpWriter *writer; @@ -112,122 +114,162 @@ int main(int argc, char *argv []) { /* 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; - } - nComponents = oldInfo->getNMoleculeStamp(); - - /* 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); - } - } - - - //creat lattice - simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); - - if (simpleLat == NULL) { - std::cerr << "Error in creating lattice" << std::endl; - exit(1); - } - - 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; + Globals* simParams = oldInfo->getSimParams(); - /* 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] - int numSites = nanoparticle.getNMol(simpleLat); + /* Create nanoparticle */ + shapedLatticeSpherical nanoParticle(latticeConstant,latticeType,particleRadius); + /* Build a lattice and get lattice points for this lattice constant */ + vector sites = nanoParticle.getPoints(); + vector orientations = nanoParticle.getPointsOrt(); - /* 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; + std::cout <<"nSites: " << sites.size() << std::endl; + + /* Get number of lattice sites */ + int nSites = sites.size(); + + + + + std::vector components = simParams->getComponents(); + std::vector molFractions; + std::vector molecularMasses; + std::vector 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 block."); + painCave.isFatal = 1; + simError(); } - int totComponents = 0; - for (int i = 0;i 1.0) { + sprintf(painCave.errMsg, "One of the requested molFractions was" + " greater than one!"); + painCave.isFatal = 1; + simError(); + } + totalFraction += molFractions.at(i); + } + 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(); + } + 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 ids; + for (int i = 0; i < sites.size(); i++) ids.push_back(i); + /* Random particle is the default case*/ + if (!args_info.ShellRadius_given){ + /* 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 +277,72 @@ 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 nMol) { ifstream oldMdFile; ofstream newMdFile; const int MAXLEN = 65535; @@ -328,13 +353,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