Develop and Download Open Source Software

Powder auto-indexing software Conograph

Browse Subversion Repository

/[conograph]/Conograph/trunk/src/lattice_symmetry/LatticeFigureOfMeritToCheckSymmetry.cc

Contents of /Conograph/trunk/src/lattice_symmetry/LatticeFigureOfMeritToCheckSymmetry.cc

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++src
File size: 15409 byte(s) 
The output format for base-centered monoclinic cells was corrected.
1 /*
2 * The MIT License
3
4 BLDConograph (Bravais lattice determination module in Conograph)
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 #include <limits>
28 #include "../utility_lattice_reduction/put_Buerger_reduced_lattice.hh"
29 #include "check_equiv.hh"
30 #include "LatticeFigureOfMeritToCheckSymmetry.hh"
31
32 LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry()
33 : m_S_red( SymMat43_Double( SymMat<Double>(3), NRMat<Int4>(4,3) ) )
34 {
35 }
36
37
38 LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const Double& rhs)
39 : m_latfom(rhs),
40 m_S_red( SymMat43_Double( SymMat<Double>(3), NRMat<Int4>(4,3) ) )
41 {
42 }
43
44
45 LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const BravaisType& brat,
46 const SymMat43_Double& S)
47 : m_S_red( SymMat43_Double( SymMat<Double>(3), NRMat<Int4>(4,3) ) )
48 {
49 setLatticeConstants43(brat, S);
50 }
51
52
53
54
55 #ifdef DEBUG
56 static bool checkInitialLatticeParameters(
57 const BravaisType& brat,
58 const SymMat43_Double& S_red)
59 {
60 const SymMat<Double> dbl_S_red( S_red.first );
61
62 if( brat.enumLaueGroup() == Ci && brat.enumCentringType() == Prim )
63 {
64 assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)
65 && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1)
66 && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2)
67 && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );
68 }
69 else if( brat.enumLaueGroup() == C2h_Y && brat.enumCentringType() == Prim )
70 {
71 assert( 0.0 <= dbl_S_red(0,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(0,0)
72 && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(0,0) );
73 }
74 else if( brat.enumLaueGroup() == C2h_Z && brat.enumCentringType() == Prim )
75 {
76 assert( 0.0 <= dbl_S_red(0,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)
77 && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1) );
78 }
79 else if( brat.enumLaueGroup() == C2h_X && brat.enumCentringType() == Prim )
80 {
81 assert( 0.0 <= dbl_S_red(1,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1)
82 && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );
83 }
84 else if( brat.enumLaueGroup() == C2h_Y && brat.enumCentringType() == BaseZ )
85 {
86 assert( 0.0 <= dbl_S_red(0,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 0.9999 < dbl_S_red(0,0) );
87 }
88 else if( brat.enumLaueGroup() == C2h_Z && brat.enumCentringType() == BaseX )
89 {
90 assert( 0.0 <= dbl_S_red(0,1) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 0.9999 < dbl_S_red(1,1) );
91 }
92 else if( brat.enumLaueGroup() == C2h_X && brat.enumCentringType() == BaseY )
93 {
94 assert( 0.0 <= dbl_S_red(1,2) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 0.9999 < dbl_S_red(2,2) );
95 }
96 else if( brat.enumLaueGroup() == D2h
97 && brat.enumCentringType() != BaseX
98 && brat.enumCentringType() != BaseY
99 && brat.enumCentringType() != BaseZ )
100 {
101 assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0) );
102 }
103
104 const SymMat<Double> S_super = transform_sym_matrix(S_red.second, S_red.first);
105 assert( S_super(0,1) <= 0.0
106 && S_super(0,2) <= 0.0
107 && S_super(0,3) <= 0.0
108 && S_super(1,2) <= 0.0
109 && S_super(1,3) <= 0.0
110 && S_super(2,3) <= 0.0 );
111
112 return true;
113 }
114 #endif
115
116
117 void LatticeFigureOfMeritToCheckSymmetry::setLatticeConstants43(const BravaisType& brat, const SymMat43_Double& S)
118 {
119 m_S_red = S;
120 assert( checkInitialLatticeParameters(brat, m_S_red) );
121
122 m_latfom.setLatticeConstants43(brat, S);
123 }
124
125
126
127 static bool operator<(const SymMat<Double>& lhs, const SymMat<Double>& rhs)
128 {
129 static const Double EPS_1 = 1.0+sqrt( numeric_limits<double>::epsilon() );
130 assert( lhs.size() == 3 );
131 assert( rhs.size() == 3 );
132
133 static const Int4 ISIZE = 3;
134 for(Int4 i=0; i<ISIZE; i++)
135 {
136 if( lhs(i,i)*EPS_1 < rhs(i,i) ) return true;
137 if( rhs(i,i)*EPS_1 < lhs(i,i) ) return false;
138
139 for(Int4 j=0; j<i; j++)
140 {
141 const Double lhs_ij = lhs(i,i)+lhs(j,j)+lhs(i,j)*2.0;
142 const Double rhs_ij = rhs(i,i)+rhs(j,j)+rhs(i,j)*2.0;
143
144 if( lhs_ij*EPS_1 < rhs_ij ) return true;
145 if( rhs_ij*EPS_1 < lhs_ij ) return false;
146 }
147 }
148
149 return false;
150 }
151
152
153
154 bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsMonoclinic(const ePointGroup& epg_new,
155 const Double& resol,
156 map< SymMat<Double>, NRMat<Int4> >& ans) const
157 {
158 ans.clear();
159
160 SymMat<Double> ans0 = m_S_red.first;
161 cal_average_crystal_system(C2h_X, ans0);
162
163 SymMat<Double> S_red(3);
164 NRMat<Int4> trans_mat2;
165 if( check_equiv_m(ans0, m_S_red.first, resol ) )
166 {
167 if( epg_new == C2h_X )
168 {
169 S_red = ans0;
170 trans_mat2 = m_S_red.second;
171 putBuergerReducedMonoclinicP<Double>(1, 2, S_red, trans_mat2);
172 }
173 else if( epg_new == C2h_Y )
174 {
175 S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);
176 trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());
177 putBuergerReducedMonoclinicP<Double>(0, 2, S_red, trans_mat2);
178 }
179 else // if( epg_new == C2h_Z )
180 {
181 // (0 0 1)(0 1 0)
182 // (1 0 0)(0 0 1)
183 // (0 1 0)(1 0 0)
184 S_red = transform_sym_matrix(put_matrix_YZX(), ans0);
185 trans_mat2 = mprod(m_S_red.second, put_matrix_ZXY());
186 putBuergerReducedMonoclinicP<Double>(0, 1, S_red, trans_mat2);
187 }
188 ans.insert( SymMat43_Double( S_red, trans_mat2) );
189 }
190
191 ans0 = m_S_red.first;
192 cal_average_crystal_system(C2h_Y, ans0);
193 if( check_equiv_m(ans0, m_S_red.first, resol ) )
194 {
195 if( epg_new == C2h_X )
196 {
197 S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);
198 trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());
199 putBuergerReducedMonoclinicP<Double>(1, 2, S_red, trans_mat2);
200 }
201 else if( epg_new == C2h_Y )
202 {
203 S_red = ans0;
204 trans_mat2 = m_S_red.second;
205 putBuergerReducedMonoclinicP<Double>(0, 2, S_red, trans_mat2);
206 }
207 else // if( epg_new == C2h_Z )
208 {
209 S_red = transform_sym_matrix(put_matrix_XZY(), ans0);
210 trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());
211 putBuergerReducedMonoclinicP<Double>(0, 1, S_red, trans_mat2);
212 }
213 ans.insert( SymMat43_Double( S_red, trans_mat2) );
214 }
215
216 ans0 = m_S_red.first;
217 cal_average_crystal_system(C2h_Z, ans0);
218 if( check_equiv_m(ans0, m_S_red.first, resol ) )
219 {
220 if( epg_new == C2h_X )
221 {
222 // (0 1 0)(0 0 1)
223 // (0 0 1)(1 0 0)
224 // (1 0 0)(0 1 0)
225 S_red = transform_sym_matrix(put_matrix_ZXY(), ans0);
226 trans_mat2 = mprod(m_S_red.second, put_matrix_YZX());
227 putBuergerReducedMonoclinicP<Double>(1, 2, S_red, trans_mat2);
228 }
229 else if( epg_new == C2h_Y )
230 {
231 S_red = transform_sym_matrix(put_matrix_XZY(), ans0);
232 trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());
233 putBuergerReducedMonoclinicP<Double>(0, 2, S_red, trans_mat2);
234 }
235 else // if( epg_new == C2h_Z )
236 {
237 S_red = ans0;
238 trans_mat2 = m_S_red.second;
239 putBuergerReducedMonoclinicP<Double>(0, 1, S_red, trans_mat2);
240 }
241 ans.insert( SymMat43_Double( S_red, trans_mat2) );
242 }
243
244 return !( ans.empty() );
245 }
246
247
248 bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsOrthorhombic(const Double& resol,
249 map< SymMat<Double>, NRMat<Int4> >& ans) const
250 {
251 ans.clear();
252
253 const BravaisType& brat = m_latfom.putBravaisType();
254
255 SymMat<Double> ans0 = m_S_red.first;
256 cal_average_crystal_system(D2h, ans0);
257 if( check_equiv_m(ans0, m_S_red.first, resol ) )
258 {
259 if( brat.enumCentringType() == BaseX )
260 {
261 if( ans0(1,1) < ans0(2,2) )
262 {
263 ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_ZYX(), ans0), mprod( m_S_red.second, put_matrix_ZYX() ) ) );
264 }
265 else
266 {
267 ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );
268 }
269 }
270 else if( brat.enumCentringType() == BaseY )
271 {
272 if( ans0(0,0) < ans0(2,2) )
273 {
274 ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_ZXY(), ans0), mprod( m_S_red.second, put_matrix_YZX() ) ) );
275 }
276 else
277 {
278 ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );
279 }
280 }
281 else if( brat.enumCentringType() == BaseZ )
282 {
283 if( ans0(0,0) < ans0(1,1) )
284 {
285 ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_YXZ(), ans0), mprod( m_S_red.second, put_matrix_YXZ() ) ) );
286 }
287 else
288 {
289 ans.insert( SymMat43_Double( transform_sym_matrix(put_matrix_XYZ(), ans0), m_S_red.second ) );
290 }
291 }
292 else
293 {
294 NRMat<Int4> trans_mat = m_S_red.second;
295 putBuergerReducedOrthorhombic(brat.enumCentringType(), ans0, trans_mat);
296 ans.insert( SymMat43_Double(ans0, trans_mat ) );
297 }
298 return true;
299 }
300 return false;
301 }
302
303
304 bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsTetragonal(const Double& resol,
305 map< SymMat<Double>, NRMat<Int4> >& ans) const
306 {
307 ans.clear();
308
309 SymMat<Double> ans0 = m_S_red.first;
310 cal_average_crystal_system(D4h_X, ans0);
311 if( check_equiv_m(ans0, m_S_red.first, resol ) )
312 {
313 ans.insert( SymMat43_Double(
314 transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );
315 }
316
317 ans0 = m_S_red.first;
318 cal_average_crystal_system(D4h_Y, ans0);
319 if( check_equiv_m(ans0, m_S_red.first, resol ) )
320 {
321 ans.insert( SymMat43_Double(
322 transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );
323 }
324
325 ans0 = m_S_red.first;
326 cal_average_crystal_system(D4h_Z, ans0);
327 if( check_equiv_m(ans0, m_S_red.first, resol ) )
328 {
329 ans.insert( SymMat43_Double(ans0, m_S_red.second ) );
330 }
331
332 return !( ans.empty() );
333 }
334
335
336
337
338 bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsHexagonal(const ePointGroup& epg_new, const Double& resol,
339 map< SymMat<Double>, NRMat<Int4> >& ans) const
340 {
341 ans.clear();
342 const BravaisType& brat = m_latfom.putBravaisType();
343
344 SymMat43_Double ans2(SymMat<Double>(3), NRMat<Int4>(3,3));
345
346 if( brat.enumLaueGroup() == C2h_X )
347 {
348 ans2.first = transform_sym_matrix(put_matrix_YZX(), m_S_red.first);
349 ans2.second = mprod( m_S_red.second, put_matrix_ZXY() );
350 }
351 else if( brat.enumLaueGroup() == C2h_Y )
352 {
353 ans2.first = transform_sym_matrix(put_matrix_ZXY(), m_S_red.first);
354 ans2.second = mprod( m_S_red.second, put_matrix_YZX() );
355 }
356 else // if( brat.enumLaueGroup() == C2h_Z )
357 {
358 ans2.first = transform_sym_matrix(put_matrix_XYZ(), m_S_red.first);
359 ans2.second = m_S_red.second;
360 }
361
362 if( ans2.first(0,1) < 0.0 )
363 {
364 ans2.first(0,1) *= -1;
365 ans2.second[0][0] *= -1;
366 ans2.second[1][0] *= -1;
367 ans2.second[2][0] *= -1;
368 }
369
370 SymMat<Double> ans0 = ans2.first;
371 cal_average_crystal_system(epg_new, ans2.first);
372 if( check_equiv_m(ans2.first, ans0, resol ) )
373 {
374 ans.insert( ans2 );
375 return true;
376 }
377 else return false;
378 }
379
380
381 bool LatticeFigureOfMeritToCheckSymmetry::checkLatticeSymmetry(const ePointGroup& epg_new, const Double& resol,
382 map< SymMat<Double>, NRMat<Int4> >& ans) const
383 {
384 ans.clear();
385 const BravaisType& brat = m_latfom.putBravaisType();
386 if( epg_new == Ci || epg_new == brat.enumLaueGroup() )
387 {
388 ans.insert( m_S_red );
389 return true;
390 }
391
392 if( epg_new == C2h_X || epg_new == C2h_Y || epg_new == C2h_Z )
393 {
394 assert( brat.enumLaueGroup() == Ci );
395 assert( brat.enumCentringType() == Prim );
396
397 return checkIfLatticeIsMonoclinic(epg_new, resol, ans);
398 }
399 else if( epg_new == D4h_Z )
400 {
401 assert( brat.enumLaueGroup() == D2h );
402 assert( brat.enumCentringType() == Prim
403 || brat.enumCentringType() == Inner );
404
405 return checkIfLatticeIsTetragonal(resol, ans);
406 }
407 else if( epg_new == D2h )
408 {
409 assert( brat.enumLaueGroup() != Ci || brat.enumCentringType() == Prim );
410 assert( brat.enumLaueGroup() != C2h_Z || brat.enumCentringType() == BaseX );
411 assert( brat.enumLaueGroup() != C2h_X || brat.enumCentringType() == BaseY );
412 assert( brat.enumLaueGroup() != C2h_Y || brat.enumCentringType() == BaseZ );
413 assert( brat.enumCentringType() != Rhom_hex );
414
415 return checkIfLatticeIsOrthorhombic(resol, ans);
416 }
417 else if( epg_new == D6h )
418 {
419 assert( brat.enumCentringType() == Prim );
420 assert( brat.enumLaueGroup() == C2h_X
421 || brat.enumLaueGroup() == C2h_Y
422 || brat.enumLaueGroup() == C2h_Z );
423 return checkIfLatticeIsHexagonal(epg_new, resol, ans);
424 }
425 else
426 {
427 assert( epg_new == Oh );
428 assert( brat.enumCentringType() == Prim
429 || brat.enumCentringType() == Inner
430 || brat.enumCentringType() == Face );
431
432 SymMat43_Double ans2 = m_S_red;
433 cal_average_crystal_system(epg_new, ans2.first);
434 if( check_equiv_m(ans2.first, m_S_red.first, resol ) )
435 {
436 ans.insert( ans2 );
437 return true;
438 }
439 }
440 return !(ans.empty());
441 }
442
443
444 void LatticeFigureOfMeritToCheckSymmetry::putLatticesOfHigherSymmetry(
445 const ePointGroup& epg, const Double& resol,
446 vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result) const
447 {
448 lattice_result.clear();
449 map< SymMat<Double>, NRMat<Int4> > S_red_tray;
450 if( !this->checkLatticeSymmetry(epg, resol, S_red_tray) ) return;
451
452 const BravaisType& ebrat_original = this->putLatticeFigureOfMerit().putBravaisType();
453 const eCentringType eblat = (ebrat_original.enumBravaisType()==Monoclinic_B?
454 (epg==D31d_rho?Prim:(epg==D3d_1_hex?Rhom_hex:BaseZ)):ebrat_original.enumCentringType());
455
456 const NRMat<Int4> matrix_min_to_sell = this->putInitialForm().second;
457
458 SymMat<Double> S_super(4);
459 NRMat<Int4> trans_mat(4,3);
460
461 for(map< SymMat<Double>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
462 {
463 // S_super = it->second * it->first * Transpose(it->second) is close to
464 // Delone-reduced form of the original lattice.
465 S_super = transform_sym_matrix(it->second, it->first );
466
467 trans_mat = identity_matrix<Int4>(4);
468
469 // S_super = trans_mat * it->second * it->first * Transpose(trans_mat * it->second).
470 put_Selling_reduced_dim_less_than_4(S_super, trans_mat);
471 moveSmallerDiagonalLeftUpper(S_super, trans_mat);
472
473 lattice_result.push_back( LatticeFigureOfMeritToCheckSymmetry( BravaisType( pair<eCentringType, ePointGroup>(eblat, epg) ),
474 SymMat43_Double(it->first, mprod(trans_mat, it->second) ) ) );
475 }
476 }
Back to OSDN">Back to OSDN
 ViewVC Help
Powered by ViewVC 1.1.26  
About OSDN
Find Software
Develop Software
Help