1 /* Copyright (C) 2010 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 <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include "cell.h"
39 #include "delaunay.h"
40 #include "hall_symbol.h"
41 #include "mathfunc.h"
42 #include "niggli.h"
43 #include "pointgroup.h"
44 #include "primitive.h"
45 #include "spacegroup.h"
46 #include "spg_database.h"
47 #include "symmetry.h"
48
49 #include "debug.h"
50
51 #define REDUCE_RATE 0.95
52 #define NUM_ATTEMPT 20
53 #define INT_PREC 0.1
54
55 static double change_of_basis_monocli[36][3][3] = {
56 {{ 1, 0, 0 }, /* b first turn; two axes are flipped in second turn */
57 { 0, 1, 0 },
58 { 0, 0, 1 }},
59 {{ 0, 0, 1 }, /* b */
60 { 0,-1, 0 },
61 { 1, 0, 0 }},
62 {{ 0, 0, 1 }, /* a */
63 { 1, 0, 0 },
64 { 0, 1, 0 }},
65 {{ 1, 0, 0 }, /* c */
66 { 0, 0, 1 },
67 { 0,-1, 0 }},
68 {{ 0, 1, 0 }, /* c */
69 { 0, 0, 1 },
70 { 1, 0, 0 }},
71 {{ 0,-1, 0 }, /* a */
72 { 1, 0, 0 },
73 { 0, 0, 1 }},
74 {{-1, 0, 1 }, /* b */
75 { 0, 1, 0 },
76 {-1, 0, 0 }},
77 {{ 1, 0,-1 }, /* b */
78 { 0,-1, 0 },
79 { 0, 0,-1 }},
80 {{ 0, 1,-1 }, /* a */
81 { 1, 0, 0 },
82 { 0, 0,-1 }},
83 {{-1,-1, 0 }, /* c */
84 { 0, 0, 1 },
85 {-1, 0, 0 }},
86 {{ 1,-1, 0 }, /* c */
87 { 0, 0, 1 },
88 { 0,-1, 0 }},
89 {{ 0, 1, 1 }, /* a */
90 { 1, 0, 0 },
91 { 0, 1, 0 }},
92 {{ 0, 0,-1 }, /* b */
93 { 0, 1, 0 },
94 { 1, 0,-1 }},
95 {{-1, 0, 0 }, /* b */
96 { 0,-1, 0 },
97 {-1, 0, 1 }},
98 {{ 0,-1, 0 }, /* a */
99 { 1, 0, 0 },
100 { 0,-1, 1 }},
101 {{ 0, 1, 0 }, /* c */
102 { 0, 0, 1 },
103 { 1, 1, 0 }},
104 {{-1, 0, 0 }, /* c */
105 { 0, 0, 1 },
106 {-1, 1, 0 }},
107 {{ 0, 0,-1 }, /* a */
108 { 1, 0, 0 },
109 { 0,-1,-1 }},
110 {{ 1, 0, 0 }, /* b two axes are flipped to look for non-acute axes */
111 { 0,-1, 0 },
112 { 0, 0,-1 }},
113 {{ 0, 0,-1 }, /* b */
114 { 0, 1, 0 },
115 { 1, 0, 0 }},
116 {{ 0, 0, 1 }, /* a */
117 {-1, 0, 0 },
118 { 0,-1, 0 }},
119 {{-1, 0, 0 }, /* c */
120 { 0, 0,-1 },
121 { 0,-1, 0 }},
122 {{ 0, 1, 0 }, /* c */
123 { 0, 0,-1 },
124 {-1, 0, 0 }},
125 {{ 0, 1, 0 }, /* a */
126 {-1, 0, 0 },
127 { 0, 0, 1 }},
128 {{-1, 0,-1 }, /* b */
129 { 0,-1, 0 },
130 {-1, 0, 0 }},
131 {{ 1, 0, 1 }, /* b */
132 { 0, 1, 0 },
133 { 0, 0, 1 }},
134 {{ 0,-1,-1 }, /* a */
135 {-1, 0, 0 },
136 { 0, 0,-1 }},
137 {{ 1,-1, 0 }, /* c */
138 { 0, 0,-1 },
139 { 1, 0, 0 }},
140 {{-1,-1, 0 }, /* c */
141 { 0, 0,-1 },
142 { 0,-1, 0 }},
143 {{ 0,-1, 1 }, /* a */
144 {-1, 0, 0 },
145 { 0,-1, 0 }},
146 {{ 0, 0, 1 }, /* b */
147 { 0,-1, 0 },
148 { 1, 0, 1 }},
149 {{-1, 0, 0 }, /* b */
150 { 0, 1, 0 },
151 {-1, 0,-1 }},
152 {{ 0, 1, 0 }, /* a */
153 {-1, 0, 0 },
154 { 0, 1, 1 }},
155 {{ 0, 1, 0 }, /* c */
156 { 0, 0,-1 },
157 {-1, 1, 0 }},
158 {{ 1, 0, 0 }, /* c */
159 { 0, 0,-1 },
160 { 1, 1, 0 }},
161 {{ 0, 0,-1 }, /* a */
162 {-1, 0, 0 },
163 { 0, 1,-1 }}};
164
165 static Centering change_of_centering_monocli[36] = {
166 C_FACE, /* first turn */
167 A_FACE,
168 B_FACE,
169 B_FACE,
170 A_FACE,
171 C_FACE,
172 BASE,
173 BASE,
174 BASE,
175 BASE,
176 BASE,
177 BASE,
178 A_FACE,
179 C_FACE,
180 C_FACE,
181 A_FACE,
182 B_FACE,
183 B_FACE,
184 C_FACE, /* second turn */
185 A_FACE,
186 B_FACE,
187 B_FACE,
188 A_FACE,
189 C_FACE,
190 BASE,
191 BASE,
192 BASE,
193 BASE,
194 BASE,
195 BASE,
196 A_FACE,
197 C_FACE,
198 C_FACE,
199 A_FACE,
200 B_FACE,
201 B_FACE};
202
203 static int change_of_unique_axis_monocli[36] = {
204 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0,
205 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0};
206
207 static double change_of_basis_ortho[6][3][3] = {{{ 1, 0, 0 },
208 { 0, 1, 0 },
209 { 0, 0, 1 }},
210 {{ 0, 0, 1 },
211 { 1, 0, 0 },
212 { 0, 1, 0 }},
213 {{ 0, 1, 0 },
214 { 0, 0, 1 },
215 { 1, 0, 0 }},
216 {{ 0, 1, 0 },
217 { 1, 0, 0 },
218 { 0, 0,-1 }},
219 {{ 1, 0, 0 },
220 { 0, 0, 1 },
221 { 0,-1, 0 }},
222 {{ 0, 0, 1 },
223 { 0, 1, 0 },
224 {-1, 0, 0 }}};
225
226 static Centering change_of_centering_ortho[6] = {C_FACE,
227 B_FACE,
228 A_FACE,
229 C_FACE,
230 B_FACE,
231 A_FACE};
232 static int change_of_unique_axis_ortho[6] = {2, 1, 0, 2, 1, 0};
233
234 static double hR_to_hP[3][3] = {{ 1, 0, 1 },
235 {-1, 1, 1 },
236 { 0,-1, 1 }};
237 static double change_of_basis_501[3][3] = {{ 0, 0, 1},
238 { 0,-1, 0},
239 { 1, 0, 0}};
240
241 static int spacegroup_to_hall_number[230] = {
242 1, 2, 3, 6, 9, 18, 21, 30, 39, 57,
243 60, 63, 72, 81, 90, 108, 109, 112, 115, 116,
244 119, 122, 123, 124, 125, 128, 134, 137, 143, 149,
245 155, 161, 164, 170, 173, 176, 182, 185, 191, 197,
246 203, 209, 212, 215, 218, 221, 227, 228, 230, 233,
247 239, 245, 251, 257, 263, 266, 269, 275, 278, 284,
248 290, 292, 298, 304, 310, 313, 316, 322, 334, 335,
249 337, 338, 341, 343, 349, 350, 351, 352, 353, 354,
250 355, 356, 357, 358, 359, 361, 363, 364, 366, 367,
251 368, 369, 370, 371, 372, 373, 374, 375, 376, 377,
252 378, 379, 380, 381, 382, 383, 384, 385, 386, 387,
253 388, 389, 390, 391, 392, 393, 394, 395, 396, 397,
254 398, 399, 400, 401, 402, 404, 406, 407, 408, 410,
255 412, 413, 414, 416, 418, 419, 420, 422, 424, 425,
256 426, 428, 430, 431, 432, 433, 435, 436, 438, 439,
257 440, 441, 442, 443, 444, 446, 447, 448, 449, 450,
258 452, 454, 455, 456, 457, 458, 460, 462, 463, 464,
259 465, 466, 467, 468, 469, 470, 471, 472, 473, 474,
260 475, 476, 477, 478, 479, 480, 481, 482, 483, 484,
261 485, 486, 487, 488, 489, 490, 491, 492, 493, 494,
262 495, 497, 498, 500, 501, 502, 503, 504, 505, 506,
263 507, 508, 509, 510, 511, 512, 513, 514, 515, 516,
264 517, 518, 520, 521, 523, 524, 525, 527, 529, 530,
265 };
266
267 static double identity[3][3] = {{ 1, 0, 0 },
268 { 0, 1, 0 },
269 { 0, 0, 1 }};
270 static double monocli_i2c[3][3] = {{ 1, 0,-1 },
271 { 0, 1, 0 },
272 { 1, 0, 0 }};
273 static double monocli_a2c[3][3] = {{ 0, 0, 1 },
274 { 0,-1, 0 },
275 { 1, 0, 0 }};
276 static double rhombo_obverse[3][3] = {{ 2./3,-1./3,-1./3 },
277 { 1./3, 1./3,-2./3 },
278 { 1./3, 1./3, 1./3 }};
279 static double rhomb_reverse[3][3] = {{ 1./3,-2./3, 1./3 },
280 { 2./3,-1./3,-1./3 },
281 { 1./3, 1./3, 1./3 }};
282 static double a2c[3][3] = {{ 0, 0, 1 },
283 { 1, 0, 0 },
284 { 0, 1, 0 }};
285 static double b2c[3][3] = {{ 0, 1, 0 },
286 { 0, 0, 1 },
287 { 1, 0, 0 }};
288
289 static double A_mat[3][3] = {{ 1, 0, 0},
290 { 0, 1./2,-1./2},
291 { 0, 1./2, 1./2}};
292 static double C_mat[3][3] = {{ 1./2, 1./2, 0},
293 {-1./2, 1./2, 0},
294 { 0, 0, 1}};
295 static double R_mat[3][3] = {{ 2./3,-1./3,-1./3 },
296 { 1./3, 1./3,-2./3 },
297 { 1./3, 1./3, 1./3 }};
298 static double I_mat[3][3] = {{-1./2, 1./2, 1./2 },
299 { 1./2,-1./2, 1./2 },
300 { 1./2, 1./2,-1./2 }};
301 static double F_mat[3][3] = {{ 0, 1./2, 1./2 },
302 { 1./2, 0, 1./2 },
303 { 1./2, 1./2, 0 }};
304
305 static Spacegroup * search_spacegroup_with_symmetry(const Cell * primitive,
306 const int candidates[],
307 const int num_candidates,
308 const Symmetry *symmetry,
309 const double symprec,
310 const double angle_tolerance);
311 static Spacegroup * get_spacegroup(const int hall_number,
312 const double origin_shift[3],
313 SPGCONST double conv_lattice[3][3]);
314 static int iterative_search_hall_number(double origin_shift[3],
315 double conv_lattice[3][3],
316 const int candidates[],
317 const int num_candidates,
318 const Cell * primitive,
319 const Symmetry * symmetry,
320 const double symprec,
321 const double angle_tolerance);
322 static int change_basis_tricli(int int_transform_mat[3][3],
323 SPGCONST double conv_lattice[3][3],
324 SPGCONST double primitive_lattice[3][3],
325 const double symprec);
326 static int change_basis_monocli(int int_transform_mat[3][3],
327 SPGCONST double conv_lattice[3][3],
328 SPGCONST double primitive_lattice[3][3],
329 const double symprec);
330 static Symmetry *
331 get_initial_conventional_symmetry(const Centering centering,
332 SPGCONST double transform_mat[3][3],
333 const Symmetry * symmetry);
334 static int search_hall_number(double origin_shift[3],
335 double conv_lattice[3][3],
336 const int candidates[],
337 const int num_candidates,
338 SPGCONST double primitive_lattice[3][3],
339 const Symmetry * symmetry,
340 const double symprec);
341 static int match_hall_symbol_db(double origin_shift[3],
342 double lattice[3][3],
343 const int hall_number,
344 const int pointgroup_number,
345 const Holohedry holohedry,
346 const Centering centering,
347 const Symmetry *symmetry,
348 const double symprec);
349 static int match_hall_symbol_db_monocli(double origin_shift[3],
350 double lattice[3][3],
351 const int hall_number,
352 const int num_hall_types,
353 const Centering centering,
354 const Symmetry *symmetry,
355 const double symprec);
356 static int match_hall_symbol_db_ortho(double origin_shift[3],
357 double lattice[3][3],
358 const int hall_number,
359 const Centering centering,
360 const Symmetry *symmetry,
361 const int num_free_axes,
362 const double symprec);
363 static Symmetry * get_conventional_symmetry(SPGCONST double transform_mat[3][3],
364 const Centering centering,
365 const Symmetry *primitive_sym);
366 static Centering get_centering(double correction_mat[3][3],
367 SPGCONST int transform_mat[3][3],
368 const Laue laue);
369 static Centering get_base_center(SPGCONST int transform_mat[3][3]);
370 static int get_centering_shifts(double shift[3][3],
371 const Centering centering);
372
373
374 /* Return spacegroup.number = 0 if failed */
spa_search_spacegroup(const Cell * primitive,const int hall_number,const double symprec,const double angle_tolerance)375 Spacegroup * spa_search_spacegroup(const Cell * primitive,
376 const int hall_number,
377 const double symprec,
378 const double angle_tolerance)
379 {
380 Spacegroup *spacegroup;
381 Symmetry *symmetry;
382 int candidate[1];
383
384 debug_print("search_spacegroup (tolerance = %f):\n", symprec);
385
386 symmetry = NULL;
387 spacegroup = NULL;
388
389 if ((symmetry = sym_get_operation(primitive, symprec, angle_tolerance)) ==
390 NULL) {
391 return NULL;
392 }
393
394 if (hall_number > 0) {
395 candidate[0] = hall_number;
396 }
397
398 if (hall_number) {
399 spacegroup = search_spacegroup_with_symmetry(primitive,
400 candidate,
401 1,
402 symmetry,
403 symprec,
404 angle_tolerance);
405 } else {
406 spacegroup = search_spacegroup_with_symmetry(primitive,
407 spacegroup_to_hall_number,
408 230,
409 symmetry,
410 symprec,
411 angle_tolerance);
412 }
413
414 sym_free_symmetry(symmetry);
415 symmetry = NULL;
416
417 return spacegroup;
418 }
419
420
spa_search_spacegroup_with_symmetry(const Symmetry * symmetry,const double symprec)421 int spa_search_spacegroup_with_symmetry(const Symmetry *symmetry,
422 const double symprec)
423 {
424 int i, hall_number;
425 Spacegroup *spacegroup;
426 Cell *primitive;
427
428 spacegroup = NULL;
429
430 primitive = cel_alloc_cell(1);
431 mat_copy_matrix_d3(primitive->lattice, identity);
432 for (i = 0; i < 3; i++) {
433 primitive->position[0][i] = 0;
434 }
435 spacegroup = search_spacegroup_with_symmetry(primitive,
436 spacegroup_to_hall_number,
437 230,
438 symmetry,
439 symprec,
440 -1.0);
441 hall_number = spacegroup->hall_number;
442 free(spacegroup);
443 spacegroup = NULL;
444
445 return hall_number;
446 }
447
448 /* Return NULL if failed */
spa_transform_to_primitive(int * mapping_table,const Cell * cell,SPGCONST double trans_mat[3][3],const Centering centering,const double symprec)449 Cell * spa_transform_to_primitive(int * mapping_table,
450 const Cell * cell,
451 SPGCONST double trans_mat[3][3],
452 const Centering centering,
453 const double symprec)
454 {
455 double tmat[3][3], tmat_inv[3][3], prim_lat[3][3];
456 Cell * primitive;
457
458 primitive = NULL;
459
460 if (!mat_inverse_matrix_d3(tmat_inv, trans_mat, symprec)) {
461 goto err;
462 }
463
464 switch (centering) {
465 case PRIMITIVE:
466 mat_copy_matrix_d3(tmat, tmat_inv);
467 break;
468 case A_FACE:
469 mat_multiply_matrix_d3(tmat, tmat_inv, A_mat);
470 break;
471 case C_FACE:
472 mat_multiply_matrix_d3(tmat, tmat_inv, C_mat);
473 break;
474 case FACE:
475 mat_multiply_matrix_d3(tmat, tmat_inv, F_mat);
476 break;
477 case BODY:
478 mat_multiply_matrix_d3(tmat, tmat_inv, I_mat);
479 break;
480 case R_CENTER:
481 mat_multiply_matrix_d3(tmat, tmat_inv, R_mat);
482 break;
483 default:
484 goto err;
485 }
486
487 mat_multiply_matrix_d3(prim_lat, cell->lattice, tmat);
488 if ((primitive = cel_trim_cell(mapping_table, prim_lat, cell, symprec))
489 == NULL) {
490 warning_print("spglib: cel_trim_cell failed.");
491 warning_print(" (line %d, %s).\n", __LINE__, __FILE__);
492 }
493
494 return primitive;
495
496 err:
497 return NULL;
498 }
499
500 /* Return NULL if failed */
spa_transform_from_primitive(const Cell * primitive,const Centering centering,const double symprec)501 Cell * spa_transform_from_primitive(const Cell * primitive,
502 const Centering centering,
503 const double symprec)
504 {
505 int multi, i, j, k, num_atom;
506 int *mapping_table;
507 double tmat[3][3], inv_tmat[3][3], shift[3][3];
508 Cell *std_cell, *trimmed_cell;
509
510 mapping_table = NULL;
511 trimmed_cell = NULL;
512 std_cell = NULL;
513
514 switch (centering) {
515 case PRIMITIVE:
516 break;
517 case A_FACE:
518 mat_copy_matrix_d3(tmat, A_mat);
519 mat_inverse_matrix_d3(inv_tmat, A_mat, 0);
520 break;
521 case C_FACE:
522 mat_copy_matrix_d3(tmat, C_mat);
523 mat_inverse_matrix_d3(inv_tmat, C_mat, 0);
524 break;
525 case FACE:
526 mat_copy_matrix_d3(tmat, F_mat);
527 mat_inverse_matrix_d3(inv_tmat, F_mat, 0);
528 break;
529 case BODY:
530 mat_copy_matrix_d3(tmat, I_mat);
531 mat_inverse_matrix_d3(inv_tmat, I_mat, 0);
532 break;
533 case R_CENTER:
534 mat_copy_matrix_d3(tmat, R_mat);
535 mat_inverse_matrix_d3(inv_tmat, R_mat, 0);
536 break;
537 default:
538 goto ret;
539 }
540
541 multi = get_centering_shifts(shift, centering);
542
543 if ((mapping_table = (int*) malloc(sizeof(int) * primitive->size * multi))
544 == NULL) {
545 warning_print("spglib: Memory could not be allocated ");
546 goto ret;
547 }
548
549 if ((std_cell = cel_alloc_cell(primitive->size * multi)) == NULL) {
550 free(mapping_table);
551 mapping_table = NULL;
552 goto ret;
553 }
554
555 mat_multiply_matrix_d3(std_cell->lattice, primitive->lattice, inv_tmat);
556
557 num_atom = 0;
558 for (i = 0; i < primitive->size; i++) {
559 mat_multiply_matrix_vector_d3(std_cell->position[num_atom],
560 tmat,
561 primitive->position[i]);
562 std_cell->types[num_atom] = primitive->types[i];
563 num_atom++;
564 }
565
566 for (i = 0; i < multi - 1; i++) {
567 for (j = 0; j < primitive->size; j++) {
568 mat_copy_vector_d3(std_cell->position[num_atom],
569 std_cell->position[j]);
570 for (k = 0; k < 3; k++) {
571 std_cell->position[num_atom][k] += shift[i][k];
572 }
573 std_cell->types[num_atom] = std_cell->types[j];
574 num_atom++;
575 }
576 }
577
578 trimmed_cell = cel_trim_cell(mapping_table,
579 std_cell->lattice,
580 std_cell,
581 symprec);
582 cel_free_cell(std_cell);
583 std_cell = NULL;
584 free(mapping_table);
585 mapping_table = NULL;
586
587 ret:
588 return trimmed_cell;
589 }
590
591 /* Return spacegroup.number = 0 if failed */
search_spacegroup_with_symmetry(const Cell * primitive,const int candidates[],const int num_candidates,const Symmetry * symmetry,const double symprec,const double angle_tolerance)592 static Spacegroup * search_spacegroup_with_symmetry(const Cell * primitive,
593 const int candidates[],
594 const int num_candidates,
595 const Symmetry *symmetry,
596 const double symprec,
597 const double angle_tolerance)
598 {
599 int hall_number;
600 double conv_lattice[3][3];
601 double origin_shift[3];
602 Spacegroup *spacegroup;
603 PointSymmetry pointsym;
604
605 debug_print("search_spacegroup (tolerance = %f):\n", symprec);
606
607 spacegroup = NULL;
608
609 pointsym = ptg_get_pointsymmetry(symmetry->rot, symmetry->size);
610 if (pointsym.size < symmetry->size) {
611 warning_print("spglib: Point symmetry of primitive cell is broken. ");
612 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
613 return NULL;
614 }
615
616 hall_number = iterative_search_hall_number(origin_shift,
617 conv_lattice,
618 candidates,
619 num_candidates,
620 primitive,
621 symmetry,
622 symprec,
623 angle_tolerance);
624 spacegroup = get_spacegroup(hall_number, origin_shift, conv_lattice);
625
626 return spacegroup;
627 }
628
629 /* Return spacegroup.number = 0 if failed */
get_spacegroup(const int hall_number,const double origin_shift[3],SPGCONST double conv_lattice[3][3])630 static Spacegroup * get_spacegroup(const int hall_number,
631 const double origin_shift[3],
632 SPGCONST double conv_lattice[3][3])
633 {
634 Spacegroup *spacegroup;
635 SpacegroupType spacegroup_type;
636
637 spacegroup = NULL;
638
639 if ((spacegroup = (Spacegroup*) malloc(sizeof(Spacegroup))) == NULL) {
640 warning_print("spglib: Memory could not be allocated.");
641 return NULL;
642 }
643
644 spacegroup_type = spgdb_get_spacegroup_type(hall_number);
645
646 if (spacegroup_type.number > 0) {
647 mat_copy_matrix_d3(spacegroup->bravais_lattice, conv_lattice);
648 mat_copy_vector_d3(spacegroup->origin_shift, origin_shift);
649 spacegroup->number = spacegroup_type.number;
650 spacegroup->hall_number = hall_number;
651 spacegroup->pointgroup_number = spacegroup_type.pointgroup_number;
652 strcpy(spacegroup->schoenflies,
653 spacegroup_type.schoenflies);
654 strcpy(spacegroup->hall_symbol,
655 spacegroup_type.hall_symbol);
656 strcpy(spacegroup->international,
657 spacegroup_type.international);
658 strcpy(spacegroup->international_long,
659 spacegroup_type.international_full);
660 strcpy(spacegroup->international_short,
661 spacegroup_type.international_short);
662 strcpy(spacegroup->choice,
663 spacegroup_type.choice);
664 }
665
666 return spacegroup;
667 }
668
669 /* Return 0 if failed */
iterative_search_hall_number(double origin_shift[3],double conv_lattice[3][3],const int candidates[],const int num_candidates,const Cell * primitive,const Symmetry * symmetry,const double symprec,const double angle_tolerance)670 static int iterative_search_hall_number(double origin_shift[3],
671 double conv_lattice[3][3],
672 const int candidates[],
673 const int num_candidates,
674 const Cell * primitive,
675 const Symmetry * symmetry,
676 const double symprec,
677 const double angle_tolerance)
678 {
679 int attempt, hall_number;
680 double tolerance;
681 Symmetry * sym_reduced;
682
683 debug_print("iterative_search_hall_number:\n");
684
685 hall_number = 0;
686 sym_reduced = NULL;
687
688 hall_number = search_hall_number(origin_shift,
689 conv_lattice,
690 candidates,
691 num_candidates,
692 primitive->lattice,
693 symmetry,
694 symprec);
695
696 if (hall_number > 0) {
697 goto ret;
698 }
699
700 tolerance = symprec;
701 for (attempt = 0; attempt < NUM_ATTEMPT; attempt++) {
702
703 warning_print("spglib: Attempt %d tolerance = %f failed",
704 attempt, tolerance);
705 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
706
707 tolerance *= REDUCE_RATE;
708 sym_reduced = sym_reduce_operation(primitive,
709 symmetry,
710 tolerance,
711 angle_tolerance);
712 hall_number = search_hall_number(origin_shift,
713 conv_lattice,
714 candidates,
715 num_candidates,
716 primitive->lattice,
717 sym_reduced,
718 symprec);
719 sym_free_symmetry(sym_reduced);
720 sym_reduced = NULL;
721 if (hall_number > 0) {
722 break;
723 }
724 }
725
726 ret:
727 return hall_number;
728 }
729
730 /* Return 0 if failed */
search_hall_number(double origin_shift[3],double conv_lattice[3][3],const int candidates[],const int num_candidates,SPGCONST double primitive_lattice[3][3],const Symmetry * symmetry,const double symprec)731 static int search_hall_number(double origin_shift[3],
732 double conv_lattice[3][3],
733 const int candidates[],
734 const int num_candidates,
735 SPGCONST double primitive_lattice[3][3],
736 const Symmetry * symmetry,
737 const double symprec)
738 {
739 int i, hall_number;
740 Centering centering;
741 Pointgroup pointgroup;
742 Symmetry * conv_symmetry;
743 int int_transform_mat[3][3];
744 double correction_mat[3][3], transform_mat[3][3];
745
746 debug_print("search_hall_number:\n");
747
748 hall_number = 0;
749 conv_symmetry = NULL;
750
751 pointgroup = ptg_get_transformation_matrix(int_transform_mat,
752 symmetry->rot,
753 symmetry->size);
754 if (pointgroup.number == 0) {
755 goto err;
756 }
757
758 mat_multiply_matrix_di3(conv_lattice, primitive_lattice, int_transform_mat);
759
760 if (pointgroup.laue == LAUE1) {
761 if (! change_basis_tricli(int_transform_mat,
762 conv_lattice,
763 primitive_lattice,
764 symprec)) {
765 goto err;
766 }
767 }
768
769 if (pointgroup.laue == LAUE2M) {
770 if (! change_basis_monocli(int_transform_mat,
771 conv_lattice,
772 primitive_lattice,
773 symprec)) {
774 goto err;
775 }
776 }
777
778 if ((centering = get_centering(correction_mat,
779 int_transform_mat,
780 pointgroup.laue)) == CENTERING_ERROR) {
781 goto err;
782 }
783
784 mat_multiply_matrix_id3(transform_mat, int_transform_mat, correction_mat);
785 mat_multiply_matrix_d3(conv_lattice, primitive_lattice, transform_mat);
786
787 if ((conv_symmetry = get_initial_conventional_symmetry(centering,
788 transform_mat,
789 symmetry)) == NULL) {
790 goto err;
791 }
792
793 for (i = 0; i < num_candidates; i++) {
794 if (match_hall_symbol_db(origin_shift,
795 conv_lattice, /* <-- modified only matched */
796 candidates[i],
797 pointgroup.number,
798 pointgroup.holohedry,
799 centering,
800 conv_symmetry,
801 symprec)) {
802 hall_number = candidates[i];
803 break;
804 }
805 }
806
807 sym_free_symmetry(conv_symmetry);
808 conv_symmetry = NULL;
809
810 return hall_number;
811
812 err:
813 return 0;
814 }
815
816 /* Triclinic: Niggli cell reduction */
817 /* Return 0 if failed */
change_basis_tricli(int int_transform_mat[3][3],SPGCONST double conv_lattice[3][3],SPGCONST double primitive_lattice[3][3],const double symprec)818 static int change_basis_tricli(int int_transform_mat[3][3],
819 SPGCONST double conv_lattice[3][3],
820 SPGCONST double primitive_lattice[3][3],
821 const double symprec)
822 {
823 int i, j;
824 double niggli_cell[9];
825 double smallest_lattice[3][3], inv_lattice[3][3], transform_mat[3][3];
826
827 for (i = 0; i < 3; i++) {
828 for (j = 0; j < 3; j++) {
829 niggli_cell[i * 3 + j] = conv_lattice[i][j];
830 }
831 }
832
833 if (! niggli_reduce(niggli_cell, symprec * symprec)) {
834 return 0;
835 }
836
837 for (i = 0; i < 3; i++) {
838 for (j = 0; j < 3; j++) {
839 smallest_lattice[i][j] = niggli_cell[i * 3 + j];
840 }
841 }
842 if (mat_get_determinant_d3(smallest_lattice) < 0) {
843 for (i = 0; i < 3; i++) {
844 for (j = 0; j < 3; j++) {
845 smallest_lattice[i][j] = -smallest_lattice[i][j];
846 }
847 }
848 }
849 mat_inverse_matrix_d3(inv_lattice, primitive_lattice, 0);
850 mat_multiply_matrix_d3(transform_mat, inv_lattice, smallest_lattice);
851 mat_cast_matrix_3d_to_3i(int_transform_mat, transform_mat);
852
853 return 1;
854 }
855
856 /* Monoclinic: choose shortest a, c lattice vectors (|a| < |c|) */
857 /* Return 0 if failed */
change_basis_monocli(int int_transform_mat[3][3],SPGCONST double conv_lattice[3][3],SPGCONST double primitive_lattice[3][3],const double symprec)858 static int change_basis_monocli(int int_transform_mat[3][3],
859 SPGCONST double conv_lattice[3][3],
860 SPGCONST double primitive_lattice[3][3],
861 const double symprec)
862 {
863 double smallest_lattice[3][3], inv_lattice[3][3], transform_mat[3][3];
864
865 if (! del_delaunay_reduce_2D(smallest_lattice,
866 conv_lattice,
867 1, /* unique axis of b */
868 symprec)) {
869 return 0;
870 }
871
872 mat_inverse_matrix_d3(inv_lattice, primitive_lattice, 0);
873 mat_multiply_matrix_d3(transform_mat, inv_lattice, smallest_lattice);
874 mat_cast_matrix_3d_to_3i(int_transform_mat, transform_mat);
875 return 1;
876 }
877
878 /* Return NULL if failed */
879 static Symmetry *
get_initial_conventional_symmetry(const Centering centering,SPGCONST double transform_mat[3][3],const Symmetry * symmetry)880 get_initial_conventional_symmetry(const Centering centering,
881 SPGCONST double transform_mat[3][3],
882 const Symmetry * symmetry)
883 {
884 Symmetry * conv_symmetry;
885
886 debug_print("get_initial_conventional_symmetry\n");
887
888 conv_symmetry = NULL;
889
890 if (centering == R_CENTER) {
891 /* hP for rhombohedral */
892 conv_symmetry = get_conventional_symmetry(transform_mat,
893 PRIMITIVE,
894 symmetry);
895 } else {
896 conv_symmetry = get_conventional_symmetry(transform_mat,
897 centering,
898 symmetry);
899 }
900
901 return conv_symmetry;
902 }
903
904 /* Return 0 if failed */
match_hall_symbol_db(double origin_shift[3],double lattice[3][3],const int hall_number,const int pointgroup_number,const Holohedry holohedry,const Centering centering,const Symmetry * symmetry,const double symprec)905 static int match_hall_symbol_db(double origin_shift[3],
906 double lattice[3][3],
907 const int hall_number,
908 const int pointgroup_number,
909 const Holohedry holohedry,
910 const Centering centering,
911 const Symmetry *symmetry,
912 const double symprec)
913 {
914 int is_found, num_hall_types;
915 SpacegroupType spacegroup_type;
916 Symmetry * changed_symmetry;
917 double changed_lattice[3][3], inv_lattice[3][3], transform_mat[3][3];
918
919 changed_symmetry = NULL;
920
921 spacegroup_type = spgdb_get_spacegroup_type(hall_number);
922 num_hall_types = (spacegroup_to_hall_number[spacegroup_type.number] -
923 spacegroup_to_hall_number[spacegroup_type.number - 1]);
924
925 if (pointgroup_number != spacegroup_type.pointgroup_number) {
926 goto err;
927 }
928
929 switch (holohedry) {
930 case MONOCLI:
931 if (match_hall_symbol_db_monocli(origin_shift,
932 lattice,
933 hall_number,
934 num_hall_types,
935 centering,
936 symmetry,
937 symprec)) {return 1;}
938 break;
939
940 case ORTHO:
941 if (spacegroup_type.number == 48 ||
942 spacegroup_type.number == 50 ||
943 spacegroup_type.number == 59 ||
944 spacegroup_type.number == 68 ||
945 spacegroup_type.number == 70) { /* uncount origin shift */
946 num_hall_types /= 2;
947 }
948
949 if (num_hall_types == 1) {
950 if (match_hall_symbol_db_ortho(origin_shift,
951 lattice,
952 hall_number,
953 centering,
954 symmetry,
955 6,
956 symprec)) {return 1;}
957 break;
958 }
959
960 if (num_hall_types == 2) {
961 if (match_hall_symbol_db_ortho(origin_shift,
962 lattice,
963 hall_number,
964 centering,
965 symmetry,
966 3,
967 symprec)) {return 1;}
968 break;
969 }
970
971 if (num_hall_types == 3) {
972 mat_copy_matrix_d3(changed_lattice, lattice);
973 if (! match_hall_symbol_db_ortho
974 (origin_shift,
975 changed_lattice,
976 spacegroup_to_hall_number[spacegroup_type.number - 1],
977 centering,
978 symmetry,
979 0,
980 symprec)) {break;}
981 mat_inverse_matrix_d3(inv_lattice, lattice, 0);
982 mat_multiply_matrix_d3(transform_mat, inv_lattice, changed_lattice);
983
984 if ((changed_symmetry = get_conventional_symmetry(transform_mat,
985 PRIMITIVE,
986 symmetry)) == NULL) {
987 goto err;
988 }
989
990 is_found = match_hall_symbol_db_ortho(origin_shift,
991 changed_lattice,
992 hall_number,
993 centering,
994 changed_symmetry,
995 2,
996 symprec);
997 sym_free_symmetry(changed_symmetry);
998 changed_symmetry = NULL;
999 if (is_found) {
1000 mat_copy_matrix_d3(lattice, changed_lattice);
1001 return 1;
1002 }
1003 break;
1004 }
1005
1006 if (num_hall_types == 6) {
1007 if (match_hall_symbol_db_ortho(origin_shift,
1008 lattice,
1009 hall_number,
1010 centering,
1011 symmetry,
1012 1,
1013 symprec)) {return 1;}
1014 break;
1015 }
1016
1017 break;
1018
1019 case CUBIC:
1020 if (hal_match_hall_symbol_db(origin_shift,
1021 lattice,
1022 hall_number,
1023 centering,
1024 symmetry,
1025 symprec)) {return 1;}
1026
1027 if (hall_number == 501) { /* Try another basis for No.205 */
1028 mat_multiply_matrix_d3(changed_lattice,
1029 lattice,
1030 change_of_basis_501);
1031 if ((changed_symmetry = get_conventional_symmetry(change_of_basis_501,
1032 PRIMITIVE,
1033 symmetry)) == NULL) {
1034 goto err;
1035 }
1036
1037 is_found = hal_match_hall_symbol_db(origin_shift,
1038 changed_lattice,
1039 hall_number,
1040 PRIMITIVE,
1041 changed_symmetry,
1042 symprec);
1043 sym_free_symmetry(changed_symmetry);
1044 changed_symmetry = NULL;
1045 if (is_found) {
1046 mat_copy_matrix_d3(lattice, changed_lattice);
1047 return 1;
1048 }
1049 }
1050 break;
1051
1052 case TRIGO:
1053 if (centering == R_CENTER) {
1054 if (hall_number == 433 ||
1055 hall_number == 436 ||
1056 hall_number == 444 ||
1057 hall_number == 450 ||
1058 hall_number == 452 ||
1059 hall_number == 458 ||
1060 hall_number == 460) {
1061 mat_multiply_matrix_d3(changed_lattice, lattice, hR_to_hP);
1062 if ((changed_symmetry =
1063 get_conventional_symmetry(hR_to_hP, R_CENTER, symmetry)) == NULL) {
1064 goto err;
1065 }
1066
1067 is_found = hal_match_hall_symbol_db(origin_shift,
1068 changed_lattice,
1069 hall_number,
1070 centering,
1071 changed_symmetry,
1072 symprec);
1073 sym_free_symmetry(changed_symmetry);
1074 changed_symmetry = NULL;
1075 if (is_found) {
1076 mat_copy_matrix_d3(lattice, changed_lattice);
1077 return 1;
1078 }
1079 } else {
1080 if (hal_match_hall_symbol_db(origin_shift,
1081 lattice,
1082 hall_number,
1083 PRIMITIVE,
1084 symmetry,
1085 symprec)) {
1086 return 1;
1087 }
1088 }
1089 break;
1090 }
1091 /* Do not break for other trigonal cases */
1092 default: /* HEXA, TETRA, TRICLI and rest of TRIGO */
1093 if (hal_match_hall_symbol_db(origin_shift,
1094 lattice,
1095 hall_number,
1096 centering,
1097 symmetry,
1098 symprec)) {
1099 return 1;
1100 }
1101 break;
1102 }
1103
1104 err:
1105 return 0;
1106 }
1107
1108 /* Return 0 if failed */
match_hall_symbol_db_monocli(double origin_shift[3],double lattice[3][3],const int hall_number,const int num_hall_types,const Centering centering,const Symmetry * symmetry,const double symprec)1109 static int match_hall_symbol_db_monocli(double origin_shift[3],
1110 double lattice[3][3],
1111 const int hall_number,
1112 const int num_hall_types,
1113 const Centering centering,
1114 const Symmetry *symmetry,
1115 const double symprec)
1116 {
1117 int i, j, k, l, is_found;
1118 double vec[2][3], norms[3];
1119 Centering changed_centering;
1120 Symmetry * changed_symmetry;
1121 double changed_lattice[3][3];
1122
1123 changed_symmetry = NULL;
1124
1125 for (i = 0; i < 36; i++) {
1126 /* centring type should be P or C */
1127 if (centering == C_FACE) {
1128 changed_centering = change_of_centering_monocli[i];
1129 } else { /* suppose PRIMITIVE */
1130 changed_centering = centering;
1131 }
1132
1133 mat_multiply_matrix_d3(changed_lattice,
1134 lattice,
1135 change_of_basis_monocli[i]);
1136
1137 /* Make non-acute and length preference */
1138 l = 0;
1139 for (j = 0; j < 3; j++) {
1140 if (j == change_of_unique_axis_monocli[i]) {continue;}
1141 for (k = 0; k < 3; k++) {vec[l][k] = changed_lattice[k][j];}
1142 norms[l] = mat_norm_squared_d3(vec[l]);
1143 l++;
1144 }
1145
1146 /* discard if principal angle is acute. */
1147 if (vec[0][0] * vec[1][0] +
1148 vec[0][1] * vec[1][1] +
1149 vec[0][2] * vec[1][2] > 0) {
1150 continue;
1151 }
1152
1153 /* Choose |a| < |b| < |c| if there are freedom. */
1154 if (num_hall_types == 3) {
1155 if (norms[0] > norms[1]) {continue;}
1156 }
1157
1158
1159 if ((changed_symmetry =
1160 get_conventional_symmetry(change_of_basis_monocli[i],
1161 PRIMITIVE,
1162 symmetry)) == NULL) {
1163 goto err;
1164 }
1165
1166 is_found = hal_match_hall_symbol_db(origin_shift,
1167 changed_lattice,
1168 hall_number,
1169 changed_centering,
1170 changed_symmetry,
1171 symprec);
1172 sym_free_symmetry(changed_symmetry);
1173 changed_symmetry = NULL;
1174 if (is_found) {
1175 mat_copy_matrix_d3(lattice, changed_lattice);
1176 return 1;
1177 }
1178 }
1179
1180 err:
1181 return 0;
1182 }
1183
1184 /* Return 0 if failed */
match_hall_symbol_db_ortho(double origin_shift[3],double lattice[3][3],const int hall_number,const Centering centering,const Symmetry * symmetry,const int num_free_axes,const double symprec)1185 static int match_hall_symbol_db_ortho(double origin_shift[3],
1186 double lattice[3][3],
1187 const int hall_number,
1188 const Centering centering,
1189 const Symmetry *symmetry,
1190 const int num_free_axes,
1191 const double symprec)
1192 {
1193 int i, j, k, l, is_found;
1194 double vec[3], norms[3];
1195 Centering changed_centering;
1196 Symmetry * changed_symmetry;
1197 double changed_lattice[3][3];
1198
1199 changed_symmetry = NULL;
1200
1201 for (i = 0; i < 6; i++) {
1202 if (centering == C_FACE) {
1203 changed_centering = change_of_centering_ortho[i];
1204 } else {
1205 changed_centering = centering;
1206 }
1207
1208 mat_multiply_matrix_d3(changed_lattice,
1209 lattice,
1210 change_of_basis_ortho[i]);
1211
1212 if (num_free_axes == 2) {
1213 l = 0;
1214 for (j = 0; j < 3; j++) {
1215 if (j == change_of_unique_axis_ortho[i]) {continue;}
1216 for (k = 0; k < 3; k++) {vec[k] = changed_lattice[k][j];}
1217 norms[l] = mat_norm_squared_d3(vec);
1218 l++;
1219 }
1220 if (norms[0] > norms[1]) {continue;}
1221 }
1222
1223 if (num_free_axes == 3) {
1224 for (j = 0; j < 3; j++) {
1225 for (k = 0; k < 3; k++) {vec[k] = changed_lattice[k][j];}
1226 norms[j] = mat_norm_squared_d3(vec);
1227 }
1228 if (norms[0] > norms[1] || norms[0] > norms[2]) {continue;}
1229 }
1230
1231 if (num_free_axes == 6) {
1232 for (j = 0; j < 3; j++) {
1233 for (k = 0; k < 3; k++) {vec[k] = changed_lattice[k][j];}
1234 norms[j] = mat_norm_squared_d3(vec);
1235 }
1236 if (norms[0] > norms[1] || norms[1] > norms[2]) {continue;}
1237 }
1238
1239 if ((changed_symmetry = get_conventional_symmetry(change_of_basis_ortho[i],
1240 PRIMITIVE,
1241 symmetry)) == NULL) {
1242 goto err;
1243 }
1244
1245 is_found = hal_match_hall_symbol_db(origin_shift,
1246 changed_lattice,
1247 hall_number,
1248 changed_centering,
1249 changed_symmetry,
1250 symprec);
1251 sym_free_symmetry(changed_symmetry);
1252 changed_symmetry = NULL;
1253 if (is_found) {
1254 mat_copy_matrix_d3(lattice, changed_lattice);
1255 return 1;
1256 }
1257 }
1258
1259 err:
1260 return 0;
1261 }
1262
1263 /* Return NULL if failed */
get_conventional_symmetry(SPGCONST double transform_mat[3][3],const Centering centering,const Symmetry * primitive_sym)1264 static Symmetry * get_conventional_symmetry(SPGCONST double transform_mat[3][3],
1265 const Centering centering,
1266 const Symmetry *primitive_sym)
1267 {
1268 int i, j, k, multi, size;
1269 double inv_tmat[3][3], shift[3][3];
1270 double symmetry_rot_d3[3][3], primitive_sym_rot_d3[3][3];
1271 Symmetry *symmetry;
1272
1273 symmetry = NULL;
1274
1275 size = primitive_sym->size;
1276
1277 switch (centering) {
1278 case FACE:
1279 if ((symmetry = sym_alloc_symmetry(size * 4)) == NULL) {
1280 return NULL;
1281 }
1282 break;
1283 case R_CENTER:
1284 if ((symmetry = sym_alloc_symmetry(size * 3)) == NULL) {
1285 return NULL;
1286 }
1287 break;
1288 case BODY:
1289 case A_FACE:
1290 case B_FACE:
1291 case C_FACE:
1292 if ((symmetry = sym_alloc_symmetry(size * 2)) == NULL) {
1293 return NULL;
1294 }
1295 break;
1296 default:
1297 if ((symmetry = sym_alloc_symmetry(size)) == NULL) {
1298 return NULL;
1299 }
1300 break;
1301 }
1302
1303 for (i = 0; i < size; i++) {
1304 mat_cast_matrix_3i_to_3d(primitive_sym_rot_d3, primitive_sym->rot[i]);
1305
1306 /* C*S*C^-1: recover conventional cell symmetry operation */
1307 mat_get_similar_matrix_d3(symmetry_rot_d3,
1308 primitive_sym_rot_d3,
1309 transform_mat,
1310 0);
1311 mat_cast_matrix_3d_to_3i(symmetry->rot[i], symmetry_rot_d3);
1312
1313 /* translation in conventional cell: C = B^-1*P */
1314 mat_inverse_matrix_d3(inv_tmat, transform_mat, 0);
1315 mat_multiply_matrix_vector_d3(symmetry->trans[i],
1316 inv_tmat,
1317 primitive_sym->trans[i]);
1318 }
1319
1320 multi = 1;
1321
1322 if (centering != PRIMITIVE) {
1323 multi = get_centering_shifts(shift, centering);
1324 for (i = 0; i < multi - 1; i++) {
1325 for (j = 0; j < size; j++) {
1326 mat_copy_matrix_i3(symmetry->rot[(i+1) * size + j],
1327 symmetry->rot[j]);
1328 for (k = 0; k < 3; k++) {
1329 symmetry->trans[(i+1) * size + j][k] =
1330 symmetry->trans[j][k] + shift[i][k];
1331 }
1332 }
1333 }
1334 }
1335
1336 for (i = 0; i < multi; i++) {
1337 for (j = 0; j < size; j++) {
1338 for (k = 0; k < 3; k++) {
1339 symmetry->trans[i * size + j][k] =
1340 mat_Dmod1(symmetry->trans[i * size + j][k]);
1341 }
1342 }
1343 }
1344
1345 return symmetry;
1346 }
1347
1348 /* Return CENTERING_ERROR if failed */
get_centering(double correction_mat[3][3],SPGCONST int transform_mat[3][3],const Laue laue)1349 static Centering get_centering(double correction_mat[3][3],
1350 SPGCONST int transform_mat[3][3],
1351 const Laue laue)
1352 {
1353 int det;
1354 double trans_corr_mat[3][3];
1355 Centering centering;
1356
1357 mat_copy_matrix_d3(correction_mat, identity);
1358 det = abs(mat_get_determinant_i3(transform_mat));
1359 debug_print("laue class: %d\n", laue);
1360 debug_print("multiplicity: %d\n", det);
1361
1362 switch (det) {
1363
1364 case 1:
1365 centering = PRIMITIVE;
1366 break;
1367
1368 case 2:
1369 centering = get_base_center(transform_mat);
1370 if (centering == A_FACE) {
1371 if (laue == LAUE2M) {
1372 debug_print("Monocli A to C\n");
1373 mat_copy_matrix_d3(correction_mat, monocli_a2c);
1374 } else {
1375 mat_copy_matrix_d3(correction_mat, a2c);
1376 }
1377 centering = C_FACE;
1378 }
1379 if (centering == B_FACE) {
1380 mat_copy_matrix_d3(correction_mat, b2c);
1381 centering = C_FACE;
1382 }
1383 if (laue == LAUE2M && centering == BODY) {
1384 debug_print("Monocli I to C\n");
1385 mat_copy_matrix_d3(correction_mat, monocli_i2c);
1386 centering = C_FACE;
1387 }
1388 break;
1389
1390 case 3:
1391 /* hP (a=b) but not hR (a=b=c) */
1392 centering = R_CENTER;
1393 mat_multiply_matrix_id3(trans_corr_mat, transform_mat, rhombo_obverse);
1394 if (mat_is_int_matrix(trans_corr_mat, INT_PREC)) {
1395 mat_copy_matrix_d3(correction_mat, rhombo_obverse);
1396 debug_print("R-center observe setting\n");
1397 debug_print_matrix_d3(trans_corr_mat);
1398 }
1399 mat_multiply_matrix_id3(trans_corr_mat, transform_mat, rhomb_reverse);
1400 if (mat_is_int_matrix(trans_corr_mat, INT_PREC)) {
1401 mat_copy_matrix_d3(correction_mat, rhomb_reverse);
1402 debug_print("R-center reverse setting\n");
1403 debug_print_matrix_d3(trans_corr_mat);
1404 }
1405 break;
1406
1407 case 4:
1408 centering = FACE;
1409 break;
1410
1411 default:
1412 centering = CENTERING_ERROR;
1413 break;
1414 }
1415
1416 return centering;
1417 }
1418
get_base_center(SPGCONST int transform_mat[3][3])1419 static Centering get_base_center(SPGCONST int transform_mat[3][3])
1420 {
1421 int i;
1422 Centering centering = PRIMITIVE;
1423
1424 debug_print("lat_get_base_center\n");
1425
1426 /* C center */
1427 for (i = 0; i < 3; i++) {
1428 if (transform_mat[i][0] == 0 &&
1429 transform_mat[i][1] == 0 &&
1430 abs(transform_mat[i][2]) == 1) {
1431 centering = C_FACE;
1432 goto end;
1433 }
1434 }
1435
1436 /* A center */
1437 for (i = 0; i < 3; i++) {
1438 if (abs(transform_mat[i][0]) == 1 &&
1439 transform_mat[i][1] == 0 &&
1440 transform_mat[i][2] == 0) {
1441 centering = A_FACE;
1442 goto end;
1443 }
1444 }
1445
1446 /* B center */
1447 for (i = 0; i < 3; i++) {
1448 if (transform_mat[i][0] == 0 &&
1449 abs(transform_mat[i][1]) == 1 &&
1450 transform_mat[i][2] == 0) {
1451 centering = B_FACE;
1452 goto end;
1453 }
1454 }
1455
1456 /* body center */
1457 if ((abs(transform_mat[0][0]) +
1458 abs(transform_mat[0][1]) +
1459 abs(transform_mat[0][2]) == 2) &&
1460 (abs(transform_mat[1][0]) +
1461 abs(transform_mat[1][1]) +
1462 abs(transform_mat[1][2]) == 2) &&
1463 (abs(transform_mat[2][0]) +
1464 abs(transform_mat[2][1]) +
1465 abs(transform_mat[2][2]) == 2)) {
1466 centering = BODY;
1467 goto end;
1468 }
1469
1470 /* This should not happen. */
1471 warning_print("spglib: No centring was found (line %d, %s).\n", __LINE__, __FILE__);
1472 return PRIMITIVE;
1473
1474 end:
1475 return centering;
1476 }
1477
get_centering_shifts(double shift[3][3],const Centering centering)1478 static int get_centering_shifts(double shift[3][3],
1479 const Centering centering)
1480 {
1481 int i, j, multi;
1482
1483 multi = 1;
1484 for (i = 0; i < 3; i++) {
1485 for (j = 0; j < 3; j++) {
1486 shift[i][j] = 0;
1487 }
1488 }
1489
1490 if (centering != PRIMITIVE) {
1491 if (centering != FACE && centering != R_CENTER) {
1492 for (i = 0; i < 3; i++) { shift[0][i] = 0.5; } /* BASE */
1493 if (centering == A_FACE) { shift[0][0] = 0; }
1494 if (centering == B_FACE) { shift[0][1] = 0; }
1495 if (centering == C_FACE) { shift[0][2] = 0; }
1496 multi = 2;
1497 }
1498
1499 if (centering == R_CENTER) {
1500 shift[0][0] = 2. / 3;
1501 shift[0][1] = 1. / 3;
1502 shift[0][2] = 1. / 3;
1503 shift[1][0] = 1. / 3;
1504 shift[1][1] = 2. / 3;
1505 shift[1][2] = 2. / 3;
1506 multi = 3;
1507 }
1508
1509 if (centering == FACE) {
1510 shift[0][0] = 0;
1511 shift[0][1] = 0.5;
1512 shift[0][2] = 0.5;
1513 shift[1][0] = 0.5;
1514 shift[1][1] = 0;
1515 shift[1][2] = 0.5;
1516 shift[2][0] = 0.5;
1517 shift[2][1] = 0.5;
1518 shift[2][2] = 0;
1519 multi = 4;
1520 }
1521 }
1522
1523 return multi;
1524 }
1525