1 /* Copyright (C) 2008 Atsushi Togo */
2 /* All rights reserved. */
3
4 /* This file is part of spglib. */
5
6 /* Redistribution and use in source and binary forms, with or without */
7 /* modification, are permitted provided that the following conditions */
8 /* are met: */
9
10 /* * Redistributions of source code must retain the above copyright */
11 /* notice, this list of conditions and the following disclaimer. */
12
13 /* * Redistributions in binary form must reproduce the above copyright */
14 /* notice, this list of conditions and the following disclaimer in */
15 /* the documentation and/or other materials provided with the */
16 /* distribution. */
17
18 /* * Neither the name of the 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 <stddef.h>
38 #include "mathfunc.h"
39 #include "kpoint.h"
40 #include "kgrid.h"
41
42 #ifdef KPTWARNING
43 #include <stdio.h>
44 #define warning_print(...) fprintf(stderr,__VA_ARGS__)
45 #else
46 #define warning_print(...)
47 #endif
48
49 #define KPT_NUM_BZ_SEARCH_SPACE 125
50
51 static int bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
52 { 0, 0, 0},
53 { 0, 0, 1},
54 { 0, 0, 2},
55 { 0, 0, -2},
56 { 0, 0, -1},
57 { 0, 1, 0},
58 { 0, 1, 1},
59 { 0, 1, 2},
60 { 0, 1, -2},
61 { 0, 1, -1},
62 { 0, 2, 0},
63 { 0, 2, 1},
64 { 0, 2, 2},
65 { 0, 2, -2},
66 { 0, 2, -1},
67 { 0, -2, 0},
68 { 0, -2, 1},
69 { 0, -2, 2},
70 { 0, -2, -2},
71 { 0, -2, -1},
72 { 0, -1, 0},
73 { 0, -1, 1},
74 { 0, -1, 2},
75 { 0, -1, -2},
76 { 0, -1, -1},
77 { 1, 0, 0},
78 { 1, 0, 1},
79 { 1, 0, 2},
80 { 1, 0, -2},
81 { 1, 0, -1},
82 { 1, 1, 0},
83 { 1, 1, 1},
84 { 1, 1, 2},
85 { 1, 1, -2},
86 { 1, 1, -1},
87 { 1, 2, 0},
88 { 1, 2, 1},
89 { 1, 2, 2},
90 { 1, 2, -2},
91 { 1, 2, -1},
92 { 1, -2, 0},
93 { 1, -2, 1},
94 { 1, -2, 2},
95 { 1, -2, -2},
96 { 1, -2, -1},
97 { 1, -1, 0},
98 { 1, -1, 1},
99 { 1, -1, 2},
100 { 1, -1, -2},
101 { 1, -1, -1},
102 { 2, 0, 0},
103 { 2, 0, 1},
104 { 2, 0, 2},
105 { 2, 0, -2},
106 { 2, 0, -1},
107 { 2, 1, 0},
108 { 2, 1, 1},
109 { 2, 1, 2},
110 { 2, 1, -2},
111 { 2, 1, -1},
112 { 2, 2, 0},
113 { 2, 2, 1},
114 { 2, 2, 2},
115 { 2, 2, -2},
116 { 2, 2, -1},
117 { 2, -2, 0},
118 { 2, -2, 1},
119 { 2, -2, 2},
120 { 2, -2, -2},
121 { 2, -2, -1},
122 { 2, -1, 0},
123 { 2, -1, 1},
124 { 2, -1, 2},
125 { 2, -1, -2},
126 { 2, -1, -1},
127 {-2, 0, 0},
128 {-2, 0, 1},
129 {-2, 0, 2},
130 {-2, 0, -2},
131 {-2, 0, -1},
132 {-2, 1, 0},
133 {-2, 1, 1},
134 {-2, 1, 2},
135 {-2, 1, -2},
136 {-2, 1, -1},
137 {-2, 2, 0},
138 {-2, 2, 1},
139 {-2, 2, 2},
140 {-2, 2, -2},
141 {-2, 2, -1},
142 {-2, -2, 0},
143 {-2, -2, 1},
144 {-2, -2, 2},
145 {-2, -2, -2},
146 {-2, -2, -1},
147 {-2, -1, 0},
148 {-2, -1, 1},
149 {-2, -1, 2},
150 {-2, -1, -2},
151 {-2, -1, -1},
152 {-1, 0, 0},
153 {-1, 0, 1},
154 {-1, 0, 2},
155 {-1, 0, -2},
156 {-1, 0, -1},
157 {-1, 1, 0},
158 {-1, 1, 1},
159 {-1, 1, 2},
160 {-1, 1, -2},
161 {-1, 1, -1},
162 {-1, 2, 0},
163 {-1, 2, 1},
164 {-1, 2, 2},
165 {-1, 2, -2},
166 {-1, 2, -1},
167 {-1, -2, 0},
168 {-1, -2, 1},
169 {-1, -2, 2},
170 {-1, -2, -2},
171 {-1, -2, -1},
172 {-1, -1, 0},
173 {-1, -1, 1},
174 {-1, -1, 2},
175 {-1, -1, -2},
176 {-1, -1, -1}
177 };
178
179 static MatINT *get_point_group_reciprocal(const MatINT * rotations,
180 const int is_time_reversal);
181 static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
182 const double symprec,
183 const size_t num_q,
184 SPGCONST double qpoints[][3]);
185 static size_t get_dense_ir_reciprocal_mesh(int grid_address[][3],
186 size_t ir_mapping_table[],
187 const int mesh[3],
188 const int is_shift[3],
189 const MatINT *rot_reciprocal);
190 static size_t get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
191 size_t ir_mapping_table[],
192 const int mesh[3],
193 const int is_shift[3],
194 const MatINT *rot_reciprocal);
195 static size_t get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
196 size_t ir_mapping_table[],
197 const int mesh[3],
198 const int is_shift[3],
199 const MatINT *rot_reciprocal);
200 static size_t get_dense_num_ir(size_t ir_mapping_table[], const int mesh[3]);
201 static size_t relocate_dense_BZ_grid_address(int bz_grid_address[][3],
202 size_t bz_map[],
203 SPGCONST int grid_address[][3],
204 const int mesh[3],
205 SPGCONST double rec_lattice[3][3],
206 const int is_shift[3]);
207 static double get_tolerance_for_BZ_reduction(SPGCONST double rec_lattice[3][3],
208 const int mesh[3]);
209 static int check_mesh_symmetry(const int mesh[3],
210 const int is_shift[3],
211 const MatINT *rot_reciprocal);
212
213 /* grid_address (e.g. 4x4x4 mesh, unless GRID_ORDER_XYZ is defined) */
214 /* [[ 0 0 0] */
215 /* [ 1 0 0] */
216 /* [ 2 0 0] */
217 /* [-1 0 0] */
218 /* [ 0 1 0] */
219 /* [ 1 1 0] */
220 /* [ 2 1 0] */
221 /* [-1 1 0] */
222 /* .... ] */
223 /* */
224 /* Each value of 'map' correspnds to the index of grid_point. */
kpt_get_irreducible_reciprocal_mesh(int grid_address[][3],int ir_mapping_table[],const int mesh[3],const int is_shift[3],const MatINT * rot_reciprocal)225 int kpt_get_irreducible_reciprocal_mesh(int grid_address[][3],
226 int ir_mapping_table[],
227 const int mesh[3],
228 const int is_shift[3],
229 const MatINT *rot_reciprocal)
230 {
231 int num_ir;
232 size_t i;
233 size_t *dense_ir_mapping_table;
234
235 if ((dense_ir_mapping_table =
236 (size_t*)malloc(sizeof(size_t) * mesh[0] * mesh[1] * mesh[2])) == NULL) {
237 warning_print("spglib: Memory of unique_rot could not be allocated.");
238 return 0;
239 }
240
241 num_ir = kpt_get_dense_irreducible_reciprocal_mesh(grid_address,
242 dense_ir_mapping_table,
243 mesh,
244 is_shift,
245 rot_reciprocal);
246
247 for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
248 ir_mapping_table[i] = dense_ir_mapping_table[i];
249 }
250
251 free(dense_ir_mapping_table);
252 dense_ir_mapping_table = NULL;
253
254 return num_ir;
255 }
256
kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],size_t ir_mapping_table[],const int mesh[3],const int is_shift[3],const MatINT * rot_reciprocal)257 size_t kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
258 size_t ir_mapping_table[],
259 const int mesh[3],
260 const int is_shift[3],
261 const MatINT *rot_reciprocal)
262 {
263 size_t num_ir;
264
265 num_ir = get_dense_ir_reciprocal_mesh(grid_address,
266 ir_mapping_table,
267 mesh,
268 is_shift,
269 rot_reciprocal);
270
271 return num_ir;
272 }
273
kpt_get_stabilized_reciprocal_mesh(int grid_address[][3],int ir_mapping_table[],const int mesh[3],const int is_shift[3],const int is_time_reversal,const MatINT * rotations,const size_t num_q,SPGCONST double qpoints[][3])274 int kpt_get_stabilized_reciprocal_mesh(int grid_address[][3],
275 int ir_mapping_table[],
276 const int mesh[3],
277 const int is_shift[3],
278 const int is_time_reversal,
279 const MatINT * rotations,
280 const size_t num_q,
281 SPGCONST double qpoints[][3])
282 {
283 int num_ir;
284 size_t i;
285 size_t *dense_ir_mapping_table;
286
287 if ((dense_ir_mapping_table =
288 (size_t*)malloc(sizeof(size_t) * mesh[0] * mesh[1] * mesh[2])) == NULL) {
289 warning_print("spglib: Memory of unique_rot could not be allocated.");
290 return 0;
291 }
292
293 num_ir = kpt_get_dense_stabilized_reciprocal_mesh(grid_address,
294 dense_ir_mapping_table,
295 mesh,
296 is_shift,
297 is_time_reversal,
298 rotations,
299 num_q,
300 qpoints);
301
302 for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
303 ir_mapping_table[i] = dense_ir_mapping_table[i];
304 }
305
306 free(dense_ir_mapping_table);
307 dense_ir_mapping_table = NULL;
308
309 return num_ir;
310 }
311
kpt_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],size_t ir_mapping_table[],const int mesh[3],const int is_shift[3],const int is_time_reversal,const MatINT * rotations,const size_t num_q,SPGCONST double qpoints[][3])312 size_t kpt_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],
313 size_t ir_mapping_table[],
314 const int mesh[3],
315 const int is_shift[3],
316 const int is_time_reversal,
317 const MatINT * rotations,
318 const size_t num_q,
319 SPGCONST double qpoints[][3])
320 {
321 size_t num_ir;
322 MatINT *rot_reciprocal, *rot_reciprocal_q;
323 double tolerance;
324
325 rot_reciprocal = NULL;
326 rot_reciprocal_q = NULL;
327
328 rot_reciprocal = get_point_group_reciprocal(rotations, is_time_reversal);
329 tolerance = 0.01 / (mesh[0] + mesh[1] + mesh[2]);
330 rot_reciprocal_q = get_point_group_reciprocal_with_q(rot_reciprocal,
331 tolerance,
332 num_q,
333 qpoints);
334
335 num_ir = get_dense_ir_reciprocal_mesh(grid_address,
336 ir_mapping_table,
337 mesh,
338 is_shift,
339 rot_reciprocal_q);
340
341 mat_free_MatINT(rot_reciprocal_q);
342 rot_reciprocal_q = NULL;
343 mat_free_MatINT(rot_reciprocal);
344 rot_reciprocal = NULL;
345 return num_ir;
346 }
347
348 void
kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],const int address_orig[3],SPGCONST int (* rot_reciprocal)[3][3],const int num_rot,const int mesh[3],const int is_shift[3])349 kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
350 const int address_orig[3],
351 SPGCONST int (*rot_reciprocal)[3][3],
352 const int num_rot,
353 const int mesh[3],
354 const int is_shift[3])
355 {
356 int i;
357 int address_double_orig[3], address_double[3];
358
359 for (i = 0; i < 3; i++) {
360 address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
361 }
362 for (i = 0; i < num_rot; i++) {
363 mat_multiply_matrix_vector_i3(address_double,
364 rot_reciprocal[i],
365 address_double_orig);
366 rot_grid_points[i] = kgd_get_dense_grid_point_double_mesh(address_double, mesh);
367 }
368 }
369
370 void
kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],const int address_orig[3],SPGCONST int (* rot_reciprocal)[3][3],const int num_rot,const int mesh[3],const int is_shift[3],const size_t bz_map[])371 kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
372 const int address_orig[3],
373 SPGCONST int (*rot_reciprocal)[3][3],
374 const int num_rot,
375 const int mesh[3],
376 const int is_shift[3],
377 const size_t bz_map[])
378 {
379 int i;
380 int address_double_orig[3], address_double[3], bzmesh[3];
381
382 for (i = 0; i < 3; i++) {
383 bzmesh[i] = mesh[i] * 2;
384 address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
385 }
386 for (i = 0; i < num_rot; i++) {
387 mat_multiply_matrix_vector_i3(address_double,
388 rot_reciprocal[i],
389 address_double_orig);
390 rot_grid_points[i] =
391 bz_map[kgd_get_dense_grid_point_double_mesh(address_double, bzmesh)];
392 }
393 }
394
kpt_relocate_BZ_grid_address(int bz_grid_address[][3],int bz_map[],SPGCONST int grid_address[][3],const int mesh[3],SPGCONST double rec_lattice[3][3],const int is_shift[3])395 int kpt_relocate_BZ_grid_address(int bz_grid_address[][3],
396 int bz_map[],
397 SPGCONST int grid_address[][3],
398 const int mesh[3],
399 SPGCONST double rec_lattice[3][3],
400 const int is_shift[3])
401 {
402 int i, num_bz_map, num_bzgp;
403 size_t *dense_bz_map;
404
405 num_bz_map = mesh[0] * mesh[1] * mesh[2] * 8;
406
407 if ((dense_bz_map =
408 (size_t*)malloc(sizeof(size_t) * num_bz_map)) == NULL) {
409 warning_print("spglib: Memory of unique_rot could not be allocated.");
410 return 0;
411 }
412
413 num_bzgp = kpt_relocate_dense_BZ_grid_address(bz_grid_address,
414 dense_bz_map,
415 grid_address,
416 mesh,
417 rec_lattice,
418 is_shift);
419
420 for (i = 0; i < num_bz_map; i++) {
421 if (dense_bz_map[i] == num_bz_map) {
422 bz_map[i] = -1;
423 } else {
424 bz_map[i] = dense_bz_map[i];
425 }
426 }
427
428 free(dense_bz_map);
429 dense_bz_map = NULL;
430
431 return num_bzgp;
432 }
433
kpt_relocate_dense_BZ_grid_address(int bz_grid_address[][3],size_t bz_map[],SPGCONST int grid_address[][3],const int mesh[3],SPGCONST double rec_lattice[3][3],const int is_shift[3])434 size_t kpt_relocate_dense_BZ_grid_address(int bz_grid_address[][3],
435 size_t bz_map[],
436 SPGCONST int grid_address[][3],
437 const int mesh[3],
438 SPGCONST double rec_lattice[3][3],
439 const int is_shift[3])
440 {
441 return relocate_dense_BZ_grid_address(bz_grid_address,
442 bz_map,
443 grid_address,
444 mesh,
445 rec_lattice,
446 is_shift);
447 }
448
kpt_get_point_group_reciprocal(const MatINT * rotations,const int is_time_reversal)449 MatINT *kpt_get_point_group_reciprocal(const MatINT * rotations,
450 const int is_time_reversal)
451 {
452 return get_point_group_reciprocal(rotations, is_time_reversal);
453 }
454
kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,const double symprec,const size_t num_q,SPGCONST double qpoints[][3])455 MatINT *kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
456 const double symprec,
457 const size_t num_q,
458 SPGCONST double qpoints[][3])
459 {
460 return get_point_group_reciprocal_with_q(rot_reciprocal,
461 symprec,
462 num_q,
463 qpoints);
464 }
465
466 /* Return NULL if failed */
get_point_group_reciprocal(const MatINT * rotations,const int is_time_reversal)467 static MatINT *get_point_group_reciprocal(const MatINT * rotations,
468 const int is_time_reversal)
469 {
470 int i, j, num_rot;
471 MatINT *rot_reciprocal, *rot_return;
472 int *unique_rot;
473 SPGCONST int inversion[3][3] = {
474 {-1, 0, 0 },
475 { 0,-1, 0 },
476 { 0, 0,-1 }
477 };
478
479 rot_reciprocal = NULL;
480 rot_return = NULL;
481 unique_rot = NULL;
482
483 if (is_time_reversal) {
484 if ((rot_reciprocal = mat_alloc_MatINT(rotations->size * 2)) == NULL) {
485 return NULL;
486 }
487 } else {
488 if ((rot_reciprocal = mat_alloc_MatINT(rotations->size)) == NULL) {
489 return NULL;
490 }
491 }
492
493 if ((unique_rot = (int*)malloc(sizeof(int) * rot_reciprocal->size)) == NULL) {
494 warning_print("spglib: Memory of unique_rot could not be allocated.");
495 mat_free_MatINT(rot_reciprocal);
496 rot_reciprocal = NULL;
497 return NULL;
498 }
499
500 for (i = 0; i < rot_reciprocal->size; i++) {
501 unique_rot[i] = -1;
502 }
503
504 for (i = 0; i < rotations->size; i++) {
505 mat_transpose_matrix_i3(rot_reciprocal->mat[i], rotations->mat[i]);
506
507 if (is_time_reversal) {
508 mat_multiply_matrix_i3(rot_reciprocal->mat[rotations->size+i],
509 inversion,
510 rot_reciprocal->mat[i]);
511 }
512 }
513
514 num_rot = 0;
515 for (i = 0; i < rot_reciprocal->size; i++) {
516 for (j = 0; j < num_rot; j++) {
517 if (mat_check_identity_matrix_i3(rot_reciprocal->mat[unique_rot[j]],
518 rot_reciprocal->mat[i])) {
519 goto escape;
520 }
521 }
522 unique_rot[num_rot] = i;
523 num_rot++;
524 escape:
525 ;
526 }
527
528 if ((rot_return = mat_alloc_MatINT(num_rot)) != NULL) {
529 for (i = 0; i < num_rot; i++) {
530 mat_copy_matrix_i3(rot_return->mat[i], rot_reciprocal->mat[unique_rot[i]]);
531 }
532 }
533
534 free(unique_rot);
535 unique_rot = NULL;
536 mat_free_MatINT(rot_reciprocal);
537 rot_reciprocal = NULL;
538
539 return rot_return;
540 }
541
542 /* Return NULL if failed */
get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,const double symprec,const size_t num_q,SPGCONST double qpoints[][3])543 static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
544 const double symprec,
545 const size_t num_q,
546 SPGCONST double qpoints[][3])
547 {
548 int i, j, k, l, is_all_ok, num_rot;
549 int *ir_rot;
550 double q_rot[3], diff[3];
551 MatINT * rot_reciprocal_q;
552
553 ir_rot = NULL;
554 rot_reciprocal_q = NULL;
555 is_all_ok = 0;
556 num_rot = 0;
557
558 if ((ir_rot = (int*)malloc(sizeof(int) * rot_reciprocal->size)) == NULL) {
559 warning_print("spglib: Memory of ir_rot could not be allocated.");
560 return NULL;
561 }
562
563 for (i = 0; i < rot_reciprocal->size; i++) {
564 ir_rot[i] = -1;
565 }
566 for (i = 0; i < rot_reciprocal->size; i++) {
567 for (j = 0; j < num_q; j++) {
568 is_all_ok = 0;
569 mat_multiply_matrix_vector_id3(q_rot,
570 rot_reciprocal->mat[i],
571 qpoints[j]);
572
573 for (k = 0; k < num_q; k++) {
574 for (l = 0; l < 3; l++) {
575 diff[l] = q_rot[l] - qpoints[k][l];
576 diff[l] -= mat_Nint(diff[l]);
577 }
578
579 if (mat_Dabs(diff[0]) < symprec &&
580 mat_Dabs(diff[1]) < symprec &&
581 mat_Dabs(diff[2]) < symprec) {
582 is_all_ok = 1;
583 break;
584 }
585 }
586
587 if (! is_all_ok) {
588 break;
589 }
590 }
591
592 if (is_all_ok) {
593 ir_rot[num_rot] = i;
594 num_rot++;
595 }
596 }
597
598 if ((rot_reciprocal_q = mat_alloc_MatINT(num_rot)) != NULL) {
599 for (i = 0; i < num_rot; i++) {
600 mat_copy_matrix_i3(rot_reciprocal_q->mat[i],
601 rot_reciprocal->mat[ir_rot[i]]);
602 }
603 }
604
605 free(ir_rot);
606 ir_rot = NULL;
607
608 return rot_reciprocal_q;
609 }
610
get_dense_ir_reciprocal_mesh(int grid_address[][3],size_t ir_mapping_table[],const int mesh[3],const int is_shift[3],const MatINT * rot_reciprocal)611 static size_t get_dense_ir_reciprocal_mesh(int grid_address[][3],
612 size_t ir_mapping_table[],
613 const int mesh[3],
614 const int is_shift[3],
615 const MatINT *rot_reciprocal)
616 {
617 if (check_mesh_symmetry(mesh, is_shift, rot_reciprocal)) {
618 return get_dense_ir_reciprocal_mesh_normal(grid_address,
619 ir_mapping_table,
620 mesh,
621 is_shift,
622 rot_reciprocal);
623 } else {
624 return get_dense_ir_reciprocal_mesh_distortion(grid_address,
625 ir_mapping_table,
626 mesh,
627 is_shift,
628 rot_reciprocal);
629 }
630 }
631
get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],size_t ir_mapping_table[],const int mesh[3],const int is_shift[3],const MatINT * rot_reciprocal)632 static size_t get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
633 size_t ir_mapping_table[],
634 const int mesh[3],
635 const int is_shift[3],
636 const MatINT *rot_reciprocal)
637 {
638 /* In the following loop, mesh is doubled. */
639 /* Even and odd mesh numbers correspond to */
640 /* is_shift[i] are 0 or 1, respectively. */
641 /* is_shift = [0,0,0] gives Gamma center mesh. */
642 /* grid: reducible grid points */
643 /* ir_mapping_table: the mapping from each point to ir-point. */
644
645 size_t i, grid_point_rot;
646 int j;
647 int address_double[3], address_double_rot[3];
648
649 kgd_get_all_grid_addresses(grid_address, mesh);
650
651 #pragma omp parallel for private(j, grid_point_rot, address_double, address_double_rot)
652 for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
653 kgd_get_grid_address_double_mesh(address_double,
654 grid_address[i],
655 mesh,
656 is_shift);
657 ir_mapping_table[i] = i;
658 for (j = 0; j < rot_reciprocal->size; j++) {
659 mat_multiply_matrix_vector_i3(address_double_rot,
660 rot_reciprocal->mat[j],
661 address_double);
662 grid_point_rot = kgd_get_dense_grid_point_double_mesh(address_double_rot, mesh);
663 if (grid_point_rot < ir_mapping_table[i]) {
664 #ifdef _OPENMP
665 ir_mapping_table[i] = grid_point_rot;
666 #else
667 ir_mapping_table[i] = ir_mapping_table[grid_point_rot];
668 break;
669 #endif
670 }
671 }
672 }
673
674 return get_dense_num_ir(ir_mapping_table, mesh);
675 }
676
677 static size_t
get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],size_t ir_mapping_table[],const int mesh[3],const int is_shift[3],const MatINT * rot_reciprocal)678 get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
679 size_t ir_mapping_table[],
680 const int mesh[3],
681 const int is_shift[3],
682 const MatINT *rot_reciprocal)
683 {
684 size_t i, grid_point_rot;
685 int j, k, indivisible;
686 int address_double[3], address_double_rot[3];
687 long long_address_double[3], long_address_double_rot[3], divisor[3];
688
689 /* divisor, long_address_double, and long_address_double_rot have */
690 /* long integer type to treat dense mesh. */
691
692 kgd_get_all_grid_addresses(grid_address, mesh);
693
694 for (j = 0; j < 3; j++) {
695 divisor[j] = mesh[(j + 1) % 3] * mesh[(j + 2) % 3];
696 }
697
698 #pragma omp parallel for private(j, k, grid_point_rot, address_double, address_double_rot, long_address_double, long_address_double_rot)
699 for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
700 kgd_get_grid_address_double_mesh(address_double,
701 grid_address[i],
702 mesh,
703 is_shift);
704 for (j = 0; j < 3; j++) {
705 long_address_double[j] = address_double[j] * divisor[j];
706 }
707 ir_mapping_table[i] = i;
708 for (j = 0; j < rot_reciprocal->size; j++) {
709
710 /* Equivalent to mat_multiply_matrix_vector_i3 except for data type */
711 for (k = 0; k < 3; k++) {
712 long_address_double_rot[k] =
713 rot_reciprocal->mat[j][k][0] * long_address_double[0] +
714 rot_reciprocal->mat[j][k][1] * long_address_double[1] +
715 rot_reciprocal->mat[j][k][2] * long_address_double[2];
716 }
717
718 for (k = 0; k < 3; k++) {
719 indivisible = long_address_double_rot[k] % divisor[k];
720 if (indivisible) {break;}
721 address_double_rot[k] = long_address_double_rot[k] / divisor[k];
722 if ((address_double_rot[k] % 2 != 0 && is_shift[k] == 0) ||
723 (address_double_rot[k] % 2 == 0 && is_shift[k] == 1)) {
724 indivisible = 1;
725 break;
726 }
727 }
728 if (indivisible) {continue;}
729 grid_point_rot =
730 kgd_get_dense_grid_point_double_mesh(address_double_rot, mesh);
731 if (grid_point_rot < ir_mapping_table[i]) {
732 #ifdef _OPENMP
733 ir_mapping_table[i] = grid_point_rot;
734 #else
735 ir_mapping_table[i] = ir_mapping_table[grid_point_rot];
736 break;
737 #endif
738 }
739 }
740 }
741
742 return get_dense_num_ir(ir_mapping_table, mesh);
743 }
744
get_dense_num_ir(size_t ir_mapping_table[],const int mesh[3])745 static size_t get_dense_num_ir(size_t ir_mapping_table[], const int mesh[3])
746 {
747 size_t i, num_ir;
748
749 num_ir = 0;
750
751 #pragma omp parallel for reduction(+:num_ir)
752 for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
753 if (ir_mapping_table[i] == i) {
754 num_ir++;
755 }
756 }
757
758 #ifdef _OPENMP
759 for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
760 ir_mapping_table[i] = ir_mapping_table[ir_mapping_table[i]];
761 }
762 #endif
763
764 return num_ir;
765 }
766
relocate_dense_BZ_grid_address(int bz_grid_address[][3],size_t bz_map[],SPGCONST int grid_address[][3],const int mesh[3],SPGCONST double rec_lattice[3][3],const int is_shift[3])767 static size_t relocate_dense_BZ_grid_address(int bz_grid_address[][3],
768 size_t bz_map[],
769 SPGCONST int grid_address[][3],
770 const int mesh[3],
771 SPGCONST double rec_lattice[3][3],
772 const int is_shift[3])
773 {
774 double tolerance, min_distance;
775 double q_vector[3], distance[KPT_NUM_BZ_SEARCH_SPACE];
776 int bzmesh[3], bz_address_double[3];
777 size_t i, boundary_num_gp, total_num_gp, bzgp, gp, num_bzmesh;
778 int j, k, min_index;
779
780 tolerance = get_tolerance_for_BZ_reduction(rec_lattice, mesh);
781 for (j = 0; j < 3; j++) {
782 bzmesh[j] = mesh[j] * 2;
783 }
784
785 num_bzmesh = bzmesh[0] * bzmesh[1] * (size_t)(bzmesh[2]);
786 for (i = 0; i < num_bzmesh; i++) {
787 bz_map[i] = num_bzmesh;
788 }
789
790 boundary_num_gp = 0;
791 total_num_gp = mesh[0] * mesh[1] * (size_t)(mesh[2]);
792
793 /* Multithreading doesn't work for this loop since gp calculated */
794 /* with boundary_num_gp is unstable to store bz_grid_address. */
795 for (i = 0; i < total_num_gp; i++) {
796 for (j = 0; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
797 for (k = 0; k < 3; k++) {
798 q_vector[k] =
799 ((grid_address[i][k] + bz_search_space[j][k] * mesh[k]) * 2 +
800 is_shift[k]) / ((double)mesh[k]) / 2;
801 }
802 mat_multiply_matrix_vector_d3(q_vector, rec_lattice, q_vector);
803 distance[j] = mat_norm_squared_d3(q_vector);
804 }
805 min_distance = distance[0];
806 min_index = 0;
807 for (j = 1; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
808 if (distance[j] < min_distance) {
809 min_distance = distance[j];
810 min_index = j;
811 }
812 }
813
814 for (j = 0; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
815 if (distance[j] < min_distance + tolerance) {
816 if (j == min_index) {
817 gp = i;
818 } else {
819 gp = boundary_num_gp + total_num_gp;
820 }
821
822 for (k = 0; k < 3; k++) {
823 bz_grid_address[gp][k] =
824 grid_address[i][k] + bz_search_space[j][k] * mesh[k];
825 bz_address_double[k] = bz_grid_address[gp][k] * 2 + is_shift[k];
826 }
827 bzgp = kgd_get_dense_grid_point_double_mesh(bz_address_double, bzmesh);
828 bz_map[bzgp] = gp;
829 if (j != min_index) {
830 boundary_num_gp++;
831 }
832 }
833 }
834 }
835
836 return boundary_num_gp + total_num_gp;
837 }
838
get_tolerance_for_BZ_reduction(SPGCONST double rec_lattice[3][3],const int mesh[3])839 static double get_tolerance_for_BZ_reduction(SPGCONST double rec_lattice[3][3],
840 const int mesh[3])
841 {
842 int i, j;
843 double tolerance;
844 double length[3];
845
846 for (i = 0; i < 3; i++) {
847 length[i] = 0;
848 for (j = 0; j < 3; j++) {
849 length[i] += rec_lattice[j][i] * rec_lattice[j][i];
850 }
851 length[i] /= mesh[i] * mesh[i];
852 }
853 tolerance = length[0];
854 for (i = 1; i < 3; i++) {
855 if (tolerance < length[i]) {
856 tolerance = length[i];
857 }
858 }
859 tolerance *= 0.01;
860
861 return tolerance;
862 }
863
check_mesh_symmetry(const int mesh[3],const int is_shift[3],const MatINT * rot_reciprocal)864 static int check_mesh_symmetry(const int mesh[3],
865 const int is_shift[3],
866 const MatINT *rot_reciprocal)
867 {
868 int i, j, k, sum;
869 int eq[3];
870
871 eq[0] = 0; /* a=b */
872 eq[1] = 0; /* b=c */
873 eq[2] = 0; /* c=a */
874
875 /* Check 3 and 6 fold rotations and non-convensional choice of unit cells */
876 for (i = 0; i < rot_reciprocal->size; i++) {
877 sum = 0;
878 for (j = 0; j < 3; j++) {
879 for (k = 0; k < 3; k++) {
880 sum += abs(rot_reciprocal->mat[i][j][k]);
881 }
882 }
883 if (sum > 3) {
884 return 0;
885 }
886 }
887
888 for (i = 0; i < rot_reciprocal->size; i++) {
889 if (rot_reciprocal->mat[i][0][0] == 0 &&
890 rot_reciprocal->mat[i][1][0] == 1 &&
891 rot_reciprocal->mat[i][2][0] == 0) {eq[0] = 1;}
892 if (rot_reciprocal->mat[i][0][0] == 0 &&
893 rot_reciprocal->mat[i][1][0] == 1 &&
894 rot_reciprocal->mat[i][2][0] == 0) {eq[1] = 1;}
895 if (rot_reciprocal->mat[i][0][0] == 0 &&
896 rot_reciprocal->mat[i][1][0] == 0 &&
897 rot_reciprocal->mat[i][2][0] == 1) {eq[2] = 1;}
898 }
899
900
901 return (((eq[0] && mesh[0] == mesh[1] && is_shift[0] == is_shift[1]) || (!eq[0])) &&
902 ((eq[1] && mesh[1] == mesh[2] && is_shift[1] == is_shift[2]) || (!eq[1])) &&
903 ((eq[2] && mesh[2] == mesh[0] && is_shift[2] == is_shift[0]) || (!eq[2])));
904 }
905