Line data Source code
1 : !!****m* ABINIT/m_fstab
2 : !! NAME
3 : !! m_fstab
4 : !!
5 : !! FUNCTION
6 : !! Tools for the management of a set of Fermi surface k-points.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG, MVer)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_fstab
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_xmpi
27 : use m_errors
28 : use m_krank
29 : use m_htetra
30 : use m_ebands
31 : use m_crystal
32 : use m_dtset
33 :
34 : use m_time, only : cwtime, cwtime_report
35 : use m_fstrings, only : itoa, sjoin, ktoa
36 : use m_numeric_tools, only : bisect
37 : use m_matrix, only : matr3inv
38 : use m_special_funcs, only : gaussian
39 : use m_kpts, only : kpts_timrev_from_kptopt, smpbz, kpts_map
40 :
41 : implicit none
42 :
43 : private
44 : !!***
45 :
46 : !!****t* m_fstab/fstab_t
47 : !! NAME
48 : !! fstab_t
49 : !!
50 : !! FUNCTION
51 : !! Tables with the correspondence between k-points on the Fermi surface (FS) and k-points
52 : !! in the IBZ (i.e. the k-points found in ebands_t).
53 : !! We use `nsppol` fstab_t objects to account for spin polarization and possibly different
54 : !! number of bands crossing the Fermi level.
55 : !!
56 : !! SOURCE
57 :
58 : type,public :: fstab_t
59 :
60 : integer :: spin = -1
61 : ! Spin index
62 :
63 : integer :: nkfs = -1
64 : ! Number of k-points on the Fermi-surface in the full BZ.
65 :
66 : integer :: nktot = -1
67 : !integer :: nkbz = -1
68 : ! Total number of k-points in the initial mesh. Used to compute integrals in the BZ.
69 :
70 : integer :: nkibz = -1
71 : ! Number of points in the IBZ
72 :
73 : integer :: bmin = -1, bmax = -1
74 : ! Min and max band index included in the calculation.
75 : ! Note that these values are obtained by taking the union over the k-points in the FS
76 : ! For each k-point, we usually have a different number of states crossing eF given by bstart_cnt_ibz.
77 :
78 : integer :: maxnb = -1
79 : ! Max number of bands on the FS i.e.: bmax - bmin + 1.
80 :
81 : integer :: eph_intmeth = 1
82 : ! Integration method.
83 : ! 1 for gaussian (including adaptive broadening if eph_fsmear is negative.
84 : ! |2| for tetrahedra.
85 : ! 2 for the optimized tetrahedron method.
86 : ! -2 for the linear tetrahedron method.
87 :
88 : integer :: nene = -1
89 : ! Number of chemical potential values used for inelastic integration.
90 :
91 : real(dp) :: eph_fsmear
92 : ! Gaussian broadening. Negative value activates adaptive gaussian broadening.
93 : ! See https://journals.aps.org/prb/pdf/10.1103/PhysRevB.92.075405
94 :
95 : real(dp) :: min_smear = tol9
96 : ! Used for the adaptive gaussian broadening: use min_smear if the broadening computed from the group velocity
97 : ! is smaller than this value to avoid divergences in the gaussian.
98 :
99 : real(dp) :: enemin
100 : ! Minimal chemical potential value used for inelastic integration.
101 :
102 : real(dp) :: deltaene
103 : ! Chemical potential increment for inelastic integration.
104 :
105 : type(krank_t) :: krank
106 : ! rank/inverse_rank pair for the k-points on the FS (kpts).
107 :
108 : integer,allocatable :: indkk_fs(:,:)
109 : ! (6, nkfs)
110 : ! Table giving the correspondence between a point in the FS-BZ and the IBZ:
111 : !
112 : ! indkk_fs(1,:) Mapping FS-BZ --> k-points in the IBZ (taken from ebands_t)
113 : ! indkk_fs(2,:) The index of the symmetry S such that kfs = tim_sign * S(k_ibz) + G0
114 : ! indkk_fs(3:5,:) The reduced components of G0.
115 : ! indkk_fs(6,:) 1 if time-reversal was used to generate the k-point, 0 otherwise
116 : !
117 : ! NB: The table is generated using the symrel convention and can therefore be used to
118 : ! symmetrize wavefunctions in k-space.
119 :
120 : integer,allocatable :: bstart_cnt_ibz(:,:)
121 : ! (2, nkibz)
122 : ! The indices of the bands within the energy window (depends on fsk)
123 : ! Note that we use the k-point index in the IBZ.
124 : !
125 : ! bstcnt(1, :) The index of the first band inside the energy window (start)
126 : ! bstcnt(2, :) Number of bands on the FS (count)
127 :
128 : real(dp) :: klatt(3, 3)
129 : ! Reciprocal of lattice vectors for full kpoint grid. Used by init_tetra
130 :
131 : real(dp) :: kmesh_cartvec(3,3)
132 : ! vectors defining the k-mesh (stored as column vector in Cartesian coords.
133 : ! Used to implement the adaptive gaussian broadening.
134 :
135 : real(dp),allocatable :: kpts(:,:)
136 : ! (3, nkfs)
137 : ! Reduced coordinates of the BZ k-points on the Fermi surface.
138 :
139 : real(dp),allocatable :: vk(:,:), vkq(:,:)
140 : ! (3, mnb)
141 : ! Velocities in cartesian coordinates. Used to implement the adaptive gaussian broadening
142 : ! Values are filled by the caller (e.g. phgamma) inside the loop over k-points.
143 :
144 : real(dp),allocatable :: tetra_wtk(:,:)
145 : ! (maxnb, nkibz)
146 : ! Weights for FS integration with tetrahedron method
147 : ! Note that the weights are dimensioned with nkibz
148 : ! (1, :) corresponds to %bmin
149 :
150 : real(dp),allocatable :: tetra_wtk_ene(:,:,:)
151 : ! (maxnb, nkibz, nene)
152 : ! Weights for FS integration with tetrahedron method for all chemical potentials
153 : ! Note that the weights are dimensioned with nkibz
154 : ! (1, :) corresponds to %bmin
155 :
156 : real(dp),allocatable :: dbldelta_tetra_weights_kfs(:,:,:)
157 : ! (maxnb, maxnb, nkfs)
158 : ! (1, 1, :) corresponds to %bmin
159 :
160 : contains
161 :
162 : procedure :: free => fstab_free
163 : ! Free memory.
164 :
165 : procedure :: findkg0 => fstab_findkg0
166 : ! Find the index of the k-point on the FS
167 :
168 : procedure :: get_dbldelta_weights => fstab_get_dbldelta_weights
169 : ! Compute weights for the integration of the double-delta.
170 :
171 : end type fstab_t
172 :
173 : public :: fstab_init ! Initialize the object.
174 : public :: fstab_print ! Print the object
175 : !!***
176 :
177 : !----------------------------------------------------------------------
178 :
179 : contains !============================================================
180 : !!***
181 :
182 : !!****f* m_fstab/fstab_free
183 : !! NAME
184 : !! fstab_free
185 : !!
186 : !! FUNCTION
187 : !! Free memory
188 : !!
189 : !! SOURCE
190 :
191 5 : subroutine fstab_free(fstab)
192 :
193 : !Arguments ------------------------------------
194 : class(fstab_t),intent(inout) :: fstab
195 : ! ************************************************************************
196 :
197 : ! integer
198 5 : ABI_SFREE(fstab%indkk_fs)
199 5 : ABI_SFREE(fstab%bstart_cnt_ibz)
200 :
201 : ! real
202 5 : ABI_SFREE(fstab%kpts)
203 5 : ABI_SFREE(fstab%vk)
204 5 : ABI_SFREE(fstab%vkq)
205 5 : ABI_SFREE(fstab%tetra_wtk)
206 5 : ABI_SFREE(fstab%tetra_wtk_ene)
207 5 : ABI_SFREE(fstab%dbldelta_tetra_weights_kfs)
208 :
209 : ! types
210 5 : call fstab%krank%free()
211 :
212 5 : end subroutine fstab_free
213 : !!***
214 :
215 : !----------------------------------------------------------------------
216 :
217 : !!****f* m_fstab/fstab_init
218 : !! NAME
219 : !! fstab_init
220 : !!
221 : !! FUNCTION
222 : !! Initialize the tables for the FS integration.
223 : !!
224 : !! INPUTS
225 : !! ebands<ebands_t>=The object describing the band structure.
226 : !! cryst<crystal_t>=Info on the crystalline structure.
227 : !! dtset:
228 : !! eph_fsewin=Energy window in Hartree. Only states in [efermi-fsewin, efermi+fsewin] are included.
229 : !! eph_intmeth=Flag selecting the integration method.
230 : !! kptrlatt(3,3)=k-point lattice specification
231 : !! nshiftk= number of shift vectors.
232 : !! shiftk(3,nshiftk)=shift vectors for k point generation
233 : !! comm=MPI communicator.
234 : !!
235 : !! OUTPUT
236 : !! fstab(nsppol)=Tables with the correspondence between points of the Fermi surface (FS)
237 : !! and the k-points in ebands_t.
238 : !! tetra: Tetrahedron object.
239 : !!
240 : !! TODO
241 : !! Use a different algorithm to select k-points if tetra. First compute tetra weights
242 : !! then k-points contributing to FS integrals are selected according to some threshold.
243 : !!
244 : !! SOURCE
245 :
246 15 : subroutine fstab_init(fstab, ebands, cryst, dtset, tetra, comm)
247 :
248 : !Arguments ------------------------------------
249 : !scalars
250 : type(ebands_t),intent(in) :: ebands
251 : type(crystal_t),intent(in) :: cryst
252 : type(dataset_type),intent(in) :: dtset
253 : type(htetra_t),intent(out) :: tetra
254 : integer,intent(in) :: comm
255 : !arrays
256 : type(fstab_t),target,intent(out) :: fstab(ebands%nsppol)
257 :
258 : !Local variables-------------------------------
259 : !scalars
260 : integer,parameter :: option0 = 0, brav1 = 1, bcorr0 = 0
261 : integer :: nkfs,spin,band,nband_k,i1,i2,ib,blow,ik_bz,ik_ibz,nkibz
262 : integer :: ik,mkpt,nkbz,ierr, nene,ifermi
263 : real(dp),parameter :: max_occ1 = one
264 : real(dp) :: elow,ehigh,ebis,enemin,enemax,deltaene,cpu,wall,gflops
265 : logical :: in_win
266 : character(len=80) :: errstr
267 : character(len=5000) :: msg
268 : type(fstab_t),pointer :: fs
269 5 : type(krank_t) :: krank
270 : !arrays
271 : integer :: kptrlatt(3,3)
272 5 : integer,allocatable :: full2ebands(:,:),bz2ibz(:), fs2bz(:),indkk(:,:) !,fs2ibz(:)
273 : real(dp) :: rlatt(3,3), klatt(3,3)
274 5 : real(dp),allocatable :: kbz(:,:), tmp_eigen(:),bdelta(:,:),btheta(:,:)
275 : ! *************************************************************************
276 :
277 5 : call cwtime(cpu, wall, gflops, "start")
278 :
279 80 : if (any(cryst%symrel(:,:,1) /= identity_3d) .and. any(abs(cryst%tnons(:,1)) > tol10) ) then
280 0 : ABI_ERROR('The first symmetry is not the identity operator!')
281 : end if
282 :
283 5 : nkibz = ebands%nkpt
284 65 : kptrlatt = dtset%kptrlatt
285 : !call kpts_ibz_from_kptrlatt(cryst, kptrlatt, ebands%kptopt, dtset%nshiftk, dtset%shiftk, &
286 : ! nkibz, kibz, wtk, nkbz, kbz, &
287 : ! new_kptrlatt, new_shiftk) ! Optional
288 :
289 : ! Call smpbz to get the full grid of k-points `kbz`
290 : ! brav1=1 is able to treat all bravais lattices (same option used in getkgrid)
291 : mkpt= kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3) &
292 : +kptrlatt(1,2)*kptrlatt(2,3)*kptrlatt(3,1) &
293 : +kptrlatt(1,3)*kptrlatt(2,1)*kptrlatt(3,2) &
294 : -kptrlatt(1,2)*kptrlatt(2,1)*kptrlatt(3,3) &
295 : -kptrlatt(1,3)*kptrlatt(2,2)*kptrlatt(3,1) &
296 5 : -kptrlatt(1,1)*kptrlatt(2,3)*kptrlatt(3,2)
297 :
298 15 : ABI_MALLOC(kbz, (3, mkpt))
299 :
300 5 : call smpbz(brav1, std_out, kptrlatt, mkpt, nkbz, dtset%nshiftk, option0, dtset%shiftk, kbz)
301 :
302 : ! Find correspondence BZ --> IBZ
303 : ! Note that we use symrel so these tables can be used to symmetrize wavefunctions.
304 15 : ABI_MALLOC(indkk, (6, nkbz))
305 :
306 5 : call krank%from_kptrlatt(ebands%nkpt, ebands%kptns, kptrlatt, compute_invrank=.False.)
307 :
308 5 : if (kpts_map("symrel", ebands%kptopt, cryst, krank, nkbz, kbz, indkk) /= 0) then
309 : write(msg, '(10a)' ) &
310 0 : 'The WFK file cannot be used to start the present calculation ',ch10, &
311 0 : 'It was asked that the wavefunctions be accurate, but',ch10, &
312 0 : 'at least one of the k points could not be generated from a symmetrical one.',ch10, &
313 0 : 'Action: check your WFK file and k-point input variables',ch10, &
314 0 : ' (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
315 0 : ABI_ERROR(msg)
316 : end if
317 :
318 5 : call krank%free()
319 5 : call cwtime_report(" fstab_init%krank", cpu, wall, gflops)
320 :
321 17930 : ABI_ICALLOC(full2ebands, (6, nkbz))
322 :
323 2565 : do ik_bz=1,nkbz
324 2560 : full2ebands(1, ik_bz) = indkk(1, ik_bz) ! ik_ibz
325 2560 : full2ebands(2, ik_bz) = indkk(2, ik_bz) ! isym
326 10240 : full2ebands(3:5, ik_bz) = indkk(3:5, ik_bz) ! g0
327 2565 : full2ebands(6, ik_bz) = indkk(6, ik_bz) ! itimrev
328 : end do
329 5 : ABI_FREE(indkk)
330 :
331 : ! Select only the k-points in the BZ that are sufficiently close to the FS.
332 : ! FIXME: Do not know why but lambda depends on eph_fsewin if gaussian.
333 5 : ABI_CHECK(dtset%eph_fsewin > tol12, "dtset%eph_fsewin < tol12")
334 5 : elow = ebands%fermie - dtset%eph_fsewin
335 5 : ehigh = ebands%fermie + dtset%eph_fsewin
336 5 : ebis = elow - abs(elow) * 0.001_dp
337 :
338 : ! Allocate workspace arrays.
339 : !ABI_MALLOC(fs2ibz, (nkbz))
340 15 : ABI_MALLOC(fs2bz, (nkbz))
341 :
342 10 : do spin=1,ebands%nsppol
343 5 : fs => fstab(spin)
344 5 : fs%spin = spin
345 15 : ABI_MALLOC(fs%bstart_cnt_ibz, (2, nkibz))
346 440 : fs%bstart_cnt_ibz = -1
347 :
348 : ! Find k-points on the FS associated to this spin.
349 5 : nkfs = 0
350 2565 : do ik_bz=1,nkbz
351 2560 : ik_ibz = full2ebands(1, ik_bz)
352 2560 : nband_k = ebands%nband(ik_ibz + (spin-1)*nkibz)
353 :
354 2560 : blow = bisect(ebands%eig(:nband_k, ik_ibz, spin), ebis)
355 2560 : if (blow == 0) blow = 1
356 : !if (blow == nband_k .or. blow == 0) cycle ! out of range
357 : !write(std_out,*)"here with blow: ", blow,nband_k
358 : !write(std_out,*)"eig_blow, eig_max, elow, ehigh:", &
359 : ! ebands%eig(blow, ik_ibz, spin), ebands%eig(nband_k, ik_ibz, spin), elow,ehigh
360 :
361 : in_win = .False.; i1 = huge(1); i2 = -1
362 14415 : do band=blow,nband_k
363 : !if (ebands%eig(band, ik_ibz, spin) > ehigh) exit
364 : !write(std_out,*)band, ebands%eig(band, ik_ibz, spin) >= elow, ebands%eig(band, ik_ibz, spin) <= ehigh
365 14415 : if (ebands%eig(band, ik_ibz, spin) >= elow .and. ebands%eig(band, ik_ibz, spin) <= ehigh) then
366 690 : in_win = .True.; i1 = min(i1, band); i2 = max(i2, band)
367 : end if
368 : end do
369 :
370 2565 : if (in_win) then
371 : ! Add this k-point and the corresponding bands.
372 : !write(std_out,*)"in win"
373 690 : nkfs = nkfs + 1
374 : !fs2ibz(nkfs) = ik_ibz
375 690 : fs2bz(nkfs) = ik_bz
376 750 : if (any(fs%bstart_cnt_ibz(:, ik_ibz) /= [-1, -1])) then
377 1980 : ABI_CHECK(all(fs%bstart_cnt_ibz(:, ik_ibz) == [i1, i2-i1+1]), "bstart_cnt_ibz!")
378 : end if
379 2070 : fs%bstart_cnt_ibz(:, ik_ibz) = [i1, i2-i1+1]
380 : end if
381 : end do ! ik_bz
382 :
383 : ! Build fstab_t for this spin.
384 5 : fs%nkibz = nkibz; fs%nkfs = nkfs; fs%nktot = nkbz
385 15 : ABI_MALLOC(fs%kpts, (3, nkfs))
386 15 : ABI_MALLOC(fs%indkk_fs, (6, nkfs))
387 695 : do ik=1,nkfs
388 : !ik_ibz = fs2ibz(ik)
389 690 : ik_bz = fs2bz(ik)
390 2760 : fs%kpts(:,ik) = kbz(:, ik_bz)
391 4835 : fs%indkk_fs(:, ik) = full2ebands(:, ik_bz)
392 : end do
393 :
394 : ! Define band indices enclosing states on the FS.
395 : ! Note that we need all k-points for a given band when computing weights with tetrahedron.
396 : ! This means that we have to be careful when selecting the weight associated to a given pair
397 : ! (band_kq, kq), (band_k, k).
398 : ! Then we have to rearrange the weights
399 5 : fs%bmin = huge(1); fs%bmax = -huge(1)
400 150 : do ik_ibz=1,nkibz
401 145 : if (fs%bstart_cnt_ibz(1, ik_ibz) /= -1) then
402 30 : fs%bmin = min(fs%bmin, fs%bstart_cnt_ibz(1,ik_ibz))
403 : end if
404 150 : if (fs%bstart_cnt_ibz(2, ik_ibz) /= -1) then
405 30 : fs%bmax = max(fs%bmax, fs%bstart_cnt_ibz(1,ik_ibz) + fs%bstart_cnt_ibz(2,ik_ibz) - 1)
406 : end if
407 : end do
408 :
409 : !write(std_out,*)"bmin, bmax for tetra: ",fs%bmin, fs%bmax
410 5 : ABI_CHECK(fs%bmin /= huge(1) .and. fs%bmax /= -huge(1), "No point on the Fermi surface!")
411 5 : fs%maxnb = fs%bmax - fs%bmin + 1
412 :
413 : ! DEBUG: use same number of bands for each k-point on the FS.
414 : !do ik_ibz=1,nkibz
415 : ! if (fs%bstart_cnt_ibz(1, ik_ibz) /= -1) then
416 : ! fs%bstart_cnt_ibz(1, ik_ibz) = fs%bmin
417 : ! fs%bstart_cnt_ibz(2, ik_ibz) = fs%maxnb
418 : ! end if
419 : !end do
420 :
421 55 : ABI_CALLOC(fs%vk, (3, fs%maxnb))
422 50 : ABI_CALLOC(fs%vkq, (3, fs%maxnb))
423 :
424 10 : call fs%krank%init(nkfs, fs%kpts)
425 : end do ! spin
426 :
427 5 : call cwtime_report(" fstab_init%fs_build:", cpu, wall, gflops)
428 :
429 : ! fix window around fermie for tetrahedron or gaussian weight calculation
430 : ! this is spin independent
431 5 : nene = 100 ! TODO: make this variable and maybe temperature dependent???
432 5 : deltaene = two * dtset%eph_fsewin / dble(nene-1)
433 5 : ifermi = int(nene / 2)
434 5 : enemin = ebands%fermie - dble(ifermi-1)*deltaene
435 5 : enemax = enemin + dble(nene-1)*deltaene
436 :
437 65 : rlatt = kptrlatt
438 5 : call matr3inv(rlatt, klatt)
439 :
440 : ! Setup FS integration
441 10 : do spin=1,ebands%nsppol
442 5 : fs => fstab(spin)
443 5 : fs%nene = nene
444 5 : fs%enemin = enemin
445 5 : fs%deltaene = deltaene
446 5 : fs%eph_intmeth = dtset%eph_intmeth
447 5 : fs%eph_fsmear = dtset%eph_fsmear
448 :
449 65 : fs%klatt = klatt
450 20 : fs%kmesh_cartvec(:, 1) = cryst%gprimd(:,1)*klatt(1,1) + cryst%gprimd(:,2)*klatt(2,1) + cryst%gprimd(:,3)*klatt(3,1)
451 20 : fs%kmesh_cartvec(:, 2) = cryst%gprimd(:,1)*klatt(1,2) + cryst%gprimd(:,2)*klatt(2,2) + cryst%gprimd(:,3)*klatt(3,2)
452 25 : fs%kmesh_cartvec(:, 3) = cryst%gprimd(:,1)*klatt(1,3) + cryst%gprimd(:,2)*klatt(2,3) + cryst%gprimd(:,3)*klatt(3,3)
453 : ! TODO: It seems that two_pi is not needed here!
454 : !fs%kmesh_cartvec = two_pi * fs%kmesh_cartvec
455 : !do i1=1,3
456 : ! write(std_out, *)"klatt:", klatt(:, i1)
457 : ! write(std_out, *)"gprimd:", cryst%gprimd(:, i1)
458 : ! write(std_out, *)"cartvec:", fs%kmesh_cartvec(:, i1)
459 : !end do
460 : end do
461 :
462 10 : ABI_MALLOC(bz2ibz, (nkbz))
463 2570 : bz2ibz = full2ebands(1, :)
464 5 : call tetra%init(bz2ibz, cryst%gprimd, klatt, kbz, nkbz, ebands%kptns, nkibz, ierr, errstr, comm)
465 5 : ABI_CHECK(ierr == 0, errstr)
466 5 : ABI_FREE(bz2ibz)
467 :
468 5 : if (abs(dtset%eph_intmeth) == 2) then
469 : ! TODO: compute weights on the fly to reduce memory? nene should be set to zero if not used!
470 6 : ABI_MALLOC(tmp_eigen, (nkibz))
471 6 : ABI_MALLOC(btheta, (nene, nkibz))
472 4 : ABI_MALLOC(bdelta, (nene, nkibz))
473 :
474 4 : do spin=1,ebands%nsppol
475 2 : fs => fstab(spin)
476 :
477 : ! Allocate tables used to store tetrahedron weights.
478 1942 : ABI_CALLOC(fs%dbldelta_tetra_weights_kfs, (fs%maxnb, fs%maxnb, fs%nkfs))
479 182 : ABI_CALLOC(fs%tetra_wtk, (fs%maxnb, nkibz))
480 17610 : ABI_CALLOC(fs%tetra_wtk_ene, (fs%maxnb, nkibz, fs%nene))
481 :
482 8 : do band=fs%bmin,fs%bmax
483 : ! Get the contribution of this band
484 124 : tmp_eigen = ebands%eig(band, :nkibz, spin)
485 :
486 : ! Calculate general integration weights at each irred kpoint
487 : ! as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
488 : call tetra%blochl_weights(tmp_eigen, enemin, enemax, max_occ1, fs%nene, nkibz, &
489 4 : bcorr0, btheta, bdelta, xmpi_comm_self)
490 :
491 : ! Save weights in the correct position.
492 4 : ib = band - fs%bmin + 1
493 122 : do ik_ibz=1,nkibz
494 116 : fs%tetra_wtk(ib, ik_ibz) = bdelta(ifermi, ik_ibz) * nkibz
495 11720 : fs%tetra_wtk_ene(ib, ik_ibz, 1:fs%nene) = bdelta(1:fs%nene, ik_ibz) * nkibz
496 : end do
497 : end do ! band
498 : end do ! spin
499 :
500 2 : ABI_FREE(tmp_eigen)
501 2 : ABI_FREE(btheta)
502 2 : ABI_FREE(bdelta)
503 : end if
504 :
505 : !ABI_FREE(fs2ibz)
506 5 : ABI_FREE(fs2bz)
507 5 : ABI_FREE(kbz)
508 5 : ABI_FREE(full2ebands)
509 :
510 5 : call cwtime_report(" fstab_init%fs_weights:", cpu, wall, gflops)
511 15 : if (xmpi_comm_rank(comm) == 0) call fstab_print(fstab, [std_out, ab_out])
512 :
513 25 : end subroutine fstab_init
514 : !!***
515 :
516 : !----------------------------------------------------------------------
517 :
518 : !!****f* m_fstab/fstab_findkg0
519 : !! NAME
520 : !! fstab_findkg0
521 : !!
522 : !! FUNCTION
523 : !! Return the index `ik_fs` of the k-point `kpt` in the FS-BZ. Return -1 if not found.
524 : !!
525 : !! INPUTS
526 : !! kpt(3)=K-point in reduced coordinates
527 : !!
528 : !! OUTPUT
529 : !! g0=Reciprocal lattice vector such that kpt = fstab%kpts(:, ik_fs) + g0
530 : !!
531 : !! SOURCE
532 :
533 9792 : integer function fstab_findkg0(fstab, kpt, g0) result(ik_fs)
534 :
535 : !Arguments ------------------------------------
536 : !scalars
537 : class(fstab_t),intent(in) :: fstab
538 : !arrays
539 : integer,intent(out) :: g0(3)
540 : real(dp),intent(in) :: kpt(3)
541 : ! *************************************************************************
542 :
543 9792 : ik_fs = fstab%krank%get_index(kpt)
544 9792 : if (ik_fs /= -1) then
545 17600 : g0 = nint(kpt - fstab%kpts(:, ik_fs))
546 : else
547 21568 : g0 = huge(1)
548 : end if
549 :
550 9792 : end function fstab_findkg0
551 : !!***
552 :
553 : !----------------------------------------------------------------------
554 :
555 : !!****f* m_fstab/fstab_get_dbldelta_weights
556 : !! NAME
557 : !! fstab_get_dbldelta_weights
558 : !!
559 : !! FUNCTION
560 : !! Return the weights for the integration of the double-delta on the Fermi-surface
561 : !!
562 : !! INPUTS
563 : !! ebands<ebands_type>=GS band structure.
564 : !! ik_ibz=Index of the k-point in the IBZ
565 : !! spin=Spin index
566 : !! nesting=Used in tetra. Set to 1 if tetra Weights cannot be computed with Tetrahedron due to nesting.
567 : !! In this case, we fallback to adaptive gaussian.
568 : !!
569 : !! OUTPUT
570 : !! wtk(fs%maxnb, fs%maxnb)=Weights for FS integration.
571 : !!
572 : !! SOURCE
573 :
574 2864 : subroutine fstab_get_dbldelta_weights(fs, ebands, ik_fs, ik_ibz, ikq_ibz, spin, nesting, wtk)
575 :
576 : !Arguments ------------------------------------
577 : !scalars
578 : class(fstab_t),intent(in) :: fs
579 : class(ebands_t),intent(in) :: ebands
580 : integer,intent(in) :: ik_fs, ik_ibz, ikq_ibz, spin, nesting
581 : !arrays
582 : real(dp),intent(out) :: wtk(fs%maxnb, fs%maxnb)
583 :
584 : !Local variables-------------------------------
585 : !scalars
586 : integer :: bstart_k, nband_k, bstart_kq, nband_kq, ib1, band1, ib2, band2, ii
587 : logical :: use_adaptive
588 : real(dp) :: g1, g2, sigma, abc(3)
589 : ! *************************************************************************
590 :
591 2864 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz)
592 2864 : ABI_CHECK(nband_k >= 1 .and. nband_k <= fs%maxnb, "Wrong nband_k")
593 :
594 2864 : bstart_kq = fs%bstart_cnt_ibz(1, ikq_ibz); nband_kq = fs%bstart_cnt_ibz(2, ikq_ibz)
595 2864 : ABI_CHECK(nband_kq >= 1 .and. nband_kq <= fs%maxnb, "Wrong nband_kq")
596 :
597 20048 : wtk = zero
598 2864 : if (fs%eph_intmeth == 1 .or. nesting /= 0) then
599 : ! Gaussian method: constant or adaptive method from group velocities if eph_fsmear is negative.
600 1536 : sigma = fs%eph_fsmear
601 1536 : use_adaptive = fs%eph_fsmear < zero .or. abs(fs%eph_intmeth) == 2
602 3072 : do ib2=1,nband_k
603 1536 : band2 = ib2 + bstart_k - 1
604 1536 : if (use_adaptive) then
605 4464 : do ii=1,3
606 14508 : abc(ii) = abs(dot_product(fs%vk(:, ib2), fs%kmesh_cartvec(:,ii)))
607 : end do
608 5580 : sigma = max(maxval(abc), fs%min_smear)
609 : !write(std_out, *)"sigma:", sigma * Ha_eV
610 : end if
611 1536 : g2 = gaussian(ebands%eig(band2, ik_ibz, spin) - ebands%fermie, sigma)
612 4608 : do ib1=1,nband_kq
613 1536 : band1 = ib1 + bstart_kq - 1
614 1536 : if (use_adaptive) then
615 4464 : do ii=1,3
616 14508 : abc(ii) = abs(dot_product(fs%vkq(:, ib1), fs%kmesh_cartvec(:,ii)))
617 : end do
618 5580 : sigma = max(maxval(abc), fs%min_smear)
619 : end if
620 1536 : g1 = gaussian(ebands%eig(band1, ikq_ibz, spin) - ebands%fermie, sigma)
621 3072 : wtk(ib1, ib2) = (g1 * g2) / fs%nktot
622 : end do
623 : end do
624 :
625 1328 : else if (abs(fs%eph_intmeth) == 2) then
626 : ! Tetrahedron method. Copy weights in the correct position.
627 2656 : do ib2=1,nband_k
628 1328 : band2 = ib2 + bstart_k - fs%bmin
629 3984 : do ib1=1,nband_kq
630 1328 : band1 = ib1 + bstart_kq - fs%bmin
631 : ! libtetrabz_dbldelta seems to report weights in this order.
632 2656 : wtk(ib1, ib2) = fs%dbldelta_tetra_weights_kfs(band1, band2, ik_fs)
633 : end do
634 : end do
635 :
636 : else
637 0 : ABI_ERROR(sjoin("Wrong integration method:", itoa(fs%eph_intmeth)))
638 : end if
639 :
640 2864 : end subroutine fstab_get_dbldelta_weights
641 : !!***
642 :
643 : !----------------------------------------------------------------------
644 :
645 : !!****f* m_fstab/fstab_print
646 : !! NAME
647 : !! fstab_print
648 : !!
649 : !! FUNCTION
650 : !! Print info on the object.
651 : !!
652 : !! INPUTS
653 : !! units=unit numbers for output
654 : !! [header]=Header string.
655 : !! [prtvol]=verbosity level
656 : !!
657 : !! SOURCE
658 :
659 10 : subroutine fstab_print(fstab, units, header, prtvol)
660 :
661 : !Arguments ------------------------------------
662 : !scalars
663 : class(fstab_t),target,intent(in) :: fstab(:)
664 : integer,intent(in) :: units(:)
665 : character(len=*),optional,intent(in) :: header
666 : integer,optional,intent(in) :: prtvol
667 :
668 : !Local variables-------------------------------
669 : !scalars
670 : integer :: my_prtvol,spin
671 : character(len=5000) :: msg
672 : ! *************************************************************************
673 :
674 5 : my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
675 :
676 5 : msg = ' ==== Fermi surface info ==== '
677 5 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
678 5 : call wrtout(units, msg)
679 :
680 8 : select case (fstab(1)%eph_intmeth)
681 : case (1)
682 3 : if (fstab(1)%eph_fsmear > zero) then
683 1 : write(msg,"(a,f5.1,a)")" FS integration done with gaussian method and broadening:", &
684 2 : fstab(1)%eph_fsmear * Ha_eV, " (meV)"
685 : else
686 2 : write(msg,"(a)")" FS integration done with adaptive gaussian method"
687 : end if
688 : case (2)
689 2 : write(msg,"(a)")" FS integration done with tetrahedron method"
690 : case (-2)
691 0 : write(msg,"(a)")" FS integration done with optimized tetrahedron method"
692 : case default
693 5 : ABI_ERROR(sjoin("Invalid value for eph_intmeth:", itoa(fstab(1)%eph_intmeth)))
694 : end select
695 :
696 5 : call wrtout(units, msg)
697 :
698 5 : write(msg,"(a,i0)")" Total number of k-points in the full mesh: ",fstab(1)%nktot
699 5 : call wrtout(units, msg)
700 : !write(msg,"(a,f5.1)")" Energy window: ",fstab(1)%eph_fsewin * Ha_eV, " (eV)
701 : !call wrtout(units, msg)
702 :
703 10 : do spin=1,size(fstab)
704 5 : associate (fs => fstab(spin))
705 5 : write(msg,"(a,i0)")" For spin: ",spin
706 5 : call wrtout(units, msg)
707 : write(msg,"(a,i0,a,f5.1,a)") &
708 5 : " Number of BZ k-points close to the Fermi surface: ",fs%nkfs," [", (100.0_dp * fs%nkfs) / fs%nktot, " %]"
709 5 : call wrtout(units, msg)
710 5 : write(msg,"(a,i0)")" Maximum number of bands crossing the Fermi level: ",fs%maxnb
711 5 : call wrtout(units, msg)
712 150 : write(msg,"(2(a,i0))")" min band: ", minval(fs%bstart_cnt_ibz(1,:), mask=fs%bstart_cnt_ibz(1,:) /= -1)
713 5 : call wrtout(units, msg)
714 5 : write(msg,"(2(a,i0))")" Max band: ", maxval(fs%bstart_cnt_ibz(1,:) + fs%bstart_cnt_ibz(2,:) - 1, &
715 155 : mask=fs%bstart_cnt_ibz(1,:) /= -1)
716 10 : call wrtout(units, msg)
717 : end associate
718 : end do
719 :
720 5 : end subroutine fstab_print
721 : !!***
722 :
723 0 : end module m_fstab
724 : !!***
|
