Branch data Line data Source code
1 : : !!****m* ABINIT/m_tetrahedron
2 : : !! NAME
3 : : !! m_tetrahedron
4 : : !!
5 : : !! FUNCTION
6 : : !! module for tetrahedron interpolation of DOS and similar quantities
7 : : !! depends on sort_tetra and on m_kpt_rank
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2010-2026 ABINIT group (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 carefully the case of degenerate tetrahedron
17 : : !! 2) Change API so that we can pass the energy mesh instead of omega_min and omega_max
18 : : !! 3) Add table ik_ibz --> tetra_list to avoid cycling inside big loop over ntetra
19 : : !! 4) Add options to get only delta and/or theta ?
20 : : !!
21 : : !! SOURCE
22 : :
23 : : #if defined HAVE_CONFIG_H
24 : : #include "config.h"
25 : : #endif
26 : :
27 : : #include "libtetra.h"
28 : :
29 : : module m_tetrahedron
30 : :
31 : : ! make sure stdout is defined, as libtetra.h needs it
32 : : use, intrinsic :: iso_fortran_env, only : stdin=>input_unit, &
33 : : stdout=>output_unit, &
34 : : stderr=>error_unit
35 : : USE_MEMORY_PROFILING
36 : : USE_MSG_HANDLING
37 : : use m_krank
38 : : USE_MPI
39 : : #ifdef HAVE_LIBTETRA_ABINIT
40 : : use m_io_tools, only : open_file
41 : : use m_xmpi
42 : : use defs_basis, only : dp, zero, one, tol14, pi
43 : : #elif
44 : : integer, parameter :: dp = kind(1.0d0)
45 : : real(dp), parameter :: tol14 = 1.d-14, zero = 0.d0, one = 1.d0
46 : : real(dp), parameter :: pi=3.141592653589793238462643383279502884197_dp
47 : : #endif
48 : :
49 : : implicit none
50 : :
51 : : #if defined HAVE_MPI1
52 : : include 'mpif.h'
53 : : #endif
54 : :
55 : : private
56 : : !!***
57 : :
58 : : !!****t* m_tetrahedron/t_tetrahedron
59 : : !! NAME
60 : : !! t_tetrahedron
61 : : !!
62 : : !! FUNCTION
63 : : !! tetrahedron geometry object
64 : : !!
65 : : !! SOURCE
66 : :
67 : : type, public :: t_tetrahedron
68 : :
69 : : integer :: ntetra = 0
70 : : ! Number of tetrahedra
71 : :
72 : : real(dp) :: vv
73 : : ! volume of the tetrahedra
74 : :
75 : : real(dp) :: klatt(3, 3)
76 : : ! reciprocal of lattice vectors for full kpoint grid
77 : :
78 : : integer,allocatable :: tetra_full(:,:,:)
79 : : !(4,2,ntetra)
80 : : ! For each tetra
81 : : ! (:,1,itetra) indices of the vertex in IBZ (symmetrical image)
82 : : ! (:,2,itetra) indices of the vertexes in the BZ
83 : :
84 : : integer,allocatable :: tetra_mult(:)
85 : : !(ntetra)
86 : : ! multiplicity of each irred tetrahedron
87 : :
88 : : integer,allocatable :: tetra_wrap(:,:,:)
89 : : !(3,4,ntetra)
90 : : ! flag to wrap tetrahedron summit into IBZ
91 : :
92 : : integer,allocatable :: ibz_tetra_count(:)
93 : : ! ibz_tetra_mapping(nkpt_ibz)
94 : : ! Number of tetrahedra associated to a point in the IBZ.
95 : :
96 : : integer,allocatable :: ibz_tetra_mapping(:,:)
97 : : ! ibz_tetra_mapping(nkpt_ibz, maxval(tetra%ibz_tetra_count)))
98 : : ! map ikbz to tetra index.
99 : :
100 : : end type t_tetrahedron
101 : :
102 : : public :: init_tetra ! Initialize the object
103 : : ! See also the high-level interface tetra_from_kptrlatt provided by m_kpts.
104 : : public :: get_tetra_weight ! Calculate integration weights and their derivatives. shape (nkpt, nene).
105 : : public :: tetra_blochl_weights ! Same as in get_tetra_weight but weights have shape (nene, nkpt).
106 : : public :: get_dbl_tetra_weight ! Calculate integration weights for double tetrahedron integration of delta functions.
107 : : ! (NB these correspond to the derivative terms in normal tetrahedron).
108 : : public :: destroy_tetra ! Free memory.
109 : : public :: tetra_write ! Write text file (XML format) with tetra info.
110 : : public :: tetralib_has_mpi ! Return True if the library has been compiled with MPI support.
111 : : public :: tetra_get_onewk ! Calculate integration weights and their derivatives for a single k-point in the IBZ.
112 : : public :: tetra_get_onewk_wvals ! Similar to tetra_get_onewk_wvalsa but receives arbitrary list of frequency points.
113 : : public :: tetra_get_onetetra_wvals ! Get weights for one tetrahedra with arbitrary list of frequency points
114 : : !!***
115 : :
116 : : contains
117 : : !!***
118 : :
119 : : !----------------------------------------------------------------------
120 : :
121 : : !!****f* m_tetrahedron/destroy_tetra
122 : : !! NAME
123 : : !! destroy_tetra
124 : : !!
125 : : !! FUNCTION
126 : : !! deallocate tetrahedra pointers if needed
127 : : !!
128 : : !! SOURCE
129 : :
130 : 11 : subroutine destroy_tetra (tetra)
131 : :
132 : : type(t_tetrahedron), intent(inout) :: tetra
133 : :
134 [ + + ]: 11 : if (allocated(tetra%tetra_full)) then
135 : 4 : TETRA_DEALLOCATE(tetra%tetra_full)
136 : : end if
137 [ + + ]: 11 : if (allocated(tetra%tetra_mult)) then
138 : 4 : TETRA_DEALLOCATE(tetra%tetra_mult)
139 : : end if
140 [ + + ]: 11 : if (allocated(tetra%tetra_wrap)) then
141 : 4 : TETRA_DEALLOCATE(tetra%tetra_wrap)
142 : : end if
143 [ + + ]: 11 : if (allocated(tetra%ibz_tetra_count)) then
144 : 4 : TETRA_DEALLOCATE(tetra%ibz_tetra_count)
145 : : end if
146 [ + + ]: 11 : if (allocated(tetra%ibz_tetra_mapping)) then
147 : 4 : TETRA_DEALLOCATE(tetra%ibz_tetra_mapping)
148 : : end if
149 : :
150 : 11 : end subroutine destroy_tetra
151 : : !!***
152 : :
153 : : !----------------------------------------------------------------------
154 : :
155 : : !!****f* m_tetrahedron/init_tetra
156 : : !! NAME
157 : : !! init_tetra
158 : : !!
159 : : !! FUNCTION
160 : : !! get tetrahedra characterized by apexes
161 : : !!
162 : : !! INPUTS
163 : : !! indkpt(nkpt_fullbz)=indexes of irred kpoints equivalent to kpt_fullbz
164 : : !! gprimd(3,3) = reciprocal space vectors
165 : : !! klatt(3,3)=reciprocal of lattice vectors for full kpoint grid
166 : : !! kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone
167 : : !! nkpt_fullbz=number of kpoints in full brillouin zone
168 : : !! comm= MPI communicator
169 : : !!
170 : : !! OUTPUT
171 : : !! tetra%tetra_full(4,2,ntetra)=for each tetrahedron,
172 : : !! the different instances of the tetrahedron (fullbz kpoints)
173 : : !! tetra%tetra_mult(ntetra) = store multiplicity of each irred tetrahedron
174 : : !! tetra%tetra_wrap(3,4,ntetra) = store flag to wrap tetrahedron summit into IBZ
175 : : !! tetra%ntetra = final number of irred tetra (dimensions of tetra_* remain larger)
176 : : !! tetra%vv = tetrahedron volume divided by full BZ volume
177 : : !!
178 : : !! SOURCE
179 : :
180 [ - + - + : 4 : subroutine init_tetra(indkpt, gprimd, klatt, kpt_fullbz, nkpt_fullbz, tetra, ierr, errorstring, comm)
- + - + ]
181 : :
182 : : !Arguments ------------------------------------
183 : : !scalars
184 : : integer,intent(in) :: nkpt_fullbz, comm
185 : : integer, intent(out) :: ierr
186 : : character(len=80), intent(out) :: errorstring
187 : : type(t_tetrahedron),intent(out) :: tetra
188 : : !arrays
189 : : integer,intent(in) :: indkpt(nkpt_fullbz)
190 : : real(dp) ,intent(in) :: gprimd(3,3),klatt(3,3),kpt_fullbz(3,nkpt_fullbz)
191 : :
192 : : !Local variables-------------------------------
193 : : !scalars
194 : : integer :: ialltetra,ikpt2,ikpt_full,isummit,itetra,jalltetra,jsummit
195 : : integer :: ii,jj,ikibz,nkpt_ibz, my_rank, nprocs
196 : : integer :: symrankkpt,mtetra,itmp,ntetra_irred
197 : : real(dp) :: shift1,shift2,shift3, rcvol,hashfactor
198 : : !real :: cpu_start, cpu_stop
199 : 4 : type(krank_t) :: krank
200 : : !arrays
201 : : integer :: ind_ibz(4), tetra_shifts(3,4,6) ! 3 dimensions, 4 summits, and 6 tetrahedra / kpoint box
202 : : real(dp) :: k1(3),k2(3),k3(3)
203 : 4 : integer,allocatable :: tetra_full_(:,:,:)
204 : 4 : integer,allocatable :: tetra_mult_(:)
205 : 4 : integer,allocatable :: tetra_wrap_(:,:,:)
206 : 4 : integer, allocatable :: reforder(:)
207 : 4 : integer, allocatable :: irred_itetra(:)
208 [ - + ]: 4 : real(dp), allocatable :: tetra_hash(:)
209 : :
210 : : ! *********************************************************************
211 : :
212 : : !call cpu_time(cpu_start)
213 : :
214 : 4 : my_rank = 0; nprocs = 1
215 : : #ifdef HAVE_LIBTETRA_ABINIT
216 : 4 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
217 : : #endif
218 : :
219 : 4 : ierr = 0
220 : 4 : errorstring = ""
221 : : !jmb
222 : 4 : shift1 = zero
223 : 4 : shift2 = zero
224 : 4 : shift3 = zero
225 : :
226 [ + + + + ]: 52 : tetra%klatt = klatt
227 : :
228 : 4 : mtetra = 6 * nkpt_fullbz
229 [ - + + - : 12 : TETRA_ALLOCATE(tetra_full_, (4,2,mtetra))
+ - + - -
+ - + ]
230 [ - + + - : 12 : TETRA_ALLOCATE(tetra_mult_, (mtetra))
+ - + - -
+ - + ]
231 [ - + + - : 12 : TETRA_ALLOCATE(tetra_wrap_, (3,4,mtetra))
+ - + - -
+ - + ]
232 : :
233 [ + + ]: 5188 : tetra_mult_ = 1
234 [ + + + + : 57028 : tetra_full_ = 0
+ + ]
235 [ + + + + : 88132 : tetra_wrap_ = 0
+ + ]
236 : :
237 : : ! tetra_shifts(:,1,1) = (/0,0,0/)
238 : : ! tetra_shifts(:,2,1) = (/0,1,0/)
239 : : ! tetra_shifts(:,3,1) = (/0,1,1/)
240 : : ! tetra_shifts(:,4,1) = (/1,1,0/)
241 : : ! tetra_shifts(:,1,2) = (/0,0,0/)
242 : : ! tetra_shifts(:,2,2) = (/0,1,1/)
243 : : ! tetra_shifts(:,3,2) = (/1,1,0/)
244 : : ! tetra_shifts(:,4,2) = (/1,1,1/)
245 : : ! tetra_shifts(:,1,3) = (/0,0,0/)
246 : : ! tetra_shifts(:,2,3) = (/1,0,0/)
247 : : ! tetra_shifts(:,3,3) = (/1,1,0/)
248 : : ! tetra_shifts(:,4,3) = (/1,1,1/)
249 : : ! tetra_shifts(:,1,4) = (/0,0,0/)
250 : : ! tetra_shifts(:,2,4) = (/0,0,1/)
251 : : ! tetra_shifts(:,3,4) = (/1,0,0/)
252 : : ! tetra_shifts(:,4,4) = (/1,1,1/)
253 : : ! tetra_shifts(:,1,5) = (/0,0,1/)
254 : : ! tetra_shifts(:,2,5) = (/1,0,0/)
255 : : ! tetra_shifts(:,3,5) = (/1,0,1/)
256 : : ! tetra_shifts(:,4,5) = (/1,1,1/)
257 : : ! tetra_shifts(:,1,6) = (/0,0,0/)
258 : : ! tetra_shifts(:,2,6) = (/0,0,1/)
259 : : ! tetra_shifts(:,3,6) = (/0,1,1/)
260 : : ! tetra_shifts(:,4,6) = (/1,1,1/)
261 : :
262 : : ! bxu, the following division scheme is according to Bloechl's paper
263 [ + + ]: 16 : tetra_shifts(:,1,1) = (/0,0,0/)
264 [ + + ]: 16 : tetra_shifts(:,2,1) = (/1,0,0/)
265 [ + + ]: 16 : tetra_shifts(:,3,1) = (/0,1,0/)
266 [ + + ]: 16 : tetra_shifts(:,4,1) = (/1,0,1/)
267 [ + + ]: 16 : tetra_shifts(:,1,2) = (/1,0,0/)
268 [ + + ]: 16 : tetra_shifts(:,2,2) = (/1,1,0/)
269 [ + + ]: 16 : tetra_shifts(:,3,2) = (/0,1,0/)
270 [ + + ]: 16 : tetra_shifts(:,4,2) = (/1,0,1/)
271 [ + + ]: 16 : tetra_shifts(:,1,3) = (/0,1,0/)
272 [ + + ]: 16 : tetra_shifts(:,2,3) = (/1,1,0/)
273 [ + + ]: 16 : tetra_shifts(:,3,3) = (/1,0,1/)
274 [ + + ]: 16 : tetra_shifts(:,4,3) = (/1,1,1/)
275 [ + + ]: 16 : tetra_shifts(:,1,4) = (/0,0,0/)
276 [ + + ]: 16 : tetra_shifts(:,2,4) = (/0,1,0/)
277 [ + + ]: 16 : tetra_shifts(:,3,4) = (/0,0,1/)
278 [ + + ]: 16 : tetra_shifts(:,4,4) = (/1,0,1/)
279 [ + + ]: 16 : tetra_shifts(:,1,5) = (/0,0,1/)
280 [ + + ]: 16 : tetra_shifts(:,2,5) = (/1,0,1/)
281 [ + + ]: 16 : tetra_shifts(:,3,5) = (/0,1,0/)
282 [ + + ]: 16 : tetra_shifts(:,4,5) = (/0,1,1/)
283 [ + + ]: 16 : tetra_shifts(:,1,6) = (/0,1,0/)
284 [ + + ]: 16 : tetra_shifts(:,2,6) = (/1,0,1/)
285 [ + + ]: 16 : tetra_shifts(:,3,6) = (/0,1,1/)
286 [ + + ]: 16 : tetra_shifts(:,4,6) = (/1,1,1/)
287 : :
288 : : ! Make full k-point rank arrays
289 : : ! TODO: Lot of memory allocated here if dense mesh e.g ~ 300 ** 3
290 : 4 : call krank%init(nkpt_fullbz, kpt_fullbz)
291 : :
292 : 4 : ialltetra = 1
293 [ + + ]: 868 : do ikpt_full=1,nkpt_fullbz
294 [ + + ]: 6052 : do itetra=1,6
295 : : !ialltetra = itetra + (ikpt_full -1) * 6
296 : : !if (mod(ialltetra, nprocs) /= my_rank) cycle ! MPI parallelism.
297 [ + + ]: 25920 : do isummit=1,4
298 : : k1(:) = kpt_fullbz(:,ikpt_full) &
299 : : + tetra_shifts(1,isummit,itetra)*klatt(:,1) &
300 : : + tetra_shifts(2,isummit,itetra)*klatt(:,2) &
301 [ + + ]: 82944 : + tetra_shifts(3,isummit,itetra)*klatt(:,3)
302 : :
303 : : ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full !
304 : 20736 : symrankkpt = krank%get_rank(k1)
305 : 20736 : ikpt2 = krank%invrank(symrankkpt)
306 [ - + ]: 20736 : if (ikpt2 < 1) then
307 : 0 : errorstring = 'Error in ranking k-points - exiting with un-initialized tetrahedra.'
308 : 0 : ierr = 2
309 : 0 : call krank%free()
310 [ # # # # ]: 0 : TETRA_ALLOCATE(tetra%tetra_full, (4,2,1))
311 [ # # # # ]: 0 : TETRA_ALLOCATE(tetra%tetra_mult, (1))
312 [ # # # # ]: 0 : TETRA_ALLOCATE(tetra%tetra_wrap, (3,4,1))
313 [ # # ]: 0 : TETRA_DEALLOCATE(tetra_full_)
314 [ # # ]: 0 : TETRA_DEALLOCATE(tetra_mult_)
315 [ # # ]: 0 : TETRA_DEALLOCATE(tetra_wrap_)
316 : 0 : return
317 : : end if
318 : :
319 : : ! Store irreducible kpoint equivalent to kpt_fullbz(:,ikpt2)
320 : 20736 : tetra_full_(isummit,1,ialltetra) = indkpt(ikpt2)
321 : 20736 : tetra_full_(isummit,2,ialltetra) = ikpt2
322 : 20736 : shift1 = k1(1)-kpt_fullbz(1,ikpt2)
323 : 20736 : shift2 = k1(2)-kpt_fullbz(2,ikpt2)
324 : 20736 : shift3 = k1(3)-kpt_fullbz(3,ikpt2)
325 [ + + ]: 20736 : if (shift1>0.5d0) then
326 : 1728 : tetra_wrap_(1,isummit,ialltetra) = 1
327 [ - + ]: 19008 : else if (shift1<-0.5d0) then
328 : 0 : tetra_wrap_(1,isummit,ialltetra) = -1
329 : : end if
330 [ + + ]: 20736 : if (shift2>0.5d0) then
331 : 1728 : tetra_wrap_(2,isummit,ialltetra) = 1
332 [ - + ]: 19008 : else if (shift2<-0.5d0) then
333 : 0 : tetra_wrap_(2,isummit,ialltetra) = -1
334 : : end if
335 [ + + ]: 20736 : if (shift3>0.5d0) then
336 : 1728 : tetra_wrap_(3,isummit,ialltetra) = 1
337 [ + - ]: 19008 : else if (shift3<-0.5d0) then
338 : 0 : tetra_wrap_(3,isummit,ialltetra) = -1
339 : : end if
340 : :
341 : : ! sort itetra summits
342 : : ! TODO: replace with sort_int
343 [ + + ]: 57024 : do jsummit=isummit,2,-1
344 [ + + ]: 51840 : if ( tetra_full_(jsummit,1,ialltetra) < tetra_full_(jsummit-1,1,ialltetra) ) then
345 : 6912 : itmp = tetra_full_(jsummit,1,ialltetra)
346 : 6912 : tetra_full_(jsummit,1,ialltetra) = tetra_full_(jsummit-1,1,ialltetra)
347 : 6912 : tetra_full_(jsummit-1,1,ialltetra) = itmp
348 : 6912 : itmp = tetra_full_(jsummit,2,ialltetra)
349 : 6912 : tetra_full_(jsummit,2,ialltetra) = tetra_full_(jsummit-1,2,ialltetra)
350 : 6912 : tetra_full_(jsummit-1,2,ialltetra) = itmp
351 : : ! keep fullbz_kpt tetrahedra points in same order
352 : 6912 : itmp = tetra_wrap_(1,jsummit,ialltetra)
353 : 6912 : tetra_wrap_(1,jsummit,ialltetra) = tetra_wrap_(1,jsummit-1,ialltetra)
354 : 6912 : tetra_wrap_(1,jsummit-1,ialltetra) = itmp
355 : 6912 : itmp = tetra_wrap_(2,jsummit,ialltetra)
356 : 6912 : tetra_wrap_(2,jsummit,ialltetra) = tetra_wrap_(2,jsummit-1,ialltetra)
357 : 6912 : tetra_wrap_(2,jsummit-1,ialltetra) = itmp
358 : 6912 : itmp = tetra_wrap_(1,jsummit,ialltetra)
359 : 6912 : tetra_wrap_(3,jsummit,ialltetra) = tetra_wrap_(3,jsummit-1,ialltetra)
360 : 6912 : tetra_wrap_(3,jsummit-1,ialltetra) = itmp
361 : : end if
362 : : end do ! jsummit
363 : :
364 : : end do ! isummit
365 : :
366 [ - + ]: 5184 : if (ialltetra > mtetra) then
367 : : write (errorstring, '(3a,i0,a,i0)' ) &
368 : 0 : 'init_tetra: BUG - ',&
369 : 0 : ' ialltetra > mtetra ',&
370 : 0 : ' ialltetra= ',ialltetra,', mtetra= ',mtetra
371 : 0 : ierr = 1
372 : 0 : return
373 : : end if
374 : 6048 : ialltetra = ialltetra+1
375 : : end do ! itetra
376 : : end do ! ikpt_full
377 : :
378 : : !call cpu_time(cpu_stop)
379 : : !write(*,*)"tetra_init ikpt_loop:", cpu_stop - cpu_start
380 : : !cpu_start = cpu_stop
381 : :
382 : 4 : call krank%free()
383 : :
384 : : rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
385 : : & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
386 : 4 : & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
387 : :
388 : : ! Volume of all tetrahedra should be the same as that of tetra 1
389 : : ! this is the volume of 1 tetrahedron, should be coherent with notation in Lehmann & Taut
390 [ + + ]: 16 : k1(:) = gprimd(:,1)*klatt(1,1) + gprimd(:,2)*klatt(2,1) + gprimd(:,3)*klatt(3,1)
391 [ + + ]: 16 : k2(:) = gprimd(:,1)*klatt(1,2) + gprimd(:,2)*klatt(2,2) + gprimd(:,3)*klatt(3,2)
392 [ + + ]: 16 : k3(:) = gprimd(:,1)*klatt(1,3) + gprimd(:,2)*klatt(2,3) + gprimd(:,3)*klatt(3,3)
393 : : tetra%vv = abs (k1(1)*(k2(2)*k3(3)-k2(3)*k3(2)) &
394 : : & -k1(2)*(k2(1)*k3(3)-k2(3)*k3(1)) &
395 : 4 : & +k1(3)*(k2(1)*k3(2)-k2(2)*k3(1))) / 6.d0 / rcvol
396 : :
397 : : ! eliminate equivalent tetrahedra by symmetry and account for them in multiplicity tetra_mult
398 : 4 : tetra%ntetra = mtetra
399 : :
400 : : ! FIXME: could we replace this with a ranking algorithm to avoid the O(tetra%ntetra^2) step? For example:
401 : : ! get tetrahedron rank - problem too many combinations in principle = nkpt_irred^4 - only a few used in practice
402 : : ! sort ranks and keep indices
403 : :
404 : : ! make hash table = tetra_full_(1)*nkptirred**3+tetra_full_(2)*nkptirred**2+tetra_full_(3)*nkptirred**1+tetra_full_(4)
405 : :
406 : 4 : hashfactor = 100.d0 ! *acos(-1.d0) ! 100 pi should be far from an integer...
407 [ - + + - : 12 : TETRA_ALLOCATE(tetra_hash, (tetra%ntetra))
+ - + - -
+ - + ]
408 [ - + + - : 12 : TETRA_ALLOCATE(reforder, (tetra%ntetra))
+ - + - -
+ - + ]
409 : :
410 : : !MG: In principle the order of the indices should not matter.
411 [ + + ]: 5188 : do ialltetra=1, tetra%ntetra
412 : : tetra_hash(ialltetra) = tetra_full_(1,1,ialltetra)*hashfactor**3+&
413 : : & tetra_full_(2,1,ialltetra)*hashfactor**2+&
414 : : & tetra_full_(3,1,ialltetra)*hashfactor**1+&
415 : 5184 : & tetra_full_(4,1,ialltetra)
416 : 5188 : reforder(ialltetra) = ialltetra
417 : : end do
418 : :
419 : 4 : call sort_tetra(tetra%ntetra, tetra_hash, reforder, tol14)
420 : : ! Most of the wall-time is spent in the preamble of this routine (up to this point).
421 : : ! sort_tetra is not easy to parallelize...
422 : :
423 : : ! determine number of tetra after reduction
424 [ - + + - : 12 : TETRA_ALLOCATE(irred_itetra, (tetra%ntetra))
+ - + - -
+ - + ]
425 : 4 : jalltetra = 1
426 : 4 : irred_itetra(1) = 1
427 [ + + ]: 5184 : do ialltetra=2, tetra%ntetra
428 [ + - ]: 5180 : if (abs(tetra_hash(ialltetra)-tetra_hash(ialltetra-1)) > tol14) then
429 : : ! found a new series
430 : 5180 : jalltetra = jalltetra + 1
431 : : end if
432 : 5184 : irred_itetra(ialltetra) = jalltetra
433 : : end do
434 : :
435 : : ! reset number of tetra
436 : 4 : ntetra_irred = jalltetra
437 : :
438 : : ! allocate definitive tetra arrays and transfer to new arrays
439 [ - + + - : 12 : TETRA_ALLOCATE(tetra%tetra_full, (4,2,ntetra_irred))
+ - - + -
+ - + ]
440 [ - + + - : 12 : TETRA_ALLOCATE(tetra%tetra_mult, (ntetra_irred))
+ - - + -
+ - + ]
441 [ - + + - : 12 : TETRA_ALLOCATE(tetra%tetra_wrap, (3,4,ntetra_irred))
+ - - + -
+ - + ]
442 : :
443 : : ! eliminate equal rank tetrahedra and accumulate multiplicity into first one
444 [ + + + + : 57028 : tetra%tetra_full = 0
+ + ]
445 [ + + ]: 5188 : tetra%tetra_mult = 0
446 [ + + + + : 88132 : tetra%tetra_wrap = 0
+ + ]
447 : 4 : jalltetra = 1
448 [ + + + + ]: 44 : tetra%tetra_full(:,:,1) = tetra_full_(:,:,reforder(1))
449 : 4 : tetra%tetra_mult(1) = 1
450 [ + + + + ]: 68 : tetra%tetra_wrap(:,:,1) = tetra_wrap_(:,:,reforder(1))
451 [ + + ]: 5184 : do ialltetra=2, tetra%ntetra
452 : : ! TODO: check if tolerance is adapted
453 [ + - ]: 5184 : if (abs(tetra_hash(ialltetra)-tetra_hash(ialltetra-1)) > tol14) then
454 : : ! found a new series
455 : 5180 : jalltetra = jalltetra + 1
456 [ + + + + ]: 56980 : tetra%tetra_full(:,:,jalltetra) = tetra_full_(:,:,reforder(ialltetra))
457 [ + + + + ]: 88060 : tetra%tetra_wrap(:,:,jalltetra) = tetra_wrap_(:,:,reforder(ialltetra))
458 : 5180 : tetra%tetra_mult(jalltetra) = 1
459 : : else
460 : : ! TODO: add real check that the tetra are equivalent...
461 : : ! otherwise increment jalltetra here as well, generate new series?
462 : 0 : tetra%tetra_mult(jalltetra) = tetra%tetra_mult(jalltetra) + tetra_mult_(reforder(ialltetra))
463 : : !tetra_mult_(reforder(ialltetra)) = 0
464 : : end if
465 : : end do
466 : :
467 : : ! reset of ntetra for final version after checks and debu
468 : 4 : tetra%ntetra = ntetra_irred
469 : :
470 [ - + ]: 4 : TETRA_DEALLOCATE(tetra_hash)
471 : 4 : TETRA_DEALLOCATE(reforder)
472 : 4 : TETRA_DEALLOCATE(irred_itetra)
473 [ - + ]: 4 : TETRA_DEALLOCATE(tetra_full_)
474 [ - + ]: 4 : TETRA_DEALLOCATE(tetra_mult_)
475 [ - + ]: 4 : TETRA_DEALLOCATE(tetra_wrap_)
476 : :
477 : : ! Create mapping between the irreducible k-points
478 : : ! and all the tetrahedron contributing with some weight
479 [ + + ]: 868 : nkpt_ibz = maxval(indkpt)
480 : :
481 : : ! 1. First we count what is the maximum number of distinct tetrahedra that each k-point contains
482 [ - + + - : 12 : TETRA_ALLOCATE(tetra%ibz_tetra_count,(nkpt_ibz))
+ - + - -
+ - + -
+ ]
483 [ + + ]: 868 : tetra%ibz_tetra_count(:) = 0
484 : :
485 : : ! Count max tetra contributing
486 [ + + ]: 5188 : do ii=1,tetra%ntetra
487 : : ! Here we need the original ordering to reference the correct irred kpoints
488 [ + + ]: 25920 : ind_ibz(:) = tetra%tetra_full(:,1,ii)
489 : : ! count max tetra contributing
490 [ + + ]: 25924 : do jj=1,4
491 : 20736 : ikibz = ind_ibz(jj)
492 [ - + ]: 20736 : if (ikibz > nkpt_ibz) cycle
493 : 25920 : tetra%ibz_tetra_count(ikibz) = tetra%ibz_tetra_count(ikibz) + 1
494 : : end do
495 : : end do
496 : :
497 : : ! 2. Then we build mapping of ikbz to tetra
498 [ - + + - : 880 : TETRA_ALLOCATE(tetra%ibz_tetra_mapping,(nkpt_ibz,maxval(tetra%ibz_tetra_count)))
+ + - + +
- + - + -
- + - + -
+ ]
499 [ + + ]: 868 : tetra%ibz_tetra_count(:) = 0
500 [ + + ]: 5188 : do ii=1,tetra%ntetra
501 : : ! Here we need the original ordering to reference the correct irred kpoints
502 [ + + ]: 25920 : ind_ibz(:) = tetra%tetra_full(:,1,ii)
503 : : ! Use the counter to move pointer and then fill index
504 [ + + ]: 25924 : do jj=1,4
505 : 20736 : ikibz = ind_ibz(jj)
506 [ - + ]: 20736 : if (ikibz > nkpt_ibz) cycle
507 : : ! avoid putting the same index twice
508 [ + + ]: 20736 : if (tetra%ibz_tetra_count(ikibz) > 0) then
509 [ - + ]: 19872 : if (tetra%ibz_tetra_mapping(ikibz,tetra%ibz_tetra_count(ikibz)) == ii) cycle
510 : : end if
511 : 20736 : tetra%ibz_tetra_count(ikibz) = tetra%ibz_tetra_count(ikibz) + 1
512 : 25920 : tetra%ibz_tetra_mapping(ikibz,tetra%ibz_tetra_count(ikibz)) = ii
513 : : end do
514 : : end do
515 : :
516 : : !call cpu_time(cpu_stop)
517 : : !write(*,*)"tetra_init 2nd part:", cpu_stop - cpu_start
518 : : !cpu_start = cpu_stop
519 : :
520 [ - + - + : 4 : end subroutine init_tetra
- + - + -
+ - + ]
521 : : !!***
522 : :
523 : : !----------------------------------------------------------------------
524 : :
525 : : !!****f* m_tetrahedron/tetra_write
526 : : !! NAME
527 : : !! tetra_write
528 : : !!
529 : : !! FUNCTION
530 : : !! Write text file with tetra info.
531 : : !!
532 : : !! INPUTS
533 : : !! tetra<t_tetrahedron>=tetrahedron geometry object
534 : : !! nkibz=Number of k-points in the IBZ used to generate tetra
535 : : !! kibz(3,nkibz)=Reduced coordinates of the IBZ
536 : : !! path=Name of output file
537 : : !!
538 : : !! OUTPUT
539 : : !! Output is written to file.
540 : : !!
541 : : !! SOURCE
542 : :
543 : 0 : subroutine tetra_write(tetra, nkibz, kibz, path)
544 : :
545 : : !Arguments ------------------------------------
546 : : !scalars
547 : : integer,intent(in) :: nkibz
548 : : character(len=*),intent(in) :: path
549 : : type(t_tetrahedron),intent(in) :: tetra
550 : : !arrays
551 : : real(dp),intent(in) :: kibz(3,nkibz)
552 : :
553 : : !Local variables-------------------------------
554 : : integer,parameter :: version=1
555 : : integer :: ik,it,unt
556 : : #ifdef HAVE_LIBTETRA_ABINIT
557 : : character(len=500) :: msg
558 : : #endif
559 : :
560 : : ! *********************************************************************
561 : :
562 : : #ifdef HAVE_LIBTETRA_ABINIT
563 [ # # ]: 0 : if (open_file(file=trim(path),iomsg=msg,newunit=unt,form="formatted",status="unknown",action="write")/=0) then
564 : 0 : TETRA_ERROR(msg)
565 : : end if
566 : : #else
567 : : open(file=trim(path),newunit=unt,form="formatted",status="unknown",action="write")
568 : : #endif
569 : :
570 : 0 : write(unt,*)version, " # version number"
571 : :
572 : : ! Write IBZ
573 : 0 : write(unt,*)nkibz, " # number of k-points in the IBZ"
574 : 0 : write(unt,"(a)")"<irreducible_zone>"
575 [ # # ]: 0 : do ik=1,nkibz
576 : 0 : write(unt,"(3es22.12)") kibz(:,ik)
577 : : end do
578 : 0 : write(unt,"(a)")"</irreducible_zone>"
579 : :
580 : : ! Write tetra info
581 : 0 : write(unt,"(i0,a)")tetra%ntetra, " # number of tetrahedra"
582 : 0 : write(unt,"(es22.12,a)")tetra%vv, " # tetrahedron volume"
583 : :
584 : 0 : write(unt,"(a)")"<tetra_full>"
585 [ # # ]: 0 : do it=1,tetra%ntetra
586 : 0 : write(unt,"(8(i0,1x))")tetra%tetra_full(:,:,it)
587 : : end do
588 : 0 : write(unt,"(a)")"</tetra_full>"
589 : :
590 : 0 : write(unt,"(a)")"<tetra_mult>"
591 [ # # ]: 0 : do it=1,tetra%ntetra
592 : 0 : write(unt,"(i0)")tetra%tetra_mult(it)
593 : : end do
594 : 0 : write(unt,"(a)")"</tetra_mult>"
595 : :
596 : 0 : write(unt,"(a)")"<tetra_wrap>"
597 [ # # ]: 0 : do it=1,tetra%ntetra
598 : 0 : write(unt,"(12(i0,1x))")tetra%tetra_wrap(:,:,it)
599 : : end do
600 : 0 : write(unt,"(a)")"</tetra_wrap>"
601 : :
602 : 0 : close(unt)
603 : :
604 : 0 : end subroutine tetra_write
605 : : !!***
606 : :
607 : : !----------------------------------------------------------------------
608 : :
609 : : !!****f* m_tetrahedron/get_tetra_weight
610 : : !! NAME
611 : : !! get_tetra_weight
612 : : !!
613 : : !! FUNCTION
614 : : !! calculate integration weights and their derivatives from Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
615 : : !!
616 : : !! INPUTS
617 : : !! eigen_in(nkpt)=eigenenergies for each k point
618 : : !! enemin=minimal energy for DOS
619 : : !! enemax=maximal energy for DOS
620 : : !! max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)
621 : : !! nene=number of energies for DOS
622 : : !! nkpt=number of irreducible kpoints
623 : : !! tetra<t_tetrahedron>
624 : : !! %ntetra=number of tetrahedra
625 : : !! %tetra_full(4,2,ntetra)=for each irred tetrahedron, the list of k point vertices
626 : : !! 1 -> irred kpoint 2 -> fullkpt
627 : : !! %tetra_mult(ntetra)=for each irred tetrahedron, its multiplicity
628 : : !! %vv = ratio of volume of one tetrahedron in reciprocal space to full BZ volume
629 : : !! bcorr=1 to include Blochl correction else 0.
630 : : !! comm=MPI communicator
631 : : !!
632 : : !! OUTPUT
633 : : !! tweight(nkpt,nene) = integration weights for each irred kpoint from all adjacent tetrahedra
634 : : !! dtweightde(nkpt,nene) = derivative of tweight wrt energy
635 : : !!
636 : : !! SOURCE
637 : :
638 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
639 : : ! THIS FUNCTION IS DEPRECATED, USE tetra_blochl_weights
640 : : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
641 : 18 : subroutine get_tetra_weight(eigen_in,enemin,enemax,max_occ,nene,nkpt,tetra,&
642 : 18 : bcorr,tweight,dtweightde,comm)
643 : :
644 : : !Arguments ------------------------------------
645 : : !scalars
646 : : integer,intent(in) :: nene,nkpt,bcorr,comm
647 : : type(t_tetrahedron), intent(in) :: tetra
648 : : real(dp) ,intent(in) :: enemax,enemin,max_occ
649 : : !arrays
650 : : real(dp) ,intent(in) :: eigen_in(nkpt)
651 : : real(dp) ,intent(out) :: dtweightde(nkpt,nene),tweight(nkpt,nene)
652 : :
653 : : !Local variables-------------------------------
654 : 18 : real(dp), allocatable :: dtweightde_ek(:, :), tweight_ek(:, :)
655 : :
656 : : ! *********************************************************************
657 : :
658 [ - + + - : 72 : TETRA_ALLOCATE(dtweightde_ek, (nene, nkpt))
- + + - +
- + - - +
- + ]
659 [ - + + - : 54 : TETRA_ALLOCATE(tweight_ek, (nene, nkpt))
- + + - -
+ - + ]
660 : :
661 : 18 : call tetra_blochl_weights(tetra,eigen_in,enemin,enemax,max_occ,nene,nkpt,bcorr,tweight_ek,dtweightde_ek,comm)
662 : :
663 : : ! transpose: otherwise the data access is crap and the code slows by an order of magnitude
664 [ + + + + ]: 782520 : tweight = transpose(tweight_ek)
665 [ + + + + ]: 782520 : dtweightde = transpose(dtweightde_ek)
666 : :
667 [ - + ]: 18 : TETRA_DEALLOCATE(dtweightde_ek)
668 [ - + ]: 18 : TETRA_DEALLOCATE(tweight_ek)
669 : :
670 : 18 : end subroutine get_tetra_weight
671 : : !!***
672 : :
673 : : !----------------------------------------------------------------------
674 : :
675 : : !!****f* m_tetrahedron/tetra_blochl_weights
676 : : !! NAME
677 : : !! tetra_blochl_weights
678 : : !!
679 : : !! FUNCTION
680 : : !! calculate integration weights and their derivatives from Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
681 : : !! Same API as get_tetra_weight but output weights here have shape (nene, nkpt)
682 : : !!
683 : : !! SOURCE
684 : :
685 : 18 : subroutine tetra_blochl_weights(tetra,eigen_in,enemin,enemax,max_occ,nene,nkpt,&
686 : 18 : bcorr,tweight_t,dtweightde_t,comm)
687 : :
688 : : !Arguments ------------------------------------
689 : : !scalars
690 : : integer,intent(in) :: nene,nkpt,bcorr,comm
691 : : type(t_tetrahedron), intent(in) :: tetra
692 : : real(dp) ,intent(in) :: enemax,enemin,max_occ
693 : : !arrays
694 : : real(dp) ,intent(in) :: eigen_in(nkpt)
695 : : real(dp) ,intent(out) :: dtweightde_t(nene,nkpt),tweight_t(nene,nkpt)
696 : :
697 : : !Local variables-------------------------------
698 : : !scalars
699 : : integer :: itetra,nprocs,my_start,my_stop,ierr,ii
700 : : !arrays
701 : : integer :: ind_ibz(4)
702 : : real(dp) :: eigen_1tetra(4)
703 : 18 : real(dp), allocatable :: tweight_tmp(:,:),dtweightde_tmp(:,:),buffer(:,:)
704 : :
705 : : ! *********************************************************************
706 : :
707 [ - + + - : 54 : TETRA_ALLOCATE(tweight_tmp, (nene, 4))
+ - + - -
+ - + ]
708 [ - + + - : 36 : TETRA_ALLOCATE(dtweightde_tmp, (nene, 4))
- + + - ]
709 [ + + + + : 1565604 : tweight_t = zero; dtweightde_t = zero
+ + + + ]
710 : :
711 : 18 : call split_work(tetra%ntetra, comm, nprocs, my_start, my_stop, ierr)
712 [ - + ]: 18 : if (ierr /= 0) TETRA_ERROR("Error in MPI layer")
713 : :
714 : : ! for each tetrahedron
715 [ + + ]: 23346 : do itetra=my_start,my_stop
716 [ + + + + ]: 18810144 : tweight_tmp = zero
717 [ + + + + ]: 18810144 : dtweightde_tmp = zero
718 : :
719 : : ! Here we need the original ordering to reference the correct irred kpoints
720 [ + + ]: 116640 : ind_ibz(:) = tetra%tetra_full(:,1,itetra)
721 : :
722 : : ! Sort energies before calling get_onetetra_
723 [ + + ]: 116640 : eigen_1tetra(:) = eigen_in(ind_ibz(:))
724 : 23328 : call sort_tetra(4, eigen_1tetra, ind_ibz, tol14)
725 : :
726 : 23328 : call get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,nene,bcorr,tweight_tmp,dtweightde_tmp)
727 : :
728 : : ! NOTE: the following blas calls are not working systematically, or do not give speed ups, strange...
729 : : !if (nene > 100) then
730 : : ! do ii=1,4
731 : : ! call daxpy (nene, 1.d0, tweight_tmp(:,ii), 1, tweight_t(:,ind_ibz(ii)), 1)
732 : : ! end do
733 : : ! do ii=1,4
734 : : ! call daxpy (nene, 1.d0, dtweightde_tmp(:,ii), 1, dtweightde_t(:,ind_ibz(ii)), 1)
735 : : ! end do
736 : : !else
737 [ + + ]: 116640 : do ii=1,4
738 [ + + ]: 18810144 : tweight_t(:,ind_ibz(ii)) = tweight_t(:,ind_ibz(ii)) + tweight_tmp(:,ii)
739 : : end do
740 [ + + ]: 116658 : do ii=1,4
741 [ + + ]: 18810144 : dtweightde_t(:,ind_ibz(ii)) = dtweightde_t(:,ind_ibz(ii)) + dtweightde_tmp(:,ii)
742 : : end do
743 : : !end if
744 : : end do ! itetra
745 : :
746 [ - + ]: 18 : TETRA_DEALLOCATE(tweight_tmp)
747 [ - + ]: 18 : TETRA_DEALLOCATE(dtweightde_tmp)
748 : :
749 [ - + ]: 18 : if (nprocs > 1) then
750 : : #ifdef HAVE_MPI
751 [ # # # # : 0 : TETRA_ALLOCATE(buffer, (nene, nkpt))
# # # # #
# # # # #
# # ]
752 : 0 : call MPI_ALLREDUCE(tweight_t,buffer,nene*nkpt,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
753 [ # # # # ]: 0 : tweight_t = buffer
754 : :
755 : 0 : call MPI_ALLREDUCE(dtweightde_t,buffer,nene*nkpt,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
756 [ # # # # ]: 0 : dtweightde_t = buffer
757 [ # # ]: 0 : TETRA_DEALLOCATE(buffer)
758 : : #endif
759 : : end if
760 : :
761 : 18 : end subroutine tetra_blochl_weights
762 : : !!***
763 : :
764 : : !----------------------------------------------------------------------
765 : :
766 : : !!****f* m_tetrahedron/get_dbl_tetra_weight
767 : : !! NAME
768 : : !! get_dbl_tetra_weight
769 : : !!
770 : : !! FUNCTION
771 : : !! calculate integration weights and their derivatives
772 : : !! for double tetrahedron method from Allen Phys Stat Sol B 120 529 (1983) [[cite:Allen1983b]]
773 : : !! the k-points and tetrahedra must be the same for both grids, of course,
774 : : !! but the range of energies is arbitrary
775 : : !!
776 : : !! Omega is called eigen1 here
777 : : !! E is called eigen2 here
778 : : !! indexing goes from 1 to 4 for the tetrahedron corners, in order of increasing eigen1
779 : : !! in Allen, from 0 to 3...
780 : : !!
781 : : !! INPUTS
782 : : !! eigen1_in(nkpt)=eigenenergies for each k point
783 : : !! eigen2_in(nkpt)=eigenenergies for each k point
784 : : !! enemin1=minimal energy for DOS in energy 1
785 : : !! enemax1=maximal energy for DOS
786 : : !! enemin2=minimal energy for DOS in energy 2
787 : : !! enemax2=maximal energy for DOS
788 : : !! max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)
789 : : !! nene1=number of energies for DOS in energy 1
790 : : !! nene2=number of energies for DOS in energy 2
791 : : !! nkpt=number of irreducible kpoints
792 : : !! tetra%ntetra=number of tetra
793 : : !! tetra%tetra_full(4,2,ntetra)=for each irred tetrahedron, the list of k point vertices
794 : : !! 1 -> irred kpoint 2 -> fullkpt
795 : : !! tetra%tetra_mult(ntetra)=for each irred tetrahedron, its multiplicity
796 : : !! tetra%vv = ratio of volume of one tetrahedron in reciprocal space to full BZ volume
797 : : !! ierr = error code on exit
798 : : !!
799 : : !! OUTPUT
800 : : !! tweight(nkpt,nene1,nene2) = integration weights for each irred kpoint from all adjacent tetrahedra
801 : : !! dtweightde(nkpt,nene1,nene2) = derivative of tweight wrt energy
802 : : !!
803 : : !! SOURCE
804 : :
805 : 0 : subroutine get_dbl_tetra_weight(eigen1_in,eigen2_in,enemin1,enemax1,enemin2,enemax2,&
806 : 0 : & max_occ,nene1,nene2,nkpt,tetra,tweight,dtweightde, ierr)
807 : :
808 : : !Arguments ------------------------------------
809 : : !scalars
810 : : integer,intent(in) :: nene1,nene2,nkpt
811 : : integer,intent(out) :: ierr
812 : : type(t_tetrahedron), intent(in) :: tetra
813 : : real(dp),intent(in) :: enemax1,enemin1
814 : : real(dp),intent(in) :: enemax2,enemin2
815 : : real(dp),intent(in) :: max_occ
816 : : !arrays
817 : : real(dp),intent(in) :: eigen1_in(nkpt)
818 : : real(dp),intent(in) :: eigen2_in(nkpt)
819 : : real(dp),intent(out) :: dtweightde(nkpt,nene1,nene2),tweight(nkpt,nene1,nene2)
820 : :
821 : : !Local variables-------------------------------
822 : : ! needed for gaussian replacement of Dirac functions
823 : : ! the three coefficients of the DOS as quadratic form,
824 : : ! in the interval [eig(ikpt-1), eig(ikpt)]
825 : : ! for ikpt = 1 we add a point below eigen(1) which doesnt
826 : : ! contribute to the DOS in any tetrahedron
827 : : !scalars
828 : : integer :: ieps1,ieps2,itetra
829 : : integer :: nn1_1,nn1_2,nn1_3,nn1_4
830 : : integer :: nn2_1,nn2_2,nn2_3
831 : : integer :: ind_a(3), ind_b(3), ind_c(3)
832 : : real(dp) :: deltaene1,eps1
833 : : real(dp) :: deltaene2,eps2
834 : : ! real(dp) :: gau_prefactor,gau_width,gau_width2
835 : : real(dp) :: epsilon1(4,4)
836 : : real(dp) :: epsilon2(4,4)
837 : : real(dp) :: inv_epsilon1(4,4)
838 : : real(dp) :: aa(3),bb(3),cc(3)
839 : : real(dp) :: delaa(3),delbb(3),delcc(3)
840 : : real(dp) :: delaa0,delbb0,delcc0
841 : : real(dp) :: inv_delaa(3),inv_delbb(3),inv_delcc(3)
842 : : real(dp) :: deleps1, deleps2
843 : : real(dp) :: inv_deleps1
844 : : real(dp) :: dccde1, dccde1_pre
845 : : real(dp) :: volconst,volconst_mult
846 : : real(dp) :: ii0, ii1, ii3
847 : : !arrays
848 : : integer :: ind_k(4)
849 : 0 : real(dp), allocatable :: tweight_tmp(:,:,:)
850 : 0 : real(dp), allocatable :: dtweightde_tmp(:,:,:)
851 : : real(dp) :: eigen1_1tetra(4)
852 : : real(dp) :: eigen2_1tetra(4)
853 : :
854 : : ! *********************************************************************
855 : :
856 : 0 : ierr = 0
857 [ # # # # ]: 0 : if (nene1 <= 1 .or. nene2 <= 1) then
858 : : !'get_dbl_tetra_weight: nene must be at least 2'
859 : 0 : ierr = 1
860 : : return
861 : : else
862 : 0 : deltaene1 = (enemax1-enemin1) / (nene1-1)
863 : 0 : deltaene2 = (enemax2-enemin2) / (nene2-1)
864 : : end if
865 : :
866 [ # # # # : 0 : TETRA_ALLOCATE(tweight_tmp, (4, nene2, nene1))
# # # # #
# # # #
# ]
867 [ # # # # : 0 : TETRA_ALLOCATE(dtweightde_tmp, (4, nene2, nene1))
# # # # #
# # # ]
868 : :
869 : : !print *, "warning: for the moment, heaviside weights are 0. The delta function / DOS weights are the only ones calculated "
870 : :
871 : 0 : volconst = tetra%vv/4.d0
872 : :
873 : : ! for each tetrahedron
874 [ # # ]: 0 : do itetra=1,tetra%ntetra
875 : : ! these are for 1 tetrahedron only.
876 [ # # # # : 0 : tweight_tmp = zero
# # ]
877 [ # # # # : 0 : dtweightde_tmp = zero
# # ]
878 : :
879 : 0 : volconst_mult = max_occ*volconst*dble(tetra%tetra_mult(itetra))
880 : :
881 : : ! Here we need the original ordering to reference the correct irred kpoints
882 : : ! ind_k refers to the index in the full k list of the summits of the present tetrahedra
883 : : ! we can forget the order of the summits within the tetrahedron, because eigen1 fixes that
884 : : ! order with its increasing value
885 : 0 : ind_k(1) = tetra%tetra_full(1,1,itetra)
886 : 0 : ind_k(2) = tetra%tetra_full(2,1,itetra)
887 : 0 : ind_k(3) = tetra%tetra_full(3,1,itetra)
888 : 0 : ind_k(4) = tetra%tetra_full(4,1,itetra)
889 : 0 : eigen1_1tetra(1) = eigen1_in(ind_k(1))
890 : 0 : eigen1_1tetra(2) = eigen1_in(ind_k(2))
891 : 0 : eigen1_1tetra(3) = eigen1_in(ind_k(3))
892 : 0 : eigen1_1tetra(4) = eigen1_in(ind_k(4))
893 : 0 : call sort_tetra(4,eigen1_1tetra,ind_k,tol14)
894 : :
895 : : ! re-sort eigen2 values according to order chosen for eigen1. Eigen2 are _not_ in increasing order!
896 : 0 : eigen2_1tetra(1) = eigen2_in(ind_k(1))
897 : 0 : eigen2_1tetra(2) = eigen2_in(ind_k(2))
898 : 0 : eigen2_1tetra(3) = eigen2_in(ind_k(3))
899 : 0 : eigen2_1tetra(4) = eigen2_in(ind_k(4))
900 : :
901 : : ! the epsilons are energy differences for the two eigenvalue sets
902 : 0 : epsilon1 = zero
903 : 0 : epsilon2 = zero
904 [ # # ]: 0 : do ieps1 = 1, 4
905 [ # # ]: 0 : do ieps2 = ieps1+1, 4
906 : 0 : epsilon1(ieps1,ieps2) = eigen1_1tetra(ieps1)-eigen1_1tetra(ieps2)
907 : 0 : epsilon1(ieps2,ieps1) = -epsilon1(ieps1,ieps2)
908 : 0 : epsilon2(ieps1,ieps2) = eigen2_1tetra(ieps1)-eigen2_1tetra(ieps2)
909 : 0 : epsilon2(ieps2,ieps1) = -epsilon2(ieps1,ieps2)
910 : : end do
911 : : end do
912 : :
913 : : ! we precalculate the inverses to avoid doing tons of divisions in the energy loops below
914 : : ! Allen formulae only require the inverses of the differences of eigen1 + the a b c below
915 : 0 : inv_epsilon1 = zero
916 [ # # ]: 0 : do ieps1 = 1, 4
917 [ # # ]: 0 : do ieps2 = ieps1+1, 4
918 [ # # ]: 0 : if (abs(epsilon1(ieps1,ieps2)) > tol14) then
919 : 0 : inv_epsilon1(ieps1,ieps2) = 1.d0 / epsilon1(ieps1,ieps2)
920 : 0 : inv_epsilon1(ieps2,ieps1) = -inv_epsilon1(ieps1,ieps2)
921 : : end if
922 : : end do
923 : : end do
924 : :
925 : : ! these bounds determine the intervals for Omega in Allen paper, and cases A, B, C
926 : 0 : nn1_1 = int((eigen1_1tetra(1)-enemin1)/deltaene1)+1
927 : 0 : nn1_2 = int((eigen1_1tetra(2)-enemin1)/deltaene1)+1
928 : 0 : nn1_3 = int((eigen1_1tetra(3)-enemin1)/deltaene1)+1
929 : 0 : nn1_4 = int((eigen1_1tetra(4)-enemin1)/deltaene1)+1
930 : :
931 : 0 : nn1_1 = max(1,nn1_1)
932 : 0 : nn1_1 = min(nn1_1,nene1)
933 : 0 : nn1_2 = max(1,nn1_2)
934 : 0 : nn1_2 = min(nn1_2,nene1)
935 : 0 : nn1_3 = max(1,nn1_3)
936 : 0 : nn1_3 = min(nn1_3,nene1)
937 : 0 : nn1_4 = max(1,nn1_4)
938 : 0 : nn1_4 = min(nn1_4,nene1)
939 : :
940 : : ! calculate Allen a_i b_i and c_i parameters
941 : : ! sort the a_i b_i c_i
942 : : !
943 : : ! NOTE: indices here go from 1 to 4 instead of 0 to 3 as in Allen...
944 : 0 : aa(1) = epsilon2(2,1) * inv_epsilon1(2,1)
945 : 0 : aa(2) = epsilon2(3,1) * inv_epsilon1(3,1)
946 : 0 : aa(3) = epsilon2(4,1) * inv_epsilon1(4,1)
947 : 0 : ind_a = (/2,3,4/)
948 : 0 : call sort_tetra(3,aa,ind_a,tol14)
949 : : ! aa are now in order a_s a_m a_l !!! Preserve the hash function ind_a to order the positions of k below
950 : 0 : delaa(1) = aa(2)-aa(1)
951 : 0 : delaa(2) = aa(3)-aa(1)
952 : 0 : delaa(3) = aa(3)-aa(2)
953 : 0 : inv_delaa = zero
954 [ # # ]: 0 : if(delaa(1)> tol14) inv_delaa(1)= 1.0d0 / delaa(1)
955 [ # # ]: 0 : if(delaa(2)> tol14) inv_delaa(2)= 1.0d0 / delaa(2)
956 [ # # ]: 0 : if(delaa(3)> tol14) inv_delaa(3)= 1.0d0 / delaa(3)
957 : :
958 : 0 : bb(1) = epsilon2(1,2) * inv_epsilon1(1,2)
959 : 0 : bb(2) = epsilon2(3,2) * inv_epsilon1(3,2)
960 : 0 : bb(3) = epsilon2(4,2) * inv_epsilon1(4,2)
961 : 0 : ind_b = (/1,3,4/)
962 : 0 : call sort_tetra(3,bb,ind_b,tol14)
963 : 0 : delbb(1) = bb(2)-bb(1)
964 : 0 : delbb(2) = bb(3)-bb(1)
965 : 0 : delbb(3) = bb(3)-bb(2)
966 : 0 : inv_delbb = zero
967 [ # # ]: 0 : if(delbb(1)> tol14) inv_delbb(1)= 1.0d0 / delbb(1)
968 [ # # ]: 0 : if(delbb(2)> tol14) inv_delbb(2)= 1.0d0 / delbb(2)
969 [ # # ]: 0 : if(delbb(3)> tol14) inv_delbb(3)= 1.0d0 / delbb(3)
970 : :
971 : 0 : cc(1) = epsilon2(1,4) * inv_epsilon1(1,4)
972 : 0 : cc(2) = epsilon2(2,4) * inv_epsilon1(2,4)
973 : 0 : cc(3) = epsilon2(3,4) * inv_epsilon1(3,4)
974 : 0 : ind_c = (/1,2,3/)
975 : 0 : call sort_tetra(3,cc,ind_c,tol14)
976 : 0 : delcc(1) = cc(2)-cc(1)
977 : 0 : delcc(2) = cc(3)-cc(1)
978 : 0 : delcc(3) = cc(3)-cc(2)
979 : 0 : inv_delcc = zero
980 [ # # ]: 0 : if(delcc(1)> tol14) inv_delcc(1)= 1.0d0 / delcc(1)
981 [ # # ]: 0 : if(delcc(2)> tol14) inv_delcc(2)= 1.0d0 / delcc(2)
982 [ # # ]: 0 : if(delcc(3)> tol14) inv_delcc(3)= 1.0d0 / delcc(3)
983 : :
984 : : !----------------------------------------------------------------------
985 : : ! start main loop A B C over eps1
986 : : !----------------------------------------------------------------------
987 : :
988 : : !
989 : : ! interval enemin1 < eps1 < e1 nothing to do
990 : : !
991 : : !
992 : : ! interval e1 < eps1 < e3 CASE A in Allen + first term in B
993 : : !
994 : : ! NB: eps1 is not updated inside the loop, only between the loops
995 : 0 : eps1 = enemin1+nn1_1*deltaene1
996 : 0 : deleps1 = eps1-eigen1_1tetra(1) ! this is Omega - omega_0
997 : 0 : dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(2,1)*inv_epsilon1(3,1)*inv_epsilon1(4,1)
998 : :
999 : : ! note we go to nn1_3
1000 [ # # ]: 0 : do ieps1=nn1_1+1,nn1_3
1001 : :
1002 : 0 : dccde1 = dccde1_pre * deleps1 ! this is f_0(Omega)*6*v
1003 : :
1004 : : ! at fixed ieps1 we can find the pivot indices for the ieps2 loop
1005 : 0 : nn2_1 = int((eigen2_1tetra(1)+deleps1*aa(1) -enemin2)/deltaene2)+1
1006 : 0 : nn2_2 = int((eigen2_1tetra(1)+deleps1*aa(2) -enemin2)/deltaene2)+1
1007 : 0 : nn2_3 = int((eigen2_1tetra(1)+deleps1*aa(3) -enemin2)/deltaene2)+1
1008 : :
1009 : 0 : nn2_1 = max(1,nn2_1)
1010 : 0 : nn2_1 = min(nn2_1,nene2)
1011 : 0 : nn2_2 = max(1,nn2_2)
1012 : 0 : nn2_2 = min(nn2_2,nene2)
1013 : 0 : nn2_3 = max(1,nn2_3)
1014 : 0 : nn2_3 = min(nn2_3,nene2)
1015 : :
1016 : 0 : inv_deleps1 = 1.0d0 / deleps1
1017 : :
1018 : 0 : eps2 = enemin2+nn2_1*deltaene2 ! this is E
1019 : 0 : deleps2 = eps2 - eigen2_1tetra(1) ! this is E-epsilon_0
1020 : :
1021 : : !-----------------------------------------------------------------------
1022 : : ! This is case AI
1023 : : !-----------------------------------------------------------------------
1024 [ # # ]: 0 : do ieps2 = nn2_1+1, nn2_2
1025 : : ! calculate running value of del "a" = a-a_s: first term should really mix eps1 and eps2
1026 : 0 : delaa0 = deleps2*inv_deleps1 - aa(1) ! a - a_s
1027 : :
1028 : 0 : ii0 = dccde1*delaa0*inv_delaa(1)*inv_delaa(2) ! this is I_0(Omega E)
1029 : :
1030 : : dtweightde_tmp(1,ieps2,ieps1) = dtweightde_tmp(1,ieps2,ieps1) + &
1031 : : & ii0*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_a(1),1)* &
1032 : : & (-2.0d0 + delaa0*inv_delaa(1)*epsilon1(ind_a(2),ind_a(1))*inv_epsilon1(ind_a(2),1) &
1033 : 0 : & + delaa0*inv_delaa(2)*epsilon1(ind_a(3),ind_a(1))*inv_epsilon1(ind_a(3),1)))
1034 : : dtweightde_tmp(ind_a(1),ieps2,ieps1) = dtweightde_tmp(ind_a(1),ieps2,ieps1) + &
1035 : 0 : & ii0*0.5d0*deleps1*inv_epsilon1(ind_a(1),1)*(2.0d0 - delaa0*inv_delaa(1) - delaa0*inv_delaa(2))
1036 : : dtweightde_tmp(ind_a(2),ieps2,ieps1) = dtweightde_tmp(ind_a(2),ieps2,ieps1) + &
1037 : 0 : & ii0*0.5d0*delaa0*inv_delaa(1)*deleps1*inv_epsilon1(ind_a(2),1)
1038 : : dtweightde_tmp(ind_a(3),ieps2,ieps1) = dtweightde_tmp(ind_a(3),ieps2,ieps1) + &
1039 : 0 : & ii0*0.5d0*delaa0*inv_delaa(2)*deleps1*inv_epsilon1(ind_a(3),1)
1040 : 0 : deleps2 = deleps2 + deltaene2
1041 : : end do
1042 : :
1043 : :
1044 : 0 : eps2 = enemin2+nn2_2*deltaene2 ! this is E
1045 : 0 : deleps2 = eps2 - eigen2_1tetra(1) ! E-E_0
1046 : :
1047 : : !-----------------------------------------------------------------------
1048 : : ! This is case AII
1049 : : !-----------------------------------------------------------------------
1050 [ # # ]: 0 : do ieps2 = nn2_2+1, nn2_3
1051 : : ! calculate running value of del "a" = a_l-a: first term should really mix eps1 and eps2
1052 : 0 : delaa0 = aa(3) - deleps2*inv_deleps1 ! a_l - a
1053 : :
1054 : 0 : ii0 = dccde1*delaa0*inv_delaa(3)*inv_delaa(2) ! this is I_0(Omega E)
1055 : :
1056 : : dtweightde_tmp(1,ieps2,ieps1) = dtweightde_tmp(1,ieps2,ieps1) + &
1057 : : & ii0*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_a(3),1)* &
1058 : : & (-2.0d0 + delaa0*inv_delaa(3)*epsilon1(ind_a(2),ind_a(3))*inv_epsilon1(ind_a(2),1) &
1059 : 0 : & + delaa0*inv_delaa(2)*epsilon1(ind_a(1),ind_a(3))*inv_epsilon1(ind_a(1),1)))
1060 : : dtweightde_tmp(ind_a(3),ieps2,ieps1) = dtweightde_tmp(ind_a(3),ieps2,ieps1) + &
1061 : 0 : & ii0*0.5d0*deleps1*inv_epsilon1(ind_a(3),1)*(2.0d0 - delaa0*inv_delaa(3) - delaa0*inv_delaa(2))
1062 : : dtweightde_tmp(ind_a(2),ieps2,ieps1) = dtweightde_tmp(ind_a(2),ieps2,ieps1) + &
1063 : 0 : & ii0*0.5d0*delaa0*inv_delaa(3)*deleps1*inv_epsilon1(ind_a(2),1)
1064 : : dtweightde_tmp(ind_a(1),ieps2,ieps1) = dtweightde_tmp(ind_a(1),ieps2,ieps1) + &
1065 : 0 : & ii0*0.5d0*delaa0*inv_delaa(2)*deleps1*inv_epsilon1(ind_a(1),1)
1066 : :
1067 : 0 : deleps2 = deleps2 + deltaene2
1068 : : end do
1069 : 0 : deleps1 = deleps1 + deltaene1
1070 : : end do
1071 : : !
1072 : : ! interval e2 < eps < e3
1073 : : !
1074 : 0 : eps1 = eps1 + (nn1_2-nn1_1)*deltaene1
1075 : :
1076 : 0 : deleps1 = eps1-eigen1_1tetra(2) ! Omega - omega_1
1077 : :
1078 : 0 : dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(2,1)*inv_epsilon1(3,2)*inv_epsilon1(4,2) ! f1 function
1079 [ # # ]: 0 : do ieps1=nn1_2+1,nn1_3
1080 : :
1081 : 0 : dccde1 = dccde1_pre * deleps1 ! f2(Omega) * 6 * v
1082 : :
1083 : : ! at fixed ieps1 we can find the pivot indices for the ieps2 loop
1084 : 0 : nn2_1 = int((eigen2_1tetra(2)+deleps1*bb(1) -enemin2)/deltaene2)+1
1085 : 0 : nn2_2 = int((eigen2_1tetra(2)+deleps1*bb(2) -enemin2)/deltaene2)+1
1086 : 0 : nn2_3 = int((eigen2_1tetra(2)+deleps1*bb(3) -enemin2)/deltaene2)+1
1087 : :
1088 : 0 : nn2_1 = max(1,nn2_1)
1089 : 0 : nn2_1 = min(nn2_1,nene2)
1090 : 0 : nn2_2 = max(1,nn2_2)
1091 : 0 : nn2_2 = min(nn2_2,nene2)
1092 : 0 : nn2_3 = max(1,nn2_3)
1093 : 0 : nn2_3 = min(nn2_3,nene2)
1094 : :
1095 : 0 : inv_deleps1 = 1.0d0 / deleps1
1096 : :
1097 : 0 : eps2 = enemin2+nn2_1*deltaene2 ! starting value for E
1098 : 0 : deleps2 = eps2 - eigen2_1tetra(2) ! E - epsilon_1
1099 : :
1100 : : !-----------------------------------------------------------------------
1101 : : ! This is case BI
1102 : : !-----------------------------------------------------------------------
1103 [ # # ]: 0 : do ieps2 = nn2_1+1, nn2_2
1104 : : ! calculate running value of del "b" = b-b_s: first term should really mix eps1 and eps2
1105 : 0 : delbb0 = deleps2*inv_deleps1 - bb(1)
1106 : :
1107 : 0 : ii1 = dccde1*delbb0*inv_delbb(1)*inv_delbb(2) ! this is I_1(Omega E)
1108 : :
1109 : : ! note negative sign here - we are correcting the I0 a0 term already calculated above
1110 : : dtweightde_tmp(2,ieps2,ieps1) = dtweightde_tmp(2,ieps2,ieps1) - &
1111 : : & ii1*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_b(1),2)* &
1112 : : & (-2.0d0 + delbb0*inv_delbb(1)*epsilon1(ind_b(2),ind_b(1))*inv_epsilon1(ind_b(2),2) &
1113 : 0 : & + delbb0*inv_delbb(2)*epsilon1(ind_b(3),ind_b(1))*inv_epsilon1(ind_b(3),2)))
1114 : : dtweightde_tmp(ind_b(1),ieps2,ieps1) = dtweightde_tmp(ind_b(1),ieps2,ieps1) - &
1115 : 0 : & ii1*0.5d0*deleps1*inv_epsilon1(ind_b(1),2)*(2.0d0 - delbb0*inv_delbb(1) - delbb0*inv_delbb(2))
1116 : : dtweightde_tmp(ind_b(2),ieps2,ieps1) = dtweightde_tmp(ind_b(2),ieps2,ieps1) - &
1117 : 0 : & ii1*0.5d0*delbb0*inv_delbb(1)*deleps1*inv_epsilon1(ind_b(2),2)
1118 : : dtweightde_tmp(ind_b(3),ieps2,ieps1) = dtweightde_tmp(ind_b(3),ieps2,ieps1) - &
1119 : 0 : & ii1*0.5d0*delbb0*inv_delbb(2)*deleps1*inv_epsilon1(ind_b(3),2)
1120 : 0 : deleps2 = deleps2 + deltaene2
1121 : : end do
1122 : :
1123 : 0 : eps2 = enemin2+nn2_2*deltaene2
1124 : 0 : deleps2 = eps2 - eigen2_1tetra(2)
1125 : :
1126 : : !-----------------------------------------------------------------------
1127 : : ! This is case BII
1128 : : !-----------------------------------------------------------------------
1129 [ # # ]: 0 : do ieps2 = nn2_2+1, nn2_3
1130 : : ! calculate running value of del "b" = b_l-b: first term should really mix eps1 and eps2
1131 : 0 : delbb0 = bb(3) - deleps2*inv_deleps1
1132 : :
1133 : 0 : ii1 = dccde1*delbb0*inv_delbb(3)*inv_delbb(2) ! this is I_1(Omega E)
1134 : :
1135 : : ! note negative sign here - we are correcting the I0 a0 term already calculated above
1136 : : dtweightde_tmp(2,ieps2,ieps1) = dtweightde_tmp(2,ieps2,ieps1) - &
1137 : : & ii1*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_b(3),2)* &
1138 : : & (-2.0d0 + delbb0*inv_delbb(3)*epsilon1(ind_b(2),ind_b(3))*inv_epsilon1(ind_b(2),2) &
1139 : 0 : & + delbb0*inv_delbb(2)*epsilon1(ind_b(1),ind_b(3))*inv_epsilon1(ind_b(1),2)))
1140 : : dtweightde_tmp(ind_b(3),ieps2,ieps1) = dtweightde_tmp(ind_b(3),ieps2,ieps1) - &
1141 : 0 : & ii1*0.5d0*deleps1*inv_epsilon1(ind_b(3),2)*(2.0d0 - delbb0*inv_delbb(3) - delbb0*inv_delbb(2))
1142 : : dtweightde_tmp(ind_b(2),ieps2,ieps1) = dtweightde_tmp(ind_b(2),ieps2,ieps1) - &
1143 : 0 : & ii1*0.5d0*delbb0*inv_delbb(3)*deleps1*inv_epsilon1(ind_b(2),2)
1144 : : dtweightde_tmp(ind_b(1),ieps2,ieps1) = dtweightde_tmp(ind_b(1),ieps2,ieps1) - &
1145 : 0 : & ii1*0.5d0*delbb0*inv_delbb(2)*deleps1*inv_epsilon1(ind_b(1),2)
1146 : :
1147 : 0 : deleps2 = deleps2 + deltaene2
1148 : : end do
1149 : :
1150 : 0 : deleps1 = deleps1 + deltaene1
1151 : : end do
1152 : :
1153 : : !
1154 : : ! interval e3 < eps < e4
1155 : : !
1156 : 0 : eps1 = eps1 + (nn1_3-nn1_2)*deltaene1
1157 : 0 : deleps1 = eps1-eigen1_1tetra(4)
1158 : 0 : dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(4,1)*inv_epsilon1(4,2)*inv_epsilon1(4,3)
1159 [ # # ]: 0 : do ieps1=nn1_3+1,nn1_4
1160 : : ! note - sign from definition of f3
1161 : 0 : dccde1 = -dccde1_pre * deleps1 ! f3(Omega) * 6 * v
1162 : :
1163 : : ! at fixed ieps1 we can find the pivot indices for the ieps2 loop
1164 : : ! NB: order is inverted for cc because deleps1 is defined negative (Omega is always less than omega_3)
1165 : 0 : nn2_1 = int((eigen2_1tetra(4)+deleps1*cc(3) -enemin2)/deltaene2)+1
1166 : 0 : nn2_2 = int((eigen2_1tetra(4)+deleps1*cc(2) -enemin2)/deltaene2)+1
1167 : 0 : nn2_3 = int((eigen2_1tetra(4)+deleps1*cc(1) -enemin2)/deltaene2)+1
1168 : :
1169 : 0 : nn2_1 = max(1,nn2_1)
1170 : 0 : nn2_1 = min(nn2_1,nene2)
1171 : 0 : nn2_2 = max(1,nn2_2)
1172 : 0 : nn2_2 = min(nn2_2,nene2)
1173 : 0 : nn2_3 = max(1,nn2_3)
1174 : 0 : nn2_3 = min(nn2_3,nene2)
1175 : 0 : inv_deleps1 = 1.0d0 / deleps1
1176 : :
1177 : 0 : eps2 = enemin2+nn2_1*deltaene2 ! starting value for E
1178 : 0 : deleps2 = eps2 - eigen2_1tetra(4) ! E - epsilon_3
1179 : :
1180 : : !-----------------------------------------------------------------------
1181 : : ! This is case CII
1182 : : !-----------------------------------------------------------------------
1183 [ # # ]: 0 : do ieps2 = nn2_1+1, nn2_2
1184 : : ! calculate running value of del "c" = c_l-c: first term should really mix eps1 and eps2
1185 : 0 : delcc0 = cc(3) - deleps2*inv_deleps1
1186 : :
1187 : 0 : ii3 = dccde1*delcc0*inv_delcc(3)*inv_delcc(2) ! this is I_3(Omega E)
1188 : :
1189 : : dtweightde_tmp(4,ieps2,ieps1) = dtweightde_tmp(4,ieps2,ieps1) + &
1190 : : & ii3*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_c(3),4)* &
1191 : : & (-2.0d0 + delcc0*inv_delcc(3)*epsilon1(ind_c(2),ind_c(3))*inv_epsilon1(ind_c(2),4) &
1192 : 0 : & + delcc0*inv_delcc(2)*epsilon1(ind_c(1),ind_c(3))*inv_epsilon1(ind_c(1),4)))
1193 : : dtweightde_tmp(ind_c(3),ieps2,ieps1) = dtweightde_tmp(ind_c(3),ieps2,ieps1) + &
1194 : 0 : & ii3*0.5d0*deleps1*inv_epsilon1(ind_c(3),4)*(2.0d0 - delcc0*inv_delcc(3) - delcc0*inv_delcc(2))
1195 : : dtweightde_tmp(ind_c(2),ieps2,ieps1) = dtweightde_tmp(ind_c(2),ieps2,ieps1) + &
1196 : 0 : & ii3*0.5d0*delcc0*inv_delcc(3)*deleps1*inv_epsilon1(ind_c(2),4)
1197 : : dtweightde_tmp(ind_c(1),ieps2,ieps1) = dtweightde_tmp(ind_c(1),ieps2,ieps1) + &
1198 : 0 : & ii3*0.5d0*delcc0*inv_delcc(2)*deleps1*inv_epsilon1(ind_c(1),4)
1199 : :
1200 : 0 : deleps2 = deleps2 + deltaene2
1201 : : end do
1202 : :
1203 : :
1204 : 0 : eps2 = enemin2+nn2_2*deltaene2
1205 : 0 : deleps2 = eps2 - eigen2_1tetra(4)
1206 : :
1207 : : !-----------------------------------------------------------------------
1208 : : ! This is case CI
1209 : : !-----------------------------------------------------------------------
1210 [ # # ]: 0 : do ieps2 = nn2_2+1, nn2_3
1211 : : ! calculate running value of del "c" = c-c_s: first term should really mix eps1 and eps2
1212 : 0 : delcc0 = deleps2*inv_deleps1 - cc(1) ! c - c_s
1213 : :
1214 : 0 : ii3 = dccde1*delcc0*inv_delcc(1)*inv_delcc(2) ! this is I_3(Omega E)
1215 : :
1216 : : dtweightde_tmp(4,ieps2,ieps1) = dtweightde_tmp(4,ieps2,ieps1) + &
1217 : : & ii3*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_c(1),4)* &
1218 : : & (-2.0d0 + delcc0*inv_delcc(1)*epsilon1(ind_c(2),ind_c(1))*inv_epsilon1(ind_c(2),4) &
1219 : 0 : & + delcc0*inv_delcc(2)*epsilon1(ind_c(3),ind_c(1))*inv_epsilon1(ind_c(3),4)))
1220 : : dtweightde_tmp(ind_c(1),ieps2,ieps1) = dtweightde_tmp(ind_c(1),ieps2,ieps1) + &
1221 : 0 : & ii3*0.5d0*deleps1*inv_epsilon1(ind_c(1),4)*(2.0d0 - delcc0*inv_delcc(1) - delcc0*inv_delcc(2))
1222 : : dtweightde_tmp(ind_c(2),ieps2,ieps1) = dtweightde_tmp(ind_c(2),ieps2,ieps1) + &
1223 : 0 : & ii3*0.5d0*delcc0*inv_delcc(1)*deleps1*inv_epsilon1(ind_c(2),4)
1224 : : dtweightde_tmp(ind_c(3),ieps2,ieps1) = dtweightde_tmp(ind_c(3),ieps2,ieps1) + &
1225 : 0 : & ii3*0.5d0*delcc0*inv_delcc(2)*deleps1*inv_epsilon1(ind_c(3),4)
1226 : 0 : deleps2 = deleps2 + deltaene2
1227 : : end do
1228 : :
1229 : 0 : deleps1 = deleps1 + deltaene1
1230 : : end do
1231 : :
1232 : 0 : eps1 = eps1 + (nn1_4-nn1_3)*deltaene1
1233 : : !
1234 : : !
1235 : : ! interval e4 < eps < enemax
1236 : : !
1237 : 0 : do ieps1=nn1_4+1,nene1
1238 : : ! dtweightde unchanged by this tetrahedron
1239 : : end do
1240 : :
1241 : : ! if we have a fully degenerate tetrahedron,
1242 : : ! 1) the tweight is a Heaviside (step) function, which is correct above, but
1243 : : ! 2) the dtweightde should contain a Dirac function: add a Gaussian here
1244 : :
1245 : : ! TODO: add treatment in double tetra case
1246 : : ! end degenerate tetrahedron if
1247 : :
1248 : : ! the following blas calls are not working systematically, or do not give speed ups, strange...
1249 : : !call daxpy (nene, 1.d0, dtweightde_tmp(:,1), 1, dtweightde_t(:,ind_ibz(1)), 1)
1250 : : !call daxpy (nene, 1.d0, dtweightde_tmp(:,2), 1, dtweightde_t(:,ind_ibz(2)), 1)
1251 : : !call daxpy (nene, 1.d0, dtweightde_tmp(:,3), 1, dtweightde_t(:,ind_ibz(3)), 1)
1252 : : !call daxpy (nene, 1.d0, dtweightde_tmp(:,4), 1, dtweightde_t(:,ind_ibz(4)), 1)
1253 : :
1254 [ # # ]: 0 : do ieps2 = 1, nene2
1255 [ # # ]: 0 : dtweightde(ind_k(1),:,ieps2) = dtweightde(ind_k(1),:,ieps2) + dtweightde_tmp(1,ieps2,:)
1256 [ # # ]: 0 : dtweightde(ind_k(2),:,ieps2) = dtweightde(ind_k(2),:,ieps2) + dtweightde_tmp(2,ieps2,:)
1257 [ # # ]: 0 : dtweightde(ind_k(3),:,ieps2) = dtweightde(ind_k(3),:,ieps2) + dtweightde_tmp(3,ieps2,:)
1258 [ # # ]: 0 : dtweightde(ind_k(4),:,ieps2) = dtweightde(ind_k(4),:,ieps2) + dtweightde_tmp(4,ieps2,:)
1259 : : !tweight(nkpt,nene1,nene2)
1260 : : end do
1261 : :
1262 : : end do ! itetra
1263 : :
1264 : : ! transpose: otherwise the data access is crap and the code slows by an order of magnitude
1265 [ # # ]: 0 : TETRA_DEALLOCATE(tweight_tmp)
1266 [ # # ]: 0 : TETRA_DEALLOCATE(dtweightde_tmp)
1267 : :
1268 : : end subroutine get_dbl_tetra_weight
1269 : : !!***
1270 : :
1271 : : !!****f* m_tetrahedron/sort_tetra
1272 : : !! NAME
1273 : : !! sort_tetra
1274 : : !!
1275 : : !! FUNCTION
1276 : : !! Sort double precision array list(n) into ascending numerical order using Heapsort
1277 : : !! algorithm, while making corresponding rearrangement of the integer
1278 : : !! array iperm. Consider that two double precision numbers
1279 : : !! within tolerance tol are equal.
1280 : : !!
1281 : : !! INPUTS
1282 : : !! n intent(in) dimension of the list
1283 : : !! tol intent(in) numbers within tolerance are equal
1284 : : !! list(n) intent(inout) list of double precision numbers to be sorted
1285 : : !! iperm(n) intent(inout) iperm(i)=i (very important)
1286 : : !!
1287 : : !! OUTPUT
1288 : : !! list(n) sorted list
1289 : : !! iperm(n) index of permutation given the right ascending order
1290 : : !!
1291 : : !! SOURCE
1292 : :
1293 : :
1294 : 23332 : subroutine sort_tetra(n,list,iperm,tol)
1295 : :
1296 : : integer, intent(in) :: n
1297 : : integer, intent(inout) :: iperm(n)
1298 : : real(dp), intent(inout) :: list(n)
1299 : : real(dp), intent(in) :: tol
1300 : :
1301 : : integer :: l,ir,iap,i,j
1302 : : real(dp) :: ap
1303 : : character(len=500) :: msg
1304 : :
1305 [ - + ]: 23332 : if (n==1) then
1306 : : ! Accomodate case of array of length 1: already sorted!
1307 : 0 : return
1308 [ - + ]: 23332 : else if (n<1) then
1309 : : ! Should not call with n<1
1310 : 0 : write(msg,1000) n
1311 : : 1000 format(/,' sort_tetra has been called with array length n=',i12,/, &
1312 : : & ' having a value less than 1. This is not allowed.')
1313 : 0 : TETRA_ERROR(msg)
1314 : :
1315 : : else ! n>1
1316 : :
1317 : : ! Conduct the usual sort
1318 : 23332 : l=n/2+1
1319 : 23332 : ir=n
1320 : :
1321 : 101080 : do ! Infinite do-loop
1322 [ + + ]: 124412 : if (l>1) then
1323 : 49248 : l=l-1
1324 : 49248 : ap=list(l)
1325 : 49248 : iap=iperm(l)
1326 : :
1327 : : else ! l<=1
1328 : 75164 : ap=list(ir)
1329 : 75164 : iap=iperm(ir)
1330 : 75164 : list(ir)=list(1)
1331 : 75164 : iperm(ir)=iperm(1)
1332 : 75164 : ir=ir-1
1333 : :
1334 [ + + ]: 75164 : if (ir==1) then
1335 : 23332 : list(1)=ap
1336 : 23332 : iperm(1)=iap
1337 : 23332 : exit ! This is the end of this algorithm
1338 : : end if
1339 : : end if ! l>1
1340 : :
1341 : 101080 : i=l
1342 : 101080 : j=l+l
1343 : :
1344 [ + + ]: 252290 : do while (j<=ir)
1345 [ + + ]: 151210 : if (j<ir) then
1346 [ + + + + : 92732 : if ( list(j)<list(j+1)-tol .or. &
+ + ]
1347 : 151210 : & (list(j)<list(j+1)+tol.and.iperm(j)<iperm(j+1))) j=j+1
1348 : : endif
1349 [ + + + + : 252290 : if (ap<list(j)-tol .or. (ap<list(j)+tol.and.iap<iperm(j))) then
+ + ]
1350 : 119156 : list(i)=list(j)
1351 : 119156 : iperm(i)=iperm(j)
1352 : 119156 : i=j
1353 : 119156 : j=j+j
1354 : : else
1355 : 32054 : j=ir+1
1356 : : end if
1357 : : enddo
1358 : :
1359 : 101080 : list(i)=ap
1360 : 101080 : iperm(i)=iap
1361 : :
1362 : : enddo ! End infinite do-loop
1363 : :
1364 : : end if ! n>1
1365 : :
1366 : : end subroutine sort_tetra
1367 : : !!***
1368 : :
1369 : : !----------------------------------------------------------------------
1370 : :
1371 : : !!****f* m_tetrahedron/tetralib_has_mpi
1372 : : !! NAME
1373 : : !! tetralib_has_mpi
1374 : : !!
1375 : : !! FUNCTION
1376 : : !! Return True if library has been compiled with MPI support
1377 : : !!
1378 : : !! SOURCE
1379 : :
1380 : 0 : logical function tetralib_has_mpi() result(ans)
1381 : :
1382 : : ans = .False.
1383 : : #ifdef HAVE_MPI
1384 : 0 : ans = .True.
1385 : : #endif
1386 : :
1387 : 0 : end function tetralib_has_mpi
1388 : : !!***
1389 : :
1390 : : !----------------------------------------------------------------------
1391 : :
1392 : : !!****f* m_tetrahedron/split_work
1393 : : !! NAME
1394 : : !! split_work
1395 : : !!
1396 : : !! FUNCTION
1397 : : !! Splits the number of tasks, ntasks, among nprocs processors. Used for the MPI parallelization of simple loops.
1398 : : !!
1399 : : !! INPUTS
1400 : : !! ntasks=number of tasks
1401 : : !! comm=MPI communicator.
1402 : : !!
1403 : : !! OUTPUT
1404 : : !! nprocs=Number of MPI processes in the communicator.
1405 : : !! my_start,my_stop= indices defining the initial and final task for this processor
1406 : : !! ierr=Exit status.
1407 : : !!
1408 : : !! NOTES
1409 : : !! If nprocs>ntasks then :
1410 : : !! my_start=ntasks+1
1411 : : !! my_stop=ntask
1412 : : !!
1413 : : !! In this particular case, loops of the form
1414 : : !!
1415 : : !! do ii=my_start,my_stop
1416 : : !! ...
1417 : : !! end do
1418 : : !!
1419 : : !! are not executed. Moreover allocation such as foo(my_start:my_stop) will generate a zero-sized array.
1420 : : !!
1421 : : !! SOURCE
1422 : :
1423 : 18 : subroutine split_work(ntasks,comm,nprocs,my_start,my_stop,ierr)
1424 : :
1425 : : !Arguments ------------------------------------
1426 : : integer,intent(in) :: ntasks,comm
1427 : : integer,intent(out) :: nprocs,my_start,my_stop,ierr
1428 : :
1429 : : !Local variables-------------------------------
1430 : : integer :: res,my_rank,block_p1,block,mpierr
1431 : :
1432 : : ! *************************************************************************
1433 : :
1434 : 18 : nprocs = 1; my_start = 1; my_stop = ntasks; ierr = 1
1435 : : #ifdef HAVE_MPI
1436 [ + - ]: 18 : call MPI_COMM_SIZE(comm,nprocs,mpierr); if (mpierr /= MPI_SUCCESS) return
1437 [ + - ]: 18 : call MPI_COMM_RANK(comm,my_rank,mpierr); if (mpierr /= MPI_SUCCESS) return
1438 : :
1439 : 18 : block = ntasks/nprocs
1440 : 18 : res = MOD(ntasks,nprocs)
1441 : 18 : block_p1= block+1
1442 : :
1443 [ - + ]: 18 : if (my_rank<res) then
1444 : 0 : my_start = my_rank *block_p1+1
1445 : 0 : my_stop = (my_rank+1)*block_p1
1446 : : else
1447 : 18 : my_start = res*block_p1 + (my_rank-res )*block + 1
1448 : 18 : my_stop = res*block_p1 + (my_rank-res+1)*block
1449 : : end if
1450 : : #endif
1451 : 18 : ierr = 0
1452 : :
1453 : : end subroutine split_work
1454 : : !!***
1455 : :
1456 : : !----------------------------------------------------------------------
1457 : :
1458 : : !!****f* m_tetrahedron/get_onetetra_
1459 : : !! NAME
1460 : : !! get_onetetra_
1461 : : !!
1462 : : !! FUNCTION
1463 : : !! Private function to calculate the contributions to the weights due to a single tetrahedron.
1464 : : !! Extracted from get_tetra_weight
1465 : : !!
1466 : : !! SOURCE
1467 : :
1468 : 23328 : pure subroutine get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,nene,bcorr, &
1469 : 23328 : & tweight_tmp,dtweightde_tmp)
1470 : :
1471 : : !Arguments ------------------------------------
1472 : : !scalars
1473 : : integer,intent(in) :: nene,bcorr,itetra
1474 : : type(t_tetrahedron), intent(in) :: tetra
1475 : : real(dp) ,intent(in) :: enemax,enemin,max_occ
1476 : : !arrays
1477 : : ! MGTODO: This layout is not optimal (lots of cache thrashing, I will optimize it later on)
1478 : : real(dp), intent(out) :: tweight_tmp(nene, 4)
1479 : : real(dp), intent(out) :: dtweightde_tmp(nene, 4)
1480 : : real(dp),intent(in) :: eigen_1tetra(4)
1481 : :
1482 : : !Local variables-------------------------------
1483 : : ! needed for gaussian replacement of Dirac functions
1484 : : ! the three coefficients of the DOS as quadratic form,
1485 : : ! in the interval [eig(ikpt-1), eig(ikpt)]
1486 : : ! for ikpt = 1 we add a point below eigen(1) which doesnt
1487 : : ! contribute to the DOS in any tetrahedron
1488 : : !scalars
1489 : : integer :: ieps,nn1,nn2,nn3,nn4
1490 : : real(dp) :: cc,cc1,cc2,cc3,dcc1de,dcc2de,dcc3de,dccde,deltaene,eps
1491 : : real(dp) :: epsilon21,epsilon31,epsilon32,epsilon41,epsilon42,epsilon43
1492 : : real(dp) :: gau_prefactor,gau_width,gau_width2,inv_epsilon21,inv_epsilon31,gval
1493 : : real(dp) :: inv_epsilon32,inv_epsilon41,inv_epsilon42,inv_epsilon43
1494 : : real(dp) :: deleps1, deleps2, deleps3, deleps4
1495 : : real(dp) :: invepsum, cc_pre, dccde_pre
1496 : : real(dp) :: cc1_pre, cc2_pre, cc3_pre
1497 : : real(dp) :: cc_tmp, dccde_tmp
1498 : : real(dp) :: dcc1de_pre, dcc2de_pre, dcc3de_pre
1499 : : real(dp) :: tmp,volconst,volconst_mult
1500 : :
1501 : : ! *********************************************************************
1502 : :
1503 : 23328 : volconst = tetra%vv/4.d0
1504 : :
1505 : 23328 : deltaene = (enemax-enemin) / (nene-1)
1506 : :
1507 : : ! This is output
1508 [ + + + + : 37596960 : tweight_tmp = zero; dtweightde_tmp = zero
+ + + + ]
1509 : :
1510 : 23328 : volconst_mult = max_occ*volconst*dble(tetra%tetra_mult(itetra))
1511 : :
1512 : : ! all notations are from Blochl PRB 49 16223 [[cite:Bloechl1994a]] Appendix B
1513 : 23328 : epsilon21 = eigen_1tetra(2)-eigen_1tetra(1)
1514 : 23328 : epsilon31 = eigen_1tetra(3)-eigen_1tetra(1)
1515 : 23328 : epsilon41 = eigen_1tetra(4)-eigen_1tetra(1)
1516 : 23328 : epsilon32 = eigen_1tetra(3)-eigen_1tetra(2)
1517 : 23328 : epsilon42 = eigen_1tetra(4)-eigen_1tetra(2)
1518 : 23328 : epsilon43 = eigen_1tetra(4)-eigen_1tetra(3)
1519 [ + + ]: 23328 : inv_epsilon21 = zero; if (epsilon21 > tol14) inv_epsilon21 = 1.d0 / epsilon21
1520 [ + - ]: 23328 : inv_epsilon31 = zero; if (epsilon31 > tol14) inv_epsilon31 = 1.d0 / epsilon31
1521 [ + - ]: 23328 : inv_epsilon41 = zero; if (epsilon41 > tol14) inv_epsilon41 = 1.d0 / epsilon41
1522 [ + + ]: 23328 : inv_epsilon32 = zero; if (epsilon32 > tol14) inv_epsilon32 = 1.d0 / epsilon32
1523 [ + - ]: 23328 : inv_epsilon42 = zero; if (epsilon42 > tol14) inv_epsilon42 = 1.d0 / epsilon42
1524 [ + + ]: 23328 : inv_epsilon43 = zero; if (epsilon43 > tol14) inv_epsilon43 = 1.d0 / epsilon43
1525 : :
1526 : 23328 : nn1 = int((eigen_1tetra(1)-enemin)/deltaene)+1
1527 : 23328 : nn2 = int((eigen_1tetra(2)-enemin)/deltaene)+1
1528 : 23328 : nn3 = int((eigen_1tetra(3)-enemin)/deltaene)+1
1529 : 23328 : nn4 = int((eigen_1tetra(4)-enemin)/deltaene)+1
1530 : :
1531 : 23328 : nn1 = max(1,nn1)
1532 : 23328 : nn1 = min(nn1,nene)
1533 : 23328 : nn2 = max(1,nn2)
1534 : 23328 : nn2 = min(nn2,nene)
1535 : 23328 : nn3 = max(1,nn3)
1536 : 23328 : nn3 = min(nn3,nene)
1537 : 23328 : nn4 = max(1,nn4)
1538 : 23328 : nn4 = min(nn4,nene)
1539 : :
1540 : 23328 : eps = enemin+nn1*deltaene
1541 : : !
1542 : : !interval enemin < eps < e1 nothing to do
1543 : : !
1544 : : !
1545 : : !interval e1 < eps < e2
1546 : : !
1547 : 23328 : deleps1 = eps-eigen_1tetra(1)
1548 : 23328 : cc_pre = volconst_mult*inv_epsilon21*inv_epsilon31*inv_epsilon41
1549 : 23328 : invepsum = inv_epsilon21+inv_epsilon31+inv_epsilon41
1550 : 23328 : dccde_pre = 3.d0*volconst_mult*inv_epsilon21*inv_epsilon31*inv_epsilon41
1551 [ + + ]: 302416 : do ieps=nn1+1,nn2
1552 : 279088 : cc = cc_pre * deleps1*deleps1*deleps1
1553 : 279088 : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + cc*(4.d0-deleps1*invepsum)
1554 : 279088 : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + cc*deleps1*inv_epsilon21
1555 : 279088 : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + cc*deleps1*inv_epsilon31
1556 : 279088 : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + cc*deleps1*inv_epsilon41
1557 : :
1558 : 279088 : dccde = dccde_pre * deleps1*deleps1
1559 : 279088 : dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + dccde*(4.d0 - deleps1*invepsum) -cc*invepsum
1560 : 279088 : dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + (dccde*deleps1 + cc) * inv_epsilon21
1561 : 279088 : dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + (dccde*deleps1 + cc) * inv_epsilon31
1562 : 279088 : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + (dccde*deleps1 + cc) * inv_epsilon41
1563 : :
1564 [ - + ]: 279088 : if (bcorr == 1) then
1565 : : ! bxu, correction terms based on Bloechl's paper
1566 : : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1567 : 0 : & 4.d0*dccde_pre*deleps1*deleps1*(epsilon21+epsilon31+epsilon41)/40.d0
1568 : : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1569 : 0 : & 4.d0*dccde_pre*deleps1*deleps1*(-epsilon21+epsilon32+epsilon42)/40.d0
1570 : : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1571 : 0 : & 4.d0*dccde_pre*deleps1*deleps1*(-epsilon31-epsilon32+epsilon43)/40.d0
1572 : : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1573 : 0 : & 4.d0*dccde_pre*deleps1*deleps1*(-epsilon41-epsilon42-epsilon43)/40.d0
1574 : :
1575 : : dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + &
1576 : 0 : & 8.d0*dccde_pre*deleps1*(epsilon21+epsilon31+epsilon41)/40.d0
1577 : : dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + &
1578 : 0 : & 8.d0*dccde_pre*deleps1*(-epsilon21+epsilon32+epsilon42)/40.d0
1579 : : dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + &
1580 : 0 : & 8.d0*dccde_pre*deleps1*(-epsilon31-epsilon32+epsilon43)/40.d0
1581 : : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + &
1582 : 0 : & 8.d0*dccde_pre*deleps1*(-epsilon41-epsilon42-epsilon43)/40.d0
1583 : : end if
1584 : :
1585 : 302416 : deleps1 = deleps1 + deltaene
1586 : : end do
1587 : :
1588 : 23328 : eps = eps + (nn2-nn1)*deltaene
1589 : : !
1590 : : ! interval e2 < eps < e3
1591 : : !
1592 : 23328 : deleps1 = eps-eigen_1tetra(1)
1593 : 23328 : deleps2 = eps-eigen_1tetra(2)
1594 : 23328 : deleps3 = eigen_1tetra(3)-eps
1595 : 23328 : deleps4 = eigen_1tetra(4)-eps
1596 : :
1597 : 23328 : cc1_pre = volconst_mult*inv_epsilon31*inv_epsilon41
1598 : 23328 : cc2_pre = volconst_mult*inv_epsilon41*inv_epsilon32*inv_epsilon31
1599 : 23328 : cc3_pre = volconst_mult*inv_epsilon42*inv_epsilon32*inv_epsilon41
1600 : :
1601 : 23328 : dcc1de_pre = 2.d0*cc1_pre
1602 : 23328 : dcc2de_pre = cc2_pre
1603 : 23328 : dcc3de_pre = cc3_pre
1604 [ + + ]: 228264 : do ieps=nn2+1,nn3
1605 : 204936 : cc1 = cc1_pre * deleps1*deleps1
1606 : 204936 : cc2 = cc2_pre * deleps1*deleps2*deleps3
1607 : 204936 : cc3 = cc3_pre * deleps2*deleps2*deleps4
1608 : :
1609 : : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1610 : 204936 : & cc1 + (cc1+cc2)*deleps3*inv_epsilon31 + (cc1+cc2+cc3)*deleps4*inv_epsilon41
1611 : : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1612 : 204936 : & cc1+cc2+cc3+(cc2+cc3)*deleps3*inv_epsilon32 + cc3*deleps4*inv_epsilon42
1613 : : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1614 : 204936 : & (cc1+cc2)*deleps1*inv_epsilon31 + (cc2+cc3)*deleps2*inv_epsilon32
1615 : : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1616 : 204936 : & (cc1+cc2+cc3)*deleps1*inv_epsilon41 + cc3*deleps2*inv_epsilon42
1617 : :
1618 : :
1619 : 204936 : dcc1de = dcc1de_pre * deleps1
1620 : 204936 : dcc2de = dcc2de_pre * (-deleps1*deleps2 +deleps1*deleps3 +deleps2*deleps3)
1621 : 204936 : dcc3de = dcc3de_pre * (2.d0*deleps2*deleps4 -deleps2*deleps2)
1622 : :
1623 : : dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) &
1624 : : & + dcc1de &
1625 : : & + ((dcc1de+dcc2de)*deleps3 -(cc1+cc2)) * inv_epsilon31 &
1626 : 204936 : & + ((dcc1de+dcc2de+dcc3de)*deleps4 -(cc1+cc2+cc3)) * inv_epsilon41
1627 : : dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) &
1628 : : & + dcc1de+dcc2de+dcc3de &
1629 : : & + ((dcc2de+dcc3de)*deleps3 -(cc2+cc3) ) * inv_epsilon32 &
1630 : 204936 : & + (dcc3de*deleps4 -cc3 ) * inv_epsilon42
1631 : : dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) &
1632 : : & + ((dcc1de+dcc2de)*deleps1 + (cc1+cc2) ) * inv_epsilon31 &
1633 : 204936 : & + ((dcc2de+dcc3de)*deleps2 + (cc2+cc3) ) * inv_epsilon32
1634 : : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) &
1635 : : & + ((dcc1de+dcc2de+dcc3de)*deleps1 + (cc1+cc2+cc3) ) * inv_epsilon41 &
1636 : 204936 : & + (dcc3de*deleps2 + cc3) * inv_epsilon42
1637 : :
1638 [ - + ]: 204936 : if (bcorr == 1) then
1639 : : ! bxu, correction terms based on Bloechl's paper
1640 : : ! The correction terms may cause the dtweightde become negative
1641 : : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1642 : : & 4.d0*cc1_pre* &
1643 : : & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1644 : 0 : & (epsilon21+epsilon31+epsilon41)/40.d0
1645 : : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1646 : : & 4.d0*cc1_pre* &
1647 : : & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1648 : 0 : & (-epsilon21+epsilon32+epsilon42)/40.d0
1649 : : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1650 : : & 4.d0*cc1_pre* &
1651 : : & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1652 : 0 : & (-epsilon31-epsilon32+epsilon43)/40.d0
1653 : : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1654 : : & 4.d0*cc1_pre* &
1655 : : & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1656 : 0 : & (-epsilon41-epsilon42-epsilon43)/40.d0
1657 : :
1658 : : dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + &
1659 : : & 4.d0*cc1_pre* &
1660 : : & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1661 : 0 : & (epsilon21+epsilon31+epsilon41)/40.d0
1662 : : dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + &
1663 : : & 4.d0*cc1_pre* &
1664 : : & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1665 : 0 : & (-epsilon21+epsilon32+epsilon42)/40.d0
1666 : : dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + &
1667 : : & 4.d0*cc1_pre* &
1668 : : & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1669 : 0 : & (-epsilon31-epsilon32+epsilon43)/40.d0
1670 : : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + &
1671 : : & 4.d0*cc1_pre* &
1672 : : & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1673 : 0 : & (-epsilon41-epsilon42-epsilon43)/40.d0
1674 : : end if
1675 : :
1676 : 204936 : deleps1 = deleps1 + deltaene
1677 : 204936 : deleps2 = deleps2 + deltaene
1678 : 204936 : deleps3 = deleps3 - deltaene
1679 : 228264 : deleps4 = deleps4 - deltaene
1680 : : end do
1681 : :
1682 : 23328 : eps = eps + (nn3-nn2)*deltaene
1683 : : !
1684 : : ! interval e3 < eps < e4
1685 : : !
1686 : 23328 : deleps4 = eigen_1tetra(4)-eps
1687 : 23328 : cc_pre = volconst_mult*inv_epsilon41*inv_epsilon42*inv_epsilon43
1688 : 23328 : invepsum = inv_epsilon41+inv_epsilon42+inv_epsilon43
1689 : 23328 : dccde_pre = -3.d0*cc_pre
1690 [ + + ]: 244400 : do ieps=nn3+1,nn4
1691 : 221072 : cc = cc_pre * deleps4*deleps4*deleps4
1692 : 221072 : cc_tmp = cc * deleps4
1693 : 221072 : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + volconst_mult - cc_tmp*inv_epsilon41
1694 : 221072 : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + volconst_mult - cc_tmp*inv_epsilon42
1695 : 221072 : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + volconst_mult - cc_tmp*inv_epsilon43
1696 : 221072 : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + volconst_mult - cc*4.d0 + cc_tmp*invepsum
1697 : :
1698 : 221072 : dccde = dccde_pre * deleps4*deleps4
1699 : 221072 : dccde_tmp = -dccde*deleps4 + cc
1700 : 221072 : dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + dccde_tmp * inv_epsilon41
1701 : 221072 : dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + dccde_tmp * inv_epsilon42
1702 : 221072 : dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + dccde_tmp * inv_epsilon43
1703 : 221072 : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) - dccde*4.d0 - dccde_tmp*invepsum
1704 : :
1705 [ - + ]: 221072 : if (bcorr == 1) then
1706 : : ! bxu, correction terms based on Bloechl's paper
1707 : : ! The correction terms may cause the dtweightde become negative
1708 : : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1709 : 0 : & 12.d0*cc_pre*deleps4*deleps4*(epsilon21+epsilon31+epsilon41)/40.d0
1710 : : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1711 : 0 : & 12.d0*cc_pre*deleps4*deleps4*(-epsilon21+epsilon32+epsilon42)/40.d0
1712 : : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1713 : 0 : & 12.d0*cc_pre*deleps4*deleps4*(-epsilon31-epsilon32+epsilon43)/40.d0
1714 : : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1715 : 0 : & 12.d0*cc_pre*deleps4*deleps4*(-epsilon41-epsilon42-epsilon43)/40.d0
1716 : :
1717 : : dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) - &
1718 : 0 : & 24.d0*cc_pre*deleps4*(epsilon21+epsilon31+epsilon41)/40.d0
1719 : : dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) - &
1720 : 0 : & 24.d0*cc_pre*deleps4*(-epsilon21+epsilon32+epsilon42)/40.d0
1721 : : dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) - &
1722 : 0 : & 24.d0*cc_pre*deleps4*(-epsilon31-epsilon32+epsilon43)/40.d0
1723 : : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) - &
1724 : 0 : & 24.d0*cc_pre*deleps4*(-epsilon41-epsilon42-epsilon43)/40.d0
1725 : : end if
1726 : :
1727 : 244400 : deleps4 = deleps4 - deltaene
1728 : : end do
1729 : 23328 : eps = eps + (nn4-nn3)*deltaene
1730 : : !
1731 : : !
1732 : : ! interval e4 < eps < enemax
1733 : : !
1734 [ + + ]: 1429904 : do ieps=nn4+1,nene
1735 : 1406576 : tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + volconst_mult
1736 : 1406576 : tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + volconst_mult
1737 : 1406576 : tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + volconst_mult
1738 : 1429904 : tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + volconst_mult
1739 : : ! dtweightde unchanged by this tetrahedron
1740 : : end do
1741 : :
1742 : : !
1743 : : ! if we have a fully degenerate tetrahedron,
1744 : : ! 1) the tweight is a Heaviside (step) function, which is correct above, but
1745 : : ! 2) the dtweightde should contain a Dirac function: add a Gaussian here
1746 : : !
1747 [ - + ]: 23328 : if (epsilon41 < tol14) then
1748 : :
1749 : : ! to ensure the gaussian will integrate properly:
1750 : : ! WARNING: this smearing could be problematic if too large
1751 : : ! and doesnt integrate well if its too small
1752 : 0 : gau_width = 10.0d0*deltaene
1753 : 0 : gau_width2 = 1.0 / gau_width / gau_width
1754 : 0 : gau_prefactor = volconst_mult / gau_width / sqrt(pi)
1755 : : !
1756 : : ! average position since bracket for epsilon41 is relatively large
1757 : 0 : cc = (eigen_1tetra(1)+eigen_1tetra(2)+eigen_1tetra(3)+eigen_1tetra(4))/4.d0
1758 : 0 : eps = enemin
1759 [ # # ]: 0 : do ieps=1,nene
1760 : 0 : tmp = eps - cc
1761 : 0 : gval = gau_prefactor*exp(-tmp*tmp*gau_width2)
1762 : : ! MG TODO: I think this is not correct, because we have divided by 4 so
1763 : : ! the other points should be accumulated as well.
1764 : : ! There are however changes in the unit tests if I activate these lines...
1765 : : !dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + gval
1766 : : !dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + gval
1767 : : !dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + gval
1768 : 0 : dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + gval
1769 : 0 : eps = eps + deltaene
1770 : : end do
1771 : : end if ! end degenerate tetrahedron if
1772 : :
1773 : 23328 : end subroutine get_onetetra_
1774 : : !!***
1775 : :
1776 : : !----------------------------------------------------------------------
1777 : :
1778 : : !!****f* m_tetrahedron/tetra_get_onewk
1779 : : !! NAME
1780 : : !! tetra_get_onewk
1781 : : !!
1782 : : !! FUNCTION
1783 : : !! Calculate integration weights and their derivatives for a single k-point in the IBZ.
1784 : : !!
1785 : : !! INPUTS
1786 : : !! tetra<t_tetrahedron>=Object with tables for tetrahedron method.
1787 : : !! ik_ibz=Index of the k-point in the IBZ array
1788 : : !! bcorr=1 to include Blochl correction else 0.
1789 : : !! nene=number of energies for DOS
1790 : : !! nibz=number of irreducible kpoints
1791 : : !! eigen_ibz(nkibz)=eigenenergies for each k point
1792 : : !! enemin=minimal energy for DOS
1793 : : !! enemax=maximal energy for DOS
1794 : : !! max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)
1795 : : !!
1796 : : !! OUTPUT
1797 : : !! weights(nene,2) = integration weights for
1798 : : !! Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
1799 : : !! for a given (band, k-point, spin).
1800 : : !!
1801 : : !! SOURCE
1802 : :
1803 : 0 : subroutine tetra_get_onewk(tetra,ik_ibz,bcorr,nene,nkibz,eig_ibz,enemin,enemax,max_occ,weights)
1804 : :
1805 : : !Arguments ------------------------------------
1806 : : !scalars
1807 : : integer,intent(in) :: ik_ibz,nene,nkibz,bcorr
1808 : : type(t_tetrahedron), intent(in) :: tetra
1809 : : real(dp) ,intent(in) :: enemin,enemax,max_occ
1810 : : !arrays
1811 : : real(dp),intent(in) :: eig_ibz(nkibz)
1812 : : real(dp),intent(out) :: weights(nene,2)
1813 : :
1814 : : !Local variables-------------------------------
1815 : : !scalars
1816 : : integer :: itetra,ii
1817 : : !arrays
1818 : : integer :: ind_ibz(4)
1819 : 0 : real(dp) :: tweight_tmp(nene,4),dtweightde_tmp(nene,4),eigen_1tetra(4)
1820 : :
1821 : : ! *********************************************************************
1822 : :
1823 [ # # # # ]: 0 : weights = zero
1824 : :
1825 : : ! For each tetrahedron
1826 [ # # ]: 0 : do itetra=1,tetra%ntetra
1827 : :
1828 : : ! Here we need the original ordering to reference the correct irred kpoints
1829 [ # # ]: 0 : ind_ibz(:) = tetra%tetra_full(:,1,itetra)
1830 : : ! Cycle if this tetra does not contribute to this k-point.
1831 [ # # # # : 0 : if (all(ind_ibz /= ik_ibz)) cycle
# # ]
1832 : :
1833 : : ! Sort energies before calling get_onetetra_
1834 [ # # ]: 0 : eigen_1tetra(:) = eig_ibz(ind_ibz(:))
1835 : 0 : call sort_tetra(4, eigen_1tetra, ind_ibz, tol14)
1836 : :
1837 : : call get_onetetra_(tetra, itetra, eigen_1tetra, enemin, enemax, max_occ, nene, bcorr, &
1838 : 0 : tweight_tmp, dtweightde_tmp)
1839 : :
1840 : : ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz)
1841 [ # # ]: 0 : do ii=1,4
1842 [ # # ]: 0 : if (ind_ibz(ii) == ik_ibz) then
1843 [ # # ]: 0 : weights(:,1) = weights(:,1) + dtweightde_tmp(:,ii)
1844 [ # # ]: 0 : weights(:,2) = weights(:,2) + tweight_tmp(:,ii)
1845 : : end if
1846 : : end do
1847 : : end do ! itetra
1848 : :
1849 : 0 : end subroutine tetra_get_onewk
1850 : : !!***
1851 : :
1852 : : !----------------------------------------------------------------------
1853 : :
1854 : : !!****f* m_tetrahedron/tetra_get_onewk_wvals
1855 : : !! NAME
1856 : : !! tetra_get_onewk_wvals
1857 : : !!
1858 : : !! FUNCTION
1859 : : !! Calculate integration weights and their derivatives for a single k-point in the IBZ.
1860 : : !!
1861 : : !! INPUTS
1862 : : !! tetra<t_tetrahedron>=Object with tables for tetrahedron method.
1863 : : !! ik_ibz=Index of the k-point in the IBZ array
1864 : : !! bcorr=1 to include Blochl correction else 0.
1865 : : !! nw=number of energies in wvals
1866 : : !! nibz=number of irreducible kpoints
1867 : : !! wvals(nw)=Frequency points.
1868 : : !! eigen_ibz(nkibz)=eigenenergies for each k point
1869 : : !! [wtol]: If present, frequency points that differ by less that wtol are treated as equivalent.
1870 : : !! and the tetrahedron integration is performed only once per frequency point.
1871 : : !!
1872 : : !! OUTPUT
1873 : : !! weights(nw,2) = integration weights for
1874 : : !! Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
1875 : : !! for a given (band, k-point, spin).
1876 : : !!
1877 : : !! SOURCE
1878 : :
1879 : 0 : subroutine tetra_get_onewk_wvals(tetra, ik_ibz, bcorr, nw, wvals, nkibz, eig_ibz, weights, wtol)
1880 : :
1881 : : !Arguments ------------------------------------
1882 : : !scalars
1883 : : integer,intent(in) :: ik_ibz,nw,nkibz,bcorr
1884 : : real(dp), optional, intent(in) :: wtol
1885 : : type(t_tetrahedron), intent(in) :: tetra
1886 : : !arrays
1887 : : real(dp),intent(in) :: wvals(nw)
1888 : : real(dp),intent(in) :: eig_ibz(nkibz)
1889 : : real(dp),intent(out) :: weights(nw, 2)
1890 : :
1891 : : !Local variables-------------------------------
1892 : : !scalars
1893 : : !integer,save :: done = 0
1894 : : integer,parameter :: nene=3
1895 : : integer :: itetra,ii,jj,iw,ie
1896 : : logical :: samew
1897 : : real(dp),parameter :: max_occ1 = one
1898 : : real(dp) :: enemin, enemax
1899 : : !arrays
1900 : : integer :: ind_ibz(4)
1901 : : real(dp) :: theta_tmp(nene,4), delta_tmp(nene,4), eigen_1tetra(4)
1902 : :
1903 : : ! *********************************************************************
1904 : :
1905 [ # # # # ]: 0 : weights = zero
1906 : :
1907 : : ! For each tetrahedron
1908 [ # # ]: 0 : do jj=1,tetra%ibz_tetra_count(ik_ibz)
1909 : 0 : itetra = tetra%ibz_tetra_mapping(ik_ibz,jj)
1910 : :
1911 : : ! Here we need the original ordering to reference the correct irred kpoints
1912 [ # # ]: 0 : ind_ibz(:) = tetra%tetra_full(:,1,itetra)
1913 : :
1914 : : ! Sort energies before calling get_onetetra_
1915 [ # # ]: 0 : eigen_1tetra(:) = eig_ibz(ind_ibz(:))
1916 : 0 : call sort_tetra(4, eigen_1tetra, ind_ibz, tol14)
1917 : :
1918 [ # # ]: 0 : do iw=1,nw
1919 : 0 : samew = .False.
1920 [ # # ]: 0 : if (present(wtol)) then
1921 [ # # ]: 0 : if (iw > 1) samew = abs(wvals(iw) - wvals(iw - 1)) < wtol
1922 : : end if
1923 [ # # ]: 0 : if (.not. samew) then
1924 : 0 : enemin = wvals(iw) - 0.01; enemax = wvals(iw) + 0.01
1925 : 0 : ie = nene / 2 + 1
1926 : : call get_onetetra_(tetra, itetra, eigen_1tetra, enemin, enemax, max_occ1, nene, bcorr, &
1927 : 0 : theta_tmp, delta_tmp)
1928 : : end if
1929 : :
1930 : : ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz)
1931 [ # # ]: 0 : do ii=1,4
1932 [ # # ]: 0 : if (ind_ibz(ii) == ik_ibz) then
1933 : 0 : weights(iw, 1) = weights(iw, 1) + delta_tmp(ie, ii)
1934 : 0 : weights(iw, 2) = weights(iw, 2) + theta_tmp(ie, ii)
1935 : : end if
1936 : : end do
1937 : : end do ! iw
1938 : : end do ! itetra
1939 : :
1940 : 0 : end subroutine tetra_get_onewk_wvals
1941 : : !!***
1942 : :
1943 : : !----------------------------------------------------------------------
1944 : :
1945 : : !!****f* m_tetrahedron/tetra_get_onetetra_wvals
1946 : : !! NAME
1947 : : !! tetra_get_onetetra_wvals
1948 : : !!
1949 : : !! FUNCTION
1950 : : !! Calculate integration weights and their derivatives for a single k-point in the IBZ.
1951 : : !!
1952 : : !! INPUTS
1953 : : !! tetra<t_tetrahedron>=Object with tables for tetrahedron method.
1954 : : !! ik_ibz=Index of the k-point in the IBZ array
1955 : : !! bcorr=1 to include Blochl correction else 0.
1956 : : !! nw=number of energies in wvals
1957 : : !! nibz=number of irreducible kpoints
1958 : : !! wvals(nw)=Frequency points.
1959 : : !! eigen_ibz(nkibz)=eigenenergies for each k point
1960 : : !! [wtol]: If present, frequency points that differ by less that wtol are treated as equivalent.
1961 : : !! and the tetrahedron integration is performed only once per frequency point.
1962 : : !!
1963 : : !! OUTPUT
1964 : : !! weights(nw,2) = integration weights for
1965 : : !! Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
1966 : : !! for a given (band, k-point, spin).
1967 : : !!
1968 : : !! SOURCE
1969 : :
1970 : 0 : subroutine tetra_get_onetetra_wvals(tetra, itetra, eigen_1tetra, bcorr, nw, wvals, weights, wtol)
1971 : :
1972 : : !Arguments ------------------------------------
1973 : : !scalars
1974 : : integer,intent(in) :: nw,bcorr
1975 : : real(dp), optional, intent(in) :: wtol
1976 : : type(t_tetrahedron), intent(in) :: tetra
1977 : : !arrays
1978 : : real(dp),intent(in) :: wvals(nw)
1979 : : real(dp),intent(out) :: weights(nw, 2, 4)
1980 : :
1981 : : !Local variables-------------------------------
1982 : : !scalars
1983 : : !integer,save :: done = 0
1984 : : integer,parameter :: nene3=3
1985 : : integer :: itetra,ii,idx,iw,ie
1986 : : integer :: ind(4)
1987 : : logical :: samew
1988 : : real(dp),parameter :: max_occ1 = one
1989 : : real(dp) :: enemin, enemax
1990 : : !arrays
1991 : : real(dp) :: theta_tmp(nene3,4), delta_tmp(nene3,4), eigen_1tetra(4)
1992 : :
1993 : : ! *********************************************************************
1994 : :
1995 : 0 : ind = [1,2,3,4]
1996 : 0 : call sort_tetra(4, eigen_1tetra, ind, tol14)
1997 [ # # # # : 0 : weights = 0
# # ]
1998 : :
1999 : : !for all the frequencies
2000 [ # # ]: 0 : do iw=1,nw
2001 : 0 : samew = .False.
2002 [ # # ]: 0 : if (present(wtol)) then
2003 [ # # ]: 0 : if (iw > 1) samew = abs(wvals(iw) - wvals(iw - 1)) < wtol
2004 : : end if
2005 [ # # ]: 0 : if (.not. samew) then
2006 : 0 : enemin = wvals(iw) - 0.01
2007 : 0 : enemax = wvals(iw) + 0.01
2008 : 0 : ie = nene3 / 2 + 1
2009 : : call get_onetetra_(tetra, itetra, eigen_1tetra, enemin, enemax, max_occ1, nene3, bcorr, &
2010 : 0 : theta_tmp, delta_tmp)
2011 : : end if
2012 : :
2013 : : ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz)
2014 [ # # ]: 0 : do ii=1,4
2015 : 0 : idx = ind(ii)
2016 : 0 : weights(iw, 1, idx) = weights(iw, 1, idx) + delta_tmp(ie, ii)
2017 : 0 : weights(iw, 2, idx) = weights(iw, 2, idx) + theta_tmp(ie, ii)
2018 : : end do
2019 : : end do !iw
2020 : :
2021 : 0 : end subroutine tetra_get_onetetra_wvals
2022 : : !!***
2023 : :
2024 [ # # # # : 0 : end module m_tetrahedron
# # # # #
# # # # #
# # # # #
# # # ]
2025 : : !!***
|
