47 |
|
#include "brains/SimInfo.hpp" |
48 |
|
#include "brains/SimCreator.hpp" |
49 |
|
#include "io/DumpWriter.hpp" |
50 |
< |
#include "math/SquareMatrix3.hpp" |
50 |
> |
#include "clusters/Icosahedron.hpp" |
51 |
> |
#include "clusters/Decahedron.hpp" |
52 |
|
|
53 |
|
using namespace OpenMD; |
54 |
|
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 |
– |
vector<pair<int, int> > Edges; |
63 |
– |
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 |
– |
const RealType HGR = ( sqrt(5.0) + 1.0 ) / 4.0; // half of the golden ratio |
103 |
– |
|
104 |
– |
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 |
– |
|
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 |
– |
} |
55 |
|
|
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 |
– |
|
56 |
|
void createMdFile(const std::string&oldMdFileName, |
57 |
|
const std::string&newMdFileName, |
58 |
|
int nMol) { |
93 |
|
MoLocator* locator; |
94 |
|
RealType latticeConstant; |
95 |
|
int nShells; |
300 |
– |
|
301 |
– |
DumpWriter *writer; |
96 |
|
|
97 |
< |
init(); |
97 |
> |
DumpWriter *writer; |
98 |
|
|
99 |
|
// Parse Command Line Arguments |
100 |
|
if (cmdline_parser(argc, argv, &args_info) != 0) |
130 |
|
|
131 |
|
if (args_info.latticeConstant_given) { |
132 |
|
latticeConstant = args_info.latticeConstant_arg; |
133 |
< |
} else { |
340 |
< |
|
341 |
< |
int count=0; |
342 |
< |
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" |
345 |
< |
"\t%d shells is %d.", nShells, count); |
135 |
> |
"\tgiven."); |
136 |
|
painCave.isFatal = 1; |
347 |
– |
painCave.severity = OPENMD_INFO; |
137 |
|
cmdline_parser_print_help(); |
138 |
|
simError(); |
139 |
|
} |
141 |
|
/* parse md file and set up the system */ |
142 |
|
SimCreator oldCreator; |
143 |
|
SimInfo* oldInfo = oldCreator.createSim(inputFileName, false); |
355 |
– |
|
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 |
|
|