applications/randomBuilder/randomBuilder.cpp

/* Copyright (c) 2006 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).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
 *
 *  randomBuilder.cpp
 *
 *  Created by Charles F. Vardeman II on 10 Apr 2006.
 *  @author  Charles F. Vardeman II
 *  @version $Id$
 *
 */


#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <string>
#include <map>
#include <fstream>

#include "applications/randomBuilder/randomBuilderCmd.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,
                  std::vector<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)

  std::vector<Component*> components = simParams->getComponents();
  std::vector<RealType> molFractions;
  std::vector<RealType> molecularMasses;
  std::vector<int> nMol;
  int nComponents = components.size();

  if (nComponents == 1) {
    molFractions.push_back(1.0);    
  } else {
    if (args_info.molFraction_given == nComponents) {
      for (int i = 0; i < nComponents; i++) {
        molFractions.push_back(args_info.molFraction_arg[i]);
      }
    } else if (args_info.molFraction_given == nComponents-1) {
      RealType remainingFraction = 1.0;
      for (int i = 0; i < nComponents-1; i++) {
        molFractions.push_back(args_info.molFraction_arg[i]);
        remainingFraction -= molFractions[i];
      }
      molFractions.push_back(remainingFraction);
    } else {    
      sprintf(painCave.errMsg, "randomBuilder can't figure out molFractions "
              "for all of the components in the <MetaData> block.");
      painCave.isFatal = 1;
      simError();
    }
  }

  // do some sanity checking:
  
  RealType totalFraction = 0.0;

  for (int i = 0; i < nComponents; i++) {
    if (molFractions.at(i) < 0.0) {
      sprintf(painCave.errMsg, "One of the requested molFractions was"
              " less than zero!");
      painCave.isFatal = 1;
      simError();
    }
    if (molFractions.at(i) > 1.0) {
      sprintf(painCave.errMsg, "One of the requested molFractions was"
              " greater than one!");
      painCave.isFatal = 1;
      simError();
    }
    totalFraction += molFractions.at(i);
  }
  if (abs(totalFraction - 1.0) > 1e-6) {
    sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
    painCave.isFatal = 1;
    simError();
  }

  int remaining = nSites;
  for (int i=0; i < nComponents-1; i++) {    
    nMol.push_back(int((RealType)nSites * molFractions.at(i)));
    remaining -= nMol.at(i);
  }
  nMol.push_back(remaining);

  // recompute actual mol fractions and perform final sanity check:

  int totalMolecules = 0;
  RealType totalMass = 0.0;
  for (int i=0; i < nComponents; i++) {
    molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
    totalMolecules += nMol.at(i);
    molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
                                         oldInfo->getForceField()));
    totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
  }
  RealType avgMass = totalMass / (RealType) totalMolecules;

  if (totalMolecules != nSites) {
    sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
            "to the number of lattice sites!");
    painCave.isFatal = 1;
    simError();
  }
     
  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, nMol);

  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

  // Randomize a vector of ints:

  vector<int> ids;
  for (int i = 0; i < sites.size(); i++) ids.push_back(i);
  std::random_shuffle(ids.begin(), ids.end());

  Molecule* mol;
  int l = 0;
  for (int i = 0; i < nComponents; i++){
    locator = new MoLocator(newInfo->getMoleculeStamp(i), 
                            newInfo->getForceField());
    for (int n = 0; n < nMol.at(i); n++) {
      mol = newInfo->getMoleculeByGlobalIndex(l);
      locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
      l++;
    }
  }
  
  // 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, 
                  std::vector<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);
 
  int i = 0;
  while (!oldMdFile.eof()) {
    
    //correct molecule number
    if (strstr(buffer, "nMol") != NULL) {
      if(i<nMol.size()){
        sprintf(buffer, "\tnMol = %i;", nMol.at(i));
        newMdFile << buffer << std::endl;
        i++;
      }
    } else
      newMdFile << buffer << std::endl;
    
    oldMdFile.getline(buffer, MAXLEN);
  }
  
  oldMdFile.close();
  newMdFile.close();

  if (i != nMol.size()) {
    sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
            "\tstatements in component blocks.  Make sure that all\n"
            "\tcomponents in the template file have nMol=1");
    painCave.isFatal = 1;
    simError();
  }

}

Revision:	1725
Committed:	Sat May 26 18:13:43 2012 UTC (12 years, 11 months ago) by gezelter
File size:	11554 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	/* Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
2	*
3	* The University of Notre Dame grants you ("Licensee") a
4	* non-exclusive, royalty free, license to use, modify and
5	* redistribute this software in source and binary code form, provided
6	* that the following conditions are met:
7	*
8	* 1. Redistributions of source code must retain the above copyright
9	* notice, this list of conditions and the following disclaimer.
10	*
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the
14	* distribution.
15	*
16	* This software is provided "AS IS," without a warranty of any
17	* kind. All express or implied conditions, representations and
18	* warranties, including any implied warranty of merchantability,
19	* fitness for a particular purpose or non-infringement, are hereby
20	* excluded. The University of Notre Dame and its licensors shall not
21	* be liable for any damages suffered by licensee as a result of
22	* using, modifying or distributing the software or its
23	* derivatives. In no event will the University of Notre Dame or its
24	* licensors be liable for any lost revenue, profit or data, or for
25	* direct, indirect, special, consequential, incidental or punitive
26	* damages, however caused and regardless of the theory of liability,
27	* arising out of the use of or inability to use software, even if the
28	* University of Notre Dame has been advised of the possibility of
29	* such damages.
30	*
31	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32	* research, please cite the appropriate papers when you publish your
33	* work. Good starting points are:
34	*
35	* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36	* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
38	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
39	* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
40	*
41	* randomBuilder.cpp
42	*
43	* Created by Charles F. Vardeman II on 10 Apr 2006.
44	* @author Charles F. Vardeman II
45	* @version $Id$
46	*
47	*/
48
49
50	#include <cstdlib>
51	#include <cstdio>
52	#include <cstring>
53	#include <cmath>
54	#include <iostream>
55	#include <string>
56	#include <map>
57	#include <fstream>
58
59	#include "applications/randomBuilder/randomBuilderCmd.h"
60	#include "lattice/LatticeFactory.hpp"
61	#include "utils/MoLocator.hpp"
62	#include "lattice/Lattice.hpp"
63	#include "brains/Register.hpp"
64	#include "brains/SimInfo.hpp"
65	#include "brains/SimCreator.hpp"
66	#include "io/DumpWriter.hpp"
67	#include "math/Vector3.hpp"
68	#include "math/SquareMatrix3.hpp"
69	#include "utils/StringUtils.hpp"
70
71	using namespace std;
72	using namespace OpenMD;
73
74	void createMdFile(const std::string&oldMdFileName,
75	const std::string&newMdFileName,
76	std::vector<int> nMol);
77
78	int main(int argc, char *argv []) {
79
80	registerLattice();
81
82	gengetopt_args_info args_info;
83	std::string latticeType;
84	std::string inputFileName;
85	std::string outputFileName;
86	Lattice *simpleLat;
87	RealType latticeConstant;
88	std::vector<RealType> lc;
89	const RealType rhoConvertConst = 1.661;
90	RealType density;
91	int nx, ny, nz;
92	Mat3x3d hmat;
93	MoLocator *locator;
94	std::vector<Vector3d> latticePos;
95	std::vector<Vector3d> latticeOrt;
96	int nMolPerCell;
97	DumpWriter *writer;
98
99	// parse command line arguments
100	if (cmdline_parser(argc, argv, &args_info) != 0)
101	exit(1);
102
103	density = args_info.density_arg;
104
105	//get lattice type
106	latticeType = "FCC";
107
108	simpleLat = LatticeFactory::getInstance()->createLattice(latticeType);
109
110	if (simpleLat == NULL) {
111	sprintf(painCave.errMsg, "Lattice Factory can not create %s lattice\n",
112	latticeType.c_str());
113	painCave.isFatal = 1;
114	simError();
115	}
116	nMolPerCell = simpleLat->getNumSitesPerCell();
117
118	//get the number of unit cells in each direction:
119
120	nx = args_info.nx_arg;
121
122	if (nx <= 0) {
123	sprintf(painCave.errMsg, "The number of unit cells in the x direction "
124	"must be greater than 0.");
125	painCave.isFatal = 1;
126	simError();
127	}
128
129	ny = args_info.ny_arg;
130
131	if (ny <= 0) {
132	sprintf(painCave.errMsg, "The number of unit cells in the y direction "
133	"must be greater than 0.");
134	painCave.isFatal = 1;
135	simError();
136	}
137
138	nz = args_info.nz_arg;
139
140	if (nz <= 0) {
141	sprintf(painCave.errMsg, "The number of unit cells in the z direction "
142	"must be greater than 0.");
143	painCave.isFatal = 1;
144	simError();
145	}
146
147	int nSites = nMolPerCell * nx * ny * nz;
148
149	//get input file name
150	if (args_info.inputs_num)
151	inputFileName = args_info.inputs[0];
152	else {
153	sprintf(painCave.errMsg, "No input .md file name was specified "
154	"on the command line");
155	painCave.isFatal = 1;
156	simError();
157	}
158
159	//parse md file and set up the system
160
161	SimCreator oldCreator;
162	SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
163	Globals* simParams = oldInfo->getSimParams();
164
165	// Calculate lattice constant (in Angstroms)
166
167	std::vector<Component*> components = simParams->getComponents();
168	std::vector<RealType> molFractions;
169	std::vector<RealType> molecularMasses;
170	std::vector<int> nMol;
171	int nComponents = components.size();
172
173	if (nComponents == 1) {
174	molFractions.push_back(1.0);
175	} else {
176	if (args_info.molFraction_given == nComponents) {
177	for (int i = 0; i < nComponents; i++) {
178	molFractions.push_back(args_info.molFraction_arg[i]);
179	}
180	} else if (args_info.molFraction_given == nComponents-1) {
181	RealType remainingFraction = 1.0;
182	for (int i = 0; i < nComponents-1; i++) {
183	molFractions.push_back(args_info.molFraction_arg[i]);
184	remainingFraction -= molFractions[i];
185	}
186	molFractions.push_back(remainingFraction);
187	} else {
188	sprintf(painCave.errMsg, "randomBuilder can't figure out molFractions "
189	"for all of the components in the <MetaData> block.");
190	painCave.isFatal = 1;
191	simError();
192	}
193	}
194
195	// do some sanity checking:
196
197	RealType totalFraction = 0.0;
198
199	for (int i = 0; i < nComponents; i++) {
200	if (molFractions.at(i) < 0.0) {
201	sprintf(painCave.errMsg, "One of the requested molFractions was"
202	" less than zero!");
203	painCave.isFatal = 1;
204	simError();
205	}
206	if (molFractions.at(i) > 1.0) {
207	sprintf(painCave.errMsg, "One of the requested molFractions was"
208	" greater than one!");
209	painCave.isFatal = 1;
210	simError();
211	}
212	totalFraction += molFractions.at(i);
213	}
214	if (abs(totalFraction - 1.0) > 1e-6) {
215	sprintf(painCave.errMsg, "The sum of molFractions was not close enough to 1.0");
216	painCave.isFatal = 1;
217	simError();
218	}
219
220	int remaining = nSites;
221	for (int i=0; i < nComponents-1; i++) {
222	nMol.push_back(int((RealType)nSites * molFractions.at(i)));
223	remaining -= nMol.at(i);
224	}
225	nMol.push_back(remaining);
226
227	// recompute actual mol fractions and perform final sanity check:
228
229	int totalMolecules = 0;
230	RealType totalMass = 0.0;
231	for (int i=0; i < nComponents; i++) {
232	molFractions[i] = (RealType)(nMol.at(i))/(RealType)nSites;
233	totalMolecules += nMol.at(i);
234	molecularMasses.push_back(getMolMass(oldInfo->getMoleculeStamp(i),
235	oldInfo->getForceField()));
236	totalMass += (RealType)(nMol.at(i)) * molecularMasses.at(i);
237	}
238	RealType avgMass = totalMass / (RealType) totalMolecules;
239
240	if (totalMolecules != nSites) {
241	sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
242	"to the number of lattice sites!");
243	painCave.isFatal = 1;
244	simError();
245	}
246
247	latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
248	(RealType)(1.0 / 3.0));
249
250	// Set the lattice constant
251
252	lc.push_back(latticeConstant);
253	simpleLat->setLatticeConstant(lc);
254
255	// Calculate the lattice sites and fill the lattice vector.
256
257	// Get the standard orientations of the cell sites
258
259	latticeOrt = simpleLat->getLatticePointsOrt();
260
261	vector<Vector3d> sites;
262	vector<Vector3d> orientations;
263
264	for(int i = 0; i < nx; i++) {
265	for(int j = 0; j < ny; j++) {
266	for(int k = 0; k < nz; k++) {
267
268	// Get the position of the cell sites
269
270	simpleLat->getLatticePointsPos(latticePos, i, j, k);
271
272	for(int l = 0; l < nMolPerCell; l++) {
273	sites.push_back(latticePos[l]);
274	orientations.push_back(latticeOrt[l]);
275	}
276	}
277	}
278	}
279
280	outputFileName = args_info.output_arg;
281
282	// create a new .md file on the fly which corrects the number of molecules
283
284	createMdFile(inputFileName, outputFileName, nMol);
285
286	if (oldInfo != NULL)
287	delete oldInfo;
288
289	// We need to read in the new SimInfo object, then Parse the
290	// md file and set up the system
291
292	SimCreator newCreator;
293	SimInfo* newInfo = newCreator.createSim(outputFileName, false);
294
295	// fill Hmat
296
297	hmat(0, 0) = nx * latticeConstant;
298	hmat(0, 1) = 0.0;
299	hmat(0, 2) = 0.0;
300
301	hmat(1, 0) = 0.0;
302	hmat(1, 1) = ny * latticeConstant;
303	hmat(1, 2) = 0.0;
304
305	hmat(2, 0) = 0.0;
306	hmat(2, 1) = 0.0;
307	hmat(2, 2) = nz * latticeConstant;
308
309	// Set Hmat
310
311	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
312
313	// place the molecules
314
315	// Randomize a vector of ints:
316
317	vector<int> ids;
318	for (int i = 0; i < sites.size(); i++) ids.push_back(i);
319	std::random_shuffle(ids.begin(), ids.end());
320
321	Molecule* mol;
322	int l = 0;
323	for (int i = 0; i < nComponents; i++){
324	locator = new MoLocator(newInfo->getMoleculeStamp(i),
325	newInfo->getForceField());
326	for (int n = 0; n < nMol.at(i); n++) {
327	mol = newInfo->getMoleculeByGlobalIndex(l);
328	locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
329	l++;
330	}
331	}
332
333	// Create DumpWriter and write out the coordinates
334
335	writer = new DumpWriter(newInfo, outputFileName);
336
337	if (writer == NULL) {
338	sprintf(painCave.errMsg, "error in creating DumpWriter");
339	painCave.isFatal = 1;
340	simError();
341	}
342
343	writer->writeDump();
344
345	// deleting the writer will put the closing at the end of the dump file.
346
347	delete writer;
348
349	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
350	"generated.\n", outputFileName.c_str());
351	painCave.isFatal = 0;
352	simError();
353	return 0;
354	}
355
356	void createMdFile(const std::string&oldMdFileName,
357	const std::string&newMdFileName,
358	std::vector<int> nMol) {
359	ifstream oldMdFile;
360	ofstream newMdFile;
361	const int MAXLEN = 65535;
362	char buffer[MAXLEN];
363
364	//create new .md file based on old .md file
365
366	oldMdFile.open(oldMdFileName.c_str());
367	newMdFile.open(newMdFileName.c_str());
368
369	oldMdFile.getline(buffer, MAXLEN);
370
371	int i = 0;
372	while (!oldMdFile.eof()) {
373
374	//correct molecule number
375	if (strstr(buffer, "nMol") != NULL) {
376	if(i<nMol.size()){
377	sprintf(buffer, "\tnMol = %i;", nMol.at(i));
378	newMdFile << buffer << std::endl;
379	i++;
380	}
381	} else
382	newMdFile << buffer << std::endl;
383
384	oldMdFile.getline(buffer, MAXLEN);
385	}
386
387	oldMdFile.close();
388	newMdFile.close();
389
390	if (i != nMol.size()) {
391	sprintf(painCave.errMsg, "Couldn't replace the correct number of nMol\n"
392	"\tstatements in component blocks. Make sure that all\n"
393	"\tcomponents in the template file have nMol=1");
394	painCave.isFatal = 1;
395	simError();
396	}
397
398	}
399