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