ViewVC Help

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

root/OpenMD/trunk/src/applications/simpleBuilder/simpleBuilder.cpp

Comparing trunk/src/applications/simpleBuilder/simpleBuilder.cpp (file contents):
Revision 488 by tim, Tue Apr 12 22:42:13 2005 UTC vs.
Revision 1978 by gezelter, Thu Mar 13 13:03:11 2014 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, 234107 (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 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
71	<	registerForceFields();
72	<	registerLattice();
73	>	registerLattice();
74
75	<	gengetopt_args_info args_info;
76	<	std::string latticeType;
77	<	std::string inputFileName;
78	<	std::string outPrefix;
79	<	std::string outMdFileName;
80	<	std::string outInitFileName;
81	<	Lattice *simpleLat;
82	<	int numMol;
83	<	double latticeConstant;
84	<	std::vector<double> lc;
85	<	double mass;
86	<	const double rhoConvertConst = 1.661;
87	<	double density;
88	<	int nx,
89	<	ny,
90	<	nz;
90	<	Mat3x3d hmat;
91	<	MoLocator *locator;
92	<	std::vector<Vector3d> latticePos;
93	<	std::vector<Vector3d> latticeOrt;
94	<	int numMolPerCell;
95	<	int curMolIndex;
96	<	DumpWriter *writer;
75	>	gengetopt_args_info args_info;
76	>	std::string latticeType;
77	>	std::string inputFileName;
78	>	std::string outputFileName;
79	>	Lattice *simpleLat;
80	>	RealType latticeConstant;
81	>	std::vector<RealType> lc;
82	>	const RealType rhoConvertConst = 1.661;
83	>	RealType density;
84	>	int nx, ny, nz;
85	>	Mat3x3d hmat;
86	>	MoLocator *locator;
87	>	std::vector<Vector3d> latticePos;
88	>	std::vector<Vector3d> latticeOrt;
89	>	int nMolPerCell;
90	>	DumpWriter *writer;
91
92	<	// parse command line arguments
93	<	if (cmdline_parser(argc, argv, &args_info) != 0)
94	<	exit(1);
92	>	// parse command line arguments
93	>	if (cmdline_parser(argc, argv, &args_info) != 0)
94	>	exit(1);
95
96	<	density = args_info.density_arg;
96	>	density = args_info.density_arg;
97
98	<	//get lattice type
99	<	latticeType = UpperCase(args_info.latticetype_arg);
98	>	//get lattice type
99	>	latticeType = "FCC";
100
101	<	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
101	>	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
102
103	<	if (simpleLat == NULL) {
104	<	sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n",
105	<	latticeType.c_str());
106	<	painCave.isFatal = 1;
107	<	simError();
108	<	}
103	>	if (simpleLat == NULL) {
104	>	sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n",
105	>	latticeType.c_str());
106	>	painCave.isFatal = 1;
107	>	simError();
108	>	}
109	>	nMolPerCell = simpleLat->getNumSitesPerCell();
110
111	<	//get the number of unit cell
117	<	nx = args_info.nx_arg;
111	>	//get the number of unit cells in each direction:
112
113	<	if (nx <= 0) {
120	<	std::cerr << "The number of unit cell in h direction must be greater than 0" << std::endl;
121	<	exit(1);
122	<	}
113	>	nx = args_info.nx_arg;
114
115	<	ny = args_info.ny_arg;
115	>	if (nx <= 0) {
116	>	sprintf(painCave.errMsg, "The number of unit cells in the x direction "
117	>	"must be greater than 0.");
118	>	painCave.isFatal = 1;
119	>	simError();
120	>	}
121
122	<	if (ny <= 0) {
127	<	std::cerr << "The number of unit cell in l direction must be greater than 0" << std::endl;
128	<	exit(1);
129	<	}
122	>	ny = args_info.ny_arg;
123
124	<	nz = args_info.nz_arg;
124	>	if (ny <= 0) {
125	>	sprintf(painCave.errMsg, "The number of unit cells in the y direction "
126	>	"must be greater than 0.");
127	>	painCave.isFatal = 1;
128	>	simError();
129	>	}
130
131	<	if (nz <= 0) {
134	<	std::cerr << "The number of unit cell in k direction must be greater than 0" << std::endl;
135	<	exit(1);
136	<	}
131	>	nz = args_info.nz_arg;
132
133	<	//get input file name
134	<	if (args_info.inputs_num)
135	<	inputFileName = args_info.inputs[0];
136	<	else {
137	<	std::cerr << "You must specify a input file name.\n" << std::endl;
138	<	cmdline_parser_print_help();
144	<	exit(1);
145	<	}
133	>	if (nz <= 0) {
134	>	sprintf(painCave.errMsg, "The number of unit cells in the z direction "
135	>	"must be greater than 0.");
136	>	painCave.isFatal = 1;
137	>	simError();
138	>	}
139
140	<	//parse md file and set up the system
148	<	SimCreator oldCreator;
149	<	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
140	>	int nSites = nMolPerCell * nx * ny * nz;
141
142	<	if (oldInfo->getNMoleculeStamp()>= 2) {
143	<	std::cerr << "can not build the system with more than two components"
144	<	<< std::endl;
145	<	exit(1);
146	<	}
142	>	//get input file name
143	>	if (args_info.inputs_num)
144	>	inputFileName = args_info.inputs[0];
145	>	else {
146	>	sprintf(painCave.errMsg, "No input .md file name was specified "
147	>	"on the command line");
148	>	painCave.isFatal = 1;
149	>	simError();
150	>	}
151
152	<	//get mass of molecule.
152	>	//parse md file and set up the system
153
154	<	mass = getMolMass(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
154	>	SimCreator oldCreator;
155	>	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
156	>	Globals* simParams = oldInfo->getSimParams();
157
158	<	//creat lattice
162	<	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
158	>	// Calculate lattice constant (in Angstroms)
159
160	<	if (simpleLat == NULL) {
161	<	std::cerr << "Error in creating lattice" << std::endl;
166	<	exit(1);
167	<	}
160	>	RealType avgMass = MoLocator::getMolMass(oldInfo->getMoleculeStamp(0),
161	>	oldInfo->getForceField());
162
163	<	numMolPerCell = simpleLat->getNumSitesPerCell();
163	>	latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
164	>	(RealType)(1.0 / 3.0));
165	>
166	>	// Set the lattice constant
167	>
168	>	lc.push_back(latticeConstant);
169	>	simpleLat->setLatticeConstant(lc);
170
171	<	//calculate lattice constant (in Angstrom)
172	<	latticeConstant = pow(rhoConvertConst * numMolPerCell * mass / density,
173	<	1.0 / 3.0);
171	>	// Calculate the lattice sites and fill the lattice vector.
172
173	<	//set lattice constant
176	<	lc.push_back(latticeConstant);
177	<	simpleLat->setLatticeConstant(lc);
173	>	// Get the standard orientations of the cell sites
174
175	<	//calculate the total number of molecules
180	<	numMol = nx * ny * nz * numMolPerCell;
175	>	latticeOrt = simpleLat->getLatticePointsOrt();
176
177	<	if (oldInfo->getNGlobalMolecules() != numMol) {
178	<	outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
179	<	outMdFileName = outPrefix + ".md";
177	>	vector<Vector3d> sites;
178	>	vector<Vector3d> orientations;
179	>
180	>	for(int i = 0; i < nx; i++) {
181	>	for(int j = 0; j < ny; j++) {
182	>	for(int k = 0; k < nz; k++) {
183
184	<	//creat new .md file on fly which corrects the number of molecule
185	<	createMdFile(inputFileName, outMdFileName, numMol);
186	<	std::cerr
187	<	<< "SimpleBuilder Error: the number of molecule and the density are not matched"
188	<	<< std::endl;
189	<	std::cerr << "A new .md file: " << outMdFileName
190	<	<< " is generated, use it to rerun the simpleBuilder" << std::endl;
191	<	exit(1);
184	>	// Get the position of the cell sites
185	>
186	>	simpleLat->getLatticePointsPos(latticePos, i, j, k);
187	>
188	>	for(int l = 0; l < nMolPerCell; l++) {
189	>	sites.push_back(latticePos[l]);
190	>	orientations.push_back(latticeOrt[l]);
191	>	}
192	>	}
193		}
194	+	}
195	+
196	+	outputFileName = args_info.output_arg;
197	+
198	+	// create a new .md file on the fly which corrects the number of molecules
199
200	<	//determine the output file names
197	<	if (args_info.output_given)
198	<	outInitFileName = args_info.output_arg;
199	<	else
200	<	outInitFileName = getPrefix(inputFileName.c_str()) + ".in";
201	<
202	<	//creat Molocator
203	<	locator = new MoLocator(oldInfo->getMoleculeStamp(0), oldInfo->getForceField());
200	>	createMdFile(inputFileName, outputFileName, nSites);
201
202	<	//fill Hmat
206	<	hmat(0, 0)= nx * latticeConstant;
207	<	hmat(0, 1) = 0.0;
208	<	hmat(0, 2) = 0.0;
202	>	delete oldInfo;
203
204	<	hmat(1, 0) = 0.0;
205	<	hmat(1, 1) = ny * latticeConstant;
212	<	hmat(1, 2) = 0.0;
204	>	// We need to read in the new SimInfo object, then Parse the
205	>	// md file and set up the system
206
207	<	hmat(2, 0) = 0.0;
208	<	hmat(2, 1) = 0.0;
216	<	hmat(2, 2) = nz * latticeConstant;
207	>	SimCreator newCreator;
208	>	SimInfo* newInfo = newCreator.createSim(outputFileName, false);
209
210	<	//set Hmat
219	<	oldInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
210	>	// fill Hmat
211
212	<	//place the molecules
212	>	hmat(0, 0) = nx * latticeConstant;
213	>	hmat(0, 1) = 0.0;
214	>	hmat(0, 2) = 0.0;
215
216	<	curMolIndex = 0;
216	>	hmat(1, 0) = 0.0;
217	>	hmat(1, 1) = ny * latticeConstant;
218	>	hmat(1, 2) = 0.0;
219
220	<	//get the orientation of the cell sites
221	<	//for the same type of molecule in same lattice, it will not change
222	<	latticeOrt = simpleLat->getLatticePointsOrt();
220	>	hmat(2, 0) = 0.0;
221	>	hmat(2, 1) = 0.0;
222	>	hmat(2, 2) = nz * latticeConstant;
223
224	<	Molecule* mol;
230	<	SimInfo::MoleculeIterator mi;
231	<	mol = oldInfo->beginMolecule(mi);
232	<	for(int i = 0; i < nx; i++) {
233	<	for(int j = 0; j < ny; j++) {
234	<	for(int k = 0; k < nz; k++) {
224	>	// Set Hmat
225
226	<	//get the position of the cell sites
237	<	simpleLat->getLatticePointsPos(latticePos, i, j, k);
226	>	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
227
228	<	for(int l = 0; l < numMolPerCell; l++) {
229	<	if (mol != NULL) {
230	<	locator->placeMol(latticePos[l], latticeOrt[l], mol);
231	<	} else {
232	<	std::cerr << std::endl;
233	<	}
234	<	mol = oldInfo->nextMolecule(mi);
235	<	}
236	<	}
237	<	}
238	<	}
228	>	// place the molecules
229	>
230	>	Molecule* mol;
231	>	locator = new MoLocator(newInfo->getMoleculeStamp(0),
232	>	newInfo->getForceField());
233	>	for (int n = 0; n < nSites; n++) {
234	>	mol = newInfo->getMoleculeByGlobalIndex(n);
235	>	locator->placeMol(sites[n], orientations[n], mol);
236	>	}
237	>
238	>	// Create DumpWriter and write out the coordinates
239
240	<	//create dumpwriter and write out the coordinates
241	<	oldInfo->setFinalConfigFileName(outInitFileName);
242	<	writer = new DumpWriter(oldInfo);
240	>	writer = new DumpWriter(newInfo, outputFileName);
241	>
242	>	if (writer == NULL) {
243	>	sprintf(painCave.errMsg, "error in creating DumpWriter");
244	>	painCave.isFatal = 1;
245	>	simError();
246	>	}
247
248	<	if (writer == NULL) {
256	<	std::cerr << "error in creating DumpWriter" << std::endl;
257	<	exit(1);
258	<	}
248	>	writer->writeDump();
249
250	<	writer->writeDump();
261	<	std::cout << "new initial configuration file: " << outInitFileName
262	<	<< " is generated." << std::endl;
250	>	// deleting the writer will put the closing at the end of the dump file.
251
252	<	//delete objects
252	>	delete writer;
253
254	<	//delete oldInfo and oldSimSetup
255	<	if (oldInfo != NULL)
256	<	delete oldInfo;
257	<
258	<	if (writer != NULL)
259	<	delete writer;
272	<
273	<	delete simpleLat;
274	<
275	<	return 0;
254	>	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
255	>	"generated.\n", outputFileName.c_str());
256	>	painCave.isFatal = 0;
257	>	painCave.severity = OPENMD_INFO;
258	>	simError();
259	>	return 0;
260		}
261
262	<	void createMdFile(const std::string&oldMdFileName, const std::string&newMdFileName,
263	<	int numMol) {
264	<	ifstream oldMdFile;
265	<	ofstream newMdFile;
266	<	const int MAXLEN = 65535;
267	<	char buffer[MAXLEN];
262	>	void createMdFile(const std::string&oldMdFileName,
263	>	const std::string&newMdFileName,
264	>	int nMol) {
265	>	ifstream oldMdFile;
266	>	ofstream newMdFile;
267	>	const int MAXLEN = 65535;
268	>	char buffer[MAXLEN];
269
270	<	//create new .md file based on old .md file
271	<	oldMdFile.open(oldMdFileName.c_str());
272	<	newMdFile.open(newMdFileName.c_str());
270	>	//create new .md file based on old .md file
271	>	oldMdFile.open(oldMdFileName.c_str());
272	>	newMdFile.open(newMdFileName.c_str());
273
274	<	oldMdFile.getline(buffer, MAXLEN);
274	>	oldMdFile.getline(buffer, MAXLEN);
275
276	<	while (!oldMdFile.eof()) {
276	>	while (!oldMdFile.eof()) {
277
278	<	//correct molecule number
279	<	if (strstr(buffer, "nMol") != NULL) {
280	<	sprintf(buffer, "\t\tnMol = %d;", numMol);
281	<	newMdFile << buffer << std::endl;
282	<	} else
283	<	newMdFile << buffer << std::endl;
278	>	//correct molecule number
279	>	if (strstr(buffer, "nMol") != NULL) {
280	>	sprintf(buffer, "\t\tnMol = %d;", nMol);
281	>	newMdFile << buffer << std::endl;
282	>	} else
283	>	newMdFile << buffer << std::endl;
284
285	<	oldMdFile.getline(buffer, MAXLEN);
286	<	}
285	>	oldMdFile.getline(buffer, MAXLEN);
286	>	}
287
288	<	oldMdFile.close();
289	<	newMdFile.close();
288	>	oldMdFile.close();
289	>	newMdFile.close();
290		}
291

Comparing trunk/src/applications/simpleBuilder/simpleBuilder.cpp (property svn:keywords):
Revision 488 by tim, Tue Apr 12 22:42:13 2005 UTC vs.
Revision 1978 by gezelter, Thu Mar 13 13:03:11 2014 UTC

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

Diff Legend

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