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 : !!***
|
