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 (13 years, 6 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.
#	Content
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, 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>
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	#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"
63
64	using namespace std;
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	registerLattice();
74
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);
95
96	density = args_info.density_arg;
97
98	//get lattice type
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 cells in each direction:
112
113	nx = args_info.nx_arg;
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
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
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	int nSites = nMolPerCell * nx * ny * nz;
141
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	//parse md file and set up the system
153
154	SimCreator oldCreator;
155	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
156	Globals* simParams = oldInfo->getSimParams();
157
158	// Calculate lattice constant (in Angstroms)
159
160	RealType avgMass = getMolMass(oldInfo->getMoleculeStamp(0),
161	oldInfo->getForceField());
162
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 the lattice sites and fill the lattice vector.
172
173	// Get the standard orientations of the cell sites
174
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	outputFileName = args_info.output_arg;
197
198	// create a new .md file on the fly which corrects the number of molecules
199
200	createMdFile(inputFileName, outputFileName, nSites);
201
202	if (oldInfo != NULL)
203	delete oldInfo;
204
205	// We need to read in the new SimInfo object, then Parse the
206	// md file and set up the system
207
208	SimCreator newCreator;
209	SimInfo* newInfo = newCreator.createSim(outputFileName, false);
210
211	// fill Hmat
212
213	hmat(0, 0) = nx * latticeConstant;
214	hmat(0, 1) = 0.0;
215	hmat(0, 2) = 0.0;
216
217	hmat(1, 0) = 0.0;
218	hmat(1, 1) = ny * latticeConstant;
219	hmat(1, 2) = 0.0;
220
221	hmat(2, 0) = 0.0;
222	hmat(2, 1) = 0.0;
223	hmat(2, 2) = nz * latticeConstant;
224
225	// Set Hmat
226
227	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
228
229	// place the molecules
230
231	Molecule* mol;
232	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	}
238
239	// Create DumpWriter and write out the coordinates
240
241	writer = new DumpWriter(newInfo, outputFileName);
242
243	if (writer == NULL) {
244	sprintf(painCave.errMsg, "error in creating DumpWriter");
245	painCave.isFatal = 1;
246	simError();
247	}
248
249	writer->writeDump();
250
251	// deleting the writer will put the closing at the end of the dump file.
252
253	delete writer;
254
255	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
256	"generated.\n", outputFileName.c_str());
257	painCave.isFatal = 0;
258	simError();
259	return 0;
260	}
261
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());
273
274	oldMdFile.getline(buffer, MAXLEN);
275
276	while (!oldMdFile.eof()) {
277
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	}
287
288	oldMdFile.close();
289	newMdFile.close();
290	}
291