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 12 by tim, Tue Sep 28 23:24:25 2004 UTC vs.
Revision 1978 by gezelter, Thu Mar 13 13:03:11 2014 UTC

#	Line 1 | Line 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
5	+	* non-exclusive, royalty free, license to use, modify and
6	+	* redistribute this software in source and binary code form, provided
7	+	* that the following conditions are met:
8	+	*
9	+	* 1. Redistributions of source code must retain the above copyright
10	+	*    notice, this list of conditions and the following disclaimer.
11	+	*
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.
16	+	*
17	+	* This software is provided "AS IS," without a warranty of any
18	+	* kind. All express or implied conditions, representations and
19	+	* warranties, including any implied warranty of merchantability,
20	+	* fitness for a particular purpose or non-infringement, are hereby
21	+	* excluded.  The University of Notre Dame and its licensors shall not
22	+	* be liable for any damages suffered by licensee as a result of
23	+	* using, modifying or distributing the software or its
24	+	* derivatives. In no event will the University of Notre Dame or its
25	+	* licensors be liable for any lost revenue, profit or data, or for
26	+	* direct, indirect, special, consequential, incidental or punitive
27	+	* damages, however caused and regardless of the theory of liability,
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>
44		#include <cstdio>
45		#include <cstring>
#	Line 7 | Line 49
49		#include <map>
50		#include <fstream>
51
10	–	#include "io/Globals.hpp"
11	–	#include "brains/SimInfo.hpp"
12	–	#include "brains/SimSetup.hpp"
52		#include "applications/simpleBuilder/simpleBuilderCmd.h"
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"
59	+	#include "io/DumpWriter.hpp"
60	+	#include "math/Vector3.hpp"
61	+	#include "math/SquareMatrix3.hpp"
62		#include "utils/StringUtils.hpp"
15	–	#include "applications/simpleBuilder/LatticeFactory.hpp"
16	–	#include "applications/simpleBuilder/Vector3d.hpp"
17	–	#include "applications/simpleBuilder/MoLocator.hpp"
18	–	#include "applications/simpleBuilder/Lattice.hpp"
63
64		using namespace std;
65	+	using namespace OpenMD;
66
67	<	void createMdFile(const string& oldMdFileName, const string& newMdFileName, int numMol);
68	<	double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF);
67	>	void createMdFile(const std::string&oldMdFileName,
68	>	const std::string&newMdFileName,
69	>	int nMol);
70
71	<	int main( int argc, char* argv[]){
71	>	int main(int argc, char *argv []) {
72
73	+	registerLattice();
74	+
75		gengetopt_args_info args_info;
76	<	string latticeType;
77	<	string inputFileName;
78	<	string outPrefix;
79	<	string outMdFileName;
80	<	string outInitFileName;
81	<	SimInfo* oldInfo;
82	<	SimSetup* oldSimSetup;
83	<	BaseLattice* simpleLat;
36	<	int numMol;
37	<	double latticeConstant;
38	<	vector<double> lc;
39	<	double mass;
40	<	const double rhoConvertConst = 1.661;
41	<	double density;
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	<	double Hmat[3][3];
85	>	Mat3x3d hmat;
86		MoLocator *locator;
87	<	vector<Vector3d> latticePos;
88	<	vector<Vector3d> latticeOrt;
89	<	int numMolPerCell;
90	<	int curMolIndex;
91	<	DumpWriter* writer;
50	<
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) ;
95	<
93	>	if (cmdline_parser(argc, argv, &args_info) != 0)
94	>	exit(1);
95	>
96		density = args_info.density_arg;
97
98		//get lattice type
99	<	latticeType = UpperCase(args_info.latticetype_arg);
100	<	if(!LatticeFactory::getInstance()->hasLatticeCreator(latticeType)){
101	<	cerr << latticeType << " is an invalid lattice type" << endl;
102	<	cerr << LatticeFactory::getInstance()->toString() << endl;
103	<	exit(1);
99	>	latticeType = "FCC";
100	>
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		}
109	+	nMolPerCell = simpleLat->getNumSitesPerCell();
110
111	<	//get the number of unit cell
111	>	//get the number of unit cells in each direction:
112	>
113		nx = args_info.nx_arg;
114	<	if(nx <= 0){
115	<	cerr << "The number of unit cell in h direction must be greater than 0" << endl;
116	<	exit(1);
114	>
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		ny = args_info.ny_arg;
123	<	if(ny <= 0){
124	<	cerr << "The number of unit cell in l direction must be greater than 0" << endl;
125	<	exit(1);
123	>
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		nz = args_info.nz_arg;
132	<	if(nz <= 0){
133	<	cerr << "The number of unit cell in k direction must be greater than 0" << endl;
134	<	exit(1);
132	>
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	<
139	>
140	>	int nSites = nMolPerCell * nx * ny * nz;
141	>
142		//get input file name
143	<	if (args_info.inputs_num)
143	>	if (args_info.inputs_num)
144		inputFileName = args_info.inputs[0];
145	<	else {
146	<	cerr <<"You must specify a input file name.\n" << endl;
147	<	cmdline_parser_print_help();
148	<	exit(1);
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
93	–
152		//parse md file and set up the system
95	–	oldInfo = new SimInfo;
96	–	if(oldInfo == NULL){
97	–	cerr << "error in creating SimInfo" << endl;
98	–	exit(1);
99	–	}
153
154	<	oldSimSetup = new SimSetup();
155	<	if(oldSimSetup == NULL){
156	<	cerr << "error in creating SimSetup" << endl;
104	<	exit(1);
105	<	}
154	>	SimCreator oldCreator;
155	>	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
156	>	Globals* simParams = oldInfo->getSimParams();
157
158	<	oldSimSetup->suspendInit();
108	<	oldSimSetup->setSimInfo(oldInfo );
109	<	oldSimSetup->parseFile(&inputFileName[0] );
110	<	oldSimSetup->createSim();
111	<
112	<	if(oldInfo->nComponents >=2){
113	<	cerr << "can not build the system with more than two components" << endl;
114	<	exit(1);
115	<	}
116	<
117	<	//get mass of molecule.
118	<	//Due to the design of OOPSE, given atom type, we have to query forcefiled to get the mass
119	<	mass = getMolMass(oldInfo->compStamps[0], oldSimSetup->getForceField());
120	<
121	<	//creat lattice
122	<	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
123	<	if(simpleLat == NULL){
124	<	cerr << "Error in creating lattice" << endl;
125	<	exit(1);
126	<	}
158	>	// Calculate lattice constant (in Angstroms)
159
160	<	numMolPerCell = simpleLat->getNumSitesPerCell();
160	>	RealType avgMass = MoLocator::getMolMass(oldInfo->getMoleculeStamp(0),
161	>	oldInfo->getForceField());
162	>
163	>	latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
164	>	(RealType)(1.0 / 3.0));
165
166	<	//calculate lattice constant (in Angstrom)
131	<	latticeConstant = pow(rhoConvertConst * numMolPerCell * mass /density, 1.0/3.0);
166	>	// Set the lattice constant
167
133	–	//set lattice constant
168		lc.push_back(latticeConstant);
169		simpleLat->setLatticeConstant(lc);
136	–
137	–	//calculate the total number of molecules
138	–	numMol = nx * ny * nz * numMolPerCell;
170
171	<	if (oldInfo->n_mol != numMol){
171	>	// Calculate the lattice sites and fill the lattice vector.
172
173	<	outPrefix = getPrefix(inputFileName.c_str()) + "_" + latticeType;
143	<	outMdFileName = outPrefix + ".md";
173	>	// Get the standard orientations of the cell sites
174
175	<	//creat new .md file on fly which corrects the number of molecule
176	<	createMdFile(inputFileName, outMdFileName, numMol);
177	<	cerr << "SimpleBuilder Error: the number of molecule and the density are not matched" <<endl;
178	<	cerr << "A new .md file: " << outMdFileName << " is generated, use it to rerun the simpleBuilder" << endl;
179	<	exit(1);
175	>	latticeOrt = simpleLat->getLatticePointsOrt();
176	>
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	>	// 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	<	//determine the output file names
153	<	if (args_info.output_given)
154	<	outInitFileName = args_info.output_arg;
155	<	else
156	<	outInitFileName = getPrefix(inputFileName.c_str())  + ".in";
196	>	outputFileName = args_info.output_arg;
197
198	<
159	<	//allocat memory for storing pos, vel and etc
160	<	oldInfo->getConfiguration()->createArrays(oldInfo->n_atoms);
161	<	for (int i = 0; i < oldInfo->n_atoms; i++)
162	<	oldInfo->atoms[i]->setCoords();
198	>	// create a new .md file on the fly which corrects the number of molecules
199
200	<	//creat Molocator
165	<	locator = new MoLocator(oldInfo->compStamps[0], oldSimSetup->getForceField());
200	>	createMdFile(inputFileName, outputFileName, nSites);
201
202	<	//fill Hmat
168	<	Hmat[0][0] = nx * latticeConstant;
169	<	Hmat[0][1] = 0.0;
170	<	Hmat[0][2] = 0.0;
202	>	delete oldInfo;
203
204	<	Hmat[1][0] = 0.0;
205	<	Hmat[1][1] = ny * latticeConstant;
174	<	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;
178	<	Hmat[2][2] = nz * latticeConstant ;
207	>	SimCreator newCreator;
208	>	SimInfo* newInfo = newCreator.createSim(outputFileName, false);
209
210	<	//set Hmat
181	<	oldInfo->setBoxM(Hmat);
182	<
183	<	//place the molecules
210	>	// fill Hmat
211
212	<	curMolIndex = 0;
212	>	hmat(0, 0) = nx * latticeConstant;
213	>	hmat(0, 1) = 0.0;
214	>	hmat(0, 2) = 0.0;
215
216	<	//get the orientation of the cell sites
217	<	//for the same type of molecule in same lattice, it will not change
218	<	latticeOrt = simpleLat->getLatticePointsOrt();
216	>	hmat(1, 0) = 0.0;
217	>	hmat(1, 1) = ny * latticeConstant;
218	>	hmat(1, 2) = 0.0;
219
220	<	for(int i =0; i < nx; i++){
221	<	for(int j=0; j < ny; j++){
222	<	for(int k = 0; k < nz; k++){
220	>	hmat(2, 0) = 0.0;
221	>	hmat(2, 1) = 0.0;
222	>	hmat(2, 2) = nz * latticeConstant;
223
224	<	//get the position of the cell sites
196	<	simpleLat->getLatticePointsPos(latticePos, i, j, k);
224	>	// Set Hmat
225
226	<	for(int l = 0; l < numMolPerCell; l++)
199	<	locator->placeMol(latticePos[l], latticeOrt[l], &(oldInfo->molecules[curMolIndex++]));
200	<	}
201	<	}
202	<	}
226	>	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
227
228	<	//create dumpwriter and write out the coordinates
229	<	oldInfo->finalName = outInitFileName;
230	<	writer = new DumpWriter( oldInfo );
231	<	if(writer == NULL){
232	<	cerr << "error in creating DumpWriter" << endl;
233	<	exit(1);
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	<	writer->writeFinal(0);
238	<	cout << "new initial configuration file: " << outInitFileName <<" is generated." << endl;
213	<	//delete objects
237	>
238	>	// Create DumpWriter and write out the coordinates
239
240	<	//delete oldInfo and oldSimSetup
216	<	if(oldInfo != NULL)
217	<	delete oldInfo;
240	>	writer = new DumpWriter(newInfo, outputFileName);
241
242	<	if(oldSimSetup != NULL)
243	<	delete oldSimSetup;
244	<
245	<	if (writer != NULL)
246	<	delete writer;
242	>	if (writer == NULL) {
243	>	sprintf(painCave.errMsg, "error in creating DumpWriter");
244	>	painCave.isFatal = 1;
245	>	simError();
246	>	}
247	>
248	>	writer->writeDump();
249	>
250	>	// deleting the writer will put the closing at the end of the dump file.
251	>
252	>	delete writer;
253	>
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 string& oldMdFileName, const string& newMdFileName, int numMol){
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;
#	Line 235 | Line 272 | void createMdFile(const string& oldMdFileName, const s
272		newMdFile.open(newMdFileName.c_str());
273
274		oldMdFile.getline(buffer, MAXLEN);
238	–	while(!oldMdFile.eof()){
275
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 << endl;
282	<	}
283	<	else
246	<	newMdFile << buffer << endl;
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		}
287
288		oldMdFile.close();
289		newMdFile.close();
253	–
290		}
291
256	–	double getMolMass(MoleculeStamp* molStamp, ForceFields* myFF){
257	–	int nAtoms;
258	–	AtomStamp* currAtomStamp;
259	–	double totMass;
260	–
261	–	totMass = 0;
262	–	nAtoms = molStamp->getNAtoms();
263	–
264	–	for(size_t i=0; i<nAtoms; i++){
265	–	currAtomStamp = molStamp->getAtom(i);
266	–	totMass += myFF->getAtomTypeMass(currAtomStamp->getType());
267	–	}
268	–
269	–	return totMass;
270	–	}

Comparing trunk/src/applications/simpleBuilder/simpleBuilder.cpp (property svn:keywords):
Revision 12 by tim, Tue Sep 28 23:24:25 2004 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