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 spglib 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 <math.h>
39 #include "cell.h"
40 #include "delaunay.h"
41 #include "hall_symbol.h"
42 #include "mathfunc.h"
43 #include "niggli.h"
44 #include "pointgroup.h"
45 #include "primitive.h"
46 #include "spacegroup.h"
47 #include "spg_database.h"
48 #include "symmetry.h"
49
50 #include "debug.h"
51
52 #define REDUCE_RATE 0.95
53 #define NUM_ATTEMPT 100
54 #define INT_PREC 0.1
55 #define ZERO_PREC 1e-10
56
57 static double change_of_basis_monocli[36][3][3] = {
58 {{ 1, 0, 0 }, /* b first turn; two axes are flipped in second turn */
59 { 0, 1, 0 },
60 { 0, 0, 1 }},
61 {{ 0, 0, 1 }, /* b */
62 { 0,-1, 0 },
63 { 1, 0, 0 }},
64 {{ 0, 0, 1 }, /* a */
65 { 1, 0, 0 },
66 { 0, 1, 0 }},
67 {{ 1, 0, 0 }, /* c */
68 { 0, 0, 1 },
69 { 0,-1, 0 }},
70 {{ 0, 1, 0 }, /* c */
71 { 0, 0, 1 },
72 { 1, 0, 0 }},
73 {{ 0,-1, 0 }, /* a */
74 { 1, 0, 0 },
75 { 0, 0, 1 }},
76 {{-1, 0, 1 }, /* b */
77 { 0, 1, 0 },
78 {-1, 0, 0 }},
79 {{ 1, 0,-1 }, /* b */
80 { 0,-1, 0 },
81 { 0, 0,-1 }},
82 {{ 0, 1,-1 }, /* a */
83 { 1, 0, 0 },
84 { 0, 0,-1 }},
85 {{-1,-1, 0 }, /* c */
86 { 0, 0, 1 },
87 {-1, 0, 0 }},
88 {{ 1,-1, 0 }, /* c */
89 { 0, 0, 1 },
90 { 0,-1, 0 }},
91 {{ 0, 1, 1 }, /* a */
92 { 1, 0, 0 },
93 { 0, 1, 0 }},
94 {{ 0, 0,-1 }, /* b */
95 { 0, 1, 0 },
96 { 1, 0,-1 }},
97 {{-1, 0, 0 }, /* b */
98 { 0,-1, 0 },
99 {-1, 0, 1 }},
100 {{ 0,-1, 0 }, /* a */
101 { 1, 0, 0 },
102 { 0,-1, 1 }},
103 {{ 0, 1, 0 }, /* c */
104 { 0, 0, 1 },
105 { 1, 1, 0 }},
106 {{-1, 0, 0 }, /* c */
107 { 0, 0, 1 },
108 {-1, 1, 0 }},
109 {{ 0, 0,-1 }, /* a */
110 { 1, 0, 0 },
111 { 0,-1,-1 }},
112 {{ 1, 0, 0 }, /* b two axes are flipped to look for non-acute axes */
113 { 0,-1, 0 },
114 { 0, 0,-1 }},
115 {{ 0, 0,-1 }, /* b */
116 { 0, 1, 0 },
117 { 1, 0, 0 }},
118 {{ 0, 0, 1 }, /* a */
119 {-1, 0, 0 },
120 { 0,-1, 0 }},
121 {{-1, 0, 0 }, /* c */
122 { 0, 0,-1 },
123 { 0,-1, 0 }},
124 {{ 0, 1, 0 }, /* c */
125 { 0, 0,-1 },
126 {-1, 0, 0 }},
127 {{ 0, 1, 0 }, /* a */
128 {-1, 0, 0 },
129 { 0, 0, 1 }},
130 {{-1, 0,-1 }, /* b */
131 { 0,-1, 0 },
132 {-1, 0, 0 }},
133 {{ 1, 0, 1 }, /* b */
134 { 0, 1, 0 },
135 { 0, 0, 1 }},
136 {{ 0,-1,-1 }, /* a */
137 {-1, 0, 0 },
138 { 0, 0,-1 }},
139 {{ 1,-1, 0 }, /* c */
140 { 0, 0,-1 },
141 { 1, 0, 0 }},
142 {{-1,-1, 0 }, /* c */
143 { 0, 0,-1 },
144 { 0,-1, 0 }},
145 {{ 0,-1, 1 }, /* a */
146 {-1, 0, 0 },
147 { 0,-1, 0 }},
148 {{ 0, 0, 1 }, /* b */
149 { 0,-1, 0 },
150 { 1, 0, 1 }},
151 {{-1, 0, 0 }, /* b */
152 { 0, 1, 0 },
153 {-1, 0,-1 }},
154 {{ 0, 1, 0 }, /* a */
155 {-1, 0, 0 },
156 { 0, 1, 1 }},
157 {{ 0, 1, 0 }, /* c */
158 { 0, 0,-1 },
159 {-1, 1, 0 }},
160 {{ 1, 0, 0 }, /* c */
161 { 0, 0,-1 },
162 { 1, 1, 0 }},
163 {{ 0, 0,-1 }, /* a */
164 {-1, 0, 0 },
165 { 0, 1,-1 }}};
166
167 static Centering change_of_centering_monocli[36] = {
168 C_FACE, /* first turn */
169 A_FACE,
170 B_FACE,
171 B_FACE,
172 A_FACE,
173 C_FACE,
174 A_FACE,
175 C_FACE,
176 C_FACE,
177 A_FACE,
178 B_FACE,
179 B_FACE,
180 BODY,
181 BODY,
182 BODY,
183 BODY,
184 BODY,
185 BODY,
186 C_FACE, /* second turn */
187 A_FACE,
188 B_FACE,
189 B_FACE,
190 A_FACE,
191 C_FACE,
192 A_FACE,
193 C_FACE,
194 C_FACE,
195 A_FACE,
196 B_FACE,
197 B_FACE,
198 BODY,
199 BODY,
200 BODY,
201 BODY,
202 BODY,
203 BODY};
204
205 static int change_of_unique_axis_monocli[36] = {
206 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0,
207 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0, 1, 1, 0, 2, 2, 0};
208
209 static double change_of_basis_ortho[6][3][3] = {{{ 1, 0, 0 },
210 { 0, 1, 0 },
211 { 0, 0, 1 }},
212 {{ 0, 0, 1 },
213 { 1, 0, 0 },
214 { 0, 1, 0 }},
215 {{ 0, 1, 0 },
216 { 0, 0, 1 },
217 { 1, 0, 0 }},
218 {{ 0, 1, 0 },
219 { 1, 0, 0 },
220 { 0, 0,-1 }},
221 {{ 1, 0, 0 },
222 { 0, 0, 1 },
223 { 0,-1, 0 }},
224 {{ 0, 0, 1 },
225 { 0, 1, 0 },
226 {-1, 0, 0 }}};
227
228 static Centering change_of_centering_ortho[6] = {C_FACE,
229 B_FACE,
230 A_FACE,
231 C_FACE,
232 B_FACE,
233 A_FACE};
234 static int change_of_unique_axis_ortho[6] = {2, 1, 0, 2, 1, 0};
235
236 /* n_l : the index of L(g) in N_\epsilon(g) of SPG No.16-74 */
237 /* See ITA: Affine normalizer or highest symmetry Euclidean normalizer */
238 /* Previous implementation below was not correct for 67, 68, 73, 74, */
239 /* Cmma, Ccca, Ibca, Imma */
240 /* 6 / ((Number of hall symbols of each spg-type) / x) */
241 /* where x=1 and x=2 with without and with centring. */
242 /* 6, 2, 2, 6, 2, */
243 /* 2, 6, 6, 6, 2, 1, 2, 1, 1, 1, */
244 /* 1, 2, 1, 2, 2, 1, 2, 1, 1, 1, */
245 /* 1, 2, 2, 2, 2, 1, 6, 6, 2, 2, */
246 /* 1, 1, 1, 1, 2, 2, 1, 2, 2, 1, */
247 /* 3, 1, 1, 1, 2, 2, 1, 1, 6, 6, */
248 /* 6, 2, 3, 1 */
249 static int num_axis_choices_ortho[59] = {
250 6, 2, 2, 6, 2, /* 16-20 */
251 2, 6, 6, 6, 2, 1, 2, 1, 1, 1, /* 21-30 */
252 1, 2, 1, 2, 2, 1, 2, 1, 1, 1, /* 31-40 */
253 1, 2, 2, 2, 2, 1, 6, 6, 2, 2, /* 41-50 */
254 1, 1, 1, 1, 2, 2, 1, 2, 2, 1, /* 51-60 */
255 3, 1, 1, 1, 2, 2, 2, 2, 6, 6, /* 61-70 */
256 6, 2, 6, 2}; /* 71-74 */
257
258 /* hR_to_hP is the inverse of rhombo_obverse and used until the */
259 /* commit e0398ba9. But now this is included in */
260 /* change_of_basis_rhombo in a way of */
261 /* dot(change_of_basis_rhombo_hex,hR_to_hP). */
262 /* static double hR_to_hP[3][3] = {{ 1, 0, 1 }, */
263 /* {-1, 1, 1 }, */
264 /* { 0,-1, 1 }}; */
265
266 static double change_of_basis_rhombo[6][3][3] = {{{ 1, 0, 0 },
267 { 0, 1, 0 },
268 { 0, 0, 1 }},
269 {{ 0, 0, 1 },
270 { 1, 0, 0 },
271 { 0, 1, 0 }},
272 {{ 0, 1, 0 },
273 { 0, 0, 1 },
274 { 1, 0, 0 }},
275 {{ 0, 0,-1 },
276 { 0,-1, 0 },
277 {-1, 0, 0 }},
278 {{ 0,-1, 0 },
279 {-1, 0, 0 },
280 { 0, 0,-1 }},
281 {{-1, 0, 0 },
282 { 0, 0,-1 },
283 { 0,-1, 0 }}};
284
285 static double change_of_basis_rhombo_hex[6][3][3] = {{{ 1, 0, 1},
286 {-1, 1, 1},
287 { 0, -1, 1}},
288 {{ 0, -1, 1},
289 { 1, 0, 1},
290 {-1, 1, 1}},
291 {{-1, 1, 1},
292 { 0, -1, 1},
293 { 1, 0, 1}},
294 {{ 0, 1, -1},
295 { 1, -1, -1},
296 {-1, 0, -1}},
297 {{ 1, -1, -1},
298 {-1, 0, -1},
299 { 0, 1, -1}},
300 {{-1, 0, -1},
301 { 0, 1, -1},
302 { 1, -1, -1}}};
303
304 static double change_of_basis_C4[8][3][3] = {{{ 1, 0, 0 },
305 { 0, 1, 0 },
306 { 0, 0, 1 }},
307 {{ 0,-1, 0 },
308 { 1, 0, 0 },
309 { 0, 0, 1 }},
310 {{-1, 0, 0 },
311 { 0,-1, 0 },
312 { 0, 0, 1 }},
313 {{ 0, 1, 0 },
314 {-1, 0, 0 },
315 { 0, 0, 1 }},
316 {{ 0, 1, 0 },
317 { 1, 0, 0 },
318 { 0, 0,-1 }},
319 {{-1, 0, 0 },
320 { 0, 1, 0 },
321 { 0, 0,-1 }},
322 {{ 0,-1, 0 },
323 {-1, 0, 0 },
324 { 0, 0,-1 }},
325 {{ 1, 0, 0 },
326 { 0,-1, 0 },
327 { 0, 0,-1 }}};
328
329 static double change_of_basis_C6[12][3][3] = {{{ 1, 0, 0 }, /* 1 */
330 { 0, 1, 0 },
331 { 0, 0, 1 }},
332 {{ 1,-1, 0 }, /* 2 */
333 { 1, 0, 0 },
334 { 0, 0, 1 }},
335 {{ 0,-1, 0 }, /* 3 */
336 { 1,-1, 0 },
337 { 0, 0, 1 }},
338 {{-1, 0, 0 }, /* 4 */
339 { 0,-1, 0 },
340 { 0, 0, 1 }},
341 {{-1, 1, 0 }, /* 5 */
342 {-1, 0, 0 },
343 { 0, 0, 1 }},
344 {{ 0, 1, 0 }, /* 6 */
345 {-1, 1, 0 },
346 { 0, 0, 1 }},
347 {{ 0, 1, 0 }, /* 7 */
348 { 1, 0, 0 },
349 { 0, 0,-1 }},
350 {{-1, 1, 0 }, /* 8 */
351 { 0, 1, 0 },
352 { 0, 0,-1 }},
353 {{-1, 0, 0 }, /* 9 */
354 {-1, 1, 0 },
355 { 0, 0,-1 }},
356 {{ 0,-1, 0 }, /* 10 */
357 {-1, 0, 0 },
358 { 0, 0,-1 }},
359 {{ 1,-1, 0 }, /* 11 */
360 { 0,-1, 0 },
361 { 0, 0,-1 }},
362 {{ 1, 0, 0 }, /* 12 */
363 { 1,-1, 0 },
364 { 0, 0,-1 }}};
365
366 /* Removed after commit 67997654 */
367 /* static double change_of_basis_501[3][3] = {{ 0, 0, 1}, */
368 /* { 0,-1, 0}, */
369 /* { 1, 0, 0}}; */
370
371 static int spacegroup_to_hall_number[230] = {
372 1, 2, 3, 6, 9, 18, 21, 30, 39, 57,
373 60, 63, 72, 81, 90, 108, 109, 112, 115, 116,
374 119, 122, 123, 124, 125, 128, 134, 137, 143, 149,
375 155, 161, 164, 170, 173, 176, 182, 185, 191, 197,
376 203, 209, 212, 215, 218, 221, 227, 228, 230, 233,
377 239, 245, 251, 257, 263, 266, 269, 275, 278, 284,
378 290, 292, 298, 304, 310, 313, 316, 322, 334, 335,
379 337, 338, 341, 343, 349, 350, 351, 352, 353, 354,
380 355, 356, 357, 358, 359, 361, 363, 364, 366, 367,
381 368, 369, 370, 371, 372, 373, 374, 375, 376, 377,
382 378, 379, 380, 381, 382, 383, 384, 385, 386, 387,
383 388, 389, 390, 391, 392, 393, 394, 395, 396, 397,
384 398, 399, 400, 401, 402, 404, 406, 407, 408, 410,
385 412, 413, 414, 416, 418, 419, 420, 422, 424, 425,
386 426, 428, 430, 431, 432, 433, 435, 436, 438, 439,
387 440, 441, 442, 443, 444, 446, 447, 448, 449, 450,
388 452, 454, 455, 456, 457, 458, 460, 462, 463, 464,
389 465, 466, 467, 468, 469, 470, 471, 472, 473, 474,
390 475, 476, 477, 478, 479, 480, 481, 482, 483, 484,
391 485, 486, 487, 488, 489, 490, 491, 492, 493, 494,
392 495, 497, 498, 500, 501, 502, 503, 504, 505, 506,
393 507, 508, 509, 510, 511, 512, 513, 514, 515, 516,
394 517, 518, 520, 521, 523, 524, 525, 527, 529, 530,
395 };
396
397 static double identity[3][3] = {{ 1, 0, 0 },
398 { 0, 1, 0 },
399 { 0, 0, 1 }};
400 static double monocli_i2c[3][3] = {{ 1, 0,-1 },
401 { 0, 1, 0 },
402 { 1, 0, 0 }};
403 static double monocli_a2c[3][3] = {{ 0, 0, 1 },
404 { 0,-1, 0 },
405 { 1, 0, 0 }};
406 static double rhombo_obverse[3][3] = {{ 2./3,-1./3,-1./3 },
407 { 1./3, 1./3,-2./3 },
408 { 1./3, 1./3, 1./3 }};
409 static double rhomb_reverse[3][3] = {{ 1./3,-2./3, 1./3 },
410 { 2./3,-1./3,-1./3 },
411 { 1./3, 1./3, 1./3 }};
412 static double a2c[3][3] = {{ 0, 0, 1 },
413 { 1, 0, 0 },
414 { 0, 1, 0 }};
415 static double b2c[3][3] = {{ 0, 1, 0 },
416 { 0, 0, 1 },
417 { 1, 0, 0 }};
418
419 static double A_mat[3][3] = {{ 1, 0, 0},
420 { 0, 1./2,-1./2},
421 { 0, 1./2, 1./2}};
422 static double C_mat[3][3] = {{ 1./2, 1./2, 0},
423 {-1./2, 1./2, 0},
424 { 0, 0, 1}};
425 static double R_mat[3][3] = {{ 2./3,-1./3,-1./3 },
426 { 1./3, 1./3,-2./3 },
427 { 1./3, 1./3, 1./3 }};
428 static double I_mat[3][3] = {{-1./2, 1./2, 1./2 },
429 { 1./2,-1./2, 1./2 },
430 { 1./2, 1./2,-1./2 }};
431 static double F_mat[3][3] = {{ 0, 1./2, 1./2 },
432 { 1./2, 0, 1./2 },
433 { 1./2, 1./2, 0 }};
434
435 static Spacegroup * search_spacegroup_with_symmetry(const Primitive * primitive,
436 const int candidates[],
437 const int num_candidates,
438 const Symmetry *symmetry,
439 const double symprec,
440 const double angle_tolerance);
441 static Spacegroup * get_spacegroup(const int hall_number,
442 const double origin_shift[3],
443 SPGCONST double conv_lattice[3][3]);
444 static int iterative_search_hall_number(double origin_shift[3],
445 double conv_lattice[3][3],
446 const int candidates[],
447 const int num_candidates,
448 const Primitive * primitive,
449 const Symmetry * symmetry,
450 const double symprec,
451 const double angle_tolerance);
452 static int change_basis_tricli(int tmat_int[3][3],
453 SPGCONST double conv_lattice[3][3],
454 SPGCONST double primitive_lattice[3][3],
455 const double symprec);
456 static int change_basis_monocli(int tmat_int[3][3],
457 SPGCONST double conv_lattice[3][3],
458 SPGCONST double primitive_lattice[3][3],
459 const double symprec);
460 static Symmetry *
461 get_initial_conventional_symmetry(const Centering centering,
462 SPGCONST double tmat[3][3],
463 const Symmetry * symmetry);
464 static int search_hall_number(double origin_shift[3],
465 double conv_lattice[3][3],
466 const int candidates[],
467 const int num_candidates,
468 const Primitive *primitive,
469 const Symmetry * symmetry,
470 const double symprec);
471 static int match_hall_symbol_db(double origin_shift[3],
472 double conv_lattice[3][3],
473 SPGCONST double (*orig_lattice)[3],
474 const int hall_number,
475 const int pointgroup_number,
476 const Holohedry holohedry,
477 const Centering centering,
478 const Symmetry *symmetry,
479 const double symprec);
480 static int match_hall_symbol_db_monocli(double origin_shift[3],
481 double conv_lattice[3][3],
482 SPGCONST double (*orig_lattice)[3],
483 const int hall_number,
484 const int space_group_number,
485 const Centering centering,
486 const Symmetry *conv_symmetry,
487 const double symprec);
488 static int match_hall_symbol_db_monocli_in_loop(double origin_shift[3],
489 double conv_lattice[3][3],
490 double norms_squared[2],
491 const int i,
492 SPGCONST double (*orig_lattice)[3],
493 const int check_norms,
494 const int hall_number,
495 const Centering centering,
496 const Symmetry *conv_symmetry,
497 const double symprec);
498 static int match_hall_symbol_db_ortho(double origin_shift[3],
499 double conv_lattice[3][3],
500 SPGCONST double (*orig_lattice)[3],
501 const int hall_number,
502 const Centering centering,
503 const Symmetry *symmetry,
504 const int num_free_axes,
505 const double symprec);
506 static int match_hall_symbol_db_ortho_in_loop(double origin_shift[3],
507 double lattice[3][3],
508 SPGCONST double (*orig_lattice)[3],
509 const int i,
510 const int hall_number,
511 const Centering centering,
512 const Symmetry *symmetry,
513 const int num_free_axes,
514 const double symprec);
515 static int match_hall_symbol_db_cubic(double origin_shift[3],
516 double conv_lattice[3][3],
517 SPGCONST double (*orig_lattice)[3],
518 const int hall_number,
519 const Centering centering,
520 const Symmetry *conv_symmetry,
521 const double symprec);
522 static int match_hall_symbol_db_cubic_in_loop(double origin_shift[3],
523 double conv_lattice[3][3],
524 SPGCONST double (*orig_lattice)[3],
525 const int i,
526 const int hall_number,
527 const Centering centering,
528 const Symmetry *conv_symmetry,
529 const double symprec);
530 static int match_hall_symbol_db_rhombo(double origin_shift[3],
531 double conv_lattice[3][3],
532 SPGCONST double (*orig_lattice)[3],
533 const int hall_number,
534 const Symmetry *conv_symmetry,
535 const double symprec);
536 static int match_hall_symbol_db_others(double origin_shift[3],
537 double conv_lattice[3][3],
538 SPGCONST double (*orig_lattice)[3],
539 const int hall_number,
540 const Centering centering,
541 const Holohedry holohedry,
542 const Symmetry *conv_symmetry,
543 const double symprec);
544 static int match_hall_symbol_db_change_of_basis_loop(
545 double origin_shift[3],
546 double conv_lattice[3][3],
547 SPGCONST double (*orig_lattice)[3],
548 SPGCONST double (*change_of_basis)[3][3],
549 const int num_change_of_basis,
550 const int hall_number,
551 const Centering centering,
552 const Symmetry *conv_symmetry,
553 const double symprec);
554 static Symmetry * get_conventional_symmetry(SPGCONST double tmat[3][3],
555 const Centering centering,
556 const Symmetry *primitive_sym);
557 static Centering get_centering(double correction_mat[3][3],
558 SPGCONST int tmat[3][3],
559 const Laue laue);
560 static Centering get_base_center(SPGCONST int tmat[3][3]);
561 static int get_centering_shifts(double shift[3][3],
562 const Centering centering);
563 static int is_equivalent_lattice(double tmat[3][3],
564 const int allow_flip,
565 SPGCONST double lattice[3][3],
566 SPGCONST double orig_lattice[3][3],
567 const double symprec);
568
569 /* Return spacegroup.number = 0 if failed */
spa_search_spacegroup(const Primitive * primitive,const int hall_number,const double symprec,const double angle_tolerance)570 Spacegroup * spa_search_spacegroup(const Primitive * primitive,
571 const int hall_number,
572 const double symprec,
573 const double angle_tolerance)
574 {
575 Spacegroup *spacegroup;
576 Symmetry *symmetry;
577 int candidate[1];
578
579 debug_print("search_spacegroup (tolerance = %f):\n", symprec);
580
581 symmetry = NULL;
582 spacegroup = NULL;
583
584 if ((symmetry =
585 sym_get_operation(primitive->cell, symprec, angle_tolerance)) == NULL) {
586 return NULL;
587 }
588
589 if (hall_number > 0) {
590 candidate[0] = hall_number;
591 }
592
593 if (hall_number) {
594 spacegroup = search_spacegroup_with_symmetry(primitive,
595 candidate,
596 1,
597 symmetry,
598 symprec,
599 angle_tolerance);
600 } else {
601 spacegroup = search_spacegroup_with_symmetry(primitive,
602 spacegroup_to_hall_number,
603 230,
604 symmetry,
605 symprec,
606 angle_tolerance);
607 }
608
609 sym_free_symmetry(symmetry);
610 symmetry = NULL;
611
612 return spacegroup;
613 }
614
615 /* Retrun 0 if failed */
spa_search_spacegroup_with_symmetry(const Symmetry * symmetry,const double symprec)616 int spa_search_spacegroup_with_symmetry(const Symmetry *symmetry,
617 const double symprec)
618 {
619 int i, hall_number;
620 Spacegroup *spacegroup;
621 Primitive *primitive;
622
623 spacegroup = NULL;
624
625 if ((primitive = prm_alloc_primitive(1)) == NULL) {
626 return 0;
627 }
628 if ((primitive->cell = cel_alloc_cell(1)) == NULL) {
629 return 0;
630 }
631 mat_copy_matrix_d3(primitive->cell->lattice, identity);
632 for (i = 0; i < 3; i++) {
633 primitive->cell->position[0][i] = 0;
634 }
635 spacegroup = search_spacegroup_with_symmetry(primitive,
636 spacegroup_to_hall_number,
637 230,
638 symmetry,
639 symprec,
640 -1.0);
641 prm_free_primitive(primitive);
642 primitive = NULL;
643
644 if (spacegroup != NULL) {
645 hall_number = spacegroup->hall_number;
646 free(spacegroup);
647 spacegroup = NULL;
648 return hall_number;
649 } else {
650 return 0;
651 }
652 }
653
654 /* 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)655 Cell * spa_transform_to_primitive(int * mapping_table,
656 const Cell * cell,
657 SPGCONST double trans_mat[3][3],
658 const Centering centering,
659 const double symprec)
660 {
661 double tmat[3][3], tmat_inv[3][3], prim_lat[3][3];
662 Cell * primitive;
663
664 primitive = NULL;
665
666 if (!mat_inverse_matrix_d3(tmat_inv, trans_mat, symprec)) {
667 goto err;
668 }
669
670 switch (centering) {
671 case PRIMITIVE:
672 mat_copy_matrix_d3(tmat, tmat_inv);
673 break;
674 case A_FACE:
675 mat_multiply_matrix_d3(tmat, tmat_inv, A_mat);
676 break;
677 case C_FACE:
678 mat_multiply_matrix_d3(tmat, tmat_inv, C_mat);
679 break;
680 case FACE:
681 mat_multiply_matrix_d3(tmat, tmat_inv, F_mat);
682 break;
683 case BODY:
684 mat_multiply_matrix_d3(tmat, tmat_inv, I_mat);
685 break;
686 case R_CENTER:
687 mat_multiply_matrix_d3(tmat, tmat_inv, R_mat);
688 break;
689 default:
690 goto err;
691 }
692
693 mat_multiply_matrix_d3(prim_lat, cell->lattice, tmat);
694 if ((primitive = cel_trim_cell(mapping_table, prim_lat, cell, symprec))
695 == NULL) {
696 warning_print("spglib: cel_trim_cell failed.");
697 warning_print(" (line %d, %s).\n", __LINE__, __FILE__);
698 }
699
700 return primitive;
701
702 err:
703 return NULL;
704 }
705
706 /* Return NULL if failed */
spa_transform_from_primitive(const Cell * primitive,const Centering centering,const double symprec)707 Cell * spa_transform_from_primitive(const Cell * primitive,
708 const Centering centering,
709 const double symprec)
710 {
711 int multi, i, j, k, num_atom;
712 int *mapping_table;
713 double tmat[3][3], inv_tmat[3][3], shift[3][3];
714 Cell *std_cell, *trimmed_cell;
715
716 mapping_table = NULL;
717 trimmed_cell = NULL;
718 std_cell = NULL;
719
720 switch (centering) {
721 case PRIMITIVE:
722 break;
723 case A_FACE:
724 mat_copy_matrix_d3(tmat, A_mat);
725 mat_inverse_matrix_d3(inv_tmat, A_mat, 0);
726 break;
727 case C_FACE:
728 mat_copy_matrix_d3(tmat, C_mat);
729 mat_inverse_matrix_d3(inv_tmat, C_mat, 0);
730 break;
731 case FACE:
732 mat_copy_matrix_d3(tmat, F_mat);
733 mat_inverse_matrix_d3(inv_tmat, F_mat, 0);
734 break;
735 case BODY:
736 mat_copy_matrix_d3(tmat, I_mat);
737 mat_inverse_matrix_d3(inv_tmat, I_mat, 0);
738 break;
739 case R_CENTER:
740 mat_copy_matrix_d3(tmat, R_mat);
741 mat_inverse_matrix_d3(inv_tmat, R_mat, 0);
742 break;
743 default:
744 goto ret;
745 }
746
747 multi = get_centering_shifts(shift, centering);
748
749 if ((mapping_table = (int*) malloc(sizeof(int) * primitive->size * multi))
750 == NULL) {
751 warning_print("spglib: Memory could not be allocated ");
752 goto ret;
753 }
754
755 if ((std_cell = cel_alloc_cell(primitive->size * multi)) == NULL) {
756 free(mapping_table);
757 mapping_table = NULL;
758 goto ret;
759 }
760
761 mat_multiply_matrix_d3(std_cell->lattice, primitive->lattice, inv_tmat);
762
763 num_atom = 0;
764 for (i = 0; i < primitive->size; i++) {
765 mat_multiply_matrix_vector_d3(std_cell->position[num_atom],
766 tmat,
767 primitive->position[i]);
768 std_cell->types[num_atom] = primitive->types[i];
769 num_atom++;
770 }
771
772 for (i = 0; i < multi - 1; i++) {
773 for (j = 0; j < primitive->size; j++) {
774 mat_copy_vector_d3(std_cell->position[num_atom],
775 std_cell->position[j]);
776 for (k = 0; k < 3; k++) {
777 std_cell->position[num_atom][k] += shift[i][k];
778 }
779 std_cell->types[num_atom] = std_cell->types[j];
780 num_atom++;
781 }
782 }
783
784 trimmed_cell = cel_trim_cell(mapping_table,
785 std_cell->lattice,
786 std_cell,
787 symprec);
788 cel_free_cell(std_cell);
789 std_cell = NULL;
790 free(mapping_table);
791 mapping_table = NULL;
792
793 ret:
794 return trimmed_cell;
795 }
796
797 /* Return NULL if failed */
search_spacegroup_with_symmetry(const Primitive * primitive,const int candidates[],const int num_candidates,const Symmetry * symmetry,const double symprec,const double angle_tolerance)798 static Spacegroup * search_spacegroup_with_symmetry(const Primitive * primitive,
799 const int candidates[],
800 const int num_candidates,
801 const Symmetry *symmetry,
802 const double symprec,
803 const double angle_tolerance)
804 {
805 int hall_number;
806 double conv_lattice[3][3];
807 double origin_shift[3];
808 Spacegroup *spacegroup;
809 PointSymmetry pointsym;
810
811 debug_print("search_spacegroup (tolerance = %f):\n", symprec);
812
813 spacegroup = NULL;
814 origin_shift[0] = 0;
815 origin_shift[1] = 0;
816 origin_shift[2] = 0;
817
818 pointsym = ptg_get_pointsymmetry(symmetry->rot, symmetry->size);
819 if (pointsym.size < symmetry->size) {
820 warning_print("spglib: Point symmetry of primitive cell is broken. ");
821 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
822 return NULL;
823 }
824
825 hall_number = iterative_search_hall_number(origin_shift,
826 conv_lattice,
827 candidates,
828 num_candidates,
829 primitive,
830 symmetry,
831 symprec,
832 angle_tolerance);
833 if (hall_number == 0) {
834 return NULL;
835 }
836
837 spacegroup = get_spacegroup(hall_number, origin_shift, conv_lattice);
838
839 return spacegroup;
840 }
841
842 /* Return spacegroup.number = 0 if failed */
get_spacegroup(const int hall_number,const double origin_shift[3],SPGCONST double conv_lattice[3][3])843 static Spacegroup * get_spacegroup(const int hall_number,
844 const double origin_shift[3],
845 SPGCONST double conv_lattice[3][3])
846 {
847 Spacegroup *spacegroup;
848 SpacegroupType spacegroup_type;
849
850 spacegroup = NULL;
851
852 if ((spacegroup = (Spacegroup*) malloc(sizeof(Spacegroup))) == NULL) {
853 warning_print("spglib: Memory could not be allocated.");
854 return NULL;
855 }
856
857 if (0 < hall_number && hall_number < 531) {
858 spacegroup_type = spgdb_get_spacegroup_type(hall_number);
859 mat_copy_matrix_d3(spacegroup->bravais_lattice, conv_lattice);
860 mat_copy_vector_d3(spacegroup->origin_shift, origin_shift);
861 spacegroup->number = spacegroup_type.number;
862 spacegroup->hall_number = hall_number;
863 spacegroup->pointgroup_number = spacegroup_type.pointgroup_number;
864 memcpy(spacegroup->schoenflies, spacegroup_type.schoenflies, 7);
865 memcpy(spacegroup->hall_symbol, spacegroup_type.hall_symbol, 17);
866 memcpy(spacegroup->international, spacegroup_type.international, 32);
867 memcpy(spacegroup->international_long, spacegroup_type.international_full, 20);
868 memcpy(spacegroup->international_short,
869 spacegroup_type.international_short, 11);
870 memcpy(spacegroup->choice, spacegroup_type.choice, 6);
871 }
872
873 return spacegroup;
874 }
875
876 /* 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 Primitive * primitive,const Symmetry * symmetry,const double symprec,const double angle_tolerance)877 static int iterative_search_hall_number(double origin_shift[3],
878 double conv_lattice[3][3],
879 const int candidates[],
880 const int num_candidates,
881 const Primitive * primitive,
882 const Symmetry * symmetry,
883 const double symprec,
884 const double angle_tolerance)
885 {
886 int attempt, hall_number;
887 double tolerance;
888 Symmetry * sym_reduced;
889
890 debug_print("iterative_search_hall_number:\n");
891
892 hall_number = 0;
893 sym_reduced = NULL;
894
895 hall_number = search_hall_number(origin_shift,
896 conv_lattice,
897 candidates,
898 num_candidates,
899 primitive,
900 symmetry,
901 symprec);
902
903 if (hall_number > 0) {
904 goto ret;
905 }
906
907 tolerance = symprec;
908 for (attempt = 0; attempt < NUM_ATTEMPT; attempt++) {
909
910 warning_print("spglib: Attempt %d tolerance = %e failed",
911 attempt, tolerance);
912 warning_print("(line %d, %s).\n", __LINE__, __FILE__);
913
914 tolerance *= REDUCE_RATE;
915 sym_reduced = sym_reduce_operation(primitive->cell,
916 symmetry,
917 tolerance,
918 angle_tolerance);
919 hall_number = search_hall_number(origin_shift,
920 conv_lattice,
921 candidates,
922 num_candidates,
923 primitive,
924 sym_reduced,
925 symprec);
926 sym_free_symmetry(sym_reduced);
927 sym_reduced = NULL;
928 if (hall_number > 0) {
929 break;
930 }
931 }
932
933 ret:
934 return hall_number;
935 }
936
937 /* Return 0 if failed */
search_hall_number(double origin_shift[3],double conv_lattice[3][3],const int candidates[],const int num_candidates,const Primitive * primitive,const Symmetry * symmetry,const double symprec)938 static int search_hall_number(double origin_shift[3],
939 double conv_lattice[3][3],
940 const int candidates[],
941 const int num_candidates,
942 const Primitive * primitive,
943 const Symmetry * symmetry,
944 const double symprec)
945 {
946 int i, hall_number;
947 Centering centering;
948 Pointgroup pointgroup;
949 Symmetry * conv_symmetry;
950 int tmat_int[3][3];
951 double correction_mat[3][3], tmat[3][3], conv_lattice_tmp[3][3];
952
953 debug_print("search_hall_number:\n");
954
955 hall_number = 0;
956 conv_symmetry = NULL;
957
958 /* primitive->cell->lattice is set right handed. */
959 /* tmat_int is set right handed. */
960
961 /* For rhombohedral system, basis vectors for hP */
962 /* (hexagonal lattice basis) are provided by tmat_int, */
963 /* but may be either obverse or reverse setting. */
964 pointgroup = ptg_get_transformation_matrix(tmat_int,
965 symmetry->rot,
966 symmetry->size);
967
968 debug_print("[line %d, %s]\n", __LINE__, __FILE__);
969 debug_print("initial tranformation matrix\n");
970 debug_print_matrix_i3(tmat_int);
971
972 if (pointgroup.number == 0) {
973 goto err;
974 }
975
976 /* For LAUE1 and LAUE2M, update tmat_int by making smallest lattice */
977 if (pointgroup.laue == LAUE1 || pointgroup.laue == LAUE2M) {
978 mat_multiply_matrix_di3(
979 conv_lattice_tmp, primitive->cell->lattice, tmat_int);
980
981 if (pointgroup.laue == LAUE1) {
982 if (! change_basis_tricli(tmat_int,
983 conv_lattice_tmp,
984 primitive->cell->lattice,
985 symprec)) {
986 goto err;
987 }
988 }
989
990 /* The angle between non-unique axes may be acute or non-acute. */
991 if (pointgroup.laue == LAUE2M) {
992 if (! change_basis_monocli(tmat_int,
993 conv_lattice_tmp,
994 primitive->cell->lattice,
995 symprec)) {
996 goto err;
997 }
998 }
999 }
1000
1001 /* For rhombohedral system, correction matrix to hR */
1002 /* (a=b=c primitive lattice basis) is returned judging */
1003 /* either obverse or reverse setting. centering=R_CENTER. */
1004 if ((centering = get_centering(correction_mat,
1005 tmat_int,
1006 pointgroup.laue)) == CENTERING_ERROR) {
1007 goto err;
1008 }
1009
1010 mat_multiply_matrix_id3(tmat, tmat_int, correction_mat);
1011 mat_multiply_matrix_d3(conv_lattice, primitive->cell->lattice, tmat);
1012
1013 debug_print("[line %d, %s]\n", __LINE__, __FILE__);
1014 debug_print("tranformation matrix\n");
1015 debug_print_matrix_d3(tmat);
1016
1017 /* For rhombohedral system, symmetry for a=b=c primitive lattice */
1018 /* basis is returned although centering == R_CENTER. */
1019 if ((conv_symmetry = get_initial_conventional_symmetry(centering,
1020 tmat,
1021 symmetry)) == NULL) {
1022 goto err;
1023 }
1024
1025 for (i = 0; i < num_candidates; i++) {
1026 /* Check if hall_number is that of rhombohedral. */
1027 if (match_hall_symbol_db(origin_shift,
1028 conv_lattice, /* <-- modified only matched */
1029 primitive->orig_lattice,
1030 candidates[i],
1031 pointgroup.number,
1032 pointgroup.holohedry,
1033 centering,
1034 conv_symmetry,
1035 symprec)) {
1036
1037 debug_print("[line %d, %s]\n", __LINE__, __FILE__);
1038 debug_print("origin shift\n");
1039 debug_print_vector_d3(origin_shift);
1040
1041 hall_number = candidates[i];
1042 break;
1043 }
1044 }
1045
1046 sym_free_symmetry(conv_symmetry);
1047 conv_symmetry = NULL;
1048
1049 return hall_number;
1050
1051 err:
1052 return 0;
1053 }
1054
1055 /* Triclinic: Niggli cell reduction */
1056 /* Return 0 if failed */
change_basis_tricli(int tmat_int[3][3],SPGCONST double conv_lattice[3][3],SPGCONST double primitive_lattice[3][3],const double symprec)1057 static int change_basis_tricli(int tmat_int[3][3],
1058 SPGCONST double conv_lattice[3][3],
1059 SPGCONST double primitive_lattice[3][3],
1060 const double symprec)
1061 {
1062 int i, j;
1063 double niggli_cell[9];
1064 double smallest_lattice[3][3], inv_lattice[3][3], tmat[3][3];
1065
1066 for (i = 0; i < 3; i++) {
1067 for (j = 0; j < 3; j++) {
1068 niggli_cell[i * 3 + j] = conv_lattice[i][j];
1069 }
1070 }
1071
1072 if (! niggli_reduce(niggli_cell, symprec * symprec)) {
1073 return 0;
1074 }
1075
1076 for (i = 0; i < 3; i++) {
1077 for (j = 0; j < 3; j++) {
1078 smallest_lattice[i][j] = niggli_cell[i * 3 + j];
1079 }
1080 }
1081 if (mat_get_determinant_d3(smallest_lattice) < 0) {
1082 for (i = 0; i < 3; i++) {
1083 for (j = 0; j < 3; j++) {
1084 smallest_lattice[i][j] = -smallest_lattice[i][j];
1085 }
1086 }
1087 }
1088 mat_inverse_matrix_d3(inv_lattice, primitive_lattice, 0);
1089 mat_multiply_matrix_d3(tmat, inv_lattice, smallest_lattice);
1090 mat_cast_matrix_3d_to_3i(tmat_int, tmat);
1091
1092 return 1;
1093 }
1094
1095 /* Monoclinic: choose shortest a, c lattice vectors (|a| < |c|) */
1096 /* Return 0 if failed */
change_basis_monocli(int tmat_int[3][3],SPGCONST double conv_lattice[3][3],SPGCONST double primitive_lattice[3][3],const double symprec)1097 static int change_basis_monocli(int tmat_int[3][3],
1098 SPGCONST double conv_lattice[3][3],
1099 SPGCONST double primitive_lattice[3][3],
1100 const double symprec)
1101 {
1102 double smallest_lattice[3][3], inv_lattice[3][3], tmat[3][3];
1103
1104 if (! del_delaunay_reduce_2D(smallest_lattice,
1105 conv_lattice,
1106 1, /* unique axis of b */
1107 symprec)) {
1108 return 0;
1109 }
1110
1111 mat_inverse_matrix_d3(inv_lattice, primitive_lattice, 0);
1112 mat_multiply_matrix_d3(tmat, inv_lattice, smallest_lattice);
1113 mat_cast_matrix_3d_to_3i(tmat_int, tmat);
1114 return 1;
1115 }
1116
1117 /* Return NULL if failed */
1118 static Symmetry *
get_initial_conventional_symmetry(const Centering centering,SPGCONST double tmat[3][3],const Symmetry * symmetry)1119 get_initial_conventional_symmetry(const Centering centering,
1120 SPGCONST double tmat[3][3],
1121 const Symmetry * symmetry)
1122 {
1123 Symmetry * conv_symmetry;
1124
1125 debug_print("get_initial_conventional_symmetry\n");
1126
1127 conv_symmetry = NULL;
1128
1129 if (centering == R_CENTER) {
1130 /* hP for rhombohedral */
1131 conv_symmetry = get_conventional_symmetry(tmat,
1132 PRIMITIVE,
1133 symmetry);
1134 } else {
1135 conv_symmetry = get_conventional_symmetry(tmat,
1136 centering,
1137 symmetry);
1138 }
1139
1140 return conv_symmetry;
1141 }
1142
1143 /* Return 0 if failed */
match_hall_symbol_db(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int hall_number,const int pointgroup_number,const Holohedry holohedry,const Centering centering,const Symmetry * symmetry,const double symprec)1144 static int match_hall_symbol_db(double origin_shift[3],
1145 double conv_lattice[3][3],
1146 SPGCONST double (*orig_lattice)[3],
1147 const int hall_number,
1148 const int pointgroup_number,
1149 const Holohedry holohedry,
1150 const Centering centering,
1151 const Symmetry *symmetry,
1152 const double symprec)
1153 {
1154 int is_found;
1155 SpacegroupType spacegroup_type;
1156 Symmetry * changed_symmetry;
1157 double changed_lattice[3][3], inv_lattice[3][3], tmat[3][3];
1158
1159 changed_symmetry = NULL;
1160
1161 spacegroup_type = spgdb_get_spacegroup_type(hall_number);
1162
1163 if (pointgroup_number != spacegroup_type.pointgroup_number) {
1164 goto err;
1165 }
1166
1167 switch (holohedry) {
1168 case MONOCLI:
1169 if (match_hall_symbol_db_monocli(origin_shift,
1170 conv_lattice,
1171 orig_lattice,
1172 hall_number,
1173 spacegroup_type.number,
1174 centering,
1175 symmetry,
1176 symprec)) {return 1;}
1177 break;
1178
1179 case ORTHO:
1180 if (num_axis_choices_ortho[spacegroup_type.number - 16] == 6) {
1181 if (match_hall_symbol_db_ortho(origin_shift,
1182 conv_lattice,
1183 orig_lattice,
1184 hall_number,
1185 centering,
1186 symmetry,
1187 6,
1188 symprec)) {return 1;}
1189 break;
1190 }
1191
1192 if (num_axis_choices_ortho[spacegroup_type.number - 16] == 3) {
1193 if (match_hall_symbol_db_ortho(origin_shift,
1194 conv_lattice,
1195 orig_lattice,
1196 hall_number,
1197 centering,
1198 symmetry,
1199 3,
1200 symprec)) {return 1;}
1201 break;
1202 }
1203
1204 /* Switching two axes */
1205 /* Two steps: */
1206 /* 1. Finding principal axis for the representative hall symbol */
1207 /* of the specified hall symbol without checking basis vector */
1208 /* lengths preference. */
1209 /* 2. Finding transformation matrix and origin shift for the */
1210 /* specified hall symbol. */
1211 if (num_axis_choices_ortho[spacegroup_type.number - 16] == 2) {
1212 mat_copy_matrix_d3(changed_lattice, conv_lattice);
1213 if (! match_hall_symbol_db_ortho(
1214 origin_shift,
1215 changed_lattice,
1216 orig_lattice,
1217 spacegroup_to_hall_number[spacegroup_type.number - 1],
1218 centering,
1219 symmetry,
1220 0,
1221 symprec)) {break;}
1222 mat_inverse_matrix_d3(inv_lattice, conv_lattice, 0);
1223 mat_multiply_matrix_d3(tmat, inv_lattice, changed_lattice);
1224
1225 if ((changed_symmetry = get_conventional_symmetry(tmat,
1226 PRIMITIVE,
1227 symmetry)) == NULL) {
1228 goto err;
1229 }
1230
1231 is_found = match_hall_symbol_db_ortho(origin_shift,
1232 changed_lattice,
1233 orig_lattice,
1234 hall_number,
1235 centering,
1236 changed_symmetry,
1237 2,
1238 symprec);
1239 sym_free_symmetry(changed_symmetry);
1240 changed_symmetry = NULL;
1241 if (is_found) {
1242 mat_copy_matrix_d3(conv_lattice, changed_lattice);
1243 return 1;
1244 }
1245 break;
1246 }
1247
1248 if (num_axis_choices_ortho[spacegroup_type.number - 16] == 1) {
1249 if (match_hall_symbol_db_ortho(origin_shift,
1250 conv_lattice,
1251 orig_lattice,
1252 hall_number,
1253 centering,
1254 symmetry,
1255 1,
1256 symprec)) {return 1;}
1257 break;
1258 }
1259
1260 break;
1261
1262 case CUBIC:
1263 if (match_hall_symbol_db_cubic(origin_shift,
1264 conv_lattice,
1265 orig_lattice,
1266 hall_number,
1267 centering,
1268 symmetry,
1269 symprec)) {return 1;}
1270 break;
1271
1272 case TRIGO:
1273 if ((centering == R_CENTER) && (
1274 hall_number == 433 ||
1275 hall_number == 434 ||
1276 hall_number == 436 ||
1277 hall_number == 437 ||
1278 hall_number == 444 ||
1279 hall_number == 445 ||
1280 hall_number == 450 ||
1281 hall_number == 451 ||
1282 hall_number == 452 ||
1283 hall_number == 453 ||
1284 hall_number == 458 ||
1285 hall_number == 459 ||
1286 hall_number == 460 ||
1287 hall_number == 461)) {
1288 /* Rhombohedral. symmetry is for a=b=c basis. */
1289 if (match_hall_symbol_db_rhombo(origin_shift,
1290 conv_lattice,
1291 orig_lattice,
1292 hall_number,
1293 symmetry,
1294 symprec)) {return 1;}
1295 break;
1296 }
1297 /* Do not break for other trigonal cases */
1298 default: /* HEXA, TETRA, TRICLI and rest of TRIGO */
1299 if (match_hall_symbol_db_others(origin_shift,
1300 conv_lattice,
1301 orig_lattice,
1302 hall_number,
1303 centering,
1304 holohedry,
1305 symmetry,
1306 symprec)) {return 1;}
1307 break;
1308 }
1309
1310 err:
1311 return 0;
1312 }
1313
1314 /* Return 0 if failed */
match_hall_symbol_db_monocli(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int hall_number,const int space_group_number,const Centering centering,const Symmetry * conv_symmetry,const double symprec)1315 static int match_hall_symbol_db_monocli(double origin_shift[3],
1316 double conv_lattice[3][3],
1317 SPGCONST double (*orig_lattice)[3],
1318 const int hall_number,
1319 const int space_group_number,
1320 const Centering centering,
1321 const Symmetry *conv_symmetry,
1322 const double symprec)
1323 {
1324 int i, check_norms, i_shortest, is_found_any;
1325 int is_found[36];
1326 double shortest_norm_sum, norm_sum;
1327 double norms_squared[36][2];
1328 double all_conv_lattices[36][3][3], all_origin_shifts[36][3];
1329
1330 /* E. Parthe and L. M. Gelato */
1331 /* "The best unit cell for monoclinic structure..." (1983) */
1332 if (space_group_number == 3 ||
1333 space_group_number == 4 ||
1334 space_group_number == 6 ||
1335 space_group_number == 10 ||
1336 space_group_number == 11) {
1337 /* |a| < |c| for unique axis b. (This is as written in the paper 1983.) */
1338 /* |a| < |b| for unique axis c. (This is spgilb definition.) */
1339 /* |b| < |c| for unique axis a. (This is spgilb definition.) */
1340 check_norms = 1;
1341 } else {
1342 check_norms = 0;
1343 }
1344
1345 for (i = 0; i < 36; i++) {
1346 mat_copy_matrix_d3(all_conv_lattices[i], conv_lattice);
1347 /* all_origin_shifts[i] is possibly overwritten. */
1348 /* is_found == 0: Not found */
1349 /* is_found == 1: Found. conv_lattice may not be similar to orig_lattice */
1350 /* is_found == 2: Found. conv_lattice is similar to orig_lattice */
1351 is_found[i] = match_hall_symbol_db_monocli_in_loop(all_origin_shifts[i],
1352 all_conv_lattices[i],
1353 norms_squared[i],
1354 i,
1355 orig_lattice,
1356 check_norms,
1357 hall_number,
1358 centering,
1359 conv_symmetry,
1360 symprec);
1361 }
1362
1363 is_found_any = 0;
1364 for (i = 0; i < 36; i++) {
1365 if (is_found[i]) {
1366 i_shortest = i;
1367 shortest_norm_sum = sqrt(norms_squared[i][0]) + sqrt(norms_squared[i][1]);
1368 is_found_any = 1;
1369 break;
1370 }
1371 }
1372
1373 if (! is_found_any) {
1374 goto err;
1375 }
1376
1377 /* Find shortest non-unique axes lengths */
1378 for (i = 0; i < 36; i++) {
1379 if (is_found[i]) {
1380 norm_sum = sqrt(norms_squared[i][0]) + sqrt(norms_squared[i][1]);
1381 if (shortest_norm_sum > norm_sum) {
1382 shortest_norm_sum = norm_sum;
1383 }
1384 }
1385 }
1386
1387 /* Prefers is_found == 2, i.e., similar to orig_lattice */
1388 i_shortest = -1;
1389 for (i = 0; i < 36; i++) {
1390 if (is_found[i]) {
1391 norm_sum = sqrt(norms_squared[i][0]) + sqrt(norms_squared[i][1]);
1392 if (mat_Dabs(norm_sum - shortest_norm_sum) < symprec) {
1393 if (is_found[i] == 2) {
1394 i_shortest = i;
1395 break;
1396 }
1397 if (i_shortest < 0) {
1398 i_shortest = i;
1399 }
1400 }
1401 }
1402 }
1403
1404 mat_copy_vector_d3(origin_shift, all_origin_shifts[i_shortest]);
1405 mat_copy_matrix_d3(conv_lattice, all_conv_lattices[i_shortest]);
1406 return 1;
1407
1408 err:
1409 return 0;
1410 }
1411
match_hall_symbol_db_monocli_in_loop(double origin_shift[3],double conv_lattice[3][3],double norms_squared[2],const int i,SPGCONST double (* orig_lattice)[3],const int check_norms,const int hall_number,const Centering centering,const Symmetry * conv_symmetry,const double symprec)1412 static int match_hall_symbol_db_monocli_in_loop(double origin_shift[3],
1413 double conv_lattice[3][3],
1414 double norms_squared[2],
1415 const int i,
1416 SPGCONST double (*orig_lattice)[3],
1417 const int check_norms,
1418 const int hall_number,
1419 const Centering centering,
1420 const Symmetry *conv_symmetry,
1421 const double symprec)
1422 {
1423 int j, k, l, is_found, retval, unique_axis;
1424 double vec[2][3];
1425 double tmat[3][3], change_of_basis[3][3];
1426 Centering changed_centering;
1427 Symmetry * changed_symmetry;
1428
1429 /* centring type should be P or C */
1430 if (centering == C_FACE) {
1431 changed_centering = change_of_centering_monocli[i];
1432 } else { /* suppose PRIMITIVE */
1433 changed_centering = centering;
1434 }
1435
1436 mat_copy_matrix_d3(change_of_basis, change_of_basis_monocli[i]);
1437 mat_multiply_matrix_d3(conv_lattice, conv_lattice, change_of_basis);
1438 unique_axis = change_of_unique_axis_monocli[i];
1439
1440 /* Make non-acute and length preference */
1441 /* norms_squared[0] and norms_squared[1] are the norms of the two */
1442 /* non unique axes among a,b,c. */
1443 l = 0;
1444 for (j = 0; j < 3; j++) {
1445 if (j == unique_axis) {continue;}
1446 for (k = 0; k < 3; k++) {
1447 vec[l][k] = conv_lattice[k][j];
1448 }
1449 norms_squared[l] = mat_norm_squared_d3(vec[l]);
1450 l++;
1451 }
1452
1453 /* discard if principal angle is acute. */
1454 if (vec[0][0] * vec[1][0] +
1455 vec[0][1] * vec[1][1] +
1456 vec[0][2] * vec[1][2] > ZERO_PREC) {
1457 goto cont;
1458 }
1459
1460 /* Choose |a| < |b| < |c| among two non-principles axes */
1461 /* if there are freedom. */
1462 if (check_norms) {
1463 if (norms_squared[0] > norms_squared[1] + ZERO_PREC) {
1464 goto cont;
1465 }
1466 }
1467
1468 /* When orig_lattice is given not NULL, try to find similar (a, b, c) */
1469 /* choises to the input (a, b, c) by flipping a, b, c axes. */
1470 /* Here flipping means a -> -a, and so on. */
1471 /* Note that flipped (a,b,c) that match those of input should not */
1472 /* change centring for monoclinic case. */
1473 retval = 1;
1474 if (orig_lattice != NULL) {
1475 if (mat_get_determinant_d3(orig_lattice) > symprec) {
1476 /* This (mode=1) effectively checks C2 rotation along unique axis. */
1477 if (is_equivalent_lattice(
1478 tmat, 1, conv_lattice, orig_lattice, symprec)) {
1479 if (tmat[(unique_axis + 1) % 3][(unique_axis + 1) % 3]
1480 * tmat[(unique_axis + 2) % 3][(unique_axis + 2) % 3] > ZERO_PREC) {
1481 mat_multiply_matrix_d3(conv_lattice, conv_lattice, tmat);
1482 mat_multiply_matrix_d3(change_of_basis, change_of_basis, tmat);
1483 retval = 2;
1484 }
1485 }
1486 }
1487 }
1488
1489 if ((changed_symmetry = get_conventional_symmetry(change_of_basis,
1490 PRIMITIVE,
1491 conv_symmetry)) == NULL) {
1492 goto cont;
1493 }
1494
1495 is_found = hal_match_hall_symbol_db(origin_shift,
1496 conv_lattice,
1497 hall_number,
1498 changed_centering,
1499 changed_symmetry,
1500 symprec);
1501 sym_free_symmetry(changed_symmetry);
1502 changed_symmetry = NULL;
1503
1504 if (is_found) {
1505 return retval;
1506 }
1507
1508 cont:
1509 return 0;
1510 }
1511
1512 /* Return 0 if failed */
match_hall_symbol_db_ortho(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int hall_number,const Centering centering,const Symmetry * conv_symmetry,const int num_free_axes,const double symprec)1513 static int match_hall_symbol_db_ortho(double origin_shift[3],
1514 double conv_lattice[3][3],
1515 SPGCONST double (*orig_lattice)[3],
1516 const int hall_number,
1517 const Centering centering,
1518 const Symmetry *conv_symmetry,
1519 const int num_free_axes,
1520 const double symprec)
1521 {
1522 int i;
1523
1524 /* Try to find the best (a, b, c) by flipping axes. */
1525 if (orig_lattice != NULL) {
1526 if (mat_get_determinant_d3(orig_lattice) > symprec) {
1527 for (i = 0; i < 6; i++) {
1528 if (match_hall_symbol_db_ortho_in_loop(origin_shift,
1529 conv_lattice,
1530 orig_lattice,
1531 i,
1532 hall_number,
1533 centering,
1534 conv_symmetry,
1535 num_free_axes,
1536 symprec)) {
1537 return 1;
1538 }
1539 }
1540 }
1541 }
1542
1543 /* If flipping didn't work, usual search is made. */
1544 for (i = 0; i < 6; i++) {
1545 if (match_hall_symbol_db_ortho_in_loop(origin_shift,
1546 conv_lattice,
1547 NULL,
1548 i,
1549 hall_number,
1550 centering,
1551 conv_symmetry,
1552 num_free_axes,
1553 symprec)) {
1554 return 1;
1555 }
1556 }
1557
1558 return 0;
1559 }
1560
match_hall_symbol_db_ortho_in_loop(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int i,const int hall_number,const Centering centering,const Symmetry * symmetry,const int num_free_axes,const double symprec)1561 static int match_hall_symbol_db_ortho_in_loop(double origin_shift[3],
1562 double conv_lattice[3][3],
1563 SPGCONST double (*orig_lattice)[3],
1564 const int i,
1565 const int hall_number,
1566 const Centering centering,
1567 const Symmetry *symmetry,
1568 const int num_free_axes,
1569 const double symprec)
1570 {
1571 int j, k, l, is_found;
1572 double vec[3], norms[3];
1573 Centering changed_centering;
1574 Symmetry * changed_symmetry;
1575 double changed_lattice[3][3], tmat[3][3], change_of_basis[3][3];
1576
1577 changed_symmetry = NULL;
1578
1579 if (centering == C_FACE) {
1580 changed_centering = change_of_centering_ortho[i];
1581 } else {
1582 changed_centering = centering;
1583 }
1584
1585 mat_multiply_matrix_d3(changed_lattice,
1586 conv_lattice,
1587 change_of_basis_ortho[i]);
1588
1589 /* When orig_lattice is given not NULL, try to find similar (a, b, c) */
1590 /* choises to the input (a, b, c) by flipping a, b, c axes. */
1591 /* Here flipping means a -> -a, and so on. */
1592 /* Note that flipping of axes doesn't change centring. */
1593 mat_copy_matrix_d3(change_of_basis, change_of_basis_ortho[i]);
1594 if (orig_lattice != NULL) {
1595 if (is_equivalent_lattice(
1596 tmat, 1, changed_lattice, orig_lattice, symprec)) {
1597 mat_multiply_matrix_d3(changed_lattice, changed_lattice, tmat);
1598 mat_multiply_matrix_d3(change_of_basis, change_of_basis, tmat);
1599 } else {
1600 goto cont; /* This is necessary to run through all */
1601 /* change_of_basis_ortho. */
1602 }
1603 }
1604
1605 if (num_free_axes == 2) {
1606 l = 0;
1607 for (j = 0; j < 3; j++) {
1608 if (j == change_of_unique_axis_ortho[i]) {continue;}
1609 for (k = 0; k < 3; k++) {vec[k] = changed_lattice[k][j];}
1610 norms[l] = mat_norm_squared_d3(vec);
1611 l++;
1612 }
1613 if (norms[0] > norms[1] + ZERO_PREC) {goto cont;}
1614 }
1615
1616 if (num_free_axes == 3) {
1617 for (j = 0; j < 3; j++) {
1618 for (k = 0; k < 3; k++) {vec[k] = changed_lattice[k][j];}
1619 norms[j] = mat_norm_squared_d3(vec);
1620 }
1621 if ((norms[0] > norms[1] + ZERO_PREC) ||
1622 (norms[0] > norms[2] + ZERO_PREC)) {goto cont;}
1623 }
1624
1625 if (num_free_axes == 6) {
1626 for (j = 0; j < 3; j++) {
1627 for (k = 0; k < 3; k++) {vec[k] = changed_lattice[k][j];}
1628 norms[j] = mat_norm_squared_d3(vec);
1629 }
1630 if ((norms[0] > norms[1] + ZERO_PREC) ||
1631 (norms[1] > norms[2] + ZERO_PREC)) {goto cont;}
1632 }
1633
1634 if ((changed_symmetry = get_conventional_symmetry(change_of_basis,
1635 PRIMITIVE,
1636 symmetry)) == NULL) {
1637 goto cont;
1638 }
1639
1640 is_found = hal_match_hall_symbol_db(origin_shift,
1641 changed_lattice,
1642 hall_number,
1643 changed_centering,
1644 changed_symmetry,
1645 symprec);
1646 sym_free_symmetry(changed_symmetry);
1647 changed_symmetry = NULL;
1648 if (is_found) {
1649 mat_copy_matrix_d3(conv_lattice, changed_lattice);
1650 return 1;
1651 }
1652
1653 cont:
1654 return 0;
1655 }
1656
match_hall_symbol_db_cubic(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int hall_number,const Centering centering,const Symmetry * conv_symmetry,const double symprec)1657 static int match_hall_symbol_db_cubic(double origin_shift[3],
1658 double conv_lattice[3][3],
1659 SPGCONST double (*orig_lattice)[3],
1660 const int hall_number,
1661 const Centering centering,
1662 const Symmetry *conv_symmetry,
1663 const double symprec)
1664 {
1665 int i;
1666
1667 /* Special treatment for No. 205 (501) is included in this change of */
1668 /* basis. To see old code, commit hash 67997654 and change_of_basis_501. */
1669 if (orig_lattice != NULL) {
1670 if (mat_get_determinant_d3(orig_lattice) > symprec) {
1671 for (i = 0; i < 6; i++) {
1672 if (match_hall_symbol_db_cubic_in_loop(origin_shift,
1673 conv_lattice,
1674 orig_lattice,
1675 i,
1676 hall_number,
1677 centering,
1678 conv_symmetry,
1679 symprec)) {
1680 return 1;
1681 }
1682 }
1683 }
1684 }
1685
1686 for (i = 0; i < 6; i++) {
1687 if (match_hall_symbol_db_cubic_in_loop(origin_shift,
1688 conv_lattice,
1689 NULL,
1690 i,
1691 hall_number,
1692 centering,
1693 conv_symmetry,
1694 symprec)) {
1695 return 1;
1696 }
1697 }
1698 return 0;
1699 }
1700
match_hall_symbol_db_cubic_in_loop(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int i,const int hall_number,const Centering centering,const Symmetry * conv_symmetry,const double symprec)1701 static int match_hall_symbol_db_cubic_in_loop(double origin_shift[3],
1702 double conv_lattice[3][3],
1703 SPGCONST double (*orig_lattice)[3],
1704 const int i,
1705 const int hall_number,
1706 const Centering centering,
1707 const Symmetry *conv_symmetry,
1708 const double symprec)
1709 {
1710 int is_found;
1711 double changed_lattice[3][3], tmat[3][3], change_of_basis[3][3];
1712 Symmetry *changed_symmetry;
1713
1714 changed_symmetry = NULL;
1715
1716 mat_copy_matrix_d3(change_of_basis, change_of_basis_ortho[i]);
1717 mat_multiply_matrix_d3(changed_lattice, conv_lattice, change_of_basis);
1718
1719 if (orig_lattice != NULL) {
1720 if (is_equivalent_lattice(
1721 tmat, 1, changed_lattice, orig_lattice, symprec)) {
1722 mat_multiply_matrix_d3(changed_lattice, changed_lattice, tmat);
1723 mat_multiply_matrix_d3(change_of_basis, change_of_basis, tmat);
1724 } else {
1725 goto cont; /* This is necessary to run through all */
1726 /* change_of_basis_ortho. */
1727 }
1728 }
1729
1730 if ((changed_symmetry = get_conventional_symmetry(change_of_basis,
1731 PRIMITIVE,
1732 conv_symmetry)) == NULL) {
1733 goto cont;
1734 }
1735
1736 is_found = hal_match_hall_symbol_db(origin_shift,
1737 changed_lattice,
1738 hall_number,
1739 centering,
1740 changed_symmetry,
1741 symprec);
1742
1743 sym_free_symmetry(changed_symmetry);
1744 changed_symmetry = NULL;
1745
1746 if (is_found) {
1747 mat_copy_matrix_d3(conv_lattice, changed_lattice);
1748 return 1;
1749 }
1750
1751 cont:
1752 return 0;
1753 }
1754
match_hall_symbol_db_rhombo(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int hall_number,const Symmetry * conv_symmetry,const double symprec)1755 static int match_hall_symbol_db_rhombo(double origin_shift[3],
1756 double conv_lattice[3][3],
1757 SPGCONST double (*orig_lattice)[3],
1758 const int hall_number,
1759 const Symmetry *conv_symmetry,
1760 const double symprec)
1761 {
1762 int is_found;
1763
1764 if (hall_number == 433 ||
1765 hall_number == 436 ||
1766 hall_number == 444 ||
1767 hall_number == 450 ||
1768 hall_number == 452 ||
1769 hall_number == 458 ||
1770 hall_number == 460) { /* Hexagonal lattice */
1771 is_found = match_hall_symbol_db_change_of_basis_loop(
1772 origin_shift,
1773 conv_lattice,
1774 orig_lattice,
1775 change_of_basis_rhombo_hex,
1776 6,
1777 hall_number,
1778 R_CENTER,
1779 conv_symmetry,
1780 symprec);
1781 if (is_found) {
1782 /* mat_copy_matrix_d3(conv_lattice, changed_lattice); */
1783 return 1;
1784 } else {
1785 return 0;
1786 }
1787 } else { /* Rhombohedral (a=b=c) lattice */
1788 return match_hall_symbol_db_change_of_basis_loop(
1789 origin_shift,
1790 conv_lattice,
1791 orig_lattice,
1792 change_of_basis_rhombo,
1793 6,
1794 hall_number,
1795 PRIMITIVE,
1796 conv_symmetry,
1797 symprec);
1798 }
1799
1800 return 0;
1801 }
1802
1803 /* HEXA, TETRA, TRICLI and TRIGO but not Rhombohedral a=b=c */
match_hall_symbol_db_others(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],const int hall_number,const Centering centering,const Holohedry holohedry,const Symmetry * conv_symmetry,const double symprec)1804 static int match_hall_symbol_db_others(double origin_shift[3],
1805 double conv_lattice[3][3],
1806 SPGCONST double (*orig_lattice)[3],
1807 const int hall_number,
1808 const Centering centering,
1809 const Holohedry holohedry,
1810 const Symmetry *conv_symmetry,
1811 const double symprec)
1812 {
1813 /* TRICLI: No check. */
1814 if (holohedry == TRICLI) {
1815 return hal_match_hall_symbol_db(origin_shift,
1816 conv_lattice,
1817 hall_number,
1818 centering,
1819 conv_symmetry,
1820 symprec);
1821 }
1822
1823 /* TETRA: Check over 4 fold rotations with/without flipping c-axis */
1824 if (holohedry == TETRA) {
1825 return match_hall_symbol_db_change_of_basis_loop(
1826 origin_shift,
1827 conv_lattice,
1828 orig_lattice,
1829 change_of_basis_C4,
1830 8,
1831 hall_number,
1832 centering,
1833 conv_symmetry,
1834 symprec);
1835 }
1836
1837 /* HEXA and TRIGO but not Rhombohedral: */
1838 /* Check over 6 fold rotations with/without flipping c-axis */
1839 return match_hall_symbol_db_change_of_basis_loop(
1840 origin_shift,
1841 conv_lattice,
1842 orig_lattice,
1843 change_of_basis_C6,
1844 12,
1845 hall_number,
1846 centering,
1847 conv_symmetry,
1848 symprec);
1849 }
1850
match_hall_symbol_db_change_of_basis_loop(double origin_shift[3],double conv_lattice[3][3],SPGCONST double (* orig_lattice)[3],SPGCONST double (* change_of_basis)[3][3],const int num_change_of_basis,const int hall_number,const Centering centering,const Symmetry * conv_symmetry,const double symprec)1851 static int match_hall_symbol_db_change_of_basis_loop(
1852 double origin_shift[3],
1853 double conv_lattice[3][3],
1854 SPGCONST double (*orig_lattice)[3],
1855 SPGCONST double (*change_of_basis)[3][3],
1856 const int num_change_of_basis,
1857 const int hall_number,
1858 const Centering centering,
1859 const Symmetry *conv_symmetry,
1860 const double symprec)
1861 {
1862 int i, is_found;
1863 double changed_lattice[3][3], tmat[3][3];
1864 Symmetry *changed_symmetry;
1865 Centering centering_for_symmetry;
1866
1867 changed_symmetry = NULL;
1868
1869 if (centering == R_CENTER) {
1870 centering_for_symmetry = R_CENTER;
1871 } else {
1872 centering_for_symmetry = PRIMITIVE;
1873 }
1874
1875 /* Check of similarity of basis vectors to those of input */
1876 if (orig_lattice != NULL) {
1877 if (mat_get_determinant_d3(orig_lattice) > symprec) {
1878 for (i = 0; i < num_change_of_basis; i++) {
1879 if ((changed_symmetry = get_conventional_symmetry(
1880 change_of_basis[i], centering_for_symmetry, conv_symmetry))
1881 == NULL) {
1882 continue;
1883 }
1884 mat_multiply_matrix_d3(
1885 changed_lattice, conv_lattice, change_of_basis[i]);
1886 if (is_equivalent_lattice(
1887 tmat, 0, changed_lattice, orig_lattice, symprec)) {
1888 /* Here R_CENTER means hP (hexagonal) setting. */
1889 is_found = hal_match_hall_symbol_db(origin_shift,
1890 changed_lattice,
1891 hall_number,
1892 centering,
1893 changed_symmetry,
1894 symprec);
1895 } else {
1896 is_found = 0;
1897 }
1898 sym_free_symmetry(changed_symmetry);
1899 changed_symmetry = NULL;
1900 if (is_found) {
1901 mat_copy_matrix_d3(conv_lattice, changed_lattice);
1902 return 1;
1903 }
1904 }
1905 }
1906 }
1907
1908 /* No check of similarity of basis vectors to those of input */
1909 for (i = 0; i < num_change_of_basis; i++) {
1910 if ((changed_symmetry = get_conventional_symmetry(
1911 change_of_basis[i], centering_for_symmetry, conv_symmetry)) == NULL) {
1912 continue;
1913 }
1914 mat_multiply_matrix_d3(changed_lattice, conv_lattice, change_of_basis[i]);
1915 is_found = hal_match_hall_symbol_db(origin_shift,
1916 changed_lattice,
1917 hall_number,
1918 centering,
1919 changed_symmetry,
1920 symprec);
1921 sym_free_symmetry(changed_symmetry);
1922 changed_symmetry = NULL;
1923 if (is_found) {
1924 mat_copy_matrix_d3(conv_lattice, changed_lattice);
1925 return 1;
1926 }
1927 }
1928
1929 return 0;
1930 }
1931
1932 /* Return NULL if failed */
get_conventional_symmetry(SPGCONST double tmat[3][3],const Centering centering,const Symmetry * primitive_sym)1933 static Symmetry * get_conventional_symmetry(SPGCONST double tmat[3][3],
1934 const Centering centering,
1935 const Symmetry *primitive_sym)
1936 {
1937 int i, j, k, multi, size;
1938 double inv_tmat[3][3], shift[3][3];
1939 double symmetry_rot_d3[3][3], primitive_sym_rot_d3[3][3];
1940 Symmetry *symmetry;
1941
1942 symmetry = NULL;
1943
1944 size = primitive_sym->size;
1945
1946 switch (centering) {
1947 case FACE:
1948 if ((symmetry = sym_alloc_symmetry(size * 4)) == NULL) {
1949 return NULL;
1950 }
1951 break;
1952 case R_CENTER:
1953 if ((symmetry = sym_alloc_symmetry(size * 3)) == NULL) {
1954 return NULL;
1955 }
1956 break;
1957 case BODY:
1958 case A_FACE:
1959 case B_FACE:
1960 case C_FACE:
1961 if ((symmetry = sym_alloc_symmetry(size * 2)) == NULL) {
1962 return NULL;
1963 }
1964 break;
1965 default:
1966 if ((symmetry = sym_alloc_symmetry(size)) == NULL) {
1967 return NULL;
1968 }
1969 break;
1970 }
1971
1972 for (i = 0; i < size; i++) {
1973 mat_cast_matrix_3i_to_3d(primitive_sym_rot_d3, primitive_sym->rot[i]);
1974
1975 /* C*S*C^-1: recover conventional cell symmetry operation */
1976 mat_get_similar_matrix_d3(symmetry_rot_d3,
1977 primitive_sym_rot_d3,
1978 tmat,
1979 0);
1980 mat_cast_matrix_3d_to_3i(symmetry->rot[i], symmetry_rot_d3);
1981
1982 /* translation in conventional cell: C = B^-1*P */
1983 mat_inverse_matrix_d3(inv_tmat, tmat, 0);
1984 mat_multiply_matrix_vector_d3(symmetry->trans[i],
1985 inv_tmat,
1986 primitive_sym->trans[i]);
1987 }
1988
1989 multi = 1;
1990
1991 if (centering != PRIMITIVE) {
1992 multi = get_centering_shifts(shift, centering);
1993 for (i = 0; i < multi - 1; i++) {
1994 for (j = 0; j < size; j++) {
1995 mat_copy_matrix_i3(symmetry->rot[(i+1) * size + j],
1996 symmetry->rot[j]);
1997 for (k = 0; k < 3; k++) {
1998 symmetry->trans[(i+1) * size + j][k] =
1999 symmetry->trans[j][k] + shift[i][k];
2000 }
2001 }
2002 }
2003 }
2004
2005 for (i = 0; i < multi; i++) {
2006 for (j = 0; j < size; j++) {
2007 for (k = 0; k < 3; k++) {
2008 symmetry->trans[i * size + j][k] =
2009 mat_Dmod1(symmetry->trans[i * size + j][k]);
2010 }
2011 }
2012 }
2013
2014 return symmetry;
2015 }
2016
2017 /* Return CENTERING_ERROR if failed */
get_centering(double correction_mat[3][3],SPGCONST int tmat[3][3],const Laue laue)2018 static Centering get_centering(double correction_mat[3][3],
2019 SPGCONST int tmat[3][3],
2020 const Laue laue)
2021 {
2022 int det;
2023 double trans_corr_mat[3][3];
2024 Centering centering;
2025
2026 mat_copy_matrix_d3(correction_mat, identity);
2027 det = abs(mat_get_determinant_i3(tmat));
2028 debug_print("laue class: %d\n", laue);
2029 debug_print("multiplicity: %d\n", det);
2030
2031 switch (det) {
2032
2033 case 1:
2034 centering = PRIMITIVE;
2035 break;
2036
2037 case 2:
2038 centering = get_base_center(tmat);
2039 if (centering == A_FACE) {
2040 if (laue == LAUE2M) {
2041 debug_print("Monocli A to C\n");
2042 mat_copy_matrix_d3(correction_mat, monocli_a2c);
2043 } else {
2044 mat_copy_matrix_d3(correction_mat, a2c);
2045 }
2046 centering = C_FACE;
2047 }
2048 if (centering == B_FACE) {
2049 mat_copy_matrix_d3(correction_mat, b2c);
2050 centering = C_FACE;
2051 }
2052 if (laue == LAUE2M && centering == BODY) {
2053 debug_print("Monocli I to C\n");
2054 mat_copy_matrix_d3(correction_mat, monocli_i2c);
2055 centering = C_FACE;
2056 }
2057 break;
2058
2059 case 3:
2060 /* hP (a=b) but not hR (a=b=c) */
2061 centering = R_CENTER;
2062 mat_multiply_matrix_id3(trans_corr_mat, tmat, rhombo_obverse);
2063 if (mat_is_int_matrix(trans_corr_mat, INT_PREC)) {
2064 mat_copy_matrix_d3(correction_mat, rhombo_obverse);
2065 debug_print("R-center observe setting\n");
2066 debug_print_matrix_d3(trans_corr_mat);
2067 }
2068 mat_multiply_matrix_id3(trans_corr_mat, tmat, rhomb_reverse);
2069 if (mat_is_int_matrix(trans_corr_mat, INT_PREC)) {
2070 mat_copy_matrix_d3(correction_mat, rhomb_reverse);
2071 debug_print("R-center reverse setting\n");
2072 debug_print_matrix_d3(trans_corr_mat);
2073 }
2074 break;
2075
2076 case 4:
2077 centering = FACE;
2078 break;
2079
2080 default:
2081 centering = CENTERING_ERROR;
2082 break;
2083 }
2084
2085 return centering;
2086 }
2087
get_base_center(SPGCONST int tmat[3][3])2088 static Centering get_base_center(SPGCONST int tmat[3][3])
2089 {
2090 int i;
2091 Centering centering = PRIMITIVE;
2092
2093 debug_print("lat_get_base_center\n");
2094
2095 /* C center */
2096 for (i = 0; i < 3; i++) {
2097 if (tmat[i][0] == 0 &&
2098 tmat[i][1] == 0 &&
2099 abs(tmat[i][2]) == 1) {
2100 centering = C_FACE;
2101 goto end;
2102 }
2103 }
2104
2105 /* A center */
2106 for (i = 0; i < 3; i++) {
2107 if (abs(tmat[i][0]) == 1 &&
2108 tmat[i][1] == 0 &&
2109 tmat[i][2] == 0) {
2110 centering = A_FACE;
2111 goto end;
2112 }
2113 }
2114
2115 /* B center */
2116 for (i = 0; i < 3; i++) {
2117 if (tmat[i][0] == 0 &&
2118 abs(tmat[i][1]) == 1 &&
2119 tmat[i][2] == 0) {
2120 centering = B_FACE;
2121 goto end;
2122 }
2123 }
2124
2125 /* body center */
2126 if ((abs(tmat[0][0]) +
2127 abs(tmat[0][1]) +
2128 abs(tmat[0][2]) == 2) &&
2129 (abs(tmat[1][0]) +
2130 abs(tmat[1][1]) +
2131 abs(tmat[1][2]) == 2) &&
2132 (abs(tmat[2][0]) +
2133 abs(tmat[2][1]) +
2134 abs(tmat[2][2]) == 2)) {
2135 centering = BODY;
2136 goto end;
2137 }
2138
2139 /* This should not happen. */
2140 warning_print("spglib: No centring was found (line %d, %s).\n", __LINE__, __FILE__);
2141 return PRIMITIVE;
2142
2143 end:
2144 return centering;
2145 }
2146
get_centering_shifts(double shift[3][3],const Centering centering)2147 static int get_centering_shifts(double shift[3][3],
2148 const Centering centering)
2149 {
2150 int i, j, multi;
2151
2152 multi = 1;
2153 for (i = 0; i < 3; i++) {
2154 for (j = 0; j < 3; j++) {
2155 shift[i][j] = 0;
2156 }
2157 }
2158
2159 if (centering != PRIMITIVE) {
2160 if (centering != FACE && centering != R_CENTER) {
2161 for (i = 0; i < 3; i++) { shift[0][i] = 0.5; } /* BODY */
2162 if (centering == A_FACE) { shift[0][0] = 0; }
2163 if (centering == B_FACE) { shift[0][1] = 0; }
2164 if (centering == C_FACE) { shift[0][2] = 0; }
2165 multi = 2;
2166 }
2167
2168 if (centering == R_CENTER) {
2169 shift[0][0] = 2. / 3;
2170 shift[0][1] = 1. / 3;
2171 shift[0][2] = 1. / 3;
2172 shift[1][0] = 1. / 3;
2173 shift[1][1] = 2. / 3;
2174 shift[1][2] = 2. / 3;
2175 multi = 3;
2176 }
2177
2178 if (centering == FACE) {
2179 shift[0][0] = 0;
2180 shift[0][1] = 0.5;
2181 shift[0][2] = 0.5;
2182 shift[1][0] = 0.5;
2183 shift[1][1] = 0;
2184 shift[1][2] = 0.5;
2185 shift[2][0] = 0.5;
2186 shift[2][1] = 0.5;
2187 shift[2][2] = 0;
2188 multi = 4;
2189 }
2190 }
2191
2192 return multi;
2193 }
2194
2195 /* allow_flip: Flip axes when abs(P)=I but P!=I. */
2196 /* lattice is overwritten. */
is_equivalent_lattice(double tmat[3][3],const int mode,SPGCONST double lattice[3][3],SPGCONST double orig_lattice[3][3],const double symprec)2197 static int is_equivalent_lattice(double tmat[3][3],
2198 const int mode,
2199 SPGCONST double lattice[3][3],
2200 SPGCONST double orig_lattice[3][3],
2201 const double symprec)
2202 {
2203 int i, j;
2204 double inv_lat[3][3], tmat_abs[3][3], metric[3][3], orig_metric[3][3];
2205 int tmat_int[3][3];
2206
2207 if (mat_Dabs(mat_get_determinant_d3(lattice) -
2208 mat_get_determinant_d3(orig_lattice)) > symprec) {
2209 goto fail;
2210 }
2211
2212 if (!mat_inverse_matrix_d3(inv_lat, lattice, symprec)) {
2213 goto fail;
2214 }
2215
2216 mat_multiply_matrix_d3(tmat, inv_lat, orig_lattice);
2217
2218 switch (mode) {
2219 case 0: /* Check identity of all elements */
2220 if (mat_check_identity_matrix_d3(identity, tmat, symprec)) {
2221 return 1;
2222 }
2223 break;
2224
2225 case 1: /* Check identity of all elements allowing axes flips */
2226 for (i = 0; i < 3; i++) {
2227 for (j = 0; j < 3; j++) {
2228 tmat_abs[i][j] = mat_Dabs(tmat[i][j]);
2229 }
2230 }
2231
2232 if (mat_check_identity_matrix_d3(identity, tmat_abs, symprec)) {
2233 return 1;
2234 }
2235 break;
2236
2237 case 2: /* Check metric tensors */
2238 mat_cast_matrix_3d_to_3i(tmat_int, tmat);
2239 if (! mat_check_identity_matrix_id3(tmat_int, tmat, symprec)) {
2240 break;
2241 }
2242 if (mat_get_determinant_i3(tmat_int) != 1) {
2243 break;
2244 }
2245 mat_get_metric(orig_metric, orig_lattice);
2246 mat_get_metric(metric, lattice);
2247 if (mat_check_identity_matrix_d3(orig_metric, metric, symprec)) {
2248 return 1;
2249 }
2250 break;
2251 }
2252
2253 fail:
2254 return 0;
2255 }
2256