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