1 /* Copyright (C) 2015 Atsushi Togo */
2 /* All rights reserved. */
3
4 /* These codes were originally parts of spglib, but only develped */
5 /* and used for phono3py. Therefore these were moved from spglib to */
6 /* phono3py. This file is part of phonopy. */
7
8 /* Redistribution and use in source and binary forms, with or without */
9 /* modification, are permitted provided that the following conditions */
10 /* are met: */
11
12 /* * Redistributions of source code must retain the above copyright */
13 /* notice, this list of conditions and the following disclaimer. */
14
15 /* * Redistributions in binary form must reproduce the above copyright */
16 /* notice, this list of conditions and the following disclaimer in */
17 /* the documentation and/or other materials provided with the */
18 /* distribution. */
19
20 /* * Neither the name of the phonopy project nor the names of its */
21 /* contributors may be used to endorse or promote products derived */
22 /* from this software without specific prior written permission. */
23
24 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
25 /* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
26 /* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
27 /* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
28 /* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
29 /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
30 /* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
31 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
32 /* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
33 /* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
34 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
35 /* POSSIBILITY OF SUCH DAMAGE. */
36
37 #include <stddef.h>
38 #include <mathfunc.h>
39 #include "triplet.h"
40 #include "triplet_iw.h"
41 #include "triplet_kpoint.h"
42
43 static size_t get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
44 size_t *map_q,
45 int (*grid_address)[3],
46 const int grid_point,
47 const int mesh[3],
48 const int is_time_reversal,
49 const int num_rot,
50 TPLCONST int (*rotations)[3][3],
51 const int swappable);
52
tpl_get_BZ_triplets_at_q(size_t (* triplets)[3],const size_t grid_point,TPLCONST int (* bz_grid_address)[3],const size_t * bz_map,const size_t * map_triplets,const size_t num_map_triplets,const int mesh[3])53 size_t tpl_get_BZ_triplets_at_q(size_t (*triplets)[3],
54 const size_t grid_point,
55 TPLCONST int (*bz_grid_address)[3],
56 const size_t *bz_map,
57 const size_t *map_triplets,
58 const size_t num_map_triplets,
59 const int mesh[3])
60 {
61 return tpk_get_BZ_triplets_at_q(triplets,
62 grid_point,
63 bz_grid_address,
64 bz_map,
65 map_triplets,
66 num_map_triplets,
67 mesh);
68 }
69
tpl_get_triplets_reciprocal_mesh_at_q(size_t * map_triplets,size_t * map_q,int (* grid_address)[3],const size_t grid_point,const int mesh[3],const int is_time_reversal,const int num_rot,TPLCONST int (* rotations)[3][3],const int swappable)70 size_t tpl_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
71 size_t *map_q,
72 int (*grid_address)[3],
73 const size_t grid_point,
74 const int mesh[3],
75 const int is_time_reversal,
76 const int num_rot,
77 TPLCONST int (*rotations)[3][3],
78 const int swappable)
79 {
80 return get_triplets_reciprocal_mesh_at_q(map_triplets,
81 map_q,
82 grid_address,
83 grid_point,
84 mesh,
85 is_time_reversal,
86 num_rot,
87 rotations,
88 swappable);
89 }
90
tpl_get_integration_weight(double * iw,char * iw_zero,const double * frequency_points,const size_t num_band0,TPLCONST int relative_grid_address[24][4][3],const int mesh[3],TPLCONST size_t (* triplets)[3],const size_t num_triplets,TPLCONST int (* bz_grid_address)[3],const size_t * bz_map,const double * frequencies1,const size_t num_band1,const double * frequencies2,const size_t num_band2,const size_t tp_type,const int openmp_per_triplets,const int openmp_per_bands)91 void tpl_get_integration_weight(double *iw,
92 char *iw_zero,
93 const double *frequency_points,
94 const size_t num_band0,
95 TPLCONST int relative_grid_address[24][4][3],
96 const int mesh[3],
97 TPLCONST size_t (*triplets)[3],
98 const size_t num_triplets,
99 TPLCONST int (*bz_grid_address)[3],
100 const size_t *bz_map,
101 const double *frequencies1,
102 const size_t num_band1,
103 const double *frequencies2,
104 const size_t num_band2,
105 const size_t tp_type,
106 const int openmp_per_triplets,
107 const int openmp_per_bands)
108 {
109 size_t i, num_band_prod;
110 int tp_relative_grid_address[2][24][4][3];
111
112 tpl_set_relative_grid_address(tp_relative_grid_address,
113 relative_grid_address,
114 tp_type);
115 num_band_prod = num_band0 * num_band1 * num_band2;
116
117 #pragma omp parallel for if (openmp_per_triplets)
118 for (i = 0; i < num_triplets; i++) {
119 tpi_get_integration_weight(iw + i * num_band_prod,
120 iw_zero + i * num_band_prod,
121 frequency_points, /* f0 */
122 num_band0,
123 tp_relative_grid_address,
124 mesh,
125 triplets[i],
126 num_triplets,
127 bz_grid_address,
128 bz_map,
129 frequencies1, /* f1 */
130 num_band1,
131 frequencies2, /* f2 */
132 num_band2,
133 tp_type,
134 openmp_per_bands);
135 }
136 }
137
138
tpl_get_integration_weight_with_sigma(double * iw,char * iw_zero,const double sigma,const double sigma_cutoff,const double * frequency_points,const size_t num_band0,TPLCONST size_t (* triplets)[3],const size_t num_triplets,const double * frequencies,const size_t num_band,const size_t tp_type)139 void tpl_get_integration_weight_with_sigma(double *iw,
140 char *iw_zero,
141 const double sigma,
142 const double sigma_cutoff,
143 const double *frequency_points,
144 const size_t num_band0,
145 TPLCONST size_t (*triplets)[3],
146 const size_t num_triplets,
147 const double *frequencies,
148 const size_t num_band,
149 const size_t tp_type)
150 {
151 size_t i, num_band_prod, const_adrs_shift;
152 double cutoff;
153
154 cutoff = sigma * sigma_cutoff;
155 num_band_prod = num_band0 * num_band * num_band;
156 const_adrs_shift = num_triplets * num_band0 * num_band * num_band;
157
158 #pragma omp parallel for
159 for (i = 0; i < num_triplets; i++) {
160 tpi_get_integration_weight_with_sigma(
161 iw + i * num_band_prod,
162 iw_zero + i * num_band_prod,
163 sigma,
164 cutoff,
165 frequency_points,
166 num_band0,
167 triplets[i],
168 const_adrs_shift,
169 frequencies,
170 num_band,
171 tp_type,
172 0);
173 }
174 }
175
176
tpl_is_N(const size_t triplet[3],const int * grid_address)177 int tpl_is_N(const size_t triplet[3], const int *grid_address)
178 {
179 int i, j, sum_q, is_N;
180
181 is_N = 1;
182 for (i = 0; i < 3; i++) {
183 sum_q = 0;
184 for (j = 0; j < 3; j++) { /* 1st, 2nd, 3rd triplet */
185 sum_q += grid_address[triplet[j] * 3 + i];
186 }
187 if (sum_q) {
188 is_N = 0;
189 break;
190 }
191 }
192 return is_N;
193 }
194
tpl_set_relative_grid_address(int tp_relative_grid_address[2][24][4][3],TPLCONST int relative_grid_address[24][4][3],const size_t tp_type)195 void tpl_set_relative_grid_address(
196 int tp_relative_grid_address[2][24][4][3],
197 TPLCONST int relative_grid_address[24][4][3],
198 const size_t tp_type)
199 {
200 size_t i, j, k, l;
201 int signs[2];
202
203 signs[0] = 1;
204 signs[1] = 1;
205 if ((tp_type == 2) || (tp_type == 3)) {
206 /* q1+q2+q3=G */
207 /* To set q2+1, q3-1 is needed to keep G */
208 signs[1] = -1;
209 }
210 /* tp_type == 4, q+k_i-k_f=G */
211
212 for (i = 0; i < 2; i++) {
213 for (j = 0; j < 24; j++) {
214 for (k = 0; k < 4; k++) {
215 for (l = 0; l < 3; l++) {
216 tp_relative_grid_address[i][j][k][l] =
217 relative_grid_address[j][k][l] * signs[i];
218 }
219 }
220 }
221 }
222 }
223
get_triplets_reciprocal_mesh_at_q(size_t * map_triplets,size_t * map_q,int (* grid_address)[3],const int grid_point,const int mesh[3],const int is_time_reversal,const int num_rot,TPLCONST int (* rotations)[3][3],const int swappable)224 static size_t get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
225 size_t *map_q,
226 int (*grid_address)[3],
227 const int grid_point,
228 const int mesh[3],
229 const int is_time_reversal,
230 const int num_rot,
231 TPLCONST int (*rotations)[3][3],
232 const int swappable)
233 {
234 MatINT *rot_real;
235 int i;
236 size_t num_ir;
237
238 rot_real = mat_alloc_MatINT(num_rot);
239 for (i = 0; i < num_rot; i++) {
240 mat_copy_matrix_i3(rot_real->mat[i], rotations[i]);
241 }
242
243 num_ir = tpk_get_ir_triplets_at_q(map_triplets,
244 map_q,
245 grid_address,
246 grid_point,
247 mesh,
248 is_time_reversal,
249 rot_real,
250 swappable);
251
252 mat_free_MatINT(rot_real);
253
254 return num_ir;
255 }
256