Branch data Line data Source code
1 : : !!****m* ABINIT/m_htetra
2 : : !! NAME
3 : : !! m_htetra
4 : : !!
5 : : !! FUNCTION
6 : : !! Module for tetrahedron integration of DOS and similar quantities
7 : : !! Uses some functions from a previous implementation by MJV
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2010-2026 ABINIT group (HM,MJV)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! TODO
16 : : !! 1) Test more carefully the case of degenerate tethraedron
17 : : !! 2) Add options to get only delta and/or theta ?
18 : : !!
19 : : !! SOURCE
20 : :
21 : : #if defined HAVE_CONFIG_H
22 : : #include "config.h"
23 : : #endif
24 : :
25 : : #include "abi_common.h"
26 : :
27 : : module m_htetra
28 : :
29 : : use defs_basis
30 : : use m_abicore
31 : : use m_krank
32 : : use m_xmpi
33 : : use m_errors
34 : :
35 : : use m_fstrings, only : sjoin, itoa, ftoa
36 : : use m_numeric_tools, only : linspace
37 : : use m_simtet, only : sim0onei, SIM0TWOI
38 : :
39 : : implicit none
40 : :
41 : : private
42 : : !!***
43 : :
44 : : integer, parameter :: TETRA_SIZE = 6
45 : : integer, parameter :: TETRA_STEP = 6
46 : :
47 : : !!****t* m_htetra/t_htetra_bucket
48 : : !! NAME
49 : : !! t_htetra_bucket
50 : : !!
51 : : !! FUNCTION
52 : : !! Store a bunch of tetrahedra
53 : : !!
54 : : !! SOURCE
55 : :
56 : : type :: htetra_bucket
57 : :
58 : : integer,allocatable :: indexes(:,:)
59 : :
60 : : end type htetra_bucket
61 : : !!***
62 : :
63 : : !!****t* m_htetra/htetra_t
64 : : !! NAME
65 : : !! htetra_t
66 : : !!
67 : : !! FUNCTION
68 : : !! tetrahedron geometry object
69 : : !!
70 : : !! SOURCE
71 : :
72 : : type, public :: htetra_t
73 : :
74 : : integer :: opt
75 : : ! Option for the generation of tetrahedra
76 : :
77 : : integer :: nkibz
78 : : ! Number of points in the irreducible Brillouin zone
79 : :
80 : : integer :: nkbz
81 : : ! Number of points in the full Brillouin zone
82 : :
83 : : integer :: nbuckets
84 : : ! Number of buckets for the hash table
85 : :
86 : : integer :: nunique_tetra
87 : : ! Number of unique tetrahedron
88 : :
89 : : integer :: nibz_tetra
90 : : ! Number of ibz tetrahedron
91 : :
92 : : integer,allocatable :: tetra_total(:)
93 : : ! (%nbkibz)
94 : : ! Equivalent tetrahedra per kpoint (number of irred tetra times multiplicity)
95 : :
96 : : integer,allocatable :: tetra_count(:)
97 : : ! (%nbkibz)
98 : : ! Inequivalent tetrahedra per kpoint (number of irred tetra)
99 : :
100 : : integer,allocatable :: ibz_multiplicity(:)
101 : : ! (%nbkibz)
102 : : ! Multiplicity of each k-point
103 : :
104 : : real(dp) :: vv
105 : : ! volume of the tetrahedra
106 : :
107 : : real(dp) :: klatt(3, 3)
108 : : ! reciprocal of equivalent supercell lattice vectors for full kpoint grid
109 : :
110 : : type(htetra_bucket),allocatable :: ibz(:)
111 : : ! indexes of the tetrahedra for each k-point
112 : :
113 : : type(htetra_bucket),allocatable :: unique_tetra(:)
114 : : ! indexes of the unique tetrahedra
115 : :
116 : : contains
117 : :
118 : : procedure :: init => htetra_init ! Initialize the object
119 : :
120 : : procedure :: free => htetra_free
121 : : ! Free memory
122 : :
123 : : procedure :: print => htetra_print
124 : : ! Print information about tetrahedron object
125 : :
126 : : procedure :: get_onewk => htetra_get_onewk
127 : : ! Calculate integration weights and their derivatives for a single k-point in the IBZ.
128 : :
129 : : procedure :: get_onewk_wvals => htetra_get_onewk_wvals
130 : : ! Similar to tetra_get_onewk but receives arbitrary list of frequency points.
131 : :
132 : : procedure :: get_onewk_wvals_zinv => htetra_get_onewk_wvals_zinv
133 : : ! Calculate integration weights for 1/(z-E(k)) for a single k-point in the IBZ.
134 : :
135 : : procedure :: weights_wvals_zinv => htetra_weights_wvals_zinv
136 : : ! Same as above but return the weight on all the kpoints by looping over tetrahedra
137 : :
138 : : procedure :: wvals_weights => htetra_wvals_weights
139 : : ! Compute delta and theta on a list of energies for all kpoints
140 : :
141 : : procedure :: wvals_weights_delta => htetra_wvals_weights_delta
142 : : ! Compute delta on a list of energies for all kpoints
143 : :
144 : : procedure :: blochl_weights => htetra_blochl_weights
145 : : ! And interface to help to facilitate the transition to the new tetrahedron implementation
146 : :
147 : : end type htetra_t
148 : : !!***
149 : :
150 : : !!***
151 : :
152 : : contains
153 : : !!***
154 : :
155 : : !----------------------------------------------------------------------
156 : :
157 : : !!****f* m_htetra/htetra_init
158 : : !! NAME
159 : : !! htetra_init
160 : : !!
161 : : !! FUNCTION
162 : : !! get tetrahedra characterized by apexes
163 : : !!
164 : : !! INPUTS
165 : : !! bz2ibz(nkpt_fullbz)=indexes of irred kpoints equivalent to kpt_fullbz
166 : : !! gprimd(3,3) = reciprocal space vectors
167 : : !! klatt(3,3)=reciprocal of lattice vectors for full kpoint grid
168 : : !! kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone
169 : : !! nkpt_fullbz=number of kpoints in full brillouin zone
170 : : !! options=1.generate 24 tetrahedra per k-point
171 : : !! faster but gives different results depending on the IBZ, small error for large grids
172 : : !! 2.generate tetrahedra on the FBZ and map to IBZ
173 : : !! slower but same results for IBZ and FBZ.
174 : : !! comm= MPI communicator
175 : : !! [opt]= 1 for Togo's version, 2 for Blochl's version (default)
176 : : !!
177 : : !! OUTPUT
178 : : !! tetra%ibz(4,24,nkibz)=for each k-point, the indexes in the IBZ
179 : : !! tetra%vv = tetrahedron volume divided by full BZ volume
180 : : !!
181 : : !! SOURCE
182 : :
183 : 486 : subroutine htetra_init(tetra, bz2ibz, gprimd, klatt, kpt_fullbz, nkpt_fullbz, kpt_ibz, nkpt_ibz, &
184 : : ierr, errorstring, comm, &
185 : : opt) ! optional
186 : :
187 : : !Arguments ------------------------------------
188 : : !scalars
189 : : class(htetra_t),intent(out),target :: tetra
190 : : integer,intent(in) :: nkpt_fullbz, nkpt_ibz, comm
191 : : integer,optional,intent(in) :: opt
192 : : integer,intent(out) :: ierr
193 : : character(len=80),intent(out) :: errorstring
194 : : !arrays
195 : : integer,intent(in) :: bz2ibz(nkpt_fullbz)
196 : : real(dp),intent(in) :: gprimd(3,3),klatt(3,3),kpt_fullbz(3,nkpt_fullbz),kpt_ibz(3,nkpt_ibz)
197 : :
198 : : !Local variables-------------------------------
199 : : !scalars
200 : : !type(octree_t) :: oct
201 : 243 : type(krank_t) :: krank
202 : : integer :: ikpt2,isummit,itetra,jtetra
203 : : integer :: ikibz,ikbz,idiag,ihash,min_idiag,my_rank,nprocs
204 : : integer :: max_ntetra, ntetra
205 : : real(dp) :: rcvol,length,min_length
206 : : !arrays
207 [ + - ]: 243 : integer,allocatable,target :: indexes(:,:), tetra_hash_count(:)
208 : : integer :: tetra_ibz(4)
209 : : integer :: tetra_shifts(3,4,24,4) ! 3 dimensions, 4 summits, 24 tetrahedra, 4 main diagonals
210 : : integer :: tetra_shifts_6(3,4,6,1) ! 3 dimensions, 4 summits, 6 tetrahedra, 4 main diagonals
211 : : integer :: main_diagonals(3,4), tetra_mibz(0:4)
212 : : real(dp) :: k1(3),k2(3),k3(3),diag(3)
213 : : ! *********************************************************************
214 : :
215 : : ! Use the shifts from kpclib developed by Atsushi Togo
216 : : ! This part is produced by a python script
217 : : ! This implementation is based on spglib and kpclib by Atsushi Togo
218 : : ! after a discussion with hin on the APS 2019 where he provided
219 : : ! details of his implementation.
220 : : ! Note that we don't use it in production as we found that is approach, although faster than the original
221 : : ! one proposed by Blochl (and implemented by MJV) does not preserve symmetries that is calculations done on the full BZ
222 : : ! and the IBZ do not produce the same result. The diff, however, decreases if the sampling is densified.
223 : :
224 [ + + ]: 972 : tetra_shifts(:, 1, 1,1) = [ 0, 0, 0]
225 [ + + ]: 972 : tetra_shifts(:, 2, 1,1) = [ 1, 0, 0]
226 [ + + ]: 972 : tetra_shifts(:, 3, 1,1) = [ 1, 1, 0]
227 [ + + ]: 972 : tetra_shifts(:, 4, 1,1) = [ 1, 1, 1]
228 [ + + ]: 972 : tetra_shifts(:, 1, 2,1) = [ 0, 0, 0]
229 [ + + ]: 972 : tetra_shifts(:, 2, 2,1) = [ 1, 0, 0]
230 [ + + ]: 972 : tetra_shifts(:, 3, 2,1) = [ 1, 0, 1]
231 [ + + ]: 972 : tetra_shifts(:, 4, 2,1) = [ 1, 1, 1]
232 [ + + ]: 972 : tetra_shifts(:, 1, 3,1) = [ 0, 0, 0]
233 [ + + ]: 972 : tetra_shifts(:, 2, 3,1) = [ 0, 1, 0]
234 [ + + ]: 972 : tetra_shifts(:, 3, 3,1) = [ 1, 1, 0]
235 [ + + ]: 972 : tetra_shifts(:, 4, 3,1) = [ 1, 1, 1]
236 [ + + ]: 972 : tetra_shifts(:, 1, 4,1) = [ 0, 0, 0]
237 [ + + ]: 972 : tetra_shifts(:, 2, 4,1) = [ 0, 1, 0]
238 [ + + ]: 972 : tetra_shifts(:, 3, 4,1) = [ 0, 1, 1]
239 [ + + ]: 972 : tetra_shifts(:, 4, 4,1) = [ 1, 1, 1]
240 [ + + ]: 972 : tetra_shifts(:, 1, 5,1) = [ 0, 0, 0]
241 [ + + ]: 972 : tetra_shifts(:, 2, 5,1) = [ 0, 0, 1]
242 [ + + ]: 972 : tetra_shifts(:, 3, 5,1) = [ 1, 0, 1]
243 [ + + ]: 972 : tetra_shifts(:, 4, 5,1) = [ 1, 1, 1]
244 [ + + ]: 972 : tetra_shifts(:, 1, 6,1) = [ 0, 0, 0]
245 [ + + ]: 972 : tetra_shifts(:, 2, 6,1) = [ 0, 0, 1]
246 [ + + ]: 972 : tetra_shifts(:, 3, 6,1) = [ 0, 1, 1]
247 [ + + ]: 972 : tetra_shifts(:, 4, 6,1) = [ 1, 1, 1]
248 [ + + ]: 972 : tetra_shifts(:, 1, 7,1) = [ 0, 0, 0]
249 [ + + ]: 972 : tetra_shifts(:, 2, 7,1) = [ 0, 1, 0]
250 [ + + ]: 972 : tetra_shifts(:, 3, 7,1) = [ 0, 1, 1]
251 [ + + ]: 972 : tetra_shifts(:, 4, 7,1) = [ -1, 0, 0]
252 [ + + ]: 972 : tetra_shifts(:, 1, 8,1) = [ 0, 0, 0]
253 [ + + ]: 972 : tetra_shifts(:, 2, 8,1) = [ 0, 0, 1]
254 [ + + ]: 972 : tetra_shifts(:, 3, 8,1) = [ 0, 1, 1]
255 [ + + ]: 972 : tetra_shifts(:, 4, 8,1) = [ -1, 0, 0]
256 [ + + ]: 972 : tetra_shifts(:, 1, 9,1) = [ 0, 0, 0]
257 [ + + ]: 972 : tetra_shifts(:, 2, 9,1) = [ 1, 0, 0]
258 [ + + ]: 972 : tetra_shifts(:, 3, 9,1) = [ 1, 0, 1]
259 [ + + ]: 972 : tetra_shifts(:, 4, 9,1) = [ 0, -1, 0]
260 [ + + ]: 972 : tetra_shifts(:, 1,10,1) = [ 0, 0, 0]
261 [ + + ]: 972 : tetra_shifts(:, 2,10,1) = [ 0, 0, 1]
262 [ + + ]: 972 : tetra_shifts(:, 3,10,1) = [ 1, 0, 1]
263 [ + + ]: 972 : tetra_shifts(:, 4,10,1) = [ 0, -1, 0]
264 [ + + ]: 972 : tetra_shifts(:, 1,11,1) = [ 0, 0, 0]
265 [ + + ]: 972 : tetra_shifts(:, 2,11,1) = [ 0, 0, 1]
266 [ + + ]: 972 : tetra_shifts(:, 3,11,1) = [ -1, -1, 0]
267 [ + + ]: 972 : tetra_shifts(:, 4,11,1) = [ 0, -1, 0]
268 [ + + ]: 972 : tetra_shifts(:, 1,12,1) = [ 0, 0, 0]
269 [ + + ]: 972 : tetra_shifts(:, 2,12,1) = [ 0, 0, 1]
270 [ + + ]: 972 : tetra_shifts(:, 3,12,1) = [ -1, -1, 0]
271 [ + + ]: 972 : tetra_shifts(:, 4,12,1) = [ -1, 0, 0]
272 [ + + ]: 972 : tetra_shifts(:, 1,13,1) = [ 0, 0, 0]
273 [ + + ]: 972 : tetra_shifts(:, 2,13,1) = [ 1, 0, 0]
274 [ + + ]: 972 : tetra_shifts(:, 3,13,1) = [ 1, 1, 0]
275 [ + + ]: 972 : tetra_shifts(:, 4,13,1) = [ 0, 0, -1]
276 [ + + ]: 972 : tetra_shifts(:, 1,14,1) = [ 0, 0, 0]
277 [ + + ]: 972 : tetra_shifts(:, 2,14,1) = [ 0, 1, 0]
278 [ + + ]: 972 : tetra_shifts(:, 3,14,1) = [ 1, 1, 0]
279 [ + + ]: 972 : tetra_shifts(:, 4,14,1) = [ 0, 0, -1]
280 [ + + ]: 972 : tetra_shifts(:, 1,15,1) = [ 0, 0, 0]
281 [ + + ]: 972 : tetra_shifts(:, 2,15,1) = [ 0, 1, 0]
282 [ + + ]: 972 : tetra_shifts(:, 3,15,1) = [ -1, 0, -1]
283 [ + + ]: 972 : tetra_shifts(:, 4,15,1) = [ 0, 0, -1]
284 [ + + ]: 972 : tetra_shifts(:, 1,16,1) = [ 0, 0, 0]
285 [ + + ]: 972 : tetra_shifts(:, 2,16,1) = [ 0, 1, 0]
286 [ + + ]: 972 : tetra_shifts(:, 3,16,1) = [ -1, 0, -1]
287 [ + + ]: 972 : tetra_shifts(:, 4,16,1) = [ -1, 0, 0]
288 [ + + ]: 972 : tetra_shifts(:, 1,17,1) = [ 0, 0, 0]
289 [ + + ]: 972 : tetra_shifts(:, 2,17,1) = [ 1, 0, 0]
290 [ + + ]: 972 : tetra_shifts(:, 3,17,1) = [ 0, -1, -1]
291 [ + + ]: 972 : tetra_shifts(:, 4,17,1) = [ 0, 0, -1]
292 [ + + ]: 972 : tetra_shifts(:, 1,18,1) = [ 0, 0, 0]
293 [ + + ]: 972 : tetra_shifts(:, 2,18,1) = [ 1, 0, 0]
294 [ + + ]: 972 : tetra_shifts(:, 3,18,1) = [ 0, -1, -1]
295 [ + + ]: 972 : tetra_shifts(:, 4,18,1) = [ 0, -1, 0]
296 [ + + ]: 972 : tetra_shifts(:, 1,19,1) = [ 0, 0, 0]
297 [ + + ]: 972 : tetra_shifts(:, 2,19,1) = [ -1, -1, -1]
298 [ + + ]: 972 : tetra_shifts(:, 3,19,1) = [ 0, -1, -1]
299 [ + + ]: 972 : tetra_shifts(:, 4,19,1) = [ 0, 0, -1]
300 [ + + ]: 972 : tetra_shifts(:, 1,20,1) = [ 0, 0, 0]
301 [ + + ]: 972 : tetra_shifts(:, 2,20,1) = [ -1, -1, -1]
302 [ + + ]: 972 : tetra_shifts(:, 3,20,1) = [ 0, -1, -1]
303 [ + + ]: 972 : tetra_shifts(:, 4,20,1) = [ 0, -1, 0]
304 [ + + ]: 972 : tetra_shifts(:, 1,21,1) = [ 0, 0, 0]
305 [ + + ]: 972 : tetra_shifts(:, 2,21,1) = [ -1, -1, -1]
306 [ + + ]: 972 : tetra_shifts(:, 3,21,1) = [ -1, 0, -1]
307 [ + + ]: 972 : tetra_shifts(:, 4,21,1) = [ 0, 0, -1]
308 [ + + ]: 972 : tetra_shifts(:, 1,22,1) = [ 0, 0, 0]
309 [ + + ]: 972 : tetra_shifts(:, 2,22,1) = [ -1, -1, -1]
310 [ + + ]: 972 : tetra_shifts(:, 3,22,1) = [ -1, 0, -1]
311 [ + + ]: 972 : tetra_shifts(:, 4,22,1) = [ -1, 0, 0]
312 [ + + ]: 972 : tetra_shifts(:, 1,23,1) = [ 0, 0, 0]
313 [ + + ]: 972 : tetra_shifts(:, 2,23,1) = [ -1, -1, -1]
314 [ + + ]: 972 : tetra_shifts(:, 3,23,1) = [ -1, -1, 0]
315 [ + + ]: 972 : tetra_shifts(:, 4,23,1) = [ 0, -1, 0]
316 [ + + ]: 972 : tetra_shifts(:, 1,24,1) = [ 0, 0, 0]
317 [ + + ]: 972 : tetra_shifts(:, 2,24,1) = [ -1, -1, -1]
318 [ + + ]: 972 : tetra_shifts(:, 3,24,1) = [ -1, -1, 0]
319 [ + + ]: 972 : tetra_shifts(:, 4,24,1) = [ -1, 0, 0]
320 [ + + ]: 972 : tetra_shifts(:, 1, 1,2) = [ 0, 0, 0]
321 [ + + ]: 972 : tetra_shifts(:, 2, 1,2) = [ 1, 0, 0]
322 [ + + ]: 972 : tetra_shifts(:, 3, 1,2) = [ 0, 1, 0]
323 [ + + ]: 972 : tetra_shifts(:, 4, 1,2) = [ 0, 1, 1]
324 [ + + ]: 972 : tetra_shifts(:, 1, 2,2) = [ 0, 0, 0]
325 [ + + ]: 972 : tetra_shifts(:, 2, 2,2) = [ 1, 0, 0]
326 [ + + ]: 972 : tetra_shifts(:, 3, 2,2) = [ 0, 0, 1]
327 [ + + ]: 972 : tetra_shifts(:, 4, 2,2) = [ 0, 1, 1]
328 [ + + ]: 972 : tetra_shifts(:, 1, 3,2) = [ 0, 0, 0]
329 [ + + ]: 972 : tetra_shifts(:, 2, 3,2) = [ -1, 1, 0]
330 [ + + ]: 972 : tetra_shifts(:, 3, 3,2) = [ -1, 1, 1]
331 [ + + ]: 972 : tetra_shifts(:, 4, 3,2) = [ -1, 0, 0]
332 [ + + ]: 972 : tetra_shifts(:, 1, 4,2) = [ 0, 0, 0]
333 [ + + ]: 972 : tetra_shifts(:, 2, 4,2) = [ -1, 0, 1]
334 [ + + ]: 972 : tetra_shifts(:, 3, 4,2) = [ -1, 1, 1]
335 [ + + ]: 972 : tetra_shifts(:, 4, 4,2) = [ -1, 0, 0]
336 [ + + ]: 972 : tetra_shifts(:, 1, 5,2) = [ 0, 0, 0]
337 [ + + ]: 972 : tetra_shifts(:, 2, 5,2) = [ -1, 1, 0]
338 [ + + ]: 972 : tetra_shifts(:, 3, 5,2) = [ 0, 1, 0]
339 [ + + ]: 972 : tetra_shifts(:, 4, 5,2) = [ -1, 1, 1]
340 [ + + ]: 972 : tetra_shifts(:, 1, 6,2) = [ 0, 0, 0]
341 [ + + ]: 972 : tetra_shifts(:, 2, 6,2) = [ 0, 1, 0]
342 [ + + ]: 972 : tetra_shifts(:, 3, 6,2) = [ -1, 1, 1]
343 [ + + ]: 972 : tetra_shifts(:, 4, 6,2) = [ 0, 1, 1]
344 [ + + ]: 972 : tetra_shifts(:, 1, 7,2) = [ 0, 0, 0]
345 [ + + ]: 972 : tetra_shifts(:, 2, 7,2) = [ -1, 0, 1]
346 [ + + ]: 972 : tetra_shifts(:, 3, 7,2) = [ 0, 0, 1]
347 [ + + ]: 972 : tetra_shifts(:, 4, 7,2) = [ -1, 1, 1]
348 [ + + ]: 972 : tetra_shifts(:, 1, 8,2) = [ 0, 0, 0]
349 [ + + ]: 972 : tetra_shifts(:, 2, 8,2) = [ 0, 0, 1]
350 [ + + ]: 972 : tetra_shifts(:, 3, 8,2) = [ -1, 1, 1]
351 [ + + ]: 972 : tetra_shifts(:, 4, 8,2) = [ 0, 1, 1]
352 [ + + ]: 972 : tetra_shifts(:, 1, 9,2) = [ 0, 0, 0]
353 [ + + ]: 972 : tetra_shifts(:, 2, 9,2) = [ 0, 0, 1]
354 [ + + ]: 972 : tetra_shifts(:, 3, 9,2) = [ 0, -1, 0]
355 [ + + ]: 972 : tetra_shifts(:, 4, 9,2) = [ 1, -1, 0]
356 [ + + ]: 972 : tetra_shifts(:, 1,10,2) = [ 0, 0, 0]
357 [ + + ]: 972 : tetra_shifts(:, 2,10,2) = [ 1, 0, 0]
358 [ + + ]: 972 : tetra_shifts(:, 3,10,2) = [ 0, 0, 1]
359 [ + + ]: 972 : tetra_shifts(:, 4,10,2) = [ 1, -1, 0]
360 [ + + ]: 972 : tetra_shifts(:, 1,11,2) = [ 0, 0, 0]
361 [ + + ]: 972 : tetra_shifts(:, 2,11,2) = [ -1, 0, 1]
362 [ + + ]: 972 : tetra_shifts(:, 3,11,2) = [ 0, -1, 0]
363 [ + + ]: 972 : tetra_shifts(:, 4,11,2) = [ -1, 0, 0]
364 [ + + ]: 972 : tetra_shifts(:, 1,12,2) = [ 0, 0, 0]
365 [ + + ]: 972 : tetra_shifts(:, 2,12,2) = [ -1, 0, 1]
366 [ + + ]: 972 : tetra_shifts(:, 3,12,2) = [ 0, 0, 1]
367 [ + + ]: 972 : tetra_shifts(:, 4,12,2) = [ 0, -1, 0]
368 [ + + ]: 972 : tetra_shifts(:, 1,13,2) = [ 0, 0, 0]
369 [ + + ]: 972 : tetra_shifts(:, 2,13,2) = [ 0, 1, 0]
370 [ + + ]: 972 : tetra_shifts(:, 3,13,2) = [ 0, 0, -1]
371 [ + + ]: 972 : tetra_shifts(:, 4,13,2) = [ 1, 0, -1]
372 [ + + ]: 972 : tetra_shifts(:, 1,14,2) = [ 0, 0, 0]
373 [ + + ]: 972 : tetra_shifts(:, 2,14,2) = [ 1, 0, 0]
374 [ + + ]: 972 : tetra_shifts(:, 3,14,2) = [ 0, 1, 0]
375 [ + + ]: 972 : tetra_shifts(:, 4,14,2) = [ 1, 0, -1]
376 [ + + ]: 972 : tetra_shifts(:, 1,15,2) = [ 0, 0, 0]
377 [ + + ]: 972 : tetra_shifts(:, 2,15,2) = [ -1, 1, 0]
378 [ + + ]: 972 : tetra_shifts(:, 3,15,2) = [ 0, 0, -1]
379 [ + + ]: 972 : tetra_shifts(:, 4,15,2) = [ -1, 0, 0]
380 [ + + ]: 972 : tetra_shifts(:, 1,16,2) = [ 0, 0, 0]
381 [ + + ]: 972 : tetra_shifts(:, 2,16,2) = [ -1, 1, 0]
382 [ + + ]: 972 : tetra_shifts(:, 3,16,2) = [ 0, 1, 0]
383 [ + + ]: 972 : tetra_shifts(:, 4,16,2) = [ 0, 0, -1]
384 [ + + ]: 972 : tetra_shifts(:, 1,17,2) = [ 0, 0, 0]
385 [ + + ]: 972 : tetra_shifts(:, 2,17,2) = [ 0, -1, -1]
386 [ + + ]: 972 : tetra_shifts(:, 3,17,2) = [ 1, -1, -1]
387 [ + + ]: 972 : tetra_shifts(:, 4,17,2) = [ 0, 0, -1]
388 [ + + ]: 972 : tetra_shifts(:, 1,18,2) = [ 0, 0, 0]
389 [ + + ]: 972 : tetra_shifts(:, 2,18,2) = [ 0, -1, -1]
390 [ + + ]: 972 : tetra_shifts(:, 3,18,2) = [ 1, -1, -1]
391 [ + + ]: 972 : tetra_shifts(:, 4,18,2) = [ 0, -1, 0]
392 [ + + ]: 972 : tetra_shifts(:, 1,19,2) = [ 0, 0, 0]
393 [ + + ]: 972 : tetra_shifts(:, 2,19,2) = [ 1, -1, -1]
394 [ + + ]: 972 : tetra_shifts(:, 3,19,2) = [ 0, 0, -1]
395 [ + + ]: 972 : tetra_shifts(:, 4,19,2) = [ 1, 0, -1]
396 [ + + ]: 972 : tetra_shifts(:, 1,20,2) = [ 0, 0, 0]
397 [ + + ]: 972 : tetra_shifts(:, 2,20,2) = [ 1, 0, 0]
398 [ + + ]: 972 : tetra_shifts(:, 3,20,2) = [ 1, -1, -1]
399 [ + + ]: 972 : tetra_shifts(:, 4,20,2) = [ 1, 0, -1]
400 [ + + ]: 972 : tetra_shifts(:, 1,21,2) = [ 0, 0, 0]
401 [ + + ]: 972 : tetra_shifts(:, 2,21,2) = [ 1, -1, -1]
402 [ + + ]: 972 : tetra_shifts(:, 3,21,2) = [ 0, -1, 0]
403 [ + + ]: 972 : tetra_shifts(:, 4,21,2) = [ 1, -1, 0]
404 [ + + ]: 972 : tetra_shifts(:, 1,22,2) = [ 0, 0, 0]
405 [ + + ]: 972 : tetra_shifts(:, 2,22,2) = [ 1, 0, 0]
406 [ + + ]: 972 : tetra_shifts(:, 3,22,2) = [ 1, -1, -1]
407 [ + + ]: 972 : tetra_shifts(:, 4,22,2) = [ 1, -1, 0]
408 [ + + ]: 972 : tetra_shifts(:, 1,23,2) = [ 0, 0, 0]
409 [ + + ]: 972 : tetra_shifts(:, 2,23,2) = [ 0, -1, -1]
410 [ + + ]: 972 : tetra_shifts(:, 3,23,2) = [ 0, 0, -1]
411 [ + + ]: 972 : tetra_shifts(:, 4,23,2) = [ -1, 0, 0]
412 [ + + ]: 972 : tetra_shifts(:, 1,24,2) = [ 0, 0, 0]
413 [ + + ]: 972 : tetra_shifts(:, 2,24,2) = [ 0, -1, -1]
414 [ + + ]: 972 : tetra_shifts(:, 3,24,2) = [ 0, -1, 0]
415 [ + + ]: 972 : tetra_shifts(:, 4,24,2) = [ -1, 0, 0]
416 [ + + ]: 972 : tetra_shifts(:, 1, 1,3) = [ 0, 0, 0]
417 [ + + ]: 972 : tetra_shifts(:, 2, 1,3) = [ 1, 0, 0]
418 [ + + ]: 972 : tetra_shifts(:, 3, 1,3) = [ 0, 1, 0]
419 [ + + ]: 972 : tetra_shifts(:, 4, 1,3) = [ 1, 0, 1]
420 [ + + ]: 972 : tetra_shifts(:, 1, 2,3) = [ 0, 0, 0]
421 [ + + ]: 972 : tetra_shifts(:, 2, 2,3) = [ 0, 1, 0]
422 [ + + ]: 972 : tetra_shifts(:, 3, 2,3) = [ 0, 0, 1]
423 [ + + ]: 972 : tetra_shifts(:, 4, 2,3) = [ 1, 0, 1]
424 [ + + ]: 972 : tetra_shifts(:, 1, 3,3) = [ 0, 0, 0]
425 [ + + ]: 972 : tetra_shifts(:, 2, 3,3) = [ -1, 1, 0]
426 [ + + ]: 972 : tetra_shifts(:, 3, 3,3) = [ 0, 0, 1]
427 [ + + ]: 972 : tetra_shifts(:, 4, 3,3) = [ -1, 0, 0]
428 [ + + ]: 972 : tetra_shifts(:, 1, 4,3) = [ 0, 0, 0]
429 [ + + ]: 972 : tetra_shifts(:, 2, 4,3) = [ -1, 1, 0]
430 [ + + ]: 972 : tetra_shifts(:, 3, 4,3) = [ 0, 1, 0]
431 [ + + ]: 972 : tetra_shifts(:, 4, 4,3) = [ 0, 0, 1]
432 [ + + ]: 972 : tetra_shifts(:, 1, 5,3) = [ 0, 0, 0]
433 [ + + ]: 972 : tetra_shifts(:, 2, 5,3) = [ 1, -1, 1]
434 [ + + ]: 972 : tetra_shifts(:, 3, 5,3) = [ 0, -1, 0]
435 [ + + ]: 972 : tetra_shifts(:, 4, 5,3) = [ 1, -1, 0]
436 [ + + ]: 972 : tetra_shifts(:, 1, 6,3) = [ 0, 0, 0]
437 [ + + ]: 972 : tetra_shifts(:, 2, 6,3) = [ 0, -1, 1]
438 [ + + ]: 972 : tetra_shifts(:, 3, 6,3) = [ 1, -1, 1]
439 [ + + ]: 972 : tetra_shifts(:, 4, 6,3) = [ 0, -1, 0]
440 [ + + ]: 972 : tetra_shifts(:, 1, 7,3) = [ 0, 0, 0]
441 [ + + ]: 972 : tetra_shifts(:, 2, 7,3) = [ 1, 0, 0]
442 [ + + ]: 972 : tetra_shifts(:, 3, 7,3) = [ 1, -1, 1]
443 [ + + ]: 972 : tetra_shifts(:, 4, 7,3) = [ 1, -1, 0]
444 [ + + ]: 972 : tetra_shifts(:, 1, 8,3) = [ 0, 0, 0]
445 [ + + ]: 972 : tetra_shifts(:, 2, 8,3) = [ 1, 0, 0]
446 [ + + ]: 972 : tetra_shifts(:, 3, 8,3) = [ 1, -1, 1]
447 [ + + ]: 972 : tetra_shifts(:, 4, 8,3) = [ 1, 0, 1]
448 [ + + ]: 972 : tetra_shifts(:, 1, 9,3) = [ 0, 0, 0]
449 [ + + ]: 972 : tetra_shifts(:, 2, 9,3) = [ 0, -1, 1]
450 [ + + ]: 972 : tetra_shifts(:, 3, 9,3) = [ 1, -1, 1]
451 [ + + ]: 972 : tetra_shifts(:, 4, 9,3) = [ 0, 0, 1]
452 [ + + ]: 972 : tetra_shifts(:, 1,10,3) = [ 0, 0, 0]
453 [ + + ]: 972 : tetra_shifts(:, 2,10,3) = [ 1, -1, 1]
454 [ + + ]: 972 : tetra_shifts(:, 3,10,3) = [ 0, 0, 1]
455 [ + + ]: 972 : tetra_shifts(:, 4,10,3) = [ 1, 0, 1]
456 [ + + ]: 972 : tetra_shifts(:, 1,11,3) = [ 0, 0, 0]
457 [ + + ]: 972 : tetra_shifts(:, 2,11,3) = [ 0, -1, 1]
458 [ + + ]: 972 : tetra_shifts(:, 3,11,3) = [ 0, -1, 0]
459 [ + + ]: 972 : tetra_shifts(:, 4,11,3) = [ -1, 0, 0]
460 [ + + ]: 972 : tetra_shifts(:, 1,12,3) = [ 0, 0, 0]
461 [ + + ]: 972 : tetra_shifts(:, 2,12,3) = [ 0, -1, 1]
462 [ + + ]: 972 : tetra_shifts(:, 3,12,3) = [ 0, 0, 1]
463 [ + + ]: 972 : tetra_shifts(:, 4,12,3) = [ -1, 0, 0]
464 [ + + ]: 972 : tetra_shifts(:, 1,13,3) = [ 0, 0, 0]
465 [ + + ]: 972 : tetra_shifts(:, 2,13,3) = [ 1, 0, 0]
466 [ + + ]: 972 : tetra_shifts(:, 3,13,3) = [ 0, 0, -1]
467 [ + + ]: 972 : tetra_shifts(:, 4,13,3) = [ 0, 1, -1]
468 [ + + ]: 972 : tetra_shifts(:, 1,14,3) = [ 0, 0, 0]
469 [ + + ]: 972 : tetra_shifts(:, 2,14,3) = [ 1, 0, 0]
470 [ + + ]: 972 : tetra_shifts(:, 3,14,3) = [ 0, 1, 0]
471 [ + + ]: 972 : tetra_shifts(:, 4,14,3) = [ 0, 1, -1]
472 [ + + ]: 972 : tetra_shifts(:, 1,15,3) = [ 0, 0, 0]
473 [ + + ]: 972 : tetra_shifts(:, 2,15,3) = [ -1, 0, -1]
474 [ + + ]: 972 : tetra_shifts(:, 3,15,3) = [ 0, 0, -1]
475 [ + + ]: 972 : tetra_shifts(:, 4,15,3) = [ -1, 1, -1]
476 [ + + ]: 972 : tetra_shifts(:, 1,16,3) = [ 0, 0, 0]
477 [ + + ]: 972 : tetra_shifts(:, 2,16,3) = [ -1, 0, -1]
478 [ + + ]: 972 : tetra_shifts(:, 3,16,3) = [ -1, 1, -1]
479 [ + + ]: 972 : tetra_shifts(:, 4,16,3) = [ -1, 0, 0]
480 [ + + ]: 972 : tetra_shifts(:, 1,17,3) = [ 0, 0, 0]
481 [ + + ]: 972 : tetra_shifts(:, 2,17,3) = [ 0, 0, -1]
482 [ + + ]: 972 : tetra_shifts(:, 3,17,3) = [ -1, 1, -1]
483 [ + + ]: 972 : tetra_shifts(:, 4,17,3) = [ 0, 1, -1]
484 [ + + ]: 972 : tetra_shifts(:, 1,18,3) = [ 0, 0, 0]
485 [ + + ]: 972 : tetra_shifts(:, 2,18,3) = [ 0, 1, 0]
486 [ + + ]: 972 : tetra_shifts(:, 3,18,3) = [ -1, 1, -1]
487 [ + + ]: 972 : tetra_shifts(:, 4,18,3) = [ 0, 1, -1]
488 [ + + ]: 972 : tetra_shifts(:, 1,19,3) = [ 0, 0, 0]
489 [ + + ]: 972 : tetra_shifts(:, 2,19,3) = [ -1, 1, 0]
490 [ + + ]: 972 : tetra_shifts(:, 3,19,3) = [ -1, 1, -1]
491 [ + + ]: 972 : tetra_shifts(:, 4,19,3) = [ -1, 0, 0]
492 [ + + ]: 972 : tetra_shifts(:, 1,20,3) = [ 0, 0, 0]
493 [ + + ]: 972 : tetra_shifts(:, 2,20,3) = [ -1, 1, 0]
494 [ + + ]: 972 : tetra_shifts(:, 3,20,3) = [ 0, 1, 0]
495 [ + + ]: 972 : tetra_shifts(:, 4,20,3) = [ -1, 1, -1]
496 [ + + ]: 972 : tetra_shifts(:, 1,21,3) = [ 0, 0, 0]
497 [ + + ]: 972 : tetra_shifts(:, 2,21,3) = [ 0, 0, -1]
498 [ + + ]: 972 : tetra_shifts(:, 3,21,3) = [ 0, -1, 0]
499 [ + + ]: 972 : tetra_shifts(:, 4,21,3) = [ 1, -1, 0]
500 [ + + ]: 972 : tetra_shifts(:, 1,22,3) = [ 0, 0, 0]
501 [ + + ]: 972 : tetra_shifts(:, 2,22,3) = [ 1, 0, 0]
502 [ + + ]: 972 : tetra_shifts(:, 3,22,3) = [ 0, 0, -1]
503 [ + + ]: 972 : tetra_shifts(:, 4,22,3) = [ 1, -1, 0]
504 [ + + ]: 972 : tetra_shifts(:, 1,23,3) = [ 0, 0, 0]
505 [ + + ]: 972 : tetra_shifts(:, 2,23,3) = [ -1, 0, -1]
506 [ + + ]: 972 : tetra_shifts(:, 3,23,3) = [ 0, 0, -1]
507 [ + + ]: 972 : tetra_shifts(:, 4,23,3) = [ 0, -1, 0]
508 [ + + ]: 972 : tetra_shifts(:, 1,24,3) = [ 0, 0, 0]
509 [ + + ]: 972 : tetra_shifts(:, 2,24,3) = [ -1, 0, -1]
510 [ + + ]: 972 : tetra_shifts(:, 3,24,3) = [ 0, -1, 0]
511 [ + + ]: 972 : tetra_shifts(:, 4,24,3) = [ -1, 0, 0]
512 [ + + ]: 972 : tetra_shifts(:, 1, 1,4) = [ 0, 0, 0]
513 [ + + ]: 972 : tetra_shifts(:, 2, 1,4) = [ 1, 0, 0]
514 [ + + ]: 972 : tetra_shifts(:, 3, 1,4) = [ 1, 1, 0]
515 [ + + ]: 972 : tetra_shifts(:, 4, 1,4) = [ 0, 0, 1]
516 [ + + ]: 972 : tetra_shifts(:, 1, 2,4) = [ 0, 0, 0]
517 [ + + ]: 972 : tetra_shifts(:, 2, 2,4) = [ 0, 1, 0]
518 [ + + ]: 972 : tetra_shifts(:, 3, 2,4) = [ 1, 1, 0]
519 [ + + ]: 972 : tetra_shifts(:, 4, 2,4) = [ 0, 0, 1]
520 [ + + ]: 972 : tetra_shifts(:, 1, 3,4) = [ 0, 0, 0]
521 [ + + ]: 972 : tetra_shifts(:, 2, 3,4) = [ 0, 1, 0]
522 [ + + ]: 972 : tetra_shifts(:, 3, 3,4) = [ -1, 0, 1]
523 [ + + ]: 972 : tetra_shifts(:, 4, 3,4) = [ -1, 0, 0]
524 [ + + ]: 972 : tetra_shifts(:, 1, 4,4) = [ 0, 0, 0]
525 [ + + ]: 972 : tetra_shifts(:, 2, 4,4) = [ 0, 1, 0]
526 [ + + ]: 972 : tetra_shifts(:, 3, 4,4) = [ -1, 0, 1]
527 [ + + ]: 972 : tetra_shifts(:, 4, 4,4) = [ 0, 0, 1]
528 [ + + ]: 972 : tetra_shifts(:, 1, 5,4) = [ 0, 0, 0]
529 [ + + ]: 972 : tetra_shifts(:, 2, 5,4) = [ 1, 0, 0]
530 [ + + ]: 972 : tetra_shifts(:, 3, 5,4) = [ 0, -1, 1]
531 [ + + ]: 972 : tetra_shifts(:, 4, 5,4) = [ 0, -1, 0]
532 [ + + ]: 972 : tetra_shifts(:, 1, 6,4) = [ 0, 0, 0]
533 [ + + ]: 972 : tetra_shifts(:, 2, 6,4) = [ 1, 0, 0]
534 [ + + ]: 972 : tetra_shifts(:, 3, 6,4) = [ 0, -1, 1]
535 [ + + ]: 972 : tetra_shifts(:, 4, 6,4) = [ 0, 0, 1]
536 [ + + ]: 972 : tetra_shifts(:, 1, 7,4) = [ 0, 0, 0]
537 [ + + ]: 972 : tetra_shifts(:, 2, 7,4) = [ -1, -1, 1]
538 [ + + ]: 972 : tetra_shifts(:, 3, 7,4) = [ -1, -1, 0]
539 [ + + ]: 972 : tetra_shifts(:, 4, 7,4) = [ 0, -1, 0]
540 [ + + ]: 972 : tetra_shifts(:, 1, 8,4) = [ 0, 0, 0]
541 [ + + ]: 972 : tetra_shifts(:, 2, 8,4) = [ -1, -1, 1]
542 [ + + ]: 972 : tetra_shifts(:, 3, 8,4) = [ -1, -1, 0]
543 [ + + ]: 972 : tetra_shifts(:, 4, 8,4) = [ -1, 0, 0]
544 [ + + ]: 972 : tetra_shifts(:, 1, 9,4) = [ 0, 0, 0]
545 [ + + ]: 972 : tetra_shifts(:, 2, 9,4) = [ -1, -1, 1]
546 [ + + ]: 972 : tetra_shifts(:, 3, 9,4) = [ 0, -1, 1]
547 [ + + ]: 972 : tetra_shifts(:, 4, 9,4) = [ 0, -1, 0]
548 [ + + ]: 972 : tetra_shifts(:, 1,10,4) = [ 0, 0, 0]
549 [ + + ]: 972 : tetra_shifts(:, 2,10,4) = [ -1, -1, 1]
550 [ + + ]: 972 : tetra_shifts(:, 3,10,4) = [ -1, 0, 1]
551 [ + + ]: 972 : tetra_shifts(:, 4,10,4) = [ -1, 0, 0]
552 [ + + ]: 972 : tetra_shifts(:, 1,11,4) = [ 0, 0, 0]
553 [ + + ]: 972 : tetra_shifts(:, 2,11,4) = [ -1, -1, 1]
554 [ + + ]: 972 : tetra_shifts(:, 3,11,4) = [ 0, -1, 1]
555 [ + + ]: 972 : tetra_shifts(:, 4,11,4) = [ 0, 0, 1]
556 [ + + ]: 972 : tetra_shifts(:, 1,12,4) = [ 0, 0, 0]
557 [ + + ]: 972 : tetra_shifts(:, 2,12,4) = [ -1, -1, 1]
558 [ + + ]: 972 : tetra_shifts(:, 3,12,4) = [ -1, 0, 1]
559 [ + + ]: 972 : tetra_shifts(:, 4,12,4) = [ 0, 0, 1]
560 [ + + ]: 972 : tetra_shifts(:, 1,13,4) = [ 0, 0, 0]
561 [ + + ]: 972 : tetra_shifts(:, 2,13,4) = [ 0, 0, -1]
562 [ + + ]: 972 : tetra_shifts(:, 3,13,4) = [ 1, 0, -1]
563 [ + + ]: 972 : tetra_shifts(:, 4,13,4) = [ 1, 1, -1]
564 [ + + ]: 972 : tetra_shifts(:, 1,14,4) = [ 0, 0, 0]
565 [ + + ]: 972 : tetra_shifts(:, 2,14,4) = [ 0, 0, -1]
566 [ + + ]: 972 : tetra_shifts(:, 3,14,4) = [ 0, 1, -1]
567 [ + + ]: 972 : tetra_shifts(:, 4,14,4) = [ 1, 1, -1]
568 [ + + ]: 972 : tetra_shifts(:, 1,15,4) = [ 0, 0, 0]
569 [ + + ]: 972 : tetra_shifts(:, 2,15,4) = [ 1, 0, 0]
570 [ + + ]: 972 : tetra_shifts(:, 3,15,4) = [ 1, 0, -1]
571 [ + + ]: 972 : tetra_shifts(:, 4,15,4) = [ 1, 1, -1]
572 [ + + ]: 972 : tetra_shifts(:, 1,16,4) = [ 0, 0, 0]
573 [ + + ]: 972 : tetra_shifts(:, 2,16,4) = [ 0, 1, 0]
574 [ + + ]: 972 : tetra_shifts(:, 3,16,4) = [ 0, 1, -1]
575 [ + + ]: 972 : tetra_shifts(:, 4,16,4) = [ 1, 1, -1]
576 [ + + ]: 972 : tetra_shifts(:, 1,17,4) = [ 0, 0, 0]
577 [ + + ]: 972 : tetra_shifts(:, 2,17,4) = [ 1, 0, 0]
578 [ + + ]: 972 : tetra_shifts(:, 3,17,4) = [ 1, 1, 0]
579 [ + + ]: 972 : tetra_shifts(:, 4,17,4) = [ 1, 1, -1]
580 [ + + ]: 972 : tetra_shifts(:, 1,18,4) = [ 0, 0, 0]
581 [ + + ]: 972 : tetra_shifts(:, 2,18,4) = [ 0, 1, 0]
582 [ + + ]: 972 : tetra_shifts(:, 3,18,4) = [ 1, 1, 0]
583 [ + + ]: 972 : tetra_shifts(:, 4,18,4) = [ 1, 1, -1]
584 [ + + ]: 972 : tetra_shifts(:, 1,19,4) = [ 0, 0, 0]
585 [ + + ]: 972 : tetra_shifts(:, 2,19,4) = [ 0, 0, -1]
586 [ + + ]: 972 : tetra_shifts(:, 3,19,4) = [ 0, 1, -1]
587 [ + + ]: 972 : tetra_shifts(:, 4,19,4) = [ -1, 0, 0]
588 [ + + ]: 972 : tetra_shifts(:, 1,20,4) = [ 0, 0, 0]
589 [ + + ]: 972 : tetra_shifts(:, 2,20,4) = [ 0, 1, 0]
590 [ + + ]: 972 : tetra_shifts(:, 3,20,4) = [ 0, 1, -1]
591 [ + + ]: 972 : tetra_shifts(:, 4,20,4) = [ -1, 0, 0]
592 [ + + ]: 972 : tetra_shifts(:, 1,21,4) = [ 0, 0, 0]
593 [ + + ]: 972 : tetra_shifts(:, 2,21,4) = [ 0, 0, -1]
594 [ + + ]: 972 : tetra_shifts(:, 3,21,4) = [ 1, 0, -1]
595 [ + + ]: 972 : tetra_shifts(:, 4,21,4) = [ 0, -1, 0]
596 [ + + ]: 972 : tetra_shifts(:, 1,22,4) = [ 0, 0, 0]
597 [ + + ]: 972 : tetra_shifts(:, 2,22,4) = [ 1, 0, 0]
598 [ + + ]: 972 : tetra_shifts(:, 3,22,4) = [ 1, 0, -1]
599 [ + + ]: 972 : tetra_shifts(:, 4,22,4) = [ 0, -1, 0]
600 [ + + ]: 972 : tetra_shifts(:, 1,23,4) = [ 0, 0, 0]
601 [ + + ]: 972 : tetra_shifts(:, 2,23,4) = [ 0, 0, -1]
602 [ + + ]: 972 : tetra_shifts(:, 3,23,4) = [ -1, -1, 0]
603 [ + + ]: 972 : tetra_shifts(:, 4,23,4) = [ 0, -1, 0]
604 [ + + ]: 972 : tetra_shifts(:, 1,24,4) = [ 0, 0, 0]
605 [ + + ]: 972 : tetra_shifts(:, 2,24,4) = [ 0, 0, -1]
606 [ + + ]: 972 : tetra_shifts(:, 3,24,4) = [ -1, -1, 0]
607 [ + + ]: 972 : tetra_shifts(:, 4,24,4) = [ -1, 0, 0]
608 : :
609 : : ! These shifts are taken from previous tetrahedron implementation by MJV and BXU
610 : : ! TODO: implement shifts for the other diagonals
611 [ + + ]: 972 : tetra_shifts_6(:,1,1,1) = [0,0,0]
612 [ + + ]: 972 : tetra_shifts_6(:,2,1,1) = [1,0,0]
613 [ + + ]: 972 : tetra_shifts_6(:,3,1,1) = [0,1,0]
614 [ + + ]: 972 : tetra_shifts_6(:,4,1,1) = [1,0,1]
615 [ + + ]: 972 : tetra_shifts_6(:,1,2,1) = [1,0,0]
616 [ + + ]: 972 : tetra_shifts_6(:,2,2,1) = [1,1,0]
617 [ + + ]: 972 : tetra_shifts_6(:,3,2,1) = [0,1,0]
618 [ + + ]: 972 : tetra_shifts_6(:,4,2,1) = [1,0,1]
619 [ + + ]: 972 : tetra_shifts_6(:,1,3,1) = [0,1,0]
620 [ + + ]: 972 : tetra_shifts_6(:,2,3,1) = [1,1,0]
621 [ + + ]: 972 : tetra_shifts_6(:,3,3,1) = [1,0,1]
622 [ + + ]: 972 : tetra_shifts_6(:,4,3,1) = [1,1,1]
623 [ + + ]: 972 : tetra_shifts_6(:,1,4,1) = [0,0,0]
624 [ + + ]: 972 : tetra_shifts_6(:,2,4,1) = [0,1,0]
625 [ + + ]: 972 : tetra_shifts_6(:,3,4,1) = [0,0,1]
626 [ + + ]: 972 : tetra_shifts_6(:,4,4,1) = [1,0,1]
627 [ + + ]: 972 : tetra_shifts_6(:,1,5,1) = [0,0,1]
628 [ + + ]: 972 : tetra_shifts_6(:,2,5,1) = [1,0,1]
629 [ + + ]: 972 : tetra_shifts_6(:,3,5,1) = [0,1,0]
630 [ + + ]: 972 : tetra_shifts_6(:,4,5,1) = [0,1,1]
631 [ + + ]: 972 : tetra_shifts_6(:,1,6,1) = [0,1,0]
632 [ + + ]: 972 : tetra_shifts_6(:,2,6,1) = [1,0,1]
633 [ + + ]: 972 : tetra_shifts_6(:,3,6,1) = [0,1,1]
634 [ + + ]: 972 : tetra_shifts_6(:,4,6,1) = [1,1,1]
635 : :
636 [ + + ]: 972 : main_diagonals(:,1) = [ 1, 1, 1] ! 0-7
637 [ + + ]: 972 : main_diagonals(:,2) = [-1, 1, 1] ! 1-6
638 [ + + ]: 972 : main_diagonals(:,3) = [ 1,-1, 1] ! 2-5
639 [ + + ]: 972 : main_diagonals(:,4) = [ 1, 1,-1] ! 3-4
640 : :
641 : 243 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
642 : 243 : tetra%nkibz = nkpt_ibz
643 : 243 : tetra%nkbz = nkpt_fullbz
644 : : ! HM: this value should be important for performance
645 : : ! more buckets means faster queries for unique tetrahedra
646 : : ! but more memory due to the initial size TETRA_SIZE of the buckets
647 : : ! when changing the number of buckets one should also change the hash function
648 : : ! to distribute the tetrahedra in the buckets as uniformly as possible
649 : : ! the simplest hash (not the best!) is:
650 : : ! ihash = mod(sum(tetra_ibz),nbuckts)
651 : : ! the value of sum(tetra_ibz) is between 1 and 4*nkibz so I use nkibz nbuckets
652 : : ! a larger number of buckets should speed up finding the irreducible tetrahedra
653 : : ! but with more memory allocated
654 : 243 : tetra%nbuckets = nkpt_ibz
655 : 243 : ierr = 0
656 [ - + + - : 225265 : ABI_CALLOC(tetra_hash_count,(tetra%nbuckets))
+ - + - -
+ - + +
+ ]
657 : :
658 : : ! Determine the smallest diagonal in k-space
659 : : min_length = huge(min_length)
660 [ + + ]: 1215 : do idiag = 1,4
661 : : diag(:) = gprimd(:,1)*main_diagonals(1,idiag)+&
662 : : gprimd(:,2)*main_diagonals(2,idiag)+&
663 [ + + ]: 3888 : gprimd(:,3)*main_diagonals(3,idiag)
664 : 972 : length = sqrt(diag(1)*diag(1) + diag(2)*diag(2) + diag(3)*diag(3))
665 [ + + ]: 1215 : if (length < min_length) then
666 : 246 : min_length = length
667 : 246 : min_idiag = idiag
668 : : end if
669 : : end do
670 : :
671 : : ! HM TODO: Avoid krank and map the k-point grid to indexes
672 : : ! Make full k-point rank arrays
673 : : !oct = octree_init(kpt_fullbz,2**4,[-one,-one,-one],[two,two,two])
674 : 243 : call krank%init(nkpt_fullbz, kpt_fullbz)
675 : : !
676 : : ! HM (13/04/2019): I implement two different versions:
677 : : ! 1. I only use 24 tetrahedra around the IBZ k-point
678 : : ! following the approach of A. Togo (phonopy, spglib, kspclib).
679 : : ! 2. I generate tetrahedra on the full Brillouin zone
680 : : ! and keep track of how many are contributing to the IBZ point and the multiplicities,
681 : : ! these can be more than 24 tetrahedra.
682 : : ! This is equivalent to what Matthieu implemented but avoids large memory allocations.
683 : : !
684 : : ! The two implementations differ specially when using low k-point sampling
685 : : ! (the second yields the same results using IBZ or FBZ, the first one not).
686 : : ! For large sampling the two approaches yield similar results, with the first
687 : : ! one using less memory, faster to generate and compute
688 : : !
689 [ - + + - : 225265 : ABI_MALLOC(tetra%unique_tetra,(tetra%nbuckets))
+ - + - -
+ - + - +
+ + ]
690 [ + + ]: 224779 : do ihash=1,tetra%nbuckets
691 [ - + - + : 8308075 : ABI_CALLOC(tetra%unique_tetra(ihash)%indexes,(0:4,TETRA_SIZE))
+ + + + ]
692 : : end do
693 [ - + ]: 243 : tetra%opt = 2; if (present(opt)) tetra%opt = opt
694 : :
695 : 243 : select case(tetra%opt)
696 : : case (1)
697 : : ! For each k-point in the IBZ store 24 tetrahedra each referring to 4 k-points
698 [ # # ]: 0 : do ikibz=1,tetra%nkibz
699 : : !if (mod(ikibz,nprocs) /= my_rank) cycle
700 [ # # ]: 0 : k1 = kpt_ibz(:,ikibz)
701 [ # # ]: 0 : tetra_loop1: do itetra=1,24
702 [ # # ]: 0 : do isummit=1,4
703 : : ! Find the index of the neighbouring k-points in the BZ
704 : : k2 = k1 + tetra_shifts(1,isummit,itetra,min_idiag)*klatt(:,1) + &
705 : : tetra_shifts(2,isummit,itetra,min_idiag)*klatt(:,2) + &
706 [ # # ]: 0 : tetra_shifts(3,isummit,itetra,min_idiag)*klatt(:,3)
707 : : ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full !
708 : : !ikpt2 = octree_find(oct,k2,dist)
709 : : !ikpt2 = octree_find_nearest_pbc(oct,k2,dist,shift)
710 : : !if (dist>tol12) call exit(1)
711 : 0 : ikpt2 = krank%get_index(k2)
712 : : ! Find the index of those points in the BZ and IBZ
713 : 0 : tetra_ibz(isummit) = bz2ibz(ikpt2)
714 : : end do
715 : : ! Sort index of irr k-point edges (need this so the comparison works)
716 : 0 : call sort_4tetra_int(tetra_ibz)
717 : :
718 : : ! Store only unique tetrahedra
719 : : ! Compute a very simple hash for each tetrahedron
720 : 0 : ihash = compute_hash(tetra,tetra_ibz) !mod(sum(tetra_ibz),tetra%nbuckets)+1
721 : : ! Loop over all tetrahedrons that contain this ikibz as first element
722 [ # # ]: 0 : do jtetra=1,tetra_hash_count(ihash)
723 : : ! if tetrahedron already exists add multiplicity
724 [ # # ]: 0 : if (tetra%unique_tetra(ihash)%indexes(1,jtetra)/=tetra_ibz(1)) cycle
725 [ # # ]: 0 : if (tetra%unique_tetra(ihash)%indexes(2,jtetra)/=tetra_ibz(2)) cycle
726 [ # # ]: 0 : if (tetra%unique_tetra(ihash)%indexes(3,jtetra)/=tetra_ibz(3)) cycle
727 [ # # ]: 0 : if (tetra%unique_tetra(ihash)%indexes(4,jtetra)/=tetra_ibz(4)) cycle
728 : 0 : tetra%unique_tetra(ihash)%indexes(0,jtetra) = tetra%unique_tetra(ihash)%indexes(0,jtetra)+1
729 : 0 : cycle tetra_loop1
730 : : end do
731 : : ! Otherwise store new tetrahedron
732 : 0 : tetra_hash_count(ihash) = tetra_hash_count(ihash)+1
733 : 0 : max_ntetra = size(tetra%unique_tetra(ihash)%indexes,2)
734 : : ! The contents don't fit the array so I have to resize it
735 [ # # ]: 0 : if (tetra_hash_count(ihash)>max_ntetra) then
736 [ # # # # : 0 : ABI_MALLOC(indexes,(0:4,max_ntetra+TETRA_STEP))
# # # # #
# # # #
# ]
737 [ # # # # ]: 0 : indexes(0:4,:max_ntetra) = tetra%unique_tetra(ihash)%indexes
738 [ # # # # ]: 0 : indexes(:,max_ntetra+1:) = 0
739 [ # # ]: 0 : ABI_MOVE_ALLOC(indexes,tetra%unique_tetra(ihash)%indexes)
740 : : end if
741 [ # # ]: 0 : tetra%unique_tetra(ihash)%indexes(1:,tetra_hash_count(ihash)) = tetra_ibz(:)
742 : 0 : tetra%unique_tetra(ihash)%indexes(0, tetra_hash_count(ihash)) = 1
743 : : end do tetra_loop1
744 : : end do
745 : : ! HM: The multiplicity of the tetrahedrons computed so far is wrong because we are using IBZ
746 : : ! I compute the k-point multiplicity so I can fix this later on
747 : : ! Only needed in blochl_weights* interface for looping over tetrahedra.
748 : : ! in the onewk routines the weight is known outside
749 [ # # # # : 0 : ABI_CALLOC(tetra%ibz_multiplicity,(tetra%nkibz))
# # # # #
# # # # #
# # ]
750 [ # # ]: 0 : do ikbz=1,nkpt_fullbz
751 : 0 : ikibz = bz2ibz(ikbz)
752 : 0 : tetra%ibz_multiplicity(ikibz) = tetra%ibz_multiplicity(ikibz) + 1
753 : : end do
754 : :
755 : : case(2)
756 : 243 : min_idiag = 1
757 : : ! For each k-point in the BZ generate the 6 tetrahedra that tessellate a microzone
758 [ + + ]: 1409584 : do ikbz=1,tetra%nkbz
759 [ + + ]: 5637364 : k1 = kpt_fullbz(:,ikbz)
760 [ + + ]: 9865630 : tetra_loop2: do itetra=1,6
761 : : ! Determine tetrahedron
762 [ + + ]: 42280230 : do isummit=1,4
763 : : ! Find the index of the neighbouring k-points in the BZ
764 : : k2 = k1 + tetra_shifts_6(1,isummit,itetra,min_idiag)*klatt(:,1) + &
765 : : tetra_shifts_6(2,isummit,itetra,min_idiag)*klatt(:,2) + &
766 [ + + ]: 135296736 : tetra_shifts_6(3,isummit,itetra,min_idiag)*klatt(:,3)
767 : : ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full !
768 : : !ikpt2 = octree_find(oct,k2,dist)
769 : : !ikpt2 = octree_find_nearest_pbc(oct,k2,dist,shift)
770 : : !if (dist>tol12) call exit(1)
771 : 33824184 : ikpt2 = krank%get_index(k2)
772 : : ! Find the index of those points in the BZ and IBZ
773 : 42280230 : tetra_ibz(isummit) = bz2ibz(ikpt2)
774 : : end do
775 : : ! Sort index of irr k-point edges (need this so the comparison works)
776 : 8456046 : call sort_4tetra_int(tetra_ibz)
777 : :
778 : : ! Store only unique tetrahedra
779 : : ! Compute a very simple hash for each tetrahedron
780 : 8456046 : ihash = compute_hash(tetra,tetra_ibz) !mod(sum(tetra_ibz),tetra%nbuckets)+1
781 : : ! Loop over all tetrahedrons that contain this ikibz as first element
782 [ + + ]: 156269209 : do jtetra=1,tetra_hash_count(ihash)
783 : : ! if tetrahedron already exists add multiplicity
784 [ + + ]: 151530110 : if (tetra%unique_tetra(ihash)%indexes(1,jtetra)/=tetra_ibz(1)) cycle
785 [ + + ]: 4805680 : if (tetra%unique_tetra(ihash)%indexes(2,jtetra)/=tetra_ibz(2)) cycle
786 [ + + ]: 3930585 : if (tetra%unique_tetra(ihash)%indexes(3,jtetra)/=tetra_ibz(3)) cycle
787 [ - + ]: 3716947 : if (tetra%unique_tetra(ihash)%indexes(4,jtetra)/=tetra_ibz(4)) cycle
788 : 3716947 : tetra%unique_tetra(ihash)%indexes(0,jtetra) = tetra%unique_tetra(ihash)%indexes(0,jtetra)+1
789 : 156269209 : cycle tetra_loop2
790 : : end do
791 : : ! Otherwise store new tetrahedron
792 : 4739099 : tetra_hash_count(ihash) = tetra_hash_count(ihash)+1
793 : 4739099 : max_ntetra = size(tetra%unique_tetra(ihash)%indexes,2)
794 : : ! The contents don't fit the array so I have to resize it
795 [ + + ]: 4739099 : if (tetra_hash_count(ihash)>max_ntetra) then
796 [ + - - + : 1975152 : ABI_MALLOC(indexes,(0:4,max_ntetra+TETRA_STEP))
+ - + - +
- - + -
+ ]
797 [ + + + + ]: 71822824 : indexes(0:4,:max_ntetra) = tetra%unique_tetra(ihash)%indexes
798 [ + + + + ]: 24360208 : indexes(:,max_ntetra+1:) = 0
799 [ + - ]: 658384 : ABI_MOVE_ALLOC(indexes,tetra%unique_tetra(ihash)%indexes)
800 : : end if
801 [ + + ]: 23695495 : tetra%unique_tetra(ihash)%indexes(1:,tetra_hash_count(ihash)) = tetra_ibz(:)
802 : 6148440 : tetra%unique_tetra(ihash)%indexes(0, tetra_hash_count(ihash)) = 1
803 : : end do tetra_loop2
804 : : end do
805 : : case default
806 : 0 : ierr = 1
807 : 0 : write(errorstring,*) 'Invalid option for the generation of tetrahedra,',ch10,&
808 : 0 : 'possible options are:',ch10,&
809 : 0 : '1. Generate 24 tetrahedra per k-point',ch10,&
810 : 0 : '2. Generate tetrahedra in the FBZ a map to IBZ (default)'
811 [ - + - ]: 243 : return
812 : : end select
813 : :
814 : : !ierr = octree_free(oct)
815 [ - + ]: 243 : ABI_FREE(tetra_hash_count)
816 : 243 : call krank%free()
817 : :
818 : : ! Do some maintenance: free unused memory and count unique tetrahedra per IBZ point
819 : 243 : tetra%nunique_tetra = 0
820 [ + + ]: 224779 : do ihash=1,tetra%nbuckets
821 : : ! Count tetrahedra in this bucket
822 [ + + + + ]: 5522056 : ntetra = count(tetra%unique_tetra(ihash)%indexes(0,:)>0)
823 : 224536 : tetra%nunique_tetra = tetra%nunique_tetra + ntetra
824 : : ! Allocate array with right size
825 [ + + + - : 673608 : ABI_MALLOC(indexes,(0:4,ntetra))
+ - + + -
+ - + ]
826 [ + - + - : 28883666 : indexes = tetra%unique_tetra(ihash)%indexes(:,:ntetra)
- + - - +
+ + + ]
827 [ + - ]: 224779 : ABI_MOVE_ALLOC(indexes, tetra%unique_tetra(ihash)%indexes)
828 : : end do
829 : :
830 : : ! Sum the multiplicity
831 [ - + + - : 729 : ABI_MALLOC(tetra%tetra_count,(tetra%nkibz))
+ - + - -
+ - + -
+ ]
832 [ - + + - : 729 : ABI_MALLOC(tetra%tetra_total,(tetra%nkibz))
+ - + - -
+ - + -
+ ]
833 [ + + ]: 224779 : tetra%tetra_count = 0
834 [ + + ]: 224779 : tetra%tetra_total = 0
835 [ + + ]: 224779 : do ihash=1,tetra%nbuckets
836 : 224536 : ntetra = size(tetra%unique_tetra(ihash)%indexes,2)
837 [ + + ]: 4963878 : do itetra=1,ntetra
838 [ + + ]: 28434594 : tetra_mibz = tetra%unique_tetra(ihash)%indexes(:,itetra)
839 [ + + ]: 23920031 : do isummit=1,4
840 : 18956396 : ikibz = tetra_mibz(isummit)
841 : 18956396 : tetra%tetra_total(ikibz) = tetra%tetra_total(ikibz) + tetra_mibz(0)
842 : 23695495 : tetra%tetra_count(ikibz) = tetra%tetra_count(ikibz) + 1
843 : : end do
844 : : end do
845 : : end do
846 [ + + ]: 224779 : tetra%nibz_tetra = sum(tetra%tetra_count)
847 : :
848 : : ! HM: This was being allocated here, however this is only used when we loop over kpoints
849 : : ! I will only allocate this memory if the htetra_get_onewk_* routines are called (lazy evaluation)
850 : : !call htetra_init_mapping_ibz(tetra)
851 : :
852 : : ! Calculate the volume of the tetrahedra
853 : : rcvol = abs(gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3))- &
854 : : gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3))+ &
855 : 243 : gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
856 : :
857 : : ! Volume of all tetrahedra should be the same as that of tetra 1
858 : : ! this is the volume of 1 tetrahedron, should be coherent with notation in Lehmann & Taut
859 [ + + ]: 972 : k1(:) = gprimd(:,1)*klatt(1,1) + gprimd(:,2)*klatt(2,1) + gprimd(:,3)*klatt(3,1)
860 [ + + ]: 972 : k2(:) = gprimd(:,1)*klatt(1,2) + gprimd(:,2)*klatt(2,2) + gprimd(:,3)*klatt(3,2)
861 [ + + ]: 972 : k3(:) = gprimd(:,1)*klatt(1,3) + gprimd(:,2)*klatt(2,3) + gprimd(:,3)*klatt(3,3)
862 : : tetra%vv = abs(k1(1)*(k2(2)*k3(3)-k2(3)*k3(2))- &
863 : : k1(2)*(k2(1)*k3(3)-k2(3)*k3(1))+ &
864 [ - - - + : 243 : k1(3)*(k2(1)*k3(2)-k2(2)*k3(1))) / 6.d0 / rcvol
- + - + -
- - - -
- ]
865 : :
866 : : contains
867 : 8456046 : integer function compute_hash(tetra,t) result(ihash)
868 : : class(htetra_t),intent(in) :: tetra
869 : : integer,intent(in) :: t(4)
870 [ - - + + ]: 42280230 : ihash = mod(sum(t), tetra%nbuckets) + 1
871 : : ! TODO: should use a more general hash function that supports more buckets
872 : : ! Something like:
873 : : ! id = t(1)*nk3+t(2)*nk2+t(3)*nk1+t(4)
874 : : ! where nk is the number of points in the IBZ
875 : : ! Computing this leads to overflow so should use
876 : : ! mod computation operations
877 : : ! (A + B) mod C = (A mod C + B mod C) mod C
878 : : ! (A * B) mod C = (A mod C * B mod C) mod C
879 : : ! A^B mod C = ( (A mod C)^B ) mod C
880 : : end function compute_hash
881 : :
882 : : end subroutine htetra_init
883 : : !!***
884 : :
885 : : !----------------------------------------------------------------------
886 : :
887 : : !!****f* m_htetra/htetra_init_mapping_ibz
888 : : !! NAME
889 : : !! htetra_init_mapping_ibz
890 : : !!
891 : : !! FUNCTION
892 : : !! The mapping to the IBZ has its own allocation routine.
893 : : !! I will only allocate this memory if the htetra_get_onewk_* routines are called (lazy evaluation)
894 : : !!
895 : : !! SOURCE
896 : :
897 : 157 : subroutine htetra_init_mapping_ibz(tetra)
898 : :
899 : : class(htetra_t),intent(inout) :: tetra
900 : :
901 : : integer :: ikibz, itetra, isummit, ihash, ntetra
902 : 157 : integer :: tetra_count(tetra%nkibz),tetra_mibz(0:4)
903 : : real(dp) :: mem_mb
904 : :
905 : : ! Only execute the following if not yet allocated
906 [ - + ]: 157 : if (allocated(tetra%ibz)) return
907 : :
908 : : ! Allocate IBZ to tetrahedron mapping
909 [ - + + - : 21562 : ABI_MALLOC(tetra%ibz, (tetra%nkibz))
+ - + - -
+ - + - +
+ + ]
910 [ + + ]: 314 : mem_mb = ABI_MEM_MB(tetra%ibz)
911 [ + + ]: 21248 : do ikibz=1,tetra%nkibz
912 [ - + + - : 63273 : ABI_MALLOC(tetra%ibz(ikibz)%indexes, (2, tetra%tetra_count(ikibz)))
+ - + - -
+ - + -
+ ]
913 : 21248 : mem_mb = mem_mb + 2 * tetra%tetra_count(ikibz) * 4 * b2Mb
914 : : end do
915 : :
916 : 157 : call wrtout(std_out, sjoin(" Allocating tetra%ibz%indexes with memory:", ftoa(mem_mb, fmt="f8.1"), " [Mb] <<< MEM"))
917 : :
918 : : ! Create mapping from IBZ to unique tetrahedra
919 [ + + ]: 21248 : tetra_count = 0
920 [ + + ]: 21248 : do ihash=1,tetra%nbuckets
921 : 21091 : ntetra = size(tetra%unique_tetra(ihash)%indexes, dim=2)
922 [ + + ]: 869866 : do itetra=1,ntetra
923 [ + + ]: 5091708 : tetra_mibz = tetra%unique_tetra(ihash)%indexes(:,itetra)
924 [ + + ]: 4264181 : do isummit=1,4
925 : 3394472 : ikibz = tetra_mibz(isummit)
926 : 3394472 : tetra_count(ikibz) = tetra_count(ikibz) + 1
927 : 3394472 : tetra%ibz(ikibz)%indexes(1, tetra_count(ikibz)) = ihash
928 : 4243090 : tetra%ibz(ikibz)%indexes(2, tetra_count(ikibz)) = itetra
929 : : end do
930 : : end do
931 : : end do
932 : :
933 : : end subroutine htetra_init_mapping_ibz
934 : : !!***
935 : :
936 : : !----------------------------------------------------------------------
937 : :
938 : : !!****f* m_htetra/htetra_get_ibz
939 : : !! NAME
940 : : !! htetra_get_ibz
941 : : !!
942 : : !! FUNCTION
943 : : !! Get the itetra tetrahedron contributing to the ikibz k-point
944 : : !!
945 : : !! SOURCE
946 : :
947 : 15710160 : pure subroutine htetra_get_ibz(tetra, ikibz, itetra, tetra_mibz)
948 : :
949 : : class(htetra_t), intent(in) :: tetra
950 : : integer,intent(in) :: ikibz, itetra
951 : : integer,intent(out) :: tetra_mibz(0:4)
952 : : integer :: ihash, jtetra
953 : :
954 : 15710160 : ihash = tetra%ibz(ikibz)%indexes(1,itetra)
955 : 15710160 : jtetra = tetra%ibz(ikibz)%indexes(2,itetra)
956 [ + + ]: 94260960 : tetra_mibz = tetra%unique_tetra(ihash)%indexes(:,jtetra)
957 : :
958 : 15710160 : end subroutine htetra_get_ibz
959 : : !!***
960 : :
961 : : !----------------------------------------------------------------------
962 : :
963 : : !!****f* m_htetra/htetra_print
964 : : !! NAME
965 : : !! htetra_print
966 : : !!
967 : : !! FUNCTION
968 : : !! write information about the tetrahedra object
969 : : !!
970 : : !! SOURCE
971 : :
972 : 62 : subroutine htetra_print(self, unit)
973 : :
974 : : class(htetra_t), intent(in) :: self
975 : : integer,intent(in) :: unit
976 : :
977 : : real(dp) :: total_size, unique_tetra_size, ibz_pointer_size
978 : :
979 [ - + ]: 62 : if (unit == dev_null) return
980 : :
981 : 62 : unique_tetra_size = self%nunique_tetra * 5* four / 1024 ** 2
982 : 62 : total_size = unique_tetra_size
983 : : !write(unit,'(a,i0)') ' htetra unique_tetra:', self%nunique_tetra
984 : : !write(unit,'(a,f12.1,a)') ' htetra unique_tetra_size ', unique_tetra_size, ' [Mb] <<< MEM'
985 [ - + ]: 62 : if (allocated(self%ibz)) then
986 : 0 : ibz_pointer_size = self%nibz_tetra*2*four / 1024 ** 2
987 : 0 : write(unit,'(a,i0)') ' ibz_tetra: ', self%nibz_tetra
988 : 0 : write(unit,'(a,f12.1,a)') ' ibz_tetra_size: ', ibz_pointer_size, ' [Mb] <<< MEM'
989 : 0 : total_size = total_size + ibz_pointer_size
990 : : end if
991 : : ! integer arrays
992 : 62 : total_size = total_size + 3 * (self%nkibz * 4) / 1024**2
993 : 62 : write(unit,'(a,f12.1,a)') ' htetra total size:', total_size, ' [Mb] <<< MEM'
994 : :
995 : : end subroutine htetra_print
996 : : !!***
997 : :
998 : : !----------------------------------------------------------------------
999 : :
1000 : : !!****f* m_htetra/htetra_free
1001 : : !! NAME
1002 : : !! htetra_free
1003 : : !!
1004 : : !! FUNCTION
1005 : : !! deallocate tetrahedra pointers if needed
1006 : : !!
1007 : : !! SOURCE
1008 : :
1009 : 374 : subroutine htetra_free(tetra)
1010 : :
1011 : : class(htetra_t), intent(inout) :: tetra
1012 : : integer :: ikibz,ihash
1013 : :
1014 [ + + ]: 374 : ABI_SFREE(tetra%tetra_count)
1015 [ + + ]: 374 : ABI_SFREE(tetra%tetra_total)
1016 [ - + ]: 374 : ABI_SFREE(tetra%ibz_multiplicity)
1017 : :
1018 [ + + ]: 374 : if (allocated(tetra%unique_tetra)) then
1019 [ + + ]: 224779 : do ihash=1,tetra%nbuckets
1020 [ + - ]: 224779 : ABI_SFREE(tetra%unique_tetra(ihash)%indexes)
1021 : : end do
1022 [ + - + + : 224779 : ABI_FREE(tetra%unique_tetra)
- + - + ]
1023 : : end if
1024 : :
1025 [ + + ]: 374 : if (allocated(tetra%ibz)) then
1026 [ + + ]: 21248 : do ikibz=1,tetra%nkibz
1027 [ + - ]: 21248 : ABI_SFREE(tetra%ibz(ikibz)%indexes)
1028 : : end do
1029 [ + - + + : 21248 : ABI_FREE(tetra%ibz)
- + - + ]
1030 : : end if
1031 : :
1032 : 374 : end subroutine htetra_free
1033 : : !!***
1034 : :
1035 : : !----------------------------------------------------------------------
1036 : :
1037 : : !!****f* m_htetra/get_onetetra_blochl
1038 : : !! NAME
1039 : : !! get_onetetra_blochl
1040 : : !!
1041 : : !! FUNCTION
1042 : : !! Private function to calculate the contributions to the weights due to a single tetrahedron.
1043 : : !! Extracted from get_tetra_weight
1044 : : !!
1045 : : !! SOURCE
1046 : :
1047 : 138582926 : pure subroutine get_onetetra_blochl(eig, energies, nene, bcorr, tweight, dweight)
1048 : :
1049 : : !Arguments ------------------------------------
1050 : : !scalars
1051 : : integer,intent(in) :: nene,bcorr
1052 : : !arrays
1053 : : real(dp),intent(in) :: eig(4), energies(nene)
1054 : : real(dp), intent(out) :: tweight(4, nene), dweight(4, nene)
1055 : :
1056 : : !Local variables-------------------------------
1057 : : integer :: ieps
1058 : : real(dp) :: cc,cc1,cc2,cc3
1059 : : real(dp) :: dcc1de,dcc2de,dcc3de,dccde,eps
1060 : : real(dp) :: e21,e31,e32,e41,e42,e43
1061 : : real(dp) :: inv_e32,inv_e41,inv_e42,inv_e43,inv_e21,inv_e31
1062 : : real(dp) :: deleps1,deleps2,deleps3,deleps4
1063 : : real(dp) :: e1,e2,e3,e4
1064 : : real(dp) :: invepsum, cc_pre, dccde_pre
1065 : : real(dp) :: cc1_pre, cc2_pre, cc3_pre
1066 : : real(dp) :: dccde_tmp
1067 : : real(dp) :: bcorr_fact
1068 : : ! *********************************************************************
1069 : :
1070 : : ! This is output
1071 [ + + + + : >13491*10^7 : tweight = zero; dweight = zero
+ + + + ]
1072 : :
1073 : : ! all notations are from Blochl PRB 49 16223 [[cite:Bloechl1994a]] Appendix B
1074 : 138582926 : e1 = eig(1)
1075 : 138582926 : e2 = eig(2)
1076 : 138582926 : e3 = eig(3)
1077 : 138582926 : e4 = eig(4)
1078 : 138582926 : e21 = e2-e1
1079 : 138582926 : e31 = e3-e1
1080 : 138582926 : e41 = e4-e1
1081 : 138582926 : e32 = e3-e2
1082 : 138582926 : e42 = e4-e2
1083 : 138582926 : e43 = e4-e3
1084 [ + + ]: 138582926 : inv_e21 = zero; if (e21 > tol14) inv_e21 = 1.d0 / e21
1085 [ + + ]: 138582926 : inv_e31 = zero; if (e31 > tol14) inv_e31 = 1.d0 / e31
1086 [ + + ]: 138582926 : inv_e41 = zero; if (e41 > tol14) inv_e41 = 1.d0 / e41
1087 [ + + ]: 138582926 : inv_e32 = zero; if (e32 > tol14) inv_e32 = 1.d0 / e32
1088 [ + + ]: 138582926 : inv_e42 = zero; if (e42 > tol14) inv_e42 = 1.d0 / e42
1089 [ + + ]: 138582926 : inv_e43 = zero; if (e43 > tol14) inv_e43 = 1.d0 / e43
1090 : :
1091 [ + + ]: 7878293144 : do ieps=1,nene
1092 : 7770985705 : eps = energies(ieps)
1093 : :
1094 : : !
1095 : : ! eps < e1 nothing to do
1096 : : !
1097 [ + + ]: 7770985705 : if (eps < e1) cycle
1098 : :
1099 : : !
1100 : : ! e1 < eps < e2
1101 : : !
1102 [ + + ]: 510470751 : if (eps < e2) then
1103 : 172325526 : deleps1 = eps-e1
1104 : 172325526 : invepsum = inv_e21+inv_e31+inv_e41
1105 : :
1106 : : ! Heaviside
1107 : 172325526 : cc = inv_e21*inv_e31*inv_e41*deleps1**3
1108 : 172325526 : tweight(1,ieps) = cc*(4.d0-deleps1*invepsum)
1109 : 172325526 : tweight(2,ieps) = cc*deleps1*inv_e21
1110 : 172325526 : tweight(3,ieps) = cc*deleps1*inv_e31
1111 : 172325526 : tweight(4,ieps) = cc*deleps1*inv_e41
1112 : :
1113 : : ! Delta
1114 : 172325526 : dccde_pre = 3.d0*inv_e21*inv_e31*inv_e41
1115 : 172325526 : dccde = dccde_pre*deleps1**2
1116 : 172325526 : dweight(1,ieps) = dccde*(4.d0-deleps1*invepsum)-cc*invepsum
1117 : 172325526 : dweight(2,ieps) = (dccde*deleps1+cc) * inv_e21
1118 : 172325526 : dweight(3,ieps) = (dccde*deleps1+cc) * inv_e31
1119 : 172325526 : dweight(4,ieps) = (dccde*deleps1+cc) * inv_e41
1120 : :
1121 [ - + ]: 172325526 : if (bcorr == 1) then
1122 : : ! bxu, correction terms based on Bloechl's paper
1123 : 0 : bcorr_fact = 4.d0/40.d0*dccde_pre*deleps1*deleps1
1124 : 0 : tweight(1,ieps) = tweight(1,ieps) + bcorr_fact*( e21+e31+e41)
1125 : 0 : tweight(2,ieps) = tweight(2,ieps) + bcorr_fact*(-e21+e32+e42)
1126 : 0 : tweight(3,ieps) = tweight(3,ieps) + bcorr_fact*(-e31-e32+e43)
1127 : 0 : tweight(4,ieps) = tweight(4,ieps) + bcorr_fact*(-e41-e42-e43)
1128 : :
1129 : 0 : bcorr_fact = 8.d0/40.d0*dccde_pre*deleps1
1130 : 0 : dweight(1,ieps) = dweight(1,ieps) + bcorr_fact*( e21+e31+e41)
1131 : 0 : dweight(2,ieps) = dweight(2,ieps) + bcorr_fact*(-e21+e32+e42)
1132 : 0 : dweight(3,ieps) = dweight(3,ieps) + bcorr_fact*(-e31-e32+e43)
1133 : 0 : dweight(4,ieps) = dweight(4,ieps) + bcorr_fact*(-e41-e42-e43)
1134 : : end if
1135 : :
1136 : : cycle
1137 : : endif
1138 : :
1139 : : !
1140 : : ! e2 < eps < e3
1141 : : !
1142 [ + + ]: 338145225 : if (eps < e3) then
1143 : 144551565 : deleps1 = eps-e1
1144 : 144551565 : deleps2 = eps-e2
1145 : 144551565 : deleps3 = e3-eps
1146 : 144551565 : deleps4 = e4-eps
1147 : 144551565 : cc1_pre = inv_e31*inv_e41
1148 : 144551565 : cc2_pre = inv_e41*inv_e32*inv_e31
1149 : 144551565 : cc3_pre = inv_e42*inv_e32*inv_e41
1150 : :
1151 : : ! Heaviside
1152 : 144551565 : cc1 = cc1_pre*deleps1*deleps1
1153 : 144551565 : cc2 = cc2_pre*deleps1*deleps2*deleps3
1154 : 144551565 : cc3 = cc3_pre*deleps2*deleps2*deleps4
1155 : :
1156 : : tweight(1,ieps) = (cc1)+&
1157 : : (cc1+cc2)*deleps3*inv_e31+&
1158 : 144551565 : (cc1+cc2+cc3)*deleps4*inv_e41
1159 : : tweight(2,ieps) = (cc1+cc2+cc3)+&
1160 : : (cc2+cc3)*deleps3*inv_e32+&
1161 : 144551565 : (cc3)*deleps4*inv_e42
1162 : : tweight(3,ieps) = (cc1+cc2)*deleps1*inv_e31+&
1163 : 144551565 : (cc2+cc3)*deleps2*inv_e32
1164 : : tweight(4,ieps) = (cc1+cc2+cc3)*deleps1*inv_e41+&
1165 : 144551565 : (cc3)*deleps2*inv_e42
1166 : :
1167 : : ! Delta
1168 : 144551565 : dcc1de = cc1_pre*(2.d0*deleps1)
1169 : 144551565 : dcc2de = cc2_pre*( -deleps1*deleps2+deleps1*deleps3+deleps2*deleps3)
1170 : 144551565 : dcc3de = cc3_pre*(2.d0*deleps2*deleps4-deleps2*deleps2)
1171 : : dweight(1,ieps) = dcc1de+&
1172 : : ((dcc1de+dcc2de)*deleps3-(cc1+cc2))*inv_e31+&
1173 : 144551565 : ((dcc1de+dcc2de+dcc3de)*deleps4-(cc1+cc2+cc3))*inv_e41
1174 : : dweight(2,ieps) = (dcc1de+dcc2de+dcc3de)+&
1175 : : ((dcc2de+dcc3de)*deleps3-(cc2+cc3))*inv_e32+&
1176 : 144551565 : (dcc3de*deleps4-cc3)*inv_e42
1177 : : dweight(3,ieps) = ((dcc1de+dcc2de)*deleps1+(cc1+cc2))*inv_e31+&
1178 : 144551565 : ((dcc2de+dcc3de)*deleps2+(cc2+cc3))*inv_e32
1179 : : dweight(4,ieps) = ((dcc1de+dcc2de+dcc3de)*deleps1+(cc1+cc2+cc3))*inv_e41+&
1180 : 144551565 : (dcc3de*deleps2+cc3)*inv_e42
1181 : :
1182 [ - + ]: 144551565 : if (bcorr == 1) then
1183 : : ! bxu, correction terms based on Bloechl's paper
1184 : : ! The correction terms may cause the dweight become negative
1185 : 0 : bcorr_fact = 4.d0/40.d0*cc1_pre*(3.d0*e21+6.d0*deleps2-3.d0*(e31+e42)*deleps2*deleps2*inv_e32*inv_e42)
1186 : 0 : tweight(1,ieps) = tweight(1,ieps) + bcorr_fact*( e21+e31+e41)
1187 : 0 : tweight(2,ieps) = tweight(2,ieps) + bcorr_fact*(-e21+e32+e42)
1188 : 0 : tweight(3,ieps) = tweight(3,ieps) + bcorr_fact*(-e31-e32+e43)
1189 : 0 : tweight(4,ieps) = tweight(4,ieps) + bcorr_fact*(-e41-e42-e43)
1190 : :
1191 : 0 : bcorr_fact = 4.d0/40.d0*cc1_pre*(6.d0-6.d0*(e31+e42)*deleps2*inv_e32*inv_e42)
1192 : 0 : dweight(1,ieps) = dweight(1,ieps) + bcorr_fact*( e21+e31+e41)
1193 : 0 : dweight(2,ieps) = dweight(2,ieps) + bcorr_fact*(-e21+e32+e42)
1194 : 0 : dweight(3,ieps) = dweight(3,ieps) + bcorr_fact*(-e31-e32+e43)
1195 : 0 : dweight(4,ieps) = dweight(4,ieps) + bcorr_fact*(-e41-e42-e43)
1196 : : end if
1197 : :
1198 : : cycle
1199 : : endif
1200 : :
1201 : : !
1202 : : ! e3 < eps < e4
1203 : : !
1204 [ + + ]: 193593660 : if (eps < e4) then
1205 : 162310083 : deleps4 = e4-eps
1206 : 162310083 : invepsum = inv_e41+inv_e42+inv_e43
1207 : :
1208 : : ! Heaviside
1209 : 162310083 : cc_pre = inv_e41*inv_e42*inv_e43
1210 : 162310083 : cc = cc_pre*deleps4**3
1211 : 162310083 : tweight(1,ieps) = one - deleps4*cc*inv_e41
1212 : 162310083 : tweight(2,ieps) = one - deleps4*cc*inv_e42
1213 : 162310083 : tweight(3,ieps) = one - deleps4*cc*inv_e43
1214 : 162310083 : tweight(4,ieps) = one - cc*(4.d0-deleps4*invepsum)
1215 : :
1216 : : ! Delta
1217 : 162310083 : dccde = -3.d0*cc_pre*deleps4**2
1218 : 162310083 : dccde_tmp = dccde*deleps4 - cc
1219 : 162310083 : dweight(1,ieps) = -dccde_tmp * inv_e41
1220 : 162310083 : dweight(2,ieps) = -dccde_tmp * inv_e42
1221 : 162310083 : dweight(3,ieps) = -dccde_tmp * inv_e43
1222 : 162310083 : dweight(4,ieps) = -4.d0*dccde + dccde_tmp*invepsum
1223 : :
1224 [ - + ]: 162310083 : if (bcorr == 1) then
1225 : : ! bxu, correction terms based on Bloechl's paper
1226 : : ! The correction terms may cause the dweight become negative
1227 : 0 : bcorr_fact = 12.d0/40.d0*cc_pre*deleps4*deleps4
1228 : 0 : tweight(1,ieps) = tweight(1,ieps) + bcorr_fact*( e21+e31+e41)
1229 : 0 : tweight(2,ieps) = tweight(2,ieps) + bcorr_fact*(-e21+e32+e42)
1230 : 0 : tweight(3,ieps) = tweight(3,ieps) + bcorr_fact*(-e31-e32+e43)
1231 : 0 : tweight(4,ieps) = tweight(4,ieps) + bcorr_fact*(-e41-e42-e43)
1232 : :
1233 : 0 : bcorr_fact = - 24.d0/40.d0*cc_pre*deleps4
1234 : 0 : dweight(1,ieps) = dweight(1,ieps) + bcorr_fact*( e21+e31+e41)
1235 : 0 : dweight(2,ieps) = dweight(2,ieps) + bcorr_fact*(-e21+e32+e42)
1236 : 0 : dweight(3,ieps) = dweight(3,ieps) + bcorr_fact*(-e31-e32+e43)
1237 : 0 : dweight(4,ieps) = dweight(4,ieps) + bcorr_fact*(-e41-e42-e43)
1238 : : end if
1239 : :
1240 : : cycle
1241 : : endif
1242 : :
1243 : : !
1244 : : ! e4 < eps
1245 : : !
1246 [ + + ]: 138591016 : if (e4 < eps) then
1247 : :
1248 : : ! Heaviside
1249 [ + + + + ]: 28719032867 : tweight(:,ieps:) = one
1250 : :
1251 : : ! Delta unchanged by this tetrahedron
1252 : : exit
1253 : : end if
1254 : :
1255 : : ! if we have a fully degenerate tetrahedron,
1256 : : ! 1) the tweight is a Heaviside (step) function, which is correct above, but
1257 : : ! 2) the dweight should contain a Dirac function
1258 : : !
1259 : : end do
1260 : :
1261 : 138582926 : end subroutine get_onetetra_blochl
1262 : : !!***
1263 : :
1264 : : !!****f* m_htetra/get_onetetra_lambinvigneron
1265 : : !! NAME
1266 : : !! get_onetetra_lambinvigneron
1267 : : !!
1268 : : !! FUNCTION
1269 : : !! Compute the complex weights according to: P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984)
1270 : : !! This routine is adapted from tdep where it was implemented
1271 : : !! by Olle Hellman, all credits go to him
1272 : : !!
1273 : : !! INPUTS
1274 : : !!
1275 : : !! OUTPUT
1276 : : !!
1277 : : !! SOURCE
1278 : :
1279 : : !pure
1280 : 0 : subroutine get_onetetra_lambinvigneron(eig, z, cw)
1281 : :
1282 : : ! dispersion values at the corners of the tetrahedron
1283 : : real(dp), intent(in) :: eig(4)
1284 : : ! energy to evaluate the weights at
1285 : : complex(dp), intent(in) :: z
1286 : : ! complex weights
1287 : : complex(dp), intent(out) :: cw(4)
1288 : : complex(dp) :: EZ1,EZ2,EZ3,EZ4
1289 : : real(dp) :: tol = tol14
1290 : : !real(dp) :: tol = tol10
1291 : : !real(dp) :: tol = tol6
1292 : : real(dp) :: Emin,Emax,Zdist
1293 : : real(dp) :: E12,E13,E14,E23,E24,E34
1294 : : real(dp) :: a,b,c,d,e,f
1295 : : complex(dp) zmE(4), verli(4) !, verm(4)
1296 : : !integer :: ii, jj
1297 : :
1298 : 0 : cw = zero
1299 : :
1300 : : ! Min and max energy
1301 : 0 : Emin=eig(1)
1302 : 0 : Emax=eig(4)
1303 : :
1304 : : ! First the complex energy differences
1305 : 0 : zmE = z - eig(:)
1306 : 0 : EZ1=z-eig(1)
1307 : 0 : EZ2=z-eig(2)
1308 : 0 : EZ3=z-eig(3)
1309 : 0 : EZ4=z-eig(4)
1310 : : ! Smallest distance |z-Ei|, to determine whether I should switch to the
1311 : : ! asymptotic behavior for numerical stability.
1312 : : Zdist=huge(Zdist)
1313 : : Zdist=min(Zdist,abs(EZ2))
1314 : : Zdist=min(Zdist,abs(EZ3))
1315 : : Zdist=min(Zdist,abs(EZ4))
1316 : : !@TODO add asymptotic thing with continued fractions
1317 : :
1318 : : ! Then the energy differences, for the coefficients. Must always be positive, I hope.
1319 : 0 : E12=eig(2)-eig(1)
1320 : 0 : E13=eig(3)-eig(1)
1321 : 0 : E14=eig(4)-eig(1)
1322 : 0 : E23=eig(3)-eig(2)
1323 : 0 : E24=eig(4)-eig(2)
1324 : 0 : E34=eig(4)-eig(3)
1325 : 0 : a=zero; if ( E12 .gt. tol ) a=one/E12
1326 : 0 : b=zero; if ( E13 .gt. tol ) b=one/E13
1327 : 0 : c=zero; if ( E14 .gt. tol ) c=one/E14
1328 : 0 : d=zero; if ( E23 .gt. tol ) d=one/E23
1329 : 0 : e=zero; if ( E24 .gt. tol ) e=one/E24
1330 : 0 : f=zero; if ( E34 .gt. tol ) f=one/E34
1331 : :
1332 : : ! Now get the actual weights
1333 : : ! e1=e2=e3=e4
1334 [ # # ]: 0 : if ( E12+E23+E34 .lt. tol ) then
1335 : : #if 0
1336 : : !print *, "e1=e2=e3=e4"
1337 : : cw(1)=0.25_dp/EZ1
1338 : : cw(2)=0.25_dp/EZ2
1339 : : cw(3)=0.25_dp/EZ3
1340 : : cw(4)=0.25_dp/EZ4
1341 : : #else
1342 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1343 : : #endif
1344 : : ! e2=e3=e4 ! diff wrt simteta
1345 [ # # ]: 0 : elseif ( E23+E34 .lt. tol ) then
1346 : : #if 0
1347 : : !print *, "e2=e3=e4"
1348 : : cw(1)=-a - (3*a**2*EZ2)*half + 3*a**3*EZ1*EZ2 + 3*a**4*EZ1**2*EZ2*Log(EZ2/EZ1)
1349 : : cw(2)=-a - (3*a**2*EZ2)*half + 3*a**3*EZ1*EZ2 + 3*a**4*EZ1**2*EZ2*Log(EZ2/EZ1)
1350 : : cw(3)=-a - (3*a**2*EZ2)*half + 3*a**3*EZ1*EZ2 + 3*a**4*EZ1**2*EZ2*Log(EZ2/EZ1)
1351 : : cw(4)=-a*third + (a**2*EZ1)*half - a**3*EZ1**2 + a**4*EZ1**3*Log(EZ2/EZ1)
1352 : : #else
1353 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1354 : : #endif
1355 : :
1356 : : ! e1=e2=e3
1357 [ # # ]: 0 : elseif ( E12+E23 .lt. tol ) then
1358 : : #if 0
1359 : : !print *, "e1=e2=e3" ! diff wrt simteta
1360 : : cw(1)=f*third - (EZ4*f**2)*half + EZ4**2*f**3 + EZ4**3*f**4*Log(EZ4/EZ3)
1361 : : cw(2)=f*third - (EZ4*f**2)*half + EZ4**2*f**3 + EZ4**3*f**4*Log(EZ4/EZ3)
1362 : : cw(3)=f*third - (EZ4*f**2)*half + EZ4**2*f**3 + EZ4**3*f**4*Log(EZ4/EZ3)
1363 : : cw(4)=-f + (3*EZ3*f**2)*half - 3*EZ3*EZ4*f**3 + 3*EZ3*EZ4**2*f**4*Log(EZ3/EZ4)
1364 : : #else
1365 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1366 : : #endif
1367 : : ! e1=e2 e3=e4
1368 [ # # ]: 0 : elseif ( E12+E34 .lt. tol ) then
1369 : : #if 0
1370 : : !print *, "e1=e2 e3=e4" ! FIXME This is Buggy, does not work even with parabolic dispersion
1371 : : cw(1)=-d - (3*d**2*EZ2)*half + 3*d**3*EZ2*EZ3 + 3*d**4*EZ2*EZ3**2*Log(EZ2/EZ3)
1372 : : cw(2)=-d - (3*d**2*EZ2)*half + 3*d**3*EZ2*EZ3 + 3*d**4*EZ2*EZ3**2*Log(EZ2/EZ3)
1373 : : cw(3)=d - (3*d**2*EZ3)*half - 3*d**3*EZ2*EZ3 + 3*d**4*EZ2**2*EZ3*Log(EZ3/EZ2)
1374 : : cw(4)=d - (3*d**2*EZ3)*half - 3*d**3*EZ2*EZ3 + 3*d**4*EZ2**2*EZ3*Log(EZ3/EZ2)
1375 : : #else
1376 : : !cw(1) = nine * EZ3**2 * EZ2 / E23**4 * log(EZ2/EZ3) * EZ2 * (EZ3 -E23)/E23**3 - one/E23
1377 : : !cw(2) = cw(1)
1378 : : !cw(3) = nine * EZ2**2 * EZ3 / E23**4 * log(EZ3/EZ2) * EZ3 * (EZ2 +E23)/E23**3 + one/E23
1379 : : !cw(4) = cw(3)
1380 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1381 : : !cw = zero
1382 : : #endif
1383 : :
1384 : : ! e3=e4
1385 [ # # ]: 0 : elseif ( E34 .lt. tol ) then
1386 : : !print *, "e3=e4"
1387 : : #if 0
1388 : : cw(1)=-(a*b**2*EZ1**2*(-1 + (a*EZ2 + 2*b*EZ3)*Log(EZ1))) + a**2*d**2*EZ2**3*Log(EZ2) - &
1389 : : b**2*d*EZ3**2*(1 + (2*b*EZ1 + d*EZ2)*Log(EZ3))
1390 : : cw(2)=a**2*b**2*EZ1**3*Log(EZ1) - a*d**2*EZ2**2*(1 + (a*EZ1 - 2*d*EZ3)*Log(EZ2)) - &
1391 : : b*d**2*EZ3**2*(1 + (b*EZ1 + 2*d*EZ2)*Log(EZ3))
1392 : : cw(3)=a*b**3*EZ1**3*Log(EZ1) - a*d**3*EZ2**3*Log(EZ2) + b*d*EZ3*(half + b*EZ1 + d*EZ2 + &
1393 : : (b**2*EZ1**2 + b*d*EZ1*EZ2 + d**2*EZ2**2)*Log(EZ3))
1394 : : cw(4)=a*b**3*EZ1**3*Log(EZ1) - a*d**3*EZ2**3*Log(EZ2) + b*d*EZ3*(half + b*EZ1 + d*EZ2 + &
1395 : : (b**2*EZ1**2 + b*d*EZ1*EZ2 + d**2*EZ2**2)*Log(EZ3))
1396 : : #else
1397 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1398 : : #endif
1399 : : ! e2=e3
1400 [ # # ]: 0 : elseif ( E23 .lt. tol ) then
1401 : : #if 0
1402 : : !print *, "e2=e3"
1403 : : cw(1)=-(a**2*c*EZ1**2*(-1 + (2*a*EZ2 + c*EZ4)*Log(EZ1))) + &
1404 : : a**2*e*EZ2**2*(1 + (2*a*EZ1 - e*EZ4)*Log(EZ2)) + c**2*e**2*EZ4**3*Log(EZ4)
1405 : : cw(2)=a**3*c*EZ1**3*Log(EZ1) - &
1406 : : a*e*EZ2*(half + a*EZ1 - e*EZ4 + (a**2*EZ1**2 - a*e*EZ1*EZ4 + e**2*EZ4**2)*Log(EZ2)) + c*e**3*EZ4**3*Log(EZ4)
1407 : : cw(3)=a**3*c*EZ1**3*Log(EZ1) - &
1408 : : a*e*EZ2*(half + a*EZ1 - e*EZ4 + (a**2*EZ1**2 - a*e*EZ1*EZ4 + e**2*EZ4**2)*Log(EZ2)) + c*e**3*EZ4**3*Log(EZ4)
1409 : : cw(4)=a**2*c**2*EZ1**3*Log(EZ1) - &
1410 : : a*e**2*EZ2**2*(1 + (a*EZ1 - 2*e*EZ4)*Log(EZ2)) - c*e**2*EZ4**2*(1 + (c*EZ1 + 2*e*EZ2)*Log(EZ4))
1411 : : #else
1412 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1413 : : #endif
1414 : :
1415 : : ! e1=e2
1416 [ # # ]: 0 : elseif ( E12 .lt. tol ) then
1417 : : !print *, "e1=e2"
1418 : : #if 0
1419 : : cw(1)=b*c*EZ1*(half - b*EZ3 - c*EZ4 + (b**2*EZ3**2 + b*c*EZ3*EZ4 + c**2*EZ4**2)*Log(EZ1)) - &
1420 : : b**3*EZ3**3*f*Log(EZ3) + c**3*EZ4**3*f*Log(EZ4)
1421 : : cw(2)=b*c*EZ1*(half - b*EZ3 - c*EZ4 + (b**2*EZ3**2 + b*c*EZ3*EZ4 + c**2*EZ4**2)*Log(EZ1)) - &
1422 : : b**3*EZ3**3*f*Log(EZ3) + c**3*EZ4**3*f*Log(EZ4)
1423 : : cw(3)=-(b**2*c*EZ1**2*(-1 + (2*b*EZ3 + c*EZ4)*Log(EZ1))) + &
1424 : : b**2*EZ3**2*f*(1 + (2*b*EZ1 - EZ4*f)*Log(EZ3)) + c**2*EZ4**3*f**2*Log(EZ4)
1425 : : cw(4)=-(b*c**2*EZ1**2*(-1 + (b*EZ3 + 2*c*EZ4)*Log(EZ1))) + &
1426 : : b**2*EZ3**3*f**2*Log(EZ3) - c**2*EZ4**2*f*(1 + (2*c*EZ1 + EZ3*f)*Log(EZ4))
1427 : : #else
1428 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1429 : : #endif
1430 : : ! e1<e2<e3<e4
1431 : : else
1432 : : !print *, "e1<e2<e3<e4"
1433 : : #if 0
1434 : : cw(1)=a**2*d*e*EZ2**3*Log(EZ2/EZ1) - b**2*d*EZ3**3*f*Log(EZ3/EZ1) + c*(a*b*EZ1**2 + c*e*EZ4**3*f*Log(EZ4/EZ1))
1435 : : cw(2)=a**2*b*c*EZ1**3*Log(EZ1/EZ2) - b*d**2*EZ3**3*f*Log(EZ3/EZ2) + e*(-(a*d*EZ2**2) + c*e*EZ4**3*f*Log(EZ4/EZ2))
1436 : : cw(3)=a*b**2*c*EZ1**3*Log(EZ1/EZ3) - a*d**2*e*EZ2**3*Log(EZ2/EZ3) + f*(b*d*EZ3**2 + c*e*EZ4**3*f*Log(EZ4/EZ3))
1437 : : cw(4)=a*b*c**2*EZ1**3*Log(EZ1/EZ4) - a*d*e**2*EZ2**3*Log(EZ2/EZ4) + f*(-(c*e*EZ4**2) + b*d*EZ3**3*f*Log(EZ3/EZ4))
1438 : : #else
1439 : : !do ii=1,4
1440 : : ! cw(ii) = zme(ii) ** 2 / prod_wo(ii)
1441 : : ! do jj=1,4
1442 : : ! if (jj == ii) cycle
1443 : : ! cw(ii) = cw(ii) + (zme(jj)**3 / prod_wo(jj) * log(zme(jj) / zme(ii)) / (eig(ii) - eig(jj)))
1444 : : ! end do
1445 : : !end do
1446 : :
1447 [ # # ]: 0 : call SIM0TWOI(cw, VERLI, z-eig)
1448 : : #endif
1449 : : endif
1450 : :
1451 : : ! HM:check this
1452 : : !cw = cw * two
1453 : :
1454 : : contains
1455 : : pure real(dp) function prod_wo(ii)
1456 : : integer,intent(in) :: ii
1457 : : integer :: kk
1458 : :
1459 : : prod_wo = one
1460 : : do kk=1,4
1461 : : if (kk == ii) cycle
1462 : : prod_wo = prod_wo * (eig(kk) - eig(ii))
1463 : : end do
1464 : : end function prod_wo
1465 : :
1466 : : end subroutine get_onetetra_lambinvigneron
1467 : : !!***
1468 : :
1469 : : !!****f* m_htetra/get_ontetratra_lambinvigneron_imag
1470 : : !! NAME
1471 : : !! get_ontetratra_lambinvigneron_imag
1472 : : !!
1473 : : !! FUNCTION
1474 : : !! Compute the complex weights according to:
1475 : : !! P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984)
1476 : : !! This routine is adapted from tdep where it was implemented
1477 : : !! by Olle Hellman, all credits go to him
1478 : : !!
1479 : : !! INPUTS
1480 : : !!
1481 : : !! OUTPUT
1482 : : !!
1483 : : !! SOURCE
1484 : :
1485 : 0 : pure subroutine get_onetetetra_lambinvigneron_imag(eig, energies, nene, wt)
1486 : :
1487 : : ! dispersion values at the corners of the tetrahedron
1488 : : real(dp), intent(in), dimension(4) :: eig
1489 : : ! number of energies
1490 : : integer, intent(in) :: nene
1491 : : ! energy to evaluate the weights at
1492 : : real(dp), intent(in) :: energies(nene)
1493 : : ! integration weights
1494 : : real(dp), intent(out) :: wt(4,nene)
1495 : :
1496 : : integer :: ie
1497 : : real(dp) :: z
1498 : : real(dp) :: EZ1,EZ2,EZ3,EZ4
1499 : : real(dp) :: Emin,Emax
1500 : : real(dp) :: E12,E13,E14,E23,E24,E34
1501 : : real(dp) :: a,b,c,d,e,f,ff0,ff1,ff2,ff3,gg0,gg1,gg2,gg3,hh0,hh1,hh2,hh3,ii0,ii1,ii2,ii3
1502 : :
1503 [ # # # # ]: 0 : wt = zero
1504 : : Emin = eig(1)
1505 : : Emax = eig(4)
1506 : :
1507 [ # # ]: 0 : do ie=1,nene
1508 : 0 : z = energies(ie)
1509 [ # # ]: 0 : if (z<eig(1)) then ! e<e1<e2<e3<e4
1510 : : cycle
1511 [ # # ]: 0 : else if (z .lt. eig(2)) then ! e1<e<e2<e3<e4
1512 : 0 : EZ1=z-eig(1)
1513 : 0 : EZ2=z-eig(2)
1514 : 0 : EZ3=z-eig(3)
1515 : 0 : EZ4=z-eig(4)
1516 : 0 : E12=eig(2)-eig(1)
1517 : 0 : E13=eig(3)-eig(1)
1518 : 0 : E14=eig(4)-eig(1)
1519 : 0 : a=one/E12
1520 : 0 : b=one/E13
1521 : 0 : c=one/E14
1522 : 0 : wt(1,ie)=a*b*c*EZ1**2*(-a*EZ2 - b*EZ3 - c*EZ4)
1523 : 0 : wt(2,ie)=a**2*b*c*EZ1**3
1524 : 0 : wt(3,ie)=a*b**2*c*EZ1**3
1525 : 0 : wt(4,ie)=a*b*c**2*EZ1**3
1526 : 0 : cycle
1527 [ # # ]: 0 : else if (z .lt. eig(3)) then ! e1<e2<e<e3<e4
1528 : 0 : EZ1=z-eig(1)
1529 : 0 : EZ2=z-eig(2)
1530 : 0 : EZ3=z-eig(3)
1531 : 0 : EZ4=z-eig(4)
1532 : 0 : E13=eig(3)-eig(1)
1533 : 0 : E14=eig(4)-eig(1)
1534 : 0 : E23=eig(3)-eig(2)
1535 : 0 : E24=eig(4)-eig(2)
1536 : 0 : b=one/E13
1537 : 0 : c=one/E14
1538 : 0 : d=one/E23
1539 : 0 : e=one/E24
1540 : 0 : ff0=-b**2*EZ3
1541 : 0 : ff2=-c**2*EZ4
1542 : 0 : gg0=-d**2*EZ3
1543 : 0 : gg2=-e**2*EZ4
1544 : 0 : hh0=d**2*EZ2
1545 : 0 : hh2=b**2*EZ1
1546 : 0 : ii0=e**2*EZ2
1547 : 0 : ii2=c**2*EZ1
1548 : 0 : ff1=-c*d*EZ1*EZ3-d*e*EZ2*EZ3-c*e*EZ1*EZ4
1549 : 0 : ff3=-b*d*EZ1*EZ3-b*e*EZ1*EZ4-d*e*EZ2*EZ4
1550 : 0 : gg1=-b*c*EZ1*EZ3-b*e*EZ2*EZ3-c*e*EZ2*EZ4
1551 : 0 : gg3=-b*d*EZ2*EZ3-b*c*EZ1*EZ4-c*d*EZ2*EZ4
1552 : 0 : hh1=-b*c*EZ1*EZ3-b*e*EZ2*EZ3-c*e*EZ2*EZ4
1553 : 0 : hh3=-c*d*EZ1*EZ3-d*e*EZ2*EZ3-c*e*EZ1*EZ4
1554 : 0 : ii1=-b*d*EZ2*EZ3-b*c*EZ1*EZ4-c*d*EZ2*EZ4
1555 : 0 : ii3=-b*d*EZ1*EZ3-b*e*EZ1*EZ4-d*e*EZ2*EZ4
1556 : 0 : wt(1,ie)=half*(ff0*ff1+ff2*ff3)
1557 : 0 : wt(2,ie)=half*(gg0*gg1+gg2*gg3)
1558 : 0 : wt(3,ie)=half*(hh0*hh1+hh2*hh3)
1559 : 0 : wt(4,ie)=half*(ii0*ii1+ii2*ii3)
1560 : 0 : cycle
1561 [ # # ]: 0 : else if (z .lt. eig(4)) then ! e1<e2<e3<e<e4
1562 : 0 : EZ1=z-eig(1)
1563 : 0 : EZ2=z-eig(2)
1564 : 0 : EZ3=z-eig(3)
1565 : 0 : EZ4=z-eig(4)
1566 : 0 : E14=eig(4)-eig(1)
1567 : 0 : E24=eig(4)-eig(2)
1568 : 0 : E34=eig(4)-eig(3)
1569 : 0 : c=one/E14
1570 : 0 : e=one/E24
1571 : 0 : f=one/E34
1572 : 0 : wt(1,ie)=-(c**2*e*EZ4**3*f)
1573 : 0 : wt(2,ie)=-(c*e**2*EZ4**3*f)
1574 : 0 : wt(3,ie)=-(c*e*EZ4**3*f**2)
1575 : 0 : wt(4,ie)=c*e*EZ4**2*f*(c*EZ1 + e*EZ2 + EZ3*f)
1576 : 0 : cycle
1577 : : else
1578 : 0 : exit
1579 : : end if
1580 : : end do
1581 [ # # # # ]: 0 : wt = wt*4.0_dp
1582 : :
1583 : 0 : end subroutine get_onetetetra_lambinvigneron_imag
1584 : : !!***
1585 : :
1586 : : !----------------------------------------------------------------------
1587 : :
1588 : : !!****f* m_htetra/htetra_get_onewk_wvals
1589 : : !! NAME
1590 : : !! htetra_get_onewk_wvals
1591 : : !!
1592 : : !! FUNCTION
1593 : : !! Calculate integration weights and their derivatives for a single k-point in the IBZ.
1594 : : !!
1595 : : !! INPUTS
1596 : : !! ik_ibz=Index of the k-point in the IBZ array
1597 : : !! bcorr=1 to include Blochl correction else 0.
1598 : : !! nw=number of energies in wvals
1599 : : !! nibz=number of irreducible kpoints
1600 : : !! wvals(nw)=Frequency points.
1601 : : !! eigen_ibz(nkibz)=eigenenergies for each k point
1602 : : !!
1603 : : !! OUTPUT
1604 : : !! weights(nw,2) = integration weights for
1605 : : !! Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
1606 : : !! for a given (band, k-point, spin).
1607 : : !!
1608 : : !! SOURCE
1609 : :
1610 : 100721 : subroutine htetra_get_onewk_wvals(tetra, ik_ibz, opt, nw, wvals, max_occ, nkibz, eig_ibz, weights)
1611 : :
1612 : : !Arguments ------------------------------------
1613 : : !scalars
1614 : : class(htetra_t), intent(inout) :: tetra
1615 : : integer,intent(in) :: ik_ibz,nw,nkibz,opt
1616 : : real(dp),intent(in) :: max_occ
1617 : : !arrays
1618 : : real(dp),intent(in) :: wvals(nw), eig_ibz(nkibz)
1619 : : real(dp),intent(out) :: weights(nw, 2)
1620 : :
1621 : : !Local variables-------------------------------
1622 : : !scalars
1623 : : integer :: itetra,isummit,tetra_count,tetra_total
1624 : : real(dp) :: tweight
1625 : : !arrays
1626 : : integer :: ind_ibz(4),tetra_mibz(0:4)
1627 : 100721 : real(dp) :: eig(4), tweight_tmp(4,nw),dweight_tmp(4,nw)
1628 : : ! *********************************************************************
1629 : :
1630 [ + + + + ]: 224806779 : weights = zero
1631 : : ! lazy evaluation of the mapping from k-points to tetrahedra
1632 [ + + ]: 100721 : if (.not.allocated(tetra%ibz)) call htetra_init_mapping_ibz(tetra)
1633 : :
1634 : : ! For each tetrahedron that belongs to this k-point
1635 : 100721 : tetra_count = tetra%tetra_count(ik_ibz)
1636 : 100721 : tetra_total = tetra%tetra_total(ik_ibz)
1637 [ + + ]: 15810881 : do itetra=1,tetra_count
1638 : :
1639 : 15710160 : call htetra_get_ibz(tetra, ik_ibz, itetra, tetra_mibz)
1640 : 15710160 : tweight = one*tetra_mibz(0) / tetra_total
1641 [ + + ]: 78550800 : do isummit=1,4
1642 : : ! Get mapping of each summit to eig_ibz
1643 : 62840640 : ind_ibz(isummit) = tetra_mibz(isummit)
1644 : 78550800 : eig(isummit) = eig_ibz(ind_ibz(isummit))
1645 : : end do
1646 : :
1647 : : ! Sort energies before calling get_onetetra_blochl
1648 : 15710160 : call sort_4tetra(eig, ind_ibz)
1649 : :
1650 : : ! HM: Here we should only compute what we will use!
1651 : 31420320 : select case (opt)
1652 : : case (0:1)
1653 : 15710160 : call get_onetetra_blochl(eig, wvals, nw, opt, tweight_tmp, dweight_tmp)
1654 : : case (2)
1655 : 0 : call get_onetetetra_lambinvigneron_imag(eig, wvals, nw, dweight_tmp)
1656 [ + - - - : 15710160 : tweight_tmp = zero
- - - ]
1657 : : end select
1658 : :
1659 : : ! Accumulate contributions to ik_ibz (there might be multiple vertices that map onto ik_ibz)
1660 [ + - ]: 38922106 : do isummit=1,4
1661 [ + + ]: 38821385 : if (ind_ibz(isummit) /= ik_ibz) cycle
1662 [ + + ]: 13286034336 : weights(:,1) = weights(:,1) + dweight_tmp(isummit,:)*tweight*max_occ
1663 [ + + ]: 13286034336 : weights(:,2) = weights(:,2) + tweight_tmp(isummit,:)*tweight*max_occ
1664 : : ! HM: This exit is important, avoids summing the same contribution more than once
1665 : 23111225 : exit
1666 : : end do
1667 : : end do ! itetra
1668 : :
1669 : 100721 : end subroutine htetra_get_onewk_wvals
1670 : : !!***
1671 : :
1672 : : !----------------------------------------------------------------------
1673 : :
1674 : : !!****f* m_htetra/tetra_get_onewk
1675 : : !! NAME
1676 : : !! tetra_get_onewk
1677 : : !!
1678 : : !! FUNCTION
1679 : : !! Calculate integration weights and their derivatives for a single k-point in the IBZ.
1680 : : !! Same as above but different calling arguments.
1681 : : !! IBZ Weights are not included
1682 : : !! HM: The above is preferred but I keep this one to ease the transition
1683 : : !!
1684 : : !! INPUTS
1685 : : !!
1686 : : !! OUTPUT
1687 : : !!
1688 : : !! SOURCE
1689 : :
1690 : 12326 : subroutine htetra_get_onewk(tetra, ik_ibz, bcorr, nw, nkibz, eig_ibz, enemin, enemax, max_occ, weights)
1691 : :
1692 : : !Arguments ------------------------------------
1693 : : !scalars
1694 : : class(htetra_t), intent(inout) :: tetra
1695 : : integer,intent(in) :: ik_ibz,nw,nkibz,bcorr
1696 : : real(dp) ,intent(in) :: enemin,enemax,max_occ
1697 : : !arrays
1698 : : real(dp),intent(in) :: eig_ibz(nkibz)
1699 : : real(dp),intent(out) :: weights(nw,2)
1700 : :
1701 : : !Local variables-------------------------------
1702 : : !scalars
1703 : 12326 : real(dp) :: wvals(nw)
1704 : : ! *********************************************************************
1705 : :
1706 [ + + + + ]: 63510798 : weights = zero
1707 : 12326 : wvals = linspace(enemin, enemax, nw)
1708 : 12326 : call htetra_get_onewk_wvals(tetra, ik_ibz, bcorr, nw, wvals, max_occ, nkibz, eig_ibz, weights)
1709 : :
1710 : 12326 : end subroutine htetra_get_onewk
1711 : : !!***
1712 : :
1713 : : !----------------------------------------------------------------------
1714 : :
1715 : : !!****f* m_htetra/htetra_get_onewk_wvals_zinv
1716 : : !! NAME
1717 : : !! htetra_get_onewk_wvals_zinv
1718 : : !!
1719 : : !! FUNCTION
1720 : : !! Calculate integration weights for 1/(z-E(k)) for a single k-point in the IBZ.
1721 : : !! Using either the implementation from:
1722 : : !! S. Kaprzyk, Computer Physics Communications 183, 347 (2012).
1723 : : !! or (TODO)
1724 : : !! P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984).
1725 : : !!
1726 : : !! INPUTS
1727 : : !! ik_ibz=Index of the k-point in the IBZ array
1728 : : !! bcorr=1 to include Blochl correction else 0.
1729 : : !! nw=number of energies in wvals
1730 : : !! nibz=number of irreducible kpoints
1731 : : !! wvals(nw)=Frequency points.
1732 : : !! eigen_ibz(nkibz)=eigenenergies for each k point
1733 : : !! opt: 1 for S. Kaprzyk routines, 2 for Lambin.
1734 : : !!
1735 : : !! OUTPUT
1736 : : !! weights(nw,2) = integration weights for
1737 : : !! Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
1738 : : !! for a given (band, k-point, spin).
1739 : : !! [erange(2)]: if present, weights are computed with an approximated asyntotic expression if
1740 : : !! real(z) is outside of this interval and with tetra if inside.
1741 : : !!
1742 : : !! SOURCE
1743 : :
1744 : 0 : subroutine htetra_get_onewk_wvals_zinv(tetra, ik_ibz, nz, zvals, max_occ, nkibz, eig_ibz, opt, cweights, erange)
1745 : :
1746 : : !Arguments ------------------------------------
1747 : : !scalars
1748 : : class(htetra_t), intent(inout) :: tetra
1749 : : integer,intent(in) :: ik_ibz,nz,nkibz,opt
1750 : : real(dp) ,intent(in) :: max_occ
1751 : : !arrays
1752 : : complex(dp),intent(in) :: zvals(nz)
1753 : : real(dp),optional,intent(in) :: erange(2)
1754 : : real(dp),intent(in) :: eig_ibz(nkibz)
1755 : : complex(dp),intent(out) :: cweights(nz)
1756 : :
1757 : : !Local variables-------------------------------
1758 : : !scalars
1759 : : integer :: itetra,isummit,tetra_total,tetra_count,iz
1760 : : real(dp) :: tweight
1761 : : !arrays
1762 : : integer :: ind_ibz(4),tetra_mibz(0:4)
1763 : : real(dp) :: eig(4), my_erange(2)
1764 : : complex(dp) :: verm(4), cw(4), verli(4)
1765 : : ! *********************************************************************
1766 : :
1767 [ # # ]: 0 : cweights = zero
1768 : : ! lazy evaluation of the mapping from k-points to tetrahedra
1769 [ # # ]: 0 : if (.not.allocated(tetra%ibz)) call htetra_init_mapping_ibz(tetra)
1770 : :
1771 [ # # ]: 0 : if (all(opt /= [1, 2])) then
1772 : 0 : ABI_ERROR(sjoin("Invalid opt:", itoa(opt)))
1773 : : end if
1774 : :
1775 [ # # ]: 0 : my_erange = [-huge(one), huge(one)]; if (present(erange)) my_erange = erange
1776 : :
1777 : : ! For each tetrahedron that belongs to this k-point
1778 : 0 : tetra_count = tetra%tetra_count(ik_ibz)
1779 : 0 : tetra_total = tetra%tetra_total(ik_ibz)
1780 [ # # ]: 0 : do itetra=1,tetra_count
1781 : :
1782 : 0 : call htetra_get_ibz(tetra, ik_ibz, itetra, tetra_mibz)
1783 : 0 : tweight = one * tetra_mibz(0) / tetra_total
1784 [ # # ]: 0 : do isummit=1,4
1785 : : ! Get mapping of each summit to eig_ibz
1786 : 0 : ind_ibz(isummit) = tetra_mibz(isummit)
1787 : 0 : eig(isummit) = eig_ibz(ind_ibz(isummit))
1788 : : end do
1789 : :
1790 : : ! Loop over frequencies
1791 [ # # ]: 0 : do iz=1,nz
1792 : :
1793 [ # # # # ]: 0 : if (real(zvals(iz)) >= my_erange(1) .and. real(zvals(iz)) <= my_erange(2)) then
1794 : 0 : select case(opt)
1795 : : case (1)
1796 [ # # ]: 0 : verm = zvals(iz) - eig
1797 : 0 : call SIM0TWOI(cw, VERLI, VERM)
1798 : : case (2)
1799 [ # # # ]: 0 : call get_onetetra_lambinvigneron(eig, zvals(iz), cw)
1800 : : end select
1801 : : else
1802 [ # # ]: 0 : cw = (one / (zvals(iz) - eig)) / four !* tetra%vv
1803 : : end if
1804 : :
1805 [ # # ]: 0 : do isummit=1,4
1806 [ # # ]: 0 : if (ind_ibz(isummit) /= ik_ibz) cycle
1807 : 0 : cweights(iz) = cweights(iz) + cw(isummit) * tweight * max_occ
1808 : : ! HM: This exit is important, avoids summing the same contribution more than once
1809 : 0 : exit
1810 : : end do
1811 : : end do
1812 : : end do ! itetra
1813 : :
1814 : 0 : end subroutine htetra_get_onewk_wvals_zinv
1815 : : !!***
1816 : :
1817 : : !----------------------------------------------------------------------
1818 : :
1819 : : !!****f* m_htetra/htetra_get_delta_mask
1820 : : !! NAME
1821 : : !! htetra_get_delta_mask
1822 : : !!
1823 : : !! FUNCTION
1824 : : !! Get a mask for the kpoints where the delta is finite
1825 : : !!
1826 : :
1827 : : subroutine htetra_get_delta_mask(tetra, eig_ibz, wvals, nw, nkpt, kmask, comm)
1828 : :
1829 : : !Arguments
1830 : : class(htetra_t), intent(in) :: tetra
1831 : : integer,intent(in) :: nw,nkpt,comm
1832 : : real(dp),intent(in) :: wvals(nw)
1833 : : real(dp),intent(in) :: eig_ibz(nkpt)
1834 : : integer,intent(out) :: kmask(nkpt)
1835 : :
1836 : : !Local variables-------------------------------
1837 : : integer :: ik_ibz,nprocs,my_rank,ierr, contrib
1838 : : integer :: tetra_count, itetra, isummit, ihash
1839 : : real(dp) :: emin,emax
1840 : : integer :: ind_ibz(4)
1841 : : real(dp) :: eig(4)
1842 : :
1843 : : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1844 : : kmask = 0
1845 : : ! For each bucket of tetrahedra
1846 : : do ihash=1,tetra%nbuckets
1847 : : if (mod(ihash,nprocs) /= my_rank) cycle
1848 : :
1849 : : ! For each tetrahedron
1850 : : tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2)
1851 : : do itetra=1,tetra_count
1852 : :
1853 : : ! Get mapping of each summit to eig_ibz
1854 : : do isummit=1,4
1855 : : ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit, itetra)
1856 : : eig(isummit) = eig_ibz(ind_ibz(isummit))
1857 : : end do
1858 : :
1859 : : ! Determine the energy range of the tetrahedra
1860 : : emin = minval(eig)
1861 : : emax = maxval(eig)
1862 : :
1863 : : ! Check if any value in wvals is between emin and emax
1864 : : contrib = 0; if (any(emin < wvals .and. wvals < emax)) contrib = 1
1865 : :
1866 : : ! Compute the union
1867 : : do isummit=1,4
1868 : : ik_ibz = ind_ibz(isummit)
1869 : : kmask(ik_ibz) = kmask(ik_ibz) + contrib
1870 : : end do
1871 : : end do ! itetra
1872 : : end do
1873 : :
1874 : : call xmpi_sum(kmask, comm, ierr)
1875 : :
1876 : : end subroutine htetra_get_delta_mask
1877 : : !!***
1878 : :
1879 : : !----------------------------------------------------------------------
1880 : :
1881 : : !!****f* m_htetra/htetra_wvals_weights
1882 : : !! NAME
1883 : : !! htetra_wvals_weights
1884 : : !!
1885 : : !! FUNCTION
1886 : : !! Emulates the behaviour of the previous tetrahedron implementation but
1887 : : !! taking a list of energies as input.
1888 : : !!
1889 : : !! HM: I find that in many routines its better to change the implementation
1890 : : !! and accumulate the tetrahedron weights in the same way as the
1891 : : !! gaussian smearing weights using htetra_get_onewk_wvals. However this requires
1892 : : !! some refactoring of the code. I provide this routine to make it easier
1893 : : !! to transition to the new tetrahedron implementation without refactoring.
1894 : : !! Looping over tetrahedra (i.e. using tetra_blochl_weights) is currently faster
1895 : : !! than looping over k-points.
1896 : : !!
1897 : : !! MG: Note, however, that tetra_blochl_weights requires more memory as
1898 : : !! one has to allocate dweight(nw,nkpt),tweight(nw,nkpt) and the size of the arrays increases
1899 : : !! quickly with the k-mesh and the number of frequencies (propto mband)
1900 : : !!
1901 : : !! INPUTS
1902 : : !!
1903 : : !! OUTPUT
1904 : : !!
1905 : : !! SOURCE
1906 : :
1907 : 118 : subroutine htetra_wvals_weights(tetra, eig_ibz, nw, wvals, max_occ, nkpt, opt, tweight, dweight, comm)
1908 : :
1909 : : !Arguments ------------------------------------
1910 : : !scalars
1911 : : integer,intent(in) :: nw,nkpt,opt,comm
1912 : : class(htetra_t), intent(in) :: tetra
1913 : : real(dp) ,intent(in) :: max_occ
1914 : : !arrays
1915 : : real(dp),intent(in) :: eig_ibz(nkpt)
1916 : : real(dp),intent(out) :: dweight(nw,nkpt),tweight(nw,nkpt)
1917 : :
1918 : : !Local variables-------------------------------
1919 : : !scalars
1920 : : integer :: ik_ibz,multiplicity,nprocs,my_rank,ierr
1921 : : integer :: tetra_count, itetra, isummit, ihash
1922 : : !arrays
1923 : : integer :: ind_ibz(4)
1924 : 118 : real(dp) :: eig(4), wvals(nw), dweight_tmp(4,nw),tweight_tmp(4,nw)
1925 : : ! *********************************************************************
1926 : :
1927 [ + + + + : 10409956 : tweight = zero; dweight = zero
+ + + + ]
1928 : 118 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1929 : :
1930 : : ! For each bucket of tetrahedra
1931 [ + + ]: 5994 : do ihash=1,tetra%nbuckets
1932 [ + - ]: 5876 : if (mod(ihash, nprocs) /= my_rank) cycle
1933 : :
1934 : : ! For each tetrahedron
1935 : 5876 : tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2)
1936 [ + + ]: 105152 : do itetra=1,tetra_count
1937 : :
1938 : : ! Get mapping of each summit to eig_ibz
1939 [ + + ]: 495790 : do isummit=1,4
1940 : 396632 : ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit, itetra)
1941 : 495790 : eig(isummit) = eig_ibz(ind_ibz(isummit))
1942 : : end do
1943 : :
1944 : : ! Sort energies before calling get_onetetra_blochl
1945 : 99158 : call sort_4tetra(eig, ind_ibz)
1946 : :
1947 : : ! Get tetrahedron weights
1948 : 198316 : select case (opt)
1949 : : case (0:1)
1950 : 99158 : call get_onetetra_blochl(eig, wvals, nw, opt, tweight_tmp, dweight_tmp)
1951 : : case (2)
1952 : 0 : call get_onetetetra_lambinvigneron_imag(eig, wvals, nw, dweight_tmp)
1953 [ + - - - : 99158 : tweight_tmp = zero
- - - ]
1954 : : end select
1955 : :
1956 : : ! Accumulate the contributions
1957 : 99158 : multiplicity = tetra%unique_tetra(ihash)%indexes(0, itetra)
1958 [ + + ]: 501666 : do isummit=1,4
1959 : 396632 : ik_ibz = ind_ibz(isummit)
1960 [ + + ]: 337052272 : dweight(:,ik_ibz) = dweight(:,ik_ibz) + dweight_tmp(isummit,:) * multiplicity * max_occ
1961 [ + + ]: 337151430 : tweight(:,ik_ibz) = tweight(:,ik_ibz) + tweight_tmp(isummit,:) * multiplicity * max_occ
1962 : : end do
1963 : : end do ! itetra
1964 : : end do
1965 : :
1966 : : ! Rescale weights
1967 : 118 : select case(tetra%opt)
1968 : : case (1)
1969 [ # # ]: 0 : do ik_ibz=1,tetra%nkibz
1970 [ # # ]: 0 : dweight(:,ik_ibz) = dweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%tetra_total(ik_ibz) / tetra%nkbz
1971 [ # # ]: 0 : tweight(:,ik_ibz) = tweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%tetra_total(ik_ibz) / tetra%nkbz
1972 : : end do
1973 : : case (2)
1974 [ + + + + ]: 5204978 : dweight = dweight*tetra%vv / 4.0_dp
1975 [ - + - + : 5205096 : tweight = tweight*tetra%vv / 4.0_dp
+ + + ]
1976 : : end select
1977 : :
1978 : 118 : call xmpi_sum(dweight, comm, ierr)
1979 : 118 : call xmpi_sum(tweight, comm, ierr)
1980 : :
1981 : 118 : end subroutine htetra_wvals_weights
1982 : : !!***
1983 : :
1984 : : !----------------------------------------------------------------------
1985 : :
1986 : : !!****f* m_htetra/htetra_wvals_weights_delta
1987 : : !! NAME
1988 : : !! htetra_wvals_weights_delta
1989 : : !!
1990 : : !! FUNCTION
1991 : : !! Same as above but computing only delta for performance and memory
1992 : : !! HM: Should find a clean way to avoid copy paste routine
1993 : : !!
1994 : : !! SOURCE
1995 : :
1996 : 30216 : subroutine htetra_wvals_weights_delta(tetra, eig_ibz, nw, wvals, max_occ, nkpt, opt, dweight, comm)
1997 : :
1998 : : !Arguments ------------------------------------
1999 : : !scalars
2000 : : class(htetra_t), intent(in) :: tetra
2001 : : integer,intent(in) :: nw, nkpt, opt, comm
2002 : : real(dp),intent(in) :: max_occ
2003 : : !arrays
2004 : : real(dp),intent(in) :: eig_ibz(nkpt), wvals(nw)
2005 : : real(dp),intent(out) :: dweight(nw, nkpt)
2006 : :
2007 : : !Local variables-------------------------------
2008 : : !scalars
2009 : : integer :: ik_ibz,multiplicity,nprocs,my_rank,ierr
2010 : : integer :: tetra_count, itetra, isummit, ihash
2011 : : !arrays
2012 : : integer :: ind_ibz(4)
2013 : 30216 : real(dp) :: eig(4), dweight_tmp(4,nw),tweight_tmp(4,nw)
2014 : : ! *********************************************************************
2015 : :
2016 [ + + + + ]: 13066728 : dweight = zero
2017 : 30216 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2018 : :
2019 : : ! For each bucket of tetrahedra
2020 [ + + ]: 6548472 : do ihash=1,tetra%nbuckets
2021 [ + - ]: 6518256 : if (mod(ihash, nprocs) /= my_rank) cycle
2022 : :
2023 : : ! For each tetrahedron
2024 : 6518256 : tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2)
2025 [ + + ]: 129322080 : do itetra=1,tetra_count
2026 : :
2027 : : ! Get mapping of each summit to eig_ibz
2028 [ + + ]: 613868040 : do isummit=1,4
2029 : 491094432 : ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit,itetra)
2030 : 613868040 : eig(isummit) = eig_ibz(ind_ibz(isummit))
2031 : : end do
2032 : :
2033 : : ! Sort energies before calling get_onetetra_blochl
2034 : 122773608 : call sort_4tetra(eig, ind_ibz)
2035 : :
2036 : : ! Get tetrahedron weights
2037 : 245547216 : select case (opt)
2038 : : case (0:1)
2039 : 122773608 : call get_onetetra_blochl(eig, wvals, nw, opt, tweight_tmp, dweight_tmp)
2040 : : case (2)
2041 [ + - - ]: 122773608 : call get_onetetetra_lambinvigneron_imag(eig, wvals, nw, dweight_tmp)
2042 : : end select
2043 : :
2044 : : ! Accumulate the contributions
2045 : 122773608 : multiplicity = tetra%unique_tetra(ihash)%indexes(0,itetra)
2046 [ + + ]: 620386296 : do isummit=1,4
2047 : 491094432 : ik_ibz = ind_ibz(isummit)
2048 [ + + ]: 1104962472 : dweight(:,ik_ibz) = dweight(:,ik_ibz) + dweight_tmp(isummit,:)*multiplicity*max_occ
2049 : : end do
2050 : : end do ! itetra
2051 : : end do
2052 : :
2053 : : ! Rescale weights
2054 : 30216 : select case(tetra%opt)
2055 : : case (1)
2056 [ # # ]: 0 : do ik_ibz=1,tetra%nkibz
2057 [ # # ]: 0 : dweight(:,ik_ibz) = dweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%tetra_total(ik_ibz) / tetra%nkbz
2058 : : end do
2059 : : case (2)
2060 [ - + - + : 13066728 : dweight = dweight * tetra%vv / 4.0_dp
+ + + ]
2061 : : end select
2062 : :
2063 : 30216 : call xmpi_sum(dweight, comm, ierr)
2064 : :
2065 : 30216 : end subroutine htetra_wvals_weights_delta
2066 : : !!***
2067 : :
2068 : : !----------------------------------------------------------------------
2069 : :
2070 : : !!****f* m_htetra/htetra_blochl_weights
2071 : : !! NAME
2072 : : !! htetra_blochl_weights
2073 : : !!
2074 : : !! FUNCTION
2075 : : !! Emulates the behaviour of the previous tetrahedron implementation.
2076 : : !! IBZ weights are included.
2077 : : !!
2078 : : !! INPUTS
2079 : : !!
2080 : : !! OUTPUT
2081 : : !!
2082 : : !! SOURCE
2083 : :
2084 : 118 : subroutine htetra_blochl_weights(tetra, eig_ibz, enemin, enemax, max_occ, nw, nkpt, bcorr, tweight, dweight, comm)
2085 : :
2086 : : !Arguments ------------------------------------
2087 : : !scalars
2088 : : class(htetra_t), intent(in) :: tetra
2089 : : integer,intent(in) :: nw,nkpt, bcorr, comm
2090 : : real(dp),intent(in) :: enemax, enemin, max_occ
2091 : : !arrays
2092 : : real(dp),intent(in) :: eig_ibz(nkpt)
2093 : : real(dp),intent(out) :: dweight(nw,nkpt), tweight(nw,nkpt)
2094 : :
2095 : : !Local variables-------------------------------
2096 : 236 : real(dp) :: wvals(nw)
2097 : : ! *********************************************************************
2098 : :
2099 : 118 : wvals = linspace(enemin, enemax, nw)
2100 : 118 : call htetra_wvals_weights(tetra,eig_ibz,nw,wvals,max_occ,nkpt,bcorr,tweight,dweight,comm)
2101 : :
2102 : 118 : end subroutine htetra_blochl_weights
2103 : : !!***
2104 : :
2105 : : !----------------------------------------------------------------------
2106 : :
2107 : : !!****f* m_htetra/htetra_blochl_weights_wvals_zinv
2108 : : !! NAME
2109 : : !! htetra_blochl_weights_wvals_zinv
2110 : : !!
2111 : : !! FUNCTION
2112 : : !! The same as htetra_get_onewk_wvals_zinv but looping over tetrahedra
2113 : : !! which is more efficient
2114 : : !!
2115 : : !! INPUTS
2116 : : !! nz: Number of frequencies
2117 : : !! zvals(nw): z-values
2118 : : !! max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)
2119 : : !! nkpt=number of irreducible kpoints
2120 : : !! zinv_opt:
2121 : : !! 1 for S. Kaprzyk routines,
2122 : : !! 2 for Lambin-Vigneron.
2123 : : !! [erange(2)]: if present, weights are computed with a standard quadrature method if
2124 : : !! real(z) is outside of this interval and with tetra if inside.
2125 : : !! comm=MPI communicator
2126 : : !!
2127 : : !! OUTPUT
2128 : : !!
2129 : : !! SOURCE
2130 : :
2131 : 18048 : subroutine htetra_weights_wvals_zinv(tetra, eig_ibz, nz, zvals, max_occ, nkpt, zinv_opt, cweight, comm, erange)
2132 : :
2133 : : !Arguments ------------------------------------
2134 : : !scalars
2135 : : integer,intent(in) :: nz, nkpt, zinv_opt, comm
2136 : : class(htetra_t), intent(in) :: tetra
2137 : : real(dp) ,intent(in) :: max_occ
2138 : : !arrays
2139 : : real(dp),intent(in) :: eig_ibz(nkpt)
2140 : : real(dp),optional,intent(in) :: erange(2)
2141 : : complex(dp),intent(in) :: zvals(nz)
2142 : : complex(dp),intent(out) :: cweight(nz, nkpt)
2143 : :
2144 : : !Local variables-------------------------------
2145 : : !scalars
2146 : : integer :: ik_ibz, iz, multiplicity, nprocs, my_rank, ierr, ii, jj, kk, esumk
2147 : : integer :: tetra_count, itetra, isummit, ihash
2148 : : !arrays
2149 : : integer :: ind_ibz(4)
2150 : : real(dp) :: eig(4), my_erange(2)
2151 : : complex(dp) :: cw(4), verli(4), verm(4), aw(4), bw(4) !, cw_lw(4)
2152 : 36096 : real(dp) :: rwg(nz, 4)
2153 : : ! *********************************************************************
2154 : :
2155 : 18048 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2156 [ + + + + ]: 8541696 : cweight = zero
2157 : :
2158 [ + - ]: 18048 : my_erange = [-huge(one), huge(one)]; if (present(erange)) my_erange = erange
2159 : :
2160 : : ! For each bucket of tetrahedra
2161 [ + + ]: 207552 : do ihash=1,tetra%nbuckets
2162 [ + - ]: 189504 : if (mod(ihash, nprocs) /= my_rank) cycle
2163 : :
2164 : : ! For each tetrahedron that belongs to this k-point
2165 : 189504 : tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2)
2166 [ + + ]: 1461888 : do itetra=1,tetra_count
2167 : :
2168 : : ! Get mapping of each summit to eig_ibz
2169 [ + + ]: 6271680 : do isummit=1,4
2170 : 5017344 : ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit, itetra)
2171 : 6271680 : eig(isummit) = eig_ibz(ind_ibz(isummit))
2172 : : end do
2173 : :
2174 : : ! Get multiplicity
2175 : 1254336 : multiplicity = tetra%unique_tetra(ihash)%indexes(0, itetra)
2176 : :
2177 : : ! Sort energies before calling get_onetetra_lambinvigneron
2178 : : ! SIM0TWOI does not require sorted energies but since we call the routine
2179 : : ! to fix get_onetetra_lambinvigneron we need to sort here.
2180 : 1254336 : call sort_4tetra(eig, ind_ibz)
2181 : :
2182 [ + - ]: 1443840 : if (zinv_opt == 1) then
2183 : : ! Loop over frequencies
2184 [ + + ]: 56418432 : do iz=1,nz
2185 : :
2186 : : ! Get tetrahedron weights
2187 [ - + - + ]: 55164096 : if (real(zvals(iz)) >= my_erange(1) .and. real(zvals(iz)) <= my_erange(2)) then
2188 : :
2189 : 55164096 : select case(zinv_opt)
2190 : : case (1)
2191 [ + + ]: 275820480 : verm = zvals(iz) - eig
2192 : 55164096 : call SIM0TWOI(cw, VERLI, VERM)
2193 : :
2194 : : !call get_onetetra_lambinvigneron(eig, zvals(iz), cw_lw)
2195 : : !if (any(abs(real(cw_lw(:)) / (real(cw(:)))) > 1.1)) then
2196 : : ! do ierr=1,4
2197 : : ! !write(std_out, *) "simte vs lw:", cw(ierr), cw_lw(ierr)
2198 : : ! write(std_out, *) "cw_lw / (simtet)", cw_lw(ierr) / (cw(ierr))
2199 : : ! end do
2200 : : !end if
2201 : :
2202 : : case (2)
2203 : : call get_onetetra_lambinvigneron(eig, zvals(iz), cw)
2204 : : end select
2205 : :
2206 : : else
2207 : : ! Use asymptotic expansion of integral for large z.
2208 [ # # # # ]: 0 : aw = (eig + sum(eig)) / five
2209 [ # # ]: 0 : do ii=1,4
2210 : 0 : bw(ii) = zero
2211 [ # # ]: 0 : do jj=1,4
2212 [ # # ]: 0 : if (jj == ii) cycle
2213 : : esumk = zero
2214 [ # # ]: 0 : do kk=1,4
2215 [ # # ]: 0 : if (kk == jj) cycle
2216 : 0 : esumk = esumk + (eig(kk) - eig(jj)) ** 2
2217 : : end do
2218 : 0 : bw(ii) = bw(ii) + three * (eig(jj) - eig(ii)) ** 2 + esumk
2219 : : end do
2220 : : end do
2221 [ # # ]: 0 : bw = bw / 300_dp
2222 [ # # ]: 0 : cw = one / (zvals(iz) - aw - bw / zvals(iz)) / four
2223 : : !This for naive integration
2224 : : !cw = (one / (zvals(iz) - eig)) / four
2225 : : end if
2226 : :
2227 : : ! Accumulate contributions
2228 [ + + ]: 277074816 : do isummit=1,4
2229 : 220656384 : ik_ibz = ind_ibz(isummit)
2230 : 275820480 : cweight(iz, ik_ibz) = cweight(iz, ik_ibz) + cw(isummit) * multiplicity * max_occ
2231 : : end do
2232 : :
2233 : : end do ! iz
2234 : :
2235 : : else
2236 [ # # # # : 0 : call get_onetetra_ppart_lv(nz, real(zvals), eig, rwg)
# # ]
2237 : :
2238 [ # # ]: 0 : do iz=1,nz
2239 : : ! Accumulate contributions
2240 [ # # ]: 0 : do isummit=1,4
2241 : 0 : ik_ibz = ind_ibz(isummit)
2242 : 0 : cweight(iz, ik_ibz) = cweight(iz, ik_ibz) + rwg(iz, isummit) * multiplicity * max_occ
2243 : : end do
2244 : : end do
2245 : : endif
2246 : :
2247 : : end do ! itetra
2248 : : end do
2249 : :
2250 : : ! Rescale weights
2251 : 18048 : select case(tetra%opt)
2252 : : case (1)
2253 [ # # ]: 0 : do ik_ibz=1,tetra%nkibz
2254 [ # # ]: 0 : cweight(:,ik_ibz) = cweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%nkbz / tetra%tetra_total(ik_ibz)
2255 : : end do
2256 : : case (2)
2257 [ - + - + : 8541696 : cweight = cweight * tetra%vv
+ + + ]
2258 : : end select
2259 : :
2260 : 18048 : call xmpi_sum(cweight, comm, ierr)
2261 : :
2262 : 18048 : end subroutine htetra_weights_wvals_zinv
2263 : : !!***
2264 : :
2265 : : !----------------------------------------------------------------------
2266 : :
2267 : : !!****f* m_htetra/sort_4tetra
2268 : : !! NAME
2269 : : !! sort_4tetra
2270 : : !!
2271 : : !! FUNCTION
2272 : : !! Sort double precision array list(4) into ascending numerical order
2273 : : !! while making corresponding rearrangement of the integer array iperm.
2274 : : !!
2275 : : !! Taken from: https://stackoverflow.com/questions/6145364/sort-4-number-with-few-comparisons
2276 : : !!
2277 : : !! INPUTS
2278 : : !! list(4) intent(inout) list of double precision numbers to be sorted
2279 : : !! perm(4) intent(inout) iperm(i)=i (very important)
2280 : : !!
2281 : : !! OUTPUT
2282 : : !! list(4) sorted list
2283 : : !! perm(4) index of permutation given the right ascending order
2284 : : !!
2285 : : !! SOURCE
2286 : :
2287 : 139837262 : pure subroutine sort_4tetra(list, perm)
2288 : :
2289 : : integer, intent(inout) :: perm(4)
2290 : : real(dp), intent(inout) :: list(4)
2291 : :
2292 : : !Local variables-------------------------------
2293 : : integer :: ia,ib,ic,id,ilow1,ilow2,ihigh1,ihigh2
2294 : : integer :: ilowest,ihighest,imiddle1,imiddle2
2295 : : real(dp) :: va,vb,vc,vd,vlow1,vlow2,vhigh1,vhigh2
2296 : : real(dp) :: vlowest,vhighest,vmiddle1,vmiddle2
2297 : :
2298 : 139837262 : va = list(1); ia = perm(1)
2299 : 139837262 : vb = list(2); ib = perm(2)
2300 : 139837262 : vc = list(3); ic = perm(3)
2301 : 139837262 : vd = list(4); id = perm(4)
2302 : :
2303 [ + + ]: 139837262 : if (va < vb) then
2304 : : vlow1 = va; vhigh1 = vb
2305 : : ilow1 = ia; ihigh1 = ib
2306 : : else
2307 : 72940132 : vlow1 = vb; vhigh1 = va
2308 : 72940132 : ilow1 = ib; ihigh1 = ia
2309 : : endif
2310 : :
2311 [ + + ]: 139837262 : if (vc < vd) then
2312 : : vlow2 = vc; vhigh2 = vd
2313 : : ilow2 = ic; ihigh2 = id
2314 : : else
2315 : 73080767 : vlow2 = vd; vhigh2 = vc
2316 : 73080767 : ilow2 = id; ihigh2 = ic
2317 : : endif
2318 : :
2319 [ + + ]: 139837262 : if (vlow1 < vlow2) then
2320 : : vlowest = vlow1; vmiddle1 = vlow2
2321 : : ilowest = ilow1; imiddle1 = ilow2
2322 : : else
2323 : 71288294 : vlowest = vlow2; vmiddle1 = vlow1
2324 : 71288294 : ilowest = ilow2; imiddle1 = ilow1
2325 : : endif
2326 : :
2327 [ + + ]: 139837262 : if (vhigh1 > vhigh2) then
2328 : : vhighest = vhigh1; vmiddle2 = vhigh2
2329 : : ihighest = ihigh1; imiddle2 = ihigh2
2330 : : else
2331 : 69068729 : vhighest = vhigh2; vmiddle2 = vhigh1
2332 : 69068729 : ihighest = ihigh2; imiddle2 = ihigh1
2333 : : endif
2334 : :
2335 [ + + ]: 139837262 : if (vmiddle1 < vmiddle2) then
2336 [ + + ]: 365584920 : list = [vlowest, vmiddle1, vmiddle2, vhighest]
2337 [ + + ]: 365584920 : perm = [ilowest, imiddle1, imiddle2, ihighest]
2338 : : else
2339 [ + + ]: 333601390 : list = [vlowest, vmiddle2, vmiddle1, vhighest]
2340 [ + + ]: 333601390 : perm = [ilowest, imiddle2, imiddle1, ihighest]
2341 : : endif
2342 : :
2343 : 139837262 : end subroutine sort_4tetra
2344 : : !!***
2345 : :
2346 : : !----------------------------------------------------------------------
2347 : :
2348 : : !!****f* m_numeric_tools/sort_4tetra_int
2349 : : !! NAME
2350 : : !! sort_4tetra_int
2351 : : !!
2352 : : !! FUNCTION
2353 : : !!
2354 : : !! INPUTS
2355 : : !!
2356 : : !! OUTPUT
2357 : : !!
2358 : : !! SOURCE
2359 : :
2360 : 8456046 : pure subroutine sort_4tetra_int(list)
2361 : :
2362 : : integer, intent(inout) :: list(4)
2363 : :
2364 : : !Local variables-------------------------------
2365 : : integer :: va,vb,vc,vd, vlow1,vlow2,vhigh1,vhigh2
2366 : : integer :: vlowest,vhighest, vmiddle1,vmiddle2
2367 : :
2368 : 8456046 : va = list(1)
2369 : 8456046 : vb = list(2)
2370 : 8456046 : vc = list(3)
2371 : 8456046 : vd = list(4)
2372 : :
2373 [ + + ]: 8456046 : if (va < vb) then
2374 : : vlow1 = va; vhigh1 = vb
2375 : : else
2376 : 4340218 : vlow1 = vb; vhigh1 = va
2377 : : endif
2378 : :
2379 [ + + ]: 8456046 : if (vc < vd) then
2380 : : vlow2 = vc; vhigh2 = vd
2381 : : else
2382 : 4257565 : vlow2 = vd; vhigh2 = vc
2383 : : endif
2384 : :
2385 [ + + ]: 8456046 : if (vlow1 < vlow2) then
2386 : : vlowest = vlow1; vmiddle1 = vlow2
2387 : : else
2388 : 4421792 : vlowest = vlow2; vmiddle1 = vlow1
2389 : : endif
2390 : :
2391 [ + + ]: 8456046 : if (vhigh1 > vhigh2) then
2392 : : vhighest = vhigh1; vmiddle2 = vhigh2
2393 : : else
2394 : 4589119 : vhighest = vhigh2; vmiddle2 = vhigh1
2395 : : endif
2396 : :
2397 [ + + ]: 8456046 : if (vmiddle1 < vmiddle2) then
2398 [ + + ]: 24813400 : list = [vlowest, vmiddle1, vmiddle2, vhighest]
2399 : : else
2400 [ + + ]: 17466830 : list = [vlowest, vmiddle2, vmiddle1, vhighest]
2401 : : endif
2402 : :
2403 : 8456046 : end subroutine sort_4tetra_int
2404 : : !!***
2405 : :
2406 : : !!****f* m_htetra/get_onetetra_ppart_lv
2407 : : !! NAME
2408 : : !! get_onetetra_ppart_lv
2409 : : !!
2410 : : !! FUNCTION
2411 : : !! Compute the complex weights according to: P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984)
2412 : : !!
2413 : : !! INPUTS
2414 : : !! nw
2415 : : !! wvals: energy to evaluate the weights at
2416 : : !! eig: eigenvalues at the corners of the tetrahedron
2417 : : !!
2418 : : !! OUTPUT
2419 : : !! rwg(nw, 4)
2420 : : !!
2421 : : !! SOURCE
2422 : :
2423 : 0 : pure subroutine get_onetetra_ppart_lv(nw, wvals, eig, rwg)
2424 : :
2425 : : integer,intent(in) :: nw
2426 : : real(dp), intent(in) :: wvals(nw), eig(4)
2427 : : real(dp), intent(out) :: rwg(nw, 4)
2428 : :
2429 : : !Local variables-------------------------------
2430 : : integer :: ii, iw !jj,
2431 : : real(dp),parameter :: tol = tol14
2432 : : !real(dp),parameter :: tol = tol20
2433 : : !real(dp),parameter :: tol = tol30
2434 : : real(dp) :: D12,D13,D14,D23,D24,D34
2435 : : real(dp) :: e10, e20, e30, e31, e32, e21
2436 : : real(dp) :: inv_e10, inv_e21, inv_e20, inv_e30, inv_e31, inv_e32
2437 : 0 : real(dp) :: E0(nw), E1(nw), E2(nw), E3(nw)
2438 : : ! *********************************************************************
2439 : :
2440 : : ! Then the energy differences, for the coefficients. Must always be positive, I hope.
2441 : 0 : D12 = eig(2) - eig(1)
2442 : 0 : D13 = eig(3) - eig(1)
2443 : 0 : D14 = eig(4) - eig(1)
2444 : 0 : D23 = eig(3) - eig(2)
2445 : 0 : D24 = eig(4) - eig(2)
2446 : 0 : D34 = eig(4) - eig(3)
2447 : :
2448 : : ! Now get the actual weights
2449 : : ! Notations
2450 : : ! eij = e_i - e_j
2451 : : ! Ej = E - e_j
2452 : :
2453 : 0 : e10 = huge(one); e20 = huge(one); e30 = huge(one);
2454 : 0 : e31 = huge(one); e32 = huge(one); e21 = huge(one)
2455 : 0 : inv_e10 = huge(one); inv_e21 = huge(one); inv_e20 = huge(one)
2456 : 0 : inv_e30 = huge(one); inv_e31 = huge(one); inv_e32 = huge(one)
2457 : :
2458 [ # # ]: 0 : E0 = wvals(:) - eig(1)
2459 [ # # ]: 0 : E1 = wvals(:) - eig(2)
2460 [ # # ]: 0 : E2 = wvals(:) - eig(3)
2461 [ # # ]: 0 : E3 = wvals(:) - eig(4)
2462 : :
2463 : : #if 0
2464 : : where (abs(E0) < tol1)
2465 : : E0 = tol1
2466 : : end where
2467 : : where (abs(E1) < tol1)
2468 : : E1 = tol1
2469 : : end where
2470 : : where (abs(E2) < tol1)
2471 : : E2 = tol1
2472 : : end where
2473 : : where (abs(E3) < tol1)
2474 : : E3 = tol1
2475 : : end where
2476 : : #endif
2477 : :
2478 : : ! e1=e2=e3=e4
2479 [ # # ]: 0 : if (D12 + D23 + D34 < tol) then
2480 [ # # ]: 0 : do ii=1,4
2481 : : !rwg(:, ii) = 0.25_dp / (wvals - eig(ii))
2482 [ # # ]: 0 : rwg(:, ii) = 0.25_dp / E0
2483 : : end do
2484 : :
2485 : : ! e2=e3=e4
2486 [ # # ]: 0 : else if (D23 + D34 < tol) then
2487 : 0 : e10 = eig(2) - eig(1); inv_e10 = one / e10
2488 : :
2489 [ # # ]: 0 : do iw=1,nw
2490 : : rwg(iw, 1) = &
2491 : : three * E0(iw)**2 * E1(iw) * inv_e10**4 * log(abs(E1(iw) / E0(iw))) &
2492 : 0 : + 1.5_dp * E1(iw) * (two * E0(iw) + e10) * inv_e10**3 + inv_e10
2493 : : end do
2494 : :
2495 [ # # ]: 0 : do iw=1,nw
2496 : : rwg(iw, 2) = &
2497 : : E0(iw) ** 3 * inv_e10**4 * log(abs(E0(iw) / E1(iw))) &
2498 : 0 : - (six * E0(iw)**2 + three * E0(iw) * e10 + two * e10**2) * inv_e10**3 / six
2499 : : end do
2500 [ # # ]: 0 : rwg(:,3) = rwg(:,2)
2501 [ # # ]: 0 : rwg(:,4) = rwg(:,2)
2502 : :
2503 : : !rwg = zero
2504 : :
2505 : : ! e1=e2=e3
2506 [ # # ]: 0 : else if (D12 + D23 < tol) then
2507 : :
2508 : 0 : e30 = eig(4) - eig(1); inv_e30 = one / e30
2509 : :
2510 [ # # ]: 0 : do iw=1,nw
2511 : : rwg(iw, 1) = &
2512 : : E3(iw)**3 * inv_e30**4 * log(abs(E3(iw) / E0(iw))) &
2513 : 0 : + (six * E3(iw)**2 - three * E3(iw) * e30 + two * e30**2) * inv_e30**3 / six
2514 : : end do
2515 [ # # ]: 0 : rwg(:,2) = rwg(:,1)
2516 [ # # ]: 0 : rwg(:,3) = rwg(:,1)
2517 : :
2518 [ # # ]: 0 : do iw=1,nw
2519 : : rwg(iw, 4) = &
2520 : : three * E0(iw) * E3(iw)**2 * inv_e30**4 * log(abs(E0(iw) / E3(iw))) &
2521 : 0 : - 1.5_dp * E0(iw) * (two * E3(iw) - e30) * inv_e30**3 - inv_e30
2522 : : end do
2523 : :
2524 : : !rwg = zero
2525 : :
2526 : : ! e1=e2 < e3=e4
2527 [ # # ]: 0 : else if (D12 + D34 < tol) then
2528 : :
2529 : 0 : e20 = eig(3) - eig(1); inv_e20 = one / e20
2530 : :
2531 [ # # ]: 0 : do iw=1,nw
2532 : : rwg(iw, 1) = &
2533 : : three * E0(iw) * E2(iw)**2 * inv_e20**4 * log(abs(E0(iw) / E2(iw))) &
2534 : 0 : - 1.5_dp * E0(iw) * (two * E2(iw) - e20) * inv_e20**3 - inv_e20
2535 : : end do
2536 [ # # ]: 0 : rwg(:,2) = rwg(:,1)
2537 : :
2538 [ # # ]: 0 : do iw=1,nw
2539 : : rwg(iw, 3) = &
2540 : : three * E0(iw)** 2 * E2(iw) * inv_e20**4 * log(abs(E2(iw) / E0(iw))) &
2541 : 0 : + 1.5_dp * E2(iw) * (two * E0(iw) + e20) * inv_e20**3 + inv_e20
2542 : : end do
2543 [ # # ]: 0 : rwg(:,4) = rwg(:,3)
2544 : :
2545 : : !rwg = zero
2546 : :
2547 : : ! e3=e4
2548 [ # # ]: 0 : else if (D34 < tol) then
2549 : :
2550 : 0 : e10 = eig(2) - eig(1); inv_e10 = one / e10
2551 : 0 : e20 = eig(3) - eig(1); inv_e20 = one / e20
2552 : 0 : e21 = eig(3) - eig(2); inv_e21 = one / e21
2553 : :
2554 [ # # ]: 0 : do iw=1,nw
2555 : : rwg(iw, 1) = &
2556 : : E0(iw)**2 * inv_e20**2 * inv_e10 &
2557 : : * (one + (-two * E2(iw) * inv_e20 - E1(iw) * inv_e10) * log(abs(E0(iw)))) &
2558 : : - E2(iw)**2 * inv_e20**2 * inv_e21 &
2559 : : * (one + (two * E0(iw) * inv_e20 + E1(iw) * inv_e21) * log(abs(E2(iw)))) &
2560 : 0 : + E1(iw)**3 * inv_e10**2 * inv_e21**2 * log(abs(E1(iw)))
2561 : : end do
2562 : :
2563 [ # # ]: 0 : do iw=1,nw
2564 : : rwg(iw, 2) = &
2565 : : -E1(iw)**2 * inv_e21**2 * inv_e10 &
2566 : : * (one + (-two * E2(iw) * inv_e21 + E0(iw) * inv_e10) * log(abs(E1(iw)))) &
2567 : : - E2(iw)**2 * inv_e21**2 * inv_e20 &
2568 : : * (one + (two * E1(iw) * inv_e21 + E0(iw) * inv_e20) * log(abs(E2(iw)))) &
2569 : 0 : + E0(iw)**3 * inv_e10**2 * inv_e20**2 * log(abs(E0(iw)))
2570 : : end do
2571 : :
2572 [ # # ]: 0 : do iw=1,nw
2573 : : rwg(iw, 3) = &
2574 : : +E0(iw)**3 * inv_e10 * inv_e20**3 * log(abs(E0(iw))) &
2575 : : -E1(iw)**3 * inv_e10 * inv_e21**3 * log(abs(E1(iw))) &
2576 : : +E2(iw) * inv_e20 * inv_e21 &
2577 : : * (half + E0(iw) * inv_e20 + E1(iw) * inv_e21 &
2578 : : +(E0(iw)**2 * inv_e20**2 + E1(iw)**2 * inv_e21**2 + &
2579 : 0 : E0(iw) * E1(iw) * inv_e20 * inv_e21) * log(abs(E2(iw))))
2580 : : end do
2581 : :
2582 [ # # ]: 0 : rwg(:,4) = rwg(:,3)
2583 : :
2584 : : !rwg = zero
2585 : :
2586 : : ! e2=e3
2587 [ # # ]: 0 : else if (D23 < tol) then
2588 : :
2589 : 0 : e10 = eig(2) - eig(1); inv_e10 = one / e10
2590 : 0 : e30 = eig(4) - eig(1); inv_e30 = one / e30
2591 : 0 : e31 = eig(4) - eig(2); inv_e31 = one / e31
2592 : :
2593 [ # # ]: 0 : do iw=1,nw
2594 : : rwg(iw, 1) = &
2595 : : E0(iw)**2 * inv_e10**2 * inv_e30 &
2596 : : * (one - (two * E1(iw) * inv_e10 + E3(iw) * inv_e30) * log(abs(E0(iw)))) &
2597 : : + E1(iw)**2 * inv_e10**2 * inv_e31 &
2598 : : * (one + (+two * E0(iw) * inv_e10 - E3(iw) * inv_e31) * log(abs(E1(iw)))) &
2599 : 0 : + E3(iw)**3 * inv_e30**2 * inv_e31**2 * log(abs(E3(iw)))
2600 : : end do
2601 : :
2602 [ # # ]: 0 : do iw=1,nw
2603 : : rwg(iw, 2) = &
2604 : : E0(iw)**3 * inv_e30 * inv_e10**3 * log(abs(E0(iw))) &
2605 : : + E3(iw)**3 * inv_e30 * inv_e31**3 * log(abs(E3(iw))) &
2606 : : - E1(iw) * inv_e10 * inv_e31 &
2607 : : * (half + E0(iw) * inv_e10 - E3(iw) * inv_e31 + &
2608 : 0 : (E0(iw)**2 * inv_e10**2 + E3(iw)**2 * inv_e31**2 - E0(iw) * E3(iw) * inv_e10 * inv_e31) * log(abs(E1(iw))))
2609 : : end do
2610 [ # # ]: 0 : rwg(:,3) = rwg(:,2)
2611 : :
2612 [ # # ]: 0 : do iw=1,nw
2613 : : rwg(iw, 4) = &
2614 : : -E3(iw)**2 * inv_e31**2 * inv_e30 &
2615 : : * (one + (two * E1(iw) * inv_e31 + E0(iw) * inv_e30) * log(abs(E3(iw)))) &
2616 : : - E1(iw)**2 * inv_e31**2 * inv_e10 &
2617 : : * (one + (-two * E3(iw) * inv_e31 + E0(iw) * inv_e10) * log(abs(E1(iw)))) &
2618 : 0 : + E0(iw)**3 * inv_e30**2 * inv_e10**2 * log(abs(E0(iw)))
2619 : : end do
2620 : :
2621 : : !rwg = zero
2622 : :
2623 : : ! e1=e2
2624 [ # # ]: 0 : else if (D12 < tol) then
2625 : :
2626 : 0 : e20 = eig(3) - eig(1); inv_e20 = one / e20
2627 : 0 : e30 = eig(4) - eig(1); inv_e30 = one / e30
2628 : 0 : e32 = eig(4) - eig(3); inv_e32 = one / e32
2629 : :
2630 [ # # ]: 0 : do iw=1,nw
2631 : : rwg(iw, 1) = &
2632 : : -E2(iw)**3 * inv_e32 * inv_e20**3 * log(abs(E2(iw))) &
2633 : : + E3(iw)**3 * inv_e32 * inv_e30**3 * log(abs(E3(iw))) &
2634 : : + E0(iw) * inv_e20 * inv_e30 &
2635 : : * (half - E2(iw) * inv_e20 - E3(iw) * inv_e30 + &
2636 : 0 : (E2(iw)**2 * inv_e20**2 + E3(iw)**2 * inv_e30**2 + E2(iw) * E3(iw) * inv_e20 * inv_e30) * log(abs(E0(iw))))
2637 : : end do
2638 [ # # ]: 0 : rwg(:,2) = rwg(:,1)
2639 : :
2640 [ # # ]: 0 : do iw=1,nw
2641 : : rwg(iw, 3) = &
2642 : : E2(iw)**2 * inv_e20**2 * inv_e32 &
2643 : : * (one + (two * E0(iw) * inv_e20 - E3(iw) * inv_e32) * log(abs(E2(iw)))) &
2644 : : + E0(iw)**2 * inv_e20**2 * inv_e30 &
2645 : : * (one - (two * E2(iw) * inv_e20 + E3(iw) * inv_e30) * log(abs(E0(iw)))) &
2646 : 0 : + (E3(iw)**3 * inv_e32**2 * inv_e30**2 * log(abs(E3(iw))))
2647 : : end do
2648 : : ! This was wrong due to a misplaced parantes
2649 : : !rwg(:,3) = zero
2650 : :
2651 [ # # ]: 0 : do iw=1,nw
2652 : : rwg(iw, 4) = &
2653 : : -E3(iw)**2 * inv_e30**2 * inv_e32 &
2654 : : * (one + (two * E0(iw) * inv_e30 + E2(iw) * inv_e32) * log(abs(E3(iw)))) &
2655 : : + E0(iw)**2 * inv_e30**2 * inv_e20 &
2656 : : * (one - (two * E3(iw) * inv_e30 + E2(iw) * inv_e20) * log(abs(E0(iw)))) &
2657 : 0 : + (E2(iw)**3 * inv_e32**2 * inv_e20**2 * log(abs(E2(iw))))
2658 : : end do
2659 : :
2660 : : !rwg = zero
2661 : :
2662 : : ! e1<e2<e3<e4
2663 : : else
2664 : :
2665 : 0 : e10 = eig(2) - eig(1); inv_e10 = one / e10
2666 : 0 : e20 = eig(3) - eig(1); inv_e20 = one / e20
2667 : 0 : e21 = eig(3) - eig(2); inv_e21 = one / e21
2668 : 0 : e30 = eig(4) - eig(1); inv_e30 = one / e30
2669 : 0 : e31 = eig(4) - eig(2); inv_e31 = one / e31
2670 : 0 : e32 = eig(4) - eig(3); inv_e32 = one / e32
2671 : :
2672 [ # # ]: 0 : do iw=1,nw
2673 : : rwg(iw, 1) = &
2674 : : E0(iw)**2 * inv_e10 * inv_e20 * inv_e30 &
2675 : : * (one - (E1(iw) * inv_e10 + E2(iw) * inv_e20 + E3(iw) * inv_e30) * log(abs(E0(iw)))) &
2676 : : + E1(iw)**3 * inv_e10**2 * inv_e21 * inv_e31 * log(abs(E1(iw))) &
2677 : : - E2(iw)**3 * inv_e20**2 * inv_e21 * inv_e32 * log(abs(E2(iw))) &
2678 : 0 : + E3(iw)**3 * inv_e30**2 * inv_e31 * inv_e32 * log(abs(E3(iw)))
2679 : : end do
2680 : :
2681 [ # # ]: 0 : do iw=1,nw
2682 : : rwg(iw, 2) = &
2683 : : -E1(iw)**2 * inv_e10 * inv_e21 * inv_e31 &
2684 : : * (one + (E0(iw) * inv_e10 - E2(iw) * inv_e21 - E3(iw) * inv_e31) * log(abs(E1(iw)))) &
2685 : : + E0(iw)**3 * inv_e10**2 * inv_e20 * inv_e30 * log(abs(E0(iw))) &
2686 : : - E2(iw)**3 * inv_e20 * inv_e21**2 * inv_e32 * log(abs(E2(iw))) &
2687 : 0 : + E3(iw)**3 * inv_e30 * inv_e31**2 * inv_e32 * log(abs(E3(iw)))
2688 : : end do
2689 : :
2690 [ # # ]: 0 : do iw=1,nw
2691 : : rwg(iw, 3) = &
2692 : : E2(iw)**2 * inv_e20 * inv_e21 * inv_e32 &
2693 : : * (one + (E0(iw) * inv_e20 + E1(iw) * inv_e21 - E3(iw) * inv_e32) * log(abs(E2(iw)))) &
2694 : : + E0(iw)**3 * inv_e10 * inv_e20**2 * inv_e30 * log(abs(E0(iw))) &
2695 : : - E1(iw)**3 * inv_e10 * inv_e21**2 * inv_e31 * log(abs(E1(iw))) &
2696 : 0 : + E3(iw)**3 * inv_e30 * inv_e31 * inv_e32**2 * log(abs(E3(iw)))
2697 : : end do
2698 : :
2699 [ # # ]: 0 : do iw=1,nw
2700 : : rwg(iw, 4) = &
2701 : : -E3(iw)**2 * inv_e30 * inv_e31 * inv_e32 &
2702 : : * (one + (E0(iw) * inv_e30 + E1(iw) * inv_e31 + E2(iw) * inv_e32) * log(abs(E3(iw)))) &
2703 : : + E0(iw)**3 * inv_e10 * inv_e20 * inv_e30**2 * log(abs(E0(iw))) &
2704 : : - E1(iw)**3 * inv_e10 * inv_e21 * inv_e31**2 * log(abs(E1(iw))) &
2705 : 0 : + E2(iw)**3 * inv_e20 * inv_e21 * inv_e32**2 * log(abs(E2(iw)))
2706 : : end do
2707 : : !rwg = zero
2708 : :
2709 : : end if
2710 : :
2711 : 0 : end subroutine get_onetetra_ppart_lv
2712 : : !!***
2713 : :
2714 [ - + - - : 729 : end module m_htetra
+ + - + +
- - + - +
+ - - + -
+ - - - -
- + - + -
- - - -
+ ]
2715 : : !!***
|
