--- trunk/src/applications/simpleBuilder/simpleBuilder.cpp 2004/09/28 23:24:25 12 +++ trunk/src/applications/simpleBuilder/simpleBuilder.cpp 2005/04/15 22:04:00 507 @@ -1,3 +1,44 @@ +/* + * 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. 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 + * notice, this list of conditions and the following disclaimer. + * + * 3. 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. + */ + #include #include #include @@ -7,179 +48,176 @@ #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 oopse; +void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName, + int numMol); -void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol); -double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF); +int main(int argc, char *argv []) { -int main( int argc, char* argv[]){ - + //register force fields + registerForceFields(); + registerLattice(); + gengetopt_args_info args_info; - string latticeType; - string inputFileName; - string outPrefix; - string outMdFileName; - string outInitFileName; - SimInfo* oldInfo; - SimSetup* oldSimSetup; - BaseLattice* simpleLat; + std::string latticeType; + std::string inputFileName; + std::string outPrefix; + std::string outMdFileName; + std::string outInitFileName; + Lattice *simpleLat; int numMol; double latticeConstant; - vector lc; + std::vector lc; double mass; const double rhoConvertConst = 1.661; double density; - int nx, ny, nz; - double Hmat[3][3]; + int nx, + ny, + nz; + Mat3x3d hmat; MoLocator *locator; - vector latticePos; - vector latticeOrt; + std::vector latticePos; + std::vector latticeOrt; int numMolPerCell; int curMolIndex; - DumpWriter* writer; - + 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); + + 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(); } //get the number of unit cell nx = args_info.nx_arg; - if(nx <= 0){ - cerr << "The number of unit cell in h direction must be greater than 0" << endl; + + if (nx <= 0) { + std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl; exit(1); } ny = args_info.ny_arg; - if(ny <= 0){ - cerr << "The number of unit cell in l direction must be greater than 0" << endl; + + if (ny <= 0) { + std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl; exit(1); } nz = args_info.nz_arg; - if(nz <= 0){ - cerr << "The number of unit cell in k direction must be greater than 0" << endl; + + if (nz <= 0) { + std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl; exit(1); } - + //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; + else { + std::cerr << "You must specify a input file name.\n" << std::endl; cmdline_parser_print_help(); exit(1); } - //parse md file and set up the system - oldInfo = new SimInfo; - if(oldInfo == NULL){ - cerr << "error in creating SimInfo" << endl; - exit(1); - } + SimCreator oldCreator; + SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); - oldSimSetup = new SimSetup(); - if(oldSimSetup == NULL){ - cerr << "error in creating SimSetup" << endl; - exit(1); + if (oldInfo->getNMoleculeStamp()>= 2) { + std::cerr << "can not build the system with more than two components" + << std::endl; + exit(1); } - 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()); - + + mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField()); + //creat lattice - simpleLat = LatticeFactory::getInstance()->createLattice(latticeType); - if(simpleLat == NULL){ - cerr << "Error in creating lattice" << endl; - exit(1); - } + 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 = pow(rhoConvertConst * numMolPerCell * mass / density, + 1.0 / 3.0); + //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){ - + if (oldInfo->getNGlobalMolecules() != numMol) { outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType; outMdFileName = outPrefix + ".md"; //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" <getConfiguration()->createArrays(oldInfo->n_atoms); - for (int i = 0; i < oldInfo->n_atoms; i++) - oldInfo->atoms[i]->setCoords(); - + outInitFileName = getPrefix(inputFileName.c_str()) + ".in"; + //creat Molocator - locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->getForceField()); + locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField()); //fill Hmat - Hmat[0][0] = nx * latticeConstant; - Hmat[0][1] = 0.0; - Hmat[0][2] = 0.0; + hmat(0, 0)= nx * latticeConstant; + hmat(0, 1) = 0.0; + hmat(0, 2) = 0.0; - Hmat[1][0] = 0.0; - Hmat[1][1] = ny * latticeConstant; - Hmat[1][2] = 0.0; + hmat(1, 0) = 0.0; + hmat(1, 1) = ny * latticeConstant; + hmat(1, 2) = 0.0; - Hmat[2][0] = 0.0; - Hmat[2][1] = 0.0; - Hmat[2][2] = nz * latticeConstant ; + hmat(2, 0) = 0.0; + hmat(2, 1) = 0.0; + hmat(2, 2) = nz * latticeConstant; //set Hmat - oldInfo->setBoxM(Hmat); - + oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat); + //place the molecules curMolIndex = 0; @@ -188,43 +226,57 @@ int main( int argc, char* argv[]){ //for the same type of molecule in same lattice, it will not change latticeOrt = simpleLat->getLatticePointsOrt(); - for(int i =0; i < nx; i++){ - for(int j=0; j < ny; j++){ - for(int k = 0; k < nz; k++){ + 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++) { - //get the position of the cell sites - simpleLat->getLatticePointsPos(latticePos, i, j, k); + //get the position of the cell sites + simpleLat->getLatticePointsPos(latticePos, i, j, k); - for(int l = 0; l < numMolPerCell; l++) - locator->placeMol(latticePos[l], latticeOrt[l], &(oldInfo->molecules[curMolIndex++])); - } + 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); + } + } } } //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); + oldInfo->setFinalConfigFileName(outInitFileName); + writer = new DumpWriter(oldInfo); + + if (writer == NULL) { + std::cerr << "error in creating DumpWriter" << std::endl; + exit(1); } - writer->writeFinal(0); - cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl; + + writer->writeDump(); + std::cout << "new initial configuration file: " << outInitFileName + << " is generated." << std::endl; + //delete objects //delete oldInfo and oldSimSetup - if(oldInfo != NULL) - delete oldInfo; - - if(oldSimSetup != NULL) - delete oldSimSetup; - + if (oldInfo != NULL) + delete oldInfo; + if (writer != NULL) delete writer; + + delete simpleLat; + return 0; } -void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol){ +void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName, + int numMol) { ifstream oldMdFile; ofstream newMdFile; const int MAXLEN = 65535; @@ -235,36 +287,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){ + if (strstr(buffer, "nMol") != NULL) { sprintf(buffer, "\t\tnMol = %d;", numMol); - newMdFile << buffer << endl; - } - else - newMdFile << buffer << endl; + 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; -}