Line data Source code
1 : !!****m* ABINIT/m_sigma_driver
2 : !! NAME
3 : !! m_sigma_driver
4 : !!
5 : !! FUNCTION
6 : !! Calculate the matrix elements of the self-energy operator.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (MG, GMR, VO, LR, RWG, MT)
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_sigma_driver
23 :
24 : use defs_basis
25 : use m_gwdefs
26 : use defs_wvltypes
27 : use m_xmpi
28 : use m_xomp
29 : use m_errors
30 : use m_abicore
31 : use m_abi_mixing
32 : use m_kxc
33 : use m_distribfft
34 : use netcdf
35 : use m_nctk
36 : use libxc_functionals
37 : use m_dtfil
38 : use m_cgtools
39 :
40 : use defs_datatypes, only : pseudopotential_type
41 : use defs_abitypes, only : MPI_type
42 : use m_time, only : timab
43 : use m_numeric_tools, only : imax_loc
44 : use m_fstrings, only : strcat, sjoin, itoa, basename, ktoa, ltoa
45 : use m_hide_blas, only : xdotc
46 : use m_io_tools, only : open_file, file_exists, iomode_from_fname
47 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
48 : use m_pstat, only : pstat_proc
49 : use m_dtset, only : dataset_type
50 : use m_hdr, only : hdr_type
51 : use m_crystal, only : crystal_t
52 : use m_geometry, only : normv, mkrdim, metric
53 : use m_fftcore, only : print_ngfft
54 : use m_fft_mesh, only : get_gfft, setmesh
55 : use m_fft, only : fourdp
56 : use m_ioarr, only : fftdatar_write, read_rhor
57 : use m_ebands, only : ebands_t, gaps_t
58 : use m_energies, only : energies_type
59 : use m_bz_mesh, only : kmesh_t, littlegroup_t, littlegroup_free, isamek, get_ng0sh
60 : use m_gsphere, only : gsphere_t, merge_and_sort_kg, setshells
61 : use m_kg, only : getph, getcut
62 : use m_xcdata, only : get_xclevel
63 : use m_wfd, only : wfdgw_t, wfdgw_copy, test_charge, wave_t
64 : use m_vcoul, only : vcoul_t
65 : use m_qparticles, only : wrqps, rdqps, rdgw, show_QP, updt_m_ks_to_qp
66 : use m_screening, only : epsm1_t
67 : use m_ppmodel, only : ppmodel_t
68 : use m_sigma, only : sigma_t, write_sigma_header
69 : use m_dyson_solver, only : solve_dyson
70 : use m_esymm, only : esymm_t, esymm_free
71 : use m_melemts, only : melflags_t, melements_t
72 : use m_pawang, only : pawang_type
73 : use m_pawrad, only : pawrad_type
74 : use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize
75 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
76 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print
77 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
78 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, &
79 : pawrhoij_inquire_dim, pawrhoij_symrhoij, pawrhoij_unpack
80 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, paw_overlap
81 : use m_pawdij, only : pawdij, symdij_all
82 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
83 : use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
84 : use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free
85 : use m_paw_slater, only : paw_mkdijexc_core, paw_dijhf
86 : use m_paw_dmft, only : paw_dmft_type
87 : use m_paw_sphharm, only : setsym_ylm
88 : use m_paw_mkrho, only : denfgr
89 : use m_paw_nhat, only : nhatgrid, pawmknhat
90 : use m_paw_tools, only : chkpawovlp, pawprt
91 : use m_paw_denpot, only : pawdenpot
92 : use m_paw_init, only : pawinit, paw_gencond
93 : use m_classify_bands,only : classify_bands
94 : use m_wfk, only : wfk_read_eigenvalues
95 : use m_io_kss, only : make_gvec_kss
96 : use m_vhxc_me, only : calc_vhxc_me
97 : use m_cohsex, only : cohsex_me
98 : use m_sigx, only : calc_sigx_me
99 : use m_sigc, only : calc_sigc_me
100 : use m_setvtr, only : setvtr
101 : use m_mkrho, only : prtrhomxmn
102 : use m_pspini, only : pspini
103 : use m_calc_ucrpa, only : calc_ucrpa
104 : use m_prep_calc_ucrpa,only : prep_calc_ucrpa
105 : use m_paw_correlations,only : pawpuxinit
106 : use m_spacepar, only : hartre
107 : use m_gwrdm, only : calc_rdmx,calc_rdmc,natoccs,update_hdr_bst,print_tot_occ,get_chkprdm,&
108 : print_chkprdm,change_matrix,print_total_energy,print_band_energies,quadrature_sigma_cw
109 : use m_plowannier, only : operwan_realspace_type,plowannier_type,init_plowannier,get_plowannier,&
110 : fullbz_plowannier,init_operwan_realspace,reduce_operwan_realspace,&
111 : destroy_operwan_realspace,destroy_plowannier,zero_operwan_realspace
112 : use minimax_grids, only : gx_minimax_grid
113 :
114 : implicit none
115 :
116 : private
117 : !!***
118 :
119 : public :: sigma
120 : !!***
121 :
122 : contains
123 : !!***
124 :
125 : !!****f* m_sigma_driver/sigma
126 : !! NAME
127 : !! sigma
128 : !!
129 : !! FUNCTION
130 : !! Calculate the matrix elements of the self-energy operator.
131 : !!
132 : !! INPUTS
133 : !! acell(3)=length scales of primitive translations (bohr)
134 : !! codvsn=code version
135 : !! Dtfil<type(datafiles_type)>=variables related to files
136 : !! Dtset<type(dataset_type)>=all input variables for this dataset
137 : !! Pawang<type(pawang_type)>=paw angular mesh and related data
138 : !! Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data
139 : !! Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
140 : !! Psps<type(pseudopotential_type)>=variables related to pseudopotentials
141 : !! Before entering the first time in sigma, a significant part of Psps has been initialized :
142 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
143 : !! and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in
144 : !! the call to pspini. The next time the code enters screening, Psps might be identical to the
145 : !! one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90.
146 : !! rprim(3,3)=dimensionless real space primitive translations
147 : !!
148 : !! OUTPUT
149 : !! Output is written on the main abinit output file. Some results are stored in external files
150 : !!
151 : !! NOTES
152 : !!
153 : !! ON THE USE OF FFT GRIDS:
154 : !! =================
155 : !! In case of PAW:
156 : !! ---------------
157 : !! Two FFT grids are used:
158 : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
159 : !! It is defined by nfft, ngfft, mgfft, ...
160 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
161 : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
162 : !! It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
163 : !! In case of norm-conserving:
164 : !! ---------------------------
165 : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
166 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
167 : !!
168 : !! SOURCE
169 :
170 201 : subroutine sigma(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim)
171 :
172 : !Arguments ------------------------------------
173 : !scalars
174 : character(len=8),intent(in) :: codvsn
175 : type(Datafiles_type),intent(in) :: Dtfil
176 : type(Dataset_type),intent(inout) :: Dtset
177 : type(Pawang_type),intent(inout) :: Pawang
178 : type(Pseudopotential_type),intent(inout) :: Psps
179 : !arrays
180 : real(dp),intent(in) :: acell(3),rprim(3,3)
181 : type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw)
182 : type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw)
183 :
184 : !Local variables-------------------------------
185 : !scalars
186 : integer,parameter :: tim_fourdp5 = 5, master = 0, cplex1 = 1, ipert0 = 0, idir0 = 0, optrhoij1 = 1, ndat1 = 1
187 : integer :: approx_type,b1gw,b2gw,cplex,cplex_dij,cplex_rhoij !,band
188 : integer :: dim_kxcg,gwcalctyp,gnt_option,has_dijU,has_dijso,iab,bmin,bmax,irr_idx1,irr_idx2
189 : integer :: iat,ib,ib1,ib2,ic,id_required,ider,ii,ik,ierr,ount
190 : integer :: ik_bz,ikcalc,ik_ibz,ikxc,npw_k,omp_ncpus,pwx,ibz
191 : integer :: isp,is_idx,istep,itypat,itypatcor,izero,jj,first_band,last_band
192 : integer :: ks_iv,lcor,lmn2_size_max,mband,my_nband
193 : integer :: mgfftf,mod10,moved_atm_inside,moved_rhor,n3xccc !,mgfft
194 : integer :: nbsc,ndij,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot
195 : integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene
196 : integer :: optcut,optgr0,optgr1,optgr2,option,option_test,option_dij,optrad,psp_gencond
197 : integer :: my_rank,rhoxsp_method,comm,use_aerhor,use_umklp,usexcnhat
198 : integer :: ioe0j,spin,io,jb,nomega_sigc, gw1rdm,x1rdm
199 : integer :: temp_unt,ncid
200 : integer :: work_size,nstates_per_proc,my_nbks !, jb_qp,ib_ks,ks_irr
201 : real(dp) :: compch_fft,compch_sph,r_s,rhoav,alpha
202 : real(dp) :: drude_plsmf,my_plsmf,ecore,ecut_eff,ecutdg_eff,ehartree
203 : real(dp) :: etot_sd,etot_mbb,evextnl_energy,ex_energy,gsqcutc_eff,gsqcutf_eff,gsqcut_shp,norm,old_fermie
204 : real(dp) :: eh_energy,ekin_energy,evext_energy,den_int,coef_hyb,exc_mbb_energy,tol_empty
205 : real(dp) :: ucvol,vxcavg,vxcavg_qp,el_temp
206 : real(dp) :: gwc_gsq,gwx_gsq,gw_gsq, gsqcut,boxcut,ecutf
207 : real(dp) :: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem
208 : complex(dp) :: max_degw,cdummy
209 : logical :: rdm_update,readchkprdm,prtchkprdm
210 : logical :: use_paw_aeur,dbg_mode,pole_screening,call_pawinit,is_dfpt=.false.
211 : character(len=500) :: msg
212 : character(len=fnlen) :: wfk_fname,pawden_fname,gw1rdm_fname
213 5226 : type(kmesh_t) :: Kmesh,Qmesh
214 402 : type(ebands_t) :: ks_ebands, qp_ebands
215 12864 : type(vcoul_t) :: Vcp, Vcp_ks, Vcp_full
216 10452 : type(crystal_t) :: Cryst
217 : type(Energies_type) :: KS_energies,QP_energies
218 804 : type(epsm1_t) :: epsm1
219 201 : type(gsphere_t) :: Gsph_Max,Gsph_x,Gsph_c
220 201 : type(hdr_type) :: Hdr_wfk,Hdr_sigma,Hdr_rhor
221 : type(melflags_t) :: KS_mflags,QP_mflags
222 603 : type(melements_t) :: KS_me, QP_me, GW1RDM_me
223 201 : type(MPI_type) :: MPI_enreg_seq
224 201 : type(paw_dmft_type) :: Paw_dmft
225 : type(pawfgr_type) :: Pawfgr
226 201 : type(ppmodel_t) :: PPm
227 201 : type(sigparams_t) :: Sigp
228 201 : type(sigma_t) :: Sr
229 201 : type(wfdgw_t),target :: Wfd, Wfdf, Wfd_nato_master
230 : type(wfdgw_t),pointer :: Wfd_nato_all
231 : type(wave_t),pointer :: wave
232 201 : type(wvl_data) :: Wvl
233 : !arrays
234 : integer :: gwc_ngfft(18),ngfftc(18),ngfftf(18),gwx_ngfft(18), units(2)
235 201 : integer,allocatable :: sigmak_todo(:)
236 201 : integer,allocatable :: nq_spl(:),nlmn_atm(:),my_spins(:)
237 402 : integer,allocatable :: tmp_gfft(:,:),ks_vbik(:,:),nband(:,:),l_size_atm(:),qp_vbik(:,:)
238 402 : integer,allocatable :: tmp_kstab(:,:,:),ks_irreptab(:,:,:),qp_irreptab(:,:,:),my_band_list(:)
239 : real(dp),parameter :: k0(3) = zero
240 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2)
241 603 : real(dp),allocatable :: weights(:),nat_occs(:,:),gw_rhor(:,:),gw_rhog(:,:),gw_vhartr(:)
242 402 : real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:)
243 201 : real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:)
244 201 : real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:), ks_taur(:,:)
245 402 : real(dp),allocatable :: kxc(:,:),qp_kxc(:,:),ph1d(:,:),ph1df(:,:)
246 402 : real(dp),allocatable :: prev_rhor(:,:),prev_taur(:,:),qp_nhat(:,:)
247 201 : real(dp),allocatable :: qp_nhatgr(:,:,:),qp_rhog(:,:),qp_rhor_paw(:,:)
248 201 : real(dp),allocatable :: qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:)
249 201 : real(dp),allocatable :: qp_rhor(:,:),qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:)
250 402 : real(dp),allocatable :: qp_taur(:,:),igwene(:,:,:)
251 402 : real(dp),allocatable :: vpsp(:),xccc3d(:),dijexc_core(:,:,:),dij_hf(:,:,:)
252 201 : real(dp),allocatable :: nl_bks(:,:,:)
253 201 : real(dp),allocatable :: ks_aepaw_rhor(:,:) !,ks_n_one_rhor(:,:),ks_nt_one_rhor(:,:), osoc_bks(:, :, :)
254 : complex(dp) :: ovlp(2)
255 201 : complex(dp),allocatable :: ctmp(:,:),hbare(:,:,:,:)
256 201 : complex(dp),target,allocatable :: sigcme(:,:,:,:,:)
257 402 : complex(dp),allocatable :: hdft(:,:,:,:),htmp(:,:,:,:),uks2qp(:,:)
258 603 : complex(dp),allocatable :: xrdm_k_full(:,:,:), rdm_k(:,:), pot_k(:,:), nateigv(:,:,:,:), old_ks_purex(:,:), new_hartr(:,:)
259 402 : complex(gwp),allocatable :: kxcg(:,:),fxc_ADA(:,:,:)
260 201 : complex(gwp),contiguous, pointer :: ug1(:)
261 402 : complex(dp),allocatable :: sigcme_k(:,:,:,:), rhot1_q_m(:,:,:,:,:,:,:), M1_q_m(:,:,:,:,:,:,:)
262 201 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:),bmask(:), bdm_mask(:,:,:),bdm2_mask(:,:,:)
263 201 : type(esymm_t),target,allocatable :: KS_sym(:,:)
264 201 : type(esymm_t),pointer :: QP_sym(:,:)
265 201 : type(pawcprj_type),allocatable :: Cp1(:,:) !,Cp2(:,:)
266 201 : type(littlegroup_t),allocatable :: Ltg_k(:)
267 402 : type(Paw_an_type),allocatable :: KS_paw_an(:),QP_paw_an(:)
268 402 : type(Paw_ij_type),allocatable :: KS_paw_ij(:),QP_paw_ij(:)
269 201 : type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
270 402 : type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:),QP_pawrhoij(:),prev_Pawrhoij(:),tmp_pawrhoij(:)
271 201 : type(pawpwff_t),allocatable :: Paw_pwff(:)
272 201 : type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
273 201 : type(plowannier_type) :: wanbz,wanibz,wanibz_in
274 201 : type(operwan_realspace_type), allocatable :: rhot1(:,:)
275 : !************************************************************************
276 :
277 : DBG_ENTER('COLL')
278 :
279 201 : call timab(401,1,tsec) ! sigma(Total)
280 201 : call timab(402,1,tsec) ! sigma(Init1)
281 603 : units = [std_out, ab_out]
282 :
283 : write(msg,'(7a)')&
284 201 : ' SIGMA: Calculation of the GW corrections ',ch10,ch10,&
285 201 : ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,&
286 402 : ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.'
287 201 : call wrtout(units, msg)
288 :
289 201 : if(dtset%ucrpa>0) then
290 0 : write(msg,'(6a)')ch10,' cRPA Calculation: Calculation of the screened Coulomb interaction (ucrpa/=0) ',ch10
291 0 : call wrtout(units, msg)
292 : end if
293 :
294 : #if defined HAVE_GW_DPC
295 : if (gwp /= 8) then
296 : write(msg,'(6a)')ch10,&
297 : 'Number of bytes for double precision complex /=8 ',ch10,&
298 : 'Cannot continue due to kind mismatch in BLAS library ',ch10,&
299 : 'Some BLAS interfaces are not generated by abilint '
300 : ABI_ERROR(msg)
301 : end if
302 201 : write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwp,ch10
303 : #else
304 : write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwp,ch10
305 : #endif
306 201 : call wrtout(units, msg)
307 :
308 201 : tol_empty=0.01 ! Initialize the tolerance used to decide if a band is empty (passed to m_sigx.F90)
309 201 : gwcalctyp=Dtset%gwcalctyp
310 201 : gw1rdm=Dtset%gw1rdm ! Input variable to decide if updates to the 1-RDM must be performed
311 201 : x1rdm=Dtset%x1rdm ! Input variable to use pure exchange correction on the 1-RDM ( Sigma_x - Vxc )
312 201 : rdm_update=(gwcalctyp==21 .and. gw1rdm>0) ! Input variable to decide whether to update GW density matrix
313 201 : readchkprdm=(Dtset%irdchkprdm==1) ! Input variable to decide if checkpoint files must be read
314 201 : prtchkprdm=(Dtset%prtchkprdm==1) ! Input variable to decide if checkpoint files must be written
315 :
316 201 : mod10 =MOD(Dtset%gwcalctyp,10)
317 :
318 : ! Perform some additional checks for hybrid functional calculations
319 201 : if (mod10 == 5) then
320 : !if (Dtset%ixc_sigma<0 .and. .not.libxc_functionals_check()) then
321 : ! ABI_ERROR('Hybrid functional calculations require the compilation with LIBXC library')
322 : !end if
323 : !XG 20171116 : I do not agree with this condition, as one might like to do a one-shot hybrid functional calculation
324 : !on top of a LDA/GGA calculation ... give the power (and risks) to the user !
325 : !if(gwcalctyp<10) then
326 : ! ABI_ERROR('gwcalctyp requires the update of energies and/or wavefunctions when performing hybrid XC calculations')
327 : !end if
328 35 : if (Dtset%usepaw == 1) then
329 0 : ABI_ERROR('PAW version of hybrid functional calculations is not implemented')
330 : end if
331 : end if
332 :
333 : !=== Initialize MPI variables, and parallelization level ===
334 : ! gwpara: 1--> parallelism over k-points, 2--> parallelism over bands.
335 : ! In case of gwpara==1 memory is not parallelized.
336 : ! If gwpara==2, bands are divided among processors but each proc has all the states where GW corrections are required.
337 201 : comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
338 :
339 201 : if (my_rank == master) then
340 173 : wfk_fname = dtfil%fnamewffk
341 173 : if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then
342 0 : ABI_ERROR(msg)
343 : end if
344 : end if
345 201 : call xmpi_bcast(wfk_fname, master, comm, ierr)
346 :
347 : ! === Some variables need to be initialized/nullify at start ===
348 201 : usexcnhat = 0
349 201 : call KS_energies%init()
350 201 : call mkrdim(acell, rprim, rprimd)
351 201 : call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
352 : !
353 : ! === Define FFT grid(s) sizes ===
354 : ! Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the
355 : ! (usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90.
356 : ! See also NOTES in the comments at the beginning of this file.
357 : ! NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn.
358 :
359 : call pawfgr_init(Pawfgr, Dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, &
360 201 : gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=gmet, k0=k0)
361 :
362 : ! Fake MPI_type for the sequential part.
363 201 : call initmpi_seq(MPI_enreg_seq)
364 201 : call MPI_enreg_seq%distribfft%init_seq('c',ngfftc(2),ngfftc(3),'all')
365 201 : call MPI_enreg_seq%distribfft%init_seq('f',ngfftf(2),ngfftf(3),'all')
366 :
367 402 : call print_ngfft([std_out], ngfftf, header='Dense FFT mesh used for densities and potentials')
368 804 : nfftf_tot=PRODUCT(ngfftf(1:3))
369 :
370 : ! ===========================================
371 : ! === Open and read pseudopotential files ===
372 : ! ===========================================
373 201 : call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm)
374 :
375 201 : call timab(402,2,tsec) ! Init1
376 : !
377 : ! ==================================================
378 : ! ==== Initialize Sigp, epsm1 and basic objects ====
379 : ! ==================================================
380 : ! Sigp is completely initialized here.
381 : ! epsm1 is only initialized with dimensions, (SCR|SUSC) file is read in epsm1%mkdump
382 201 : call timab(403,1,tsec) ! setup_sigma
383 :
384 : call setup_sigma(codvsn,wfk_fname,acell,rprim,Dtset,Dtfil,Psps,Pawtab,&
385 201 : gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_sigma,Cryst,Kmesh,Qmesh,ks_ebands,Gsph_Max,Gsph_x,Gsph_c,Vcp,epsm1,Sigp,comm)
386 :
387 201 : call pstat_proc%print(_PSTAT_ARGS_)
388 :
389 201 : call timab(403,2,tsec) ! setup_sigma
390 201 : call timab(402,1,tsec) ! Init1
391 :
392 201 : if (nprocs > Sigp%nbnds) then
393 0 : write(msg,"(2(a,i0))")"The number of MPI procs: ", nprocs, " is greater than nband: ", sigp%nbnds
394 0 : ABI_ERROR(msg)
395 : end if
396 :
397 201 : pole_screening = .FALSE.
398 201 : if (epsm1%fform==2002) then
399 0 : pole_screening = .TRUE.
400 0 : ABI_WARNING(' EXPERIMENTAL - Using a pole-fit screening!')
401 : end if
402 :
403 402 : call print_ngfft([std_out], gwc_ngfft, header='FFT mesh for oscillator strengths used for Sigma_c')
404 402 : call print_ngfft([std_out], gwx_ngfft, header='FFT mesh for oscillator strengths used for Sigma_x')
405 :
406 201 : b1gw = Sigp%minbdgw; b2gw = Sigp%maxbdgw
407 :
408 804 : gwc_nfftot=PRODUCT(gwc_ngfft(1:3))
409 : gwc_nfft =gwc_nfftot !no FFT //
410 :
411 : gwx_nfftot=PRODUCT(gwx_ngfft(1:3))
412 : gwx_nfft =gwx_nfftot !no FFT //
413 :
414 : ! TRYING TO RECREATE AN "ABINIT ENVIRONMENT"
415 201 : KS_energies%e_corepsp=ecore/Cryst%ucvol
416 :
417 : ! === Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ====
418 : ! * ks_vbk gives the (valence|last Fermi band) index for each k and spin.
419 : ! * spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals
420 804 : ABI_MALLOC(ks_vbik, (ks_ebands%nkpt, ks_ebands%nsppol))
421 603 : ABI_MALLOC(qp_vbik, (ks_ebands%nkpt, ks_ebands%nsppol))
422 :
423 : !call ks_ebands%update_occ(Dtset%spinmagntarget,prtvol=0)
424 201 : ks_vbik(:,:) = ks_ebands%get_valence_idx()
425 :
426 : ! ============================
427 : ! ==== PAW initialization ====
428 : ! ============================
429 201 : if (dtset%usepaw == 1) then
430 5 : call chkpawovlp(cryst%natom, cryst%ntypat, dtset%pawovlp, pawtab, cryst%rmet, cryst%typat, cryst%xred)
431 :
432 5 : cplex_dij = dtset%nspinor; cplex = 1; ndij = 1
433 :
434 46 : ABI_MALLOC(ks_pawrhoij, (cryst%natom))
435 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij, nspden_rhoij=nspden_rhoij, &
436 5 : nspden=dtset%nspden, spnorb=dtset%pawspnorb, cpxocc=dtset%pawcpxocc)
437 5 : call pawrhoij_alloc(ks_pawrhoij, cplex_rhoij, nspden_rhoij, dtset%nspinor, dtset%nsppol, cryst%typat, pawtab=pawtab)
438 :
439 : ! Test if we have to call pawinit
440 5 : gnt_option = 1; if (dtset%pawxcdev == 2 .or. (dtset%pawxcdev == 1 .and. dtset%positron /= 0)) gnt_option = 2
441 5 : call paw_gencond(dtset, gnt_option, "test", call_pawinit)
442 :
443 5 : if (psp_gencond == 1 .or. call_pawinit) then
444 0 : call timab(553, 1, tsec)
445 0 : gsqcut_shp = two * abs(dtset%diecut) * dtset%dilatmx**2 / pi**2
446 : call pawinit(dtset%effmass_free, gnt_option, gsqcut_shp, zero, dtset%pawlcutd, dtset%pawlmix, &
447 : psps%mpsang, dtset%pawnphi, cryst%nsym, dtset%pawntheta, pawang, pawrad, &
448 0 : dtset%pawspnorb, pawtab, dtset%pawxcdev, dtset%ixc, dtset%usepotzero)
449 0 : call timab(553,2,tsec)
450 :
451 : ! Update internal values
452 0 : call paw_gencond(dtset, gnt_option, "save", call_pawinit)
453 : else
454 5 : if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij = 2
455 5 : if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla = 2
456 : end if
457 :
458 13 : psps%n1xccc = maxval(pawtab(1:cryst%ntypat)%usetcore)
459 :
460 : ! Initialize optional flags in Pawtab to zero
461 : ! Cannot be done in Pawinit since the routine is called only if some parts are changed
462 13 : pawtab(:)%has_nabla = 0
463 13 : pawtab(:)%lamb_shielding = zero
464 :
465 5 : call setsym_ylm(gprimd, pawang%l_max-1, cryst%nsym, dtset%pawprtvol, cryst%rprimd, cryst%symrec, pawang%zarot)
466 :
467 : ! Initialize and compute data for DFT+U
468 5 : Paw_dmft%use_dmft=Dtset%usedmft
469 : call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,&
470 : is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,&
471 5 : Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
472 :
473 5 : if (my_rank == master) call pawtab_print(Pawtab)
474 :
475 : ! Get Pawrhoij from the header of the WFK file.
476 5 : call pawrhoij_copy(Hdr_wfk%pawrhoij, KS_Pawrhoij)
477 :
478 : ! Evaluate form factor of radial part of phi.phj - tphi.tphj.
479 : ! The q-grid must contain the FFT mesh used for sigma_c and the G-sphere for the exchange part.
480 : ! We use the FFT mesh for sigma_c since COHSEX and the extrapolar method require oscillator
481 : ! strengths on the FFT mesh.
482 15 : ABI_MALLOC(tmp_gfft,(3, gwc_nfftot))
483 5 : call get_gfft(gwc_ngfft, k0, gmet, gwc_gsq, tmp_gfft)
484 5 : ABI_FREE(tmp_gfft)
485 :
486 : ! Set up q-grid, make qmax 20% larger than largest expected.
487 15 : ABI_MALLOC(nq_spl, (Psps%ntypat))
488 15 : ABI_MALLOC(qmax, (Psps%ntypat))
489 5 : gwx_gsq = Dtset%ecutsigx / (two*pi**2)
490 5 : gw_gsq = max(gwx_gsq, gwc_gsq)
491 13 : qmax = sqrt(gw_gsq)*1.2d0
492 13 : nq_spl = Psps%mqgrid_ff
493 : ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax
494 :
495 5 : rhoxsp_method = 1 ! Arnaud-Alouani (default in sigma)
496 : !rhoxsp_method = 2 ! Shiskin-Kresse
497 5 : if (dtset%pawoptosc /= 0) rhoxsp_method = dtset%pawoptosc
498 :
499 83 : ABI_MALLOC(paw_pwff, (psps%ntypat))
500 5 : call pawpwff_init(paw_pwff, rhoxsp_method, nq_spl, qmax, gmet, pawrad, pawtab, psps)
501 :
502 5 : ABI_FREE(nq_spl)
503 5 : ABI_FREE(qmax)
504 :
505 : ! Variables/arrays related to the fine FFT grid
506 169497 : ABI_CALLOC(ks_nhat, (nfftf, Dtset%nspden))
507 :
508 46 : ABI_MALLOC(pawfgrtab, (cryst%natom))
509 5 : call pawtab_get_lsize(pawtab, l_size_atm, cryst%natom, cryst%typat)
510 :
511 5 : cplex = 1
512 5 : call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat)
513 5 : ABI_FREE(l_size_atm)
514 5 : compch_fft=greatest_real
515 13 : usexcnhat = MAXVAL(Pawtab(:)%usexcnhat)
516 : ! * 0 if Vloc in atomic data is Vbare (Blochl's formulation)
517 : ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse's formulation)
518 5 : call wrtout(std_out, sjoin(' using usexcnhat: ', itoa(usexcnhat)))
519 : !
520 : ! Identify parts of the rectangular grid where the density has to be calculated ===
521 5 : optcut = 0; optgr0 = Dtset%pawstgylm; optgr1 = 0; optgr2 = 0; optrad = 1 - Dtset%pawstgylm
522 5 : if (Dtset%pawcross==1) optrad=1
523 5 : if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
524 :
525 : call nhatgrid(Cryst%atindx1, gmet, Cryst%natom, Cryst%natom, Cryst%nattyp, ngfftf, Cryst%ntypat,&
526 5 : optcut, optgr0, optgr1, optgr2, optrad, Pawfgrtab, Pawtab, Cryst%rprimd, Cryst%typat, Cryst%ucvol, Cryst%xred)
527 :
528 20 : call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol)
529 :
530 : else
531 196 : ABI_MALLOC(Paw_pwff, (0))
532 196 : ABI_MALLOC(Pawfgrtab, (0))
533 : end if ! End of PAW Initialization
534 :
535 : ! Consistency check and additional stuff done only for GW with PAW.
536 201 : ABI_MALLOC(Paw_onsite, (0))
537 :
538 201 : if (Dtset%usepaw == 1) then
539 5 : if (Dtset%ecutwfn < Dtset%ecut) then
540 : write(msg,"(3a)")&
541 0 : " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,&
542 0 : " an excessive truncation of the planewave basis set can lead to unphysical results."
543 0 : ABI_WARNING(msg)
544 : end if
545 :
546 5 : ABI_CHECK(Dtset%useexexch == 0, "LEXX not yet implemented in GW")
547 5 : ABI_CHECK(Paw_dmft%use_dmft == 0, "DMFT + GW not available")
548 :
549 : ! Optionally read core orbitals from file and calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for the HF decoupling.
550 5 : if (Sigp%use_sigxcore == 1) then
551 2 : lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
552 5 : ABI_MALLOC(dijexc_core,(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat))
553 :
554 1 : call paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,Dtset%prtvol,Psps%filpsp)
555 : end if ! HF decoupling
556 :
557 5 : if (Dtset%pawcross==1) then
558 0 : ABI_SFREE(Paw_onsite)
559 0 : ABI_MALLOC(Paw_onsite,(Cryst%natom))
560 : call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat, &
561 0 : Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab)
562 : end if
563 : end if
564 :
565 : ! Allocate these arrays anyway, since they are passed to subroutines.
566 201 : if (.not.allocated(ks_nhat)) then
567 392 : ABI_MALLOC(ks_nhat, (nfftf, 0))
568 : end if
569 201 : if (.not.allocated(dijexc_core)) then
570 200 : ABI_MALLOC(dijexc_core, (1, 1, 0))
571 : end if
572 :
573 : ! ==================================================
574 : ! ==== Read KS band structure from the KSS file ====
575 : ! ==================================================
576 : !
577 : ! Initialize Wfd, allocate wavefunctions and precalculate tables to do the FFT using the coarse gwc_ngfft.
578 201 : mband=Sigp%nbnds
579 1005 : ABI_MALLOC(bks_mask, (mband, Kmesh%nibz, Sigp%nsppol))
580 804 : ABI_MALLOC(keep_ur , (mband, Kmesh%nibz, Sigp%nsppol))
581 61415 : keep_ur=.FALSE.; bks_mask=.FALSE.
582 :
583 201 : if (rdm_update) then
584 20 : ABI_MALLOC(bdm_mask, (mband, Kmesh%nibz, Sigp%nsppol))
585 280 : bdm_mask=.FALSE.
586 5 : if (my_rank==master) then
587 280 : bdm_mask=.TRUE.
588 : end if
589 20 : ABI_MALLOC(bdm2_mask, (mband,Kmesh%nibz,Sigp%nsppol))
590 280 : bdm2_mask=.FALSE.
591 : end if
592 804 : ABI_MALLOC(nband,(Kmesh%nibz, Sigp%nsppol))
593 1638 : nband=mband
594 :
595 : ! autoparal section
596 201 : if (dtset%max_ncpus/=0) then
597 0 : ount =ab_out
598 : ! Temporary table needed to estimate memory
599 0 : ABI_MALLOC(nlmn_atm,(Cryst%natom))
600 0 : if (Dtset%usepaw==1) then
601 0 : do iat=1,Cryst%natom
602 0 : nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size
603 : end do
604 : end if
605 :
606 0 : write(ount,'(a)')"--- !Autoparal"
607 0 : write(ount,"(a)")'#Autoparal section for Sigma runs.'
608 0 : write(ount,"(a)") "info:"
609 0 : write(ount,"(a,i0)")" autoparal: ",dtset%autoparal
610 0 : write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus
611 0 : write(ount,"(a,i0)")" gwpara: ",dtset%gwpara
612 0 : write(ount,"(a,i0)")" nkpt: ",dtset%nkpt
613 0 : write(ount,"(a,i0)")" nsppol: ",dtset%nsppol
614 0 : write(ount,"(a,i0)")" nspinor: ",dtset%nspinor
615 0 : write(ount,"(a,i0)")" nbnds: ",Sigp%nbnds
616 :
617 0 : work_size = mband * Kmesh%nibz * Sigp%nsppol
618 :
619 : ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI.
620 0 : nonscal_mem = (two*gwp*epsm1%npwe**2*epsm1%nomega*(epsm1%mqmem+1)*b2Mb) * 1.1_dp
621 :
622 : ! List of configurations.
623 : ! Assuming an OpenMP implementation with perfect speedup!
624 0 : write(ount,"(a)")"configurations:"
625 0 : do ii=1,dtset%max_ncpus
626 0 : nstates_per_proc = 0
627 0 : eff = HUGE(one)
628 0 : max_wfsmem_mb = zero
629 :
630 0 : do my_rank=0,ii-1
631 0 : call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,ii,my_spins,bks_mask,keep_ur,ierr)
632 0 : ABI_FREE(my_spins)
633 0 : if (ierr /= 0) exit
634 0 : my_nbks = COUNT(bks_mask)
635 0 : nstates_per_proc = MAX(nstates_per_proc, my_nbks)
636 0 : eff = MIN(eff, (one * work_size) / (ii * nstates_per_proc))
637 :
638 : ! Memory needed for Fourier components ug.
639 0 : ug_mem = two*gwp*Dtset%nspinor*Sigp%npwwfn*my_nbks*b2Mb
640 : ! Memory needed for real space ur (use gwc_nfft, instead of gwx_nfft)
641 0 : ur_mem = two*gwp*Dtset%nspinor*gwc_nfft*COUNT(keep_ur)*b2Mb
642 : ! Memory needed for PAW projections Cprj
643 0 : cprj_mem = zero
644 0 : if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb
645 0 : max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem)
646 : end do
647 0 : if (ierr /= 0) cycle
648 :
649 : ! Add the non-scalable part and increase by 10% to account for other datastructures.
650 0 : mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
651 0 : do omp_ncpus=1,xomp_get_max_threads()
652 0 : write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus
653 0 : write(ount,"(a,i0)")" mpi_ncpus: ",ii
654 0 : write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus
655 0 : write(ount,"(a,f12.9)")" efficiency: ",eff
656 0 : write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb
657 : end do
658 : end do
659 0 : write(ount,'(a)')"..."
660 :
661 0 : ABI_FREE(nlmn_atm)
662 0 : ABI_ERROR_NODUMP("aborting now")
663 :
664 : else
665 201 : call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr)
666 201 : ABI_CHECK(ierr==0, "Error in sigma_bksmask")
667 : end if
668 :
669 : ! Each core stores the wavefunctions where GW corrections are required.
670 406 : do isp=1,SIZE(my_spins)
671 205 : spin = my_spins(isp)
672 1053 : do ikcalc=1,Sigp%nkptgw
673 647 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
674 647 : ii=Sigp%minbnd(ikcalc,spin); jj=Sigp%maxbnd(ikcalc,spin)
675 6098 : bks_mask(ii:jj,ik_ibz,spin) = .TRUE.
676 6154 : if (MODULO(Dtset%gwmem,10)==1) keep_ur(ii:jj,ik_ibz,spin)=.TRUE.
677 : end do
678 : end do
679 :
680 201 : ABI_FREE(my_spins)
681 :
682 : call wfd%init(Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,&
683 : Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
684 201 : Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,use_fnl_dir0der0=(Dtset%gw1rdm==2))
685 :
686 : ! MRM: also initialize the Wfd_nato_master for GW 1-RDM if required.
687 : ! Warning, this should be replaced by copy but copy fails due to bands being allocated in different manners.
688 : ! FIXME: Do it in the future!
689 201 : if (rdm_update) then
690 285 : bdm2_mask=bks_mask ! As bks_mask is going to be removed, save it in bdm2_mask to use it in Evext_nl
691 : call Wfd_nato_master%init(Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bdm_mask,&
692 : Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
693 5 : Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,xmpi_comm_self)!comm) ! MPI_COMM_SELF
694 5 : call Wfd_nato_master%read_wfk(wfk_fname,iomode_from_fname(wfk_fname))
695 : end if
696 :
697 201 : call timab(402,2,tsec) ! sigma(Init1)
698 201 : call timab(404,1,tsec) ! rdkss
699 :
700 201 : call wfd%read_wfk(wfk_fname, iomode_from_fname(wfk_fname))
701 : ! This test has been disabled (too expensive!)
702 : if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL")
703 :
704 201 : if (Dtset%pawcross==1 .or. dtset%userie == 456) then
705 : call Wfdf%init(Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,&
706 : Dtset%nspden,Dtset%nspinor,dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,&
707 0 : Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm)
708 0 : if (dtset%userie == 456) then
709 0 : call wrtout(std_out, "Reading states from supercell WFK file")
710 0 : call wfdf%read_wfk("SC_WFK", iomode_from_fname("SC_WFK"))
711 : else
712 0 : call wfdgw_copy(Wfd, Wfdf)
713 : end if
714 0 : call wfdf%change_ngfft(Cryst, Psps, ngfftf)
715 0 : if (dtset%userie == 456) call wfdf%print(units, "UNPERTURBED WFDF for GWTPT")
716 : end if
717 :
718 201 : call pstat_proc%print(_PSTAT_ARGS_)
719 201 : call timab(404,2,tsec) ! rdkss
720 201 : call timab(405,1,tsec) ! Init2
721 :
722 201 : ABI_FREE(bks_mask)
723 201 : ABI_FREE(nband)
724 201 : ABI_FREE(keep_ur)
725 :
726 : ! ==============================================================
727 : ! ==== Find little group of the k-points for GW corrections ====
728 : ! ==============================================================
729 : ! The little group is used only if symsigma == 1
730 : ! If use_umklp == 1 then symmetries requiring an umklapp to preserve k_gw are included as well.
731 1242 : ABI_MALLOC(Ltg_k, (Sigp%nkptgw))
732 201 : use_umklp = 1
733 840 : do ikcalc=1,Sigp%nkptgw
734 840 : if (Sigp%symsigma /= 0) call Ltg_k(ikcalc)%init(Sigp%kptgw(:,ikcalc), Qmesh%nbz, Qmesh%bz, Cryst, use_umklp, npwe=0)
735 : end do
736 :
737 : ! Compute structure factor phases and large sphere cut-off
738 603 : ABI_MALLOC(ph1d,(2, 3 * (2 * Dtset%mgfft + 1) * Cryst%natom))
739 603 : ABI_MALLOC(ph1df,(2, 3 * (2 * mgfftf + 1) * Cryst%natom))
740 :
741 201 : call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred)
742 :
743 201 : if (Psps%usepaw == 1.and. Pawfgr%usefinegrid == 1) then
744 2 : call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred)
745 : else
746 170362 : ph1df(:,:)=ph1d(:,:)
747 : end if
748 :
749 : !===================================================================================
750 : !==== Classify the GW wavefunctions according to the irreducible representation ====
751 : !===================================================================================
752 : !* Warning still under development.
753 : !* Only for SCGW.
754 : !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
755 :
756 2241 : ABI_MALLOC(KS_sym,(Wfd%nkibz,Wfd%nsppol))
757 :
758 201 : if (Sigp%symsigma==1.and.gwcalctyp>=20) then
759 : ! call check_zarot(Gsph_c%ng,Cryst,gwc_ngfft,Gsph_c%gvec,Psps,Pawang,Gsph_c%rottb,Gsph_c%rottbm1)
760 0 : use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric
761 0 : do spin=1,Wfd%nsppol
762 0 : do ikcalc=1,Sigp%nkptgw
763 0 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
764 0 : first_band = Sigp%minbnd(ikcalc,spin)
765 0 : last_band = Sigp%maxbnd(ikcalc,spin)
766 : call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,ks_ebands,Pawtab,Pawrad,Pawang,Psps,&
767 0 : Dtset%tolsym,KS_sym(ik_ibz,spin))
768 : end do
769 : end do
770 : ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
771 0 : call sigma_tables(Sigp,Kmesh, esymm=KS_sym)
772 : end if
773 :
774 201 : call timab(405,2,tsec) ! Init2
775 201 : call timab(406,1,tsec) ! make_vhxc
776 :
777 : !===========================
778 : !=== COMPUTE THE DENSITY ===
779 : !===========================
780 : ! Evaluate the planewave part (complete charge in case of NC pseudos).
781 804 : ABI_MALLOC(ks_rhor, (nfftf, dtset%nspden))
782 804 : ABI_MALLOC(ks_taur, (nfftf, dtset%nspden * dtset%usekden))
783 :
784 201 : call wfd%mkrho(cryst, psps, ks_ebands, ngfftf, nfftf, ks_rhor)
785 :
786 201 : if ((rdm_update .and. Dtset%prtden /= 0) .and. Wfd%my_rank == master) then
787 : ! Print initial (KS) density file as read (useful to compare DEN files, cubes, etc.)
788 5 : gw1rdm_fname = trim(dtfil%fnameabo_ks_den) ! and used on Sigma grids
789 : call fftdatar_write("density",gw1rdm_fname,dtset%iomode,hdr_sigma,&
790 5 : Cryst,ngfftf,cplex1,nfftf,dtset%nspden,ks_rhor,mpi_enreg_seq,ebands=ks_ebands)
791 : end if
792 :
793 201 : if (Dtset%usekden == 1) call wfd%mkrho(cryst, psps, ks_ebands, ngfftf, nfftf, ks_taur, optcalc=1)
794 :
795 : !========================================
796 : !==== Additional computation for PAW ====
797 : !========================================
798 201 : nhatgrdim = 0
799 201 : if (Dtset%usepaw==1) then
800 : ! Get electronic temperature from dtset
801 5 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
802 : ! Calculate the compensation charge nhat.
803 5 : if (Dtset%xclevel==2) nhatgrdim = usexcnhat * Dtset%pawnhatxc
804 5 : cplex = 1; ider = 2 * nhatgrdim; izero = 0
805 5 : if (nhatgrdim > 0) then
806 0 : ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
807 : end if
808 5 : if (nhatgrdim == 0) then
809 5 : ABI_MALLOC(ks_nhatgr,(0,0,0))
810 : end if
811 :
812 : call pawmknhat(compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,&
813 : Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
814 : Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
815 5 : Cryst%ucvol,dtset%usewvl,Cryst%xred)
816 :
817 : ! === Evaluate onsite energies, potentials, densities ===
818 : ! * Initialize variables/arrays related to the PAW spheres.
819 : ! * Initialize also lmselect (index of non-zero LM-moments of densities).
820 46 : ABI_MALLOC(KS_paw_ij, (Cryst%natom))
821 5 : has_dijso = Dtset%pawspnorb; has_dijU = merge(0, 1, Dtset%usepawu == 0)
822 :
823 5 : call paw_ij_nullify(KS_paw_ij)
824 : call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
825 : Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
826 : has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
827 5 : has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
828 :
829 5 : nkxc1 = 0
830 46 : ABI_MALLOC(KS_paw_an, (Cryst%natom))
831 5 : call paw_an_nullify(KS_paw_an)
832 : call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
833 5 : cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1)
834 :
835 : ! Calculate onsite vxc with and without core charge.
836 5 : nzlmopt=-1; option=0; compch_sph=greatest_real
837 : call pawdenpot(compch_sph,el_temp,Cryst%gprimd,ipert0,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
838 : Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_energies%paw,KS_paw_ij,&
839 : Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
840 5 : Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,Cryst%xred,Cryst%ucvol,Psps%znuclpsp,Dtset%spinaxis)
841 :
842 : else
843 196 : ABI_MALLOC(ks_nhatgr, (0, 0, 0))
844 196 : ABI_MALLOC(KS_paw_ij, (0))
845 196 : ABI_MALLOC(KS_paw_an, (0))
846 : end if ! PAW
847 :
848 : call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,&
849 201 : Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
850 :
851 : ! For PAW, add the compensation charge on the FFT mesh, then get rho(G).
852 169683 : if (Dtset%usepaw==1) ks_rhor = ks_rhor + ks_nhat
853 :
854 201 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol)
855 201 : if (Dtset%usekden==1) call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=ucvol)
856 :
857 : ! FFT n(r) --> n(g)
858 603 : ABI_MALLOC(ks_rhog, (2, nfftf))
859 201 : call fourdp(1, ks_rhog, ks_rhor(:,1),-1, MPI_enreg_seq, nfftf, 1, ngfftf, tim_fourdp5)
860 :
861 : ! The following steps have been gathered in the setvtr routine:
862 : ! - get Ewald energy and Ewald forces
863 : ! - compute local ionic pseudopotential vpsp
864 : ! - eventually compute 3D core electron density xccc3d
865 : ! - eventually compute vxc and vhartr
866 : ! - set up ks_vtrial
867 : !
868 : !*******************************************************************
869 : !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
870 : !*******************************************************************
871 :
872 201 : ngrvdw = 0
873 201 : ABI_MALLOC(grvdw, (3, ngrvdw))
874 603 : ABI_MALLOC(grchempottn, (3, Cryst%natom))
875 402 : ABI_MALLOC(grewtn, (3, Cryst%natom))
876 201 : nkxc = 0
877 201 : if (Dtset%nspden == 1) nkxc = 2
878 201 : if (Dtset%nspden >= 2) nkxc = 3 ! check GGA and spinor, quite a messy part!!!
879 : ! In case of MGGA, fxc and kxc are not available and we dont need them in sigma (for now ...)
880 201 : if (Dtset%ixc < 0 .and. libxc_functionals_ismgga()) nkxc = 0
881 201 : if (nkxc /= 0) then
882 804 : ABI_MALLOC(kxc, (nfftf, nkxc))
883 : end if
884 :
885 201 : n3xccc = 0; if (Psps%n1xccc /= 0) n3xccc = nfftf
886 603 : ABI_MALLOC(xccc3d, (n3xccc))
887 603 : ABI_MALLOC(ks_vhartr, (nfftf))
888 804 : ABI_MALLOC(ks_vtrial, (nfftf, Dtset%nspden))
889 402 : ABI_MALLOC(vpsp, (nfftf))
890 603 : ABI_MALLOC(ks_vxc, (nfftf, Dtset%nspden))
891 :
892 201 : optene = 4; moved_atm_inside = 0; moved_rhor = 0; istep = 1
893 :
894 : call setvtr(Cryst%atindx1,Dtset,KS_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
895 : istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
896 : Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
897 : optene,Pawang,Pawrad,KS_Pawrhoij,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,&
898 201 : Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred,taur=ks_taur)
899 :
900 : !============================
901 : !==== Compute KS PAW Dij ====
902 : !============================
903 201 : if (Dtset%usepaw == 1) then
904 5 : call timab(561,1,tsec)
905 :
906 : ! Calculate the unsymmetrized Dij.
907 : call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,&
908 : Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
909 : Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,&
910 : Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
911 : k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),&
912 : ks_vtrial,ks_vxc,Cryst%xred,Dtset%znucl,&
913 5 : nucdipmom=Dtset%nucdipmom,spinaxis=Dtset%spinaxis)
914 :
915 : ! Symmetrize KS Dij
916 : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,&
917 : Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,&
918 5 : Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
919 :
920 : ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
921 5 : call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab)
922 5 : call timab(561,2,tsec)
923 : end if
924 :
925 201 : call timab(406,2,tsec) ! make_vhxc
926 :
927 : !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW ===
928 : ! * This part is parallelized within wfd%comm since each node has all GW wavefunctions.
929 : ! * Note that vH matrix elements are calculated using the true uncutted interaction.
930 201 : call timab(407,1,tsec) ! vHxc_me
931 :
932 201 : call KS_mflags%reset()
933 201 : if (rdm_update) then
934 5 : KS_mflags%has_hbare=1
935 5 : KS_mflags%has_kinetic=1
936 : end if
937 201 : KS_mflags%has_vhartree=1
938 201 : KS_mflags%has_vxc =1
939 201 : KS_mflags%has_vxcval =1
940 201 : if (Dtset%usepawu /= 0 ) KS_mflags%has_vu = 1
941 201 : if (Dtset%useexexch /= 0) KS_mflags%has_lexexch= 1
942 201 : if (Sigp%use_sigxcore == 1) KS_mflags%has_sxcore = 1
943 201 : if (gwcalctyp<10 ) KS_mflags%only_diago = 1 ! off-diagonal elements only for SC on wavefunctions.
944 :
945 : if (.FALSE.) then ! quick and dirty hack to test HF contribution.
946 : ABI_WARNING("testing on-site HF")
947 : lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
948 : ABI_MALLOC(dij_hf,(cplex_dij*lmn2_size_max,ndij,Cryst%natom))
949 : call paw_dijhf(ndij,cplex_dij,1,lmn2_size_max,Cryst%natom,Cryst%ntypat,Pawtab,Pawrad,Pawang,&
950 : KS_Pawrhoij,dij_hf,Dtset%prtvol)
951 :
952 : do iat=1,Cryst%natom
953 : itypat = Cryst%typat(iat)
954 : ii = Pawtab(itypat)%lmn2_size
955 : KS_Paw_ij(iat)%dijxc(:,:) = dij_hf(1:cplex_dij*ii,:,iat)
956 : KS_Paw_ij(iat)%dijxc_hat(:,:) = zero
957 : end do
958 : ABI_FREE(dij_hf)
959 :
960 : !option_dij=3
961 : !call symdij(Cryst%gprimd,Cryst%indsym,ipert0,&
962 : !& Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
963 : !& KS_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
964 : !& Cryst%symafm,Cryst%symrec)
965 : end if
966 :
967 804 : ABI_MALLOC(tmp_kstab, (2, Wfd%nkibz, Wfd%nsppol))
968 4102 : tmp_kstab=0
969 406 : do spin=1,Sigp%nsppol
970 1053 : do ikcalc=1,Sigp%nkptgw ! No spin dependent!
971 647 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
972 647 : tmp_kstab(1, ik_ibz, spin) = Sigp%minbnd(ikcalc, spin)
973 852 : tmp_kstab(2, ik_ibz, spin) = Sigp%maxbnd(ikcalc, spin)
974 : end do
975 : end do
976 :
977 201 : if (dtset%userie == 456) then
978 0 : call wrtout(std_out, "Computing KS vxc matrix elements using supercell WFK file")
979 : call calc_vhxc_me(Wfdf, KS_mflags, KS_me, Cryst, Dtset, nfftf, ngfftf, &
980 : ks_vtrial, ks_vhartr, ks_vxc, Psps, Pawtab, KS_paw_an, Pawang, Pawfgrtab, KS_paw_ij, dijexc_core, &
981 0 : ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur)
982 : else
983 :
984 : call calc_vhxc_me(Wfd, KS_mflags, KS_me, Cryst, Dtset, nfftf, ngfftf, &
985 : ks_vtrial, ks_vhartr, ks_vxc, Psps, Pawtab, KS_paw_an, Pawang, Pawfgrtab, KS_paw_ij, dijexc_core, &
986 201 : ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur)
987 : end if
988 201 : ABI_FREE(tmp_kstab)
989 :
990 : ! Set KS matrix elements connecting different irreps to zero. Do not touch unknown bands!.
991 201 : if (gwcalctyp >= 20 .and. Sigp%symsigma > 0) then
992 0 : bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
993 0 : ABI_MALLOC(ks_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
994 0 : ks_irreptab=0
995 0 : do spin=1,Sigp%nsppol
996 0 : do ikcalc=1,Sigp%nkptgw
997 0 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
998 0 : first_band = Sigp%minbnd(ikcalc,spin)
999 0 : last_band = Sigp%maxbnd(ikcalc,spin)
1000 0 : if (.not. KS_sym(ik_ibz,spin)%failed()) then
1001 0 : ks_irreptab(first_band:last_band,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(first_band:last_band)
1002 : !ks_irreptab(bmin:bmax,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(bmin:bmax)
1003 : end if
1004 : end do
1005 : end do
1006 0 : call KS_me%zero(ks_irreptab)
1007 0 : ABI_FREE(ks_irreptab)
1008 : end if
1009 :
1010 201 : call KS_me%print(header="Matrix elements in the KS basis set", prtvol=Dtset%prtvol)
1011 :
1012 : ! If possible, calculate the EXX energy between the frozen core
1013 : ! and the valence electrons using KS wavefunctions
1014 : ! MG: Be careful here, since ex_energy is meaningful only if all occupied states are calculated.
1015 201 : if (KS_mflags%has_sxcore ==1) then
1016 : ! TODO
1017 : !ex_energy = mels_get_exene_core(KS_me,kmesh,ks_ebands)
1018 1 : ex_energy=zero
1019 2 : do spin=1,Sigp%nsppol
1020 21 : do ik=1,Kmesh%nibz
1021 153 : do ib=b1gw,b2gw
1022 152 : if (Sigp%nsig_ab==1) then
1023 133 : ex_energy = ex_energy + half*ks_ebands%occ(ib,ik,spin)*Kmesh%wt(ik)*KS_me%sxcore(ib,ib,ik,spin)
1024 : else
1025 0 : ex_energy = ex_energy + half*ks_ebands%occ(ib,ik,spin)*Kmesh%wt(ik)*SUM(KS_me%sxcore(ib,ib,ik,:))
1026 : end if
1027 : end do
1028 : end do
1029 : end do
1030 1 : write(msg,'(a,2(es16.6,a))')' CORE Exchange energy with KS wavefunctions: ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
1031 1 : call wrtout(std_out, msg)
1032 : end if
1033 :
1034 201 : call timab(407,2,tsec) ! vHxc_me
1035 201 : call timab(408,1,tsec) ! hqp_init
1036 :
1037 : ! Do not break this coding!
1038 : ! When gwcalctyp>10, the order of the bands can be exchanged after
1039 : ! the diagonalization. Therefore, we have to correctly assign the matrix elements to the corresponding
1040 : ! bands and we cannot skip the following even though it looks useless.
1041 201 : if (gwcalctyp >= 10) call wrtout(std_out, ch10//' *************** KS Energies *******************')
1042 :
1043 : !=== qp_ebands stores energies and occ. used for the calculation ===
1044 : ! * Initialize qp_ebands with KS values.
1045 : ! * In case of SC update qp_ebands using the QPS file.
1046 201 : call ks_ebands%copy(qp_ebands)
1047 :
1048 804 : ABI_MALLOC(qp_rhor, (nfftf, Dtset%nspden))
1049 804 : ABI_MALLOC(qp_taur, (nfftf, Dtset%nspden * Dtset%usekden))
1050 201 : QP_sym => KS_sym
1051 :
1052 201 : if (gwcalctyp<10) then
1053 : ! one-shot GW, just do a copy of the KS density.
1054 1095499 : qp_rhor=ks_rhor
1055 125 : if(Dtset%usekden==1)qp_taur=ks_taur
1056 : QP_sym => KS_sym
1057 : else
1058 : ! Self-consistent GW.
1059 : ! Read the unitary matrix and the QP energies of the previous step from the QPS file.
1060 76 : call QP_energies%init()
1061 76 : QP_energies%e_corepsp=ecore/Cryst%ucvol
1062 :
1063 : ! m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>
1064 : ! Initialize the QP amplitudes with KS wavefunctions.
1065 456 : ABI_MALLOC(Sr%m_ks_to_qp, (Sigp%nbnds, Sigp%nbnds, Kmesh%nibz, Sigp%nsppol))
1066 73347 : Sr%m_ks_to_qp=czero
1067 986 : do ib=1,Sigp%nbnds
1068 7162 : Sr%m_ks_to_qp(ib,ib,:,:)=cone
1069 : end do
1070 :
1071 : ! Now read m_ks_to_qp and update the energies in qp_ebands.
1072 : ! TODO switch on the renormalization of n in sigma.
1073 304 : ABI_MALLOC(prev_rhor, (nfftf, Dtset%nspden))
1074 304 : ABI_MALLOC(prev_taur, (nfftf, Dtset%nspden*Dtset%usekden))
1075 228 : ABI_MALLOC(prev_Pawrhoij, (Cryst%natom*Psps%usepaw))
1076 :
1077 : call rdqps(qp_ebands,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,&
1078 76 : nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,Sr%m_ks_to_qp,prev_rhor,prev_Pawrhoij)
1079 :
1080 : !Find the irreps associated to the QP amplitudes starting from the analogous table for the KS states.
1081 : !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1082 : !allocate(qp_irreptab(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1083 : !qp_irreptab=0
1084 : !!qp_irreptab=ks_irreptab
1085 :
1086 : !do jb_qp=bmin,bmax
1087 : !do ib_ks=bmin,bmax
1088 : !if (ABS(Sr%m_ks_to_qp(ib_ks,jb_qp,ik_ibz,spin)) > tol12) then ! jb_qp has same the same character as ib_ks.
1089 : !ks_irr = ks_irreptab(ib_ks,ib_ks,ik_ibz,spin)
1090 : !qp_irreptab(jb_qp,jb_qp,ik_ibz,spin) = ks_irr
1091 : !do ii=bmin,bmax
1092 : !if (ks_irr == ks_irreptab(ii,ib_ks,ik_ibz,spin)) then
1093 : !qp_irreptab(jb_qp,ii,ik_ibz,spin) = ks_irr
1094 : !end if
1095 : !end do
1096 : !end if
1097 : !end do
1098 : !end do
1099 :
1100 1163604 : if (nscf==0) prev_rhor=ks_rhor
1101 76 : if (nscf==0 .and. Dtset%usekden==1) prev_taur=ks_taur
1102 :
1103 76 : if (nscf>0.and.gwcalctyp>=20.and. wfd%my_rank == master) then
1104 : ! Print the unitary transformation on std_out.
1105 28 : call show_QP(qp_ebands,Sr%m_ks_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp)
1106 : end if
1107 :
1108 : ! Compute QP wfg as linear combination of KS states.
1109 : ! * Wfd%ug is modified inside calc_wf_qp
1110 : ! * For PAW, update also the on-site projections.
1111 : ! * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz
1112 : ! TODO here we should use nbsc instead of nbnds
1113 :
1114 76 : call wfd%rotate(Cryst, Sr%m_ks_to_qp)
1115 :
1116 : ! Reinit the storage mode of Wfd as ug have been changed ===
1117 : ! Update also the wavefunctions for GW corrections on each processor
1118 76 : call wfd%reset_ur_cprj()
1119 :
1120 : ! This test has been disabled (too expensive!)
1121 : if (.False.) call wfd%test_ortho(Cryst, Pawtab, unit=std_out)
1122 :
1123 : ! Compute QP occupation numbers.
1124 76 : call wrtout(std_out,'sigma: calculating QP occupation numbers:')
1125 76 : call qp_ebands%update_occ(Dtset%spinmagntarget, prtvol=0)
1126 76 : qp_vbik(:,:) = qp_ebands%get_valence_idx()
1127 :
1128 : ! Calculate the irreducible representations of the new QP amplitdues.
1129 76 : if (Sigp%symsigma==1.and.gwcalctyp>=20) then
1130 0 : ABI_MALLOC(QP_sym,(Wfd%nkibz,Wfd%nsppol))
1131 0 : use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric
1132 0 : do spin=1,Wfd%nsppol
1133 0 : do ikcalc=1,Sigp%nkptgw
1134 0 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1135 : ! Quick fix for SCGW+symm TODO fix properly!
1136 0 : first_band = Sigp%minbnd(ikcalc,spin)
1137 0 : last_band = Sigp%maxbnd(ikcalc,spin)
1138 : ! first_band = MINVAL(Sigp%minbnd(:,spin))
1139 : ! last_band = MAXVAL(Sigp%maxbnd(:,spin))
1140 : ! call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,qp_ebands,Pawtab,Pawrad,Pawang,Psps,&
1141 : call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,qp_ebands,Pawtab,Pawrad,Pawang,Psps,&
1142 0 : & Dtset%tolsym,QP_sym(ik_ibz,spin))
1143 : end do
1144 : end do
1145 :
1146 : ! Recreate the Sig_ij tables taking advantage of the classification of the bands.
1147 0 : call sigma_tables(Sigp, Kmesh, esymm=QP_sym)
1148 : end if
1149 :
1150 : ! Compute QP density using the updated wfg.
1151 76 : call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, qp_rhor)
1152 76 : if (Dtset%usekden==1) call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, qp_taur, optcalc=1)
1153 :
1154 : ! ========================================
1155 : ! ==== QP self-consistent GW with PAW ====
1156 : ! ========================================
1157 76 : if (Dtset%usepaw==1) then
1158 0 : ABI_MALLOC(qp_nhat,(nfftf,Dtset%nspden))
1159 0 : nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat
1160 0 : ABI_MALLOC(qp_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
1161 :
1162 0 : ABI_MALLOC(QP_pawrhoij,(Cryst%natom))
1163 0 : ABI_MALLOC(QP_paw_ij,(Cryst%natom))
1164 0 : ABI_MALLOC(QP_paw_an,(Cryst%natom))
1165 :
1166 : ! Calculate new QP quantities: nhat, nhatgr, rho_ij, paw_ij, and paw_an.
1167 : call paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,qp_ebands,&
1168 : Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij,&
1169 0 : QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim,qp_nhatgr,compch_sph,compch_fft)
1170 : else
1171 76 : ABI_MALLOC(qp_nhat, (0, 0))
1172 76 : ABI_MALLOC(qp_nhatgr, (0, 0, 0))
1173 76 : ABI_MALLOC(QP_pawrhoij, (0))
1174 76 : ABI_MALLOC(QP_paw_ij, (0))
1175 76 : ABI_MALLOC(QP_paw_an, (0))
1176 : end if
1177 :
1178 : ! here I should renormalize the density
1179 : call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,qp_rhor,Cryst%ucvol,&
1180 76 : Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
1181 :
1182 76 : if (Dtset%usepaw==1) qp_rhor(:,:)=qp_rhor(:,:)+qp_nhat(:,:) ! Add the "hat" term.
1183 :
1184 76 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_rhor,ucvol=ucvol)
1185 76 : if (Dtset%usekden==1) then
1186 0 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_taur,optrhor=1,ucvol=ucvol)
1187 : end if
1188 :
1189 : ! Simple mixing of the PW density to damp oscillations in the Hartree potential.
1190 76 : if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then
1191 12 : write(msg,'(2a,f6.3)')ch10,' sigma: mixing QP densities using rhoqpmix= ',Dtset%rhoqpmix
1192 12 : call wrtout(std_out, msg)
1193 92260 : qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor)
1194 12 : if(Dtset%usekden==1) qp_taur = prev_taur + Dtset%rhoqpmix*(qp_taur-prev_taur) ! mix taur.
1195 : end if
1196 :
1197 76 : ABI_FREE(prev_rhor)
1198 76 : ABI_FREE(prev_taur)
1199 76 : if (Psps%usepaw==1.and.nscf>0) call pawrhoij_free(prev_pawrhoij)
1200 76 : ABI_FREE(prev_pawrhoij)
1201 :
1202 228 : ABI_MALLOC(qp_rhog,(2,nfftf))
1203 76 : call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp5)
1204 :
1205 : ! ===========================================
1206 : ! ==== Optional output of the QP density ====
1207 : ! ===========================================
1208 76 : if (Dtset%prtden/=0 .and. wfd%my_rank == master) then
1209 : call fftdatar_write("qp_rhor",dtfil%fnameabo_qp_den,dtset%iomode,hdr_sigma,&
1210 60 : cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor,mpi_enreg_seq,ebands=qp_ebands)
1211 : end if
1212 :
1213 : ! ===========================================
1214 : ! === Optional output of the full QP density
1215 : ! ===========================================
1216 76 : if (Wfd%usepaw==1.and.Dtset%prtden==2) then
1217 0 : ABI_MALLOC(qp_rhor_paw ,(nfftf,Wfd%nspden))
1218 0 : ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden))
1219 0 : ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden))
1220 :
1221 : call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,&
1222 : Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,QP_pawrhoij,Pawtab,Dtset%prtvol,&
1223 0 : qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
1224 :
1225 0 : ABI_FREE(qp_rhor_n_one)
1226 0 : ABI_FREE(qp_rhor_nt_one)
1227 :
1228 0 : if (Dtset%prtvol > 9) then
1229 : ! Print a normalisation check
1230 0 : norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3))
1231 0 : write(msg,'(a,F8.4)') ' QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm
1232 0 : call wrtout(std_out, msg)
1233 : end if
1234 :
1235 : ! Write the density to file
1236 0 : if (my_rank==master) then
1237 : call fftdatar_write("qp_pawrhor",dtfil%fnameabo_qp_pawden,dtset%iomode,hdr_sigma,&
1238 0 : cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor_paw,mpi_enreg_seq,ebands=qp_ebands)
1239 : end if
1240 0 : ABI_FREE(qp_rhor_paw)
1241 : end if
1242 :
1243 76 : nkxc=0
1244 76 : if (Dtset%nspden==1) nkxc=2
1245 76 : if (Dtset%nspden>=2) nkxc=3 !check GGA and spinor that is messy !!!
1246 : !In case of MGGA, fxc and kxc are not available and we dont need them for the screening part (for now ...)
1247 76 : if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0
1248 76 : if (nkxc/=0) then
1249 304 : ABI_MALLOC(qp_kxc,(nfftf,nkxc))
1250 : end if
1251 :
1252 : ! **** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
1253 76 : n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf
1254 228 : ABI_MALLOC(qp_vhartr, (nfftf))
1255 304 : ABI_MALLOC(qp_vtrial, (nfftf,Dtset%nspden))
1256 228 : ABI_MALLOC(qp_vxc, (nfftf,Dtset%nspden))
1257 :
1258 76 : optene=4; moved_atm_inside=0; moved_rhor=0; istep=1
1259 :
1260 : call setvtr(Cryst%atindx1,Dtset,QP_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
1261 : istep,qp_kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
1262 : Cryst%nattyp,nfftf,ngfftf,ngrvdw,qp_nhat,qp_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
1263 : optene,Pawang,Pawrad,QP_pawrhoij,Pawtab,ph1df,Psps,qp_rhog,qp_rhor,Cryst%rmet,&
1264 : Cryst%rprimd,strsxc,Cryst%ucvol,usexcnhat,qp_vhartr,vpsp,qp_vtrial,qp_vxc,vxcavg_qp,Wvl,&
1265 76 : xccc3d,Cryst%xred,taur=qp_taur)
1266 :
1267 76 : ABI_SFREE(qp_kxc)
1268 :
1269 76 : if (Dtset%usepaw==1) then
1270 0 : call timab(561,1,tsec)
1271 :
1272 : ! Compute QP Dij
1273 : call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,&
1274 : Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
1275 : Dtset%nspden,Cryst%ntypat,QP_paw_an,QP_paw_ij,Pawang,Pawfgrtab,&
1276 : Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
1277 : k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),&
1278 : qp_vtrial,qp_vxc,Cryst%xred,Dtset%znucl,&
1279 0 : nucdipmom=Dtset%nucdipmom,spinaxis=Dtset%spinaxis)
1280 :
1281 : ! Symmetrize total Dij
1282 0 : option_dij=0
1283 :
1284 : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,&
1285 : Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
1286 : QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,&
1287 0 : Cryst%symafm,Cryst%symrec)
1288 :
1289 : ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij.
1290 0 : call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab)
1291 0 : call timab(561,2,tsec)
1292 : end if
1293 :
1294 1967517 : ehartree=half*SUM(qp_rhor(:,1)*qp_vhartr(:))/DBLE(nfftf)*Cryst%ucvol
1295 :
1296 6156 : write(msg,'(a,80a)')ch10,('-',ii=1,80)
1297 76 : call wrtout(ab_out,msg)
1298 76 : write(msg,'(5a,f9.4,3a,es21.14,2a,es21.14)')ch10,&
1299 76 : ' QP results after the unitary transformation in the KS subspace: ',ch10,ch10,&
1300 76 : ' Number of electrons = ',qp_rhog(1,1)*Cryst%ucvol,ch10,ch10,&
1301 76 : ' QP Band energy [Ha] = ',qp_ebands%get_bandenergy(),ch10,&
1302 152 : ' QP Hartree energy [Ha] = ',ehartree
1303 76 : call wrtout(ab_out,msg)
1304 6156 : write(msg,'(a,80a)')ch10,('-',ii=1,80)
1305 152 : call wrtout(ab_out,msg)
1306 :
1307 : ! TODO Since plasmonpole model 2-3-4 depend on the Fourier components of the density
1308 : ! in case of self-consistency we might calculate here the ppm coefficients using qp_rhor
1309 : end if ! gwcalctyp>=10
1310 :
1311 : ! KS hamiltonian: hdft(b1,b1,k,s)= <b1,k,s|H_s|b1,k,s>
1312 106123 : ABI_CALLOC(hdft, (b1gw:b2gw, b1gw:b2gw, Kmesh%nibz, Sigp%nsppol*Sigp%nsig_ab))
1313 :
1314 201 : if (Dtset%nspinor == 1) then
1315 392 : do spin=1,Sigp%nsppol
1316 1593 : do ik=1,Kmesh%nibz
1317 9781 : do ib=b1gw,b2gw
1318 9583 : hdft(ib, ib, ik, spin) = ks_ebands%eig(ib, ik, spin)
1319 : end do
1320 : end do
1321 : end do
1322 : else
1323 : ! Spinorial case
1324 : ! * Note that here vxc contains the contribution of the core.
1325 : ! * Scale ovlp if orthonormalization is not satisfied as npwwfn might be < npwvec.
1326 : ! TODO add spin-orbit case and gwpara 2
1327 21 : ABI_MALLOC(my_band_list, (wfd%mband))
1328 14 : ABI_MALLOC(bmask, (wfd%mband))
1329 281 : bmask = .False.; bmask(b1gw:b2gw) = .True.
1330 :
1331 7 : if (Wfd%usepaw == 1) then
1332 0 : ABI_MALLOC(Cp1,(Wfd%natom, Wfd%nspinor))
1333 0 : call pawcprj_alloc(Cp1, 0, Wfd%nlmn_atm)
1334 : end if
1335 :
1336 14 : do spin=1,Sigp%nsppol
1337 45 : do ik_ibz=1,Kmesh%nibz
1338 : ! Distribute bands in [b1gw, b2gw] range
1339 31 : call wfd%distribute_bands(ik_ibz, spin, my_nband, my_band_list, bmask=bmask)
1340 31 : if (my_nband == 0) cycle
1341 31 : npw_k = Wfd%npwarr(ik_ibz)
1342 440 : do ii=1,my_nband ! ib=b1gw,b2gw in sequential
1343 402 : ib = my_band_list(ii)
1344 402 : ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave, msg) == 0, msg)
1345 402 : ug1 => wave%ug
1346 402 : cdummy = xdotc(npw_k*Wfd%nspinor,ug1,1,ug1,1)
1347 402 : ovlp(1) = REAL(cdummy); ovlp(2) = AIMAG(cdummy)
1348 :
1349 402 : if (Psps%usepaw==1) then
1350 0 : call wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.)
1351 0 : ovlp = ovlp + paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab)
1352 : end if
1353 : ! write(std_out,*)ovlp(1),ovlp(2)
1354 402 : norm=DBLE(ovlp(1)+ovlp(2))
1355 402 : ovlp(1)=DBLE(ovlp(1)/norm); ovlp(2)=DBLE(ovlp(2)/norm) ! ovlp(2)=cone-ovlp(1)
1356 402 : hdft(ib,ib,ik_ibz,1) = ks_ebands%eig(ib,ik_ibz,1)*ovlp(1) - KS_me%vxc(ib,ib,ik_ibz,3)
1357 402 : hdft(ib,ib,ik_ibz,2) = ks_ebands%eig(ib,ik_ibz,1)*ovlp(2) - KS_me%vxc(ib,ib,ik_ibz,4)
1358 402 : hdft(ib,ib,ik_ibz,3) = KS_me%vxc(ib,ib,ik_ibz,3)
1359 433 : hdft(ib,ib,ik_ibz,4) = KS_me%vxc(ib,ib,ik_ibz,4)
1360 : end do
1361 : end do
1362 : end do
1363 :
1364 7 : call xmpi_sum(hdft, wfd%comm, ierr)
1365 :
1366 7 : ABI_FREE(my_band_list)
1367 7 : ABI_FREE(bmask)
1368 14 : if (Wfd%usepaw==1) then
1369 0 : call pawcprj_free(Cp1)
1370 0 : ABI_FREE(Cp1)
1371 : end if
1372 : end if
1373 :
1374 : ! Initialize Sigma results. TODO it is better if we use ragged arrays indexed by the k-point
1375 201 : call Sr%init(Sigp, Kmesh%nibz, Dtset%usepawu)
1376 :
1377 : ! Setup bare Hamiltonian := T + v_{loc} + v_{nl} + v_H.
1378 : !
1379 : ! * The representation depends whether we are updating the wfs or not.
1380 : ! * ks_vUme is zero unless we are using DFT+U as starting point, see calc_vHxc_braket
1381 : ! * Note that vH matrix elements are calculated using the true uncutted interaction
1382 : ! This should be changed if the cutoff is also used in the GS run.
1383 :
1384 201 : if (gwcalctyp < 10) then
1385 : ! For one-shot GW use the KS representation.
1386 52716 : Sr%hhartree = hdft - KS_me%vxcval
1387 :
1388 : ! Additional stuff for PAW
1389 : ! * DFT +U Hamiltonian
1390 : ! * LEXX.
1391 : ! * Core contribution estimated using Fock exchange.
1392 125 : if (Dtset%usepaw==1) then
1393 1090 : if (Sigp%use_sigxcore == 1) Sr%hhartree = hdft - (KS_me%vxc - KS_me%sxcore)
1394 5 : if (Dtset%usepawu /= 0) Sr%hhartree = Sr%hhartree - KS_me%vu
1395 5 : if (Dtset%useexexch /= 0) then
1396 0 : ABI_ERROR("useexexch > 0 not implemented")
1397 0 : Sr%hhartree = Sr%hhartree - KS_me%vlexx
1398 : end if
1399 : end if
1400 :
1401 : else
1402 : ! Self-consistent on energies and|or wavefunctions.
1403 : ! * For NC get the bare Hamiltonian $H_{bare}= T+v_{loc}+ v_{nl}$ in the KS representation
1404 : ! * For PAW, calculate the matrix elements of h0, store also the new Dij in QP_Paw_ij.
1405 : ! * h0 is defined as T+vH[tn+nhat+tnZc] + vxc[tnc] + dij_eff and
1406 : ! dij_eff = dij^0 + dij^hartree + dij^xc-dij^xc_val + dijhat - dijhat_val.
1407 : ! In the above expression tn, tnhat are QP quantities.
1408 76 : if (Dtset%usepaw==0) then
1409 456 : ABI_MALLOC(hbare, (b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab))
1410 52603 : hbare = hdft - KS_me%vhartree - KS_me%vxcval
1411 :
1412 : ! Change basis from KS to QP, hbare is overwritten: A_{QP} = U^\dagger A_{KS} U
1413 380 : ABI_MALLOC(htmp, (b1gw:b2gw,b1gw:b2gw, Kmesh%nibz, Sigp%nsppol*Sigp%nsig_ab))
1414 304 : ABI_MALLOC(ctmp, (b1gw:b2gw, b1gw:b2gw))
1415 228 : ABI_MALLOC(uks2qp, (b1gw:b2gw, b1gw:b2gw))
1416 105054 : htmp = hbare; hbare = czero
1417 :
1418 154 : do spin=1,Sigp%nsppol
1419 607 : do ik=1,Kmesh%nibz
1420 52373 : uks2qp(:,:) = Sr%m_ks_to_qp(b1gw:b2gw,b1gw:b2gw,ik,spin)
1421 984 : do iab=1,Sigp%nsig_ab
1422 453 : is_idx=spin; if (Sigp%nsig_ab>1) is_idx=iab
1423 2820233 : ctmp(:,:)=MATMUL(htmp(:,:,ik,is_idx),uks2qp(:,:))
1424 2820686 : hbare(:,:,ik,is_idx)=MATMUL(TRANSPOSE(CONJG(uks2qp)),ctmp)
1425 : end do
1426 : end do
1427 : end do
1428 76 : ABI_FREE(htmp)
1429 76 : ABI_FREE(ctmp)
1430 76 : ABI_FREE(uks2qp)
1431 : end if ! usepaw==0
1432 :
1433 : ! Calculate the QP matrix elements
1434 : ! This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions.
1435 : ! For PAW, we have to construct the new bare Hamiltonian.
1436 76 : call wrtout(std_out,ch10//' *************** QP Energies *******************')
1437 :
1438 76 : call QP_mflags%reset()
1439 : ! if (gwcalctyp<20) QP_mflags%only_diago=1 ! For e-only, no need of off-diagonal elements.
1440 76 : QP_mflags%has_vhartree=1
1441 76 : if (Dtset%usepaw==1) then
1442 0 : QP_mflags%has_kinetic =1
1443 0 : QP_mflags%has_hbare =1
1444 : end if
1445 : !QP_mflags%has_vxc =1
1446 : !QP_mflags%has_vxcval =1
1447 : !if (Sigp%gwcalctyp >100) QP_mflags%has_vxcval_hybrid=1
1448 76 : if (mod10==5 .and. &
1449 : (Dtset%ixc_sigma==-402 .or. Dtset%ixc_sigma==-406 .or. Dtset%ixc_sigma==-427 .or. Dtset%ixc_sigma==-428 .or. &
1450 : Dtset%ixc_sigma==-456 .or. Dtset%ixc_sigma==41 .or. Dtset%ixc_sigma==42)) then
1451 23 : QP_mflags%has_vxcval_hybrid = 1
1452 : end if
1453 :
1454 : !if (Sigp%use_sigxcore==1) QP_mflags%has_sxcore =1
1455 : !if (Dtset%usepawu/=0) QP_mflags%has_vu =1
1456 : !if (Dtset%useexexch/=0) QP_mflags%has_lexexch=1
1457 :
1458 304 : ABI_MALLOC(tmp_kstab, (2, Wfd%nkibz, Wfd%nsppol))
1459 1513 : tmp_kstab=0
1460 154 : do spin=1,Sigp%nsppol
1461 553 : do ikcalc=1,Sigp%nkptgw ! No spin dependent!
1462 399 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1463 399 : tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin)
1464 477 : tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin)
1465 : end do
1466 : end do
1467 :
1468 : call calc_vhxc_me(Wfd,QP_mflags,QP_me,Cryst,Dtset,nfftf,ngfftf,&
1469 : qp_vtrial,qp_vhartr,qp_vxc,Psps,Pawtab,QP_paw_an,Pawang,Pawfgrtab,QP_paw_ij,dijexc_core,&
1470 76 : qp_rhor,usexcnhat,qp_nhat,qp_nhatgr,nhatgrdim,tmp_kstab,taur=qp_taur)
1471 76 : ABI_FREE(tmp_kstab)
1472 :
1473 76 : if (gwcalctyp>=20 .and. Sigp%symsigma>0) then
1474 0 : bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1475 0 : ABI_MALLOC(qp_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol))
1476 0 : qp_irreptab=0
1477 0 : do spin=1,Sigp%nsppol
1478 0 : do ikcalc=1,Sigp%nkptgw
1479 0 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1480 0 : first_band = Sigp%minbnd(ikcalc,spin)
1481 0 : last_band = Sigp%maxbnd(ikcalc,spin)
1482 0 : if (.not. QP_sym(ik_ibz,spin)%failed()) then
1483 0 : qp_irreptab(first_band:last_band,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(first_band:last_band)
1484 : !qp_irreptab(bmin:bmax,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(bmin:bmax)
1485 : end if
1486 : end do
1487 : end do
1488 0 : call QP_me%zero(qp_irreptab)
1489 0 : ABI_FREE(qp_irreptab)
1490 : end if
1491 :
1492 76 : call QP_me%print(header="Matrix elements in the QP basis set", prtvol=Dtset%prtvol)
1493 :
1494 : ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij.
1495 76 : if (Dtset%usepaw==1) then
1496 0 : call wrtout(std_out," *** After calc_vHxc_braket *** ")
1497 : ! TODO finalize the implementation of this routine.
1498 0 : call paw_ij_print(QP_Paw_ij,unit=std_out,pawprtvol=Dtset%pawprtvol,pawspnorb=Dtset%pawspnorb,mode_paral="COLL")
1499 0 : call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab)
1500 : end if
1501 :
1502 76 : if (Dtset%usepaw==0) then
1503 : ! GA: We have an odd bug here. I have to unroll this loop, otherwise it
1504 : ! might cause segfault when running on several nodes.
1505 : !
1506 : ! Sr%hhartree = hbare + QP_me%vhartree
1507 76 : if (QP_mflags%has_vxcval_hybrid == 0) then
1508 108 : do spin=1,Sigp%nsppol*Sr%nsig_ab
1509 380 : do ikcalc=1,Sr%nkibz
1510 2691 : do ib1=b1gw,b2gw
1511 24256 : do ib2=b1gw,b2gw
1512 23984 : Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + QP_me%vhartree(ib2,ib1,ikcalc,spin)
1513 : end do
1514 : end do
1515 : end do
1516 : end do
1517 : else
1518 46 : do spin=1,Sigp%nsppol*Sr%nsig_ab
1519 227 : do ikcalc=1,Sr%nkibz
1520 2300 : do ib1=b1gw,b2gw
1521 28117 : do ib2=b1gw,b2gw
1522 : Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + &
1523 27936 : QP_me%vhartree(ib2,ib1,ikcalc,spin) + QP_me%vxcval_hybrid(ib2,ib1,ikcalc,spin)
1524 : end do
1525 : end do
1526 : end do
1527 : end do
1528 : end if
1529 : else
1530 0 : Sr%hhartree=QP_me%hbare
1531 : end if
1532 :
1533 76 : if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then
1534 : ! bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw
1535 0 : do spin=1,Sigp%nsppol
1536 0 : do ik_ibz=1,Kmesh%nibz
1537 0 : if (.not. QP_sym(ik_ibz,spin)%failed()) then
1538 0 : bmin=Sigp%minbnd(ik_ibz,spin); bmax=Sigp%minbnd(ik_ibz,spin)
1539 0 : do ib2=bmin,bmax
1540 0 : irr_idx2 = QP_sym(ik_ibz,spin)%b2irrep(ib2)
1541 0 : do ib1=bmin,bmax
1542 0 : irr_idx1 = QP_sym(ik_ibz,spin)%b2irrep(ib1)
1543 0 : if (irr_idx1/=irr_idx2 .and. ALL((/irr_idx1,irr_idx2/)/=0) ) Sr%hhartree(ib1,ib2,ik_ibz,spin) = czero
1544 : end do
1545 : end do
1546 : end if
1547 : end do
1548 : end do
1549 : end if
1550 :
1551 76 : ABI_FREE(qp_rhog)
1552 76 : ABI_FREE(qp_vhartr)
1553 76 : ABI_FREE(qp_vxc)
1554 76 : ABI_FREE(qp_nhat)
1555 76 : ABI_FREE(qp_nhatgr)
1556 76 : call QP_me%free()
1557 : end if ! gwcalctyp<10
1558 :
1559 : ! Free some memory
1560 201 : ABI_SFREE(hbare)
1561 201 : ABI_SFREE(hdft)
1562 :
1563 : ! Prepare the storage of QP amplitudes and energies
1564 : ! Initialize with KS wavefunctions and energies.
1565 : ! FIXME: This array should be allocated only if self-consistent
1566 73472 : if (sr%needs_eigvec_qp) Sr%eigvec_qp=czero
1567 30808 : Sr%en_qp_diago=zero
1568 5828 : do ib=1,Sigp%nbnds
1569 40474 : Sr%en_qp_diago(ib,:,:) = ks_ebands%eig(ib,:,:)
1570 12004 : if (sr%needs_eigvec_qp) Sr%eigvec_qp(ib,ib,:,:) = cone
1571 : end do
1572 :
1573 : ! Store <n,k,s|V_xc[n_val]|n,k,s> and <n,k,s|V_U|n,k,s> ===
1574 : ! Note that we store the matrix elements of V_xc in the KS basis set, not in the QP basis set
1575 : ! Matrix elements of V_U are zero unless we are using DFT+U as starting point
1576 1629 : do ib=b1gw,b2gw
1577 13142 : Sr%vxcme(ib,:,:)=KS_me%vxcval(ib,ib,:,:)
1578 1629 : if (Dtset%usepawu/=0) Sr%vUme (ib,:,:)=KS_me%vu(ib,ib,:,:)
1579 : end do
1580 :
1581 : ! Initial guess for the GW energies
1582 : ! Save the energies of the previous iteration.
1583 406 : do spin=1,Sigp%nsppol
1584 1638 : do ik=1,Kmesh%nibz
1585 30402 : do ib=1,Sigp%nbnds
1586 29170 : Sr%e0 (ib,ik,spin) = qp_ebands%eig(ib,ik,spin)
1587 30402 : Sr%egw(ib,ik,spin) = qp_ebands%eig(ib,ik,spin)
1588 : end do
1589 1232 : Sr%e0gap(ik,spin) = zero
1590 1232 : ks_iv = ks_vbik(ik, spin)
1591 1437 : if (Sigp%nbnds>=ks_iv+1) Sr%e0gap(ik,spin)=Sr%e0(ks_iv+1,ik,spin)-Sr%e0(ks_iv,ik,spin)
1592 : end do
1593 : end do
1594 : !
1595 : !=== If required apply a scissor operator or update the energies ===
1596 : !TODO check if other Sr entries have to be updated
1597 : !moreover this part should be done only in case of semiconductors
1598 : !FIXME To me it makes more sense if we apply the scissor to ks_ebands but I have to RECHECK csigme
1599 201 : if (ABS(Sigp%mbpt_sciss)>tol6) then
1600 2 : write(msg,'(6a,f10.5,2a)')ch10,&
1601 2 : ' sigma : performing a first self-consistency',ch10,&
1602 2 : ' update of the energies in G by a scissor operator',ch10, &
1603 4 : ' applying a scissor operator of ',Sigp%mbpt_sciss*Ha_eV,' [eV] ',ch10
1604 2 : call wrtout(units, msg)
1605 4 : do spin=1,Sigp%nsppol
1606 8 : do ik=1,Kmesh%nibz
1607 4 : ks_iv=ks_vbik(ik,spin)
1608 6 : if (Sigp%nbnds>=ks_iv+1) then
1609 28 : Sr%egw (ks_iv+1:Sigp%nbnds,ik,spin) = Sr%egw (ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss
1610 28 : qp_ebands%eig(ks_iv+1:Sigp%nbnds,ik,spin) = qp_ebands%eig(ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss
1611 : end if
1612 : end do
1613 : end do
1614 : !call apply_scissor(qp_ebands,Sigp%mbpt_sciss)
1615 : else if (.FALSE.) then
1616 : write(msg,'(4a)')ch10,&
1617 : ' sigma : performing a first self-consistency',ch10,&
1618 : ' update of the energies in G by a previous GW calculation'
1619 : call wrtout(units, msg)
1620 : ! TODO Recheck this part, is not clear to me!
1621 : ABI_MALLOC(igwene,(qp_ebands%mband, qp_ebands%nkpt, qp_ebands%nsppol))
1622 : call rdgw(qp_ebands, '__in.gw__', igwene, extrapolate=.TRUE.)
1623 : ABI_FREE(igwene)
1624 : Sr%egw=qp_ebands%eig
1625 : !
1626 : ! * Recalculate the new fermi level.
1627 : call qp_ebands%update_occ(Dtset%spinmagntarget, prtvol=0)
1628 : end if
1629 :
1630 : ! In case of AC refer all the energies wrt to the fermi level
1631 : ! Be careful because results from ppmodel cannot be used for AC
1632 : ! FIXME check ks_energy or qp_energy (in case of SCGW?)
1633 :
1634 201 : if (mod10 == SIG_GW_AC) then
1635 : ! All these quantities will be passed to csigme
1636 : ! if I skipped the self-consistent part then here I have to use fermi
1637 942 : qp_ebands%eig = qp_ebands%eig - qp_ebands%fermie
1638 942 : Sr%egw = Sr%egw - qp_ebands%fermie
1639 942 : Sr%e0 = Sr%e0 - qp_ebands%fermie
1640 9 : old_fermie = qp_ebands%fermie
1641 : ! TODO Recheck fermi
1642 : ! Clean EVERYTHING in particulare the treatment of E fermi
1643 9 : qp_ebands%fermie = zero
1644 : end if
1645 :
1646 : ! Setup frequencies around the KS\QP eigenvalues to compute Sigma derivatives (notice the spin) ===
1647 : ! TODO it is better using an odd Sr%nomega4sd so that the KS\QP eigenvalue is in the middle
1648 201 : ioe0j=Sr%nomega4sd/2+1
1649 406 : do spin=1,Sigp%nsppol
1650 1053 : do ikcalc=1,Sigp%nkptgw
1651 647 : ib1=Sigp%minbnd(ikcalc,spin)
1652 647 : ib2=Sigp%maxbnd(ikcalc,spin)
1653 647 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc))
1654 6303 : do jb=ib1,ib2
1655 25845 : do io=1,Sr%nomega4sd
1656 25198 : Sr%omega4sd(jb,ik_ibz,io,spin)=Sr%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j)
1657 : end do
1658 : end do
1659 : end do
1660 : end do
1661 :
1662 201 : call pstat_proc%print(_PSTAT_ARGS_)
1663 201 : call timab(408,2,tsec) ! hqp_init
1664 201 : call timab(409,1,tsec) ! getW
1665 :
1666 : ! Get epsilon^{-1} either from the _SCR or the _SUSC file and store it in epsm1%epsm1
1667 : ! If epsm1%mqmem==0, allocate and read a single q-slice inside csigme.
1668 : ! TODO epsm1%nomega should be initialized so that only the frequencies really needed are stored in memory
1669 : ! TODO The same piece of code is present in screening.
1670 201 : if (sigp%needs_w()) then
1671 :
1672 327 : select case (dtset%gwgamma)
1673 : case (0)
1674 161 : id_required = 4; ikxc = 0; approx_type = 0; option_test = 0; dim_kxcg = 0
1675 322 : ABI_MALLOC(kxcg, (nfftf_tot,dim_kxcg))
1676 :
1677 : case (1, 2)
1678 : ! ALDA TDDFT kernel vertex
1679 1 : ABI_CHECK(epsm1%ID == 0, "epsm1%ID should be 0")
1680 :
1681 1 : if (Dtset%usepaw==1) then
1682 : ! If we have PAW, we need the full density on the fine grid
1683 0 : ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1684 0 : if (Dtset%getpawden==0 .and. Dtset%irdpawden==0) then
1685 0 : ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1686 : end if
1687 0 : call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1688 0 : if (.not. file_exists(dtfil%filpawdensin)) then
1689 0 : ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1690 : end if
1691 :
1692 0 : ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1693 : call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1694 0 : ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1695 :
1696 0 : call Hdr_rhor%free()
1697 0 : call pawrhoij_free(tmp_pawrhoij)
1698 0 : ABI_FREE(tmp_pawrhoij)
1699 : end if ! usepaw==1
1700 :
1701 1 : id_required=4; ikxc=7; approx_type=1; dim_kxcg=1
1702 1 : if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1703 1 : if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1704 3 : ABI_MALLOC(kxcg, (nfftf_tot, dim_kxcg))
1705 :
1706 1 : dbg_mode=.FALSE.
1707 1 : if (Dtset%usepaw==1) then
1708 : ! Use PAW all-electron density
1709 : call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_aepaw_rhor,&
1710 0 : epsm1%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode)
1711 : else
1712 : ! Norm-conserving
1713 : call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_rhor,&
1714 1 : epsm1%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode)
1715 : end if
1716 :
1717 : case (3, 4)
1718 : ! ADA non-local kernel vertex
1719 0 : ABI_CHECK(epsm1%ID==0, "epsm1%ID should be 0")
1720 0 : ABI_CHECK(Sigp%nsppol==1, "ADA vertex for GWGamma not available yet for spin-polarised cases")
1721 :
1722 0 : if (Dtset%usepaw==1) then
1723 : ! If we have PAW, we need the full density on the fine grid
1724 0 : ABI_MALLOC(ks_aepaw_rhor,(nfftf, Sigp%nsppol))
1725 0 : if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1726 0 : ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1727 : end if
1728 0 : call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1729 0 : if (.not. file_exists(dtfil%filpawdensin)) then
1730 0 : ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1731 : end if
1732 :
1733 0 : ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1734 :
1735 : call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1736 0 : ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1737 :
1738 0 : call Hdr_rhor%free()
1739 0 : call pawrhoij_free(tmp_pawrhoij)
1740 0 : ABI_FREE(tmp_pawrhoij)
1741 : end if ! Dtset%usepaw==1
1742 :
1743 0 : id_required=4; ikxc=7; approx_type=2;
1744 0 : if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma
1745 0 : if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only
1746 0 : ABI_MALLOC(fxc_ADA, (epsm1%npwe, epsm1%npwe, epsm1%nqibz))
1747 : ! Use userrd to set kappa
1748 0 : if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp
1749 : ! Set correct value of kappa (should be scaled with alpha*r_s where)
1750 : ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3)
1751 0 : rhoav = (drude_plsmf*drude_plsmf)/four_pi
1752 0 : r_s = (three/(four_pi*rhoav))**third
1753 0 : alpha = (four*ninth*piinv)**third
1754 0 : Dtset%userrd = Dtset%userrd/(alpha*r_s)
1755 :
1756 0 : dbg_mode=.TRUE.
1757 0 : if (Dtset%usepaw==1) then
1758 : ! Use PAW all-electron density
1759 : call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,&
1760 : ks_aepaw_rhor,epsm1%npwe,epsm1%nqibz,epsm1%qibz,&
1761 0 : fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode)
1762 : else
1763 : ! Norm conserving
1764 : call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,&
1765 : ks_rhor,epsm1%npwe,epsm1%nqibz,epsm1%qibz,&
1766 0 : fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode)
1767 : end if
1768 :
1769 0 : dim_kxcg = 0
1770 0 : ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1771 :
1772 : case (-3, -4, -5, -6, -7, -8)
1773 : ! bootstrap kernels
1774 3 : ABI_CHECK(epsm1%ID==0,"epsm1%ID should be 0")
1775 :
1776 3 : if (Dtset%usepaw==1) then
1777 : ! If we have PAW, we need the full density on the fine grid
1778 0 : ABI_MALLOC(ks_aepaw_rhor, (nfftf,Wfd%nspden))
1779 0 : if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1780 0 : ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1781 : end if
1782 0 : call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1783 0 : if (.not. file_exists(dtfil%filpawdensin)) then
1784 0 : ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1785 : end if
1786 :
1787 0 : ABI_MALLOC(tmp_pawrhoij, (cryst%natom*wfd%usepaw))
1788 :
1789 : call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1790 0 : ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1791 :
1792 0 : call Hdr_rhor%free()
1793 0 : call pawrhoij_free(tmp_pawrhoij)
1794 0 : ABI_FREE(tmp_pawrhoij)
1795 : end if ! Dtset%usepaw==1
1796 :
1797 3 : id_required=4; ikxc=7; dim_kxcg=0
1798 :
1799 3 : if (dtset%gwgamma>-5) then
1800 2 : approx_type=4 ! full fxc(G,G')
1801 1 : else if (dtset%gwgamma>-7) then
1802 1 : approx_type=5 ! fxc(0,0) one-shot
1803 : else
1804 0 : approx_type=6 ! rpa-type bootstrap
1805 : end if
1806 :
1807 3 : option_test=MOD(Dtset%gwgamma,2)
1808 : ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1809 : ! 0 -> TESTPARTICLE, vertex in chi0 only
1810 6 : ABI_MALLOC(kxcg, (nfftf_tot,dim_kxcg))
1811 :
1812 : case (-11)
1813 : ! LR+ALDA kernel
1814 1 : ABI_CHECK(epsm1%ID==0,"epsm1%ID should be 0")
1815 :
1816 1 : if (Dtset%usepaw==1) then
1817 : ! If we have PAW, we need the full density on the fine grid
1818 0 : ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden))
1819 0 : if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then
1820 0 : ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!")
1821 : end if
1822 0 : call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin))
1823 0 : if (.not. file_exists(dtfil%filpawdensin)) then
1824 0 : ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin))
1825 : end if
1826 :
1827 0 : ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw))
1828 :
1829 : call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1830 0 : ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm)
1831 :
1832 0 : call hdr_rhor%free()
1833 0 : call pawrhoij_free(tmp_pawrhoij)
1834 0 : ABI_FREE(tmp_pawrhoij)
1835 : end if ! Dtset%usepaw==1
1836 :
1837 1 : id_required=4; ikxc=7; dim_kxcg=1
1838 1 : approx_type=7
1839 1 : option_test=1
1840 : ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma
1841 3 : ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg))
1842 :
1843 : case default
1844 166 : ABI_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma)))
1845 : end select
1846 :
1847 : ! Set plasma frequency
1848 166 : my_plsmf = drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf = Dtset%ppmfrq
1849 166 : Dtset%ppmfrq = my_plsmf
1850 :
1851 166 : if (Dtset%gwgamma < 3) then
1852 : call epsm1%mkdump(Vcp,epsm1%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,&
1853 : approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,&
1854 166 : nfftf_tot,ngfftf,comm)
1855 : else
1856 : call epsm1%mkdump(Vcp,epsm1%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,&
1857 : approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,&
1858 0 : nfftf_tot,ngfftf,comm,fxc_ADA=fxc_ADA)
1859 : end if
1860 166 : ABI_SFREE(kxcg)
1861 166 : ABI_SFREE(fxc_ADA)
1862 166 : ABI_SFREE(ks_aepaw_rhor)
1863 : end if
1864 : !
1865 : ! ================================================
1866 : ! ==== Calculate plasmonpole model parameters ====
1867 : ! ================================================
1868 : ! TODO Maybe its better if we use mqmem as input variable
1869 201 : use_aerhor=0
1870 402 : ABI_MALLOC(ks_aepaw_rhor, (nfftf,Wfd%nspden*use_aerhor))
1871 :
1872 201 : if (sigp%needs_ppm()) then
1873 : ! If epsm1 is MPI-shared, we have to start the RMA epoch. Note that epsm1%epsm1 is read-only.
1874 130 : if (epsm1%use_mpi_shared_win) then
1875 24 : call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_win, ierr)
1876 24 : ABI_CHECK_MPI(ierr, "")
1877 : end if
1878 :
1879 130 : my_plsmf = drude_plsmf; if (Dtset%ppmfrq > tol6) my_plsmf = Dtset%ppmfrq
1880 130 : call PPm%init(epsm1%mqmem, epsm1%nqibz, epsm1%npwe, Sigp%ppmodel, my_plsmf, Dtset%gw_invalid_freq)
1881 : ! PPm%force_plsmf= force_ppmfrq ! this line to change the plasma frequency in the HL expression.
1882 :
1883 130 : if (Wfd%usepaw==1 .and. Ppm%userho==1) then
1884 : ! For PAW and ppmodel 2-3-4 we need the AE rho(G) without compensation charge.
1885 : ! It would be possible to calculate rho(G) using Paw_pwff, though. It should be faster but
1886 : ! results will depend on the expression used for the matrix elements. This approach is safer.
1887 0 : use_aerhor=1
1888 0 : ABI_FREE(ks_aepaw_rhor)
1889 0 : ABI_MALLOC(ks_aepaw_rhor, (nfftf,Wfd%nspden))
1890 :
1891 : ! Check if input density file is available, otherwise compute
1892 0 : pawden_fname = strcat(Dtfil%filnam_ds(3), '_PAWDEN')
1893 0 : call wrtout(std_out,sjoin('Checking for existence of:',pawden_fname))
1894 0 : if (file_exists(pawden_fname)) then
1895 : ! Read density from file
1896 0 : ABI_MALLOC(tmp_pawrhoij, (cryst%natom*wfd%usepaw))
1897 :
1898 : call read_rhor(pawden_fname, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, &
1899 0 : ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm)
1900 :
1901 0 : call Hdr_rhor%free()
1902 0 : call pawrhoij_free(tmp_pawrhoij)
1903 0 : ABI_FREE(tmp_pawrhoij)
1904 : else
1905 : ! Have to calculate PAW AW rhor from scratch
1906 0 : ABI_MALLOC(qp_rhor_n_one, (pawfgr%nfft,Dtset%nspden))
1907 0 : ABI_MALLOC(qp_rhor_nt_one, (pawfgr%nfft,Dtset%nspden))
1908 :
1909 : ! FIXME
1910 0 : ABI_WARNING(" denfgr in sigma seems to produce wrong results")
1911 0 : write(std_out,*)" input tilde ks_rhor integrates: ",SUM(ks_rhor(:,1))*Cryst%ucvol/nfftf
1912 :
1913 : call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,ks_nhat,Wfd%nspinor,&
1914 : Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_Pawrhoij,Pawtab,Dtset%prtvol, &
1915 0 : ks_rhor,ks_aepaw_rhor,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
1916 :
1917 0 : ABI_FREE(qp_rhor_n_one)
1918 0 : ABI_FREE(qp_rhor_nt_one)
1919 : end if
1920 :
1921 0 : write(msg,'(a,f8.4)')' sigma: PAW AE density used for PPmodel integrates to: ',SUM(ks_aepaw_rhor(:,1))*Cryst%ucvol/nfftf
1922 0 : call wrtout(std_out, msg)
1923 :
1924 0 : if (epsm1%mqmem /= 0) then
1925 : ! Calculate ppmodel parameters for all q-points.
1926 0 : call PPm%setup(Cryst,Qmesh,epsm1%npwe,epsm1%nomega,epsm1%omega,epsm1%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_aepaw_rhor(:,1))
1927 : end if
1928 :
1929 : else
1930 : ! NC or PAW with PPmodel 1.
1931 130 : if (epsm1%mqmem /= 0) then
1932 : ! Calculate ppmodel parameters for all q-points
1933 113 : call PPm%setup(Cryst, Qmesh, epsm1%npwe, epsm1%nomega, epsm1%omega, epsm1%epsm1, nfftf, Gsph_c%gvec, ngfftf, ks_rhor(:,1))
1934 : end if
1935 : end if ! PAW or NC PPm and/or needs density
1936 :
1937 : ! If epsm1 is MPI-shared, we have to close the RMA epoch.
1938 130 : if (epsm1%use_mpi_shared_win) then
1939 24 : call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr)
1940 24 : ABI_CHECK_MPI(ierr, "")
1941 : end if
1942 : end if ! sigma_needs_ppm
1943 :
1944 201 : call pstat_proc%print(_PSTAT_ARGS_)
1945 201 : call timab(409,2,tsec) ! getW
1946 :
1947 201 : if (wfd%my_rank == master) then
1948 : ! Write info on the run on ab_out, then open files to store final results.
1949 173 : call ks_ebands%report_gap(header='KS Band Gaps',unit=ab_out)
1950 173 : if(dtset%ucrpa == 0) call write_sigma_header(Sigp,epsm1,Cryst,Kmesh,Qmesh)
1951 :
1952 : ! unt_gw: File with GW corrections.
1953 : ! unt_sig: Self-energy as a function of frequency.
1954 : ! unt_sgr: Derivative wrt omega of the Self-energy.
1955 : ! unt_sigc: Sigma_c(eik) MRM
1956 : ! unt_sgm: Sigma on the Matsubara axis (imag axis)
1957 :
1958 173 : if (open_file(Dtfil%fnameabo_gw, msg, unit=unt_gw, status='unknown', form='formatted') /= 0) then
1959 0 : ABI_ERROR(msg)
1960 : end if
1961 173 : write(unt_gw,"(a)")"# QP energies E in eV"
1962 173 : write(unt_gw,"(a)")"# Format:"
1963 173 : write(unt_gw,"(a)")"# kpoint"
1964 173 : write(unt_gw,"(a)")"# number of bands computed"
1965 173 : write(unt_gw,"(2a)")"# band index, Re(E), E-E0, Im(E)", ch10
1966 173 : write(unt_gw,"(2(i0,1x),a)")Sigp%nkptgw,Sigp%nsppol, "# nkptgw, nsppol"
1967 :
1968 173 : if (open_file(Dtfil%fnameabo_gwdiag, msg, unit=unt_gwdiag, status='unknown', form='formatted') /= 0) then
1969 0 : ABI_ERROR(msg)
1970 : end if
1971 173 : write(unt_gwdiag,*)Sigp%nkptgw,Sigp%nsppol
1972 :
1973 173 : if (open_file(Dtfil%fnameabo_sig, msg, unit=unt_sig, status='unknown', form='formatted') /= 0) then
1974 0 : ABI_ERROR(msg)
1975 : end if
1976 173 : write(unt_sig,"(a)")"# Sigma_xc and spectral function A along the real frequency axis in eV units"
1977 173 : write(unt_sig,"(a)")"# Format:"
1978 173 : write(unt_sig,"(a)")"# kpoint"
1979 173 : write(unt_sig,"(a)")"# min_band(k) max_band(k)"
1980 173 : write(unt_sig,"(a)")"# For each frequency w:"
1981 173 : write(unt_sig,"(2a)")"# w, {Re(Sigma_b(w)), Im(Sigma_b(w), A_b(w) for b in [min_band, max_band]}",ch10
1982 :
1983 173 : if (open_file(Dtfil%fnameabo_sgr, msg, unit=unt_sgr, status='unknown', form='formatted') /= 0) then
1984 0 : ABI_ERROR(msg)
1985 : end if
1986 173 : write(unt_sgr,"(a)")"# Derivatives of Sigma_c(omega) wrt omega in eV units"
1987 :
1988 : ! Sigma_c(eik) MRM
1989 173 : if (open_file(trim(Dtfil%fnameabo_sgr)//'_SIGC', msg, unit=unt_sigc, status='unknown', form='formatted') /= 0) then
1990 0 : ABI_ERROR(msg)
1991 : end if
1992 :
1993 173 : if (mod10 == SIG_GW_AC) then
1994 : ! Sigma along the imaginary axis.
1995 9 : if (open_file(Dtfil%fnameabo_sgm, msg, unit=unt_sgm, status='unknown', form='formatted') /= 0) then
1996 0 : ABI_ERROR(msg)
1997 : end if
1998 9 : write(unt_sgm,"(a)")"# Sigma_xc along the imaginary frequency axis in eV units"
1999 : end if
2000 : end if
2001 :
2002 : !=======================================================================
2003 : !==== Calculate self-energy and output the results for each k-point ====
2004 : !=======================================================================
2005 : ! Here it would be possible to calculate the QP correction for the same k-point using a PPmodel
2006 : ! in the first iteration just to improve the initial guess and CD or AC in the second step. Useful
2007 : ! if the KS level is far from the expected QP results. Previously it was possible to do such
2008 : ! calculation by simply listing the same k-point twice in kptgw. Now this trick is not allowed anymore.
2009 : ! Everything, indeed, should be done in a clean and transparent way inside csigme.
2010 :
2011 402 : call wfd%print([std_out])
2012 201 : call wrtout(std_out, sigma_type_from_key(mod10))
2013 :
2014 201 : if (gwcalctyp < 10) then
2015 125 : msg = " Perturbative Calculation"
2016 125 : if (gwcalctyp == 1) msg = " Newton Raphson method "
2017 76 : else if (gwcalctyp < 20) then
2018 11 : msg = " Self-Consistent on Energies only"
2019 65 : else if (gwcalctyp == 21) then
2020 5 : msg = " GW 1-RDM Corrections"
2021 : else
2022 60 : msg = " Self-Consistent on Energies and Wavefunctions"
2023 : end if
2024 201 : if (Dtset%ucrpa == 0) call wrtout(std_out,msg)
2025 :
2026 : !=================================================
2027 : !==== Calculate the matrix elements of Sigma =====
2028 : !=================================================
2029 :
2030 201 : nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
2031 :
2032 : ! min and max band indices for GW corrections.
2033 201 : ib1=Sigp%minbdgw; ib2=Sigp%maxbdgw
2034 :
2035 : !MG TODO: I don't like the fact that ib1 and ib2 are redefined here because this
2036 : ! prevents me from refactoring the code. In particular I want to store the self-energy
2037 : ! results inside the sigma_results datatypes hence one needs to know all the dimensions
2038 : ! at the beginning of the execution (e.g. in setup_sigma) so that one can easily allocate the arrays in the type.
2039 201 : if (Dtset%ucrpa >= 1) then
2040 : ! Read the band
2041 0 : if (dtset%plowan_compute<10)then
2042 0 : if (open_file("forlb.ovlp",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then
2043 0 : ABI_ERROR(msg)
2044 : end if
2045 0 : rewind(temp_unt)
2046 0 : read(temp_unt,*)
2047 0 : read(temp_unt,*)
2048 0 : read(temp_unt,*) msg, ib1, ib2
2049 0 : close(temp_unt)
2050 : else
2051 0 : if (open_file("data.plowann",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then
2052 0 : ABI_ERROR(msg)
2053 : end if
2054 0 : rewind(temp_unt)
2055 0 : read(temp_unt,*)
2056 0 : read(temp_unt,*)
2057 0 : read(temp_unt,'(a7,2i4)') msg, ib1, ib2
2058 0 : close(temp_unt)
2059 : endif
2060 : endif
2061 :
2062 : ! Do not store it for rdm_update because it might be too large!
2063 201 : if (.not. rdm_update) then
2064 285117 : ABI_CALLOC(sigcme, (nomega_sigc,ib1:ib2,ib1:ib2,Sigp%nkptgw,Sigp%nsppol*Sigp%nsig_ab))
2065 : endif
2066 :
2067 : ! if (.False. .and. psps%usepaw == 0 .and. wfd%nspinor == 1 .and. any(dtset%so_psp /= 0)) then
2068 : ! call wrtout(std_out, "Computing SOC contribution with first-order perturbation theory")
2069 : ! ABI_MALLOC(bks_mask, (wfd%mband, wfd%nkibz, wfd%nsppol))
2070 : ! bks_mask = .False.
2071 : ! do spin=1,wfd%nsppol
2072 : ! do ikcalc=1,Sigp%nkptgw
2073 : ! ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
2074 : ! ii=Sigp%minbnd(ikcalc, spin); jj=Sigp%maxbnd(ikcalc, spin)
2075 : ! bks_mask(ii:jj, ik_ibz, spin) = .True.
2076 : ! end do
2077 : ! end do
2078 : !
2079 : ! call wfd_get_socpert(wfd, cryst, psps, pawtab, bks_mask, osoc_bks)
2080 : ! ABI_FREE(bks_mask)
2081 : ! ABI_FREE(osoc_bks)
2082 : ! end if
2083 :
2084 : !==========================================================
2085 : !==== Exchange part using the dense gwx_ngfft FFT mesh ====
2086 : !==========================================================
2087 : !TODO : load distribution is not optimal if band parallelism is used.
2088 : !but this problem was also affecting the old implementation.
2089 201 : call timab(421,1,tsec) ! calc_sigx_me
2090 :
2091 : ! This routine shows how the wavefunctions are distributed.
2092 : !call wfd_show_bkstab(Wfd,unit=std_out)
2093 :
2094 201 : if (Dtset%ucrpa>=1) then
2095 : ! Calculation of the Rho_n for calc_ucrpa
2096 0 : call wrtout(std_out, sjoin("begin of Ucrpa calc for a nb of kpoints of: ",itoa(Sigp%nkptgw)))
2097 : ! Wannier basis: rhot1_q_m will need an atom index in the near future.
2098 0 : ic=0
2099 0 : do itypat=1,cryst%ntypat
2100 0 : if(pawtab(itypat)%lpawu.ne.-1) then
2101 0 : ndim=2*dtset%lpawu(itypat)+1
2102 0 : ic=ic+1
2103 0 : itypatcor=itypat
2104 0 : lcor=dtset%lpawu(itypat)
2105 : end if
2106 : end do
2107 0 : if(ic>1) then
2108 0 : ABI_ERROR("number of correlated species is larger than one")
2109 : end if
2110 0 : ABI_MALLOC(rhot1_q_m, (cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz))
2111 0 : ABI_MALLOC(M1_q_m, (cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz))
2112 :
2113 0 : M1_q_m=czero
2114 0 : rhot1_q_m=czero
2115 :
2116 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2117 : !Initialization of wan objects and getting psichies
2118 : !Allocation of rhot1
2119 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2120 0 : if (dtset%plowan_compute>=10)then
2121 0 : call wrtout(units, " cRPA calculations using wannier weights from data.plowann")
2122 : call init_plowannier(dtset%plowan_bandf,dtset%plowan_bandi,dtset%plowan_compute,dtset%plowan_iatom,&
2123 : dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,dtset%plowan_nbl,dtset%plowan_nt,&
2124 : dtset%plowan_projcalc,dtset%acell_orig,dtset%kptns,sum(dtset%plowan_nbl),dtset%nimage,dtset%nkpt,dtset%nspinor,&
2125 0 : dtset%nsppol,dtset%wtk,dtset%dmft_t2g,wanibz_in)
2126 0 : call get_plowannier(wanibz_in,wanibz,dtset)
2127 0 : call fullbz_plowannier(dtset,kmesh,cryst,pawang,wanibz,wanbz)
2128 0 : ABI_MALLOC(rhot1,(sigp%npwx,Qmesh%nibz))
2129 0 : do pwx=1,sigp%npwx
2130 0 : do ibz=1,Qmesh%nibz
2131 0 : call init_operwan_realspace(wanbz,rhot1(pwx,ibz))
2132 0 : call zero_operwan_realspace(wanbz,rhot1(pwx,ibz))
2133 : end do
2134 : end do
2135 : endif
2136 :
2137 : !do ikcalc=1,Sigp%nkptgw
2138 : !if(cryst%nsym==1) nkcalc=Kmesh%nbz
2139 : !if(cryst%nsym>1) nkcalc=Kmesh%nibz
2140 0 : nkcalc=Kmesh%nbz
2141 : !if(1==1)then!DEBUG
2142 : !open(67,file="test.rhot1",status="REPLACE")
2143 0 : do ikcalc=1,nkcalc ! for the oscillator strength, spins are identical without SOC
2144 : ! if(Sigp%nkptgw/=Kmesh%nbz) then
2145 : ! write(msg,'(6a)')ch10,&
2146 : !& ' nkptgw and nbz differs: this is not allowed to compute U in cRPA '
2147 : ! ABI_ERROR(msg)
2148 : ! endif
2149 0 : write(std_out,*) " ikcalc",ikcalc,size(Sigp%kptgw2bz),nkcalc
2150 :
2151 0 : if(mod(ikcalc-1,nprocs)==Wfd%my_rank) then
2152 0 : if(cryst%nsym==1) ik_ibz=Kmesh%tab(ikcalc) ! Index of the irred k-point for GW
2153 0 : if(cryst%nsym>1) ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2154 :
2155 : call prep_calc_ucrpa(ik_ibz,ikcalc,itypatcor,ib1,ib2,Cryst,qp_ebands,Sigp,Gsph_x,Vcp,&
2156 : & Kmesh,Qmesh,lcor,M1_q_m,Pawtab,Pawang,Paw_pwff,&
2157 : & Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,gwx_ngfft,&
2158 0 : & ngfftf,Dtset%prtvol,Dtset%pawcross,Dtset%plowan_compute,rhot1_q_m,wanbz,rhot1)
2159 : end if
2160 : end do
2161 0 : if (dtset%plowan_compute<10)then
2162 0 : call xmpi_sum(rhot1_q_m,Wfd%comm,ierr)
2163 0 : call xmpi_sum(M1_q_m,Wfd%comm,ierr)
2164 0 : M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol
2165 0 : rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol
2166 : !do pwx=1,sigp%npwx; do ibz=1,Qmesh%nibz; do ispinor1=1,dtset%nspinor; do ispinor2=1,dtset%nspinor; do iatom1=1,cryst%nattyp(itypatcor)
2167 : ! do im1=1,2*lcor+1; do im2=1,2*lcor+1
2168 : ! write(67,*)ibz,im1,im2,rhot1_q_m(iatom1,ispinor1,ispinor2,im1,im2,pwx,ibz)
2169 : !enddo; enddo; enddo; enddo; enddo; enddo; enddo
2170 : else
2171 0 : call reduce_operwan_realspace(wanbz,rhot1,sigp%npwx,Qmesh%nibz,Wfd%comm,Kmesh%nbz,Wfd%nsppol)
2172 : !call cwtime(cpu,wall,gflops,"start") !reduction of rhot1
2173 : ! dim=0
2174 : ! do pwx=1,sigp%npwx; do ibz=1,Qmesh%nibz
2175 : ! do spin=1,wanbz%nsppol; do ispinor1=1,wanbz%nspinor; ispinor2=1,wanbz%nspinor
2176 : ! do iatom1=1,wanbz%natom_wan; do iatom2=1,wanbz%natom_wan
2177 : ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1); do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2178 : ! do il1=1,wanbz%nbl_atom_wan(iatom1); do il2=1,wanbz%nbl_atom_wan(iatom2)
2179 : ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2180 : ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2181 : ! dim=dim+1
2182 : ! enddo!im2
2183 : ! enddo!im1
2184 : ! enddo!il2
2185 : ! enddo!il1
2186 : ! enddo!pos2
2187 : ! enddo!pos1
2188 : ! enddo!iatom2
2189 : ! enddo!iatom1
2190 : ! enddo!ispinor2
2191 : ! enddo!ispinor1
2192 : ! enddo!spin
2193 : ! enddo!ibz
2194 : ! enddo!pwx
2195 : ! ABI_MALLOC(buffer,(dim))
2196 : ! nnn=0
2197 : ! do pwx=1,sigp%npwx; do ibz=1,Qmesh%nibz; do spin=1,wanbz%nsppol; do ispinor1=1,wanbz%nspinor; do ispinor2=1,wanbz%nspinor
2198 : ! do iatom1=1,wanbz%natom_wan; do iatom2=1,wanbz%natom_wan
2199 : ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1); do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2200 : ! do il1=1,wanbz%nbl_atom_wan(iatom1); do il2=1,wanbz%nbl_atom_wan(iatom2)
2201 : ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2202 : ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2203 : ! nnn=nnn+1
2204 : ! buffer(nnn)=rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)
2205 : ! enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
2206 : ! call xmpi_sum(buffer,Wfd%comm,ierr)
2207 : ! buffer=buffer/Kmesh%nbz/Wfd%nsppol
2208 : ! nnn=0
2209 : ! do pwx=1,sigp%npwx; do ibz=1,Qmesh%nibz; do spin=1,wanbz%nsppol; do ispinor1=1,wanbz%nspinor; do ispinor2=1,wanbz%nspinor
2210 : ! do iatom1=1,wanbz%natom_wan; do iatom2=1,wanbz%natom_wan
2211 : ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1); do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2212 : ! do il1=1,wanbz%nbl_atom_wan(iatom1); do il2=1,wanbz%nbl_atom_wan(iatom2)
2213 : ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2214 : ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2215 : ! nnn=nnn+1
2216 : ! rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)=buffer(nnn)
2217 : ! !write(67,*)ibz,im1,im2,rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)
2218 : ! enddo; enddo; enddo;enddo; enddo; enddo;enddo; enddo; enddo; enddo; enddo; enddo; enddo
2219 : ! ABI_FREE(buffer)
2220 : endif
2221 :
2222 : !close(67)
2223 : ! call xmpi_barrier(Wfd%comm)
2224 : ! call cwtime(cpu,wall,gflops,"stop")!reduction of rhot1
2225 : ! write(6,*)cpu,wall,gflops
2226 : ! if(Cryst%nsym==1) then
2227 : ! M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol
2228 : ! rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol
2229 : ! endif
2230 : ! Calculation of U in cRPA: need to treat with a different cutoff the
2231 : ! bare coulomb and the screened coulomb interaction.
2232 : ! else !DEBUG
2233 : ! write(6,*)"DEBUGGING calc_ucrpa"
2234 : ! open(67,file="test.rhot1",status="OLD")
2235 : ! rewind(67)
2236 : ! do pwx=1,sigp%npwx
2237 : ! do ibz=1,Qmesh%nibz
2238 : ! do spin=1,wanbz%nsppol
2239 : ! do ispinor1=1,wanbz%nspinor
2240 : ! do ispinor2=1,wanbz%nspinor
2241 : ! do iatom1=1,wanbz%natom_wan
2242 : ! do iatom2=1,wanbz%natom_wan
2243 : ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
2244 : ! do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
2245 : ! do il1=1,wanbz%nbl_atom_wan(iatom1)
2246 : ! do il2=1,wanbz%nbl_atom_wan(iatom2)
2247 : ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
2248 : ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
2249 : ! read(67,*)dummy,dummy,dummy,xx
2250 : ! rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)=xx
2251 : ! enddo!im2
2252 : ! enddo!im1
2253 : ! enddo!il2
2254 : ! enddo!il1
2255 : ! enddo!pos2
2256 : ! enddo!pos1
2257 : ! enddo!iatom2
2258 : ! enddo!iatom1
2259 : ! enddo!ispinor2
2260 : ! enddo!ispinor1
2261 : ! enddo!spin
2262 : ! enddo!ibz
2263 : ! enddo!pwx
2264 : ! close(67)
2265 : ! endif!DEBUG
2266 : call calc_ucrpa(itypatcor,cryst,Kmesh,lcor,M1_q_m,Qmesh,epsm1%npwe,sigp%npwx,&
2267 : Cryst%nsym,Sigp%nomegasr,Sigp%minomega_r,Sigp%maxomega_r,ib1,ib2,&
2268 0 : 'Gsum',Cryst%ucvol,Wfd,epsm1%fname,dtset%plowan_compute,rhot1,wanbz)
2269 :
2270 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2271 : !Deallocation of wan and rhot1
2272 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
2273 0 : if (dtset%plowan_compute >=10) then
2274 0 : do pwx=1,sigp%npwx
2275 0 : do ibz=1,Qmesh%nibz
2276 0 : call destroy_operwan_realspace(wanbz,rhot1(pwx,ibz))
2277 : end do
2278 : end do
2279 0 : ABI_FREE(rhot1)
2280 0 : call destroy_plowannier(wanbz)
2281 : endif
2282 0 : ABI_FREE(rhot1_q_m)
2283 0 : ABI_FREE(M1_q_m)
2284 :
2285 : else
2286 : ! MRM: sigmak_todo is an array indicating whether the k-point must be computed for GW correction.
2287 : ! This array is set to DO IT (to 1) for any GW calculation
2288 : ! and it will only change to NOT DO IT (to 0) in case that a density matrix update calc. is used and we are reading checkpoint files
2289 603 : ABI_MALLOC(sigmak_todo,(Wfd%nkibz))
2290 1423 : sigmak_todo(:)=1
2291 :
2292 : ! ======================================================
2293 : ! ==== Initialize arrays for density matrix updates ====
2294 : ! ======================================================
2295 201 : if (rdm_update) then
2296 : ! Note: all subroutines of 70_gw/m_gwrdm.F90 are implemented assuming nsppol == 1
2297 5 : ABI_CHECK(dtset%nsppol == 1, "1-RDM GW correction only implemented for restricted closed-shell calculations!")
2298 :
2299 2225 : ABI_CALLOC(nateigv, (Wfd%mband, Wfd%mband, Wfd%nkibz, Sigp%nsppol))
2300 290 : ABI_CALLOC(nat_occs, (Wfd%mband, Wfd%nkibz))
2301 2119 : ABI_CALLOC(xrdm_k_full, (b1gw:b2gw, b1gw:b2gw, Wfd%nkibz))
2302 :
2303 5 : write(msg,'(a34,2i9)')' Bands used for the GW 1RDM arrays',b1gw,b2gw
2304 5 : call wrtout(units, msg)
2305 :
2306 35 : do ik_ibz=1,Wfd%nkibz
2307 264 : do ib=b1gw,b2gw
2308 264 : xrdm_k_full(ib,ib,ik_ibz) = qp_ebands%occ(ib,ik_ibz,1)
2309 : end do
2310 275 : do ib=1,Wfd%mband
2311 : ! Copy initial occ numbers (in principle 2 or 0 from KS-DFT)
2312 240 : nat_occs(ib,ik_ibz) = qp_ebands%occ(ib,ik_ibz,1)
2313 : ! Set to identity matrix
2314 270 : nateigv(ib,ib,ik_ibz,1) = cone
2315 : end do
2316 : end do
2317 5 : if (readchkprdm) then
2318 1 : gw1rdm_fname = trim(dtfil%fnameabi_chkp_rdm)
2319 1 : call get_chkprdm(Wfd,Kmesh,Sigp,qp_ebands,nat_occs,nateigv,sigmak_todo,my_rank,gw1rdm_fname)
2320 : end if
2321 5 : call xmpi_barrier(Wfd%comm)
2322 :
2323 : ! Prepare arrays for the imaginary freq. integration (quadrature) of Sigma_c(iw)
2324 15 : ABI_MALLOC(weights, (Sigp%nomegasi))
2325 5 : call quadrature_sigma_cw(Sigp,Sr,weights)
2326 : end if
2327 :
2328 : ! ===================================
2329 : ! ==== Static part (Sigma_x-Vxc) ====
2330 : ! ===================================
2331 840 : do ikcalc=1,Sigp%nkptgw
2332 : ! Index of the irred k-point for GW
2333 639 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
2334 639 : call pstat_proc%print(_PSTAT_ARGS_)
2335 :
2336 : ! Do not compute MELS if the k-point was read from the checkpoint file
2337 : ! this IF only affects GW density matrix update!
2338 :
2339 840 : if (sigmak_todo(ik_ibz) == 1) then
2340 : ! min and max band indices for GW corrections (for this k-point)
2341 1280 : ib1 = MINVAL(Sigp%minbnd(ikcalc,:))
2342 1280 : ib2 = MAXVAL(Sigp%maxbnd(ikcalc,:))
2343 : call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,qp_ebands,dtset, Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,Ltg_k(ikcalc),&
2344 : Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2345 636 : gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty)
2346 :
2347 636 : if (rdm_update) then
2348 : ! Only recompute exchange update to the 1-RDM if the point read was broken or not precomputed
2349 :
2350 : ! Compute Sigma_x - Vxc or DELTA Sigma_x - Vxc
2351 : ! where DELTA Sigma_x = Sigma_x - hyb_parameter Vx^exact for hyb Functionals.
2352 : ! NB: Only restricted closed-shell calcs are implemented here
2353 1790 : ABI_CALLOC(pot_k, (ib1:ib2, ib1:ib2))
2354 1763 : ABI_CALLOC(rdm_k, (ib1:ib2, ib1:ib2))
2355 1709 : pot_k(ib1:ib2,ib1:ib2) = Sr%x_mat(ib1:ib2,ib1:ib2,ik_ibz,1) - KS_me%vxcval(ib1:ib2,ib1:ib2,ik_ibz,1)
2356 :
2357 27 : call calc_rdmx(ib1, ib2, ik_ibz, pot_k, rdm_k, qp_ebands)
2358 :
2359 : ! Update the full 1RDM with the exchange corrected one for this k-point
2360 1709 : xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz) = xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz) + rdm_k(ib1:ib2,ib1:ib2)
2361 :
2362 : ! Compute NAT ORBS for exchange corrected 1-RDM
2363 223 : do ib=ib1,ib2
2364 223 : rdm_k(ib,ib) = rdm_k(ib,ib) + qp_ebands%occ(ib,ik_ibz,1) ! Only restricted closed-shell calcs
2365 : end do
2366 27 : call natoccs(ib1, ib2, rdm_k, nateigv, nat_occs, qp_ebands, ik_ibz, iinfo=0) ! Only restricted closed-shell calcs
2367 27 : ABI_FREE(pot_k)
2368 27 : ABI_FREE(rdm_k)
2369 : end if
2370 : else
2371 3 : write(msg,'(a1)') ' '
2372 3 : call wrtout(std_out,msg)
2373 3 : write(msg,'(a64,i5)') ' Skipping the calc. of Sigma_x - ( Vxc + hyb*Fock) for k-point: ',ik_ibz
2374 3 : call wrtout(std_out,msg)
2375 3 : write(msg,'(a1)') ' '
2376 3 : call wrtout(std_out,msg)
2377 : end if
2378 : end do ! ikcalc
2379 :
2380 : ! for the time being, do not remove this barrier!
2381 201 : call xmpi_barrier(Wfd%comm)
2382 201 : call pstat_proc%print(_PSTAT_ARGS_)
2383 201 : call timab(421,2,tsec) ! calc_sigx_me
2384 :
2385 : ! ==========================================================
2386 : ! ==== Correlation part using the coarse gwc_ngfft mesh ====
2387 : ! ==========================================================
2388 201 : if (mod10/=SIG_HF) then
2389 635 : do ikcalc=1,Sigp%nkptgw
2390 469 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2391 944 : ib1 = MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2392 944 : ib2 = MAXVAL(Sigp%maxbnd(ikcalc,:))
2393 :
2394 274770 : ABI_CALLOC(sigcme_k, (nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab))
2395 :
2396 469 : if (any(mod10 == [SIG_SEX, SIG_COHSEX])) then
2397 : ! Calculate static COHSEX or SEX using the coarse gwc_ngfft mesh.
2398 : call cohsex_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,dtset, Cryst,qp_ebands,Sigp,Sr,epsm1,Gsph_c,Vcp,Kmesh,Qmesh,&
2399 : Ltg_k(ikcalc),Pawtab,Pawang,Paw_pwff,Psps,Wfd,QP_sym,&
2400 5 : gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_k)
2401 : else
2402 : ! Compute correlated part using the coarse gwc_ngfft mesh.
2403 464 : if (x1rdm/=1 .and. sigmak_todo(ik_ibz)==1) then
2404 : ! Do not compute correlation MELS if the k-point was read from the checkpoint file
2405 : ! this IF only affects GW density matrix update
2406 : call calc_sigc_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Dtset,dtfil, Cryst,qp_ebands, &
2407 : Sigp,Sr,epsm1,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,&
2408 : Ltg_k(ikcalc),PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2409 455 : gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_k)
2410 : else
2411 9 : write(msg,'(a1)') ' '
2412 9 : call wrtout(std_out, msg)
2413 9 : write(msg,'(a44,i5)') ' Skipping the calc. of Sigma_c for k-point: ',ik_ibz
2414 9 : call wrtout(std_out, msg)
2415 9 : write(msg,'(a1)') ' '
2416 9 : call wrtout(std_out, msg)
2417 : end if
2418 : end if
2419 :
2420 : ! MRM: compute density matrix numerical correction (freq. integration G0 Sigma_c(iw) G0).
2421 469 : if (rdm_update) then
2422 30 : if (sigmak_todo(ik_ibz) == 1) then
2423 : ! Only recompute correlation update to the 1-RDM if the point read was broken or not precomputed
2424 27 : if (x1rdm /= 1) then
2425 1334 : ABI_CALLOC(rdm_k, (ib1:ib2,ib1:ib2))
2426 21 : call calc_rdmc(ib1, ib2, ik_ibz, Sr%omega_i, weights, sigcme_k, qp_ebands, rdm_k)
2427 : ! Update the full 1RDM with the GW corrected one for this k-point
2428 : ! Only restricted closed-shell calcs
2429 1271 : rdm_k(ib1:ib2,ib1:ib2) = xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz) + rdm_k(ib1:ib2,ib1:ib2)
2430 : ! Compute nat orbs and occ numbers at k-point ik_ibz
2431 21 : call natoccs(ib1, ib2, rdm_k, nateigv, nat_occs, qp_ebands, ik_ibz, iinfo=1)
2432 21 : ABI_FREE(rdm_k)
2433 : endif
2434 : ! Print the checkpoint file if required
2435 27 : if (prtchkprdm) then
2436 6 : gw1rdm_fname = trim(dtfil%fnameabo_chkp_rdm)
2437 6 : call print_chkprdm(Wfd,nat_occs,nateigv,ik_ibz,my_rank,gw1rdm_fname)
2438 : end if
2439 : end if
2440 : else
2441 213867 : sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) = sigcme_k
2442 : end if
2443 635 : ABI_FREE(sigcme_k)
2444 : end do
2445 166 : call xmpi_barrier(Wfd%comm)
2446 166 : if (rdm_update) then
2447 5 : ABI_FREE(xrdm_k_full)
2448 5 : ABI_FREE(weights)
2449 : end if
2450 : end if
2451 :
2452 201 : call xmpi_barrier(Wfd%comm)
2453 201 : ABI_FREE(sigmak_todo)
2454 :
2455 : ! 1) Print WFK and DEN files.
2456 : ! 2) Build band corrections and compute new energies.
2457 201 : if (rdm_update) then
2458 19025 : ABI_CALLOC(gw_rhor, (nfftf, dtset%nspden))
2459 : !
2460 : ! NRM WARNING: only the master has bands on Wfd_nato_master so it prints everything and computes gw_rhor
2461 : !
2462 : ! All procs. update the qp_ebands and the Hdr_sigma
2463 5 : call update_hdr_bst(Wfd_nato_master, nat_occs, b1gw, b2gw, qp_ebands, Hdr_sigma, Dtset%ngfft(1:3))
2464 :
2465 : ! Compute unit cell (averaged) occ = \sum _k weight_k occ_k
2466 5 : call print_tot_occ(qp_ebands)
2467 :
2468 5 : if (my_rank == master) then
2469 5 : call Wfd_nato_master%rotate(cryst, nateigv, bmask=bdm_mask) ! Let it use bdm_mask and build NOs
2470 5 : call Wfd_nato_master%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, gw_rhor) ! Construct the density
2471 5 : if (dtset%prtwf == 1) then
2472 : ! Print WFK file, here qp_ebands contains nat. orb. occs.
2473 5 : call Wfd_nato_master%write_wfk(Hdr_sigma, qp_ebands, dtfil%fnameabo_wfk, wfknocheck=.True.)
2474 : end if
2475 5 : if (dtset%prtden == 1) then
2476 : ! Print DEN file
2477 : call fftdatar_write("density",dtfil%fnameabo_den,dtset%iomode,Hdr_sigma,&
2478 5 : Cryst,ngfftf,cplex1,nfftf,dtset%nspden,gw_rhor,mpi_enreg_seq,ebands=qp_ebands)
2479 : end if
2480 : end if
2481 5 : call xmpi_bcast(gw_rhor,master, Wfd%comm, ierr)
2482 :
2483 : ! We no longer need Wfd_nato_master.
2484 5 : call Wfd_nato_master%free()
2485 5 : call xmpi_barrier(Wfd%comm)
2486 :
2487 5 : ABI_FREE(bdm_mask) ! The master already used bdm_mask
2488 5 : ABI_FREE(nat_occs) ! Occs were already placed in qp_ebands
2489 5 : call epsm1%free() ! We no longer need epsm1 for GW@KS-DFT 1RDM but we may need space on the RAM memory
2490 :
2491 5 : if (gw1rdm==2 .and. Sigp%nkptgw==Wfd%nkibz) then
2492 : ! Compute energies only if all k-points are available
2493 : ! We need the hole 1-RDM to build Fock[GW.1RDM]!
2494 232 : ABI_CALLOC(old_ks_purex, (b1gw:b2gw, Sigp%nkptgw))
2495 228 : ABI_CALLOC(new_hartr, (b1gw:b2gw, Sigp%nkptgw))
2496 33012 : ABI_CALLOC(gw_rhog, (2, nfftf))
2497 11012 : ABI_CALLOC(gw_vhartr, (nfftf))
2498 : !
2499 : ! A) Compute Evext = int rho(r) vext(r) dr -> simply dot product on the FFT grid
2500 : ! Only restricted closed-shell calcs
2501 : !
2502 11004 : den_int = sum(gw_rhor(:,1)) * ucvol / nfftf
2503 11004 : evext_energy = sum(gw_rhor(:,1) * vpsp(:)) * ucvol / nfftf
2504 : !
2505 : ! B) Coulomb <KS_i|Vh[NO]|KS_j>
2506 : !
2507 : ! FFT to build gw_rhog
2508 4 : call fourdp(1, gw_rhog, gw_rhor(:,1), -1, MPI_enreg_seq, nfftf, ndat1, ngfftf, tim_fourdp5)
2509 4 : ecutf = dtset%ecut
2510 4 : if (psps%usepaw == 1) then
2511 0 : ecutf = dtset%pawecutdg
2512 0 : call wrtout(std_out, ch10//' FFT (fine) grid used in PAW GW update:')
2513 : end if
2514 4 : call getcut(boxcut, ecutf, gmet, gsqcut, dtset%iboxcut, std_out, k0, ngfftf)
2515 : call hartre(1, gsqcut, dtset%icutcoul, psps%usepaw, MPI_enreg_seq, nfftf, ngfftf, dtset%nkpt, dtset%rcut, &
2516 4 : gw_rhog, cryst%rprimd, dtset%vcutgeo, gw_vhartr)
2517 :
2518 : ! Build Vhartree -> gw_vhartr
2519 92 : ABI_ICALLOC(tmp_kstab, (2,Wfd%nkibz,Wfd%nsppol))
2520 8 : do spin=1,Sigp%nsppol
2521 32 : do ikcalc=1,Sigp%nkptgw
2522 24 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
2523 24 : tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin)
2524 28 : tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin)
2525 : end do
2526 : end do
2527 :
2528 : ! Build matrix elements from gw_vhartr -> GW1RDM_me
2529 : call calc_vhxc_me(Wfd,KS_mflags,GW1RDM_me,Cryst,Dtset,nfftf,ngfftf,&
2530 : ks_vtrial,gw_vhartr,ks_vxc,Psps,Pawtab,KS_paw_an,Pawang,Pawfgrtab,KS_paw_ij,dijexc_core,&
2531 4 : gw_rhor,usexcnhat,ks_nhat,ks_nhatgr,nhatgrdim,tmp_kstab,taur=ks_taur)
2532 :
2533 4 : call xmpi_barrier(Wfd%comm)
2534 4 : ABI_FREE(tmp_kstab)
2535 4 : ABI_FREE(gw_rhor)
2536 4 : ABI_FREE(gw_rhog)
2537 4 : ABI_FREE(gw_vhartr)
2538 :
2539 : ! Save new <i|Hartree[NO]|i> in KS basis for Delta eik
2540 28 : do ikcalc=1,Sigp%nkptgw
2541 24 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
2542 48 : ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2543 48 : ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2544 220 : do ib=b1gw,b2gw
2545 216 : new_hartr(ib,ikcalc)=GW1RDM_me%vhartree(ib,ib,ik_ibz,1)
2546 : end do
2547 : end do
2548 : !
2549 : ! C) Exchange <KS_i|hyb*K^RANGE-SEP?[KS]|KS_j> and save it in old_ks_purex
2550 : !
2551 : ! Build Vcp_ks and Vcp_full
2552 4 : call setup_vcp(Vcp_ks,Vcp_full,Dtset,Gsph_x,Gsph_c,Cryst,Qmesh,Kmesh,coef_hyb,comm)
2553 4 : call xmpi_barrier(Wfd%comm)
2554 4 : if (coef_hyb > tol8) then
2555 0 : do ikcalc=1,Sigp%nkptgw
2556 0 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
2557 0 : ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2558 0 : ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2559 : ! Notice we need KS occs and band energy diffs
2560 : call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,ks_ebands,dtset, Sigp,Sr,Gsph_x,Vcp_ks,Kmesh,Qmesh,Ltg_k(ikcalc),&
2561 : Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,&
2562 0 : gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty)
2563 :
2564 : ! Build <KS_i|RS?_Hyb?_Sigma_x[KS]|KS_j> matrix (Use Wfd)
2565 0 : call xmpi_barrier(Wfd%comm)
2566 0 : do ib=b1gw,b2gw
2567 : ! Save old alpha*<i|K[KS]|i> from the GS calc. for Delta eik
2568 0 : old_ks_purex(ib,ikcalc)=Sr%x_mat(ib,ib,ik_ibz,1)
2569 : end do
2570 : end do
2571 : endif
2572 4 : call xmpi_barrier(Wfd%comm)
2573 4 : call Vcp_ks%free()
2574 :
2575 : ! Use the Wfd with 'new name' (Wfd_nato_all) because it will contain the GW 1RDM nat. orbs. (bands)
2576 4 : ABI_COMMENT("From now on, the Wfd bands will contain the nat. orbs. ones")
2577 4 : Wfd_nato_all => Wfd
2578 : ! Let rotate build the NOs in Wfd_nato_all (KS->NO)
2579 4 : call Wfd_nato_all%rotate(Cryst, nateigv)
2580 4 : call xmpi_barrier(Wfd%comm)
2581 : !
2582 : ! D) Non-local <NO_i|Vnl|NO_i> terms [saved on nl_bks(band,k,spin)]
2583 : !
2584 4 : call Wfd_nato_all%get_nl_me(Cryst,Psps,Pawtab,bdm2_mask,nl_bks)
2585 4 : evextnl_energy=zero
2586 8 : do spin=1,Sigp%nsppol
2587 32 : do ikcalc=1,Sigp%nkptgw
2588 220 : do ib=Sr%b1gw,Sr%b2gw
2589 216 : evextnl_energy=evextnl_energy+Kmesh%wt(ikcalc)*qp_ebands%occ(ib,ikcalc,spin)*nl_bks(ib,ikcalc,spin)
2590 : end do
2591 : end do
2592 : end do
2593 4 : call xmpi_barrier(Wfd%comm)
2594 4 : ABI_FREE(nl_bks)
2595 4 : ABI_FREE(bdm2_mask)
2596 : !
2597 : ! E) Exchange <NO_i|K[NO]|NO_j>
2598 : !
2599 : ! Get the matrix elements and save them on Sr%x_mat
2600 4 : tol_empty=0.0001 ! Allow lower occ numbers
2601 28 : do ikcalc=1,Sigp%nkptgw
2602 24 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
2603 48 : ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2604 48 : ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2605 :
2606 : ! Build <NO_i|Sigma_x[NO]|NO_j> matrix
2607 : call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,qp_ebands,dtset, Sigp,Sr,Gsph_x,Vcp_full,Kmesh,Qmesh,Ltg_k(ikcalc),&
2608 : Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd_nato_all,Wfdf,QP_sym,&
2609 28 : gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty)
2610 : end do
2611 4 : tol_empty=0.01 ! Recover standard value for tolerance on the occ numbers
2612 4 : call xmpi_barrier(Wfd%comm)
2613 4 : call Vcp_full%free()
2614 : ! Save the new total exchange energy Ex = Ex[GW.1RDM]
2615 4 : ex_energy = Sr%get_exene(Kmesh, qp_ebands)
2616 : ! Save the new total exchange-correlation MBB energy Exc = Exc^MBB[GW.1RDM]
2617 4 : exc_mbb_energy = Sr%get_excene(Kmesh,qp_ebands)
2618 :
2619 : ! Transform <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j>,
2620 : ! <KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j>,
2621 : ! and <KS_i|T|KS_j> -> <NO_i|T|NO_j>
2622 : !
2623 4 : call change_matrix(Sigp,Sr,GW1RDM_me,Kmesh,nateigv)
2624 4 : call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data
2625 :
2626 : ! Print Delta eik and band (state) energies
2627 : !
2628 4 : call print_band_energies(b1gw,b2gw,Sr,Sigp,KS_me,Kmesh,ks_ebands,new_hartr,old_ks_purex)
2629 4 : call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data
2630 : !
2631 : ! Print the updated total energy and all energy components
2632 : !
2633 4 : eh_energy = Sr%get_haene(GW1RDM_me,Kmesh,qp_ebands)
2634 4 : ekin_energy = Sr%get_kiene(GW1RDM_me,Kmesh,qp_ebands)
2635 : ! SD 2-RDM
2636 4 : etot_sd=ekin_energy+evext_energy+evextnl_energy+QP_energies%e_corepsp+QP_energies%e_ewald+eh_energy+ex_energy
2637 : ! MBB 2-RDM
2638 4 : etot_mbb=ekin_energy+evext_energy+evextnl_energy+QP_energies%e_corepsp+QP_energies%e_ewald+eh_energy+exc_mbb_energy
2639 : call print_total_energy(ekin_energy,evext_energy,evextnl_energy,QP_energies%e_corepsp,eh_energy,ex_energy,&
2640 4 : exc_mbb_energy,QP_energies%e_ewald,etot_sd,etot_mbb,den_int)
2641 : !
2642 : ! Clean GW1RDM_me and allocated arrays
2643 : !
2644 4 : call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data
2645 4 : call GW1RDM_me%free() ! Deallocate GW1RD_me
2646 4 : ABI_FREE(old_ks_purex)
2647 8 : ABI_FREE(new_hartr)
2648 : else
2649 1 : ABI_FREE(gw_rhor)
2650 1 : ABI_FREE(bdm2_mask)
2651 : end if
2652 5 : ABI_FREE(nateigv)
2653 : endif
2654 :
2655 201 : call xmpi_barrier(Wfd%comm)
2656 :
2657 : ! MRM: skip the rest for rdm_update=true (density matrix update)
2658 201 : if (.not.rdm_update) then
2659 : ! =====================================================
2660 : ! ==== Solve Dyson equation storing results in Sr% ====
2661 : ! =====================================================
2662 : ! * Use perturbative approach or AC to find QP corrections.
2663 : ! * If qp-GW, diagonalize also H0+Sigma in the KS basis set to get the
2664 : ! new QP amplitudes and energies (Sr%eigvec_qp and Sr%en_qp_diago.
2665 : ! TODO AC with spinor not implemented yet.
2666 : ! TODO Diagonalization of Sigma+hhartre with AC is wrong.
2667 : !
2668 196 : call timab(425,1,tsec) ! solve_dyson
2669 805 : do ikcalc=1,Sigp%nkptgw
2670 609 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW
2671 1226 : ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
2672 1226 : ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
2673 :
2674 3654 : ABI_MALLOC(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab))
2675 274882 : sigcme_k=sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:)
2676 :
2677 : call solve_dyson(ikcalc, ib1, ib2, nomega_sigc, dtset, Sigp, Kmesh, sigcme_k, qp_ebands%eig, &
2678 609 : Sr, ks_me, Dtfil, Wfd%comm)
2679 609 : ABI_FREE(sigcme_k)
2680 : !
2681 : ! Calculate direct gap for each spin and print out final results.
2682 : ! We use the valence index of the KS system because we still do not know
2683 : ! all the QP corrections. Ideally one should use the QP valence index
2684 1226 : do spin=1,Sigp%nsppol
2685 617 : if (Sigp%maxbnd(ikcalc,spin) >= ks_vbik(ik_ibz,spin)+1 .and. &
2686 609 : Sigp%minbnd(ikcalc,spin) <= ks_vbik(ik_ibz,spin) ) then
2687 597 : ks_iv=ks_vbik(ik_ibz,spin)
2688 597 : Sr%egwgap (ik_ibz,spin)= Sr%egw(ks_iv+1,ik_ibz,spin) - Sr%egw(ks_iv,ik_ibz,spin)
2689 597 : Sr%degwgap(ik_ibz,spin)= Sr%degw(ks_iv+1,ik_ibz,spin) - Sr%degw(ks_iv,ik_ibz,spin)
2690 : else
2691 : ! The "gap" cannot be computed
2692 20 : Sr%e0gap(ik_ibz,spin)=zero; Sr%egwgap (ik_ibz,spin)=zero; Sr%degwgap(ik_ibz,spin)=zero
2693 : end if
2694 : end do
2695 :
2696 805 : if (wfd%my_rank == master) call Sr%write_results(ikcalc, ik_ibz, Sigp, ks_ebands)
2697 : end do !ikcalc
2698 :
2699 196 : call timab(425,2,tsec) ! solve_dyson
2700 196 : call timab(426,1,tsec) ! finalize
2701 :
2702 : ! Update the energies in qp_ebands
2703 : ! If QPSCGW, use diagonalized eigenvalues otherwise perturbative results.
2704 196 : if (gwcalctyp>=10) then
2705 941 : do ib=1,Sigp%nbnds
2706 6837 : qp_ebands%eig(ib,:,:)=Sr%en_qp_diago(ib,:,:)
2707 : end do
2708 : else
2709 25070 : qp_ebands%eig=Sr%egw
2710 : end if
2711 :
2712 : ! ================================================================================
2713 : ! ==== This part is done only if all k-points in the IBZ have been calculated ====
2714 : ! ================================================================================
2715 196 : if (Sigp%nkptgw==Kmesh%nibz) then
2716 :
2717 : ! Recalculate new occupations and Fermi level.
2718 86 : call qp_ebands%update_occ(Dtset%spinmagntarget,prtvol=Dtset%prtvol)
2719 86 : qp_vbik(:,:) = qp_ebands%get_valence_idx()
2720 :
2721 86 : write(msg,'(2a,3x,2(es16.6,a))')ch10,' New Fermi energy : ',qp_ebands%fermie,' Ha ,',qp_ebands%fermie*Ha_eV,' eV'
2722 86 : call wrtout(units, msg)
2723 :
2724 : ! === If all k-points and all occupied bands are calculated, output EXX ===
2725 : ! FIXME here be careful about the check on ks_vbik in case of metals
2726 : ! if (my_rank==master.and.Sigp%nkptgw==Kmesh%nibz.and.ALL(Sigp%minbnd(:)==1).and.ALL(Sigp%maxbnd(:)>=MAXVAL(nbv(:)))) then
2727 : ! if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=MAXVAL(MAXVAL(ks_vbik(:,:),DIM=1))) ) then
2728 1164 : if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=ks_vbik) ) then ! FIXME here the two arrays use a different indexing.
2729 :
2730 80 : ex_energy = Sr%get_exene(Kmesh,qp_ebands)
2731 80 : write(msg,'(a,2(es16.6,a))')' New Exchange energy : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
2732 80 : call wrtout(units, msg)
2733 : end if
2734 :
2735 : ! Report the QP gaps (Fundamental and direct)
2736 86 : call qp_ebands%report_gap(header='QP Band Gaps',unit=ab_out)
2737 :
2738 : ! Band structure interpolation from QP energies computed on the k-mesh.
2739 1146 : if (nint(dtset%einterp(1)) /= 0 .and. all(sigp%minbdgw == sigp%minbnd) .and. all(sigp%maxbdgw == sigp%maxbnd)) then
2740 3 : call qp_ebands%interpolate_kpath(dtset, cryst, [sigp%minbdgw, sigp%maxbdgw], dtfil%filnam_ds(4), comm)
2741 : end if
2742 : end if ! Sigp%nkptgw==Kmesh%nibz
2743 : !
2744 : ! Write SCF data in case of self-consistent calculation ===
2745 : ! Save Sr%en_qp_diago, Sr%eigvec_qp and m_ks_to_qp in the _QPS file.
2746 : ! Note that in the first iteration qp_rhor contains KS rhor, then the mixed rhor.
2747 196 : if (gwcalctyp>=10) then
2748 : ! Calculate the new m_ks_to_qp
2749 71 : call updt_m_ks_to_qp(Sigp,Kmesh,nscf,Sr,Sr%m_ks_to_qp)
2750 :
2751 71 : if (wfd%my_rank == master) then
2752 : call wrqps(Dtfil%fnameabo_qps,Sigp,Cryst,Kmesh,Psps,Pawtab,QP_Pawrhoij,&
2753 55 : Dtset%nspden,nscf,nfftf,ngfftf,Sr,qp_ebands,Sr%m_ks_to_qp,qp_rhor)
2754 : end if
2755 :
2756 : ! Report the MAX variation for each kptgw and spin
2757 71 : call wrtout(ab_out,ch10//' Convergence of QP corrections ')
2758 144 : do spin=1,Sigp%nsppol
2759 73 : write(msg,'(a,i2,a)')' >>>>> For spin ',spin,' <<<<< '
2760 73 : call wrtout(ab_out, msg)
2761 513 : do ikcalc=1,Sigp%nkptgw
2762 369 : ib1 = Sigp%minbnd(ikcalc,spin); ib2 = Sigp%maxbnd(ikcalc,spin)
2763 369 : ik_bz = Sigp%kptgw2bz(ikcalc); ik_ibz = Kmesh%tab(ik_bz)
2764 4143 : ii = imax_loc(ABS(Sr%degw(ib1:ib2,ik_ibz,spin)))
2765 369 : max_degw = Sr%degw(ii,ik_ibz,spin)
2766 369 : write(msg,('(a,i3,a,2f8.3,a,i3)'))'. kptgw no:',ikcalc,'; Maximum DeltaE = (',max_degw*Ha_eV,') for band index:',ii
2767 442 : call wrtout(ab_out, msg)
2768 : end do
2769 : end do
2770 : end if
2771 :
2772 : ! ============================================
2773 : ! ==== Save the GW results in NETCDF file ====
2774 : ! ============================================
2775 196 : if (wfd%my_rank == master) then
2776 168 : NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), '_SIGRES.nc'), xmpi_comm_self))
2777 336 : NCF_CHECK(nctk_defnwrite_ivars(ncid, ["sigres_version"], [1]))
2778 168 : NCF_CHECK(cryst%ncwrite(ncid))
2779 168 : NCF_CHECK(ks_ebands%ncwrite(ncid))
2780 168 : NCF_CHECK(Sr%ncwrite(Sigp, epsm1, ncid)) ! WARNING!! If gw1rdm>0 then epsm1 is no longer present!!
2781 : ! Add qp_rhor. Note that qp_rhor == ks_rhor if wavefunctions are not updated.
2782 : !ncerr = nctk_write_datar("qp_rhor",path,ngfft,cplex,nfft,nspden,&
2783 : ! comm_fft,fftn3_distrib,ffti3_local,datar,action)
2784 168 : NCF_CHECK(nf90_close(ncid))
2785 :
2786 : !if (Sigp%nkptgw==Wfd%nkibz) then
2787 : ! call write_qpdata(dtset, sr)
2788 : !end if
2789 : end if
2790 : end if ! MRM skipped if GW density matrix update
2791 : end if ! ucrpa
2792 :
2793 : !----------------------------- END OF THE CALCULATION ------------------------
2794 : !
2795 : !=====================
2796 : !==== Close Files ====
2797 : !=====================
2798 201 : if (wfd%my_rank == master) then
2799 173 : close(unt_gw); close(unt_gwdiag); close(unt_sig); close(unt_sgr); close(unt_sigc)
2800 173 : if (mod10==SIG_GW_AC) close(unt_sgm)
2801 : end if
2802 : !
2803 : !===============================================
2804 : !==== Free arrays and local data structures ====
2805 : !===============================================
2806 201 : ABI_FREE(ks_vbik)
2807 201 : ABI_FREE(qp_vbik)
2808 201 : ABI_FREE(ph1d)
2809 201 : ABI_FREE(ph1df)
2810 201 : ABI_FREE(qp_rhor)
2811 201 : ABI_FREE(ks_rhor)
2812 201 : ABI_FREE(ks_rhog)
2813 201 : ABI_FREE(qp_taur)
2814 201 : ABI_FREE(ks_taur)
2815 201 : ABI_FREE(ks_vhartr)
2816 201 : ABI_FREE(ks_vtrial)
2817 201 : ABI_FREE(vpsp)
2818 201 : ABI_FREE(ks_vxc)
2819 201 : ABI_FREE(xccc3d)
2820 201 : ABI_FREE(grchempottn)
2821 201 : ABI_FREE(grewtn)
2822 201 : ABI_FREE(grvdw)
2823 201 : ABI_SFREE(sigcme)
2824 201 : ABI_SFREE(kxc)
2825 201 : ABI_SFREE(qp_vtrial)
2826 201 : ABI_FREE(ks_nhat)
2827 201 : ABI_FREE(ks_nhatgr)
2828 201 : ABI_FREE(dijexc_core)
2829 201 : ABI_FREE(ks_aepaw_rhor)
2830 201 : call pawfgr_destroy(Pawfgr)
2831 :
2832 : ! Deallocation for PAW.
2833 201 : if (Dtset%usepaw==1) then
2834 5 : call pawrhoij_free(KS_Pawrhoij)
2835 36 : ABI_FREE(KS_Pawrhoij)
2836 5 : call pawfgrtab_free(Pawfgrtab)
2837 5 : call paw_ij_free(KS_paw_ij)
2838 5 : call paw_an_free(KS_paw_an)
2839 5 : call pawpwff_free(Paw_pwff)
2840 5 : if (gwcalctyp>=10) then
2841 0 : call pawrhoij_free(QP_pawrhoij)
2842 0 : call paw_ij_free(QP_paw_ij)
2843 0 : call paw_an_free(QP_paw_an)
2844 : end if
2845 5 : if (Dtset%pawcross==1) then
2846 0 : call paw_pwaves_lmn_free(Paw_onsite)
2847 0 : call wfdf%free()
2848 : end if
2849 : end if
2850 201 : ABI_FREE(Paw_onsite)
2851 232 : ABI_FREE(Pawfgrtab)
2852 209 : ABI_FREE(Paw_pwff)
2853 232 : ABI_FREE(KS_paw_ij)
2854 232 : ABI_FREE(KS_paw_an)
2855 201 : if (gwcalctyp>=10) then
2856 76 : ABI_FREE(QP_pawrhoij)
2857 76 : ABI_FREE(QP_paw_an)
2858 76 : ABI_FREE(QP_paw_ij)
2859 : end if
2860 :
2861 201 : call wfd%free(); call destroy_mpi_enreg(MPI_enreg_seq)
2862 201 : call Kmesh%free(); call Qmesh%free(); call Gsph_Max%free(); call Gsph_x%free(); call Gsph_c%free()
2863 201 : call Vcp%free(); call cryst%free(); call Sr%free()
2864 201 : if (.not.rdm_update) call epsm1%free()
2865 201 : call PPm%free(); call Hdr_sigma%free(); call Hdr_wfk%free(); call ks_ebands%free();call qp_ebands%free(); call KS_me%free()
2866 201 : call littlegroup_free(Ltg_k)
2867 840 : ABI_FREE(Ltg_k)
2868 201 : call esymm_free(KS_sym)
2869 1433 : ABI_FREE(KS_sym)
2870 :
2871 201 : if (Sigp%symsigma == 1 .and. gwcalctyp>=20) then
2872 0 : call esymm_free(QP_sym)
2873 0 : ABI_FREE(QP_sym)
2874 : end if
2875 :
2876 201 : call Sigp%free()
2877 :
2878 201 : call timab(426,2,tsec) ! finalize
2879 201 : call timab(401,2,tsec)
2880 :
2881 : DBG_EXIT('COLL')
2882 :
2883 1005 : end subroutine sigma
2884 : !!***
2885 :
2886 : !!****f* ABINIT/setup_sigma
2887 : !! NAME
2888 : !! setup_sigma
2889 : !!
2890 : !! FUNCTION
2891 : !! Initialize the data type containing parameters for a sigma calculation.
2892 : !!
2893 : !! INPUTS
2894 : !! acell(3)=length scales of primitive translations (bohr)
2895 : !! wfk_fname=Name of the WFK file.
2896 : !! Dtset<type(dataset_type)>=all input variables for this dataset
2897 : !! Dtfil<type(datafiles_type)>=variables related to files
2898 : !! rprim(3,3)=dimensionless real space primitive translations
2899 : !! Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the WFK file
2900 : !!
2901 : !! OUTPUT
2902 : !! Sigp<sigparams_t>=Parameters governing the self-energy calculation.
2903 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
2904 : !! Qmesh <kmesh_t>=Structure describing the q-point sampling.
2905 : !! Cryst<crystal_t>=Info on unit cell and symmetries.
2906 : !! Gsph_Max<gsphere_t>=Info on the G-sphere
2907 : !! Gsph_c<gsphere_t>=Info on the G-sphere for W and Sigma_c
2908 : !! Gsph_x<gsphere_t>=Info on the G-sphere for and Sigma_x
2909 : !! Hdr_wfk<hdr_type>=The header of the WFK file
2910 : !! Hdr_out<hdr_type>=The header to be used for the results of sigma calculations.
2911 : !! Vcp<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique.
2912 : !! epsm1<epsm1_t>=Datatype storing data used to construct the screening (partially Initialized in OUTPUT)
2913 : !! ks_ebands<ebands_t>=The KS energies and occupation factors.
2914 : !! gwc_ngfft(18), gwx_ngfft(18)= FFT meshes for the oscillator strengths used for the correlated and the
2915 : !! exchange part of the self-energy, respectively.
2916 : !! comm=MPI communicator.
2917 : !!
2918 : !! SOURCE
2919 :
2920 21105 : subroutine setup_sigma(codvsn,wfk_fname,acell,rprim,Dtset,Dtfil,Psps,Pawtab,&
2921 : gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,ks_ebands,Gsph_Max,Gsph_x,Gsph_c,Vcp,epsm1,Sigp,comm)
2922 :
2923 : !Arguments ------------------------------------
2924 : !scalars
2925 : integer,intent(in) :: comm
2926 : character(len=8),intent(in) :: codvsn
2927 : character(len=*),intent(in) :: wfk_fname
2928 : type(Datafiles_type),intent(in) :: Dtfil
2929 : type(Dataset_type),intent(inout) :: Dtset
2930 : type(Pseudopotential_type),intent(in) :: Psps
2931 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw)
2932 : type(sigparams_t),intent(out) :: Sigp
2933 : type(epsm1_t),intent(out) :: epsm1
2934 : type(ebands_t),intent(out) :: ks_ebands
2935 : type(kmesh_t),intent(out) :: Kmesh,Qmesh
2936 : type(crystal_t),intent(out) :: Cryst
2937 : type(gsphere_t),intent(out) :: Gsph_Max,Gsph_x,Gsph_c
2938 : type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out
2939 : type(vcoul_t),intent(out) :: Vcp
2940 : !arrays
2941 : integer,intent(out) :: gwc_ngfft(18),gwx_ngfft(18)
2942 : real(dp),intent(in) :: acell(3),rprim(3,3)
2943 :
2944 : !Local variables-------------------------------
2945 : !scalars
2946 : integer,parameter :: pertcase0 = 0, master = 0
2947 : integer :: bantot,enforce_sym,gwcalctyp,ib,ibtot,icutcoul_eff,ii,ikcalc,ikibz,io,isppol,itypat,jj,method
2948 : integer :: mod10,mqmem,mband,ng_kss,nsheps,ikcalc2bz,ierr,gap_err,ng, nsppol
2949 : integer :: gwc_nfftot,gwx_nfftot,nqlwl,test_npwkss,my_rank,nprocs,ik,nk_found,ifo,timrev,usefock_ixc
2950 : integer :: iqbz,isym,iq_ibz,itim,ic,pinv,ig1,ng_sigx,spin,gw_qprange,ivcoul_init,nvcoul_init,xclevel_ixc
2951 : real(dp),parameter :: OMEGASIMIN=0.01d0
2952 : real(dp) :: domegas,domegasi,ucvol,rcut, drude_plasmon_freq, wmax
2953 : logical :: ltest,remove_inv,changed,found
2954 : character(len=500) :: msg, iw_mesh_type
2955 : character(len=fnlen) :: fname,fcore,string
2956 201 : type(wvl_internal_type) :: wvl
2957 201 : type(gaps_t) :: ks_gaps
2958 : !arrays
2959 : integer :: ng0sh_opt(3),G0(3),q_umklp(3),kpos(6), units(2)
2960 201 : integer,allocatable :: npwarr(:),val_indices(:,:)
2961 201 : integer,pointer :: gvec_kss(:,:),gsphere_sigx_p(:,:)
2962 201 : integer,pointer :: test_gvec_kss(:,:)
2963 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3),q_bz(3),gamma_point(3,1)
2964 201 : real(dp),pointer :: energies_p(:,:,:)
2965 201 : real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:)
2966 201 : type(Pawrhoij_type),allocatable :: Pawrhoij(:)
2967 4422 : type(vcoul_t) :: Vcp_ks
2968 : integer :: nbsum
2969 : real(dp) :: te_min = -one, te_max = one
2970 : real(dp) :: ft_max_error(3) = -one
2971 : real(dp) :: cosft_duality_error = -one
2972 201 : real(dp),allocatable :: tau_mesh(:), tau_wgs(:), iw_mesh(:), iw_wgs(:)
2973 201 : real(dp),allocatable :: cosft_wt(:,:), cosft_tw(:,:), sinft_wt(:,:)
2974 : ! *************************************************************************
2975 :
2976 : DBG_ENTER('COLL')
2977 603 : units = [std_out, ab_out]
2978 201 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
2979 :
2980 201 : nsppol = dtset%nsppol
2981 1433 : ABI_CHECK(ALL(Dtset%nband(1:Dtset%nkpt*nsppol) == Dtset%nband(1)), 'Dtset%nband(:) must be constant')
2982 :
2983 : ! Basic parameters
2984 201 : Sigp%ppmodel = Dtset%ppmodel
2985 201 : Sigp%gwcalctyp = Dtset%gwcalctyp
2986 201 : Sigp%nbnds = Dtset%nband(1)
2987 201 : Sigp%symsigma = Dtset%symsigma
2988 201 : Sigp%zcut = Dtset%zcut
2989 201 : Sigp%mbpt_sciss = Dtset%mbpt_sciss
2990 :
2991 201 : timrev= 2 ! This information is not reported in the header
2992 : ! 1 => do not use time-reversal symmetry
2993 : ! 2 => take advantage of time-reversal symmetry
2994 201 : if (any(dtset%kptopt == [3, 4])) timrev = 1
2995 :
2996 : ! For HF, SEX or COHSEX use Hybertsen-Louie PPM (only $\omega=0$) ===
2997 : ! Use fake screening for HF.
2998 : ! FIXME Why, we should not redefine Sigp%ppmodel
2999 201 : gwcalctyp = Sigp%gwcalctyp
3000 201 : mod10 = MOD(Sigp%gwcalctyp, 10)
3001 201 : if (any(mod10 == [5, 6, 7])) Sigp%ppmodel=2
3002 163 : if (mod10<5 .and. MOD(Sigp%gwcalctyp,1)/=1) then !
3003 : ! One shot GW (PPM or contour deformation).
3004 124 : if (Dtset%nomegasrd==1) then
3005 : ! avoid division by zero!
3006 0 : Sigp%nomegasrd =1
3007 0 : Sigp%maxomega4sd=zero
3008 0 : Sigp%deltae =zero
3009 : else
3010 124 : Sigp%nomegasrd = Dtset%nomegasrd
3011 124 : Sigp%maxomega4sd = Dtset%omegasrdmax
3012 124 : Sigp%deltae = (2*Sigp%maxomega4sd)/(Sigp%nomegasrd-1)
3013 : endif
3014 : else
3015 : ! For AC no need to evaluate derivative by finite differences.
3016 77 : Sigp%nomegasrd =1
3017 77 : Sigp%maxomega4sd=zero
3018 77 : Sigp%deltae =zero
3019 : end if
3020 :
3021 : ! Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume
3022 201 : call mkrdim(acell,rprim,rprimd)
3023 201 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
3024 :
3025 201 : Sigp%npwwfn = Dtset%npwwfn
3026 201 : Sigp%npwx = Dtset%npwsigx
3027 :
3028 : ! Read parameters of the WFK, verifify them and retrieve all G-vectors.
3029 201 : if (dtset%userie == 456) then
3030 0 : call wfk_read_eigenvalues("SC_WFK",energies_p,Hdr_wfk,comm)
3031 : else
3032 201 : call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm)
3033 : end if
3034 1433 : mband = MAXVAL(Hdr_wfk%nband)
3035 :
3036 201 : remove_inv = .FALSE.
3037 201 : if (dtset%userie /= 456) then
3038 201 : call hdr_wfk%vs_dtset(dtset)
3039 : end if
3040 :
3041 201 : test_npwkss = 0
3042 : call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,&
3043 201 : gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr)
3044 201 : ABI_CHECK(ierr==0, "Fatal error in make_gvec_kss")
3045 :
3046 603 : ABI_MALLOC(gvec_kss,(3, test_npwkss))
3047 1093537 : gvec_kss = test_gvec_kss
3048 201 : ng_kss = test_npwkss
3049 :
3050 201 : ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2))
3051 201 : ierr = 0
3052 136868 : do ig1=1,ng
3053 546869 : if (ANY(gvec_kss(:,ig1)/=test_gvec_kss(:,ig1))) then
3054 0 : ierr=ierr+1
3055 0 : write(std_out,*)" gvec_kss ",ig1,"/",ng,gvec_kss(:,ig1),test_gvec_kss(:,ig1)
3056 : end if
3057 : end do
3058 201 : ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss")
3059 201 : ABI_FREE(test_gvec_kss)
3060 :
3061 : ! Get important dimensions from the WFK header
3062 201 : Sigp%nsppol = Hdr_wfk%nsppol
3063 201 : Sigp%nspinor = Hdr_wfk%nspinor
3064 201 : Sigp%nsig_ab = Hdr_wfk%nspinor**2 ! TODO Is it useful calculating only diagonal terms?
3065 :
3066 201 : if (Sigp%nbnds > mband) then
3067 : write(msg,'(2a,2(a,i0))') &
3068 0 : 'Number of bands stored WFK file is less than required. ',ch10,&
3069 0 : "WFK mband: ", mband, ", self-energy nband: ", sigp%nbnds
3070 0 : ABI_ERROR(msg)
3071 : end if
3072 :
3073 : ! Check input
3074 201 : if (any(Sigp%ppmodel== [3, 4])) then
3075 6 : if (gwcalctyp >= 10) then
3076 0 : ABI_ERROR(sjoin('The ppmodel chosen and gwcalctyp: ', itoa(Dtset%gwcalctyp),' are not compatible.'))
3077 : end if
3078 6 : if (Sigp%nspinor==2) then
3079 0 : ABI_ERROR(sjoin('The ppmodel chosen and nspinor: ', itoa(Sigp%nspinor), ' are not compatible.'))
3080 : end if
3081 : end if
3082 :
3083 : ! Create crystal_t data type
3084 201 : cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv)
3085 201 : call cryst%print()
3086 :
3087 201 : if (Sigp%npwwfn > ng_kss) then ! cannot use more G"s for the wfs than those stored on file
3088 2 : Sigp%npwwfn =ng_kss
3089 2 : Dtset%npwwfn =ng_kss
3090 : write(msg,'(2a,(a,i0,a))')&
3091 2 : 'Number of G-vectors for WFS found in the KSS file is less than required',ch10,&
3092 4 : 'calculation will proceed with npwwfn = ',Sigp%npwwfn,ch10
3093 2 : ABI_WARNING(msg)
3094 : end if
3095 :
3096 201 : if (Sigp%npwx>ng_kss) then
3097 : ! Have to recalculate the (large) sphere for Sigma_x.
3098 0 : pinv=1; if (remove_inv.and.Cryst%timrev==2) pinv=-1
3099 0 : gamma_point(:,1) = (/zero,zero,zero/); nullify(gsphere_sigx_p)
3100 :
3101 : call merge_and_sort_kg(1,gamma_point,Dtset%ecutsigx,Cryst%nsym,pinv,Cryst%symrel,&
3102 0 : Cryst%gprimd,gsphere_sigx_p,Dtset%prtvol)
3103 :
3104 0 : ng_sigx=SIZE(gsphere_sigx_p,DIM=2)
3105 0 : Sigp%npwx = ng_sigx
3106 0 : Dtset%npwsigx = ng_sigx
3107 :
3108 : write(msg,'(2a,(a,i0,a))')&
3109 0 : 'Number of G-vectors for Sigma_x found in the KSS file is less than required',ch10,&
3110 0 : 'calculation will proceed with npwsigx = ',Sigp%npwx,ch10
3111 0 : ABI_WARNING(msg)
3112 :
3113 0 : ltest = (Sigp%npwx >= ng_kss)
3114 0 : ABI_CHECK(ltest,"Sigp%npwx<ng_kss!")
3115 :
3116 : ! Fill gvec_kss with larger sphere.
3117 0 : ABI_FREE(gvec_kss)
3118 0 : ABI_MALLOC(gvec_kss,(3,Sigp%npwx))
3119 0 : gvec_kss = gsphere_sigx_p
3120 0 : ABI_FREE(gsphere_sigx_p)
3121 : end if
3122 :
3123 : ! Set up of the k-points and tables in the whole BZ
3124 : ! TODO Recheck symmorphy and inversion
3125 201 : call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.)
3126 : !call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.)
3127 :
3128 : ! Some required information are not filled up inside kmesh_init
3129 : ! So doing it here, even though it is not clean
3130 2613 : Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:)
3131 201 : Kmesh%nshift =Dtset%nshiftk
3132 603 : ABI_MALLOC(Kmesh%shift,(3,Kmesh%nshift))
3133 1005 : Kmesh%shift(:,:) =Dtset%shiftk(:,1:Dtset%nshiftk)
3134 :
3135 201 : call Kmesh%print(units, header="K-mesh for the wavefunctions", prtvol=Dtset%prtvol)
3136 :
3137 : ! Initialize the band structure datatype
3138 : ! Copy WFK energies and occupations up to Sigp%nbnds==Dtset%nband(:)
3139 1433 : bantot = SUM(Dtset%nband(1:Dtset%nkpt*nsppol))
3140 603 : ABI_MALLOC(doccde,(bantot))
3141 402 : ABI_MALLOC(eigen,(bantot))
3142 402 : ABI_MALLOC(occfact,(bantot))
3143 87711 : doccde(:)=zero; eigen(:)=zero; occfact(:)=zero
3144 :
3145 201 : jj=0; ibtot=0
3146 406 : do isppol=1,nsppol
3147 1638 : do ikibz=1,Dtset%nkpt
3148 46196 : do ib=1,Hdr_wfk%nband(ikibz+(isppol-1)*Dtset%nkpt)
3149 44759 : ibtot=ibtot+1
3150 45991 : if (ib<=Sigp%nbnds) then
3151 29170 : jj=jj+1
3152 29170 : occfact(jj)=Hdr_wfk%occ(ibtot)
3153 29170 : eigen (jj)=energies_p(ib,ikibz,isppol)
3154 : end if
3155 : end do
3156 : end do
3157 : end do
3158 201 : ABI_FREE(energies_p)
3159 :
3160 : ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of
3161 : ! symmetry operations in the old GW code (symmorphy and inversion)
3162 1423 : ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6))
3163 201 : ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt')
3164 :
3165 603 : ABI_MALLOC(npwarr,(Dtset%nkpt))
3166 1423 : npwarr(:)=Sigp%npwwfn
3167 :
3168 : call ks_ebands%init(bantot, Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
3169 : doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,&
3170 : Kmesh%nibz,npwarr,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,&
3171 : dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig,&
3172 201 : dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
3173 :
3174 201 : ABI_FREE(doccde)
3175 201 : ABI_FREE(eigen)
3176 201 : ABI_FREE(npwarr)
3177 :
3178 : ! Calculate KS occupation numbers and ks_vbk(nkibz, nsppol)
3179 : ! ks_vbk gives the (valence|last Fermi band) index for each k and spin.
3180 : ! spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals
3181 201 : call ks_ebands%update_occ(Dtset%spinmagntarget, prtvol=0)
3182 :
3183 201 : ks_gaps = ks_ebands%get_gaps(gap_err)
3184 402 : call ks_gaps%print([std_out])
3185 201 : call ks_ebands%report_gap(unit=std_out)
3186 :
3187 804 : ABI_MALLOC(val_indices,(ks_ebands%nkpt, nsppol))
3188 201 : val_indices = ks_ebands%get_valence_idx()
3189 :
3190 : ! Create Sigma header
3191 : ! TODO Fix problems with symmorphy and k-points
3192 201 : call Hdr_out%init(ks_ebands,codvsn,Dtset,Pawtab,pertcase0,Psps,wvl)
3193 :
3194 : ! Get Pawrhoij from the header of the WFK file
3195 634 : ABI_MALLOC(Pawrhoij, (Cryst%natom*Dtset%usepaw))
3196 201 : if (Dtset%usepaw==1) then
3197 5 : call pawrhoij_alloc(Pawrhoij, 1, Dtset%nspden, Dtset%nspinor, nsppol, Cryst%typat, pawtab=Pawtab)
3198 5 : call pawrhoij_copy(Hdr_wfk%Pawrhoij, Pawrhoij)
3199 : end if
3200 :
3201 201 : call hdr_out%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1))
3202 201 : ABI_FREE(occfact)
3203 201 : call pawrhoij_free(Pawrhoij)
3204 232 : ABI_FREE(Pawrhoij)
3205 :
3206 : ! ===========================================================
3207 : ! ==== Setup of k-points and bands for the GW corrections ====
3208 : ! ===========================================================
3209 : ! * maxbdgw and minbdgw are the Max and min band index for GW corrections over k-points.
3210 : ! They are used to dimension the wavefunctions and to calculate the matrix elements.
3211 : !
3212 201 : if (dtset%nkptgw == 0) then
3213 : !
3214 : ! Use qp_range to select the interesting k-points and the corresponding bands.
3215 : !
3216 : ! 0 --> Compute the QP corrections only for the fundamental and the direct gap.
3217 : ! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num`
3218 : ! bands above and below the Fermi level.
3219 : ! -num --> Compute the QP corrections for all the k-points in the irreducible zone.
3220 : ! Include all occupied states and `num` empty states.
3221 :
3222 30 : call wrtout(std_out, "nkptgw == 0 ==> Automatic selection of k-points and bands for the corrections.")
3223 :
3224 30 : if (gap_err /= 0 .and. dtset%gw_qprange == 0) then
3225 1 : msg = "Problem while computing fundamental/direct gap (likely metal). Will replace gw_qprange=0 with gw_qprange=1"
3226 1 : ABI_WARNING(msg)
3227 1 : dtset%gw_qprange = 1
3228 : end if
3229 :
3230 30 : gw_qprange = dtset%gw_qprange
3231 :
3232 30 : if (dtset%ucrpa > 0) then
3233 0 : dtset%nkptgw = Kmesh%nbz
3234 0 : Sigp%nkptgw = dtset%nkptgw
3235 0 : ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw))
3236 0 : ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3237 0 : ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3238 0 : Sigp%kptgw(:,:) = Kmesh%bz(:,:)
3239 0 : Sigp%minbnd = 1
3240 0 : Sigp%maxbnd = Sigp%nbnds
3241 :
3242 30 : else if (gw_qprange /= 0) then
3243 : ! Include all the k-points in the IBZ.
3244 : ! Note that kptgw == ebands%kptns so we can use a single ik index in the loop over k-points
3245 : ! No need to map kptgw onto ebands%kptns.
3246 25 : dtset%nkptgw = Kmesh%nibz
3247 25 : Sigp%nkptgw = dtset%nkptgw
3248 75 : ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw))
3249 100 : ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3250 75 : ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3251 649 : Sigp%kptgw(:,:) = Kmesh%ibz(:,:)
3252 206 : Sigp%minbnd = 1
3253 206 : Sigp%maxbnd = Sigp%nbnds
3254 :
3255 25 : if (gw_qprange > 0) then
3256 : ! All k-points: Add buffer of bands above and below the Fermi level.
3257 2 : do spin=1,nsppol
3258 8 : do ik=1,Sigp%nkptgw
3259 6 : Sigp%minbnd(ik, spin) = MAX(val_indices(ik,spin) - gw_qprange, 1)
3260 7 : Sigp%maxbnd(ik, spin) = MIN(val_indices(ik,spin) + gw_qprange + 1, Sigp%nbnds)
3261 : end do
3262 : end do
3263 :
3264 : else
3265 : ! All k-points: include all occupied states and -gw_qprange empty states.
3266 198 : Sigp%minbnd = 1
3267 48 : do spin=1,nsppol
3268 198 : do ik=1,Sigp%nkptgw
3269 174 : Sigp%maxbnd(ik, spin) = MIN(val_indices(ik,spin) - gw_qprange, Sigp%nbnds)
3270 : end do
3271 : end do
3272 : end if
3273 :
3274 : else
3275 : ! gw_qprange is not specified in the input.
3276 : ! Include the direct and the fundamental KS gap.
3277 : ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore
3278 : ! we have compute the union of the k-points where the fundamental and the direct KS gaps are located.
3279 : !
3280 : ! Find the list of `interesting` kpoints.
3281 5 : ABI_CHECK(gap_err == 0, "gw_qprange 0 cannot be used because I cannot find the gap (gap_err !=0)")
3282 5 : nk_found = 1; kpos(1) = ks_gaps%fo_kpos(1,1)
3283 :
3284 10 : do spin=1,nsppol
3285 25 : do ifo=1,3
3286 15 : ik = ks_gaps%fo_kpos(ifo, spin)
3287 15 : found = .FALSE.; jj = 0
3288 30 : do while (.not. found .and. jj < nk_found)
3289 15 : jj = jj + 1; found = (kpos(jj) == ik)
3290 : end do
3291 20 : if (.not. found) then
3292 5 : nk_found = nk_found + 1; kpos(nk_found) = ik
3293 : end if
3294 : end do
3295 : end do
3296 :
3297 : ! Now we can define the list of k-points and the bands range.
3298 5 : dtset%nkptgw = nk_found
3299 5 : Sigp%nkptgw = dtset%nkptgw
3300 :
3301 15 : ABI_MALLOC(Sigp%kptgw,(3, Sigp%nkptgw))
3302 20 : ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3303 15 : ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3304 :
3305 15 : do ii=1,Sigp%nkptgw
3306 10 : ik = kpos(ii)
3307 40 : Sigp%kptgw(:,ii) = Kmesh%ibz(:,ik)
3308 25 : do spin=1,nsppol
3309 10 : Sigp%minbnd(ii,spin) = val_indices(ik, spin)
3310 20 : Sigp%maxbnd(ii,spin) = val_indices(ik, spin) + 1
3311 : end do
3312 : end do
3313 : end if
3314 :
3315 : else
3316 : ! Treat only the k-points and bands specified in the input file.
3317 171 : Sigp%nkptgw = dtset%nkptgw
3318 513 : ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw))
3319 684 : ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol))
3320 513 : ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol))
3321 :
3322 346 : do spin=1,nsppol
3323 656 : Sigp%minbnd(:,spin)= dtset%bdgw(1,:,spin)
3324 827 : Sigp%maxbnd(:,spin)= dtset%bdgw(2,:,spin)
3325 : end do
3326 :
3327 684 : do ii=1,3
3328 2103 : do ikcalc=1,Sigp%nkptgw
3329 1932 : Sigp%kptgw(ii,ikcalc) = Dtset%kptgw(ii,ikcalc)
3330 : end do
3331 : end do
3332 :
3333 346 : do spin=1,nsppol
3334 827 : do ikcalc=1,Sigp%nkptgw
3335 656 : if (Dtset%bdgw(2,ikcalc,spin) > Sigp%nbnds) then
3336 : write(msg,'(a,2i0,2(a,i0),2a,i0)')&
3337 0 : "For (k,s) ",ikcalc, spin," bdgw= ",Dtset%bdgw(2,ikcalc,spin), " > nbnds=",Sigp%nbnds,ch10,&
3338 0 : "Calculation will continue with bdgw =",Sigp%nbnds
3339 0 : ABI_COMMENT(msg)
3340 0 : Dtset%bdgw(2, ikcalc, spin) = Sigp%nbnds
3341 : end if
3342 : end do
3343 : end do
3344 :
3345 : end if
3346 :
3347 : ! Make sure that all the degenerate states are included.
3348 : ! * We will have to average the GW corrections over degenerate states if symsigma=1 is used.
3349 : ! * KS states belonging to the same irreducible representation should be included in the basis set used for SCGW.
3350 201 : if (Sigp%symsigma /=0 .or. gwcalctyp >= 10) then
3351 299 : do isppol=1,nsppol
3352 869 : do ikcalc=1,Sigp%nkptgw
3353 :
3354 668 : if (kmesh%has_IBZ_item(Sigp%kptgw(:,ikcalc), ikibz, G0)) then
3355 : call ks_ebands%enclose_degbands(ikibz,isppol, &
3356 517 : Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol),changed,dtset%symsigma_de)
3357 :
3358 517 : if (changed) then
3359 : write(msg,'(2(a,i0),2a,2(1x,i0))')&
3360 48 : "Not all the degenerate states at ikcalc= ",ikcalc,", spin= ",isppol, ch10, &
3361 96 : "were included in the bdgw set. bdgw has been changed to: ",Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol)
3362 48 : ABI_COMMENT(msg)
3363 : end if
3364 : else
3365 : write(msg,'(3a)')&
3366 0 : ' not in the list of k points treated in the preparatory SCF run.',ch10, &
3367 0 : ' Change kptgw, or shiftk of previous run.'
3368 0 : ABI_ERROR(sjoin('k-point', ktoa(Sigp%kptgw(:,ikcalc)),trim(msg)))
3369 : end if
3370 :
3371 : end do
3372 : end do
3373 : end if
3374 :
3375 1053 : Sigp%minbdgw = MINVAL(Sigp%minbnd)
3376 1053 : Sigp%maxbdgw = MAXVAL(Sigp%maxbnd)
3377 :
3378 : ! Check if there are duplicated k-point in Sigp%
3379 840 : do ii=1,Sigp%nkptgw
3380 2282 : do jj=ii+1,Sigp%nkptgw
3381 2081 : if (isamek(Sigp%kptgw(:,ii), Sigp%kptgw(:,jj), G0)) then
3382 : write(msg,'(5a)')&
3383 0 : 'kptgw contains duplicated k-points. This is not allowed since ',ch10,&
3384 0 : 'the QP corrections for this k-point will be calculated more than once. ',ch10,&
3385 0 : 'Check your input file. '
3386 0 : ABI_ERROR(msg)
3387 : end if
3388 : end do
3389 : end do
3390 :
3391 : !=== Check if the k-points are in the BZ ===
3392 : ! FB: Honestly the code is not able to treat k-points, which are not in the input list of points in the IBZ.
3393 : ! This extension should require to change the code in different places.
3394 : ! Therefore, one should by now prevent the user from calculating sigma for a k-point not in the IBZ.
3395 603 : ABI_MALLOC(Sigp%kptgw2bz, (Sigp%nkptgw))
3396 :
3397 840 : do ikcalc=1,Sigp%nkptgw
3398 840 : if (kmesh%has_BZ_item(Sigp%kptgw(:,ikcalc), ikcalc2bz, G0)) then
3399 639 : Sigp%kptgw2bz(ikcalc) = ikcalc2bz
3400 : else
3401 : write(msg,'(3a)')&
3402 0 : ' not in the list of k points treated in the preparatory SCF run.',ch10,' Change kptgw, or shiftk of previous run.'
3403 0 : ABI_ERROR(sjoin('k-point:', ktoa(Sigp%kptgw(:,ikcalc)),trim(msg)))
3404 : end if
3405 : end do
3406 :
3407 : ! Warn the user if SCGW run and not all the k-points are included.
3408 201 : if (gwcalctyp >= 10 .and. Sigp%nkptgw /= Hdr_wfk%nkpt) then
3409 3 : write(msg,'(3a,2(a,i0),2a)')ch10,&
3410 3 : " COMMENT: In a self-consistent GW run, the QP corrections should be calculated for all the k-points of the KSS file ",ch10,&
3411 3 : " but nkptgw= ",Sigp%nkptgw," and WFK nkpt= ",Hdr_wfk%nkpt,ch10,&
3412 6 : " Assuming expert user. Execution will continue. "
3413 3 : call wrtout(ab_out, msg)
3414 : end if
3415 :
3416 : ! Setup of the table used in the case of SCGW on wavefunctions to reduce the number
3417 : ! of elements <i,kgw,s|\Sigma|j,kgw,s> that have to be calculated. No use of symmetries, except for Hermiticity.
3418 201 : call sigma_tables(Sigp, Kmesh)
3419 :
3420 : ! === Read external file and initialize basic dimension of epsm1 ===
3421 : ! TODO use mqmem as input variable instead of gwmem
3422 :
3423 : ! === If required, use a matrix for $\Sigma_c$ which is smaller than that stored on file ===
3424 : ! * By default the entire matrix is read and used,
3425 : ! * Define consistently npweps and ecuteps for \Sigma_c according the input
3426 201 : if (Dtset%npweps>0.or.Dtset%ecuteps>0) then
3427 128 : if (Dtset%npweps>0) Dtset%ecuteps=zero
3428 128 : nsheps=0
3429 128 : call setshells(Dtset%ecuteps,Dtset%npweps,nsheps,Dtset%nsym,gmet,gprimd,Dtset%symrel,'eps',ucvol)
3430 : end if
3431 :
3432 201 : mqmem=0; if (Dtset%gwmem/10==1) mqmem=1
3433 :
3434 201 : if (dtset%getscr /=0 .or. dtset%irdscr/=0 .or. dtset%getscr_filepath /= ABI_NOFILE) then
3435 161 : fname = Dtfil%fnameabi_scr
3436 40 : else if (Dtset%getsuscep/=0 .or. Dtset%irdsuscep/=0) then
3437 5 : fname = Dtfil%fnameabi_sus
3438 : else
3439 35 : fname = Dtfil%fnameabi_scr
3440 : !FIXME this has to be cleaned, in tgw2_3 Dtset%get* and Dtset%ird* are not defined
3441 : !ABI_ERROR("getsuscep or irdsuscep are not defined")
3442 : end if
3443 : !
3444 : ! === Setup of q-mesh in the whole BZ ===
3445 : ! Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed,
3446 : ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed)
3447 : !
3448 201 : if (sigp%needs_w()) then
3449 166 : if (.not. file_exists(fname)) then
3450 166 : fname = nctk_ncify(fname)
3451 166 : ABI_COMMENT(sjoin("File not found. Will try netcdf file:", fname))
3452 : end if
3453 :
3454 : ! Initialize epsm1 from fname.
3455 166 : call epsm1%from_file(fname, mqmem, Dtset%npweps, comm)
3456 :
3457 166 : Sigp%npwc=epsm1%npwe
3458 166 : if (Sigp%npwc>Sigp%npwx) then
3459 14 : Sigp%npwc=Sigp%npwx
3460 14 : ABI_COMMENT("Found npw_correlation > npw_exchange, Imposing npwc=npwx")
3461 : ! There is a good reason for doing so, see csigme.F90 and the size of the arrays
3462 : ! rhotwgp and rhotwgp: we need to define a max size and we opt for Sigp%npwx.
3463 : end if
3464 :
3465 166 : if (Dtset%nfreqim_conv==0) then
3466 : ! If no extra frequencies is requested, keep the original grids.
3467 165 : epsm1%nomega_i_conv = 0
3468 1 : else if (Dtset%nfreqim_conv < 0) then
3469 : ! If negative number of frequencies is requested, multiply the number of frequencies in the file by the absolute value.
3470 0 : epsm1%nomega_i_conv = abs(Dtset%nfreqim_conv) * epsm1%nomega_i
3471 0 : Dtset%nfreqim_conv = epsm1%nomega_i_conv
3472 1 : else if (Dtset%nfreqim_conv >= epsm1%nomega_i) then
3473 : ! If the requested number of frequencies is larger than the number in the file, use input value.
3474 1 : epsm1%nomega_i_conv = Dtset%nfreqim_conv
3475 : else if (Dtset%nfreqim_conv < epsm1%nomega_i) then
3476 : ! If the requested number of frequencies is less than the number in the file, give an error
3477 0 : ABI_ERROR(sjoin("Requested number of frequencies for convolution is non-zero and less than nfreqim in the file: ", itoa(Dtset%nfreqim_conv)," < ", itoa(epsm1%nomega_i)))
3478 : end if
3479 :
3480 166 : epsm1%npwe=Sigp%npwc
3481 166 : Dtset%npweps=epsm1%npwe
3482 166 : call Qmesh%init(Cryst,epsm1%nqibz,epsm1%qibz,Dtset%kptopt)
3483 :
3484 : else
3485 35 : epsm1%npwe =1
3486 35 : Sigp%npwc =1
3487 35 : Dtset%npweps=1
3488 35 : call qmesh%find_qmesh(Cryst,Kmesh)
3489 35 : ABI_MALLOC(epsm1%gvec,(3,1))
3490 140 : epsm1%gvec(:,1) = [0, 0, 0]
3491 35 : epsm1%nomega = 0
3492 : end if
3493 :
3494 201 : call Qmesh%print(units, header="Q-mesh for screening function", prtvol=Dtset%prtvol)
3495 :
3496 10315 : do iqbz=1,Qmesh%nbz
3497 10114 : call qmesh%get_BZ_item(iqbz, q_bz, iq_ibz, isym, itim, umklp=q_umklp)
3498 : !print *, "iqbz, q_bz, iq_ibz, isym, itim, umklp"
3499 : !print *, iqbz, q_bz, iq_ibz, isym, itim, q_umklp
3500 :
3501 40657 : if (ANY(q_umklp/=0)) then
3502 0 : sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz))
3503 0 : write(std_out,*) sq,Qmesh%bz(:,iqbz)
3504 : write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')&
3505 0 : 'qpoint ',Qmesh%bz(:,iqbz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,&
3506 0 : 'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,&
3507 0 : 'however a non-zero umklapp G_o vector is required and this is not yet allowed'
3508 0 : ABI_ERROR(msg)
3509 : end if
3510 : end do
3511 : !stop
3512 : !
3513 : ! === Find optimal value for G-sphere enlargement due to oscillator matrix elements ===
3514 : ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true size bz is buggy
3515 : ! * -one is used because we loop over all the possible differences, unlike screening
3516 :
3517 201 : call get_ng0sh(Sigp%nkptgw,Sigp%kptgw,Kmesh%nbz,Kmesh%bz,Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt)
3518 201 : call wrtout(std_out, sjoin(' Optimal value for ng0sh ', ltoa(ng0sh_opt)))
3519 804 : Sigp%mG0 = ng0sh_opt
3520 :
3521 : ! G-sphere for W and Sigma_c is initialized from the SCR file.
3522 201 : call Gsph_c%init(Cryst, epsm1%npwe, gvec=epsm1%gvec)
3523 201 : call Gsph_x%init(Cryst, Sigp%npwx, gvec=gvec_kss)
3524 201 : Sigp%ecuteps = Gsph_c%ecut
3525 201 : Dtset%ecuteps = Sigp%ecuteps
3526 :
3527 : ! Make biggest G-sphere of Sigp%npwvec vectors.
3528 201 : Sigp%npwvec=MAX(Sigp%npwwfn,Sigp%npwx)
3529 201 : call Gsph_Max%init(Cryst, Sigp%npwvec, gvec=gvec_kss)
3530 :
3531 : !BEGIN DEBUG
3532 : !Make sure that the two G-spheres are equivalent.
3533 : !ierr=0
3534 : !if (sigp%needs_w()) then
3535 : ! ng = MIN(SIZE(Gsph_c%gvec,DIM=2),SIZE(gvec_kss,DIM=2))
3536 : ! do ig1=1,ng
3537 : ! if (ANY(Gsph_c%gvec(:,ig1)/=gvec_kss(:,ig1))) then
3538 : ! ierr=ierr+1
3539 : ! write(std_out,*)" Gsph_c, gvec_kss ",ig1,"/",ng,Gsph_c%gvec(:,ig1),gvec_kss(:,ig1)
3540 : ! end if
3541 : ! end do
3542 : ! ABI_CHECK(ierr==0,"Mismatch between Gsph_c and gvec_kss")
3543 : !end if
3544 : !ierr=0
3545 : !ng = MIN(SIZE(Gsph_x%gvec,DIM=2),SIZE(gvec_kss,DIM=2))
3546 : !do ig1=1,ng
3547 : ! if (ANY(Gsph_x%gvec(:,ig1)/=gvec_kss(:,ig1))) then
3548 : ! ierr=ierr+1
3549 : ! write(std_out,*)" Gsph_x, gvec_kss ",ig1,"/",ng,Gsph_x%gvec(:,ig1),gvec_kss(:,ig1)
3550 : ! end if
3551 : !end do
3552 : !ABI_CHECK(ierr==0,"Mismatch between Gsph_x and gvec_kss")
3553 : !END DEBUG
3554 :
3555 201 : ABI_FREE(gvec_kss)
3556 : !
3557 : ! === Get Fourier components of the Coulomb term for all q-points in the IBZ ===
3558 : ! * If required, use a cutoff in the interaction
3559 : ! * Pcv%vc_sqrt contains Vc^{-1/2}
3560 : ! * Setup also the analytical calculation of the q->0 component
3561 : ! FIXME recheck ngfftf since I got different charge outside the cutoff region
3562 :
3563 201 : if (Dtset%gw_nqlwl==0) then
3564 201 : nqlwl=1
3565 201 : ABI_MALLOC(qlwl,(3,nqlwl))
3566 804 : qlwl(:,1)= GW_Q0_DEFAULT
3567 : else
3568 0 : nqlwl=Dtset%gw_nqlwl
3569 0 : ABI_MALLOC(qlwl,(3,nqlwl))
3570 0 : qlwl(:,:)=Dtset%gw_qlwl(:,1:nqlwl)
3571 : end if
3572 :
3573 : ! The Coulomb interaction used here might have two terms:
3574 : ! the first term generates the usual sigma self-energy, but possibly, one should subtract
3575 : ! from it the Coulomb interaction already present in the Kohn-Sham basis,
3576 : ! if the usefock associated to ixc is one.
3577 : ! The latter excludes (in the present implementation) mod(Dtset%gwcalctyp,10)==5
3578 201 : nvcoul_init=1
3579 201 : call get_xclevel(Dtset%ixc,xclevel_ixc,usefock_ixc)
3580 :
3581 201 : if(usefock_ixc==1)then
3582 3 : nvcoul_init=2
3583 3 : if (mod(Dtset%gwcalctyp, 10) == 5)then
3584 0 : write(msg,'(4a,i3,a,i3,4a,i5)')ch10,&
3585 0 : ' The starting wavefunctions were obtained from SCF calculations in the planewave basis set',ch10,&
3586 0 : ' with ixc = ',Dtset%ixc,' associated with usefock =',usefock_ixc,ch10,&
3587 0 : ' In this case, the present implementation does not allow that the self-energy for sigma corresponds to',ch10,&
3588 0 : ' mod(gwcalctyp,10)==5, while your gwcalctyp= ',Dtset%gwcalctyp
3589 0 : ABI_ERROR(msg)
3590 : endif
3591 : endif
3592 405 : do ivcoul_init=1,nvcoul_init
3593 204 : rcut = Dtset%gw_rcut
3594 204 : icutcoul_eff=Dtset%gw_icutcoul
3595 204 : Sigp%sigma_mixing=one
3596 204 : if( mod(Dtset%gwcalctyp,10)==5 .or. ivcoul_init==2)then
3597 38 : if(abs(Dtset%hyb_mixing)>tol8)then
3598 : ! Warning: the absolute value is needed, because of the singular way used to define the default for this input variable
3599 20 : Sigp%sigma_mixing=abs(Dtset%hyb_mixing)
3600 18 : else if(abs(Dtset%hyb_mixing_sr)>tol8)then
3601 18 : Sigp%sigma_mixing=abs(Dtset%hyb_mixing_sr)
3602 18 : icutcoul_eff=5
3603 : endif
3604 38 : if(abs(rcut)<tol6 .and. abs(Dtset%hyb_range_fock)>tol8) rcut=one/Dtset%hyb_range_fock
3605 : endif
3606 :
3607 405 : if (ivcoul_init == 1) then
3608 :
3609 201 : if (Gsph_x%ng > Gsph_c%ng) then
3610 : call Vcp%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3611 134 : Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
3612 : else
3613 : call Vcp%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3614 67 : Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
3615 : end if
3616 :
3617 : else
3618 :
3619 : ! Use a temporary Vcp_ks to compute the Coulomb interaction already present
3620 : ! in the Fock part of the Kohn-Sham Hamiltonian
3621 3 : if (Gsph_x%ng > Gsph_c%ng) then
3622 : call Vcp_ks%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3623 3 : Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
3624 : else
3625 : call Vcp_ks%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,&
3626 0 : Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
3627 : end if
3628 :
3629 : ! Now compute the residual Coulomb interaction.
3630 1194 : Vcp%vc_sqrt_resid=sqrt(Vcp%vc_sqrt**2-Sigp%sigma_mixing*Vcp_ks%vc_sqrt**2)
3631 3 : Vcp%i_sz_resid=Vcp%i_sz-Sigp%sigma_mixing*Vcp_ks%i_sz
3632 :
3633 : ! The mixing factor has already been accounted for, so set it back to one
3634 3 : Sigp%sigma_mixing=one
3635 3 : call Vcp_ks%free()
3636 :
3637 : !write(std_out,'(a)')' setup_sigma : the residual Coulomb interaction has been computed'
3638 : endif
3639 :
3640 : !#else
3641 : ! call Vcp%init(Gsph_Max,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,ivcoul_init,Dtset%vcutgeo,&
3642 : !& Dtset%ecutsigx,Sigp%npwx,nqlwl,qlwl,ngfftf,comm)
3643 : !#endif
3644 :
3645 : end do
3646 :
3647 : #if 1
3648 : ! Using the random q for the optical limit is one of the reasons
3649 : ! why sigma breaks the initial energy degeneracies.
3650 201 : if (dtset%userra > 100) then
3651 0 : call wrtout(units, "I am setting Vcp%i_sz=zero, Vcp%vc_sqrt(1,1)=czero, Vcp%vcqlwl_sqrt(1,1)=czero")
3652 0 : Vcp%i_sz=zero
3653 0 : Vcp%vc_sqrt(1,1)=czero
3654 0 : Vcp%vcqlwl_sqrt(1,1)=czero
3655 : end if
3656 : #endif
3657 :
3658 201 : ABI_FREE(qlwl)
3659 :
3660 201 : Sigp%ecuteps = Dtset%ecuteps
3661 201 : Sigp%ecutwfn = Dtset%ecutwfn
3662 201 : Sigp%ecutsigx = Dtset%ecutsigx
3663 :
3664 : ! === Setup of the FFT mesh for the oscillator strengths ===
3665 : ! * Init gwc_ngfft(7:18) and gwx_ngfft(7:18) with Dtset%ngfft(7:18)
3666 : ! * Here we redefine gwc_ngfft(1:6) according to the following options:
3667 : !
3668 : ! method == 0 --> FFT grid read from fft.in (debugging purpose)
3669 : ! method == 1 --> Normal FFT mesh
3670 : ! method == 2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90)
3671 : ! method == 3 --> Doubled FFT grid, same as the the FFT for the density,
3672 : !
3673 : ! enforce_sym == 1 --> Enforce a FFT mesh compatible with all the symmetry operation and FFT library
3674 : ! enforce_sym == 0 --> Find the smallest FFT grid compatible with the library, do not care about symmetries
3675 : !
3676 3819 : gwc_ngfft(1:18) = Dtset%ngfft(1:18)
3677 3819 : gwx_ngfft(1:18) = Dtset%ngfft(1:18)
3678 :
3679 201 : method = 2
3680 201 : if (Dtset%fftgw == 00 .or. Dtset%fftgw == 01) method = 0
3681 201 : if (Dtset%fftgw == 10 .or. Dtset%fftgw == 11) method = 1
3682 201 : if (Dtset%fftgw == 20 .or. Dtset%fftgw == 21) method = 2
3683 201 : if (Dtset%fftgw == 30 .or. Dtset%fftgw == 31) method = 3
3684 201 : enforce_sym = mod(dtset%fftgw, 10)
3685 :
3686 : ! FFT mesh for sigma_x.
3687 : call setmesh(gmet, Gsph_Max%gvec, gwx_ngfft, Sigp%npwvec, Sigp%npwx, Sigp%npwwfn, &
3688 201 : gwx_nfftot, method, Sigp%mG0, Cryst, enforce_sym)
3689 :
3690 : ! FFT mesh for sigma_c.
3691 : call setmesh(gmet, Gsph_Max%gvec, gwc_ngfft, Sigp%npwvec, epsm1%npwe, Sigp%npwwfn,&
3692 201 : gwc_nfftot, method, Sigp%mG0, Cryst, enforce_sym, unit=dev_null)
3693 :
3694 : ! ======================================================================
3695 : ! ==== Check for presence of files with core orbitals, for PAW only ====
3696 : ! ======================================================================
3697 201 : Sigp%use_sigxcore=0
3698 201 : if (Dtset%usepaw==1.and.Dtset%gw_sigxcore==1) then
3699 1 : ii = 0
3700 2 : do itypat=1,Cryst%ntypat
3701 1 : string = Psps%filpsp(itypat)
3702 1 : fcore = "CORE_"//TRIM(basename(string))
3703 1 : ic = INDEX (TRIM(string), "/" , back=.TRUE.)
3704 1 : if (ic>0 .and. ic<LEN_TRIM(string)) then
3705 : ! string defines a path, prepend path to fcore
3706 1 : fcore = Psps%filpsp(itypat)(1:ic)//TRIM(fcore)
3707 : end if
3708 2 : if (file_exists(fcore)) then
3709 1 : ii = ii+1
3710 : else
3711 0 : ABI_WARNING(sjoin("HF decoupling is required but cannot find file:", fcore))
3712 : end if
3713 : end do
3714 :
3715 1 : Sigp%use_sigxcore=1
3716 1 : if (ii/=Cryst%ntypat) then
3717 0 : ABI_ERROR("Files with core orbitals not found")
3718 : end if
3719 : end if ! PAW+HF decoupling
3720 : !
3721 : ! ==============================
3722 : ! ==== Extrapolar technique ====
3723 : ! ==============================
3724 201 : Sigp%gwcomp = Dtset%gwcomp
3725 201 : Sigp%gwencomp = Dtset%gwencomp
3726 :
3727 201 : if (Sigp%gwcomp==1) then
3728 19 : write(msg,'(6a,e11.4,a)')ch10,&
3729 19 : 'Using the extrapolar approximation to accelerate convergence',ch10,&
3730 19 : 'with respect to the number of bands included',ch10,&
3731 38 : 'with gwencomp: ',Sigp%gwencomp*Ha_eV,' [eV]'
3732 19 : call wrtout(std_out, msg)
3733 : end if
3734 : !
3735 : ! ===================================
3736 : ! ==== Final compatibility tests ====
3737 : ! ===================================
3738 1433 : ltest=(ks_ebands%mband == Sigp%nbnds .and. ALL(ks_ebands%nband == Sigp%nbnds))
3739 0 : ABI_CHECK(ltest, 'BUG in definition of ks_ebands%nband')
3740 :
3741 : ! FIXME
3742 201 : if (Dtset%symsigma/=0 .and. Sigp%nomegasr/=0) then
3743 1 : if (cryst%idx_spatial_inversion() == 0) then
3744 0 : write(msg,'(5a)')' setup_sigma : BUG :',ch10,&
3745 0 : 'It is not possible to use symsigma /= 0 to calculate the spectral function ',ch10,&
3746 0 : 'when the system does not have the spatial inversion. Please use symsigma=0 '
3747 0 : ABI_WARNING(msg)
3748 : end if
3749 : end if
3750 :
3751 201 : if (mod10==SIG_GW_AC) then
3752 : !if (Sigp%gwcalctyp/=1) ABI_ERROR("Self-consistency with AC not implemented") ! MRM: let's allow it
3753 9 : if (Sigp%gwcomp==1) ABI_ERROR("AC with extrapolar technique not implemented")
3754 : end if
3755 :
3756 : ! For analytic continuation define the number of imaginary frequencies for Sigma
3757 : ! Tests show than more than 12 freqs in the Pade approximant worsen the results!
3758 192 : Sigp%nomegasi = 0
3759 :
3760 : !TODO this should not be done here but in init_sigma_t
3761 : if (mod10 == 1) then
3762 9 : Sigp%nomegasi = abs(Dtset%nomegasi)
3763 9 : iw_mesh_type = "linear"
3764 9 : if (dtset%nomegasi < 0) iw_mesh_type = "minimax"
3765 9 : Sigp%omegasimax = Dtset%omegasimax
3766 9 : Sigp%omegasimin = OMEGASIMIN
3767 27 : ABI_MALLOC(Sigp%omegasi, (Sigp%nomegasi))
3768 :
3769 7 : select case (iw_mesh_type)
3770 : case ("linear")
3771 : ! Linear mesh along the imaginary axis.
3772 7 : domegasi=Sigp%omegasimax/(Sigp%nomegasi-1)
3773 187 : do io=1,Sigp%nomegasi
3774 187 : Sigp%omegasi(io)=CMPLX(zero,(io-1)*domegasi)
3775 : end do
3776 :
3777 : case ("minimax")
3778 : ! Minimax mesh along the imaginary axis.
3779 2 : nbsum = sigp%nbnds
3780 : !call wrtout(std_out, "Computing minimax grid")
3781 :
3782 : ! Compute min/max transition energy taking into account nsppol if any.
3783 6 : te_min = minval(ks_gaps%cb_min - ks_gaps%vb_max)
3784 16 : te_max = maxval(ks_ebands%eig(nbsum,:,:) - ks_ebands%eig(1,:,:))
3785 2 : if (te_min <= tol6) then
3786 0 : te_min = tol6
3787 0 : ABI_ERROR("System is metallic or with a very small fundamental gap! Check energies in WFK file!")
3788 : end if
3789 :
3790 : call gx_minimax_grid(sigp%nomegasi, te_min, te_max, & ! in
3791 : tau_mesh, tau_wgs, iw_mesh, iw_wgs, & ! out args allocated by the routine.
3792 : cosft_wt, cosft_tw, sinft_wt, &
3793 2 : ft_max_error, cosft_duality_error, ierr)
3794 2 : ABI_CHECK(ierr == 0, "Error in gx_minimax_grid")
3795 :
3796 34 : Sigp%omegasi = j_dpc * iw_mesh
3797 :
3798 2 : ABI_FREE_NOCOUNT(tau_mesh)
3799 2 : ABI_FREE_NOCOUNT(tau_wgs)
3800 2 : ABI_FREE_NOCOUNT(iw_mesh)
3801 2 : ABI_FREE_NOCOUNT(iw_wgs)
3802 2 : ABI_FREE_NOCOUNT(cosft_wt)
3803 2 : ABI_FREE_NOCOUNT(cosft_tw)
3804 2 : ABI_FREE_NOCOUNT(sinft_wt)
3805 :
3806 : case default
3807 9 : ABI_ERROR(sjoin("Invalid iw_mesh_type:", iw_mesh_type))
3808 : end select
3809 :
3810 9 : write(msg,'(7a,i3,2(2a,f8.3),a)')ch10,&
3811 9 : ' Parameters for the analytic continuation of Sigma_c(i omega): ',ch10,&
3812 9 : ' Mesh type: = ',trim(iw_mesh_type), ch10, &
3813 9 : ' Number of imaginary frequencies = ',Sigp%nomegasi,ch10,&
3814 228 : ' Min frequency on imag axis (eV) = ',minval(aimag(sigp%omegasi)) * Ha_eV,ch10,&
3815 237 : ' Max frequency on imag axis (eV) = ',maxval(aimag(sigp%omegasi)) * Ha_eV,ch10
3816 9 : call wrtout(units, msg)
3817 :
3818 : ! MRM: do not print for 1-RDM correction
3819 9 : if (Sigp%gwcalctyp /= 21) then
3820 4 : write(msg,'(4a)')ch10,&
3821 4 : ' setup_sigma: calculating Sigma(iw)',&
3822 8 : ' at imaginary frequencies (eV) (Fermi Level set to 0) ',ch10
3823 4 : call wrtout(units, msg)
3824 54 : do io=1,Sigp%nomegasi
3825 50 : write(msg,'(2(f10.3,2x))')Sigp%omegasi(io)*Ha_eV
3826 54 : call wrtout(units, msg)
3827 : end do
3828 4 : call wrtout(units, "")
3829 : endif
3830 :
3831 9 : ltest=(Sigp%omegasimax>0.1d-4.and.Sigp%nomegasi>0)
3832 0 : ABI_CHECK(ltest,'Wrong value of omegasimax or nomegasi')
3833 : !if (Sigp%gwcalctyp/=1) then !
3834 : ! !ABI_ERROR("SC-GW with analytic continuation is not coded") ! MRM: let's allow it
3835 : !end if
3836 : end if
3837 :
3838 201 : if (Sigp%symsigma/=0.and.gwcalctyp>=20) then
3839 0 : ABI_WARNING("SC-GW with symmetries is still under development. Use at your own risk!")
3840 0 : ABI_ERROR("SC-GW requires symsigma == 0 in input. New default in Abinit9 is symsigma 1!")
3841 : end if
3842 :
3843 : ! Setup parameters for Spectral function.
3844 201 : if (Dtset%gw_customnfreqsp/=0) then
3845 1 : Sigp%nomegasr = Dtset%gw_customnfreqsp
3846 1 : ABI_WARNING('Custom grid for spectral function specified. Assuming experienced user.')
3847 1 : if (Dtset%gw_customnfreqsp/=0) then
3848 1 : Dtset%nfreqsp = Dtset%gw_customnfreqsp
3849 1 : ABI_WARNING('nfreqsp has been set to the same number as gw_customnfreqsp')
3850 : end if
3851 : else
3852 200 : Sigp%nomegasr =Dtset%nfreqsp
3853 200 : Sigp%minomega_r=Dtset%freqspmin
3854 200 : Sigp%maxomega_r=Dtset%freqspmax
3855 :
3856 : ! TODO: Mesh should be centered on e0
3857 : if (mod10 == 1 .and. sigp%nomegasr == 0 .and. Sigp%gwcalctyp /= 21) then
3858 : ! Note that in AC computing quantities on the real-axis is really cheap
3859 : ! so we can use very dense meshes without affecting performance.
3860 : ! The default for nfresp and freqspmax is zero.
3861 : ! Here we compute wr_step and nwrt so that we have e0 +- the expected plasmom frequency
3862 : drude_plasmon_freq = sqrt(four_pi * ks_ebands%nelect / cryst%ucvol)
3863 : wmax = dtset%freqspmax; if (abs(wmax) < tol6) wmax = two * drude_plasmon_freq
3864 : !sigp%nomegasr = nint(wmax / (0.05_dp * eV_Ha))
3865 : !if (mod(sigp%nomegasr, 2) == 0) sigp%nomegasr = sigp%nomegasr + 1
3866 : !sigp%maxomega_r = wmax
3867 : end if
3868 :
3869 : end if
3870 :
3871 201 : if (Sigp%nomegasr > 0) then
3872 5 : if (Dtset%gw_customnfreqsp == 0) then
3873 : ! Check
3874 4 : if (Sigp%minomega_r >= Sigp%maxomega_r) then
3875 0 : ABI_ERROR('freqspmin must be smaller than freqspmax!')
3876 : end if
3877 4 : domegas = zero
3878 4 : if(Sigp%nomegasr /= 1) domegas = (Sigp%maxomega_r-Sigp%minomega_r)/(Sigp%nomegasr-1)
3879 :
3880 : !TODO this should be moved to Sr% and done in init_sigma_t
3881 12 : ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr))
3882 804 : do io=1,Sigp%nomegasr
3883 804 : Sigp%omega_r(io) = CMPLX(Sigp%minomega_r + domegas*(io-1),zero)
3884 : end do
3885 4 : write(msg,'(4a,i8,3(2a,f8.3),a)')ch10,&
3886 4 : ' Parameters for the calculation of the spectral function : ',ch10,&
3887 4 : ' Number of points = ',Sigp%nomegasr,ch10,&
3888 4 : ' Min frequency [eV] = ',Sigp%minomega_r*Ha_eV,ch10,&
3889 4 : ' Max frequency [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,&
3890 8 : ' Frequency step [eV] = ',domegas*Ha_eV,ch10
3891 4 : call wrtout(std_out, msg)
3892 : else
3893 10 : Sigp%minomega_r = MINVAL(Dtset%gw_freqsp(:))
3894 10 : Sigp%maxomega_r = MAXVAL(Dtset%gw_freqsp(:))
3895 : !TODO this should be moved to Sr% and done in init_sigma_t
3896 3 : ABI_MALLOC(Sigp%omega_r, (Sigp%nomegasr))
3897 9 : do io=1,Sigp%nomegasr
3898 9 : Sigp%omega_r(io) = CMPLX(Dtset%gw_freqsp(io), zero)
3899 : end do
3900 1 : write(msg,'(4a,i8,2(2a,f8.3),3a)')ch10,&
3901 1 : ' Parameters for the calculation of the spectral function : ',ch10,&
3902 1 : ' Number of points = ',Sigp%nomegasr,ch10,&
3903 1 : ' Min frequency [eV] = ',Sigp%minomega_r*Ha_eV,ch10,&
3904 1 : ' Max frequency [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,&
3905 2 : ' A custom set of frequencies is used! See the input file for values.',ch10
3906 1 : call wrtout(std_out, msg)
3907 : end if
3908 : end if
3909 :
3910 201 : ABI_FREE(val_indices)
3911 201 : call ks_gaps%free()
3912 :
3913 : DBG_EXIT('COLL')
3914 :
3915 804 : end subroutine setup_sigma
3916 : !!***
3917 :
3918 : !----------------------------------------------------------------------
3919 :
3920 : !!****f* ABINIT/sigma_tables
3921 : !! NAME
3922 : !! sigma_tables
3923 : !!
3924 : !! FUNCTION
3925 : !! Build symmetry tables used to speedup self-consistent GW calculations.
3926 : !!
3927 : !! INPUTS
3928 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
3929 : !! [esymm(Kmesh%nibz,Sigp%nsppol)]= Bands_Symmetries
3930 : !!
3931 : !! SiDE EFFECTS
3932 : !! Sigp<sigparams_t>=This routine initializes the tables:
3933 : !! %Sigcij_tab
3934 : !! %Sigxij_tab
3935 : !! that are used to select the matrix elements of the self-energy that have to be calculated.
3936 : !!
3937 : !! SOURCE
3938 :
3939 402 : subroutine sigma_tables(Sigp, Kmesh, esymm)
3940 :
3941 : use m_sigtk, only : sigtk_sigma_tables
3942 :
3943 : !Arguments ------------------------------------
3944 : !scalars
3945 : type(sigparams_t),target,intent(inout) :: Sigp
3946 : type(kmesh_t),intent(in) :: Kmesh
3947 : !arrays
3948 : type(esymm_t),optional,intent(in) :: esymm(Kmesh%nibz, Sigp%nsppol)
3949 :
3950 : !Local variables-------------------------------
3951 : !scalars
3952 : integer :: ikcalc,ik_ibz
3953 : logical :: sigc_is_herm, only_diago
3954 : !arrays
3955 201 : integer,allocatable :: kcalc2ibz(:)
3956 : ! *************************************************************************
3957 :
3958 201 : only_diago = sigp%gwcalctyp < 20
3959 201 : sigc_is_herm = sigp%is_herm()
3960 :
3961 603 : ABI_MALLOC(kcalc2ibz, (sigp%nkptgw))
3962 840 : do ikcalc=1,sigp%nkptgw
3963 639 : ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc))
3964 840 : kcalc2ibz(ikcalc) = ik_ibz
3965 : end do
3966 :
3967 201 : if (present(esymm)) then
3968 : call sigtk_sigma_tables(sigp%nkptgw, kmesh%nibz, sigp%nsppol, sigp%minbnd, sigp%maxbnd, kcalc2ibz, &
3969 0 : only_diago, sigc_is_herm, sigp%sigxij_tab, sigp%sigcij_tab, esymm=esymm)
3970 : else
3971 : call sigtk_sigma_tables(sigp%nkptgw, kmesh%nibz, sigp%nsppol, sigp%minbnd, sigp%maxbnd, kcalc2ibz, &
3972 201 : only_diago, sigc_is_herm, sigp%sigxij_tab, sigp%sigcij_tab)
3973 : end if
3974 :
3975 201 : ABI_FREE(kcalc2ibz)
3976 :
3977 201 : end subroutine sigma_tables
3978 : !!***
3979 :
3980 : !----------------------------------------------------------------------
3981 :
3982 : !!****f* ABINIT/sigma_bksmask
3983 : !! NAME
3984 : !! sigma_bksmask
3985 : !!
3986 : !! FUNCTION
3987 : !! Compute tables for the distribution and the storage of the wavefunctions in the SIGMA code.
3988 : !!
3989 : !! INPUTS
3990 : !! Dtset<type(dataset_type)>=all input variables for this dataset
3991 : !! Sigp<sigparams_t>=Parameters governing the self-energy calculation.
3992 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
3993 : !! nprocs=Total number of MPI processors
3994 : !! my_rank=Rank of this this processor.
3995 : !!
3996 : !! OUTPUT
3997 : !! my_spins(:)=Pointer to NULL in input. In output: list of spins treated by this node.
3998 : !! bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state.
3999 : !! keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space.
4000 : !! ierr=Exit status.
4001 : !!
4002 : !! SOURCE
4003 :
4004 201 : subroutine sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr)
4005 :
4006 : !Arguments ------------------------------------
4007 : !scalars
4008 : integer,intent(in) :: my_rank,nprocs
4009 : integer,intent(out) :: ierr
4010 : type(Dataset_type),intent(in) :: Dtset
4011 : type(sigparams_t),intent(in) :: Sigp
4012 : type(kmesh_t),intent(in) :: Kmesh
4013 : !arrays
4014 : integer,allocatable,intent(out) :: my_spins(:)
4015 : logical,intent(out) :: bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
4016 : logical,intent(out) :: keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
4017 :
4018 : !Local variables-------------------------------
4019 : !scalars
4020 : integer :: my_nspins,my_maxb,my_minb,isp,spin,band,nsppol,rank_spin
4021 : logical :: store_ur
4022 : !arrays
4023 402 : integer :: tmp_spins(Sigp%nsppol),nprocs_spin(Sigp%nsppol)
4024 : ! *************************************************************************
4025 :
4026 201 : ierr=0; nsppol=Sigp%nsppol
4027 :
4028 : ! List of spins for each node, number of processors per each spin
4029 : ! and the MPI rank in the "spin" communicator.
4030 201 : my_nspins=nsppol
4031 603 : ABI_MALLOC(my_spins, (nsppol))
4032 1017 : my_spins= [(isp, isp=1,nsppol)]
4033 406 : nprocs_spin = nprocs; rank_spin = my_rank
4034 :
4035 201 : if (nsppol==2 .and. nprocs>1) then
4036 : ! Distribute spins (optimal distribution if nprocs is even)
4037 0 : nprocs_spin(1) = nprocs/2
4038 0 : nprocs_spin(2) = nprocs - nprocs/2
4039 0 : my_nspins=1
4040 0 : my_spins(1)=1
4041 0 : if (my_rank+1>nprocs/2) then
4042 : ! I will treat spin=2, compute shifted rank.
4043 0 : my_spins(1)=2
4044 0 : rank_spin = my_rank - nprocs/2
4045 : end if
4046 : end if
4047 :
4048 201 : store_ur = (MODULO(Dtset%gwmem,10)==1)
4049 61415 : keep_ur=.FALSE.; bks_mask=.FALSE.
4050 :
4051 243 : select case (Dtset%gwpara)
4052 : case (1)
4053 : ! Parallelization over transitions **without** memory distributions (Except for the spin).
4054 42 : my_minb=1; my_maxb=Sigp%nbnds
4055 42 : if (dtset%ucrpa>0) then
4056 0 : bks_mask(my_minb:my_maxb,:,:)=.TRUE.
4057 0 : if (store_ur) keep_ur(my_minb:my_maxb,:,:)=.TRUE.
4058 : else
4059 84 : do isp=1,my_nspins
4060 42 : spin = my_spins(isp)
4061 4934 : bks_mask(my_minb:my_maxb,:,spin)=.TRUE.
4062 4976 : if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE.
4063 : end do
4064 : end if
4065 :
4066 : case (2)
4067 : ! Distribute bands and spin (alternating planes of bands)
4068 322 : do isp=1,my_nspins
4069 163 : spin = my_spins(isp)
4070 5403 : do band=1,Sigp%nbnds
4071 5244 : if (xmpi_distrib_with_replicas(band,Sigp%nbnds,rank_spin,nprocs_spin(spin))) then
4072 : !if (MODULO(band-1,nprocs_spin(spin))==rank_spin) then
4073 26347 : bks_mask(band,:,spin)=.TRUE.
4074 25895 : if (store_ur) keep_ur(band,:,spin)=.TRUE.
4075 : end if
4076 : end do
4077 : end do
4078 :
4079 : #if 0
4080 : ! Each node store the full set of occupied states to speed up Sigma_x.
4081 : do isp=1,my_nspins
4082 : spin = my_spins(isp)
4083 : do ik_ibz=1,Kmesh%nibz
4084 : ks_iv=ks_vbik(ik_ibz,spin) ! Valence index for this (k,s)
4085 : bks_mask(1:ks_iv,:,spin)=.TRUE.
4086 : if (store_ur) keep_ur(1:ks_iv,:,spin)=.TRUE.
4087 : end do
4088 : end do
4089 : #endif
4090 :
4091 : case default
4092 0 : ABI_WARNING("Wrong value for gwpara")
4093 201 : ierr = 1
4094 : end select
4095 :
4096 : ! Return my_spins with correct size.
4097 406 : tmp_spins(1:my_nspins) = my_spins(1:my_nspins)
4098 :
4099 201 : ABI_FREE(my_spins)
4100 603 : ABI_MALLOC(my_spins, (my_nspins))
4101 607 : my_spins = tmp_spins(1:my_nspins)
4102 :
4103 201 : end subroutine sigma_bksmask
4104 : !!***
4105 :
4106 : !!****f* ABINIT/paw_qpscgw
4107 : !! NAME
4108 : !! paw_qpscgw
4109 : !!
4110 : !! FUNCTION
4111 : !! This routine is called during QP self-consistent GW calculations. It calculates the new QP on-site quantities
4112 : !! using the QP amplitudes read from the QPS file.
4113 : !!
4114 : !! INPUTS
4115 : !! Wfd<wfdgw_t>=Datatype gathering data on QP amplitudes.
4116 : !! nscf=Number of QPSCGW iterations done so far (read from the QPS file).
4117 : !! nfftf=Number of points in the fine FFT grid.
4118 : !! ngfft(18)=information on the fine FFT grid used for densities and potentials.
4119 : !! Dtset<dataset_type>=All input variables for this dataset.
4120 : !! Cryst<crystal_t>=Info on unit cell and symmetries.
4121 : !! Kmesh<kmesh_t>=Structure describing the k-point sampling.
4122 : !! Psps<Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file
4123 : !! Pawang<pawang_type>=PAW angular mesh and related data.
4124 : !! Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data
4125 : !! Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
4126 : !! Pawfgrtab(natom)<Pawfgrtab_type>= For PAW, various arrays giving data related to fine grid for a given atom.
4127 : !! prev_Pawrhoij(Cryst%natom))<Pawrhoij_type>=Previous QP rhoij used for mixing if nscf>0 and rhoqpmix/=one.
4128 : !! MPI_enreg=Information about MPI parallelization.
4129 : !! qp_ebands<ebands_t>=QP band structure.
4130 : !! QP_energies<Energies_type>=Simple datastructure to gather all part of total energy.
4131 : !! nhatgrdim= 0 if pawgrnhat array is not used ; 1 otherwise
4132 : !!
4133 : !! OUTPUT
4134 : !! QP_pawrhoij(Cryst%natom))<Pawrhoij_type>=on-site densities calculated from the QP amplitudes.
4135 : !! qp_nhat(nfftf,Dtset%nspden)=Compensation charge density calculated from the QP amplitudes.
4136 : !! qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim)=Derivatives of the QP nhat on fine rectangular grid (and derivatives).
4137 : !! qp_compch_sph=QP compensation charge integral inside spheres computed over spherical meshes.
4138 : !! qp_compch_fft=QP compensation charge inside spheres computed over fine fft grid.
4139 : !! QP_paw_ij(Cryst%natom)<Paw_ij_type>=Non-local D_ij strengths of the QP Hamiltonian.
4140 : !! QP_paw_an(Cryst%natom)<Paw_an_type>=Various arrays related to the Hamiltonian given
4141 : !! on ANgular mesh or ANgular moments.
4142 : !!
4143 : !! SOURCE
4144 :
4145 0 : subroutine paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,qp_ebands, &
4146 0 : Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij, &
4147 0 : QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim, &
4148 0 : qp_nhatgr,qp_compch_sph,qp_compch_fft)
4149 :
4150 : !Arguments ------------------------------------
4151 : !scalars
4152 : integer,intent(in) :: nfftf,nscf,nhatgrdim
4153 : real(dp),intent(out) :: qp_compch_fft,qp_compch_sph
4154 : type(kmesh_t),intent(in) :: Kmesh
4155 : type(crystal_t),intent(in) :: Cryst
4156 : type(Dataset_type),intent(in) :: Dtset
4157 : type(Pseudopotential_type),intent(in) :: Psps
4158 : type(Pawang_type),intent(in) :: Pawang
4159 : type(ebands_t),intent(in) :: qp_ebands
4160 : type(Energies_type),intent(inout) :: QP_energies
4161 : type(wfdgw_t),intent(inout) :: Wfd
4162 : !arrays
4163 : integer,intent(in) :: ngfftf(18)
4164 : real(dp),intent(out) :: qp_nhat(nfftf,Dtset%nspden)
4165 : real(dp),intent(out) :: qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim)
4166 : type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom)
4167 : type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
4168 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
4169 : type(Pawrhoij_type),intent(inout) :: prev_Pawrhoij(Cryst%natom)
4170 : type(Pawrhoij_type),intent(out) :: QP_pawrhoij(Cryst%natom)
4171 : type(Paw_ij_type),intent(out) :: QP_paw_ij(Cryst%natom)
4172 : type(Paw_an_type),intent(inout) :: QP_paw_an(Cryst%natom)
4173 :
4174 : !Local variables ------------------------------
4175 : !scalars
4176 : integer,parameter :: ipert0 = 0, idir0 = 0, optrhoij1 = 1
4177 : integer :: choice,cplex,cplex_rhoij,has_dijU,has_dijso,iat,ider
4178 : integer :: izero,nkxc1,nspden_rhoij,nzlmopt, option,usexcnhat
4179 : character(len=500) :: msg
4180 : real(dp) :: el_temp
4181 : !arrays
4182 : real(dp) :: k0(3)
4183 : !************************************************************************
4184 :
4185 : ABI_UNUSED(Kmesh%nibz)
4186 : !
4187 : ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation)
4188 : ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation)
4189 : usexcnhat=MAXVAL(Pawtab(:)%usexcnhat)
4190 : !
4191 : ! Calculate new rhoij_qp from updated Cprj_ibz, note use_rhoij_=1.
4192 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
4193 0 : nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc)
4194 : call pawrhoij_alloc(QP_pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,&
4195 0 : pawtab=Pawtab,use_rhoij_=1,use_rhoijres=1)
4196 :
4197 : ! FIXME kptop should be passed via Kmesh, in GW time reversal is always assumed.
4198 0 : call wfd%pawrhoij(Cryst,qp_ebands,Dtset%kptopt,QP_pawrhoij,Dtset%pawprtvol)
4199 : !
4200 : ! Symmetrize QP $\rho_{ij}$.
4201 0 : choice=1
4202 : call pawrhoij_symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,&
4203 : Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij1,Pawang,Dtset%pawprtvol,&
4204 0 : Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
4205 :
4206 : ! ======================
4207 : ! ==== Make QP nhat ====
4208 : ! ======================
4209 0 : cplex=1; ider=2*nhatgrdim; izero=0; k0(:)=zero
4210 :
4211 : call pawmknhat(qp_compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,&
4212 : Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
4213 : Pawfgrtab,qp_nhatgr,qp_nhat,QP_Pawrhoij,QP_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
4214 0 : Cryst%ucvol,Dtset%usewvl,Cryst%xred)
4215 :
4216 : ! Allocate quantities related to the PAW spheres for the QP Hamiltonian.
4217 : ! TODO call paw_ij_init in scfcv and respfn, fix small issues
4218 0 : has_dijso=Dtset%pawspnorb; has_dijU=merge(0,1,Dtset%usepawu==0)
4219 :
4220 0 : call paw_ij_nullify(QP_paw_ij)
4221 : call paw_ij_init(QP_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
4222 : Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
4223 : has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
4224 0 : has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1)
4225 :
4226 0 : call paw_an_nullify(QP_paw_an); nkxc1=0 ! No kernel
4227 : call paw_an_init(QP_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
4228 0 : cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1)
4229 :
4230 : ! =====================================================
4231 : ! ==== Optional mixing of the PAW onsite densities ====
4232 : ! =====================================================
4233 0 : if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then
4234 0 : write(msg,'(a,f6.3)')' sigma: mixing on-site QP rho_ij densities using rhoqpmix= ',Dtset%rhoqpmix
4235 0 : call wrtout(std_out, msg)
4236 : ! qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor)
4237 :
4238 0 : call pawrhoij_unpack(QP_Pawrhoij) ! Unpack new QP %rhoijp
4239 0 : call pawrhoij_unpack(prev_Pawrhoij) ! Unpack previous QP %rhoijp
4240 :
4241 0 : do iat=1,Cryst%natom
4242 : QP_pawrhoij(iat)%rhoij_ = prev_Pawrhoij(iat)%rhoij_ &
4243 0 : + Dtset%rhoqpmix * (QP_pawrhoij(iat)%rhoij_ - prev_pawrhoij(iat)%rhoij_)
4244 :
4245 0 : prev_pawrhoij(iat)%use_rhoij_=0
4246 0 : ABI_FREE(prev_pawrhoij(iat)%rhoij_)
4247 : end do
4248 :
4249 : ! Re-Symmetrize mixed QP $\rho_{ij}$.
4250 : choice=1
4251 : call pawrhoij_symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,&
4252 : Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij1,Pawang,Dtset%pawprtvol,&
4253 0 : Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat)
4254 : end if
4255 :
4256 0 : do iat=1,Cryst%natom
4257 0 : QP_pawrhoij(iat)%use_rhoij_=0
4258 0 : ABI_FREE(QP_pawrhoij(iat)%rhoij_)
4259 : end do
4260 : !
4261 : ! =================================================================================
4262 : ! ==== Evaluate on-site energies, potentials, densities using (mixed) QP rhoij ====
4263 : ! =================================================================================
4264 : ! * Initialize also "lmselect" (index of non-zero LM-moments of densities).
4265 :
4266 0 : nzlmopt=-1; option=0; qp_compch_sph=greatest_real
4267 :
4268 : ! Get electronic temperature from dtset
4269 0 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
4270 :
4271 : call pawdenpot(qp_compch_sph,el_temp,Cryst%gprimd,ipert0,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
4272 : Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,QP_paw_an,QP_paw_an,QP_energies%paw,&
4273 : QP_paw_ij,Pawang,Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
4274 0 : Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,Cryst%xred,Cryst%ucvol,Psps%znuclpsp,Dtset%spinaxis)
4275 :
4276 0 : end subroutine paw_qpscgw
4277 : !!***
4278 :
4279 : !!****f* ABINIT/setup_vcp
4280 : !! NAME
4281 : !! setup_vcp
4282 : !!
4283 : !! FUNCTION
4284 : !! Initialize the Vcp and compute the corresponding residue.
4285 : !!
4286 : !! INPUTS
4287 : !! Dtset<type(dataset_type)>=all input variables for this dataset
4288 : !! Gsph_c<gsphere_t>=Info on the G-sphere for W and Sigma_c
4289 : !! Gsph_x<gsphere_t>=Info on the G-sphere for and Sigma_x
4290 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
4291 : !! Qmesh <kmesh_t>=Structure describing the q-point sampling.
4292 : !! Cryst<crystal_t>=Info on unit cell and symmetries.
4293 : !! comm=Information about the xmpi_world
4294 : !!
4295 : !! OUTPUT
4296 : !! Vcp_full<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique.
4297 : !! Vcp_ks<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique.
4298 : !! coef_hyb=real variable containing the amount of GLOBAL hybridization
4299 :
4300 4 : subroutine setup_vcp(Vcp_ks,Vcp_full,Dtset,Gsph_x,Gsph_c,Cryst,Qmesh,Kmesh,coef_hyb,comm)
4301 :
4302 : !Arguments ------------------------------------
4303 : !scalars
4304 : type(Dataset_type),intent(in) :: Dtset
4305 : type(gsphere_t),intent(in) :: Gsph_x,Gsph_c
4306 : type(crystal_t),intent(in) :: Cryst
4307 : type(kmesh_t),intent(in) :: Kmesh, Qmesh
4308 : type(vcoul_t),intent(inout) :: Vcp_ks, Vcp_full
4309 : integer,intent(in) :: comm
4310 : real(dp),intent(inout) :: coef_hyb
4311 :
4312 : !Local variables-------------------------------
4313 : !scalars
4314 : integer :: usefock_ixc,nqlwl,xclevel_ixc,icsing_eff
4315 : real(dp) :: rcut
4316 : !arrays
4317 4 : real(dp),allocatable :: qlwl(:,:)
4318 : !************************************************************************
4319 :
4320 : ! Build Vcp_ks and Vcp_full
4321 4 : if (Dtset%gw_nqlwl==0) then
4322 4 : nqlwl=1
4323 4 : ABI_MALLOC(qlwl,(3,nqlwl))
4324 16 : qlwl(:,1)= GW_Q0_DEFAULT
4325 : else
4326 0 : nqlwl=Dtset%gw_nqlwl
4327 0 : ABI_MALLOC(qlwl,(3,nqlwl))
4328 0 : qlwl(:,:)=Dtset%gw_qlwl(:,1:nqlwl)
4329 : end if
4330 4 : rcut=Dtset%gw_rcut
4331 4 : icsing_eff=Dtset%gw_icutcoul
4332 : ! 1st part: Use a Vcp_full to compute the full Coulomb interaction for NOs
4333 4 : if (Gsph_x%ng > Gsph_c%ng) then
4334 : call Vcp_full%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icsing_eff,Dtset%vcutgeo,&
4335 4 : Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
4336 : else
4337 : call Vcp_full%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icsing_eff,Dtset%vcutgeo,&
4338 0 : Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
4339 : end if
4340 : ! 2nd part: Use a Vcp_ks to compute the Coulomb interaction already present in the Fock part of the Kohn-Sham Hamiltonian
4341 4 : coef_hyb=zero
4342 4 : call get_xclevel(Dtset%ixc,xclevel_ixc,usefock_ixc) ! usefock is 0 for 0.0 Fock contrib in KS-GS calc
4343 4 : if (usefock_ixc==1)then
4344 0 : if (abs(Dtset%hyb_mixing)>tol8) then
4345 0 : coef_hyb=abs(Dtset%hyb_mixing)
4346 : if(abs(coef_hyb+999.0d0)<tol8) then ! HF
4347 : coef_hyb=1.0d0
4348 : end if
4349 0 : else if(abs(Dtset%hyb_mixing_sr)>tol8)then
4350 0 : coef_hyb=abs(Dtset%hyb_mixing_sr)
4351 0 : icsing_eff=5
4352 : end if
4353 0 : if (abs(rcut)<tol6 .and. abs(Dtset%hyb_range_fock)>tol8) then
4354 0 : rcut=one/Dtset%hyb_range_fock
4355 : end if
4356 : end if
4357 4 : if (Gsph_x%ng > Gsph_c%ng) then
4358 : call Vcp_ks%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icsing_eff,Dtset%vcutgeo,&
4359 4 : Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm)
4360 : else
4361 : call Vcp_ks%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icsing_eff,Dtset%vcutgeo,&
4362 0 : Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm)
4363 : end if
4364 4 : ABI_FREE(qlwl)
4365 : ! Now compute the "residual" Coulomb interactions
4366 4 : ABI_COMMENT("From now on, Vcp_ks%vc_sqrt_resid contains the amount of exact exchange present on the GS.")
4367 2216 : Vcp_ks%vc_sqrt_resid=sqrt(coef_hyb*Vcp_ks%vc_sqrt**2)
4368 4 : Vcp_ks%i_sz_resid=coef_hyb*Vcp_ks%i_sz
4369 :
4370 4 : end subroutine setup_vcp
4371 : !!***
4372 :
4373 0 : end module m_sigma_driver
4374 : !!***
|
