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