--- trunk/src/applications/simpleBuilder/simpleBuilder.cpp 2005/04/12 22:42:13 488 +++ trunk/src/applications/simpleBuilder/simpleBuilder.cpp 2014/03/13 13:03:11 1978 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -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,16 @@ * 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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -61,246 +62,230 @@ using namespace std; #include "utils/StringUtils.hpp" using namespace std; -using namespace oopse; -void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName, - int numMol); +using namespace OpenMD; +void createMdFile(const std::string&oldMdFileName, + const std::string&newMdFileName, + int nMol); + int main(int argc, char *argv []) { - //register force fields - registerForceFields(); - registerLattice(); + registerLattice(); - gengetopt_args_info args_info; - std::string latticeType; - std::string inputFileName; - std::string outPrefix; - std::string outMdFileName; - std::string outInitFileName; - Lattice *simpleLat; - int numMol; - double latticeConstant; - std::vector lc; - double mass; - const double rhoConvertConst = 1.661; - double density; - int nx, - ny, - nz; - Mat3x3d hmat; - MoLocator *locator; - std::vector latticePos; - std::vector latticeOrt; - int numMolPerCell; - int curMolIndex; - DumpWriter *writer; + gengetopt_args_info args_info; + std::string latticeType; + std::string inputFileName; + std::string outputFileName; + Lattice *simpleLat; + RealType latticeConstant; + std::vector lc; + const RealType rhoConvertConst = 1.661; + RealType density; + int nx, ny, nz; + Mat3x3d hmat; + MoLocator *locator; + std::vector latticePos; + std::vector latticeOrt; + int nMolPerCell; + DumpWriter *writer; - // parse command line arguments - if (cmdline_parser(argc, argv, &args_info) != 0) - exit(1); + // parse command line arguments + if (cmdline_parser(argc, argv, &args_info) != 0) + exit(1); - density = args_info.density_arg; + density = args_info.density_arg; - //get lattice type - latticeType = UpperCase(args_info.latticetype_arg); + //get lattice type + latticeType = "FCC"; - simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); + simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); - if (simpleLat == NULL) { - sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n", - latticeType.c_str()); - painCave.isFatal = 1; - simError(); - } + if (simpleLat == NULL) { + sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n", + latticeType.c_str()); + painCave.isFatal = 1; + simError(); + } + nMolPerCell = simpleLat->getNumSitesPerCell(); - //get the number of unit cell - nx = args_info.nx_arg; + //get the number of unit cells in each direction: - if (nx <= 0) { - std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl; - exit(1); - } + nx = args_info.nx_arg; - ny = args_info.ny_arg; + if (nx <= 0) { + sprintf(painCave.errMsg, "The number of unit cells in the x direction " + "must be greater than 0."); + painCave.isFatal = 1; + simError(); + } - if (ny <= 0) { - std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl; - exit(1); - } + ny = args_info.ny_arg; - nz = args_info.nz_arg; + if (ny <= 0) { + sprintf(painCave.errMsg, "The number of unit cells in the y direction " + "must be greater than 0."); + painCave.isFatal = 1; + simError(); + } - if (nz <= 0) { - std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl; - exit(1); - } + nz = args_info.nz_arg; - //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; - cmdline_parser_print_help(); - exit(1); - } + if (nz <= 0) { + sprintf(painCave.errMsg, "The number of unit cells in the z direction " + "must be greater than 0."); + painCave.isFatal = 1; + simError(); + } - //parse md file and set up the system - SimCreator oldCreator; - SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); + int nSites = nMolPerCell * nx * ny * nz; - if (oldInfo->getNMoleculeStamp()>= 2) { - std::cerr << "can not build the system with more than two components" - << std::endl; - exit(1); - } + //get input file name + 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"); + painCave.isFatal = 1; + simError(); + } - //get mass of molecule. + //parse md file and set up the system - mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField()); + SimCreator oldCreator; + SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); + Globals* simParams = oldInfo->getSimParams(); - //creat lattice - simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); + // Calculate lattice constant (in Angstroms) - if (simpleLat == NULL) { - std::cerr << "Error in creating lattice" << std::endl; - exit(1); - } + RealType avgMass = MoLocator::getMolMass(oldInfo->getMoleculeStamp(0), + oldInfo->getForceField()); - numMolPerCell = simpleLat->getNumSitesPerCell(); + latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density, + (RealType)(1.0 / 3.0)); + + // Set the lattice constant + + lc.push_back(latticeConstant); + simpleLat->setLatticeConstant(lc); - //calculate lattice constant (in Angstrom) - latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density, - 1.0 / 3.0); + // Calculate the lattice sites and fill the lattice vector. - //set lattice constant - lc.push_back(latticeConstant); - simpleLat->setLatticeConstant(lc); + // Get the standard orientations of the cell sites - //calculate the total number of molecules - numMol = nx * ny * nz * numMolPerCell; + latticeOrt = simpleLat->getLatticePointsOrt(); - if (oldInfo->getNGlobalMolecules() != numMol) { - outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType; - outMdFileName = outPrefix + ".md"; + vector sites; + vector orientations; + + for(int i = 0; i < nx; i++) { + for(int j = 0; j < ny; j++) { + for(int k = 0; k < nz; k++) { - //creat new .md file on fly which corrects the number of molecule - createMdFile(inputFileName, outMdFileName, numMol); - std::cerr - << "SimpleBuilder Error: the number of molecule and the density are not matched" - << std::endl; - std::cerr << "A new .md file: " << outMdFileName - << " is generated, use it to rerun the simpleBuilder" << std::endl; - exit(1); + // Get the position of the cell sites + + simpleLat->getLatticePointsPos(latticePos, i, j, k); + + for(int l = 0; l < nMolPerCell; l++) { + sites.push_back(latticePos[l]); + orientations.push_back(latticeOrt[l]); + } + } } + } + + outputFileName = args_info.output_arg; + + // create a new .md file on the fly which corrects the number of molecules - //determine the output file names - if (args_info.output_given) - outInitFileName = args_info.output_arg; - else - outInitFileName = getPrefix(inputFileName.c_str()) + ".in"; - - //creat Molocator - locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField()); + createMdFile(inputFileName, outputFileName, nSites); - //fill Hmat - hmat(0, 0)= nx * latticeConstant; - hmat(0, 1) = 0.0; - hmat(0, 2) = 0.0; + delete oldInfo; - hmat(1, 0) = 0.0; - hmat(1, 1) = ny * latticeConstant; - hmat(1, 2) = 0.0; + // We need to read in the new SimInfo object, then Parse the + // md file and set up the system - hmat(2, 0) = 0.0; - hmat(2, 1) = 0.0; - hmat(2, 2) = nz * latticeConstant; + SimCreator newCreator; + SimInfo* newInfo = newCreator.createSim(outputFileName, false); - //set Hmat - oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat); + // fill Hmat - //place the molecules + hmat(0, 0) = nx * latticeConstant; + hmat(0, 1) = 0.0; + hmat(0, 2) = 0.0; - curMolIndex = 0; + hmat(1, 0) = 0.0; + hmat(1, 1) = ny * latticeConstant; + hmat(1, 2) = 0.0; - //get the orientation of the cell sites - //for the same type of molecule in same lattice, it will not change - latticeOrt = simpleLat->getLatticePointsOrt(); + hmat(2, 0) = 0.0; + hmat(2, 1) = 0.0; + hmat(2, 2) = nz * latticeConstant; - Molecule* mol; - SimInfo::MoleculeIterator mi; - mol = oldInfo->beginMolecule(mi); - for(int i = 0; i < nx; i++) { - for(int j = 0; j < ny; j++) { - for(int k = 0; k < nz; k++) { + // Set Hmat - //get the position of the cell sites - simpleLat->getLatticePointsPos(latticePos, i, j, k); + newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat); - for(int l = 0; l < numMolPerCell; l++) { - if (mol != NULL) { - locator->placeMol(latticePos[l], latticeOrt[l], mol); - } else { - std::cerr << std::endl; - } - mol = oldInfo->nextMolecule(mi); - } - } - } - } + // place the molecules + + Molecule* mol; + locator = new MoLocator(newInfo->getMoleculeStamp(0), + newInfo->getForceField()); + for (int n = 0; n < nSites; n++) { + mol = newInfo->getMoleculeByGlobalIndex(n); + locator->placeMol(sites[n], orientations[n], mol); + } + + // Create DumpWriter and write out the coordinates - //create dumpwriter and write out the coordinates - oldInfo->setFinalConfigFileName(outInitFileName); - writer = new DumpWriter(oldInfo); + writer = new DumpWriter(newInfo, outputFileName); + + if (writer == NULL) { + sprintf(painCave.errMsg, "error in creating DumpWriter"); + painCave.isFatal = 1; + simError(); + } - if (writer == NULL) { - std::cerr << "error in creating DumpWriter" << std::endl; - exit(1); - } + writer->writeDump(); - 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 objects + delete writer; - //delete oldInfo and oldSimSetup - if (oldInfo != NULL) - delete oldInfo; - - if (writer != NULL) - delete writer; - - delete simpleLat; - - return 0; + sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been " + "generated.\n", outputFileName.c_str()); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + return 0; } -void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName, - int numMol) { - ifstream oldMdFile; - ofstream newMdFile; - const int MAXLEN = 65535; - char buffer[MAXLEN]; +void createMdFile(const std::string&oldMdFileName, + const std::string&newMdFileName, + int nMol) { + ifstream oldMdFile; + ofstream newMdFile; + const int MAXLEN = 65535; + char buffer[MAXLEN]; - //create new .md file based on old .md file - oldMdFile.open(oldMdFileName.c_str()); - newMdFile.open(newMdFileName.c_str()); + //create new .md file based on old .md file + oldMdFile.open(oldMdFileName.c_str()); + newMdFile.open(newMdFileName.c_str()); - oldMdFile.getline(buffer, MAXLEN); + oldMdFile.getline(buffer, MAXLEN); - while (!oldMdFile.eof()) { + while (!oldMdFile.eof()) { - //correct molecule number - if (strstr(buffer, "nMol") != NULL) { - sprintf(buffer, "\t\tnMol = %d;", numMol); - newMdFile << buffer << std::endl; - } else - newMdFile << buffer << std::endl; + //correct molecule number + if (strstr(buffer, "nMol") != NULL) { + sprintf(buffer, "\t\tnMol = %d;", nMol); + newMdFile << buffer << std::endl; + } else + newMdFile << buffer << std::endl; - oldMdFile.getline(buffer, MAXLEN); - } + oldMdFile.getline(buffer, MAXLEN); + } - oldMdFile.close(); - newMdFile.close(); + oldMdFile.close(); + newMdFile.close(); }