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 <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <string>
#include <map>
#include <fstream>
#include <algorithm>

#include "config.h"
#include "icosahedralBuilderCmd.h"
#include "utils/MoLocator.hpp"
#include "utils/Tuple.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 OpenMD;
using namespace std;

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

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