| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27 #include "G4CrystalUnitCell.hh"
28
29 #include "G4PhysicalConstants.hh"
30
31 #include <cmath>
32
33 G4CrystalUnitCell::G4CrystalUnitCell(G4double
34 G4double beta, G4double gamma, G4int spacegr
35 : theSize(G4ThreeVector(sizeA, sizeB, sizeC)
36 theAngle(G4ThreeVector(alpha, beta, gamma)
37 theSpaceGroup(spacegroup)
38 {
39 nullVec = G4ThreeVector(0., 0., 0.);
40 theUnitBasis[0] = CLHEP::HepXHat;
41 theUnitBasis[1] = CLHEP::HepYHat;
42 theUnitBasis[2] = CLHEP::HepZHat;
43
44 theRecUnitBasis[0] = CLHEP::HepXHat;
45 theRecUnitBasis[1] = CLHEP::HepYHat;
46 theRecUnitBasis[2] = CLHEP::HepZHat;
47
48 cosa = std::cos(alpha), cosb = std::cos(beta
49 sina = std::sin(alpha), sinb = std::sin(beta
50
51 cosar = (cosb * cosg - cosa) / (sinb * sing)
52 cosbr = (cosa * cosg - cosb) / (sina * sing)
53 cosgr = (cosa * cosb - cosg) / (sina * sinb)
54
55 theVolume = ComputeCellVolume();
56 theRecVolume = 1. / theVolume;
57
58 theRecSize[0] = sizeB * sizeC * sina / theVo
59 theRecSize[1] = sizeC * sizeA * sinb / theVo
60 theRecSize[2] = sizeA * sizeB * sing / theVo
61
62 theRecAngle[0] = std::acos(cosar);
63 theRecAngle[1] = std::acos(cosbr);
64 theRecAngle[2] = std::acos(cosgr);
65
66 G4double x3, y3, z3;
67
68 switch (GetLatticeSystem(theSpaceGroup)) {
69 case Amorphous:
70 break;
71 case Cubic: // Cubic, C44 set
72 break;
73 case Tetragonal:
74 break;
75 case Orthorhombic:
76 break;
77 case Rhombohedral:
78 theUnitBasis[1].rotateZ(gamma - CLHEP::h
79 // Z' axis computed by hand to get both
80 // X'.Z' = cos(alpha), Y'.Z' = cos(beta)
81 x3 = cosa, y3 = (cosb - cosa * cosg) / s
82 theUnitBasis[2] = G4ThreeVector(x3, y3,
83 break;
84 case Monoclinic:
85 theUnitBasis[2].rotateX(beta - CLHEP::ha
86 break;
87 case Triclinic:
88 theUnitBasis[1].rotateZ(gamma - CLHEP::h
89 // Z' axis computed by hand to get both
90 // X'.Z' = cos(alpha), Y'.Z' = cos(beta)
91 x3 = cosa, y3 = (cosb - cosa * cosg) / s
92 theUnitBasis[2] = G4ThreeVector(x3, y3,
93 break;
94 case Hexagonal: // Tetragonal, C16=0
95 theUnitBasis[1].rotateZ(30. * CLHEP::deg
96 break;
97 default:
98 break;
99 }
100
101 for (auto i : {0, 1, 2}) {
102 theBasis[i] = theUnitBasis[i] * theSize[i]
103 theRecBasis[i] = theRecUnitBasis[i] * theR
104 }
105 }
106
107 //....oooOO0OOooo........oooOO0OOooo........oo
108
109 theLatticeSystemType G4CrystalUnitCell::GetLat
110 {
111 if (aGroup >= 1 && aGroup <= 2) {
112 return Triclinic;
113 }
114 if (aGroup >= 3 && aGroup <= 15) {
115 return Monoclinic;
116 }
117 if (aGroup >= 16 && aGroup <= 74) {
118 return Orthorhombic;
119 }
120 if (aGroup >= 75 && aGroup <= 142) {
121 return Tetragonal;
122 }
123 if (aGroup == 146 || aGroup == 148 || aGroup
124 aGroup == 166 || aGroup == 167)
125 {
126 return Rhombohedral;
127 }
128 if (aGroup >= 143 && aGroup <= 167) {
129 return Hexagonal;
130 }
131 if (aGroup >= 168 && aGroup <= 194) {
132 return Hexagonal;
133 }
134 if (aGroup >= 195 && aGroup <= 230) {
135 return Cubic;
136 }
137
138 return Amorphous;
139 }
140
141 //....oooOO0OOooo........oooOO0OOooo........oo
142
143 const G4ThreeVector& G4CrystalUnitCell::GetUni
144 {
145 return (idx >= 0 && idx < 3 ? theUnitBasis[i
146 }
147
148 //....oooOO0OOooo........oooOO0OOooo........oo
149
150 const G4ThreeVector& G4CrystalUnitCell::GetBas
151 {
152 return (idx >= 0 && idx < 3 ? theBasis[idx]
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oo
156
157 const G4ThreeVector& G4CrystalUnitCell::GetRec
158 {
159 return (idx >= 0 && idx < 3 ? theRecUnitBasi
160 }
161
162 //....oooOO0OOooo........oooOO0OOooo........oo
163
164 const G4ThreeVector& G4CrystalUnitCell::GetRec
165 {
166 return (idx >= 0 && idx < 3 ? theRecBasis[id
167 }
168
169 //....oooOO0OOooo........oooOO0OOooo........oo
170
171 G4ThreeVector G4CrystalUnitCell::GetUnitBasisT
172 {
173 // Z' axis computed by hand to get both open
174 // X'.Z' = cos(alpha), Y'.Z' = cos(beta), so
175 G4double x3 = cosa, y3 = (cosb - cosa * cosg
176 return G4ThreeVector(x3, y3, z3).unit();
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181 G4bool G4CrystalUnitCell::FillAtomicUnitPos(G4
182 {
183 // Just for testing the infrastructure
184 G4ThreeVector aaa = pos;
185 vecout.push_back(aaa);
186 vecout.emplace_back(2., 5., 3.);
187 return true;
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oo
191
192 G4bool G4CrystalUnitCell::FillAtomicPos(G4Thre
193 {
194 FillAtomicUnitPos(posin, vecout);
195 for (auto& vec : vecout) {
196 vec.setX(vec.x() * theSize[0]);
197 vec.setY(vec.y() * theSize[1]);
198 vec.setZ(vec.z() * theSize[2]);
199 }
200 return true;
201 }
202
203 //....oooOO0OOooo........oooOO0OOooo........oo
204
205 G4bool G4CrystalUnitCell::FillElReduced(G4doub
206 {
207 switch (GetLatticeSystem()) {
208 case Amorphous:
209 return FillAmorphous(Cij);
210 break;
211 case Cubic: // Cubic, C44 set
212 return FillCubic(Cij);
213 break;
214 case Tetragonal:
215 return FillTetragonal(Cij);
216 break;
217 case Orthorhombic:
218 return FillOrthorhombic(Cij);
219 break;
220 case Rhombohedral:
221 return FillRhombohedral(Cij);
222 break;
223 case Monoclinic:
224 return FillMonoclinic(Cij);
225 break;
226 case Triclinic:
227 return FillTriclinic(Cij);
228 break;
229 case Hexagonal: // Tetragonal, C16=0
230 return FillHexagonal(Cij);
231 break;
232 default:
233 break;
234 }
235
236 return false;
237 }
238
239 //....oooOO0OOooo........oooOO0OOooo........oo
240
241 G4bool G4CrystalUnitCell::FillAmorphous(G4doub
242 {
243 Cij[3][3] = 0.5 * (Cij[0][0] - Cij[0][1]);
244 return true;
245 }
246
247 //....oooOO0OOooo........oooOO0OOooo........oo
248
249 G4bool G4CrystalUnitCell::FillCubic(G4double C
250 {
251 G4double C11 = Cij[0][0], C12 = Cij[0][1], C
252
253 for (size_t i = 0; i < 6; i++) {
254 for (size_t j = i; j < 6; j++) {
255 if (i < 3 && j < 3) {
256 Cij[i][j] = (i == j) ? C11 : C12;
257 }
258 else if (i == j && i >= 3) {
259 Cij[i][i] = C44;
260 }
261 else {
262 Cij[i][j] = 0.;
263 }
264 }
265 }
266
267 ReflectElReduced(Cij);
268
269 return (C11 != 0. && C12 != 0. && C44 != 0.)
270 }
271
272 //....oooOO0OOooo........oooOO0OOooo........oo
273
274 G4bool G4CrystalUnitCell::FillTetragonal(G4dou
275 {
276 G4double C11 = Cij[0][0], C12 = Cij[0][1], C
277 G4double C33 = Cij[2][2], C44 = Cij[3][3], C
278
279 Cij[1][1] = C11; // Copy small number of in
280 Cij[1][2] = C13;
281 Cij[1][5] = -C16;
282 Cij[4][4] = C44;
283
284 ReflectElReduced(Cij);
285
286 // NOTE: Do not test for C16 != 0., to allo
287 return (C11 != 0. && C12 != 0. && C13 != 0.
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oo
291
292 G4bool G4CrystalUnitCell::FillOrthorhombic(G4d
293 {
294 // No degenerate elements; just check for al
295 ReflectElReduced(Cij);
296
297 G4bool good = true;
298 for (size_t i = 0; i < 6; i++) {
299 for (size_t j = i + 1; j < 3; j++) {
300 good &= (Cij[i][j] != 0);
301 }
302 }
303
304 return good;
305 }
306
307 //....oooOO0OOooo........oooOO0OOooo........oo
308
309 G4bool G4CrystalUnitCell::FillRhombohedral(G4d
310 {
311 G4double C11 = Cij[0][0], C12 = Cij[0][1], C
312 G4double C15 = Cij[0][4], C33 = Cij[2][2], C
313
314 Cij[1][1] = C11; // Copy small number of in
315 Cij[1][2] = C13;
316 Cij[1][3] = -C14;
317 Cij[1][4] = -C15;
318 Cij[3][5] = -C15;
319 Cij[4][4] = C44;
320 Cij[4][5] = C14;
321
322 // NOTE: C15 may be zero (c.f. rhombohedral
323 return (C11 != 0 && C12 != 0 && C13 != 0 &&
324 }
325
326 //....oooOO0OOooo........oooOO0OOooo........oo
327
328 G4bool G4CrystalUnitCell::FillMonoclinic(G4dou
329 {
330 // The monoclinic matrix has 13 independent
331 // Sanity condition is same as orthorhombic,
332
333 return (FillOrthorhombic(Cij) && Cij[0][5] !
334 Cij[3][4] != 0.);
335 }
336
337 //....oooOO0OOooo........oooOO0OOooo........oo
338
339 G4bool G4CrystalUnitCell::FillTriclinic(G4doub
340 {
341 // The triclinic matrix has the entire upper
342
343 ReflectElReduced(Cij);
344
345 G4bool good = true;
346 for (size_t i = 0; i < 6; i++) {
347 for (size_t j = i; j < 6; j++) {
348 good &= (Cij[i][j] != 0);
349 }
350 }
351
352 return good;
353 }
354
355 //....oooOO0OOooo........oooOO0OOooo........oo
356
357 G4bool G4CrystalUnitCell::FillHexagonal(G4doub
358 {
359 Cij[0][5] = 0.;
360 Cij[4][5] = 0.5 * (Cij[0][0] - Cij[0][1]);
361 return true;
362 }
363
364 //....oooOO0OOooo........oooOO0OOooo........oo
365
366 G4bool G4CrystalUnitCell::ReflectElReduced(G4d
367 {
368 for (size_t i = 1; i < 6; i++) {
369 for (size_t j = i + 1; j < 6; j++) {
370 Cij[j][i] = Cij[i][j];
371 }
372 }
373 return true;
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oo
377
378 G4double G4CrystalUnitCell::ComputeCellVolume(
379 {
380 G4double a = theSize[0], b = theSize[1], c =
381
382 switch (GetLatticeSystem()) {
383 case Amorphous:
384 return 0.;
385 break;
386 case Cubic:
387 return a * a * a;
388 break;
389 case Tetragonal:
390 return a * a * c;
391 break;
392 case Orthorhombic:
393 return a * b * c;
394 break;
395 case Rhombohedral:
396 return a * a * a * std::sqrt(1. - 3. * c
397 break;
398 case Monoclinic:
399 return a * b * c * sinb;
400 break;
401 case Triclinic:
402 return a * b * c *
403 std::sqrt(1. - cosa * cosa - cosb
404 break;
405 case Hexagonal:
406 return std::sqrt(3.0) / 2. * a * a * c;
407 break;
408 default:
409 break;
410 }
411
412 return 0.;
413 }
414
415 //....oooOO0OOooo........oooOO0OOooo........oo
416
417 G4double G4CrystalUnitCell::GetIntSp2(G4int h,
418 {
419 /* Reference:
420 Table 2.4, pag. 65
421
422 @Inbook{Ladd2003,
423 author="Ladd, Mark and Palmer, Rex",
424 title="Lattices and Space-Group Theory",
425 bookTitle="Structure Determination by X-ray
426 year="2003",
427 publisher="Springer US",
428 address="Boston, MA",
429 pages="51--116",
430 isbn="978-1-4615-0101-5",
431 doi="10.1007/978-1-4615-0101-5_2",
432 url="http://dx.doi.org/10.1007/978-1-4615-0
433 }
434 */
435
436 G4double a = theSize[0], b = theSize[1], c =
437 G4double a2 = a * a, b2 = b * b, c2 = c * c;
438 G4double h2 = h * h, k2 = k * k, l2 = l * l;
439
440 G4double cos2a, sin2a, sin2b;
441 G4double R, T;
442
443 switch (GetLatticeSystem()) {
444 case Amorphous:
445 return 0.;
446 break;
447 case Cubic:
448 return a2 / (h2 + k2 + l2);
449 break;
450 case Tetragonal:
451 return 1.0 / ((h2 + k2) / a2 + l2 / c2);
452 break;
453 case Orthorhombic:
454 return 1.0 / (h2 / a2 + k2 / b2 + l2 / c
455 break;
456 case Rhombohedral:
457 cos2a = cosa * cosa;
458 sin2a = sina * sina;
459 T = h2 + k2 + l2 + 2. * (h * k + k * l +
460 R = sin2a / (1. - 3 * cos2a + 2. * cos2a
461 return a * a / (T * R);
462 break;
463 case Monoclinic:
464 sin2b = sinb * sinb;
465 return 1. / (1. / sin2b * (h2 / a2 + l2
466 break;
467 case Triclinic:
468 return 1. / GetRecIntSp2(h, k, l);
469 break;
470 case Hexagonal:
471 return 1. / ((4. * (h2 + k2 + h * k) / (
472 break;
473 default:
474 break;
475 }
476
477 return 0.;
478 }
479
480 //....oooOO0OOooo........oooOO0OOooo........oo
481
482 G4double G4CrystalUnitCell::GetRecIntSp2(G4int
483 {
484 /* Reference:
485 Table 2.4, pag. 65
486
487 @Inbook{Ladd2003,
488 author="Ladd, Mark and Palmer, Rex",
489 title="Lattices and Space-Group Theory",
490 bookTitle="Structure Determination by X-ray
491 year="2003",
492 publisher="Springer US",
493 address="Boston, MA",
494 pages="51--116",
495 isbn="978-1-4615-0101-5",
496 doi="10.1007/978-1-4615-0101-5_2",
497 url="http://dx.doi.org/10.1007/978-1-4615-0
498 }
499 */
500
501 G4double a = theRecSize[0], b = theRecSize[1
502 G4double a2 = a * a, b2 = b * b, c2 = c * c;
503 G4double h2 = h * h, k2 = k * k, l2 = l * l;
504
505 switch (GetLatticeSystem()) {
506 case Amorphous:
507 return 0.;
508 break;
509 case Cubic:
510 return a2 * (h2 + k2 + l2);
511 break;
512 case Tetragonal:
513 return (h2 + k2) * a2 + l2 * c2;
514 break;
515 case Orthorhombic:
516 return h2 * a2 + k2 + b2 + h2 * c2;
517 break;
518 case Rhombohedral:
519 return (h2 + k2 + l2 + 2. * (h * k + k *
520 break;
521 case Monoclinic:
522 return h2 * a2 + k2 * b2 + l2 * c2 + 2.
523 break;
524 case Triclinic:
525 return h2 * a2 + k2 * b2 + l2 * c2 + 2.
526 2. * h * k * a * b * cosgr;
527 break;
528 case Hexagonal:
529 return (h2 + k2 + h * k) * a2 + l2 * c2;
530 break;
531 default:
532 break;
533 }
534
535 return 0.;
536 }
537
538 //....oooOO0OOooo........oooOO0OOooo........oo
539
540 G4double G4CrystalUnitCell::GetIntCosAng(G4int
541 {
542 /* Reference:
543 Table 2.4, pag. 65
544
545 @Inbook{Kelly2012,
546 author="Anthony A. Kelly and Kevin M. Knowl
547 title="Appendix 3 Interplanar Spacings and
548 bookTitle="Crystallography and Crystal Defe
549 year="2012",
550 publisher="John Wiley & Sons, Ltd.",
551 isbn="978-0-470-75014-8",
552 doi="10.1002/9781119961468",
553 url="http://onlinelibrary.wiley.com/book/10
554 }
555 */
556
557 G4double a = theRecSize[0], b = theRecSize[1
558 G4double a2 = a * a, b2 = b * b, c2 = c * c;
559 G4double dsp1dsp2;
560 switch (GetLatticeSystem()) {
561 case Amorphous:
562 return 0.;
563 break;
564 case Cubic:
565 return (h1 * h2 + k1 * k2 + l1 + l2) /
566 (std::sqrt(h1 * h1 + k1 * k1 + l1
567 break;
568 case Tetragonal:
569 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l
570 return 0.;
571 break;
572 case Orthorhombic:
573 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l
574 return dsp1dsp2 * (h1 * h2 * a2 + k1 * k
575 break;
576 case Rhombohedral:
577 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l
578 return dsp1dsp2 *
579 (h1 * h2 * a2 + k1 * k2 * b2 + l1
580 (h1 * l2 + h2 * l1) * a * c * c
581 break;
582 case Monoclinic:
583 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l
584 return dsp1dsp2 *
585 (h1 * h2 * a2 + k1 * k2 * b2 + l1
586 (h1 * l2 + h2 * l1) * a * c * c
587 break;
588 case Triclinic:
589 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l
590 return dsp1dsp2 *
591 (h1 * h2 * a2 + k1 * k2 * b2 + l1
592 (h1 * l2 + h2 * l1) * a * c * c
593 break;
594 case Hexagonal:
595 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l
596 return dsp1dsp2 * ((h1 * h2 + k1 * k2 +
597 break;
598 default:
599 break;
600 }
601
602 return 0.;
603 }
604
605 //....oooOO0OOooo........oooOO0OOooo........oo
606