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-2024 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
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_bethe_salpeter 25 26 use defs_basis 27 use defs_wvltypes 28 use m_bs_defs 29 use m_abicore 30 use m_xmpi 31 use m_errors 32 use m_screen 33 use m_nctk 34 use m_distribfft 35 use netcdf 36 use m_hdr 37 use m_dtset 38 use m_dtfil 39 use m_crystal 40 41 use defs_datatypes, only : pseudopotential_type, ebands_t 42 use defs_abitypes, only : MPI_type 43 use m_gwdefs, only : GW_Q0_DEFAULT 44 use m_time, only : timab 45 use m_fstrings, only : strcat, sjoin, endswith, itoa 46 use m_io_tools, only : file_exists, iomode_from_fname 47 use m_geometry, only : mkrdim, metric, normv 48 use m_hide_lapack, only : matrginv 49 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 50 use m_fftcore, only : print_ngfft 51 use m_fft_mesh, only : rotate_FFT_mesh, get_gftt, setmesh 52 use m_fft, only : fourdp 53 use m_bz_mesh, only : kmesh_t, get_ng0sh, find_qmesh, make_mesh 54 use m_double_grid, only : double_grid_t, double_grid_init, double_grid_free 55 use m_ebands, only : ebands_init, ebands_print, ebands_copy, ebands_free, & 56 ebands_update_occ, ebands_get_valence_idx, ebands_apply_scissors, ebands_report_gap 57 use m_kg, only : getph 58 use m_gsphere, only : gsphere_t 59 use m_vcoul, only : vcoul_t 60 use m_qparticles, only : rdqps, rdgw !, show_QP , rdgw 61 use m_wfd, only : wfd_init, wfdgw_t, test_charge 62 use m_wfk, only : wfk_read_eigenvalues 63 use m_energies, only : energies_type, energies_init 64 use m_io_screening, only : hscr_t, hscr_io 65 use m_haydock, only : exc_haydock_driver 66 use m_exc_diago, only : exc_diago_driver 67 use m_exc_analyze, only : exc_den 68 use m_eprenorms, only : eprenorms_t, eprenorms_free, eprenorms_from_epnc, eprenorms_bcast 69 use m_pawang, only : pawang_type 70 use m_pawrad, only : pawrad_type 71 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 72 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 73 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 74 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init 75 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free,& 76 & pawrhoij_inquire_dim, pawrhoij_symrhoij 77 use m_pawdij, only : pawdij, symdij 78 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 79 use m_paw_hr, only : pawhur_t, pawhur_free, pawhur_init 80 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 81 use m_paw_sphharm, only : setsym_ylm 82 use m_paw_denpot, only : pawdenpot 83 use m_paw_init, only : pawinit,paw_gencond 84 use m_paw_onsite, only : pawnabla_init 85 use m_paw_dmft, only : paw_dmft_type 86 use m_paw_mkrho, only : denfgr 87 use m_paw_nhat, only : nhatgrid,pawmknhat 88 use m_paw_tools, only : chkpawovlp, pawprt 89 use m_paw_correlations,only : pawpuxinit 90 use m_exc_build, only : exc_build_ham 91 use m_setvtr, only : setvtr 92 use m_mkrho, only : prtrhomxmn 93 use m_pspini, only : pspini 94 use m_drivexc, only : mkdenpos 95 96 implicit none 97 98 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
165 subroutine bethe_salpeter(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim,xred) 166 167 !Arguments ------------------------------------ 168 !scalars 169 character(len=8),intent(in) :: codvsn 170 type(datafiles_type),intent(inout) :: Dtfil 171 type(dataset_type),intent(inout) :: Dtset 172 type(pawang_type),intent(inout) :: Pawang 173 type(pseudopotential_type),intent(inout) :: Psps 174 !arrays 175 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,Dtset%natom) 176 type(pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw) 177 type(pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw) 178 179 !Local variables ------------------------------ 180 !scalars 181 integer,parameter :: tim_fourdp0=0,level=40,ipert0=0,idir0=0,cplex1=1,master=0,option1=1 182 integer :: band,cplex_rhoij,spin,ik_ibz,mqmem,iwarn 183 integer :: has_dijU,has_dijso,gnt_option 184 integer :: ik_bz,mband 185 integer :: choice 186 integer :: ider 187 integer :: usexcnhat,nfft_osc,mgfft_osc 188 integer :: isym,izero 189 integer :: optcut,optgr0,optgr1,optgr2,option,optrad,optrhoij,psp_gencond 190 integer :: ngrvdw,nhatgrdim,nkxc1,nprocs,nspden_rhoij,nzlmopt,ifft 191 integer :: my_rank,rhoxsp_method,comm 192 integer :: mgfftf,spin_opt,which_fixed 193 integer :: nscf,nbsc,nkxc,n3xccc 194 integer :: nfftf,nfftf_tot,nfftot_osc,my_minb,my_maxb 195 integer :: optene,moved_atm_inside,moved_rhor,initialized,istep,ierr 196 real(dp) :: ucvol,drude_plsmf,ecore,ecut_eff,ecutdg_eff,norm 197 real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp 198 real(dp) :: compch_fft,compch_sph,gsq_osc 199 real(dp) :: vxcavg 200 logical :: iscompatibleFFT,is_dfpt=.false.,paw_add_onsite,call_pawinit 201 character(len=500) :: msg 202 character(len=fnlen) :: wfk_fname,w_fname 203 type(Pawfgr_type) :: Pawfgr 204 type(excfiles) :: BS_files 205 type(excparam) :: BSp 206 type(paw_dmft_type) :: Paw_dmft 207 type(MPI_type) :: MPI_enreg_seq 208 type(crystal_t) :: Cryst 209 type(kmesh_t) :: Kmesh,Qmesh 210 type(gsphere_t) :: Gsph_x,Gsph_c,Gsph_x_dense,Gsph_c_dense 211 type(Hdr_type) :: Hdr_wfk,Hdr_bse 212 type(ebands_t) :: ks_ebands,qp_ebands,ks_ebands_dense,qp_ebands_dense 213 type(Energies_type) :: KS_energies 214 type(vcoul_t) :: Vcp, Vcp_dense 215 type(wfdgw_t) :: Wfd, Wfd_dense 216 type(screen_t) :: W 217 type(screen_info_t) :: W_info 218 type(wvl_data) :: wvl 219 type(kmesh_t) :: Kmesh_dense,Qmesh_dense 220 type(Hdr_type) :: Hdr_wfk_dense 221 type(double_grid_t) :: grid 222 type(eprenorms_t) :: Epren 223 !arrays 224 integer :: ngfft_osc(18),ngfftc(18),ngfftf(18),nrcell(3) 225 integer,allocatable :: ktabr(:,:),l_size_atm(:) 226 integer,allocatable :: nband(:,:),nq_spl(:),irottb(:,:) 227 integer,allocatable :: qp_vbik(:,:) 228 integer,allocatable :: gfft_osc(:,:) 229 real(dp),parameter :: k0(3)=zero 230 real(dp) :: tsec(2),gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3),rprimd(3,3),eh_rcoord(3),strsxc(6) 231 real(dp),allocatable :: ph1df(:,:),prev_rhor(:,:),ph1d(:,:) 232 real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:),ks_rhor(:,:),qp_aerhor(:,:) 233 real(dp),allocatable :: qp_rhor(:,:),qp_rhog(:,:) !,qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:) 234 real(dp),allocatable :: qp_rhor_paw(:,:),qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:),qp_nhat(:,:) 235 real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:) 236 real(dp),allocatable :: vpsp(:),xccc3d(:) 237 real(dp),allocatable :: ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:) 238 real(dp),allocatable :: kxc(:,:) !,qp_kxc(:,:) 239 complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:) 240 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:) 241 type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:) 242 type(Pawrhoij_type),allocatable :: prev_Pawrhoij(:) !QP_pawrhoij(:), 243 type(pawpwff_t),allocatable :: Paw_pwff(:) 244 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 245 type(pawhur_t),allocatable :: Hur(:) 246 type(Paw_ij_type),allocatable :: KS_paw_ij(:) 247 type(Paw_an_type),allocatable :: KS_paw_an(:) 248 249 !************************************************************************ 250 251 DBG_ENTER('COLL') 252 253 call timab(650,1,tsec) ! bse(Total) 254 call timab(651,1,tsec) ! bse(Init1) 255 256 comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 257 258 wfk_fname = dtfil%fnamewffk 259 260 if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then 261 ABI_ERROR(msg) 262 end if 263 call xmpi_bcast(wfk_fname, master, comm, ierr) 264 265 write(msg,'(8a)')& 266 ' Exciton: Calculation of dielectric properties by solving the Bethe-Salpeter equation ',ch10,& 267 ' in frequency domain and reciprocal space on a transitions basis set. ',ch10,& 268 ' Based on a program developed by L. Reining, V. Olevano, F. Sottile, ',ch10,& 269 ' S. Albrecht, and G. Onida. Incorporated in ABINIT by M. Giantomassi. ',ch10 270 call wrtout([std_out, ab_out], msg) 271 272 #ifdef HAVE_GW_DPC 273 if (gwpc/=8) then 274 write(msg,'(6a)')ch10,& 275 ' Number of bytes for double precision complex /=8 ',ch10,& 276 ' Cannot continue due to kind mismatch in BLAS library ',ch10,& 277 ' Some BLAS interfaces are not generated by abilint ' 278 ABI_ERROR(msg) 279 end if 280 write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10 281 #else 282 write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10 283 #endif 284 call wrtout([std_out, ab_out], msg) 285 286 !=== Some variables need to be initialized/nullify at start === 287 call energies_init(KS_energies) 288 usexcnhat=0 289 call mkrdim(acell,rprim,rprimd) 290 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 291 ! 292 !=== Define FFT grid(s) sizes === 293 !* Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the 294 !(usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 295 !See also NOTES in the comments at the beginning of this file. 296 !NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 297 298 !TODO Recheck getng, should use same trick as that used in screening and sigma. 299 call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 300 gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0) 301 302 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 303 nfftf_tot=PRODUCT(ngfftf(1:3)) 304 305 ! Fake MPI_type for the sequential part. 306 call initmpi_seq(MPI_enreg_seq) 307 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 308 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 309 310 ! =========================================== 311 ! === Open and read pseudopotential files === 312 ! =========================================== 313 call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm) 314 315 ! === Initialization of basic objects including the BSp structure that defines the parameters of the run === 316 call setup_bse(codvsn,acell,rprim,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,& 317 Cryst,Kmesh,Qmesh,ks_ebands,qp_ebands,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,wvl%descr) 318 319 if (BSp%use_interp) then 320 call setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh,Kmesh_dense,& 321 Qmesh_dense,ks_ebands_dense,qp_ebands_dense,Gsph_x_dense,Gsph_c_dense,& 322 Vcp_dense,Hdr_wfk_dense,grid,comm) 323 end if 324 325 !call timab(652,2,tsec) ! setup_bse 326 327 nfftot_osc=PRODUCT(ngfft_osc(1:3)) 328 nfft_osc =nfftot_osc !no FFT // 329 mgfft_osc =MAXVAL(ngfft_osc(1:3)) 330 331 call print_ngfft(ngfft_osc,header='FFT mesh used for oscillator strengths') 332 333 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 334 KS_energies%e_corepsp=ecore/Cryst%ucvol 335 336 ! 337 !============================ 338 !==== PAW initialization ==== 339 !============================ 340 if (Dtset%usepaw==1) then 341 call chkpawovlp(Cryst%natom,Cryst%ntypat,Dtset%pawovlp,Pawtab,Cryst%rmet,Cryst%typat,xred) 342 343 ABI_MALLOC(KS_Pawrhoij,(Cryst%natom)) 344 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 345 nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc) 346 call pawrhoij_alloc(KS_Pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 347 348 ! Initialize values for several basic arrays === 349 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 350 351 ! Test if we have to call pawinit 352 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 353 354 if (psp_gencond==1.or.call_pawinit) then 355 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 356 call pawinit(Dtset%effmass_free,gnt_option,gsqcut_shp,zero,Dtset%pawlcutd,Dtset%pawlmix,& 357 Psps%mpsang,Dtset%pawnphi,Cryst%nsym,Dtset%pawntheta,Pawang,Pawrad,& 358 Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,Dtset%ixc,Dtset%usepotzero) 359 360 ! Update internal values 361 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 362 else 363 if (Pawtab(1)%has_kij ==1) Pawtab(1:Cryst%ntypat)%has_kij =2 364 if (Pawtab(1)%has_nabla==1) Pawtab(1:Cryst%ntypat)%has_nabla=2 365 end if 366 Psps%n1xccc=MAXVAL(Pawtab(1:Cryst%ntypat)%usetcore) 367 368 ! Initialize optional flags in Pawtab to zero 369 ! (Cannot be done in Pawinit since the routine is called only if some parameters are changed) 370 Pawtab(:)%has_nabla = 0 371 Pawtab(:)%usepawu = 0 372 Pawtab(:)%useexexch = 0 373 Pawtab(:)%exchmix =zero 374 Pawtab(:)%lamb_shielding = zero 375 376 ! * Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit. 377 ! TODO solve problem with memory leak and clean this part as well as the associated flag 378 call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab) 379 380 call setsym_ylm(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,Cryst%rprimd,Cryst%symrec,Pawang%zarot) 381 382 ! Initialize and compute data for DFT+U 383 Paw_dmft%use_dmft=Dtset%usedmft 384 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 385 is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Dtset%nspinor,Cryst%ntypat,Dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,& 386 Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu) 387 if (Dtset%usepawu>0.or.Dtset%useexexch>0) then 388 ABI_ERROR('BS equation with DFT+U not completely coded!') 389 end if 390 if (my_rank == master) call pawtab_print(Pawtab) 391 392 ! Get Pawrhoij from the header of the WFK file. 393 call pawrhoij_copy(Hdr_wfk%pawrhoij,KS_Pawrhoij) 394 395 ! Re-symmetrize rhoij === 396 ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly 397 choice=1; optrhoij=1 398 ! call pawrhoij_symrhoij(KS_Pawrhoij,KS_Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,& 399 ! & Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,& 400 ! & Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 401 402 ! Evaluate form factor of radial part of phi.phj-tphi.tphj === 403 rhoxsp_method=1 ! Arnaud-Alouani 404 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 405 406 ABI_MALLOC(gfft_osc,(3,nfftot_osc)) 407 call get_gftt(ngfft_osc,k0,gmet,gsq_osc,gfft_osc) 408 ABI_FREE(gfft_osc) 409 410 ! Set up q grids, make qmax 20% larger than largest expected: 411 ABI_MALLOC(nq_spl,(Psps%ntypat)) 412 ABI_MALLOC(qmax,(Psps%ntypat)) 413 nq_spl = Psps%mqgrid_ff 414 qmax = SQRT(gsq_osc)*1.2d0 ! qmax=Psps%qgrid_ff(Psps%mqgrid_ff) 415 ABI_MALLOC(Paw_pwff,(Psps%ntypat)) 416 417 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 418 419 ABI_FREE(nq_spl) 420 ABI_FREE(qmax) 421 ! 422 ! Variables/arrays related to the fine FFT grid === 423 ABI_MALLOC(ks_nhat,(nfftf,Dtset%nspden)) 424 ks_nhat=zero 425 ABI_MALLOC(Pawfgrtab,(Cryst%natom)) 426 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 427 call pawfgrtab_init(Pawfgrtab,cplex1,l_size_atm,Dtset%nspden,Dtset%typat) 428 ABI_FREE(l_size_atm) 429 compch_fft=greatest_real 430 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 431 ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation) 432 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 433 write(msg,'(a,i0)')' bethe_salpeter : using usexcnhat = ',usexcnhat 434 call wrtout(std_out,msg) 435 ! 436 ! Identify parts of the rectangular grid where the density has to be calculated === 437 optcut=0; optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm 438 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 439 440 call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 441 optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 442 else 443 ABI_MALLOC(Paw_pwff,(0)) 444 end if !End of PAW Initialization 445 446 ! Consistency check and additional stuff done only for GW with PAW. 447 if (Dtset%usepaw==1) then 448 if (Dtset%ecutwfn < Dtset%ecut) then 449 write(msg,"(5a)")& 450 "WARNING - ",ch10,& 451 " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 452 " an excessive truncation of the planewave basis set can lead to unphysical results." 453 call wrtout(ab_out,msg,'COLL') 454 end if 455 456 ABI_CHECK(Dtset%usedmft==0,"DMFT + BSE not allowed") 457 ABI_CHECK(Dtset%useexexch==0,"LEXX + BSE not allowed") 458 end if 459 460 ! Allocate these arrays anyway, since they are passed to subroutines. 461 if (.not.allocated(ks_nhat)) then 462 ABI_MALLOC(ks_nhat,(nfftf,0)) 463 end if 464 465 !================================================== 466 !==== Read KS band structure from the KSS file ==== 467 !================================================== 468 469 ! Initialize wave function handler, allocate wavefunctions. 470 my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds 471 ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol)) 472 nband=mband 473 474 !At present, no memory distribution, each node has the full set of states. 475 ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol)) 476 bks_mask=.TRUE. 477 478 ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol)) 479 keep_ur=.FALSE.; if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE. 480 481 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Dtset%nsppol,bks_mask,& 482 Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_osc,& 483 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 484 485 ABI_FREE(bks_mask) 486 ABI_FREE(nband) 487 ABI_FREE(keep_ur) 488 489 call wfd%print(header="Wavefunctions used to construct the e-h basis set",mode_paral='PERS') 490 491 call timab(651,2,tsec) ! bse(Init1) 492 call timab(653,1,tsec) ! bse(rdkss) 493 494 call wfd%read_wfk(wfk_fname,iomode_from_fname(wfk_fname)) 495 496 ! This test has been disabled (too expensive!) 497 if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 498 499 call timab(653,2,tsec) ! bse(rdkss) 500 call timab(655,1,tsec) ! bse(mkrho) 501 502 !TODO : check the consistency of Wfd with Wfd_dense !!! 503 if (BSp%use_interp) then 504 ! Initialize wave function handler, allocate wavefunctions. 505 my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds 506 ABI_MALLOC(nband,(Kmesh_dense%nibz,Dtset%nsppol)) 507 nband=mband 508 509 ! At present, no memory distribution, each node has the full set of states. 510 ! albeit we allocate only the states that are used. 511 ABI_MALLOC(bks_mask,(mband,Kmesh_dense%nibz,Dtset%nsppol)) 512 bks_mask=.False. 513 do spin=1,Bsp%nsppol 514 bks_mask(Bsp%lomo_spin(spin):Bsp%humo_spin(spin),:,spin) = .True. 515 end do 516 !bks_mask=.TRUE. 517 518 ABI_MALLOC(keep_ur,(mband,Kmesh_dense%nibz,Dtset%nsppol)) 519 keep_ur=.FALSE.; if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE. 520 521 call wfd_init(Wfd_dense,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh_dense%nibz,Dtset%nsppol,& 522 bks_mask,Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,ngfft_osc,& 523 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 524 525 ABI_FREE(bks_mask) 526 ABI_FREE(nband) 527 ABI_FREE(keep_ur) 528 529 call wfd_dense%print(header="Wavefunctions on the dense K-mesh used for interpolation",mode_paral='PERS') 530 531 call wfd_dense%read_wfk(Dtfil%fnameabi_wfkfine, iomode_from_fname(dtfil%fnameabi_wfkfine)) 532 !call wfd_dense%update_bkstab() 533 534 ! This test has been disabled (too expensive!) 535 if (.False.) call wfd_dense%test_ortho(Cryst,Pawtab,unit=std_out,mode_paral="COLL") 536 end if 537 538 !=== Calculate the FFT index of $(R^{-1}(r-\tau))$ === 539 !* S=\transpose R^{-1} and k_BZ = S k_IBZ 540 !* irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 541 ABI_MALLOC(irottb,(nfftot_osc,Cryst%nsym)) 542 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfft_osc,irottb,iscompatibleFFT) 543 544 ABI_MALLOC(ktabr,(nfftot_osc,Kmesh%nbz)) 545 do ik_bz=1,Kmesh%nbz 546 isym=Kmesh%tabo(ik_bz) 547 do ifft=1,nfftot_osc 548 ktabr(ifft,ik_bz)=irottb(ifft,isym) 549 end do 550 end do 551 ABI_FREE(irottb) 552 ! 553 !=========================== 554 !=== COMPUTE THE DENSITY === 555 !=========================== 556 !* Evaluate Planewave part (complete charge in case of NC pseudos). 557 ! 558 ABI_MALLOC(ks_rhor, (nfftf, Wfd%nspden)) 559 call wfd%mkrho(cryst, psps, ks_ebands, ngfftf, nfftf, ks_rhor) 560 ! 561 !=== Additional computation for PAW === 562 nhatgrdim=0 563 if (Dtset%usepaw==1) then 564 ! 565 ! Calculate the compensation charge nhat. 566 if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc 567 ider=2*nhatgrdim; izero=0; qphon(:)=zero 568 if (nhatgrdim>0) then 569 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3)) 570 end if 571 572 call pawmknhat(compch_fft,cplex1,ider,idir0,ipert0,izero,Cryst%gprimd,& 573 Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 574 Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,qphon,Cryst%rprimd,& 575 Cryst%ucvol,Dtset%usewvl,Cryst%xred) 576 577 ! Evaluate onsite energies, potentials, densities === 578 ! * Initialize variables/arrays related to the PAW spheres. 579 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 580 ABI_MALLOC(KS_paw_ij,(Cryst%natom)) 581 call paw_ij_nullify(KS_paw_ij) 582 583 has_dijso=Dtset%pawspnorb 584 has_dijU=merge(0,1,Dtset%usepawu==0) 585 586 call paw_ij_init(KS_paw_ij,cplex1,Dtset%nspinor,Dtset%nsppol,& 587 Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 588 has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=0,has_dijxc_hat=0,has_dijxc_val=0,& 589 has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1) 590 591 ABI_MALLOC(KS_paw_an,(Cryst%natom)) 592 call paw_an_nullify(KS_paw_an) 593 594 nkxc1=0 595 call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,& 596 cplex1,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0) 597 598 ! Calculate onsite vxc with and without core charge === 599 nzlmopt=-1; option=0; compch_sph=greatest_real 600 601 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,& 602 Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 603 Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,& 604 Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,& 605 Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,& 606 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_ebands%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,Pawang,Pawrad,KS_Pawrhoij,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,& 681 & Cryst%rprimd,strsxc,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_ebands stores energies and occ. used for the calculation === 717 !* Initialize qp_ebands with KS values. 718 !* In case of SC update qp_ebands 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_ebands,qp_ebands) 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_ebands. 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_ebands,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_ebands,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_ebands,Dtset%spinmagntarget,prtvol=0) 811 ABI_MALLOC(qp_vbik,(qp_ebands%nkpt,qp_ebands%nsppol)) 812 qp_vbik(:,:) = ebands_get_valence_idx(qp_ebands) 813 ABI_FREE(qp_vbik) 814 815 call wfd%mkrho(cryst, psps, qp_ebands, 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_ebands,qp_ebands,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_ebands,qp_ebands,Wfd,Psps,Pawtab,Hur,Epren, & 935 kmesh_dense=Kmesh_dense, ks_bst_dense=ks_ebands_dense, qp_bst_dense=qp_ebands_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_ebands,qp_ebands,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_ebands) 975 call ebands_free(qp_ebands) 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_ebands_dense) 993 call ebands_free(qp_ebands_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_ebands<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_ebands,qp_ebands,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_ebands,qp_ebands 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_indices(:,:),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_ebands,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_ebands,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 1542 1543 call ebands_print(ks_ebands,"Band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol) 1544 1545 call ebands_report_gap(ks_ebands,header=" KS band structure",unit=std_out,mode_paral="COLL") 1546 1547 ABI_MALLOC(val_indices,(ks_ebands%nkpt,ks_ebands%nsppol)) 1548 val_indices = ebands_get_valence_idx(ks_ebands) 1549 1550 do spin=1,ks_ebands%nsppol 1551 val_idx(spin) = val_indices(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_indices(1,spin) /= val_indices(:,spin)) ) then 1555 ABI_ERROR("BSE code does not support metals") 1556 end if 1557 end do 1558 1559 ABI_FREE(val_indices) 1560 ! 1561 ! === Create the BSE header === 1562 call hdr_init(ks_ebands,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 call hdr_bse%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 1572 1573 ABI_FREE(occfact) 1574 1575 if (Dtset%usepaw==1) call pawrhoij_free(Pawrhoij) 1576 ABI_FREE(Pawrhoij) 1577 1578 ! Find optimal value for G-sphere enlargment due to oscillator matrix elements 1579 ! We will split k-points over processors 1580 call xmpi_split_work(Kmesh%nbz, comm, my_k1, my_k2) 1581 1582 ! If there is no work to do, just skip the computation 1583 if (my_k2-my_k1+1 <= 0) then 1584 ng0sh_opt(:)=(/zero,zero,zero/) 1585 else 1586 ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true since bz is buggy 1587 ! * -one is used because we loop over all the possibile differences, unlike screening 1588 call get_ng0sh(my_k2-my_k1+1,Kmesh%bz(:,my_k1:my_k2),Kmesh%nbz,Kmesh%bz,& 1589 & Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt) 1590 end if 1591 1592 call xmpi_max(ng0sh_opt,BSp%mg0,comm,ierr) 1593 1594 write(msg,'(a,3(i0,1x))') ' optimal value for ng0sh = ',BSp%mg0 1595 call wrtout(std_out,msg) 1596 1597 ! === Setup of the FFT mesh for the oscilator strengths === 1598 ! * ngfft_osc(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening. 1599 ! * Here we redefine ngfft_osc(1:6) according to the following options : 1600 ! 1601 ! method==0 --> FFT grid read from fft.in (debugging purpose) 1602 ! method==1 --> Normal FFT mesh 1603 ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 1604 ! method==3 --> Doubled FFT grid, same as the the FFT for the density, 1605 ! 1606 ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library 1607 ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries 1608 ! 1609 ngfft_osc(1:18)=Dtset%ngfft(1:18); method=2 1610 if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0 1611 if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1 1612 if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2 1613 if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3 1614 enforce_sym=MOD(Dtset%fftgw,10) 1615 1616 call setmesh(gmet,Gsph_x%gvec,ngfft_osc,BSp%npwvec,BSp%npweps,BSp%npwwfn,nfftot_osc,method,BSp%mg0,Cryst,enforce_sym) 1617 nfftot_osc=PRODUCT(ngfft_osc(1:3)) 1618 1619 call print_ngfft(ngfft_osc,"FFT mesh for oscillator matrix elements",std_out,"COLL",prtvol=Dtset%prtvol) 1620 ! 1621 ! BSp%homo gives the 1622 !BSp%homo = val_idx(1) 1623 ! highest occupied band for each spin 1624 BSp%homo_spin = val_idx 1625 1626 ! TODO generalize the code to account for this unlikely case. 1627 !if (Dtset%nsppol==2) then 1628 ! ABI_CHECK(BSp%homo == val_idx(2),"Different valence indices for spin up and down") 1629 !end if 1630 1631 !BSp%lumo = BSp%homo + 1 1632 !BSp%humo = BSp%nbnds 1633 !BSp%nbndv = BSp%homo - BSp%lomo + 1 1634 !BSp%nbndc = BSp%nbnds - BSp%homo 1635 1636 BSp%lumo_spin = BSp%homo_spin + 1 1637 BSp%humo_spin = BSp%nbnds 1638 BSp%nbndv_spin = BSp%homo_spin - BSp%lomo_spin + 1 1639 BSp%nbndc_spin = BSp%nbnds - BSp%homo_spin 1640 BSp%maxnbndv = MAXVAL(BSp%nbndv_spin(:)) 1641 BSp%maxnbndc = MAXVAL(BSp%nbndc_spin(:)) 1642 1643 BSp%nkbz = Kmesh%nbz 1644 1645 call ebands_copy(ks_ebands,qp_ebands) 1646 ABI_MALLOC(igwene,(qp_ebands%mband,qp_ebands%nkpt,qp_ebands%nsppol)) 1647 igwene=zero 1648 1649 call bsp_calctype2str(Bsp,msg) 1650 call wrtout(std_out,"Calculation type: "//TRIM(msg)) 1651 1652 SELECT CASE (Bsp%calc_type) 1653 CASE (BSE_HTYPE_RPA_KS) 1654 if (ABS(BSp%mbpt_sciss)>tol6) then 1655 write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies." 1656 call wrtout(std_out,msg) 1657 call ebands_apply_scissors(qp_ebands,BSp%mbpt_sciss) 1658 else 1659 write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]." 1660 call wrtout(std_out,msg) 1661 end if 1662 1663 CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw 1664 gw_fname=TRIM(Dtfil%filnam_ds(4))//'_GW' 1665 gw_fname="__in.gw__" 1666 if (.not.file_exists(gw_fname)) then 1667 msg = " File "//TRIM(gw_fname)//" not found. Aborting now" 1668 ABI_ERROR(msg) 1669 end if 1670 1671 call rdgw(qp_ebands,gw_fname,igwene,extrapolate=.FALSE.) ! here gwenergy is real 1672 1673 do isppol=1,Dtset%nsppol 1674 write(std_out,*) ' k GW energies [eV]' 1675 do ik_ibz=1,Kmesh%nibz 1676 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(qp_ebands%eig(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds) 1677 end do 1678 write(std_out,*) ' k Im GW energies [eV]' 1679 do ik_ibz=1,Kmesh%nibz 1680 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(igwene(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds) 1681 end do 1682 end do 1683 ! 1684 ! If required apply the scissors operator on top of the QP bands structure (!) 1685 if (ABS(BSp%mbpt_sciss)>tol6) then 1686 write(msg,'(a,f8.2,a)')' Applying a scissors operator ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the QP energies!" 1687 ABI_COMMENT(msg) 1688 call ebands_apply_scissors(qp_ebands,BSp%mbpt_sciss) 1689 end if 1690 1691 CASE (BSE_HTYPE_RPA_QP) 1692 ABI_ERROR("Not implemented error!") 1693 1694 CASE DEFAULT 1695 write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type 1696 ABI_ERROR(msg) 1697 END SELECT 1698 1699 call ebands_report_gap(qp_ebands,header=" QP band structure",unit=std_out,mode_paral="COLL") 1700 1701 ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index. 1702 ! FIXME: linewidths not coded. 1703 ABI_MALLOC(gw_energy,(BSp%nbnds,Kmesh%nibz,Dtset%nsppol)) 1704 gw_energy = qp_ebands%eig 1705 1706 BSp%have_complex_ene = ANY(igwene > tol16) 1707 1708 ! Compute the number of resonant transitions, nreh, for the two spin channels and initialize BSp%Trans. 1709 ABI_MALLOC(Bsp%nreh,(Bsp%nsppol)) 1710 1711 ! Possible cutoff on the transitions. 1712 BSp%ircut = Dtset%bs_eh_cutoff(1) 1713 BSp%uvcut = Dtset%bs_eh_cutoff(2) 1714 1715 call init_transitions(BSp%Trans,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz,Bsp%nbnds,Bsp%nkibz,& 1716 & BSp%nsppol,Dtset%nspinor,gw_energy,qp_ebands%occ,Kmesh%tab,minmax_tene,Bsp%nreh) 1717 1718 ! Setup of the frequency mesh for the absorption spectrum. 1719 ! If not specified, use the min-max resonant transition energy and make it 10% smaller|larger. 1720 1721 !if (ABS(Dtset%bs_freq_mesh(1)) < tol6) then 1722 ! Dtset%bs_freq_mesh(1) = MAX(minmax_tene(1) - minmax_tene(1) * 0.1, zero) 1723 !end if 1724 1725 if (ABS(Dtset%bs_freq_mesh(2)) < tol6) then 1726 Dtset%bs_freq_mesh(2) = minmax_tene(2) + minmax_tene(2) * 0.1 1727 end if 1728 1729 Bsp%omegai = Dtset%bs_freq_mesh(1) 1730 Bsp%omegae = Dtset%bs_freq_mesh(2) 1731 Bsp%domega = Dtset%bs_freq_mesh(3) 1732 BSp%broad = Dtset%zcut 1733 1734 ! The frequency mesh (including the complex imaginary shift) 1735 BSp%nomega = (BSp%omegae - BSp%omegai)/BSp%domega + 1 1736 ABI_MALLOC(BSp%omega,(BSp%nomega)) 1737 do io=1,BSp%nomega 1738 BSp%omega(io) = (BSp%omegai + (io-1)*BSp%domega) + j_dpc*BSp%broad 1739 end do 1740 1741 ABI_FREE(gw_energy) 1742 ABI_FREE(igwene) 1743 1744 do spin=1,Bsp%nsppol 1745 write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh(spin) 1746 call wrtout(std_out,msg) 1747 end do 1748 1749 if (ANY(Bsp%nreh/=Bsp%nreh(1))) then 1750 write(msg,'(a,2(i0,1x))')" BSE code with different number of transitions for the two spin channels: ",Bsp%nreh 1751 ABI_WARNING(msg) 1752 end if 1753 ! 1754 ! Create transition table vcks2t 1755 Bsp%lomo_min = MINVAL(BSp%lomo_spin) 1756 Bsp%homo_max = MAXVAL(BSp%homo_spin) 1757 Bsp%lumo_min = MINVAL(BSp%lumo_spin) 1758 Bsp%humo_max = MAXVAL(BSp%humo_spin) 1759 1760 ABI_MALLOC(Bsp%vcks2t,(BSp%lomo_min:BSp%homo_max,BSp%lumo_min:BSp%humo_max,BSp%nkbz,Dtset%nsppol)) 1761 Bsp%vcks2t = 0 1762 1763 do spin=1,BSp%nsppol 1764 do it=1,BSp%nreh(spin) 1765 BSp%vcks2t(BSp%Trans(it,spin)%v,BSp%Trans(it,spin)%c,BSp%Trans(it,spin)%k,spin) = it 1766 end do 1767 end do 1768 1769 hexc_size = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hexc_size=2*hexc_size 1770 if (Bsp%nstates<=0) then 1771 Bsp%nstates=hexc_size 1772 else 1773 if (Bsp%nstates>hexc_size) then 1774 Bsp%nstates=hexc_size 1775 write(msg,'(2(a,i0),2a)')& 1776 "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates,ch10,& 1777 "the number of excitonic states nstates has been modified" 1778 ABI_WARNING(msg) 1779 end if 1780 end if 1781 1782 msg=' Fundamental parameters for the solution of the Bethe-Salpeter equation:' 1783 call print_bs_parameters(BSp,unit=std_out,header=msg,mode_paral="COLL",prtvol=Dtset%prtvol) 1784 call print_bs_parameters(BSp,unit=ab_out, header=msg,mode_paral="COLL") 1785 1786 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 1787 write(msg,'(3a,9i2,2a,3f6.3,2a)')& 1788 "The first symmetry operation should be the Identity with zero tnons while ",ch10,& 1789 "symrec(:,:,1) = ",Cryst%symrec(:,:,1),ch10,& 1790 "tnons(:,1) = ",Cryst%tnons(:,1),ch10,& 1791 "This is not allowed, sym_rhotwgq0 should be changed." 1792 ABI_ERROR(msg) 1793 end if 1794 ! 1795 ! Prefix for generic output files. 1796 BS_files%out_basename = TRIM(Dtfil%filnam_ds(4)) 1797 ! 1798 ! Search for files to restart from. 1799 if (Dtset%gethaydock/=0 .or. Dtset%irdhaydock/=0) then 1800 BS_files%in_haydock_basename = TRIM(Dtfil%fnameabi_haydock) 1801 end if 1802 1803 test_file = Dtfil%fnameabi_bsham_reso 1804 if (file_exists(test_file)) then 1805 BS_files%in_hreso = test_file 1806 else 1807 BS_files%out_hreso = TRIM(Dtfil%filnam_ds(4))//'_BSR' 1808 end if 1809 1810 test_file = Dtfil%fnameabi_bsham_coup 1811 if (file_exists(test_file) ) then 1812 BS_files%in_hcoup = test_file 1813 else 1814 BS_files%out_hcoup = TRIM(Dtfil%filnam_ds(4))//'_BSC' 1815 end if 1816 ! 1817 ! in_eig is the name of the input file with eigenvalues and eigenvectors 1818 ! constructed from getbseig or irdbseig. out_eig is the name of the output file 1819 ! produced by this dataset. in_eig_exists checks for the presence of the input file. 1820 ! 1821 if (file_exists(Dtfil%fnameabi_bseig)) then 1822 BS_files%in_eig = Dtfil%fnameabi_bseig 1823 else 1824 BS_files%out_eig = TRIM(BS_files%out_basename)//"_BSEIG" 1825 end if 1826 1827 call print_bs_files(BS_files,unit=std_out) 1828 ! 1829 ! ========================================================== 1830 ! ==== Temperature dependence of the spectrum ============== 1831 ! ========================================================== 1832 BSp%do_ep_renorm = .FALSE. 1833 BSp%do_lifetime = .FALSE. ! Not yet implemented 1834 1835 ep_nc_fname = 'test_EP.nc' 1836 if(file_exists(ep_nc_fname)) then 1837 BSp%do_ep_renorm = .TRUE. 1838 1839 if(my_rank == master) then 1840 call eprenorms_from_epnc(Epren,ep_nc_fname) 1841 end if 1842 call eprenorms_bcast(Epren,master,comm) 1843 end if 1844 ! 1845 ! ========================================================== 1846 ! ==== Final check on the parameters of the calculation ==== 1847 ! ========================================================== 1848 if ( Bsp%use_coupling>0 .and. ALL(Bsp%algorithm /= [BSE_ALGO_DDIAGO, BSE_ALGO_HAYDOCK]) ) then 1849 ABI_ERROR("Resonant+Coupling is only available with the direct diagonalization or the haydock method.") 1850 end if 1851 1852 ! autoparal section 1853 if (dtset%max_ncpus /=0 .and. dtset%autoparal /=0 ) then 1854 ount = ab_out 1855 ! TODO: 1856 ! nsppol and calculation with coupling! 1857 1858 ! Temporary table needed to estimate memory 1859 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 1860 if (Dtset%usepaw==1) then 1861 do iat=1,Cryst%natom 1862 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 1863 end do 1864 end if 1865 1866 tot_nreh = SUM(BSp%nreh) 1867 work_size = tot_nreh * (tot_nreh + 1) / 2 1868 1869 write(ount,'(a)')"--- !Autoparal" 1870 write(ount,"(a)")'#Autoparal section for Bethe-Salpeter runs.' 1871 1872 write(ount,"(a)") "info:" 1873 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 1874 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 1875 write(ount,"(a,i0)")" nkibz: ",Bsp%nkibz 1876 write(ount,"(a,i0)")" nkbz: ",Bsp%nkbz 1877 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 1878 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 1879 write(ount,"(a,i0)")" lomo_min: ",Bsp%lomo_min 1880 write(ount,"(a,i0)")" humo_max: ",Bsp%humo_max 1881 write(ount,"(a,i0)")" tot_nreh: ",tot_nreh 1882 !write(ount,"(a,i0)")" nbnds: ",Ep%nbnds 1883 1884 ! Wavefunctions are not distributed. We read all the bands 1885 ! from 1 up to Bsp%nbnds because we have to recompute rhor 1886 ! but then we deallocate all the states that are not used for the construction of the e-h 1887 ! before allocating the EXC hamiltonian. Hence we can safely use (humo - lomo + 1) instead of Bsp%nbnds. 1888 !my_nbks = (Bsp%humo - Bsp%lomo +1) * Bsp%nkibz * Dtset%nsppol 1889 1890 ! This one overestimates the memory but it seems to be safer. 1891 my_nbks = Bsp%nbnds * Dtset%nkpt * Dtset%nsppol 1892 1893 ! Memory needed for Fourier components ug. 1894 ug_mem = two*gwpc*Dtset%nspinor*Bsp%npwwfn*my_nbks*b2Mb 1895 1896 ! Memory needed for real space ur. 1897 ur_mem = zero 1898 if (MODULO(Dtset%gwmem,10)==1) then 1899 ur_mem = two*gwpc*Dtset%nspinor*nfftot_osc*my_nbks*b2Mb 1900 end if 1901 1902 ! Memory needed for PAW projections Cprj 1903 cprj_mem = zero 1904 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 1905 1906 wfsmem_mb = ug_mem + ur_mem + cprj_mem 1907 1908 ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI: wavefunctions + W 1909 nonscal_mem = (wfsmem_mb + two*gwpc*BSp%npweps**2*b2Mb) * 1.1_dp 1910 1911 ! List of configurations. 1912 write(ount,"(a)")"configurations:" 1913 do il=1,dtset%max_ncpus 1914 if (il > work_size) cycle 1915 neh_per_proc = work_size / il 1916 neh_per_proc = neh_per_proc + MOD(work_size, il) 1917 eff = (one * work_size) / (il * neh_per_proc) 1918 1919 ! EXC matrix is distributed. 1920 mempercpu_mb = nonscal_mem + two * dpc * neh_per_proc * b2Mb 1921 1922 write(ount,"(a,i0)")" - tot_ncpus: ",il 1923 write(ount,"(a,i0)")" mpi_ncpus: ",il 1924 !write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 1925 write(ount,"(a,f12.9)")" efficiency: ",eff 1926 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 1927 end do 1928 1929 write(ount,'(a)')"..." 1930 1931 ABI_FREE(nlmn_atm) 1932 ABI_ERROR_NODUMP("aborting now") 1933 end if 1934 1935 DBG_EXIT("COLL") 1936 1937 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_ebands<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
1972 subroutine setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh, & 1973 Kmesh_dense,Qmesh_dense,ks_ebands_dense,qp_ebands_dense,Gsph_x,Gsph_c,Vcp_dense,Hdr_wfk_dense,grid,comm) 1974 1975 !Arguments ------------------------------------ 1976 !scalars 1977 integer,intent(in) :: comm 1978 type(dataset_type),intent(in) :: Dtset 1979 type(datafiles_type),intent(inout) :: Dtfil 1980 type(excparam),intent(inout) :: Bsp 1981 type(hdr_type),intent(out) :: Hdr_wfk_dense 1982 type(crystal_t),intent(in) :: Cryst 1983 type(kmesh_t),intent(in) :: Kmesh 1984 type(kmesh_t),intent(out) :: Kmesh_dense,Qmesh_dense 1985 type(ebands_t),intent(out) :: ks_ebands_dense,qp_ebands_dense 1986 type(double_grid_t),intent(out) :: grid 1987 type(vcoul_t),intent(out) :: Vcp_dense 1988 type(gsphere_t),intent(out) :: Gsph_x,Gsph_c 1989 !arrays 1990 1991 !Local variables ------------------------------ 1992 !scalars 1993 integer,parameter :: pertcase0=0,master=0 1994 integer :: bantot_dense,ib,ibtot,ik_ibz,isppol,jj 1995 integer :: nbnds_kss_dense 1996 integer :: spin,hexc_size 1997 integer :: my_rank 1998 integer :: it 1999 integer :: nprocs 2000 integer :: is1,is2,is3,is4 2001 real(dp) :: nelect_hdr_dense 2002 logical,parameter :: remove_inv=.FALSE. 2003 character(len=500) :: msg 2004 character(len=fnlen) :: wfk_fname_dense 2005 integer :: nqlwl 2006 !arrays 2007 integer,allocatable :: npwarr(:) 2008 real(dp),allocatable :: shiftk(:,:) 2009 real(dp),allocatable :: doccde(:),eigen(:),occfact(:) 2010 real(dp),pointer :: energies_p_dense(:,:,:) 2011 complex(dpc),allocatable :: gw_energy(:,:,:) 2012 integer,allocatable :: nbands_temp(:) 2013 integer :: kptrlatt_dense(3,3) 2014 real(dp),allocatable :: qlwl(:,:) 2015 real(dp) :: minmax_tene(2) 2016 2017 !************************************************************************ 2018 2019 DBG_ENTER("COLL") 2020 2021 kptrlatt_dense = zero 2022 2023 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 2024 2025 SELECT CASE(BSp%interp_mode) 2026 CASE (1,2,3,4) 2027 nbnds_kss_dense = -1 2028 wfk_fname_dense = Dtfil%fnameabi_wfkfine 2029 call wrtout(std_out,"BSE Interpolation: will read energies from: "//trim(wfk_fname_dense),"COLL") 2030 2031 if (nctk_try_fort_or_ncfile(wfk_fname_dense, msg) /= 0) then 2032 ABI_ERROR(msg) 2033 end if 2034 2035 Dtfil%fnameabi_wfkfine = wfk_fname_dense 2036 2037 call wfk_read_eigenvalues(wfk_fname_dense,energies_p_dense,Hdr_wfk_dense,comm) 2038 nbnds_kss_dense = MAXVAL(Hdr_wfk_dense%nband) 2039 CASE DEFAULT 2040 ABI_ERROR("Not yet implemented") 2041 END SELECT 2042 2043 nelect_hdr_dense = Hdr_wfk_dense%nelect 2044 2045 if (ABS(Dtset%nelect-nelect_hdr_dense)>tol6) then 2046 write(msg,'(2(a,f8.2))')& 2047 & "File contains ", nelect_hdr_dense," electrons but nelect initialized from input is ",Dtset%nelect 2048 ABI_ERROR(msg) 2049 end if 2050 2051 ! Setup of the k-point list and symmetry tables in the BZ 2052 SELECT CASE(BSp%interp_mode) 2053 CASE (1,2,3,4) 2054 if(Dtset%chksymbreak == 0) then 2055 ABI_MALLOC(shiftk,(3,Dtset%nshiftk)) 2056 kptrlatt_dense(:,1) = BSp%interp_kmult(1)*Dtset%kptrlatt(:,1) 2057 kptrlatt_dense(:,2) = BSp%interp_kmult(2)*Dtset%kptrlatt(:,2) 2058 kptrlatt_dense(:,3) = BSp%interp_kmult(3)*Dtset%kptrlatt(:,3) 2059 do jj = 1,Dtset%nshiftk 2060 shiftk(:,jj) = Bsp%interp_kmult(:)*Dtset%shiftk(:,jj) 2061 end do 2062 call make_mesh(Kmesh_dense,Cryst,Dtset%kptopt,kptrlatt_dense,Dtset%nshiftk,shiftk,break_symmetry=.TRUE.) 2063 ABI_FREE(shiftk) 2064 else 2065 !Initialize Kmesh with no wrapping inside ]-0.5;0.5] 2066 call Kmesh_dense%init(Cryst,Hdr_wfk_dense%nkpt,Hdr_wfk_dense%kptns,Dtset%kptopt) 2067 end if 2068 CASE DEFAULT 2069 ABI_ERROR("Not yet implemented") 2070 END SELECT 2071 2072 ! Init Qmesh 2073 call find_qmesh(Qmesh_dense,Cryst,Kmesh_dense) 2074 2075 call Gsph_c%init(Cryst, 0, ecut=Dtset%ecuteps) 2076 2077 call double_grid_init(Kmesh,Kmesh_dense,Dtset%kptrlatt,BSp%interp_kmult,grid) 2078 2079 BSp%nkibz_interp = Kmesh_dense%nibz !We might allow for a smaller number of points.... 2080 2081 call Kmesh_dense%print("Interpolated K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 2082 call Kmesh_dense%print("Interpolated K-mesh for the wavefunctions",ab_out, 0, "COLL") 2083 2084 if (nbnds_kss_dense < Dtset%nband(1)) then 2085 write(msg,'(2(a,i0),3a,i0)')& 2086 'Interpolated WFK file contains only ', nbnds_kss_dense,' levels instead of ',Dtset%nband(1),' required;',ch10,& 2087 'The calculation will be done with nbands= ',nbnds_kss_dense 2088 ABI_WARNING(msg) 2089 ABI_ERROR("Not supported yet !") 2090 end if 2091 2092 ABI_MALLOC(nbands_temp,(Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol)) 2093 do isppol=1,Hdr_wfk_dense%nsppol 2094 do ik_ibz=1,Hdr_wfk_dense%nkpt 2095 nbands_temp(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) = Dtset%nband(1) 2096 end do 2097 end do 2098 2099 call Gsph_c%extend(Cryst, Dtset%ecutwfn, Gsph_x) 2100 call Gsph_x%print(unit=std_out,prtvol=Dtset%prtvol) 2101 2102 nqlwl=1 2103 ABI_MALLOC(qlwl,(3,nqlwl)) 2104 qlwl(:,nqlwl)= GW_Q0_DEFAULT 2105 2106 ! Compute Coulomb term on the largest G-sphere. 2107 if (Gsph_x%ng > Gsph_c%ng ) then 2108 call Vcp_dense%init(Gsph_x,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%gw_icutcoul,& 2109 Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm) 2110 else 2111 call Vcp_dense%init(Gsph_c,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%gw_icutcoul,& 2112 Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm) 2113 end if 2114 2115 ABI_FREE(qlwl) 2116 2117 bantot_dense=SUM(Hdr_wfk_dense%nband(1:Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol)) 2118 ABI_MALLOC(doccde,(bantot_dense)) 2119 ABI_MALLOC(eigen,(bantot_dense)) 2120 ABI_MALLOC(occfact,(bantot_dense)) 2121 doccde=zero; eigen=zero; occfact=zero 2122 2123 jj=0; ibtot=0 2124 do isppol=1,Hdr_wfk_dense%nsppol 2125 do ik_ibz=1,Hdr_wfk_dense%nkpt 2126 do ib=1,Hdr_wfk_dense%nband(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) 2127 ibtot=ibtot+1 2128 if (ib<=BSP%nbnds) then 2129 jj=jj+1 2130 occfact(jj)=Hdr_wfk_dense%occ(ibtot) 2131 eigen (jj)=energies_p_dense(ib,ik_ibz,isppol) 2132 end if 2133 end do 2134 end do 2135 end do 2136 2137 ABI_FREE(energies_p_dense) 2138 2139 ABI_MALLOC(npwarr,(kmesh_dense%nibz)) 2140 npwarr=BSP%npwwfn 2141 2142 call ebands_init(bantot_dense,ks_ebands_dense,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& 2143 doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,& 2144 Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,& 2145 Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,& 2146 hdr_wfk_dense%cellcharge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, & 2147 hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk) 2148 2149 ABI_FREE(doccde) 2150 ABI_FREE(eigen) 2151 ABI_FREE(npwarr) 2152 2153 ABI_FREE(nbands_temp) 2154 2155 ABI_FREE(occfact) 2156 2157 !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when 2158 ! the occupation scheme for semiconductors is used. 2159 call ebands_update_occ(ks_ebands_dense,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 2160 2161 call ebands_print(ks_ebands_dense,"Interpolated band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol) 2162 2163 call ebands_report_gap(ks_ebands_dense,header="Interpolated KS band structure",unit=std_out,mode_paral="COLL") 2164 2165 BSp%nkbz_interp = Kmesh_dense%nbz 2166 2167 call ebands_copy(ks_ebands_dense,qp_ebands_dense) 2168 2169 SELECT CASE (Bsp%calc_type) 2170 CASE (BSE_HTYPE_RPA_KS) 2171 if (ABS(BSp%mbpt_sciss)>tol6) then 2172 write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies." 2173 call wrtout(std_out,msg) 2174 call ebands_apply_scissors(qp_ebands_dense,BSp%mbpt_sciss) 2175 else 2176 write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]." 2177 call wrtout(std_out,msg) 2178 end if 2179 ! 2180 CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw 2181 ABI_ERROR("Not yet implemented with interpolation !") 2182 CASE (BSE_HTYPE_RPA_QP) 2183 ABI_ERROR("Not implemented error!") 2184 CASE DEFAULT 2185 write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type 2186 ABI_ERROR(msg) 2187 END SELECT 2188 2189 call ebands_report_gap(qp_ebands_dense,header=" Interpolated QP band structure",unit=std_out,mode_paral="COLL") 2190 2191 ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index. 2192 ! FIXME: linewidths not coded. 2193 ABI_MALLOC(gw_energy, (BSp%nbnds,Kmesh_dense%nibz,Dtset%nsppol)) 2194 gw_energy = qp_ebands_dense%eig 2195 2196 ABI_MALLOC(Bsp%nreh_interp,(Hdr_wfk_dense%nsppol)) 2197 Bsp%nreh_interp=zero 2198 2199 call init_transitions(BSp%Trans_interp,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz_interp,Bsp%nbnds,& 2200 & Bsp%nkibz_interp,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,gw_energy,qp_ebands_dense%occ,Kmesh_dense%tab,minmax_tene,& 2201 & Bsp%nreh_interp) 2202 2203 ABI_FREE(gw_energy) 2204 2205 do spin=1,Dtset%nsppol 2206 write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh_interp(spin) 2207 call wrtout(std_out,msg) 2208 end do 2209 2210 if (ANY(Bsp%nreh_interp/=Bsp%nreh_interp(1))) then 2211 write(msg,'(a,(i0))')" BSE code does not support different number of transitions for the two spin channels",Bsp%nreh 2212 ABI_ERROR(msg) 2213 end if 2214 ! 2215 ! Create transition table vcks2t 2216 is1=BSp%lomo_min;is2=BSp%homo_max;is3=BSp%lumo_min;is4=BSp%humo_max 2217 ABI_MALLOC(Bsp%vcks2t_interp, (is1:is2,is3:is4,BSp%nkbz_interp,Dtset%nsppol)) 2218 Bsp%vcks2t_interp = 0 2219 2220 do spin=1,Dtset%nsppol 2221 do it=1,BSp%nreh_interp(spin) 2222 BSp%vcks2t_interp(BSp%Trans_interp(it,spin)%v,BSp%Trans_interp(it,spin)%c, BSp%Trans_interp(it,spin)%k,spin) = it 2223 end do 2224 end do 2225 2226 hexc_size = SUM(Bsp%nreh_interp); if (Bsp%use_coupling>0) hexc_size=2*hexc_size 2227 if (Bsp%nstates_interp<=0) then 2228 Bsp%nstates_interp=hexc_size 2229 else 2230 if (Bsp%nstates_interp>hexc_size) then 2231 Bsp%nstates_interp=hexc_size 2232 write(msg,'(2(a,i0),2a)')& 2233 "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates_interp,ch10,& 2234 "the number of excitonic states nstates has been modified" 2235 ABI_WARNING(msg) 2236 end if 2237 end if 2238 2239 DBG_EXIT("COLL") 2240 2241 end subroutine setup_bse_interp