applications/nanoparticleBuilder/icosahedralBuilder.cpp

/*
 * Copyright (c) 2013 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, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

#include "config.h"
#include "icosahedralBuilderCmd.h"
#include "utils/MoLocator.hpp"
#include "utils/Tuple.hpp"
#include "brains/SimInfo.hpp"
#include "brains/SimCreator.hpp"
#include "io/DumpWriter.hpp"
#include "math/SquareMatrix3.hpp"

using namespace OpenMD;
using namespace std;

//
// Create Mackay icosaheron structure.
//
// Heavily modified from a code created by: Yanting Wang    07/21/2003
//

vector<Vector3d> Points;
vector<pair<int, int> > Edges;
vector<tuple3<int, int, int> > Facets;
vector<Vector3d> Basis; // Basis vectors of the edges

//
// function np
//
// Calculate number of particles on the nth layer.
//

int np( int n ) {
  if( n<0 ) return -1;
  else if( n==0 ) return 1;
  else if( n==1 ) return 12;
  else if( n==2 ) return 42;
  else {
    int count = 0;    
    count += 12;   // edge particles
    count += (n-1)*30;   // side particles
    for( int i = 1; i <= n-2; i++ ) count += i*20;   // body particles
    return count;
  }
}

//
// function init
//
// Initialize some constants.
//

void init() {
  
  Basis.clear();
  Edges.clear();
  Facets.clear();
  Points.clear();
  
  //
  // Initialize Basis vectors.
  //
  const RealType HGR = ( sqrt(5.0) + 1.0 ) / 4.0;   // half of the golden ratio
  
  Basis.push_back( Vector3d(  HGR,  0.0,  0.5 ));
  Basis.push_back( Vector3d(  HGR,  0.0, -0.5 ));
  Basis.push_back( Vector3d(  0.5,  HGR,  0.0 ));
  Basis.push_back( Vector3d( -0.5,  HGR,  0.0 ));
  Basis.push_back( Vector3d(  0.0,  0.5,  HGR ));
  Basis.push_back( Vector3d(  0.0, -0.5,  HGR ));
  Basis.push_back( Vector3d(  0.5, -HGR,  0.0 ));
  Basis.push_back( Vector3d(  0.0,  0.5, -HGR ));
  Basis.push_back( Vector3d( -HGR,  0.0,  0.5 ));
  Basis.push_back( Vector3d(  0.0, -0.5, -HGR ));
  Basis.push_back( Vector3d( -HGR,  0.0, -0.5 ));
  Basis.push_back( Vector3d( -0.5, -HGR,  0.0 ));
  
  //
  // Initialize 30 edges
  //
  
  Edges.push_back(std::make_pair(0, 1));
  Edges.push_back(std::make_pair(0, 2));
  Edges.push_back(std::make_pair(0, 4));
  Edges.push_back(std::make_pair(0, 5));
  Edges.push_back(std::make_pair(0, 6));
  
  Edges.push_back(std::make_pair(10, 3));
  Edges.push_back(std::make_pair(10, 7));
  Edges.push_back(std::make_pair(10, 8));
  Edges.push_back(std::make_pair(10, 9));
  Edges.push_back(std::make_pair(10, 11));
  
  Edges.push_back(std::make_pair(1, 2));
  Edges.push_back(std::make_pair(1, 6));
  Edges.push_back(std::make_pair(1, 7));
  Edges.push_back(std::make_pair(1, 9));
  
  Edges.push_back(std::make_pair(8, 3));
  Edges.push_back(std::make_pair(8, 4));
  Edges.push_back(std::make_pair(8, 5));
  Edges.push_back(std::make_pair(8, 11));
  
  Edges.push_back(std::make_pair(2, 3));
  Edges.push_back(std::make_pair(2, 4));
  Edges.push_back(std::make_pair(2, 7));
  
  Edges.push_back(std::make_pair(11, 5));
  Edges.push_back(std::make_pair(11, 6));
  Edges.push_back(std::make_pair(11, 9));
  
  Edges.push_back(std::make_pair(6, 5));
  Edges.push_back(std::make_pair(6, 9));
  
  Edges.push_back(std::make_pair(3, 4));
  Edges.push_back(std::make_pair(3, 7));
  
  Edges.push_back(std::make_pair(7, 9));
  
  Edges.push_back(std::make_pair(5, 4));
  
  //
  // Initialize 20 facets
  //
  
  Facets.push_back(make_tuple3(0, 1, 2));
  Facets.push_back(make_tuple3(0, 2, 4));
  Facets.push_back(make_tuple3(0, 4, 5));
  Facets.push_back(make_tuple3(0, 5, 6));
  Facets.push_back(make_tuple3(0, 1, 6));
  
  Facets.push_back(make_tuple3(10, 3, 7));
  Facets.push_back(make_tuple3(10, 3, 8));
  Facets.push_back(make_tuple3(10, 8, 11));
  Facets.push_back(make_tuple3(10, 9, 11));
  Facets.push_back(make_tuple3(10, 7, 9));
                                    
  Facets.push_back(make_tuple3(1, 2, 7));
  Facets.push_back(make_tuple3(1, 7, 9));
  Facets.push_back(make_tuple3(1, 6, 9));
  
  Facets.push_back(make_tuple3(8, 5, 11));
  Facets.push_back(make_tuple3(8, 4, 5));
  Facets.push_back(make_tuple3(8, 3, 4));
  
  Facets.push_back(make_tuple3(2, 3, 7));
  Facets.push_back(make_tuple3(2, 3, 4));
  
  Facets.push_back(make_tuple3(11, 5, 6));
  Facets.push_back(make_tuple3(11, 6, 9));
}

//
// function ih
//
// Create nth layer particles.
// The distance between nearest neighbors has the unit length of 1.

vector<Vector3d> ih( int n ) {

  if( n < 0 ) return Points; 
  
  if( n==0 ) {
    // center particle only
    Points.push_back(Vector3d( 0.0, 0.0, 0.0 ));
    return Points;
  } 
  
  //
  // Generate edge particles
  //
  for( vector<Vector3d>::iterator i = Basis.begin(); i != Basis.end(); ++i ) {
    
    Points.push_back( (*i) * RealType(n) );
  }
  
  //
  // Generate side particles
  //
  
  if( n<2 ) return Points;
  
  for( vector<pair<int,int> >::iterator i=Edges.begin(); 
       i != Edges.end(); ++i ) {
    
    Vector3d e1 = Basis[ (*i).first  ] * RealType(n);
    Vector3d e2 = Basis[ (*i).second ] * RealType(n);
    
    for( int j = 1; j <= n-1; j++ ) {
      Points.push_back( e1 + (e2-e1) * RealType(j) / RealType(n));
    }      
  }
  
  //
  // Generate body particles
  //
  
  if( n<3 ) return Points;
  
  for( vector<tuple3<int,int,int> >::iterator i = Facets.begin();
       i != Facets.end(); ++i) {
    
    Vector3d e1 = Basis[ (*i).first  ] * RealType(n);
    Vector3d e2 = Basis[ (*i).second ] * RealType(n);
    Vector3d e3 = Basis[ (*i).third  ] * RealType(n);
    
    for( int j=1; j<=n-2; j++ ) {
      
      Vector3d v1 = e1 + (e2-e1) * RealType(j+1) / RealType(n);
      Vector3d v2 = e1 + (e3-e1) * RealType(j+1) / RealType(n);
      
      for( int k=1; k<=j; k++ ) {
        Points.push_back(v1 + (v2-v1) * RealType(k) / RealType(j+1));
      }
    }
  }
  return Points;
}


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);

  unsigned int i = 0;
  while (!oldMdFile.eof()) {

    //correct molecule number
    if (strstr(buffer, "nMol") != NULL) {
      sprintf(buffer, "\tnMol = %i;", nMol);
      newMdFile << buffer << std::endl;
    } else {
      newMdFile << buffer << std::endl;
    }
    
    oldMdFile.getline(buffer, MAXLEN);
  }
  
  oldMdFile.close();
  newMdFile.close();
}

int main(int argc, char *argv []) {
  
  gengetopt_args_info args_info;
  std::string inputFileName;
  std::string outputFileName;
  
  MoLocator* locator;
  RealType latticeConstant;
  int nShells;
 
  DumpWriter *writer;

  init();
  
  // Parse Command Line Arguments
  if (cmdline_parser(argc, argv, &args_info) != 0)
    exit(1);
         
  /* 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;
    cmdline_parser_print_help();
    simError();
  }

  if (args_info.shells_given) {   
    nShells = args_info.shells_arg;    
    if( nShells < 0 ) {
      sprintf(painCave.errMsg, "icosahedralBuilder:  The number of shells\n"
              "\tmust be greater than or equal to zero.");
      painCave.isFatal = 1;
      cmdline_parser_print_help();
      simError();
    }
  } else {
    sprintf(painCave.errMsg, "icosahedralBuilder:  The number of shells\n"
            "\tis required to build a Mackay Icosahedron.");
    painCave.isFatal = 1;
    cmdline_parser_print_help();
    simError();
  }

  if (args_info.latticeConstant_given) {    
    latticeConstant = args_info.latticeConstant_arg;
  } else  {
    
    int count=0;
    for( int i = 0; i <= nShells; i++ ) count += np( i );
    sprintf(painCave.errMsg, "icosahedralBuilder:  No lattice constant\n"
            "\tgiven.  Total number of atoms in a Mackay Icosahedron with\n"
            "\t%d shells is %d.", nShells, count);
    painCave.isFatal = 1;
    painCave.severity = OPENMD_INFO;
    cmdline_parser_print_help();
    simError();
  }

  /* parse md file and set up the system */
  SimCreator oldCreator;
  SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);

  Globals* simParams = oldInfo->getSimParams();

  //generate the coordinates
  for( int i = 0; i <= nShells; i++ ) ih( i );

  outputFileName = args_info.output_arg;
   
  //creat new .md file on fly which corrects the number of molecule    

  createMdFile(inputFileName, outputFileName, Points.size());
  
  if (oldInfo != NULL)
    delete oldInfo;
  
  SimCreator newCreator;
  SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
    
  // Place molecules
  Molecule* mol;
  SimInfo::MoleculeIterator mi;
  mol = NewInfo->beginMolecule(mi);

  int l = 0;

  locator = new MoLocator(NewInfo->getMoleculeStamp(0), 
                          NewInfo->getForceField());
    
  Vector3d boxMax;
  Vector3d boxMin;

  for (int n = 0; n < Points.size(); n++) {
    mol = NewInfo->getMoleculeByGlobalIndex(l);
    Vector3d location = Points[n] * latticeConstant;
    Vector3d orientation = Vector3d(0, 0, 1.0);

    if (n == 0) {
      boxMax = location;
      boxMin = location;
    } else {
      for (int i = 0; i < 3; i++) {
        boxMax[i] = max(boxMax[i], location[i]);
        boxMin[i] = min(boxMin[i], location[i]);
      }
    }
    
    locator->placeMol(location, orientation, mol);
    l++;
  }

  Mat3x3d boundingBox;
  boundingBox(0,0) = 10.0*(boxMax[0] - boxMin[0]);
  boundingBox(1,1) = 10.0*(boxMax[1] - boxMin[1]);
  boundingBox(2,2) = 10.0*(boxMax[2] - boxMin[2]);
 
  //set Hmat
  NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat( boundingBox );
  
  //create dumpwriter and write out the coordinates
  writer = new DumpWriter(NewInfo, outputFileName);
  
  if (writer == NULL) {
    sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
    painCave.isFatal = 1;
    simError();
  }
  
  writer->writeDump();

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

  delete writer;

  // cleanup a by calling sim error.....
  sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
          "generated.\n", outputFileName.c_str());
  painCave.isFatal = 0;
  painCave.severity = OPENMD_INFO;
  simError();
  return 0;
}
Revision:	1861
Committed:	Tue Apr 9 19:45:54 2013 UTC (12 years, 1 month ago) by gezelter
File size:	12518 byte(s)
Log Message:	Added a Harmonic Torsion Type, fixed some bugs in RNEMD and waterReplacer.
#	User	Rev	Content
1	gezelter	1860	/*
2			* Copyright (c) 2013 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 "config.h"
44			#include "icosahedralBuilderCmd.h"
45			#include "utils/MoLocator.hpp"
46			#include "utils/Tuple.hpp"
47			#include "brains/SimInfo.hpp"
48			#include "brains/SimCreator.hpp"
49			#include "io/DumpWriter.hpp"
50			#include "math/SquareMatrix3.hpp"
51
52			using namespace OpenMD;
53			using namespace std;
54
55			//
56			// Create Mackay icosaheron structure.
57			//
58			// Heavily modified from a code created by: Yanting Wang 07/21/2003
59			//
60
61			vector<Vector3d> Points;
62	gezelter	1861	vector<pair<int, int> > Edges;
63	gezelter	1860	vector<tuple3<int, int, int> > Facets;
64			vector<Vector3d> Basis; // Basis vectors of the edges
65
66			//
67			// function np
68			//
69			// Calculate number of particles on the nth layer.
70			//
71
72			int np( int n ) {
73			if( n<0 ) return -1;
74			else if( n==0 ) return 1;
75			else if( n==1 ) return 12;
76			else if( n==2 ) return 42;
77			else {
78			int count = 0;
79			count += 12; // edge particles
80			count += (n-1)*30; // side particles
81			for( int i = 1; i <= n-2; i++ ) count += i*20; // body particles
82			return count;
83			}
84			}
85
86			//
87			// function init
88			//
89			// Initialize some constants.
90			//
91
92			void init() {
93
94			Basis.clear();
95			Edges.clear();
96			Facets.clear();
97			Points.clear();
98
99			//
100			// Initialize Basis vectors.
101			//
102	gezelter	1861	const RealType HGR = ( sqrt(5.0) + 1.0 ) / 4.0; // half of the golden ratio
103	gezelter	1860
104	gezelter	1861	Basis.push_back( Vector3d( HGR, 0.0, 0.5 ));
105			Basis.push_back( Vector3d( HGR, 0.0, -0.5 ));
106			Basis.push_back( Vector3d( 0.5, HGR, 0.0 ));
107			Basis.push_back( Vector3d( -0.5, HGR, 0.0 ));
108			Basis.push_back( Vector3d( 0.0, 0.5, HGR ));
109			Basis.push_back( Vector3d( 0.0, -0.5, HGR ));
110			Basis.push_back( Vector3d( 0.5, -HGR, 0.0 ));
111			Basis.push_back( Vector3d( 0.0, 0.5, -HGR ));
112			Basis.push_back( Vector3d( -HGR, 0.0, 0.5 ));
113			Basis.push_back( Vector3d( 0.0, -0.5, -HGR ));
114			Basis.push_back( Vector3d( -HGR, 0.0, -0.5 ));
115			Basis.push_back( Vector3d( -0.5, -HGR, 0.0 ));
116	gezelter	1860
117			//
118			// Initialize 30 edges
119			//
120
121			Edges.push_back(std::make_pair(0, 1));
122			Edges.push_back(std::make_pair(0, 2));
123			Edges.push_back(std::make_pair(0, 4));
124			Edges.push_back(std::make_pair(0, 5));
125			Edges.push_back(std::make_pair(0, 6));
126
127			Edges.push_back(std::make_pair(10, 3));
128			Edges.push_back(std::make_pair(10, 7));
129			Edges.push_back(std::make_pair(10, 8));
130			Edges.push_back(std::make_pair(10, 9));
131			Edges.push_back(std::make_pair(10, 11));
132
133			Edges.push_back(std::make_pair(1, 2));
134			Edges.push_back(std::make_pair(1, 6));
135			Edges.push_back(std::make_pair(1, 7));
136			Edges.push_back(std::make_pair(1, 9));
137
138			Edges.push_back(std::make_pair(8, 3));
139			Edges.push_back(std::make_pair(8, 4));
140			Edges.push_back(std::make_pair(8, 5));
141			Edges.push_back(std::make_pair(8, 11));
142
143			Edges.push_back(std::make_pair(2, 3));
144			Edges.push_back(std::make_pair(2, 4));
145			Edges.push_back(std::make_pair(2, 7));
146
147			Edges.push_back(std::make_pair(11, 5));
148			Edges.push_back(std::make_pair(11, 6));
149			Edges.push_back(std::make_pair(11, 9));
150
151			Edges.push_back(std::make_pair(6, 5));
152			Edges.push_back(std::make_pair(6, 9));
153
154			Edges.push_back(std::make_pair(3, 4));
155			Edges.push_back(std::make_pair(3, 7));
156
157			Edges.push_back(std::make_pair(7, 9));
158
159			Edges.push_back(std::make_pair(5, 4));
160
161			//
162			// Initialize 20 facets
163			//
164
165			Facets.push_back(make_tuple3(0, 1, 2));
166			Facets.push_back(make_tuple3(0, 2, 4));
167			Facets.push_back(make_tuple3(0, 4, 5));
168			Facets.push_back(make_tuple3(0, 5, 6));
169			Facets.push_back(make_tuple3(0, 1, 6));
170
171			Facets.push_back(make_tuple3(10, 3, 7));
172			Facets.push_back(make_tuple3(10, 3, 8));
173			Facets.push_back(make_tuple3(10, 8, 11));
174			Facets.push_back(make_tuple3(10, 9, 11));
175			Facets.push_back(make_tuple3(10, 7, 9));
176
177			Facets.push_back(make_tuple3(1, 2, 7));
178			Facets.push_back(make_tuple3(1, 7, 9));
179			Facets.push_back(make_tuple3(1, 6, 9));
180
181			Facets.push_back(make_tuple3(8, 5, 11));
182			Facets.push_back(make_tuple3(8, 4, 5));
183			Facets.push_back(make_tuple3(8, 3, 4));
184
185			Facets.push_back(make_tuple3(2, 3, 7));
186			Facets.push_back(make_tuple3(2, 3, 4));
187
188			Facets.push_back(make_tuple3(11, 5, 6));
189			Facets.push_back(make_tuple3(11, 6, 9));
190			}
191
192			//
193			// function ih
194			//
195			// Create nth layer particles.
196			// The distance between nearest neighbors has the unit length of 1.
197
198			vector<Vector3d> ih( int n ) {
199
200			if( n < 0 ) return Points;
201
202			if( n==0 ) {
203			// center particle only
204			Points.push_back(Vector3d( 0.0, 0.0, 0.0 ));
205			return Points;
206			}
207
208			//
209			// Generate edge particles
210			//
211			for( vector<Vector3d>::iterator i = Basis.begin(); i != Basis.end(); ++i ) {
212
213			Points.push_back( (i) RealType(n) );
214			}
215
216			//
217			// Generate side particles
218			//
219
220			if( n<2 ) return Points;
221
222			for( vector<pair<int,int> >::iterator i=Edges.begin();
223			i != Edges.end(); ++i ) {
224
225			Vector3d e1 = Basis[ (i).first ] RealType(n);
226			Vector3d e2 = Basis[ (i).second ] RealType(n);
227
228			for( int j = 1; j <= n-1; j++ ) {
229			Points.push_back( e1 + (e2-e1) * RealType(j) / RealType(n));
230			}
231			}
232
233			//
234			// Generate body particles
235			//
236
237			if( n<3 ) return Points;
238
239			for( vector<tuple3<int,int,int> >::iterator i = Facets.begin();
240			i != Facets.end(); ++i) {
241
242			Vector3d e1 = Basis[ (i).first ] RealType(n);
243			Vector3d e2 = Basis[ (i).second ] RealType(n);
244			Vector3d e3 = Basis[ (i).third ] RealType(n);
245
246			for( int j=1; j<=n-2; j++ ) {
247
248			Vector3d v1 = e1 + (e2-e1) * RealType(j+1) / RealType(n);
249			Vector3d v2 = e1 + (e3-e1) * RealType(j+1) / RealType(n);
250
251			for( int k=1; k<=j; k++ ) {
252			Points.push_back(v1 + (v2-v1) * RealType(k) / RealType(j+1));
253			}
254			}
255			}
256			return Points;
257			}
258
259
260			void createMdFile(const std::string&oldMdFileName,
261			const std::string&newMdFileName,
262			int nMol) {
263			ifstream oldMdFile;
264			ofstream newMdFile;
265			const int MAXLEN = 65535;
266			char buffer[MAXLEN];
267
268			//create new .md file based on old .md file
269			oldMdFile.open(oldMdFileName.c_str());
270			newMdFile.open(newMdFileName.c_str());
271			oldMdFile.getline(buffer, MAXLEN);
272
273			unsigned int i = 0;
274			while (!oldMdFile.eof()) {
275
276			//correct molecule number
277			if (strstr(buffer, "nMol") != NULL) {
278			sprintf(buffer, "\tnMol = %i;", nMol);
279			newMdFile << buffer << std::endl;
280			} else {
281			newMdFile << buffer << std::endl;
282			}
283
284			oldMdFile.getline(buffer, MAXLEN);
285			}
286
287			oldMdFile.close();
288			newMdFile.close();
289			}
290
291			int main(int argc, char *argv []) {
292
293			gengetopt_args_info args_info;
294			std::string inputFileName;
295			std::string outputFileName;
296
297			MoLocator* locator;
298			RealType latticeConstant;
299			int nShells;
300
301			DumpWriter *writer;
302
303			init();
304
305			// Parse Command Line Arguments
306			if (cmdline_parser(argc, argv, &args_info) != 0)
307			exit(1);
308
309			/* get input file name */
310			if (args_info.inputs_num)
311			inputFileName = args_info.inputs[0];
312			else {
313			sprintf(painCave.errMsg, "No input .md file name was specified "
314			"on the command line");
315			painCave.isFatal = 1;
316			cmdline_parser_print_help();
317			simError();
318			}
319
320			if (args_info.shells_given) {
321			nShells = args_info.shells_arg;
322			if( nShells < 0 ) {
323			sprintf(painCave.errMsg, "icosahedralBuilder: The number of shells\n"
324			"\tmust be greater than or equal to zero.");
325			painCave.isFatal = 1;
326			cmdline_parser_print_help();
327			simError();
328			}
329			} else {
330			sprintf(painCave.errMsg, "icosahedralBuilder: The number of shells\n"
331			"\tis required to build a Mackay Icosahedron.");
332			painCave.isFatal = 1;
333			cmdline_parser_print_help();
334			simError();
335			}
336
337			if (args_info.latticeConstant_given) {
338			latticeConstant = args_info.latticeConstant_arg;
339			} else {
340
341			int count=0;
342			for( int i = 0; i <= nShells; i++ ) count += np( i );
343			sprintf(painCave.errMsg, "icosahedralBuilder: No lattice constant\n"
344			"\tgiven. Total number of atoms in a Mackay Icosahedron with\n"
345			"\t%d shells is %d.", nShells, count);
346			painCave.isFatal = 1;
347	gezelter	1861	painCave.severity = OPENMD_INFO;
348	gezelter	1860	cmdline_parser_print_help();
349			simError();
350			}
351
352			/* parse md file and set up the system */
353			SimCreator oldCreator;
354			SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
355
356			Globals* simParams = oldInfo->getSimParams();
357
358			//generate the coordinates
359			for( int i = 0; i <= nShells; i++ ) ih( i );
360
361			outputFileName = args_info.output_arg;
362
363			//creat new .md file on fly which corrects the number of molecule
364
365			createMdFile(inputFileName, outputFileName, Points.size());
366
367			if (oldInfo != NULL)
368			delete oldInfo;
369
370			SimCreator newCreator;
371			SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
372
373			// Place molecules
374			Molecule* mol;
375			SimInfo::MoleculeIterator mi;
376			mol = NewInfo->beginMolecule(mi);
377
378			int l = 0;
379
380			locator = new MoLocator(NewInfo->getMoleculeStamp(0),
381			NewInfo->getForceField());
382
383			Vector3d boxMax;
384			Vector3d boxMin;
385
386			for (int n = 0; n < Points.size(); n++) {
387			mol = NewInfo->getMoleculeByGlobalIndex(l);
388			Vector3d location = Points[n] * latticeConstant;
389			Vector3d orientation = Vector3d(0, 0, 1.0);
390
391			if (n == 0) {
392			boxMax = location;
393			boxMin = location;
394			} else {
395			for (int i = 0; i < 3; i++) {
396			boxMax[i] = max(boxMax[i], location[i]);
397			boxMin[i] = min(boxMin[i], location[i]);
398			}
399			}
400
401			locator->placeMol(location, orientation, mol);
402			l++;
403			}
404
405			Mat3x3d boundingBox;
406			boundingBox(0,0) = 10.0*(boxMax[0] - boxMin[0]);
407			boundingBox(1,1) = 10.0*(boxMax[1] - boxMin[1]);
408			boundingBox(2,2) = 10.0*(boxMax[2] - boxMin[2]);
409
410			//set Hmat
411			NewInfo->getSnapshotManager()->getCurrentSnapshot()->setHmat( boundingBox );
412
413			//create dumpwriter and write out the coordinates
414			writer = new DumpWriter(NewInfo, outputFileName);
415
416			if (writer == NULL) {
417			sprintf(painCave.errMsg, "Error in creating dumpwriter object ");
418			painCave.isFatal = 1;
419			simError();
420			}
421
422			writer->writeDump();
423
424			// deleting the writer will put the closing at the end of the dump file
425
426			delete writer;
427
428			// cleanup a by calling sim error.....
429			sprintf(painCave.errMsg, "A new OpenMD file called \"%s\" has been "
430			"generated.\n", outputFileName.c_str());
431			painCave.isFatal = 0;
432			painCave.severity = OPENMD_INFO;
433			simError();
434			return 0;
435			}