Branch data 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, 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 : : mkpt= kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3) &
291 : : +kptrlatt(1,2)*kptrlatt(2,3)*kptrlatt(3,1) &
292 : : +kptrlatt(1,3)*kptrlatt(2,1)*kptrlatt(3,2) &
293 : : -kptrlatt(1,2)*kptrlatt(2,1)*kptrlatt(3,3) &
294 : : -kptrlatt(1,3)*kptrlatt(2,2)*kptrlatt(3,1) &
295 : 5 : -kptrlatt(1,1)*kptrlatt(2,3)*kptrlatt(3,2)
296 : :
297 [ - + + - : 15 : ABI_MALLOC(kbz, (3, mkpt))
+ - + - -
+ - + ]
298 : :
299 : 5 : call smpbz(std_out, kptrlatt, mkpt, nkbz, dtset%nshiftk, option0, dtset%shiftk, kbz)
300 : :
301 : : ! Find correspondence BZ --> IBZ
302 : : ! Note that we use symrel so these tables can be used to symmetrize wavefunctions.
303 [ - + + - : 15 : ABI_MALLOC(indkk, (6, nkbz))
+ - + - -
+ - + ]
304 : :
305 : 5 : call krank%from_kptrlatt(ebands%nkpt, ebands%kptns, kptrlatt, compute_invrank=.False.)
306 : :
307 [ - + ]: 5 : if (kpts_map("symrel", ebands%kptopt, cryst, krank, nkbz, kbz, indkk) /= 0) then
308 : : write(msg, '(10a)' ) &
309 : 0 : 'The WFK file cannot be used to start the present calculation ',ch10, &
310 : 0 : 'It was asked that the wavefunctions be accurate, but',ch10, &
311 : 0 : 'at least one of the k points could not be generated from a symmetrical one.',ch10, &
312 : 0 : 'Action: check your WFK file and k-point input variables',ch10, &
313 : 0 : ' (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
314 : 0 : ABI_ERROR(msg)
315 : : end if
316 : :
317 : 5 : call krank%free()
318 : 5 : call cwtime_report(" fstab_init%krank", cpu, wall, gflops)
319 : :
320 [ - + + - : 17930 : ABI_ICALLOC(full2ebands, (6, nkbz))
- + - + +
+ + + ]
321 : :
322 [ + + ]: 2565 : do ik_bz=1,nkbz
323 : 2560 : full2ebands(1, ik_bz) = indkk(1, ik_bz) ! ik_ibz
324 : 2560 : full2ebands(2, ik_bz) = indkk(2, ik_bz) ! isym
325 [ + + ]: 10240 : full2ebands(3:5, ik_bz) = indkk(3:5, ik_bz) ! g0
326 : 2565 : full2ebands(6, ik_bz) = indkk(6, ik_bz) ! itimrev
327 : : end do
328 [ - + ]: 5 : ABI_FREE(indkk)
329 : :
330 : : ! Select only the k-points in the BZ that are sufficiently close to the FS.
331 : : ! FIXME: Do not know why but lambda depends on eph_fsewin if gaussian.
332 [ - + ]: 5 : ABI_CHECK(dtset%eph_fsewin > tol12, "dtset%eph_fsewin < tol12")
333 : 5 : elow = ebands%fermie - dtset%eph_fsewin
334 : 5 : ehigh = ebands%fermie + dtset%eph_fsewin
335 : 5 : ebis = elow - abs(elow) * 0.001_dp
336 : :
337 : : ! Allocate workspace arrays.
338 : : !ABI_MALLOC(fs2ibz, (nkbz))
339 [ - + + - : 15 : ABI_MALLOC(fs2bz, (nkbz))
+ - + - -
+ - + ]
340 : :
341 [ + + ]: 10 : do spin=1,ebands%nsppol
342 : 5 : fs => fstab(spin)
343 : 5 : fs%spin = spin
344 [ - + + - : 15 : ABI_MALLOC(fs%bstart_cnt_ibz, (2, nkibz))
+ - + - -
+ - + -
+ ]
345 [ + + + + ]: 440 : fs%bstart_cnt_ibz = -1
346 : :
347 : : ! Find k-points on the FS associated to this spin.
348 : 5 : nkfs = 0
349 [ + + ]: 2565 : do ik_bz=1,nkbz
350 : 2560 : ik_ibz = full2ebands(1, ik_bz)
351 : 2560 : nband_k = ebands%nband(ik_ibz + (spin-1)*nkibz)
352 : :
353 : 2560 : blow = bisect(ebands%eig(:nband_k, ik_ibz, spin), ebis)
354 [ - + ]: 2560 : if (blow == 0) blow = 1
355 : : !if (blow == nband_k .or. blow == 0) cycle ! out of range
356 : : !write(std_out,*)"here with blow: ", blow,nband_k
357 : : !write(std_out,*)"eig_blow, eig_max, elow, ehigh:", &
358 : : ! ebands%eig(blow, ik_ibz, spin), ebands%eig(nband_k, ik_ibz, spin), elow,ehigh
359 : :
360 : : in_win = .False.; i1 = huge(1); i2 = -1
361 [ + + ]: 14415 : do band=blow,nband_k
362 : : !if (ebands%eig(band, ik_ibz, spin) > ehigh) exit
363 : : !write(std_out,*)band, ebands%eig(band, ik_ibz, spin) >= elow, ebands%eig(band, ik_ibz, spin) <= ehigh
364 [ + + + + ]: 14415 : if (ebands%eig(band, ik_ibz, spin) >= elow .and. ebands%eig(band, ik_ibz, spin) <= ehigh) then
365 : 690 : in_win = .True.; i1 = min(i1, band); i2 = max(i2, band)
366 : : end if
367 : : end do
368 : :
369 [ + + ]: 2565 : if (in_win) then
370 : : ! Add this k-point and the corresponding bands.
371 : : !write(std_out,*)"in win"
372 : 690 : nkfs = nkfs + 1
373 : : !fs2ibz(nkfs) = ik_ibz
374 : 690 : fs2bz(nkfs) = ik_bz
375 [ + + + + : 750 : if (any(fs%bstart_cnt_ibz(:, ik_ibz) /= [-1, -1])) then
+ + ]
376 [ + + + - : 1980 : ABI_CHECK(all(fs%bstart_cnt_ibz(:, ik_ibz) == [i1, i2-i1+1]), "bstart_cnt_ibz!")
- + ]
377 : : end if
378 [ + + ]: 2070 : fs%bstart_cnt_ibz(:, ik_ibz) = [i1, i2-i1+1]
379 : : end if
380 : : end do ! ik_bz
381 : :
382 : : ! Build fstab_t for this spin.
383 : 5 : fs%nkibz = nkibz; fs%nkfs = nkfs; fs%nktot = nkbz
384 [ - + + - : 15 : ABI_MALLOC(fs%kpts, (3, nkfs))
+ - + - -
+ - + -
+ ]
385 [ - + + - : 15 : ABI_MALLOC(fs%indkk_fs, (6, nkfs))
+ - + - -
+ - + -
+ ]
386 [ + + ]: 695 : do ik=1,nkfs
387 : : !ik_ibz = fs2ibz(ik)
388 : 690 : ik_bz = fs2bz(ik)
389 [ + + ]: 2760 : fs%kpts(:,ik) = kbz(:, ik_bz)
390 [ + + ]: 4835 : fs%indkk_fs(:, ik) = full2ebands(:, ik_bz)
391 : : end do
392 : :
393 : : ! Define band indices enclosing states on the FS.
394 : : ! Note that we need all k-points for a given band when computing weights with tetrahedron.
395 : : ! This means that we have to be careful when selecting the weight associated to a given pair
396 : : ! (band_kq, kq), (band_k, k).
397 : : ! Then we have to rearrange the weights
398 : 5 : fs%bmin = huge(1); fs%bmax = -huge(1)
399 [ + + ]: 150 : do ik_ibz=1,nkibz
400 [ + + ]: 145 : if (fs%bstart_cnt_ibz(1, ik_ibz) /= -1) then
401 : 30 : fs%bmin = min(fs%bmin, fs%bstart_cnt_ibz(1,ik_ibz))
402 : : end if
403 [ + + ]: 150 : if (fs%bstart_cnt_ibz(2, ik_ibz) /= -1) then
404 : 30 : fs%bmax = max(fs%bmax, fs%bstart_cnt_ibz(1,ik_ibz) + fs%bstart_cnt_ibz(2,ik_ibz) - 1)
405 : : end if
406 : : end do
407 : :
408 : : !write(std_out,*)"bmin, bmax for tetra: ",fs%bmin, fs%bmax
409 [ + - - + ]: 5 : ABI_CHECK(fs%bmin /= huge(1) .and. fs%bmax /= -huge(1), "No point on the Fermi surface!")
410 : 5 : fs%maxnb = fs%bmax - fs%bmin + 1
411 : :
412 : : ! DEBUG: use same number of bands for each k-point on the FS.
413 : : !do ik_ibz=1,nkibz
414 : : ! if (fs%bstart_cnt_ibz(1, ik_ibz) /= -1) then
415 : : ! fs%bstart_cnt_ibz(1, ik_ibz) = fs%bmin
416 : : ! fs%bstart_cnt_ibz(2, ik_ibz) = fs%maxnb
417 : : ! end if
418 : : !end do
419 : :
420 [ - + + - : 55 : ABI_CALLOC(fs%vk, (3, fs%maxnb))
+ - + - -
+ - + - +
+ + + + ]
421 [ - + + - : 50 : ABI_CALLOC(fs%vkq, (3, fs%maxnb))
- + - + -
+ + + +
+ ]
422 : :
423 : 10 : call fs%krank%init(nkfs, fs%kpts)
424 : : end do ! spin
425 : :
426 : 5 : call cwtime_report(" fstab_init%fs_build:", cpu, wall, gflops)
427 : :
428 : : ! fix window around fermie for tetrahedron or gaussian weight calculation
429 : : ! this is spin independent
430 : 5 : nene = 100 ! TODO: make this variable and maybe temperature dependent???
431 : 5 : deltaene = two * dtset%eph_fsewin / dble(nene-1)
432 : 5 : ifermi = int(nene / 2)
433 : 5 : enemin = ebands%fermie - dble(ifermi-1)*deltaene
434 : 5 : enemax = enemin + dble(nene-1)*deltaene
435 : :
436 [ + + + + ]: 65 : rlatt = kptrlatt
437 : 5 : call matr3inv(rlatt, klatt)
438 : :
439 : : ! Setup FS integration
440 [ + + ]: 10 : do spin=1,ebands%nsppol
441 : 5 : fs => fstab(spin)
442 : 5 : fs%nene = nene
443 : 5 : fs%enemin = enemin
444 : 5 : fs%deltaene = deltaene
445 : 5 : fs%eph_intmeth = dtset%eph_intmeth
446 : 5 : fs%eph_fsmear = dtset%eph_fsmear
447 : :
448 [ + + + + ]: 65 : fs%klatt = klatt
449 [ + + ]: 20 : fs%kmesh_cartvec(:, 1) = cryst%gprimd(:,1)*klatt(1,1) + cryst%gprimd(:,2)*klatt(2,1) + cryst%gprimd(:,3)*klatt(3,1)
450 [ + + ]: 20 : fs%kmesh_cartvec(:, 2) = cryst%gprimd(:,1)*klatt(1,2) + cryst%gprimd(:,2)*klatt(2,2) + cryst%gprimd(:,3)*klatt(3,2)
451 [ + + ]: 25 : fs%kmesh_cartvec(:, 3) = cryst%gprimd(:,1)*klatt(1,3) + cryst%gprimd(:,2)*klatt(2,3) + cryst%gprimd(:,3)*klatt(3,3)
452 : : ! TODO: It seems that two_pi is not needed here!
453 : : !fs%kmesh_cartvec = two_pi * fs%kmesh_cartvec
454 : : !do i1=1,3
455 : : ! write(std_out, *)"klatt:", klatt(:, i1)
456 : : ! write(std_out, *)"gprimd:", cryst%gprimd(:, i1)
457 : : ! write(std_out, *)"cartvec:", fs%kmesh_cartvec(:, i1)
458 : : !end do
459 : : end do
460 : :
461 [ - + + - : 10 : ABI_MALLOC(bz2ibz, (nkbz))
- + - + ]
462 [ - + + + ]: 2570 : bz2ibz = full2ebands(1, :)
463 : 5 : call tetra%init(bz2ibz, cryst%gprimd, klatt, kbz, nkbz, ebands%kptns, nkibz, ierr, errstr, comm)
464 [ - + ]: 5 : ABI_CHECK(ierr == 0, errstr)
465 [ - + ]: 5 : ABI_FREE(bz2ibz)
466 : :
467 [ + + ]: 5 : if (abs(dtset%eph_intmeth) == 2) then
468 : : ! TODO: compute weights on the fly to reduce memory? nene should be set to zero if not used!
469 [ - + + - : 6 : ABI_MALLOC(tmp_eigen, (nkibz))
+ - + - -
+ - + ]
470 [ - + + - : 6 : ABI_MALLOC(btheta, (nene, nkibz))
+ - + - -
+ - + ]
471 [ - + + - : 4 : ABI_MALLOC(bdelta, (nene, nkibz))
- + - + ]
472 : :
473 [ + + ]: 4 : do spin=1,ebands%nsppol
474 : 2 : fs => fstab(spin)
475 : :
476 : : ! Allocate tables used to store tetrahedron weights.
477 [ - + + - : 1942 : ABI_CALLOC(fs%dbldelta_tetra_weights_kfs, (fs%maxnb, fs%maxnb, fs%nkfs))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
478 [ - + + - : 182 : ABI_CALLOC(fs%tetra_wtk, (fs%maxnb, nkibz))
- + + - +
- + - - +
- + - + +
+ + + ]
479 [ - + + - : 17610 : ABI_CALLOC(fs%tetra_wtk_ene, (fs%maxnb, nkibz, fs%nene))
- + + - -
+ + - + -
+ - - + -
+ - + + +
+ + + + ]
480 : :
481 [ + + ]: 8 : do band=fs%bmin,fs%bmax
482 : : ! Get the contribution of this band
483 [ + - - + : 124 : tmp_eigen = ebands%eig(band, :nkibz, spin)
- - + + ]
484 : :
485 : : ! Calculate general integration weights at each irred kpoint
486 : : ! as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
487 : : call tetra%blochl_weights(tmp_eigen, enemin, enemax, max_occ1, fs%nene, nkibz, &
488 : 4 : bcorr0, btheta, bdelta, xmpi_comm_self)
489 : :
490 : : ! Save weights in the correct position.
491 : 4 : ib = band - fs%bmin + 1
492 [ + + ]: 122 : do ik_ibz=1,nkibz
493 : 116 : fs%tetra_wtk(ib, ik_ibz) = bdelta(ifermi, ik_ibz) * nkibz
494 [ + + ]: 11720 : fs%tetra_wtk_ene(ib, ik_ibz, 1:fs%nene) = bdelta(1:fs%nene, ik_ibz) * nkibz
495 : : end do
496 : : end do ! band
497 : : end do ! spin
498 : :
499 [ - + ]: 2 : ABI_FREE(tmp_eigen)
500 [ - + ]: 2 : ABI_FREE(btheta)
501 [ - + ]: 2 : ABI_FREE(bdelta)
502 : : end if
503 : :
504 : : !ABI_FREE(fs2ibz)
505 [ - + ]: 5 : ABI_FREE(fs2bz)
506 [ - + ]: 5 : ABI_FREE(kbz)
507 [ - + ]: 5 : ABI_FREE(full2ebands)
508 : :
509 : 5 : call cwtime_report(" fstab_init%fs_weights:", cpu, wall, gflops)
510 [ + - + + ]: 15 : if (xmpi_comm_rank(comm) == 0) call fstab_print(fstab, [std_out, ab_out])
511 : :
512 [ - + - + : 25 : end subroutine fstab_init
- + ]
513 : : !!***
514 : :
515 : : !----------------------------------------------------------------------
516 : :
517 : : !!****f* m_fstab/fstab_findkg0
518 : : !! NAME
519 : : !! fstab_findkg0
520 : : !!
521 : : !! FUNCTION
522 : : !! Return the index `ik_fs` of the k-point `kpt` in the FS-BZ. Return -1 if not found.
523 : : !!
524 : : !! INPUTS
525 : : !! kpt(3)=K-point in reduced coordinates
526 : : !!
527 : : !! OUTPUT
528 : : !! g0=Reciprocal lattice vector such that kpt = fstab%kpts(:, ik_fs) + g0
529 : : !!
530 : : !! SOURCE
531 : :
532 : 9792 : integer function fstab_findkg0(fstab, kpt, g0) result(ik_fs)
533 : :
534 : : !Arguments ------------------------------------
535 : : !scalars
536 : : class(fstab_t),intent(in) :: fstab
537 : : !arrays
538 : : integer,intent(out) :: g0(3)
539 : : real(dp),intent(in) :: kpt(3)
540 : : ! *************************************************************************
541 : :
542 : 9792 : ik_fs = fstab%krank%get_index(kpt)
543 [ + + ]: 9792 : if (ik_fs /= -1) then
544 [ + + ]: 17600 : g0 = nint(kpt - fstab%kpts(:, ik_fs))
545 : : else
546 [ + + ]: 21568 : g0 = huge(1)
547 : : end if
548 : :
549 : 9792 : end function fstab_findkg0
550 : : !!***
551 : :
552 : : !----------------------------------------------------------------------
553 : :
554 : : !!****f* m_fstab/fstab_get_dbldelta_weights
555 : : !! NAME
556 : : !! fstab_get_dbldelta_weights
557 : : !!
558 : : !! FUNCTION
559 : : !! Return the weights for the integration of the double-delta on the Fermi-surface
560 : : !!
561 : : !! INPUTS
562 : : !! ebands<ebands_type>=GS band structure.
563 : : !! ik_ibz=Index of the k-point in the IBZ
564 : : !! spin=Spin index
565 : : !! nesting=Used in tetra. Set to 1 if tetra Weights cannot be computed with Tetrahedron due to nesting.
566 : : !! In this case, we fallback to adaptive gaussian.
567 : : !!
568 : : !! OUTPUT
569 : : !! wtk(fs%maxnb, fs%maxnb)=Weights for FS integration.
570 : : !!
571 : : !! SOURCE
572 : :
573 : 2864 : subroutine fstab_get_dbldelta_weights(fs, ebands, ik_fs, ik_ibz, ikq_ibz, spin, nesting, wtk)
574 : :
575 : : !Arguments ------------------------------------
576 : : !scalars
577 : : class(fstab_t),intent(in) :: fs
578 : : class(ebands_t),intent(in) :: ebands
579 : : integer,intent(in) :: ik_fs, ik_ibz, ikq_ibz, spin, nesting
580 : : !arrays
581 : : real(dp),intent(out) :: wtk(fs%maxnb, fs%maxnb)
582 : :
583 : : !Local variables-------------------------------
584 : : !scalars
585 : : integer :: bstart_k, nband_k, bstart_kq, nband_kq, ib1, band1, ib2, band2, ii
586 : : logical :: use_adaptive
587 : : real(dp) :: g1, g2, sigma, abc(3)
588 : : ! *************************************************************************
589 : :
590 : 2864 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz)
591 [ + - - + ]: 2864 : ABI_CHECK(nband_k >= 1 .and. nband_k <= fs%maxnb, "Wrong nband_k")
592 : :
593 : 2864 : bstart_kq = fs%bstart_cnt_ibz(1, ikq_ibz); nband_kq = fs%bstart_cnt_ibz(2, ikq_ibz)
594 [ + - - + ]: 2864 : ABI_CHECK(nband_kq >= 1 .and. nband_kq <= fs%maxnb, "Wrong nband_kq")
595 : :
596 [ + + + + ]: 20048 : wtk = zero
597 [ + + + + ]: 2864 : if (fs%eph_intmeth == 1 .or. nesting /= 0) then
598 : : ! Gaussian method: constant or adaptive method from group velocities if eph_fsmear is negative.
599 : 1536 : sigma = fs%eph_fsmear
600 [ + + + + ]: 1536 : use_adaptive = fs%eph_fsmear < zero .or. abs(fs%eph_intmeth) == 2
601 [ + + ]: 3072 : do ib2=1,nband_k
602 : 1536 : band2 = ib2 + bstart_k - 1
603 [ + + ]: 1536 : if (use_adaptive) then
604 [ + + ]: 4464 : do ii=1,3
605 [ + + ]: 14508 : abc(ii) = abs(dot_product(fs%vk(:, ib2), fs%kmesh_cartvec(:,ii)))
606 : : end do
607 [ - + + - : 5580 : sigma = max(maxval(abc), fs%min_smear)
+ + ]
608 : : !write(std_out, *)"sigma:", sigma * Ha_eV
609 : : end if
610 : 1536 : g2 = gaussian(ebands%eig(band2, ik_ibz, spin) - ebands%fermie, sigma)
611 [ + + ]: 4608 : do ib1=1,nband_kq
612 : 1536 : band1 = ib1 + bstart_kq - 1
613 [ + + ]: 1536 : if (use_adaptive) then
614 [ + + ]: 4464 : do ii=1,3
615 [ + + ]: 14508 : abc(ii) = abs(dot_product(fs%vkq(:, ib1), fs%kmesh_cartvec(:,ii)))
616 : : end do
617 [ - + + - : 5580 : sigma = max(maxval(abc), fs%min_smear)
+ + ]
618 : : end if
619 : 1536 : g1 = gaussian(ebands%eig(band1, ikq_ibz, spin) - ebands%fermie, sigma)
620 : 3072 : wtk(ib1, ib2) = (g1 * g2) / fs%nktot
621 : : end do
622 : : end do
623 : :
624 [ + - ]: 1328 : else if (abs(fs%eph_intmeth) == 2) then
625 : : ! Tetrahedron method. Copy weights in the correct position.
626 [ + + ]: 2656 : do ib2=1,nband_k
627 : 1328 : band2 = ib2 + bstart_k - fs%bmin
628 [ + + ]: 3984 : do ib1=1,nband_kq
629 : 1328 : band1 = ib1 + bstart_kq - fs%bmin
630 : : ! libtetrabz_dbldelta seems to report weights in this order.
631 : 2656 : wtk(ib1, ib2) = fs%dbldelta_tetra_weights_kfs(band1, band2, ik_fs)
632 : : end do
633 : : end do
634 : :
635 : : else
636 : 0 : ABI_ERROR(sjoin("Wrong integration method:", itoa(fs%eph_intmeth)))
637 : : end if
638 : :
639 : 2864 : end subroutine fstab_get_dbldelta_weights
640 : : !!***
641 : :
642 : : !----------------------------------------------------------------------
643 : :
644 : : !!****f* m_fstab/fstab_print
645 : : !! NAME
646 : : !! fstab_print
647 : : !!
648 : : !! FUNCTION
649 : : !! Print info on the object.
650 : : !!
651 : : !! INPUTS
652 : : !! units=unit numbers for output
653 : : !! [header]=Header string.
654 : : !! [prtvol]=verbosity level
655 : : !!
656 : : !! SOURCE
657 : :
658 : 10 : subroutine fstab_print(fstab, units, header, prtvol)
659 : :
660 : : !Arguments ------------------------------------
661 : : !scalars
662 : : class(fstab_t),target,intent(in) :: fstab(:)
663 : : integer,intent(in) :: units(:)
664 : : character(len=*),optional,intent(in) :: header
665 : : integer,optional,intent(in) :: prtvol
666 : :
667 : : !Local variables-------------------------------
668 : : !scalars
669 : : integer :: my_prtvol,spin
670 : : character(len=5000) :: msg
671 : : ! *************************************************************************
672 : :
673 : 5 : my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
674 : :
675 : 5 : msg = ' ==== Fermi surface info ==== '
676 [ - + - - : 5 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
- - ]
677 : 5 : call wrtout(units, msg)
678 : :
679 : 8 : select case (fstab(1)%eph_intmeth)
680 : : case (1)
681 [ + + ]: 3 : if (fstab(1)%eph_fsmear > zero) then
682 : 1 : write(msg,"(a,f5.1,a)")" FS integration done with gaussian method and broadening:", &
683 : 2 : fstab(1)%eph_fsmear * Ha_eV, " (meV)"
684 : : else
685 : 2 : write(msg,"(a)")" FS integration done with adaptive gaussian method"
686 : : end if
687 : : case (2)
688 : 2 : write(msg,"(a)")" FS integration done with tetrahedron method"
689 : : case (-2)
690 : 0 : write(msg,"(a)")" FS integration done with optimized tetrahedron method"
691 : : case default
692 [ + + - - ]: 5 : ABI_ERROR(sjoin("Invalid value for eph_intmeth:", itoa(fstab(1)%eph_intmeth)))
693 : : end select
694 : :
695 : 5 : call wrtout(units, msg)
696 : :
697 : 5 : write(msg,"(a,i0)")" Total number of k-points in the full mesh: ",fstab(1)%nktot
698 : 5 : call wrtout(units, msg)
699 : : !write(msg,"(a,f5.1)")" Energy window: ",fstab(1)%eph_fsewin * Ha_eV, " (eV)
700 : : !call wrtout(units, msg)
701 : :
702 [ + + ]: 10 : do spin=1,size(fstab)
703 : 5 : associate (fs => fstab(spin))
704 : 5 : write(msg,"(a,i0)")" For spin: ",spin
705 : 5 : call wrtout(units, msg)
706 : : write(msg,"(a,i0,a,f5.1,a)") &
707 : 5 : " Number of BZ k-points close to the Fermi surface: ",fs%nkfs," [", (100.0_dp * fs%nkfs) / fs%nktot, " %]"
708 : 5 : call wrtout(units, msg)
709 : 5 : write(msg,"(a,i0)")" Maximum number of bands crossing the Fermi level: ",fs%maxnb
710 : 5 : call wrtout(units, msg)
711 [ + + + + ]: 150 : write(msg,"(2(a,i0))")" min band: ", minval(fs%bstart_cnt_ibz(1,:), mask=fs%bstart_cnt_ibz(1,:) /= -1)
712 : 5 : call wrtout(units, msg)
713 : 5 : write(msg,"(2(a,i0))")" Max band: ", maxval(fs%bstart_cnt_ibz(1,:) + fs%bstart_cnt_ibz(2,:) - 1, &
714 [ + + + + ]: 155 : mask=fs%bstart_cnt_ibz(1,:) /= -1)
715 : 10 : call wrtout(units, msg)
716 : : end associate
717 : : end do
718 : :
719 [ - + ]: 5 : end subroutine fstab_print
720 : : !!***
721 : :
722 [ # # # # : 0 : end module m_fstab
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
723 : : !!***
|
