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.
#	User	Rev	Content
1	chuckv	950	/* 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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
9	chuckv	950	* notice, this list of conditions and the following disclaimer.
10			*
11	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
12	chuckv	950	* 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	gezelter	1390	* 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	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
39			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
40	chuckv	950	*
41			* randomBuilder.cpp
42			*
43			* Created by Charles F. Vardeman II on 10 Apr 2006.
44			* @author Charles F. Vardeman II
45	gezelter	1442	* @version $Id$
46	chuckv	950	*
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	gezelter	1390	using namespace OpenMD;
73	chuckv	950
74			void createMdFile(const std::string&oldMdFileName,
75			const std::string&newMdFileName,
76	gezelter	1065	std::vector<int> nMol);
77	chuckv	950
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	gezelter	1062	std::string outputFileName;
86	chuckv	950	Lattice *simpleLat;
87	tim	963	RealType latticeConstant;
88			std::vector<RealType> lc;
89			const RealType rhoConvertConst = 1.661;
90			RealType density;
91	gezelter	1062	int nx, ny, nz;
92	chuckv	950	Mat3x3d hmat;
93			MoLocator *locator;
94			std::vector<Vector3d> latticePos;
95			std::vector<Vector3d> latticeOrt;
96	gezelter	1065	int nMolPerCell;
97	chuckv	950	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	gezelter	1067	latticeType = "FCC";
107	chuckv	950
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	gezelter	1065	nMolPerCell = simpleLat->getNumSitesPerCell();
117	chuckv	950
118	gezelter	1062	//get the number of unit cells in each direction:
119
120	chuckv	950	nx = args_info.nx_arg;
121
122			if (nx <= 0) {
123	gezelter	1062	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	chuckv	950	}
128
129			ny = args_info.ny_arg;
130
131			if (ny <= 0) {
132	gezelter	1062	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	chuckv	950	}
137
138			nz = args_info.nz_arg;
139
140			if (nz <= 0) {
141	gezelter	1062	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	chuckv	950	}
146
147	gezelter	1065	int nSites = nMolPerCell * nx * ny * nz;
148
149	chuckv	950	//get input file name
150			if (args_info.inputs_num)
151			inputFileName = args_info.inputs[0];
152			else {
153	gezelter	1062	sprintf(painCave.errMsg, "No input .md file name was specified "
154			"on the command line");
155			painCave.isFatal = 1;
156			simError();
157	chuckv	950	}
158
159			//parse md file and set up the system
160	gezelter	1062
161	chuckv	950	SimCreator oldCreator;
162			SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
163			Globals* simParams = oldInfo->getSimParams();
164
165	gezelter	1065	// 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	gezelter	1062	painCave.isFatal = 1;
217			simError();
218	chuckv	950	}
219
220	gezelter	1065	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	chuckv	950
227	gezelter	1065	// recompute actual mol fractions and perform final sanity check:
228	chuckv	950
229	gezelter	1065	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	gezelter	1062
240	gezelter	1065	if (totalMolecules != nSites) {
241			sprintf(painCave.errMsg, "Computed total number of molecules is not equal "
242			"to the number of lattice sites!");
243	gezelter	1062	painCave.isFatal = 1;
244			simError();
245	chuckv	950	}
246	gezelter	1065
247			latticeConstant = pow(rhoConvertConst * nMolPerCell * avgMass / density,
248	tim	963	(RealType)(1.0 / 3.0));
249	gezelter	1065
250	gezelter	1062	// Set the lattice constant
251	gezelter	1065
252	chuckv	950	lc.push_back(latticeConstant);
253			simpleLat->setLatticeConstant(lc);
254	gezelter	1065
255	gezelter	1062	// 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	chuckv	950
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	gezelter	1062	// Get the position of the cell sites
269
270	chuckv	950	simpleLat->getLatticePointsPos(latticePos, i, j, k);
271
272	gezelter	1065	for(int l = 0; l < nMolPerCell; l++) {
273	gezelter	1062	sites.push_back(latticePos[l]);
274			orientations.push_back(latticeOrt[l]);
275	chuckv	950	}
276			}
277			}
278			}
279
280	gezelter	1062	outputFileName = args_info.output_arg;
281	chuckv	950
282	gezelter	1062	// create a new .md file on the fly which corrects the number of molecules
283	chuckv	950
284	gezelter	1065	createMdFile(inputFileName, outputFileName, nMol);
285	chuckv	950
286	gezelter	1062	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	chuckv	950	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	gezelter	1062	// Set Hmat
310	chuckv	950
311	gezelter	1062	newInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat(hmat);
312	chuckv	950
313	gezelter	1062	// place the molecules
314
315			// Randomize a vector of ints:
316	chuckv	950
317	gezelter	1062	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	chuckv	950
321			Molecule* mol;
322			int l = 0;
323	gezelter	1062	for (int i = 0; i < nComponents; i++){
324			locator = new MoLocator(newInfo->getMoleculeStamp(i),
325			newInfo->getForceField());
326	gezelter	1065	for (int n = 0; n < nMol.at(i); n++) {
327	gezelter	1062	mol = newInfo->getMoleculeByGlobalIndex(l);
328			locator->placeMol(sites[ids[l]], orientations[ids[l]], mol);
329			l++;
330			}
331	chuckv	950	}
332	gezelter	1062
333			// Create DumpWriter and write out the coordinates
334	chuckv	950
335	gezelter	1062	writer = new DumpWriter(newInfo, outputFileName);
336	chuckv	950
337			if (writer == NULL) {
338	gezelter	1062	sprintf(painCave.errMsg, "error in creating DumpWriter");
339			painCave.isFatal = 1;
340			simError();
341	chuckv	950	}
342
343	gezelter	1062	writer->writeDump();
344
345			// deleting the writer will put the closing at the end of the dump file.
346
347			delete writer;
348
349	gezelter	1390	sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
350	gezelter	1065	"generated.\n", outputFileName.c_str());
351	gezelter	1062	painCave.isFatal = 0;
352			simError();
353	chuckv	950	return 0;
354			}
355
356	gezelter	1062	void createMdFile(const std::string&oldMdFileName,
357			const std::string&newMdFileName,
358	gezelter	1065	std::vector<int> nMol) {
359	chuckv	950	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	gezelter	1062
366	chuckv	950	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	gezelter	1065	if(i<nMol.size()){
377			sprintf(buffer, "\tnMol = %i;", nMol.at(i));
378	chuckv	950	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	gezelter	1077
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	chuckv	950	}
399