1 |
#include "latticeBuilder.hpp" |
2 |
#include <cmath> |
3 |
#include <cstdlib> |
4 |
#include "simError.h" |
5 |
|
6 |
|
7 |
Lattice::Lattice(int latticeType, double latticeParameter){ |
8 |
int hasError; |
9 |
|
10 |
latticePosX = NULL; |
11 |
latticePosY = NULL; |
12 |
latticePosZ = NULL; |
13 |
|
14 |
startX = 0.0; |
15 |
startY = 0.0; |
16 |
startZ = 0.0; |
17 |
|
18 |
switch(latticeType){ |
19 |
|
20 |
case FCC_LATTICE_TYPE: |
21 |
hasError = createFccLattice(latticeParameter); |
22 |
break; |
23 |
|
24 |
case BCC_LATTICE_TYPE: |
25 |
hasError = createBccLattice(latticeParameter); |
26 |
break; |
27 |
|
28 |
case HCP_LATTICE_TYPE: |
29 |
hasError = createHcpLattice(latticeParameter); |
30 |
break; |
31 |
|
32 |
case HCPWATER_LATTICE_TYPE: |
33 |
hasError = createHcpWaterLattice(latticeParameter); |
34 |
break; |
35 |
|
36 |
case ORTHORHOMBIC_LATTICE_TYPE: |
37 |
hasError = createFccLattice(latticeParameter); |
38 |
break; |
39 |
|
40 |
|
41 |
default: |
42 |
//simerror here..... |
43 |
; |
44 |
} |
45 |
|
46 |
latticePosX = new double[nCellSites]; |
47 |
latticePosY = new double[nCellSites]; |
48 |
latticePosZ = new double[nCellSites]; |
49 |
|
50 |
return; |
51 |
} |
52 |
|
53 |
Lattice::Lattice(int latticeType, double latticeParameters[3]){ |
54 |
int hasError; |
55 |
|
56 |
latticePosX = NULL; |
57 |
latticePosY = NULL; |
58 |
latticePosZ = NULL; |
59 |
|
60 |
startX = 0.0; |
61 |
startY = 0.0; |
62 |
startZ = 0.0; |
63 |
|
64 |
switch(latticeType){ |
65 |
|
66 |
case ORTHORHOMBIC_LATTICE_TYPE: |
67 |
hasError = createOrthorhombicLattice(latticeParameters); |
68 |
break; |
69 |
|
70 |
|
71 |
default: |
72 |
//simerror here..... |
73 |
; |
74 |
} |
75 |
|
76 |
latticePosX = new double[nCellSites]; |
77 |
latticePosY = new double[nCellSites]; |
78 |
latticePosZ = new double[nCellSites]; |
79 |
|
80 |
return; |
81 |
} |
82 |
|
83 |
|
84 |
Lattice::~Lattice(void){ |
85 |
if (latticePosX != NULL) delete[] latticePosX; |
86 |
if (latticePosY != NULL) delete[] latticePosY; |
87 |
if (latticePosZ != NULL) delete[] latticePosZ; |
88 |
} |
89 |
|
90 |
|
91 |
int Lattice::createFccLattice(double latticeParameter){ |
92 |
double cell2; |
93 |
double rroot3; |
94 |
|
95 |
cellLengthX = latticeParameter; |
96 |
cellLengthY = latticeParameter; |
97 |
cellLengthZ = latticeParameter; |
98 |
|
99 |
cell2 = 0.5 * latticeParameter; |
100 |
rroot3 = 1.0 / sqrt(3.0); |
101 |
|
102 |
nCellSites = 4; |
103 |
// create new unit cells |
104 |
thisUnitCell.sx = new double[nCellSites]; |
105 |
thisUnitCell.sy = new double[nCellSites]; |
106 |
thisUnitCell.sz = new double[nCellSites]; |
107 |
thisUnitCell.s_ex = new double[nCellSites]; |
108 |
thisUnitCell.s_ey = new double[nCellSites]; |
109 |
thisUnitCell.s_ez = new double[nCellSites]; |
110 |
|
111 |
// add members to each unit cell |
112 |
// Molecule 1 |
113 |
thisUnitCell.sx[0] = 0.0; |
114 |
thisUnitCell.sy[0] = 0.0; |
115 |
thisUnitCell.sz[0] = 0.0; |
116 |
thisUnitCell.s_ex[0] = rroot3; |
117 |
thisUnitCell.s_ey[0] = rroot3; |
118 |
thisUnitCell.s_ez[0] = rroot3; |
119 |
|
120 |
// Molecule 2 |
121 |
thisUnitCell.sx[1] = 0.0; |
122 |
thisUnitCell.sy[1] = cell2; |
123 |
thisUnitCell.sz[1] = cell2; |
124 |
thisUnitCell.s_ex[1] = -rroot3; |
125 |
thisUnitCell.s_ey[1] = rroot3; |
126 |
thisUnitCell.s_ez[1] = -rroot3; |
127 |
|
128 |
// Molecule 3 |
129 |
thisUnitCell.sx[2] = cell2; |
130 |
thisUnitCell.sy[2] = cell2; |
131 |
thisUnitCell.sz[2] = 0.0; |
132 |
thisUnitCell.s_ex[2] = rroot3; |
133 |
thisUnitCell.s_ey[2] = -rroot3; |
134 |
thisUnitCell.s_ez[2] = -rroot3; |
135 |
|
136 |
// Molecule 4 |
137 |
thisUnitCell.sx[3] = cell2; |
138 |
thisUnitCell.sy[3] = 0.0; |
139 |
thisUnitCell.sz[3] = cell2; |
140 |
thisUnitCell.s_ex[3] = -rroot3; |
141 |
thisUnitCell.s_ey[3] = -rroot3; |
142 |
thisUnitCell.s_ez[3] = rroot3; |
143 |
|
144 |
return 0; |
145 |
} |
146 |
|
147 |
|
148 |
int Lattice::createOrthorhombicLattice(double latticeParameters[3]){ |
149 |
double cellx2, celly2, cellz2; |
150 |
double rroot3; |
151 |
|
152 |
cellLengthX = latticeParameters[0]; |
153 |
cellLengthY = latticeParameters[1]; |
154 |
cellLengthZ = latticeParameters[2]; |
155 |
|
156 |
cellx2 = 0.5 * latticeParameters[0]; |
157 |
celly2 = 0.5 * latticeParameters[1]; |
158 |
cellz2 = 0.5 * latticeParameters[2]; |
159 |
rroot3 = 1.0 / sqrt(3.0); |
160 |
|
161 |
nCellSites = 4; |
162 |
// create new unit cells |
163 |
thisUnitCell.sx = new double[nCellSites]; |
164 |
thisUnitCell.sy = new double[nCellSites]; |
165 |
thisUnitCell.sz = new double[nCellSites]; |
166 |
thisUnitCell.s_ex = new double[nCellSites]; |
167 |
thisUnitCell.s_ey = new double[nCellSites]; |
168 |
thisUnitCell.s_ez = new double[nCellSites]; |
169 |
|
170 |
// add members to each unit cell |
171 |
// Molecule 1 |
172 |
thisUnitCell.sx[0] = 0.0; |
173 |
thisUnitCell.sy[0] = 0.0; |
174 |
thisUnitCell.sz[0] = 0.0; |
175 |
thisUnitCell.s_ex[0] = rroot3; |
176 |
thisUnitCell.s_ey[0] = rroot3; |
177 |
thisUnitCell.s_ez[0] = rroot3; |
178 |
|
179 |
// Molecule 2 |
180 |
thisUnitCell.sx[1] = 0.0; |
181 |
thisUnitCell.sy[1] = celly2; |
182 |
thisUnitCell.sz[1] = cellz2; |
183 |
thisUnitCell.s_ex[1] = -rroot3; |
184 |
thisUnitCell.s_ey[1] = rroot3; |
185 |
thisUnitCell.s_ez[1] = -rroot3; |
186 |
|
187 |
// Molecule 3 |
188 |
thisUnitCell.sx[2] = cellx2; |
189 |
thisUnitCell.sy[2] = celly2; |
190 |
thisUnitCell.sz[2] = 0.0; |
191 |
thisUnitCell.s_ex[2] = rroot3; |
192 |
thisUnitCell.s_ey[2] = -rroot3; |
193 |
thisUnitCell.s_ez[2] = -rroot3; |
194 |
|
195 |
// Molecule 4 |
196 |
thisUnitCell.sx[3] = cellx2; |
197 |
thisUnitCell.sy[3] = 0.0; |
198 |
thisUnitCell.sz[3] = cellz2; |
199 |
thisUnitCell.s_ex[3] = -rroot3; |
200 |
thisUnitCell.s_ey[3] = -rroot3; |
201 |
thisUnitCell.s_ez[3] = rroot3; |
202 |
|
203 |
return 0; |
204 |
} |
205 |
|
206 |
|
207 |
// Body centered cubic lattice |
208 |
int Lattice::createBccLattice(double latticeParameter){ |
209 |
return 0; |
210 |
} |
211 |
|
212 |
|
213 |
// Standard HCP lattice |
214 |
int Lattice::createHcpLattice(double latticeParameter){ |
215 |
return 0; |
216 |
} |
217 |
|
218 |
// HCP contains tetrahedral sites for waters |
219 |
int Lattice::createHcpWaterLattice(double latticeParameter){ |
220 |
|
221 |
return 0; |
222 |
|
223 |
double rroot3; |
224 |
double cell; |
225 |
double cell2; |
226 |
double cell4; |
227 |
double cell16; |
228 |
double cell23; |
229 |
double cell34; |
230 |
double cell09; |
231 |
double cell59; |
232 |
|
233 |
nCellSites = 16; |
234 |
|
235 |
// create new unit cells |
236 |
thisUnitCell.sx = new double[nCellSites]; |
237 |
thisUnitCell.sy = new double[nCellSites]; |
238 |
thisUnitCell.sz = new double[nCellSites]; |
239 |
thisUnitCell.s_ex = new double[nCellSites]; |
240 |
thisUnitCell.s_ey = new double[nCellSites]; |
241 |
thisUnitCell.s_ez = new double[nCellSites]; |
242 |
|
243 |
rroot3 = 1.0 / sqrt(3.0); |
244 |
cell = latticeParameter; |
245 |
cell2 = 0.5 * cell; |
246 |
cell4 = 0.25 * cell; |
247 |
cell16 = cell / 6.0; |
248 |
cell23 = (2.0 * cell) / 3.0; |
249 |
cell34 = 0.75 * cell; |
250 |
cell09 = 0.11238 * cell; |
251 |
cell59 = 0.61238 * cell; |
252 |
|
253 |
// ** build the unit cell ** |
254 |
|
255 |
// molecule 1 |
256 |
|
257 |
thisUnitCell.sx[0] = 0.0; |
258 |
thisUnitCell.sy[0] = 0.0; |
259 |
thisUnitCell.sz[0] = 0.0; |
260 |
|
261 |
// molecule 2 |
262 |
|
263 |
thisUnitCell.sx[1] = cell4; |
264 |
thisUnitCell.sy[1] = cell2; |
265 |
thisUnitCell.sz[1] = 0.0; |
266 |
|
267 |
// molecule 3 |
268 |
|
269 |
thisUnitCell.sx[2] = cell4; |
270 |
thisUnitCell.sy[2] = cell16; |
271 |
thisUnitCell.sz[2] = cell09; |
272 |
|
273 |
// molecule 4 |
274 |
|
275 |
thisUnitCell.sx[3] = 0.0; |
276 |
thisUnitCell.sy[3] = cell23; |
277 |
thisUnitCell.sz[3] = cell09; |
278 |
|
279 |
// molecule 5 |
280 |
|
281 |
thisUnitCell.sx[4] = cell4; |
282 |
thisUnitCell.sy[4] = cell16; |
283 |
thisUnitCell.sz[4] = cell2; |
284 |
|
285 |
// molecule 6 |
286 |
|
287 |
thisUnitCell.sx[5] = 0.0; |
288 |
thisUnitCell.sy[5] = cell23; |
289 |
thisUnitCell.sz[5] = cell2; |
290 |
|
291 |
// molecule 7 |
292 |
|
293 |
thisUnitCell.sx[6] = 0.0; |
294 |
thisUnitCell.sy[6] = 0.0; |
295 |
thisUnitCell.sz[6] = cell59; |
296 |
|
297 |
|
298 |
// molecule 8 |
299 |
|
300 |
thisUnitCell.sx[7] = cell4; |
301 |
thisUnitCell.sy[7] = cell2; |
302 |
thisUnitCell.sz[7] = cell59; |
303 |
|
304 |
// molecule 9 |
305 |
|
306 |
thisUnitCell.sz[8] = cell2; |
307 |
thisUnitCell.sy[8] = 0.0; |
308 |
thisUnitCell.sz[8] = 0.0; |
309 |
|
310 |
// molecule 10 |
311 |
|
312 |
thisUnitCell.sz[9] = cell34; |
313 |
thisUnitCell.sy[9] = cell2; |
314 |
thisUnitCell.sz[9] = 0.0; |
315 |
|
316 |
// molecule 11 |
317 |
|
318 |
thisUnitCell.sz[10] = cell34; |
319 |
thisUnitCell.sy[10] = cell16; |
320 |
thisUnitCell.sz[10] = cell09; |
321 |
|
322 |
// molecule 12 |
323 |
|
324 |
thisUnitCell.sz[11] = cell2; |
325 |
thisUnitCell.sy[11] = cell23; |
326 |
thisUnitCell.sz[11] = cell09; |
327 |
|
328 |
// molecule 13 |
329 |
|
330 |
thisUnitCell.sz[12] = cell34; |
331 |
thisUnitCell.sy[12] = cell16; |
332 |
thisUnitCell.sz[12] = cell2; |
333 |
|
334 |
// molecule 14 |
335 |
|
336 |
thisUnitCell.sz[13]= cell2; |
337 |
thisUnitCell.sy[13]= cell23; |
338 |
thisUnitCell.sz[13]= cell2; |
339 |
|
340 |
// molecule 15 |
341 |
|
342 |
thisUnitCell.sz[14] = cell2; |
343 |
thisUnitCell.sy[14] = 0.0; |
344 |
thisUnitCell.sz[14] = cell59; |
345 |
|
346 |
// molecule 16 |
347 |
|
348 |
thisUnitCell.sz[15] = cell34; |
349 |
thisUnitCell.sy[15] = cell2; |
350 |
thisUnitCell.sz[15] = cell59; |
351 |
|
352 |
return 0; |
353 |
} |
354 |
|
355 |
|
356 |
//Returns lattice points when called repeatedly |
357 |
int Lattice::getLatticePoints(double** thePosX, double** thePosY, |
358 |
double** thePosZ, |
359 |
int ix, int iy, int iz){ |
360 |
|
361 |
int iref; |
362 |
|
363 |
for( iref=0;iref < nCellSites;iref++){ |
364 |
|
365 |
latticePosX[iref] = startX + thisUnitCell.sx[iref] + |
366 |
cellLengthX * (double( ix ) - 0.5); |
367 |
latticePosY[iref] = startY + thisUnitCell.sy[iref] + |
368 |
cellLengthY * (double( iy ) - 0.5); |
369 |
latticePosZ[iref] = startZ + thisUnitCell.sz[iref] + |
370 |
cellLengthZ * (double( iz ) - 0.5); |
371 |
|
372 |
} |
373 |
|
374 |
*thePosX = latticePosX; |
375 |
*thePosY = latticePosY; |
376 |
*thePosZ = latticePosZ; |
377 |
|
378 |
return nCellSites; |
379 |
} |
380 |
|