ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/nanoparticleBuilder/icosahedralBuilder.cpp
(Generate patch)

Comparing branches/development/src/applications/nanoparticleBuilder/icosahedralBuilder.cpp (file contents):
Revision 1860 by gezelter, Tue Apr 9 16:14:15 2013 UTC vs.
Revision 1876 by gezelter, Fri May 17 17:10:11 2013 UTC

# Line 40 | Line 40
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 –
43   #include "config.h"
44   #include "icosahedralBuilderCmd.h"
45   #include "utils/MoLocator.hpp"
46   #include "utils/Tuple.hpp"
57 – #include "brains/Register.hpp"
47   #include "brains/SimInfo.hpp"
48   #include "brains/SimCreator.hpp"
49   #include "io/DumpWriter.hpp"
50 < #include "math/Vector3.hpp"
51 < #include "math/SquareMatrix3.hpp"
63 < #include "utils/StringUtils.hpp"
50 > #include "clusters/Icosahedron.hpp"
51 > #include "clusters/Decahedron.hpp"
52  
53   using namespace OpenMD;
54   using namespace std;
55  
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 –
56   void createMdFile(const std::string&oldMdFileName,
57                    const std::string&newMdFileName,
58                    int nMol) {
# Line 310 | Line 93 | int main(int argc, char *argv []) {
93    MoLocator* locator;
94    RealType latticeConstant;
95    int nShells;
313 –
314 –  DumpWriter *writer;
96  
97 <  init();
97 >  DumpWriter *writer;
98    
99    // Parse Command Line Arguments
100    if (cmdline_parser(argc, argv, &args_info) != 0)
# Line 349 | Line 130 | int main(int argc, char *argv []) {
130  
131    if (args_info.latticeConstant_given) {    
132      latticeConstant = args_info.latticeConstant_arg;
133 <  } else  {
353 <    
354 <    int count=0;
355 <    for( int i = 0; i <= nShells; i++ ) count += np( i );
133 >  } else  {  
134      sprintf(painCave.errMsg, "icosahedralBuilder:  No lattice constant\n"
135 <            "\tgiven.  Total number of atoms in a Mackay Icosahedron with\n"
358 <            "\t%d shells is %d.", nShells, count);
135 >            "\tgiven.");
136      painCave.isFatal = 1;
137      cmdline_parser_print_help();
138      simError();
# Line 364 | Line 141 | int main(int argc, char *argv []) {
141    /* parse md file and set up the system */
142    SimCreator oldCreator;
143    SimInfo* oldInfo = oldCreator.createSim(inputFileName, false);
367 –
144    Globals* simParams = oldInfo->getSimParams();
145  
146 <  //generate the coordinates
147 <  for( int i = 0; i <= nShells; i++ ) ih( i );
146 >  vector<Vector3d> Points;
147 >  if (args_info.ico_given){
148 >    Icosahedron* ico = new Icosahedron();
149 >    Points = ico->getPoints(nShells);
150 >  } else if (args_info.deca_given) {
151 >    RegularDecahedron* deca = new RegularDecahedron(nShells);
152 >    Points = deca->getPoints();
153 >  } else if (args_info.ino_given) {
154 >    int columnAtoms = args_info.columnAtoms_arg;    
155 >    InoDecahedron* ino = new InoDecahedron(columnAtoms, nShells);
156 >    Points = ino->getPoints();
157 >  } else if (args_info.marks_given) {
158 >    int columnAtoms = args_info.columnAtoms_arg;    
159 >    int twinAtoms = args_info.twinAtoms_arg;
160 >    Decahedron* marks = new Decahedron(columnAtoms, nShells, twinAtoms);
161 >    Points = marks->getPoints();
162 >  }
163 >  else if (args_info.stone_given) {
164 >    int columnAtoms = args_info.columnAtoms_arg;    
165 >    int twinAtoms = args_info.twinAtoms_arg;
166 >    int truncatedPlanes = args_info.truncatedPlanes_arg;
167 >    CurlingStoneDecahedron* csd = new CurlingStoneDecahedron(columnAtoms, nShells, twinAtoms, truncatedPlanes);
168 >    Points = csd->getPoints();
169 >  }
170 >  
171 >    
172  
173    outputFileName = args_info.output_arg;
174    
# Line 376 | Line 176 | int main(int argc, char *argv []) {
176  
177    createMdFile(inputFileName, outputFileName, Points.size());
178    
179 <  if (oldInfo != NULL)
380 <    delete oldInfo;
179 >  delete oldInfo;
180    
181    SimCreator newCreator;
182    SimInfo* NewInfo = newCreator.createSim(outputFileName, false);
# Line 395 | Line 194 | int main(int argc, char *argv []) {
194    Vector3d boxMax;
195    Vector3d boxMin;
196  
197 <  for (int n = 0; n < Points.size(); n++) {
197 >  for (unsigned int n = 0; n < Points.size(); n++) {
198      mol = NewInfo->getMoleculeByGlobalIndex(l);
199      Vector3d location = Points[n] * latticeConstant;
200      Vector3d orientation = Vector3d(0, 0, 1.0);
201 <
201 >    
202      if (n == 0) {
203        boxMax = location;
204        boxMin = location;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines