--- trunk/src/applications/randomBuilder/randomBuilder.cpp 2006/10/10 02:44:13 1062 +++ trunk/src/applications/randomBuilder/randomBuilder.cpp 2009/11/25 20:02:06 1390 @@ -5,19 +5,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,12 +28,21 @@ * 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). * + * * randomBuilder.cpp * * Created by Charles F. Vardeman II on 10 Apr 2006. * @author Charles F. Vardeman II - * @version $Id: randomBuilder.cpp,v 1.4 2006-10-10 02:44:13 gezelter Exp $ + * @version $Id: randomBuilder.cpp,v 1.8 2009-11-25 20:01:58 gezelter Exp $ * */ @@ -69,11 +69,11 @@ 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, - int components, int* numMol); + std::vector nMol); int main(int argc, char *argv []) { @@ -86,10 +86,8 @@ int main(int argc, char *argv []) { std::string inputFileName; std::string outputFileName; Lattice *simpleLat; - int* numMol; RealType latticeConstant; std::vector lc; - RealType mass; const RealType rhoConvertConst = 1.661; RealType density; int nx, ny, nz; @@ -97,8 +95,7 @@ int main(int argc, char *argv []) { MoLocator *locator; std::vector latticePos; std::vector latticeOrt; - int numMolPerCell; - int curMolIndex; + int nMolPerCell; DumpWriter *writer; // parse command line arguments @@ -108,7 +105,7 @@ int main(int argc, char *argv []) { density = args_info.density_arg; //get lattice type - latticeType = UpperCase(args_info.latticetype_arg); + latticeType = "FCC"; simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); @@ -118,6 +115,7 @@ int main(int argc, char *argv []) { painCave.isFatal = 1; simError(); } + nMolPerCell = simpleLat->getNumSitesPerCell(); //get the number of unit cells in each direction: @@ -148,6 +146,8 @@ int main(int argc, char *argv []) { simError(); } + int nSites = nMolPerCell * nx * ny * nz; + //get input file name if (args_info.inputs_num) inputFileName = args_info.inputs[0]; @@ -164,44 +164,96 @@ int main(int argc, char *argv []) { SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); Globals* simParams = oldInfo->getSimParams(); - int nComponents =simParams->getNComponents(); - if (oldInfo->getNMoleculeStamp() > 2) { - sprintf(painCave.errMsg, "randomBuilder can't yet build a system with " - "more than two components."); - painCave.isFatal = 1; - simError(); - } + // Calculate lattice constant (in Angstroms) - //get mass of molecule. + std::vector components = simParams->getComponents(); + std::vector molFractions; + std::vector molecularMasses; + std::vector nMol; + int nComponents = components.size(); - mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField()); + 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, "randomBuilder can't figure out molFractions " + "for all of the components in the block."); + painCave.isFatal = 1; + simError(); + } + } - // Create the lattice + // do some sanity checking: + + RealType totalFraction = 0.0; - simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); - - if (simpleLat == NULL) { - sprintf(painCave.errMsg, "Error in creating lattice."); + 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(); } - numMolPerCell = simpleLat->getNumSitesPerCell(); + 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); - // Calculate lattice constant (in Angstroms) + // recompute actual mol fractions and perform final sanity check: - latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density, - (RealType)(1.0 / 3.0)); + 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(); + } + + latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density, + (RealType)(1.0 / 3.0)); + // Set the lattice constant - + lc.push_back(latticeConstant); simpleLat->setLatticeConstant(lc); - - // Calculate the total number of molecules - - int totMol = nx * ny * nz * numMolPerCell; - + // Calculate the lattice sites and fill the lattice vector. // Get the standard orientations of the cell sites @@ -219,37 +271,19 @@ int main(int argc, char *argv []) { simpleLat->getLatticePointsPos(latticePos, i, j, k); - for(int l = 0; l < numMolPerCell; l++) { + for(int l = 0; l < nMolPerCell; l++) { sites.push_back(latticePos[l]); orientations.push_back(latticeOrt[l]); } } } } - - int numSites = sites.size(); - - numMol = new int[nComponents]; - if (nComponents != args_info.molFraction_given && nComponents != 1){ - sprintf(painCave.errMsg, "There needs to be the same number of " - "molFraction arguments as there are components in the " - " block."); - painCave.isFatal = 1; - simError(); - } - int totComponents = 0; - for (int i = 0;i ids; @@ -293,7 +325,7 @@ int main(int argc, char *argv []) { for (int i = 0; i < nComponents; i++){ locator = new MoLocator(newInfo->getMoleculeStamp(i), newInfo->getForceField()); - for (int n = 0; n < numMol[i]; n++) { + for (int n = 0; n < nMol.at(i); n++) { mol = newInfo->getMoleculeByGlobalIndex(l); locator->placeMol(sites[ids[l]], orientations[ids[l]], mol); l++; @@ -316,8 +348,8 @@ int main(int argc, char *argv []) { delete writer; - sprintf(painCave.errMsg, "A new OOPSE MD file called \"%s\" has been " - "generated.", outputFileName.c_str()); + sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been " + "generated.\n", outputFileName.c_str()); painCave.isFatal = 0; simError(); return 0; @@ -325,7 +357,7 @@ void createMdFile(const std::string&oldMdFileName, 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; @@ -343,8 +375,8 @@ void createMdFile(const std::string&oldMdFileName, //correct molecule number if (strstr(buffer, "nMol") != NULL) { - if(i