TABLE OF CONTENTS
- ABINIT/m_bethe_salpeter
- m_bethe_salpeter/bethe_salpeter
- m_bethe_salpeter/setup_bse
- m_bethe_salpeter/setup_bse_interp
ABINIT/m_bethe_salpeter [ Modules ]
NAME
m_bethe_salpeter
FUNCTION
Main routine to calculate dielectric properties by solving the Bethe-Salpeter equation in Frequency-Reciprocal space on a transition (electron-hole) basis set.
COPYRIGHT
Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida) Copyright (C) 2009-2022 ABINIT group (MG, YG) 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
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_bethe_salpeter 26 27 use defs_basis 28 use defs_wvltypes 29 use m_bs_defs 30 use m_abicore 31 use m_xmpi 32 use m_errors 33 use m_screen 34 use m_nctk 35 use m_distribfft 36 use netcdf 37 use m_hdr 38 use m_dtset 39 use m_dtfil 40 use m_crystal 41 42 use defs_datatypes, only : pseudopotential_type, ebands_t 43 use defs_abitypes, only : MPI_type 44 use m_gwdefs, only : GW_Q0_DEFAULT 45 use m_time, only : timab 46 use m_fstrings, only : strcat, sjoin, endswith, itoa 47 use m_io_tools, only : file_exists, iomode_from_fname 48 use m_geometry, only : mkrdim, metric, normv 49 use m_hide_lapack, only : matrginv 50 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 51 use m_fftcore, only : print_ngfft 52 use m_fft_mesh, only : rotate_FFT_mesh, get_gftt, setmesh 53 use m_fft, only : fourdp 54 use m_bz_mesh, only : kmesh_t, get_ng0sh, find_qmesh, make_mesh 55 use m_double_grid, only : double_grid_t, double_grid_init, double_grid_free 56 use m_ebands, only : ebands_init, ebands_print, ebands_copy, ebands_free, & 57 ebands_update_occ, ebands_get_valence_idx, ebands_apply_scissors, ebands_report_gap 58 use m_kg, only : getph 59 use m_gsphere, only : gsphere_t 60 use m_vcoul, only : vcoul_t 61 use m_qparticles, only : rdqps, rdgw !, show_QP , rdgw 62 use m_wfd, only : wfd_init, wfdgw_t, test_charge 63 use m_wfk, only : wfk_read_eigenvalues 64 use m_energies, only : energies_type, energies_init 65 use m_io_screening, only : hscr_t, hscr_io 66 use m_haydock, only : exc_haydock_driver 67 use m_exc_diago, only : exc_diago_driver 68 use m_exc_analyze, only : exc_den 69 use m_eprenorms, only : eprenorms_t, eprenorms_free, eprenorms_from_epnc, eprenorms_bcast 70 use m_pawang, only : pawang_type 71 use m_pawrad, only : pawrad_type 72 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 73 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 74 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 75 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init 76 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free,& 77 & pawrhoij_inquire_dim, pawrhoij_symrhoij 78 use m_pawdij, only : pawdij, symdij 79 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 80 use m_paw_hr, only : pawhur_t, pawhur_free, pawhur_init 81 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 82 use m_paw_sphharm, only : setsym_ylm 83 use m_paw_denpot, only : pawdenpot 84 use m_paw_init, only : pawinit,paw_gencond 85 use m_paw_onsite, only : pawnabla_init 86 use m_paw_dmft, only : paw_dmft_type 87 use m_paw_mkrho, only : denfgr 88 use m_paw_nhat, only : nhatgrid,pawmknhat 89 use m_paw_tools, only : chkpawovlp, pawprt 90 use m_paw_correlations,only : pawpuxinit 91 use m_exc_build, only : exc_build_ham 92 use m_setvtr, only : setvtr 93 use m_mkrho, only : prtrhomxmn 94 use m_pspini, only : pspini 95 use m_drivexc, only : mkdenpos 96 97 implicit none 98 99 private
m_bethe_salpeter/bethe_salpeter [ Functions ]
[ Top ] [ m_bethe_salpeter ] [ Functions ]
NAME
bethe_salpeter
FUNCTION
Main routine to calculate dielectric properties by solving the Bethe-Salpeter equation in Frequency-Reciprocal space on a transition (electron-hole) basis set.
INPUTS
acell(3)=Length scales of primitive translations (bohr) codvsn=Code version Dtfil<datafiles_type>=Variables related to files. Dtset<dataset_type>=All input variables for this dataset. Pawang<pawang_type)>=PAW angular mesh and related data. Pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. Pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. Psps<pseudopotential_type>=Variables related to pseudopotentials. Before entering the first time in the routine, a significant part of Psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in the call to pspini. The next time the code enters bethe_salpeter, Psps might be identical to the one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90. rprim(3,3)=Dimensionless real space primitive translations. xred(3,natom)=Reduced atomic coordinates. Input files used during the calculation. KSS : Kohn Sham electronic structure file. SCR (SUSC) : Files containing the symmetrized epsilon^-1 or the irreducible RPA polarizability, respectively. Used to construct the screening W. GW file : Optional file with the GW QP corrections.
OUTPUT
Output is written on the main output file and on the following external files: * _RPA_NLF_MDF: macroscopic RPA dielectric function without non-local field effects. * _GW_NLF_MDF: macroscopic RPA dielectric function without non-local field effects calculated with GW energies or the scissors operator. * _EXC_MDF: macroscopic dielectric function with excitonic effects obtained by solving the Bethe-Salpeter problem at different level of sophistication.
NOTES
ON THE USE OF FFT GRIDS: ================= In case of PAW: --------------- Two FFT grids are used: - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis. It is defined by nfft, ngfft, mgfft, ... Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid. - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres. It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid. In case of norm-conserving: --------------------------- - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ... For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
166 subroutine bethe_salpeter(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim,xred) 167 168 !Arguments ------------------------------------ 169 !scalars 170 character(len=8),intent(in) :: codvsn 171 type(datafiles_type),intent(inout) :: Dtfil 172 type(dataset_type),intent(inout) :: Dtset 173 type(pawang_type),intent(inout) :: Pawang 174 type(pseudopotential_type),intent(inout) :: Psps 175 !arrays 176 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,Dtset%natom) 177 type(pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw) 178 type(pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw) 179 180 !Local variables ------------------------------ 181 !scalars 182 integer,parameter :: tim_fourdp0=0,level=40,ipert0=0,idir0=0,cplex1=1,master=0,option1=1 183 integer :: band,cplex_rhoij,spin,ik_ibz,mqmem,iwarn 184 integer :: has_dijU,has_dijso,gnt_option 185 integer :: ik_bz,mband 186 integer :: choice 187 integer :: ider 188 integer :: usexcnhat,nfft_osc,mgfft_osc 189 integer :: isym,izero 190 integer :: optcut,optgr0,optgr1,optgr2,option,optrad,optrhoij,psp_gencond 191 integer :: ngrvdw,nhatgrdim,nkxc1,nprocs,nspden_rhoij,nzlmopt,ifft 192 integer :: my_rank,rhoxsp_method,comm 193 integer :: mgfftf,spin_opt,which_fixed 194 integer :: nscf,nbsc,nkxc,n3xccc 195 integer :: nfftf,nfftf_tot,nfftot_osc,my_minb,my_maxb 196 integer :: optene,moved_atm_inside,moved_rhor,initialized,istep,ierr 197 real(dp) :: ucvol,drude_plsmf,ecore,ecut_eff,ecutdg_eff,norm 198 real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp 199 real(dp) :: compch_fft,compch_sph,gsq_osc 200 real(dp) :: vxcavg 201 logical :: iscompatibleFFT,is_dfpt=.false.,paw_add_onsite,call_pawinit 202 character(len=500) :: msg 203 character(len=fnlen) :: wfk_fname,w_fname 204 type(Pawfgr_type) :: Pawfgr 205 type(excfiles) :: BS_files 206 type(excparam) :: BSp 207 type(paw_dmft_type) :: Paw_dmft 208 type(MPI_type) :: MPI_enreg_seq 209 type(crystal_t) :: Cryst 210 type(kmesh_t) :: Kmesh,Qmesh 211 type(gsphere_t) :: Gsph_x,Gsph_c,Gsph_x_dense,Gsph_c_dense 212 type(Hdr_type) :: Hdr_wfk,Hdr_bse 213 type(ebands_t) :: KS_BSt,QP_BSt,KS_BSt_dense,QP_BSt_dense 214 type(Energies_type) :: KS_energies 215 type(vcoul_t) :: Vcp, Vcp_dense 216 type(wfdgw_t) :: Wfd, Wfd_dense 217 type(screen_t) :: W 218 type(screen_info_t) :: W_info 219 type(wvl_data) :: wvl 220 type(kmesh_t) :: Kmesh_dense,Qmesh_dense 221 type(Hdr_type) :: Hdr_wfk_dense 222 type(double_grid_t) :: grid 223 type(eprenorms_t) :: Epren 224 !arrays 225 integer :: ngfft_osc(18),ngfftc(18),ngfftf(18),nrcell(3) 226 integer,allocatable :: ktabr(:,:),l_size_atm(:) 227 integer,allocatable :: nband(:,:),nq_spl(:),irottb(:,:) 228 integer,allocatable :: qp_vbik(:,:) 229 integer,allocatable :: gfft_osc(:,:) 230 real(dp),parameter :: k0(3)=zero 231 real(dp) :: tsec(2),gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3),rprimd(3,3),eh_rcoord(3),strsxc(6) 232 real(dp),allocatable :: ph1df(:,:),prev_rhor(:,:),ph1d(:,:) 233 real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:),ks_rhor(:,:),qp_aerhor(:,:) 234 real(dp),allocatable :: qp_rhor(:,:),qp_rhog(:,:) !,qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:) 235 real(dp),allocatable :: qp_rhor_paw(:,:),qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:),qp_nhat(:,:) 236 real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:) 237 real(dp),allocatable :: vpsp(:),xccc3d(:) 238 real(dp),allocatable :: ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:) 239 real(dp),allocatable :: kxc(:,:) !,qp_kxc(:,:) 240 complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:) 241 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:) 242 type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:) 243 type(Pawrhoij_type),allocatable :: prev_Pawrhoij(:) !QP_pawrhoij(:), 244 type(pawpwff_t),allocatable :: Paw_pwff(:) 245 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 246 type(pawhur_t),allocatable :: Hur(:) 247 type(Paw_ij_type),allocatable :: KS_paw_ij(:) 248 type(Paw_an_type),allocatable :: KS_paw_an(:) 249 250 !************************************************************************ 251 252 DBG_ENTER('COLL') 253 254 call timab(650,1,tsec) ! bse(Total) 255 call timab(651,1,tsec) ! bse(Init1) 256 257 comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 258 259 wfk_fname = dtfil%fnamewffk 260 261 if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then 262 ABI_ERROR(msg) 263 end if 264 call xmpi_bcast(wfk_fname, master, comm, ierr) 265 266 write(msg,'(8a)')& 267 ' Exciton: Calculation of dielectric properties by solving the Bethe-Salpeter equation ',ch10,& 268 ' in frequency domain and reciprocal space on a transitions basis set. ',ch10,& 269 ' Based on a program developed by L. Reining, V. Olevano, F. Sottile, ',ch10,& 270 ' S. Albrecht, and G. Onida. Incorporated in ABINIT by M. Giantomassi. ',ch10 271 call wrtout([std_out, ab_out], msg) 272 273 #ifdef HAVE_GW_DPC 274 if (gwpc/=8) then 275 write(msg,'(6a)')ch10,& 276 ' Number of bytes for double precision complex /=8 ',ch10,& 277 ' Cannot continue due to kind mismatch in BLAS library ',ch10,& 278 ' Some BLAS interfaces are not generated by abilint ' 279 ABI_ERROR(msg) 280 end if 281 write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10 282 #else 283 write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10 284 #endif 285 call wrtout([std_out, ab_out], msg) 286 287 !=== Some variables need to be initialized/nullify at start === 288 call energies_init(KS_energies) 289 usexcnhat=0 290 call mkrdim(acell,rprim,rprimd) 291 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 292 ! 293 !=== Define FFT grid(s) sizes === 294 !* Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the 295 !(usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 296 !See also NOTES in the comments at the beginning of this file. 297 !NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 298 299 !TODO Recheck getng, should use same trick as that used in screening and sigma. 300 call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 301 gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0) 302 303 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 304 nfftf_tot=PRODUCT(ngfftf(1:3)) 305 306 ! Fake MPI_type for the sequential part. 307 call initmpi_seq(MPI_enreg_seq) 308 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 309 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 310 311 ! =========================================== 312 ! === Open and read pseudopotential files === 313 ! =========================================== 314 call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm) 315 316 ! === Initialization of basic objects including the BSp structure that defines the parameters of the run === 317 call setup_bse(codvsn,acell,rprim,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,& 318 Cryst,Kmesh,Qmesh,KS_BSt,QP_BSt,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,wvl%descr) 319 320 if (BSp%use_interp) then 321 call setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh,Kmesh_dense,& 322 Qmesh_dense,KS_BSt_dense,QP_BSt_dense,Gsph_x_dense,Gsph_c_dense,& 323 Vcp_dense,Hdr_wfk_dense,grid,comm) 324 end if 325 326 !call timab(652,2,tsec) ! setup_bse 327 328 nfftot_osc=PRODUCT(ngfft_osc(1:3)) 329 nfft_osc =nfftot_osc !no FFT // 330 mgfft_osc =MAXVAL(ngfft_osc(1:3)) 331 332 call print_ngfft(ngfft_osc,header='FFT mesh used for oscillator strengths') 333 334 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 335 KS_energies%e_corepsp=ecore/Cryst%ucvol 336 337 ! 338 !============================ 339 !==== PAW initialization ==== 340 !============================ 341 if (Dtset%usepaw==1) then 342 call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,xred) 343 344 ABI_MALLOC(KS_Pawrhoij,(Cryst%natom)) 345 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 346 nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc) 347 call pawrhoij_alloc(KS_Pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 348 349 ! Initialize values for several basic arrays === 350 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 351 352 ! Test if we have to call pawinit 353 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 354 355 if (psp_gencond==1.or.call_pawinit) then 356 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 357 call pawinit(Dtset%effmass_free,gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,& 358 Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,& 359 Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%ixc,Dtset%usepotzero) 360 361 ! Update internal values 362 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 363 else 364 if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2 365 if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2 366 end if 367 Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore) 368 369 ! Initialize optional flags in Pawtab to zero 370 ! (Cannot be done in Pawinit since the routine is called only if some parameters are changed) 371 Pawtab(:)%has_nabla = 0 372 Pawtab(:)%usepawu = 0 373 Pawtab(:)%useexexch = 0 374 Pawtab(:)%exchmix =zero 375 Pawtab(:)%lamb_shielding = zero 376 377 ! * Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit. 378 ! TODO solve problem with memory leak and clean this part as well as the associated flag 379 call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab) 380 381 call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,Cryst%rprimd,Cryst%symrec,Pawang%zarot) 382 383 ! Initialize and compute data for DFT+U 384 Paw_dmft%use_dmft=Dtset%usedmft 385 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 386 is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Dtset%nspinor,Cryst%ntypat,Dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,& 387 Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu) 388 if (Dtset%usepawu>0.or.Dtset%useexexch>0) then 389 ABI_ERROR('BS equation with DFT+U not completely coded!') 390 end if 391 if (my_rank == master) call pawtab_print(Pawtab) 392 393 ! Get Pawrhoij from the header of the WFK file. 394 call pawrhoij_copy(Hdr_wfk%pawrhoij,KS_Pawrhoij) 395 396 ! Re-symmetrize rhoij === 397 ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly 398 choice=1; optrhoij=1 399 ! call pawrhoij_symrhoij(KS_Pawrhoij,KS_Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,& 400 ! & Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,& 401 ! & Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 402 403 ! Evaluate form factor of radial part of phi.phj-tphi.tphj === 404 rhoxsp_method=1 ! Arnaud-Alouani 405 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 406 407 ABI_MALLOC(gfft_osc,(3,nfftot_osc)) 408 call get_gftt(ngfft_osc,k0,gmet,gsq_osc,gfft_osc) 409 ABI_FREE(gfft_osc) 410 411 ! Set up q grids, make qmax 20% larger than largest expected: 412 ABI_MALLOC(nq_spl,(Psps%ntypat)) 413 ABI_MALLOC(qmax,(Psps%ntypat)) 414 nq_spl = Psps%mqgrid_ff 415 qmax = SQRT(gsq_osc)*1.2d0 ! qmax=Psps%qgrid_ff(Psps%mqgrid_ff) 416 ABI_MALLOC(Paw_pwff,(Psps%ntypat)) 417 418 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 419 420 ABI_FREE(nq_spl) 421 ABI_FREE(qmax) 422 ! 423 ! Variables/arrays related to the fine FFT grid === 424 ABI_MALLOC(ks_nhat,(nfftf,Dtset%nspden)) 425 ks_nhat=zero 426 ABI_MALLOC(Pawfgrtab,(Cryst%natom)) 427 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 428 call pawfgrtab_init(Pawfgrtab,cplex1,l_size_atm,Dtset%nspden,Dtset%typat) 429 ABI_FREE(l_size_atm) 430 compch_fft=greatest_real 431 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 432 ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation) 433 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 434 write(msg,'(a,i0)')' bethe_salpeter : using usexcnhat = ',usexcnhat 435 call wrtout(std_out,msg) 436 ! 437 ! Identify parts of the rectangular grid where the density has to be calculated === 438 optcut=0; optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm 439 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 440 441 call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 442 optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 443 else 444 ABI_MALLOC(Paw_pwff,(0)) 445 end if !End of PAW Initialization 446 447 ! Consistency check and additional stuff done only for GW with PAW. 448 if (Dtset%usepaw==1) then 449 if (Dtset%ecutwfn < Dtset%ecut) then 450 write(msg,"(5a)")& 451 "WARNING - ",ch10,& 452 " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 453 " an excessive truncation of the planewave basis set can lead to unphysical results." 454 call wrtout(ab_out,msg,'COLL') 455 end if 456 457 ABI_CHECK(Dtset%usedmft==0,"DMFT + BSE not allowed") 458 ABI_CHECK(Dtset%useexexch==0,"LEXX + BSE not allowed") 459 end if 460 461 ! Allocate these arrays anyway, since they are passed to subroutines. 462 if (.not.allocated(ks_nhat)) then 463 ABI_MALLOC(ks_nhat,(nfftf,0)) 464 end if 465 466 !================================================== 467 !==== Read KS band structure from the KSS file ==== 468 !================================================== 469 470 ! Initialize wave function handler, allocate wavefunctions. 471 my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds 472 ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol)) 473 nband=mband 474 475 !At present, no memory distribution, each node has the full set of states. 476 ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol)) 477 bks_mask=.TRUE. 478 479 ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol)) 480 keep_ur=.FALSE.; if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE. 481 482 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Dtset%nsppol,bks_mask,& 483 Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_osc,& 484 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 485 486 ABI_FREE(bks_mask) 487 ABI_FREE(nband) 488 ABI_FREE(keep_ur) 489 490 call wfd%print(header="Wavefunctions used to construct the e-h basis set",mode_paral='PERS') 491 492 call timab(651,2,tsec) ! bse(Init1) 493 call timab(653,1,tsec) ! bse(rdkss) 494 495 call wfd%read_wfk(wfk_fname,iomode_from_fname(wfk_fname)) 496 497 ! This test has been disabled (too expensive!) 498 if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 499 500 call timab(653,2,tsec) ! bse(rdkss) 501 call timab(655,1,tsec) ! bse(mkrho) 502 503 !TODO : check the consistency of Wfd with Wfd_dense !!! 504 if (BSp%use_interp) then 505 ! Initialize wave function handler, allocate wavefunctions. 506 my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds 507 ABI_MALLOC(nband,(Kmesh_dense%nibz,Dtset%nsppol)) 508 nband=mband 509 510 ! At present, no memory distribution, each node has the full set of states. 511 ! albeit we allocate only the states that are used. 512 ABI_MALLOC(bks_mask,(mband,Kmesh_dense%nibz,Dtset%nsppol)) 513 bks_mask=.False. 514 do spin=1,Bsp%nsppol 515 bks_mask(Bsp%lomo_spin(spin):Bsp%humo_spin(spin),:,spin) = .True. 516 end do 517 !bks_mask=.TRUE. 518 519 ABI_MALLOC(keep_ur,(mband,Kmesh_dense%nibz,Dtset%nsppol)) 520 keep_ur=.FALSE.; if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE. 521 522 call wfd_init(Wfd_dense,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh_dense%nibz,Dtset%nsppol,& 523 bks_mask,Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,ngfft_osc,& 524 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 525 526 ABI_FREE(bks_mask) 527 ABI_FREE(nband) 528 ABI_FREE(keep_ur) 529 530 call wfd_dense%print(header="Wavefunctions on the dense K-mesh used for interpolation",mode_paral='PERS') 531 532 call wfd_dense%read_wfk(Dtfil%fnameabi_wfkfine, iomode_from_fname(dtfil%fnameabi_wfkfine)) 533 !call wfd_dense%update_bkstab() 534 535 ! This test has been disabled (too expensive!) 536 if (.False.) call wfd_dense%test_ortho(Cryst,Pawtab,unit=std_out,mode_paral="COLL") 537 end if 538 539 !=== Calculate the FFT index of $(R^{-1}(r-\tau))$ === 540 !* S=\transpose R^{-1} and k_BZ = S k_IBZ 541 !* irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 542 ABI_MALLOC(irottb,(nfftot_osc,Cryst%nsym)) 543 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfft_osc,irottb,iscompatibleFFT) 544 545 ABI_MALLOC(ktabr,(nfftot_osc,Kmesh%nbz)) 546 do ik_bz=1,Kmesh%nbz 547 isym=Kmesh%tabo(ik_bz) 548 do ifft=1,nfftot_osc 549 ktabr(ifft,ik_bz)=irottb(ifft,isym) 550 end do 551 end do 552 ABI_FREE(irottb) 553 ! 554 !=========================== 555 !=== COMPUTE THE DENSITY === 556 !=========================== 557 !* Evaluate Planewave part (complete charge in case of NC pseudos). 558 ! 559 ABI_MALLOC(ks_rhor,(nfftf,Wfd%nspden)) 560 call wfd%mkrho(Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_rhor) 561 ! 562 !=== Additional computation for PAW === 563 nhatgrdim=0 564 if (Dtset%usepaw==1) then 565 ! 566 ! Calculate the compensation charge nhat. 567 if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc 568 ider=2*nhatgrdim; izero=0; qphon(:)=zero 569 if (nhatgrdim>0) then 570 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3)) 571 end if 572 573 call pawmknhat(compch_fft,cplex1,ider,idir0,ipert0,izero,Cryst%gprimd,& 574 Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 575 Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,qphon,Cryst%rprimd,& 576 Cryst%ucvol,Dtset%usewvl,Cryst%xred) 577 578 ! Evaluate onsite energies, potentials, densities === 579 ! * Initialize variables/arrays related to the PAW spheres. 580 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 581 ABI_MALLOC(KS_paw_ij,(Cryst%natom)) 582 call paw_ij_nullify(KS_paw_ij) 583 584 has_dijso=Dtset%pawspnorb 585 has_dijU=merge(0,1,Dtset%usepawu==0) 586 587 call paw_ij_init(KS_paw_ij,cplex1,Dtset%nspinor,Dtset%nsppol,& 588 Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 589 has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=0,has_dijxc_hat=0,has_dijxc_val=0,& 590 has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1) 591 592 ABI_MALLOC(KS_paw_an,(Cryst%natom)) 593 call paw_an_nullify(KS_paw_an) 594 595 nkxc1=0 596 call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,& 597 cplex1,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0) 598 599 ! Calculate onsite vxc with and without core charge === 600 nzlmopt=-1; option=0; compch_sph=greatest_real 601 602 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,& 603 Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 604 Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,& 605 Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,& 606 Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp) 607 end if !PAW 608 609 if (.not.allocated(ks_nhatgr)) then 610 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,0)) 611 end if 612 613 call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,& 614 Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 615 616 ! === For PAW, add the compensation charge on the FFT mesh, then get rho(G) === 617 if (Dtset%usepaw==1) ks_rhor(:,:)=ks_rhor(:,:)+ks_nhat(:,:) 618 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol) 619 620 ABI_MALLOC(ks_rhog,(2,nfftf)) 621 622 call fourdp(1,ks_rhog,ks_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp0) 623 call timab(655,2,tsec) ! bse(mkrho) 624 625 ! 626 ! The following steps have been gathered in the setvtr routine: 627 ! - get Ewald energy and Ewald forces 628 ! - compute local ionic pseudopotential vpsp 629 ! - eventually compute 3D core electron density xccc3d 630 ! - eventually compute vxc and vhartr 631 ! - set up ks_vtrial 632 ! 633 ! ******************************************************************* 634 ! **** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 635 ! ******************************************************************* 636 ngrvdw=0 637 ABI_MALLOC(grvdw,(3,ngrvdw)) 638 ABI_MALLOC(grchempottn,(3,Cryst%natom)) 639 ABI_MALLOC(grewtn,(3,Cryst%natom)) 640 nkxc=0 641 !if (Wfd%nspden==1) nkxc=2 642 !if (Wfd%nspden>=2) nkxc=3 ! check GGA and spinor, quite a messy part!!! 643 ABI_MALLOC(kxc,(nfftf,nkxc)) 644 645 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 646 ABI_MALLOC(xccc3d,(n3xccc)) 647 ABI_MALLOC(ks_vhartr,(nfftf)) 648 ABI_MALLOC(ks_vtrial,(nfftf,Wfd%nspden)) 649 ABI_MALLOC(vpsp,(nfftf)) 650 ABI_MALLOC(ks_vxc,(nfftf,Wfd%nspden)) 651 652 optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1 653 ! 654 !=== Compute structure factor phases and large sphere cut-off === 655 !WARNING cannot use Dtset%mgfft, this has to be checked better 656 !mgfft=MAXVAL(ngfftc(:)) 657 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom)) 658 write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf 659 if (Dtset%mgfftdg/=mgfftf) then 660 ! write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf" 661 ! write(std_out,*)'HACKING Dtset%mgfftf' 662 ! Dtset%mgfftdg=mgfftf 663 end if 664 ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom)) 665 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom)) 666 667 call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred) 668 669 if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then 670 call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred) 671 else 672 ph1df(:,:)=ph1d(:,:) 673 end if 674 675 ABI_FREE(ph1d) 676 677 call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 678 & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 679 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 680 & optene,Pawrad,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,& 681 & Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,wvl,xccc3d,Cryst%xred) 682 683 ABI_FREE(ph1df) 684 ABI_FREE(vpsp) 685 686 ! ============================ 687 ! ==== Compute KS PAW Dij ==== 688 ! ============================ 689 if (Wfd%usepaw==1) then 690 call timab(561,1,tsec) 691 ! 692 ! Calculate the unsymmetrized Dij. 693 call pawdij(cplex1,Dtset%enunit,Cryst%gprimd,ipert0,& 694 Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 695 Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,& 696 Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 697 k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,& 698 nucdipmom=Dtset%nucdipmom) 699 700 ! Symmetrize KS Dij 701 call symdij(Cryst%gprimd,Cryst%indsym,ipert0,& 702 Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,0,KS_paw_ij,Pawang,& 703 Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 704 705 ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 706 call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab) 707 call timab(561,2,tsec) 708 end if 709 710 ABI_FREE(kxc) 711 ABI_FREE(xccc3d) 712 ABI_FREE(grchempottn) 713 ABI_FREE(grewtn) 714 ABI_FREE(grvdw) 715 716 !=== QP_BSt stores energies and occ. used for the calculation === 717 !* Initialize QP_BSt with KS values. 718 !* In case of SC update QP_BSt using the QPS file. 719 ABI_MALLOC(qp_rhor,(nfftf,Dtset%nspden)) 720 qp_rhor=ks_rhor 721 722 ! AE density used for the model dielectric function. 723 ABI_MALLOC(qp_aerhor,(nfftf,Dtset%nspden)) 724 qp_aerhor=ks_rhor 725 726 ! PAW: Compute AE rhor. Under testing 727 if (Wfd%usepaw==1 .and. BSp%mdlf_type/=0) then 728 ABI_WARNING("Entering qp_aerhor with PAW") 729 730 ABI_MALLOC(qp_rhor_paw ,(nfftf,Wfd%nspden)) 731 ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden)) 732 ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden)) 733 734 ABI_MALLOC(qp_nhat,(nfftf,Wfd%nspden)) 735 qp_nhat = ks_nhat 736 ! TODO: I pass KS_pawrhoij instead of QP_pawrhoij but in the present version there's no difference. 737 738 call denfgr(Cryst%atindx1,Cryst%gmet,Wfd%comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,& 739 & Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_pawrhoij,Pawtab,Dtset%prtvol,& 740 & qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 741 742 norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3)) 743 write(msg,'(a,F8.4)') ' QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm 744 call wrtout(std_out,msg) 745 write(std_out,*)"MAX", MAXVAL(qp_rhor_paw(:,1)) 746 write(std_out,*)"MIN", MINVAL(qp_rhor_paw(:,1)) 747 748 ABI_FREE(qp_nhat) 749 ABI_FREE(qp_rhor_n_one) 750 ABI_FREE(qp_rhor_nt_one) 751 752 ! Use ae density for the model dielectric function. 753 iwarn=0 754 call mkdenpos(iwarn,nfftf,Wfd%nspden,option1,qp_rhor_paw,dtset%xc_denpos) 755 qp_aerhor = qp_rhor_paw 756 ABI_FREE(qp_rhor_paw) 757 end if 758 759 !call copy_bandstructure(KS_BSt,QP_BSt) 760 761 if (.FALSE.) then 762 ! $ m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}> $ 763 ABI_MALLOC(m_ks_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 764 m_ks_to_qp=czero 765 do spin=1,Wfd%nsppol 766 do ik_ibz=1,Wfd%nkibz 767 do band=1,Wfd%nband(ik_ibz,spin) 768 m_ks_to_qp(band,band,ik_ibz,spin)=cone ! Initialize the QP amplitudes with KS wavefunctions. 769 end do 770 end do 771 end do 772 ! 773 ! Now read m_ks_to_qp and update the energies in QP_BSt. 774 ! TODO switch on the renormalization of n in sigma. 775 ABI_MALLOC(prev_rhor,(nfftf,Wfd%nspden)) 776 ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Wfd%usepaw)) 777 778 call rdqps(QP_BSt,Dtfil%fnameabi_qps,Wfd%usepaw,Wfd%nspden,1,nscf,& 779 nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,prev_rhor,prev_Pawrhoij) 780 781 ABI_FREE(prev_rhor) 782 if (Psps%usepaw==1.and.nscf>0) then 783 call pawrhoij_free(prev_pawrhoij) 784 end if 785 ABI_FREE(prev_pawrhoij) 786 ! 787 ! if (nscf>0.and.wfd_iam_master(Wfd)) then ! Print the unitary transformation on std_out. 788 ! call show_QP(QP_BSt,m_ks_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp) 789 ! end if 790 ! 791 ! === Compute QP wfg as linear combination of KS states === 792 ! * Wfd%ug is modified inside calc_wf_qp 793 ! * For PAW, update also the on-site projections. 794 ! * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz 795 ! TODO here we should use nbsc instead of nbnds 796 797 call wfd%rotate(Cryst,m_ks_to_qp) 798 799 ABI_FREE(m_ks_to_qp) 800 ! 801 ! === Reinit the storage mode of Wfd as ug have been changed === 802 ! * Update also the wavefunctions for GW corrections on each processor 803 call wfd%reset_ur_cprj() 804 805 call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 806 807 ! Compute QP occupation numbers. 808 call wrtout(std_out,'bethe_salpeter: calculating QP occupation numbers') 809 810 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0) 811 ABI_MALLOC(qp_vbik,(QP_BSt%nkpt,QP_BSt%nsppol)) 812 qp_vbik(:,:) = ebands_get_valence_idx(QP_BSt) 813 ABI_FREE(qp_vbik) 814 815 call wfd%mkrho(Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_rhor) 816 end if 817 818 ABI_MALLOC(qp_rhog,(2,nfftf)) 819 call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,0) 820 821 ! States up to lomo-1 are useless now since only the states close to the gap are 822 ! needed to construct the EXC Hamiltonian. Here we deallocate the wavefunctions 823 ! to make room for the excitonic Hamiltonian that is allocated in exc_build_ham. 824 ! and for the screening that is allocated below. 825 ! Hereafter bands from 1 up to lomo-1 and all bands above humo+1 should not be accessed! 826 ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 827 828 bks_mask=.FALSE. 829 do spin=1,Bsp%nsppol 830 if (Bsp%lomo_spin(spin)>1) bks_mask(1:Bsp%lomo_spin(spin)-1,:,spin)=.TRUE. 831 if (Bsp%humo_spin(spin)+1<=Wfd%mband) bks_mask(Bsp%humo_spin(spin)+1:,:,spin)=.TRUE. 832 end do 833 call wfd%wave_free(what="All", bks_mask=bks_mask) 834 ABI_FREE(bks_mask) 835 ! 836 ! ================================================================ 837 ! Build the screened interaction W in the irreducible wedge. 838 ! * W(q,G1,G2) = vc^{1/2} (q,G1) e^{-1}(q,G1,G2) vc^{1/2) (q,G2) 839 ! * Use Coulomb term for q-->0, 840 ! * Only the first small Q is used, shall we average if nqlwl>1? 841 ! ================================================================ 842 ! TODO clean this part and add an option to retrieve a single frequency to save memory. 843 call timab(654,1,tsec) ! bse(rdmkeps^-1) 844 845 call screen_nullify(W) 846 if (BSp%use_coulomb_term) then ! Init W. 847 ! Incore or out-of-core solution? 848 mqmem=0; if (Dtset%gwmem/10==1) mqmem=Qmesh%nibz 849 850 W_info%invalid_freq = Dtset%gw_invalid_freq 851 W_info%mat_type = MAT_INV_EPSILON 852 !W_info%mat_type = MAT_W 853 !W_info%vtx_family 854 !W_info%ixc 855 !W_info%use_ada 856 W_info%use_mdf = BSp%mdlf_type 857 !W_info%use_ppm 858 !W_info%vtx_test 859 !W_info%wint_method 860 ! 861 !W_info%ada_kappa 862 W_info%eps_inf = BSp%eps_inf 863 !W_info%drude_plsmf 864 865 call screen_init(W,W_Info,Cryst,Qmesh,Gsph_c,Vcp,w_fname,mqmem,Dtset%npweps,& 866 Dtset%iomode,ngfftf,nfftf_tot,Wfd%nsppol,Wfd%nspden,qp_aerhor,Wfd%prtvol,Wfd%comm) 867 end if 868 call timab(654,2,tsec) ! bse(rdmkeps^-1) 869 ! 870 ! ============================================= 871 ! ==== Build of the excitonic Hamiltonian ===== 872 ! ============================================= 873 ! 874 call timab(656,1,tsec) ! bse(mkexcham) 875 876 call exc_build_ham(BSp,BS_files,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,& 877 & Wfd,W,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff) 878 ! 879 ! Free Em1 to make room for the full excitonic Hamiltonian. 880 call screen_free(W) 881 882 call timab(656,2,tsec) ! bse(mkexcham) 883 ! 884 ! ========================================= 885 ! ==== Macroscopic dielectric function ==== 886 ! ========================================= 887 call timab(657,1,tsec) ! bse(mkexceps) 888 ! 889 ! First deallocate the internal %ur buffers to make room for the excitonic Hamiltonian. 890 call timab(658,1,tsec) ! bse(wfd_wave_free) 891 ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 892 bks_mask=.TRUE. 893 call wfd%wave_free(what="Real_space", bks_mask=bks_mask) 894 ABI_FREE(bks_mask) 895 call timab(658,2,tsec) ! bse(wfd_wave_free) 896 897 ! Compute the commutator [r,Vu] (PAW only). 898 ABI_MALLOC(HUr,(Cryst%natom*Wfd%usepaw)) 899 900 call timab(659,1,tsec) ! bse(make_pawhur_t) 901 if (Bsp%inclvkb/=0 .and. Wfd%usepaw==1 .and. Dtset%usepawu/=0) then !TODO here I need KS_Paw_ij 902 ABI_WARNING("Commutator for DFT+U not tested") 903 call pawhur_init(hur,Wfd%nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,KS_Paw_ij) 904 end if 905 call timab(659,2,tsec) ! bse(make_pawhur_t) 906 907 select case (BSp%algorithm) 908 case (BSE_ALGO_NONE) 909 ABI_COMMENT("Skipping solution of the BSE equation") 910 911 case (BSE_ALGO_DDIAGO, BSE_ALGO_CG) 912 call timab(660,1,tsec) ! bse(exc_diago_driver) 913 call exc_diago_driver(Wfd,Bsp,BS_files,KS_BSt,QP_BSt,Cryst,Kmesh,Psps, Pawtab,Hur,Hdr_bse,drude_plsmf,Epren) 914 call timab(660,2,tsec) ! bse(exc_diago_driver) 915 916 if (.FALSE.) then ! Calculate electron-hole excited state density. Not tested at all. 917 call exc_den(BSp,BS_files,ngfftf,nfftf,Kmesh,ktabr,Wfd) 918 end if 919 920 if (.FALSE.) then 921 paw_add_onsite=.FALSE.; spin_opt=1; which_fixed=1; eh_rcoord=(/zero,zero,zero/); nrcell=(/2,2,2/) 922 ! call exc_plot(Bsp,Bs_files,Wfd,Kmesh,Cryst,Psps,Pawtab,Pawrad,paw_add_onsite,spin_opt,which_fixed,eh_rcoord,nrcell,ngfftf) 923 end if 924 925 if (BSp%use_interp) then 926 ABI_ERROR("Interpolation technique not coded for diagonalization and CG") 927 end if 928 929 case (BSE_ALGO_Haydock) 930 call timab(661,1,tsec) ! bse(exc_haydock_driver) 931 932 if (BSp%use_interp) then 933 934 call exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_BSt,Wfd,Psps,Pawtab,Hur,Epren, & 935 kmesh_dense=Kmesh_dense, ks_bst_dense=KS_BSt_dense, qp_bst_dense=QP_BSt_dense,wfd_dense=Wfd_dense, & 936 vcp_dense=Vcp_dense, grid=grid) 937 else 938 call exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_BSt,Wfd,Psps,Pawtab,Hur,Epren) 939 end if 940 941 call timab(661,2,tsec) ! bse(exc_haydock_driver) 942 943 case default 944 write(msg,'(a,i0)')" Wrong BSE algorithm: ",BSp%algorithm 945 ABI_ERROR(msg) 946 end select 947 948 call timab(657,2,tsec) ! bse(mkexceps) 949 950 !===================== 951 !==== Free memory ==== 952 !===================== 953 ABI_FREE(ktabr) 954 ABI_FREE(ks_vhartr) 955 ABI_FREE(ks_vtrial) 956 ABI_FREE(ks_vxc) 957 ABI_FREE(ks_nhat) 958 ABI_FREE(ks_nhatgr) 959 ABI_FREE(ks_rhog) 960 ABI_FREE(ks_rhor) 961 ABI_FREE(qp_rhog) 962 ABI_FREE(qp_rhor) 963 ABI_FREE(qp_aerhor) 964 ! 965 ! Free local data structures. 966 call destroy_mpi_enreg(MPI_enreg_seq) 967 call cryst%free() 968 call Gsph_x%free() 969 call Gsph_c%free() 970 call Kmesh%free() 971 call Qmesh%free() 972 call Hdr_wfk%free() 973 call Hdr_bse%free() 974 call ebands_free(KS_BSt) 975 call ebands_free(QP_BSt) 976 call Vcp%free() 977 call pawhur_free(Hur) 978 ABI_FREE(Hur) 979 call bs_parameters_free(BSp) 980 call wfd%free() 981 call pawfgr_destroy(Pawfgr) 982 call eprenorms_free(Epren) 983 984 ! Free memory used for interpolation. 985 if (BSp%use_interp) then 986 call double_grid_free(grid) 987 call wfd_dense%free() 988 call Gsph_x_dense%free() 989 call Gsph_c_dense%free() 990 call Kmesh_dense%free() 991 call Qmesh_dense%free() 992 call ebands_free(KS_BSt_dense) 993 call ebands_free(QP_BSt_dense) 994 call Vcp_dense%free() 995 call Hdr_wfk_dense%free() 996 end if 997 998 ! Optional deallocation for PAW. 999 if (Dtset%usepaw==1) then 1000 call pawrhoij_free(KS_Pawrhoij) 1001 ABI_FREE(KS_Pawrhoij) 1002 call pawfgrtab_free(Pawfgrtab) 1003 ABI_FREE(Pawfgrtab) 1004 call paw_ij_free(KS_paw_ij) 1005 ABI_FREE(KS_paw_ij) 1006 call paw_an_free(KS_paw_an) 1007 ABI_FREE(KS_paw_an) 1008 call pawpwff_free(Paw_pwff) 1009 end if 1010 ABI_FREE(Paw_pwff) 1011 1012 call timab(650,2,tsec) ! bse(Total) 1013 1014 DBG_EXIT('COLL') 1015 1016 end subroutine bethe_salpeter
m_bethe_salpeter/setup_bse [ Functions ]
[ Top ] [ m_bethe_salpeter ] [ Functions ]
NAME
setup_bse
FUNCTION
This routine performs the initialization of basic objects and quantities used for Bethe-Salpeter calculations. In particular the excparam data type that defines the parameters of the calculation is completely initialized starting from the content of Dtset and the parameters read from the external WFK and SCR (SUSC) file.
INPUTS
codvsn=Code version ngfft_gw(18)=Information about 3D FFT for density and potentials, see ~abinit/doc/variables/vargs.htm#ngfft acell(3)=Length scales of primitive translations (bohr) rprim(3,3)=Dimensionless real space primitive translations. Dtset<dataset_type>=All input variables for this dataset. Some of them might be redefined here TODO Dtfil=filenames and unit numbers used in abinit. Psps <pseudopotential_type>=variables related to pseudopotentials Pawtab(Psps%ntypat*Dtset%usepaw)<pawtab_type>=PAW tabulated starting data
OUTPUT
Cryst<crystal_t>=Info on the crystalline Structure. Kmesh<kmesh_t>=Structure defining the k-sampling for the wavefunctions. Qmesh<kmesh_t>=Structure defining the q-sampling for the symmetrized inverse dielectric matrix. Gsph_x<gsphere_t=Data type gathering info on the G-sphere for wave functions and e^{-1}, KS_BSt<ebands_t>=The KS band structure (energies, occupancies, k-weights...) Vcp<vcoul_t>=Structure gathering information on the Coulomb interaction in reciprocal space, including a possible cutoff in real space. ngfft_osc(18)=Contain all needed information about the 3D FFT for the oscillator matrix elements. See ~abinit/doc/variables/vargs.htm#ngfft Bsp<excparam>=Basic parameters defining the Bethe-Salpeter run. Completely initialed in output. Hdr_wfk<Hdr_type>=The header of the WFK file. Hdr_bse<Hdr_type>=Local header initialized from the parameters used for the Bethe-Salpeter calculation. BS_files<excfiles>=Files used in the calculation. w_file=File name used to construct W. Set to ABI_NOFILE if no external file is used.
SOURCE
1057 subroutine setup_bse(codvsn,acell,rprim,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,& 1058 & Cryst,Kmesh,Qmesh,KS_BSt,QP_bst,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,Wvl) 1059 1060 !Arguments ------------------------------------ 1061 !scalars 1062 integer,intent(in) :: comm 1063 character(len=8),intent(in) :: codvsn 1064 character(len=fnlen),intent(out) :: w_fname 1065 type(dataset_type),intent(inout) :: Dtset 1066 type(datafiles_type),intent(in) :: Dtfil 1067 type(pseudopotential_type),intent(in) :: Psps 1068 type(excparam),intent(inout) :: Bsp 1069 type(hdr_type),intent(out) :: Hdr_wfk,Hdr_bse 1070 type(crystal_t),intent(out) :: Cryst 1071 type(kmesh_t),intent(out) :: Kmesh,Qmesh 1072 type(gsphere_t),intent(out) :: Gsph_x,Gsph_c 1073 type(ebands_t),intent(out) :: KS_BSt,QP_Bst 1074 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw) 1075 type(vcoul_t),intent(out) :: Vcp 1076 type(excfiles),intent(out) :: BS_files 1077 type(wvl_internal_type), intent(in) :: Wvl 1078 type(eprenorms_t),intent(out) :: Epren 1079 !arrays 1080 integer,intent(out) :: ngfft_osc(18) 1081 real(dp),intent(in) :: acell(3),rprim(3,3) 1082 1083 !Local variables ------------------------------ 1084 !scalars 1085 integer,parameter :: pertcase0=0,master=0 1086 integer(i8b) :: work_size,tot_nreh,neh_per_proc,il 1087 integer :: bantot,enforce_sym,ib,ibtot,ik_ibz,isppol,jj,method,iat,ount !ii, 1088 integer :: mband,io,nfftot_osc,spin,hexc_size,nqlwl,iq 1089 integer :: timrev,iq_bz,isym,iq_ibz,itim 1090 integer :: my_rank,nprocs,fform,npwe_file,ierr,my_k1, my_k2,my_nbks 1091 integer :: first_dig,second_dig,it 1092 real(dp) :: ucvol,qnorm 1093 real(dp):: eff,mempercpu_mb,wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem 1094 logical,parameter :: remove_inv=.FALSE. 1095 logical :: ltest,occ_from_dtset 1096 character(len=500) :: msg 1097 character(len=fnlen) :: gw_fname,test_file,wfk_fname 1098 character(len=fnlen) :: ep_nc_fname 1099 type(hscr_t) :: Hscr 1100 !arrays 1101 integer :: ng0sh_opt(3),val_idx(Dtset%nsppol) 1102 integer,allocatable :: npwarr(:),val_indeces(:,:),nlmn_atm(:) 1103 real(dp) :: qpt_bz(3),minmax_tene(2) 1104 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3) 1105 real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:) 1106 real(dp),allocatable :: igwene(:,:,:) 1107 real(dp),pointer :: energies_p(:,:,:) 1108 complex(dpc),allocatable :: gw_energy(:,:,:) 1109 type(Pawrhoij_type),allocatable :: Pawrhoij(:) 1110 1111 !************************************************************************ 1112 1113 DBG_ENTER("COLL") 1114 1115 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1116 1117 ! === Check for calculations that are not implemented === 1118 ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1)) 1119 ABI_CHECK(ltest,'Dtset%nband must be constant') 1120 ABI_CHECK(Dtset%nspinor==1,"nspinor==2 not coded") 1121 1122 ! === Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume === 1123 call mkrdim(acell,rprim,rprimd) 1124 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1125 1126 ! Read energies and header from the WFK file. 1127 wfk_fname = dtfil%fnamewffk 1128 if (.not. file_exists(wfk_fname)) then 1129 wfk_fname = nctk_ncify(wfk_fname) 1130 ABI_COMMENT(sjoin("File not found. Will try netcdf file: ", wfk_fname)) 1131 end if 1132 1133 call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm) 1134 mband = MAXVAL(Hdr_wfk%nband) 1135 1136 call hdr_wfk%vs_dtset(dtset) 1137 1138 ! === Create crystal_t data type === 1139 !remove_inv= .FALSE. !(nsym_kss/=Hdr_wfk%nsym) 1140 timrev= 2 ! This information is not reported in the header 1141 ! 1 => do not use time-reversal symmetry 1142 ! 2 => take advantage of time-reversal symmetry 1143 1144 cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv) 1145 call cryst%print() 1146 ! 1147 ! Setup of the k-point list and symmetry tables in the BZ ----------------------------------- 1148 if (Dtset%chksymbreak==0) then 1149 ABI_WARNING("Calling make_mesh") 1150 call make_mesh(Kmesh,Cryst,Dtset%kptopt,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk,break_symmetry=.TRUE.) 1151 ! TODO 1152 !Check if kibz from KSS file corresponds to the one returned by make_mesh. 1153 else 1154 call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt) 1155 end if 1156 BSp%nkibz = Kmesh%nibz !We might allow for a smaller number of points.... 1157 1158 call Kmesh%print("K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 1159 call Kmesh%print("K-mesh for the wavefunctions",ab_out, 0, "COLL") 1160 1161 nqlwl = 0 1162 w_fname = ABI_NOFILE 1163 if (Dtset%getscr/=0 .or. Dtset%irdscr/=0 .or. dtset%getscr_filepath /= ABI_NOFILE) then 1164 w_fname=Dtfil%fnameabi_scr 1165 else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then 1166 w_fname=Dtfil%fnameabi_sus 1167 ABI_ERROR("(get|ird)suscep not implemented") 1168 end if 1169 1170 if (w_fname /= ABI_NOFILE) then ! Read dimensions from the external file 1171 1172 if (.not. file_exists(w_fname)) then 1173 w_fname = nctk_ncify(w_fname) 1174 ABI_COMMENT(sjoin("File not found. Will try netcdf file: ", w_fname)) 1175 end if 1176 1177 if (my_rank==master) then 1178 ! Master reads npw and nqlwl from SCR file. 1179 call wrtout(std_out,sjoin('Testing file: ', w_fname)) 1180 1181 call hscr%from_file(w_fname, fform, xmpi_comm_self) 1182 ! Echo the header. 1183 if (Dtset%prtvol>0) call Hscr%print() 1184 1185 npwe_file = Hscr%npwe ! Have to change %npweps if it was larger than dim on disk. 1186 nqlwl = Hscr%nqlwl 1187 1188 if (Dtset%npweps>npwe_file) then 1189 write(msg,'(2(a,i0),2a,i0)')& 1190 "The number of G-vectors stored on file (",npwe_file,") is smaller than Dtset%npweps = ",Dtset%npweps,ch10,& 1191 "Calculation will proceed with the maximum available set, npwe_file = ",npwe_file 1192 ABI_WARNING(msg) 1193 Dtset%npweps = npwe_file 1194 else if (Dtset%npweps<npwe_file) then 1195 write(msg,'(2(a,i0),2a,i0)')& 1196 "The number of G-vectors stored on file (",npwe_file,") is larger than Dtset%npweps = ",Dtset%npweps,ch10,& 1197 "Calculation will proceed with Dtset%npweps = ",Dtset%npweps 1198 ABI_COMMENT(msg) 1199 end if 1200 end if 1201 1202 call Hscr%bcast(master, my_rank, comm) 1203 call xmpi_bcast(Dtset%npweps,master,comm,ierr) 1204 call xmpi_bcast(nqlwl,master,comm,ierr) 1205 1206 if (nqlwl>0) then 1207 ABI_MALLOC(qlwl,(3,nqlwl)) 1208 qlwl = Hscr%qlwl 1209 end if 1210 ! 1211 ! Init Qmesh from the SCR file. 1212 call Qmesh%init(Cryst,Hscr%nqibz,Hscr%qibz,Dtset%kptopt) 1213 1214 ! The G-sphere for W and Sigma_c is initialized from the gvectors found in the SCR file. 1215 call Gsph_c%init(Cryst, Dtset%npweps, gvec=Hscr%gvec) 1216 1217 call Hscr%free() 1218 else 1219 ! Init Qmesh from the K-mesh reported in the WFK file. 1220 call find_qmesh(Qmesh,Cryst,Kmesh) 1221 1222 ! The G-sphere for W and Sigma_c is initialized from ecuteps. 1223 call Gsph_c%init(Cryst, 0, ecut=Dtset%ecuteps) 1224 Dtset%npweps = Gsph_c%ng 1225 end if 1226 1227 BSp%npweps = Dtset%npweps 1228 BSp%ecuteps = Dtset%ecuteps 1229 1230 if (nqlwl==0) then 1231 nqlwl=1 1232 ABI_MALLOC(qlwl,(3,nqlwl)) 1233 qlwl(:,nqlwl)= GW_Q0_DEFAULT 1234 write(msg,'(3a,i2,a,3f9.6)')& 1235 "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,& 1236 "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1) 1237 ABI_COMMENT(msg) 1238 end if 1239 write(std_out,*)"nqlwl and qlwl for Coulomb singularity and e^-1",nqlwl,qlwl 1240 1241 ! === Setup of q-mesh in the whole BZ === 1242 ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 1243 ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed) 1244 ! 1245 call Qmesh%print("Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL") 1246 call Qmesh%print("Q-mesh for the screening function",ab_out ,0 ,"COLL") 1247 1248 do iq_bz=1,Qmesh%nbz 1249 call qmesh%get_BZ_item(iq_bz,qpt_bz,iq_ibz,isym,itim) 1250 sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz)) 1251 if (ANY(ABS(Qmesh%bz(:,iq_bz)-sq )>1.0d-4)) then 1252 write(std_out,*) sq,Qmesh%bz(:,iq_bz) 1253 write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')& 1254 'qpoint ',Qmesh%bz(:,iq_bz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,& 1255 'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,& 1256 'however a non zero umklapp G_o vector is required and this is not yet allowed' 1257 ABI_ERROR(msg) 1258 end if 1259 end do 1260 1261 BSp%algorithm = Dtset%bs_algorithm 1262 BSp%nstates = Dtset%bs_nstates 1263 Bsp%nsppol = Dtset%nsppol 1264 Bsp%hayd_term = Dtset%bs_hayd_term 1265 1266 ! Define the algorithm for solving the BSE. 1267 if (BSp%algorithm == BSE_ALGO_HAYDOCK) then 1268 BSp%niter = Dtset%bs_haydock_niter 1269 BSp%haydock_tol = Dtset%bs_haydock_tol 1270 1271 else if (BSp%algorithm == BSE_ALGO_CG) then 1272 ! FIXME For the time being use an hardcoded value. 1273 ! TODO change name in Dtset% 1274 BSp%niter = Dtset%nstep !100 1275 BSp%cg_tolwfr = Dtset%tolwfr 1276 BSp%nline = Dtset%nline 1277 BSp%nbdbuf = Dtset%nbdbuf 1278 BSp%nstates = Dtset%bs_nstates 1279 ABI_WARNING("Check CG setup") 1280 else 1281 !BSp%niter = 0 1282 !BSp%tol_iter = HUGE(one) 1283 end if 1284 ! 1285 ! Shall we include Local field effects? 1286 SELECT CASE (Dtset%bs_exchange_term) 1287 CASE (0,1) 1288 BSp%exchange_term = Dtset%bs_exchange_term 1289 CASE DEFAULT 1290 write(msg,'(a,i0)')" Wrong bs_exchange_term: ",Dtset%bs_exchange_term 1291 ABI_ERROR(msg) 1292 END SELECT 1293 ! 1294 ! Treatment of the off-diagonal coupling block. 1295 SELECT CASE (Dtset%bs_coupling) 1296 CASE (0) 1297 BSp%use_coupling = 0 1298 msg = 'RESONANT ONLY CALCULATION' 1299 CASE (1) 1300 BSp%use_coupling = 1 1301 msg = ' RESONANT+COUPLING CALCULATION ' 1302 CASE DEFAULT 1303 write(msg,'(a,i0)')" Wrong bs_coupling: ",Dtset%bs_coupling 1304 ABI_ERROR(msg) 1305 END SELECT 1306 call wrtout(std_out,msg) 1307 1308 BSp%use_diagonal_Wgg = .FALSE. 1309 Bsp%use_coulomb_term = .TRUE. 1310 BSp%eps_inf=zero 1311 Bsp%mdlf_type=0 1312 1313 first_dig =MOD(Dtset%bs_coulomb_term,10) 1314 second_dig=Dtset%bs_coulomb_term/10 1315 1316 Bsp%wtype = second_dig 1317 SELECT CASE (second_dig) 1318 CASE (BSE_WTYPE_NONE) 1319 call wrtout(std_out,"Coulomb term won't be calculated") 1320 Bsp%use_coulomb_term = .FALSE. 1321 1322 CASE (BSE_WTYPE_FROM_SCR) 1323 call wrtout(std_out,"W is read from an external SCR file") 1324 Bsp%use_coulomb_term = .TRUE. 1325 1326 CASE (BSE_WTYPE_FROM_MDL) 1327 call wrtout(std_out,"W is approximated with the model dielectric function") 1328 Bsp%use_coulomb_term = .TRUE. 1329 BSp%mdlf_type = MDL_BECHSTEDT 1330 BSp%eps_inf = Dtset%mdf_epsinf 1331 ABI_CHECK(Bsp%eps_inf > zero, "mdf_epsinf <= 0") 1332 1333 CASE DEFAULT 1334 write(msg,'(a,i0)')" Wrong second digit in bs_coulomb_term: ",Dtset%bs_coulomb_term 1335 ABI_ERROR(msg) 1336 END SELECT 1337 ! 1338 ! Diagonal approximation or full matrix? 1339 BSp%use_diagonal_Wgg = .TRUE. 1340 if (Bsp%wtype /= BSE_WTYPE_NONE) then 1341 SELECT CASE (first_dig) 1342 CASE (0) 1343 call wrtout(std_out,"Using diagonal approximation W_GG") 1344 BSp%use_diagonal_Wgg = .TRUE. 1345 CASE (1) 1346 call wrtout(std_out,"Using full W_GG' matrix ") 1347 BSp%use_diagonal_Wgg = .FALSE. 1348 CASE DEFAULT 1349 write(msg,'(a,i0)')" Wrong first digit in bs_coulomb_term: ",Dtset%bs_coulomb_term 1350 ABI_ERROR(msg) 1351 END SELECT 1352 end if 1353 1354 !TODO move the initialization of the parameters for the interpolation in setup_bse_interp 1355 1356 BSp%use_interp = .FALSE. 1357 BSp%interp_mode = BSE_INTERP_YG 1358 BSp%interp_kmult(1:3) = 0 1359 BSp%prep_interp = .FALSE. 1360 BSp%sum_overlaps = .TRUE. ! Sum over the overlaps 1361 1362 ! Printing ncham 1363 BSp%prt_ncham = .FALSE. 1364 1365 ! Deactivate Interpolation Technique by default 1366 ! if (.FALSE.) then 1367 1368 ! Reading parameters from the input file 1369 BSp%use_interp = (dtset%bs_interp_mode /= 0) 1370 BSp%prep_interp = (dtset%bs_interp_prep == 1) 1371 1372 SELECT CASE (dtset%bs_interp_mode) 1373 CASE (0) 1374 ! No interpolation, do not print anything ! 1375 CASE (1) 1376 call wrtout(std_out,"Using interpolation technique with energies and wavefunctions from dense WFK") 1377 CASE (2) 1378 call wrtout(std_out,"Interpolation technique with energies and wfn on dense WFK + treatment ABC of divergence") 1379 CASE (3) 1380 call wrtout(std_out,"Interpolation technique + divergence ABC along diagonal") 1381 CASE (4) 1382 call wrtout(std_out,"Using interpolation technique mode 1 with full computation of hamiltonian") 1383 CASE DEFAULT 1384 ABI_ERROR(sjoin("Wrong interpolation mode for bs_interp_mode:", itoa(dtset%bs_interp_mode))) 1385 END SELECT 1386 1387 ! Read from dtset 1388 if(BSp%use_interp) then 1389 BSp%interp_method = dtset%bs_interp_method 1390 BSp%rl_nb = dtset%bs_interp_rl_nb 1391 BSp%interp_m3_width = dtset%bs_interp_m3_width 1392 BSp%interp_kmult(1:3) = dtset%bs_interp_kmult(1:3) 1393 BSp%interp_mode = dtset%bs_interp_mode 1394 end if 1395 1396 ! Dimensions and parameters of the calculation. 1397 ! TODO one should add npwx as well 1398 !BSp%npweps=Dtset%npweps 1399 !BSp%npwwfn=Dtset%npwwfn 1400 1401 ABI_MALLOC(Bsp%lomo_spin, (Bsp%nsppol)) 1402 ABI_MALLOC(Bsp%homo_spin, (Bsp%nsppol)) 1403 ABI_MALLOC(Bsp%lumo_spin, (Bsp%nsppol)) 1404 ABI_MALLOC(Bsp%humo_spin, (Bsp%nsppol)) 1405 ABI_MALLOC(Bsp%nbndv_spin, (Bsp%nsppol)) 1406 ABI_MALLOC(Bsp%nbndc_spin, (Bsp%nsppol)) 1407 1408 ! FIXME use bs_loband(nsppol) 1409 Bsp%lomo_spin = Dtset%bs_loband 1410 write(std_out,*)"bs_loband",Dtset%bs_loband 1411 !if (Bsp%nsppol == 2) Bsp%lomo_spin(2) = Dtset%bs_loband 1412 1413 ! Check lomo correct only for unpolarized semiconductors 1414 !if (Dtset%nsppol == 1 .and. Bsp%lomo > Dtset%nelect/2) then 1415 ! write(msg,'(a,i0,a,f8.3)') " Bsp%lomo = ",Bsp%lomo," cannot be greater than nelect/2 = ",Dtset%nelect/2 1416 ! ABI_ERROR(msg) 1417 !end if 1418 ! 1419 ! ============================================== 1420 ! ==== Setup of the q for the optical limit ==== 1421 ! ============================================== 1422 Bsp%inclvkb = Dtset%inclvkb 1423 1424 if (Dtset%gw_nqlwl == 0) then 1425 ! Predefined list of 6 q-versors (b vectors and Cart axis) 1426 call cryst%get_redcart_qdirs(Bsp%nq, Bsp%q) 1427 else 1428 BSp%nq = Dtset%gw_nqlwl 1429 ABI_MALLOC(BSp%q, (3,BSp%nq)) 1430 BSp%q = Dtset%gw_qlwl 1431 do iq=1,BSp%nq ! normalization 1432 qnorm = normv(BSp%q(:,iq), Cryst%gmet,"G") 1433 BSp%q(:,iq) = BSp%q(:,iq) / qnorm 1434 end do 1435 end if 1436 1437 ! ====================================================== 1438 ! === Define the flags defining the calculation type === 1439 ! ====================================================== 1440 Bsp%calc_type = Dtset%bs_calctype 1441 1442 BSp%mbpt_sciss = zero ! Shall we use the scissors operator to open the gap? 1443 if (ABS(Dtset%mbpt_sciss)>tol6) BSp%mbpt_sciss = Dtset%mbpt_sciss 1444 1445 !now test input parameters from input and WFK file and assume some defaults 1446 ! 1447 ! TODO Add the possibility of using a randomly shifted k-mesh with nsym>1. 1448 ! so that densities and potentials are correctly symmetrized but 1449 ! the list of the k-point in the IBZ is not expanded. 1450 1451 if (mband < Dtset%nband(1)) then 1452 write(msg,'(2(a,i0),3a,i0)')& 1453 'WFK file contains only ', mband,' levels instead of ',Dtset%nband(1),' required;',ch10,& 1454 'The calculation will be done with nbands= ',mband 1455 ABI_WARNING(msg) 1456 Dtset%nband(:) = mband 1457 end if 1458 1459 BSp%nbnds = Dtset%nband(1) ! TODO Note the change in the meaning of input variables 1460 1461 if (BSp%nbnds<=Dtset%nelect/2) then 1462 write(msg,'(2a,a,i0,a,f8.2)')& 1463 'BSp%nbnds cannot be smaller than homo ',ch10,& 1464 'while BSp%nbnds = ',BSp%nbnds,' and Dtset%nelect = ',Dtset%nelect 1465 ABI_ERROR(msg) 1466 end if 1467 1468 !TODO add new dim for exchange part and consider the possibility of having npwsigx > npwwfn (see setup_sigma). 1469 1470 ! === Build enlarged G-sphere for the exchange part === 1471 call Gsph_c%extend(Cryst, Dtset%ecutwfn, Gsph_x) 1472 call Gsph_x%print(unit=std_out,prtvol=Dtset%prtvol) 1473 1474 ! NPWVEC as the biggest between npweps and npwwfn. MG RECHECK this part. 1475 !BSp%npwwfn = Dtset%npwwfn 1476 Bsp%npwwfn = Gsph_x%ng ! FIXME temporary hack 1477 BSp%npwvec=MAX(BSp%npwwfn,BSp%npweps) 1478 Bsp%ecutwfn = Dtset%ecutwfn 1479 1480 ! Compute Coulomb term on the largest G-sphere. 1481 if (Gsph_x%ng > Gsph_c%ng ) then 1482 call Vcp%init(Gsph_x,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,& 1483 nqlwl,qlwl,comm) 1484 else 1485 call Vcp%init(Gsph_c,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,& 1486 nqlwl,qlwl,comm) 1487 end if 1488 1489 ABI_FREE(qlwl) 1490 1491 bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)) 1492 ABI_MALLOC(doccde,(bantot)) 1493 ABI_MALLOC(eigen,(bantot)) 1494 ABI_MALLOC(occfact,(bantot)) 1495 doccde=zero; eigen=zero; occfact=zero 1496 1497 ! Get occupation from input if occopt == 2 1498 occ_from_dtset = (Dtset%occopt == 2) 1499 1500 jj=0; ibtot=0 1501 do isppol=1,Dtset%nsppol 1502 do ik_ibz=1,Dtset%nkpt 1503 do ib=1,Hdr_wfk%nband(ik_ibz+(isppol-1)*Dtset%nkpt) 1504 ibtot=ibtot+1 1505 if (ib<=BSP%nbnds) then 1506 jj=jj+1 1507 eigen (jj)=energies_p(ib,ik_ibz,isppol) 1508 if (occ_from_dtset) then 1509 !Not occupations must be the same for different images 1510 occfact(jj)=Dtset%occ_orig(ibtot,1) 1511 else 1512 occfact(jj)=Hdr_wfk%occ(ibtot) 1513 end if 1514 end if 1515 end do 1516 end do 1517 end do 1518 1519 ABI_FREE(energies_p) 1520 ! 1521 ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of 1522 ! symmetry operations in the old GW code (symmorphy and inversion) 1523 ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6)) 1524 ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt') 1525 1526 ABI_MALLOC(npwarr,(Dtset%nkpt)) 1527 npwarr=BSP%npwwfn 1528 1529 call ebands_init(bantot,KS_BSt,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& 1530 doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,& 1531 Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,& 1532 dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1533 dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1534 1535 ABI_FREE(doccde) 1536 ABI_FREE(eigen) 1537 ABI_FREE(npwarr) 1538 1539 !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when 1540 ! the occupation scheme for semiconductors is used. 1541 call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 1542 1543 call ebands_print(KS_BSt,"Band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol) 1544 1545 call ebands_report_gap(KS_BSt,header=" KS band structure",unit=std_out,mode_paral="COLL") 1546 1547 ABI_MALLOC(val_indeces,(KS_BSt%nkpt,KS_BSt%nsppol)) 1548 val_indeces = ebands_get_valence_idx(KS_BSt) 1549 1550 do spin=1,KS_BSt%nsppol 1551 val_idx(spin) = val_indeces(1,spin) 1552 write(msg,'(a,i2,a,i0)')" For spin : ",spin," val_idx ",val_idx(spin) 1553 call wrtout(std_out,msg) 1554 if ( ANY(val_indeces(1,spin) /= val_indeces(:,spin)) ) then 1555 ABI_ERROR("BSE code does not support metals") 1556 end if 1557 end do 1558 1559 ABI_FREE(val_indeces) 1560 ! 1561 ! === Create the BSE header === 1562 call hdr_init(KS_BSt,codvsn,Dtset,Hdr_bse,Pawtab,pertcase0,Psps,wvl) 1563 1564 ! === Get Pawrhoij from the header of the WFK file === 1565 ABI_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw)) 1566 if (Dtset%usepaw==1) then 1567 call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 1568 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 1569 end if 1570 1571 ! CP modified 1572 !call hdr_bse%update(bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 1573 call hdr_bse%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 1574 ! End CP modified 1575 1576 ABI_FREE(occfact) 1577 1578 if (Dtset%usepaw==1) call pawrhoij_free(Pawrhoij) 1579 ABI_FREE(Pawrhoij) 1580 1581 ! Find optimal value for G-sphere enlargment due to oscillator matrix elements 1582 ! We will split k-points over processors 1583 call xmpi_split_work(Kmesh%nbz, comm, my_k1, my_k2) 1584 1585 ! If there is no work to do, just skip the computation 1586 if (my_k2-my_k1+1 <= 0) then 1587 ng0sh_opt(:)=(/zero,zero,zero/) 1588 else 1589 ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true since bz is buggy 1590 ! * -one is used because we loop over all the possibile differences, unlike screening 1591 call get_ng0sh(my_k2-my_k1+1,Kmesh%bz(:,my_k1:my_k2),Kmesh%nbz,Kmesh%bz,& 1592 & Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt) 1593 end if 1594 1595 call xmpi_max(ng0sh_opt,BSp%mg0,comm,ierr) 1596 1597 write(msg,'(a,3(i0,1x))') ' optimal value for ng0sh = ',BSp%mg0 1598 call wrtout(std_out,msg) 1599 1600 ! === Setup of the FFT mesh for the oscilator strengths === 1601 ! * ngfft_osc(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening. 1602 ! * Here we redefine ngfft_osc(1:6) according to the following options : 1603 ! 1604 ! method==0 --> FFT grid read from fft.in (debugging purpose) 1605 ! method==1 --> Normal FFT mesh 1606 ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 1607 ! method==3 --> Doubled FFT grid, same as the the FFT for the density, 1608 ! 1609 ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library 1610 ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries 1611 ! 1612 ngfft_osc(1:18)=Dtset%ngfft(1:18); method=2 1613 if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0 1614 if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1 1615 if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2 1616 if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3 1617 enforce_sym=MOD(Dtset%fftgw,10) 1618 1619 call setmesh(gmet,Gsph_x%gvec,ngfft_osc,BSp%npwvec,BSp%npweps,BSp%npwwfn,nfftot_osc,method,BSp%mg0,Cryst,enforce_sym) 1620 nfftot_osc=PRODUCT(ngfft_osc(1:3)) 1621 1622 call print_ngfft(ngfft_osc,"FFT mesh for oscillator matrix elements",std_out,"COLL",prtvol=Dtset%prtvol) 1623 ! 1624 ! BSp%homo gives the 1625 !BSp%homo = val_idx(1) 1626 ! highest occupied band for each spin 1627 BSp%homo_spin = val_idx 1628 1629 ! TODO generalize the code to account for this unlikely case. 1630 !if (Dtset%nsppol==2) then 1631 ! ABI_CHECK(BSp%homo == val_idx(2),"Different valence indeces for spin up and down") 1632 !end if 1633 1634 !BSp%lumo = BSp%homo + 1 1635 !BSp%humo = BSp%nbnds 1636 !BSp%nbndv = BSp%homo - BSp%lomo + 1 1637 !BSp%nbndc = BSp%nbnds - BSp%homo 1638 1639 BSp%lumo_spin = BSp%homo_spin + 1 1640 BSp%humo_spin = BSp%nbnds 1641 BSp%nbndv_spin = BSp%homo_spin - BSp%lomo_spin + 1 1642 BSp%nbndc_spin = BSp%nbnds - BSp%homo_spin 1643 BSp%maxnbndv = MAXVAL(BSp%nbndv_spin(:)) 1644 BSp%maxnbndc = MAXVAL(BSp%nbndc_spin(:)) 1645 1646 BSp%nkbz = Kmesh%nbz 1647 1648 call ebands_copy(KS_BSt,QP_bst) 1649 ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol)) 1650 igwene=zero 1651 1652 call bsp_calctype2str(Bsp,msg) 1653 call wrtout(std_out,"Calculation type: "//TRIM(msg)) 1654 1655 SELECT CASE (Bsp%calc_type) 1656 CASE (BSE_HTYPE_RPA_KS) 1657 if (ABS(BSp%mbpt_sciss)>tol6) then 1658 write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies." 1659 call wrtout(std_out,msg) 1660 call ebands_apply_scissors(QP_BSt,BSp%mbpt_sciss) 1661 else 1662 write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]." 1663 call wrtout(std_out,msg) 1664 end if 1665 1666 CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw 1667 gw_fname=TRIM(Dtfil%filnam_ds(4))//'_GW' 1668 gw_fname="__in.gw__" 1669 if (.not.file_exists(gw_fname)) then 1670 msg = " File "//TRIM(gw_fname)//" not found. Aborting now" 1671 ABI_ERROR(msg) 1672 end if 1673 1674 call rdgw(QP_Bst,gw_fname,igwene,extrapolate=.FALSE.) ! here gwenergy is real 1675 1676 do isppol=1,Dtset%nsppol 1677 write(std_out,*) ' k GW energies [eV]' 1678 do ik_ibz=1,Kmesh%nibz 1679 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(QP_bst%eig(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds) 1680 end do 1681 write(std_out,*) ' k Im GW energies [eV]' 1682 do ik_ibz=1,Kmesh%nibz 1683 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(igwene(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds) 1684 end do 1685 end do 1686 ! 1687 ! If required apply the scissors operator on top of the QP bands structure (!) 1688 if (ABS(BSp%mbpt_sciss)>tol6) then 1689 write(msg,'(a,f8.2,a)')' Applying a scissors operator ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the QP energies!" 1690 ABI_COMMENT(msg) 1691 call ebands_apply_scissors(QP_BSt,BSp%mbpt_sciss) 1692 end if 1693 1694 CASE (BSE_HTYPE_RPA_QP) 1695 ABI_ERROR("Not implemented error!") 1696 1697 CASE DEFAULT 1698 write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type 1699 ABI_ERROR(msg) 1700 END SELECT 1701 1702 call ebands_report_gap(QP_BSt,header=" QP band structure",unit=std_out,mode_paral="COLL") 1703 1704 ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index. 1705 ! FIXME: linewidths not coded. 1706 ABI_MALLOC(gw_energy,(BSp%nbnds,Kmesh%nibz,Dtset%nsppol)) 1707 gw_energy = QP_BSt%eig 1708 1709 BSp%have_complex_ene = ANY(igwene > tol16) 1710 1711 ! Compute the number of resonant transitions, nreh, for the two spin channels and initialize BSp%Trans. 1712 ABI_MALLOC(Bsp%nreh,(Bsp%nsppol)) 1713 1714 ! Possible cutoff on the transitions. 1715 BSp%ircut = Dtset%bs_eh_cutoff(1) 1716 BSp%uvcut = Dtset%bs_eh_cutoff(2) 1717 1718 call init_transitions(BSp%Trans,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz,Bsp%nbnds,Bsp%nkibz,& 1719 & BSp%nsppol,Dtset%nspinor,gw_energy,QP_BSt%occ,Kmesh%tab,minmax_tene,Bsp%nreh) 1720 1721 ! Setup of the frequency mesh for the absorption spectrum. 1722 ! If not specified, use the min-max resonant transition energy and make it 10% smaller|larger. 1723 1724 !if (ABS(Dtset%bs_freq_mesh(1)) < tol6) then 1725 ! Dtset%bs_freq_mesh(1) = MAX(minmax_tene(1) - minmax_tene(1) * 0.1, zero) 1726 !end if 1727 1728 if (ABS(Dtset%bs_freq_mesh(2)) < tol6) then 1729 Dtset%bs_freq_mesh(2) = minmax_tene(2) + minmax_tene(2) * 0.1 1730 end if 1731 1732 Bsp%omegai = Dtset%bs_freq_mesh(1) 1733 Bsp%omegae = Dtset%bs_freq_mesh(2) 1734 Bsp%domega = Dtset%bs_freq_mesh(3) 1735 BSp%broad = Dtset%zcut 1736 1737 ! The frequency mesh (including the complex imaginary shift) 1738 BSp%nomega = (BSp%omegae - BSp%omegai)/BSp%domega + 1 1739 ABI_MALLOC(BSp%omega,(BSp%nomega)) 1740 do io=1,BSp%nomega 1741 BSp%omega(io) = (BSp%omegai + (io-1)*BSp%domega) + j_dpc*BSp%broad 1742 end do 1743 1744 ABI_FREE(gw_energy) 1745 ABI_FREE(igwene) 1746 1747 do spin=1,Bsp%nsppol 1748 write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh(spin) 1749 call wrtout(std_out,msg) 1750 end do 1751 1752 if (ANY(Bsp%nreh/=Bsp%nreh(1))) then 1753 write(msg,'(a,2(i0,1x))')" BSE code with different number of transitions for the two spin channels: ",Bsp%nreh 1754 ABI_WARNING(msg) 1755 end if 1756 ! 1757 ! Create transition table vcks2t 1758 Bsp%lomo_min = MINVAL(BSp%lomo_spin) 1759 Bsp%homo_max = MAXVAL(BSp%homo_spin) 1760 Bsp%lumo_min = MINVAL(BSp%lumo_spin) 1761 Bsp%humo_max = MAXVAL(BSp%humo_spin) 1762 1763 ABI_MALLOC(Bsp%vcks2t,(BSp%lomo_min:BSp%homo_max,BSp%lumo_min:BSp%humo_max,BSp%nkbz,Dtset%nsppol)) 1764 Bsp%vcks2t = 0 1765 1766 do spin=1,BSp%nsppol 1767 do it=1,BSp%nreh(spin) 1768 BSp%vcks2t(BSp%Trans(it,spin)%v,BSp%Trans(it,spin)%c,BSp%Trans(it,spin)%k,spin) = it 1769 end do 1770 end do 1771 1772 hexc_size = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hexc_size=2*hexc_size 1773 if (Bsp%nstates<=0) then 1774 Bsp%nstates=hexc_size 1775 else 1776 if (Bsp%nstates>hexc_size) then 1777 Bsp%nstates=hexc_size 1778 write(msg,'(2(a,i0),2a)')& 1779 "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates,ch10,& 1780 "the number of excitonic states nstates has been modified" 1781 ABI_WARNING(msg) 1782 end if 1783 end if 1784 1785 msg=' Fundamental parameters for the solution of the Bethe-Salpeter equation:' 1786 call print_bs_parameters(BSp,unit=std_out,header=msg,mode_paral="COLL",prtvol=Dtset%prtvol) 1787 call print_bs_parameters(BSp,unit=ab_out, header=msg,mode_paral="COLL") 1788 1789 if (ANY (Cryst%symrec(:,:,1) /= RESHAPE ( (/1,0,0,0,1,0,0,0,1/),(/3,3/) )) .or. ANY( ABS(Cryst%tnons(:,1)) > tol6) ) then 1790 write(msg,'(3a,9i2,2a,3f6.3,2a)')& 1791 "The first symmetry operation should be the Identity with zero tnons while ",ch10,& 1792 "symrec(:,:,1) = ",Cryst%symrec(:,:,1),ch10,& 1793 "tnons(:,1) = ",Cryst%tnons(:,1),ch10,& 1794 "This is not allowed, sym_rhotwgq0 should be changed." 1795 ABI_ERROR(msg) 1796 end if 1797 ! 1798 ! Prefix for generic output files. 1799 BS_files%out_basename = TRIM(Dtfil%filnam_ds(4)) 1800 ! 1801 ! Search for files to restart from. 1802 if (Dtset%gethaydock/=0 .or. Dtset%irdhaydock/=0) then 1803 BS_files%in_haydock_basename = TRIM(Dtfil%fnameabi_haydock) 1804 end if 1805 1806 test_file = Dtfil%fnameabi_bsham_reso 1807 if (file_exists(test_file)) then 1808 BS_files%in_hreso = test_file 1809 else 1810 BS_files%out_hreso = TRIM(Dtfil%filnam_ds(4))//'_BSR' 1811 end if 1812 1813 test_file = Dtfil%fnameabi_bsham_coup 1814 if (file_exists(test_file) ) then 1815 BS_files%in_hcoup = test_file 1816 else 1817 BS_files%out_hcoup = TRIM(Dtfil%filnam_ds(4))//'_BSC' 1818 end if 1819 ! 1820 ! in_eig is the name of the input file with eigenvalues and eigenvectors 1821 ! constructed from getbseig or irdbseig. out_eig is the name of the output file 1822 ! produced by this dataset. in_eig_exists checks for the presence of the input file. 1823 ! 1824 if (file_exists(Dtfil%fnameabi_bseig)) then 1825 BS_files%in_eig = Dtfil%fnameabi_bseig 1826 else 1827 BS_files%out_eig = TRIM(BS_files%out_basename)//"_BSEIG" 1828 end if 1829 1830 call print_bs_files(BS_files,unit=std_out) 1831 ! 1832 ! ========================================================== 1833 ! ==== Temperature dependence of the spectrum ============== 1834 ! ========================================================== 1835 BSp%do_ep_renorm = .FALSE. 1836 BSp%do_lifetime = .FALSE. ! Not yet implemented 1837 1838 ep_nc_fname = 'test_EP.nc' 1839 if(file_exists(ep_nc_fname)) then 1840 BSp%do_ep_renorm = .TRUE. 1841 1842 if(my_rank == master) then 1843 call eprenorms_from_epnc(Epren,ep_nc_fname) 1844 end if 1845 call eprenorms_bcast(Epren,master,comm) 1846 end if 1847 ! 1848 ! ========================================================== 1849 ! ==== Final check on the parameters of the calculation ==== 1850 ! ========================================================== 1851 if ( Bsp%use_coupling>0 .and. ALL(Bsp%algorithm /= [BSE_ALGO_DDIAGO, BSE_ALGO_HAYDOCK]) ) then 1852 ABI_ERROR("Resonant+Coupling is only available with the direct diagonalization or the haydock method.") 1853 end if 1854 1855 ! autoparal section 1856 if (dtset%max_ncpus /=0 .and. dtset%autoparal /=0 ) then 1857 ount = ab_out 1858 ! TODO: 1859 ! nsppol and calculation with coupling! 1860 1861 ! Temporary table needed to estimate memory 1862 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 1863 if (Dtset%usepaw==1) then 1864 do iat=1,Cryst%natom 1865 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 1866 end do 1867 end if 1868 1869 tot_nreh = SUM(BSp%nreh) 1870 work_size = tot_nreh * (tot_nreh + 1) / 2 1871 1872 write(ount,'(a)')"--- !Autoparal" 1873 write(ount,"(a)")'#Autoparal section for Bethe-Salpeter runs.' 1874 1875 write(ount,"(a)") "info:" 1876 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 1877 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 1878 write(ount,"(a,i0)")" nkibz: ",Bsp%nkibz 1879 write(ount,"(a,i0)")" nkbz: ",Bsp%nkbz 1880 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 1881 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 1882 write(ount,"(a,i0)")" lomo_min: ",Bsp%lomo_min 1883 write(ount,"(a,i0)")" humo_max: ",Bsp%humo_max 1884 write(ount,"(a,i0)")" tot_nreh: ",tot_nreh 1885 !write(ount,"(a,i0)")" nbnds: ",Ep%nbnds 1886 1887 ! Wavefunctions are not distributed. We read all the bands 1888 ! from 1 up to Bsp%nbnds because we have to recompute rhor 1889 ! but then we deallocate all the states that are not used for the construction of the e-h 1890 ! before allocating the EXC hamiltonian. Hence we can safely use (humo - lomo + 1) instead of Bsp%nbnds. 1891 !my_nbks = (Bsp%humo - Bsp%lomo +1) * Bsp%nkibz * Dtset%nsppol 1892 1893 ! This one overestimates the memory but it seems to be safer. 1894 my_nbks = Bsp%nbnds * Dtset%nkpt * Dtset%nsppol 1895 1896 ! Memory needed for Fourier components ug. 1897 ug_mem = two*gwpc*Dtset%nspinor*Bsp%npwwfn*my_nbks*b2Mb 1898 1899 ! Memory needed for real space ur. 1900 ur_mem = zero 1901 if (MODULO(Dtset%gwmem,10)==1) then 1902 ur_mem = two*gwpc*Dtset%nspinor*nfftot_osc*my_nbks*b2Mb 1903 end if 1904 1905 ! Memory needed for PAW projections Cprj 1906 cprj_mem = zero 1907 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 1908 1909 wfsmem_mb = ug_mem + ur_mem + cprj_mem 1910 1911 ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI: wavefunctions + W 1912 nonscal_mem = (wfsmem_mb + two*gwpc*BSp%npweps**2*b2Mb) * 1.1_dp 1913 1914 ! List of configurations. 1915 write(ount,"(a)")"configurations:" 1916 do il=1,dtset%max_ncpus 1917 if (il > work_size) cycle 1918 neh_per_proc = work_size / il 1919 neh_per_proc = neh_per_proc + MOD(work_size, il) 1920 eff = (one * work_size) / (il * neh_per_proc) 1921 1922 ! EXC matrix is distributed. 1923 mempercpu_mb = nonscal_mem + two * dpc * neh_per_proc * b2Mb 1924 1925 write(ount,"(a,i0)")" - tot_ncpus: ",il 1926 write(ount,"(a,i0)")" mpi_ncpus: ",il 1927 !write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 1928 write(ount,"(a,f12.9)")" efficiency: ",eff 1929 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 1930 end do 1931 1932 write(ount,'(a)')"..." 1933 1934 ABI_FREE(nlmn_atm) 1935 ABI_ERROR_NODUMP("aborting now") 1936 end if 1937 1938 DBG_EXIT("COLL") 1939 1940 end subroutine setup_bse
m_bethe_salpeter/setup_bse_interp [ Functions ]
[ Top ] [ m_bethe_salpeter ] [ Functions ]
NAME
setup_bse_interp
FUNCTION
INPUTS
ngfft_gw(18)=Information about 3D FFT for density and potentials, see ~abinit/doc/variables/vargs.htm#ngfft acell(3)=Length scales of primitive translations (bohr) rprim(3,3)=Dimensionless real space primitive translations. Dtset<dataset_type>=All input variables for this dataset. Some of them might be redefined here TODO Dtfil=filenames and unit numbers used in abinit. fnameabi_wfkfile is changed is Fortran file is not found but a netcdf version with similar name is available.
OUTPUT
Cryst<crystal_structure>=Info on the crystalline Structure. Kmesh<BZ_mesh_type>=Structure defining the k-sampling for the wavefunctions. Qmesh<BZ_mesh_type>=Structure defining the q-sampling for the symmetrized inverse dielectric matrix. Gsph_x<gsphere_t=Data type gathering info on the G-sphere for wave functions and e^{-1}, KS_BSt<Bandstructure_type>=The KS band structure (energies, occupancies, k-weights...) Vcp<vcoul_t>=Structure gathering information on the Coulomb interaction in reciprocal space, including a possible cutoff in real space. ngfft_osc(18)=Contain all needed information about the 3D FFT for the oscillator matrix elements. See ~abinit/doc/variables/vargs.htm#ngfft Bsp<excparam>=Basic parameters defining the Bethe-Salpeter run. Completely initialed in output. Hdr_wfk<Hdr_type>=The header of the WFK file. Hdr_bse<Hdr_type>=Local header initialized from the parameters used for the Bethe-Salpeter calculation. w_file=File name used to construct W. Set to ABI_NOFILE if no external file is used.
SOURCE
1975 subroutine setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh, & 1976 Kmesh_dense,Qmesh_dense,KS_BSt_dense,QP_bst_dense,Gsph_x,Gsph_c,Vcp_dense,Hdr_wfk_dense,grid,comm) 1977 1978 !Arguments ------------------------------------ 1979 !scalars 1980 integer,intent(in) :: comm 1981 type(dataset_type),intent(in) :: Dtset 1982 type(datafiles_type),intent(inout) :: Dtfil 1983 type(excparam),intent(inout) :: Bsp 1984 type(hdr_type),intent(out) :: Hdr_wfk_dense 1985 type(crystal_t),intent(in) :: Cryst 1986 type(kmesh_t),intent(in) :: Kmesh 1987 type(kmesh_t),intent(out) :: Kmesh_dense,Qmesh_dense 1988 type(ebands_t),intent(out) :: KS_BSt_dense,QP_Bst_dense 1989 type(double_grid_t),intent(out) :: grid 1990 type(vcoul_t),intent(out) :: Vcp_dense 1991 type(gsphere_t),intent(out) :: Gsph_x,Gsph_c 1992 !arrays 1993 1994 !Local variables ------------------------------ 1995 !scalars 1996 integer,parameter :: pertcase0=0,master=0 1997 integer :: bantot_dense,ib,ibtot,ik_ibz,isppol,jj 1998 integer :: nbnds_kss_dense 1999 integer :: spin,hexc_size 2000 integer :: my_rank 2001 integer :: it 2002 integer :: nprocs 2003 integer :: is1,is2,is3,is4 2004 real(dp) :: nelect_hdr_dense 2005 logical,parameter :: remove_inv=.FALSE. 2006 character(len=500) :: msg 2007 character(len=fnlen) :: wfk_fname_dense 2008 integer :: nqlwl 2009 !arrays 2010 integer,allocatable :: npwarr(:) 2011 real(dp),allocatable :: shiftk(:,:) 2012 real(dp),allocatable :: doccde(:),eigen(:),occfact(:) 2013 real(dp),pointer :: energies_p_dense(:,:,:) 2014 complex(dpc),allocatable :: gw_energy(:,:,:) 2015 integer,allocatable :: nbands_temp(:) 2016 integer :: kptrlatt_dense(3,3) 2017 real(dp),allocatable :: qlwl(:,:) 2018 real(dp) :: minmax_tene(2) 2019 2020 !************************************************************************ 2021 2022 DBG_ENTER("COLL") 2023 2024 kptrlatt_dense = zero 2025 2026 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 2027 2028 SELECT CASE(BSp%interp_mode) 2029 CASE (1,2,3,4) 2030 nbnds_kss_dense = -1 2031 wfk_fname_dense = Dtfil%fnameabi_wfkfine 2032 call wrtout(std_out,"BSE Interpolation: will read energies from: "//trim(wfk_fname_dense),"COLL") 2033 2034 if (nctk_try_fort_or_ncfile(wfk_fname_dense, msg) /= 0) then 2035 ABI_ERROR(msg) 2036 end if 2037 2038 Dtfil%fnameabi_wfkfine = wfk_fname_dense 2039 2040 call wfk_read_eigenvalues(wfk_fname_dense,energies_p_dense,Hdr_wfk_dense,comm) 2041 nbnds_kss_dense = MAXVAL(Hdr_wfk_dense%nband) 2042 CASE DEFAULT 2043 ABI_ERROR("Not yet implemented") 2044 END SELECT 2045 2046 nelect_hdr_dense = Hdr_wfk_dense%nelect 2047 2048 if (ABS(Dtset%nelect-nelect_hdr_dense)>tol6) then 2049 write(msg,'(2(a,f8.2))')& 2050 & "File contains ", nelect_hdr_dense," electrons but nelect initialized from input is ",Dtset%nelect 2051 ABI_ERROR(msg) 2052 end if 2053 2054 ! Setup of the k-point list and symmetry tables in the BZ 2055 SELECT CASE(BSp%interp_mode) 2056 CASE (1,2,3,4) 2057 if(Dtset%chksymbreak == 0) then 2058 ABI_MALLOC(shiftk,(3,Dtset%nshiftk)) 2059 kptrlatt_dense(:,1) = BSp%interp_kmult(1)*Dtset%kptrlatt(:,1) 2060 kptrlatt_dense(:,2) = BSp%interp_kmult(2)*Dtset%kptrlatt(:,2) 2061 kptrlatt_dense(:,3) = BSp%interp_kmult(3)*Dtset%kptrlatt(:,3) 2062 do jj = 1,Dtset%nshiftk 2063 shiftk(:,jj) = Bsp%interp_kmult(:)*Dtset%shiftk(:,jj) 2064 end do 2065 call make_mesh(Kmesh_dense,Cryst,Dtset%kptopt,kptrlatt_dense,Dtset%nshiftk,shiftk,break_symmetry=.TRUE.) 2066 ABI_FREE(shiftk) 2067 else 2068 !Initialize Kmesh with no wrapping inside ]-0.5;0.5] 2069 call Kmesh_dense%init(Cryst,Hdr_wfk_dense%nkpt,Hdr_wfk_dense%kptns,Dtset%kptopt) 2070 end if 2071 CASE DEFAULT 2072 ABI_ERROR("Not yet implemented") 2073 END SELECT 2074 2075 ! Init Qmesh 2076 call find_qmesh(Qmesh_dense,Cryst,Kmesh_dense) 2077 2078 call Gsph_c%init(Cryst, 0, ecut=Dtset%ecuteps) 2079 2080 call double_grid_init(Kmesh,Kmesh_dense,Dtset%kptrlatt,BSp%interp_kmult,grid) 2081 2082 BSp%nkibz_interp = Kmesh_dense%nibz !We might allow for a smaller number of points.... 2083 2084 call Kmesh_dense%print("Interpolated K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 2085 call Kmesh_dense%print("Interpolated K-mesh for the wavefunctions",ab_out, 0, "COLL") 2086 2087 if (nbnds_kss_dense < Dtset%nband(1)) then 2088 write(msg,'(2(a,i0),3a,i0)')& 2089 'Interpolated WFK file contains only ', nbnds_kss_dense,' levels instead of ',Dtset%nband(1),' required;',ch10,& 2090 'The calculation will be done with nbands= ',nbnds_kss_dense 2091 ABI_WARNING(msg) 2092 ABI_ERROR("Not supported yet !") 2093 end if 2094 2095 ABI_MALLOC(nbands_temp,(Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol)) 2096 do isppol=1,Hdr_wfk_dense%nsppol 2097 do ik_ibz=1,Hdr_wfk_dense%nkpt 2098 nbands_temp(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) = Dtset%nband(1) 2099 end do 2100 end do 2101 2102 call Gsph_c%extend(Cryst, Dtset%ecutwfn, Gsph_x) 2103 call Gsph_x%print(unit=std_out,prtvol=Dtset%prtvol) 2104 2105 nqlwl=1 2106 ABI_MALLOC(qlwl,(3,nqlwl)) 2107 qlwl(:,nqlwl)= GW_Q0_DEFAULT 2108 2109 ! Compute Coulomb term on the largest G-sphere. 2110 if (Gsph_x%ng > Gsph_c%ng ) then 2111 call Vcp_dense%init(Gsph_x,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%gw_icutcoul,& 2112 Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm) 2113 else 2114 call Vcp_dense%init(Gsph_c,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%gw_icutcoul,& 2115 Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm) 2116 end if 2117 2118 ABI_FREE(qlwl) 2119 2120 bantot_dense=SUM(Hdr_wfk_dense%nband(1:Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol)) 2121 ABI_MALLOC(doccde,(bantot_dense)) 2122 ABI_MALLOC(eigen,(bantot_dense)) 2123 ABI_MALLOC(occfact,(bantot_dense)) 2124 doccde=zero; eigen=zero; occfact=zero 2125 2126 jj=0; ibtot=0 2127 do isppol=1,Hdr_wfk_dense%nsppol 2128 do ik_ibz=1,Hdr_wfk_dense%nkpt 2129 do ib=1,Hdr_wfk_dense%nband(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) 2130 ibtot=ibtot+1 2131 if (ib<=BSP%nbnds) then 2132 jj=jj+1 2133 occfact(jj)=Hdr_wfk_dense%occ(ibtot) 2134 eigen (jj)=energies_p_dense(ib,ik_ibz,isppol) 2135 end if 2136 end do 2137 end do 2138 end do 2139 2140 ABI_FREE(energies_p_dense) 2141 2142 ABI_MALLOC(npwarr,(kmesh_dense%nibz)) 2143 npwarr=BSP%npwwfn 2144 2145 ! CP modified 2146 !call ebands_init(bantot_dense,KS_BSt_dense,Dtset%nelect,doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,& 2147 !& Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,& 2148 !& Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,& 2149 !& hdr_wfk_dense%cellcharge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, & 2150 !& hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk) 2151 call ebands_init(bantot_dense,KS_BSt_dense,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& ! CP added 2152 & doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,& 2153 & Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,& 2154 & Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,& 2155 & hdr_wfk_dense%cellcharge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, & 2156 & hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk) 2157 ! End CP modified 2158 2159 ABI_FREE(doccde) 2160 ABI_FREE(eigen) 2161 ABI_FREE(npwarr) 2162 2163 ABI_FREE(nbands_temp) 2164 2165 ABI_FREE(occfact) 2166 2167 !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when 2168 ! the occupation scheme for semiconductors is used. 2169 call ebands_update_occ(KS_BSt_dense,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 2170 2171 call ebands_print(KS_BSt_dense,"Interpolated band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol) 2172 2173 call ebands_report_gap(KS_BSt_dense,header="Interpolated KS band structure",unit=std_out,mode_paral="COLL") 2174 2175 BSp%nkbz_interp = Kmesh_dense%nbz 2176 2177 call ebands_copy(KS_BSt_dense,QP_bst_dense) 2178 2179 SELECT CASE (Bsp%calc_type) 2180 CASE (BSE_HTYPE_RPA_KS) 2181 if (ABS(BSp%mbpt_sciss)>tol6) then 2182 write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies." 2183 call wrtout(std_out,msg) 2184 call ebands_apply_scissors(QP_BSt_dense,BSp%mbpt_sciss) 2185 else 2186 write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]." 2187 call wrtout(std_out,msg) 2188 end if 2189 ! 2190 CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw 2191 ABI_ERROR("Not yet implemented with interpolation !") 2192 CASE (BSE_HTYPE_RPA_QP) 2193 ABI_ERROR("Not implemented error!") 2194 CASE DEFAULT 2195 write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type 2196 ABI_ERROR(msg) 2197 END SELECT 2198 2199 call ebands_report_gap(QP_BSt_dense,header=" Interpolated QP band structure",unit=std_out,mode_paral="COLL") 2200 2201 ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index. 2202 ! FIXME: linewidths not coded. 2203 ABI_MALLOC(gw_energy, (BSp%nbnds,Kmesh_dense%nibz,Dtset%nsppol)) 2204 gw_energy = QP_BSt_dense%eig 2205 2206 ABI_MALLOC(Bsp%nreh_interp,(Hdr_wfk_dense%nsppol)) 2207 Bsp%nreh_interp=zero 2208 2209 call init_transitions(BSp%Trans_interp,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz_interp,Bsp%nbnds,& 2210 & Bsp%nkibz_interp,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,gw_energy,QP_BSt_dense%occ,Kmesh_dense%tab,minmax_tene,& 2211 & Bsp%nreh_interp) 2212 2213 ABI_FREE(gw_energy) 2214 2215 do spin=1,Dtset%nsppol 2216 write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh_interp(spin) 2217 call wrtout(std_out,msg) 2218 end do 2219 2220 if (ANY(Bsp%nreh_interp/=Bsp%nreh_interp(1))) then 2221 write(msg,'(a,(i0))')" BSE code does not support different number of transitions for the two spin channels",Bsp%nreh 2222 ABI_ERROR(msg) 2223 end if 2224 ! 2225 ! Create transition table vcks2t 2226 is1=BSp%lomo_min;is2=BSp%homo_max;is3=BSp%lumo_min;is4=BSp%humo_max 2227 ABI_MALLOC(Bsp%vcks2t_interp, (is1:is2,is3:is4,BSp%nkbz_interp,Dtset%nsppol)) 2228 Bsp%vcks2t_interp = 0 2229 2230 do spin=1,Dtset%nsppol 2231 do it=1,BSp%nreh_interp(spin) 2232 BSp%vcks2t_interp(BSp%Trans_interp(it,spin)%v,BSp%Trans_interp(it,spin)%c, BSp%Trans_interp(it,spin)%k,spin) = it 2233 end do 2234 end do 2235 2236 hexc_size = SUM(Bsp%nreh_interp); if (Bsp%use_coupling>0) hexc_size=2*hexc_size 2237 if (Bsp%nstates_interp<=0) then 2238 Bsp%nstates_interp=hexc_size 2239 else 2240 if (Bsp%nstates_interp>hexc_size) then 2241 Bsp%nstates_interp=hexc_size 2242 write(msg,'(2(a,i0),2a)')& 2243 "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates_interp,ch10,& 2244 "the number of excitonic states nstates has been modified" 2245 ABI_WARNING(msg) 2246 end if 2247 end if 2248 2249 DBG_EXIT("COLL") 2250 2251 end subroutine setup_bse_interp