Line data Source code
1 : !!****m* ABINIT/m_gwpt
2 : !! NAME
3 : !! m_gwpt
4 : !!
5 : !! FUNCTION
6 : !! Compute electron-phonon matrix elements within the GWPT formalism.
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_gwpt
23 :
24 : use, intrinsic :: iso_c_binding
25 : use defs_basis
26 : use m_abicore
27 : USE_MPI
28 : use m_xmpi
29 : use m_mpinfo
30 : use m_errors
31 : use m_clib
32 : use m_hide_blas
33 : use m_copy
34 : use m_fft
35 : use m_hamiltonian
36 : use m_pawcprj
37 : use m_sort
38 : use m_sigtk
39 : use m_ephtk
40 : use netcdf
41 : use m_nctk
42 : use m_mkffnl
43 : use m_xcdata
44 :
45 : use m_gwdefs, only : czero_gw
46 : use defs_abitypes, only : mpi_type
47 : use defs_datatypes, only : pseudopotential_type
48 : use m_array, only : array2_gwpc_t
49 : use m_dtfil, only : datafiles_type
50 : use m_dtset, only : dataset_type
51 : use m_gwdefs, only : GW_Q0_DEFAULT
52 : use m_time, only : cwtime, cwtime_report, timab, sec2str, timab
53 : use m_fstrings, only : itoa, ftoa, sjoin, ktoa, ltoa, strcat
54 : use m_numeric_tools, only : arth, c2r, r2c, get_diag, linfit, iseven, simpson_cplx, print_arr, inrange !, bins_t
55 : use m_io_tools, only : iomode_from_fname
56 : use m_fftcore, only : ngfft_seq, sphereboundary, print_ngfft
57 : use m_fft_mesh, only : setmesh
58 : use m_cgtk, only : cgtk_rotate, cgtk_change_gsphere
59 : use m_crystal, only : crystal_t
60 : use m_kpts, only : kpts_ibz_from_kptrlatt, kpts_timrev_from_kptopt, kpts_map
61 : use m_kg, only : getph, mkkin
62 : use m_bz_mesh, only : isamek, kmesh_t
63 : use m_gsphere, only : gsphere_t
64 : use m_getgh1c, only : getgh1c, rf_transgrid_and_pack
65 : use m_ioarr, only : read_rhor
66 : use m_hdr, only : hdr_type
67 : use m_dvdb, only : dvdb_t
68 : use m_ifc, only : ifc_type
69 : use m_wfd, only : wfd_t
70 : use m_pawang, only : pawang_type
71 : use m_pawrad, only : pawrad_type
72 : use m_pawtab, only : pawtab_type
73 : use m_pawrhoij, only : pawrhoij_type
74 : use m_pawfgr, only : pawfgr_type
75 : use m_dfpt_cgwf, only : stern_t
76 : use m_io_screening, only : hscr_t, get_hscr_qmesh_gsph, read_screening
77 : use m_vcoul, only : vcoul_t
78 : use m_gstore, only : gstore_t, gqk_t, gstore_check_restart, &
79 : GSTORE_KQ_MISSING, GSTORE_KQ_COMPUTED, GSTORE_GTYPE_GWPT, gstore_symmetrize
80 : use m_rhotoxc, only : rhotoxc
81 : use m_drivexc, only : check_kxc
82 : use m_occ, only : get_fact_spin_tol_empty
83 : use m_ebands, only : ebands_t
84 : use m_pstat, only : pstat_proc
85 : use m_screen, only : em1_symmetrize_op
86 : use m_ppmodel, only : ppmodel_t
87 : use m_lgroup, only : lgroup_t
88 : use m_initylmg, only : initylmg_k
89 :
90 : implicit none
91 :
92 : private
93 : !!***
94 :
95 : #ifdef HAVE_MPI1
96 : include 'mpif.h'
97 : #endif
98 :
99 : public :: gwpt_run ! Main entry point to compute GWPT e-ph matrix elements
100 :
101 : !----------------------------------------------------------------------
102 :
103 : contains !=====================================================
104 : !!***
105 :
106 : !----------------------------------------------------------------------
107 :
108 : !!****f* m_gwpt/gwpt_run
109 : !! NAME
110 : !! gwpt_run
111 : !!
112 : !! FUNCTION
113 : !! Compute e-ph matrix elements with the GWPT formalism.
114 : !!
115 : !! INPUTS
116 : !! wfk0_path=String with the path to the GS unperturbed WFK file.
117 : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
118 : !! dtset<dataset_type>=All input variables for this dataset.
119 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
120 : !! dvdb=Database with the DFPT SCF potentials.
121 : !! drhovdb=Database with the DFPT SCF densities.
122 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
123 : !! wfk_hdr=Header of the WFK file.
124 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
125 : !! pawang<pawang_type)>=PAW angular mesh and related data.
126 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
127 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
128 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
129 : !! comm=MPI communicator.
130 : !!
131 : !! NOTES
132 : !!
133 : !! 1) Conventions used for g-spheres and the periodic part of the KS states:
134 : !!
135 : !! _kq --> k + q
136 : !! _kmp --> k - p
137 : !! _kqmp --> k + q - p
138 : !!
139 : !! 2) The routines used to symmetrize wavefunctions and DFPT scattering potentials
140 : !! expect in input symmetry tables generated using different conventions.
141 : !! For the wavefunctions, we use the symrel convention while for the scattering potentials we use the symrec convention.
142 : !! We encode this in the name of the variable using e.g. mapc_qq for the symrec convention (C) and mapl_k convention (L)
143 : !!
144 : !! 3) The DFPT routines operate on double-precision wavefunctions stored in arrays with real/imag part e.g. cg(1:2,npw_k)
145 : !! while the GW routines operate on complex arrays of kind=gwp where gwp is defined at configure-time.
146 : !! The default value of gwp is single-precision.
147 : !! We use the following conventions for the buffers used to store the wavefunctions:
148 : !!
149 : !! cg_kq, cr_kq
150 : !! cg1_kqmp, cr1_kqmp
151 : !!
152 : !! OUTPUT
153 : !! GSTORE.nc file.
154 : !!
155 : !! TODO
156 : !! Problems to be addressed:
157 : !!
158 : !! 1)
159 : !! Sigma is usually split into Sigma_c(w) and Sigma_x where Sigma_x is the static Fock operator
160 : !! evaluated with KS orbitals. The advantage of such partitioning is that Sigma_x = iGv
161 : !! can be computed by summing over occupied states only. Sigma_x requires more G-vectors to converge
162 : !! as the bare Coulomb interaction goes as 1/|q+G|^2 that is not integrable in 3D but this "expensive"
163 : !! operations are needed only inside a sum over bands that is restricted to occupied states.
164 : !! On the other hand, Sigma_c(w) is way more expensive as we have to sum a large number of empty states
165 : !! while taking the w-dependence of the screening into account.
166 : !! Fortunately, all the operations can be restricted to a small G-sphere of kinetic energy ecuteps that can be handled
167 : !! with a coarser FFT mesh.
168 : !! Another distinct advantage of such splitting is that one can handle the divergence in v(q,G) for |q+G| --> 0
169 : !! using well know techniques from GW and the anisotropic behavior of eps-1(q) for q --> 0 in low-dimensional systems.
170 : !! The disavantage is that one needs to compute the GWPT e-ph matrix in two steps, first Sigma_c and then Sigma_x,
171 : !! so certain operations such as the k-point mapping, and the computation of the form factors are performed twice
172 : !! Note, however, that MG believes that Sigma_x is a much better approximation than v_xc when one is interested
173 : !! in the e-ph matrix elements connecting low-energy states such as band edges to high-energy states.
174 : !!
175 : !! 2)
176 : !! We need to solve the NSCF Sternheimer for q and -q. In principle one can solve the equation only at q
177 : !! and then use spatial inversion or TR to get the solution at -q but this requires solving the Sternheimer
178 : !! for all the pp wavevectors in the BZ (or better in the IBZ_{q,k,alpha}. The use of symmetries is rendered complicated
179 : !! by the parallelism over pp but perhaps one can precompute \Delta psi with all MPI procs and write the results to temporary files.
180 : !!
181 : !! 3) Check gstore_brange less than nband
182 : !!
183 : !! Debugging options:
184 : !! useria = 1 # write gvals=gxc to GSTORE.nc
185 : !! useria = 0 (default) # write gvals=gks-gxc+gsigx+gsigc to GSTORE.nc
186 : !! useria = -1 # write gvals=gsigx (when userid=0) + gsigc (when useric=0) to GSTORE.nc
187 : !! useria = 888 # use xmpi_comm_self when reopening GSTORE.nc to avoid parallel NetCDF deadlocks
188 : !!
189 : !! userib = 0 (default) # calculate all k and q
190 : !! userib = 1 # filter k and q, only calculate k=Lambda, and q=L (for comparesion with finite difference GW)
191 : !!
192 : !! useric = 0 (default) # include the correlation part of Sigma
193 : !! useric = 1 # exclude the correlation part of Sigma
194 : !!
195 : !! userid = 0 (default) # include the exchange part of Sigma
196 : !! userid = 1 # exclude the exchange part of Sigma
197 : !!
198 : !! userie = 1 : Exclude pp = Gamma
199 : !! userie = 0 : Use all pp of the given p-grid
200 : !! userie = -1 : Use only pp = Gamma
201 : !!
202 : !! SOURCE
203 :
204 3 : subroutine gwpt_run(wfk0_path, dtfil, ngfft, ngfftf, dtset, cryst, ebands, dvdb, drhodb, ifc, wfk_hdr, &
205 3 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
206 :
207 : !Arguments ------------------------------------
208 : !scalars
209 : character(len=*),intent(in) :: wfk0_path
210 : integer,intent(in) :: comm
211 : type(datafiles_type),intent(in) :: dtfil
212 : type(dataset_type),intent(inout) :: dtset
213 : type(crystal_t),intent(in) :: cryst
214 : type(ebands_t),target,intent(in) :: ebands
215 : type(dvdb_t),intent(inout) :: dvdb, drhodb
216 : type(pawang_type),intent(in) :: pawang
217 : type(pseudopotential_type),intent(in) :: psps
218 : type(pawfgr_type),intent(in) :: pawfgr
219 : type(ifc_type),intent(in) :: ifc
220 : type(hdr_type),intent(in) :: wfk_hdr
221 : type(mpi_type),intent(inout) :: mpi_enreg
222 : !arrays
223 : integer,intent(in) :: ngfft(18), ngfftf(18)
224 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
225 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
226 :
227 : !Local variables ------------------------------
228 : !scalars
229 : integer,parameter :: LOG_MODQ = 1, LOG_MODK = 4, LOG_MODP = 4, ENOUGH_STERN = 5
230 : integer,parameter :: tim_getgh1c1 = 1, berryopt0 = 0, ider0 = 0, idir0 = 0, istwfk1 = 1, cplex1 = 1, pawread0 = 0
231 : integer,parameter :: master = 0, ndat1 = 1, with_cplex0 = 0, n3xccc0 = 0, optder0 = 0
232 : integer :: band, band_me, nband_me, stern_comm, nkpt, my_rank, nsppol, iq_ibz, iq_bz, my_npert
233 : integer :: nb_k, nb_kq, bstart_k, bstop_k, bstart_kq, bstop_kq, matblk, method, enforce_sym
234 : integer :: cplex,drho_cplex,nkxc,nk3xc,option,usexcnhat,db_iqpt,natom,natom3,ipc,nspinor,nproc
235 : integer :: ib_sum, ii, u1_band !,u1c_ib_k, jj, iw !ib_kq, band_ks, ib_k, ibsum_kq, u1_master, ip
236 : integer :: my_is, spin, idir,ipert, ig, max_npw_xc, min_npw_xc, npw_x, npw_c, nw_nk, nw_mkq
237 : integer :: isym_q, trev_q, ip_ibz
238 : integer :: ik_ibz, isym_k, trev_k, npw_k, istwf_k, npw_k_ibz, istwf_k_ibz, ik_glob, ik_bz
239 : integer :: ikq_ibz, isym_kq, trev_kq, npw_kq, istwf_kq, npw_kq_ibz, istwf_kq_ibz
240 : integer :: ikmp_ibz, isym_kmp, trev_kmp, npw_kmp, istwf_kmp
241 : integer :: ikqmp_ibz, isym_kqmp, trev_kqmp, npw_kqmp, istwf_kqmp, npw_kqmp_ibz, istwf_kqmp_ibz, mpw,ierr,nqbz,ncerr !,spad
242 : integer :: n1,n2,n3,n4,n5,n6,nspden, mqmem, im_kq, m_kq, in_k, n_k, restart, root_ncid, spin_ncid, usecprj
243 : integer :: nfft,nfftf,mgfft,mgfftf,nkpg_k,nkpg_kq,nkpg_kqmp,nkpg_kmp,imyp, cnt, nvloc, iw_nk, iw_mkq, ndone, nmiss
244 : integer :: my_ipp, ipp_bz, ipp_ibz, isym_pp, itim_pp, comm_rpt, nqlwl, scr_iomode, stern_qq_ierr, stern_mq_ierr
245 : integer :: qptopt, my_iq, my_ik, qbuf_size, iqbuf_cnt, timrev_k, timrev_q
246 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_qq, wall_qq, gflops_qq, cpu_kk, wall_kk, gflops_kk, cpu_pp, wall_pp, gflops_pp
247 : real(dp) :: drude_plsmf, my_plsmf !, rtmp
248 : real(dp) :: fact_spin, theta_mu_minus_e0i, tol_empty, tol_empty_in, e_nk, e_mkq, ebsum_kmp, ebsum_kqmp ! e0i
249 : real(dp) :: eta_sum, delta_E, inv_delta_E
250 3 : real(dp),contiguous, pointer :: qp_ene(:,:,:), qp_occ(:,:,:)
251 : real(dp) :: weight_q,bigexc,bigsxc,vxcavg ! ediff, eshift, q0rad, bz_vol
252 : logical :: isirr_k, isirr_kq, isirr_kmp, isirr_kqmp, qq_is_gamma, pp_is_gamma, isirr_q
253 : logical :: stern_use_cache, need_ftinterp, symmetrize, use_lgk
254 : logical :: print_time_qq, print_time_kk, print_time_pp, non_magnetic_xc, need_x_kmp, need_x_kqmp, test_sigma
255 : complex(dp) :: ieta !, idelta_sum
256 3 : type(wfd_t) :: wfd
257 3 : type(gs_hamiltonian_type) :: gs_ham_kqmp, gs_ham_kmp
258 3 : type(rf_hamiltonian_type) :: rf_ham_kqmp, rf_ham_kmp
259 312 : type(crystal_t) :: pot_cryst, den_cryst
260 6 : type(hdr_type) :: pot_hdr, den_hdr
261 3 : type(stern_t) :: stern_kmp, stern_kqmp
262 78 : type(kmesh_t) :: pp_mesh, kmesh
263 3 : type(gsphere_t),target :: gsph_x, gsph_c
264 12 : type(hscr_t),target :: hscr
265 66 : type(vcoul_t) :: vcp
266 12 : type(gstore_t),target :: gstore
267 : type(gqk_t),pointer :: gqk
268 3 : type(lgroup_t) :: lg_myq
269 : type(xcdata_type) :: xcdata
270 3 : type(ppmodel_t) :: ppm
271 : character(len=fnlen) :: screen_filepath, gstore_filepath
272 : character(len=5000) :: msg, qq_bz_string, kk_string, qkp_string, pp_string
273 : character(len=500) :: init_mode
274 : !arrays
275 6 : integer :: nbsum, my_bsum_start(dtset%nsppol), my_bsum_stop(dtset%nsppol), my_nbsum(dtset%nsppol)
276 : integer :: g0_k(3), g0_q(3), g0_kq(3), g0_kmp(3), g0_kqmp(3), units(2), work_ngfft(18), gmax(3)
277 : integer :: mapl_k(6), mapl_kq(6), mapl_kqmp(6), mapl_kmp(6), mapc_qq(6), mapc_qq2dvdb(6)
278 9 : integer :: my_pp_start_spin(dtset%nsppol), my_pp_stop_spin(dtset%nsppol), my_npp(dtset%nsppol)
279 3 : integer(i1b),allocatable :: itreatq_dvdb(:)
280 6 : integer,allocatable :: kg_k(:,:), kg_kq(:,:), kg_kmp(:,:), kg_kqmp(:,:), my_pp_inds(:)
281 3 : integer,allocatable :: gbound_k(:,:), gbound_kq(:,:), gbound_kmp(:,:), gbound_kqmp(:,:), gbound_c(:,:), gbound_x(:,:)
282 9 : integer,allocatable :: nband(:,:), wfd_istwfk(:), qibz2dvdb(:) ! count_bk(:,:),
283 6 : integer,allocatable :: iq_buf(:,:), done_qbz_spin(:,:), state_kq(:,:)
284 3 : integer(i1b),allocatable :: itreat_qibz(:)
285 3 : integer, contiguous, pointer :: kg_c(:,:), kg_x(:,:)
286 : !real(dp) :: eig0nk !, cpu, wall, gflops !, cpu_q, wall_q, gflops_q, cpu_all, wall_all, gflops_all
287 : complex(gwp) :: ctmp_gwpc, xdot_tmp
288 : complex(dp) :: ctmp_dp
289 : !type(bins_t) :: bins
290 : !arrays
291 6 : real(dp) :: fermie1_idir_ipert(3,cryst%natom), ylmgr_dum(1,1,1), dum_nhat(0), dum_xccc3d(0), tsec(2)
292 : real(dp) :: kk(3),kq(3),kk_ibz(3),kq_ibz(3), kqmp(3), kmp(3), pp(3), kmp_ibz(3), kqmp_ibz(3), qq_ibz(3), qq_bz(3)
293 6 : real(dp),allocatable :: qlwl(:,:), kpg_k(:,:),kpg_kq(:,:),kpg_kmp(:,:),kpg_kqmp(:,:)
294 9 : real(dp),allocatable :: ffnl_kmp(:,:,:,:),ffnl_kqmp(:,:,:,:), kinpw_kqmp(:), kinpw_kmp(:), ph3d_kqmp(:,:,:), ph3d_kmp(:,:,:)
295 3 : real(dp),allocatable, target :: vxc1_qq(:,:,:,:)
296 3 : real(dp),allocatable :: gxc_atm(:,:,:,:), gks_atm(:,:,:,:), gks_atm2(:,:,:,:), gsig_atm(:,:,:,:)
297 9 : real(dp),allocatable :: cg_work(:,:), ug_k(:,:), ug_kq(:,:), ph1d(:,:)
298 6 : real(dp),allocatable :: vlocal(:,:,:,:), vlocal1_qq(:,:,:,:,:), v1scf_qq(:,:,:,:), vlocal1_mqq(:,:,:,:,:), v1scf_mq(:,:,:,:)
299 9 : real(dp),allocatable :: ylm_kmp(:,:), ylm_kqmp(:,:), vtrial(:,:), work(:,:,:,:), rhor(:,:), vxc(:,:), kxc(:,:)
300 3 : real(dp),allocatable :: omegame0i_nk(:), omegame0i_mkq(:), omegas_nk(:), omegas_mkq(:)
301 3 : real(dp),allocatable :: my_gbuf(:,:,:,:,:,:), my_gbuf_ks(:,:,:,:,:,:)
302 9 : real(dp),allocatable :: cg_kmp(:,:), cg_kqmp(:,:), cg1_kqmp(:,:), cg1_kmp(:,:), full_cg1_kqmp(:,:), full_cg1_kmp(:,:), vxc_nk(:,:)
303 3 : complex(dp), contiguous, pointer :: cvxc1_qq_ptr(:,:,:)
304 6 : complex(gwp),allocatable :: ur_star_kmp(:), ur_star_kqmp(:), cwork_ur(:), rhotwg_c(:), rhotwg_x(:), vc_sqrt_gx(:)
305 6 : complex(gwp),allocatable :: full_ur1_kqmp(:), full_ur1_star_kmp(:), sigcme_nk(:), sigcme_mkq(:), ur_nk(:,:), ur_mkq(:,:)
306 3 : complex(gwp),allocatable :: vec_gwc_nk(:,:,:), vec_gwc_mkq(:,:,:), vec_gx_nk(:,:), vec_gx_mkq(:,:)
307 3 : complex(gwp),allocatable :: vec_coh_nk(:,:), vec_coh_mkq(:,:), wc0_pbz(:,:)
308 9 : complex(gwp),allocatable :: botsq_pbz(:,:), otq_pbz(:,:), dmeig_pbz(:,:), epsm1_ggw(:,:,:)
309 3 : complex(gwp),allocatable :: trans_botsq_pbz(:,:), trans_otq_pbz(:,:), trans_dmeig_pbz(:,:), sigx_nk(:,:), sigce0_nk(:,:)
310 6 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
311 3 : type(pawcprj_type),allocatable :: cwaveprj0(:,:), cwaveprj(:,:)
312 6 : type(pawrhoij_type),allocatable :: pot_pawrhoij(:), den_pawrhoij(:)
313 3 : type(lgroup_t),allocatable :: lg_myk(:)
314 3 : type(array2_gwpc_t),allocatable :: wc0_pibz(:)
315 : !************************************************************************
316 :
317 3 : if (psps%usepaw == 1) then
318 0 : ABI_ERROR("PAW not implemented")
319 : ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
320 : end if
321 :
322 3 : ABI_CHECK_IEQ(dtset%nspinor, 1, "GWPT with nspinor 2 not coded")
323 3 : ABI_CHECK_IEQ(dtset%nsppol, 1, "GWPT with nsppol 2 not tested")
324 :
325 3 : test_sigma = .False.
326 : !test_sigma = my_rank == 2
327 :
328 : !dtset%useric = 1 ! exclude the correlation part of Sigma
329 : !dtset%userid = 1 ! exclude the exchange part of Sigma
330 :
331 3 : call cwtime(cpu_all, wall_all, gflops_all, "start")
332 :
333 : ! Copy important dimensions
334 9 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm); units = [std_out, ab_out]
335 3 : natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor; nspden = dtset%nspden
336 3 : ieta = +j_dpc * dtset%zcut
337 : !idelta_sum = +j_dpc * dtset%elph2_imagden
338 3 : eta_sum = dtset%elph2_imagden
339 :
340 : ! Set tolerance used to decide if a band is empty
341 : ! and normalization of theta_mu_minus_esum. If nsppol == 2, qp_occ $\in [0,1]$
342 3 : tol_empty_in = 0.01
343 3 : call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
344 :
345 3 : qp_ene => ebands%eig; qp_occ => ebands%occ
346 :
347 : ! Check if a previous GSTORE.nc file is present to restart the calculation if dtset%eph_restart == 1,
348 : ! and use done_qbz_spin mask to cycle the loops below if restart /= 0.
349 3 : gstore_filepath = strcat(dtfil%filnam_ds(4), "_GSTORE.nc")
350 :
351 : !gstore_filepath = dtfil%filgstorein
352 3 : call gstore_check_restart(gstore_filepath, dtset, nqbz, done_qbz_spin, restart, comm)
353 :
354 3 : if (restart == 0) then
355 : ! Build new gstore object from dtset input variables.
356 3 : call gstore%init(gstore_filepath, dtset, dtfil, wfk_hdr, cryst, ebands, ifc, comm, gtype=GSTORE_GTYPE_GWPT)
357 12 : ABI_REMALLOC(done_qbz_spin, (gstore%nqbz, nsppol))
358 86 : done_qbz_spin = 0
359 : else
360 : ! Init gstore from pre-existent file. gstore_gname and read_dw are not relevant here.
361 : call gstore%from_ncpath(gstore_filepath, with_cplex0, dtset, dtfil, cryst, ebands, ifc, &
362 0 : "atom", dtset%gstore_gname, .False., comm)
363 : end if
364 :
365 3 : if (nqbz /= 0) then
366 0 : ABI_CHECK_IEQ(nqbz, gstore%nqbz, "nqbz from file differs from gstore%nqbz")
367 : else
368 3 : nqbz = gstore%nqbz
369 : end if
370 :
371 : ! Open GSTORE.nc file and go to data mode.
372 : ! On some systems, parallel NetCDF access can deadlock when restarting GWPT.
373 : ! useria == 888 activates independent I/O as a temporary portability workaround.
374 3 : if (dtset%useria == 888) then
375 0 : NCF_CHECK(nctk_open_modify(root_ncid, gstore%path, xmpi_comm_self))
376 : else
377 3 : NCF_CHECK(nctk_open_modify(root_ncid, gstore%path, comm))
378 : end if
379 3 : NCF_CHECK(nctk_set_datamode(root_ncid))
380 :
381 3 : call gstore%get_missing_qbz_spin(done_qbz_spin, ndone, nmiss)
382 : !call wrtout(units, sjoin("- Number of q-points/spin completed:", itoa(count(done_qbz_spin == 1)))
383 :
384 : ! TODO: Check restart capabilities. It seems they don't work properly.
385 : ! NB: Write phonon data here as we are not guaranteed to have all the IBZ q-points
386 : ! inside the loop over my_iq if filtering has been used.
387 : ! Make sure internal table with gstore_done_qbz_spin is properly filled.
388 3 : if (ndone == 0) then
389 3 : call gstore%compute_and_write_ph(root_ncid)
390 : else
391 0 : call wrtout(units, sjoin("- Restarting GSTORE calculation from file:", gstore%path))
392 0 : call wrtout(units, sjoin(" Found: ", itoa(ndone), " (qpt, spin) entries already computed"))
393 : end if
394 :
395 : ! ================
396 : ! HANDLE SCREENING
397 : ! ================
398 : ! Read g-sphere and pp_mesh from SCR file and init gsph_c for the correlated part.
399 3 : screen_filepath = dtfil%fnameabi_scr
400 3 : ABI_CHECK(dtfil%fnameabi_scr /= ABI_NOFILE, "SCR file must be specified")
401 3 : call get_hscr_qmesh_gsph(screen_filepath , dtset, cryst, hscr, pp_mesh, gsph_c, qlwl, comm)
402 3 : call hscr%print(units, dtset%prtvol, header="Header of the SCR file")
403 :
404 3 : nqlwl = size(qlwl, dim=2)
405 3 : if (nqlwl == 0) then
406 0 : nqlwl=1
407 0 : ABI_MALLOC(qlwl,(3,nqlwl))
408 0 : qlwl(:,nqlwl)= GW_Q0_DEFAULT
409 : write(msg,'(3a,i0,a,3f9.6)')&
410 0 : "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,&
411 0 : "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1)
412 0 : ABI_COMMENT(msg)
413 : end if
414 :
415 : ! Init g-sphere for the exchange part from ecutsigx.
416 3 : call gsph_c%extend(cryst, dtset%ecutsigx, gsph_x)
417 :
418 3 : call gsph_c%print(units, dtset%prtvol, header="Header of the correlation g-sphere")
419 3 : call gsph_x%print(units, dtset%prtvol, header="Header of the exchange g-sphere")
420 :
421 : ! TODO:
422 : ! Here we sort the pp_mesh by stars so that we can split the pp wavevectors in blocks and therefore
423 : ! reduce the number of wavevectors in the IBZ that must be stored in memory.
424 : !call pp_mesh%pack_by_stars()
425 :
426 : ! Distribute the sum over pp wavevectors inside pp_sum_comm using block distribution.
427 9 : my_pp_start_spin = -1; my_pp_stop_spin = 0
428 6 : do my_is=1,gstore%my_nspins
429 3 : spin = gstore%my_spins(my_is); gqk => gstore%gqk(my_is)
430 3 : call xmpi_split_block(pp_mesh%nbz, gqk%pp_sum_comm%value, my_npp(spin), my_pp_inds)
431 3 : if (my_npp(spin) > 0) then
432 3 : my_pp_start_spin(spin) = my_pp_inds(1); my_pp_stop_spin(spin) = my_pp_inds(my_npp(spin))
433 : end if
434 6 : ABI_SFREE(my_pp_inds)
435 : end do ! my_is
436 :
437 : ! TODO: Use same FFT mesh as GW code but make sure that all potentials are properly interpolated.
438 :
439 : ! === Setup of the FFT mesh for the oscillator strengths ===
440 : ! * Init gwc_ngfft(7:18) and gwx_ngfft(7:18) with Dtset%ngfft(7:18)
441 : ! * Here we redefine gwc_ngfft(1:6) according to the following options:
442 : !
443 : ! method == 0 --> FFT grid read from fft.in (debugging purpose)
444 : ! method == 1 --> Normal FFT mesh
445 : ! method == 2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
446 : ! method == 3 --> Doubled FFT grid, same as the the FFT for the density,
447 : !
448 : ! enforce_sym == 1 --> Enforce a FFT mesh compatible with all the symmetry operation and FFT library
449 : ! enforce_sym == 0 --> Find the smallest FFT grid compatible with the library, do not care about symmetries
450 : !
451 : !gwc_ngfft(1:18) = dtset%ngfft(1:18)
452 : !gwx_ngfft(1:18) = dtset%ngfft(1:18)
453 :
454 3 : method = 2
455 3 : if (dtset%fftgw == 00 .or. dtset%fftgw == 01) method = 0
456 : if (dtset%fftgw == 10 .or. dtset%fftgw == 11) method = 1
457 : if (dtset%fftgw == 20 .or. dtset%fftgw == 21) method = 2
458 3 : if (dtset%fftgw == 30 .or. dtset%fftgw == 31) method = 3
459 3 : enforce_sym = mod(dtset%fftgw, 10)
460 :
461 : #if 0
462 : ! FFT mesh for sigma_x.
463 : call setmesh(cryst%gmet, Gsph_Max%gvec, gwx_ngfft, Sigp%npwvec, Sigp%npwx, Sigp%npwwfn, &
464 : gwx_nfftot, method, Sigp%mG0, cryst, enforce_sym)
465 :
466 : ! FFT mesh for sigma_c.
467 : call setmesh(cryst%gmet, Gsph_Max%gvec, gwc_ngfft, Sigp%npwvec, epsm1%npwe, Sigp%npwwfn,&
468 : gwc_nfftot, method, Sigp%mG0, cryst, enforce_sym, unit=dev_null)
469 : #endif
470 :
471 : ! Initialize Coulomb term on the IBZ of the pp_mesh. Use largest G-sphere.
472 3 : call kmesh%init(cryst, wfk_hdr%nkpt, wfk_hdr%kptns, dtset%kptopt)
473 :
474 : ! Some required information are not filled up inside kmesh_init
475 : ! So doing it here, even though it is not clean.
476 39 : Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:)
477 3 : Kmesh%nshift =Dtset%nshiftk
478 9 : ABI_MALLOC(Kmesh%shift,(3,Kmesh%nshift))
479 15 : Kmesh%shift(:,:) =Dtset%shiftk(:,1:Dtset%nshiftk)
480 :
481 3 : npw_x = gsph_x%ng; npw_c = gsph_c%ng
482 3 : min_npw_xc = min(npw_x, npw_c)
483 3 : max_npw_xc = max(npw_x, npw_c)
484 3 : if (gsph_x%ng >= gsph_c%ng) then
485 : call vcp%init(gsph_x, cryst, pp_mesh, kmesh, dtset%gw_rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, gsph_x%ng, &
486 3 : nqlwl, qlwl, comm)
487 : else
488 : call vcp%init(gsph_c, cryst, pp_mesh, kmesh, dtset%gw_rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, gsph_c%ng, &
489 0 : nqlwl, qlwl, comm)
490 : end if
491 3 : ABI_FREE(qlwl)
492 :
493 3 : if (my_rank == master) then
494 3 : call kmesh%print(units, header="K-mesh for wavefunctions", prtvol=dtset%prtvol)
495 3 : call pp_mesh%print(units, header="P-mesh for Screening", prtvol=dtset%prtvol)
496 3 : call gsph_x%print(units, dtset%prtvol, header="G-sphere for exchange")
497 3 : call gsph_c%print(units, dtset%prtvol, header="G-sphere for correlation")
498 3 : call vcp%print(units, prtvol=dtset%prtvol)
499 : end if
500 3 : call kmesh%free()
501 :
502 3 : ABI_CHECK_IGE(npw_x, 1, "npw_x <= 1")
503 3 : ABI_CHECK_IGE(npw_c, 1, "npw_c <= 1")
504 3 : ABI_CHECK_IGEQ(npw_x, npw_c, "npw_x < npw_c is not supported")
505 :
506 : ! Initialize the wave function descriptor.
507 : ! Each node has all k-points and spins and bands between my_bsum_start and my_bsum_stop
508 : ! TODO: One can exploit qq, kk and pp parallelism to find the wavevectors in the IBZ
509 : ! that will be needed in the loops below and allocate only these wavevectors so that memory scales.
510 :
511 12 : nbsum = dtset%mband; my_bsum_start = 1; my_bsum_stop = nbsum; my_nbsum = my_bsum_stop - my_bsum_start + 1
512 : ! FIXME: This term is needed in metals.
513 27 : fermie1_idir_ipert = zero
514 :
515 3 : nkpt = wfk_hdr%nkpt
516 12 : ABI_MALLOC(nband, (nkpt, nsppol))
517 15 : ABI_MALLOC(bks_mask, (dtset%mband, nkpt, nsppol))
518 12 : ABI_MALLOC(keep_ur, (dtset%mband, nkpt, nsppol))
519 493 : nband = dtset%mband; bks_mask = .False.; keep_ur = .False.
520 :
521 : ! Distribute wavefunctions according to the set of kk, qq and pp wavevectors treated by this MPI proc.
522 : ! Also, compute mpw and gmax including the additional pp
523 : ! This is the maximum number of PWs for all possible k+q treated.
524 : call gstore%fill_bks_mask_pp_mesh(dtset%ecut, dtset%mband, nkpt, nsppol, my_pp_start_spin, my_pp_stop_spin, pp_mesh, &
525 3 : my_bsum_start, my_bsum_stop, bks_mask, mpw, gmax)
526 :
527 : ! Init work_ngfft
528 12 : gmax = gmax + 4 ! FIXME: this is to account for umklapp, should also consider Gamma-only and istwfk
529 12 : gmax = 2*gmax + 1
530 : !gmax = 12*gmax + 1
531 :
532 3 : call ngfft_seq(work_ngfft, gmax)
533 : !print *, "gmax:", gmax
534 : !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
535 15 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
536 :
537 : !if (dtset%userie == 124) then
538 : ! ! Debugging section have all states on each MPI rank.
539 : ! bks_mask = .True.; call wrtout(std_out, " Storing all bands for debugging purposes.")
540 : !end if
541 :
542 : ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts)
543 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
544 9 : ABI_MALLOC(wfd_istwfk, (nkpt))
545 22 : wfd_istwfk = 1
546 :
547 : call wfd%init(cryst, pawtab, psps, keep_ur, dtset%mband, nband, nkpt, nsppol, bks_mask,&
548 : nspden, nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,&
549 3 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
550 :
551 3 : ABI_FREE(nband)
552 3 : ABI_FREE(bks_mask)
553 3 : ABI_FREE(keep_ur)
554 3 : ABI_FREE(wfd_istwfk)
555 3 : call pstat_proc%print(_PSTAT_ARGS_)
556 :
557 : ! Read wavefunctions.
558 3 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
559 :
560 : ! FFT meshes from input file, not necessarily equal to the ones found in the external files.
561 : ! NB: ur arrays are always allocated with nfft and not with product(ngfft(4:6)).
562 12 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
563 21 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
564 3 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6)
565 :
566 6 : call print_ngfft([std_out], ngfft, header="FFT mesh")
567 :
568 : ! Set the FFT mesh
569 3 : call wfd%change_ngfft(cryst, psps, ngfft)
570 3 : call wfd%print(units, header="Wavefunctions for GWPT calculation.")
571 :
572 3 : usecprj = dtset%usepaw
573 12 : ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
574 12 : ABI_MALLOC(cwaveprj, (natom, nspinor*usecprj))
575 12 : ABI_MALLOC(gbound_k, (2*mgfft+8, 2))
576 6 : ABI_MALLOC(gbound_kq, (2*mgfft+8, 2))
577 6 : ABI_MALLOC(gbound_kmp, (2*mgfft+8, 2))
578 6 : ABI_MALLOC(gbound_kqmp, (2*mgfft+8, 2))
579 6 : ABI_MALLOC(gbound_c, (2*mgfft+8, 2))
580 6 : ABI_MALLOC(gbound_x, (2*mgfft+8, 2))
581 9 : ABI_MALLOC(cg_work, (2, mpw*nspinor))
582 9 : ABI_MALLOC(full_ur1_kqmp, (nfft*nspinor))
583 6 : ABI_MALLOC(full_ur1_star_kmp, (nfft*nspinor))
584 :
585 : ! ============================
586 : ! Compute v_nk matrix elements
587 : ! ============================
588 : ! Diagonal elements of velocity operator in cartesian coordinates for all kk in the IBZ.
589 : ! Use ndone to understand if velocities have been already compured in a previous run.
590 :
591 : ! TODO: A similar section of code is found in gstore%compute. Should have common routine.
592 :
593 3 : if (gstore%with_vk /= 0 .and. ndone == 0) then
594 3 : call gstore%compute_and_write_vk(mpw, wfd, ebands, psps, pawtab, root_ncid)
595 : end if
596 :
597 : ! Radius of sphere with volume equivalent to the micro zone.
598 : !q0rad = two_pi * (three / (four_pi * cryst%ucvol * gstore%nqbz)) ** third
599 : !bz_vol = two_pi**3 / cryst%ucvol
600 :
601 : ! Open the DVDB file with first-order potentials and drhodb with the first-order densities.
602 3 : call dvdb%open_read(ngfftf, xmpi_comm_self)
603 3 : call drhodb%open_read(ngfftf, xmpi_comm_self)
604 3 : ABI_CHECK(dvdb%has_fields("pot1", msg), msg)
605 3 : ABI_CHECK(drhodb%has_fields("den1", msg), msg)
606 :
607 : ! Make sure that dvdb and drhodb have the same q-points.
608 3 : ABI_CHECK_IEQ(dvdb%nqpt, drhodb%nqpt, "Different number of q-points in DVDB and DRHODB")
609 3 : ierr = 0
610 22 : do ii=1,dvdb%nqpt
611 79 : if (any(abs(dvdb%qpts(:, ii) - drhodb%qpts(:, ii)) > tol12)) then
612 0 : ierr = ierr + 1; call wrtout(std_out, sjoin(ktoa(dvdb%qpts(:, ii)), " /= ", ktoa(drhodb%qpts(:, ii))))
613 : end if
614 : end do
615 3 : ABI_CHECK(ierr == 0, "Found different q-points in DVDB and DRHODB. See messages above!")
616 :
617 : ! Check if the q-points are present in the DVDB
618 : ! qmap_symrec gives the mapping gstore%ibz --> dvdb%ibz
619 : !call dvdb%need_ftinterp(gstore%nqibz, gstore%qibz, gstore%qptopt, qmap_symrec, need_ftinterp)
620 : !ABI_FREE(qmap_symrec)
621 : !need_ftinterp = .True.
622 :
623 : !if (.not. need_ftinterp .and. dtset%eph_need_ftinterp /= 0) then
624 : ! ABI_WARNING("Enforcing FT interpolation for q-points even if it's not strictly needed.")
625 : ! need_ftinterp = .True.
626 : !end if
627 :
628 : !if (need_ftinterp) then
629 : ! call wrtout(units, " Cannot find all q-points in the DVDB --> Activating Fourier interpolation.")
630 : ! call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt, 1, dtset%ddb_shiftq, nfftf, ngfftf, xmpi_comm_self)
631 : !else
632 : ! call wrtout(units, " DVDB file contains all q-points along the path --> Reading DFPT potentials from file.")
633 : !end if
634 :
635 : ! Activate parallelism over perturbations at the level of the DVDB, my_npert is output
636 3 : call gstore%set_perts_distrib(cryst, dvdb, my_npert)
637 3 : call gstore%set_perts_distrib(cryst, drhodb, my_npert)
638 : !print *, "Treating my_npert", my_npert
639 :
640 : ! This part is taken from dfpt_vtorho
641 : !==== Initialize most of the Hamiltonian (and derivative) ====
642 : ! 1) Allocate all arrays and initialize quantities that do not depend on k and spin.
643 : ! 2) Perform the setup needed for the non-local factors:
644 : !
645 : ! Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_ham_kqmp.
646 : ! PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
647 :
648 : ! Get one-dimensional structure factor information on the coarse grid.
649 9 : ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
650 3 : call getph(cryst%atindx, natom, n1, n2, n3, ph1d, cryst%xred)
651 :
652 : call gs_ham_kqmp%init(psps, pawtab, nspinor, nsppol, nspden, natom, &
653 : dtset%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg, &
654 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab, &
655 3 : usecprj=usecprj, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option)
656 :
657 : call gs_ham_kmp%init(psps, pawtab, nspinor, nsppol, nspden, natom, &
658 : dtset%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg, &
659 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab, &
660 3 : usecprj=usecprj, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option)
661 :
662 3 : matblk = gs_ham_kmp%matblk
663 :
664 : ! Allocate workspace arrays.
665 : ! vtrial and vlocal are required for Sternheimer (H0). DFPT routines do not need it.
666 : ! Note nvloc in vlocal (we will select one/four spin components afterwards)
667 3 : nvloc = gs_ham_kqmp%nvloc
668 14039 : ABI_CALLOC(vtrial, (nfftf, nspden))
669 16703 : ABI_CALLOC(vlocal, (n4, n5, n6, nvloc))
670 :
671 : ! Read the GS potential (vtrial) from input POT file
672 : ! In principle one may store GS vtrial in the DVDB but getpot_filepath is simpler to implement.
673 3 : call wrtout(units, sjoin(" Reading KS GS potential for Sternheimer from: ", dtfil%filpotin))
674 : call read_rhor(dtfil%filpotin, cplex1, nspden, nfftf, ngfftf, pawread0, mpi_enreg, vtrial, pot_hdr, pot_pawrhoij, comm, &
675 3 : allow_interp=.True., want_varname="vtrial")
676 3 : pot_cryst = pot_hdr%get_crystal()
677 3 : if (cryst%compare(pot_cryst, header=" Comparing input crystal with POT crystal") /= 0) then
678 0 : ABI_ERROR("Crystal structure from WFK and POT do not agree! Check messages above!")
679 : end if
680 3 : call pot_cryst%free(); call pot_hdr%free()
681 :
682 : ! Find correspondence IBZ --> set of q-points in DVDB.
683 : ! need_ftinterp selects whether DFPT potentials should be read from the DVDB or Fourier-interpolated on the fly.
684 : ! Activate FT interpolation automatically if required q-points in the IBZ are not found in the DVDB.
685 :
686 : ! qibz2dvdb gives the mapping gstore%ibz --> dvdb%ibz
687 : ! TODO: Make sure we have the same ibz in rho1%ibz
688 3 : need_ftinterp = .False.
689 9 : ABI_MALLOC(qibz2dvdb, (gstore%nqibz))
690 3 : if (dvdb%find_qpts(gstore%nqibz, gstore%qibz, qibz2dvdb, comm) /= 0) then
691 0 : call wrtout(units, " Cannot find eph_ngqpt_fine q-points in DVDB --> Activating Fourier interpolation.")
692 0 : need_ftinterp = .True.
693 : else
694 3 : call wrtout(units, " DVDB file contains all q-points in the IBZ --> Reading DFPT potentials from file.")
695 3 : need_ftinterp = .False.
696 : end if
697 : !need_ftinterp = .True.
698 :
699 : ! Distribute DFPT potentials (IBZ q-points) inside qpt_comm.
700 : ! Note that we distribute IBZ instead of the full BZ or the IBZ_k inside the loop over ikcalc.
701 : ! This means that the load won't be equally distributed but memory will scale with qpt_comm%nproc.
702 : ! To reduce load imbalance, we sort the qibz points by norm and use cyclic distribution inside qpt_comm
703 :
704 : ! itreat_qibz(nqibz)
705 : ! Table used to distribute potentials over q-points in the IBZ.
706 : ! The loop over qpts in the IBZ(k) is MPI distributed inside qpt_comm according to this table.
707 : ! 0 if this IBZ point is not treated by this proc.
708 : ! 1 if this IBZ is treated.
709 :
710 : ! TODO: Recheck this part
711 28 : ABI_ICALLOC(itreat_qibz, (gstore%nqibz))
712 22 : itreat_qibz = 1
713 : !call sort_rpts(gstore%nqibz, gstore%qibz, cryst%gmet, iperm)
714 : !do ii=1,gstore%nqibz
715 : ! iq_ibz = iperm(ii)
716 : ! do my_spin=1,gstore%my_nspins
717 : ! gqk => gstore%gqk(my_spin)
718 : ! if (mod(ii, gqk%qpt_comm%nproc) == gqk%qpt_comm%me) itreat_qibz(iq_ibz) = 1
719 : ! end do
720 : !end do
721 : !ABI_FREE(iperm)
722 :
723 22 : call wrtout(std_out, sjoin("P Number of q-points in the IBZ treated by this proc: " ,itoa(count(itreat_qibz == 1))))
724 :
725 3 : if (need_ftinterp) then
726 : ! Use ddb_ngqpt q-mesh to compute the real-space representation of DFPT v1scf_qq potentials to prepare Fourier interpolation.
727 : ! R-points are distributed inside comm_rpt
728 : ! Note that when R-points are distributed inside qpt_comm we cannot interpolate potentials on-the-fly
729 : ! inside the loop over q-points.
730 : ! In this case, indeed, the interpolation must be done in gwpt_setup_qloop once we know the q-points contributing
731 : ! to the integral and the potentials must be cached.
732 : !FIXME: qpt_comm is buggy.
733 : !if (gwpt%imag_only) comm_rpt = xmpi_comm_self
734 : !comm_rpt = gwpt%bsum_comm%value
735 0 : comm_rpt = xmpi_comm_self
736 0 : qptopt = ebands%kptopt; if (dtset%qptopt /= 0) qptopt = dtset%qptopt
737 0 : call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt)
738 0 : call drhodb%ftinterp_setup(dtset%ddb_ngqpt, qptopt, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt)
739 : end if
740 :
741 6 : call dvdb%print([std_out], "DVDB file", dtset%prtvol)
742 6 : call drhodb%print([std_out], "DRHODB file", dtset%prtvol)
743 3 : call pstat_proc%print(_PSTAT_ARGS_)
744 :
745 3 : if (.not. need_ftinterp) then
746 : ! Need to translate itreat_qibz into itreatq_dvdb.
747 28 : ABI_ICALLOC(itreatq_dvdb, (dvdb%nqpt))
748 22 : do iq_ibz=1,gstore%nqibz
749 19 : if (itreat_qibz(iq_ibz) == 0) cycle
750 19 : db_iqpt = qibz2dvdb(iq_ibz)
751 19 : ABI_CHECK(db_iqpt /= -1, sjoin("Could not find IBZ q-point:", ktoa(gstore%qibz(:, iq_ibz)), "in the DVDB file."))
752 22 : itreatq_dvdb(db_iqpt) = 1
753 : end do
754 3 : ABI_FREE(itreatq_dvdb)
755 : end if
756 :
757 3 : ABI_FREE(itreat_qibz)
758 3 : ABI_FREE(qibz2dvdb)
759 :
760 : ! Read the GS density needed to compute vxc1_qq.
761 14039 : ABI_CALLOC(rhor, (nfftf, nspden))
762 : call read_rhor(dtfil%fildensin, cplex1, nspden, nfftf, ngfftf, pawread0, mpi_enreg, rhor, den_hdr, den_pawrhoij, comm, &
763 3 : allow_interp=.True., want_varname="density")
764 3 : den_cryst = den_hdr%get_crystal()
765 3 : if (cryst%compare(den_cryst, header=" Comparing input crystal with DEN crystal") /= 0) then
766 0 : ABI_ERROR("Crystal structures from WFK and DEN do not agree! Check messages above!")
767 : end if
768 3 : call den_cryst%free(); call den_hdr%free()
769 :
770 : ! Initialize plasmon-pole object.
771 3 : mqmem = pp_mesh%nibz
772 3 : drude_plsmf = sqrt(four_pi * ebands%nelect / cryst%ucvol)
773 3 : my_plsmf = drude_plsmf; if (dtset%ppmfrq > tol6) my_plsmf = dtset%ppmfrq
774 3 : call ppm%init(mqmem, pp_mesh%nibz, npw_c, dtset%ppmodel, my_plsmf, dtset%gw_invalid_freq)
775 3 : if (my_rank == master) call ppm%print(units)
776 :
777 : ! Read symmetrized em1 from file and build ppmodel parameters.
778 : ! TODO: MPI-shared memory + compute only my set of pp-vectors in ppm%new_setup
779 3 : scr_iomode = iomode_from_fname(screen_filepath)
780 15 : ABI_MALLOC(epsm1_ggw, (npw_c, npw_c, hscr%nomega))
781 :
782 3 : if (dtset%gwcomp == 2) then
783 : ! Allocate memory to store static screening in the IBZ and workspace array for the value at pp_bz in the BZ.
784 : ! TODO: Store only the ip_ibz needed by this MPI rank to save memory.
785 0 : call wrtout(units, " Activating COH remainder technique")
786 0 : ABI_MALLOC(wc0_pibz, (pp_mesh%nibz))
787 0 : ABI_MALLOC(wc0_pbz, (npw_c, npw_c))
788 : end if
789 :
790 22 : do iq_ibz=1,pp_mesh%nibz
791 : call read_screening("inverse_dielectric_function", screen_filepath, &
792 19 : npw_c, 1, hscr%nomega, epsm1_ggw, scr_iomode, comm, iqiA=iq_ibz)
793 :
794 19 : call ppm%new_setup(iq_ibz, cryst, pp_mesh, npw_c, hscr%nomega, hscr%omega, epsm1_ggw, nfftf, gsph_c%gvec, ngfftf, rhor(:,1))
795 :
796 22 : if (dtset%gwcomp == 2) then
797 : ! Keep static limit of Wc in memory.
798 0 : do ii=1,npw_c
799 0 : epsm1_ggw(ii,ii,1) = epsm1_ggw(ii,ii,1) - one
800 : end do
801 0 : ABI_MALLOC(wc0_pibz(iq_ibz)%vals, (npw_c, npw_c))
802 0 : wc0_pibz(iq_ibz)%vals = epsm1_ggw(:,:,1)
803 : end if
804 : end do ! iq_ibz
805 :
806 3 : ABI_FREE(epsm1_ggw)
807 3 : call hscr%free()
808 :
809 : ! Allocate g-vectors centered on k, k+q, k-p, and k+q-p.
810 9 : ABI_MALLOC(kg_k, (3, mpw))
811 6 : ABI_MALLOC(kg_kq, (3, mpw))
812 6 : ABI_MALLOC(kg_kmp, (3, mpw))
813 6 : ABI_MALLOC(kg_kqmp, (3, mpw))
814 :
815 : ! GS wavefunctions
816 9 : ABI_MALLOC(ur_star_kmp, (nfft*nspinor))
817 6 : ABI_MALLOC(ur_star_kqmp, (nfft*nspinor))
818 6 : ABI_MALLOC(cwork_ur, (nfft*nspinor))
819 9 : ABI_MALLOC(cg_kmp, (2, mpw*nspinor))
820 6 : ABI_MALLOC(cg_kqmp, (2, mpw*nspinor))
821 : ! First order change (full term including the active space).
822 6 : ABI_MALLOC(cg1_kqmp, (2, mpw*nspinor))
823 6 : ABI_MALLOC(cg1_kmp, (2, mpw*nspinor))
824 :
825 : !stern_use_cache = merge(.True., .False., dtset%eph_stern == 1)
826 3 : stern_use_cache = .False.
827 :
828 6 : if (my_rank == master) call gstore%print([std_out])
829 3 : call pstat_proc%print(_PSTAT_ARGS_)
830 :
831 : ! TODO
832 : !if (dtset%gstore_iv1p_comm /= 0) then
833 : ! call gstore%compute_and_write_commutator(mpw, gmax, ngfft, ngfftf, dtset, cryst, pawfgr, psps, &
834 : ! wfd, mpi_enreg, kg_k, ebands, dvdb, gs_ham_kq, root_ncid)
835 : !end if
836 :
837 : ! This parameter defines the size of the q-buffer used to store the g(k, q) e-ph matrix elements
838 : ! for all the k-point treated by this MPI rank.
839 : ! Increasing the buffer size increases the memory requirements
840 : ! but it leads to better performance as the number of IO operations is decreased.
841 : ! TODO: Should compute it on the basis of my_nkpt and my_nqpt
842 3 : qbuf_size = 1
843 : !qbuf_size = 16
844 3 : call wrtout(std_out, sjoin(" Begin computation of GWPT e-ph matrix elements with qbuf_size:", itoa(qbuf_size)), pre_newlines=1)
845 :
846 : ! If True, only k-points in the IBZ and q-points in the IBZ_k are computed.
847 : ! Matrix elements in full BZs are then reconstructed by symmetry at the end of the run by calling
848 : ! gstore_symmetrize.
849 3 : symmetrize = (dtset%gstore_kzone == "bz" .and. dtset%gstore_qzone == "bz" .and. dtset%gstore_sym > 0)
850 :
851 : ! Compute kxc needed for vxc1.
852 : ! A similar piece of code is used in m_respfn_driver.
853 : ! option 2 for xc and kxc (no paramagnetic part if xcdata%nspden=1).
854 : ! Note usage of dum_xccc3d to ignore non-linear core correction.
855 3 : nkxc = 2*min(dtset%nspden,2)-1; if (dtset%xclevel==2) nkxc = 12*min(dtset%nspden,2)-5
856 3 : call xcdata_init(xcdata, dtset=dtset)
857 3 : non_magnetic_xc = (dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
858 :
859 3 : call check_kxc(dtset%ixc, dtset%optdriver)
860 12 : ABI_MALLOC(kxc, (nfft, nkxc))
861 12 : ABI_MALLOC(vxc, (nfft, dtset%nspden))
862 :
863 3 : nk3xc=1; option=2; usexcnhat=0
864 : call rhotoxc(bigexc, bigsxc, kxc, mpi_enreg, nfft, ngfft, &
865 : dum_nhat, 0, dum_nhat, 0, nkxc, nk3xc, non_magnetic_xc, n3xccc0, option, rhor, &
866 3 : cryst%rprimd, usexcnhat, vxc, vxcavg, dum_xccc3d, xcdata)
867 :
868 : ! Here we decide if the q-points can be reduced to the IBZ(k)
869 3 : use_lgk = (dtset%gstore_use_lgk /= 0 .or. dtset%gstore_sym == 2)
870 : if (use_lgk) then
871 3 : call wrtout(units, " Only q-points in the IBZ_k will be computed.")
872 0 : else if (dtset%gstore_use_lgq /= 0) then
873 0 : call wrtout(units, " Only k-points in the IBZ_q will be computed.")
874 : else
875 0 : call wrtout(units, " Little group operations won't be used")
876 : end if
877 :
878 : ! ===================================================
879 : ! Loop over MPI distributed spins in Sigma (gqk%comm)
880 : ! ===================================================
881 3 : call pstat_proc%print(_PSTAT_ARGS_)
882 :
883 3 : stern_qq_ierr = 0; stern_mq_ierr = 0
884 :
885 6 : do my_is=1,gstore%my_nspins
886 3 : spin = gstore%my_spins(my_is); gqk => gstore%gqk(my_is); my_npert = gqk%my_npert
887 : ABI_CHECK_IEQ(my_npert, gqk%my_npert, "my_npert")
888 3 : ABI_CHECK_IEQ(gqk%bstart_kq, 1, "bstart_qk should be 1")
889 3 : ABI_CHECK_IGEQ(nbsum, gqk%bstop_kq, "nband must be greater than the max band in the e-ph matrix elements")
890 :
891 3 : NCF_CHECK(nf90_inq_ncid(root_ncid, strcat("gqk", "_spin", itoa(spin)), spin_ncid))
892 3 : NCF_CHECK(nctk_prepare_mpiio(root_ncid, "gstore_done_qbz_spin"))
893 3 : NCF_CHECK(nctk_prepare_mpiio(spin_ncid, "gvals"))
894 3 : NCF_CHECK(nctk_prepare_mpiio(spin_ncid, "gvals_ks"))
895 :
896 : ! Note the possibility of specifying different number of states for the incoming and the intermediate states.
897 3 : nb_k = gqk%nb_k; bstart_k = gqk%bstart_k; bstop_k = gqk%bstop_k
898 3 : nb_kq = gqk%nb_kq; bstart_kq = gqk%bstart_kq; bstop_kq = gqk%bstop_kq
899 :
900 9 : ABI_MALLOC(iq_buf, (2, qbuf_size))
901 12 : ABI_MALLOC(state_kq, (gqk%my_nk, qbuf_size))
902 11 : state_kq = GSTORE_KQ_MISSING
903 15 : ABI_MALLOC(gsig_atm, (2, nb_kq, nb_k, natom3))
904 12 : ABI_MALLOC(gks_atm, (2, nb_kq, nb_k, natom3))
905 12 : ABI_MALLOC(gks_atm2, (2, nb_kq, nb_k, natom3))
906 12 : ABI_MALLOC(gxc_atm, (2, nb_kq, nb_k, natom3))
907 :
908 : ! Arrays used to compare gwpt with gw
909 : if (test_sigma) then
910 : ABI_CALLOC(vxc_nk, (nb_k, gqk%glob_nk))
911 : ABI_CALLOC(sigx_nk, (nb_k, gqk%glob_nk))
912 : ABI_CALLOC(sigce0_nk, (nb_k, gqk%glob_nk))
913 : end if
914 :
915 12 : ABI_MALLOC(ur_nk, (nfft*nspinor, bstart_k:bstop_k))
916 12 : ABI_MALLOC(ur_mkq, (nfft*nspinor, bstart_kq:bstop_kq))
917 :
918 : ! Inside the loops, we compute gsig_atm(2, nb_kq, nb_k, natom3)
919 24 : ABI_MALLOC(my_gbuf, (gqk%cplex, nb_kq, nb_k, natom3, gqk%my_nk, qbuf_size))
920 24 : ABI_MALLOC(my_gbuf_ks, (gqk%cplex, nb_kq, nb_k, natom3, gqk%my_nk, qbuf_size))
921 :
922 : ! Allocate memory to deal with frequencies in Sigma(w).
923 :
924 5 : select case (dtset%gwpt_wmode)
925 : case (1)
926 : ! Prepare list of omegas: first e_nk then e_mkq for all m indices.
927 : ! Note that nw_nk depends on bstart_kq and bstop_kq so nw_nk depends on bands at _kq.
928 2 : nw_nk = 1 + (bstop_kq - bstart_kq + 1)
929 2 : nw_mkq = 1 + (bstop_k - bstart_k + 1)
930 2 : call wrtout(units, " GWPT frequency convolution is evaluated at E_nk and E_mkq$ and the average is taken.")
931 : case (2)
932 1 : nw_nk = 1
933 1 : nw_mkq = (bstop_k - bstart_k + 1)
934 1 : call wrtout(units, " GWPT frequency convolution is evaluated at E_nk.")
935 : case default
936 3 : ABI_ERROR(sjoin("Invalid gwpt_wmode:", itoa(dtset%gwpt_wmode)))
937 : end select
938 :
939 9 : ABI_MALLOC(omegame0i_nk, (nw_nk))
940 9 : ABI_MALLOC(omegame0i_mkq, (nw_mkq))
941 6 : ABI_MALLOC(omegas_nk, (nw_nk))
942 6 : ABI_MALLOC(omegas_mkq, (nw_mkq))
943 9 : ABI_MALLOC(sigcme_nk, (nw_nk))
944 9 : ABI_MALLOC(sigcme_mkq, (nw_mkq))
945 :
946 : ! Correlated part
947 15 : ABI_MALLOC(vec_gwc_nk, (npw_c*nspinor, nw_nk, bstart_k:bstop_k))
948 15 : ABI_MALLOC(vec_gwc_mkq, (npw_c*nspinor, nw_mkq, bstart_kq:bstop_kq))
949 : ! Exchange part
950 12 : ABI_MALLOC(vec_gx_nk, (npw_x*nspinor, bstart_k:bstop_k))
951 12 : ABI_MALLOC(vec_gx_mkq, (npw_x*nspinor, bstart_kq:bstop_kq))
952 :
953 3 : if (dtset%gwcomp == 2) then
954 0 : ABI_MALLOC(vec_coh_nk, (npw_c, bstart_k:bstop_k))
955 0 : ABI_MALLOC(vec_coh_mkq, (npw_c, bstart_kq:bstop_kq))
956 : end if
957 :
958 : ! Compute the little group of the k-point so that we can compute g(k,q) only for q in the IBZ_k
959 3 : if (use_lgk) then
960 3 : timrev_k = kpts_timrev_from_kptopt(ebands%kptopt)
961 14 : ABI_MALLOC(lg_myk, (gqk%my_nk))
962 8 : do my_ik=1,gqk%my_nk
963 20 : kk = gqk%my_kpts(:, my_ik)
964 8 : call lg_myk(my_ik)%init(cryst, kk, timrev_k, gstore%nqbz, gstore%qbz, gstore%nqibz, gstore%qibz, xmpi_comm_self)
965 : end do
966 : end if
967 :
968 : ! ============================================================
969 : ! Loop over MPI distributed q-points in Sigma_q (gqk%qpt_comm)
970 : ! ============================================================
971 : ! the loop of my_iq can be cycled, so it is good to set iqbuf_cnt to 0 before entering the loop
972 3 : iqbuf_cnt = 0
973 83 : do my_iq=1,gqk%my_nq
974 80 : call gqk%myqpt(my_iq, gstore, weight_q, qq_bz)
975 :
976 115 : if (dtset%userib /= 0 .and. (any(abs(qq_bz - [0.5, 0.0, 0.0]) > tol14))) cycle
977 :
978 80 : iq_bz = gqk%my_q2bz(my_iq); qq_is_gamma = sum(qq_bz**2) < tol14; qq_bz_string = ktoa(qq_bz)
979 :
980 : ! Handle possible restart.
981 80 : if (done_qbz_spin(iq_bz, spin) == 1) then
982 0 : call wrtout(std_out, sjoin(" iq_bz:", itoa(iq_bz), ", spin: ", itoa(spin), " already computed --> skipping iteration"), do_flush=.True.)
983 0 : cycle
984 : end if
985 :
986 : ! Compute the little group of the q-point so that we can compute g(k,q) only for k in the IBZ_q
987 80 : if (dtset%gstore_use_lgq /= 0) then
988 0 : timrev_q = kpts_timrev_from_kptopt(gstore%qptopt)
989 0 : call lg_myq%init(cryst, qq_bz, timrev_q, gstore%nkbz, gstore%kbz, gstore%nkibz, gstore%kibz, xmpi_comm_self)
990 : end if
991 :
992 : ! Note symrec conventions here as needed to symmetrize the DFPT potentials.
993 80 : iq_ibz = gqk%my_q2ibz(1, my_iq) ; isym_q = gqk%my_q2ibz(2, my_iq)
994 80 : trev_q = gqk%my_q2ibz(6, my_iq) ; g0_q = gqk%my_q2ibz(3:5, my_iq)
995 : isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0))
996 : !tsign_q = 1; if (trev_q == 1) tsign_q = -1
997 320 : qq_ibz = gstore%qibz(:, iq_ibz)
998 560 : mapc_qq = gqk%my_q2ibz(:, my_iq)
999 :
1000 80 : print_time_qq = my_rank == 0 .and. (my_iq <= LOG_MODQ .or. mod(my_iq, LOG_MODQ) == 0)
1001 80 : if (print_time_qq) then
1002 80 : call cwtime(cpu_qq, wall_qq, gflops_qq, "start")
1003 80 : call inds2str(0, sjoin(" Computing g^Sigma(k, q) for qq_bz:", qq_bz_string), my_iq, gqk%my_nq, gqk%glob_nq, msg)
1004 80 : call wrtout(std_out, sjoin(msg, ", and spin:", itoa(spin)), pre_newlines=1)
1005 : !print *, "iq_ibz:", iq_ibz, "qq_bz:", qq_bz, "qq_ibz:", qq_ibz
1006 : end if
1007 :
1008 : ! ==================================================
1009 : ! Get DFPT potentials and densities for this q-point
1010 : ! ==================================================
1011 : ! After this branch we know `cplex` and we have allocated:
1012 : !
1013 : ! v1scf_qq(cplex, nfftf, nspden, my_npert))
1014 : ! vxc1_qq(cplex, nfft, nspden, my_npert)
1015 : !
1016 : ! Important: vxc1_qq does not include the contribution due to the model core charge (if any).
1017 :
1018 80 : if (need_ftinterp) then
1019 : ! Use Fourier interpolation to get DFPT potentials and DFPT densities for this qpt.
1020 0 : call dvdb%get_ftqbz(qq_bz, cplex, nfftf, ngfftf, v1scf_qq, gqk%pert_comm%value)
1021 :
1022 : call drhodb%get_vxc1_ftqbz(dtset, cryst, qq_bz, drho_cplex, nfftf, ngfftf, nkxc, kxc, &
1023 0 : vxc1_qq, non_magnetic_xc, usexcnhat, gqk%pert_comm%value)
1024 : else
1025 : ! Read and reconstruct the dvscf potentials and the densities for this qq_bz and my_npert perturbations.
1026 80 : db_iqpt = dvdb%findq(qq_ibz)
1027 80 : ABI_CHECK(db_iqpt /= -1, sjoin("Could not find symmetric of q-point:", ktoa(qq_bz), "in DVDB file."))
1028 : ! The first entry in mapc_qq2dvdb gives the index in dvdb%qpts. The other entries in mapc_qq are OK as they refer to symmetries.
1029 80 : mapc_qq2dvdb = mapc_qq; mapc_qq2dvdb(1) = db_iqpt
1030 80 : call dvdb%readsym_qbz(cryst, qq_bz, mapc_qq2dvdb, cplex, nfftf, ngfftf, v1scf_qq, gqk%pert_comm%value)
1031 :
1032 80 : db_iqpt = drhodb%findq(qq_ibz)
1033 80 : ABI_CHECK(db_iqpt /= -1, sjoin("Could not find symmetric of q-point:", ktoa(qq_bz), "in DRHODB file."))
1034 80 : mapc_qq2dvdb = mapc_qq; mapc_qq2dvdb(1) = db_iqpt
1035 : call drhodb%read_vxc1_qbz(dtset, cryst, qq_bz, mapc_qq2dvdb, drho_cplex, nfftf, ngfftf, nkxc, kxc, &
1036 80 : vxc1_qq, non_magnetic_xc, usexcnhat, gqk%pert_comm%value)
1037 : end if
1038 :
1039 80 : ABI_CHECK_IEQ(cplex, drho_cplex, "Different values of cplex for v1 and rho1!")
1040 311 : cvxc1_qq_ptr => null(); if (cplex == 2) call c_f_pointer(c_loc(vxc1_qq), cvxc1_qq_ptr, [nfft, nspden, my_npert])
1041 :
1042 : ! Allocate vlocal1_qq with correct cplex and nvloc.
1043 560 : ABI_MALLOC(vlocal1_qq, (cplex*n4, n5, n6, nvloc, my_npert))
1044 480 : ABI_MALLOC(vlocal1_mqq, (cplex*n4, n5, n6, nvloc, my_npert))
1045 :
1046 : ! Build DFPT potential at -qq by taking CC.
1047 10364596 : v1scf_mq = v1scf_qq; if (cplex == 2) v1scf_mq(2,:,:,:) = -v1scf_mq(2,:,:,:)
1048 :
1049 : ! =============================================================
1050 : ! Loop over k-points in the e-ph matrix elements (gqk%kpt_comm)
1051 : ! =============================================================
1052 176 : do my_ik=1,gqk%my_nk
1053 384 : kk = gqk%my_kpts(:, my_ik); kk_string = ktoa(kk); ik_glob = gqk%my_k2glob(my_ik)
1054 :
1055 : ! NB: All procs in gqk%pert_comm and gqk%bsum_comm and gqk%pp_sum_comm enter this section.
1056 : ! Set entry to zero BEFORE the cycle instructions below, otherwise cycled (filtered)
1057 : ! k-points may keep stale data in my_gbuf and write them to disk.
1058 96 : iqbuf_cnt = 1 + mod(my_iq - 1, qbuf_size)
1059 288 : iq_buf(:, iqbuf_cnt) = [my_iq, iq_bz]
1060 90912 : my_gbuf(:,:,:,:, my_ik, iqbuf_cnt) = zero
1061 90912 : my_gbuf_ks(:,:,:,:, my_ik, iqbuf_cnt) = zero
1062 :
1063 : ! The k-point and the symmetries relating the BZ k-point to the IBZ.
1064 96 : ik_ibz = gqk%my_k2ibz(1, my_ik) ; isym_k = gqk%my_k2ibz(2, my_ik)
1065 384 : trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5,my_ik)
1066 384 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
1067 672 : mapl_k = gqk%my_k2ibz(:, my_ik)
1068 :
1069 384 : kk_ibz = ebands%kptns(:,ik_ibz)
1070 :
1071 : ! If we are going to reconstruct g(k,q) in the full BZ by symmetry at the end of the run,
1072 : ! only k-points in the IBZ need to be computed here.
1073 96 : if (symmetrize .and. .not. isirr_k) then
1074 0 : state_kq(my_ik, iqbuf_cnt) = GSTORE_KQ_MISSING; cycle
1075 : end if
1076 :
1077 96 : if (dtset%userib /= 0) then
1078 0 : if (any(abs(gqk%my_kpts(:, my_ik) - [0.25, 0.0, 0.0]) > tol14) .and. &
1079 : any(abs(gqk%my_kpts(:, my_ik) - [-0.25, 0.0, 0.0]) > tol14)) then
1080 0 : state_kq(my_ik, iqbuf_cnt) = GSTORE_KQ_MISSING; cycle
1081 : end if
1082 : end if
1083 :
1084 : ! Here we skip points if little group tricks are activated.
1085 96 : if (use_lgk) then
1086 96 : if (lg_myk(my_ik)%findq_ibzk(qq_bz) == -1) then
1087 65 : state_kq(my_ik, iqbuf_cnt) = GSTORE_KQ_MISSING; cycle
1088 : end if
1089 : end if
1090 :
1091 31 : if (dtset%gstore_use_lgq /= 0) then
1092 0 : if (lg_myq%findq_ibzk(kk) == -1) then
1093 0 : state_kq(my_ik, iqbuf_cnt) = GSTORE_KQ_MISSING; cycle
1094 : end if
1095 : end if
1096 :
1097 36805 : gks_atm = zero
1098 36805 : gks_atm2 = zero
1099 :
1100 31 : call inds2str(0, sjoin(" Computing g^Sigma(k,q) for kpt:", kk_string), my_ik, gqk%my_nk, gqk%glob_nk, msg)
1101 31 : call wrtout(std_out, sjoin(msg, ", for spin:", itoa(spin)), pre_newlines=1)
1102 31 : call timab(1940, 1, tsec)
1103 :
1104 31 : istwf_k_ibz = wfd%istwfk(ik_ibz); npw_k_ibz = wfd%npwarr(ik_ibz)
1105 :
1106 31 : print_time_kk = my_rank == 0 .and. (my_ik <= LOG_MODK .or. mod(my_ik, LOG_MODK) == 0)
1107 31 : if (print_time_kk) call cwtime(cpu_kk, wall_kk, gflops_kk, "start")
1108 :
1109 : ! Get npw_k, kg_k for kk
1110 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kk, ik_ibz, isirr_k, dtset%nloalg, & ! in
1111 31 : istwf_k, npw_k, kg_k, nkpg_k, kpg_k, gbound_k) ! out
1112 31 : ABI_CHECK_ILEQ(npw_k, mpw, "npw_k > mpw!")
1113 :
1114 : ! Find k + q in the extended zone and extract symmetry info.
1115 : ! Be careful here because there are two umklapp vectors to be considered as:
1116 : !
1117 : ! k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz
1118 : !
1119 124 : kq = kk + qq_bz
1120 :
1121 31 : if (kpts_map("symrel", ebands%kptopt, cryst, gstore%krank_ibz, 1, kq, mapl_kq) /= 0) then
1122 0 : write(msg, '(4a)' )"k-mesh is not closed!",ch10, "k+q could not be generated from a symmetrical one.",trim(ltoa(kq))
1123 0 : ABI_ERROR(msg)
1124 : end if
1125 124 : ikq_ibz = mapl_kq(1); isym_kq = mapl_kq(2); trev_kq = mapl_kq(6); g0_kq = mapl_kq(3:5)
1126 99 : isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0))
1127 124 : kq_ibz = ebands%kptns(:, ikq_ibz)
1128 31 : istwf_kq_ibz = wfd%istwfk(ikq_ibz); npw_kq_ibz = wfd%npwarr(ikq_ibz)
1129 :
1130 : ! Get npw_kq, kg_kq for k+q.
1131 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kq, ikq_ibz, isirr_kq, dtset%nloalg, & ! in
1132 31 : istwf_kq, npw_kq, kg_kq, nkpg_kq, kpg_kq, gbound_kq) ! out
1133 31 : ABI_CHECK_ILEQ(npw_kq, mpw, "npw_kq > mpw!")
1134 :
1135 104952 : ABI_CALLOC(ug_k, (2, npw_k*nspinor))
1136 106581 : ABI_CALLOC(ug_kq, (2, npw_kq*nspinor))
1137 :
1138 : ! Precompute ur_nk and ur_mkq for all m and n band indices treated.
1139 : ! TODO: Can distribute operations inside gqk%pert_comm
1140 :
1141 201 : do n_k=bstart_k, bstop_k
1142 170 : in_k = n_k - bstart_k + 1 !; if (gqk%pert_ppsum_comm%skip(in_k)) cycle ! MPI parallelism inside pert_comm
1143 : call wfd%rotate_cg(n_k, ndat1, spin, kk_ibz, npw_k, kg_k, istwf_k, &
1144 170 : cryst, mapl_k, gbound_k, work_ngfft, work, ug_k, urs_kbz=ur_nk(:,n_k))
1145 31 : if (test_sigma) vxc_nk(in_k, ik_glob) = dot_product(ur_nk(:,n_k), vxc(:, spin) * ur_nk(:,n_k)) / nfftf
1146 : end do
1147 : !call xmpi_sum(ur_nk, gqk%pert_ppsum_comm%value, ierr)
1148 :
1149 387 : do m_kq=bstart_kq, bstop_kq
1150 356 : im_kq = m_kq - bstart_kq + 1 !; if (gqk%pert_ppsum_comm%skip(im_kq)) cycle ! MPI parallelism inside pert_comm
1151 : call wfd%rotate_cg(m_kq, ndat1, spin, kq_ibz, npw_kq, kg_kq, istwf_kq, &
1152 387 : cryst, mapl_kq, gbound_kq, work_ngfft, work, ug_kq, urs_kbz=ur_mkq(:,m_kq))
1153 : end do
1154 : !call xmpi_sum(ur_mkq, gqk%pert_ppsum_comm%value, ierr)
1155 :
1156 : ! ===========================
1157 : ! Compute <m,k+q|vxc1_qq|n,k>
1158 : ! ===========================
1159 :
1160 36805 : gxc_atm = czero; cnt = 0
1161 387 : do m_kq=bstart_kq, bstop_kq
1162 356 : im_kq = m_kq - bstart_kq + 1
1163 2363 : do n_k=bstart_k, bstop_k
1164 1976 : in_k = n_k - bstart_k + 1
1165 1976 : cnt = cnt + 1
1166 1976 : if (gqk%pp_sum_comm%skip(cnt)) cycle ! MPI parallelism inside pp_sum_comm
1167 14188 : do imyp=1,gqk%my_npert
1168 11856 : if (cplex == 1) then
1169 8479536 : ctmp_gwpc = sum(GWPC_CONJG(ur_mkq(:,m_kq)) * ur_nk(:,n_k) * vxc1_qq(1,:,spin,imyp)) / nfftf
1170 : else
1171 43427616 : ctmp_gwpc = sum(GWPC_CONJG(ur_mkq(:,m_kq)) * ur_nk(:,n_k) * cvxc1_qq_ptr(:,spin,imyp)) / nfftf
1172 : end if
1173 11856 : ipc = gqk%my_pertcases(imyp)
1174 37544 : gxc_atm(:, im_kq, in_k, ipc) = [real(ctmp_gwpc), aimag(ctmp_gwpc)]
1175 : end do ! imyp
1176 : end do ! n_k
1177 : end do ! m_kq
1178 :
1179 : ! TODO: this is an all_gatherv but oh well.
1180 31 : call xmpi_sum(gxc_atm, gqk%pert_ppsum_comm%value, ierr)
1181 31 : call timab(1940, 2, tsec)
1182 :
1183 : ! ==========================================================
1184 : ! MPI sum over the pp momenta in the full BZ gqk%pp_sum_comm
1185 : ! ==========================================================
1186 : !
1187 : ! Be careful here because pp should run over the list of wavevectors in the screening matrix!
1188 : ! as pp_mesh%bz is not necessarily equivalent to the k-mesh for the wavefunctions.
1189 : ! Moreover we have to use the ipp_bz index to symmetrize W(pp_bz) from W(pp_ibz).
1190 : !
1191 : ! TODO: Should order nbz in shells so that one can reduce the memory required
1192 : ! to store W(pp) if pp_parallelism is activated.
1193 36805 : gsig_atm = zero
1194 :
1195 727 : do ipp_bz=my_pp_start_spin(spin), my_pp_stop_spin(spin)
1196 : ! All procs in gqk%pert_comm and gqk%bsum_com enter this section.
1197 696 : call timab(1941, 1, tsec)
1198 :
1199 696 : my_ipp = ipp_bz - my_pp_start_spin(spin) + 1
1200 696 : print_time_pp = my_rank == 0 .and. (my_ipp <= LOG_MODP .or. mod(my_ipp, LOG_MODP) == 0)
1201 267 : if (print_time_pp) call cwtime(cpu_pp, wall_pp, gflops_pp, "start")
1202 :
1203 4872 : pp = pp_mesh%bz(:,ipp_bz); pp_is_gamma = sum(pp**2) < tol14; pp_string = ktoa(pp)
1204 :
1205 : ! Debug, include only pp=Gamma
1206 696 : if (dtset%userie > 0) then
1207 0 : if (.not. pp_is_gamma) cycle
1208 : end if
1209 :
1210 696 : qkp_string = sjoin("While treating qq_bz: ", qq_bz_string, "kpt:", kk_string, "pp:", pp_string, ch10)
1211 :
1212 : ! Symmetry tables and g-sphere centered on k-p.
1213 2784 : kmp = kk - pp
1214 696 : if (kpts_map("symrel", ebands%kptopt, cryst, gstore%krank_ibz, 1, kmp, mapl_kmp) /= 0) then
1215 0 : write(msg, '(4a)' )"k-mesh is not closed!",ch10, &
1216 0 : "k-p could not be generated from a symmetrical one.",trim(ltoa(kmp))
1217 0 : ABI_ERROR(msg)
1218 : end if
1219 2784 : ikmp_ibz = mapl_kmp(1); isym_kmp = mapl_kmp(2); trev_kmp = mapl_kmp(6); g0_kmp = mapl_kmp(3:5)
1220 2349 : isirr_kmp = (isym_kmp == 1 .and. trev_kmp == 0 .and. all(g0_kmp == 0))
1221 2784 : kmp_ibz = ebands%kptns(:, ikmp_ibz)
1222 :
1223 : ! Get npw_kmp, kg_kmp for k-p.
1224 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kmp, ikmp_ibz, isirr_kmp, dtset%nloalg, & ! in
1225 696 : istwf_kmp, npw_kmp, kg_kmp, nkpg_kmp, kpg_kmp, gbound_kmp) ! out
1226 696 : ABI_CHECK_ILEQ(npw_kmp, mpw, "npw_kmp > mpw!")
1227 :
1228 2784 : ABI_MALLOC(ph3d_kmp, (2, npw_kmp, matblk))
1229 2088 : ABI_MALLOC(kinpw_kmp, (npw_kmp))
1230 696 : call mkkin(dtset%ecut, dtset%ecutsm, dtset%effmass_free, cryst%gmet, kg_kmp, kinpw_kmp, kmp, npw_kmp, 0, 0)
1231 :
1232 : ! Compute nonlocal form factors ffnl_kmp at (k-p+G).
1233 3480 : ABI_MALLOC(ffnl_kmp, (npw_kmp, 1, psps%lmnmax, psps%ntypat))
1234 :
1235 2784 : ABI_MALLOC(ylm_kmp, (npw_kmp, psps%mpsang**2 * psps%useylm))
1236 696 : if (psps%useylm == 1) then
1237 0 : call initylmg_k(npw_kmp, psps%mpsang, optder0, cryst%rprimd, cryst%gprimd, kmp, kg_kmp, ylm_kmp, ylmgr_dum)
1238 : end if
1239 :
1240 : call mkffnl_objs(cryst, psps, 1, ffnl_kmp, ider0, idir0, kg_kmp, kpg_kmp, kmp, nkpg_kmp, &
1241 696 : npw_kmp, ylm_kmp, ylmgr_dum) !, comm=gqk%pert_comm%value, request=ffnl_kmp_request)
1242 :
1243 : ! Symmetry tables and g-sphere centered on k+q-p.
1244 2784 : kqmp = kq - pp
1245 696 : if (kpts_map("symrel", ebands%kptopt, cryst, gstore%krank_ibz, 1, kqmp, mapl_kqmp) /= 0) then
1246 0 : ABI_ERROR(sjoin("k-mesh not closed!",ch10, "k+q-p could not be generated from a symmetrical one:",trim(ltoa(kqmp))))
1247 : end if
1248 2784 : ikqmp_ibz = mapl_kqmp(1); isym_kqmp = mapl_kqmp(2); trev_kqmp = mapl_kqmp(6); g0_kqmp = mapl_kqmp(3:5)
1249 1832 : isirr_kqmp = (isym_kqmp == 1 .and. trev_kqmp == 0 .and. all(g0_kqmp == 0))
1250 2784 : kqmp_ibz = ebands%kptns(:, ikqmp_ibz)
1251 696 : istwf_kqmp_ibz = wfd%istwfk(ikqmp_ibz); npw_kqmp_ibz = wfd%npwarr(ikqmp_ibz)
1252 :
1253 : ! Get npw_kqmp, kg_kqmp for k+q-p.
1254 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kqmp, ikqmp_ibz, isirr_kqmp, dtset%nloalg, & ! in
1255 696 : istwf_kqmp, npw_kqmp, kg_kqmp, nkpg_kqmp, kpg_kqmp, gbound_kqmp) ! out
1256 696 : ABI_CHECK_ILEQ(npw_kqmp, mpw, "npw_kqmp > mpw!")
1257 :
1258 2088 : ABI_MALLOC(full_cg1_kqmp, (2, npw_kqmp*nspinor))
1259 2088 : ABI_MALLOC(full_cg1_kmp, (2, npw_kmp*nspinor))
1260 :
1261 2784 : ABI_MALLOC(ph3d_kqmp, (2, npw_kqmp, matblk))
1262 2088 : ABI_MALLOC(kinpw_kqmp, (npw_kqmp))
1263 696 : call mkkin(dtset%ecut, dtset%ecutsm, dtset%effmass_free, cryst%gmet, kg_kqmp, kinpw_kqmp, kqmp, npw_kqmp, 0, 0)
1264 :
1265 : ! Compute nonlocal form factors ffnl_kqmp at (k+q-p+G).
1266 3480 : ABI_MALLOC(ffnl_kqmp, (npw_kqmp, 1, psps%lmnmax, psps%ntypat))
1267 :
1268 2784 : ABI_MALLOC(ylm_kqmp, (npw_kqmp, psps%mpsang**2 * psps%useylm))
1269 696 : if (psps%useylm == 1) then
1270 0 : call initylmg_k(npw_kqmp, psps%mpsang, optder0, cryst%rprimd, cryst%gprimd, kqmp, kg_kqmp, ylm_kqmp, ylmgr_dum)
1271 : end if
1272 :
1273 : call mkffnl_objs(cryst, psps, 1, ffnl_kqmp, ider0, idir0, kg_kqmp, kpg_kqmp, kqmp, nkpg_kqmp, &
1274 696 : npw_kqmp, ylm_kqmp, ylmgr_dum) ! , comm=gqk%pert_comm%value, request=ffnl_kqmp_request)
1275 :
1276 : ! =============================================
1277 : ! Load the k dependent parts of the Hamiltonian
1278 : ! =============================================
1279 :
1280 : ! (k+q-p, k-p)
1281 : call gs_ham_kqmp%load_k(kpt_k=kmp, npw_k=npw_kmp, istwf_k=istwf_kmp, kg_k=kg_kmp, kpg_k=kpg_kmp, kinpw_k=kinpw_kmp, &
1282 696 : ph3d_k=ph3d_kmp, ffnl_k=ffnl_kmp, compute_ph3d=.true., compute_gbound=.true.)
1283 :
1284 : call gs_ham_kqmp%load_kprime(kpt_kp=kqmp, npw_kp=npw_kqmp, istwf_kp=istwf_kqmp, kg_kp=kg_kqmp, kpg_kp=kpg_kqmp, &
1285 696 : kinpw_kp=kinpw_kqmp, ph3d_kp=ph3d_kqmp, ffnl_kp=ffnl_kqmp, compute_ph3d=.true., compute_gbound=.true.)
1286 :
1287 : ! (k-p, k+q-p)
1288 : call gs_ham_kmp%load_k(kpt_k=kqmp, npw_k=npw_kqmp, istwf_k=istwf_kqmp, kg_k=kg_kqmp, kpg_k=kpg_kqmp, kinpw_k=kinpw_kqmp, &
1289 696 : ph3d_k=ph3d_kqmp, ffnl_k=ffnl_kqmp, compute_ph3d=.true., compute_gbound=.true.)
1290 :
1291 : call gs_ham_kmp%load_kprime(kpt_kp=kmp, npw_kp=npw_kmp, istwf_kp=istwf_kmp, kg_kp=kg_kmp, kpg_kp=kpg_kmp, &
1292 696 : kinpw_kp=kinpw_kmp, ph3d_kp=ph3d_kmp, ffnl_kp=ffnl_kmp, compute_ph3d=.true., compute_gbound=.true.)
1293 :
1294 : ! ====================================
1295 : ! This is the g-sphere for W_{gg'}(pp)
1296 : ! ====================================
1297 : ! Note that in this case, the spheres are always Gamma-centered i.e. it does not depend on the pp wavevector.
1298 696 : kg_c => gsph_c%gvec(:, 1:npw_c)
1299 696 : kg_x => gsph_x%gvec(:, 1:npw_x)
1300 696 : call sphereboundary(gbound_c, istwfk1, kg_c, mgfft, npw_c)
1301 696 : call sphereboundary(gbound_x, istwfk1, kg_x, mgfft, npw_x)
1302 :
1303 2088 : ABI_MALLOC(rhotwg_c, (npw_c*nspinor))
1304 2088 : ABI_MALLOC(rhotwg_x, (npw_x*nspinor))
1305 2088 : ABI_MALLOC(vc_sqrt_gx, (npw_x))
1306 :
1307 : ! We need two stern_t objects to compute the first order change of the wavefunctions at k-p and k+q-p.
1308 : ! Clearly, we should not duplicate the work when pp = 0.
1309 : ! When pp == 0, we also get the gks matrix elements after stern_solve.
1310 : ! Alternatively, one can solve the Sternheimer in the IBZ(kappa, alpha), store the results on disk
1311 : ! and then use symmetries to reconstruct delta_u in the full BZ on the fly assuming spatial inversion or TR.
1312 : ! Also, one should handle more carefully the integration in g_sigma around pp = Gamma in the case of semiconductors.
1313 :
1314 696 : nband_me = nbsum; stern_comm = xmpi_comm_self
1315 :
1316 : ! =======================
1317 : ! Delta_{qka} psi_{n'k-p}
1318 : ! =======================
1319 : ! (k+q-p, k-p)
1320 : call stern_kmp%init(dtset, npw_kmp, npw_kqmp, nspinor, nbsum, nband_me, fermie1_idir_ipert, &
1321 696 : stern_use_cache, work_ngfft, mpi_enreg, stern_comm)
1322 :
1323 : ! Use symmetries to reconstruct GS wavefunctions at k+q-p and store them in stern_kmp%cgq.
1324 8024 : do ib_sum=my_bsum_start(spin), my_bsum_stop(spin)
1325 : call wfd%rotate_cg(ib_sum, ndat1, spin, kqmp_ibz, npw_kqmp, kg_kqmp, istwf_kqmp, &
1326 7328 : cryst, mapl_kqmp, gbound_kqmp, work_ngfft, work, cg_kqmp)
1327 :
1328 : ! NB: cg_kqmp is dimensioned with mpw --> have to slice cg_kqmp
1329 27033572 : stern_kmp%cgq(:,:,ib_sum) = cg_kqmp(:,1:npw_kqmp*nspinor)
1330 : end do ! ib_sum
1331 :
1332 : ! ==========================
1333 : ! Delta_{-qka} psi_{n'k+q-p}
1334 : ! ==========================
1335 : ! (k-p, k+q-p)
1336 : call stern_kqmp%init(dtset, npw_kqmp, npw_kmp, nspinor, nbsum, nband_me, fermie1_idir_ipert, &
1337 696 : stern_use_cache, work_ngfft, mpi_enreg, stern_comm)
1338 :
1339 : ! Use symmetries to reconstruct GS wavefunctions at k-p and store them in stern_kqmp%cgq.
1340 8024 : do ib_sum=my_bsum_start(spin), my_bsum_stop(spin)
1341 : call wfd%rotate_cg(ib_sum, ndat1, spin, kmp_ibz, npw_kmp, kg_kmp, istwf_kmp, &
1342 7328 : cryst, mapl_kmp, gbound_kmp, work_ngfft, work, cg_kmp)
1343 :
1344 : ! NB: cg_kmp is dimensioned with mpw --> have to slice cg_kmp
1345 27033572 : stern_kqmp%cgq(:,:,ib_sum) = cg_kmp(:,1:npw_kmp*nspinor)
1346 : end do ! ib_sum
1347 :
1348 : ! Get PPM parameters at pp_bz to applying W_{gg'}(pp_bz).
1349 : ! Note: Sq_ibz = q_bz + G0 with non-zero G0 is not supported here
1350 696 : call ppm%get_qbz(gsph_c, pp_mesh, ipp_bz, botsq_pbz, otq_pbz, dmeig_pbz)
1351 :
1352 : ! Need transpose of PPM matrices when summing over the G index.
1353 696 : if (ppm%dm2_botsq /= 0) then
1354 696 : call alloc_copy(botsq_pbz, trans_botsq_pbz)
1355 696 : call sqmat_itranspose(npw_c, trans_botsq_pbz)
1356 : end if
1357 696 : if (ppm%dm2_otq /= 0) then
1358 696 : call alloc_copy(otq_pbz, trans_otq_pbz)
1359 696 : call sqmat_itranspose(npw_c, trans_otq_pbz)
1360 : end if
1361 696 : if (ppm%dm_eig /= 0) then
1362 696 : call alloc_copy(dmeig_pbz, trans_dmeig_pbz)
1363 696 : call sqmat_itranspose(ppm%dm_eig, trans_dmeig_pbz)
1364 : end if
1365 :
1366 : ! Find the corresponding irred pp-point in the pp_mesh.
1367 696 : call pp_mesh%get_bz_item(ipp_bz, pp, ipp_ibz, isym_pp, itim_pp)
1368 :
1369 : ! Get Fourier components of the Coulomb interaction in the BZ
1370 : ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2
1371 : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
1372 : ! NOTE: vc_sqrt_gx is dimensioned with npw_x --> use rottb table from gsph_x.
1373 396144 : do ig=1,npw_x
1374 396144 : vc_sqrt_gx(gsph_x%rottb(ig, itim_pp, isym_pp)) = vcp%vc_sqrt(ig, ipp_ibz)
1375 : end do
1376 :
1377 696 : if (dtset%gwcomp == 2) then
1378 : ! Compute static limit at pp_bz from the symmetrical image in the IBZ
1379 0 : call em1_symmetrize_op(ipp_bz, npw_c, 1, gsph_c, pp_mesh, wc0_pibz(ipp_ibz)%vals, wc0_pbz)
1380 : end if
1381 696 : call timab(1941, 2, tsec)
1382 :
1383 : ! ===========================================
1384 : ! Sum over bands (n' index in equations)
1385 : ! All procs in gqk%pert_comm enter this part.
1386 : ! ===========================================
1387 :
1388 8024 : do ib_sum=my_bsum_start(spin), my_bsum_stop(spin)
1389 7328 : call timab(1942, 1, tsec)
1390 : ! Get u_{n',k-p}(r)^*, stored in ur_star_kmp.
1391 : call wfd%rotate_cg(ib_sum, ndat1, spin, kmp_ibz, npw_kmp, kg_kmp, istwf_kmp, &
1392 7328 : cryst, mapl_kmp, gbound_kmp, work_ngfft, work, cg_kmp, urs_kbz=ur_star_kmp)
1393 :
1394 38911136 : ur_star_kmp = GWPC_CONJG(ur_star_kmp)
1395 7328 : ebsum_kmp = ebands%eig(ib_sum, ikmp_ibz, spin)
1396 7328 : ebsum_kqmp = ebands%eig(ib_sum, ikqmp_ibz, spin)
1397 :
1398 : ! =====================================
1399 : ! Precompute oscillator matrix elements
1400 : ! =====================================
1401 : ! These terms do not depend on (idir, ipert) and can be reused in the loop over perturbations below.
1402 : ! If the n' bands in the sum are distributed, one should transmit the (m, n) indices.
1403 :
1404 7328 : theta_mu_minus_e0i = fact_spin * qp_occ(ib_sum, ikmp_ibz, spin)
1405 7328 : need_x_kmp = (abs(theta_mu_minus_e0i / fact_spin) >= tol_empty) ! allow negative occ numbers
1406 : !need_x_kmp = .True.
1407 : !print *, "kmp, ib_sum, theta_mu_minus_e0i", ib_sum, theta_mu_minus_e0i
1408 :
1409 : ! Contract immediately over g' with the frequency convolution:
1410 : !
1411 : ! \int de' Wc_{gg'}(pp, e') / (omega - e_{bsum, kmp) - e')
1412 : !
1413 : ! vec_gwc_nk(:,:,n_k) stores:
1414 : !
1415 : ! sum_g' \int de' Wc_{gg'}(pp, e') / (omega - e_{bsum, kmp) - e') <bsum,k-p|e^{-i(p+g')}r|n,k>
1416 : !
1417 7328 : if (gqk%pert_comm%nproc > 1) then
1418 0 : vec_gwc_nk = zero
1419 0 : vec_gx_nk = zero
1420 0 : if (dtset%gwcomp == 2) vec_coh_nk = zero
1421 : end if
1422 :
1423 41056 : do n_k=bstart_k, bstop_k
1424 33728 : in_k = n_k - bstart_k + 1; if (gqk%pert_comm%skip(in_k)) cycle ! MPI parallelism inside pert_comm
1425 :
1426 : ! Compute <bsum,k-p|e^{-i(p+G')}r|n,k> * vc_sqrt(p,G')
1427 173770624 : cwork_ur = ur_star_kmp * ur_nk(:,n_k)
1428 :
1429 33728 : if (need_x_kmp) then
1430 20800 : call fft_ur(npw_x, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_x, gbound_x, cwork_ur, rhotwg_x)
1431 20800 : call sigtk_multiply_by_vc_sqrt("N", npw_x, nspinor, ndat1, vc_sqrt_gx, rhotwg_x)
1432 20800 : if (pp_is_gamma) then
1433 808 : rhotwg_x(1) = czero_gw; if (ib_sum == n_k) rhotwg_x(1) = cmplx(sqrt(vcp%i_sz),0.0_gwp)
1434 : end if
1435 :
1436 12579456 : vec_gx_nk(:,n_k) = rhotwg_x(1:npw_x*nspinor)
1437 20800 : if (dtset%userid /= 0) vec_gx_nk(:,n_k) = zero
1438 : ! FIXME: This is wrong if nspinor == 2
1439 553600 : rhotwg_c(:) = rhotwg_x(1:npw_c*nspinor)
1440 :
1441 : if (qq_is_gamma .and. test_sigma) then
1442 : sigx_nk(in_k, ik_glob) = sigx_nk(in_k, ik_glob) + dot_product(rhotwg_x, rhotwg_x)
1443 : end if
1444 :
1445 : else
1446 12928 : call fft_ur(npw_c, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_c, gbound_c, cwork_ur, rhotwg_c)
1447 12928 : call sigtk_multiply_by_vc_sqrt("N", npw_c, nspinor, ndat1, vc_sqrt_gx, rhotwg_c)
1448 12928 : if (pp_is_gamma) then
1449 1168 : rhotwg_c(1) = czero_gw; if (ib_sum == n_k) rhotwg_c(1) = cmplx(sqrt(vcp%i_sz),0.0_gwp)
1450 : end if
1451 : end if
1452 :
1453 : ! Prepare list of omegas: first e_nk then e_mkq for all m indices.
1454 33728 : omegas_nk(1) = qp_ene(n_k, ik_ibz, spin); cnt = 1
1455 :
1456 33728 : if (dtset%gwpt_wmode == 1) then
1457 172224 : do m_kq=bstart_kq, bstop_kq
1458 172224 : cnt = cnt + 1; omegas_nk(cnt) = qp_ene(m_kq, ikq_ibz, spin)
1459 : end do
1460 : end if
1461 260160 : omegame0i_nk = omegas_nk - qp_ene(ib_sum, ikmp_ibz, spin)
1462 :
1463 : ! Note that the i/two_pi factor in Sigma(w) is included in calc_sigc
1464 11728192 : vec_gwc_nk(:,:,n_k) = zero
1465 : call ppm%calc_sigc(nspinor, npw_c, nw_nk, rhotwg_c, botsq_pbz, otq_pbz, &
1466 : omegame0i_nk, dtset%zcut, theta_mu_minus_e0i, dmeig_pbz, npw_c, &
1467 33728 : vec_gwc_nk(:,:,n_k), sigcme_nk)
1468 :
1469 33728 : if (dtset%gwcomp == 2) then
1470 0 : vec_coh_nk(:, n_k) = matmul(wc0_pbz, rhotwg_c)
1471 : !call xgemv("N", npw_c, npw_c, cone_gw, wc0_pbz, npw_c, rhotwg_c, 1, czero_gw, vec_coh_nk(:, n_k), 1)
1472 : end if
1473 :
1474 : if (qq_is_gamma .and. test_sigma) then
1475 : sigce0_nk(in_k, ik_glob) = sigce0_nk(in_k, ik_glob) + dot_product(rhotwg_c, vec_gwc_nk(:,1,n_k))
1476 : end if
1477 :
1478 41056 : if (dtset%useric /= 0) vec_gwc_nk(:,:,n_k) = zero
1479 : end do ! n_k
1480 :
1481 : ! TODO: this is an all_gatherv but oh well.
1482 7328 : if (gqk%pert_comm%nproc > 1) then
1483 0 : call xmpi_sum(vec_gx_nk, gqk%pert_comm%value, ierr)
1484 0 : call xmpi_sum(vec_gwc_nk, gqk%pert_comm%value, ierr)
1485 0 : if (dtset%gwcomp == 2) call xmpi_sum(vec_coh_nk, gqk%pert_comm%value, ierr)
1486 : end if
1487 :
1488 7328 : call timab(1942, 2, tsec)
1489 7328 : call timab(1943, 1, tsec)
1490 :
1491 : ! Get u_{n',k+q-p}(r), stored in ur_star_kqmp
1492 : call wfd%rotate_cg(ib_sum, ndat1, spin, kqmp_ibz, npw_kqmp, kg_kqmp, istwf_kqmp, &
1493 7328 : cryst, mapl_kqmp, gbound_kqmp, work_ngfft, work, cg_kqmp, urs_kbz=ur_star_kqmp)
1494 :
1495 38911136 : ur_star_kqmp = GWPC_CONJG(ur_star_kqmp)
1496 7328 : theta_mu_minus_e0i = fact_spin * qp_occ(ib_sum, ikqmp_ibz, spin)
1497 :
1498 7328 : need_x_kqmp = (abs(theta_mu_minus_e0i / fact_spin) >= tol_empty) ! allow negative occ numbers
1499 : !print *, "kqmp, ib_sum, theta_mu_minus_e0i", ib_sum, theta_mu_minus_e0i
1500 : !need_x_kqmp = .True.
1501 :
1502 : ! Contract immediately over g with the frequency convolution:
1503 : !
1504 : ! \int de' Wc_{gg'}(pp, e') / (omega - e_{bsum, kqmp) - e')
1505 : !
1506 : ! Store results in vec_gwc_mkq(:,:,m_kq).
1507 : !
1508 : ! sum_g \int de' Wc_{gg'}(pp, e') / (omega - e_{bsum, kqmp) - e') <m,k+q|e^{+i(p+g)}r|bsum,k+q-p>
1509 : !
1510 7328 : if (gqk%pert_comm%nproc > 1) then
1511 0 : vec_gwc_mkq = zero
1512 0 : if (need_x_kqmp) vec_gx_mkq = zero
1513 0 : if (dtset%gwcomp == 2) vec_coh_mkq = zero
1514 : end if
1515 :
1516 85024 : do m_kq=bstart_kq, bstop_kq
1517 77696 : im_kq = m_kq - bstart_kq + 1; if (gqk%pert_comm%skip(im_kq)) cycle ! MPI parallelism inside pert_comm
1518 :
1519 : ! <m,k+q|e^{i(p+G)}r|bsum,k+q-p> * vc_sqrt(p,G).
1520 : ! Exchange bra and ket and take the CC of the FFT in sigtk_multiply_by_vc_sqrt.
1521 407281408 : cwork_ur = ur_star_kqmp * ur_mkq(:,m_kq)
1522 :
1523 77696 : if (need_x_kqmp) then
1524 49792 : call fft_ur(npw_x, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_x, gbound_x, cwork_ur, rhotwg_x)
1525 49792 : call sigtk_multiply_by_vc_sqrt("C", npw_x, nspinor, ndat1, vc_sqrt_gx, rhotwg_x)
1526 49792 : if (pp_is_gamma) then
1527 1744 : rhotwg_x(1) = czero_gw; if (ib_sum == m_kq) rhotwg_x(1) = cmplx(sqrt(vcp%i_sz),0.0_gwp)
1528 : end if
1529 :
1530 31302912 : vec_gx_mkq(:,m_kq) = rhotwg_x(1:npw_x*nspinor)
1531 49792 : if (dtset%userid /= 0) vec_gx_mkq(:,m_kq) = zero
1532 1238272 : rhotwg_c(:) = rhotwg_x(1:npw_c*nspinor)
1533 :
1534 : else
1535 27904 : call fft_ur(npw_c, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_c, gbound_c, cwork_ur, rhotwg_c)
1536 27904 : call sigtk_multiply_by_vc_sqrt("C", npw_c, nspinor, ndat1, vc_sqrt_gx, rhotwg_c)
1537 27904 : if (pp_is_gamma) then
1538 2368 : rhotwg_c(1) = czero_gw; if (ib_sum == m_kq) rhotwg_c(1) = cmplx(sqrt(vcp%i_sz),0.0_gwp)
1539 : end if
1540 : end if
1541 :
1542 77696 : if (dtset%gwpt_wmode == 1) then
1543 : ! Prepare list of omegas: first e_mkq then e_nk for all n indices.
1544 26496 : omegas_mkq(1) = qp_ene(m_kq, ikq_ibz, spin); cnt = 1
1545 185472 : do n_k=bstart_k, bstop_k
1546 185472 : cnt = cnt + 1; omegas_mkq(cnt) = qp_ene(n_k, ik_ibz, spin)
1547 : end do
1548 : else
1549 : ! Here we compute the convolution at e_nk for all n_k bands at fixed m_kq
1550 51200 : cnt = 0
1551 256000 : do n_k=bstart_k, bstop_k
1552 256000 : cnt = cnt + 1; omegas_mkq(cnt) = qp_ene(n_k, ik_ibz, spin)
1553 : end do
1554 : end if
1555 :
1556 545664 : omegame0i_mkq = omegas_mkq - qp_ene(ib_sum, ikqmp_ibz, spin)
1557 :
1558 : ! Here we sum over G instead of G' so we have to pass the transpose of the PPM matrix elements.
1559 : ! TODO: Generalize ppm%calc_sigc with BLAS-like API.
1560 15595648 : vec_gwc_mkq(:,:,m_kq) = zero
1561 : call ppm%calc_sigc(nspinor, npw_c, nw_mkq, rhotwg_c, trans_botsq_pbz, trans_otq_pbz, &
1562 : omegame0i_mkq, dtset%zcut, theta_mu_minus_e0i, trans_dmeig_pbz, npw_c, &
1563 77696 : vec_gwc_mkq(:,:,m_kq), sigcme_mkq)
1564 :
1565 77696 : if (dtset%gwcomp == 2) then
1566 0 : vec_coh_mkq(:, m_kq) = matmul(conjg(rhotwg_c), wc0_pbz)
1567 : !call xgemv("H", npw_c, npw_c, cone_gw, wc0_pbz, npw_c, rhotwg_c, 1, czero_gw, vec_coh_mkq(:, m_kq), 1)
1568 : end if
1569 :
1570 85024 : if (dtset%useric /= 0) vec_gwc_mkq(:,:,m_kq) = zero
1571 : end do ! m_kq
1572 :
1573 : ! TODO: this is an all_gatherv but oh well.
1574 7328 : if (gqk%pert_comm%nproc > 1) then
1575 0 : call xmpi_sum(vec_gwc_mkq, gqk%pert_comm%value, ierr)
1576 0 : if (need_x_kqmp) call xmpi_sum(vec_gx_mkq, gqk%pert_comm%value, ierr)
1577 0 : if (dtset%gwcomp == 2) call xmpi_sum(vec_coh_mkq, gqk%pert_comm%value, ierr)
1578 : end if
1579 :
1580 7328 : call timab(1943, 2, tsec)
1581 :
1582 : ! ========================================
1583 : ! Loop over my set of atomic perturbations
1584 : ! ========================================
1585 : ! For each perturbation:
1586 : ! - setup H1 from vlocal1_qq or vlocal1_mqq.
1587 : ! For each band in band_sum:
1588 : ! - Solve the Sternheimer non-self-consistently and get the KS e-ph matrix elements.
1589 : ! - Build the full first-order wavefunction including the active subspace.
1590 :
1591 51296 : do imyp=1,gqk%my_npert
1592 43968 : call timab(1944, 1, tsec)
1593 : ! NB: Only one proc enters this section. No MPI parallelism is allowed here.
1594 43968 : idir = dvdb%my_pinfo(1, imyp); ipert = dvdb%my_pinfo(2, imyp); ipc = dvdb%my_pinfo(3, imyp)
1595 : !print *, "For kk, ", kk, "pp:", pp, "idir, ipert", idir, ipert
1596 :
1597 : ! Set up local potential vlocal1_qq with proper dimensioning, from vtrial1 taking into account the spin
1598 : ! and prepare application of the NL part. Each MPI rank prepares its own potential.
1599 : call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, nvloc, &
1600 43968 : pawfgr, mpi_enreg, vtrial, v1scf_qq(:,:,:,imyp), vlocal, vlocal1_qq(:,:,:,:,imyp))
1601 :
1602 43968 : call gs_ham_kqmp%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.)
1603 :
1604 43968 : call rf_ham_kqmp%init(cplex, gs_ham_kqmp, ipert, has_e1kbsc=.true.)
1605 43968 : call rf_ham_kqmp%load_spin(spin, vlocal1=vlocal1_qq(:,:,:,:,imyp), with_nonlocal=.true.)
1606 :
1607 : ! ======================
1608 : ! NSCF Sternheimer at qq
1609 : ! ======================
1610 : ! Compute Delta_{q,idir,ipert} \psi_{bsum, k-p}
1611 : call wfd%rotate_cg(ib_sum, ndat1, spin, kmp_ibz, npw_kmp, kg_kmp, istwf_kmp, &
1612 43968 : cryst, mapl_kmp, gbound_kmp, work_ngfft, work, cg_kmp)
1613 :
1614 510144 : stern_kmp%bands_treated_now(:) = 0; stern_kmp%bands_treated_now(ib_sum) = 1
1615 510144 : stern_kmp%rank_band = 0; u1_band = ib_sum; band_me = ib_sum
1616 :
1617 : ! (k+q-p, k-p)
1618 43968 : init_mode = "None"
1619 : call stern_kmp%solve(u1_band, band_me, idir, ipert, qq_bz, gs_ham_kqmp, rf_ham_kqmp, &
1620 : ebands%eig(:,ikmp_ibz,spin), ebands%eig(:,ikqmp_ibz,spin), &
1621 : cg_kmp, cwaveprj0, cg1_kqmp, cwaveprj, msg, ierr, &
1622 43968 : full_cg1=full_cg1_kqmp, full_ur1=full_ur1_kqmp, init_mode=init_mode)
1623 :
1624 : ! The last bands may fail to converge with resid=-two. In this case we ignore the contribution.
1625 43968 : if (ierr /= 0) then
1626 1446230 : full_cg1_kqmp = zero; full_ur1_kqmp = zero; stern_qq_ierr = stern_qq_ierr + 1
1627 182 : if (stern_qq_ierr <= ENOUGH_STERN) then
1628 15 : ABI_WARNING(sjoin("Stern at +q", qkp_string, msg))
1629 : end if
1630 : end if
1631 :
1632 : ! Store KS e-ph matrix elements for this ipc perturbation.
1633 43968 : if (pp_is_gamma .and. (ib_sum >= gqk%bstart_k .and. ib_sum <= gqk%bstop_k)) then
1634 1020 : ii = ib_sum - bstart_k + 1
1635 36588 : gks_atm(:,:,ii,ipc) = stern_kmp%eig1_k(:, bstart_kq:bstop_kq, ib_sum)
1636 : end if
1637 :
1638 : ! Compute <m,k+q|e^{i(p+G)r}|Delta_q psi_{bsum,k-p}>
1639 : ! Exchange bra and ket and take the CC of the FFT.
1640 233466816 : full_ur1_kqmp = GWPC_CONJG(full_ur1_kqmp)
1641 :
1642 510144 : do m_kq=bstart_kq, bstop_kq
1643 466176 : im_kq = m_kq - bstart_kq + 1
1644 2443688448 : cwork_ur = full_ur1_kqmp * ur_mkq(:,m_kq)
1645 :
1646 466176 : ctmp_dp = r2c(stern_kmp%eig1_k(:, m_kq, ib_sum))
1647 466176 : e_mkq = ebands%eig(m_kq, ikq_ibz, spin)
1648 466176 : delta_E = ebsum_kmp - e_mkq
1649 466176 : inv_delta_E = delta_E / ( delta_E ** 2 + eta_sum ** 2)
1650 466176 : ctmp_dp = ctmp_dp * inv_delta_E * sqrt(vcp%i_sz)
1651 : !ctmp_dp = zero
1652 : !if (qq_is_gamma) ctmp_dp = zero
1653 :
1654 466176 : if (need_x_kqmp) then
1655 298752 : call fft_ur(npw_x, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_x, gbound_x, cwork_ur, rhotwg_x)
1656 : !if (pp_is_gamma) print *, "rhotwg_x(1)", rhotwg_x(1), "vc(0)", vc_sqrt_gx(1)
1657 298752 : ctmp_gwpc = rhotwg_x(1)
1658 298752 : call sigtk_multiply_by_vc_sqrt("C", npw_x, nspinor, ndat1, vc_sqrt_gx, rhotwg_x)
1659 298752 : if (pp_is_gamma) then
1660 10464 : rhotwg_x(1) = ctmp_gwpc * cmplx(sqrt(vcp%i_sz),0.0_gwp)
1661 : !print *, "ctmp_dp +qq", ctmp_gwpc
1662 : !print *, "rhotwg_x(1)", rhotwg_x(1)
1663 : !rhotwg_x(1) = ctmp_dp
1664 : end if
1665 7429632 : rhotwg_c(:) = rhotwg_x(1:npw_c*nspinor)
1666 : else
1667 167424 : call fft_ur(npw_c, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_c, gbound_c, cwork_ur, rhotwg_c)
1668 167424 : ctmp_gwpc = rhotwg_c(1)
1669 167424 : call sigtk_multiply_by_vc_sqrt("C", npw_c, nspinor, ndat1, vc_sqrt_gx, rhotwg_c)
1670 167424 : if (pp_is_gamma) then
1671 14208 : rhotwg_c(1) = ctmp_gwpc * cmplx(sqrt(vcp%i_sz),0.0_gwp)
1672 : !rhotwg_c(1) = ctmp_dp
1673 : end if
1674 : end if
1675 :
1676 2692800 : do n_k=bstart_k, bstop_k
1677 2182656 : in_k = n_k - bstart_k + 1
1678 :
1679 2182656 : if (dtset%gwpt_wmode == 1) then
1680 : ! +2 because omegas_nk(1) = e_nk and omegas_nk(2:cnt) = e_mkq for all m_kq bands
1681 953856 : iw_mkq = m_kq - bstart_kq + 2
1682 : ! Take the average at e_nk and e_mkq
1683 62954496 : ctmp_gwpc = half * sum(rhotwg_c(:) * (vec_gwc_nk(:,1,n_k) + vec_gwc_nk(:,iw_mkq,n_k)))
1684 : else
1685 : ! Use the value at e_nk
1686 19660800 : ctmp_gwpc = sum(rhotwg_c(:) * (vec_gwc_nk(:,1,n_k)))
1687 : end if
1688 :
1689 2182656 : if (need_x_kqmp) then
1690 1300992 : xdot_tmp = - xdotu(npw_x*nspinor, rhotwg_x, 1, vec_gx_nk(:,n_k), 1)
1691 1300992 : ctmp_gwpc = ctmp_gwpc + xdot_tmp ! * theta_mu_minus_e0i ! theta_mu_minus_e0i is only needed for metals
1692 : end if
1693 :
1694 2182656 : if (dtset%gwcomp == 2) then
1695 0 : ctmp_gwpc = ctmp_gwpc - quarter * dot_product(conjg(rhotwg_c), vec_coh_nk(:, n_k))
1696 : end if
1697 :
1698 2648832 : if (dtset%userie < 0) then
1699 0 : if (pp_is_gamma) then
1700 : gsig_atm(1, im_kq, in_k, ipc) = gsig_atm(1, im_kq, in_k, ipc)
1701 : gsig_atm(2, im_kq, in_k, ipc) = gsig_atm(2, im_kq, in_k, ipc)
1702 : else
1703 0 : gsig_atm(1, im_kq, in_k, ipc) = gsig_atm(1, im_kq, in_k, ipc) + real(ctmp_gwpc)
1704 0 : gsig_atm(2, im_kq, in_k, ipc) = gsig_atm(2, im_kq, in_k, ipc) + aimag(ctmp_gwpc)
1705 : end if
1706 : else
1707 2182656 : gsig_atm(1, im_kq, in_k, ipc) = gsig_atm(1, im_kq, in_k, ipc) + real(ctmp_gwpc)
1708 2182656 : gsig_atm(2, im_kq, in_k, ipc) = gsig_atm(2, im_kq, in_k, ipc) + aimag(ctmp_gwpc)
1709 : end if
1710 :
1711 : end do ! n_k
1712 : end do ! m_kq
1713 :
1714 : !if (.not. qq_is_gamma) then
1715 : ! ==========================
1716 : ! Same operations but for -q
1717 : ! ==========================
1718 : ! Set up local potential vlocal1_qq with proper dimensioning, from vtrial1 taking into account the spin
1719 : ! and prepare application of the NL part. Each MPI rank prepares its own potential.
1720 : call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, nvloc, &
1721 43968 : pawfgr, mpi_enreg, vtrial, v1scf_mq(:,:,:,imyp), vlocal, vlocal1_mqq(:,:,:,:,imyp))
1722 :
1723 43968 : call gs_ham_kmp%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.)
1724 :
1725 43968 : call rf_ham_kmp%init(cplex, gs_ham_kmp, ipert, has_e1kbsc=.true.)
1726 43968 : call rf_ham_kmp%load_spin(spin, vlocal1=vlocal1_mqq(:,:,:,:,imyp), with_nonlocal=.true.)
1727 :
1728 : ! ======================
1729 : ! NSCF Sternheimer at -q
1730 : ! ======================
1731 : ! Compute Delta_{-q,idir,ipert} \psi_{bsum, k+q-p}
1732 : call wfd%rotate_cg(ib_sum, ndat1, spin, kqmp_ibz, npw_kqmp, kg_kqmp, istwf_kqmp, &
1733 43968 : cryst, mapl_kqmp, gbound_kqmp, work_ngfft, work, cg_kqmp)
1734 :
1735 510144 : stern_kqmp%bands_treated_now(:) = 0; stern_kqmp%bands_treated_now(ib_sum) = 1
1736 510144 : stern_kqmp%rank_band = 0; u1_band = ib_sum; band_me = ib_sum
1737 :
1738 : !init_mode = "input"
1739 43968 : init_mode = "None"
1740 43968 : if (init_mode == "input") then
1741 : call cgtk_change_gsphere(nspinor, &
1742 : npw_kqmp, istwfk1, kg_kqmp, cg1_kqmp, &
1743 0 : npw_kmp, istwfk1, kg_kmp, cg1_kmp, work_ngfft, work)
1744 : end if
1745 :
1746 : ! (k-p, k+q-p)
1747 : call stern_kqmp%solve(u1_band, band_me, idir, ipert, -qq_bz, gs_ham_kmp, rf_ham_kmp, &
1748 : ebands%eig(:,ikqmp_ibz,spin), ebands%eig(:,ikmp_ibz,spin), &
1749 : cg_kqmp, cwaveprj0, cg1_kmp, cwaveprj, msg, ierr, &
1750 175872 : full_cg1=full_cg1_kmp, full_ur1=full_ur1_star_kmp, init_mode=init_mode)
1751 :
1752 233466816 : full_ur1_star_kmp = GWPC_CONJG(full_ur1_star_kmp)
1753 :
1754 : ! The last bands may fail to converge with resid=-two. In this case we ignore the contribution.
1755 43968 : if (ierr /= 0) then
1756 1446230 : full_cg1_kmp = zero; full_ur1_star_kmp = zero; stern_mq_ierr = stern_mq_ierr + 1
1757 182 : if (stern_mq_ierr <= ENOUGH_STERN) then
1758 15 : ABI_WARNING(sjoin("Stern at -q:", qkp_string, msg))
1759 : end if
1760 : end if
1761 :
1762 : ! For debug, gks_atm2 and gks_atm should be consistent
1763 43968 : if (pp_is_gamma .and. (ib_sum >= gqk%bstart_k .and. ib_sum <= gqk%bstop_k)) then
1764 1020 : in_k = ib_sum - bstart_k + 1
1765 36588 : gks_atm2(:,:,in_k,ipc) = stern_kqmp%eig1_k(:, bstart_kq:bstop_kq, ib_sum)
1766 : end if
1767 :
1768 43968 : call timab(1944, 2, tsec)
1769 43968 : call timab(1945, 1, tsec)
1770 246336 : do n_k=bstart_k, bstop_k
1771 202368 : in_k = n_k - bstart_k + 1
1772 :
1773 : ! <Delta_{-q} psi_{bsum,k+q-p}|e^{-i(p+G')r}|n,k>
1774 1042623744 : cwork_ur = full_ur1_star_kmp * ur_nk(:,n_k)
1775 :
1776 : ! TODO: Recheck equations, may need a complex conjugate
1777 202368 : ctmp_dp = r2c(stern_kqmp%eig1_k(:, n_k, ib_sum))
1778 202368 : e_nk = ebands%eig(n_k, ik_ibz, spin)
1779 202368 : delta_E = ebsum_kqmp - e_nk
1780 202368 : inv_delta_E = delta_E / ( delta_E ** 2 + eta_sum ** 2)
1781 202368 : ctmp_dp = ctmp_dp * inv_delta_E * sqrt(vcp%i_sz)
1782 : !ctmp_dp = zero
1783 : !if (qq_is_gamma) ctmp_dp = zero
1784 :
1785 202368 : if (need_x_kmp) then
1786 124800 : call fft_ur(npw_x, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_x, gbound_x, cwork_ur, rhotwg_x)
1787 124800 : ctmp_gwpc = rhotwg_x(1)
1788 124800 : call sigtk_multiply_by_vc_sqrt("N", npw_x, nspinor, ndat1, vc_sqrt_gx, rhotwg_x)
1789 124800 : if (pp_is_gamma) then
1790 4848 : rhotwg_x(1) = ctmp_gwpc * cmplx(sqrt(vcp%i_sz),0.0_gwp)
1791 : !print *, "ctmp_dp -qq", ctmp_dp
1792 : !print *, "rhotwg_x(1)", rhotwg_x(1)
1793 : !rhotwg_x(1) = ctmp_dp
1794 : end if
1795 :
1796 3321600 : rhotwg_c(:) = rhotwg_x(1:npw_c*nspinor)
1797 : else
1798 77568 : call fft_ur(npw_c, nfft, nspinor, ndat1, mgfft, ngfft, istwfk1, kg_c, gbound_c, cwork_ur, rhotwg_c)
1799 77568 : ctmp_gwpc = rhotwg_c(1)
1800 77568 : call sigtk_multiply_by_vc_sqrt("N", npw_c, nspinor, ndat1, vc_sqrt_gx, rhotwg_c)
1801 77568 : if (pp_is_gamma) then
1802 7008 : rhotwg_c(1) = ctmp_gwpc * cmplx(sqrt(vcp%i_sz),0.0_gwp)
1803 : !rhotwg_c(1) = ctmp_dp
1804 : end if
1805 : end if
1806 :
1807 2428992 : do m_kq=bstart_kq, bstop_kq
1808 2182656 : im_kq = m_kq - bstart_kq + 1
1809 :
1810 2182656 : if (dtset%gwpt_wmode == 1) then
1811 : ! +2 because omegas_mkq(1) = e_mkq and omegas_mkq(2:cnt) = e_nk for all n_k bands
1812 953856 : iw_nk = n_k - bstart_k + 2
1813 : ! Take the average at e_nk and e_mkq
1814 62954496 : ctmp_gwpc = half * sum(rhotwg_c(:) * (vec_gwc_mkq(:,1,m_kq) + vec_gwc_mkq(:,iw_nk,m_kq)))
1815 : else
1816 : ! Use the value at e_nk. Note in_k index
1817 19660800 : ctmp_gwpc = sum(rhotwg_c(:) * vec_gwc_mkq(:,in_k,m_kq))
1818 : end if
1819 :
1820 2182656 : if (need_x_kmp) then
1821 : ! TODO recheck
1822 1300992 : xdot_tmp = - xdotu(npw_x*nspinor, vec_gx_mkq(:,m_kq), 1, rhotwg_x, 1)
1823 1300992 : ctmp_gwpc = ctmp_gwpc + xdot_tmp ! * theta_mu_minus_e0i ! theta_mu_minus_e0i is only needed for metals
1824 : end if
1825 :
1826 2182656 : if (dtset%gwcomp == 2) then
1827 0 : ctmp_gwpc = ctmp_gwpc - quarter * dot_product(conjg(rhotwg_c), vec_coh_mkq(:, m_kq))
1828 : end if
1829 :
1830 2385024 : if (dtset%userie < 0) then
1831 0 : if (pp_is_gamma) then
1832 : gsig_atm(1, im_kq, in_k, ipc) = gsig_atm(1, im_kq, in_k, ipc)
1833 : gsig_atm(2, im_kq, in_k, ipc) = gsig_atm(2, im_kq, in_k, ipc)
1834 : else
1835 0 : gsig_atm(1, im_kq, in_k, ipc) = gsig_atm(1, im_kq, in_k, ipc) + real(ctmp_gwpc)
1836 0 : gsig_atm(2, im_kq, in_k, ipc) = gsig_atm(2, im_kq, in_k, ipc) + aimag(ctmp_gwpc)
1837 : end if
1838 : else
1839 2182656 : gsig_atm(1, im_kq, in_k, ipc) = gsig_atm(1, im_kq, in_k, ipc) + real(ctmp_gwpc)
1840 2182656 : gsig_atm(2, im_kq, in_k, ipc) = gsig_atm(2, im_kq, in_k, ipc) + aimag(ctmp_gwpc)
1841 : end if
1842 :
1843 : end do ! m_kq
1844 : end do ! n_k
1845 51296 : call timab(1945, 2, tsec)
1846 : !end if ! .not qq_is_gamma.
1847 : end do ! imyp (my perturbations)
1848 :
1849 8024 : call rf_ham_kqmp%free(); call rf_ham_kmp%free()
1850 : end do ! ib_sum (sum over bands)
1851 :
1852 : ! Show progress of all pp points.
1853 696 : if (print_time_pp) then
1854 267 : call inds2str(1, " My pp-point:", my_ipp, my_npp(spin), pp_mesh%nbz, msg)
1855 267 : call cwtime_report(msg, cpu_pp, wall_pp, gflops_pp);
1856 267 : if (my_ipp == LOG_MODP) call wrtout(std_out, "...", do_flush=.True.)
1857 : end if
1858 :
1859 696 : ABI_FREE(kpg_kmp)
1860 696 : ABI_FREE(kpg_kqmp)
1861 696 : ABI_FREE(ffnl_kmp)
1862 696 : ABI_FREE(ffnl_kqmp)
1863 696 : ABI_FREE(full_cg1_kqmp)
1864 696 : ABI_FREE(full_cg1_kmp)
1865 696 : ABI_FREE(rhotwg_c)
1866 696 : ABI_FREE(rhotwg_x)
1867 696 : ABI_FREE(vc_sqrt_gx)
1868 696 : ABI_FREE(kinpw_kmp)
1869 696 : ABI_FREE(kinpw_kqmp)
1870 696 : ABI_FREE(ph3d_kmp)
1871 696 : ABI_FREE(ph3d_kqmp)
1872 :
1873 696 : ABI_SFREE(ylm_kmp)
1874 696 : ABI_SFREE(ylm_kqmp)
1875 696 : ABI_SFREE(botsq_pbz)
1876 696 : ABI_SFREE(otq_pbz)
1877 696 : ABI_SFREE(dmeig_pbz)
1878 696 : ABI_SFREE(trans_botsq_pbz)
1879 696 : ABI_SFREE(trans_otq_pbz)
1880 696 : ABI_SFREE(trans_dmeig_pbz)
1881 :
1882 727 : call stern_kmp%free(); call stern_kqmp%free()
1883 : end do ! ipp_bz
1884 :
1885 31 : ABI_FREE(kpg_k)
1886 31 : ABI_FREE(kpg_kq)
1887 31 : ABI_FREE(ug_k)
1888 31 : ABI_FREE(ug_kq)
1889 :
1890 : ! Here we are outside of the loop over pp_sum, band_sum and perturbations.
1891 : ! Collect gsig_atm and gks_atm inside pert_ppsum_comm so that all procs can write the data.
1892 : !call xmpi_sum_master(gsig_atm, master, gqk%pert_ppsum_bsum_comm%value, ierr)
1893 : !call xmpi_sum_master(gks_atm , master, gqk%pert_ppsum_bsum_comm%value, ierr)
1894 : !call xmpi_sum_master(gks_atm2, master, gqk%pert_ppsum_bsum_comm%value, ierr)
1895 : ! SC: not sure why, but xmpi_sum_master sometimes causes gvals to be filled
1896 : ! with zeros and none in the GSTORE.nc file when running in parallel (observed
1897 : ! on my desktop, lemaitre4, and lucia).
1898 31 : call xmpi_sum(gsig_atm, gqk%pert_ppsum_bsum_comm%value, ierr)
1899 31 : call xmpi_sum(gks_atm, gqk%pert_ppsum_bsum_comm%value, ierr)
1900 31 : call xmpi_sum(gks_atm2, gqk%pert_ppsum_bsum_comm%value, ierr)
1901 :
1902 36805 : gsig_atm = gsig_atm / (cryst%ucvol * pp_mesh%nbz)
1903 :
1904 31 : if (dtset%useria == 0) then
1905 36836 : gsig_atm = gsig_atm + gks_atm - gxc_atm
1906 0 : else if (dtset%useria > 0) then
1907 0 : gsig_atm = gxc_atm
1908 : else if (dtset%useria < 0) then
1909 : gsig_atm = gsig_atm
1910 : end if
1911 :
1912 : ! Save e-ph matrix elements in the buffer.
1913 36805 : my_gbuf(:,:,:,:, my_ik, iqbuf_cnt) = gsig_atm
1914 36805 : my_gbuf_ks(:,:,:,:, my_ik, iqbuf_cnt) = gks_atm
1915 31 : state_kq(my_ik, iqbuf_cnt) = GSTORE_KQ_COMPUTED
1916 :
1917 173 : if (print_time_kk) then
1918 31 : call inds2str(3, "My k-point", my_ik, gqk%my_nk, gqk%glob_nk, msg)
1919 31 : call cwtime_report(msg, cpu_kk, wall_kk, gflops_kk); if (my_ik == LOG_MODK) call wrtout(std_out, "...", do_flush=.True.)
1920 : end if
1921 : end do ! my_ik
1922 :
1923 : ! Dump buffer inside the loop over my_iq (and outside of my_ik) otherwise when restarting GWPT some my_ik points will be missing
1924 80 : if (iqbuf_cnt == qbuf_size) call dump_my_gbuf()
1925 :
1926 80 : ABI_FREE(v1scf_qq)
1927 80 : ABI_FREE(vlocal1_qq)
1928 80 : ABI_FREE(vlocal1_mqq)
1929 80 : ABI_FREE(vxc1_qq)
1930 :
1931 80 : if (print_time_qq .and. wall_qq > tol3) then
1932 80 : call inds2str(2, "My q-point", my_iq, gqk%my_nq, gqk%glob_nq, msg)
1933 80 : call cwtime_report(msg, cpu_qq, wall_qq, gflops_qq); if (my_iq == LOG_MODQ) call wrtout(std_out, "...", do_flush=.True.)
1934 : end if
1935 83 : call lg_myq%free()
1936 : end do ! my_iq
1937 :
1938 : ! Dump the remainder.
1939 3 : if (iqbuf_cnt /= 0) call dump_my_gbuf()
1940 :
1941 3 : ABI_FREE(ur_nk)
1942 3 : ABI_FREE(ur_mkq)
1943 3 : ABI_FREE(iq_buf)
1944 3 : ABI_FREE(state_kq)
1945 3 : ABI_FREE(my_gbuf)
1946 3 : ABI_FREE(my_gbuf_ks)
1947 3 : ABI_FREE(gsig_atm)
1948 3 : ABI_FREE(gks_atm)
1949 3 : ABI_FREE(gks_atm2)
1950 3 : ABI_FREE(gxc_atm)
1951 3 : ABI_FREE(omegas_nk)
1952 3 : ABI_FREE(omegas_mkq)
1953 3 : ABI_FREE(omegame0i_nk)
1954 3 : ABI_FREE(omegame0i_mkq)
1955 3 : ABI_FREE(vec_gwc_nk)
1956 3 : ABI_FREE(vec_gwc_mkq)
1957 3 : ABI_FREE(vec_gx_nk)
1958 3 : ABI_FREE(vec_gx_mkq)
1959 3 : ABI_FREE(sigcme_nk)
1960 3 : ABI_FREE(sigcme_mkq)
1961 3 : ABI_SFREE(vec_coh_nk)
1962 3 : ABI_SFREE(vec_coh_mkq)
1963 :
1964 : ! Collect self-energy matrix elements.
1965 : if (test_sigma) then
1966 : call xmpi_sum(vxc_nk, gqk%kpt_comm%value, ierr)
1967 : call xmpi_sum(sigx_nk, gqk%kpt_comm%value, ierr)
1968 : call xmpi_sum(sigce0_nk, gqk%kpt_comm%value, ierr)
1969 : sigx_nk = -sigx_nk * (one / (cryst%ucvol * pp_mesh%nbz))
1970 : sigce0_nk = sigce0_nk * (one / (cryst%ucvol * pp_mesh%nbz))
1971 :
1972 : if (gqk%comm%me == master) then
1973 : write(ab_out, "(2a)") ch10, " Sigma^x_nk and Sigma^c_nk(E0) in eV:"
1974 : do ik_glob=1, gqk%glob_nk
1975 : ik_bz = gstore%kglob2bz(ik_glob, spin)
1976 : ik_ibz = gstore%kbz2ibz(1, ik_bz)
1977 : write(ab_out, "(2a)") "Band E0 <VxcDFT> SigX SigC(E0) for k-point:", trim(ktoa(gstore%kbz(:, ik_bz)))
1978 : do band=gqk%bstart_k, gqk%bstop_k
1979 : in_k = band - gqk%bstart_k + 1
1980 : write(ab_out, "(i5, 4(f8.3))") &
1981 : band, ebands%eig(band, ik_ibz, spin) * Ha_eV, &
1982 : vxc_nk(in_k, ik_glob) * Ha_eV, &
1983 : real(sigx_nk(in_k, ik_glob)) * Ha_eV, &
1984 : real(sigce0_nk(in_k, ik_glob)) * Ha_eV
1985 : end do
1986 : end do
1987 : end if
1988 :
1989 : ABI_SFREE(vxc_nk)
1990 : ABI_SFREE(sigx_nk)
1991 : ABI_SFREE(sigce0_nk)
1992 : end if ! test_sigma
1993 :
1994 6 : if (use_lgk) then
1995 8 : do my_ik=1,gqk%my_nk
1996 8 : call lg_myk(my_ik)%free()
1997 : end do
1998 8 : ABI_FREE(lg_myk)
1999 : end if
2000 : end do ! my_is
2001 :
2002 3 : call xmpi_sum(stern_qq_ierr, comm, ierr)
2003 3 : call xmpi_sum(stern_mq_ierr, comm, ierr)
2004 3 : call wrtout(units, sjoin(" Total number of failures in Sternheimer solvers at +q:", itoa(stern_qq_ierr)))
2005 3 : call wrtout(units, sjoin(" Total number of failures in Sternheimer solvers at -q:", itoa(stern_mq_ierr)))
2006 :
2007 3 : call cwtime_report(" gwpt_eph full calculation", cpu_all, wall_all, gflops_all, end_str=ch10)
2008 :
2009 : ! Set gstore_completed to 1 so that we can easily check if restarted is needed.
2010 : !if (my_rank == master) then
2011 3 : NCF_CHECK(nf90_put_var(root_ncid, root_vid("gstore_completed"), 1))
2012 : !end if
2013 : ! nf90_sync can deadlock on lemaitre4 for some MPI decompositions.
2014 3 : if (dtset%useria /= 888) then
2015 3 : NCF_CHECK(nf90_sync(root_ncid))
2016 : end if
2017 3 : NCF_CHECK(nf90_close(root_ncid))
2018 3 : call xmpi_barrier(comm)
2019 :
2020 : ! Reconstruct g(k,q) and g_KS(k,q) matrix elements in the full BZ by symmetry.
2021 3 : if (symmetrize) then
2022 0 : call gstore_symmetrize(gstore%path, wfk0_path, ngfft, dtset, dtfil, cryst, psps, pawtab, ebands, ifc, comm)
2023 : end if
2024 :
2025 : ! Output some of the results to ab_out for testing purposes
2026 3 : call gstore%print_for_abitests(dtset, ebands, .True., with_ks=.True.)
2027 :
2028 : ! Free memory
2029 3 : ABI_FREE(kg_k)
2030 3 : ABI_FREE(kg_kq)
2031 3 : ABI_FREE(kg_kmp)
2032 3 : ABI_FREE(kg_kqmp)
2033 3 : ABI_FREE(cg_work)
2034 3 : ABI_FREE(ur_star_kmp)
2035 3 : ABI_FREE(ur_star_kqmp)
2036 3 : ABI_FREE(full_ur1_kqmp)
2037 3 : ABI_FREE(full_ur1_star_kmp)
2038 3 : ABI_FREE(cwork_ur)
2039 3 : ABI_FREE(cg_kmp)
2040 3 : ABI_FREE(cg_kqmp)
2041 3 : ABI_FREE(cg1_kqmp)
2042 3 : ABI_FREE(cg1_kmp)
2043 3 : ABI_FREE(vtrial)
2044 3 : ABI_FREE(work)
2045 3 : ABI_FREE(ph1d)
2046 3 : ABI_FREE(vlocal)
2047 3 : ABI_FREE(gbound_k)
2048 3 : ABI_FREE(gbound_kq)
2049 3 : ABI_FREE(gbound_kmp)
2050 3 : ABI_FREE(gbound_kqmp)
2051 3 : ABI_FREE(gbound_c)
2052 3 : ABI_FREE(gbound_x)
2053 3 : ABI_FREE(done_qbz_spin)
2054 3 : ABI_FREE(rhor)
2055 3 : ABI_FREE(kxc)
2056 3 : ABI_FREE(vxc)
2057 :
2058 3 : if (dtset%gwcomp == 2) then
2059 0 : ABI_FREE(wc0_pbz)
2060 0 : do ip_ibz=1,pp_mesh%nibz
2061 0 : call wc0_pibz(ip_ibz)%free()
2062 : end do
2063 0 : ABI_FREE(wc0_pibz)
2064 : end if
2065 :
2066 3 : call gs_ham_kqmp%free(); call gs_ham_kmp%free(); call wfd%free(); call vcp%free(); call ppm%free()
2067 3 : call pp_mesh%free(); call gsph_c%free(); call gsph_x%free(); call gstore%free()
2068 3 : call pawcprj_free(cwaveprj0); call pawcprj_free(cwaveprj)
2069 3 : ABI_FREE(cwaveprj0)
2070 3 : ABI_FREE(cwaveprj)
2071 :
2072 3 : call xmpi_barrier(comm) ! This to make sure that the parallel output of GSTORE is completed
2073 27 : call cwtime_report(" gwpt_run: MPI barrier before returning.", cpu_all, wall_all, gflops_all, end_str=ch10, comm=comm)
2074 :
2075 : contains
2076 :
2077 489 : subroutine inds2str(level, prefix, my_ik, my_nk, nk_tot, out_str)
2078 : integer,intent(in) :: level, my_ik, my_nk, nk_tot
2079 : character(len=*),intent(in) :: prefix
2080 : character(len=*),intent(out) :: out_str
2081 :
2082 489 : out_str = sjoin(prefix, itoa(my_ik), "/", itoa(my_nk), "[", itoa(nk_tot), "]")
2083 2569 : out_str = repeat(' ', 4 * level) // trim(out_str)
2084 489 : end subroutine inds2str
2085 :
2086 83 : subroutine dump_my_gbuf()
2087 :
2088 : ! This function is called inside the double loop over (my_is, my_iq) or when we exit
2089 : ! from the my_iq loop to dump the remainder that is still in the q-buffer,
2090 : ! All the MPI procs in the (kpt_comm x pert_comm) grid shall call this contained routine
2091 : ! as we have side-effects i.e. iqbuf_cnt set to 0.
2092 :
2093 : ! On disk we have the global arrays:
2094 : !
2095 : ! nctkarr_t("gvals", "dp", "two, nb_kq, nb_k, natom3, glob_nk, glob_nq")
2096 : !
2097 : ! while the local MPI buffers are dimensioned as follows:
2098 : !
2099 : ! my_gbuf(2, nb_kq, nb_k, natom3, gqk%my_nk, qbuf_size)
2100 :
2101 : ! If parallelism over perturbation is activated, only the procs treating the first perturbation
2102 : ! i.e. the procs treating different k-points for this q are involved in IO
2103 : ! as all the local buffers store results for all natom3 perturbations.
2104 :
2105 : ! NOTE: A similar routine is used in m_gstore. The two implementations should be kept in synch.
2106 :
2107 : integer :: ii, iq_bz, iq_glob, my_iq
2108 :
2109 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2110 : ! FIXME: Recheck this part as we have way more levels of parallelism in GWPT
2111 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2112 : !if (gqk%coords_qkpb_sumbp(3) /= 0) goto 10 ! Yes, I'm very proud of this GOTO.
2113 : !if (gqk%pert_ppsum_bsum_comm%me /= 0) goto 10 ! Yes, I'm very proud of this GOTO.
2114 : ! SC: I comment out the above two GOTOs to avoid the deadlock issue on my desktop, lemaitre4 and lucia
2115 : !return
2116 :
2117 : !iq_buf(:, iqbuf_cnt) = [my_iq, iq_bz]
2118 83 : my_iq = iq_buf(1, 1)
2119 83 : iq_glob = my_iq + gqk%my_qstart - 1
2120 :
2121 : !print *, "in dump_my_gbuf with start: ", [1, 1, 1, 1, gqk%my_kstart, iq_glob]
2122 : !print *, " count; ", [2, gqk%nb_kq, gqk%nb_k, gqk%natom3, gqk%my_nk, iqbuf_cnt]
2123 : !print *, "my_gbuf", my_gbuf(:,:,:,natom3,1,1)
2124 :
2125 : ! Output g^Sigma
2126 : ! NB: this is an individual IO operation
2127 : ncerr = nf90_put_var(spin_ncid, spin_vid("gvals"), my_gbuf, &
2128 : start=[1, 1, 1, 1, gqk%my_kstart, iq_glob], &
2129 1079 : count=[2, gqk%nb_kq, gqk%nb_k, gqk%natom3, gqk%my_nk, iqbuf_cnt])
2130 83 : NCF_CHECK(ncerr)
2131 :
2132 : ! Output g^KS
2133 : ncerr = nf90_put_var(spin_ncid, spin_vid("gvals_ks"), my_gbuf_ks, &
2134 : start=[1, 1, 1, 1, gqk%my_kstart, iq_glob], &
2135 1079 : count=[2, gqk%nb_kq, gqk%nb_k, gqk%natom3, gqk%my_nk, iqbuf_cnt])
2136 83 : NCF_CHECK(ncerr)
2137 :
2138 : ncerr = nf90_put_var(root_ncid, root_vid("gstore_glob_state_kqs"), state_kq(:, 1:iqbuf_cnt), &
2139 : start=[gqk%my_kstart, iq_glob, spin], &
2140 581 : count=[gqk%my_nk, iqbuf_cnt, 1])
2141 83 : NCF_CHECK(ncerr)
2142 :
2143 : ! Only one proc sets the entry in done_qbz_spin to 1 for all the q-points in the buffer.
2144 : !if (all(gqk%coords_qkpb_sumbp(2:3) == [0, 0])) then
2145 166 : do ii=1,iqbuf_cnt
2146 83 : iq_bz = iq_buf(2, ii)
2147 332 : NCF_CHECK(nf90_put_var(root_ncid, root_vid("gstore_done_qbz_spin"), 1, start=[iq_bz, spin]))
2148 : end do
2149 : !end if
2150 :
2151 : ! Zero the counter before returning
2152 : !10 iqbuf_cnt = 0
2153 267 : state_kq = GSTORE_KQ_MISSING
2154 :
2155 : ! These syncs can deadlock on lemaitre4 for some MPI decompositions.
2156 83 : if (dtset%useria /= 888) then
2157 83 : NCF_CHECK(nf90_sync(spin_ncid))
2158 83 : NCF_CHECK(nf90_sync(root_ncid))
2159 : end if
2160 :
2161 83 : end subroutine dump_my_gbuf
2162 :
2163 169 : integer function root_vid(var_name)
2164 : character(len=*),intent(in) :: var_name
2165 169 : root_vid = nctk_idname(root_ncid, var_name)
2166 : end function root_vid
2167 :
2168 166 : integer function spin_vid(var_name)
2169 : character(len=*),intent(in) :: var_name
2170 166 : spin_vid = nctk_idname(spin_ncid, var_name)
2171 : end function spin_vid
2172 :
2173 : end subroutine gwpt_run
2174 : !!***
2175 :
2176 : end module m_gwpt
2177 : !!***
|
