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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 33 - (hide 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 rtomiyasu 25 /*
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 rtomiyasu 33 static const T zerro = 0;
318 rtomiyasu 25
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 rtomiyasu 33 assert( zerro < S_red(i,i) );
344 rtomiyasu 25 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 rtomiyasu 33 assert( zerro < S_red(j,j) );
366 rtomiyasu 25 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 rtomiyasu 33 static const T zerro = 0;
409 rtomiyasu 25
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 rtomiyasu 33 assert( zerro < S_red( ibase_axis, ibase_axis ) );
435 rtomiyasu 25 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 rtomiyasu 33 assert( zerro < S_red( ir, ir ) );
457 rtomiyasu 25 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