Line data Source code
1 : !!****m* ABINIT/m_ephwg
2 : !! NAME
3 : !! m_ephwg
4 : !!
5 : !! FUNCTION
6 : !! Tools and objects to compute the weights used for the BZ integration of EPH quantities.
7 : !! More specifically the integration of quantities such as the imaginary part of the self-energy
8 : !! involving delta functions. Different approaches are available.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2008-2026 ABINIT group (MG, HM)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_ephwg
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_errors
29 : use m_xmpi
30 : use m_copy
31 : use m_dtset
32 : use m_htetra
33 : use m_nctk
34 : use netcdf
35 : use m_crystal
36 : use m_ifc
37 : use m_lgroup
38 : use m_ebands
39 : use m_eph_double_grid
40 : use m_krank
41 :
42 : use m_time, only : cwtime, cwtime_report
43 : use m_matrix, only : matr3inv
44 : use m_numeric_tools, only : arth, inrange, wrap2_pmhalf
45 : use m_special_funcs, only : gaussian
46 : use m_fstrings, only : strcat, ltoa, itoa, ftoa, ktoa, sjoin
47 : use m_simtet, only : sim0onei, SIM0TWOI
48 : use m_kpts, only : kpts_timrev_from_kptopt, kpts_ibz_from_kptrlatt, kpts_map
49 : use m_occ, only : occ_fd, occ_be
50 :
51 : implicit none
52 :
53 : private
54 : !!***
55 :
56 : !----------------------------------------------------------------------
57 :
58 : !!****t* m_ephwg/ephwg_t
59 : !! NAME
60 : !! ephwg_t
61 : !!
62 : !! FUNCTION
63 : !! Stores electron eigevalues and phonon frequencies in the IBZ (assume same mesh for e and ph).
64 : !! Provides tools to compute (e_{k+q} - w{q}) in the IBZ(k)
65 : !! and integrate the delta functions for phonon emission/absorption with the tetrahedron method.
66 : !!
67 : !! SOURCE
68 :
69 : type, public :: ephwg_t
70 :
71 : integer :: natom3
72 : ! 3 * natom
73 :
74 : integer :: nsppol
75 : ! Number of independent spin polarizations.
76 :
77 : integer :: nbcount
78 : ! Number of bands treated.
79 :
80 : integer :: bstart
81 : ! The fist band (global index) starts at bstart.
82 : ! Used to select bands around the Fermi level.
83 :
84 : integer :: kptopt
85 : ! Option for k-point generation.
86 :
87 : integer :: timrev
88 : ! 1 if the use of time-reversal is allowed; 0 otherwise
89 :
90 : integer :: nibz, nbz
91 : ! Number of q-points in IBZ and full BZ.
92 :
93 : integer :: nq_k
94 : ! Number of points in IBZ(k) i.e. the irreducible wedge
95 : ! defined by the operations of the little group of k.
96 :
97 : !real(dp) :: max_phfrq
98 : ! Max Phonon frequency, computed from phfrq_ibz
99 :
100 : integer :: kptrlatt(3,3)
101 : ! Value of kptrlatt after inkpts. So one shift
102 :
103 : integer, allocatable :: kq2ibz(:)
104 : ! kq2ibz(nq_k)
105 : ! Mapping (k + q) --> initial IBZ array
106 :
107 : real(dp),allocatable :: ibz(:,:)
108 : ! ibz(3, nibz)
109 : ! The initial IBZ.
110 :
111 : real(dp),allocatable :: bz(:,:)
112 : ! bz(3, nbz)
113 : ! points in full BZ.
114 :
115 : real(dp) :: klatt(3, 3)
116 : ! Reciprocal of lattice vectors for full kpoint grid. Used by init_tetra
117 :
118 : integer,allocatable :: lgk2ibz(:)
119 : ! lgk2ibz(nq_k)
120 : ! Mapping Little-group IBZ_k --> initial IBZ
121 : ! TODO: This should be generalized to have the symmetry indices as well so
122 : ! that we can use it in sigmaph but then we have to implement similar algo for double grid.
123 :
124 : real(dp),allocatable :: phfrq_ibz(:,:)
125 : ! (nibz, natom3)
126 : ! Phonon frequencies in the IBZ
127 :
128 : real(dp),allocatable :: eigkbs_ibz(:, :, :)
129 : ! (nibz, nbcount, nsppol)
130 : ! Electron eigenvalues in the IBZ for nbcount states
131 : ! (not necessarily equal to global nband, see also bstart and bcount)
132 :
133 : type(crystal_t), pointer :: cryst => null()
134 : ! Pointer to input structure (does not own memory)
135 :
136 : type(lgroup_t) :: lgk
137 : ! Little group of the k-point
138 :
139 : type(htetra_t) :: tetra_k
140 : ! Used to evaluate delta(w - e_{k+q} +/- phw_q) with tetrahedron method.
141 :
142 : contains
143 :
144 : procedure :: setup_kpoint => ephwg_setup_kpoint
145 : ! Prepare tetrahedron method for given external k-point.
146 :
147 : procedure :: double_grid_setup_kpoint => ephwg_double_grid_setup_kpoint
148 : ! Prepare tetrahedron method for given external k-point using double grid routines.
149 :
150 : procedure :: report_stats => ephwg_report_stats
151 : ! Report how much memory is being used by this object
152 :
153 : procedure :: get_deltas => ephwg_get_deltas
154 : ! Compute weights for $ \int \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
155 :
156 : procedure :: get_deltas_wvals => ephwg_get_deltas_wvals
157 : ! Compute weights for $ \int \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
158 :
159 : procedure :: get_deltas_qibzk => ephwg_get_deltas_qibzk
160 : ! Compute weights for $ \delta(\omega - \omega_{q\nu}} $ using IBZ(k) ordering
161 :
162 : procedure :: get_zinv_weights => ephwg_get_zinv_weights
163 : ! Compute weights for $ \int 1 / (\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
164 :
165 : procedure :: free => ephwg_free
166 : ! Free memory
167 : end type ephwg_t
168 :
169 : public :: ephwg_new ! Basic Constructor
170 : public :: ephwg_from_ebands ! Build object from ebands_t
171 :
172 : contains
173 : !!***
174 :
175 : !----------------------------------------------------------------------
176 :
177 : !!****f* m_ephwg/ephwg_new
178 : !! NAME
179 : !! ephwg_new
180 : !!
181 : !! FUNCTION
182 : !! Initialize the object from the electronic eigenvalues given in the IBZ.
183 : !!
184 : !! INPUTS
185 : !! cryst<cryst_t>=Crystalline structure.
186 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
187 : !! bstart=Index of the first band to be included.
188 : !! nbcount=Number of bands included
189 : !! kptopt=Option for the k-point generation.
190 : !! kptrlatt(3,3)=k-point lattice specification
191 : !! nshiftk= number of shift vectors.
192 : !! shiftk(3,nshiftk)=shift vectors for k point generation
193 : !! nkibz=Number of points in the IBZ
194 : !! kibz(3,nkibz)=Reduced coordinates of the k-points in the IBZ.
195 : !! nsppol=Number of independent spin polarizations.
196 : !! eig_ibz(nbcount, nkibz, nsppol) = Electron eigenvalues for nbcount states
197 : !!
198 : !! OUTPUT
199 : !!
200 : !! SOURCE
201 :
202 23 : type(ephwg_t) function ephwg_new( &
203 23 : cryst, ifc, bstart, nbcount, kptopt, kptrlatt, nshiftk, shiftk, nkibz, kibz, nsppol, eig_ibz, comm) result(new)
204 :
205 : !Arguments ------------------------------------
206 : !scalars
207 : integer,intent(in) :: kptopt, nshiftk, nkibz, bstart, nbcount, nsppol, comm
208 : type(crystal_t),target,intent(in) :: cryst
209 : type(ifc_type),intent(in) :: ifc
210 : !arrays
211 : integer,intent(in) :: kptrlatt(3,3)
212 : real(dp),intent(in) :: shiftk(3, nshiftk), kibz(3, nkibz)
213 : real(dp),intent(in) :: eig_ibz(nbcount, nkibz, nsppol)
214 :
215 : !Local variables-------------------------------
216 : !scalars
217 : integer :: nprocs, my_rank, ik, ierr, out_nkibz
218 : real(dp) :: cpu, wall, gflops
219 : !arrays
220 : real(dp) :: rlatt(3,3)
221 : integer :: out_kptrlatt(3,3)
222 46 : real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom), phfrq(3*cryst%natom)
223 23 : real(dp),allocatable :: out_kibz(:,:), out_wtk(:)
224 : !----------------------------------------------------------------------
225 :
226 23 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
227 :
228 23 : new%natom3 = ifc%natom * 3
229 23 : new%nsppol = nsppol
230 23 : new%nbcount = nbcount
231 23 : new%bstart = bstart
232 23 : new%kptopt = kptopt
233 23 : new%timrev = kpts_timrev_from_kptopt(new%kptopt)
234 23 : new%nibz = nkibz
235 23 : new%cryst => cryst
236 23 : call alloc_copy(kibz, new%ibz)
237 :
238 23 : call cwtime(cpu, wall, gflops, "start")
239 :
240 : ! Get full BZ (new%nbz, new%bz) and new kptrlatt for tetra.
241 : call kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, out_nkibz, out_kibz, out_wtk, new%nbz, new%bz, &
242 23 : new_kptrlatt=out_kptrlatt)
243 23 : call cwtime_report(" ephwg_new: kpts_ibz_from_kptrlatt", cpu, wall, gflops)
244 :
245 299 : new%kptrlatt = out_kptrlatt
246 299 : rlatt = out_kptrlatt; call matr3inv(rlatt, new%klatt)
247 :
248 23 : ABI_CHECK(size(out_kibz, dim=2) == new%nibz, "mismatch in nkibz!")
249 23 : ABI_FREE(out_kibz)
250 23 : ABI_FREE(out_wtk)
251 :
252 : ! Copy eigenvalues in IBZ. Change shape for better performance in other routines.
253 115 : ABI_MALLOC(new%eigkbs_ibz, (new%nibz, new%nbcount, new%nsppol))
254 6087 : do ik=1,new%nibz
255 20350 : new%eigkbs_ibz(ik, :, :) = eig_ibz(:, ik, :)
256 : end do
257 :
258 : ! Fourier interpolate phonon frequencies on the same mesh.
259 36176 : ABI_CALLOC(new%phfrq_ibz, (new%nibz, new%natom3))
260 :
261 6087 : do ik=1,new%nibz
262 6064 : if (mod(ik, nprocs) /= my_rank) cycle ! mpi-parallelism
263 6064 : call ifc%fourq(cryst, new%ibz(:, ik), phfrq, displ_cart)
264 42036 : new%phfrq_ibz(ik, :) = phfrq
265 : end do
266 :
267 : ! Collect results on each rank
268 23 : call xmpi_sum(new%phfrq_ibz, comm, ierr)
269 :
270 : !new%max_phfrq = maxval(phfrq_ibz)
271 :
272 23 : call cwtime_report(" ephwg_new: ifc_fourq", cpu, wall, gflops)
273 :
274 92 : end function ephwg_new
275 : !!***
276 :
277 : !----------------------------------------------------------------------
278 :
279 : !!****f* m_ephwg/ephwg_from_ebands
280 : !! NAME
281 : !! ephwg_from_ebands
282 : !!
283 : !! FUNCTION
284 : !! Convenience constructor to initialize the object from an ebands_t object
285 :
286 23 : type(ephwg_t) function ephwg_from_ebands(cryst, ifc, ebands, bstart, nbcount, comm) result(new)
287 :
288 : !Arguments ------------------------------------
289 : !scalars
290 : integer,intent(in) :: bstart, nbcount, comm
291 : type(crystal_t),intent(in) :: cryst
292 : type(ifc_type),intent(in) :: ifc
293 : type(ebands_t),intent(in) :: ebands
294 :
295 : !Local variables-------------------------------
296 23 : real(dp),allocatable :: eig_ibz(:, :, :)
297 : !----------------------------------------------------------------------
298 :
299 12 : if (bstart == 1 .and. nbcount == ebands%mband) then
300 : new = ephwg_new(cryst, ifc, bstart, nbcount, ebands%kptopt, ebands%kptrlatt, ebands%nshiftk, ebands%shiftk, ebands%nkpt, &
301 1 : ebands%kptns, ebands%nsppol, ebands%eig, comm)
302 : else
303 66 : ABI_CHECK(inrange(bstart, [1, ebands%mband]), "Wrong bstart")
304 66 : ABI_CHECK(inrange(bstart + nbcount - 1, [1, ebands%mband]), "Wrong nbcount")
305 : ! Copy submatrix of eigenvalues
306 110 : ABI_MALLOC(eig_ibz, (nbcount, ebands%nkpt, ebands%nsppol))
307 13889 : eig_ibz = ebands%eig(bstart:bstart+nbcount-1, : , :)
308 : new = ephwg_new(cryst, ifc, bstart, nbcount, ebands%kptopt, ebands%kptrlatt, ebands%nshiftk, ebands%shiftk, ebands%nkpt, &
309 22 : ebands%kptns, ebands%nsppol, eig_ibz, comm)
310 22 : ABI_FREE(eig_ibz)
311 : end if
312 :
313 23 : end function ephwg_from_ebands
314 : !!***
315 :
316 : !!****f* m_ephwg/ephwg_setup_kpoint
317 : !! NAME
318 : !! ephwg_setup_kpoint
319 : !!
320 : !! FUNCTION
321 : !! Set internal tables and object required to compute integration weights for a given k-point.
322 : !!
323 : !! INPUTS
324 : !! kpoint(3): k-point in reduced coordinates.
325 : !! prtvol: Verbosity level
326 : !! comm: MPI communicator
327 : !!
328 : !! OUTPUT
329 : !!
330 : !! SOURCE
331 :
332 62 : subroutine ephwg_setup_kpoint(self, kpoint, prtvol, comm, skip_mapping)
333 :
334 : !Arguments ------------------------------------
335 : !scalars
336 : class(ephwg_t),target,intent(inout) :: self
337 : integer,intent(in) :: prtvol, comm
338 : logical,optional,intent(in) :: skip_mapping
339 : !arrays
340 : real(dp),intent(in) :: kpoint(3)
341 :
342 : !Local variables-------------------------------
343 : !scalars
344 : integer :: ierr,ii
345 : logical :: do_mapping
346 : real(dp) :: cpu, wall, gflops
347 : character(len=80) :: errorstring
348 : !character(len=500) :: msg
349 : type(crystal_t),pointer :: cryst
350 62 : type(krank_t) :: krank
351 : !arrays
352 62 : integer,allocatable :: indkk(:,:)
353 : !----------------------------------------------------------------------
354 :
355 62 : do_mapping = .true.; if (present(skip_mapping)) do_mapping = .not. skip_mapping
356 62 : cryst => self%cryst
357 62 : call cwtime(cpu, wall, gflops, "start")
358 :
359 : ! Get little group of the (external) kpoint.
360 62 : call self%lgk%free()
361 62 : call self%lgk%init(self%cryst, kpoint, self%timrev, self%nbz, self%bz, self%nibz, self%ibz, comm)
362 :
363 68 : if (prtvol > 0) call self%lgk%print([std_out])
364 62 : self%nq_k = self%lgk%nibz
365 :
366 62 : call cwtime_report(" lgroup_new", cpu, wall, gflops)
367 :
368 62 : if (do_mapping) then
369 : ! TODO: Use symrec conventions although this means that we cannot reuse these tables
370 : ! to symmetrize wavefunctions and potentials that require S-1 i.e. the symrel convention.
371 :
372 : ! Get mapping IBZ_k --> initial IBZ (self%lgk%ibz --> self%ibz)
373 0 : ABI_MALLOC(indkk, (6, self%nq_k))
374 :
375 0 : call krank%from_kptrlatt(self%nibz, self%ibz, self%kptrlatt, compute_invrank=.False.)
376 :
377 0 : if (kpts_map("symrel", self%kptopt, cryst, krank, self%nq_k, self%lgk%ibz, indkk) /= 0) then
378 0 : ABI_ERROR("At least one of the points in IBZ(k) could not be generated from a symmetrical one.")
379 : end if
380 :
381 0 : call krank%free()
382 :
383 0 : ABI_SFREE(self%lgk2ibz)
384 0 : call alloc_copy(indkk(1, :), self%lgk2ibz)
385 0 : ABI_FREE(indkk)
386 0 : call cwtime_report(" listkk1", cpu, wall, gflops)
387 :
388 : ! Get mapping (k + q) --> initial IBZ.
389 0 : do ii=1,self%nq_k
390 0 : self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) + kpoint
391 : end do
392 0 : ABI_MALLOC(indkk, (6, self%nq_k))
393 :
394 0 : call krank%from_kptrlatt(self%nibz, self%ibz, self%kptrlatt, compute_invrank=.False.)
395 :
396 0 : if (kpts_map("symrel", self%kptopt, cryst, krank, self%nq_k, self%lgk%ibz, indkk) /= 0) then
397 0 : ABI_ERROR("At least one of the points in IBZ(k) + q could not be generated from a symmetrical one.")
398 : end if
399 0 : call krank%free()
400 :
401 0 : call cwtime_report(" listkk2", cpu, wall, gflops)
402 :
403 0 : ABI_SFREE(self%kq2ibz)
404 0 : call alloc_copy(indkk(1, :), self%kq2ibz)
405 0 : ABI_FREE(indkk)
406 :
407 : ! Revert changes
408 0 : do ii=1,self%nq_k
409 0 : self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) - kpoint
410 : end do
411 : end if
412 :
413 : ! Get mapping BZ --> IBZ_k (self%bz --> self%lgrp%ibz) required for tetrahedron method
414 186 : ABI_MALLOC(indkk, (self%nbz, 1))
415 429950 : indkk(:, 1) = self%lgk%bz2ibz_smap(1, :)
416 :
417 : ! Build tetrahedron object using IBZ(k) as the effective IBZ
418 : ! This means that input data for tetra routines must be provided in lgk%kibz_q
419 62 : call self%tetra_k%free()
420 : call self%tetra_k%init(indkk(:, 1), cryst%gprimd, self%klatt, self%bz, self%nbz, &
421 62 : self%lgk%ibz, self%nq_k, ierr, errorstring, comm)
422 : !call tetra_write(self%tetra_k, self%lgk%nibz, self%lgk%ibz, strcat("tetrak_", ktoa(kpoint)))
423 62 : ABI_CHECK(ierr == 0, errorstring)
424 :
425 62 : if (xmpi_comm_rank(comm) == 0) call self%tetra_k%print(std_out)
426 62 : ABI_FREE(indkk)
427 :
428 62 : call cwtime_report(" init_tetra", cpu, wall, gflops)
429 :
430 186 : end subroutine ephwg_setup_kpoint
431 : !!***
432 :
433 : !!****f* m_ephwg/ephwg_double_grid_setup_kpoint
434 : !! NAME
435 : !! ephwg_setup_kpoint
436 : !!
437 : !! FUNCTION
438 : !! Set internal tables and object required to compute integration weights for a given k-point
439 : !! using the double grid routines to map the different k-points.
440 : !! This version should be more efficient than its counterpart without the double grid.
441 : !!
442 : !! INPUTS
443 : !! kpoint(3): k-point in reduced coordinates.
444 : !! prtvol: Verbosity level
445 : !! comm: MPI communicator
446 : !!
447 : !! OUTPUT
448 : !!
449 : !! SOURCE
450 :
451 19 : subroutine ephwg_double_grid_setup_kpoint(self, eph_doublegrid, kpoint, prtvol, comm)
452 :
453 : !Arguments ------------------------------------
454 : !scalars
455 : class(ephwg_t),target,intent(inout) :: self
456 : type(eph_double_grid_t),intent(inout) :: eph_doublegrid
457 : integer,intent(in) :: prtvol, comm
458 : !arrays
459 : real(dp),intent(in) :: kpoint(3)
460 :
461 : !Local variables-------------------------------
462 : !scalars
463 : integer,parameter :: timrev0 = 0
464 : integer :: ierr,ii,ik_idx
465 : character(len=80) :: errorstring
466 : !character(len=500) :: msg
467 : type(crystal_t),pointer :: cryst
468 : !arrays
469 19 : integer,allocatable :: lgkibz2bz(:) !indkk(:,:),
470 19 : integer,allocatable :: bz2lgkibz(:), bz2lgkibzkq(:) !, bz2bz(:), mapping(:,:)
471 : !real(dp) :: kpt(3), wrap_kpt(3), shift
472 : !----------------------------------------------------------------------
473 :
474 19 : cryst => self%cryst
475 :
476 : ! Get little group of the (external) kpoint.
477 19 : call self%lgk%free()
478 19 : call self%lgk%init(self%cryst, kpoint, self%timrev, self%nbz, self%bz, self%nibz, self%ibz, comm)
479 19 : if (prtvol > 0) call self%lgk%print([std_out])
480 19 : self%nq_k = self%lgk%nibz
481 :
482 : ! get dg%bz --> self%lgrp%ibz
483 57 : ABI_REMALLOC(eph_doublegrid%bz2lgkibz, (eph_doublegrid%dense_nbz))
484 :
485 : ! Old version using all crystal symmetries
486 : !call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
487 : ! cryst%symrel, cryst%nsym, &
488 : ! eph_doublegrid%bz2lgkibz, has_timrev=1)
489 : call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
490 : self%lgk%symrec_lg, self%lgk%nsym_lg, &
491 19 : eph_doublegrid%bz2lgkibz, timrev0, use_symrec=.true.)
492 :
493 : ! self%lgrp%ibz --> dg%bz
494 80998 : ABI_ICALLOC(lgkibz2bz, (self%lgk%nibz))
495 339987 : do ii=1,self%nbz
496 339968 : ik_idx = eph_doublegrid%bz2lgkibz(ii)
497 339987 : lgkibz2bz(ik_idx) = ii
498 : enddo
499 :
500 : ! get self%lgrp%ibz --> dg%bz --> self%ibz
501 57 : ABI_REMALLOC(self%lgk2ibz, (self%nq_k))
502 80960 : do ii=1,self%nq_k
503 80941 : ik_idx = lgkibz2bz(ii)
504 80960 : self%lgk2ibz(ii) = eph_doublegrid%bz2ibz_dense(ik_idx)
505 : enddo
506 :
507 : ! calculate k+q
508 80960 : do ii=1,self%nq_k
509 323783 : self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) + kpoint
510 : end do
511 :
512 : ! get dg%bz --> lgk%ibz (k+q)
513 57 : ABI_MALLOC(bz2lgkibzkq, (eph_doublegrid%dense_nbz))
514 : ! Old version using all crystal symmetries
515 : !call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
516 : ! cryst%symrel, cryst%nsym, &
517 : ! bz2lgkibzkq, has_timrev=1)
518 : call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
519 : self%lgk%symrec_lg, self%lgk%nsym_lg, &
520 19 : bz2lgkibzkq, timrev0, use_symrec=.true.)
521 :
522 : ! self%lgrp%ibz (k+q) --> dg%bz
523 339987 : do ii=1,self%nbz
524 339968 : ik_idx = bz2lgkibzkq(ii)
525 339987 : lgkibz2bz(ik_idx) = ii
526 : enddo
527 19 : ABI_FREE(bz2lgkibzkq)
528 :
529 : ! get self%lgrp%ibz (k+q) --> dg%bz --> self%ibz
530 57 : ABI_REMALLOC(self%kq2ibz, (self%nq_k))
531 80960 : do ii=1,self%nq_k
532 80941 : ik_idx = lgkibz2bz(ii)
533 80960 : self%kq2ibz(ii) = eph_doublegrid%bz2ibz_dense(ik_idx)
534 : enddo
535 19 : ABI_FREE(lgkibz2bz)
536 :
537 : ! revert change
538 80960 : do ii=1,self%nq_k
539 323783 : self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) - kpoint
540 : end do
541 :
542 : ! get self%bz --> dg%bz --> self%lgrp%ibz
543 57 : ABI_MALLOC(bz2lgkibz, (self%nbz))
544 :
545 339987 : do ii=1,self%nbz
546 : ! get self%bz --> dg%bz
547 339968 : ik_idx = eph_doublegrid%get_index(self%bz(:,ii),2)
548 : ! dg%bz --> self%lgrp%ibz
549 339987 : bz2lgkibz(ii) = eph_doublegrid%bz2lgkibz(ik_idx)
550 : end do
551 :
552 : ! Build tetrahedron object using IBZ(k) as the effective IBZ
553 : ! This means that input data for tetra routines must be provided in lgk%kibz_q
554 19 : call self%tetra_k%free()
555 : call self%tetra_k%init(bz2lgkibz, cryst%gprimd, self%klatt, self%bz, self%nbz, &
556 19 : self%lgk%ibz, self%nq_k, ierr, errorstring, comm)
557 19 : if (ierr /= 0) then
558 0 : ABI_ERROR(errorstring)
559 : end if
560 19 : ABI_FREE(bz2lgkibz)
561 :
562 19 : end subroutine ephwg_double_grid_setup_kpoint
563 : !!***
564 :
565 : !----------------------------------------------------------------------
566 :
567 : !!****f* m_ephwg/ephwg_report_stats
568 : !! NAME
569 : !! ephwg_report_stats
570 : !!
571 : !! FUNCTION
572 : !!
573 : !! INPUTS
574 : !!
575 : !! OUTPUT
576 : !!
577 : !! SOURCE
578 :
579 81 : subroutine ephwg_report_stats(self)
580 :
581 : !Arguments ------------------------------------
582 : !scalars
583 : class(ephwg_t),intent(in) :: self
584 :
585 : !Variables
586 : real(dp) :: mem_tot
587 : !----------------------------------------------------------------------
588 :
589 : ! IBZ qpoints
590 81 : mem_tot = 3 * self%nibz * dp
591 : ! BZ qpoints
592 81 : mem_tot = mem_tot + 3 * self%nbz * dp
593 : ! lgk2ibz and kq2ibz
594 81 : mem_tot = mem_tot + self%nq_k * 2 * 4
595 : ! phonon frequencies
596 81 : mem_tot = mem_tot + self%nibz * self%natom3 * dp
597 : ! eigenvalues
598 81 : mem_tot = mem_tot + self%nibz * self%nbcount * self%nsppol * dp
599 :
600 81 : write(std_out,"(a,f8.1,a)") " Memory allocated for ephwg weights:", mem_tot * b2Mb, " [Mb] <<< MEM"
601 :
602 81 : end subroutine ephwg_report_stats
603 : !!***
604 :
605 : !----------------------------------------------------------------------
606 :
607 : !!****f* m_ephwg/ephwg_get_deltas
608 : !! NAME
609 : !! ephwg_get_deltas
610 : !!
611 : !! FUNCTION
612 : !! Compute weights for $ \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
613 : !! for a given (band, spin) and phonon mode nu.
614 : !!
615 : !! INPUTS
616 : !! band=band index (global index i.e. unshifted)
617 : !! spin=Spin index
618 : !! nu=Phonon branch.
619 : !! nene=number of energies for DOS
620 : !! eminmax=min and energy in delta (linear mesh)
621 : !! bcorr=1 to include Blochl correction else 0.
622 : !! comm=MPI communicator
623 : !! [broad]=Gaussian broadening
624 : !!
625 : !! OUTPUT
626 : !! deltaw_pm(nene, nq_k, 2) (plus, minus) including the weights for BZ integration.
627 : ! These arrays have the same order as the little group used in sigmaph.
628 : !!
629 : !! SOURCE
630 :
631 0 : subroutine ephwg_get_deltas(self, band, spin, nu, nene, eminmax, bcorr, deltaw_pm, comm, &
632 : broad) ! optional
633 :
634 : !Arguments ------------------------------------
635 : !scalars
636 : integer,intent(in) :: band, spin, nu, nene, bcorr, comm
637 : class(ephwg_t),intent(in) :: self
638 : real(dp),optional,intent(in) :: broad
639 : !arrays
640 : real(dp),intent(in) :: eminmax(2)
641 : real(dp),intent(out) :: deltaw_pm(nene, self%nq_k, 2)
642 :
643 : !Local variables-------------------------------
644 : !scalars
645 : integer :: iq,iq_ibz,ikpq_ibz,ib,ie
646 : real(dp),parameter :: max_occ1 = one
647 : real(dp) :: omega_step
648 : !arrays
649 0 : real(dp) :: wme0(nene)
650 0 : real(dp),allocatable :: thetaw(:,:), pme_k(:,:)
651 : !----------------------------------------------------------------------
652 :
653 0 : ib = band - self%bstart + 1
654 :
655 0 : ABI_MALLOC(thetaw, (nene, self%nq_k))
656 0 : ABI_MALLOC(pme_k, (self%nq_k, 2))
657 :
658 : ! Fill array for e_{k+q, b} +- w_{q,nu)
659 0 : do iq=1,self%nq_k
660 0 : iq_ibz = self%lgk2ibz(iq) ! IBZ_k --> IBZ
661 0 : ikpq_ibz = self%kq2ibz(iq) ! k + q --> IBZ
662 0 : pme_k(iq, 1) = self%eigkbs_ibz(ikpq_ibz, ib, spin) - self%phfrq_ibz(iq_ibz, nu)
663 0 : pme_k(iq, 2) = self%eigkbs_ibz(ikpq_ibz, ib, spin) + self%phfrq_ibz(iq_ibz, nu)
664 : end do
665 :
666 0 : if (present(broad)) then
667 0 : omega_step = (eminmax(2) - eminmax(1)) / (nene - 1)
668 : ! Use thetaw as workspace array
669 0 : thetaw(:, 1) = arth(eminmax(1), omega_step, nene)
670 0 : do iq=1,self%nq_k
671 0 : do ie=1,2
672 0 : wme0 = thetaw(:, 1) - pme_k(iq, ie)
673 0 : deltaw_pm(:, iq, ie) = gaussian(wme0, broad)
674 : end do
675 : end do
676 :
677 : ! Multiply by weights
678 0 : do ie=1,nene
679 0 : deltaw_pm(ie, :, 1) = deltaw_pm(ie, :, 1) * self%lgk%weights
680 0 : deltaw_pm(ie, :, 2) = deltaw_pm(ie, :, 2) * self%lgk%weights
681 : end do
682 :
683 : else
684 : ! TODO Add routine to compute only delta
685 : call self%tetra_k%blochl_weights(pme_k(:,1), eminmax(1), eminmax(2), max_occ1, nene, self%nq_k, &
686 0 : bcorr, thetaw, deltaw_pm(:,:,1), comm)
687 : call self%tetra_k%blochl_weights(pme_k(:,2), eminmax(1), eminmax(2), max_occ1, nene, self%nq_k, &
688 0 : bcorr, thetaw, deltaw_pm(:,:,2), comm)
689 : end if
690 :
691 0 : ABI_FREE(thetaw)
692 0 : ABI_FREE(pme_k)
693 :
694 0 : end subroutine ephwg_get_deltas
695 : !!***
696 :
697 : !----------------------------------------------------------------------
698 :
699 : !!****f* m_ephwg/ephwg_get_deltas_wvals
700 : !! NAME
701 : !! ephwg_get_deltas_wvals
702 : !!
703 : !! FUNCTION
704 : !! Compute weights for $ \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
705 : !! for a given (band, spin) and phonon mode nu.
706 : !!
707 : !! INPUTS
708 : !! band=band index (global index i.e. unshifted)
709 : !! spin=Spin index
710 : !! nu=Phonon branch.
711 : !! nene=number of energies for DOS
712 : !! eminmax=min and energy in delta (linear mesh)
713 : !! bcorr=1 to include Blochl correction else 0.
714 : !! comm=MPI communicator
715 : !! [broad]=Gaussian broadening
716 : !!
717 : !! OUTPUT
718 : !! deltaw_pm(nene, nq_k, 2) (plus, minus) including the weights for BZ integration.
719 : ! These arrays have the same order as the little group used in sigmaph.
720 : !!
721 : !! SOURCE
722 :
723 15108 : subroutine ephwg_get_deltas_wvals(self, band, spin, nu, neig, eig, bcorr, deltaw_pm, comm, &
724 : broad) ! optional
725 :
726 : !Arguments ------------------------------------
727 : !scalars
728 : integer,intent(in) :: band, spin, nu, neig, bcorr, comm
729 : real(dp),intent(in) :: eig(neig)
730 : class(ephwg_t),intent(in) :: self
731 : real(dp),optional,intent(in) :: broad
732 : !arrays
733 : real(dp),intent(out) :: deltaw_pm(neig,self%nq_k, 2)
734 :
735 : !Local variables-------------------------------
736 : !scalars
737 : real(dp),parameter :: max_occ1 = one
738 : integer :: iq, iq_ibz, ikpq_ibz, ib, nprocs, my_rank
739 30216 : real(dp) :: wme0(neig)
740 : !arrays
741 15108 : real(dp),allocatable :: pme_k(:,:)
742 : !----------------------------------------------------------------------
743 :
744 15108 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
745 15108 : ib = band - self%bstart + 1
746 13081836 : deltaw_pm = zero
747 :
748 45324 : ABI_MALLOC(pme_k, (self%nq_k, 2))
749 :
750 : ! Fill array for e_{k+q, b} +- w_{q,nu)
751 3274236 : do iq=1,self%nq_k
752 3259128 : iq_ibz = self%lgk2ibz(iq) ! IBZ_k --> IBZ
753 3259128 : ikpq_ibz = self%kq2ibz(iq) ! k + q --> IBZ
754 3259128 : pme_k(iq, 1) = self%eigkbs_ibz(ikpq_ibz, ib, spin) - self%phfrq_ibz(iq_ibz, nu)
755 3274236 : pme_k(iq, 2) = self%eigkbs_ibz(ikpq_ibz, ib, spin) + self%phfrq_ibz(iq_ibz, nu)
756 : end do
757 :
758 : ! Compute the tetrahedron or gaussian weights
759 15108 : if (present(broad)) then
760 0 : do iq_ibz=1,self%nq_k
761 0 : if (mod(iq_ibz, nprocs) /= my_rank) cycle ! MPI parallelism
762 0 : wme0 = eig - pme_k(iq_ibz, 1)
763 0 : deltaw_pm(:,iq_ibz,1) = gaussian(wme0, broad) * self%lgk%weights(iq_ibz)
764 0 : wme0 = eig - pme_k(iq_ibz, 2)
765 0 : deltaw_pm(:,iq_ibz,2) = gaussian(wme0, broad) * self%lgk%weights(iq_ibz)
766 : end do
767 : else
768 15108 : call self%tetra_k%wvals_weights_delta(pme_k(:, 1), neig, eig, max_occ1, self%nq_k, bcorr, deltaw_pm(:,:,1), comm)
769 15108 : call self%tetra_k%wvals_weights_delta(pme_k(:, 2), neig, eig, max_occ1, self%nq_k, bcorr, deltaw_pm(:,:,2), comm)
770 : end if
771 :
772 15108 : ABI_FREE(pme_k)
773 :
774 15108 : end subroutine ephwg_get_deltas_wvals
775 : !!***
776 :
777 : !----------------------------------------------------------------------
778 :
779 : !!****f* m_ephwg/ephwg_get_deltas_qibzk
780 : !! NAME
781 : !! ephwg_get_deltas_qibzk
782 : !!
783 : !! FUNCTION
784 : !! Compute weights for $ \delta(\omega - \omega_{q\nu} $ for given nu, using q-opints in the IBZ(k)
785 : !!
786 : !! INPUTS
787 : !! nu=Phonon branch.
788 : !! nene=number of energies for DOS
789 : !! eminmax=min and energy in delta (linear mesh)
790 : !! bcorr=1 to include Blochl correction else 0.
791 : !! comm=MPI communicator
792 : !! [with_qweights]= .False. if q-point weights should not be included in dt_weights
793 : !!
794 : !! OUTPUT
795 : !! dt_weights(nene, nq_k, 2) weights for BZ integration (delta and theta function)
796 : ! These arrays have the same order as the q-points in the little group of the k-point.
797 : !!
798 : !! SOURCE
799 :
800 114 : subroutine ephwg_get_deltas_qibzk(self, nu, nene, eminmax, bcorr, dt_weights, comm, with_qweights)
801 :
802 : !Arguments ------------------------------------
803 : !scalars
804 : integer,intent(in) :: nu, nene, bcorr, comm
805 : class(ephwg_t),intent(in) :: self
806 : logical,optional,intent(in) :: with_qweights
807 : !arrays
808 : real(dp),intent(in) :: eminmax(2)
809 : real(dp),intent(out) :: dt_weights(nene, self%nq_k, 2)
810 :
811 : !Local variables-------------------------------
812 : !scalars
813 : integer :: iq, iq_ibz, ie, ii
814 : real(dp),parameter :: max_occ1 = one
815 : !arrays
816 114 : real(dp),allocatable :: eigen_in(:)
817 : !----------------------------------------------------------------------
818 :
819 342 : ABI_MALLOC(eigen_in, (self%nq_k))
820 :
821 : ! Fill eigen_in
822 5874 : do iq=1,self%nq_k
823 5760 : iq_ibz = self%lgk2ibz(iq) ! IBZ_k --> IBZ
824 5874 : eigen_in(iq) = self%phfrq_ibz(iq_ibz, nu)
825 : end do
826 :
827 : call self%tetra_k%blochl_weights(eigen_in, eminmax(1), eminmax(2), max_occ1, nene, self%nq_k, &
828 114 : bcorr, dt_weights(:,:,2), dt_weights(:,:,1), comm)
829 :
830 114 : if (present(with_qweights)) then
831 114 : if (.not. with_qweights) then
832 0 : do ii=1,2
833 0 : do ie=1,nene
834 0 : dt_weights(ie, :, ii) = dt_weights(ie, :, ii) * self%lgk%weights
835 : end do
836 : end do
837 : end if
838 : end if
839 :
840 114 : ABI_FREE(eigen_in)
841 :
842 114 : end subroutine ephwg_get_deltas_qibzk
843 : !!***
844 :
845 : !----------------------------------------------------------------------
846 :
847 : !!****f* m_ephwg/ephwg_get_zinv_weights
848 : !! NAME
849 : !! ephwg_get_zinv_weights
850 : !!
851 : !! FUNCTION
852 : !! Compute weights for a given (kpoint, qpoint, spin) for all phonon modes.
853 : !!
854 : !! INPUTS
855 : !! nz: Number of frequencies
856 : !! nbcalc=Number of bands in self-energy matrix elements.
857 : !! zvals(nw): z-values
858 : !! iband_sum = band index in self-energy sum. (global index i.e. unshifted)
859 : !! spin=Spin index
860 : !! nu=Phonon branch index
861 : !! zinv_opt:
862 : !! 1 for S. Kaprzyk routines,
863 : !! 2 for Lambin-Vigneron.
864 : !! comm=MPI communicator
865 : !! [use_bzsum]= By default the weights are multiplied by the Nstar(q) / Nq where
866 : !! Nstar(q) is the number of points in the star of the q-point (using the symmetries of the little group of k)
867 : !! If use_bzsum is set to True, the Nstar(q) coefficient is removed so that the caller can
868 : !! integrate over the BZ without using symmetries.
869 : !! [erange(2)]: if present, weights are computed with an approximated asyntotic expression if
870 : !! real(z) is outside of this interval and with tetra if inside.
871 :
872 : !! OUTPUT
873 : !! cweights(nz, 2, nbcalc, %nq_k) (plus, minus)
874 : !! include weights for BZ integration.
875 : !!
876 : !! SOURCE
877 :
878 1128 : subroutine ephwg_get_zinv_weights(self, nz, nbcalc, zvals, iband_sum, spin, nu, zinv_opt, cweights, comm, use_bzsum, erange)
879 :
880 : !Arguments ------------------------------------
881 : !scalars
882 : integer,intent(in) :: iband_sum, spin, nu, nz, nbcalc, zinv_opt, comm
883 : class(ephwg_t),intent(in) :: self
884 : logical, optional, intent(in) :: use_bzsum
885 : !arrays
886 : complex(dp),intent(in) :: zvals(nz, nbcalc)
887 : complex(dp),intent(out) :: cweights(nz, 2, nbcalc, self%nq_k)
888 : real(dp),optional,intent(in) :: erange(2)
889 :
890 : !Local variables-------------------------------
891 : !scalars
892 : integer,parameter :: master = 0
893 : integer :: iq_ibz, ikpq_ibz, ib, ii, iq, nprocs, my_rank !, ierr
894 : real(dp),parameter :: max_occ1 = one
895 : !real(dp) :: emin, emax
896 : logical :: use_bzsum_
897 : !arrays
898 : real(dp) :: my_erange(2)
899 1128 : real(dp),allocatable :: pme_k(:,:)
900 1128 : complex(dp),allocatable :: cweights_tmp(:,:)
901 : !----------------------------------------------------------------------
902 :
903 1128 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
904 :
905 1128 : use_bzsum_ = .False.; if (present(use_bzsum)) use_bzsum_ = use_bzsum
906 1128 : my_erange = [-huge(one), huge(one)]; if (present(erange)) my_erange = erange
907 :
908 : ! Allocate array for e_{k+q, b} +- w_{q,nu)
909 3384 : ABI_MALLOC(pme_k, (self%nq_k, 2))
910 :
911 1128 : ib = iband_sum - self%bstart + 1
912 12972 : do iq=1,self%nq_k
913 11844 : iq_ibz = self%lgk2ibz(iq) ! IBZ_k --> IBZ
914 11844 : ikpq_ibz = self%kq2ibz(iq) ! k + q --> IBZ
915 11844 : pme_k(iq, 1) = self%eigkbs_ibz(ikpq_ibz, ib, spin) - self%phfrq_ibz(iq_ibz, nu)
916 12972 : pme_k(iq, 2) = self%eigkbs_ibz(ikpq_ibz, ib, spin) + self%phfrq_ibz(iq_ibz, nu)
917 : end do
918 :
919 : ! As this part is quite demanding, especially when nz is large, use input z-mesh
920 : ! when we are inside the window in which the denominator can blow up (+- some tolerance)
921 : ! Outside the window, downsample the mesh use to compute the weights and spline the results.
922 : !if (zinv_opt == 2) then
923 : ! emin = minval(self%eigkbs_ibz(:, ib, spin))
924 : ! emax = maxval(self%eigkbs_ibz(:, ib, spin))
925 : ! !my_erange = [emin - half * abs(emin), emax + half * abs(emax)]
926 : ! my_erange = [emin - tol2 * abs(emin), emax + tol2 * abs(emax)]
927 : !end if
928 :
929 8631372 : cweights = zero
930 4512 : ABI_MALLOC(cweights_tmp, (nz, self%nq_k))
931 :
932 10152 : do ib=1,nbcalc
933 28200 : do ii=1,2
934 : call self%tetra_k%weights_wvals_zinv(pme_k(:, ii), nz, zvals(:, ib), max_occ1, self%nq_k, zinv_opt, &
935 18048 : cweights_tmp, comm, erange=my_erange)
936 216576 : do iq=1,self%nq_k
937 8541696 : cweights(:, ii, ib, iq) = cweights_tmp(:, iq)
938 : end do
939 : end do
940 : end do
941 :
942 1128 : ABI_FREE(cweights_tmp)
943 :
944 : ! Rescale weights so that the caller can sum over the full BZ.
945 : !if (use_bzsum_) cweights = cweights / ( self%lgk%weights(iqlk) * self%nbz )
946 :
947 1128 : ABI_FREE(pme_k)
948 : !call xmpi_sum(cweights, comm, ierr)
949 :
950 1128 : end subroutine ephwg_get_zinv_weights
951 : !!***
952 :
953 : !!****f* m_ephwg/ephwg_free
954 : !! NAME
955 : !! ephwg_free
956 : !!
957 : !! FUNCTION
958 : !! Deallocate memory
959 : !!
960 : !! INPUTS
961 : !!
962 : !! OUTPUT
963 : !!
964 : !! SOURCE
965 :
966 131 : subroutine ephwg_free(self)
967 :
968 : !Arguments ------------------------------------
969 : class(ephwg_t),intent(inout) :: self
970 : !----------------------------------------------------------------------
971 :
972 : ! integer
973 131 : ABI_SFREE(self%kq2ibz)
974 :
975 : ! Real
976 131 : ABI_SFREE(self%ibz)
977 131 : ABI_SFREE(self%bz)
978 131 : ABI_SFREE(self%lgk2ibz)
979 131 : ABI_SFREE(self%phfrq_ibz)
980 131 : ABI_SFREE(self%eigkbs_ibz)
981 :
982 : ! types
983 131 : call self%tetra_k%free()
984 131 : call self%lgk%free()
985 :
986 : ! nullify pointers
987 131 : self%cryst => null()
988 :
989 131 : end subroutine ephwg_free
990 : !!***
991 :
992 0 : end module m_ephwg
993 : !!***
|
