1 /* pointgroup.c */
2 /* Copyright (C) 2008 Atsushi Togo */
3
4 #include <string.h>
5 #include <stdio.h>
6 #include "lattice.h"
7 #include "pointgroup.h"
8 #include "symmetry.h"
9 #include "mathfunc.h"
10
11 #include "debug.h"
12
13 #define NUM_ROT_AXES 73
14
15 typedef struct {
16 int table[10];
17 char symbol[6];
18 Holohedry holohedry;
19 Laue laue;
20 } PointgroupType;
21
22 static PointgroupType pointgroup_data[32] = {
23 {
24 {0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
25 "1 ",
26 TRICLI,
27 LAUE1,
28 },
29 {
30 {0, 0, 0, 0, 1, 1, 0, 0, 0, 0},
31 "-1 ",
32 TRICLI,
33 LAUE1,
34 },
35 {
36 {0, 0, 0, 0, 0, 1, 1, 0, 0, 0},
37 "2 ",
38 MONOCLI,
39 LAUE2M,
40 },
41 {
42 {0, 0, 0, 1, 0, 1, 0, 0, 0, 0},
43 "m ",
44 MONOCLI,
45 LAUE2M,
46 },
47 {
48 {0, 0, 0, 1, 1, 1, 1, 0, 0, 0},
49 "2/m ",
50 MONOCLI,
51 LAUE2M,
52 },
53 {
54 {0, 0, 0, 0, 0, 1, 3, 0, 0, 0},
55 "222 ",
56 ORTHO,
57 LAUEMMM,
58 },
59 {
60 {0, 0, 0, 2, 0, 1, 1, 0, 0, 0},
61 "mm2 ",
62 ORTHO,
63 LAUEMMM,
64 },
65 {
66 {0, 0, 0, 3, 1, 1, 3, 0, 0, 0},
67 "mmm ",
68 ORTHO,
69 LAUEMMM,
70 },
71 {
72 {0, 0, 0, 0, 0, 1, 1, 0, 2, 0},
73 "4 ",
74 TETRA,
75 LAUE4M,
76 },
77 {
78 {0, 2, 0, 0, 0, 1, 1, 0, 0, 0},
79 "-4 ",
80 TETRA,
81 LAUE4M,
82 },
83 {
84 {0, 2, 0, 1, 1, 1, 1, 0, 2, 0},
85 "4/m ",
86 TETRA,
87 LAUE4M,
88 },
89 {
90 {0, 0, 0, 0, 0, 1, 5, 0, 2, 0},
91 "422 ",
92 TETRA,
93 LAUE4MMM,
94 },
95 {
96 {0, 0, 0, 4, 0, 1, 1, 0, 2, 0},
97 "4mm ",
98 TETRA,
99 LAUE4MMM,
100 },
101 {
102 {0, 2, 0, 2, 0, 1, 3, 0, 0, 0},
103 "-42m ",
104 TETRA,
105 LAUE4MMM,
106 },
107 {
108 {0, 2, 0, 5, 1, 1, 5, 0, 2, 0},
109 "4/mmm",
110 TETRA,
111 LAUE4MMM,
112 },
113 {
114 {0, 0, 0, 0, 0, 1, 0, 2, 0, 0},
115 "3 ",
116 TRIGO,
117 LAUE3,
118 },
119 {
120 {0, 0, 2, 0, 1, 1, 0, 2, 0, 0},
121 "-3 ",
122 TRIGO,
123 LAUE3,
124 },
125 {
126 {0, 0, 0, 0, 0, 1, 3, 2, 0, 0},
127 "32 ",
128 TRIGO,
129 LAUE3M,
130 },
131 {
132 {0, 0, 0, 3, 0, 1, 0, 2, 0, 0},
133 "3m ",
134 TRIGO,
135 LAUE3M,
136 },
137 {
138 {0, 0, 2, 3, 1, 1, 3, 2, 0, 0},
139 "-3m ",
140 TRIGO,
141 LAUE3M,
142 },
143 {
144 {0, 0, 0, 0, 0, 1, 1, 2, 0, 2},
145 "6 ",
146 HEXA,
147 LAUE6M,
148 },
149 {
150 {2, 0, 0, 1, 0, 1, 0, 2, 0, 0},
151 "-6 ",
152 HEXA,
153 LAUE6M,
154 },
155 {
156 {2, 0, 2, 1, 1, 1, 1, 2, 0, 2},
157 "6/m ",
158 HEXA,
159 LAUE6M,
160 },
161 {
162 {0, 0, 0, 0, 0, 1, 7, 2, 0, 2},
163 "622 ",
164 HEXA,
165 LAUE6MMM,
166 },
167 {
168 {0, 0, 0, 6, 0, 1, 1, 2, 0, 2},
169 "6mm ",
170 HEXA,
171 LAUE6MMM,
172 },
173 {
174 {2, 0, 0, 4, 0, 1, 3, 2, 0, 0},
175 "-62m ",
176 HEXA,
177 LAUE6MMM,
178 },
179 {
180 {2, 0, 2, 7, 1, 1, 7, 2, 0, 2},
181 "6/mmm",
182 HEXA,
183 LAUE6MMM,
184 },
185 {
186 {0, 0, 0, 0, 0, 1, 3, 8, 0, 0},
187 "23 ",
188 CUBIC,
189 LAUEM3,
190 },
191 {
192 {0, 0, 8, 3, 1, 1, 3, 8, 0, 0},
193 "m-3 ",
194 CUBIC,
195 LAUEM3,
196 },
197 {
198 {0, 0, 0, 0, 0, 1, 9, 8, 6, 0},
199 "432 ",
200 CUBIC,
201 LAUEM3M,
202 },
203 {
204 {0, 6, 0, 6, 0, 1, 3, 8, 0, 0},
205 "-43m ",
206 CUBIC,
207 LAUEM3M,
208 },
209 {
210 {0, 6, 8, 9, 1, 1, 9, 8, 6, 0},
211 "m-3m ",
212 CUBIC,
213 LAUEM3M,
214 }
215 };
216
217 static int identity[3][3] = {
218 { 1, 0, 0},
219 { 0, 1, 0},
220 { 0, 0, 1},
221 };
222
223 static int inversion[3][3] = {
224 {-1, 0, 0},
225 { 0,-1, 0},
226 { 0, 0,-1},
227 };
228
229 static int rot_axes[][3] = {
230 { 1, 0, 0},
231 { 0, 1, 0},
232 { 0, 0, 1},
233 { 0, 1, 1},
234 { 1, 0, 1},
235 { 1, 1, 0},
236 { 0, 1,-1},
237 {-1, 0, 1},
238 { 1,-1, 0},
239 { 1, 1, 1}, /* 10 */
240 {-1, 1, 1},
241 { 1,-1, 1},
242 { 1, 1,-1},
243 { 0, 1, 2},
244 { 2, 0, 1},
245 { 1, 2, 0},
246 { 0, 2, 1},
247 { 1, 0, 2},
248 { 2, 1, 0},
249 { 0,-1, 2}, /* 20 */
250 { 2, 0,-1},
251 {-1, 2, 0},
252 { 0,-2, 1},
253 { 1, 0,-2},
254 {-2, 1, 0},
255 { 2, 1, 1},
256 { 1, 2, 1},
257 { 1, 1, 2},
258 { 2,-1,-1},
259 {-1, 2,-1}, /* 30 */
260 {-1,-1, 2},
261 { 2, 1,-1},
262 {-1, 2, 1},
263 { 1,-1, 2},
264 { 2,-1, 1}, /* 35 */
265 { 1, 2,-1},
266 {-1, 1, 2},
267 { 3, 1, 2},
268 { 2, 3, 1},
269 { 1, 2, 3}, /* 40 */
270 { 3, 2, 1},
271 { 1, 3, 2},
272 { 2, 1, 3},
273 { 3,-1, 2},
274 { 2, 3,-1}, /* 45 */
275 {-1, 2, 3},
276 { 3,-2, 1},
277 { 1, 3,-2},
278 {-2, 1, 3},
279 { 3,-1,-2}, /* 50 */
280 {-2, 3,-1},
281 {-1,-2, 3},
282 { 3,-2,-1},
283 {-1, 3,-2},
284 {-2,-1, 3}, /* 55 */
285 { 3, 1,-2},
286 {-2, 3, 1},
287 { 1,-2, 3},
288 { 3, 2,-1},
289 {-1, 3, 2}, /* 60 */
290 { 2,-1, 3},
291 { 1, 1, 3},
292 {-1, 1, 3},
293 { 1,-1, 3},
294 {-1,-1, 3}, /* 65 */
295 { 1, 3, 1},
296 {-1, 3, 1},
297 { 1, 3,-1},
298 {-1, 3,-1},
299 { 3, 1, 1}, /* 70 */
300 { 3, 1,-1},
301 { 3,-1, 1},
302 { 3,-1,-1},
303 };
304
305 static int get_pointgroup_number( const Symmetry * symmetry );
306 static int get_pointgroup_class_table( int table[10],
307 const Symmetry * symmetry );
308 static int get_rotation_type( SPGCONST int rot[3][3] );
309 static int get_rotation_axis( SPGCONST int rot[3][3] );
310 static int get_orthogonal_axis( int ortho_axes[],
311 SPGCONST int proper_rot[3][3],
312 const int rot_order );
313 static int laue2m( int axes[3],
314 const Symmetry * symmetry );
315
316 #ifdef DEBUG
317 static int lauemmm( int axes[3],
318 const Symmetry * symmetry );
319 static int laue4m( int axes[3],
320 const Symmetry * symmetry );
321 static int laue4mmm( int axes[3],
322 const Symmetry * symmetry );
323 static int laue3( int axes[3],
324 const Symmetry * symmetry );
325 static int laue3m( int axes[3],
326 const Symmetry * symmetry );
327 static int lauem3m( int axes[3],
328 const Symmetry * symmetry );
329 #endif
330
331 static int laue_one_axis( int axes[3],
332 const Symmetry * symmetry,
333 const int rot_order );
334 static int lauennn( int axes[3],
335 const Symmetry * symmetry,
336 const int rot_order );
337 static int get_axes( int axes[3],
338 const Laue laue,
339 const Symmetry * symmetry );
340 static void get_proper_rotation( int prop_rot[3][3],
341 SPGCONST int rot[3][3] );
342 static void get_transform_matrix( int mat[3][3],
343 const int axes[3] );
344 static int is_exist_axis( const int axis_vec[3], const int axis_index );
345
346
ptg_get_pointgroup_number(const Symmetry * symmetry)347 int ptg_get_pointgroup_number( const Symmetry * symmetry )
348 {
349 return get_pointgroup_number( symmetry );
350 }
351
ptg_get_pointgroup(const int pointgroup_number)352 Pointgroup ptg_get_pointgroup( const int pointgroup_number )
353 {
354 Pointgroup pointgroup;
355 PointgroupType pointgroup_type;
356
357 pointgroup_type = pointgroup_data[ pointgroup_number ];
358 strcpy( pointgroup.symbol, pointgroup_type.symbol );
359 pointgroup.holohedry = pointgroup_type.holohedry;
360 pointgroup.laue = pointgroup_type.laue;
361
362 debug_print("Point group: %s\n", pointgroup_type.symbol );
363
364 return pointgroup;
365 }
366
367 /* pointgroup is modified. */
ptg_get_transformation_matrix(Pointgroup * pointgroup,const Symmetry * symmetry)368 void ptg_get_transformation_matrix( Pointgroup * pointgroup,
369 const Symmetry * symmetry )
370 {
371 int axes[3];
372 int transform_mat[3][3];
373
374 get_axes( axes, pointgroup->laue, symmetry );
375 get_transform_matrix( transform_mat, axes );
376 mat_copy_matrix_i3( pointgroup->transform_mat, transform_mat );
377 }
378
get_pointgroup_number(const Symmetry * symmetry)379 static int get_pointgroup_number( const Symmetry * symmetry )
380 {
381 int i, j, pg_num, counter;
382 int table[10];
383 PointgroupType pointgroup_type;
384
385 /* Get list of rotation part of symmetry operations */
386 if ( ! get_pointgroup_class_table( table, symmetry ) ) {
387 pg_num = -1;
388 goto end;
389 }
390
391 pg_num = -1;
392 for ( i = 0; i < 32; i++ ) {
393 counter = 0;
394 pointgroup_type = pointgroup_data[ i ];
395 for ( j = 0; j < 10; j++ ) {
396 if ( pointgroup_type.table[j] == table[j] ) { counter++; }
397 }
398 if ( counter == 10 ) {
399 pg_num = i;
400 break;
401 }
402 }
403
404 end:
405 return pg_num;
406 }
407
get_pointgroup_class_table(int table[10],const Symmetry * symmetry)408 static int get_pointgroup_class_table( int table[10],
409 const Symmetry * symmetry )
410 {
411 /* Look-up table */
412 /* Operation -6 -4 -3 -2 -1 1 2 3 4 6 */
413 /* Trace - 2 -1 0 1 -3 3 -1 0 1 2 */
414 /* Determinant -1 -1 -1 -1 -1 1 1 1 1 1 */
415
416 /* table[0] = -6 axis */
417 /* table[1] = -4 axis */
418 /* table[2] = -3 axis */
419 /* table[3] = -2 axis */
420 /* table[4] = -1 axis */
421 /* table[5] = 1 axis */
422 /* table[6] = 2 axis */
423 /* table[7] = 3 axis */
424 /* table[8] = 4 axis */
425 /* table[9] = 6 axis */
426
427 int i, rot_type;
428
429 for ( i = 0; i < 10; i++ ) { table[i] = 0; }
430 for ( i = 0; i < symmetry->size; i++ ) {
431 rot_type = get_rotation_type( symmetry->rot[i] );
432 if ( rot_type == -1 ) {
433 goto err;
434 } else {
435 table[rot_type]++;
436 }
437 }
438
439 return 1;
440
441 err:
442 warning_print("spglib: No point group symbol found ");
443 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
444 return 0;
445 }
446
get_rotation_type(SPGCONST int rot[3][3])447 static int get_rotation_type( SPGCONST int rot[3][3] )
448 {
449 int rot_type;
450
451 if ( mat_get_determinant_i3( rot ) == -1 ) {
452 switch ( mat_get_trace_i3( rot ) ) {
453 case -2: /* -6 */
454 rot_type = 0;
455 break;
456 case -1: /* -4 */
457 rot_type = 1;
458 break;
459 case 0: /* -3 */
460 rot_type = 2;
461 break;
462 case 1: /* -2 */
463 rot_type = 3;
464 break;
465 case -3: /* -1 */
466 rot_type = 4;
467 break;
468 default:
469 rot_type = -1;
470 break;
471 }
472 } else {
473 switch ( mat_get_trace_i3( rot ) ) {
474 case 3: /* 1 */
475 rot_type = 5;
476 break;
477 case -1: /* 2 */
478 rot_type = 6;
479 break;
480 case 0: /* 3 */
481 rot_type = 7;
482 break;
483 case 1: /* 4 */
484 rot_type = 8;
485 break;
486 case 2: /* 6 */
487 rot_type = 9;
488 break;
489 default:
490 rot_type = -1;
491 break;
492 }
493 }
494
495 return rot_type;
496 }
497
get_axes(int axes[3],const Laue laue,const Symmetry * symmetry)498 static int get_axes( int axes[3],
499 const Laue laue,
500 const Symmetry * symmetry )
501 {
502 switch (laue) {
503 case LAUE1:
504 axes[0] = 0;
505 axes[1] = 1;
506 axes[2] = 2;
507 break;
508 case LAUE2M:
509 laue2m( axes, symmetry );
510 break;
511 case LAUEMMM:
512 lauennn( axes, symmetry, 2 );
513 break;
514 case LAUE4M:
515 laue_one_axis( axes, symmetry, 4 );
516 break;
517 case LAUE4MMM:
518 laue_one_axis( axes, symmetry, 4 );
519 break;
520 case LAUE3:
521 laue_one_axis( axes, symmetry, 3 );
522 break;
523 case LAUE3M:
524 laue_one_axis( axes, symmetry, 3 );
525 break;
526 case LAUE6M:
527 laue_one_axis( axes, symmetry, 3 );
528 break;
529 case LAUE6MMM:
530 laue_one_axis( axes, symmetry, 3 );
531 break;
532 case LAUEM3:
533 lauennn( axes, symmetry, 2 );
534 break;
535 case LAUEM3M:
536 lauennn( axes, symmetry, 4 );
537 break;
538 default:
539 break;
540 }
541
542 return 1;
543 }
544
laue2m(int axes[3],const Symmetry * symmetry)545 static int laue2m( int axes[3],
546 const Symmetry * symmetry )
547 {
548 int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
549 int prop_rot[3][3], t_mat[3][3];
550 int ortho_axes[NUM_ROT_AXES];
551
552 for ( i = 0; i < symmetry->size; i++ ) {
553 get_proper_rotation( prop_rot, symmetry->rot[i] );
554
555 /* Search two-fold rotation */
556 if ( ! ( mat_get_trace_i3( prop_rot ) == -1 ) ) {
557 continue;
558 }
559
560 /* The first axis */
561 axes[1] = get_rotation_axis( prop_rot );
562 break;
563 }
564
565 /* The second axis */
566 num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, 2 );
567 if ( ! num_ortho_axis ) { goto err; }
568
569 min_norm = 8;
570 is_found = 0;
571 for ( i = 0; i < num_ortho_axis; i++ ) {
572 norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
573 if ( norm < min_norm ) {
574 min_norm = norm;
575 axes[0] = ortho_axes[i];
576 is_found = 1;
577 }
578 }
579 if ( ! is_found ) { goto err; }
580
581 /* The third axis */
582 min_norm = 8;
583 is_found = 0;
584 for ( i = 0; i < num_ortho_axis; i++ ) {
585 norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
586 if ( norm < min_norm && ( ! ( ortho_axes[i] == axes[0] ) ) ) {
587 min_norm = norm;
588 axes[2] = ortho_axes[i];
589 is_found = 1;
590 }
591 }
592 if ( ! is_found ) { goto err; }
593
594 get_transform_matrix( t_mat, axes );
595 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
596 tmpval = axes[0];
597 axes[0] = axes[2];
598 axes[2] = tmpval;
599 }
600
601 return 1;
602
603 err:
604 return 0;
605 }
606
607 #ifdef DEBUG
lauemmm(int axes[3],const Symmetry * symmetry)608 static int lauemmm( int axes[3],
609 const Symmetry * symmetry )
610 {
611 int i, count, axis, tmpval;
612 int prop_rot[3][3], t_mat[3][3];
613
614 for ( i = 0; i < 3; i++ ) { axes[i] = -1; }
615 count = 0;
616 for ( i = 0; i < symmetry->size; i++ ) {
617 get_proper_rotation( prop_rot, symmetry->rot[i] );
618
619 /* Search two-fold rotation */
620 if ( ! ( mat_get_trace_i3( prop_rot ) == -1 ) ) {
621 continue;
622 }
623
624 axis = get_rotation_axis( prop_rot );
625 if ( ! ( ( axis == axes[0] ) ||
626 ( axis == axes[1] ) ||
627 ( axis == axes[2] ) ) ) {
628 axes[count] = axis;
629 count++;
630 }
631 }
632
633 get_transform_matrix( t_mat, axes );
634 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
635 tmpval = axes[0];
636 axes[0] = axes[1];
637 axes[1] = tmpval;
638 }
639
640 return 1;
641 }
642
laue4m(int axes[3],const Symmetry * symmetry)643 static int laue4m( int axes[3],
644 const Symmetry * symmetry )
645 {
646 int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
647 int axis_vec[3];
648 int prop_rot[3][3], t_mat[3][3];
649 int ortho_axes[NUM_ROT_AXES];
650
651 for ( i = 0; i < symmetry->size; i++ ) {
652 get_proper_rotation( prop_rot, symmetry->rot[i] );
653
654 /* Search foud-fold rotation */
655 if ( mat_get_trace_i3( prop_rot ) == 1 ) {
656 /* The first axis */
657 axes[2] = get_rotation_axis( prop_rot );
658 break;
659 }
660 }
661
662 /* The second axis */
663 num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, 4 );
664 if ( ! num_ortho_axis ) { goto err; }
665
666 min_norm = 8;
667 is_found = 0;
668 for ( i = 0; i < num_ortho_axis; i++ ) {
669 norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
670 if ( norm < min_norm ) {
671 min_norm = norm;
672 axes[0] = ortho_axes[i];
673 is_found = 1;
674 }
675 }
676 if ( ! is_found ) { goto err; }
677
678 /* The third axis */
679 mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
680 is_found = 0;
681 for ( i = 0; i < NUM_ROT_AXES; i++ ) {
682 if ( is_exist_axis( axis_vec, i ) ) {
683 is_found = 1;
684 axes[1] = i;
685 break;
686 }
687 }
688 if ( ! is_found ) { goto err; }
689
690 get_transform_matrix( t_mat, axes );
691 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
692 tmpval = axes[0];
693 axes[0] = axes[1];
694 axes[1] = tmpval;
695 }
696
697 return 1;
698
699 err:
700 return 0;
701 }
702
laue4mmm(int axes[3],const Symmetry * symmetry)703 static int laue4mmm( int axes[3],
704 const Symmetry * symmetry )
705 {
706 int i, is_found, tmpval, axis;
707 int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
708 int axis_vec[3];
709
710 for ( i = 0; i < symmetry->size; i++ ) {
711 get_proper_rotation( prop_rot, symmetry->rot[i] );
712
713 /* Search foud-fold rotation */
714 if ( mat_get_trace_i3( prop_rot ) == 1 ) {
715 /* The first axis */
716 axes[2] = get_rotation_axis( prop_rot );
717 break;
718 }
719 }
720
721 is_found = 0;
722 for ( i = 0; i < symmetry->size; i++ ) {
723 get_proper_rotation( prop_rot2, symmetry->rot[i] );
724
725 /* Search two-fold rotation */
726 if ( ! ( mat_get_trace_i3( prop_rot2 ) == -1 ) ) {
727 continue;
728 }
729
730 /* The second axis */
731 axis = get_rotation_axis( prop_rot2 );
732 if ( ! ( axis == axes[2] ) ) {
733 axes[0] = axis;
734 is_found = 1;
735 break;
736 }
737 }
738 if ( ! is_found ) { goto err; }
739
740
741 /* The third axis */
742 mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
743 is_found = 0;
744 for ( i = 0; i < NUM_ROT_AXES; i++ ) {
745 if ( is_exist_axis( axis_vec, i ) ) {
746 is_found = 1;
747 axes[1] = i;
748 break;
749 }
750 }
751 if ( ! is_found ) { goto err; }
752
753 get_transform_matrix( t_mat, axes );
754 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
755 tmpval = axes[0];
756 axes[0] = axes[1];
757 axes[1] = tmpval;
758 }
759
760 return 1;
761
762 err:
763 return 0;
764 }
765
laue3(int axes[3],const Symmetry * symmetry)766 static int laue3( int axes[3],
767 const Symmetry * symmetry )
768 {
769 int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
770 int prop_rot[3][3], t_mat[3][3];
771 int axis_vec[3];
772 int ortho_axes[NUM_ROT_AXES];
773
774 for ( i = 0; i < symmetry->size; i++ ) {
775 get_proper_rotation( prop_rot, symmetry->rot[i] );
776
777 /* Search thee-fold rotation */
778 if ( mat_get_trace_i3( prop_rot ) == 0 ) {
779 /* The first axis */
780 axes[2] = get_rotation_axis( prop_rot );
781 break;
782 }
783 }
784
785 /* The second axis */
786 num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, 3 );
787 if ( ! num_ortho_axis ) { goto err; }
788 min_norm = 8;
789 is_found = 0;
790 for ( i = 0; i < num_ortho_axis; i++ ) {
791 norm = mat_norm_squared_i3( rot_axes[ortho_axes[i]] );
792 if ( norm < min_norm ) {
793 min_norm = norm;
794 axes[0] = ortho_axes[i];
795 is_found = 1;
796 }
797 }
798 if ( ! is_found ) { goto err; }
799
800 /* The third axis */
801 mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
802 is_found = 0;
803 for ( i = 0; i < NUM_ROT_AXES; i++ ) {
804 is_found = is_exist_axis( axis_vec, i );
805 if ( is_found == 1 ) {
806 axes[1] = i;
807 break;
808 }
809 if ( is_found == -1 ) {
810 axes[1] = i + NUM_ROT_AXES;
811 break;
812 }
813 }
814 if ( ! is_found ) { goto err; }
815
816 get_transform_matrix( t_mat, axes );
817 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
818 tmpval = axes[0];
819 axes[0] = axes[1];
820 axes[1] = tmpval;
821 }
822
823 return 1;
824
825 err:
826 return 0;
827 }
828
laue3m(int axes[3],const Symmetry * symmetry)829 static int laue3m( int axes[3],
830 const Symmetry * symmetry )
831 {
832 int i, is_found, tmpval, axis;
833 int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
834 int axis_vec[3];
835
836 for ( i = 0; i < symmetry->size; i++ ) {
837 get_proper_rotation( prop_rot, symmetry->rot[i] );
838
839 /* Search three-fold rotation */
840 if ( mat_get_trace_i3( prop_rot ) == 0 ) {
841 /* The first axis */
842 axes[2] = get_rotation_axis( prop_rot );
843 debug_print("laue3m prop_rot\n");
844 debug_print_matrix_i3( prop_rot );
845 break;
846 }
847 }
848
849 is_found = 0;
850 for ( i = 0; i < symmetry->size; i++ ) {
851 get_proper_rotation( prop_rot2, symmetry->rot[i] );
852
853 /* Search two-fold rotation */
854 if ( ! ( mat_get_trace_i3( prop_rot2 ) == -1 ) ) {
855 continue;
856 }
857
858 /* The second axis */
859 axis = get_rotation_axis( prop_rot2 );
860 if ( ! ( axis == axes[2] ) ) {
861 axes[0] = axis;
862 is_found = 1;
863 break;
864 }
865 }
866 if ( ! is_found ) { goto err; }
867
868 /* The third axis */
869 mat_multiply_matrix_vector_i3( axis_vec, prop_rot, rot_axes[axes[0]] );
870 is_found = 0;
871 for ( i = 0; i < NUM_ROT_AXES; i++ ) {
872 is_found = is_exist_axis( axis_vec, i );
873 if ( is_found == 1 ) {
874 axes[1] = i;
875 break;
876 }
877 if ( is_found == -1 ) {
878 axes[1] = i + NUM_ROT_AXES;
879 break;
880 }
881 }
882 if ( ! is_found ) { goto err; }
883
884 get_transform_matrix( t_mat, axes );
885 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
886 tmpval = axes[0];
887 axes[0] = axes[1];
888 axes[1] = tmpval;
889 }
890
891 return 1;
892
893 err:
894 return 0;
895 }
896
lauem3m(int axes[3],const Symmetry * symmetry)897 static int lauem3m( int axes[3],
898 const Symmetry * symmetry )
899 {
900 int i, count, axis, tmpval;
901 int prop_rot[3][3], t_mat[3][3];
902
903 for ( i = 0; i < 3; i++ ) { axes[i] = -1; }
904 count = 0;
905 for ( i = 0; i < symmetry->size; i++ ) {
906 get_proper_rotation( prop_rot, symmetry->rot[i] );
907
908 /* Search four-fold rotation */
909 if ( ! ( mat_get_trace_i3( prop_rot ) == 1 ) ) {
910 continue;
911 }
912
913 axis = get_rotation_axis( prop_rot );
914 if ( ! ( ( axis == axes[0] ) ||
915 ( axis == axes[1] ) ||
916 ( axis == axes[2] ) ) ) {
917 axes[count] = axis;
918 count++;
919 }
920 }
921
922 get_transform_matrix( t_mat, axes );
923 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
924 tmpval = axes[0];
925 axes[0] = axes[1];
926 axes[1] = tmpval;
927 }
928
929 return 1;
930 }
931 #endif
932
laue_one_axis(int axes[3],const Symmetry * symmetry,const int rot_order)933 static int laue_one_axis( int axes[3],
934 const Symmetry * symmetry,
935 const int rot_order )
936 {
937 int i, j, num_ortho_axis, det, min_det, is_found, tmpval;
938 int axis_vec[3], tmp_axes[3];
939 int prop_rot[3][3], t_mat[3][3];
940 int ortho_axes[NUM_ROT_AXES];
941
942 for ( i = 0; i < symmetry->size; i++ ) {
943 get_proper_rotation( prop_rot, symmetry->rot[i] );
944
945 /* Search foud-fold rotation */
946 if ( rot_order == 4 ) {
947 if ( mat_get_trace_i3( prop_rot ) == 1 ) {
948 /* The first axis */
949 axes[2] = get_rotation_axis( prop_rot );
950 break;
951 }
952 }
953
954 /* Search three-fold rotation */
955 if ( rot_order == 3 ) {
956 if ( mat_get_trace_i3( prop_rot ) == 0 ) {
957 /* The first axis */
958 axes[2] = get_rotation_axis( prop_rot );
959 break;
960 }
961 }
962 }
963
964 /* Candidates of the second axis */
965 num_ortho_axis = get_orthogonal_axis( ortho_axes, prop_rot, rot_order );
966 if ( ! num_ortho_axis ) { goto err; }
967
968 tmp_axes[2] = axes[2];
969 min_det = 4;
970 is_found = 0;
971 for ( i = 0; i < num_ortho_axis; i++ ) {
972 tmp_axes[0] = ortho_axes[i];
973 mat_multiply_matrix_vector_i3( axis_vec, prop_rot,
974 rot_axes[tmp_axes[0]] );
975 for ( j = 0; j < num_ortho_axis; j++ ) {
976 is_found = is_exist_axis( axis_vec, ortho_axes[j] );
977 if ( is_found == 1 ) {
978 tmp_axes[1] = ortho_axes[j];
979 break;
980 }
981 if ( is_found == -1 ) {
982 tmp_axes[1] = ortho_axes[j] + NUM_ROT_AXES;
983 break;
984 }
985 }
986
987 get_transform_matrix( t_mat, tmp_axes );
988 det = mat_get_determinant_i3( t_mat );
989 if ( det < 0 ) { det = -det; }
990 if ( det < min_det ) {
991 min_det = det;
992 axes[0] = tmp_axes[0];
993 axes[1] = tmp_axes[1];
994 }
995 }
996 if ( ! is_found ) { goto err; }
997
998 get_transform_matrix( t_mat, axes );
999 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
1000 tmpval = axes[0];
1001 axes[0] = axes[1];
1002 axes[1] = tmpval;
1003 }
1004
1005 return 1;
1006
1007 err:
1008 return 0;
1009 }
1010
lauennn(int axes[3],const Symmetry * symmetry,const int rot_order)1011 static int lauennn( int axes[3],
1012 const Symmetry * symmetry,
1013 const int rot_order )
1014 {
1015 int i, count, axis, tmpval;
1016 int prop_rot[3][3], t_mat[3][3];
1017
1018 for ( i = 0; i < 3; i++ ) { axes[i] = -1; }
1019 count = 0;
1020 for ( i = 0; i < symmetry->size; i++ ) {
1021 get_proper_rotation( prop_rot, symmetry->rot[i] );
1022
1023 /* Search two- or four-fold rotation */
1024 if ( ( mat_get_trace_i3( prop_rot ) == -1 && rot_order == 2) ||
1025 ( mat_get_trace_i3( prop_rot ) == 1 && rot_order == 4 ) ) {
1026 axis = get_rotation_axis( prop_rot );
1027 if ( ! ( ( axis == axes[0] ) ||
1028 ( axis == axes[1] ) ||
1029 ( axis == axes[2] ) ) ) {
1030 axes[count] = axis;
1031 count++;
1032 }
1033 }
1034 }
1035
1036 get_transform_matrix( t_mat, axes );
1037 if ( mat_get_determinant_i3( t_mat ) < 0 ) {
1038 tmpval = axes[0];
1039 axes[0] = axes[1];
1040 axes[1] = tmpval;
1041 }
1042
1043 return 1;
1044 }
1045
get_rotation_axis(SPGCONST int proper_rot[3][3])1046 static int get_rotation_axis( SPGCONST int proper_rot[3][3] )
1047 {
1048 int i, axis = -1;
1049 int vec[3];
1050
1051 /* No specific axis for I and -I */
1052 if ( mat_check_identity_matrix_i3( proper_rot, identity ) ) {
1053 goto end;
1054 }
1055
1056 /* Look for eigenvector = rotation axis */
1057 for ( i = 0; i < NUM_ROT_AXES; i++ ) {
1058 mat_multiply_matrix_vector_i3( vec, proper_rot, rot_axes[i] );
1059 if ( vec[0] == rot_axes[i][0] &&
1060 vec[1] == rot_axes[i][1] &&
1061 vec[2] == rot_axes[i][2] ) {
1062 axis = i;
1063 break;
1064 }
1065 }
1066
1067 end:
1068 #ifdef DEBUG
1069 if ( axis == -1 ) {
1070 printf("rotation axis cound not found.\n");
1071 }
1072 #endif
1073
1074 return axis;
1075 }
1076
get_orthogonal_axis(int ortho_axes[],SPGCONST int proper_rot[3][3],const int rot_order)1077 static int get_orthogonal_axis( int ortho_axes[],
1078 SPGCONST int proper_rot[3][3],
1079 const int rot_order )
1080 {
1081 int i, num_ortho_axis;
1082 int vec[3];
1083 int sum_rot[3][3], rot[3][3];
1084
1085 num_ortho_axis = 0;
1086
1087 mat_copy_matrix_i3( sum_rot, identity );
1088 mat_copy_matrix_i3( rot, identity );
1089 for ( i = 0; i < rot_order-1; i++ ) {
1090 mat_multiply_matrix_i3( rot, proper_rot, rot );
1091 mat_add_matrix_i3( sum_rot, rot, sum_rot );
1092 }
1093
1094 for ( i = 0; i < NUM_ROT_AXES; i++ ) {
1095 mat_multiply_matrix_vector_i3( vec, sum_rot, rot_axes[i] );
1096 if ( vec[0] == 0 &&
1097 vec[1] == 0 &&
1098 vec[2] == 0 ) {
1099 ortho_axes[ num_ortho_axis ] = i;
1100 num_ortho_axis++;
1101 }
1102 }
1103
1104 return num_ortho_axis;
1105 }
1106
get_proper_rotation(int prop_rot[3][3],SPGCONST int rot[3][3])1107 static void get_proper_rotation( int prop_rot[3][3],
1108 SPGCONST int rot[3][3] )
1109 {
1110 if ( mat_get_determinant_i3( rot ) == -1 ) {
1111 mat_multiply_matrix_i3( prop_rot, inversion, rot );
1112 } else {
1113 mat_copy_matrix_i3( prop_rot, rot );
1114 }
1115 }
1116
get_transform_matrix(int mat[3][3],const int axes[3])1117 static void get_transform_matrix( int mat[3][3],
1118 const int axes[3] )
1119 {
1120 int i, j, s[3];
1121
1122 for ( i = 0; i < 3; i++ ) {
1123 if ( axes[i] < NUM_ROT_AXES ) {
1124 s[i] = 1;
1125 } else {
1126 s[i] = -1;
1127 }
1128 }
1129 for ( i = 0; i < 3; i++ ) {
1130 for ( j = 0; j < 3; j++ ) {
1131 mat[i][j] = s[j] * rot_axes[axes[j]%NUM_ROT_AXES][i];
1132 }
1133 }
1134 }
1135
is_exist_axis(const int axis_vec[3],const int axis_index)1136 static int is_exist_axis( const int axis_vec[3], const int axis_index )
1137 {
1138 if ( ( axis_vec[0] == rot_axes[axis_index][0] ) &&
1139 ( axis_vec[1] == rot_axes[axis_index][1] ) &&
1140 ( axis_vec[2] == rot_axes[axis_index][2] ) ) { return 1; }
1141 if ( ( axis_vec[0] == -rot_axes[axis_index][0] ) &&
1142 ( axis_vec[1] == -rot_axes[axis_index][1] ) &&
1143 ( axis_vec[2] == -rot_axes[axis_index][2] ) ) { return -1; }
1144 return 0;
1145 }
1146
1147