Develop and Download Open Source Software

Powder auto-indexing software Conograph

Browse Subversion Repository

/[conograph]/Conograph/trunk/src/utility_lattice_reduction/put_Buerger_reduced_lattice.hh

Contents of /Conograph/trunk/src/utility_lattice_reduction/put_Buerger_reduced_lattice.hh

Parent Directory Parent Directory | Revision Log Revision Log

Revision 33 - (show annotations) (download) (as text)
Wed Sep 7 04:38:51 2016 UTC (7 years, 6 months ago) by rtomiyasu
File MIME type: text/x-c++hdr
File size: 14127 byte(s) 
The output format for base-centered monoclinic cells was corrected.
1 /*
2 * The MIT License
3
4 Conograph (powder auto-indexing program)
5
6 Copyright (c) <2012> <Ryoko Oishi-Tomiyasu, KEK>
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this software and associated documentation files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 THE SOFTWARE.
25 *
26 */
27 #ifndef PUT_Buerger_REDUCED_LATTICE_HH_
28 #define PUT_Buerger_REDUCED_LATTICE_HH_
29
30 #include "../utility_data_structure/SymMat.hh"
31 #include "../utility_data_structure/VCData.hh"
32 #include "../utility_data_structure/nrutil_nr.hh"
33 #include "../utility_lattice_reduction/put_Selling_reduced_lattice.hh"
34 #include "../point_group/enumPointGroup.hh"
35 #include "../point_group/point_gp_data.hh"
36 #include "../centring_type/enumCentringType.hh"
37 #include "../bravais_type/BravaisType.hh"
38
39
40 template<class T>
41 static void arrangeNondiagonalSign(SymMat<T>& inv_S_red, NRMat<Int4>& trans_mat)
42 {
43 static const T zerro = 0;
44
45 NRMat<Int4> mat2(3,3,0);
46 mat2[0][0] = 1;
47 mat2[1][1] = 1;
48 mat2[2][2] = 1;
49
50 if( zerro < inv_S_red(0,1) )
51 {
52 if( zerro < inv_S_red(0,2) )
53 {
54 if( zerro < inv_S_red(1,2) ) return;
55 else // if( inv_S_red(1,2) <= zerro )
56 {
57 mat2[0][0] = -1;
58 inv_S_red(0,1) *= -1;
59 inv_S_red(0,2) *= -1;
60 }
61 }
62 else // if( inv_S_red(0,2) <= zerro )
63 {
64 if( zerro <= inv_S_red(1,2) )
65 {
66 mat2[1][1] = -1;
67 inv_S_red(0,1) *= -1;
68 inv_S_red(1,2) *= -1;
69 }
70 else // if( inv_S_red(1,2) < zerro )
71 {
72 if( zerro <= inv_S_red(0,2) ) // zerro == inv_S_red(0,2)
73 {
74 mat2[0][0] = -1;
75 inv_S_red(0,1) *= -1;
76 }
77 else // inv_S_red(0,2) < zerro
78 {
79 mat2[2][2] = -1;
80 inv_S_red(0,2) *= -1;
81 inv_S_red(1,2) *= -1;
82 }
83 }
84 }
85 }
86 else //if( inv_S_red(0,1) <= zerro )
87 {
88 if( zerro <= inv_S_red(0,2) )
89 {
90 if( zerro <= inv_S_red(1,2) )
91 {
92 mat2[2][2] = -1;
93 inv_S_red(0,2) *= -1;
94 inv_S_red(1,2) *= -1;
95 }
96 else // if( inv_S_red(1,2) < zerro )
97 {
98 if( inv_S_red(0,2) <= zerro ) return;
99 else if( zerro <= inv_S_red(0,1) ) // zerro == inv_S_red(0,1).
100 {
101 mat2[0][0] = -1;
102 inv_S_red(0,2) *= -1;
103 }
104 else // inv_S_red(0,1) < zerro && inv_S_red(0,2) > zerro.
105 {
106 mat2[1][1] = -1;
107 inv_S_red(0,1) *= -1;
108 inv_S_red(1,2) *= -1;
109 }
110 }
111 }
112 else // if( inv_S_red(0,2) < zerro )
113 {
114 if( zerro < inv_S_red(1,2) )
115 {
116 if( zerro <= inv_S_red(0,1) ) // inv_S_red(0,1) = zerro.
117 {
118 mat2[1][1] = -1;
119 inv_S_red(1,2) *= -1;
120 }
121 else // inv_S_red(0,1) < zerro.
122 {
123 mat2[0][0] = -1;
124 inv_S_red(0,1) *= -1;
125 inv_S_red(0,2) *= -1;
126 }
127 }
128 else // if( inv_S_red(1,2) <= zerro )
129 return;
130 }
131 }
132 trans_mat = mprod(mat2, trans_mat);
133 }
134
135
136 // inv_S_red = trans_mat2*put_transform_matrix_34() * inv_S_super * Transpose(trans_mat2*put_transform_matrix_34()).
137 template<class T>
138 void putBuergerReducedMatrix(const SymMat<T>& inv_S_super, const bool& inv_flag, SymMat<T>& inv_S_red,
139 NRMat<Int4>& trans_mat2)
140 {
141 static const T zerro = 0;
142
143 assert( inv_S_super.size() == 4 );
144 assert( inv_S_red.size() == 3 );
145 assert( inv_S_super(0,0) <= inv_S_super(1,1) );
146 assert( inv_S_super(1,1) <= inv_S_super(2,2) );
147 assert( inv_S_super(2,2) <= inv_S_super(3,3) );
148 assert( inv_S_super(0,1) <= zerro
149 && inv_S_super(0,2) <= zerro
150 && inv_S_super(0,3) <= zerro
151 && inv_S_super(1,2) <= zerro
152 && inv_S_super(1,3) <= zerro
153 && inv_S_super(2,3) <= zerro );
154
155 trans_mat2 = NRMat<Int4>(3,3,0);
156 trans_mat2[0][0] = 1;
157 trans_mat2[1][1] = 1;
158 trans_mat2[2][2] = 1;
159
160 if( inv_S_super(0,0) + inv_S_super(0,1)*2 < zerro )
161 {
162 trans_mat2[1][0] = 1;
163 }
164 else if( inv_S_super(0,0) + inv_S_super(0,2)*2 < zerro )
165 {
166 trans_mat2[2][0] = 1;
167 }
168 else if( inv_S_super(1,1) + inv_S_super(1,2)*2 < zerro )
169 {
170 trans_mat2[2][1] = 1;
171 }
172 else if( inv_S_super(0,0) + inv_S_super(1,1) + (inv_S_super(0,1) + inv_S_super(0,2) + inv_S_super(1,2) )*2 < zerro )
173 {
174 trans_mat2[2][0] = -1;
175 trans_mat2[2][1] = -1;
176 trans_mat2[2][2] = -1;
177 }
178 inv_S_red = transform_sym_matrix(trans_mat2, put_sym_matrix_sizeNplus1toN(inv_S_super));
179 if( inv_flag ) moveSmallerDiagonalLeftUpper(inv_S_red, trans_mat2);
180 else moveLargerDiagonalLeftUpper(inv_S_red, trans_mat2);
181 arrangeNondiagonalSign(inv_S_red, trans_mat2);
182 }
183
184
185
186 template <class T>
187 void cal_average_crystal_system(const ePointGroup& epg, SymMat<T>& ans)
188 {
189 if(epg == Ci) return;
190 else if(epg == C2h_X)
191 {
192 ans(0,1) = 0;
193 ans(0,2) = 0;
194 }
195 else if(epg == C2h_Y)
196 {
197 ans(0,1) = 0;
198 ans(1,2) = 0;
199 }
200 else if(epg == C2h_Z)
201 {
202 ans(0,2) = 0;
203 ans(1,2) = 0;
204 }
205 else if(epg == D2h)
206 {
207 ans(0,1) = 0;
208 ans(0,2) = 0;
209 ans(1,2) = 0;
210 }
211 else if(epg == D4h_X)
212 {
213 ans(0,1) = 0;
214 ans(0,2) = 0;
215 ans(1,2) = 0;
216 ans(1,1) = (ans(1,1)+ans(2,2))/2;
217 ans(2,2) = ans(1,1);
218 }
219 else if(epg == D4h_Y)
220 {
221 ans(0,1) = 0;
222 ans(0,2) = 0;
223 ans(1,2) = 0;
224 ans(0,0) = (ans(0,0)+ans(2,2))/2;
225 ans(2,2) = ans(0,0);
226 }
227 else if(epg == D4h_Z)
228 {
229 ans(0,1) = 0;
230 ans(0,2) = 0;
231 ans(1,2) = 0;
232 ans(0,0) = (ans(0,0)+ans(1,1))/2;
233 ans(1,1) = ans(0,0);
234 }
235 else if(epg == D31d_rho)
236 {
237 ans(0,0) = (ans(0,0)+ans(1,1)+ans(2,2))/3;
238 ans(1,1) = ans(0,0);
239 ans(2,2) = ans(0,0);
240 ans(0,1) = (ans(0,1)+ans(0,2)+ans(1,2))/3;
241 ans(0,2) = ans(0,1);
242 ans(1,2) = ans(0,1);
243 }
244 else if(epg == D3d_1_hex || epg == D6h)
245 {
246 ans(0,0) = (ans(0,0)+ans(1,1))/2;
247 ans(1,1) = ans(0,0);
248 ans(0,1) = ans(0,0)/2;
249 ans(0,2) = 0;
250 ans(1,2) = 0;
251 }
252 else if(epg == Oh)
253 {
254 ans(0,0) = (ans(0,0)+ans(1,1)+ans(2,2))/3;
255 ans(1,1) = ans(0,0);
256 ans(2,2) = ans(0,0);
257 ans(0,1) = 0;
258 ans(0,2) = 0;
259 ans(1,2) = 0;
260 }
261 else
262 {
263 assert( false );
264 }
265 };
266
267 inline void putBuergerReducedMatrix(const SymMat<Double>& inv_S_super, SymMat<Double>& inv_S_red, NRMat<Int4>& trans_mat2)
268 {
269 putBuergerReducedMatrix(inv_S_super, true, inv_S_red, trans_mat2);
270
271 #ifdef DEBUG
272 assert( ( inv_S_red(0,1) <= 0.0 && inv_S_red(0,2) <= 0.0 && inv_S_red(1,2) <= 0.0 )
273 || ( 0.0 < inv_S_red(0,1) && 0.0 < inv_S_red(0,2) && 0.0 < inv_S_red(1,2) ) );
274 assert( inv_S_red(1,1)*0.9999 < inv_S_red(2,2) );
275 assert( inv_S_red(0,0)*0.9999 < inv_S_red(1,1) );
276 assert( inv_S_red(0,1) * (-1.9999) < inv_S_red(0,0)
277 && inv_S_red(0,1) * 1.9999 < inv_S_red(0,0)
278 && inv_S_red(0,2) * (-1.9999) < inv_S_red(0,0)
279 && inv_S_red(0,2) * 1.9999 < inv_S_red(0,0)
280 && inv_S_red(1,2) * (-1.9999) < inv_S_red(1,1)
281 && inv_S_red(1,2) * 1.9999 < inv_S_red(1,1) );
282 assert( 0.0 < inv_S_red(0,0) + inv_S_red(1,1) + ( inv_S_red(0,1) + inv_S_red(0,2) + inv_S_red(1,2) ) * 1.9999
283 && 0.0 < inv_S_red(0,0) + inv_S_red(1,1) + ( inv_S_red(0,1) - inv_S_red(0,2) - inv_S_red(1,2) ) * 1.9999
284 && 0.0 < inv_S_red(0,0) + inv_S_red(1,1) - ( inv_S_red(0,1) - inv_S_red(0,2) + inv_S_red(1,2) ) * 1.9999
285 && 0.0 < inv_S_red(0,0) + inv_S_red(1,1) - ( inv_S_red(0,1) + inv_S_red(0,2) - inv_S_red(1,2) ) * 1.9999 );
286 #endif
287 }
288
289
290 inline void putBuergerReducedMatrix(const SymMat<VCData>& S_super, SymMat<VCData>& S_red, NRMat<Int4>& trans_mat2)
291 {
292 putBuergerReducedMatrix(S_super, false, S_red, trans_mat2);
293
294 #ifdef DEBUG
295 assert( ( S_red(0,1).Value() <= 0.0 && S_red(0,2).Value() <= 0.0 && S_red(1,2).Value() <= 0.0 )
296 || ( 0.0 < S_red(0,1).Value() && 0.0 < S_red(0,2).Value() && 0.0 < S_red(1,2).Value() ) );
297 assert( S_red(1,1).Value()*0.9999 < S_red(0,0).Value() );
298 assert( S_red(2,2).Value()*0.9999 < S_red(1,1).Value() );
299 assert( S_red(0,1).Value() * (-1.9999) < S_red(1,1).Value()
300 && S_red(0,1).Value() * 1.9999 < S_red(1,1).Value()
301 && S_red(0,2).Value() * (-1.9999) < S_red(2,2).Value()
302 && S_red(0,2).Value() * 1.9999 < S_red(2,2).Value()
303 && S_red(1,2).Value() * (-1.9999) < S_red(2,2).Value()
304 && S_red(1,2).Value() * 1.9999 < S_red(2,2).Value() );
305 assert( 0.0 < S_red(1,1).Value() + S_red(2,2).Value() + ( S_red(0,1) + S_red(0,2) + S_red(1,2) ).Value() * 1.9999
306 && 0.0 < S_red(1,1).Value() + S_red(2,2).Value() + ( S_red(0,1) - S_red(0,2) - S_red(1,2) ).Value() * 1.9999
307 && 0.0 < S_red(1,1).Value() + S_red(2,2).Value() - ( S_red(0,1) - S_red(0,2) + S_red(1,2) ).Value() * 1.9999
308 && 0.0 < S_red(1,1).Value() + S_red(2,2).Value() - ( S_red(0,1) + S_red(0,2) - S_red(1,2) ).Value() * 1.9999 );
309 #endif
310 }
311
312
313 template<class T>
314 void putBuergerReducedMonoclinicP(const Int4& i, const Int4& j,
315 SymMat<T>& S_red, NRMat<Int4>& trans_mat2)
316 {
317 static const T zerro = 0;
318
319 assert(S_red.size()==3);
320 assert(trans_mat2.ncols()==3);
321 assert(i < j);
322 const Int4 irow = trans_mat2.nrows();
323
324 if( S_red(i,j) < zerro )
325 {
326 S_red(i,j) *= -1;
327 for(Int4 l=0; l<irow; l++)
328 {
329 trans_mat2[l][i] *= -1;
330 }
331 }
332
333 do{
334 if( S_red(j,j) < S_red(i,j) * 2 )
335 {
336 // S_red(i,j).Value() <= S_red(j,j).Value() * m
337 // i : -1 m 0
338 // j : 0 1 0
339 // : 0 0 1
340 const Int8 m = iceil( S_red(i,j) / S_red(j,j) );
341
342 S_red(i,i) += ( S_red(j,j) * m - S_red(i,j) * 2 ) * m;
343 assert( zerro < S_red(i,i) );
344 S_red(i,j) = S_red(j,j) * m - S_red(i,j);
345
346 for(Int4 l=0; l<irow; l++)
347 {
348 trans_mat2[l][j] += trans_mat2[l][i]*m;
349 }
350 for(Int4 l=0; l<irow; l++)
351 {
352 trans_mat2[l][i] *= -1;
353 }
354 }
355
356 if( S_red(i,i) < S_red(i,j) * 2 )
357 {
358 // S_red(i,j).Value() <= S_red(i,i).Value() * n
359 // i : 1 0 0
360 // j : n -1 0
361 // : 0 0 1
362 const Int8 n = iceil( S_red(i,j) / S_red(i,i) );
363
364 S_red(j,j) += ( S_red(i,i) * n - S_red(i,j) * 2 ) * n;
365 assert( zerro < S_red(j,j) );
366 S_red(i,j) = S_red(i,i) * n - S_red(i,j);
367
368 for(Int4 l=0; l<irow; l++)
369 {
370 trans_mat2[l][i] += trans_mat2[l][j]*n;
371 }
372 for(Int4 l=0; l<irow; l++)
373 {
374 trans_mat2[l][j] *= -1;
375 }
376 }
377 }
378 while( S_red(j,j) < S_red(i,j) * 2 || S_red(i,i) < S_red(i,j) * 2 );
379
380 if( S_red(i,i) < S_red(j,j) )
381 {
382 const Int4 k = put_complement_set3(i, j);
383 swap(S_red(i,i), S_red(j,j));
384 swap(S_red(i,k), S_red(j,k));
385 for(Int4 l=0; l<irow; l++)
386 {
387 swap(trans_mat2[l][i], trans_mat2[l][j]);
388 }
389 }
390 }
391
392
393 inline Double operator/(const VCData& lhs, const VCData& rhs)
394 {
395 assert(rhs.Value() != 0.0);
396 return lhs.Value() / rhs.Value();
397 }
398
399 template<class T>
400 void putBuergerReducedMonoclinicB(
401 const BravaisType& monoclinic_b_type,
402 SymMat<T>& S_red, NRMat<Int4>& trans_mat2)
403 {
404 const Int4 ibase_axis = monoclinic_b_type.enumBASEaxis();
405 const Int4 iabc_axis = monoclinic_b_type.enumABCaxis();
406 const Int4 ir = put_complement_set3(iabc_axis, ibase_axis);
407
408 static const T zerro = 0;
409
410 assert(S_red.size()==3);
411 assert(trans_mat2.nrows()==0 || trans_mat2.ncols()==3);
412 const Int4 irow = trans_mat2.nrows();
413
414 if( S_red( ibase_axis, ir ) < zerro )
415 {
416 S_red( ibase_axis, ir ) *= -1;
417 for(Int4 l=0; l<irow; l++)
418 {
419 trans_mat2[l][ir] *= -1;
420 }
421 }
422
423 do{
424 if( S_red( ir, ir ) < S_red( ibase_axis, ir ) )
425 {
426 // S_red( ibase_axis, ir ).Value() <= S_red( ir, ir ).Value() * m2
427 // ir : 1 0 0
428 // ibase_axis : m2 -1 0
429 // : 0 0 1
430 const Int8 m1 = iceil( ( S_red( ibase_axis, ir ) / S_red( ir, ir ) ) * 0.5 );
431 const Int8 m2 = m1*2;
432
433 S_red( ibase_axis, ibase_axis ) += ( S_red( ir, ir ) * m2 - S_red( ibase_axis, ir ) * 2 ) * m2;
434 assert( zerro < S_red( ibase_axis, ibase_axis ) );
435 S_red( ibase_axis, ir ) = S_red( ir, ir ) * m2 - S_red( ibase_axis, ir );
436
437 for(Int4 l=0; l<irow; l++)
438 {
439 trans_mat2[l][ir] += trans_mat2[l][ibase_axis]*m2;
440 }
441 for(Int4 l=0; l<irow; l++)
442 {
443 trans_mat2[l][ibase_axis] *= -1;
444 }
445 }
446
447 if( S_red( ibase_axis, ibase_axis ) < S_red( ibase_axis, ir ) * 2 )
448 {
449 // S_red( ibase_axis, ir ).Value() <= S_red( ibase_axis, ibase_axis ).Value() * n
450 // ir :-1 n 0
451 // ibase_axis : 0 1 0
452 // : 0 0 1
453 const Int4 n = iceil( S_red( ibase_axis, ir ) / S_red( ibase_axis, ibase_axis ) );
454
455 S_red( ir, ir ) += ( S_red( ibase_axis, ibase_axis ) * n - S_red( ibase_axis, ir ) * 2 ) * n;
456 assert( zerro < S_red( ir, ir ) );
457 S_red( ibase_axis, ir ) = S_red( ibase_axis, ibase_axis ) * n - S_red( ibase_axis, ir );
458
459 for(Int4 l=0; l<irow; l++)
460 {
461 trans_mat2[l][ibase_axis] += trans_mat2[l][ir]*n;
462 }
463 for(Int4 l=0; l<irow; l++)
464 {
465 trans_mat2[l][ir] *= -1;
466 }
467 }
468 }
469 while( S_red( ir, ir ) < S_red( ibase_axis, ir ) || S_red( ibase_axis, ibase_axis ) < S_red( ibase_axis, ir ) * 2 );
470 }
471
472
473 template<class T>
474 void putBuergerReducedOrthorhombic(const eCentringType& brat,
475 SymMat<T>& S_red, NRMat<Int4>& trans_mat)
476 {
477 assert( brat != BaseX );
478 assert( brat != BaseY );
479
480 if( brat == BaseZ )
481 {
482 if( S_red(0,0) < S_red(1,1) )
483 {
484 S_red = transform_sym_matrix(put_matrix_YXZ(), S_red);
485 trans_mat = mprod( trans_mat, put_matrix_YXZ() );
486 }
487 }
488 else
489 {
490 NRMat<Int4> trans_mat2 = identity_matrix<Int4>(3);
491 moveLargerDiagonalLeftUpper(S_red, trans_mat2);
492 trans_mat = mprod(trans_mat, transpose(trans_mat2)); // inverse(trans_mat2) = transpose(trans_mat2).
493 }
494 }
495
496 #endif /*PUT_MINKOWSKI_REDUCED_LATTICE_HH_*/
Back to OSDN">Back to OSDN
 ViewVC Help
Powered by ViewVC 1.1.26  
About OSDN
Find Software
Develop Software
Help