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