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