Line data Source code
1 : !!****m* ABINIT/m_sigtk
2 : !! NAME
3 : !! m_sigtk
4 : !!
5 : !! FUNCTION
6 : !! Helper functions common to electron self-energy calculations. Provides tools to:
7 : !! Define list of k-points and bands in sel-energy matrix elements from input variables.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_sigtk
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_ebands
29 : use m_crystal
30 : use m_xmpi
31 : use netcdf
32 : use m_nctk
33 : use m_hdr
34 : use m_dtset
35 : use m_krank
36 :
37 : use m_build_info, only : abinit_version
38 : use m_fstrings, only : sjoin, ltoa, strcat, itoa, ftoa
39 : use m_io_tools, only : open_file
40 : use defs_datatypes, only : pseudopotential_type
41 : use defs_wvltypes, only : wvl_internal_type
42 : use m_gwdefs, only : sigijtab_t, sigijtab_free
43 : use m_esymm, only : esymm_t
44 : use m_pawtab, only : pawtab_type
45 : use m_kpts, only : kpts_ibz_from_kptrlatt, kpts_timrev_from_kptopt, kpts_map
46 :
47 : implicit none
48 :
49 : private
50 :
51 : public :: sigtk_kcalc_from_nkptgw
52 : public :: sigtk_kcalc_from_qprange
53 : public :: sigtk_kcalc_from_gaps
54 : public :: sigtk_kcalc_from_erange
55 : public :: sigtk_kpts_in_erange
56 : public :: sigtk_sigma_tables
57 : public :: sigtk_multiply_by_vc_sqrt
58 : public :: sigtk_dw_tpp_red
59 : !!***
60 :
61 : ! Tables for degenerated KS states.
62 : type, public :: bids_t
63 : integer, allocatable :: vals(:)
64 : end type bids_t
65 :
66 : type, public :: degtab_t
67 : type(bids_t), allocatable :: bids(:)
68 : contains
69 : procedure :: free => degtab_free
70 : end type degtab_t
71 :
72 : public :: degtab_array_free ! Free array of degtab_t objects.
73 : !!***
74 :
75 : contains !=====================================================
76 : !!***
77 :
78 : !!****f* m_sigtk/sigtk_kcalc_from_nkptgw
79 : !! NAME
80 : !! sigtk_kcalc_from_nkptgw
81 : !!
82 : !! FUNCTION
83 : !! Initialize list of k-points and bands for self-energy matrix elements from nkptgw.
84 : !!
85 : !! INPUT
86 : !! dtset<dataset_type>=All input variables for this dataset.
87 : !! mband: Max number of bands.
88 : !!
89 : !! OUTPUT
90 : !! nkcalc: Number of k-points in self-energy matrix elements.
91 : !! kcalc(3, nkcalc): List of k-points where the self-energy is computed.
92 : !! bstart_ks(nkcalc, nsppol): Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
93 : !! nbcalc_ks(nkcalc, nsppol): Number of bands included in self-energy matrix elements for each k-point in kcalc.
94 : !!
95 : !! SOURCE
96 :
97 27 : subroutine sigtk_kcalc_from_nkptgw(dtset, mband, nkcalc, kcalc, bstart_ks, nbcalc_ks)
98 :
99 : !Arguments ------------------------------------
100 : type(dataset_type),intent(in) :: dtset
101 : integer,intent(in) :: mband
102 : integer,intent(out) :: nkcalc
103 : !arrays
104 : real(dp),allocatable,intent(out) :: kcalc(:,:)
105 : integer,allocatable,intent(out) :: bstart_ks(:,:)
106 : integer,allocatable,intent(out) :: nbcalc_ks(:,:)
107 :
108 : !Local variables ------------------------------
109 : !scalars
110 : integer :: spin, ierr, ikcalc
111 : character(len=500) :: msg
112 : ! *************************************************************************
113 :
114 27 : call wrtout(std_out, " Generating list of k-points for self-energy from kptgw and bdgw.")
115 :
116 27 : nkcalc = dtset%nkptgw
117 81 : ABI_MALLOC(kcalc, (3, nkcalc))
118 108 : ABI_MALLOC(bstart_ks, (nkcalc, dtset%nsppol))
119 81 : ABI_MALLOC(nbcalc_ks, (nkcalc, dtset%nsppol))
120 :
121 222 : kcalc = dtset%kptgw(:,1:nkcalc)
122 54 : do spin=1,dtset%nsppol
123 69 : bstart_ks(:,spin) = dtset%bdgw(1,1:nkcalc,spin)
124 96 : nbcalc_ks(:,spin) = dtset%bdgw(2,1:nkcalc,spin) - dtset%bdgw(1,1:nkcalc,spin) + 1
125 : end do
126 :
127 : ! Consistency check on bdgw and mband
128 27 : ierr = 0
129 54 : do spin=1,dtset%nsppol
130 96 : do ikcalc=1,nkcalc
131 69 : if (dtset%bdgw(2,ikcalc,spin) > mband) then
132 0 : ierr = ierr + 1
133 : write(msg,'(a,2(i0,1x),2(a,i0))')&
134 0 : "For (k, s) ",ikcalc,spin," bdgw= ",dtset%bdgw(2,ikcalc,spin), " > mband = ",mband
135 0 : ABI_WARNING(msg)
136 : end if
137 : end do
138 : end do
139 27 : ABI_CHECK(ierr == 0, "Not enough bands in WFK file. See messages above. Aborting now.")
140 :
141 27 : end subroutine sigtk_kcalc_from_nkptgw
142 : !!***
143 :
144 : !!****f* m_sigtk/sigtk_kcalc_from_qprange
145 : !! NAME
146 : !! sigtk_kcalc_from_qprange
147 : !!
148 : !! FUNCTION
149 : !! Use qprange to select the interesting k-points and the corresponding bands.
150 : !!
151 : !! 0 --> Compute the QP corrections only for the fundamental and the direct gap.
152 : !! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num`
153 : !! bands above and below the Fermi level.
154 : !! -num --> Compute the QP corrections for all the k-points in the irreducible zone.
155 : !! Include all occupied states and `num` empty states.
156 : !!
157 : !! INPUT
158 : !! dtset<dataset_type>=All input variables for this dataset.
159 : !! cryst<crystal_t>=Crystalline structure
160 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
161 : !! qprange: See above description
162 : !!
163 : !! OUTPUT
164 : !! nkcalc: Number of k-points in self-energy matrix elements.
165 : !! kcalc(3, nkcalc): List of k-points where the self-energy is computed.
166 : !! bstart_ks(nkcalc, nsppol): Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
167 : !! nbcalc_ks(nkcalc, nsppol): Number of bands included in self-energy matrix elements for each k-point in kcalc.
168 : !!
169 : !! SOURCE
170 :
171 11 : subroutine sigtk_kcalc_from_qprange(dtset, cryst, ebands, qprange, nkcalc, kcalc, bstart_ks, nbcalc_ks)
172 :
173 : !Arguments ------------------------------------
174 : type(dataset_type),intent(in) :: dtset
175 : type(crystal_t),intent(in) :: cryst
176 : type(ebands_t),intent(in) :: ebands
177 : integer,intent(in) :: qprange
178 : integer,intent(out) :: nkcalc
179 : !arrays
180 : real(dp),allocatable,intent(out) :: kcalc(:,:)
181 : integer,allocatable,intent(out) :: bstart_ks(:,:), nbcalc_ks(:,:)
182 :
183 : !Local variables ------------------------------
184 : !scalars
185 : integer :: spin, ik, bstop, mband, sigma_nkbz
186 : !arrays
187 22 : integer :: kptrlatt(3,3), val_indices(ebands%nkpt, ebands%nsppol)
188 11 : real(dp),allocatable :: sigma_wtk(:),sigma_kbz(:,:)
189 : ! *************************************************************************
190 :
191 11 : mband = ebands%mband
192 :
193 11 : val_indices = ebands%get_valence_idx()
194 :
195 32 : if (any(dtset%sigma_ngkpt /= 0)) then
196 4 : call wrtout(std_out, " Generating list of k-points for self-energy from sigma_ngkpt and qprange.")
197 4 : ABI_CHECK(qprange /= 0, "qprange must be != 0")
198 : ! Get %kcalc from sigma_ngkpt
199 4 : kptrlatt = 0
200 4 : kptrlatt(1,1) = dtset%sigma_ngkpt(1); kptrlatt(2,2) = dtset%sigma_ngkpt(2); kptrlatt(3,3) = dtset%sigma_ngkpt(3)
201 : call kpts_ibz_from_kptrlatt(cryst, kptrlatt, dtset%kptopt, dtset%sigma_nshiftk, dtset%sigma_shiftk, &
202 4 : nkcalc, kcalc, sigma_wtk, sigma_nkbz, sigma_kbz)
203 4 : ABI_FREE(sigma_kbz)
204 4 : ABI_FREE(sigma_wtk)
205 : else
206 : ! Include all the k-points in the IBZ.
207 : ! Note that kcalc == ebands%kptns so we can use a single ik index in the loop over k-points.
208 : ! No need to map kcalc onto ebands%kptns.
209 7 : call wrtout(std_out, " nkptgw set to 0 ==> Include all k-points in the IBZ for Sigma_nk.")
210 7 : nkcalc = ebands%nkpt
211 21 : ABI_MALLOC(kcalc, (3, nkcalc))
212 238 : kcalc = ebands%kptns
213 : end if
214 :
215 44 : ABI_MALLOC(bstart_ks, (nkcalc, dtset%nsppol))
216 33 : ABI_MALLOC(nbcalc_ks, (nkcalc, dtset%nsppol))
217 :
218 11 : if (qprange > 0) then
219 11 : call wrtout(std_out, " Using buffer of bands above and below the Fermi level.")
220 22 : do spin=1,dtset%nsppol
221 90 : do ik=1,nkcalc
222 68 : bstart_ks(ik,spin) = max(val_indices(ik,spin) - qprange, 1)
223 68 : bstop = min(val_indices(ik,spin) + qprange, mband)
224 79 : nbcalc_ks(ik,spin) = bstop - bstart_ks(ik,spin) + 1
225 : end do
226 : end do
227 :
228 : else
229 0 : call wrtout(std_out, " Including all occupied states and -qprange empty states.")
230 0 : bstart_ks = 1
231 0 : do spin=1,dtset%nsppol
232 0 : do ik=1,nkcalc
233 0 : nbcalc_ks(ik,spin) = min(val_indices(ik,spin) - qprange, mband)
234 : end do
235 : end do
236 : end if
237 :
238 11 : end subroutine sigtk_kcalc_from_qprange
239 : !!***
240 :
241 : !!****f* m_sigtk/sigtk_kcalc_from_gaps
242 : !! NAME
243 : !! sigtk_kcalc_from_gaps
244 : !!
245 : !! FUNCTION
246 : !! Select list of k-points and bands for self-energy matrix elements so that fundamental and direct gaps are included.
247 : !!
248 : !! INPUT
249 : !! dtset<dataset_type>=All input variables for this dataset.
250 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
251 : !! gaps<gaps_t>=Store location of gaps.
252 : !!
253 : !! OUTPUT
254 : !! nkcalc: Number of k-points in self-energy matrix elements.
255 : !! kcalc(3, nkcalc): List of k-points where the self-energy is computed.
256 : !! bstart_ks(nkcalc, nsppol): Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
257 : !! nbcalc_ks(nkcalc, nsppol): Number of bands included in self-energy matrix elements for each k-point in kcalc.
258 : !!
259 : !! SOURCE
260 :
261 14 : subroutine sigtk_kcalc_from_gaps(dtset, ebands, gaps, nkcalc, kcalc, bstart_ks, nbcalc_ks)
262 :
263 : !Arguments ------------------------------------
264 : type(dataset_type),intent(in) :: dtset
265 : type(ebands_t),intent(in) :: ebands
266 : type(gaps_t) :: gaps
267 : integer,intent(out) :: nkcalc
268 : !arrays
269 : real(dp),allocatable,intent(out) :: kcalc(:,:)
270 : integer,allocatable,intent(out) :: bstart_ks(:,:)
271 : integer,allocatable,intent(out) :: nbcalc_ks(:,:)
272 :
273 : !Local variables ------------------------------
274 : !scalars
275 : integer :: spin, nsppol, ii, ik_ibz, nk_found, ifo, jj, ib_min, ib_max
276 : logical :: found, changed
277 : !arrays
278 28 : integer :: val_indices(ebands%nkpt, ebands%nsppol), kpos(6)
279 : ! *************************************************************************
280 :
281 : ABI_UNUSED((/dtset%natom/))
282 :
283 14 : call wrtout(std_out, " Including direct and fundamental KS gap in Sigma_nk")
284 28 : ABI_CHECK(maxval(gaps%ierr) == 0, "qprange 0 cannot be used because I cannot find the gap (gap_err !=0)")
285 :
286 14 : nsppol = ebands%nsppol
287 14 : val_indices = ebands%get_valence_idx()
288 :
289 : ! Include the direct and the fundamental KS gap.
290 : ! The problem here is that kptgw and nkptgw do not depend on the spin and therefore
291 : ! we have compute the union of the k-points where the fundamental and the direct gaps are located.
292 14 : nk_found = 1; kpos(1) = gaps%fo_kpos(1,1)
293 :
294 : ! Find the list of `interesting` kpoints.
295 28 : do spin=1,nsppol
296 70 : do ifo=1,3
297 42 : ik_ibz = gaps%fo_kpos(ifo, spin)
298 42 : found = .False.; jj = 0
299 84 : do while (.not. found .and. jj < nk_found)
300 42 : jj = jj + 1; found = (kpos(jj) == ik_ibz)
301 : end do
302 56 : if (.not. found) then
303 10 : nk_found = nk_found + 1; kpos(nk_found) = ik_ibz
304 : end if
305 : end do
306 : end do
307 :
308 : ! Now we can define the list of k-points and the bands range.
309 14 : nkcalc = nk_found
310 42 : ABI_MALLOC(kcalc, (3, nkcalc))
311 56 : ABI_MALLOC(bstart_ks, (nkcalc, nsppol))
312 42 : ABI_MALLOC(nbcalc_ks, (nkcalc, nsppol))
313 :
314 38 : do ii=1,nkcalc
315 24 : ik_ibz = kpos(ii)
316 96 : kcalc(:,ii) = ebands%kptns(:,ik_ibz)
317 62 : do spin=1,nsppol
318 : ! Enlarge initial band range to include degenerate states.
319 24 : ib_min = val_indices(ik_ibz, spin)
320 24 : ib_max = ib_min + 1
321 24 : call ebands%enclose_degbands(ik_ibz, spin, ib_min, ib_max, changed, dtset%symsigma_de)
322 24 : bstart_ks(ii,spin) = ib_min
323 48 : nbcalc_ks(ii,spin) = ib_max - ib_min + 1
324 : end do
325 : end do
326 :
327 14 : end subroutine sigtk_kcalc_from_gaps
328 : !!***
329 :
330 : !!****f* m_sigtk/sigtk_kcalc_from_erange
331 : !! NAME
332 : !! sigtk_kcalc_from_erange
333 : !!
334 : !! FUNCTION
335 : !! Select list of k-points and bands for self-energy matrix elements on the basis of their positions
336 : !! wrt to the (band edges|fermi level) and the value of sigma_erange.
337 : !! Useful when computing electron-lifetimes for transport calculations.
338 : !!
339 : !! INPUT
340 : !! dtset<dataset_type>=All input variables for this dataset.
341 : !! cryst<crystal_t>=Crystalline structure
342 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
343 : !! gaps<gaps_t>=Store location of gaps.
344 : !! comm: MPI communicator.
345 : !!
346 : !! OUTPUT
347 : !! nkcalc: Number of k-points in self-energy matrix elements.
348 : !! kcalc(3, nkcalc): List of k-points where the self-energy is computed.
349 : !! bstart_ks(nkcalc, nsppol): Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
350 : !! nbcalc_ks(nkcalc, nsppol): Number of bands included in self-energy matrix elements for each k-point in kcalc.
351 : !!
352 : !! SOURCE
353 :
354 14 : subroutine sigtk_kcalc_from_erange(dtset, cryst, ebands, gaps, nkcalc, kcalc, bstart_ks, nbcalc_ks, comm)
355 :
356 : !Arguments ------------------------------------
357 : type(dataset_type),intent(in) :: dtset
358 : type(crystal_t),intent(in) :: cryst
359 : type(ebands_t),intent(in) :: ebands
360 : type(gaps_t),intent(in) :: gaps
361 : integer,intent(in) :: comm
362 : integer,intent(out) :: nkcalc
363 : !arrays
364 : real(dp),allocatable,intent(out) :: kcalc(:,:)
365 : integer,allocatable,intent(out) :: bstart_ks(:,:)
366 : integer,allocatable,intent(out) :: nbcalc_ks(:,:)
367 :
368 : !Local variables ------------------------------
369 : !scalars
370 : integer,parameter :: master = 0
371 : integer :: spin, ik, band, ii, ic, nsppol, tmp_nkpt, sigma_nkbz, my_rank
372 : logical :: found
373 : real(dp) :: cmin, vmax, ee
374 : logical :: assume_gap
375 : character(len=500) :: msg
376 14 : type(krank_t) :: krank
377 : !arrays
378 28 : integer :: kptrlatt(3,3), units(1), kpos(ebands%nkpt)
379 14 : integer,allocatable :: ib_work(:,:,:), sigmak2ebands(:), indkk(:,:)
380 14 : real(dp),allocatable :: sigma_wtk(:),sigma_kbz(:,:),tmp_kcalc(:,:)
381 : ! *************************************************************************
382 :
383 14 : my_rank = xmpi_comm_rank(comm) !; nprocs = xmpi_comm_size(comm)
384 28 : units = [std_out]
385 16 : assume_gap = .not. all(dtset%sigma_erange < zero)
386 :
387 14 : if (my_rank == master) then
388 14 : write(std_out, "(a)")" Selecting k-points and bands according to their position wrt the band edges (sigma_erange)."
389 42 : write(std_out, "(a, 2(f6.3, 1x), a)")" sigma_erange: ", dtset%sigma_erange(:) * Ha_eV, " (eV)"
390 14 : if (assume_gap) then
391 26 : call gaps%print([std_out])
392 26 : ABI_CHECK(maxval(gaps%ierr) == 0, "sigma_erange 0 cannot be used because I cannot find the gap (gap_err !=0)")
393 : end if
394 : end if
395 :
396 50 : if (any(dtset%sigma_ngkpt /= 0)) then
397 2 : call wrtout(std_out, sjoin(" Generating initial list of k-points from sigma_nkpt:", ltoa(dtset%sigma_ngkpt)))
398 : ! Get tentative tmp_nkpt and tmp_kcalc from sigma_ngkpt.
399 2 : kptrlatt = 0
400 2 : kptrlatt(1,1) = dtset%sigma_ngkpt(1); kptrlatt(2,2) = dtset%sigma_ngkpt(2); kptrlatt(3,3) = dtset%sigma_ngkpt(3)
401 : call kpts_ibz_from_kptrlatt(cryst, kptrlatt, dtset%kptopt, dtset%sigma_nshiftk, dtset%sigma_shiftk, &
402 2 : tmp_nkpt, tmp_kcalc, sigma_wtk, sigma_nkbz, sigma_kbz)
403 :
404 2 : ABI_FREE(sigma_kbz)
405 2 : ABI_FREE(sigma_wtk)
406 :
407 : ! Map tmp_kcalc to ebands%kpts
408 :
409 6 : ABI_MALLOC(indkk, (6, tmp_nkpt))
410 :
411 2 : call krank%from_kptrlatt(ebands%nkpt, ebands%kptns, ebands%kptrlatt, compute_invrank=.False.)
412 :
413 2 : if (kpts_map("symrec", ebands%kptopt, cryst, krank, tmp_nkpt, tmp_kcalc, indkk) /= 0) then
414 : write(msg, '(3a)' )&
415 0 : "At least one of the k-points could not be generated from a symmetrical one in the WFK.",ch10,&
416 0 : 'Action: check your WFK file and the value of sigma_nkpt, sigma_shiftk in the input file.'
417 0 : ABI_ERROR(msg)
418 : end if
419 :
420 2 : call krank%free()
421 :
422 6 : ABI_MALLOC(sigmak2ebands, (tmp_nkpt))
423 10 : sigmak2ebands = indkk(1, :)
424 2 : ABI_FREE(tmp_kcalc)
425 2 : ABI_FREE(indkk)
426 :
427 : else
428 : ! Include all the k-points in the IBZ in the initial list.
429 12 : call wrtout(std_out, " Generating initial list of k-points from input ebands%kptns.")
430 12 : tmp_nkpt = ebands%nkpt
431 : ! Trivial map
432 36 : ABI_MALLOC(sigmak2ebands, (tmp_nkpt))
433 10506 : sigmak2ebands = [(ii, ii=1, ebands%nkpt)]
434 : end if
435 :
436 14 : nsppol = ebands%nsppol
437 56 : ABI_MALLOC(ib_work, (2, tmp_nkpt, nsppol))
438 :
439 28 : do spin=1,nsppol
440 :
441 14 : if (assume_gap) then
442 : ! Get CBM and VBM with some tolerance
443 13 : vmax = gaps%vb_max(spin) + tol2 * eV_Ha
444 13 : cmin = gaps%cb_min(spin) - tol2 * eV_Ha
445 : else
446 1 : vmax = ebands%fermie
447 1 : cmin = ebands%fermie
448 : end if
449 :
450 3528 : do ii=1,tmp_nkpt
451 : ! Index of k-point in ebands.
452 3500 : ik = sigmak2ebands(ii)
453 : ! Will use this initial values to understand if k-point is in energy window.
454 3500 : ib_work(1, ii, spin) = huge(1)
455 3500 : ib_work(2, ii, spin) = -huge(1)
456 33688 : do band=1,ebands%nband(ik + (spin-1) * ebands%nkpt)
457 30174 : ee = ebands%eig(band, ik, spin)
458 30174 : if (abs(dtset%sigma_erange(1)) > zero) then
459 3582 : if (ee <= vmax .and. vmax - ee <= abs(dtset%sigma_erange(1))) then
460 38 : ib_work(1, ii, spin) = min(ib_work(1, ii, spin), band)
461 38 : ib_work(2, ii, spin) = max(ib_work(2, ii, spin), band)
462 : !write(std_out, *), "Adding valence band", band, " with ee [eV]: ", ee * Ha_eV
463 : end if
464 : end if
465 33674 : if (abs(dtset%sigma_erange(2)) > zero) then
466 30174 : if (ee >= cmin .and. ee - cmin <= abs(dtset%sigma_erange(2))) then
467 44 : ib_work(1, ii, spin) = min(ib_work(1, ii, spin), band)
468 44 : ib_work(2, ii, spin) = max(ib_work(2, ii, spin), band)
469 : !write(std_out, *)"Adding conduction band", band, " with ee [eV]: ", ee * Ha_eV
470 : end if
471 : end if
472 : end do
473 : end do
474 : end do
475 :
476 : ! Now we can define the list of k-points and the bands range.
477 : ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore
478 : ! we have to compute the union of the k-points.
479 14 : nkcalc = 0
480 3514 : do ii=1,tmp_nkpt
481 6940 : found = .False.
482 6940 : do spin=1,nsppol
483 6940 : if (ib_work(1, ii, spin) <= ib_work(2, ii, spin)) then
484 : found = .True.; exit
485 : end if
486 : end do
487 3514 : if (found) then
488 60 : nkcalc = nkcalc + 1
489 60 : kpos(nkcalc) = ii
490 : end if
491 : end do
492 :
493 42 : ABI_MALLOC(kcalc, (3, nkcalc))
494 56 : ABI_MALLOC(bstart_ks, (nkcalc, nsppol))
495 42 : ABI_MALLOC(nbcalc_ks, (nkcalc, nsppol))
496 :
497 74 : do ic=1,nkcalc
498 : ! Index in the ib_work array
499 60 : ii = kpos(ic)
500 : ! Index in ebands.
501 60 : ik = sigmak2ebands(ii)
502 240 : kcalc(:,ic) = ebands%kptns(:,ik)
503 134 : do spin=1,nsppol
504 60 : bstart_ks(ic, spin) = 0
505 60 : nbcalc_ks(ic, spin) = 0
506 60 : if (ib_work(1, ii, spin) <= ib_work(2, ii, spin)) then
507 60 : bstart_ks(ic, spin) = ib_work(1, ii, spin)
508 60 : nbcalc_ks(ic, spin) = ib_work(2, ii, spin) - ib_work(1, ii, spin) + 1
509 : end if
510 120 : if (nbcalc_ks(ic, spin) == 0) then
511 0 : ABI_WARNING("Spin-polarized case with nbcalc_ks == 0, don't know if code can handle it!")
512 : end if
513 : end do
514 : end do
515 :
516 14 : if (my_rank == master) then
517 : ! Write info about k-points used in the calculation.
518 : write(msg, "(a, i0, a, 2(f6.3, 1x), a)") &
519 42 : " Found ", nkcalc, " k-points within sigma_erange: ", dtset%sigma_erange(:) * Ha_eV, " (eV)"
520 14 : call wrtout(units, msg)
521 50 : if (any(dtset%sigma_ngkpt /= 0)) then
522 2 : call wrtout(units, sjoin(" These k-points belong to the sigma_ngkpt k-mesh:", ltoa(dtset%sigma_ngkpt)))
523 : end if
524 162 : write(msg, "(2(a, i0))")" min(nbcalc_ks): ", minval(nbcalc_ks), " Max(nbcalc_ks): ", maxval(nbcalc_ks)
525 14 : call wrtout(units, msg)
526 : end if
527 :
528 14 : ABI_FREE(ib_work)
529 14 : ABI_FREE(sigmak2ebands)
530 :
531 14 : end subroutine sigtk_kcalc_from_erange
532 : !!***
533 :
534 : !!****f* m_sigmaph/sigtk_kpts_in_erange
535 : !! NAME
536 : !! sigtk_kpts_in_erange
537 : !!
538 : !! FUNCTION
539 : !! Use star functions interpolation and [[einterp]] to interpolate KS energies onto dense k-mesh
540 : !! defined by [[sigma_ngkpt]] and [[sigma_shiftk]].
541 : !! find k-points inside (electron/hole) pockets according to the values specified by [[sigma_erange]].
542 : !! write kerange.nc file with the tables required by abinit to automate nscf band structure calculations
543 : !! mainly used to prepare eph calculations in which only selected k-points are nededed (imaginary part of self-energies).
544 : !!
545 : !! INPUTS
546 : !! dtset <dataset_type>=all input variables for this dataset
547 : !! cryst<crystal_t>=Crystalline structure
548 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
549 : !! psps <pseudopotential_type>=all the information about psps
550 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
551 : !! prefix=Prefix for output file.
552 : !! comm: MPI communicator.
553 : !!
554 : !! SOURCE
555 :
556 2 : subroutine sigtk_kpts_in_erange(dtset, cryst, ebands, psps, pawtab, prefix, comm)
557 :
558 : !Arguments ------------------------------------
559 : !scalars
560 : type(dataset_type),intent(in) :: dtset
561 : type(crystal_t),intent(in) :: cryst
562 : type(ebands_t),intent(in) :: ebands
563 : type(pseudopotential_type),intent(in) :: psps
564 : character(len=*),intent(in) :: prefix
565 : integer,intent(in) :: comm
566 : !arrays
567 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
568 :
569 : !Local variables ------------------------------
570 : !scalars
571 : integer,parameter :: master = 0, pertcase0 = 0, image1 = 1
572 : integer :: ii, my_rank, nprocs, spin, ikf_ibz, band, nkpt_inerange, gap_err, unt, ncid, cnt, ncerr
573 : logical :: assume_gap
574 : real(dp) :: ee, cmin, vmax
575 : character(len=500) :: msg
576 : character(len=fnlen) :: path
577 2 : type(ebands_t) :: fine_ebands
578 2 : type(gaps_t) :: gaps, fine_gaps
579 2 : type(wvl_internal_type) :: dummy_wvl
580 2 : type(hdr_type) :: fine_hdr
581 : !arrays
582 : integer :: fine_kptrlatt(3,3), band_block(2), units(2)
583 2 : integer,allocatable :: kshe_mask(:,:,:), krange2ibz(:)
584 : real(dp) :: params(4)
585 : ! *************************************************************************
586 :
587 2 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
588 :
589 : ! (-num, -num) activate treatment of metals with energy window around Efermi.
590 4 : assume_gap = .not. all(dtset%sigma_erange < zero)
591 6 : units = [std_out, ab_out]
592 :
593 2 : if (my_rank == master) then
594 2 : call wrtout(units, sjoin(ch10, repeat("=", 92)))
595 2 : call wrtout(units, " Using SKW interpolation to interpolate KS energies onto dense k-mesh.")
596 2 : call wrtout(units, sjoin(" defined by sigma_ngkpt:", trim(ltoa(dtset%sigma_ngkpt))))
597 2 : ABI_CHECK(allocated(dtset%sigma_shiftk), "sigma_nshiftk must be specified in input.")
598 2 : write(std_out, "(2a)") " and sigma_shiftk shifts:"
599 4 : do ii=1,dtset%nshiftk
600 4 : call wrtout(units, sjoin(itoa(ii), ltoa(dtset%sigma_shiftk(:, ii))))
601 : end do
602 :
603 2 : if (assume_gap) then
604 1 : call wrtout(units, " Finding k-points inside (electron/hole) pockets (assuming semiconductor).")
605 : else
606 1 : call wrtout(units, " Finding k-points inside energy window around Fermi level (assuming metal).")
607 : end if
608 6 : write(msg, "(a, 2(f6.3, 1x), a)")" Using sigma_erange: ", dtset%sigma_erange(:) * Ha_eV, " (eV)"
609 2 : call wrtout(units, msg)
610 2 : call wrtout(units, sjoin(" SKW parameters (einterp): ", ltoa(dtset%einterp)))
611 2 : call wrtout(units, sjoin(repeat("=", 92), ch10))
612 : !call ebands%print([std_out], header, prtvol=dtset%prtvol)
613 :
614 : ! Consistency check.
615 2 : if (all(dtset%sigma_erange == zero)) then
616 0 : ABI_ERROR("sigma_erange must be specified in input when calling sigtk_kpts_in_erange.")
617 : end if
618 2 : if (all(dtset%sigma_ngkpt == 0)) then
619 0 : ABI_ERROR("sigma_ngkpt must be specified in input when calling sigtk_kpts_in_erange.")
620 : end if
621 : end if
622 :
623 2 : if (assume_gap) then
624 : ! Compute gaps using input ebands.
625 1 : gaps = ebands%get_gaps(gap_err)
626 1 : if (gap_err /= 0) then
627 0 : ABI_ERROR("Cannot compute fundamental and direct gap (likely metal).")
628 : end if
629 :
630 1 : if (my_rank == master) call gaps%print(units, header="Gaps from input WFK")
631 1 : call gaps%free()
632 : else
633 1 : call wrtout(units, sjoin("Using Fermi level:", ftoa(ebands%fermie * Ha_eV, fmt="f6.2"), " (eV)"))
634 : end if
635 :
636 : ! Interpolate band energies with star functions.
637 : ! In the EPH code, we will need eigens in the IBZ to compute efermi not just energies inside pockets.
638 2 : fine_kptrlatt = 0
639 8 : do ii=1,3
640 8 : fine_kptrlatt(ii, ii) = dtset%sigma_ngkpt(ii)
641 : end do
642 6 : band_block = [1, ebands%mband]
643 10 : params = 0; params(1) = 1; params(2) = 5; if (nint(dtset%einterp(1)) == 1) params = dtset%einterp
644 :
645 : fine_ebands = ebands%interp_kmesh(cryst, params, fine_kptrlatt, &
646 2 : dtset%sigma_nshiftk, dtset%sigma_shiftk, band_block, comm)
647 176 : fine_ebands%istwfk = 1
648 :
649 2 : call fine_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol)
650 4 : call fine_ebands%print([std_out], header="FINE EBANDS", prtvol=dtset%prtvol)
651 :
652 2 : if (assume_gap) then
653 : ! Compute gaps using fine_ebands.
654 1 : fine_gaps = fine_ebands%get_gaps(gap_err)
655 1 : if (gap_err /= 0) then
656 0 : ABI_ERROR("Cannot compute fundamental and direct gap (likely metal).")
657 : end if
658 :
659 1 : if (my_rank == master) call fine_gaps%print(units, header="Gaps from SKW interpolated eigenvalues")
660 : end if
661 :
662 : ! Build new header with fine k-mesh (note kptrlatt_orig == kptrlatt)
663 : call fine_hdr%init_lowlvl(fine_ebands, psps, pawtab, dummy_wvl, abinit_version, pertcase0, &
664 : dtset%natom, dtset%nsym, dtset%nspden, dtset%ecut, dtset%pawecutdg, dtset%ecutsm, dtset%dilatmx, &
665 : dtset%intxc, dtset%ixc, dtset%stmbias, dtset%usewvl, dtset%pawcpxocc, dtset%pawspnorb, dtset%ngfft, dtset%ngfftdg, &
666 : dtset%so_psp, dtset%qptn, cryst%rprimd, cryst%xred, cryst%symrel, cryst%tnons, cryst%symafm, cryst%typat, &
667 : dtset%amu_orig(:, image1), dtset%icoulomb, &
668 : dtset%kptopt, dtset%nelect, dtset%ne_qFD, dtset%nh_qFD, dtset%ivalence, dtset%cellcharge(1), &
669 2 : fine_kptrlatt, fine_kptrlatt, dtset%sigma_nshiftk, dtset%sigma_nshiftk, dtset%sigma_shiftk, dtset%sigma_shiftk)
670 :
671 : ! Find k-points inside sigma_erange energy window.
672 : ! Set entry to the number of states inside the pocket at (ikpt, spin)
673 : ! (last index discerns between hole and electron pockets)
674 366 : ABI_ICALLOC(kshe_mask, (fine_ebands%nkpt, ebands%nsppol, 2))
675 :
676 4 : do spin=1,ebands%nsppol
677 : ! Get CBM and VBM with some tolerance.
678 2 : if (assume_gap) then
679 1 : vmax = fine_gaps%vb_max(spin) + tol2 * eV_Ha
680 1 : cmin = fine_gaps%cb_min(spin) - tol2 * eV_Ha
681 : else
682 : ! Note that we use the Fermi level from ebands instead of fine_ebands.
683 1 : vmax = ebands%fermie
684 1 : cmin = ebands%fermie
685 : end if
686 :
687 178 : do ikf_ibz=1,fine_ebands%nkpt
688 3111 : do band=1,ebands%mband
689 2952 : ee = fine_ebands%eig(band, ikf_ibz, spin)
690 : ! Check whether the interpolated eigenvalue is inside the sigma_erange window.
691 2952 : if (abs(dtset%sigma_erange(1)) > zero) then
692 2952 : if (ee <= vmax .and. vmax - ee <= abs(dtset%sigma_erange(1))) then
693 11 : kshe_mask(ikf_ibz, spin, 1) = kshe_mask(ikf_ibz, spin, 1) + 1; exit
694 : end if
695 : end if
696 3098 : if (abs(dtset%sigma_erange(2)) > zero) then
697 2941 : if (ee >= cmin .and. ee - cmin <= abs(dtset%sigma_erange(2))) then
698 6 : kshe_mask(ikf_ibz, spin, 2) = kshe_mask(ikf_ibz, spin, 2) + 1; exit
699 : end if
700 : end if
701 : end do
702 : end do
703 : end do
704 :
705 : ! Build list of k-points inside pockets. Use over dimensioned array.
706 358 : cnt = count(kshe_mask /= 0)
707 6 : ABI_MALLOC(krange2ibz, (cnt))
708 176 : cnt = 0
709 176 : do ikf_ibz=1,fine_ebands%nkpt
710 816 : if (any(kshe_mask(ikf_ibz,:,:) /= 0)) then
711 17 : cnt = cnt + 1; krange2ibz(cnt) = ikf_ibz
712 : end if
713 : end do
714 2 : nkpt_inerange = cnt
715 :
716 : ! Possible extensions that may be implemented at this level:
717 : ! 1. Find image points in the BZ?
718 : ! 2. Compute tetra and q-points for EPH calculation or use +/- wmax window and heuristic approach in sigmaph at runtime?
719 : ! 3. Compute SKW 1st and 2nd derivatives needed to treat Frohlich?
720 :
721 : ! Write output files with k-point list.
722 2 : if (my_rank == master .and. len_trim(prefix) /= 0) then
723 2 : write(std_out, "(a,i0,a,f5.1,a)")" Found: ", nkpt_inerange, " kpoints in sigma_erange energy windows. (nkeff / nkibz): ", &
724 4 : (100.0_dp * nkpt_inerange) / fine_ebands%nkpt, " [%]"
725 :
726 : ! Write text file with Abinit input variables (mainly for testing purposes).
727 2 : path = strcat(prefix, "_KERANGE")
728 2 : if (open_file(path, msg, newunit=unt, form="formatted") /= 0) then
729 0 : ABI_ERROR(msg)
730 : end if
731 2 : write(unt, "(a)")"kptopt 0"
732 2 : write(unt, "(a, i0)")"nkpt ", nkpt_inerange
733 2 : write(unt, "(a)")"kpt"
734 19 : do ii=1,nkpt_inerange
735 19 : write(unt, "(3(es16.8,1x))") fine_ebands%kptns(:, krange2ibz(ii))
736 : end do
737 2 : write(unt, "(a, i0)")"wtk"
738 19 : do ii=1,nkpt_inerange
739 19 : write(unt, "(es16.8)") fine_ebands%wtk(krange2ibz(ii))
740 : end do
741 2 : close(unt)
742 :
743 : ! Write netcdf file used to perform NSCF run and EPH calculations with eph_task = -4.
744 2 : path = strcat(prefix, "_KERANGE.nc")
745 2 : NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
746 : ! Write crystalline structure, fine_hdr and fine_ebands defined on the fine k-mesh.
747 : ! fine_ebands will be used to compare with the ab-initio NSCF eigenvalues.
748 : !
749 : ! TODO: The size of the KERANGE.nc quickly increases with the k-mesh.
750 : ! It is ~700 Mb for a ~ 300^3 grid due to occ and eigens
751 : ! But these quantities are now used in inkpts so it may be possible to avoid writing them to disk.
752 : !
753 2 : NCF_CHECK(fine_hdr%ncwrite(ncid, fform_from_ext("KERANGE.nc"), nc_define=.True.))
754 2 : NCF_CHECK(cryst%ncwrite(ncid))
755 2 : NCF_CHECK(fine_ebands%ncwrite(ncid))
756 4 : NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("nkpt_inerange", nkpt_inerange)], defmode=.True.))
757 : ! Define extra arrays.
758 : ncerr = nctk_def_arrays(ncid, [ &
759 : nctkarr_t("kshe_mask", "int", "number_of_kpoints, number_of_spins, two"), &
760 : nctkarr_t("krange2ibz", "int", "nkpt_inerange"), &
761 : nctkarr_t("sigma_erange", "dp", "two"), &
762 : nctkarr_t("einterp", "dp", "four") &
763 10 : ], defmode=.True.)
764 2 : NCF_CHECK(ncerr)
765 : ! Write extra arrays.
766 2 : NCF_CHECK(nctk_set_datamode(ncid))
767 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kshe_mask"), kshe_mask))
768 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "krange2ibz"), krange2ibz(1:nkpt_inerange)))
769 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma_erange"), dtset%sigma_erange))
770 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "einterp"), params))
771 2 : NCF_CHECK(nf90_close(ncid))
772 : end if
773 :
774 2 : ABI_FREE(kshe_mask)
775 2 : ABI_FREE(krange2ibz)
776 :
777 2 : call fine_gaps%free(); call fine_ebands%free(); call fine_hdr%free()
778 :
779 2 : end subroutine sigtk_kpts_in_erange
780 : !!***
781 :
782 12 : subroutine degtab_free(degtab)
783 : class(degtab_t),intent(inout) :: degtab
784 : integer :: ii
785 46 : do ii=1,size(degtab%bids)
786 46 : ABI_SFREE(degtab%bids(ii)%vals)
787 : end do
788 46 : ABI_FREE(degtab%bids)
789 12 : end subroutine degtab_free
790 :
791 56 : subroutine degtab_array_free(degtab)
792 : class(degtab_t),intent(inout) :: degtab(:,:)
793 :
794 : integer :: jj, ii, ideg
795 :
796 112 : do jj=1,size(degtab, dim=2)
797 284 : do ii=1,size(degtab, dim=1)
798 172 : if (.not. allocated(degtab(ii, jj)%bids)) cycle
799 630 : do ideg=1,size(degtab(ii, jj)%bids)
800 630 : ABI_SFREE(degtab(ii, jj)%bids(ideg)%vals)
801 : end do
802 686 : ABI_SFREE(degtab(ii, jj)%bids)
803 : end do
804 : end do
805 :
806 56 : end subroutine degtab_array_free
807 : !!***
808 :
809 : !----------------------------------------------------------------------
810 :
811 : !!****f* m_sigtk/sigtk_sigma_tables
812 : !! NAME
813 : !! sigtk_sigma_tables
814 : !!
815 : !! FUNCTION
816 : !! Build tables with the band indices used to compute the matrix elements of sigma_x and sigma_c
817 : !! taking into account the kind of self-energies and symmetries from esymm.
818 : !!
819 : !! INPUTS
820 : !! nkcalc: Number of k-points to compute.
821 : !! nkibz: Number of k-points in the IBZ.
822 : !! nsppol: Number of spins.
823 : !! bstart_ks, bstop_ks: First and last band for each (ikcalc, spin).
824 : !! kcalc2ibz: Mapping kcalc --> IBZ.
825 : !! only_diago: True if only diagonal matrix elements are wanted.
826 : !! sigc_is_herm: True is Sigma_c is Hermitian.
827 : !! [esymm]: Band symmetries
828 : !!
829 : !! OUTPUT
830 : !! sigxij_tab, sigcij_tab
831 : !!
832 : !! SOURCE
833 :
834 402 : subroutine sigtk_sigma_tables(nkcalc, nkibz, nsppol, bstart_ks, bstop_ks, kcalc2ibz, &
835 0 : only_diago, sigc_is_herm, sigxij_tab, sigcij_tab, esymm)
836 :
837 : !Arguments ------------------------------------
838 : integer,intent(in) :: nkcalc, nkibz, nsppol
839 : integer,intent(in) :: bstart_ks(nkcalc, nsppol), bstop_ks(nkcalc, nsppol)
840 : logical,intent(in) :: only_diago, sigc_is_herm
841 : integer,intent(in) :: kcalc2ibz(nkcalc)
842 : type(sigijtab_t),allocatable,intent(inout) :: Sigxij_tab(:,:), Sigcij_tab(:,:)
843 : type(esymm_t),optional,intent(in) :: esymm(nkibz, nsppol)
844 :
845 : !Local variables-------------------------------
846 : !scalars
847 : integer :: spin,ikcalc,ik_ibz,bmin,bmax,bcol,brow
848 : integer :: ii,idx_x,idx_c,irr_idx1,irr_idx2
849 : !arrays
850 201 : integer,allocatable :: sigc_bidx(:), sigx_bidx(:)
851 402 : logical :: use_sym_at(nkibz, nsppol)
852 : ! *************************************************************************
853 :
854 201 : if (allocated(Sigxij_tab)) then
855 0 : call sigijtab_free(Sigxij_tab)
856 0 : ABI_FREE(Sigxij_tab)
857 : end if
858 201 : if (allocated(Sigcij_tab)) then
859 0 : call sigijtab_free(Sigcij_tab)
860 0 : ABI_FREE(Sigcij_tab)
861 : end if
862 :
863 1656 : ABI_MALLOC(Sigcij_tab, (nkcalc, nsppol))
864 1455 : ABI_MALLOC(Sigxij_tab, (nkcalc, nsppol))
865 :
866 1638 : use_sym_at = .FALSE.
867 201 : if (present(esymm)) then
868 : ! Create the Sig_ij tables taking advantage of the classification of the bands.
869 0 : do spin=1,nsppol
870 0 : do ikcalc=1,nkcalc
871 0 : ik_ibz = kcalc2ibz(ikcalc)
872 0 : use_sym_at(ik_ibz, spin) = .not. esymm(ik_ibz, spin)%failed()
873 : end do
874 : end do
875 : end if
876 :
877 406 : do spin=1,nsppol
878 1053 : do ikcalc=1,nkcalc
879 647 : ik_ibz = kcalc2ibz(ikcalc)
880 :
881 852 : if (use_sym_at(ik_ibz, spin)) then
882 0 : if (only_diago) then
883 0 : ABI_ERROR("You should not be here!")
884 : end if
885 :
886 0 : bmin = bstart_ks(ikcalc, spin); bmax = bstop_ks(ikcalc, spin)
887 0 : ABI_MALLOC(Sigxij_tab(ikcalc, spin)%col, (bmin:bmax))
888 0 : ABI_MALLOC(Sigcij_tab(ikcalc, spin)%col, (bmin:bmax))
889 :
890 0 : do bcol=bmin,bmax
891 0 : ABI_MALLOC(sigc_bidx, (bmax - bmin + 1))
892 0 : ABI_MALLOC(sigx_bidx, (bmax - bmin + 1))
893 :
894 0 : if (esymm(ik_ibz,spin)%err_status /= 0) then
895 : ! Band classification failed.
896 0 : sigc_bidx = [(ii, ii=bmin, bmax)]
897 0 : idx_c = bmax - bmin + 1
898 0 : sigx_bidx = [(ii,ii=bmin,bcol)] ! Hermitian
899 0 : idx_x = bcol - bmin + 1
900 : else
901 0 : irr_idx2 = esymm(ik_ibz,spin)%b2irrep(bcol)
902 0 : idx_c = 0
903 0 : do brow=bmin,bmax
904 0 : irr_idx1 = esymm(ik_ibz,spin)%b2irrep(brow)
905 0 : if (sigc_is_herm .and. bcol < brow) CYCLE ! Only the upper triangle for HF, SEX, or COHSEX.
906 0 : if (irr_idx1 == irr_idx2) then ! same character, add this row to the list.
907 0 : idx_c = idx_c + 1
908 0 : sigc_bidx(idx_c) = brow
909 : end if
910 : end do
911 : idx_x = 0
912 0 : do brow=bmin,bcol
913 0 : irr_idx1 = esymm(ik_ibz,spin)%b2irrep(brow)
914 : if (bcol<brow) CYCLE ! Sig_x is always Hermitian.
915 0 : if (irr_idx1 == irr_idx2) then ! same character, add this row to the list.
916 0 : idx_x = idx_x +1
917 0 : sigx_bidx(idx_x) = brow
918 : end if
919 : end do
920 : end if
921 :
922 : ! Table for Sigma_x matrix elements taking into account symmetries of the bands.
923 0 : ABI_MALLOC(Sigxij_tab(ikcalc, spin)%col(bcol)%bidx, (idx_x))
924 :
925 0 : Sigxij_tab(ikcalc, spin)%col(bcol)%size1 = idx_x
926 0 : Sigxij_tab(ikcalc, spin)%col(bcol)%bidx(:) = sigx_bidx(1:idx_x)
927 : !write(std_out,*)" Sigxij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol
928 : !write(std_out,*)" size: ",idx_x,(Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(ii),ii=1,idx_x)
929 : !
930 : ! Table for Sigma_c matrix elements taking into account symmetries of the bands.
931 0 : ABI_MALLOC(Sigcij_tab(ikcalc, spin)%col(bcol)%bidx, (idx_c))
932 :
933 0 : Sigcij_tab(ikcalc, spin)%col(bcol)%size1= idx_c
934 0 : Sigcij_tab(ikcalc, spin)%col(bcol)%bidx(:) = sigc_bidx(1:idx_c)
935 : !write(std_out,*)" Sigcij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol
936 : !write(std_out,*)" size: ",idx_c,(Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(ii), ii=1,idx_c)
937 :
938 0 : ABI_FREE(sigx_bidx)
939 0 : ABI_FREE(sigc_bidx)
940 : end do ! bcol
941 :
942 : else
943 : ! Symmetries cannot be used for this (k,s).
944 647 : bmin = bstart_ks(ikcalc, spin); bmax = bstop_ks(ikcalc, spin)
945 7392 : ABI_MALLOC(Sigcij_tab (ikcalc, spin)%col, (bmin:bmax))
946 6745 : ABI_MALLOC(Sigxij_tab (ikcalc, spin)%col, (bmin:bmax))
947 :
948 647 : if (only_diago) then
949 : ! QP wavefunctions == KS, therefore only diagonal elements are calculated.
950 1837 : do bcol=bmin,bmax
951 1571 : ABI_MALLOC(Sigcij_tab(ikcalc, spin)%col(bcol)%bidx, (1:1))
952 1571 : Sigcij_tab(ikcalc, spin)%col(bcol)%size1= 1
953 1571 : Sigcij_tab(ikcalc, spin)%col(bcol)%bidx(1) = bcol
954 :
955 1571 : ABI_MALLOC(Sigxij_tab(ikcalc, spin)%col(bcol)%bidx, (1:1))
956 1571 : Sigxij_tab(ikcalc, spin)%col(bcol)%size1 = 1
957 1837 : Sigxij_tab(ikcalc, spin)%col(bcol)%bidx(1) = bcol
958 : end do
959 : else
960 : ! Use QP wavefunctions, Sigma_ij matrix is sparse but we have to classify the states in sigma.
961 : ! The only thing we can do here is filling the entire matrix taking advantage of Hermiticity (if any).
962 4261 : do bcol=bmin,bmax
963 11640 : ABI_MALLOC(Sigxij_tab(ikcalc, spin)%col(bcol)%bidx, (bcol-bmin+1))
964 3880 : Sigxij_tab(ikcalc, spin)%col(bcol)%size1= bcol-bmin+1
965 50520 : Sigxij_tab(ikcalc, spin)%col(bcol)%bidx(:) = [(ii, ii=bmin,bcol)] ! Sigma_x is Hermitian.
966 : !write(std_out,*)"Sigxij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol,Sigxij_tab(ikcalc,spin)%col(bcol)%bidx(:)
967 :
968 11640 : ABI_MALLOC(sigc_bidx, (bmax-bmin+1))
969 46640 : idx_c = 0
970 46640 : do brow=bmin,bmax
971 42760 : if (sigc_is_herm .and. bcol < brow) CYCLE ! Only the upper triangle of Sigc_ij is needed (SEX, COHSEX).
972 30860 : idx_c = idx_c +1
973 46640 : sigc_bidx(idx_c) = brow
974 : end do
975 11640 : ABI_MALLOC(Sigcij_tab(ikcalc, spin)%col(bcol)%bidx,(idx_c))
976 3880 : Sigcij_tab(ikcalc, spin)%col(bcol)%size1= idx_c
977 34740 : Sigcij_tab(ikcalc, spin)%col(bcol)%bidx(:) = sigc_bidx(1:idx_c)
978 4261 : ABI_FREE(sigc_bidx)
979 : !write(std_out,*)"Sigcij_tab: ikcalc, spin, bcol ",ikcalc,spin,bcol,Sigcij_tab(ikcalc,spin)%col(bcol)%bidx(:)
980 : end do
981 : end if
982 : end if
983 :
984 : end do !ikcalc
985 : end do !spin
986 :
987 201 : end subroutine sigtk_sigma_tables
988 : !!***
989 :
990 : !!****f* m_sigtk/sigtk_multiply_by_vc_sqrt
991 : !! NAME
992 : !! sigtk_multiply_by_vc_sqrt
993 : !!
994 : !! FUNCTION
995 : !! Multiply rhotwg vector by the square root of the Coulomb term taking into account nspinor.
996 : !!
997 : !! INPUTS
998 : !! trans="C" to take the complex conjugate of rhotwg. "N" to use rhotwg directly.
999 : !! npw=Number of PWs.
1000 : !! nspinor: Number of spinor components.
1001 : !! ndat=Number of bands in rhotwh.
1002 : !! vc_sqrt: square root of the Coulomb interaction vc(q,g).
1003 : !!
1004 : !! SIDE EFFECTS
1005 : !! rhotgw:
1006 : !! In input: <k+q|e^{-i(q+g)r|k>
1007 : !! In output: <k+q|e^{-i(q+g)r|k> * vc_sqrt(q, g)
1008 : !!
1009 : !! SOURCE
1010 :
1011 779968 : subroutine sigtk_multiply_by_vc_sqrt(trans, npw, nspinor, ndat, vc_sqrt, rhotwg)
1012 :
1013 : character(len=1),intent(in) :: trans
1014 : integer,intent(in) :: npw, nspinor, ndat
1015 : complex(gwp),intent(in) :: vc_sqrt(npw)
1016 : complex(gwp),intent(inout) :: rhotwg(npw*nspinor, ndat)
1017 :
1018 : !Local variables ------------------------------
1019 : integer :: ii, spad, idat
1020 : !************************************************************************
1021 :
1022 : select case (trans)
1023 : case ("N")
1024 472192 : do idat=1, ndat
1025 708288 : do ii=1,nspinor
1026 236096 : spad = (ii-1) * npw
1027 92831424 : rhotwg(spad+1:spad+npw, idat) = rhotwg(spad+1:spad+npw, idat) * vc_sqrt(1:npw)
1028 : end do
1029 : end do
1030 :
1031 : case ("C")
1032 : ! Take the complex conjugate of rhotwg.
1033 1087744 : do idat=1, ndat
1034 1631616 : do ii=1,nspinor
1035 543872 : spad = (ii-1) * npw
1036 228971904 : rhotwg(spad+1:spad+npw, idat) = GWPC_CONJG(rhotwg(spad+1:spad+npw, idat)) * vc_sqrt(1:npw)
1037 : end do
1038 : end do
1039 :
1040 : case default
1041 779968 : ABI_ERROR(sjoin("Invalid trans", trans))
1042 : end select
1043 :
1044 779968 : end subroutine sigtk_multiply_by_vc_sqrt
1045 : !!***
1046 :
1047 : !!****f* m_sigtk/sigtk_dw_tpp_red
1048 : !! NAME
1049 : !! sigtk_dw_tpp_red
1050 : !!
1051 : !! FUNCTION
1052 : !! Compute T_pp'(q,nu) matrix in reduced coordinates.
1053 : !!
1054 : !! INPUTS
1055 : !!
1056 : !! OUTPUTS
1057 : !!
1058 : !! SOURCE
1059 :
1060 7118 : pure subroutine sigtk_dw_tpp_red(natom, displ_red, tpp_red)
1061 :
1062 : integer,intent(in) :: natom
1063 : real(dp),intent(in) :: displ_red(2, 3, natom)
1064 : complex(dp),intent(out) :: tpp_red(3*natom,3*natom)
1065 :
1066 : !Local variables ------------------------------
1067 : integer :: ip1, ip2, idir1, idir2, ipert1, ipert2
1068 : complex(dp) :: dka, dkap, dkpa, dkpap
1069 : !************************************************************************
1070 :
1071 49826 : do ip2=1,natom*3
1072 42708 : idir2 = mod(ip2-1, 3) + 1; ipert2 = (ip2 - idir2) / 3 + 1
1073 306074 : do ip1=1,natom*3
1074 256248 : idir1 = mod(ip1-1, 3) + 1; ipert1 = (ip1 - idir1) / 3 + 1
1075 : ! (k,a) (k,a')* + (k',a) (k',a')*
1076 256248 : dka = dcmplx(displ_red(1, idir1, ipert1), displ_red(2, idir1, ipert1))
1077 256248 : dkap = dcmplx(displ_red(1, idir2, ipert1), displ_red(2, idir2, ipert1))
1078 256248 : dkpa = dcmplx(displ_red(1, idir1, ipert2), displ_red(2, idir1, ipert2))
1079 256248 : dkpap = dcmplx(displ_red(1, idir2, ipert2), displ_red(2, idir2, ipert2))
1080 298956 : tpp_red(ip1, ip2) = dka * dconjg(dkap) + dkpa * dconjg(dkpap)
1081 : end do
1082 : end do
1083 :
1084 7118 : end subroutine sigtk_dw_tpp_red
1085 : !!***
1086 :
1087 0 : end module m_sigtk
1088 : !!***
|
