ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/branches/development/src/applications/simpleBuilder/simpleBuilder.cpp

Comparing:
trunk/src/applications/simpleBuilder/simpleBuilder.cpp (file contents), Revision 376 by tim, Thu Feb 24 20:55:07 2005 UTC vs.
branches/development/src/applications/simpleBuilder/simpleBuilder.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

#	Line 1 | Line 1
1	<	/*
1	>	/*
2		* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	+	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include <cstdlib>
#	Line 49 | Line 50
50		#include <fstream>
51
52		#include "applications/simpleBuilder/simpleBuilderCmd.h"
53	<	#include "applications/simpleBuilder/LatticeFactory.hpp"
54	<	#include "applications/simpleBuilder/MoLocator.hpp"
55	<	#include "applications/simpleBuilder/Lattice.hpp"
53	>	#include "lattice/LatticeFactory.hpp"
54	>	#include "utils/MoLocator.hpp"
55	>	#include "lattice/Lattice.hpp"
56		#include "brains/Register.hpp"
57		#include "brains/SimInfo.hpp"
58		#include "brains/SimCreator.hpp"
#	Line 61 | Line 62 | using namespace std;
62		#include "utils/StringUtils.hpp"
63
64		using namespace std;
65	<	using namespace oopse;
65	<	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
66	<	int numMol);
65	>	using namespace OpenMD;
66
67	+	void createMdFile(const std::string&oldMdFileName,
68	+	const std::string&newMdFileName,
69	+	int nMol);
70	+
71		int main(int argc, char *argv []) {
72
73	<	//register force fields
74	<	registerForceFields();
73	>	// register force fields
74	>	registerForceFields();
75	>	registerLattice();
76
77	<	gengetopt_args_info args_info;
78	<	std::string latticeType;
79	<	std::string inputFileName;
80	<	std::string outPrefix;
81	<	std::string outMdFileName;
82	<	std::string outInitFileName;
83	<	BaseLattice *simpleLat;
84	<	int numMol;
85	<	double latticeConstant;
86	<	std::vector<double> lc;
87	<	double mass;
88	<	const double rhoConvertConst = 1.661;
89	<	double density;
90	<	int nx,
91	<	ny,
92	<	nz;
89	<	Mat3x3d hmat;
90	<	MoLocator *locator;
91	<	std::vector<Vector3d> latticePos;
92	<	std::vector<Vector3d> latticeOrt;
93	<	int numMolPerCell;
94	<	int curMolIndex;
95	<	DumpWriter *writer;
77	>	gengetopt_args_info args_info;
78	>	std::string latticeType;
79	>	std::string inputFileName;
80	>	std::string outputFileName;
81	>	Lattice *simpleLat;
82	>	RealType latticeConstant;
83	>	std::vector<RealType> lc;
84	>	const RealType rhoConvertConst = 1.661;
85	>	RealType density;
86	>	int nx, ny, nz;
87	>	Mat3x3d hmat;
88	>	MoLocator *locator;
89	>	std::vector<Vector3d> latticePos;
90	>	std::vector<Vector3d> latticeOrt;
91	>	int nMolPerCell;
92	>	DumpWriter *writer;
93
94	<	// parse command line arguments
95	<	if (cmdline_parser(argc, argv, &args_info) != 0)
96	<	exit(1);
94	>	// parse command line arguments
95	>	if (cmdline_parser(argc, argv, &args_info) != 0)
96	>	exit(1);
97
98	<	density = args_info.density_arg;
98	>	density = args_info.density_arg;
99
100	<	//get lattice type
101	<	latticeType = UpperCase(args_info.latticetype_arg);
100	>	//get lattice type
101	>	latticeType = "FCC";
102
103	<	if (!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)) {
104	<	std::cerr << latticeType << " is an invalid lattice type" << std::endl;
105	<	std::cerr << LatticeFactory::getInstance()->toString() << std::endl;
106	<	exit(1);
107	<	}
103	>	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
104	>
105	>	if (simpleLat == NULL) {
106	>	sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n",
107	>	latticeType.c_str());
108	>	painCave.isFatal = 1;
109	>	simError();
110	>	}
111	>	nMolPerCell = simpleLat->getNumSitesPerCell();
112
113	<	//get the number of unit cell
113	<	nx = args_info.nx_arg;
113	>	//get the number of unit cells in each direction:
114
115	<	if (nx <= 0) {
116	<	std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl;
117	<	exit(1);
118	<	}
115	>	nx = args_info.nx_arg;
116
117	<	ny = args_info.ny_arg;
117	>	if (nx <= 0) {
118	>	sprintf(painCave.errMsg, "The number of unit cells in the x direction "
119	>	"must be greater than 0.");
120	>	painCave.isFatal = 1;
121	>	simError();
122	>	}
123
124	<	if (ny <= 0) {
123	<	std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl;
124	<	exit(1);
125	<	}
124	>	ny = args_info.ny_arg;
125
126	<	nz = args_info.nz_arg;
126	>	if (ny <= 0) {
127	>	sprintf(painCave.errMsg, "The number of unit cells in the y direction "
128	>	"must be greater than 0.");
129	>	painCave.isFatal = 1;
130	>	simError();
131	>	}
132
133	<	if (nz <= 0) {
130	<	std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl;
131	<	exit(1);
132	<	}
133	>	nz = args_info.nz_arg;
134
135	<	//get input file name
136	<	if (args_info.inputs_num)
137	<	inputFileName = args_info.inputs[0];
138	<	else {
139	<	std::cerr << "You must specify a input file name.\n" << std::endl;
140	<	cmdline_parser_print_help();
140	<	exit(1);
141	<	}
135	>	if (nz <= 0) {
136	>	sprintf(painCave.errMsg, "The number of unit cells in the z direction "
137	>	"must be greater than 0.");
138	>	painCave.isFatal = 1;
139	>	simError();
140	>	}
141
142	<	//parse md file and set up the system
144	<	SimCreator oldCreator;
145	<	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
142	>	int nSites = nMolPerCell * nx * ny * nz;
143
144	<	if (oldInfo->getNMoleculeStamp()>= 2) {
145	<	std::cerr << "can not build the system with more than two components"
146	<	<< std::endl;
147	<	exit(1);
148	<	}
144	>	//get input file name
145	>	if (args_info.inputs_num)
146	>	inputFileName = args_info.inputs[0];
147	>	else {
148	>	sprintf(painCave.errMsg, "No input .md file name was specified "
149	>	"on the command line");
150	>	painCave.isFatal = 1;
151	>	simError();
152	>	}
153
154	<	//get mass of molecule.
154	>	//parse md file and set up the system
155
156	<	mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
156	>	SimCreator oldCreator;
157	>	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
158	>	Globals* simParams = oldInfo->getSimParams();
159
160	<	//creat lattice
158	<	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
160	>	// Calculate lattice constant (in Angstroms)
161
162	<	if (simpleLat == NULL) {
163	<	std::cerr << "Error in creating lattice" << std::endl;
162	<	exit(1);
163	<	}
162	>	RealType avgMass = getMolMass(oldInfo->getMoleculeStamp(0),
163	>	oldInfo->getForceField());
164
165	<	numMolPerCell = simpleLat->getNumSitesPerCell();
165	>	latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
166	>	(RealType)(1.0 / 3.0));
167	>
168	>	// Set the lattice constant
169	>
170	>	lc.push_back(latticeConstant);
171	>	simpleLat->setLatticeConstant(lc);
172
173	<	//calculate lattice constant (in Angstrom)
168	<	latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
169	<	1.0 / 3.0);
173	>	// Calculate the lattice sites and fill the lattice vector.
174
175	<	//set lattice constant
172	<	lc.push_back(latticeConstant);
173	<	simpleLat->setLatticeConstant(lc);
175	>	// Get the standard orientations of the cell sites
176
177	<	//calculate the total number of molecules
176	<	numMol = nx * ny * nz * numMolPerCell;
177	>	latticeOrt = simpleLat->getLatticePointsOrt();
178
179	<	if (oldInfo->getNGlobalMolecules() != numMol) {
180	<	outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
181	<	outMdFileName = outPrefix + ".md";
179	>	vector<Vector3d> sites;
180	>	vector<Vector3d> orientations;
181	>
182	>	for(int i = 0; i < nx; i++) {
183	>	for(int j = 0; j < ny; j++) {
184	>	for(int k = 0; k < nz; k++) {
185
186	<	//creat new .md file on fly which corrects the number of molecule
187	<	createMdFile(inputFileName, outMdFileName, numMol);
188	<	std::cerr
189	<	<< "SimpleBuilder Error: the number of molecule and the density are not matched"
190	<	<< std::endl;
191	<	std::cerr << "A new .md file: " << outMdFileName
192	<	<< " is generated, use it to rerun the simpleBuilder" << std::endl;
193	<	exit(1);
186	>	// Get the position of the cell sites
187	>
188	>	simpleLat->getLatticePointsPos(latticePos, i, j, k);
189	>
190	>	for(int l = 0; l < nMolPerCell; l++) {
191	>	sites.push_back(latticePos[l]);
192	>	orientations.push_back(latticeOrt[l]);
193	>	}
194	>	}
195		}
196	+	}
197	+
198	+	outputFileName = args_info.output_arg;
199	+
200	+	// create a new .md file on the fly which corrects the number of molecules
201
202	<	//determine the output file names
193	<	if (args_info.output_given)
194	<	outInitFileName = args_info.output_arg;
195	<	else
196	<	outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
197	<
198	<	//creat Molocator
199	<	locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
202	>	createMdFile(inputFileName, outputFileName, nSites);
203
204	<	//fill Hmat
205	<	hmat(0, 0)= nx * latticeConstant;
203	<	hmat(0, 1) = 0.0;
204	<	hmat(0, 2) = 0.0;
204	>	if (oldInfo != NULL)
205	>	delete oldInfo;
206
207	<	hmat(1, 0) = 0.0;
208	<	hmat(1, 1) = ny * latticeConstant;
208	<	hmat(1, 2) = 0.0;
207	>	// We need to read in the new SimInfo object, then Parse the
208	>	// md file and set up the system
209
210	<	hmat(2, 0) = 0.0;
211	<	hmat(2, 1) = 0.0;
212	<	hmat(2, 2) = nz * latticeConstant;
210	>	SimCreator newCreator;
211	>	SimInfo* newInfo = newCreator.createSim(outputFileName, false);
212
213	<	//set Hmat
215	<	oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
213	>	// fill Hmat
214
215	<	//place the molecules
215	>	hmat(0, 0) = nx * latticeConstant;
216	>	hmat(0, 1) = 0.0;
217	>	hmat(0, 2) = 0.0;
218
219	<	curMolIndex = 0;
219	>	hmat(1, 0) = 0.0;
220	>	hmat(1, 1) = ny * latticeConstant;
221	>	hmat(1, 2) = 0.0;
222
223	<	//get the orientation of the cell sites
224	<	//for the same type of molecule in same lattice, it will not change
225	<	latticeOrt = simpleLat->getLatticePointsOrt();
223	>	hmat(2, 0) = 0.0;
224	>	hmat(2, 1) = 0.0;
225	>	hmat(2, 2) = nz * latticeConstant;
226
227	<	Molecule* mol;
226	<	SimInfo::MoleculeIterator mi;
227	<	mol = oldInfo->beginMolecule(mi);
228	<	for(int i = 0; i < nx; i++) {
229	<	for(int j = 0; j < ny; j++) {
230	<	for(int k = 0; k < nz; k++) {
227	>	// Set Hmat
228
229	<	//get the position of the cell sites
233	<	simpleLat->getLatticePointsPos(latticePos, i, j, k);
229	>	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
230
231	<	for(int l = 0; l < numMolPerCell; l++) {
232	<	if (mol != NULL) {
233	<	locator->placeMol(latticePos[l], latticeOrt[l], mol);
234	<	} else {
235	<	std::cerr << std::endl;
236	<	}
237	<	mol = oldInfo->nextMolecule(mi);
238	<	}
239	<	}
240	<	}
241	<	}
231	>	// place the molecules
232	>
233	>	Molecule* mol;
234	>	locator = new MoLocator(newInfo->getMoleculeStamp(0),
235	>	newInfo->getForceField());
236	>	for (int n = 0; n < nSites; n++) {
237	>	mol = newInfo->getMoleculeByGlobalIndex(n);
238	>	locator->placeMol(sites[n], orientations[n], mol);
239	>	}
240	>
241	>	// Create DumpWriter and write out the coordinates
242
243	<	//create dumpwriter and write out the coordinates
244	<	oldInfo->setFinalConfigFileName(outInitFileName);
245	<	writer = new DumpWriter(oldInfo);
243	>	writer = new DumpWriter(newInfo, outputFileName);
244	>
245	>	if (writer == NULL) {
246	>	sprintf(painCave.errMsg, "error in creating DumpWriter");
247	>	painCave.isFatal = 1;
248	>	simError();
249	>	}
250
251	<	if (writer == NULL) {
252	<	std::cerr << "error in creating DumpWriter" << std::endl;
253	<	exit(1);
254	<	}
251	>	writer->writeDump();
252
253	<	writer->writeDump();
257	<	std::cout << "new initial configuration file: " << outInitFileName
258	<	<< " is generated." << std::endl;
253	>	// deleting the writer will put the closing at the end of the dump file.
254
255	<	//delete objects
255	>	delete writer;
256
257	<	//delete oldInfo and oldSimSetup
258	<	if (oldInfo != NULL)
259	<	delete oldInfo;
260	<
261	<	if (writer != NULL)
267	<	delete writer;
268	<
269	<	return 0;
257	>	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
258	>	"generated.\n", outputFileName.c_str());
259	>	painCave.isFatal = 0;
260	>	simError();
261	>	return 0;
262		}
263
264	<	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
265	<	int numMol) {
266	<	ifstream oldMdFile;
267	<	ofstream newMdFile;
268	<	const int MAXLEN = 65535;
269	<	char buffer[MAXLEN];
264	>	void createMdFile(const std::string&oldMdFileName,
265	>	const std::string&newMdFileName,
266	>	int nMol) {
267	>	ifstream oldMdFile;
268	>	ofstream newMdFile;
269	>	const int MAXLEN = 65535;
270	>	char buffer[MAXLEN];
271
272	<	//create new .md file based on old .md file
273	<	oldMdFile.open(oldMdFileName.c_str());
274	<	newMdFile.open(newMdFileName.c_str());
272	>	//create new .md file based on old .md file
273	>	oldMdFile.open(oldMdFileName.c_str());
274	>	newMdFile.open(newMdFileName.c_str());
275
276	<	oldMdFile.getline(buffer, MAXLEN);
276	>	oldMdFile.getline(buffer, MAXLEN);
277
278	<	while (!oldMdFile.eof()) {
278	>	while (!oldMdFile.eof()) {
279
280	<	//correct molecule number
281	<	if (strstr(buffer, "nMol") != NULL) {
282	<	sprintf(buffer, "\t\tnMol = %d;", numMol);
283	<	newMdFile << buffer << std::endl;
284	<	} else
285	<	newMdFile << buffer << std::endl;
280	>	//correct molecule number
281	>	if (strstr(buffer, "nMol") != NULL) {
282	>	sprintf(buffer, "\t\tnMol = %d;", nMol);
283	>	newMdFile << buffer << std::endl;
284	>	} else
285	>	newMdFile << buffer << std::endl;
286
287	<	oldMdFile.getline(buffer, MAXLEN);
288	<	}
287	>	oldMdFile.getline(buffer, MAXLEN);
288	>	}
289
290	<	oldMdFile.close();
291	<	newMdFile.close();
290	>	oldMdFile.close();
291	>	newMdFile.close();
292		}
293

Comparing:
trunk/src/applications/simpleBuilder/simpleBuilder.cpp (property svn:keywords), Revision 376 by tim, Thu Feb 24 20:55:07 2005 UTC vs.
branches/development/src/applications/simpleBuilder/simpleBuilder.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines