TABLE OF CONTENTS
ABINIT/cchi0 [ Functions ]
NAME
cchi0
FUNCTION
Main calculation of the independent-particle susceptibility chi0 for qpoint != 0
INPUTS
use_tr=If .TRUE. valence states are stored in Wfs_val and only resonant transitions are calculated (time reversal is assumed) Dtset <type(dataset_type)>=all input variables in this dataset Cryst<crystal_t>= data type gathering info on symmetries and unit cell %natom=number of atoms %nsym=number of symmetries %xred(3,natom)=reduced coordinated of atoms %typat(natom)=type of each atom %rprimd(3,3)=dimensional primitive translations in real space (bohr) %timrev= 2 if time reversal can be used, 1 otherwise qpoint(3)=reciprocal space coordinates of the q wavevector Ep<type(em1params_t_type)>= Parameters related to the calculation of the inverse dielectric matrix. %nbnds=number of bands summed over %npwe=number of planewaves for the irreducible polarizability X^0_GGp %npwvec=maximum number of G vectors used to define the dimension of some arrays e.g igfft %nsppol=1 for unpolarized, 2 for spin-polarized %nomega=total number of frequencies in X^0 (both real and imaginary) %nomegasf=number of real frequencies used to sample the imaginary part of X^0 (spectral method) %spmeth=1 if we use the spectral method, 0 for standard Adler-Wiser expression %spsmear=gaussian broadening used to approximate the delta distribution %zcut=small imaginary shift to avoid poles in X^0 Psps <type(pseudopotential_type)>=variables related to pseudopotentials Kmesh <kmesh_t>= datatype gathering parameters related to the k-point sampling %nibz=number of k-points in the IBZ %nbz=number of k-points in the BZ %bz(3,nbz)=reduced coordinates for k-points in the full Brillouin zone %ibz(3,nibz)=reduced coordinates for k-points in the irreducible wedge %tab(nbz)=mapping between a kpt in the BZ (array bz) and the irred point in the array ibz %tabi(nbz)= -1 if inversion is needed to obtain this particular kpt in the BZ, 1 means identity %tabo(nbz)= for each point in the BZ, the index of the symmetry operation S in reciprocal space which rotates k_IBZ onto \pm k_BZ (depending on tabi) %tabp(nbz)= For each k_BZ, it gives the phase factors associated to non-symmorphic operations, i.e e^{-i 2 \pi k_IBZ \cdot R{^-1}t} == e{-i 2\pi k_BZ cdot t} where : \transpose R{-1}=S and (S k_IBZ) = \pm k_BZ (depending on ktabi) %tabr(nfftot,nbz) For each point r on the real mesh and for each k-point in the BZ, tabr gives the index of (R^-1 (r-t)) in the FFT array where R=\transpose S^{-1} and k_BZ=S k_IBZ. t is the fractional translation associated to R Gsph_epsG0<gsphere_t data type> The G-sphere used to describe chi0/eps. (including umklapp G0 vectors) %ng=number of G vectors for chi0 %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G %phmGt(Ep%npwe,nsym)=phase factors e^{-iG \cdot t} needed to symmetrize oscillator matrix elements and epsilon %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1) %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$). nbvw=number of bands in the arrays wfrv ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths (see input variable) nfftot_gw=Total number of points in the GW FFT grid Ltg_q<Little group>=Data type gathering information on the little group of the q-points. Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data Pawang<pawang_type> angular mesh discretization and related data: qp_ebands<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities) %mband=MAX number of bands over k-points and spin (==Ep%nbnds) %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave. Wfd<wfdgw_t>=Object used to access the wavefunctions
OUTPUT
chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix at wavevector qpoint and each frequeny defined by Ep%omega and Ep%nomega
SOURCE
1139 subroutine cchi0(use_tr,Dtset,Cryst,qpoint,Ep,Psps,Kmesh,qp_ebands,Gsph_epsG0,& 1140 Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftot_gw,ngfftf,nfftf_tot,& 1141 chi0,ktabr,ktabrf,Ltg_q,chi0_sumrule,Wfd,Wfdf,wan) 1142 1143 !Arguments ------------------------------------ 1144 !scalars 1145 integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot 1146 logical,intent(in) :: use_tr 1147 type(ebands_t),target,intent(in) :: qp_ebands 1148 type(kmesh_t),intent(in) :: Kmesh 1149 type(crystal_t),intent(in) :: Cryst 1150 type(Dataset_type),intent(in) :: Dtset 1151 type(em1params_t),intent(in) :: Ep 1152 type(gsphere_t),intent(in) :: Gsph_epsG0 1153 type(littlegroup_t),intent(in) :: Ltg_q 1154 type(Pawang_type),intent(in) :: Pawang 1155 type(Pseudopotential_type),intent(in) :: Psps 1156 type(wfdgw_t),target,intent(inout) :: Wfd,Wfdf 1157 !arrays 1158 integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz) 1159 integer,intent(in) :: ngfft_gw(18),ngfftf(18) 1160 real(dp),intent(in) :: qpoint(3) 1161 real(dp),intent(out) :: chi0_sumrule(Ep%npwe) 1162 complex(gwpc),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega) 1163 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 1164 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 1165 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 1166 type(plowannier_type),intent(inout) :: wan 1167 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom) 1168 1169 !Local variables ------------------------------ 1170 !scalars 1171 integer,parameter :: tim_fourdp1 = 1, two_poles = 2, one_pole = 1, ndat1 = 1 1172 integer :: bandinf,bandsup,dim_rtwg,band1,band2,ierr,band1c,band2c 1173 integer :: ig1,ig2,iat1,iat2,iat,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz 1174 integer :: io,iomegal,iomegar,ispinor1,ispinor2,isym_k,itypatcor,nfft,il1,il2 1175 integer :: isym_kmq,itim_k,itim_kmq,m1,m2,my_wl,my_wr,size_chi0 1176 integer :: nfound,nkpt_summed,nspinor,nsppol,mband 1177 integer :: comm,gw_mgfft,use_padfft,gw_fftalga,lcor,mgfftf,use_padfftf 1178 integer :: my_nbbp,my_nbbpks,spin,nbmax,dummy 1179 real(dp) :: cpu_time,wall_time,gflops 1180 real(dp) :: deltaeGW_b1kmq_b2k,deltaeGW_enhigh_b2k,deltaf_b1kmq_b2k 1181 real(dp) :: e_b1_kmq,en_high,fac,fac2,fac3,f_b1_kmq,factor,max_rest,min_rest,my_max_rest 1182 real(dp) :: my_min_rest,numerator,spin_fact,weight,wl,wr 1183 real(dp) :: gw_gsq,memreq 1184 complex(dpc) :: ph_mkmqt,ph_mkt 1185 complex(gwpc) :: local_czero_gw 1186 logical :: qzero,isirred_k,isirred_kmq,luwindow,is_metallic 1187 character(len=500) :: msg,allup 1188 type(gsphere_t) :: Gsph_FFT 1189 !arrays 1190 integer :: G0(3),umklp_k(3),umklp_kmq(3) 1191 integer :: ucrpa_bands(2) 1192 integer :: wtk_ltg(Kmesh%nbz),got(Wfd%nproc) 1193 integer,allocatable :: tabr_k(:),tabr_kmq(:),tabrf_k(:),tabrf_kmq(:) 1194 integer,allocatable :: igfftepsG0(:),gspfft_igfft(:),igfftepsG0f(:) 1195 integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:) 1196 integer,allocatable :: bbp_ks_distrb(:,:,:,:) 1197 real(dp) :: kbz(3),kmq_bz(3),spinrot_k(4),spinrot_kmq(4),q0(3),tsec(2) 1198 real(dp),ABI_CONTIGUOUS pointer :: qp_eig(:,:,:),qp_occ(:,:,:) 1199 real(dp),allocatable :: omegasf(:) 1200 complex(dpc),allocatable :: green_enhigh_w(:),green_w(:),kkweight(:,:) 1201 complex(gwpc),allocatable :: sf_chi0(:,:,:),rhotwg(:) 1202 complex(gwpc),allocatable :: ur1_kmq_ibz(:),ur2_k_ibz(:),wfwfg(:) 1203 complex(gwpc),allocatable :: usr1_kmq(:),ur2_k(:) 1204 complex(gwpc),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:) 1205 complex(gwpc),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:) 1206 complex(dpc), allocatable :: coeffW_BZ(:,:,:,:,:,:) 1207 logical,allocatable :: bbp_mask(:,:) 1208 type(pawcprj_type),allocatable :: Cprj1_kmq(:,:),Cprj2_k(:,:) 1209 type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:) 1210 !************************************************************************ 1211 1212 DBG_ENTER("COLL") 1213 1214 call timab(331,1,tsec) ! cchi0 1215 call cwtime(cpu_time,wall_time,gflops,"start") 1216 1217 nsppol = Wfd%nsppol; nspinor = Wfd%nspinor 1218 is_metallic = ebands_has_metal_scheme(qp_ebands) 1219 1220 ucrpa_bands(1)=dtset%ucrpa_bands(1) 1221 ucrpa_bands(2)=dtset%ucrpa_bands(2) 1222 luwindow=.false. 1223 if(abs(dtset%ucrpa_window(1)+1_dp)>tol8.or.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) then 1224 luwindow=.true. 1225 endif 1226 !write(6,*)"ucrpa_bands",ucrpa_bands; write(6,*)"ucrpa_window",dtset%ucrpa_window; write(6,*)"luwindow",luwindow 1227 1228 ! For cRPA calculation of U: read forlb.ovlp 1229 if(dtset%ucrpa>=1 .AND. dtset%plowan_compute <10) then 1230 call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,& 1231 nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands) 1232 endif 1233 ! End of reading forlb.ovlp 1234 1235 if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,ngfft_gw) 1236 gw_mgfft = MAXVAL(ngfft_gw(1:3)) 1237 gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10) 1238 1239 if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3)) 1240 1241 ! == Copy some values === 1242 comm = Wfd%comm 1243 mband = Wfd%mband 1244 nfft = Wfd%nfft 1245 ABI_CHECK(Wfd%nfftot==nfftot_gw,"Wrong nfftot_gw") 1246 1247 dim_rtwg=1 !; if (nspinor==2) dim_rtwg=2 ! can reduce size depending on Ep%nI and Ep%nj 1248 size_chi0 = Ep%npwe*Ep%nI*Ep%npwe*Ep%nJ*Ep%nomega 1249 1250 qp_eig => qp_ebands%eig; qp_occ => qp_ebands%occ 1251 1252 ! Initialize the completeness correction 1253 if (Ep%gwcomp==1) then 1254 en_high=MAXVAL(qp_eig(Ep%nbnds,:,:)) + Ep%gwencomp 1255 write(msg,'(a,f8.2,a)')' Using completeness correction with the energy ',en_high*Ha_eV,' [eV]' 1256 call wrtout(std_out, msg) 1257 1258 ! Allocation of wfwfg and green_enhigh_w moved inside openmp loop 1259 ! Init the largest G-sphere contained in the FFT box for the wavefunctions. 1260 call Gsph_FFT%in_fftbox(Cryst,Wfd%ngfft) 1261 1262 !call Gsph_FFT%print(unit=std_out,prtvol=10) 1263 1264 ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng)) 1265 ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2)) 1266 1267 ! Mapping between G-sphere and FFT box. 1268 call Gsph_FFT%fft_tabs((/0,0,0/),Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft) 1269 ABI_FREE(dummy_gbound) 1270 1271 if (Psps%usepaw==1) then ! * Prepare the onsite contributions on the GW FFT mesh. 1272 ABI_MALLOC(gw_gfft,(3,nfft)) 1273 q0=zero 1274 call get_gftt(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! Get the set of plane waves in the FFT Box. 1275 ABI_MALLOC(Pwij_fft,(Psps%ntypat)) 1276 call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 1277 end if 1278 end if 1279 1280 ! Setup weights (2 for spin unpolarized system, 1 for polarized). 1281 ! spin_fact is used to normalize the occupation factors to one. Consider also the AFM case. 1282 select case (nsppol) 1283 case (1) 1284 weight = two / Kmesh%nbz; spin_fact = half 1285 if (Wfd%nspden==2) then 1286 weight = one / Kmesh%nbz; spin_fact = half 1287 end if 1288 if (nspinor==2) then 1289 weight = one / Kmesh%nbz; spin_fact = one 1290 end if 1291 case (2) 1292 weight = one / Kmesh%nbz; spin_fact = one 1293 case default 1294 ABI_BUG("Wrong nsppol") 1295 end select 1296 1297 ! Weight for points in the IBZ_q. 1298 wtk_ltg(:) = 1 1299 if (Ep%symchi==1) then 1300 do ik_bz=1,Ltg_q%nbz 1301 wtk_ltg(ik_bz)=0 1302 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only k-points in the IBZ_q. 1303 wtk_ltg(ik_bz)=SUM(Ltg_q%wtksym(:,:,ik_bz)) 1304 end do 1305 end if 1306 1307 write(msg,'(a,i2,2a,i2)')& 1308 ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,& 1309 ' Using symmetries to sum only over the IBZ_q = ',Ep%symchi 1310 call wrtout(std_out, msg) 1311 1312 if (use_tr) then 1313 ! Special care has to be taken in metals and/or spin dependent systems 1314 ! as Wfs_val might contain unoccupied states. 1315 call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.') 1316 else 1317 call wrtout(std_out,' Using slow algorithm without time reversal symmetry.') 1318 end if 1319 1320 ! TODO this table can be calculated for each k-point 1321 my_nbbpks=0; allup="All"; got=0 1322 ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol)) 1323 call wrtout(std_out, sjoin(' Memory needed for bbp_ks_distrb: ', & 1324 ftoa(four*mband**2*Kmesh%nbz*nsppol*b2Mb, fmt="f8.1"), ' [Mb] <<< MEM')) 1325 1326 ABI_MALLOC(bbp_mask,(mband, mband)) 1327 1328 do spin=1,nsppol 1329 do ik_bz=1,Kmesh%nbz 1330 if (Ep%symchi == 1) then 1331 if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only IBZ_q 1332 end if 1333 1334 ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz. 1335 call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k) 1336 1337 ! Get index of k-q in the BZ, stop if not found as the weight=one/nkbz is not correct. 1338 call kmesh%get_BZ_diff(kbz,qpoint,ikmq_bz,g0,nfound) 1339 ABI_CHECK(nfound==1,"Check kmesh") 1340 1341 ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz. 1342 call kmesh%get_BZ_item(ikmq_bz,kmq_bz,ikmq_ibz,isym_kmq,itim_kmq) 1343 1344 call chi0_bbp_mask(ikmq_ibz, ik_ibz, spin, spin_fact, use_tr, & 1345 ep%gwcomp, ep%spmeth, ep%nbnds, mband, qp_ebands, bbp_mask) 1346 1347 call wfd%distribute_kb_kpbp(ikmq_ibz,ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got=got,bbp_mask=bbp_mask) 1348 my_nbbpks = my_nbbpks + my_nbbp 1349 end do 1350 end do 1351 1352 ABI_FREE(bbp_mask) 1353 1354 write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0." 1355 call wrtout(std_out, msg) 1356 1357 if (Psps%usepaw==1) then 1358 ABI_MALLOC(Pwij,(Psps%ntypat)) 1359 call pawpwij_init(Pwij,Ep%npwepG0,qpoint,Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 1360 ! Allocate statements moved to inside openmp loop 1361 end if 1362 1363 SELECT CASE (Ep%spmeth) 1364 CASE (0) 1365 call wrtout(std_out,' Calculating chi0(q,omega,G,G")') 1366 ! Allocation of green_w moved inside openmp loop 1367 1368 CASE (1, 2) 1369 call wrtout(std_out,' Calculating Im chi0(q,omega,G,G")') 1370 1371 ! Find Max and min resonant transitions for this q, report also treated by this proc. 1372 call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,& 1373 & max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_eig,qp_occ,qpoint,bbp_ks_distrb) 1374 ! 1375 ! Calculate frequency dependent weights for Hilbert transform. 1376 ABI_MALLOC(omegasf,(Ep%nomegasf)) 1377 ABI_MALLOC(kkweight,(Ep%nomegasf,Ep%nomega)) 1378 !my_wl=1; my_wr=Ep%nomegasf 1379 call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,& 1380 & 0,Ep%zcut,zero,my_wl,my_wr,kkweight) 1381 1382 if (.not.use_tr) then 1383 ABI_BUG('spectral method requires time-reversal') 1384 end if 1385 1386 memreq = two*gwpc*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb 1387 write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwpc*Ep%npwe**2*b2Mb,' [Mb]' 1388 call wrtout(std_out,msg) 1389 write(msg,'(a,f10.4,a)')' memory required by sf_chi0: ',memreq,' [Gb]' 1390 call wrtout(std_out,msg) 1391 if (memreq > two) then 1392 ABI_WARNING(' Memory required for sf_chi0 is larger than 2.0 Gb!') 1393 end if 1394 ABI_MALLOC_OR_DIE(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr) 1395 sf_chi0=czero_gw 1396 1397 CASE DEFAULT 1398 ABI_BUG("Wrong spmeth") 1399 END SELECT 1400 1401 nkpt_summed=Kmesh%nbz 1402 if (Ep%symchi == 1) then 1403 nkpt_summed = Ltg_q%nibz_ltg 1404 call Ltg_q%print(std_out, Dtset%prtvol) 1405 end if 1406 1407 write(msg,'(a,i0,a)')' Calculation status: ',nkpt_summed,' k-points to be completed' 1408 call wrtout(std_out, msg) 1409 1410 ! ============================================ 1411 ! === Begin big fat loop over transitions === 1412 ! ============================================ 1413 chi0=czero_gw; chi0_sumrule=zero 1414 1415 ! === Loop on spin to calculate trace $\chi_{up,up}+\chi_{down,down}$ === 1416 ! Only $\chi_{up,up} for AFM. 1417 do spin=1,nsppol 1418 if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE 1419 1420 ! Allocation of arrays that are private to loop 1421 if (Ep%gwcomp==1) then 1422 ABI_MALLOC(wfwfg,(nfft*nspinor**2)) 1423 end if 1424 if (Ep%gwcomp==1) then 1425 ABI_MALLOC(green_enhigh_w,(Ep%nomega)) 1426 end if 1427 if (Ep%spmeth==0) then 1428 ABI_MALLOC(green_w,(Ep%nomega)) 1429 end if 1430 if (Psps%usepaw==1) then 1431 ABI_MALLOC(Cprj2_k ,(Cryst%natom,nspinor)) 1432 call pawcprj_alloc(Cprj2_k, 0,Wfd%nlmn_atm) 1433 ABI_MALLOC(Cprj1_kmq,(Cryst%natom,nspinor)) 1434 call pawcprj_alloc(Cprj1_kmq,0,Wfd%nlmn_atm) 1435 if (Dtset%pawcross==1) then 1436 ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor)) 1437 ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor)) 1438 ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor)) 1439 ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor)) 1440 ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor)) 1441 ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor)) 1442 ABI_MALLOC(igfftepsG0f,(Ep%npwepG0)) 1443 ABI_MALLOC(tabrf_k,(nfftf_tot)) 1444 ABI_MALLOC(tabrf_kmq,(nfftf_tot)) 1445 end if 1446 end if 1447 1448 ABI_MALLOC(rhotwg,(Ep%npwepG0*dim_rtwg)) 1449 ABI_MALLOC(tabr_k,(nfft)) 1450 ABI_MALLOC(tabr_kmq,(nfft)) 1451 ABI_MALLOC(ur1_kmq_ibz,(nfft*nspinor)) 1452 ABI_MALLOC(ur2_k_ibz,(nfft*nspinor)) 1453 ABI_MALLOC(usr1_kmq,(nfft*nspinor)) 1454 ABI_MALLOC(ur2_k, (nfft*nspinor)) 1455 ABI_MALLOC(igfftepsG0,(Ep%npwepG0)) 1456 1457 ! Loop over k-points in the BZ. 1458 do ik_bz=1,Kmesh%nbz 1459 1460 if (Ep%symchi==1) then 1461 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 1462 end if 1463 1464 if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 1465 1466 write(msg,'(2(a,i4),a,i2,a,i3)')' ik= ',ik_bz,'/',Kmesh%nbz,' spin= ',spin,' done by mpi rank:',Wfd%my_rank 1467 call wrtout(std_out,msg) 1468 1469 ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz. 1470 call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt,umklp_k,isirred_k) 1471 1472 call kmesh%get_BZ_diff(kbz,qpoint,ikmq_bz,G0,nfound) 1473 if (nfound==0) then 1474 ABI_ERROR("Cannot find kbz - qpoint in Kmesh") 1475 end if 1476 1477 ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz. 1478 call kmesh%get_BZ_item(ikmq_bz,kmq_bz,ikmq_ibz,isym_kmq,itim_kmq,ph_mkmqt,umklp_kmq,isirred_kmq) 1479 1480 !BEGIN DEBUG 1481 !if (ANY(umklp_k /=0)) then 1482 ! write(msg,'(a,3i2)')" umklp_k /= 0 ",umklp_k 1483 ! ABI_ERROR(msg) 1484 !end if 1485 !if (ANY( g0 /= -umklp_kmq + umklp_k) ) then 1486 !if (ANY( g0 /= -umklp_kmq ) ) then 1487 ! write(msg,'(a,3(1x,3i2))')" g0 /= -umklp_kmq + umklp_k ",g0, umklp_kmq, umklp_k 1488 ! ABI_ERROR(msg) 1489 !end if 1490 !g0 = -umklp_k + umklp_kmq 1491 !g0 = +umklp_k - umklp_kmq 1492 !if (ANY (ABS(g0) > Ep%mg0) ) then 1493 ! write(msg,'(a,6(1x,i0))')" ABS(g0) > Ep%mg0 ",g0,Ep%mg0 1494 ! ABI_ERROR(msg) 1495 !end if 1496 !END DEBUG 1497 1498 ! Copy tables for rotated FFT points 1499 tabr_k(:) =ktabr(:,ik_bz) 1500 spinrot_k(:)=Cryst%spinrot(:,isym_k) 1501 1502 tabr_kmq(:)=ktabr(:,ikmq_bz) 1503 spinrot_kmq(:)=Cryst%spinrot(:,isym_kmq) 1504 1505 if (Dtset%pawcross==1) then 1506 tabrf_k(:) =ktabrf(:,ik_bz) 1507 tabrf_kmq(:)=ktabrf(:,ikmq_bz) 1508 end if 1509 ! 1510 ! Tables for the FFT of the oscillators. 1511 ! a) FFT index of G-G0. 1512 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 1513 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2)) 1514 call Gsph_epsG0%fft_tabs(g0,gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igfftepsG0) 1515 if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 1516 #ifdef FC_IBM 1517 ! XLF does not deserve this optimization (problem with [v67mbpt][t03]) 1518 use_padfft = 0 1519 #endif 1520 if (use_padfft==0) then 1521 ABI_FREE(gw_gbound) 1522 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft)) 1523 end if 1524 1525 if (Dtset%pawcross==1) then 1526 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 1527 call Gsph_epsG0%fft_tabs(g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f) 1528 if (ANY(gw_fftalga == [2, 4])) use_padfftf=0 1529 if (use_padfftf==0) then 1530 ABI_FREE(gboundf) 1531 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 1532 end if 1533 end if 1534 1535 nbmax=Ep%nbnds 1536 do band1=1,nbmax ! Loop over "conduction" states. 1537 if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 1538 1539 call wfd%get_ur(band1,ikmq_ibz,spin,ur1_kmq_ibz) 1540 1541 if (Psps%usepaw==1) then 1542 call wfd%get_cprj(band1,ikmq_ibz,spin,Cryst,Cprj1_kmq,sorted=.FALSE.) 1543 call paw_symcprj(ikmq_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_kmq) 1544 if (Dtset%pawcross==1) then 1545 call wfdf%paw_get_aeur(band1,ikmq_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1) 1546 end if 1547 end if 1548 1549 e_b1_kmq=qp_eig(band1,ikmq_ibz,spin) 1550 f_b1_kmq= qp_occ(band1,ikmq_ibz,spin) 1551 1552 do band2=1,nbmax ! Loop over "valence" states. 1553 !debug if (.not.luwindow.AND.dtset%ucrpa==1.AND.band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 1554 !debug& .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) CYClE 1555 !write(6,*) "ik,band1,band2",ik_bz,band1,band2 1556 if (luwindow.AND.dtset%ucrpa==1 & 1557 .AND.((qp_ebands%eig(band1,ik_ibz ,spin)-qp_ebands%fermie)<=dtset%ucrpa_window(2)) & 1558 .AND.((qp_ebands%eig(band1,ik_ibz ,spin)-qp_ebands%fermie)>=dtset%ucrpa_window(1)) & 1559 .AND.((qp_ebands%eig(band2,ikmq_ibz,spin)-qp_ebands%fermie)<=dtset%ucrpa_window(2)) & 1560 .AND.((qp_ebands%eig(band2,ikmq_ibz,spin)-qp_ebands%fermie)>=dtset%ucrpa_window(1))) CYCLE 1561 1562 if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE 1563 1564 deltaf_b1kmq_b2k=spin_fact*(f_b1_kmq-qp_occ(band2,ik_ibz,spin)) 1565 1566 if (Ep%gwcomp==0) then ! Skip negligible transitions. 1567 if (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC) CYCLE 1568 else 1569 ! When the completeness correction is used, 1570 ! we need to also consider transitions with vanishing deltaf 1571 !if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC) CYCLE 1572 ! 1573 ! Rangel This is to compute chi correctly when using the extrapolar method 1574 if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. band1<band2)) CYCLE 1575 end if 1576 1577 deltaeGW_b1kmq_b2k=e_b1_kmq-qp_eig(band2,ik_ibz,spin) 1578 1579 call wfd%get_ur(band2,ik_ibz,spin,ur2_k_ibz) 1580 1581 if (Psps%usepaw==1) then 1582 call wfd%get_cprj(band2,ik_ibz,spin,Cryst,Cprj2_k,sorted=.FALSE.) 1583 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_k) 1584 if (Dtset%pawcross==1) then 1585 call wfdf%paw_get_aeur(band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae2,ur_ae_onsite2,ur_ps_onsite2) 1586 end if 1587 end if 1588 1589 SELECT CASE (Ep%spmeth) 1590 CASE (0) 1591 ! Standard Adler-Wiser expression. 1592 ! Add the small imaginary of the Time-Ordered RF only for non-zero real omega ! FIXME What about metals? 1593 if (.not.use_tr) then 1594 ! Have to sum over all possible resonant and anti-resonant transitions. 1595 do io=1,Ep%nomega 1596 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole) 1597 end do 1598 1599 else 1600 if (Ep%gwcomp==0) then ! cannot be completely skipped in case of completeness correction 1601 if (band1<band2) CYCLE ! Here we GAIN a factor ~2 1602 end if 1603 1604 do io=1,Ep%nomega 1605 !Rangel: In metals, the intra-band transitions term does not contain the antiresonant part 1606 !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0) 1607 if (band1==band2) then 1608 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole) 1609 else 1610 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,two_poles) 1611 end if 1612 1613 if (Ep%gwcomp==1) then ! Calculate the completeness correction 1614 numerator= -spin_fact*qp_occ(band2,ik_ibz,spin) 1615 deltaeGW_enhigh_b2k = en_high-qp_eig(band2,ik_ibz,spin) 1616 1617 if (REAL(Ep%omega(io))<GW_TOL_W0) then ! Completeness correction is NOT valid for real frequencies 1618 green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2k,Ep%zcut,GW_TOL_W0,two_poles) 1619 else 1620 green_enhigh_w(io) = local_czero_gw 1621 end if 1622 ! 1623 !Rangel Correction for metals 1624 !if (deltaf_b1kmq_b2k<0.d0) then 1625 if (band1>=band2 .and. ABS(deltaf_b1kmq_b2k) > GW_TOL_DOCC ) then 1626 green_w(io)= green_w(io) - green_enhigh_w(io) 1627 else ! Disregard green_w, since it is already accounted for through the time-reversal 1628 green_w(io)= - green_enhigh_w(io) 1629 end if 1630 end if !gwcomp==1 1631 end do !io 1632 1633 if (Ep%gwcomp==1.and.band1==band2) then 1634 ! Add the "delta part" of the extrapolar method. TODO doesnt work for spinor 1635 call calc_wfwfg(tabr_k,itim_k,spinrot_k,nfft,nspinor,ngfft_gw,ur2_k_ibz,ur2_k_ibz,wfwfg) 1636 1637 if (Psps%usepaw==1) then 1638 call paw_rho_tw_g(nfft,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,gw_gfft,& 1639 & Cprj2_k,Cprj2_k,Pwij_fft,wfwfg) 1640 1641 ! Add PAW cross term 1642 if (Dtset%pawcross==1) then 1643 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 1644 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,& 1645 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k ,tabrf_k ,ph_mkt ,spinrot_k,dim_rtwg,wfwfg) 1646 end if 1647 end if 1648 1649 qzero=.FALSE. 1650 call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,& 1651 & nfft,ngfft_gw,gspfft_igfft,gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0) 1652 1653 end if 1654 end if ! use_tr 1655 1656 CASE (1, 2) 1657 ! Spectral method, WARNING time-reversal here is always assumed! 1658 if (deltaeGW_b1kmq_b2k<0) CYCLE 1659 call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1kmq_b2k,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth) 1660 END SELECT 1661 1662 ! Form rho-twiddle(r)=u^*_{b1,kmq_bz}(r) u_{b2,kbz}(r) and its FFT transform. 1663 call rho_tw_g(nspinor,Ep%npwepG0,nfft,ndat1,ngfft_gw,1,use_padfft,igfftepsG0,gw_gbound,& 1664 & ur1_kmq_ibz,itim_kmq,tabr_kmq,ph_mkmqt,spinrot_kmq,& 1665 & ur2_k_ibz, itim_k ,tabr_k ,ph_mkt ,spinrot_k,dim_rtwg,rhotwg) 1666 1667 if (Psps%usepaw==1) then 1668 ! Add PAW on-site contribution, projectors are already in the BZ. 1669 call paw_rho_tw_g(Ep%npwepG0,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,& 1670 & Cprj1_kmq,Cprj2_k,Pwij,rhotwg) 1671 1672 ! Add PAW cross term 1673 if (Dtset%pawcross==1) then 1674 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 1675 & ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,& 1676 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k ,tabrf_k ,ph_mkt ,spinrot_k,dim_rtwg,rhotwg) 1677 end if 1678 end if 1679 1680 SELECT CASE (Ep%spmeth) 1681 1682 CASE (0) ! Adler-Wiser. 1683 !debug if(dtset%ucrpa==2) then 1684 if(dtset%ucrpa>=1.and..not.luwindow) then 1685 fac=one 1686 fac2=one 1687 fac3=one 1688 m1=-1 1689 m2=-1 1690 call flush_unit(std_out) 1691 if(dtset%ucrpa<=2) then 1692 if ( band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 1693 & .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 1694 if (dtset%plowan_compute >=10) then 1695 band1c=band1-wan%bandi_wan+1 1696 band2c=band2-wan%bandi_wan+1 1697 do iat1=1, wan%natom_wan 1698 do iat2=1, wan%natom_wan 1699 do ispinor1=1,wan%nspinor 1700 do ispinor2=1,wan%nspinor 1701 do il1=1,wan%nbl_atom_wan(iat1) 1702 do il2=1,wan%nbl_atom_wan(iat2) 1703 do m1=1,2*wan%latom_wan(iat1)%lcalc(il1)+1 1704 do m2=1,2*wan%latom_wan(iat2)%lcalc(il2)+1 1705 fac=fac - real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1)*& 1706 &conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*& 1707 &wan%psichi(ikmq_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)*& 1708 &conjg(wan%psichi(ikmq_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2))) 1709 enddo !m2 1710 enddo !m1 1711 enddo !il2 1712 enddo !il1 1713 enddo !ispinor2 1714 enddo !isspinor1 1715 enddo !iat2 1716 enddo !iat1 1717 else !plowan_compute>=10 1718 do iat=1, cryst%nattyp(itypatcor) 1719 do ispinor1=1,nspinor 1720 do ispinor2=1,nspinor 1721 do m1=1,2*lcor+1 1722 do m2=1,2*lcor+1 1723 fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 1724 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* & 1725 & coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2)*& 1726 & conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2))) 1727 enddo !m2 1728 enddo !m1 1729 enddo !ispinor2 1730 enddo !ispinor1 1731 enddo !iat 1732 endif !plowan_compute>=10 1733 if(dtset%ucrpa==1) fac=zero 1734 endif 1735 else if (dtset%ucrpa>=3) then 1736 if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then 1737 do iat=1, cryst%nattyp(itypatcor) 1738 do ispinor1=1,nspinor 1739 do m1=1,2*lcor+1 1740 fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 1741 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))) 1742 enddo 1743 enddo 1744 enddo 1745 if(dtset%ucrpa==4) fac2=zero 1746 endif 1747 if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 1748 do iat=1, cryst%nattyp(itypatcor) 1749 do ispinor1=1,nspinor 1750 do m1=1,2*lcor+1 1751 fac3=fac3-real(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1)*& 1752 & conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1))) 1753 enddo 1754 enddo 1755 enddo 1756 if(dtset%ucrpa==4) fac3=zero 1757 endif 1758 fac=real(fac2*fac3) 1759 endif 1760 1761 ! if(dtset%prtvol>=10) write(6,'(6i3,e15.5,a)') ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q=/0" 1762 ! if(dtset%prtvol>=10.and.abs(fac-one)>0.00001) & 1763 !& write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0" 1764 ! if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0" 1765 green_w=green_w*fac 1766 endif 1767 1768 call assemblychi0_sym(is_metallic,ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0) 1769 1770 CASE (1, 2) 1771 ! Spectral method (not yet adapted for nspinor=2) 1772 call assemblychi0sf(ik_bz,Ep%symchi,Ltg_q,Ep%npwepG0,Ep%npwe,rhotwg,Gsph_epsG0,& 1773 & deltaf_b1kmq_b2k,my_wl,iomegal,wl,my_wr,iomegar,wr,Ep%nomegasf,sf_chi0) 1774 1775 CASE DEFAULT 1776 ABI_BUG("Wrong spmeth") 1777 END SELECT 1778 1779 ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. [[cite:Mahan2000]] 1780 ! TODO Does not work with spinor 1781 factor=spin_fact*qp_occ(band2,ik_ibz,spin) 1782 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1kmq_b2k,& 1783 & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 1784 1785 ! Include also the completeness correction in the sum rule 1786 if (Ep%gwcomp==1) then 1787 factor=-spin_fact*qp_occ(band2,ik_ibz,spin) 1788 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2k,& 1789 & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 1790 if (band1==Ep%nbnds) then 1791 chi0_sumrule(:)=chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2k 1792 end if 1793 end if 1794 1795 end do !band2 1796 end do !band1 1797 1798 ABI_FREE(gw_gbound) 1799 if (Dtset%pawcross==1) then 1800 ABI_FREE(gboundf) 1801 end if 1802 end do !ik_bz 1803 1804 ! Deallocation of arrays private to the spin loop. 1805 ABI_FREE(igfftepsG0) 1806 ABI_FREE(ur1_kmq_ibz) 1807 ABI_FREE(ur2_k_ibz) 1808 ABI_FREE(usr1_kmq) 1809 ABI_FREE(ur2_k) 1810 ABI_FREE(rhotwg) 1811 ABI_FREE(tabr_k) 1812 ABI_FREE(tabr_kmq) 1813 1814 ABI_SFREE(green_w) 1815 ABI_SFREE(wfwfg) 1816 ABI_SFREE(green_enhigh_w) 1817 if (Psps%usepaw==1) then 1818 call pawcprj_free(Cprj2_k) 1819 ABI_FREE(Cprj2_k) 1820 call pawcprj_free(Cprj1_kmq) 1821 ABI_FREE(Cprj1_kmq) 1822 if (Dtset%pawcross==1) then 1823 ABI_FREE(ur_ae1) 1824 ABI_FREE(ur_ae_onsite1) 1825 ABI_FREE(ur_ps_onsite1) 1826 ABI_FREE(ur_ae2) 1827 ABI_FREE(ur_ae_onsite2) 1828 ABI_FREE(ur_ps_onsite2) 1829 ABI_FREE(tabrf_k) 1830 ABI_FREE(tabrf_kmq) 1831 ABI_FREE(gboundf) 1832 ABI_FREE(igfftepsG0f) 1833 end if 1834 end if 1835 end do !spin 1836 1837 ! After big loop over transitions, now MPI 1838 ! Master took care of the contribution in case of metallic|spin polarized systems. 1839 SELECT CASE (Ep%spmeth) 1840 CASE (0) 1841 ! Adler-Wiser 1842 ! Collective sum of the contributions of each node. 1843 ! Looping on frequencies to avoid problems with the size of the MPI packet 1844 do io=1,Ep%nomega 1845 call xmpi_sum(chi0(:,:,io),comm,ierr) 1846 end do 1847 1848 CASE (1, 2) 1849 ! Spectral method. 1850 call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth) 1851 1852 ! Deallocate here before xmpi_sum 1853 ABI_SFREE(sf_chi0) 1854 1855 ! Collective sum of the contributions. 1856 ! Looping over frequencies to avoid problems with the size of the MPI packet 1857 do io=1,Ep%nomega 1858 call xmpi_sum(chi0(:,:,io),comm,ierr) 1859 end do 1860 1861 CASE DEFAULT 1862 ABI_BUG("Wrong spmeth") 1863 END SELECT 1864 1865 ! Divide by the volume 1866 !$OMP PARALLEL WORKSHARE 1867 chi0=chi0*weight/Cryst%ucvol 1868 !$OMP END PARALLEL WORKSHARE 1869 1870 ! === Collect the sum rule === 1871 ! * The pi factor comes from Im[1/(x-ieta)] = pi delta(x) 1872 call xmpi_sum(chi0_sumrule,comm,ierr) 1873 chi0_sumrule=chi0_sumrule*pi*weight/Cryst%ucvol 1874 ! 1875 ! ************************************************* 1876 ! **** Now each node has chi0(q,G,Gp,Ep%omega) **** 1877 ! ************************************************* 1878 1879 ! Impose Hermiticity (valid only for zero or purely imaginary frequencies) 1880 ! MG what about metals, where we have poles around zero? 1881 ! FB because of the intraband term, chi0 is never hermitian in case of metals 1882 ! FIXME: as of today, hermitianity is also enforced for metallic systems 1883 !if (.not. is_metallic) then 1884 do io=1,Ep%nomega 1885 if (ABS(REAL(Ep%omega(io))) <0.00001) then 1886 do ig2=1,Ep%npwe 1887 do ig1=1,ig2-1 1888 chi0(ig2,ig1,io) = GWPC_CONJG(chi0(ig1,ig2,io)) 1889 end do 1890 end do 1891 end if 1892 end do 1893 !endif 1894 1895 ! === Symmetrize chi0 in case of AFM system === 1896 ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$. 1897 ! Works only in case of magnetic group Shubnikov type IV. 1898 if (Cryst%use_antiferro) then 1899 call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0=chi0) 1900 end if 1901 1902 ! ===================== 1903 ! ==== Free memory ==== 1904 ! ===================== 1905 ABI_FREE(bbp_ks_distrb) 1906 1907 ABI_SFREE(gw_gfft) 1908 ABI_SFREE(kkweight) 1909 ABI_SFREE(omegasf) 1910 ABI_SFREE(gspfft_igfft) 1911 1912 call Gsph_FFT%free() 1913 1914 ! deallocation for PAW. 1915 if (Psps%usepaw==1) then 1916 call pawpwij_free(Pwij) 1917 ABI_FREE(Pwij) 1918 if (allocated(Pwij_fft)) then 1919 call pawpwij_free(Pwij_fft) 1920 ABI_FREE(Pwij_fft) 1921 end if 1922 end if 1923 1924 if(dtset%ucrpa>=1 .AND. dtset%plowan_compute<10) then 1925 ABI_FREE(coeffW_BZ) 1926 endif 1927 1928 call timab(331,2,tsec) 1929 call cwtime(cpu_time,wall_time,gflops,"stop") 1930 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 1931 1932 DBG_EXIT("COLL") 1933 1934 end subroutine cchi0
ABINIT/cchi0q0 [ Functions ]
NAME
cchi0q0
FUNCTION
Calculate chi0 in the limit q --> 0
INPUTS
use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry) Dtset <type(dataset_type)>=all input variables in this dataset Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp) %ng=number of G vectors %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G in the array gvec %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0 %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$). %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1) Ep%inclvkb=flag to include (or not) the grad of Vkb Ltg_q= little group datatype nbvw=number of bands in the arrays wfrv,wfgv Kmesh<kmesh_t> The k-point mesh %kbz(3,nbz)=k-point coordinates, full Brillouin zone %tab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding irreducible point (kIBZ), where kBZ= (IS) kIBZ and I is either the inversion or the identity %tabi(nbzx)= for each point in the BZ defines whether inversion has to be considered in the relation kBZ=(IS) kIBZ (1 => only S; -1 => -S) %tabo(nbzx)= the symmetry operation S that takes kIBZ to each kBZ %tabp(nbzx)= phase factor associated to tnons e^{-i 2 \pi k\cdot R{^-1}t} ktabr(nfftot_gw,Kmesh%nbz) index of R^-(r-t) in the FFT array, where k_BZ = (IS) k_IBZ and S = \transpose R^{-1} Ep%nbnds=number of bands ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths. Ep%nomega=number of frequencies Cryst<crystal_t>= data type gathering info on symmetries and unit cell %natom=number of atoms %nsym=number of symmetry operations %symrec(3,3,nsym)=symmetry operations in reciprocal space %typat(natom)=type of each atom %xred(3,natom)=reduced coordinated of atoms %rprimd(3,3)=dimensional primitive translations in real space (bohr) %timrev=2 if time-reversal symmetry can be used, 1 otherwise Ep%npwe=number of planewaves for sigma exchange (input variable) nfftot_gw=Total number of points in the GW FFT grid Ep%omega(Ep%nomega)=frequencies Psps <type(pseudopotential_type)>=variables related to pseudopotentials %mpsang=1+maximum angular momentum for nonlocal pseudopotential Pawang<pawang_type> angular mesh discretization and related data: Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels qp_ebands<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities) %mband=MAX number of bands over k-points and spin (==Ep%nbnds) %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes ks_ebands<ebands_t>=KS energies and occupations. %eig(mband,nkpt,nsppol)=KS energies Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave. Wfd<wfdgw_t>=Object used to access the wavefunctions
OUTPUT
chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq, and frequencies defined by Ep%omega chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)= Lower wings chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)= Upper wings chi0_head(3,3,Ep%nomega)=Head of chi0.
NOTES
*) The terms "head", "wings" and "body" of chi(G,Gp) refer to G=Gp=0, either G or Gp=0, and neither=0 respectively *) Symmetry conventions: 1) symmetry in real space is defined as: R_t f(r) = f(R^-1(r-t)) 2) S=\transpose R^-1 3) kbz=S kibz The wavefunctions for the k-point in the BZ are (assuming nondegenerate states): u(G,b, Sk) = u ( S^-1G,b,k)* e^{-i(Sk+G)*t) u(G,b,-Sk) = u*(-S^-1G,b,k)* e^{ i(Sk-G)*t) u(r,b, Sk) = u (R^-1(r-t),b,k) e^{-iSk*t} u(r,b,-Sk) = u*(R^-1(r-t),b,k) e^{ iSK*t} The gradient of Vnl(K,Kp) for the k-point in the BZ should be: gradvnl(SG,SGp,Sk)=S gradvnl(G,Gp,kibz)
TODO
Check npwepG0 before activating umklapps
SOURCE
167 subroutine cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,qp_ebands,ks_ebands,Gsph_epsG0,& 168 Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,& 169 nfftot_gw,ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q,chi0_sumrule,Wfd,Wfdf,wan) 170 171 !Arguments ------------------------------------ 172 !scalars 173 integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot 174 logical,intent(in) :: use_tr 175 type(ebands_t),target,intent(in) :: qp_ebands,ks_ebands 176 type(crystal_t),intent(in) :: Cryst 177 type(Dataset_type),intent(in) :: Dtset 178 type(littlegroup_t),intent(in) :: Ltg_q 179 type(em1params_t),intent(in) :: Ep 180 type(kmesh_t),intent(in) :: Kmesh 181 type(gsphere_t),intent(in) :: Gsph_epsG0 182 type(Pseudopotential_type),intent(in) :: Psps 183 type(Pawang_type),intent(in) :: Pawang 184 type(wfdgw_t),target,intent(inout) :: Wfd,Wfdf 185 !arrays 186 integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz) 187 integer,intent(in) :: ngfft_gw(18),ngfftf(18) 188 real(dp),intent(out) :: chi0_sumrule(Ep%npwe) 189 complex(gwpc),intent(out) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega) 190 complex(dpc),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3) 191 complex(dpc),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3) 192 complex(dpc),intent(out) :: chi0_head(3,3,Ep%nomega) 193 type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw) 194 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 195 type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw) 196 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 197 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 198 type(plowannier_type),intent(inout) :: wan 199 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom) 200 201 !Local variables ------------------------------ 202 !scalars 203 integer,parameter :: tim_fourdp=1,enough=10,two_poles=2,one_pole=1,ndat1=1 204 integer :: bandinf,bandsup,lcor,nspinor,npw_k,istwf_k,mband,nfft,band1c,band2c 205 integer :: band1,band2,iat1,iat2,iat,ig,ig1,ig2,itim_k,ik_bz,ik_ibz,io,iqlwl,ispinor1,ispinor2,isym_k,il1,il2 206 integer :: itypatcor,m1,m2,nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,use_padfftf,mgfftf 207 integer :: my_nbbp,my_nbbpks,spin,nsppol,iq,nq 208 integer :: comm,ierr,my_wl,my_wr,iomegal,iomegar,gw_mgfft,dummy 209 real(dp) :: cpu_time,wall_time,gflops 210 real(dp) :: fac,fac1,fac2,fac3,fac4,spin_fact,deltaf_b1b2,weight,factor 211 real(dp) :: max_rest,min_rest,my_max_rest,my_min_rest, qlen 212 real(dp) :: en_high,deltaeGW_enhigh_b2,wl,wr,numerator,deltaeGW_b1b2,gw_gsq,memreq 213 complex(dpc) :: deltaeKS_b1b2 214 logical :: qzero,luwindow,is_metallic 215 character(len=500) :: msg_tmp,msg,allup 216 type(gsphere_t) :: Gsph_FFT 217 type(wave_t),pointer :: wave1, wave2 218 !arrays 219 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 220 integer :: ucrpa_bands(2), got(Wfd%nproc) 221 integer :: wtk_ltg(Kmesh%nbz) 222 integer,allocatable :: tabr_k(:),tabrf_k(:), igffteps0(:),gspfft_igfft(:),igfftepsG0f(:) 223 integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:), bbp_ks_distrb(:,:,:,:) 224 real(dp) :: kbz(3),spinrot_kbz(4),q0(3) 225 real(dp),ABI_CONTIGUOUS pointer :: ks_eig(:,:,:),qp_eig(:,:,:),qp_occ(:,:,:) 226 real(dp),allocatable :: omegasf(:), qdirs(:,:) 227 complex(gwpc) :: rhotwx(3,Wfd%nspinor**2) 228 complex(gwpc),allocatable :: rhotwg(:) 229 complex(dpc),allocatable :: green_w(:),green_enhigh_w(:) 230 complex(dpc),allocatable :: sf_lwing(:,:,:),sf_uwing(:,:,:),sf_head(:,:,:) 231 complex(dpc) :: chq(3), wng(3) 232 complex(dpc) :: ph_mkt 233 complex(dpc),allocatable :: kkweight(:,:) 234 complex(gwpc),allocatable :: ur1_kibz(:),ur2_kibz(:), usr1_k(:),ur2_k(:), wfwfg(:), sf_chi0(:,:,:) 235 complex(gwpc),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:) 236 complex(gwpc),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:) 237 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:),ug2(:) 238 complex(dpc), allocatable :: coeffW_BZ(:,:,:,:,:,:), head_qvals(:) 239 logical :: gradk_not_done(Kmesh%nibz) 240 logical,allocatable :: bbp_mask(:,:) 241 type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj2_bz(:,:), Cprj1_ibz(:,:),Cprj2_ibz(:,:) 242 type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:) 243 type(pawhur_t),allocatable :: Hur(:) 244 type(vkbr_t),allocatable :: vkbr(:) 245 246 !************************************************************************ 247 248 DBG_ENTER("COLL") 249 250 call cwtime(cpu_time, wall_time, gflops, "start") 251 252 ! Change FFT mesh if needed 253 if (ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3))) call wfd%change_ngfft(Cryst,Psps,ngfft_gw) 254 255 gw_mgfft = MAXVAL(ngfft_gw(1:3)); gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10) 256 if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3)) 257 258 ! Copy important variables. 259 comm = Wfd%comm; nsppol = Wfd%nsppol; nspinor = Wfd%nspinor; mband = Wfd%mband; nfft = Wfd%nfft 260 ABI_CHECK(Wfd%nfftot == nfftot_gw, "Wrong nfftot_gw") 261 dim_rtwg = 1 !; if (nspinor==2) dim_rtwg=2 ! Can reduce size depending on Ep%nI and Ep%nj 262 263 is_metallic = ebands_has_metal_scheme(qp_ebands) 264 ucrpa_bands(1)=dtset%ucrpa_bands(1) 265 ucrpa_bands(2)=dtset%ucrpa_bands(2) 266 luwindow=.false. 267 if (abs(dtset%ucrpa_window(1)+1_dp)>tol8.or.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) luwindow=.true. 268 269 ! For cRPA calculation of U: read forlb.ovlp 270 if(dtset%ucrpa>=1 .AND. dtset%plowan_compute<10) then 271 call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,& 272 nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands) 273 endif 274 275 ks_eig => ks_ebands%eig 276 qp_eig => qp_ebands%eig; qp_occ => qp_ebands%occ 277 278 chi0_lwing = czero; chi0_uwing = czero; chi0_head = czero 279 280 if (Psps%usepaw == 0) then 281 if (Ep%inclvkb /= 0) then 282 ! Include the term <n,k|[Vnl,iqr]|n"k>' for q -> 0. 283 ABI_CHECK(nspinor == 1, "nspinor with inclvkb not coded") 284 else 285 ABI_WARNING('Neglecting <n,k|[Vnl,iqr]|m,k>') 286 end if 287 288 else 289 ! For PAW+DFT+U, precalculate <\phi_i|[Hu,r]|phi_j\> 290 ABI_MALLOC(HUr, (Cryst%natom)) 291 if (Dtset%usepawu /= 0) then 292 call pawhur_init(hur,nsppol,Dtset%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij) 293 end if 294 end if 295 296 ! Initialize the completeness correction. 297 ABI_MALLOC(green_enhigh_w, (Ep%nomega)) 298 green_enhigh_w = czero 299 300 if (Ep%gwcomp == 1) then 301 en_high = MAXVAL(qp_eig(Ep%nbnds,:,:))+Ep%gwencomp 302 write(msg,'(a,f8.2,a)')' Using completeness correction with energy ',en_high*Ha_eV,' [eV] ' 303 call wrtout(std_out, msg) 304 ABI_MALLOC(wfwfg,(nfft*nspinor**2)) 305 306 ! Init the largest G-sphere contained in the FFT box for the wavefunctions. 307 call Gsph_FFT%in_fftbox(Cryst,Wfd%ngfft) 308 call Gsph_FFT%print(unit=std_out, prtvol=10) 309 310 ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng)) 311 ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2)) 312 313 ! Mapping between G-sphere and FFT box. 314 call Gsph_FFT%fft_tabs([0, 0, 0],Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft) 315 ABI_FREE(dummy_gbound) 316 317 if (Psps%usepaw==1) then 318 ! Prepare the onsite contributions on the GW FFT mesh. 319 ABI_MALLOC(gw_gfft,(3,nfft)) 320 q0=zero 321 call get_gftt(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! The set of plane waves in the FFT Box. 322 ABI_MALLOC(Pwij_fft,(Psps%ntypat)) 323 call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 324 end if 325 end if 326 327 ! Setup weight (2 for spin unpolarized systems, 1 for polarized). 328 ! spin_fact is used to normalize the occupation factors to one. 329 ! Consider also the AFM case. 330 select case (nsppol) 331 case (1) 332 weight = two / Kmesh%nbz; spin_fact = half 333 if (Wfd%nspden == 2) then 334 weight = one / Kmesh%nbz; spin_fact = half 335 end if 336 if (nspinor == 2) then 337 weight = one / Kmesh%nbz; spin_fact = one 338 end if 339 case (2) 340 weight = one / Kmesh%nbz; spin_fact = one 341 case default 342 ABI_BUG(sjoin("Wrong nsppol:", itoa(nsppol))) 343 end select 344 345 ! k-weights for points in the IBZ_q 346 wtk_ltg(:) = 1 347 if (Ep%symchi == 1) then 348 do ik_bz=1,Ltg_q%nbz 349 wtk_ltg(ik_bz) = 0 350 if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only k in IBZ_q 351 wtk_ltg(ik_bz) = sum(Ltg_q%wtksym(:,:,ik_bz)) 352 end do 353 end if 354 355 write(msg,'(a,i3,a)')' Q-points for long wave-length limit. # ',Ep%nqlwl,ch10 356 do iqlwl=1,Ep%nqlwl 357 write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',Ep%qlwl(:,iqlwl),ch10 358 msg=TRIM(msg)//msg_tmp 359 end do 360 call wrtout(std_out, msg) 361 362 write(msg,'(a,i2,2a,i2)')& 363 ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,& 364 ' Using symmetries to sum only over the IBZ_q = ',Ep%symchi 365 call wrtout(std_out, msg) 366 367 if (use_tr) then 368 call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.') 369 else 370 call wrtout(std_out,' Using slow algorithm without time reversal symmetry.') 371 end if 372 373 ! Evaluate oscillator matrix elements btw partial waves. Note q=Gamma 374 if (Psps%usepaw == 1) then 375 ABI_MALLOC(Pwij,(Psps%ntypat)) 376 call pawpwij_init(Pwij,Ep%npwepG0, [zero, zero, zero],Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 377 378 ABI_MALLOC(Cprj1_bz,(Cryst%natom,nspinor)) 379 call pawcprj_alloc(Cprj1_bz,0,Wfd%nlmn_atm) 380 ABI_MALLOC(Cprj2_bz,(Cryst%natom,nspinor)) 381 call pawcprj_alloc(Cprj2_bz,0,Wfd%nlmn_atm) 382 383 ABI_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor)) 384 call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm) 385 ABI_MALLOC(Cprj2_ibz,(Cryst%natom,nspinor)) 386 call pawcprj_alloc(Cprj2_ibz,0,Wfd%nlmn_atm) 387 if (Dtset%pawcross==1) then 388 ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor)) 389 ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor)) 390 ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor)) 391 ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor)) 392 ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor)) 393 ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor)) 394 ABI_MALLOC(igfftepsG0f,(Ep%npwepG0)) 395 ABI_MALLOC(tabrf_k,(nfftf_tot)) 396 end if 397 end if 398 399 ABI_MALLOC(rhotwg,(Ep%npwe*dim_rtwg)) 400 ABI_MALLOC(tabr_k,(nfft)) 401 ABI_MALLOC(ur1_kibz,(nfft*nspinor)) 402 ABI_MALLOC(ur2_kibz,(nfft*nspinor)) 403 404 ABI_MALLOC(usr1_k,(nfft*nspinor)) 405 ABI_MALLOC(ur2_k,(nfft*nspinor)) 406 ! 407 ! Tables for the FFT of the oscillators. 408 ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere). 409 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 410 411 ABI_MALLOC(igffteps0,(Gsph_epsG0%ng)) 412 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2)) 413 call Gsph_epsG0%fft_tabs([0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0) 414 if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 415 #ifdef FC_IBM 416 ! XLF does not deserve this optimization (problem with [v67mbpt][t03]) 417 use_padfft = 0 418 #endif 419 if (use_padfft==0) then 420 ABI_FREE(gw_gbound) 421 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft)) 422 end if 423 if (Dtset%pawcross==1) then 424 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 425 call Gsph_epsG0%fft_tabs((/0,0,0/),mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f) 426 if ( ANY(gw_fftalga == (/2,4/)) ) use_padfftf=0 427 if (use_padfftf==0) then 428 ABI_FREE(gboundf) 429 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 430 end if 431 end if 432 433 ! TODO this table can be calculated for each k-point 434 my_nbbpks=0; allup="All"; got=0 435 ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol)) 436 call wrtout(std_out, sjoin(' Memory needed for bbp_ks_distrb: ', & 437 ftoa(four*mband**2*Kmesh%nbz*nsppol*b2Mb, fmt="f8.1"), ' [Mb] <<< MEM')) 438 439 ABI_MALLOC(bbp_mask, (mband, mband)) 440 441 do spin=1,nsppol 442 do ik_bz=1,Kmesh%nbz 443 444 if (Ep%symchi==1) then 445 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 446 end if 447 ik_ibz=Kmesh%tab(ik_bz) 448 449 call chi0_bbp_mask(ik_ibz, ik_ibz, spin, spin_fact, use_tr, & 450 ep%gwcomp, ep%spmeth, ep%nbnds, mband, qp_ebands, bbp_mask) 451 452 call wfd%distribute_bbp(ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got=got,bbp_mask=bbp_mask) 453 my_nbbpks = my_nbbpks + my_nbbp 454 end do 455 end do 456 457 ABI_FREE(bbp_mask) 458 459 write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0q0." 460 call wrtout(std_out, msg) 461 462 SELECT CASE (Ep%spmeth) 463 CASE (0) 464 call wrtout(std_out,' Calculating chi0(q=(0,0,0),omega,G,G")') 465 ABI_MALLOC(green_w, (Ep%nomega)) 466 467 CASE (1, 2) 468 call wrtout(std_out,' Calculating Im chi0(q=(0,0,0),omega,G,G")') 469 ! 470 ! === Find max and min resonant transitions for this q, report values for this processor === 471 call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,& 472 max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_eig,qp_occ, [zero,zero,zero], bbp_ks_distrb) 473 474 ! === Calculate frequency dependent weights for Kramers Kronig transform === 475 ABI_MALLOC(omegasf, (Ep%nomegasf)) 476 ABI_MALLOC(kkweight, (Ep%nomegasf,Ep%nomega)) 477 !my_wl=1; my_wr=Ep%nomegasf 478 call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,& 479 0,Ep%zcut,zero,my_wl,my_wr,kkweight) 480 481 ABI_CHECK(use_tr, 'Hilbert transform requires time-reversal') 482 483 ! allocate head and wings of the spectral function. 484 ABI_MALLOC(sf_head,(3,3,my_wl:my_wr)) 485 ABI_MALLOC(sf_lwing,(Ep%npwe,my_wl:my_wr,3)) 486 ABI_MALLOC(sf_uwing,(Ep%npwe,my_wl:my_wr,3)) 487 sf_head=czero; sf_lwing=czero; sf_uwing=czero 488 489 memreq = two*gwpc*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb 490 write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwpc*Ep%npwe**2*b2Mb,' [Mb]' 491 call wrtout(std_out, msg) 492 write(msg,'(a,f10.4,a)')' memory required by sf_chi0q0: ',memreq,' [Gb]' 493 call wrtout(std_out, msg) 494 if (memreq > two) then 495 ABI_WARNING(' Memory required for sf_chi0q0 is larger than 2.0 Gb!') 496 end if 497 ABI_MALLOC_OR_DIE(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr) 498 sf_chi0=czero_gw 499 500 CASE DEFAULT 501 ABI_BUG("Wrong spmeth") 502 END SELECT 503 504 nkpt_summed = Kmesh%nbz 505 if (Ep%symchi /= 0) then 506 nkpt_summed = Ltg_q%nibz_ltg 507 call Ltg_q%print(std_out, Dtset%prtvol) 508 end if 509 call wrtout(std_out, sjoin(' Calculation status: ', itoa(nkpt_summed), ' k-points to be completed')) 510 511 ABI_MALLOC(vkbr, (Kmesh%nibz)) 512 gradk_not_done = .TRUE. 513 514 ! ============================================ 515 ! === Begin big fat loop over transitions ==== 516 ! ============================================ 517 chi0 = czero_gw; chi0_sumrule = zero 518 519 ! Loop on spin to calculate $\chi_{\up,\up} + \chi_{\down,\down}$ 520 do spin=1,nsppol 521 if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE 522 523 ! Loop over k-points in the BZ. 524 do ik_bz=1,Kmesh%nbz 525 if (Ep%symchi == 1) then 526 if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only IBZ_q 527 end if 528 529 if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 530 531 write(msg,'(2(a,i4),a,i2,a,i3)')' ik= ',ik_bz,'/',Kmesh%nbz,' spin=',spin,' done by mpi rank:',Wfd%my_rank 532 call wrtout(std_out, msg) 533 534 ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz. 535 call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt) 536 tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points 537 spinrot_kbz(:)=Cryst%spinrot(:,isym_k) 538 if (Dtset%pawcross==1) tabrf_k(:) = ktabrf(:,ik_bz) 539 540 istwf_k = Wfd%istwfk(ik_ibz) 541 npw_k = Wfd%npwarr(ik_ibz) 542 kg_k => Wfd%Kdata(ik_ibz)%kg_k 543 544 if (psps%usepaw == 0 .and. Ep%inclvkb /= 0 .and. gradk_not_done(ik_ibz)) then 545 ! Include term <n,k|[Vnl,iqr]|n"k>' for q -> 0. 546 call vkbr_init(vkbr(ik_ibz), Cryst, Psps, Ep%inclvkb, istwf_k, npw_k, Kmesh%ibz(:,ik_ibz), kg_k) 547 gradk_not_done(ik_ibz) = .FALSE. 548 end if 549 550 ! Loop over "conduction" states. 551 do band1=1,Ep%nbnds 552 if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE 553 554 ABI_CHECK(wfd%get_wave_ptr(band1, ik_ibz, spin, wave1, msg) == 0, msg) 555 ug1 => wave1%ug 556 call wfd%get_ur(band1,ik_ibz,spin,ur1_kibz) 557 558 if (Psps%usepaw==1) then 559 call wfd%get_cprj(band1,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.) 560 call pawcprj_copy(Cprj1_ibz,Cprj1_bz) 561 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz) 562 if (Dtset%pawcross==1) then 563 call wfdf%paw_get_aeur(band1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1) 564 end if 565 end if 566 567 ! Loop over "valence" states. 568 do band2=1,Ep%nbnds 569 570 ! write(std_out,*) "ik,band1,band2",ik_bz,band1,band2 571 ! ----------------- cRPA for U 572 !debug if (.not.luwindow.AND.dtset%ucrpa==1.AND.band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 573 !debug& .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) CYCLE 574 if (luwindow.AND.dtset%ucrpa==1 & 575 .AND.((ks_ebands%eig(band1,ik_ibz,spin)-ks_ebands%fermie)<=dtset%ucrpa_window(2)) & 576 .AND.((ks_ebands%eig(band1,ik_ibz,spin)-ks_ebands%fermie)>=dtset%ucrpa_window(1)) & 577 .AND.((ks_ebands%eig(band2,ik_ibz,spin)-ks_ebands%fermie)<=dtset%ucrpa_window(2)) & 578 .AND.((ks_ebands%eig(band2,ik_ibz,spin)-ks_ebands%fermie)>=dtset%ucrpa_window(1))) CYCLE 579 !----------------- cRPA for U 580 581 if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE 582 583 deltaeKS_b1b2 = ks_eig(band1, ik_ibz, spin) - ks_eig(band2, ik_ibz, spin) 584 deltaf_b1b2 = spin_fact * (qp_occ(band1, ik_ibz, spin) - qp_occ(band2, ik_ibz, spin)) 585 deltaeGW_b1b2 = qp_eig(band1, ik_ibz, spin) - qp_eig(band2, ik_ibz, spin) 586 587 if (Ep%gwcomp == 0) then 588 ! Skip negligible transitions. 589 if (abs(deltaf_b1b2) < GW_TOL_DOCC) CYCLE 590 else 591 ! when the completeness trick is used, we need to also consider transitions with vanishing deltaf 592 ! Rangel Correction for metals 593 if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. ( ABS(deltaf_b1b2)< GW_TOL_DOCC .or. band1 < band2)) CYCLE 594 end if 595 596 ABI_CHECK(wfd%get_wave_ptr(band2, ik_ibz, spin, wave2, msg) == 0, msg) 597 ug2 => wave2%ug 598 call wfd%get_ur(band2,ik_ibz,spin,ur2_kibz) 599 600 if (Psps%usepaw==1) then 601 call wfd%get_cprj(band2,ik_ibz,spin,Cryst,Cprj2_ibz,sorted=.FALSE.) 602 call pawcprj_copy(Cprj2_ibz,Cprj2_bz) 603 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_bz) 604 if (Dtset%pawcross==1) then 605 call wfdf%paw_get_aeur(band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 606 ur_ae2,ur_ae_onsite2,ur_ps_onsite2) 607 end if 608 end if 609 610 SELECT CASE (Ep%spmeth) 611 CASE (0) 612 ! Adler-Wiser expression. 613 ! Add small imaginary of the Time-Ordered response function but only for non-zero real omega 614 ! FIXME What about metals? 615 616 if (.not. use_tr) then 617 ! Adler-Wiser without time-reversal. 618 do io=1,Ep%nomega 619 green_w(io) = g0g0w(Ep%omega(io), deltaf_b1b2, deltaeGW_b1b2, Ep%zcut, GW_TOL_W0, one_pole) 620 end do 621 622 else 623 if (Ep%gwcomp == 0) then ! cannot be completely skipped in case of completeness correction 624 if (band1 < band2) CYCLE ! Here we GAIN a factor ~2 625 end if 626 627 do io=1,Ep%nomega 628 ! Rangel: In metals, the intra-band transitions term does not contain the antiresonant part 629 ! if(abs(deltaeGW_b1b2)>GW_TOL_W0) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0) 630 if (band1 == band2) green_w(io) = g0g0w(Ep%omega(io), deltaf_b1b2, deltaeGW_b1b2, Ep%zcut, GW_TOL_W0, one_pole) 631 if (band1 /= band2) green_w(io) = g0g0w(Ep%omega(io), deltaf_b1b2, deltaeGW_b1b2, Ep%zcut, GW_TOL_W0, two_poles) 632 633 if (Ep%gwcomp == 1) then 634 ! Calculate the completeness correction 635 numerator= -spin_fact * qp_occ(band2,ik_ibz,spin) 636 deltaeGW_enhigh_b2 = en_high - qp_eig(band2,ik_ibz,spin) 637 ! Completeness correction is NOT valid for real frequencies 638 if (REAL(Ep%omega(io)) < GW_TOL_W0) then 639 green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2,Ep%zcut,GW_TOL_W0,two_poles) 640 else 641 green_enhigh_w(io) = czero_gw 642 endif 643 ! 644 ! Rangel Correction for metals 645 if (band1 >= band2 .and. abs(deltaf_b1b2) > GW_TOL_DOCC) then 646 green_w(io)= green_w(io) - green_enhigh_w(io) 647 else 648 ! Disregard green_w, since it is already accounted for through the time-reversal 649 green_w(io)= - green_enhigh_w(io) 650 end if 651 end if !gwcomp == 1 652 end do !io 653 654 if (Ep%gwcomp == 1 .and. band1 == band2) then 655 ! Add the "delta part", symmetrization is done inside the routine. 656 call calc_wfwfg(tabr_k,itim_k,spinrot_kbz,nfft,nspinor,ngfft_gw,ur2_kibz,ur2_kibz,wfwfg) 657 658 if (Psps%usepaw==1) then 659 call paw_rho_tw_g(nfft,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,gw_gfft,& 660 Cprj2_bz,Cprj2_bz,Pwij_fft,wfwfg) 661 662 ! Add PAW cross term 663 if (Dtset%pawcross==1) then 664 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 665 ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 666 ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 667 dim_rtwg,wfwfg) 668 end if 669 end if 670 671 qzero = .TRUE. 672 call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,& 673 nfft,ngfft_gw,gspfft_igfft,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0) 674 end if 675 end if ! use_tr 676 677 CASE (1, 2) 678 ! Spectral method, here time-reversal is always assumed. 679 if (deltaeGW_b1b2 < 0) CYCLE 680 call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1b2,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth) 681 END SELECT 682 683 ! FFT of u^*_{b1,k}(r) u_{b2,k}(r) and (q,G=0) limit using small q and k.p perturbation theory 684 call rho_tw_g(nspinor,Ep%npwe,nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,& 685 ur1_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 686 ur2_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 687 dim_rtwg,rhotwg) 688 689 if (psps%usepaw == 0) then 690 ! Matrix elements of i[H,r] for NC pseudopotentials. 691 rhotwx = nc_ihr_comm(vkbr(ik_ibz), cryst, psps, npw_k, nspinor, istwf_k, Ep%inclvkb, & 692 Kmesh%ibz(:,ik_ibz), ug1, ug2, kg_k) 693 694 else 695 ! 1) Add PAW onsite contribution, projectors are already in the BZ. 696 call paw_rho_tw_g(Ep%npwe,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,& 697 Cprj1_bz,Cprj2_bz,Pwij,rhotwg) 698 699 ! 2) Matrix elements of i[H,r] for PAW. 700 rhotwx = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug1,ug2,kg_k,Cprj1_ibz,Cprj2_ibz,HUr) 701 702 ! Add PAW cross term 703 if (Dtset%pawcross==1) then 704 call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,& 705 ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 706 ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,& 707 dim_rtwg,rhotwg) 708 709 ! Add cross-term contribution to the commutator 710 if (Dtset%userib/=111) then 711 call paw_cross_ihr_comm(rhotwx,nspinor,nfftf_tot,Cryst,Pawfgrtab,Paw_onsite,& 712 ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj1_ibz,Cprj2_ibz) 713 end if 714 end if 715 end if 716 717 ! Treat a possible degeneracy between v and c. 718 if (abs(deltaeKS_b1b2) > GW_TOL_W0) then 719 rhotwx = -rhotwx / deltaeKS_b1b2 720 else 721 rhotwx = czero_gw 722 end if 723 724 SELECT CASE (Ep%spmeth) 725 CASE (0) 726 ! ---------------- Ucrpa (begin) 727 if(dtset%ucrpa>=1.and..not.luwindow) then 728 fac=one 729 fac1=zero 730 fac2=zero 731 fac3=zero 732 fac4=one 733 m1=-1 734 m2=-1 735 if(dtset%ucrpa<=2) then 736 call flush_unit(std_out) 737 if ( band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)& 738 & .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 739 ! if(dtset%prtvol>=10)write(6,*)"calculation is in progress",band1,band2,ucrpa_bands(1),ucrpa_bands(2) 740 if (dtset%plowan_compute>=10) then 741 band1c=band1-wan%bandi_wan+1 742 band2c=band2-wan%bandi_wan+1 743 do iat1=1,wan%natom_wan 744 do ispinor1=1,wan%nspinor 745 do il1=1,wan%nbl_atom_wan(iat1) 746 do m1=1,2*(wan%latom_wan(iat1)%lcalc(il1))+1 747 fac1=fac1 + real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*& 748 &conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1)) 749 fac2=fac2 + real(wan%psichi(ik_bz,band2c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*& 750 &conjg(wan%psichi(ik_bz,band2c,iat1)%atom(il1)%matl(m1,spin,ispinor1)) 751 do iat2=1,wan%natom_wan 752 do ispinor2=1,wan%nspinor 753 do il2=1,wan%nbl_atom_wan(iat2) 754 do m2=1,2*(wan%latom_wan(iat2)%lcalc(il2))+1 755 fac=fac - real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1)*& 756 & conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*& 757 & wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)*& 758 & conjg(wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2))) 759 fac3=fac3+ real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*& 760 &conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*& 761 &wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)*& 762 &conjg(wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)) 763 enddo !m2 764 enddo !il2 765 enddo !ispinor2 766 enddo !iat2 767 enddo !m1 768 enddo !il1 769 enddo !ispinor1 770 enddo !iat 771 else !plowan_compute 772 do iat=1, cryst%nattyp(itypatcor) 773 do ispinor1=1,nspinor 774 do m1=1,2*lcor+1 775 fac1=fac1+ real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 776 &conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))) 777 fac2=fac2+ real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*& 778 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1))) 779 do ispinor2=1,nspinor 780 do m2=1,2*lcor+1 781 fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 782 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))*& 783 & coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*& 784 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2))) 785 fac3=fac3 + real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 786 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* & 787 & coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*& 788 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2))) 789 ! if(dtset%prtvol>=10)write(6,*) fac,fac3 790 enddo !m2 791 enddo !ispinor2 792 ! if(dtset%prtvol>=10)write(6,*) fac,fac3,fac1,fac2,fac1*fac2 793 enddo !m1 794 enddo !ispinor1 795 enddo !iat 796 endif !plowan_compute>=10 797 fac4=fac 798 ! fac=zero 799 if(dtset%ucrpa==1) fac=zero 800 ! write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0" 801 endif 802 else if (dtset%ucrpa==3) then 803 if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then 804 do iat=1, cryst%nattyp(itypatcor) 805 do ispinor1=1,nspinor 806 do m1=1,2*lcor+1 807 fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*& 808 & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))) 809 enddo 810 enddo 811 enddo 812 if(dtset%ucrpa==4) fac2=zero 813 endif 814 if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then 815 do iat=1, cryst%nattyp(itypatcor) 816 do ispinor1=1,nspinor 817 do m1=1,2*lcor+1 818 fac3=fac3-real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*& 819 & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1))) 820 enddo 821 enddo 822 enddo 823 if(dtset%ucrpa==4) fac3=zero 824 endif 825 fac=real(fac2*fac3) 826 endif 827 !if(dtset%prtvol>=10) write(6,'(6i4,e15.5,a)') ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0" 828 !if(abs(fac-one)>0.00001) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0" 829 ! if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac4," q==0" 830 green_w=green_w*fac 831 endif 832 ! ---------------- Ucrpa (end) 833 834 ! Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?) 835 call accumulate_chi0_q0(is_metallic,ik_bz,isym_k,itim_k,Ep%gwcomp,nspinor,Ep%npwepG0,Ep,& 836 Cryst,Ltg_q,Gsph_epsG0,chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing) 837 838 CASE (1, 2) 839 ! Spectral method, to be consistent here we use the KS eigenvalues. 840 call accumulate_sfchi0_q0(ik_bz,isym_k,itim_k,nspinor,Ep%symchi,Ep%npwepG0,Ep%npwe,Cryst,Ltg_q,& 841 Gsph_epsG0,deltaf_b1b2,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,Ep%nomegasf,& 842 sf_chi0,sf_head,sf_lwing,sf_uwing) 843 844 CASE DEFAULT 845 ABI_BUG("Wrong spmeth") 846 END SELECT 847 848 ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. [[cite:Mahan2000]] 849 factor = spin_fact * qp_occ(band2,ik_ibz,spin) 850 851 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1b2,& 852 Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 853 854 if (Ep%gwcomp == 1) then 855 ! Include also the completeness correction in the sum rule. 856 factor=-spin_fact*qp_occ(band2,ik_ibz,spin) 857 call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2,& 858 Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule) 859 if (band1 == Ep%nbnds) then 860 chi0_sumrule(:) = chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2 861 end if 862 end if 863 864 end do ! band2 865 end do ! band1 866 867 if (Psps%usepaw == 0 .and. Ep%inclvkb /= 0 .and. Ep%symchi == 1) then 868 call vkbr_free(vkbr(ik_ibz)) ! Not need anymore as we loop only over IBZ. 869 end if 870 end do !ik_bz 871 end do !spin 872 873 ABI_FREE(igffteps0) 874 875 call vkbr_free(vkbr) 876 ABI_FREE(vkbr) 877 878 ! === After big fat loop over transitions, now MPI === 879 ! * Master took care of the contribution in case of (metallic|spin) polarized systems. 880 select case (Ep%spmeth) 881 case (0) 882 ! Adler-Wiser expression 883 ! Sum contributions from each proc. 884 ! Looping on frequencies to avoid problems with the size of the MPI packet. 885 do io=1,Ep%nomega 886 call xmpi_sum(chi0(:,:,io),comm,ierr) 887 end do 888 889 case (1, 2) 890 ! Spectral method. 891 call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth) 892 893 ABI_SFREE(sf_chi0) 894 895 ! Sum contributions from each proc 896 ! Looping on frequencies to avoid problems with the size of the MPI packet 897 do io=1,Ep%nomega 898 call xmpi_sum(chi0(:,:,io),comm,ierr) 899 end do 900 901 call hilbert_transform_headwings(Ep%npwe,Ep%nomega,Ep%nomegasf,& 902 my_wl,my_wr,kkweight,sf_lwing,sf_uwing,sf_head,chi0_lwing,& 903 chi0_uwing,chi0_head,Ep%spmeth) 904 905 case default 906 ABI_BUG(sjoin("Wrong spmeth:", itoa(ep%spmeth))) 907 end select 908 909 ! Divide by the volume 910 !$OMP PARALLEL WORKSHARE 911 chi0 = chi0 * weight / Cryst%ucvol 912 !$OMP END PARALLEL WORKSHARE 913 914 ! Collect sum rule. pi comes from Im[1/(x-ieta)] = pi delta(x) 915 call xmpi_sum(chi0_sumrule, comm, ierr) 916 chi0_sumrule = chi0_sumrule * pi * weight / Cryst%ucvol 917 918 ! Collect head and wings. 919 call xmpi_sum(chi0_head, comm, ierr) 920 call xmpi_sum(chi0_lwing, comm, ierr) 921 call xmpi_sum(chi0_uwing, comm, ierr) 922 923 chi0_head = chi0_head * weight / cryst%ucvol 924 ! Tensor in terms of reciprocal lattice vectors. 925 do io=1,Ep%nomega 926 chi0_head(:,:,io) = matmul(chi0_head(:,:,io), cryst%gmet) * (two_pi**2) 927 end do 928 chi0_lwing = chi0_lwing * weight / cryst%ucvol 929 chi0_uwing = chi0_uwing * weight / cryst%ucvol 930 931 ! =============================================== 932 ! ==== Symmetrize chi0 in case of AFM system ==== 933 ! =============================================== 934 ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$. 935 ! Works only in the case of magnetic group Shubnikov type IV. 936 if (cryst%use_antiferro) then 937 call symmetrize_afm_chi0(Cryst, Gsph_epsG0, Ltg_q, Ep%npwe, Ep%nomega, chi0=chi0, & 938 chi0_head=chi0_head, chi0_lwing=chi0_lwing, chi0_uwing=chi0_uwing) 939 end if 940 941 ! =================================================== 942 ! ==== Construct head and wings from the tensor ===== 943 ! =================================================== 944 do io=1,Ep%nomega 945 do ig=2,Ep%npwe 946 wng = chi0_uwing(ig,io,:) 947 chi0(1,ig,io) = vdotw(Ep%qlwl(:,1), wng, Cryst%gmet,"G") 948 wng = chi0_lwing(ig,io,:) 949 chi0(ig,1,io) = vdotw(Ep%qlwl(:,1), wng, Cryst%gmet,"G") 950 end do 951 chq = matmul(chi0_head(:,:,io), Ep%qlwl(:,1)) 952 chi0(1,1,io) = vdotw(Ep%qlwl(:,1), chq, Cryst%gmet,"G") ! Use user-defined small q 953 end do 954 955 if (wfd%my_rank == 0 .and. dtset%prtvol > 20) then 956 qlen = tol3 957 call cryst%get_redcart_qdirs(nq, qdirs, qlen=qlen) 958 ABI_MALLOC(head_qvals, (nq)) 959 call wrtout([std_out, ab_out], "Head of the irreducible polarizability for q --> 0", pre_newlines=1) 960 call wrtout([std_out, ab_out], sjoin(" q0_len:", ftoa(qlen), "(Bohr^-1)")) 961 write(msg, "(*(a14))") "omega_re (eV)", "omega_im (eV)", "[100]", "[010]", "[001]", "x", "y", "z" 962 call wrtout([std_out, ab_out], msg) 963 do io=1,Ep%nomega 964 do iq=1,nq 965 chq = matmul(chi0_head(:,:,io), qdirs(:,iq)) 966 head_qvals(iq) = vdotw(qdirs(:, iq), chq, cryst%gmet, "G") 967 end do 968 write(msg, "(*(es12.5,2x))") ep%omega(io) * Ha_eV, real(head_qvals(:)) 969 call wrtout([std_out, ab_out], msg) 970 ! Write imag part to std_out 971 write(msg, "(*(es12.5,2x))") ep%omega(io) * Ha_eV, aimag(head_qvals(:)) 972 call wrtout(std_out, msg) 973 end do 974 call wrtout([std_out, ab_out], " ") 975 ABI_FREE(qdirs) 976 ABI_FREE(head_qvals) 977 end if 978 !stop 979 980 ! Impose Hermiticity (valid only for zero or purely imaginary frequencies) 981 ! MG: what about metals, where we have poles around zero? 982 ! FB: because of the intraband term, chi0 is never hermitian in case of metals 983 if (.not. is_metallic) then 984 do io=1,Ep%nomega 985 if (ABS(REAL(Ep%omega(io))) < 0.00001) then 986 do ig2=1,Ep%npwe 987 do ig1=1,ig2-1 988 chi0(ig2,ig1,io)=GWPC_CONJG(chi0(ig1,ig2,io)) 989 end do 990 end do 991 end if 992 end do 993 end if 994 995 ! ===================== 996 ! ==== Free memory ==== 997 ! ===================== 998 ABI_FREE(bbp_ks_distrb) 999 ABI_FREE(rhotwg) 1000 ABI_FREE(tabr_k) 1001 ABI_FREE(ur1_kibz) 1002 ABI_FREE(ur2_kibz) 1003 ABI_FREE(usr1_k) 1004 ABI_FREE(ur2_k) 1005 ABI_FREE(gw_gbound) 1006 1007 if (Dtset%pawcross==1) then 1008 ABI_FREE(gboundf) 1009 end if 1010 1011 ABI_SFREE(green_enhigh_w) 1012 ABI_SFREE(gw_gfft) 1013 ABI_SFREE(wfwfg) 1014 ABI_SFREE(kkweight) 1015 ABI_SFREE(omegasf) 1016 ABI_SFREE(green_w) 1017 ABI_SFREE(sf_head) 1018 ABI_SFREE(sf_lwing) 1019 ABI_SFREE(sf_uwing) 1020 ABI_SFREE(gspfft_igfft) 1021 1022 call Gsph_FFT%free() 1023 1024 if (Psps%usepaw==1) then 1025 ! deallocation for PAW. 1026 call pawcprj_free(Cprj1_bz) 1027 ABI_FREE(Cprj1_bz) 1028 call pawcprj_free(Cprj2_bz) 1029 ABI_FREE(Cprj2_bz) 1030 call pawcprj_free(Cprj1_ibz) 1031 ABI_FREE(Cprj1_ibz) 1032 call pawcprj_free(Cprj2_ibz) 1033 ABI_FREE(Cprj2_ibz) 1034 call pawpwij_free(Pwij) 1035 ABI_FREE(Pwij) 1036 if (allocated(Pwij_fft)) then 1037 call pawpwij_free(Pwij_fft) 1038 ABI_FREE(Pwij_fft) 1039 end if 1040 call pawhur_free(Hur) 1041 ABI_FREE(Hur) 1042 if (Dtset%pawcross==1) then 1043 ABI_FREE(ur_ae1) 1044 ABI_FREE(ur_ae_onsite1) 1045 ABI_FREE(ur_ps_onsite1) 1046 ABI_FREE(ur_ae2) 1047 ABI_FREE(ur_ae_onsite2) 1048 ABI_FREE(ur_ps_onsite2) 1049 ABI_FREE(tabrf_k) 1050 ABI_FREE(gboundf) 1051 ABI_FREE(igfftepsG0f) 1052 end if 1053 end if 1054 1055 if(dtset%ucrpa >= 1 .AND. dtset%plowan_compute < 10) then 1056 ABI_FREE(coeffW_BZ) 1057 endif 1058 1059 call cwtime(cpu_time, wall_time, gflops, "stop") 1060 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 1061 1062 DBG_EXIT("COLL") 1063 1064 end subroutine cchi0q0
ABINIT/chi0q0_intraband [ Functions ]
NAME
chi0q0_intraband
FUNCTION
Calculate chi0 in the limit q-->0
INPUTS
use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry) Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp) %ng=number of G vectors %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G in the array gvec %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0 %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$). %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1) Ep%nbnds=number of bands ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths. Ep%nomega=number of frequencies Cryst<crystal_t>= data type gathering info on symmetries and unit cell %natom=number of atoms %nsym=number of symmetry operations %symrec(3,3,nsym)=symmetry operations in reciprocal space %typat(natom)=type of each atom %xred(3,natom)=reduced coordinated of atoms %rprimd(3,3)=dimensional primitive translations in real space (bohr) %timrev=2 if time-reversal symmetry can be used, 1 otherwise Ep%npwe=number of planewaves for sigma exchange (input variable) Ep%nsppol=1 for unpolarized, 2 for spin-polarized Ep%omega(Ep%nomega)=frequencies Psps <type(pseudopotential_type)>=variables related to pseudopotentials %mpsang=1+maximum angular momentum for nonlocal pseudopotential Pawang<pawang_type> angular mesh discretization and related data: Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels BSt<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities) %mband=MAX number of bands over k-points and spin (==Ep%nbnds) %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
OUTPUT
chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq, and frequencies defined by Ep%omega
NOTES
*) The terms "head", "wings" and "body" of chi(G,Gp) refer to G=Gp=0, either G or Gp=0, and neither=0 respectively
TODO
Check npwepG0 before Switching on umklapp
SOURCE
1992 subroutine chi0q0_intraband(Wfd,Cryst,Ep,Psps,BSt,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,usepawu,& 1993 & ngfft_gw,chi0,chi0_head,chi0_lwing,chi0_uwing) 1994 1995 !Arguments ------------------------------------ 1996 !scalars 1997 integer,intent(in) :: usepawu 1998 logical,intent(in) :: use_tr 1999 type(ebands_t),intent(in) :: BSt 2000 type(crystal_t),intent(in) :: Cryst 2001 type(em1params_t),intent(in) :: Ep 2002 type(gsphere_t),intent(in) :: Gsph_epsG0 2003 type(Pseudopotential_type),intent(in) :: Psps 2004 type(Pawang_type),intent(in) :: Pawang 2005 type(wfdgw_t),target,intent(inout) :: Wfd 2006 !arrays 2007 integer,intent(in) :: ngfft_gw(18) 2008 complex(gwpc),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega) 2009 complex(dpc),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3) 2010 complex(dpc),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3) 2011 complex(dpc),intent(out) :: chi0_head(3,3,Ep%nomega) 2012 type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw) 2013 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 2014 type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw) 2015 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 2016 2017 !Local variables ------------------------------ 2018 !scalars 2019 integer,parameter :: tim_fourdp1=1,two_poles=2,one_pole=1,ndat1=1 2020 integer,parameter :: unitdos0=0,option1=1,NOMEGA_PRINTED=15 2021 integer :: nqlwl,nband_k,iomega,istwf_k,npw_k,my_nband,lbidx 2022 integer :: band,itim_k,ik_bz,ik_ibz,io,isym_k,spin,iqlwl !ig,ig1,ig2,my_nbbp,my_nbbpks 2023 integer :: nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,ifft 2024 integer :: kptopt,isym,nsppol,nspinor 2025 integer :: comm,ierr,gw_mgfft,use_umklp,inclvkb 2026 real(dp) :: spin_fact,deltaf_b1b2,weight 2027 real(dp) :: deltaeGW_b1b2,zcut 2028 real(dp),parameter :: dummy_dosdeltae=HUGE(zero) 2029 real(dp) :: o_entropy,o_nelect,maxocc 2030 complex(dpc) :: ph_mkt 2031 logical :: iscompatibleFFT,is_metallic !,ltest 2032 character(len=500) :: msg,msg_tmp !,allup 2033 type(kmesh_t) :: Kmesh 2034 type(littlegroup_t) :: Ltg_q 2035 type(vkbr_t) :: vkbr 2036 type(wave_t),pointer :: wave 2037 !arrays 2038 integer :: my_band_list(Wfd%mband) 2039 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 2040 integer,allocatable :: ktabr(:,:),irottb(:,:) 2041 !integer :: got(Wfd%nproc) 2042 integer,allocatable :: tabr_k(:),igffteps0(:),gw_gbound(:,:) 2043 real(dp),parameter :: q0(3)=(/zero,zero,zero/) 2044 real(dp) :: kpt(3),dedk(3),kbz(3),spinrot_kbz(4) 2045 !real(dp),ABI_CONTIGUOUS pointer :: ks_eig(:,:,:),qp_eig(:,:,:),qp_occ(:,:,:) 2046 real(dp) :: shift_ene(BSt%mband,BSt%nkpt,BSt%nsppol) 2047 real(dp) :: delta_occ(BSt%mband,BSt%nkpt,BSt%nsppol) 2048 !real(dp) :: eigen_vec(BSt%bantot) 2049 real(dp) :: o_doccde(BSt%bantot) 2050 real(dp) :: eigen_pdelta_vec(BSt%bantot),eigen_mdelta_vec(BSt%bantot) 2051 real(dp) :: o_occ_pdelta(BSt%bantot),o_occ_mdelta(BSt%bantot) 2052 real(dp) :: delta_ene(BSt%mband,BSt%nkpt,BSt%nsppol) 2053 real(dp) :: test_docc(BSt%mband,BSt%nkpt,BSt%nsppol) 2054 real(dp),allocatable :: qlwl(:,:) 2055 complex(gwpc) :: comm_kbbs(3,Wfd%nspinor**2) 2056 complex(dpc),allocatable :: ihr_comm(:,:,:,:,:) 2057 complex(gwpc),allocatable :: rhotwg(:) 2058 complex(dpc) :: green_w(Ep%nomega) 2059 complex(gwpc),allocatable :: ur1(:) 2060 complex(gwpc),ABI_CONTIGUOUS pointer :: ug(:) 2061 logical :: bmask(Wfd%mband) 2062 type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj1_ibz(:,:),Cp_bks(:,:) 2063 type(pawpwij_t),allocatable :: Pwij(:) 2064 type(pawhur_t),allocatable :: Hur(:) 2065 2066 !************************************************************************ 2067 2068 DBG_ENTER("COLL") 2069 2070 nsppol = Wfd%nsppol 2071 nspinor = Wfd%nspinor 2072 is_metallic = ebands_has_metal_scheme(BSt) 2073 2074 gw_mgfft = MAXVAL(ngfft_gw(1:3)) 2075 gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10) 2076 2077 ! Calculate <k,b1|i[H,r]|k',b2>. 2078 inclvkb=2; if (Wfd%usepaw==1) inclvkb=0 2079 ABI_MALLOC(ihr_comm,(3,nspinor**2,Wfd%mband,Wfd%nkibz,nsppol)) 2080 ihr_comm = czero 2081 2082 if (Wfd%usepaw==1) then 2083 ABI_MALLOC(Cp_bks,(Cryst%natom,nspinor)) 2084 call pawcprj_alloc(Cp_bks,0,Wfd%nlmn_atm) 2085 ABI_MALLOC(HUr,(Cryst%natom)) 2086 if (usepawu/=0) then ! For PAW+DFT+U, precalculate <\phi_i|[Hu,r]|phi_j\>. 2087 call pawhur_init(hur,nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij) 2088 end if 2089 end if 2090 2091 do spin=1,nsppol 2092 do ik_ibz=1,Wfd%nkibz 2093 npw_k = Wfd%npwarr(ik_ibz) 2094 nband_k= Wfd%nband(ik_ibz,spin) 2095 kpt = Wfd%kibz(:,ik_ibz) 2096 kg_k => Wfd%Kdata(ik_ibz)%kg_k 2097 istwf_k = Wfd%istwfk(ik_ibz) 2098 ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded") 2099 2100 ! Distribute bands. 2101 bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included. 2102 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,bmask=bmask) 2103 if (my_nband==0) CYCLE 2104 2105 if (Wfd%usepaw==0.and.inclvkb/=0) then ! Include term <n,k|[Vnl,iqr]|n"k>' for q->0. 2106 call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,kpt,kg_k) 2107 end if 2108 2109 do lbidx=1,my_nband 2110 band=my_band_list(lbidx) 2111 2112 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 2113 ug => wave%ug 2114 2115 if (Wfd%usepaw==0) then 2116 ! Matrix elements of i[H,r] for NC pseudopotentials. 2117 comm_kbbs = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug,ug,kg_k) 2118 else 2119 ! Matrix elements of i[H,r] for PAW. 2120 call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cp_bks,sorted=.FALSE.) 2121 comm_kbbs = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug,ug,kg_k,Cp_bks,Cp_bks,HUr) 2122 end if 2123 2124 ihr_comm(:,:,band,ik_ibz,spin) = comm_kbbs 2125 end do 2126 2127 call vkbr_free(vkbr) ! Not need anymore as we loop only over IBZ. 2128 end do 2129 end do 2130 ! 2131 ! Gather the commutator on each node. 2132 call xmpi_sum(ihr_comm,Wfd%comm,ierr) 2133 2134 if (Wfd%usepaw==1) then 2135 call pawcprj_free(Cp_bks) 2136 ABI_FREE(Cp_bks) 2137 call pawhur_free(Hur) 2138 ABI_FREE(Hur) 2139 end if 2140 2141 nqlwl=1 2142 ABI_MALLOC(qlwl,(3,nqlwl)) 2143 !qlwl = GW_Q0_DEFAULT(3) 2144 qlwl(:,1) = (/0.00001_dp, 0.00002_dp, 0.00003_dp/) 2145 ! 2146 write(msg,'(a,i3,a)')' Q-points for long wave-length limit in chi0q_intraband. # ',nqlwl,ch10 2147 do iqlwl=1,nqlwl 2148 write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',qlwl(:,iqlwl),ch10 2149 msg=TRIM(msg)//msg_tmp 2150 end do 2151 call wrtout(std_out, msg) 2152 ! 2153 ! delta_ene = e_{b,k-q} - e_{b,k} = -q. <b,k| i[H,r] |b,k> + O(q^2). 2154 delta_ene = zero 2155 do spin=1,nsppol 2156 do ik_ibz=1,Wfd%nkibz 2157 do band=1,Wfd%nband(ik_ibz,spin) 2158 dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin)) 2159 delta_ene(band,ik_ibz,spin) = -vdotw(qlwl(:,1),dedk,Cryst%gmet,"G") 2160 end do 2161 end do 2162 end do 2163 2164 maxocc=two/(nsppol*nspinor) 2165 2166 ! Calculate the occupations at f(e+delta/2). 2167 shift_ene = BSt%eig + half*delta_ene 2168 2169 call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_pdelta_vec) 2170 2171 if (BSt%occopt < 9) then 2172 call getnel(o_doccde,dummy_dosdeltae,eigen_pdelta_vec,o_entropy,BSt%fermie,BSt%fermie,maxocc,BSt%mband,BSt%nband,& 2173 & o_nelect,BSt%nkpt,BSt%nsppol,o_occ_pdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk,1,BSt%nband(1)) 2174 ! CP: adding 1 and BSt%nband(1) as dummy arguments since here we already test for occopt==9 2175 !write(std_out,*)"nelect1: ",o_nelect 2176 else 2177 ABI_ERROR('occopt 9 not implemented for GW calculations') 2178 end if 2179 ! 2180 ! Calculate the occupations at f(e-delta/2). 2181 shift_ene = BSt%eig - half*delta_ene 2182 2183 call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_mdelta_vec) 2184 2185 ! CP modified 2186 !call getnel(o_doccde,dummy_dosdeltae,eigen_mdelta_vec,o_entropy,BSt%fermie,maxocc,BSt%mband,BSt%nband,& 2187 !& o_nelect,BSt%nkpt,BSt%nsppol,o_occ_mdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk) 2188 !write(std_out,*)"nelect2: ",o_nelect 2189 if (BSt%occopt < 9) then 2190 call getnel(o_doccde,dummy_dosdeltae,eigen_mdelta_vec,o_entropy,BSt%fermie,BSt%fermie,maxocc,BSt%mband,BSt%nband,& 2191 & o_nelect,BSt%nkpt,BSt%nsppol,o_occ_mdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk,1,BSt%nband(1)) 2192 write(std_out,*)"nelect2: ",o_nelect 2193 else 2194 ABI_ERROR("occopt 9 not implemented for GW calculations") 2195 end if 2196 ! End CP modified 2197 2198 ! f(e-delta/2) - f(e+delta/2). 2199 o_occ_pdelta = o_occ_mdelta - o_occ_pdelta 2200 2201 call unpack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,o_occ_pdelta,delta_occ) 2202 ! 2203 ! Expand f(e-delta/2) - f(e+delta/2) up to the first order in the small q. 2204 do spin=1,nsppol 2205 do ik_ibz=1,Wfd%nkibz 2206 do band=1,Wfd%nband(ik_ibz,spin) 2207 dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin)) 2208 test_docc(band,ik_ibz,spin) = +vdotw(qlwl(:,1),dedk,Cryst%gmet,"G") * BSt%doccde(band,ik_ibz,spin) 2209 write(std_out,'(a,3(i0,1x),1x,3es16.8)')" spin,ik_ibz,band, delta_occ: ",& 2210 & spin,ik_ibz,band,delta_occ(band,ik_ibz,spin),& 2211 & test_docc(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin)-test_docc(band,ik_ibz,spin) 2212 end do 2213 end do 2214 end do 2215 2216 ! ABI_ERROR("DONE") 2217 ! do spin=1,nsppol 2218 ! do ik_ibz=1,Wfd%nkibz 2219 ! nband_k = Wfd%nband(ik_ibz,spin) 2220 ! do band=1,nband_k 2221 ! write(std_out,'(a,3i3,2es14.6)')" spin, band, ik_ibz, delta_ene, delta_occ ",& 2222 !& spin,band,ik_ibz,delta_ene(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin) 2223 ! end do 2224 ! end do 2225 ! end do 2226 2227 ABI_FREE(ihr_comm) 2228 ABI_FREE(qlwl) 2229 2230 if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,ngfft_gw) 2231 2232 ! TODO take into account the case of random k-meshes. 2233 kptopt=3 2234 call Kmesh%init(Cryst,Wfd%nkibz,Wfd%kibz,kptopt) 2235 ! 2236 !=== Get the FFT index of $ (R^{-1}(r-\tau)) $ === 2237 !* S= $\transpose R^{-1}$ and k_BZ = S k_IBZ 2238 !* irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk. 2239 ABI_MALLOC(irottb,(Wfd%nfftot,Cryst%nsym)) 2240 2241 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT) 2242 ABI_CHECK(iscompatibleFFT,"FFT mesh not compatible with symmetries") 2243 2244 ABI_MALLOC(ktabr,(Wfd%nfftot,Kmesh%nbz)) 2245 do ik_bz=1,Kmesh%nbz 2246 isym=Kmesh%tabo(ik_bz) 2247 do ifft=1,Wfd%nfftot 2248 ktabr(ifft,ik_bz)=irottb(ifft,isym) 2249 end do 2250 end do 2251 ABI_FREE(irottb) 2252 ! 2253 ! === Setup weight (2 for spin unpolarized systems, 1 for polarized) === 2254 ! spin_fact is used to normalize the occupation factors to one. 2255 ! Consider also the AFM case. 2256 SELECT CASE (nsppol) 2257 CASE (1) 2258 weight=two/Kmesh%nbz; spin_fact=half 2259 if (Wfd%nspden==2) then 2260 weight=one/Kmesh%nbz; spin_fact=half 2261 end if 2262 if (nspinor==2) then 2263 weight=one/Kmesh%nbz; spin_fact=one 2264 end if 2265 CASE (2) 2266 weight=one/Kmesh%nbz; spin_fact=one 2267 CASE DEFAULT 2268 ABI_BUG("Wrong nsppol") 2269 END SELECT 2270 2271 use_umklp = 0 2272 call Ltg_q%init(q0, Kmesh%nbz, Kmesh%bz, Cryst, use_umklp, Ep%npwepG0, gvec=Gsph_epsG0%gvec) 2273 2274 write(msg,'(a,i2)')' Using symmetries to sum only over the IBZ_q = ',Ep%symchi 2275 call wrtout(std_out, msg) 2276 ! 2277 ! Evaluate oscillator matrix elements btw partial waves. Note that q=Gamma is used. 2278 if (Psps%usepaw==1) then 2279 ABI_MALLOC(Pwij,(Psps%ntypat)) 2280 call pawpwij_init(Pwij,Ep%npwepG0, [zero,zero,zero], Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 2281 2282 ABI_MALLOC(Cprj1_bz ,(Cryst%natom,nspinor)) 2283 call pawcprj_alloc(Cprj1_bz, 0,Wfd%nlmn_atm) 2284 ABI_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor)) 2285 call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm) 2286 end if 2287 2288 ABI_MALLOC(rhotwg,(Ep%npwe*nspinor**2)) 2289 ABI_MALLOC(tabr_k,(Wfd%nfftot)) 2290 ABI_MALLOC(ur1,(Wfd%nfft*nspinor)) 2291 ! 2292 ! Tables for the FFT of the oscillators. 2293 ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere). 2294 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 2295 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2)) 2296 ABI_MALLOC(igffteps0,(Gsph_epsG0%ng)) 2297 2298 call Gsph_epsG0%fft_tabs([0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0) 2299 if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 2300 if (use_padfft==0) then 2301 ABI_FREE(gw_gbound) 2302 ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft)) 2303 end if 2304 2305 nkpt_summed=Kmesh%nbz 2306 if (Ep%symchi/=0) then 2307 nkpt_summed=Ltg_q%nibz_ltg 2308 call Ltg_q%print(std_out, Wfd%prtvol) 2309 end if 2310 ! 2311 ! ============================================ 2312 ! === Begin big fat loop over transitions ==== 2313 ! ============================================ 2314 chi0 = czero_gw 2315 chi0_head = czero_gw; chi0_lwing = czero_gw; chi0_uwing = czero_gw 2316 dim_rtwg=1; if (nspinor==2) dim_rtwg=2 !can reduce size depending on Ep%nI and Ep%nj 2317 2318 zcut = Ep%zcut 2319 zcut = 0.1/Ha_eV 2320 write(std_out,*)" using zcut ",zcut*Ha_eV," [eV]" 2321 2322 ! Loop on spin to calculate $ \chi_{\up,\up} + \chi_{\down,\down} 2323 do spin=1,nsppol 2324 ! Loop over k-points in the BZ. 2325 do ik_bz=1,Kmesh%nbz 2326 if (Ep%symchi==1) then 2327 if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q 2328 end if 2329 2330 ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz. 2331 call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt) 2332 tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points 2333 spinrot_kbz(:)=Cryst%spinrot(:,isym_k) 2334 nband_k=Wfd%nband(ik_ibz,spin) 2335 2336 ! Distribute bands. 2337 bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included. 2338 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,bmask=bmask) 2339 if (my_nband==0) CYCLE 2340 2341 write(msg,'(2(a,i4),a,i2,a,i3)')' ik = ',ik_bz,' / ',Kmesh%nbz,' spin = ',spin,' done by processor ',Wfd%my_rank 2342 call wrtout(std_out, msg) 2343 2344 do lbidx=1,my_nband 2345 ! Loop over bands treated by this node. 2346 band=my_band_list(lbidx) 2347 call wfd%get_ur(band,ik_ibz,spin,ur1) 2348 2349 if (Psps%usepaw==1) then 2350 call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.) 2351 call pawcprj_copy(Cprj1_ibz,Cprj1_bz) 2352 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz) 2353 end if 2354 2355 deltaf_b1b2 = spin_fact*delta_occ(band,ik_ibz,spin) 2356 deltaeGW_b1b2= delta_ene(band,ik_ibz,spin) 2357 2358 ! Add small imaginary of the Time-Ordered resp function but only for non-zero real omega FIXME What about metals? 2359 if (.not.use_tr) then 2360 do io=1,Ep%nomega 2361 !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,one_pole) 2362 green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,one_pole) 2363 end do 2364 else 2365 do io=1,Ep%nomega ! This expression implements time-reversal even when the input k-mesh breaks it. 2366 !green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,two_poles) 2367 green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,two_poles) 2368 end do !io 2369 end if ! use_tr 2370 2371 ! FFT of u^*_{b1,k}(r) u_{b2,k}(r). 2372 call rho_tw_g(nspinor,Ep%npwe,Wfd%nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,& 2373 & ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 2374 & ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,& 2375 & dim_rtwg,rhotwg) 2376 2377 if (Psps%usepaw==1) then 2378 ! Add PAW onsite contribution, projectors are already in the BZ. 2379 call paw_rho_tw_g(Ep%npwe,dim_rtwg,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_epsG0%gvec,& 2380 & Cprj1_bz,Cprj1_bz,Pwij,rhotwg) 2381 end if 2382 2383 ! ==== Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?) ==== 2384 call assemblychi0_sym(is_metallic,ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0) 2385 end do !band 2386 end do !ik_bz 2387 end do !spin 2388 2389 ! Collect body, head and wings within comm 2390 comm=Wfd%comm 2391 do io=1,Ep%nomega 2392 call xmpi_sum(chi0(:,:,io),comm,ierr) 2393 end do 2394 call xmpi_sum(chi0_head,comm,ierr) 2395 call xmpi_sum(chi0_lwing,comm,ierr) 2396 call xmpi_sum(chi0_uwing,comm,ierr) 2397 2398 ! Divide by the volume 2399 chi0 = chi0 * weight/Cryst%ucvol 2400 chi0_head = chi0_head * weight/Cryst%ucvol 2401 do io=1,Ep%nomega ! Tensor in the basis of the reciprocal lattice vectors. 2402 chi0_head(:,:,io) = MATMUL(chi0_head(:,:,io),Cryst%gmet) * (two_pi**2) 2403 end do 2404 chi0_lwing = chi0_lwing * weight/Cryst%ucvol 2405 chi0_uwing = chi0_uwing * weight/Cryst%ucvol 2406 ! 2407 ! =============================================== 2408 ! ==== Symmetrize chi0 in case of AFM system ==== 2409 ! =============================================== 2410 ! * Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$. 2411 ! * Works only in the case of magnetic group Shubnikov type IV. 2412 if (Cryst%use_antiferro) then 2413 call symmetrize_afm_chi0(Cryst, Gsph_epsG0, Ltg_q, Ep%npwe, Ep%nomega, chi0=chi0, & 2414 chi0_head=chi0_head, chi0_lwing=chi0_lwing, chi0_uwing=chi0_uwing) 2415 end if 2416 ! 2417 ! ================================================== 2418 ! ==== Construct head and wings from the tensor ==== 2419 ! ================================================== 2420 !do io=1,Ep%nomega 2421 ! do ig=2,Ep%npwe 2422 ! wng = chi0_uwing(ig,io,:) 2423 ! chi0(1,ig,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G") 2424 ! wng = chi0_lwing(ig,io,:) 2425 ! chi0(ig,1,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G") 2426 ! end do 2427 ! chq = MATMUL(chi0_head(:,:,io), Ep%qlwl(:,1)) 2428 ! chi0(1,1,io) = vdotw(Ep%qlwl(:,1),chq,Cryst%gmet,"G") ! Use user-defined small q 2429 !end do 2430 !call wfd_barrier(Wfd) 2431 2432 ! Impose Hermiticity (valid only for zero or purely imaginary frequencies) 2433 ! MG what about metals, where we have poles around zero? 2434 !do io=1,Ep%nomega 2435 ! if (ABS(REAL(Ep%omega(io)))<0.00001) then 2436 ! do ig2=1,Ep%npwe 2437 ! do ig1=1,ig2-1 2438 ! chi0(ig2,ig1,io)=CONJG(chi0(ig1,ig2,io)) 2439 ! end do 2440 ! end do 2441 ! end if 2442 !end do 2443 2444 do iomega=1,MIN(Ep%nomega,NOMEGA_PRINTED) 2445 write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0_intra(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]' 2446 call wrtout(std_out, msg) 2447 call print_arr(chi0(:,:,iomega),unit=std_out) 2448 end do 2449 2450 ! ===================== 2451 ! ==== Free memory ==== 2452 ! ===================== 2453 ABI_FREE(rhotwg) 2454 ABI_FREE(tabr_k) 2455 ABI_FREE(ur1) 2456 ABI_FREE(gw_gbound) 2457 ABI_FREE(ktabr) 2458 ABI_FREE(igffteps0) 2459 2460 ! deallocation for PAW. 2461 if (Psps%usepaw==1) then 2462 call pawcprj_free(Cprj1_bz) 2463 ABI_FREE(Cprj1_bz) 2464 call pawcprj_free(Cprj1_ibz) 2465 ABI_FREE(Cprj1_ibz) 2466 call pawpwij_free(Pwij) 2467 ABI_FREE(Pwij) 2468 end if 2469 2470 call Ltg_q%free() 2471 call Kmesh%free() 2472 2473 DBG_EXIT("COLL") 2474 2475 end subroutine chi0q0_intraband
ABINIT/m_chi0 [ Modules ]
NAME
m_chi0
FUNCTION
COPYRIGHT
Copyright (C) 1999-2022 ABINIT group (GMR, VO, LR, RWG, MG, RShaltaf) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_chi0 24 25 use defs_basis 26 use m_abicore 27 use m_xmpi 28 use m_errors 29 use m_hide_blas 30 use m_time 31 use m_wfd 32 use m_dtset 33 34 use defs_datatypes, only : pseudopotential_type, ebands_t 35 use m_fstrings, only : ftoa, sjoin, itoa 36 use m_gwdefs, only : GW_TOL_DOCC, GW_TOL_W0, czero_gw, em1params_t, g0g0w 37 use m_numeric_tools, only : imin_loc, print_arr 38 use m_geometry, only : normv, vdotw 39 use m_crystal, only : crystal_t 40 use m_fft_mesh, only : rotate_FFT_mesh, get_gftt 41 use m_occ, only : getnel 42 use m_ebands, only : pack_eneocc, unpack_eneocc, ebands_has_metal_scheme 43 use m_bz_mesh, only : kmesh_t, littlegroup_t 44 use m_gsphere, only : gsphere_t 45 use m_io_tools, only : flush_unit 46 use m_oscillators, only : rho_tw_g, calc_wfwfg 47 use m_vkbr, only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm 48 use m_chi0tk, only : hilbert_transform, setup_spectral, assemblychi0_sym, assemblychi0sf, symmetrize_afm_chi0, & 49 approxdelta, completechi0_deltapart, accumulate_chi0sumrule, make_transitions, & 50 chi0_bbp_mask, accumulate_chi0_q0, accumulate_sfchi0_q0, hilbert_transform_headwings 51 use m_pawang, only : pawang_type 52 use m_pawrad, only : pawrad_type 53 use m_pawtab, only : pawtab_type 54 use m_paw_ij, only : paw_ij_type 55 use m_pawfgrtab, only : pawfgrtab_type 56 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 57 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g 58 use m_paw_sym, only : paw_symcprj 59 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t 60 use m_paw_hr, only : pawhur_t, pawhur_free, pawhur_init, paw_ihr, paw_cross_ihr_comm 61 use m_read_plowannier, only : read_plowannier 62 use m_plowannier, only : plowannier_type 63 64 implicit none 65 66 private