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