TABLE OF CONTENTS
ABINIT/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.
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (M.Giantomassi, L. Reining, V. Olevano, F. Sottile, S. Albrecht, G. Onida) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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.
PARENTS
driver
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.
CHILDREN
bs_parameters_free,chkpawovlp,crystal_free,denfgr,destroy_mpi_enreg double_grid_free,ebands_free,ebands_update_occ,energies_init eprenorms_free,exc_build_ham,exc_den,exc_diago_driver exc_haydock_driver,fourdp,get_gftt,getph,gsph_free,hdr_free init_distribfft_seq,initmpi_seq,kmesh_free,metric,mkdenpos,mkrdim nhatgrid,paw_an_free,paw_an_init,paw_an_nullify,paw_gencond,paw_ij_free paw_ij_init,paw_ij_nullify,pawdenpot,pawdij,pawfgr_destroy,pawfgr_init pawfgrtab_free,pawfgrtab_init,pawhur_free,pawhur_init,pawinit,pawmknhat pawnabla_init,pawprt,pawpuxinit,pawpwff_free,pawpwff_init pawrhoij_alloc,pawrhoij_copy,pawrhoij_free,pawtab_get_lsize pawtab_print,print_ngfft,prtrhomxmn,pspini,rdqps,rotate_fft_mesh screen_free,screen_init,screen_nullify,setsymrhoij,setup_bse setup_bse_interp,setvtr,symdij,test_charge,timab,vcoul_free,wfd_free wfd_init,wfd_mkrho,wfd_print,wfd_read_wfk,wfd_reset_ur_cprj,wfd_rotate wfd_test_ortho,wfd_wave_free,wrtout,xmpi_bcast
SOURCE
87 #if defined HAVE_CONFIG_H 88 #include "config.h" 89 #endif 90 91 #include "abi_common.h" 92 93 subroutine bethe_salpeter(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim,xred) 94 95 use defs_basis 96 use defs_datatypes 97 use defs_abitypes 98 use defs_wvltypes 99 use m_bs_defs 100 use m_profiling_abi 101 use m_xmpi 102 use m_errors 103 use m_screen 104 use m_nctk 105 #ifdef HAVE_NETCDF 106 use netcdf 107 #endif 108 use m_hdr 109 110 use m_fstrings, only : strcat, sjoin, endswith 111 use m_io_tools, only : file_exists, iomode_from_fname 112 use m_geometry, only : mkrdim, metric 113 use m_mpinfo, only : destroy_mpi_enreg 114 use m_fftcore, only : print_ngfft 115 use m_fft_mesh, only : rotate_FFT_mesh, get_gftt 116 use m_crystal, only : crystal_t, crystal_free 117 use m_crystal_io, only : crystal_ncwrite 118 use m_bz_mesh, only : kmesh_t, kmesh_free, make_path 119 use m_double_grid, only : double_grid_t, double_grid_free 120 use m_ebands, only : ebands_update_occ, ebands_free, ebands_copy, get_valence_idx 121 use m_gsphere, only : gsphere_t, gsph_free 122 use m_vcoul, only : vcoul_t, vcoul_free 123 use m_qparticles, only : rdqps !, show_QP , rdgw 124 use m_wfd, only : wfd_init, wfd_free, wfd_print, wfd_t, wfd_test_ortho,& 125 & wfd_read_wfk, wfd_wave_free, wfd_rotate, wfd_reset_ur_cprj 126 use m_energies, only : energies_type, energies_init 127 use m_haydock, only : exc_haydock_driver 128 use m_exc_diago, only : exc_diago_driver 129 use m_eprenorms, only : eprenorms_t, eprenorms_free 130 use m_pawang, only : pawang_type 131 use m_pawrad, only : pawrad_type 132 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 133 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 134 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 135 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init 136 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy,& 137 & pawrhoij_free, pawrhoij_get_nspden, symrhoij 138 use m_pawdij, only : pawdij, symdij 139 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 140 use m_pawhr, only : pawhur_t, pawhur_free, pawhur_init 141 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 142 use m_paw_dmft, only : paw_dmft_type 143 use m_kg, only : getph 144 145 !This section has been created automatically by the script Abilint (TD). 146 !Do not modify the following lines by hand. 147 #undef ABI_FUNC 148 #define ABI_FUNC 'bethe_salpeter' 149 use interfaces_14_hidewrite 150 use interfaces_18_timing 151 use interfaces_41_xc_lowlevel 152 use interfaces_51_manage_mpi 153 use interfaces_53_ffts 154 use interfaces_64_psp 155 use interfaces_65_paw 156 use interfaces_67_common 157 use interfaces_69_wfdesc 158 use interfaces_71_bse 159 !End of the abilint section 160 161 implicit none 162 163 !Arguments ------------------------------------ 164 !scalars 165 character(len=6),intent(in) :: codvsn 166 type(datafiles_type),intent(inout) :: Dtfil 167 type(dataset_type),intent(inout) :: Dtset 168 type(pawang_type),intent(inout) :: Pawang 169 type(pseudopotential_type),intent(inout) :: Psps 170 !arrays 171 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,Dtset%natom) 172 type(pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw) 173 type(pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw) 174 175 !Local variables ------------------------------ 176 !scalars 177 integer,parameter :: tim_fourdp0=0,level=40,ipert0=0,idir0=0,cplex1=1,master=0,option1=1 178 integer :: band,spin,ik_ibz,mqmem,iwarn !ii, 179 integer :: has_dijU,has_dijso,gnt_option,iomode 180 integer :: ik_bz,mband 181 integer :: choice 182 integer :: ider !,ierr 183 integer :: usexcnhat,nfft_osc,mgfft_osc 184 integer :: isym,izero 185 integer :: optcut,optgr0,optgr1,optgr2,option,optrad,optrhoij,psp_gencond 186 integer :: ngrvdw,nhatgrdim,nkxc1,nprocs,nspden_rhoij,nzlmopt,ifft 187 integer :: my_rank,rhoxsp_method,comm 188 integer :: mgfftf,spin_opt,which_fixed 189 integer :: nscf,nbsc,nkxc,n3xccc 190 integer :: nfftf,nfftf_tot,nfftot_osc,my_minb,my_maxb 191 integer :: optene,moved_atm_inside,moved_rhor,initialized,istep,ierr 192 real(dp) :: ucvol,drude_plsmf,ecore,ecut_eff,ecutdg_eff,opt_ecut,norm 193 real(dp) :: gsqcutc_eff,gsqcutf_eff,gsqcut_shp 194 real(dp) :: compch_fft,compch_sph,gsq_osc 195 real(dp) :: vxcavg 196 logical :: iscompatibleFFT,paw_add_onsite,call_pawinit 197 character(len=500) :: msg 198 character(len=fnlen) :: wfk_fname,w_fname 199 type(Pawfgr_type) :: Pawfgr 200 type(excfiles) :: BS_files 201 type(excparam) :: BSp 202 type(paw_dmft_type) :: Paw_dmft 203 type(MPI_type) :: MPI_enreg_seq 204 type(crystal_t) :: Cryst 205 type(kmesh_t) :: Kmesh,Qmesh 206 type(gsphere_t) :: Gsph_x,Gsph_c 207 type(gsphere_t) :: Gsph_x_dense,Gsph_c_dense 208 type(Hdr_type) :: Hdr_wfk,Hdr_bse 209 type(ebands_t) :: KS_BSt,QP_BSt 210 type(Energies_type) :: KS_energies 211 type(vcoul_t) :: Vcp 212 type(wfd_t) :: Wfd 213 type(screen_t) :: W 214 type(screen_info_t) :: W_info 215 type(wvl_data) :: wvl 216 type(kmesh_t) :: Kmesh_dense,Qmesh_dense 217 type(Hdr_type) :: Hdr_wfk_dense 218 type(ebands_t) :: KS_BSt_dense,QP_BSt_dense 219 type(wfd_t) :: Wfd_dense 220 type(double_grid_t) :: grid 221 type(vcoul_t) :: Vcp_dense 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_lda_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 MSG_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,msg,'COLL') 271 call wrtout(ab_out,msg,'COLL') 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 MSG_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,msg,'COLL') 286 call wrtout(ab_out,msg,'COLL') 287 288 !=== Some variables need to be initialized/nullify at start === 289 call energies_init(KS_energies) 290 usexcnhat=0 291 call mkrdim(acell,rprim,rprimd) 292 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 293 ! 294 !=== Define FFT grid(s) sizes === 295 !* Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the 296 !(usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 297 !See also NOTES in the comments at the beginning of this file. 298 !NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 299 300 !TODO Recheck getng, should use same trick as that used in screening and sigma. 301 call pawfgr_init(Pawfgr,Dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 302 & gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=gmet,k0=k0) 303 304 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 305 nfftf_tot=PRODUCT(ngfftf(1:3)) 306 ! 307 ! * Fake MPI_type for the sequential part. 308 call initmpi_seq(MPI_enreg_seq) 309 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 310 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 311 ! 312 ! =========================================== 313 ! === Open and read pseudopotential files === 314 ! =========================================== 315 call pspini(Dtset,Dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,Pawrad,Pawtab,Psps,rprimd,comm_mpi=comm) 316 317 ! === Initialization of basic objects including the BSp structure that defines the parameters of the run === 318 call setup_bse(codvsn,acell,rprim,ngfftf,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,& 319 & Cryst,Kmesh,Qmesh,KS_BSt,QP_BSt,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,wvl%descr) 320 321 if (BSp%use_interp) then 322 call setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh,Kmesh_dense,& 323 & Qmesh_dense,KS_BSt_dense,QP_BSt_dense,Gsph_x_dense,Gsph_c_dense,& 324 & Vcp_dense,Hdr_wfk_dense,ngfftf,grid,comm) 325 end if 326 327 !call timab(652,2,tsec) ! setup_bse 328 329 nfftot_osc=PRODUCT(ngfft_osc(1:3)) 330 nfft_osc =nfftot_osc !no FFT // 331 mgfft_osc =MAXVAL(ngfft_osc(1:3)) 332 333 call print_ngfft(ngfft_osc,header='FFT mesh used for oscillator strengths') 334 335 !TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 336 KS_energies%e_corepsp=ecore/Cryst%ucvol 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_DT_MALLOC(KS_Pawrhoij,(Cryst%natom)) 345 nspden_rhoij=pawrhoij_get_nspden(Dtset%nspden,Dtset%nspinor,Dtset%pawspnorb) 346 call pawrhoij_alloc(KS_Pawrhoij,Dtset%pawcpxocc,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(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%xclevel,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 375 ! * Evaluate <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for the long wavelength limit. 376 ! TODO solve problem with memory leak and clean this part as well as the associated flag 377 call pawnabla_init(Psps%mpsang,Cryst%ntypat,Pawrad,Pawtab) 378 379 call setsymrhoij(gprimd,Pawang%l_max-1,Cryst%nsym,Dtset%pawprtvol,Cryst%rprimd,Cryst%symrec,Pawang%zarot) 380 381 ! Initialize and compute data for LDA+U 382 if (Dtset%usepawu>0.or.Dtset%useexexch>0) then 383 Paw_dmft%use_dmft=Dtset%usedmft 384 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 385 & Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,Cryst%ntypat,Pawang,Dtset%pawprtvol,& 386 & Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu) 387 MSG_ERROR("BS equation with LDA+U not completely coded") 388 end if 389 if (my_rank == master) call pawtab_print(Pawtab) 390 391 ! Get Pawrhoij from the header of the WFK file. 392 call pawrhoij_copy(Hdr_wfk%pawrhoij,KS_Pawrhoij) 393 394 ! Re-symmetrize symrhoij === 395 ! this call leads to a SIGFAULT, likely some pointer is not initialized correctly 396 choice=1; optrhoij=1 397 ! call symrhoij(KS_Pawrhoij,KS_Pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,& 398 ! & Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij,Pawang,Dtset%pawprtvol,Pawtab,& 399 ! & Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 400 401 ! Evaluate form factor of radial part of phi.phj-tphi.tphj === 402 rhoxsp_method=1 ! Arnaud-Alouani 403 if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc 404 405 ABI_MALLOC(gfft_osc,(3,nfftot_osc)) 406 call get_gftt(ngfft_osc,k0,gmet,gsq_osc,gfft_osc) 407 ABI_FREE(gfft_osc) 408 409 ! Set up q grids, make qmax 20% larger than largest expected: 410 ABI_MALLOC(nq_spl,(Psps%ntypat)) 411 ABI_MALLOC(qmax,(Psps%ntypat)) 412 nq_spl = Psps%mqgrid_ff 413 qmax = SQRT(gsq_osc)*1.2d0 ! qmax=Psps%qgrid_ff(Psps%mqgrid_ff) 414 ABI_DT_MALLOC(Paw_pwff,(Psps%ntypat)) 415 416 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 417 418 ABI_FREE(nq_spl) 419 ABI_FREE(qmax) 420 ! 421 ! Variables/arrays related to the fine FFT grid === 422 ABI_MALLOC(ks_nhat,(nfftf,Dtset%nspden)) 423 ks_nhat=zero 424 ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom)) 425 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 426 call pawfgrtab_init(Pawfgrtab,cplex1,l_size_atm,Dtset%nspden,Dtset%typat) 427 ABI_FREE(l_size_atm) 428 compch_fft=greatest_real 429 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 430 ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation) 431 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 432 write(msg,'(a,i2)')' bethe_salpeter : using usexcnhat = ',usexcnhat 433 call wrtout(std_out,msg,'COLL') 434 ! 435 ! Identify parts of the rectangular grid where the density has to be calculated === 436 optcut=0; optgr0=Dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-Dtset%pawstgylm 437 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 438 439 call nhatgrid(Cryst%atindx1,gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 440 optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 441 else 442 ABI_DT_MALLOC(Paw_pwff,(0)) 443 end if !End of PAW Initialization 444 445 ! Consistency check and additional stuff done only for GW with PAW. 446 if (Dtset%usepaw==1) then 447 if (Dtset%ecutwfn < Dtset%ecut) then 448 write(msg,"(5a)")& 449 & "WARNING - ",ch10,& 450 & " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 451 & " an excessive truncation of the planewave basis set can lead to unphysical results." 452 call wrtout(ab_out,msg,'COLL') 453 end if 454 455 ABI_CHECK(Dtset%usedmft==0,"DMFT + BSE not allowed") 456 ABI_CHECK(Dtset%useexexch==0,"LEXX + BSE not allowed") 457 end if 458 459 ! Allocate these arrays anyway, since they are passed to subroutines. 460 if (.not.allocated(ks_nhat)) then 461 ABI_MALLOC(ks_nhat,(nfftf,0)) 462 end if 463 464 !================================================== 465 !==== Read KS band structure from the KSS file ==== 466 !================================================== 467 468 !* Initialize wave function handler, allocate wavefunctions. 469 my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds 470 ABI_MALLOC(nband,(Kmesh%nibz,Dtset%nsppol)) 471 nband=mband 472 473 !At present, no memory distribution, each node has the full set of states. 474 ABI_MALLOC(bks_mask,(mband,Kmesh%nibz,Dtset%nsppol)) 475 bks_mask=.TRUE. 476 477 ABI_MALLOC(keep_ur,(mband,Kmesh%nibz,Dtset%nsppol)) 478 keep_ur=.FALSE. 479 if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE. 480 481 !opt_ecut=zero 482 opt_ecut=Dtset%ecutwfn 483 484 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,BSp%npwwfn,mband,nband,Kmesh%nibz,Dtset%nsppol,bks_mask,& 485 & Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,ngfft_osc,& 486 & Gsph_x%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut) 487 488 ABI_FREE(bks_mask) 489 ABI_FREE(nband) 490 ABI_FREE(keep_ur) 491 492 call wfd_print(Wfd,header="Wavefunctions used to construct the e-h basis set",mode_paral='PERS') 493 494 call timab(651,2,tsec) ! bse(Init1) 495 call timab(653,1,tsec) ! bse(rdkss) 496 497 call wfd_read_wfk(Wfd,wfk_fname,iomode_from_fname(wfk_fname)) 498 499 ! This test has been disabled (too expensive!) 500 if (.False.) call wfd_test_ortho(Wfd,Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 501 502 call timab(653,2,tsec) ! bse(rdkss) 503 call timab(655,1,tsec) ! bse(mkrho) 504 505 !TODO : check the consistency of Wfd with Wfd_dense !!! 506 if (BSp%use_interp) then 507 ! Initialize wave function handler, allocate wavefunctions. 508 my_minb=1; my_maxb=BSp%nbnds; mband=BSp%nbnds 509 ABI_MALLOC(nband,(Kmesh_dense%nibz,Dtset%nsppol)) 510 nband=mband 511 512 ! At present, no memory distribution, each node has the full set of states. 513 ! albeit we allocate only the states that are used. 514 ABI_MALLOC(bks_mask,(mband,Kmesh_dense%nibz,Dtset%nsppol)) 515 bks_mask=.False. 516 do spin=1,Bsp%nsppol 517 bks_mask(Bsp%lomo_spin(spin):Bsp%humo_spin(spin),:,spin) = .True. 518 end do 519 !bks_mask=.TRUE. 520 521 ABI_MALLOC(keep_ur,(mband,Kmesh_dense%nibz,Dtset%nsppol)) 522 keep_ur=.FALSE. 523 if (MODULO(Dtset%gwmem,10)==1) keep_ur = .TRUE. 524 525 ! opt_ecut=zero 526 opt_ecut=Dtset%ecutwfn 527 528 call wfd_init(Wfd_dense,Cryst,Pawtab,Psps,keep_ur,Dtset%paral_kgb,BSp%npwwfn,mband,nband,Kmesh_dense%nibz,Dtset%nsppol,& 529 & bks_mask,Dtset%nspden,Dtset%nspinor,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,ngfft_osc,& 530 & Gsph_x_dense%gvec,Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,opt_ecut=opt_ecut) 531 532 ABI_FREE(bks_mask) 533 ABI_FREE(nband) 534 ABI_FREE(keep_ur) 535 536 call wfd_print(Wfd_dense,header="Wavefunctions on the dense K-mesh used for interpolation",mode_paral='PERS') 537 538 iomode = iomode_from_fname(dtfil%fnameabi_wfkfine) 539 call wfd_read_wfk(Wfd_dense,Dtfil%fnameabi_wfkfine,iomode) 540 541 ! This test has been disabled (too expensive!) 542 if (.False.) call wfd_test_ortho(Wfd_dense,Cryst,Pawtab,unit=std_out,mode_paral="COLL") 543 end if 544 545 !=== Calculate the FFT index of $(R^{-1}(r-\tau))$ === 546 !* S=\transpose R^{-1} and k_BZ = S k_IBZ 547 !* irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 548 ABI_MALLOC(irottb,(nfftot_osc,Cryst%nsym)) 549 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfft_osc,irottb,iscompatibleFFT) 550 551 ABI_MALLOC(ktabr,(nfftot_osc,Kmesh%nbz)) 552 do ik_bz=1,Kmesh%nbz 553 isym=Kmesh%tabo(ik_bz) 554 do ifft=1,nfftot_osc 555 ktabr(ifft,ik_bz)=irottb(ifft,isym) 556 end do 557 end do 558 ABI_FREE(irottb) 559 ! 560 !=========================== 561 !=== COMPUTE THE DENSITY === 562 !=========================== 563 !* Evaluate Planewave part (complete charge in case of NC pseudos). 564 ! 565 ABI_MALLOC(ks_rhor,(nfftf,Wfd%nspden)) 566 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,KS_BSt,ngfftf,nfftf,ks_rhor) 567 ! 568 !=== Additional computation for PAW === 569 nhatgrdim=0 570 if (Dtset%usepaw==1) then 571 ! 572 ! Calculate the compensation charge nhat. 573 if (Dtset%xclevel==2) nhatgrdim=usexcnhat*Dtset%pawnhatxc 574 ider=2*nhatgrdim; izero=0; qphon(:)=zero 575 if (nhatgrdim>0) then 576 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3)) 577 end if 578 579 call pawmknhat(compch_fft,cplex1,ider,idir0,ipert0,izero,Cryst%gprimd,& 580 & Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 581 & Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,qphon,Cryst%rprimd,& 582 & Cryst%ucvol,Dtset%usewvl,Cryst%xred) 583 584 ! Evaluate onsite energies, potentials, densities === 585 ! * Initialize variables/arrays related to the PAW spheres. 586 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 587 ABI_DT_MALLOC(KS_paw_ij,(Cryst%natom)) 588 call paw_ij_nullify(KS_paw_ij) 589 590 has_dijso=Dtset%pawspnorb 591 has_dijU=Dtset%usepawu 592 has_dijso=Dtset%pawspnorb 593 has_dijU=Dtset%usepawu 594 595 call paw_ij_init(KS_paw_ij,cplex1,Dtset%nspinor,Dtset%nsppol,& 596 & Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 597 & has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=0,has_dijxc_hat=0,has_dijxc_val=0,& 598 & has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1) 599 600 ABI_DT_MALLOC(KS_paw_an,(Cryst%natom)) 601 call paw_an_nullify(KS_paw_an) 602 603 nkxc1=0 604 call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,Dtset%nspden,& 605 & cplex1,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=0) 606 607 ! Calculate onsite vxc with and without core charge === 608 nzlmopt=-1; option=0; compch_sph=greatest_real 609 610 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,& 611 & Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 612 & Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,& 613 & Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,& 614 & Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Cryst%ucvol,Psps%znuclpsp) 615 end if !PAW 616 617 if (.not.allocated(ks_nhatgr)) then 618 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,0)) 619 end if 620 621 call test_charge(nfftf,KS_BSt%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,& 622 & Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 623 ! 624 ! === For PAW, add the compensation charge on the FFT mesh, then get rho(G) === 625 if (Dtset%usepaw==1) ks_rhor(:,:)=ks_rhor(:,:)+ks_nhat(:,:) 626 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol) 627 628 ABI_MALLOC(ks_rhog,(2,nfftf)) 629 630 call fourdp(1,ks_rhog,ks_rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Dtset%paral_kgb,tim_fourdp0) 631 632 call timab(655,2,tsec) ! bse(mkrho) 633 ! 634 ! 635 ! The following steps have been gathered in the setvtr routine: 636 ! - get Ewald energy and Ewald forces 637 ! - compute local ionic pseudopotential vpsp 638 ! - eventually compute 3D core electron density xccc3d 639 ! - eventually compute vxc and vhartr 640 ! - set up ks_vtrial 641 ! 642 ! ******************************************************************* 643 ! **** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 644 ! ******************************************************************* 645 ngrvdw=0 646 ABI_MALLOC(grvdw,(3,ngrvdw)) 647 ABI_MALLOC(grchempottn,(3,Cryst%natom)) 648 ABI_MALLOC(grewtn,(3,Cryst%natom)) 649 nkxc=0 650 !if (Wfd%nspden==1) nkxc=2 651 !if (Wfd%nspden>=2) nkxc=3 ! check GGA and spinor, quite a messy part!!! 652 ABI_MALLOC(kxc,(nfftf,nkxc)) 653 654 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 655 ABI_MALLOC(xccc3d,(n3xccc)) 656 ABI_MALLOC(ks_vhartr,(nfftf)) 657 ABI_MALLOC(ks_vtrial,(nfftf,Wfd%nspden)) 658 ABI_MALLOC(vpsp,(nfftf)) 659 ABI_MALLOC(ks_vxc,(nfftf,Wfd%nspden)) 660 661 optene=4; moved_atm_inside=0; moved_rhor=0; initialized=1; istep=1 662 ! 663 !=== Compute structure factor phases and large sphere cut-off === 664 !WARNING cannot use Dtset%mgfft, this has to be checked better 665 !mgfft=MAXVAL(ngfftc(:)) 666 !allocate(ph1d(2,3*(2*mgfft+1)*Cryst%natom),ph1df(2,3*(2*mgfftf+1)*Cryst%natom)) 667 write(std_out,*)' CHECK ',Dtset%mgfftdg,mgfftf 668 if (Dtset%mgfftdg/=mgfftf) then 669 ! write(std_out,*)"WARNING Dtset%mgfftf /= mgfftf" 670 ! write(std_out,*)'HACKING Dtset%mgfftf' 671 ! Dtset%mgfftdg=mgfftf 672 end if 673 ABI_MALLOC(ph1d,(2,3*(2*Dtset%mgfft+1)*Cryst%natom)) 674 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*Cryst%natom)) 675 676 call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred) 677 678 if (Psps%usepaw==1.and.Pawfgr%usefinegrid==1) then 679 call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred) 680 else 681 ph1df(:,:)=ph1d(:,:) 682 end if 683 684 ABI_FREE(ph1d) 685 686 call setvtr(Cryst%atindx1,Dtset,KS_energies,Cryst%gmet,Cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 687 & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 688 & Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 689 & optene,Pawrad,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,& 690 & Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,wvl,xccc3d,Cryst%xred) 691 692 ABI_FREE(ph1df) 693 ABI_FREE(vpsp) 694 695 ! ============================ 696 ! ==== Compute KS PAW Dij ==== 697 ! ============================ 698 if (Wfd%usepaw==1) then 699 _IBM6("Another silly write for IBM6") 700 call timab(561,1,tsec) 701 ! 702 ! Calculate the unsymmetrized Dij. 703 call pawdij(cplex1,Dtset%enunit,Cryst%gprimd,ipert0,& 704 & Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 705 & Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,& 706 & Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 707 & k0,Dtset%spnorbscl,Cryst%ucvol,dtset%charge,ks_vtrial,ks_vxc,Cryst%xred,& 708 & nucdipmom=Dtset%nucdipmom) 709 ! 710 ! Symmetrize KS Dij 711 call symdij(Cryst%gprimd,Cryst%indsym,ipert0,& 712 & Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,0,KS_paw_ij,Pawang,& 713 & Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 714 ! 715 ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 716 call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab) 717 call timab(561,2,tsec) 718 end if 719 720 ABI_FREE(kxc) 721 ABI_FREE(xccc3d) 722 ABI_FREE(grchempottn) 723 ABI_FREE(grewtn) 724 ABI_FREE(grvdw) 725 726 !=== QP_BSt stores energies and occ. used for the calculation === 727 !* Initialize QP_BSt with KS values. 728 !* In case of SC update QP_BSt using the QPS file. 729 ABI_MALLOC(qp_rhor,(nfftf,Dtset%nspden)) 730 qp_rhor=ks_rhor 731 732 ! AE density used for the model dielectric function. 733 ABI_MALLOC(qp_aerhor,(nfftf,Dtset%nspden)) 734 qp_aerhor=ks_rhor 735 736 ! PAW: Compute AE rhor. Under testing 737 if (Wfd%usepaw==1 .and. BSp%mdlf_type/=0) then 738 MSG_WARNING("Entering qp_aerhor with PAW") 739 740 ABI_MALLOC(qp_rhor_paw ,(nfftf,Wfd%nspden)) 741 ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden)) 742 ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden)) 743 744 ABI_MALLOC(qp_nhat,(nfftf,Wfd%nspden)) 745 qp_nhat = ks_nhat 746 ! TODO: I pass KS_pawrhoij instead of QP_pawrhoij but in the present version there's no difference. 747 748 call denfgr(Cryst%atindx1,Cryst%gmet,Wfd%comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,& 749 & Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_pawrhoij,Pawtab,Dtset%prtvol,& 750 & qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 751 752 norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3)) 753 write(msg,'(a,F8.4)') ' QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm 754 call wrtout(std_out,msg,'COLL') 755 write(std_out,*)"MAX", MAXVAL(qp_rhor_paw(:,1)) 756 write(std_out,*)"MIN", MINVAL(qp_rhor_paw(:,1)) 757 758 ABI_FREE(qp_nhat) 759 ABI_FREE(qp_rhor_n_one) 760 ABI_FREE(qp_rhor_nt_one) 761 762 ! Use ae density for the model dielectric function. 763 iwarn=0 764 call mkdenpos(iwarn,nfftf,Wfd%nspden,option1,qp_rhor_paw,dtset%xc_denpos) 765 qp_aerhor = qp_rhor_paw 766 ABI_FREE(qp_rhor_paw) 767 end if 768 769 !call copy_bandstructure(KS_BSt,QP_BSt) 770 771 if (.FALSE.) then 772 ! $ m_lda_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}> $ 773 ABI_MALLOC(m_lda_to_qp,(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 774 m_lda_to_qp=czero 775 do spin=1,Wfd%nsppol 776 do ik_ibz=1,Wfd%nkibz 777 do band=1,Wfd%nband(ik_ibz,spin) 778 m_lda_to_qp(band,band,ik_ibz,spin)=cone ! Initialize the QP amplitudes with KS wavefunctions. 779 end do 780 end do 781 end do 782 ! 783 ! * Now read m_lda_to_qp and update the energies in QP_BSt. 784 ! TODO switch on the renormalization of n in sigma. 785 ABI_MALLOC(prev_rhor,(nfftf,Wfd%nspden)) 786 ABI_DT_MALLOC(prev_Pawrhoij,(Cryst%natom*Wfd%usepaw)) 787 788 call rdqps(QP_BSt,Dtfil%fnameabi_qps,Wfd%usepaw,Wfd%nspden,1,nscf,& 789 & nfftf,ngfftf,Cryst%ucvol,Wfd%paral_kgb,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_lda_to_qp,prev_rhor,prev_Pawrhoij) 790 791 ABI_FREE(prev_rhor) 792 if (Psps%usepaw==1.and.nscf>0) then 793 call pawrhoij_free(prev_pawrhoij) 794 end if 795 ABI_DT_FREE(prev_pawrhoij) 796 ! 797 ! if (nscf>0.and.wfd_iam_master(Wfd)) then ! Print the unitary transformation on std_out. 798 ! call show_QP(QP_BSt,m_lda_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp) 799 ! end if 800 ! 801 ! === Compute QP wfg as linear combination of KS states === 802 ! * Wfd%ug is modified inside calc_wf_qp 803 ! * For PAW, update also the on-site projections. 804 ! * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz 805 ! TODO here we should use nbsc instead of nbnds 806 807 call wfd_rotate(Wfd,Cryst,m_lda_to_qp) 808 809 ABI_FREE(m_lda_to_qp) 810 ! 811 ! === Reinit the storage mode of Wfd as ug have been changed === 812 ! * Update also the wavefunctions for GW corrections on each processor 813 call wfd_reset_ur_cprj(Wfd) 814 815 call wfd_test_ortho(Wfd,Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 816 817 ! Compute QP occupation numbers. 818 call wrtout(std_out,'bethe_salpeter: calculating QP occupation numbers','COLL') 819 820 call ebands_update_occ(QP_BSt,Dtset%spinmagntarget,prtvol=0) 821 ABI_MALLOC(qp_vbik,(QP_BSt%nkpt,QP_BSt%nsppol)) 822 qp_vbik(:,:) = get_valence_idx(QP_BSt) 823 ABI_FREE(qp_vbik) 824 825 call wfd_mkrho(Wfd,Cryst,Psps,Kmesh,QP_BSt,ngfftf,nfftf,qp_rhor) 826 end if 827 828 ABI_MALLOC(qp_rhog,(2,nfftf)) 829 call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,ngfftf,Wfd%paral_kgb,0) 830 831 ! States up to lomo-1 are useless now since only the states close to the gap are 832 ! needed to construct the EXC Hamiltonian. Here we deallocate the wavefunctions 833 ! to make room for the excitonic Hamiltonian that is allocated in exc_build_ham. 834 ! and for the screening that is allocated below. 835 ! Hereafter bands from 1 up to lomo-1 and all bands above humo+1 should not be accessed! 836 ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 837 838 bks_mask=.FALSE. 839 do spin=1,Bsp%nsppol 840 if (Bsp%lomo_spin(spin)>1) bks_mask(1:Bsp%lomo_spin(spin)-1,:,spin)=.TRUE. 841 if (Bsp%humo_spin(spin)+1<=Wfd%mband) bks_mask(Bsp%humo_spin(spin)+1:,:,spin)=.TRUE. 842 end do 843 call wfd_wave_free(Wfd,"All",bks_mask) 844 ABI_FREE(bks_mask) 845 ! 846 ! ================================================================ 847 ! Build the screened interaction W in the irreducible wedge. 848 ! * W(q,G1,G2) = vc^{1/2} (q,G1) e^{-1}(q,G1,G2) vc^{1/2) (q,G2) 849 ! * Use Coulomb term for q-->0, 850 ! * Only the first small Q is used, shall we average if nqlwl>1? 851 ! ================================================================ 852 ! TODO clean this part and add an option to retrieve a single frequency to save memory. 853 call timab(654,1,tsec) ! bse(rdmkeps^-1) 854 855 call screen_nullify(W) 856 if (BSp%use_coulomb_term) then ! Init W. 857 ! Incore or out-of-core solution? 858 mqmem=0; if (Dtset%gwmem/10==1) mqmem=Qmesh%nibz 859 860 W_info%invalid_freq = Dtset%gw_invalid_freq 861 W_info%mat_type = MAT_INV_EPSILON 862 !W_info%mat_type = MAT_W 863 !W_info%vtx_family 864 !W_info%ixc 865 !W_info%use_ada 866 W_info%use_mdf = BSp%mdlf_type 867 !W_info%use_ppm 868 !W_info%vtx_test 869 !W_info%wint_method 870 ! 871 !W_info%ada_kappa 872 W_info%eps_inf = BSp%eps_inf 873 !W_info%drude_plsmf 874 875 call screen_init(W,W_Info,Cryst,Qmesh,Gsph_c,Vcp,w_fname,mqmem,Dtset%npweps,& 876 & Dtset%iomode,ngfftf,nfftf_tot,Wfd%nsppol,Wfd%nspden,qp_aerhor,Wfd%prtvol,Wfd%comm) 877 end if 878 call timab(654,2,tsec) ! bse(rdmkeps^-1) 879 ! 880 ! ============================================= 881 ! ==== Build of the excitonic Hamiltonian ===== 882 ! ============================================= 883 ! 884 call timab(656,1,tsec) ! bse(mkexcham) 885 886 call exc_build_ham(BSp,BS_files,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,& 887 & Wfd,W,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff) 888 ! 889 ! Free Em1 to make room for the full excitonic Hamiltonian. 890 call screen_free(W) 891 892 call timab(656,2,tsec) ! bse(mkexcham) 893 ! 894 ! ========================================= 895 ! ==== Macroscopic dielectric function ==== 896 ! ========================================= 897 call timab(657,1,tsec) ! bse(mkexceps) 898 ! 899 ! First deallocate the internal %ur buffers to make room for the excitonic Hamiltonian. 900 call timab(658,1,tsec) ! bse(wfd_wave_free) 901 ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 902 bks_mask=.TRUE. 903 call wfd_wave_free(Wfd,"Real_space",bks_mask) 904 ABI_FREE(bks_mask) 905 call timab(658,2,tsec) ! bse(wfd_wave_free) 906 907 ! Compute the commutator [r,Vu] (PAW only). 908 ABI_DT_MALLOC(HUr,(Cryst%natom*Wfd%usepaw)) 909 910 call timab(659,1,tsec) ! bse(make_pawhur_t) 911 if (Bsp%inclvkb/=0 .and. Wfd%usepaw==1 .and. Dtset%usepawu/=0) then !TODO here I need KS_Paw_ij 912 MSG_WARNING("Commutator for LDA+U not tested") 913 call pawhur_init(hur,Wfd%nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,KS_Paw_ij) 914 end if 915 call timab(659,2,tsec) ! bse(make_pawhur_t) 916 917 select case (BSp%algorithm) 918 case (BSE_ALGO_NONE) 919 MSG_COMMENT("Skipping solution of the BSE equation") 920 921 case (BSE_ALGO_DDIAGO, BSE_ALGO_CG) 922 call timab(660,1,tsec) ! bse(exc_diago_driver) 923 call exc_diago_driver(Wfd,Bsp,BS_files,KS_BSt,QP_BSt,Cryst,Kmesh,Psps,& 924 & Pawtab,Hur,Hdr_bse,drude_plsmf,Epren) 925 call timab(660,2,tsec) ! bse(exc_diago_driver) 926 927 if (.FALSE.) then ! Calculate electron-hole excited state density. Not tested at all. 928 call exc_den(BSp,BS_files,ngfftf,nfftf,Kmesh,ktabr,Wfd) 929 end if 930 931 if (.FALSE.) then 932 paw_add_onsite=.FALSE.; spin_opt=1; which_fixed=1; eh_rcoord=(/zero,zero,zero/); nrcell=(/2,2,2/) 933 ! call exc_plot(Bsp,Bs_files,Wfd,Kmesh,Cryst,Psps,Pawtab,Pawrad,paw_add_onsite,spin_opt,which_fixed,eh_rcoord,nrcell,ngfftf) 934 end if 935 936 if (BSp%use_interp) then 937 MSG_ERROR("Interpolation technique not coded for diagonalization and CG") 938 end if 939 940 case (BSE_ALGO_Haydock) 941 call timab(661,1,tsec) ! bse(exc_haydock_driver) 942 943 if (BSp%use_interp) then 944 945 call exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_BSt,Wfd,Psps,Pawtab,Hur,Epren,& 946 & Kmesh_dense,KS_BSt_dense,QP_BSt_dense,Wfd_dense,Vcp_dense,grid) 947 else 948 call exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_BSt,Wfd,Psps,Pawtab,Hur,Epren) 949 end if 950 951 call timab(661,2,tsec) ! bse(exc_haydock_driver) 952 case default 953 write(msg,'(a,i0)')" Wrong BSE algorithm: ",BSp%algorithm 954 MSG_ERROR(msg) 955 end select 956 957 call timab(657,2,tsec) ! bse(mkexceps) 958 959 !===================== 960 !==== Free memory ==== 961 !===================== 962 ABI_FREE(ktabr) 963 ABI_FREE(ks_vhartr) 964 ABI_FREE(ks_vtrial) 965 ABI_FREE(ks_vxc) 966 ABI_FREE(ks_nhat) 967 ABI_FREE(ks_nhatgr) 968 ABI_FREE(ks_rhog) 969 ABI_FREE(ks_rhor) 970 ABI_FREE(qp_rhog) 971 ABI_FREE(qp_rhor) 972 ABI_FREE(qp_aerhor) 973 ! 974 !* Destroy local data structures. 975 call destroy_mpi_enreg(MPI_enreg_seq) 976 call crystal_free(Cryst) 977 call gsph_free(Gsph_x) 978 call gsph_free(Gsph_c) 979 call kmesh_free(Kmesh) 980 call kmesh_free(Qmesh) 981 call hdr_free(Hdr_wfk) 982 call hdr_free(Hdr_bse) 983 call ebands_free(KS_BSt) 984 call ebands_free(QP_BSt) 985 call vcoul_free(Vcp) 986 call pawhur_free(Hur) 987 ABI_DT_FREE(Hur) 988 call bs_parameters_free(BSp) 989 call wfd_free(Wfd) 990 call pawfgr_destroy(Pawfgr) 991 call eprenorms_free(Epren) 992 993 ! Free memory used for interpolation. 994 if (BSp%use_interp) then 995 call double_grid_free(grid) 996 call wfd_free(Wfd_dense) 997 call gsph_free(Gsph_x_dense) 998 call gsph_free(Gsph_c_dense) 999 call kmesh_free(Kmesh_dense) 1000 call kmesh_free(Qmesh_dense) 1001 call ebands_free(KS_BSt_dense) 1002 call ebands_free(QP_BSt_dense) 1003 call vcoul_free(Vcp_dense) 1004 call hdr_free(Hdr_wfk_dense) 1005 end if 1006 1007 ! Optional deallocation for PAW. 1008 if (Dtset%usepaw==1) then 1009 call pawrhoij_free(KS_Pawrhoij) 1010 ABI_DT_FREE(KS_Pawrhoij) 1011 call pawfgrtab_free(Pawfgrtab) 1012 ABI_DT_FREE(Pawfgrtab) 1013 call paw_ij_free(KS_paw_ij) 1014 ABI_DT_FREE(KS_paw_ij) 1015 call paw_an_free(KS_paw_an) 1016 ABI_DT_FREE(KS_paw_an) 1017 call pawpwff_free(Paw_pwff) 1018 end if 1019 ABI_DT_FREE(Paw_pwff) 1020 1021 call timab(650,2,tsec) ! bse(Total) 1022 1023 DBG_EXIT('COLL') 1024 1025 end subroutine bethe_salpeter