1 /* Copyright (C) 2008 Atsushi Togo */
2 /* All rights reserved. */
3
4 /* This file is part of spglib. */
5
6 /* Redistribution and use in source and binary forms, with or without */
7 /* modification, are permitted provided that the following conditions */
8 /* are met: */
9
10 /* * Redistributions of source code must retain the above copyright */
11 /* notice, this list of conditions and the following disclaimer. */
12
13 /* * Redistributions in binary form must reproduce the above copyright */
14 /* notice, this list of conditions and the following disclaimer in */
15 /* the documentation and/or other materials provided with the */
16 /* distribution. */
17
18 /* * Neither the name of the phonopy project nor the names of its */
19 /* contributors may be used to endorse or promote products derived */
20 /* from this software without specific prior written permission. */
21
22 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
23 /* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
24 /* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
25 /* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
26 /* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
27 /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
28 /* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
29 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
30 /* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
31 /* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
32 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
33 /* POSSIBILITY OF SUCH DAMAGE. */
34
35 #include <string.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include "pointgroup.h"
39 #include "symmetry.h"
40 #include "mathfunc.h"
41
42 #include "debug.h"
43
44 #define NUM_ROT_AXES 73
45
46 typedef struct {
47 int table[10];
48 char symbol[6];
49 char schoenflies[4];
50 Holohedry holohedry;
51 Laue laue;
52 } PointgroupType;
53
54 static PointgroupType pointgroup_data[33] = {
55 { /* 0 */
56 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
57 " ",
58 " ",
59 HOLOHEDRY_NONE,
60 LAUE_NONE,
61 },
62 { /* 1 */
63 {0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
64 "1 ",
65 "C1 ",
66 TRICLI,
67 LAUE1,
68 },
69 { /* 2 */
70 {0, 0, 0, 0, 1, 1, 0, 0, 0, 0},
71 "-1 ",
72 "Ci ",
73 TRICLI,
74 LAUE1,
75 },
76 { /* 3 */
77 {0, 0, 0, 0, 0, 1, 1, 0, 0, 0},
78 "2 ",
79 "C2 ",
80 MONOCLI,
81 LAUE2M,
82 },
83 { /* 4 */
84 {0, 0, 0, 1, 0, 1, 0, 0, 0, 0},
85 "m ",
86 "Cs ",
87 MONOCLI,
88 LAUE2M,
89 },
90 { /* 5 */
91 {0, 0, 0, 1, 1, 1, 1, 0, 0, 0},
92 "2/m ",
93 "C2h",
94 MONOCLI,
95 LAUE2M,
96 },
97 { /* 6 */
98 {0, 0, 0, 0, 0, 1, 3, 0, 0, 0},
99 "222 ",
100 "D2 ",
101 ORTHO,
102 LAUEMMM,
103 },
104 { /* 7 */
105 {0, 0, 0, 2, 0, 1, 1, 0, 0, 0},
106 "mm2 ",
107 "C2v",
108 ORTHO,
109 LAUEMMM,
110 },
111 { /* 8 */
112 {0, 0, 0, 3, 1, 1, 3, 0, 0, 0},
113 "mmm ",
114 "D2h",
115 ORTHO,
116 LAUEMMM,
117 },
118 { /* 9 */
119 {0, 0, 0, 0, 0, 1, 1, 0, 2, 0},
120 "4 ",
121 "C4 ",
122 TETRA,
123 LAUE4M,
124 },
125 { /* 10 */
126 {0, 2, 0, 0, 0, 1, 1, 0, 0, 0},
127 "-4 ",
128 "S4 ",
129 TETRA,
130 LAUE4M,
131 },
132 { /* 11 */
133 {0, 2, 0, 1, 1, 1, 1, 0, 2, 0},
134 "4/m ",
135 "C4h",
136 TETRA,
137 LAUE4M,
138 },
139 { /* 12 */
140 {0, 0, 0, 0, 0, 1, 5, 0, 2, 0},
141 "422 ",
142 "D4 ",
143 TETRA,
144 LAUE4MMM,
145 },
146 { /* 13 */
147 {0, 0, 0, 4, 0, 1, 1, 0, 2, 0},
148 "4mm ",
149 "C4v",
150 TETRA,
151 LAUE4MMM,
152 },
153 { /* 14 */
154 {0, 2, 0, 2, 0, 1, 3, 0, 0, 0},
155 "-42m ",
156 "D2d",
157 TETRA,
158 LAUE4MMM,
159 },
160 { /* 15 */
161 {0, 2, 0, 5, 1, 1, 5, 0, 2, 0},
162 "4/mmm",
163 "D4h",
164 TETRA,
165 LAUE4MMM,
166 },
167 { /* 16 */
168 {0, 0, 0, 0, 0, 1, 0, 2, 0, 0},
169 "3 ",
170 "C3 ",
171 TRIGO,
172 LAUE3,
173 },
174 { /* 17 */
175 {0, 0, 2, 0, 1, 1, 0, 2, 0, 0},
176 "-3 ",
177 "C3i",
178 TRIGO,
179 LAUE3,
180 },
181 { /* 18 */
182 {0, 0, 0, 0, 0, 1, 3, 2, 0, 0},
183 "32 ",
184 "D3 ",
185 TRIGO,
186 LAUE3M,
187 },
188 { /* 19 */
189 {0, 0, 0, 3, 0, 1, 0, 2, 0, 0},
190 "3m ",
191 "C3v",
192 TRIGO,
193 LAUE3M,
194 },
195 { /* 20 */
196 {0, 0, 2, 3, 1, 1, 3, 2, 0, 0},
197 "-3m ",
198 "D3d",
199 TRIGO,
200 LAUE3M,
201 },
202 { /* 21 */
203 {0, 0, 0, 0, 0, 1, 1, 2, 0, 2},
204 "6 ",
205 "C6 ",
206 HEXA,
207 LAUE6M,
208 },
209 { /* 22 */
210 {2, 0, 0, 1, 0, 1, 0, 2, 0, 0},
211 "-6 ",
212 "C3h",
213 HEXA,
214 LAUE6M,
215 },
216 { /* 23 */
217 {2, 0, 2, 1, 1, 1, 1, 2, 0, 2},
218 "6/m ",
219 "C6h",
220 HEXA,
221 LAUE6M,
222 },
223 { /* 24 */
224 {0, 0, 0, 0, 0, 1, 7, 2, 0, 2},
225 "622 ",
226 "D6 ",
227 HEXA,
228 LAUE6MMM,
229 },
230 { /* 25 */
231 {0, 0, 0, 6, 0, 1, 1, 2, 0, 2},
232 "6mm ",
233 "C6v",
234 HEXA,
235 LAUE6MMM,
236 },
237 { /* 26 */
238 {2, 0, 0, 4, 0, 1, 3, 2, 0, 0},
239 "-6m2 ",
240 "D3h",
241 HEXA,
242 LAUE6MMM,
243 },
244 { /* 27 */
245 {2, 0, 2, 7, 1, 1, 7, 2, 0, 2},
246 "6/mmm",
247 "D6h",
248 HEXA,
249 LAUE6MMM,
250 },
251 { /* 28 */
252 {0, 0, 0, 0, 0, 1, 3, 8, 0, 0},
253 "23 ",
254 "T ",
255 CUBIC,
256 LAUEM3,
257 },
258 { /* 29 */
259 {0, 0, 8, 3, 1, 1, 3, 8, 0, 0},
260 "m-3 ",
261 "Th ",
262 CUBIC,
263 LAUEM3,
264 },
265 { /* 30 */
266 {0, 0, 0, 0, 0, 1, 9, 8, 6, 0},
267 "432 ",
268 "O ",
269 CUBIC,
270 LAUEM3M,
271 },
272 { /* 31 */
273 {0, 6, 0, 6, 0, 1, 3, 8, 0, 0},
274 "-43m ",
275 "Td ",
276 CUBIC,
277 LAUEM3M,
278 },
279 { /* 32 */
280 {0, 6, 8, 9, 1, 1, 9, 8, 6, 0},
281 "m-3m ",
282 "Oh ",
283 CUBIC,
284 LAUEM3M,
285 }
286 };
287
288 static int identity[3][3] = {
289 { 1, 0, 0},
290 { 0, 1, 0},
291 { 0, 0, 1},
292 };
293
294 static int inversion[3][3] = {
295 {-1, 0, 0},
296 { 0,-1, 0},
297 { 0, 0,-1},
298 };
299
300 static int rot_axes[][3] = {
301 { 1, 0, 0},
302 { 0, 1, 0},
303 { 0, 0, 1},
304 { 0, 1, 1},
305 { 1, 0, 1},
306 { 1, 1, 0},
307 { 0, 1,-1},
308 {-1, 0, 1},
309 { 1,-1, 0},
310 { 1, 1, 1}, /* 10 */
311 {-1, 1, 1},
312 { 1,-1, 1},
313 { 1, 1,-1},
314 { 0, 1, 2},
315 { 2, 0, 1},
316 { 1, 2, 0},
317 { 0, 2, 1},
318 { 1, 0, 2},
319 { 2, 1, 0},
320 { 0,-1, 2}, /* 20 */
321 { 2, 0,-1},
322 {-1, 2, 0},
323 { 0,-2, 1},
324 { 1, 0,-2},
325 {-2, 1, 0},
326 { 2, 1, 1},
327 { 1, 2, 1},
328 { 1, 1, 2},
329 { 2,-1,-1},
330 {-1, 2,-1}, /* 30 */
331 {-1,-1, 2},
332 { 2, 1,-1},
333 {-1, 2, 1},
334 { 1,-1, 2},
335 { 2,-1, 1}, /* 35 */
336 { 1, 2,-1},
337 {-1, 1, 2},
338 { 3, 1, 2},
339 { 2, 3, 1},
340 { 1, 2, 3}, /* 40 */
341 { 3, 2, 1},
342 { 1, 3, 2},
343 { 2, 1, 3},
344 { 3,-1, 2},
345 { 2, 3,-1}, /* 45 */
346 {-1, 2, 3},
347 { 3,-2, 1},
348 { 1, 3,-2},
349 {-2, 1, 3},
350 { 3,-1,-2}, /* 50 */
351 {-2, 3,-1},
352 {-1,-2, 3},
353 { 3,-2,-1},
354 {-1, 3,-2},
355 {-2,-1, 3}, /* 55 */
356 { 3, 1,-2},
357 {-2, 3, 1},
358 { 1,-2, 3},
359 { 3, 2,-1},
360 {-1, 3, 2}, /* 60 */
361 { 2,-1, 3},
362 { 1, 1, 3},
363 {-1, 1, 3},
364 { 1,-1, 3},
365 {-1,-1, 3}, /* 65 */
366 { 1, 3, 1},
367 {-1, 3, 1},
368 { 1, 3,-1},
369 {-1, 3,-1},
370 { 3, 1, 1}, /* 70 */
371 { 3, 1,-1},
372 { 3,-1, 1},
373 { 3,-1,-1},
374 };
375
376 static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],
377 const int num_rotations);
378 static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym);
379 static int get_pointgroup_class_table(int table[10],
380 SPGCONST PointSymmetry * pointsym);
381 static int get_rotation_type(SPGCONST int rot[3][3]);
382 static int get_rotation_axis(SPGCONST int rot[3][3]);
383 static int get_orthogonal_axis(int ortho_axes[],
384 SPGCONST int proper_rot[3][3],
385 const int rot_order);
386 static int laue2m(int axes[3],
387 SPGCONST PointSymmetry * pointsym);
388
389 #ifdef SPGDEBUG
390 static int lauemmm(int axes[3],
391 SPGCONST PointSymmetry * pointsym);
392 static int laue4m(int axes[3],
393 SPGCONST PointSymmetry * pointsym);
394 static int laue4mmm(int axes[3],
395 SPGCONST PointSymmetry * pointsym);
396 static int laue3(int axes[3],
397 SPGCONST PointSymmetry * pointsym);
398 static int laue3m(int axes[3],
399 SPGCONST PointSymmetry * pointsym);
400 static int lauem3m(int axes[3],
401 SPGCONST PointSymmetry * pointsym);
402 #endif
403
404 static int laue_one_axis(int axes[3],
405 SPGCONST PointSymmetry * pointsym,
406 const int rot_order);
407 static int lauennn(int axes[3],
408 SPGCONST PointSymmetry * pointsym,
409 const int rot_order);
410 static int get_axes(int axes[3],
411 const Laue laue,
412 SPGCONST PointSymmetry * pointsym);
413 static void get_proper_rotation(int prop_rot[3][3],
414 SPGCONST int rot[3][3]);
415 static void set_transformation_matrix(int tmat[3][3],
416 const int axes[3]);
417 static int is_exist_axis(const int axis_vec[3], const int axis_index);
418 static void sort_axes(int axes[3]);
419
420 /* Retrun pointgroup.number = 0 if failed */
ptg_get_transformation_matrix(int transform_mat[3][3],SPGCONST int rotations[][3][3],const int num_rotations)421 Pointgroup ptg_get_transformation_matrix(int transform_mat[3][3],
422 SPGCONST int rotations[][3][3],
423 const int num_rotations)
424 {
425 int i, j, pg_num;
426 int axes[3];
427 PointSymmetry pointsym;
428 Pointgroup pointgroup;
429
430 debug_print("ptg_get_transformation_matrix:\n");
431
432 for (i = 0; i < 3; i++) {
433 for (j = 0; j < 3; j++) {
434 transform_mat[i][j] = 0;
435 }
436 }
437
438 pg_num = get_pointgroup_number_by_rotations(rotations, num_rotations);
439
440 if (pg_num > 0) {
441 pointgroup = ptg_get_pointgroup(pg_num);
442 pointsym = ptg_get_pointsymmetry(rotations, num_rotations);
443 get_axes(axes, pointgroup.laue, &pointsym);
444 set_transformation_matrix(transform_mat, axes);
445 } else {
446 pointgroup = ptg_get_pointgroup(0);
447 }
448
449 return pointgroup;
450 }
451
ptg_get_pointgroup(const int pointgroup_number)452 Pointgroup ptg_get_pointgroup(const int pointgroup_number)
453 {
454 int i;
455 Pointgroup pointgroup;
456 PointgroupType pointgroup_type;
457
458 pointgroup.number = pointgroup_number;
459 pointgroup_type = pointgroup_data[pointgroup_number];
460 strcpy(pointgroup.symbol, pointgroup_type.symbol);
461 strcpy(pointgroup.schoenflies, pointgroup_type.schoenflies);
462 for (i = 0; i < 5; i++) {
463 if (pointgroup.symbol[i] == ' ') {pointgroup.symbol[i] = '\0';}
464 }
465 for (i = 0; i < 3; i++) {
466 if (pointgroup.schoenflies[i] == ' ') {pointgroup.schoenflies[i] = '\0';}
467 }
468 pointgroup.holohedry = pointgroup_type.holohedry;
469 pointgroup.laue = pointgroup_type.laue;
470
471 debug_print("ptg_get_pointgroup: %s\n", pointgroup_type.symbol);
472
473 return pointgroup;
474 }
475
ptg_get_pointsymmetry(SPGCONST int rotations[][3][3],const int num_rotations)476 PointSymmetry ptg_get_pointsymmetry(SPGCONST int rotations[][3][3],
477 const int num_rotations)
478 {
479 int i, j;
480 PointSymmetry pointsym;
481
482 pointsym.size = 0;
483 for (i = 0; i < num_rotations; i++) {
484 for (j = 0; j < pointsym.size; j++) {
485 if (mat_check_identity_matrix_i3(rotations[i], pointsym.rot[j])) {
486 goto escape;
487 }
488 }
489 mat_copy_matrix_i3(pointsym.rot[pointsym.size], rotations[i]);
490 pointsym.size++;
491 escape:
492 ;
493 }
494
495 return pointsym;
496 }
497
get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],const int num_rotations)498 static int get_pointgroup_number_by_rotations(SPGCONST int rotations[][3][3],
499 const int num_rotations)
500 {
501 PointSymmetry pointsym;
502
503 pointsym = ptg_get_pointsymmetry(rotations, num_rotations);
504 return get_pointgroup_number(&pointsym);
505 }
506
get_pointgroup_number(SPGCONST PointSymmetry * pointsym)507 static int get_pointgroup_number(SPGCONST PointSymmetry * pointsym)
508 {
509 int i, j, pg_num, counter;
510 int table[10];
511 PointgroupType pointgroup_type;
512
513 debug_print("get_pointgroup_number:");
514
515
516 pg_num = 0;
517
518 /* Get list of point symmetry operations */
519 if (! get_pointgroup_class_table(table, pointsym)) {
520 goto end;
521 }
522
523 for (i = 1; i < 33; i++) {
524 counter = 0;
525 pointgroup_type = pointgroup_data[i];
526 for (j = 0; j < 10; j++) {
527 if (pointgroup_type.table[j] == table[j]) {counter++;}
528 }
529 if (counter == 10) {
530 pg_num = i;
531 break;
532 }
533 }
534
535 end:
536 debug_print(" %d\n", pg_num);
537 return pg_num;
538 }
539
get_pointgroup_class_table(int table[10],SPGCONST PointSymmetry * pointsym)540 static int get_pointgroup_class_table(int table[10],
541 SPGCONST PointSymmetry * pointsym)
542 {
543 /* Look-up table */
544 /* Operation -6 -4 -3 -2 -1 1 2 3 4 6 */
545 /* Trace - 2 -1 0 1 -3 3 -1 0 1 2 */
546 /* Determinant -1 -1 -1 -1 -1 1 1 1 1 1 */
547
548 /* table[0] = -6 axis */
549 /* table[1] = -4 axis */
550 /* table[2] = -3 axis */
551 /* table[3] = -2 axis */
552 /* table[4] = -1 axis */
553 /* table[5] = 1 axis */
554 /* table[6] = 2 axis */
555 /* table[7] = 3 axis */
556 /* table[8] = 4 axis */
557 /* table[9] = 6 axis */
558
559 int i, rot_type;
560
561 for (i = 0; i < 10; i++) { table[i] = 0; }
562 for (i = 0; i < pointsym->size; i++) {
563 rot_type = get_rotation_type(pointsym->rot[i]);
564 if (rot_type == -1) {
565 goto err;
566 } else {
567 table[rot_type]++;
568 }
569 }
570
571 return 1;
572
573 err:
574 warning_print("spglib: No point group symbol found ");
575 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
576 return 0;
577 }
578
get_rotation_type(SPGCONST int rot[3][3])579 static int get_rotation_type(SPGCONST int rot[3][3])
580 {
581 int rot_type;
582
583 if (mat_get_determinant_i3(rot) == -1) {
584 switch (mat_get_trace_i3(rot)) {
585 case -2: /* -6 */
586 rot_type = 0;
587 break;
588 case -1: /* -4 */
589 rot_type = 1;
590 break;
591 case 0: /* -3 */
592 rot_type = 2;
593 break;
594 case 1: /* -2 */
595 rot_type = 3;
596 break;
597 case -3: /* -1 */
598 rot_type = 4;
599 break;
600 default:
601 rot_type = -1;
602 break;
603 }
604 } else {
605 switch (mat_get_trace_i3(rot)) {
606 case 3: /* 1 */
607 rot_type = 5;
608 break;
609 case -1: /* 2 */
610 rot_type = 6;
611 break;
612 case 0: /* 3 */
613 rot_type = 7;
614 break;
615 case 1: /* 4 */
616 rot_type = 8;
617 break;
618 case 2: /* 6 */
619 rot_type = 9;
620 break;
621 default:
622 rot_type = -1;
623 break;
624 }
625 }
626
627 return rot_type;
628 }
629
get_axes(int axes[3],const Laue laue,SPGCONST PointSymmetry * pointsym)630 static int get_axes(int axes[3],
631 const Laue laue,
632 SPGCONST PointSymmetry * pointsym)
633 {
634 switch (laue) {
635 case LAUE1:
636 axes[0] = 0;
637 axes[1] = 1;
638 axes[2] = 2;
639 break;
640 case LAUE2M:
641 laue2m(axes, pointsym);
642 break;
643 case LAUEMMM:
644 lauennn(axes, pointsym, 2);
645 break;
646 case LAUE4M:
647 laue_one_axis(axes, pointsym, 4);
648 break;
649 case LAUE4MMM:
650 laue_one_axis(axes, pointsym, 4);
651 break;
652 case LAUE3:
653 laue_one_axis(axes, pointsym, 3);
654 break;
655 case LAUE3M:
656 laue_one_axis(axes, pointsym, 3);
657 break;
658 case LAUE6M:
659 laue_one_axis(axes, pointsym, 3);
660 break;
661 case LAUE6MMM:
662 laue_one_axis(axes, pointsym, 3);
663 break;
664 case LAUEM3:
665 lauennn(axes, pointsym, 2);
666 break;
667 case LAUEM3M:
668 lauennn(axes, pointsym, 4);
669 break;
670 default:
671 break;
672 }
673
674 return 1;
675 }
676
laue2m(int axes[3],SPGCONST PointSymmetry * pointsym)677 static int laue2m(int axes[3],
678 SPGCONST PointSymmetry * pointsym)
679 {
680 int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
681 int prop_rot[3][3], t_mat[3][3];
682 int ortho_axes[NUM_ROT_AXES];
683
684 for (i = 0; i < pointsym->size; i++) {
685 get_proper_rotation(prop_rot, pointsym->rot[i]);
686
687 /* Search two-fold rotation */
688 if (! (mat_get_trace_i3(prop_rot) == -1)) {
689 continue;
690 }
691
692 /* The first axis */
693 axes[1] = get_rotation_axis(prop_rot);
694 break;
695 }
696
697 /* The second axis */
698 num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, 2);
699 if (! num_ortho_axis) { goto err; }
700
701 min_norm = 8;
702 is_found = 0;
703 for (i = 0; i < num_ortho_axis; i++) {
704 norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
705 if (norm < min_norm) {
706 min_norm = norm;
707 axes[0] = ortho_axes[i];
708 is_found = 1;
709 }
710 }
711 if (! is_found) { goto err; }
712
713 /* The third axis */
714 min_norm = 8;
715 is_found = 0;
716 for (i = 0; i < num_ortho_axis; i++) {
717 norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
718 if (norm < min_norm && (! (ortho_axes[i] == axes[0]))) {
719 min_norm = norm;
720 axes[2] = ortho_axes[i];
721 is_found = 1;
722 }
723 }
724 if (! is_found) { goto err; }
725
726 set_transformation_matrix(t_mat, axes);
727 if (mat_get_determinant_i3(t_mat) < 0) {
728 tmpval = axes[0];
729 axes[0] = axes[2];
730 axes[2] = tmpval;
731 }
732
733 return 1;
734
735 err:
736 return 0;
737 }
738
739 #ifdef SPGDEBUG
lauemmm(int axes[3],SPGCONST PointSymmetry * pointsym)740 static int lauemmm(int axes[3],
741 SPGCONST PointSymmetry * pointsym)
742 {
743 int i, count, axis;
744 int prop_rot[3][3];
745
746
747 for (i = 0; i < 3; i++) { axes[i] = -1; }
748 count = 0;
749 for (i = 0; i < pointsym->size; i++) {
750 get_proper_rotation(prop_rot, pointsym->rot[i]);
751
752 /* Search two-fold rotation */
753 if (! (mat_get_trace_i3(prop_rot) == -1)) {
754 continue;
755 }
756
757 axis = get_rotation_axis(prop_rot);
758 if (! ((axis == axes[0]) ||
759 (axis == axes[1]) ||
760 (axis == axes[2]))) {
761 axes[count] = axis;
762 count++;
763 }
764 }
765
766 sort_axes(axes);
767
768 return 1;
769 }
770
laue4m(int axes[3],SPGCONST PointSymmetry * pointsym)771 static int laue4m(int axes[3],
772 SPGCONST PointSymmetry * pointsym)
773 {
774 int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
775 int axis_vec[3];
776 int prop_rot[3][3], t_mat[3][3];
777 int ortho_axes[NUM_ROT_AXES];
778
779 for (i = 0; i < pointsym->size; i++) {
780 get_proper_rotation(prop_rot, pointsym->rot[i]);
781
782 /* Search foud-fold rotation */
783 if ( mat_get_trace_i3(prop_rot) == 1) {
784 /* The first axis */
785 axes[2] = get_rotation_axis(prop_rot);
786 break;
787 }
788 }
789
790 /* The second axis */
791 num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, 4);
792 if (! num_ortho_axis) { goto err; }
793
794 min_norm = 8;
795 is_found = 0;
796 for (i = 0; i < num_ortho_axis; i++) {
797 norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
798 if (norm < min_norm) {
799 min_norm = norm;
800 axes[0] = ortho_axes[i];
801 is_found = 1;
802 }
803 }
804 if (! is_found) { goto err; }
805
806 /* The third axis */
807 mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
808 is_found = 0;
809 for (i = 0; i < NUM_ROT_AXES; i++) {
810 if (is_exist_axis(axis_vec, i)) {
811 is_found = 1;
812 axes[1] = i;
813 break;
814 }
815 }
816 if (! is_found) { goto err; }
817
818 set_transformation_matrix(t_mat, axes);
819 if (mat_get_determinant_i3(t_mat) < 0) {
820 tmpval = axes[0];
821 axes[0] = axes[1];
822 axes[1] = tmpval;
823 }
824
825 return 1;
826
827 err:
828 return 0;
829 }
830
laue4mmm(int axes[3],SPGCONST PointSymmetry * pointsym)831 static int laue4mmm(int axes[3],
832 SPGCONST PointSymmetry * pointsym)
833 {
834 int i, is_found, tmpval, axis;
835 int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
836 int axis_vec[3];
837
838 for (i = 0; i < pointsym->size; i++) {
839 get_proper_rotation(prop_rot, pointsym->rot[i]);
840
841 /* Search foud-fold rotation */
842 if (mat_get_trace_i3(prop_rot) == 1) {
843 /* The first axis */
844 axes[2] = get_rotation_axis(prop_rot);
845 break;
846 }
847 }
848
849 is_found = 0;
850 for (i = 0; i < pointsym->size; i++) {
851 get_proper_rotation(prop_rot2, pointsym->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
869 /* The third axis */
870 mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
871 is_found = 0;
872 for (i = 0; i < NUM_ROT_AXES; i++) {
873 if (is_exist_axis(axis_vec, i)) {
874 is_found = 1;
875 axes[1] = i;
876 break;
877 }
878 }
879 if (! is_found) { goto err; }
880
881 set_transformation_matrix(t_mat, axes);
882 if (mat_get_determinant_i3(t_mat) < 0) {
883 tmpval = axes[0];
884 axes[0] = axes[1];
885 axes[1] = tmpval;
886 }
887
888 return 1;
889
890 err:
891 return 0;
892 }
893
laue3(int axes[3],SPGCONST PointSymmetry * pointsym)894 static int laue3(int axes[3],
895 SPGCONST PointSymmetry * pointsym)
896 {
897 int i, num_ortho_axis, norm, min_norm, is_found, tmpval;
898 int prop_rot[3][3], t_mat[3][3];
899 int axis_vec[3];
900 int ortho_axes[NUM_ROT_AXES];
901
902 for (i = 0; i < pointsym->size; i++) {
903 get_proper_rotation(prop_rot, pointsym->rot[i]);
904
905 /* Search thee-fold rotation */
906 if (mat_get_trace_i3(prop_rot) == 0) {
907 /* The first axis */
908 axes[2] = get_rotation_axis(prop_rot);
909 break;
910 }
911 }
912
913 /* The second axis */
914 num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, 3);
915 if (! num_ortho_axis) { goto err; }
916 min_norm = 8;
917 is_found = 0;
918 for (i = 0; i < num_ortho_axis; i++) {
919 norm = mat_norm_squared_i3(rot_axes[ortho_axes[i]]);
920 if (norm < min_norm) {
921 min_norm = norm;
922 axes[0] = ortho_axes[i];
923 is_found = 1;
924 }
925 }
926 if (! is_found) { goto err; }
927
928 /* The third axis */
929 mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
930 is_found = 0;
931 for (i = 0; i < NUM_ROT_AXES; i++) {
932 is_found = is_exist_axis(axis_vec, i);
933 if (is_found == 1) {
934 axes[1] = i;
935 break;
936 }
937 if (is_found == -1) {
938 axes[1] = i + NUM_ROT_AXES;
939 break;
940 }
941 }
942 if (! is_found) { goto err; }
943
944 set_transformation_matrix(t_mat, axes);
945 if (mat_get_determinant_i3(t_mat) < 0) {
946 tmpval = axes[0];
947 axes[0] = axes[1];
948 axes[1] = tmpval;
949 }
950
951 return 1;
952
953 err:
954 return 0;
955 }
956
laue3m(int axes[3],SPGCONST PointSymmetry * pointsym)957 static int laue3m(int axes[3],
958 SPGCONST PointSymmetry * pointsym)
959 {
960 int i, is_found, tmpval, axis;
961 int prop_rot[3][3], prop_rot2[3][3], t_mat[3][3];
962 int axis_vec[3];
963
964 for (i = 0; i < pointsym->size; i++) {
965 get_proper_rotation(prop_rot, pointsym->rot[i]);
966
967 /* Search three-fold rotation */
968 if (mat_get_trace_i3(prop_rot) == 0) {
969 /* The first axis */
970 axes[2] = get_rotation_axis(prop_rot);
971 debug_print("laue3m prop_rot\n");
972 debug_print_matrix_i3(prop_rot);
973 break;
974 }
975 }
976
977 is_found = 0;
978 for (i = 0; i < pointsym->size; i++) {
979 get_proper_rotation(prop_rot2, pointsym->rot[i]);
980
981 /* Search two-fold rotation */
982 if (! (mat_get_trace_i3(prop_rot2) == -1)) {
983 continue;
984 }
985
986 /* The second axis */
987 axis = get_rotation_axis(prop_rot2);
988 if (! (axis == axes[2])) {
989 axes[0] = axis;
990 is_found = 1;
991 break;
992 }
993 }
994 if (! is_found) { goto err; }
995
996 /* The third axis */
997 mat_multiply_matrix_vector_i3(axis_vec, prop_rot, rot_axes[axes[0]]);
998 is_found = 0;
999 for (i = 0; i < NUM_ROT_AXES; i++) {
1000 is_found = is_exist_axis(axis_vec, i);
1001 if (is_found == 1) {
1002 axes[1] = i;
1003 break;
1004 }
1005 if (is_found == -1) {
1006 axes[1] = i + NUM_ROT_AXES;
1007 break;
1008 }
1009 }
1010 if (! is_found) { goto err; }
1011
1012 set_transformation_matrix(t_mat, axes);
1013 if (mat_get_determinant_i3(t_mat) < 0) {
1014 tmpval = axes[0];
1015 axes[0] = axes[1];
1016 axes[1] = tmpval;
1017 }
1018
1019 return 1;
1020
1021 err:
1022 return 0;
1023 }
1024
lauem3m(int axes[3],SPGCONST PointSymmetry * pointsym)1025 static int lauem3m(int axes[3],
1026 SPGCONST PointSymmetry * pointsym)
1027 {
1028 int i, count, axis;
1029 int prop_rot[3][3];
1030
1031 for (i = 0; i < 3; i++) { axes[i] = -1; }
1032 count = 0;
1033 for (i = 0; i < pointsym->size; i++) {
1034 get_proper_rotation(prop_rot, pointsym->rot[i]);
1035
1036 /* Search four-fold rotation */
1037 if (! (mat_get_trace_i3(prop_rot) == 1)) {
1038 continue;
1039 }
1040
1041 axis = get_rotation_axis(prop_rot);
1042 if (! ((axis == axes[0]) ||
1043 (axis == axes[1]) ||
1044 (axis == axes[2]))) {
1045 axes[count] = axis;
1046 count++;
1047 }
1048 }
1049
1050 sort_axes(axes);
1051
1052 return 1;
1053 }
1054 #endif
1055
laue_one_axis(int axes[3],SPGCONST PointSymmetry * pointsym,const int rot_order)1056 static int laue_one_axis(int axes[3],
1057 SPGCONST PointSymmetry * pointsym,
1058 const int rot_order)
1059 {
1060 int i, j, num_ortho_axis, det, is_found, tmpval;
1061 int axis_vec[3], tmp_axes[3];
1062 int prop_rot[3][3], t_mat[3][3];
1063 int ortho_axes[NUM_ROT_AXES];
1064
1065 debug_print("laue_one_axis with rot_order %d\n", rot_order);
1066
1067 for (i = 0; i < pointsym->size; i++) {
1068 get_proper_rotation(prop_rot, pointsym->rot[i]);
1069
1070 /* Search foud-fold rotation */
1071 if (rot_order == 4) {
1072 if (mat_get_trace_i3(prop_rot) == 1) {
1073 /* The first axis */
1074 axes[2] = get_rotation_axis(prop_rot);
1075 break;
1076 }
1077 }
1078
1079 /* Search three-fold rotation */
1080 if (rot_order == 3) {
1081 if (mat_get_trace_i3(prop_rot) == 0) {
1082 /* The first axis */
1083 axes[2] = get_rotation_axis(prop_rot);
1084 break;
1085 }
1086 }
1087 }
1088
1089 /* Candidates of the second axis */
1090 num_ortho_axis = get_orthogonal_axis(ortho_axes, prop_rot, rot_order);
1091 if (! num_ortho_axis) { goto err; }
1092
1093 tmp_axes[1] = -1;
1094 tmp_axes[2] = axes[2];
1095 for (i = 0; i < num_ortho_axis; i++) {
1096 is_found = 0;
1097 tmp_axes[0] = ortho_axes[i];
1098 mat_multiply_matrix_vector_i3(axis_vec,
1099 prop_rot,
1100 rot_axes[tmp_axes[0]]);
1101 for (j = 0; j < num_ortho_axis; j++) {
1102 is_found = is_exist_axis(axis_vec, ortho_axes[j]);
1103 if (is_found == 1) {
1104 tmp_axes[1] = ortho_axes[j];
1105 break;
1106 }
1107 if (is_found == -1) {
1108 tmp_axes[1] = ortho_axes[j] + NUM_ROT_AXES;
1109 break;
1110 }
1111 }
1112
1113 if (!is_found) { continue; }
1114
1115 set_transformation_matrix(t_mat, tmp_axes);
1116 det = abs(mat_get_determinant_i3(t_mat));
1117 if (det < 4) { /* to avoid F-center choice det=4 */
1118 axes[0] = tmp_axes[0];
1119 axes[1] = tmp_axes[1];
1120 goto end;
1121 }
1122 }
1123
1124 err: /* axes are not correctly found. */
1125 warning_print("spglib: Secondary axis is not found.");
1126 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
1127 return 0;
1128
1129 end:
1130 set_transformation_matrix(t_mat, axes);
1131 if (mat_get_determinant_i3(t_mat) < 0) {
1132 tmpval = axes[0];
1133 axes[0] = axes[1];
1134 axes[1] = tmpval;
1135 }
1136
1137 debug_print("axes[0] = %d\n", axes[0]);
1138 debug_print("axes[1] = %d\n", axes[1]);
1139 debug_print("axes[2] = %d\n", axes[2]);
1140
1141 return 1;
1142
1143 }
1144
lauennn(int axes[3],SPGCONST PointSymmetry * pointsym,const int rot_order)1145 static int lauennn(int axes[3],
1146 SPGCONST PointSymmetry * pointsym,
1147 const int rot_order)
1148 {
1149 int i, count, axis;
1150 int prop_rot[3][3];
1151
1152 for (i = 0; i < 3; i++) {
1153 axes[i] = -1;
1154 }
1155
1156 count = 0;
1157 for (i = 0; i < pointsym->size; i++) {
1158 get_proper_rotation(prop_rot, pointsym->rot[i]);
1159
1160 /* Search two- or four-fold rotation */
1161 if ((mat_get_trace_i3(prop_rot) == -1 && rot_order == 2) ||
1162 (mat_get_trace_i3(prop_rot) == 1 && rot_order == 4)) {
1163 axis = get_rotation_axis(prop_rot);
1164 if (! ((axis == axes[0]) ||
1165 (axis == axes[1]) ||
1166 (axis == axes[2]))) {
1167 axes[count] = axis;
1168 count++;
1169 }
1170 }
1171 }
1172
1173 sort_axes(axes);
1174
1175 return 1;
1176 }
1177
get_rotation_axis(SPGCONST int proper_rot[3][3])1178 static int get_rotation_axis(SPGCONST int proper_rot[3][3])
1179 {
1180 int i, axis = -1;
1181 int vec[3];
1182
1183 /* No specific axis for I and -I */
1184 if (mat_check_identity_matrix_i3(proper_rot, identity)) {
1185 goto end;
1186 }
1187
1188 /* Look for eigenvector = rotation axis */
1189 for (i = 0; i < NUM_ROT_AXES; i++) {
1190 mat_multiply_matrix_vector_i3(vec, proper_rot, rot_axes[i]);
1191 if (vec[0] == rot_axes[i][0] &&
1192 vec[1] == rot_axes[i][1] &&
1193 vec[2] == rot_axes[i][2]) {
1194 axis = i;
1195 break;
1196 }
1197 }
1198
1199 end:
1200 #ifdef SPGDEBUG
1201 if (axis == -1) {
1202 printf("rotation axis cound not found.\n");
1203 }
1204 #endif
1205
1206 return axis;
1207 }
1208
get_orthogonal_axis(int ortho_axes[],SPGCONST int proper_rot[3][3],const int rot_order)1209 static int get_orthogonal_axis(int ortho_axes[],
1210 SPGCONST int proper_rot[3][3],
1211 const int rot_order)
1212 {
1213 int i, num_ortho_axis;
1214 int vec[3];
1215 int sum_rot[3][3], rot[3][3];
1216
1217 num_ortho_axis = 0;
1218
1219 mat_copy_matrix_i3(sum_rot, identity);
1220 mat_copy_matrix_i3(rot, identity);
1221 for (i = 0; i < rot_order - 1; i++) {
1222 mat_multiply_matrix_i3(rot, proper_rot, rot);
1223 mat_add_matrix_i3(sum_rot, rot, sum_rot);
1224 }
1225
1226 for (i = 0; i < NUM_ROT_AXES; i++) {
1227 mat_multiply_matrix_vector_i3(vec, sum_rot, rot_axes[i]);
1228 if (vec[0] == 0 &&
1229 vec[1] == 0 &&
1230 vec[2] == 0) {
1231 ortho_axes[num_ortho_axis] = i;
1232 num_ortho_axis++;
1233 }
1234 }
1235
1236 return num_ortho_axis;
1237 }
1238
get_proper_rotation(int prop_rot[3][3],SPGCONST int rot[3][3])1239 static void get_proper_rotation(int prop_rot[3][3],
1240 SPGCONST int rot[3][3])
1241 {
1242 if (mat_get_determinant_i3(rot) == -1) {
1243 mat_multiply_matrix_i3(prop_rot, inversion, rot);
1244 } else {
1245 mat_copy_matrix_i3(prop_rot, rot);
1246 }
1247 }
1248
set_transformation_matrix(int tmat[3][3],const int axes[3])1249 static void set_transformation_matrix(int tmat[3][3],
1250 const int axes[3])
1251 {
1252 int i, j, s[3];
1253
1254 for (i = 0; i < 3; i++) {
1255 if (axes[i] < NUM_ROT_AXES) {
1256 s[i] = 1;
1257 } else {
1258 s[i] = -1; /* axes[i] + NUM_ROT_AXES means improper rotation. */
1259 }
1260 }
1261 for (i = 0; i < 3; i++) {
1262 for (j = 0; j < 3; j++) {
1263 tmat[i][j] = s[j] * rot_axes[axes[j]%NUM_ROT_AXES][i];
1264 }
1265 }
1266 }
1267
is_exist_axis(const int axis_vec[3],const int axis_index)1268 static int is_exist_axis(const int axis_vec[3], const int axis_index)
1269 {
1270 if ((axis_vec[0] == rot_axes[axis_index][0]) &&
1271 (axis_vec[1] == rot_axes[axis_index][1]) &&
1272 (axis_vec[2] == rot_axes[axis_index][2])) { return 1; }
1273 if ((axis_vec[0] == -rot_axes[axis_index][0]) &&
1274 (axis_vec[1] == -rot_axes[axis_index][1]) &&
1275 (axis_vec[2] == -rot_axes[axis_index][2])) { return -1; }
1276 return 0;
1277 }
1278
sort_axes(int axes[3])1279 static void sort_axes(int axes[3])
1280 {
1281 int axis;
1282 int t_mat[3][3];
1283
1284 if (axes[1] > axes[2]) {
1285 axis = axes[1];
1286 axes[1] = axes[2];
1287 axes[2] = axis;
1288 }
1289
1290 if (axes[0] > axes[1]) {
1291 axis = axes[0];
1292 axes[0] = axes[1];
1293 axes[1] = axis;
1294 }
1295
1296 if (axes[1] > axes[2]) {
1297 axis = axes[1];
1298 axes[1] = axes[2];
1299 axes[2] = axis;
1300 }
1301
1302 set_transformation_matrix(t_mat, axes);
1303 if (mat_get_determinant_i3(t_mat) < 0) {
1304 axis = axes[1];
1305 axes[1] = axes[2];
1306 axes[2] = axis;
1307 }
1308 }
1309
1310