applications/simpleBuilder/simpleBuilder.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <string>
#include <map>
#include <fstream>

#include "applications/simpleBuilder/simpleBuilderCmd.h"
#include "lattice/LatticeFactory.hpp"
#include "utils/MoLocator.hpp"
#include "lattice/Lattice.hpp"
#include "brains/Register.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"

using namespace std;
using namespace OpenMD;

void createMdFile(const std::string&oldMdFileName, 
                  const std::string&newMdFileName,
                  int nMol);

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

  registerLattice();
    
  gengetopt_args_info args_info;
  std::string latticeType;
  std::string inputFileName;
  std::string outputFileName;
  Lattice *simpleLat;
  RealType latticeConstant;
  std::vector<RealType> lc;
  const RealType rhoConvertConst = 1.661;
  RealType density;
  int nx, ny, nz;
  Mat3x3d hmat;
  MoLocator *locator;
  std::vector<Vector3d> latticePos;
  std::vector<Vector3d> latticeOrt;
  int nMolPerCell;
  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 = "FCC";

  simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
    
  if (simpleLat == NULL) {
    sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n",
            latticeType.c_str());
    painCave.isFatal = 1;
    simError();
  }
  nMolPerCell = simpleLat->getNumSitesPerCell();

  //get the number of unit cells in each direction:

  nx = args_info.nx_arg;

  if (nx <= 0) {
    sprintf(painCave.errMsg, "The number of unit cells in the x direction "
            "must be greater than 0.");
    painCave.isFatal = 1;
    simError();
  }

  ny = args_info.ny_arg;

  if (ny <= 0) {
    sprintf(painCave.errMsg, "The number of unit cells in the y direction "
            "must be greater than 0.");
    painCave.isFatal = 1;
    simError();
  }

  nz = args_info.nz_arg;

  if (nz <= 0) {
    sprintf(painCave.errMsg, "The number of unit cells in the z direction "
            "must be greater than 0.");
    painCave.isFatal = 1;
    simError();
  }

  int nSites = nMolPerCell * nx * ny * nz;

  //get input file name
  if (args_info.inputs_num)
    inputFileName = args_info.inputs[0];
  else {
    sprintf(painCave.errMsg, "No input .md file name was specified "
            "on the command line");
    painCave.isFatal = 1;
    simError();
  }

  //parse md file and set up the system

  SimCreator oldCreator;
  SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
  Globals* simParams = oldInfo->getSimParams();

  // Calculate lattice constant (in Angstroms)

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

  latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
                        (RealType)(1.0 / 3.0));
  
  // Set the lattice constant
  
  lc.push_back(latticeConstant);
  simpleLat->setLatticeConstant(lc);

  // Calculate the lattice sites and fill the lattice vector.

  // Get the standard orientations of the cell sites

  latticeOrt = simpleLat->getLatticePointsOrt();

  vector<Vector3d> sites;
  vector<Vector3d> orientations;
  
  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 < nMolPerCell; l++) {
          sites.push_back(latticePos[l]);
          orientations.push_back(latticeOrt[l]);
        }
      }
    }
  }
  
  outputFileName = args_info.output_arg;
  
  // create a new .md file on the fly which corrects the number of molecules

  createMdFile(inputFileName, outputFileName, nSites);

  if (oldInfo != NULL)
    delete oldInfo;

  // We need to read in the new SimInfo object, then Parse the 
  // md file and set up the system

  SimCreator newCreator;
  SimInfo* newInfo = newCreator.createSim(outputFileName, false);

  // 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

  newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);

  // place the molecules
  
  Molecule* mol;
  locator = new MoLocator(newInfo->getMoleculeStamp(0), 
                          newInfo->getForceField());
  for (int n = 0; n < nSites; n++) {
    mol = newInfo->getMoleculeByGlobalIndex(n);
    locator->placeMol(sites[n], orientations[n], mol);
  }
   
  // Create DumpWriter and write out the coordinates

  writer = new DumpWriter(newInfo, outputFileName);
  
  if (writer == NULL) {
    sprintf(painCave.errMsg, "error in creating DumpWriter");
    painCave.isFatal = 1;
    simError();
  }

  writer->writeDump();

  // deleting the writer will put the closing at the end of the dump file.

  delete writer;

  sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
          "generated.\n", outputFileName.c_str());
  painCave.isFatal = 0;
  simError();
  return 0;
}

void createMdFile(const std::string&oldMdFileName, 
                  const std::string&newMdFileName,
                  int nMol) {
  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;", nMol);
      newMdFile << buffer << std::endl;
    } else
      newMdFile << buffer << std::endl;

    oldMdFile.getline(buffer, MAXLEN);
  }

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

Revision:	1725
Committed:	Sat May 26 18:13:43 2012 UTC (12 years, 11 months ago) by gezelter
File size:	8295 byte(s)
Log Message:	Individual ForceField classes have been removed (they were essentially all duplicates anyway). ForceField has moved to brains, and since only one force field is in play at any time, the ForceFieldFactory and Register methods have been removed.
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* 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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* 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	gezelter	1390	*
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	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43	tim	12	#include <cstdlib>
44			#include <cstdio>
45			#include <cstring>
46			#include <cmath>
47			#include <iostream>
48			#include <string>
49			#include <map>
50			#include <fstream>
51
52			#include "applications/simpleBuilder/simpleBuilderCmd.h"
53	gezelter	481	#include "lattice/LatticeFactory.hpp"
54			#include "utils/MoLocator.hpp"
55			#include "lattice/Lattice.hpp"
56	gezelter	246	#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"
63	tim	12
64			using namespace std;
65	gezelter	1390	using namespace OpenMD;
66	tim	12
67	gezelter	1067	void createMdFile(const std::string&oldMdFileName,
68			const std::string&newMdFileName,
69			int nMol);
70
71	gezelter	246	int main(int argc, char *argv []) {
72	tim	12
73	gezelter	507	registerLattice();
74	gezelter	246
75	gezelter	507	gengetopt_args_info args_info;
76			std::string latticeType;
77			std::string inputFileName;
78	gezelter	1067	std::string outputFileName;
79	gezelter	507	Lattice *simpleLat;
80	gezelter	1067	RealType latticeConstant;
81			std::vector<RealType> lc;
82			const RealType rhoConvertConst = 1.661;
83			RealType density;
84			int nx, ny, nz;
85	gezelter	507	Mat3x3d hmat;
86			MoLocator *locator;
87			std::vector<Vector3d> latticePos;
88			std::vector<Vector3d> latticeOrt;
89	gezelter	1067	int nMolPerCell;
90	gezelter	507	DumpWriter *writer;
91	tim	12
92	gezelter	507	// parse command line arguments
93			if (cmdline_parser(argc, argv, &args_info) != 0)
94			exit(1);
95	tim	12
96	gezelter	507	density = args_info.density_arg;
97	tim	12
98	gezelter	507	//get lattice type
99	gezelter	1067	latticeType = "FCC";
100	tim	12
101	gezelter	507	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
102	tim	487
103	gezelter	507	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	gezelter	1067	nMolPerCell = simpleLat->getNumSitesPerCell();
110	tim	12
111	gezelter	1067	//get the number of unit cells in each direction:
112
113	gezelter	507	nx = args_info.nx_arg;
114	tim	12
115	gezelter	507	if (nx <= 0) {
116	gezelter	1067	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	gezelter	507	}
121	tim	12
122	gezelter	507	ny = args_info.ny_arg;
123	tim	12
124	gezelter	507	if (ny <= 0) {
125	gezelter	1067	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	gezelter	507	}
130	tim	12
131	gezelter	507	nz = args_info.nz_arg;
132	tim	12
133	gezelter	507	if (nz <= 0) {
134	gezelter	1067	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	gezelter	507	}
139	tim	12
140	gezelter	1067	int nSites = nMolPerCell * nx * ny * nz;
141
142	gezelter	507	//get input file name
143			if (args_info.inputs_num)
144			inputFileName = args_info.inputs[0];
145			else {
146	gezelter	1067	sprintf(painCave.errMsg, "No input .md file name was specified "
147			"on the command line");
148			painCave.isFatal = 1;
149			simError();
150	gezelter	507	}
151	tim	12
152	gezelter	507	//parse md file and set up the system
153	gezelter	1067
154	gezelter	507	SimCreator oldCreator;
155			SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
156	gezelter	1067	Globals* simParams = oldInfo->getSimParams();
157	tim	12
158	gezelter	1067	// Calculate lattice constant (in Angstroms)
159	tim	12
160	gezelter	1067	RealType avgMass = getMolMass(oldInfo->getMoleculeStamp(0),
161			oldInfo->getForceField());
162	tim	12
163	gezelter	1067	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	tim	12
171	gezelter	1067	// Calculate the lattice sites and fill the lattice vector.
172	tim	12
173	gezelter	1067	// Get the standard orientations of the cell sites
174	tim	12
175	gezelter	1067	latticeOrt = simpleLat->getLatticePointsOrt();
176	tim	12
177	gezelter	1067	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	tim	12
184	gezelter	1067	// 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	tim	12
200	gezelter	1067	createMdFile(inputFileName, outputFileName, nSites);
201	tim	12
202	gezelter	1067	if (oldInfo != NULL)
203			delete oldInfo;
204	gezelter	246
205	gezelter	1067	// We need to read in the new SimInfo object, then Parse the
206			// md file and set up the system
207	tim	12
208	gezelter	1067	SimCreator newCreator;
209			SimInfo* newInfo = newCreator.createSim(outputFileName, false);
210	tim	12
211	gezelter	1067	// fill Hmat
212
213			hmat(0, 0) = nx * latticeConstant;
214	gezelter	507	hmat(0, 1) = 0.0;
215			hmat(0, 2) = 0.0;
216	tim	12
217	gezelter	507	hmat(1, 0) = 0.0;
218			hmat(1, 1) = ny * latticeConstant;
219			hmat(1, 2) = 0.0;
220	tim	12
221	gezelter	507	hmat(2, 0) = 0.0;
222			hmat(2, 1) = 0.0;
223			hmat(2, 2) = nz * latticeConstant;
224	tim	12
225	gezelter	1067	// Set Hmat
226	tim	12
227	gezelter	1067	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
228	gezelter	246
229	gezelter	1067	// place the molecules
230
231	gezelter	507	Molecule* mol;
232	gezelter	1067	locator = new MoLocator(newInfo->getMoleculeStamp(0),
233			newInfo->getForceField());
234			for (int n = 0; n < nSites; n++) {
235			mol = newInfo->getMoleculeByGlobalIndex(n);
236			locator->placeMol(sites[n], orientations[n], mol);
237	gezelter	507	}
238	gezelter	1067
239			// Create DumpWriter and write out the coordinates
240	tim	12
241	gezelter	1067	writer = new DumpWriter(newInfo, outputFileName);
242
243	gezelter	507	if (writer == NULL) {
244	gezelter	1067	sprintf(painCave.errMsg, "error in creating DumpWriter");
245			painCave.isFatal = 1;
246			simError();
247	gezelter	507	}
248	tim	12
249	gezelter	507	writer->writeDump();
250	gezelter	246
251	gezelter	1067	// deleting the writer will put the closing at the end of the dump file.
252	gezelter	246
253	gezelter	1067	delete writer;
254	gezelter	246
255	gezelter	1390	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
256	gezelter	1067	"generated.\n", outputFileName.c_str());
257			painCave.isFatal = 0;
258			simError();
259	gezelter	507	return 0;
260	tim	12	}
261
262	gezelter	1067	void createMdFile(const std::string&oldMdFileName,
263			const std::string&newMdFileName,
264			int nMol) {
265	gezelter	507	ifstream oldMdFile;
266			ofstream newMdFile;
267			const int MAXLEN = 65535;
268			char buffer[MAXLEN];
269	tim	12
270	gezelter	507	//create new .md file based on old .md file
271			oldMdFile.open(oldMdFileName.c_str());
272			newMdFile.open(newMdFileName.c_str());
273	tim	12
274	gezelter	507	oldMdFile.getline(buffer, MAXLEN);
275	gezelter	246
276	gezelter	507	while (!oldMdFile.eof()) {
277	gezelter	246
278	gezelter	507	//correct molecule number
279			if (strstr(buffer, "nMol") != NULL) {
280	gezelter	1067	sprintf(buffer, "\t\tnMol = %d;", nMol);
281	gezelter	507	newMdFile << buffer << std::endl;
282			} else
283			newMdFile << buffer << std::endl;
284	gezelter	246
285	gezelter	507	oldMdFile.getline(buffer, MAXLEN);
286			}
287	gezelter	246
288	gezelter	507	oldMdFile.close();
289			newMdFile.close();
290	tim	12	}
291	gezelter	246