applications/simpleBuilder/simpleBuilder.cpp

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <string>
#include <map>
#include <fstream>

#include "applications/simpleBuilder/simpleBuilderCmd.h"
#include "applications/simpleBuilder/LatticeFactory.hpp"
#include "applications/simpleBuilder/MoLocator.hpp"
#include "applications/simpleBuilder/Lattice.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 "UseTheForce/DUFF.hpp"
#include "UseTheForce/EAM.hpp"
#include "UseTheForce/ForceFieldCreator.hpp"

using namespace std;
using namespace oopse;
void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
                  int numMol);

int main(int argc, char *argv []) {

    ForceFieldBuilder<DUFF> DUFFCreator("DUFF");
    ForceFieldBuilder<DUFF> WATERCreator("WATER");
    ForceFieldBuilder<DUFF> LJCreator("LJ");
    //in theory, EAM can also be merged
    ForceFieldBuilder<EAM> EAMCreator("EAM");
    
    gengetopt_args_info args_info;
    std::string latticeType;
    std::string inputFileName;
    std::string outPrefix;
    std::string outMdFileName;
    std::string outInitFileName;
    BaseLattice *simpleLat;
    int numMol;
    double latticeConstant;
    std::vector<double> lc;
    double mass;
    const double rhoConvertConst = 1.661;
    double density;
    int nx,
    ny,
    nz;
    Mat3x3d hmat;
    MoLocator *locator;
    std::vector<Vector3d> latticePos;
    std::vector<Vector3d> latticeOrt;
    int numMolPerCell;
    int curMolIndex;
    DumpWriter *writer;

    // parse command line arguments
    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)) {
        std::cerr << latticeType << " is an invalid lattice type" << std::endl;
        std::cerr << LatticeFactory::getInstance()->toString() << std::endl;
        exit(1);
    }

    //get the number of unit cell
    nx = args_info.nx_arg;

    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) {
        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) {
        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)
        inputFileName = args_info.inputs[0];
    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
    SimCreator oldCreator;
    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);

    if (oldInfo->getNMoleculeStamp()>= 2) {
        std::cerr << "can not build the system with more than two components"
            << std::endl;
        exit(1);
    }

    //get mass of molecule. 

    mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());

    //creat lattice
    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);

    //set lattice constant
    lc.push_back(latticeConstant);
    simpleLat->setLatticeConstant(lc);

    //calculate the total number of molecules
    numMol = nx * ny * nz * numMolPerCell;

    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);
        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);
    }

    //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());

    //fill Hmat
    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(2, 0) = 0.0;
    hmat(2, 1) = 0.0;
    hmat(2, 2) = nz * latticeConstant;

    //set Hmat
    oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);

    //place the molecules

    curMolIndex = 0;

    //get the orientation of the cell sites
    //for the same type of molecule in same lattice, it will not change
    latticeOrt = simpleLat->getLatticePointsOrt();

    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);

                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->setFinalConfigFileName(outInitFileName);
    writer = new DumpWriter(oldInfo, oldInfo->getDumpFileName());

    if (writer == NULL) {
        std::cerr << "error in creating DumpWriter" << std::endl;
        exit(1);
    }

    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 (writer != NULL)
        delete writer;

    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];

    //create new .md file based on old .md file
    oldMdFile.open(oldMdFileName.c_str());
    newMdFile.open(newMdFileName.c_str());

    oldMdFile.getline(buffer, MAXLEN);

    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;

        oldMdFile.getline(buffer, MAXLEN);
    }

    oldMdFile.close();
    newMdFile.close();
}

Revision:	1903
Committed:	Thu Jan 6 00:16:07 2005 UTC (19 years, 6 months ago) by tim
File size:	7568 byte(s)
Log Message:	simpleBuilder in progress
#	User	Rev	Content
1	tim	1770	#include <cstdlib>
2			#include <cstdio>
3			#include <cstring>
4			#include <cmath>
5			#include <iostream>
6			#include <string>
7			#include <map>
8			#include <fstream>
9
10			#include "applications/simpleBuilder/simpleBuilderCmd.h"
11			#include "applications/simpleBuilder/LatticeFactory.hpp"
12			#include "applications/simpleBuilder/MoLocator.hpp"
13			#include "applications/simpleBuilder/Lattice.hpp"
14	tim	1903	#include "brains/SimInfo.hpp"
15			#include "brains/SimCreator.hpp"
16			#include "io/DumpWriter.hpp"
17			#include "math/Vector3.hpp"
18			#include "math/SquareMatrix3.hpp"
19			#include "utils/StringUtils.hpp"
20			#include "UseTheForce/DUFF.hpp"
21			#include "UseTheForce/EAM.hpp"
22			#include "UseTheForce/ForceFieldCreator.hpp"
23	tim	1770
24	tim	1903	using namespace std;
25	tim	1891	using namespace oopse;
26	tim	1832	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
27	tim	1770	int numMol);
28
29			int main(int argc, char *argv []) {
30	tim	1903
31			ForceFieldBuilder<DUFF> DUFFCreator("DUFF");
32			ForceFieldBuilder<DUFF> WATERCreator("WATER");
33			ForceFieldBuilder<DUFF> LJCreator("LJ");
34			//in theory, EAM can also be merged
35			ForceFieldBuilder<EAM> EAMCreator("EAM");
36
37	tim	1770	gengetopt_args_info args_info;
38	tim	1903	std::string latticeType;
39			std::string inputFileName;
40			std::string outPrefix;
41			std::string outMdFileName;
42			std::string outInitFileName;
43	tim	1770	BaseLattice *simpleLat;
44			int numMol;
45			double latticeConstant;
46	tim	1903	std::vector<double> lc;
47	tim	1770	double mass;
48			const double rhoConvertConst = 1.661;
49			double density;
50			int nx,
51			ny,
52			nz;
53	tim	1903	Mat3x3d hmat;
54	tim	1770	MoLocator *locator;
55	tim	1903	std::vector<Vector3d> latticePos;
56			std::vector<Vector3d> latticeOrt;
57	tim	1770	int numMolPerCell;
58			int curMolIndex;
59			DumpWriter *writer;
60
61			// parse command line arguments
62			if (cmdline_parser(argc, argv, &args_info) != 0)
63			exit(1);
64
65			density = args_info.density_arg;
66
67			//get lattice type
68			latticeType = UpperCase(args_info.latticetype_arg);
69
70			if (!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)) {
71	tim	1830	std::cerr << latticeType << " is an invalid lattice type" << std::endl;
72			std::cerr << LatticeFactory::getInstance()->toString() << std::endl;
73	tim	1770	exit(1);
74			}
75
76			//get the number of unit cell
77			nx = args_info.nx_arg;
78
79			if (nx <= 0) {
80	tim	1903	std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl;
81	tim	1770	exit(1);
82			}
83
84			ny = args_info.ny_arg;
85
86			if (ny <= 0) {
87	tim	1903	std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl;
88	tim	1770	exit(1);
89			}
90
91			nz = args_info.nz_arg;
92
93			if (nz <= 0) {
94	tim	1903	std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl;
95	tim	1770	exit(1);
96			}
97
98			//get input file name
99			if (args_info.inputs_num)
100			inputFileName = args_info.inputs[0];
101			else {
102	tim	1830	std::cerr << "You must specify a input file name.\n" << std::endl;
103	tim	1770	cmdline_parser_print_help();
104			exit(1);
105			}
106
107			//parse md file and set up the system
108	tim	1903	SimCreator oldCreator;
109			SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
110	tim	1770
111	tim	1903	if (oldInfo->getNMoleculeStamp()>= 2) {
112	tim	1830	std::cerr << "can not build the system with more than two components"
113			<< std::endl;
114	tim	1770	exit(1);
115			}
116
117			//get mass of molecule.
118
119	tim	1903	mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
120
121	tim	1770	//creat lattice
122			simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
123
124			if (simpleLat == NULL) {
125	tim	1830	std::cerr << "Error in creating lattice" << std::endl;
126	tim	1770	exit(1);
127			}
128
129			numMolPerCell = simpleLat->getNumSitesPerCell();
130
131			//calculate lattice constant (in Angstrom)
132			latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
133			1.0 / 3.0);
134
135			//set lattice constant
136			lc.push_back(latticeConstant);
137			simpleLat->setLatticeConstant(lc);
138
139			//calculate the total number of molecules
140			numMol = nx * ny * nz * numMolPerCell;
141
142	tim	1903	if (oldInfo->getNGlobalMolecules() != numMol) {
143	tim	1770	outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
144			outMdFileName = outPrefix + ".md";
145
146			//creat new .md file on fly which corrects the number of molecule
147			createMdFile(inputFileName, outMdFileName, numMol);
148	tim	1830	std::cerr
149	tim	1770	<< "SimpleBuilder Error: the number of molecule and the density are not matched"
150	tim	1830	<< std::endl;
151			std::cerr << "A new .md file: " << outMdFileName
152			<< " is generated, use it to rerun the simpleBuilder" << std::endl;
153	tim	1770	exit(1);
154			}
155
156			//determine the output file names
157			if (args_info.output_given)
158			outInitFileName = args_info.output_arg;
159			else
160			outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
161	tim	1903
162	tim	1770	//creat Molocator
163	tim	1903	locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
164	tim	1770
165			//fill Hmat
166	tim	1903	hmat(0, 0)= nx * latticeConstant;
167			hmat(0, 1) = 0.0;
168			hmat(0, 2) = 0.0;
169	tim	1770
170	tim	1903	hmat(1, 0) = 0.0;
171			hmat(1, 1) = ny * latticeConstant;
172			hmat(1, 2) = 0.0;
173	tim	1770
174	tim	1903	hmat(2, 0) = 0.0;
175			hmat(2, 1) = 0.0;
176			hmat(2, 2) = nz * latticeConstant;
177	tim	1770
178			//set Hmat
179	tim	1903	oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
180	tim	1770
181			//place the molecules
182
183			curMolIndex = 0;
184
185			//get the orientation of the cell sites
186			//for the same type of molecule in same lattice, it will not change
187			latticeOrt = simpleLat->getLatticePointsOrt();
188
189	tim	1903	Molecule* mol;
190			SimInfo::MoleculeIterator mi;
191			mol = oldInfo->beginMolecule(mi);
192	tim	1770	for(int i = 0; i < nx; i++) {
193			for(int j = 0; j < ny; j++) {
194			for(int k = 0; k < nz; k++) {
195
196			//get the position of the cell sites
197			simpleLat->getLatticePointsPos(latticePos, i, j, k);
198
199	tim	1903	for(int l = 0; l < numMolPerCell; l++) {
200			if (mol != NULL) {
201			locator->placeMol(latticePos[l], latticeOrt[l], mol);
202			} else {
203			std::cerr << std::endl;
204			}
205			mol = oldInfo->nextMolecule(mi);
206			}
207	tim	1770	}
208			}
209			}
210
211			//create dumpwriter and write out the coordinates
212	tim	1903	oldInfo->setFinalConfigFileName(outInitFileName);
213			writer = new DumpWriter(oldInfo, oldInfo->getDumpFileName());
214	tim	1770
215			if (writer == NULL) {
216	tim	1830	std::cerr << "error in creating DumpWriter" << std::endl;
217	tim	1770	exit(1);
218			}
219
220	tim	1903	writer->writeDump();
221			std::cout << "new initial configuration file: " << outInitFileName
222	tim	1830	<< " is generated." << std::endl;
223	tim	1770
224			//delete objects
225
226			//delete oldInfo and oldSimSetup
227			if (oldInfo != NULL)
228			delete oldInfo;
229
230			if (writer != NULL)
231			delete writer;
232
233			return 0;
234			}
235
236	tim	1832	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
237	tim	1770	int numMol) {
238			ifstream oldMdFile;
239			ofstream newMdFile;
240			const int MAXLEN = 65535;
241			char buffer[MAXLEN];
242
243			//create new .md file based on old .md file
244			oldMdFile.open(oldMdFileName.c_str());
245			newMdFile.open(newMdFileName.c_str());
246
247			oldMdFile.getline(buffer, MAXLEN);
248
249			while (!oldMdFile.eof()) {
250
251			//correct molecule number
252			if (strstr(buffer, "nMol") != NULL) {
253			sprintf(buffer, "\t\tnMol = %d;", numMol);
254	tim	1830	newMdFile << buffer << std::endl;
255	tim	1770	} else
256	tim	1830	newMdFile << buffer << std::endl;
257	tim	1770
258			oldMdFile.getline(buffer, MAXLEN);
259			}
260
261			oldMdFile.close();
262			newMdFile.close();
263			}
264