--- trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp 2006/10/13 20:16:59 1069 +++ trunk/src/applications/nanoparticleBuilder/nanoparticleBuilder.cpp 2009/11/25 20:02:06 1390 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,15 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ #include @@ -64,7 +64,7 @@ using namespace std; #include "utils/StringUtils.hpp" using namespace std; -using namespace oopse; +using namespace OpenMD; void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName, std::vector numMol); @@ -78,41 +78,25 @@ int main(int argc, char *argv []) { gengetopt_args_info args_info; std::string latticeType; std::string inputFileName; - std::string outPrefix; std::string outputFileName; - std::string outInitFileName; - - - Lattice *simpleLat; MoLocator* locator; - int* numMol; int nComponents; double latticeConstant; std::vector lc; - double mass; - const double rhoConvertConst = 1.661; - double density; - double particleRadius; - - + RealType particleRadius; Mat3x3d hmat; std::vector latticePos; std::vector latticeOrt; - int numMolPerCell; - int nShells; /* Number of shells in nanoparticle*/ - DumpWriter *writer; // Parse Command Line Arguments if (cmdline_parser(argc, argv, &args_info) != 0) exit(1); - - - + /* get lattice type */ latticeType = "FCC"; @@ -120,8 +104,8 @@ int main(int argc, char *argv []) { if (args_info.inputs_num) inputFileName = args_info.inputs[0]; else { - sprintf(painCave.errMsg, "No input .md file name was specified" - "on the command line"); + sprintf(painCave.errMsg, "No input .md file name was specified " + "on the command line"); painCave.isFatal = 1; cmdline_parser_print_help(); simError(); @@ -131,46 +115,111 @@ int main(int argc, char *argv []) { SimCreator oldCreator; SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); - - latticeConstant = args_info.latticeCnst_arg; + latticeConstant = args_info.latticeConstant_arg; particleRadius = args_info.radius_arg; Globals* simParams = oldInfo->getSimParams(); - - /* create Molocators */ - locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField()); - /* Create nanoparticle */ - shapedLatticeSpherical nanoParticle(latticeConstant,latticeType,particleRadius); + shapedLatticeSpherical nanoParticle(latticeConstant, latticeType, + particleRadius); /* Build a lattice and get lattice points for this lattice constant */ - vector sites = nanoParticle.getPoints(); - vector orientations = nanoParticle.getPointsOrt(); - - std::cout <<"nSites: " << sites.size() << std::endl; - - /* Get number of lattice sites */ - int nSites = sites.size(); + vector sites = nanoParticle.getSites(); + vector orientations = nanoParticle.getOrientations(); + std::vector vacancyTargets; + vector isVacancy; + + Vector3d myLoc; + RealType myR; + for (int i = 0; i < sites.size(); i++) + isVacancy.push_back(false); + if (args_info.vacancyPercent_given) { + if (args_info.vacancyPercent_arg < 0.0 || args_info.vacancyPercent_arg > 100.0) { + sprintf(painCave.errMsg, "vacancyPercent was set to a non-sensical value."); + painCave.isFatal = 1; + simError(); + } else { + RealType vF = args_info.vacancyPercent_arg / 100.0; + RealType vIR; + RealType vOR; + if (args_info.vacancyInnerRadius_given) { + vIR = args_info.vacancyInnerRadius_arg; + } else { + vIR = 0.0; + } + if (args_info.vacancyOuterRadius_given) { + vOR = args_info.vacancyOuterRadius_arg; + } else { + vOR = particleRadius; + } + if (vIR >= 0.0 && vOR <= particleRadius && vOR >= vIR) { + + for (int i = 0; i < sites.size(); i++) { + myLoc = sites[i]; + myR = myLoc.length(); + if (myR >= vIR && myR <= vOR) { + vacancyTargets.push_back(i); + } + } + std::random_shuffle(vacancyTargets.begin(), vacancyTargets.end()); + + int nTargets = vacancyTargets.size(); + vacancyTargets.resize((int)(vF * nTargets)); + + + sprintf(painCave.errMsg, "Removing %d atoms from randomly-selected\n" + "\tsites between %lf and %lf.", (int) vacancyTargets.size(), + vIR, vOR); + painCave.isFatal = 0; + simError(); + isVacancy.clear(); + for (int i = 0; i < sites.size(); i++) { + bool vac = false; + for (int j = 0; j < vacancyTargets.size(); j++) { + if (i == vacancyTargets[j]) vac = true; + } + isVacancy.push_back(vac); + } + + } else { + sprintf(painCave.errMsg, "Something is strange about the vacancy\n" + "\tinner or outer radii. Check their values."); + painCave.isFatal = 1; + simError(); + } + } + } + /* Get number of lattice sites */ + int nSites = sites.size() - vacancyTargets.size(); + std::vector components = simParams->getComponents(); std::vector molFractions; + std::vector shellRadii; std::vector molecularMasses; std::vector nMol; + std::map componentFromSite; nComponents = components.size(); - - - if (nComponents == 1) { + if (args_info.molFraction_given && args_info.shellRadius_given) { + sprintf(painCave.errMsg, "Specify either molFraction or shellRadius " + "arguments, but not both!"); + painCave.isFatal = 1; + simError(); + } + + if (nComponents == 1) { molFractions.push_back(1.0); - } else { - if (args_info.molFraction_given == nComponents) { + shellRadii.push_back(particleRadius); + } else if (args_info.molFraction_given) { + if ((int)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) { + } else if ((int)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]); @@ -183,165 +232,217 @@ int main(int argc, char *argv []) { painCave.isFatal = 1; simError(); } - } - - RealType totalFraction = 0.0; - - /* Do some simple sanity checking*/ - - - - for (int i = 0; i < nComponents; i++) { - if (molFractions.at(i) < 0.0) { - sprintf(painCave.errMsg, "One of the requested molFractions was" - " less than zero!"); + } else if ((int)args_info.shellRadius_given) { + if ((int)args_info.shellRadius_given == nComponents) { + for (int i = 0; i < nComponents; i++) { + shellRadii.push_back(args_info.shellRadius_arg[i]); + } + } else if ((int)args_info.shellRadius_given == nComponents-1) { + for (int i = 0; i < nComponents-1; i++) { + shellRadii.push_back(args_info.shellRadius_arg[i]); + } + shellRadii.push_back(particleRadius); + } else { + sprintf(painCave.errMsg, "nanoparticleBuilder can't figure out the\n" + "\tshell radii for all of the components in the block."); painCave.isFatal = 1; simError(); } - if (molFractions.at(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"); + } else { + sprintf(painCave.errMsg, "You have a multi-component block,\n" + "\tbut have not specified either molFraction or shellRadius arguments."); 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); + + if (args_info.molFraction_given) { + RealType totalFraction = 0.0; + + /* Do some simple sanity checking*/ + + for (int i = 0; i < nComponents; i++) { + if (molFractions.at(i) < 0.0) { + sprintf(painCave.errMsg, "One of the requested molFractions was" + " less than zero!"); + painCave.isFatal = 1; + simError(); + } + if (molFractions.at(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; + for (int i=0; i < nComponents; i++) { + molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites; + totalMolecules += nMol.at(i); + } + + if (totalMolecules != nSites) { + sprintf(painCave.errMsg, "Computed total number of molecules is not equal " + "to the number of lattice sites!"); + painCave.isFatal = 1; + simError(); + } + } else { - - -// 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); + for (int i = 0; i < shellRadii.size(); i++) { + if (shellRadii.at(i) > particleRadius + 1e-6 ) { + sprintf(painCave.errMsg, "One of the shellRadius values exceeds the particle Radius."); + painCave.isFatal = 1; + simError(); + } + if (shellRadii.at(i) <= 0.0 ) { + sprintf(painCave.errMsg, "One of the shellRadius values is smaller than zero!"); + painCave.isFatal = 1; + simError(); + } + } } - 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; + vector ids; + if ((int)args_info.molFraction_given){ + sprintf(painCave.errMsg, "Creating a randomized spherical nanoparticle."); + painCave.isFatal = 0; simError(); - } - + /* Random particle is the default case*/ - - - 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 */ + for (int i = 0; i < sites.size(); i++) + if (!isVacancy[i]) ids.push_back(i); + 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){ - sprintf(painCave.errMsg, "Number of .md components " - "does not match the number of shell radius specifications"); - painCave.isFatal = 1; - simError(); - } - } + } else{ + sprintf(painCave.errMsg, "Creating a core-shell spherical nanoparticle."); + painCave.isFatal = 0; + simError(); + RealType smallestSoFar; + int myComponent = -1; + nMol.clear(); + nMol.resize(nComponents); + + for (int i = 0; i < sites.size(); i++) { + myLoc = sites[i]; + myR = myLoc.length(); + smallestSoFar = particleRadius; + if (!isVacancy[i]) { + for (int j = 0; j < nComponents; j++) { + if (myR <= shellRadii[j]) { + if (shellRadii[j] <= smallestSoFar) { + smallestSoFar = shellRadii[j]; + myComponent = j; + } + } + } + componentFromSite[i] = myComponent; + nMol[myComponent]++; + } + } + } - - outputFileName = args_info.output_arg; - - - //creat new .md file on fly which corrects the number of molecule + + //creat new .md file on fly which corrects the number of molecule createMdFile(inputFileName, outputFileName, nMol); if (oldInfo != NULL) delete oldInfo; - - // 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); - - + // Place molecules Molecule* mol; SimInfo::MoleculeIterator mi; mol = NewInfo->beginMolecule(mi); + int l = 0; + int whichSite = 0; 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++; + + if (!args_info.molFraction_given) { + for (int n = 0; n < sites.size(); n++) { + if (!isVacancy[n]) { + if (componentFromSite[n] == i) { + mol = NewInfo->getMoleculeByGlobalIndex(l); + locator->placeMol(sites[n], orientations[n], mol); + l++; + } + } + } + } else { + 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)= 2.0*particleRadius; + hmat(0, 0)= 10.0*particleRadius; hmat(0, 1) = 0.0; hmat(0, 2) = 0.0; hmat(1, 0) = 0.0; - hmat(1, 1) = 2.0*particleRadius; + hmat(1, 1) = 10.0*particleRadius; hmat(1, 2) = 0.0; hmat(2, 0) = 0.0; hmat(2, 1) = 0.0; - hmat(2, 2) = 2.0*particleRadius; + hmat(2, 2) = 10.0*particleRadius; //set Hmat NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat); //create dumpwriter and write out the coordinates - writer = new DumpWriter(NewInfo,outputFileName); + writer = new DumpWriter(NewInfo, outputFileName); if (writer == NULL) { - sprintf(painCave.errMsg, "Error in creating dumpwrite object "); - painCave.isFatal = 1; - simError(); + sprintf(painCave.errMsg, "Error in creating dumpwriter object "); + painCave.isFatal = 1; + simError(); } writer->writeDump(); - std::cout << "new initial configuration file: " << outInitFileName - << " is generated." << std::endl; - // 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 " + sprintf(painCave.errMsg, "A new OpenMD 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, +void createMdFile(const std::string&oldMdFileName, + const std::string&newMdFileName, std::vector nMol) { ifstream oldMdFile; ofstream newMdFile; @@ -351,16 +452,15 @@ void createMdFile(const std::string&oldMdFileName, con //create new .md file based on old .md file oldMdFile.open(oldMdFileName.c_str()); newMdFile.open(newMdFileName.c_str()); - oldMdFile.getline(buffer, MAXLEN); - + int i = 0; while (!oldMdFile.eof()) { - + //correct molecule number if (strstr(buffer, "nMol") != NULL) { if(i