Line data Source code
1 : !!****m* ABINIT/m_wkk
2 : !! NAME
3 : !! m_wkk
4 : !!
5 : !! FUNCTION
6 : !! Computation of the matrix elements of the screened interaction W_kk'.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt.
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_wkk
23 :
24 : !use, intrinsic :: iso_c_binding
25 : use defs_basis
26 : use m_abicore
27 : use m_xmpi
28 : use m_mpinfo
29 : use m_errors
30 : use m_fft
31 : use m_wfd
32 : use m_hdr
33 : use netcdf
34 : use m_nctk
35 : use m_dtset
36 : use m_dtfil
37 : use m_fstab
38 : use m_ephtk
39 :
40 : use defs_abitypes, only : mpi_type
41 : use defs_datatypes, only : pseudopotential_type
42 : use m_gwdefs, only : GW_Q0_DEFAULT, cone_gw, czero_gw
43 : use m_hide_blas, only : xgemm, xdotc
44 : use m_special_funcs, only : gaussian
45 : use m_time, only : cwtime, cwtime_report, timab, sec2str
46 : use m_fstrings, only : itoa, ftoa, sjoin, ktoa, ltoa, strcat
47 : use m_numeric_tools, only : arth, linspace ! print_arr
48 : use m_io_tools, only : iomode_from_fname
49 : use m_fftcore, only : ngfft_seq, sphereboundary, get_kg, print_ngfft
50 : use m_crystal, only : crystal_t
51 : use m_bz_mesh, only : kmesh_t, findqg0
52 : use m_htetra, only : htetra_t
53 : use m_gsphere, only : gsphere_t
54 : use m_pawtab, only : pawtab_type
55 : use m_io_screening, only : hscr_t, get_hscr_qmesh_gsph
56 : use m_vcoul, only : vcoul_t
57 : !use m_occ, only : get_fact_spin_tol_empty
58 : use m_ebands, only : ebands_t, edos_t
59 : use m_pstat, only : pstat_proc
60 : use m_sigtk, only : sigtk_multiply_by_vc_sqrt
61 : use m_screening, only : epsm1_t
62 : use m_ddk, only : ddkop_t
63 :
64 : implicit none
65 :
66 : private
67 : !!***
68 :
69 : public :: wkk_run ! Main entry point to compute Wkk' matrix elements.
70 :
71 : !----------------------------------------------------------------------
72 :
73 : contains !=====================================================
74 : !!***
75 :
76 : !----------------------------------------------------------------------
77 :
78 : !!****f* m_wkk/wkk_run
79 : !! NAME
80 : !! wkk_run
81 : !!
82 : !! FUNCTION
83 : !! Computation of the matrix elements of the screened interaction W_kk'
84 : !! between two Cooper pairs.
85 : !!
86 : !! INPUTS
87 : !! wfk0_path=String with the path to the GS unperturbed WFK file.
88 : !! dtfil<datafiles_type>=Variables related to files.
89 : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
90 : !! dtset<dataset_type>=All input variables for this dataset.
91 : !! cryst: Crystalline structure
92 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
93 : !! wfk_hdr=Header of the WFK file.
94 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
95 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
96 : !! comm=MPI communicator.
97 : !!
98 : !! OUTPUT
99 : !! Results are written to file in netcdf format.
100 : !!
101 : !! SOURCE
102 :
103 0 : subroutine wkk_run(wfk0_path, dtfil, ngfft, ngfftf, dtset, cryst, ebands, wfk_hdr, &
104 0 : pawtab, psps, mpi_enreg, comm)
105 :
106 : !Arguments ------------------------------------
107 : !scalars
108 : character(len=*),intent(in) :: wfk0_path
109 : type(datafiles_type),intent(in) :: dtfil
110 : type(dataset_type),intent(inout) :: dtset
111 : type(crystal_t),intent(in) :: cryst
112 : type(ebands_t),target,intent(in) :: ebands
113 : type(hdr_type),intent(in) :: wfk_hdr
114 : type(pseudopotential_type),intent(in) :: psps
115 : type(mpi_type),intent(inout) :: mpi_enreg
116 : integer,intent(in) :: comm
117 : !arrays
118 : integer,intent(in) :: ngfft(18), ngfftf(18)
119 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
120 :
121 : !Local variables ------------------------------
122 : !scalars
123 : integer,parameter :: LOG_MODQ = 1, LOG_MODK = 4, istwfk1 = 1, master = 0, ndat1 = 1
124 : integer :: id_required, approx_type, ikxc, option_test, nkxc, ig
125 : integer :: nkibz, my_rank, nsppol, iq_ibz, iq_bz, isym_qq, itim_qq
126 : integer :: nspinor,nprocs, ii, spin, npw_x, npw_c ! max_npw_xc, min_npw_xc, my_is, cplex, ib,
127 : integer :: bstart_k, bstop_k, nband_k, ncols, mpw, ierr !, ncerr
128 : integer :: bstart_kp, bstop_kp, nband_kp, ik_bz, bmin, bmax, max_nb
129 : integer :: ikp_ibz, isym_kp, trev_kp, npw_kp, istwf_kp, npw_kp_ibz, istwf_kp_ibz
130 : integer :: ik_ibz, isym_k, trev_k, npw_k, istwf_k, npw_k_ibz, istwf_k_ibz
131 : integer :: m_k, im_k, n_kp, in_kp, root_ncid, wkk_mode, ne
132 : integer :: nfft,nfftf,mgfft,mgfftf,nkpg_kp,nkpg_k,edos_intmeth, nqlwl, scr_iomode ! cnt,
133 : integer :: ikp_bz, my_ikp, my_nkp !, my_iq,
134 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_kp, wall_kp, gflops_kp ! cpu, wall, gflops,
135 : real(dp) :: edos_step, edos_broad, e_mk, e_nkp, e_min, e_max, e_step, smear_mk, smear_nkp, faq
136 : logical :: isirr_k, isirr_kp, qq_is_gamma, remove_exchange, print_time_kp
137 0 : type(wfd_t) :: wfd
138 0 : type(kmesh_t) :: qmesh, kmesh
139 0 : type(gsphere_t),target :: gsph_x, gsph_c
140 0 : type(hscr_t),target :: hscr
141 0 : type(vcoul_t) :: vcp
142 0 : type(edos_t) :: edos
143 0 : type(epsm1_t) :: epsm1
144 0 : type(ddkop_t) :: ddkop
145 0 : type(htetra_t) :: tetra
146 : character(len=fnlen) :: screen_filepath
147 : character(len=5000) :: msg !, qkp_string
148 : !arrays
149 : integer :: g0_k(3), g0_kp(3), g0_qq(3), units(2), work_ngfft(18), gmax(3), mapl_k(6), mapl_kp(6), mg0(3)
150 0 : integer,allocatable :: kg_kp(:,:), kg_k(:,:), gbound_kp(:,:), gbound_k(:,:), gbound_c(:,:), gbound_x(:,:)
151 0 : integer,allocatable :: nband(:,:), wfd_istwfk(:)
152 0 : integer, contiguous, pointer :: kg_c(:,:), kg_x(:,:)
153 : real(dp) :: kk_ibz(3), kk_bz(3), kp_ibz(3), kp_bz(3), qq_bz(3), kk_diff(3) ! qq_ibz(3)
154 : complex(gwp) :: ctmp_gwpc
155 : character(len=fnlen) :: path
156 : !arrays
157 0 : real(dp) :: n0(ebands%nsppol)
158 0 : real(dp),allocatable :: qlwl(:,:), kpg_kp(:,:), kpg_k(:,:), ug_kp(:,:,:), ug_k(:,:) !, cg_work(:,:)
159 0 : real(dp),allocatable :: work(:,:,:,:), e_mesh(:), e_args(:), wgt_mk(:), wgt_nkp(:,:) ! vcart_ibz(:,:,:,:),
160 0 : complex(gwp),allocatable :: cwork_ur(:), rhotwg_mn_x(:,:,:), rhotwg_mn_c(:,:,:), w_rhotwg_mn_c(:,:,:)
161 0 : complex(gwp),allocatable :: vc_sqrt_gx(:), ur_nkp(:,:), ur_mk(:,:), kxcg(:,:), mu_mn(:,:)
162 0 : complex(dp),allocatable :: w_ee(:,:) ! w_kkp(:,:,:,:)
163 0 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
164 0 : type(fstab_t),target,allocatable :: fstab(:)
165 : !************************************************************************
166 :
167 0 : ABI_CHECK(psps%usepaw == 0, "PAW not implemented")
168 0 : ABI_CHECK(dtset%nsppol == 1, "nsppol 2 implemented!")
169 0 : ABI_CHECK(dtset%nspinor == 1, "nspinor 2 implemented!")
170 :
171 : ABI_UNUSED(mpi_enreg%nproc)
172 :
173 0 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm); units = [std_out, ab_out]
174 0 : call cwtime(cpu_all, wall_all, gflops_all, "start")
175 :
176 0 : nsppol = ebands%nsppol; nspinor = ebands%nspinor
177 :
178 : ! Compute electron DOS.
179 0 : edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos
180 : edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear
181 0 : edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha
182 0 : edos = ebands%get_edos(cryst, edos_intmeth, edos_step, edos_broad, comm)
183 :
184 : ! Get DOS per spin channel
185 0 : n0(:) = edos%gef(1:edos%nsppol)
186 0 : if (my_rank == master) then
187 0 : call edos%print(units)
188 0 : path = strcat(dtfil%filnam_ds(4), "_EDOS")
189 0 : call edos%write(path)
190 : end if
191 :
192 : ! Find Fermi surface k-points.
193 : ! TODO: support kptopt, change setup of k-points if tetra: fist tetra weights then k-points on the Fermi surface!
194 0 : ABI_MALLOC(fstab, (nsppol))
195 0 : call fstab_init(fstab, ebands, cryst, dtset, tetra, comm)
196 0 : call tetra%free()
197 :
198 0 : bmin = huge(1); bmax = -1
199 0 : do spin=1,nsppol
200 0 : bmin = min(bmin, fstab(spin)%bmin)
201 0 : bmax = max(bmax, fstab(spin)%bmax)
202 : end do
203 0 : max_nb = bmax - bmin + 1
204 :
205 : ! Build linear mesh for W(e,e').
206 0 : spin = 1
207 : associate (fs => fstab(spin))
208 0 : e_min = minval(ebands%eig(fs%bmin,:,spin)) - 0.1_dp * eV_Ha
209 0 : e_max = maxval(ebands%eig(fs%bmax,:,spin)) + 0.1_dp * eV_Ha
210 0 : e_step = 0.002 * eV_Ha
211 0 : ne = 1 + (e_max - e_min) / e_step
212 0 : ABI_CALLOC(w_ee, (ne, ne))
213 0 : ABI_MALLOC(e_mesh, (ne))
214 0 : ABI_MALLOC(e_args, (ne))
215 0 : ABI_MALLOC(wgt_mk, (ne))
216 0 : ABI_MALLOC(wgt_nkp, (ne, max_nb))
217 0 : e_mesh = linspace(e_min, e_max, ne)
218 : end associate
219 :
220 : ! ================
221 : ! HANDLE SCREENING
222 : ! ================
223 0 : screen_filepath = dtfil%fnameabi_scr
224 0 : ABI_CHECK(dtfil%fnameabi_scr /= ABI_NOFILE, "SCR file must be specified in input!")
225 :
226 : ! Read g-sphere for correlation and qmesh from SCR file.
227 0 : call get_hscr_qmesh_gsph(screen_filepath, dtset, cryst, hscr, qmesh, gsph_c, qlwl, comm)
228 0 : call hscr%print(units, dtset%prtvol, header="Header of the SCR file")
229 :
230 0 : nqlwl = size(qlwl, dim=2)
231 0 : if (nqlwl == 0) then
232 0 : nqlwl=1
233 0 : ABI_MALLOC(qlwl,(3,nqlwl))
234 0 : qlwl(:,nqlwl)= GW_Q0_DEFAULT
235 : write(msg,'(3a,i0,a,3f9.6)')&
236 0 : "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,&
237 0 : "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1)
238 0 : ABI_COMMENT(msg)
239 : end if
240 :
241 : ! Init g-sphere for the exchange part from ecutsigx.
242 0 : call gsph_c%extend(cryst, dtset%ecutsigx, gsph_x)
243 :
244 : ! Initialize Coulomb term on the IBZ of the qmesh. Use largest G-sphere.
245 0 : call kmesh%init(cryst, wfk_hdr%nkpt, wfk_hdr%kptns, dtset%kptopt)
246 0 : faq = one / (cryst%ucvol*kmesh%nbz)
247 :
248 : ! TODO:
249 : ! Here we sort the pp_mesh by stars so that we can split the pp wavevectors in blocks and therefore
250 : ! reduce the number of wavevectors in the IBZ that must be stored in memory.
251 : !call qmesh%pack_by_stars()
252 :
253 0 : npw_x = gsph_x%ng; npw_c = gsph_c%ng !; min_npw_xc = min(npw_x, npw_c); max_npw_xc = max(npw_x, npw_c)
254 0 : if (gsph_x%ng >= gsph_c%ng) then
255 : call vcp%init(gsph_x, cryst, qmesh, kmesh, dtset%rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, gsph_x%ng, &
256 0 : nqlwl, qlwl, comm)
257 : else
258 : call vcp%init(gsph_c, cryst, qmesh, kmesh, dtset%rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, gsph_c%ng, &
259 0 : nqlwl, qlwl, comm)
260 : end if
261 :
262 0 : ABI_FREE(qlwl)
263 :
264 0 : if (my_rank == master) then
265 0 : call kmesh%print(units, header="K-mesh for wavefunctions", prtvol=dtset%prtvol)
266 0 : call gsph_x%print(units, dtset%prtvol, header="G-sphere for exchange")
267 0 : call gsph_c%print(units, dtset%prtvol, header="G-sphere for correlation")
268 0 : call vcp%print(units, prtvol=dtset%prtvol)
269 : end if
270 0 : call kmesh%free()
271 :
272 0 : ABI_CHECK_IGE(npw_x, 1, "npw_x <= 1")
273 0 : ABI_CHECK_IGE(npw_c, 1, "npw_c <= 1")
274 :
275 : ! Note that in this case, the sphere is always Gamma-centered i.e. it does not depend on the qq wavevector
276 0 : kg_c => gsph_c%gvec(:, 1:npw_c)
277 0 : kg_x => gsph_x%gvec(:, 1:npw_x)
278 :
279 : ! Initialize the wave function descriptor.
280 : ! Only wavefunctions on the FS are stored in wfd.
281 : ! Need all k-points on the FS because of k+q, spin is not distributed for the time being.
282 : ! It would be possible to reduce the memory allocated per MPI-rank via OpenMP.
283 0 : nkibz = wfk_hdr%nkpt
284 0 : ABI_MALLOC(nband, (nkibz, nsppol))
285 0 : ABI_MALLOC(bks_mask, (dtset%mband, nkibz, nsppol))
286 0 : ABI_MALLOC(keep_ur, (dtset%mband, nkibz, nsppol))
287 0 : nband = dtset%mband; bks_mask = .False.; keep_ur = .False.
288 :
289 0 : do spin=1,ebands%nsppol
290 0 : associate (fs => fstab(spin))
291 0 : do ik_bz=1,fs%nkfs
292 0 : ik_ibz = fs%indkk_fs(1, ik_bz)
293 0 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz)
294 0 : bks_mask(bstart_k:bstart_k+nband_k-1, ik_ibz, spin) = .True.
295 : end do
296 : end associate
297 : end do
298 :
299 : ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts)
300 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
301 0 : ABI_MALLOC(wfd_istwfk, (nkibz))
302 0 : wfd_istwfk = 1
303 :
304 : call wfd%init(cryst, pawtab, psps, keep_ur, dtset%mband, nband, nkibz, nsppol, bks_mask,&
305 : dtset%nspden, nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,&
306 0 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
307 :
308 0 : call pstat_proc%print(_PSTAT_ARGS_)
309 :
310 0 : ABI_FREE(nband)
311 0 : ABI_FREE(bks_mask)
312 0 : ABI_FREE(keep_ur)
313 0 : ABI_FREE(wfd_istwfk)
314 :
315 : ! Read wavefunctions.
316 0 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
317 :
318 : ! TODO
319 : ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
320 : ! We also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
321 : ! that will be used to symmetrize the wavefunctions in G-space.
322 : ! This is the maximum number of PWs for all possible k+q treated.
323 :
324 0 : call ephtk_get_mpw_gmax(nkibz, wfk_hdr%kptns, dtset%ecut, cryst%gmet, mpw, gmax, comm)
325 : ! FIXME
326 : !gmax = 2 * gmax
327 :
328 0 : call ddkop%init(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft)
329 :
330 : #if 0
331 : call cwtime(cpu, wall, gflops, "start", msg=" Computing v_nk matrix elements for all states on the FS...")
332 :
333 : ABI_CALLOC(vcart_ibz, (3, bmin:bmax, nkibz, nsppol))
334 : ABI_MALLOC(cg_work, (2, mpw * nspinor))
335 :
336 : cnt = 0
337 : do spin=1,nsppol
338 : associate (fs => fstab(spin))
339 : do ik_ibz=1,ebands%nkpt
340 : kk_ibz = ebands%kptns(:, ik_ibz)
341 : npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz)
342 : ! NB: The two checks below are global --> all procs will cycle.
343 : if (all(bks_mask(:, ik_ibz, spin) .eqv. .False.)) cycle
344 : if (npw_k == 1) cycle
345 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
346 :
347 : call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk_ibz, istwf_k, npw_k, wfd%kdata(ik_ibz)%kg_k)
348 :
349 : do band_k=fs%bmin,fs%bmax
350 : if (.not. bks_mask(band_k, ik_ibz, spin)) cycle
351 : call wfd%copy_cg(band_k, ik_ibz, spin, cg_work)
352 : eig0nk = ebands%eig(band_k, ik_ibz, spin)
353 : vk = ddkop%get_vdiag(eig0nk, istwf_k, npw_k, wfd%nspinor, cg_work, cwaveprj0)
354 : vcart_ibz(:, band_k, ik_ibz, spin) = vk
355 : end do
356 : end do
357 : end associate
358 : end do ! spin
359 :
360 : call xmpi_sum(vcart_ibz, comm, ierr)
361 : ABI_FREE(cg_work)
362 : ABI_FREE(vcart_ibz)
363 : call cwtime_report(" Velocities", cpu, wall, gflops)
364 : #endif
365 0 : call ddkop%free()
366 :
367 0 : call ngfft_seq(work_ngfft, gmax)
368 0 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
369 :
370 : ! FFT meshes from input file, not necessarily equal to the ones found in the external files.
371 : ! NB: ur arrays are always allocated with nfft and not with product(ngfft(4:6)).
372 0 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
373 0 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
374 :
375 0 : call print_ngfft([std_out], ngfft, header="FFT mesh")
376 :
377 : ! Set the FFT mesh
378 0 : call wfd%change_ngfft(cryst, psps, ngfft)
379 0 : call wfd%print(units, header="Wavefunctions for GWPT calculation.")
380 :
381 0 : ABI_MALLOC(gbound_kp, (2*mgfft+8, 2))
382 0 : ABI_MALLOC(gbound_k, (2*mgfft+8, 2))
383 0 : ABI_MALLOC(gbound_c, (2*mgfft+8, 2))
384 0 : ABI_MALLOC(gbound_x, (2*mgfft+8, 2))
385 :
386 : ! ====================================
387 : ! This is the g-sphere for W_{gg'}(qq)
388 : ! ====================================
389 : ! Init g-sphere for the exchange part from ecutsigx.
390 0 : call sphereboundary(gbound_c, istwfk1, kg_c, mgfft, npw_c)
391 0 : call sphereboundary(gbound_x, istwfk1, kg_x, mgfft, npw_x)
392 :
393 0 : ABI_MALLOC(vc_sqrt_gx, (npw_x))
394 :
395 : ! Allocate g-vectors centered on k, k'
396 0 : ABI_MALLOC(kg_k, (3, mpw))
397 0 : ABI_MALLOC(kg_kp, (3, mpw))
398 :
399 0 : spin = 1
400 : associate (fs => fstab(spin))
401 0 : ABI_MALLOC_OR_DIE(ur_nkp, (nfft*nspinor, fs%bmax-fs%bmin+1), ierr)
402 0 : ABI_MALLOC_OR_DIE(ur_mk, (nfft*nspinor, fs%bmax-fs%bmin+1), ierr)
403 0 : ABI_MALLOC(cwork_ur, (nfft*nspinor))
404 :
405 : ! Read symmetrized em1 from file
406 0 : id_required = 4; ikxc = 0; approx_type = 0; option_test = 0; nkxc = 0
407 0 : ABI_MALLOC(kxcg, (nfftf, nkxc))
408 0 : scr_iomode = iomode_from_fname(screen_filepath)
409 :
410 0 : call epsm1%from_file(screen_filepath, qmesh%nibz, npw_c, comm)
411 : call epsm1%mkdump(vcp, npw_c, kg_c, nkxc, kxcg, id_required, approx_type, &
412 0 : ikxc, option_test, dtfil%fnameabo_scr, scr_iomode, nfftf, ngfft, comm)
413 :
414 0 : call epsm1%malloc_epsm1_qbz(npw_c, epsm1%nomega)
415 0 : ABI_FREE(kxcg)
416 :
417 0 : call pstat_proc%print(_PSTAT_ARGS_)
418 :
419 0 : wkk_mode = 1
420 : ! Master creates the netcdf file used to store the results of the calculation.
421 0 : if (my_rank == master) then
422 0 : NCF_CHECK(nctk_open_create(root_ncid, strcat(dtfil%filnam_ds(4), "_WKK.nc") , xmpi_comm_self))
423 0 : NCF_CHECK(cryst%ncwrite(root_ncid))
424 0 : NCF_CHECK(ebands%ncwrite(root_ncid))
425 0 : NCF_CHECK(edos%ncwrite(root_ncid))
426 :
427 : !ncerr = nctk_def_dims(ncid, [ &
428 : ! nctkdim_t("smat_bsize1", smat_bsize1), nctkdim_t("smat_bsize2", smat_bsize2) &
429 : ! ], defmode=.True.)
430 : !NCF_CHECK(ncerr)
431 : !ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
432 : ! "gwr_completed", "scf_iteration" &
433 : !])
434 : !NCF_CHECK(ncerr)
435 : !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
436 : ! "wr_step", "gwr_boxcutmin", "cosft_duality_error", "regterm" &
437 : !])
438 : !NCF_CHECK(ncerr)
439 :
440 : !ncerr = nctk_def_arrays(ncid, [ &
441 : ! nctkarr_t("gwr_task", "char", "character_string_length"), &
442 : !])
443 : !NCF_CHECK(ncerr)
444 : !NCF_CHECK(nctk_set_datamode(ncid))
445 :
446 : !ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
447 : ! "gwr_completed", "sig_diago", "b1gw", "b2gw", "symsigma", "symchi", "scf_iteration"], &
448 : ! [0, merge(1, 0, gwr%sig_diago), gwr%b1gw, gwr%b2gw, gwr%dtset%symsigma, dtset%symchi, gwr%scf_iteration])
449 : !NCF_CHECK(ncerr)
450 :
451 : !ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
452 : ! "wr_step", "ecuteps", "ecut", "ecutwfn", "ecutsigx", "gwr_boxcutmin", &
453 : ! "cosft_duality_error", "regterm"], &
454 : ! [gwr%wr_step, dtset%ecuteps, dtset%ecut, dtset%ecutwfn, dtset%ecutsigx, dtset%gwr_boxcutmin, &
455 : ! gwr%ft_max_error(1), gwr%ft_max_error(2), gwr%ft_max_error(3), gwr%cosft_duality_error, regterm &
456 : ! ])
457 : !NCF_CHECK(ncerr)
458 : !NCF_CHECK(nf90_put_var(ncid, vid("gwr_task"), trim(dtset%gwr_task)))
459 :
460 0 : NCF_CHECK(nf90_close(root_ncid))
461 : end if
462 0 : call xmpi_barrier(comm)
463 :
464 : ! Open WKK.nc file and go to data mode.
465 0 : NCF_CHECK(nctk_open_modify(root_ncid, strcat(dtfil%filnam_ds(4), "_WKK.nc") , comm))
466 0 : NCF_CHECK(nctk_set_datamode(root_ncid))
467 :
468 0 : mg0 = [1, 1, 1]
469 0 : remove_exchange = .True.
470 :
471 0 : ABI_CALLOC(mu_mn, (bmin:bmax, bmin:bmax))
472 :
473 : ! Loop over k'-points in the energy window.
474 0 : do ikp_bz=1,fs%nkfs
475 0 : my_ikp = ikp_bz; my_nkp = fs%nkfs
476 :
477 0 : print_time_kp = my_rank == 0 .and. (ikp_bz <= LOG_MODK .or. mod(ikp_bz, LOG_MODK) == 0)
478 0 : if (print_time_kp) call cwtime(cpu_kp, wall_kp, gflops_kp, "start")
479 :
480 : ! The k-point and the symmetries relating the BZ k-point to the IBZ.
481 0 : kp_bz = fs%kpts(:, ikp_bz)
482 0 : ikp_ibz = fs%indkk_fs(1, ikp_bz) ; isym_kp = fs%indkk_fs(2, ikp_bz)
483 0 : trev_kp = fs%indkk_fs(6, ikp_bz); g0_kp = fs%indkk_fs(3:5,ikp_bz)
484 0 : isirr_kp = (isym_kp == 1 .and. trev_kp == 0 .and. all(g0_kp == 0))
485 0 : mapl_kp = fs%indkk_fs(:, ikp_bz)
486 :
487 0 : kp_ibz = ebands%kptns(:,ikp_ibz)
488 0 : istwf_kp_ibz = wfd%istwfk(ikp_ibz); npw_kp_ibz = wfd%npwarr(ikp_ibz)
489 :
490 : ! Get npw_kp, kg_kp for this kp_bz.
491 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kp_bz, ikp_ibz, isirr_kp, dtset%nloalg, & ! in
492 0 : istwf_kp, npw_kp, kg_kp, nkpg_kp, kpg_kp, gbound_kp) ! out
493 :
494 : ! Rotate from IBZ to BZ and compute ur_nkp for all n_kp bands.
495 0 : bstart_kp = fs%bstart_cnt_ibz(1, ikp_ibz); nband_kp = fs%bstart_cnt_ibz(2, ikp_ibz); bstop_kp = bstart_kp + nband_kp - 1
496 0 : ABI_CALLOC(ug_kp, (2, npw_kp*nspinor, nband_kp))
497 : call wfd%rotate_cg(bstart_kp, nband_kp, spin, kp_ibz, npw_kp, kg_kp, istwf_kp, &
498 0 : cryst, mapl_kp, gbound_kp, work_ngfft, work, ug_kp, urs_kbz=ur_nkp)
499 0 : do ii=1, nband_kp
500 : !print *, "isirr_kp:", isirr_kp
501 : !print *, "g:", sum(ug_kp(1,:,ii)**2 + ug_kp(2,:,ii)**2)
502 : !print *, "r:", sum(abs(ur_nkp(:,ii)) ** 2) / nfft ! * cryst%ucvol /
503 : end do
504 :
505 : ! Loop over k-points in the energy window.
506 0 : do ik_bz=1,fs%nkfs
507 :
508 0 : mapl_k = fs%indkk_fs(:, ik_bz)
509 0 : ik_ibz = mapl_k(1); isym_k = mapl_k(2); trev_k = mapl_k(6); g0_k = mapl_k(3:5)
510 0 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
511 0 : kk_ibz = ebands%kptns(:, ik_ibz)
512 0 : istwf_k_ibz = wfd%istwfk(ik_ibz); npw_k_ibz = wfd%npwarr(ik_ibz)
513 :
514 0 : kk_bz = fs%kpts(:, ik_bz)
515 0 : kk_ibz = ebands%kptns(:,ik_ibz)
516 :
517 : ! Get npw_k, kg_k for this k.
518 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kk_bz, ik_ibz, isirr_k, dtset%nloalg, & ! in
519 0 : istwf_k, npw_k, kg_k, nkpg_k, kpg_k, gbound_k) ! out
520 0 : ABI_FREE(kpg_k)
521 :
522 0 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz); bstop_k = bstart_k + nband_k - 1
523 0 : ABI_MALLOC(ug_k, (2, npw_k*nspinor))
524 0 : ABI_MALLOC(rhotwg_mn_x, (npw_x*nspinor, nband_k, nband_kp))
525 :
526 0 : do m_k=bstart_k, bstop_k
527 : ! Rotate from IBZ to BZ and compute ur_mk
528 0 : im_k = m_k - bstart_k + 1
529 : call wfd%rotate_cg(m_k, ndat1, spin, kk_ibz, npw_k, kg_k, istwf_k, &
530 0 : cryst, mapl_k, gbound_k, work_ngfft, work, ug_k, urs_kbz=ur_mk(:,im_k))
531 :
532 0 : do n_kp=bstart_kp, bstop_kp
533 0 : in_kp = n_kp - bstart_kp + 1
534 0 : cwork_ur = conjg(ur_nkp(:,in_kp)) * ur_mk(:,im_k)
535 0 : call fft_ur(npw_x, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_x, gbound_x, cwork_ur, rhotwg_mn_x(:, im_k, in_kp))
536 0 : call sigtk_multiply_by_vc_sqrt("N", npw_x, nspinor, ndat1, vc_sqrt_gx, rhotwg_mn_x(:, im_k, in_kp))
537 : end do
538 : end do
539 :
540 0 : ABI_FREE(ug_k)
541 :
542 : ! Identify q and G0 where q + G0 = k_GW - k_i
543 0 : kk_diff = kk_bz - kp_bz
544 :
545 0 : call findqg0(iq_bz, g0_qq, kk_diff, qmesh%nbz, qmesh%bz, mG0)
546 :
547 : ! Find the corresponding irred qq-point in qmesh.
548 0 : call qmesh%get_bz_item(iq_bz, qq_bz, iq_ibz, isym_qq, itim_qq)
549 0 : qq_is_gamma = sum(qq_bz**2) < tol14
550 :
551 : ! Find k + q in the extended zone and extract symmetry info.
552 : ! Be careful here because there are two umklapp vectors to be considered as:
553 : !
554 : ! k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz
555 : !
556 : ! Get Fourier components of the Coulomb interaction in the BZ
557 : ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2
558 : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
559 : ! NOTE: vc_sqrt_gx is dimensioned with npw_x --> use rottb table from gsph_x.
560 0 : do ig=1,npw_x
561 0 : vc_sqrt_gx(gsph_x%rottb(ig, itim_qq, isym_qq)) = vcp%vc_sqrt(ig, iq_ibz)
562 : end do
563 :
564 0 : call epsm1%rotate_iqbz(iq_bz, hscr%nomega, npw_c, gsph_c, qmesh, remove_exchange)
565 :
566 : ! Compute \sum_{gg'} (M_g^{kn,k'm})^* W_{gg'}(q) M_g'^{kn,k'm} with k' = k + q.
567 : ! TODO: nspinor
568 0 : ABI_MALLOC(rhotwg_mn_c, (npw_c*nspinor, nband_k, nband_kp))
569 0 : rhotwg_mn_c(:,:,:) = rhotwg_mn_x(1:npw_c,:,:)
570 :
571 0 : ncols = nband_k * nband_kp
572 0 : ABI_MALLOC(w_rhotwg_mn_c, (npw_c*nspinor, nband_k, nband_kp))
573 :
574 : call xgemm("N", "N", npw_c, ncols, npw_c, cone_gw, epsm1%epsm1_qbz(:,:,1), npw_c, rhotwg_mn_c(:,:,1), npw_c*nspinor, &
575 0 : czero_gw, w_rhotwg_mn_c(:,:,1), npw_c)
576 :
577 : !do in_kp=1, nband_k
578 : ! do im_k=1, nband_kq
579 : ! w_rhotwg_mn_c(:,in_kp,im_k) = matmul(epsm1%epsm1_qbz(:,:,1), rhotwg_mn_c(:,im_k,in_kp))
580 : ! end do
581 : !end do
582 : !print *, " rhotwg_mn_c", maxval(abs(rhotwg_mn_c)), maxloc(abs(rhotwg_mn_c))
583 : !print *, "w_rhotwg_mn_c:", maxval(abs(w_rhotwg_mn_c))
584 :
585 : ! Precompute delta(e - e_nkp - ef).
586 0 : do n_kp=bstart_kp, bstop_kp
587 0 : in_kp = n_kp - bstart_kp + 1
588 0 : e_nkp = ebands%eig(n_kp, ik_ibz, spin)
589 0 : e_args = e_mesh - e_nkp
590 0 : smear_nkp = 0.1_dp * eV_Ha
591 0 : wgt_nkp(:, im_k) = gaussian(e_args, smear_nkp)
592 : end do
593 :
594 0 : do n_kp=bstart_kp, bstop_kp
595 0 : in_kp = n_kp - bstart_kp + 1
596 0 : e_nkp = ebands%eig(n_kp, ikp_ibz, spin)
597 : !fs%tetra_wtk(n_kp - fs%bmin + 1, ikp_ibz)
598 0 : do m_k=bstart_k, bstop_k
599 0 : im_k = m_k - bstart_k + 1
600 0 : e_mk = ebands%eig(m_k, ik_ibz, spin)
601 : !fs%tetra_wtk(m_k - fs%bmin + 1, ik_ibz)
602 0 : smear_mk = 0.1_dp * eV_Ha
603 0 : ctmp_gwpc = xdotc(npw_c*nspinor, rhotwg_mn_c(:,im_k,in_kp), 1, w_rhotwg_mn_c(:,im_k,in_kp), 1)
604 : if (remove_exchange) then
605 0 : ctmp_gwpc = ctmp_gwpc + xdotc(npw_x*nspinor, rhotwg_mn_x(:,im_k,in_kp), 1, rhotwg_mn_x(:,im_k,in_kp), 1)
606 : end if
607 0 : ctmp_gwpc = faq * ctmp_gwpc
608 : !print *, "ctmp_gwpc:", ctmp_gwpc
609 : !mu_mn(m_k, n_kp) = mu_mn(m_k, n_kp) + &
610 : ! ctmp_gwpc * gaussian(e_mk - ebands%fermie, smear_mk) * gaussian(e_mp, smear_mk)
611 : ! Computes: A := A + alpha * x * y^T
612 : !w_ee =
613 : !call dger(ne, ne, alpha, x, incx, y, incy, w_ee, lda)
614 : end do
615 : end do
616 :
617 0 : ABI_FREE(rhotwg_mn_x)
618 0 : ABI_FREE(rhotwg_mn_c)
619 0 : ABI_FREE(w_rhotwg_mn_c)
620 : end do ! ik_bz
621 :
622 0 : if (print_time_kp) then
623 0 : call inds2str(1, "My kp-point", my_ikp, my_nkp, fs%nkfs, msg)
624 0 : call cwtime_report(msg, cpu_kp, wall_kp, gflops_kp); if (ikp_bz == LOG_MODK) call wrtout(std_out, "...", do_flush=.True.)
625 : end if
626 :
627 0 : ABI_FREE(ug_kp)
628 0 : ABI_FREE(kpg_kp)
629 : end do ! ikp_bz
630 : end associate
631 :
632 0 : ABI_FREE(ur_nkp)
633 0 : ABI_FREE(ur_mk)
634 : !end do ! my_is
635 :
636 0 : call cwtime_report(" wkk calculation", cpu_all, wall_all, gflops_all, end_str=ch10)
637 :
638 0 : call xmpi_sum(w_ee, comm, ierr)
639 0 : call xmpi_sum(mu_mn, comm, ierr)
640 :
641 : if (my_rank == master) then
642 : !do m_k=bmin, bmax
643 : ! do n_kp=bmin, bmax
644 : ! call print_arr(units, mu_mn, max_r=max_nb, max_c=max_nb)
645 : ! end do
646 : !end do
647 : !call wrtout(units, sjoin("Re mu: ", ftoa(sum(real(mu_mn)))))
648 : !call wrtout(units, sjoin("Im mu: ", ftoa(sum(aimag(mu_mn)))))
649 : end if
650 :
651 : ! Set vqk_completed to 1 so that we can easily check if restarted is needed.
652 : !if (wkk_mode /= 0) then
653 : !if (my_rank == master) then
654 : ! NCF_CHECK(nf90_put_var(root_ncid, root_vid("vqk_completed"), 1))
655 : !end if
656 0 : NCF_CHECK(nf90_close(root_ncid))
657 0 : call xmpi_barrier(comm)
658 : !end if
659 :
660 : ! Output some of the results to ab_out for testing purposes
661 : !call gstore%print_for_abitests(dtset)
662 :
663 : ! Free memory
664 0 : ABI_FREE(kg_kp)
665 0 : ABI_FREE(kg_k)
666 0 : ABI_FREE(cwork_ur)
667 0 : ABI_FREE(work)
668 0 : ABI_FREE(gbound_kp)
669 0 : ABI_FREE(gbound_k)
670 0 : ABI_FREE(gbound_c)
671 0 : ABI_FREE(gbound_x)
672 0 : ABI_FREE(vc_sqrt_gx)
673 0 : ABI_FREE(w_ee)
674 0 : ABI_FREE(e_mesh)
675 0 : ABI_FREE(e_args)
676 0 : ABI_FREE(wgt_mk)
677 0 : ABI_FREE(wgt_nkp)
678 0 : ABI_FREE(mu_mn)
679 :
680 0 : call wfd%free(); call vcp%free(); call qmesh%free(); call gsph_x%free(); call gsph_c%free(); call hscr%free(); call edos%free()
681 0 : call epsm1%free()
682 :
683 0 : do spin=1,ebands%nsppol
684 0 : call fstab(spin)%free()
685 : end do
686 0 : ABI_FREE(fstab)
687 :
688 0 : call xmpi_barrier(comm) ! This to make sure that the parallel output of WKK is completed
689 0 : call cwtime_report(" wkk_run: MPI barrier before returning.", cpu_all, wall_all, gflops_all, end_str=ch10, comm=comm)
690 :
691 : contains
692 :
693 0 : subroutine inds2str(ilevel, prefix, my_ik, my_nk, nk_tot, out_str)
694 : character(len=*),intent(in) :: prefix
695 : integer,intent(in) :: ilevel, my_ik, my_nk, nk_tot
696 : character(len=*),intent(out) :: out_str
697 :
698 0 : out_str = sjoin(prefix, itoa(my_ik), "/", itoa(my_nk), "[", itoa(nk_tot), "]")
699 0 : out_str = repeat(' ', 4 * ilevel) // trim(out_str)
700 0 : end subroutine inds2str
701 :
702 : integer function root_vid(var_name)
703 : character(len=*),intent(in) :: var_name
704 : root_vid = nctk_idname(root_ncid, var_name)
705 : end function root_vid
706 :
707 : !integer function spin_vid(var_name)
708 : ! character(len=*),intent(in) :: var_name
709 : ! spin_vid = nctk_idname(spin_ncid, var_name)
710 : !end function spin_vid
711 :
712 : end subroutine wkk_run
713 : !!***
714 :
715 : end module m_wkk
716 : !!***
|
