ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CrystalUnitCell.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4CrystalUnitCell.cc
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 
28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
29 
30 // 21-04-16, created by E.Bagli
31 
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "G4CrystalUnitCell.hh"
35 #include "G4PhysicalConstants.hh"
36 #include <cmath>
37 
38 G4CrystalUnitCell::G4CrystalUnitCell(G4double sizeA,
39  G4double sizeB,
40  G4double sizeC,
41  G4double alpha,
42  G4double beta,
43  G4double gamma,
44  G4int spacegroup):
45 theSpaceGroup(spacegroup),
46 theSize(G4ThreeVector(sizeA,sizeB,sizeC)),
47 theAngle(G4ThreeVector(alpha,beta,gamma))
48 {
49 
50  nullVec = G4ThreeVector(0.,0.,0.);
51  theUnitBasis[0] = CLHEP::HepXHat;
52  theUnitBasis[1] = CLHEP::HepYHat;
53  theUnitBasis[2] = CLHEP::HepZHat;
54 
55  theRecUnitBasis[0] = CLHEP::HepXHat;
56  theRecUnitBasis[1] = CLHEP::HepYHat;
57  theRecUnitBasis[2] = CLHEP::HepZHat;
58 
59  cosa=std::cos(alpha), cosb=std::cos(beta), cosg=std::cos(gamma);
60  sina=std::sin(alpha), sinb=std::sin(beta), sing=std::sin(gamma);
61 
62  cosar = (cosb*cosg-cosa)/(sinb*sing);
63  cosbr = (cosa*cosg-cosb)/(sina*sing);
64  cosgr = (cosa*cosb-cosg)/(sina*sinb);
65 
66  theVolume = ComputeCellVolume();
67  theRecVolume = 1. / theVolume;
68 
69  theRecSize[0] = sizeB * sizeC * sina / theVolume;
70  theRecSize[1] = sizeC * sizeA * sinb / theVolume;
71  theRecSize[2] = sizeA * sizeB * sing / theVolume;
72 
73  theRecAngle[0] = std::acos(cosar);
74  theRecAngle[1] = std::acos(cosbr);
75  theRecAngle[2] = std::acos(cosgr);
76 
77  G4double x3,y3,z3;
78 
79  switch (GetLatticeSystem(theSpaceGroup)) {
80  case Amorphous:
81  break;
82  case Cubic: // Cubic, C44 set
83  break;
84  case Tetragonal:
85  break;
86  case Orthorhombic:
87  break;
88  case Rhombohedral:
89  theUnitBasis[1].rotateZ(gamma-CLHEP::halfpi); // X-Y opening angle
90  // Z' axis computed by hand to get both opening angles right
91  // X'.Z' = cos(alpha), Y'.Z' = cos(beta), solve for Z' components
92  x3=cosa, y3=(cosb-cosa*cosg)/sing, z3=std::sqrt(1.-x3*x3-y3*y3);
93  theUnitBasis[2] = G4ThreeVector(x3, y3, z3).unit();
94  break;
95  case Monoclinic:
96  theUnitBasis[2].rotateX(beta-CLHEP::halfpi); // Z-Y opening angle
97  break;
98  case Triclinic:
99  theUnitBasis[1].rotateZ(gamma-CLHEP::halfpi); // X-Y opening angle
100  // Z' axis computed by hand to get both opening angles right
101  // X'.Z' = cos(alpha), Y'.Z' = cos(beta), solve for Z' components
102  x3=cosa, y3=(cosb-cosa*cosg)/sing, z3=std::sqrt(1.-x3*x3-y3*y3);
103  theUnitBasis[2] = G4ThreeVector(x3, y3, z3).unit();
104  break;
105  case Hexagonal: // Tetragonal, C16=0
106  theUnitBasis[1].rotateZ(30.*CLHEP::deg); // X-Y opening angle
107  break;
108  default:
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 
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 G4CrystalUnitCell::~G4CrystalUnitCell(){;}
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
148 
149 theLatticeSystemType G4CrystalUnitCell::GetLatticeSystem(G4int aGroup){
150 
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;
164 }
165 
166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167 /*
168 theBravaisLatticeType G4CrystalUnitCell::GetBravaisLattice(G4int aGroup){
169  ;
170 }
171 */
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173 
174 const G4ThreeVector& G4CrystalUnitCell::GetUnitBasis(G4int idx) const {
175  return (idx>=0 && idx<3 ? theUnitBasis[idx] : nullVec);
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179 
180 const G4ThreeVector& G4CrystalUnitCell::GetBasis(G4int idx) const {
181  return (idx>=0 && idx<3 ? theBasis[idx] : nullVec);
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 
186 const G4ThreeVector& G4CrystalUnitCell::GetRecUnitBasis(G4int idx) const {
187  return (idx>=0 && idx<3 ? theRecUnitBasis[idx] : nullVec);
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
192 const G4ThreeVector& G4CrystalUnitCell::GetRecBasis(G4int idx) const {
193  return (idx>=0 && idx<3 ? theRecBasis[idx] : nullVec);
194 }
195 
196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197 
198 G4ThreeVector G4CrystalUnitCell::GetUnitBasisTrigonal(){
199  // Z' axis computed by hand to get both opening angles right
200  // X'.Z' = cos(alpha), Y'.Z' = cos(beta), solve for Z' components
201  G4double x3=cosa, y3=(cosb-cosa*cosg)/sing, z3=std::sqrt(1.-x3*x3-y3*y3);
202  return G4ThreeVector(x3, y3, z3).unit();
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
206 
207 G4bool G4CrystalUnitCell::FillAtomicUnitPos(G4ThreeVector& pos, std::vector<G4ThreeVector>& vecout){
208  // Just for testing the infrastructure
209  G4ThreeVector aaa = pos;
210  vecout.push_back(aaa);
211  vecout.push_back(G4ThreeVector(2.,5.,3.));
212  return true;
213 }
214 
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
216 
217 G4bool G4CrystalUnitCell::FillAtomicPos(G4ThreeVector& posin, std::vector<G4ThreeVector>& vecout){
218  FillAtomicUnitPos(posin,vecout);
219  for(auto &vec:vecout){
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;
225 }
226 
227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228 
229 G4bool G4CrystalUnitCell::FillElReduced(G4double Cij[6][6]) {
230  switch (GetLatticeSystem()) {
231  case Amorphous:
232  return FillAmorphous(Cij);
233  break;
234  case Cubic: // Cubic, C44 set
235  return FillCubic(Cij);
236  break;
237  case Tetragonal:
238  return FillTetragonal(Cij);
239  break;
240  case Orthorhombic:
241  return FillOrthorhombic(Cij);
242  break;
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;
257  }
258 
259  return false;
260 }
261 
262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
263 
264 G4bool G4CrystalUnitCell::FillAmorphous(G4double Cij[6][6]) const {
265  Cij[3][3] = 0.5*(Cij[0][0]-Cij[0][1]);
266  return true;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 
271 G4bool G4CrystalUnitCell::FillCubic(G4double Cij[6][6]) const {
272  G4double C11=Cij[0][0], C12=Cij[0][1], C44=Cij[3][3];
273 
274  for (size_t i=0; i<6; i++) {
275  for (size_t j=i; j<6; j++) {
276  if (i<3 && j<3) Cij[i][j] = (i==j) ? C11 : C12;
277  else if (i==j && i>=3) Cij[i][i] = C44;
278  else Cij[i][j] = 0.;
279  }
280  }
281 
282  ReflectElReduced(Cij);
283 
284  return (C11!=0. && C12!=0. && C44!=0.);
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 G4bool G4CrystalUnitCell::FillTetragonal(G4double Cij[6][6]) const {
290  G4double C11=Cij[0][0], C12=Cij[0][1], C13=Cij[0][2], C16=Cij[0][5];
291  G4double C33=Cij[2][2], C44=Cij[3][3], C66=Cij[5][5];
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 }
303 
304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
305 
306 G4bool G4CrystalUnitCell::FillOrthorhombic(G4double Cij[6][6]) const {
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;
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 
321 G4bool G4CrystalUnitCell::FillRhombohedral(G4double Cij[6][6]) const {
322  G4double C11=Cij[0][0], C12=Cij[0][1], C13=Cij[0][2], C14=Cij[0][3];
323  G4double C15=Cij[0][4], C33=Cij[2][2], C44=Cij[3][3], C66=0.5*(C11-C12);
324 
325  Cij[1][1] = C11; // Copy small number of individual elements
326  Cij[1][2] = C13;
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.);
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
339 
340 G4bool G4CrystalUnitCell::FillMonoclinic(G4double Cij[6][6]) const {
341  // The monoclinic matrix has 13 independent elements with no degeneracies
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 }
347 
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349 
350 G4bool G4CrystalUnitCell::FillTriclinic(G4double Cij[6][6]) const {
351  // The triclinic matrix has the entire upper half filled (21 elements)
352 
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);
358  }
359 
360  return good;
361 }
362 
363 
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
365 
366 G4bool G4CrystalUnitCell::FillHexagonal(G4double Cij[6][6]) const {
367  Cij[0][5] = 0.;
368  Cij[4][5] = 0.5*(Cij[0][0] - Cij[0][1]);
369  return true;
370 }
371 
372 
373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
374 
375 G4bool G4CrystalUnitCell::ReflectElReduced(G4double Cij[6][6]) const {
376  for (size_t i=1; i<6; i++) {
377  for (size_t j=i+1; j<6; j++) {
378  Cij[j][i] = Cij[i][j];
379  }
380  }
381  return true;
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385 
386 G4double G4CrystalUnitCell::ComputeCellVolume(){
387  G4double a = theSize[0], b = theSize[1], c = theSize[2];
388 
389  switch(GetLatticeSystem())
390  {
391  case Amorphous:
392  return 0.;
393  break;
394  case Cubic:
395  return a * a * a;
396  break;
397  case Tetragonal:
398  return a * a * c;
399  break;
400  case Orthorhombic:
401  return a * b * c;
402  break;
403  case Rhombohedral:
404  return a*a*a*std::sqrt(1.-3.*cosa*cosa+2.*cosa*cosa*cosa);
405  break;
406  case Monoclinic:
407  return a*b*c*sinb;
408  break;
409  case Triclinic:
410  return a*b*c*std::sqrt(1.-cosa*cosa-cosb*cosb-cosg*cosg*2.*cosa*cosb*cosg);
411  break;
412  case Hexagonal:
413  return std::sqrt(3.0)/2.*a*a*c;
414  break;
415  default:
416  break;
417  }
418 
419  return 0.;
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
423 
424 G4double G4CrystalUnitCell::GetIntSp2(G4int h,
425  G4int k,
426  G4int l){
427 
428  /* Reference:
429  Table 2.4, pag. 65
430 
431  @Inbook{Ladd2003,
432  author="Ladd, Mark and Palmer, Rex",
433  title="Lattices and Space-Group Theory",
434  bookTitle="Structure Determination by X-ray Crystallography",
435  year="2003",
436  publisher="Springer US",
437  address="Boston, MA",
438  pages="51--116",
439  isbn="978-1-4615-0101-5",
440  doi="10.1007/978-1-4615-0101-5_2",
441  url="http://dx.doi.org/10.1007/978-1-4615-0101-5_2"
442  }
443  */
444 
445  G4double a = theSize[0], b = theSize[1], c = theSize[2];
446  G4double a2 = a*a, b2 = b*b, c2 = c*c;
447  G4double h2 = h*h, k2 = k*k, l2 = l*l;
448 
449  G4double cos2a,sin2a,sin2b;
450  G4double R,T;
451 
452  switch(GetLatticeSystem())
453  {
454  case Amorphous:
455  return 0.;
456  break;
457  case Cubic:
458  return a2 / ( h2+k2+l2 );
459  break;
460  case Tetragonal:
461  return 1.0 / ( (h2 + k2)/a2 + l2/c2 );
462  break;
463  case Orthorhombic:
464  return 1.0 / ( h2/a2 + k2/b2 + l2/c2 );
465  break;
466  case Rhombohedral:
467  cos2a=cosa*cosa; sin2a=sina*sina;
468  T = h2+k2+l2+2.*(h*k+k*l+h*l) * ((cos2a-cosa)/sin2a);
469  R = sin2a / (1. - 3*cos2a + 2.*cos2a*cosa);
470  return a*a / (T*R);
471  break;
472  case Monoclinic:
473  sin2b=sinb*sinb;
474  return 1./(1./sin2b * (h2/a2+l2/c2-2*h*l*cosb/(a*c)) + k2/b2);
475  break;
476  case Triclinic:
477  return 1./GetRecIntSp2(h,k,l);
478  break;
479  case Hexagonal:
480  return 1. / ( (4.*(h2+k2+h*k) / (3.*a2)) + l2/c2 );
481  break;
482  default:
483  break;
484  }
485 
486  return 0.;
487 }
488 
489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
490 
491 G4double G4CrystalUnitCell::GetRecIntSp2(G4int h,
492  G4int k,
493  G4int l){
494  /* Reference:
495  Table 2.4, pag. 65
496 
497  @Inbook{Ladd2003,
498  author="Ladd, Mark and Palmer, Rex",
499  title="Lattices and Space-Group Theory",
500  bookTitle="Structure Determination by X-ray Crystallography",
501  year="2003",
502  publisher="Springer US",
503  address="Boston, MA",
504  pages="51--116",
505  isbn="978-1-4615-0101-5",
506  doi="10.1007/978-1-4615-0101-5_2",
507  url="http://dx.doi.org/10.1007/978-1-4615-0101-5_2"
508  }
509  */
510 
511  G4double a = theRecSize[0], b = theRecSize[1], c = theRecSize[2];
512  G4double a2 = a*a, b2 = b*b, c2 = c*c;
513  G4double h2 = h*h, k2 = k*k, l2 = l*l;
514 
515  switch(GetLatticeSystem())
516  {
517  case Amorphous:
518  return 0.;
519  break;
520  case Cubic:
521  return a2 * (h2+k2+l2);
522  break;
523  case Tetragonal:
524  return (h2+k2)*a2 + l2*c2 ;
525  break;
526  case Orthorhombic:
527  return h2*a2 + k2+b2 + h2*c2;
528  break;
529  case Rhombohedral:
530  return (h2+k2+l2+2.*(h*k+k*l+h*l) * cosar)*a2;
531  break;
532  case Monoclinic:
533  return h2*a2+k2*b2+l2*c2+2.*h*l*a*c*cosbr;
534  break;
535  case Triclinic:
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;
537  break;
538  case Hexagonal:
539  return (h2+k2+h*k)*a2 + l2*c2;
540  break;
541  default:
542  break;
543  }
544 
545  return 0.;
546 }
547 
548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
549 
550 G4double G4CrystalUnitCell::GetIntCosAng(G4int h1,
551  G4int k1,
552  G4int l1,
553  G4int h2,
554  G4int k2,
555  G4int l2){
556 
557  /* Reference:
558  Table 2.4, pag. 65
559 
560  @Inbook{Kelly2012,
561  author="Anthony A. Kelly and Kevin M. Knowles",
562  title="Appendix 3 Interplanar Spacings and Interplanar Angles",
563  bookTitle="Crystallography and Crystal Defects, 2nd Edition",
564  year="2012",
565  publisher="John Wiley & Sons, Ltd.",
566  isbn="978-0-470-75014-8",
567  doi="10.1002/9781119961468",
568  url="http://onlinelibrary.wiley.com/book/10.1002/9781119961468"
569  }
570  */
571 
572  G4double a = theRecSize[0], b = theRecSize[1], c = theRecSize[2];
573  G4double a2 = a*a, b2 = b*b, c2 = c*c;
574  G4double dsp1dsp2;
575  switch(GetLatticeSystem())
576  {
577  case Amorphous:
578  return 0.;
579  break;
580  case Cubic:
581  return (h1*h2 + k1*k2 + l1+l2) / (std::sqrt(h1*h1 + k1*k1 + l1*l1) * std::sqrt(h2*h2 + k2*k2 + l2*l2));
582  break;
583  case Tetragonal:
584  dsp1dsp2 = std::sqrt(GetIntSp2(h1,k1,l1)*GetIntSp2(h2,k2,l2));
585  return 0. ;
586  break;
587  case Orthorhombic:
588  dsp1dsp2 = std::sqrt(GetIntSp2(h1,k1,l1)*GetIntSp2(h2,k2,l2));
589  return dsp1dsp2 * (h1*h2*a2 + k1*k2*a2 + l1*l2*c2);
590  break;
591  case Rhombohedral:
592  dsp1dsp2 = std::sqrt(GetIntSp2(h1,k1,l1)*GetIntSp2(h2,k2,l2));
593  return dsp1dsp2 * (h1*h2*a2 + k1*k2*b2 + l1*l2*c2+
594  (k1*l2+k2*l1)*b*c*cosar+
595  (h1*l2+h2*l1)*a*c*cosbr+
596  (h1*k2+h2*k1)*a*b*cosgr);
597  break;
598  case Monoclinic:
599  dsp1dsp2 = std::sqrt(GetIntSp2(h1,k1,l1)*GetIntSp2(h2,k2,l2));
600  return dsp1dsp2 * (h1*h2*a2 + k1*k2*b2 + l1*l2*c2+
601  (k1*l2+k2*l1)*b*c*cosar+
602  (h1*l2+h2*l1)*a*c*cosbr+
603  (h1*k2+h2*k1)*a*b*cosgr);
604  break;
605  case Triclinic:
606  dsp1dsp2 = std::sqrt(GetIntSp2(h1,k1,l1)*GetIntSp2(h2,k2,l2));
607  return dsp1dsp2 * (h1*h2*a2 + k1*k2*b2 + l1*l2*c2+
608  (k1*l2+k2*l1)*b*c*cosar+
609  (h1*l2+h2*l1)*a*c*cosbr+
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  }
619 
620  return 0.;
621 }
622 
623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
624