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
#	Content
1	#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	#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
24	using namespace std;
25	using namespace oopse;
26	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
27	int numMol);
28
29	int main(int argc, char *argv []) {
30
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	gengetopt_args_info args_info;
38	std::string latticeType;
39	std::string inputFileName;
40	std::string outPrefix;
41	std::string outMdFileName;
42	std::string outInitFileName;
43	BaseLattice *simpleLat;
44	int numMol;
45	double latticeConstant;
46	std::vector<double> lc;
47	double mass;
48	const double rhoConvertConst = 1.661;
49	double density;
50	int nx,
51	ny,
52	nz;
53	Mat3x3d hmat;
54	MoLocator *locator;
55	std::vector<Vector3d> latticePos;
56	std::vector<Vector3d> latticeOrt;
57	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	std::cerr << latticeType << " is an invalid lattice type" << std::endl;
72	std::cerr << LatticeFactory::getInstance()->toString() << std::endl;
73	exit(1);
74	}
75
76	//get the number of unit cell
77	nx = args_info.nx_arg;
78
79	if (nx <= 0) {
80	std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl;
81	exit(1);
82	}
83
84	ny = args_info.ny_arg;
85
86	if (ny <= 0) {
87	std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl;
88	exit(1);
89	}
90
91	nz = args_info.nz_arg;
92
93	if (nz <= 0) {
94	std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl;
95	exit(1);
96	}
97
98	//get input file name
99	if (args_info.inputs_num)
100	inputFileName = args_info.inputs[0];
101	else {
102	std::cerr << "You must specify a input file name.\n" << std::endl;
103	cmdline_parser_print_help();
104	exit(1);
105	}
106
107	//parse md file and set up the system
108	SimCreator oldCreator;
109	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
110
111	if (oldInfo->getNMoleculeStamp()>= 2) {
112	std::cerr << "can not build the system with more than two components"
113	<< std::endl;
114	exit(1);
115	}
116
117	//get mass of molecule.
118
119	mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
120
121	//creat lattice
122	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
123
124	if (simpleLat == NULL) {
125	std::cerr << "Error in creating lattice" << std::endl;
126	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	if (oldInfo->getNGlobalMolecules() != numMol) {
143	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	std::cerr
149	<< "SimpleBuilder Error: the number of molecule and the density are not matched"
150	<< std::endl;
151	std::cerr << "A new .md file: " << outMdFileName
152	<< " is generated, use it to rerun the simpleBuilder" << std::endl;
153	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
162	//creat Molocator
163	locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
164
165	//fill Hmat
166	hmat(0, 0)= nx * latticeConstant;
167	hmat(0, 1) = 0.0;
168	hmat(0, 2) = 0.0;
169
170	hmat(1, 0) = 0.0;
171	hmat(1, 1) = ny * latticeConstant;
172	hmat(1, 2) = 0.0;
173
174	hmat(2, 0) = 0.0;
175	hmat(2, 1) = 0.0;
176	hmat(2, 2) = nz * latticeConstant;
177
178	//set Hmat
179	oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
180
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	Molecule* mol;
190	SimInfo::MoleculeIterator mi;
191	mol = oldInfo->beginMolecule(mi);
192	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	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	}
208	}
209	}
210
211	//create dumpwriter and write out the coordinates
212	oldInfo->setFinalConfigFileName(outInitFileName);
213	writer = new DumpWriter(oldInfo, oldInfo->getDumpFileName());
214
215	if (writer == NULL) {
216	std::cerr << "error in creating DumpWriter" << std::endl;
217	exit(1);
218	}
219
220	writer->writeDump();
221	std::cout << "new initial configuration file: " << outInitFileName
222	<< " is generated." << std::endl;
223
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	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
237	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	newMdFile << buffer << std::endl;
255	} else
256	newMdFile << buffer << std::endl;
257
258	oldMdFile.getline(buffer, MAXLEN);
259	}
260
261	oldMdFile.close();
262	newMdFile.close();
263	}
264