--- trunk/src/applications/simpleBuilder/simpleBuilder.cpp 2004/09/28 23:24:25 12 +++ trunk/src/applications/simpleBuilder/simpleBuilder.cpp 2014/03/13 13:03:11 1978 @@ -1,3 +1,45 @@ +/* + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * + * The University of Notre Dame grants you ("Licensee") a + * non-exclusive, royalty free, license to use, modify and + * redistribute this software in source and binary code form, provided + * that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 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. + * + * This software is provided "AS IS," without a warranty of any + * kind. All express or implied conditions, representations and + * warranties, including any implied warranty of merchantability, + * fitness for a particular purpose or non-infringement, are hereby + * excluded. The University of Notre Dame and its licensors shall not + * be liable for any damages suffered by licensee as a result of + * using, modifying or distributing the software or its + * derivatives. In no event will the University of Notre Dame or its + * licensors be liable for any lost revenue, profit or data, or for + * direct, indirect, special, consequential, incidental or punitive + * damages, however caused and regardless of the theory of liability, + * 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 #include #include @@ -7,224 +49,219 @@ #include #include -#include "io/Globals.hpp" -#include "brains/SimInfo.hpp" -#include "brains/SimSetup.hpp" #include "applications/simpleBuilder/simpleBuilderCmd.h" +#include "lattice/LatticeFactory.hpp" +#include "utils/MoLocator.hpp" +#include "lattice/Lattice.hpp" +#include "brains/Register.hpp" +#include "brains/SimInfo.hpp" +#include "brains/SimCreator.hpp" +#include "io/DumpWriter.hpp" +#include "math/Vector3.hpp" +#include "math/SquareMatrix3.hpp" #include "utils/StringUtils.hpp" -#include "applications/simpleBuilder/LatticeFactory.hpp" -#include "applications/simpleBuilder/Vector3d.hpp" -#include "applications/simpleBuilder/MoLocator.hpp" -#include "applications/simpleBuilder/Lattice.hpp" using namespace std; +using namespace OpenMD; -void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol); -double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF); +void createMdFile(const std::string&oldMdFileName, + const std::string&newMdFileName, + int nMol); -int main( int argc, char* argv[]){ +int main(int argc, char *argv []) { + registerLattice(); + gengetopt_args_info args_info; - string latticeType; - string inputFileName; - string outPrefix; - string outMdFileName; - string outInitFileName; - SimInfo* oldInfo; - SimSetup* oldSimSetup; - BaseLattice* simpleLat; - int numMol; - double latticeConstant; - vector lc; - double mass; - const double rhoConvertConst = 1.661; - double density; + 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; - double Hmat[3][3]; + Mat3x3d hmat; MoLocator *locator; - vector latticePos; - vector latticeOrt; - int numMolPerCell; - int curMolIndex; - DumpWriter* writer; - + std::vector latticePos; + std::vector latticeOrt; + int nMolPerCell; + DumpWriter *writer; + // parse command line arguments - if (cmdline_parser (argc, argv, &args_info) != 0) - exit(1) ; - + if (cmdline_parser(argc, argv, &args_info) != 0) + exit(1); + density = args_info.density_arg; //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()->toString() << endl; - exit(1); + latticeType = "FCC"; + + 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(); } + nMolPerCell = simpleLat->getNumSitesPerCell(); - //get the number of unit cell + //get the number of unit cells in each direction: + 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); + + if (nx <= 0) { + sprintf(painCave.errMsg, "The number of unit cells in the x direction " + "must be greater than 0."); + painCave.isFatal = 1; + simError(); } 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); + + if (ny <= 0) { + sprintf(painCave.errMsg, "The number of unit cells in the y direction " + "must be greater than 0."); + painCave.isFatal = 1; + simError(); } 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); + + if (nz <= 0) { + sprintf(painCave.errMsg, "The number of unit cells in the z direction " + "must be greater than 0."); + painCave.isFatal = 1; + simError(); } - + + int nSites = nMolPerCell * nx * ny * nz; + //get input file name - if (args_info.inputs_num) + if (args_info.inputs_num) inputFileName = args_info.inputs[0]; - else { - cerr <<"You must specify a input file name.\n" << endl; - cmdline_parser_print_help(); - exit(1); + else { + sprintf(painCave.errMsg, "No input .md file name was specified " + "on the command line"); + painCave.isFatal = 1; + simError(); } - //parse md file and set up the system - oldInfo = new SimInfo; - if(oldInfo == NULL){ - cerr << "error in creating SimInfo" << endl; - exit(1); - } - oldSimSetup = new SimSetup(); - if(oldSimSetup == NULL){ - cerr << "error in creating SimSetup" << endl; - exit(1); - } + SimCreator oldCreator; + SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); + Globals* simParams = oldInfo->getSimParams(); - oldSimSetup->suspendInit(); - oldSimSetup->setSimInfo(oldInfo ); - oldSimSetup->parseFile(&inputFileName[0] ); - 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){ - cerr << "Error in creating lattice" << endl; - exit(1); - } + // Calculate lattice constant (in Angstroms) - numMolPerCell = simpleLat->getNumSitesPerCell(); + RealType avgMass = MoLocator::getMolMass(oldInfo->getMoleculeStamp(0), + oldInfo->getForceField()); + + latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density, + (RealType)(1.0 / 3.0)); - //calculate lattice constant (in Angstrom) - latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0); + // Set the lattice constant - //set lattice constant lc.push_back(latticeConstant); simpleLat->setLatticeConstant(lc); - - //calculate the total number of molecules - numMol = nx * ny * nz * numMolPerCell; - if (oldInfo->n_mol != numMol){ + // Calculate the lattice sites and fill the lattice vector. - outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType; - outMdFileName = outPrefix + ".md"; + // Get the standard orientations of the cell sites - //creat new .md file on fly which corrects the number of molecule - createMdFile(inputFileName, outMdFileName, numMol); - cerr << "SimpleBuilder Error: the number of molecule and the density are not matched" <getLatticePointsOrt(); + + 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++) { + + // 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]); + } + } + } } - //determine the output file names - if (args_info.output_given) - outInitFileName = args_info.output_arg; - else - outInitFileName = getPrefix(inputFileName.c_str()) + ".in"; + outputFileName = args_info.output_arg; - - //allocat memory for storing pos, vel and etc - oldInfo->getConfiguration()->createArrays(oldInfo->n_atoms); - for (int i = 0; i < oldInfo->n_atoms; i++) - oldInfo->atoms[i]->setCoords(); + // create a new .md file on the fly which corrects the number of molecules - //creat Molocator - locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->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->setBoxM(Hmat); - - //place the molecules + // fill Hmat - curMolIndex = 0; + hmat(0, 0) = nx * latticeConstant; + hmat(0, 1) = 0.0; + hmat(0, 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(1, 0) = 0.0; + hmat(1, 1) = ny * latticeConstant; + hmat(1, 2) = 0.0; - for(int i =0; i < nx; i++){ - for(int j=0; j < ny; j++){ - for(int k = 0; k < nz; k++){ + hmat(2, 0) = 0.0; + hmat(2, 1) = 0.0; + hmat(2, 2) = nz * latticeConstant; - //get the position of the cell sites - simpleLat->getLatticePointsPos(latticePos, i, j, k); + // Set Hmat - for(int l = 0; l < numMolPerCell; l++) - locator->placeMol(latticePos[l], latticeOrt[l], &(oldInfo->molecules[curMolIndex++])); - } - } - } + newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat); - //create dumpwriter and write out the coordinates - oldInfo->finalName = outInitFileName; - writer = new DumpWriter( oldInfo ); - if(writer == NULL){ - cerr << "error in creating DumpWriter" << endl; - exit(1); + // 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); } - writer->writeFinal(0); - cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl; - //delete objects + + // Create DumpWriter and write out the coordinates - //delete oldInfo and oldSimSetup - if(oldInfo != NULL) - delete oldInfo; + writer = new DumpWriter(newInfo, outputFileName); - if(oldSimSetup != NULL) - delete oldSimSetup; - - if (writer != NULL) - delete writer; + if (writer == NULL) { + sprintf(painCave.errMsg, "error in creating DumpWriter"); + painCave.isFatal = 1; + simError(); + } + + writer->writeDump(); + + // deleting the writer will put the closing at the end of the dump file. + + delete writer; + + 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 string& oldMdFileName, const string& newMdFileName, int numMol){ +void createMdFile(const std::string&oldMdFileName, + const std::string&newMdFileName, + int nMol) { ifstream oldMdFile; ofstream newMdFile; const int MAXLEN = 65535; @@ -235,36 +272,20 @@ void createMdFile(const string& oldMdFileName, const s newMdFile.open(newMdFileName.c_str()); 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 << endl; - } - else - newMdFile << buffer << endl; + 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.close(); 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; -}