| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26
27 #include "G4CrystalUnitCell.hh"
28
29 #include "G4PhysicalConstants.hh"
30
31 #include <cmath>
32
33 G4CrystalUnitCell::G4CrystalUnitCell(G4double sizeA, G4double sizeB, G4double sizeC, G4double alpha,
34 G4double beta, G4double gamma, G4int spacegroup)
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), cosg = std::cos(gamma);
49 sina = std::sin(alpha), sinb = std::sin(beta), sing = std::sin(gamma);
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 / theVolume;
59 theRecSize[1] = sizeC * sizeA * sinb / theVolume;
60 theRecSize[2] = sizeA * sizeB * sing / theVolume;
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::halfpi); // X-Y opening angle
79 // Z' axis computed by hand to get both opening angles right
80 // X'.Z' = cos(alpha), Y'.Z' = cos(beta), solve for Z' components
81 x3 = cosa, y3 = (cosb - cosa * cosg) / sing, z3 = std::sqrt(1. - x3 * x3 - y3 * y3);
82 theUnitBasis[2] = G4ThreeVector(x3, y3, z3).unit();
83 break;
84 case Monoclinic:
85 theUnitBasis[2].rotateX(beta - CLHEP::halfpi); // Z-Y opening angle
86 break;
87 case Triclinic:
88 theUnitBasis[1].rotateZ(gamma - CLHEP::halfpi); // X-Y opening angle
89 // Z' axis computed by hand to get both opening angles right
90 // X'.Z' = cos(alpha), Y'.Z' = cos(beta), solve for Z' components
91 x3 = cosa, y3 = (cosb - cosa * cosg) / sing, z3 = std::sqrt(1. - x3 * x3 - y3 * y3);
92 theUnitBasis[2] = G4ThreeVector(x3, y3, z3).unit();
93 break;
94 case Hexagonal: // Tetragonal, C16=0
95 theUnitBasis[1].rotateZ(30. * CLHEP::deg); // X-Y opening angle
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] * theRecSize[i];
104 }
105 }
106
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108
109 theLatticeSystemType G4CrystalUnitCell::GetLatticeSystem(G4int aGroup)
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 == 155 || aGroup == 160 || aGroup == 161 ||
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........oooOO0OOooo........oooOO0OOooo......
142
143 const G4ThreeVector& G4CrystalUnitCell::GetUnitBasis(G4int idx) const
144 {
145 return (idx >= 0 && idx < 3 ? theUnitBasis[idx] : nullVec);
146 }
147
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149
150 const G4ThreeVector& G4CrystalUnitCell::GetBasis(G4int idx) const
151 {
152 return (idx >= 0 && idx < 3 ? theBasis[idx] : nullVec);
153 }
154
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156
157 const G4ThreeVector& G4CrystalUnitCell::GetRecUnitBasis(G4int idx) const
158 {
159 return (idx >= 0 && idx < 3 ? theRecUnitBasis[idx] : nullVec);
160 }
161
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
163
164 const G4ThreeVector& G4CrystalUnitCell::GetRecBasis(G4int idx) const
165 {
166 return (idx >= 0 && idx < 3 ? theRecBasis[idx] : nullVec);
167 }
168
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170
171 G4ThreeVector G4CrystalUnitCell::GetUnitBasisTrigonal()
172 {
173 // Z' axis computed by hand to get both opening angles right
174 // X'.Z' = cos(alpha), Y'.Z' = cos(beta), solve for Z' components
175 G4double x3 = cosa, y3 = (cosb - cosa * cosg) / sing, z3 = std::sqrt(1. - x3 * x3 - y3 * y3);
176 return G4ThreeVector(x3, y3, z3).unit();
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181 G4bool G4CrystalUnitCell::FillAtomicUnitPos(G4ThreeVector& pos, std::vector<G4ThreeVector>& vecout)
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........oooOO0OOooo........oooOO0OOooo......
191
192 G4bool G4CrystalUnitCell::FillAtomicPos(G4ThreeVector& posin, std::vector<G4ThreeVector>& vecout)
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........oooOO0OOooo........oooOO0OOooo......
204
205 G4bool G4CrystalUnitCell::FillElReduced(G4double Cij[6][6])
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........oooOO0OOooo........oooOO0OOooo......
240
241 G4bool G4CrystalUnitCell::FillAmorphous(G4double Cij[6][6]) const
242 {
243 Cij[3][3] = 0.5 * (Cij[0][0] - Cij[0][1]);
244 return true;
245 }
246
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248
249 G4bool G4CrystalUnitCell::FillCubic(G4double Cij[6][6]) const
250 {
251 G4double C11 = Cij[0][0], C12 = Cij[0][1], C44 = Cij[3][3];
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........oooOO0OOooo........oooOO0OOooo......
273
274 G4bool G4CrystalUnitCell::FillTetragonal(G4double Cij[6][6]) const
275 {
276 G4double C11 = Cij[0][0], C12 = Cij[0][1], C13 = Cij[0][2], C16 = Cij[0][5];
277 G4double C33 = Cij[2][2], C44 = Cij[3][3], C66 = Cij[5][5];
278
279 Cij[1][1] = C11; // Copy small number of individual elements
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 allow calling from Hexagonal
287 return (C11 != 0. && C12 != 0. && C13 != 0. && C33 != 0. && C44 != 0. && C66 != 0.);
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291
292 G4bool G4CrystalUnitCell::FillOrthorhombic(G4double Cij[6][6]) const
293 {
294 // No degenerate elements; just check for all non-zero
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........oooOO0OOooo........oooOO0OOooo......
308
309 G4bool G4CrystalUnitCell::FillRhombohedral(G4double Cij[6][6]) const
310 {
311 G4double C11 = Cij[0][0], C12 = Cij[0][1], C13 = Cij[0][2], C14 = Cij[0][3];
312 G4double C15 = Cij[0][4], C33 = Cij[2][2], C44 = Cij[3][3], C66 = 0.5 * (C11 - C12);
313
314 Cij[1][1] = C11; // Copy small number of individual elements
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(I) vs. (II))
323 return (C11 != 0 && C12 != 0 && C13 != 0 && C14 != 0. && C33 != 0. && C44 != 0. && C66 != 0.);
324 }
325
326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
327
328 G4bool G4CrystalUnitCell::FillMonoclinic(G4double Cij[6][6]) const
329 {
330 // The monoclinic matrix has 13 independent elements with no degeneracies
331 // Sanity condition is same as orthorhombic, plus C45, C(1,2,3)6
332
333 return (FillOrthorhombic(Cij) && Cij[0][5] != 0. && Cij[1][5] != 0. && Cij[2][5] != 0. &&
334 Cij[3][4] != 0.);
335 }
336
337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
338
339 G4bool G4CrystalUnitCell::FillTriclinic(G4double Cij[6][6]) const
340 {
341 // The triclinic matrix has the entire upper half filled (21 elements)
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........oooOO0OOooo........oooOO0OOooo......
356
357 G4bool G4CrystalUnitCell::FillHexagonal(G4double Cij[6][6]) const
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........oooOO0OOooo........oooOO0OOooo......
365
366 G4bool G4CrystalUnitCell::ReflectElReduced(G4double Cij[6][6]) const
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........oooOO0OOooo........oooOO0OOooo......
377
378 G4double G4CrystalUnitCell::ComputeCellVolume()
379 {
380 G4double a = theSize[0], b = theSize[1], c = theSize[2];
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. * cosa * cosa + 2. * cosa * cosa * cosa);
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 * cosb - cosg * cosg * 2. * cosa * cosb * cosg);
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........oooOO0OOooo........oooOO0OOooo....
416
417 G4double G4CrystalUnitCell::GetIntSp2(G4int h, G4int k, G4int l)
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 Crystallography",
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-0101-5_2"
433 }
434 */
435
436 G4double a = theSize[0], b = theSize[1], c = theSize[2];
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 / c2);
455 break;
456 case Rhombohedral:
457 cos2a = cosa * cosa;
458 sin2a = sina * sina;
459 T = h2 + k2 + l2 + 2. * (h * k + k * l + h * l) * ((cos2a - cosa) / sin2a);
460 R = sin2a / (1. - 3 * cos2a + 2. * cos2a * cosa);
461 return a * a / (T * R);
462 break;
463 case Monoclinic:
464 sin2b = sinb * sinb;
465 return 1. / (1. / sin2b * (h2 / a2 + l2 / c2 - 2 * h * l * cosb / (a * c)) + k2 / b2);
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) / (3. * a2)) + l2 / c2);
472 break;
473 default:
474 break;
475 }
476
477 return 0.;
478 }
479
480 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
481
482 G4double G4CrystalUnitCell::GetRecIntSp2(G4int h, G4int k, G4int l)
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 Crystallography",
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-0101-5_2"
498 }
499 */
500
501 G4double a = theRecSize[0], b = theRecSize[1], c = theRecSize[2];
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 * l + h * l) * cosar) * a2;
520 break;
521 case Monoclinic:
522 return h2 * a2 + k2 * b2 + l2 * c2 + 2. * h * l * a * c * cosbr;
523 break;
524 case Triclinic:
525 return h2 * a2 + k2 * b2 + l2 * c2 + 2. * k * l * b * c * cosar + 2. * l * h * c * a * cosbr +
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........oooOO0OOooo........oooOO0OOooo....
539
540 G4double G4CrystalUnitCell::GetIntCosAng(G4int h1, G4int k1, G4int l1, G4int h2, G4int k2, G4int l2)
541 {
542 /* Reference:
543 Table 2.4, pag. 65
544
545 @Inbook{Kelly2012,
546 author="Anthony A. Kelly and Kevin M. Knowles",
547 title="Appendix 3 Interplanar Spacings and Interplanar Angles",
548 bookTitle="Crystallography and Crystal Defects, 2nd Edition",
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.1002/9781119961468"
554 }
555 */
556
557 G4double a = theRecSize[0], b = theRecSize[1], c = theRecSize[2];
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 * l1) * std::sqrt(h2 * h2 + k2 * k2 + l2 * l2));
567 break;
568 case Tetragonal:
569 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l1) * GetIntSp2(h2, k2, l2));
570 return 0.;
571 break;
572 case Orthorhombic:
573 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l1) * GetIntSp2(h2, k2, l2));
574 return dsp1dsp2 * (h1 * h2 * a2 + k1 * k2 * a2 + l1 * l2 * c2);
575 break;
576 case Rhombohedral:
577 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l1) * GetIntSp2(h2, k2, l2));
578 return dsp1dsp2 *
579 (h1 * h2 * a2 + k1 * k2 * b2 + l1 * l2 * c2 + (k1 * l2 + k2 * l1) * b * c * cosar +
580 (h1 * l2 + h2 * l1) * a * c * cosbr + (h1 * k2 + h2 * k1) * a * b * cosgr);
581 break;
582 case Monoclinic:
583 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l1) * GetIntSp2(h2, k2, l2));
584 return dsp1dsp2 *
585 (h1 * h2 * a2 + k1 * k2 * b2 + l1 * l2 * c2 + (k1 * l2 + k2 * l1) * b * c * cosar +
586 (h1 * l2 + h2 * l1) * a * c * cosbr + (h1 * k2 + h2 * k1) * a * b * cosgr);
587 break;
588 case Triclinic:
589 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l1) * GetIntSp2(h2, k2, l2));
590 return dsp1dsp2 *
591 (h1 * h2 * a2 + k1 * k2 * b2 + l1 * l2 * c2 + (k1 * l2 + k2 * l1) * b * c * cosar +
592 (h1 * l2 + h2 * l1) * a * c * cosbr + (h1 * k2 + h2 * k1) * a * b * cosgr);
593 break;
594 case Hexagonal:
595 dsp1dsp2 = std::sqrt(GetIntSp2(h1, k1, l1) * GetIntSp2(h2, k2, l2));
596 return dsp1dsp2 * ((h1 * h2 + k1 * k2 + 0.5 * (h1 * k2 + k1 * h2)) * a2 + l1 * l2 * c2);
597 break;
598 default:
599 break;
600 }
601
602 return 0.;
603 }
604
605 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
606