Branch data 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 self-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 : : !!***
|
