TABLE OF CONTENTS
ABINIT/setup_bse_interp [ Functions ]
NAME
setup_bse_interp
FUNCTION
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (Y. Gillet) 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
ngfft_gw(18)=Information about 3D FFT for density and potentials, see ~abinit/doc/variables/vargs.htm#ngfft acell(3)=Length scales of primitive translations (bohr) rprim(3,3)=Dimensionless real space primitive translations. Dtset<dataset_type>=All input variables for this dataset. Some of them might be redefined here TODO Dtfil=filenames and unit numbers used in abinit. fnameabi_wfkfile is changed is Fortran file is not found but a netcdf version with similar name is available.
OUTPUT
Cryst<crystal_structure>=Info on the crystalline Structure. Kmesh<BZ_mesh_type>=Structure defining the k-sampling for the wavefunctions. Qmesh<BZ_mesh_type>=Structure defining the q-sampling for the symmetrized inverse dielectric matrix. Gsph_x<gsphere_t=Data type gathering info on the G-sphere for wave functions and e^{-1}, KS_BSt<Bandstructure_type>=The KS band structure (energies, occupancies, k-weights...) Vcp<vcoul_t>=Structure gathering information on the Coulomb interaction in reciprocal space, including a possible cutoff in real space. ngfft_osc(18)=Contain all needed information about the 3D FFT for the oscillator matrix elements. See ~abinit/doc/variables/vargs.htm#ngfft Bsp<excparam>=Basic parameters defining the Bethe-Salpeter run. Completely initialed in output. Hdr_wfk<Hdr_type>=The header of the WFK file. Hdr_bse<Hdr_type>=Local header initialized from the parameters used for the Bethe-Salpeter calculation. w_file=File name used to construct W. Set to ABI_NOFILE if no external file is used.
PARENTS
bethe_salpeter
CHILDREN
apply_scissor,double_grid_init,ebands_copy,ebands_init,ebands_print ebands_report_gap,ebands_update_occ,find_qmesh,gsph_extend,gsph_init init_transitions,kmesh_init,kmesh_print,make_mesh,print_gsphere vcoul_init,wfk_read_eigenvalues,wrtout
SOURCE
50 #if defined HAVE_CONFIG_H 51 #include "config.h" 52 #endif 53 54 #include "abi_common.h" 55 56 subroutine setup_bse_interp(Dtset,Dtfil,BSp,Cryst,Kmesh,& 57 & Kmesh_dense,Qmesh_dense,KS_BSt_dense,QP_bst_dense,Gsph_x,Gsph_c,Vcp_dense,Hdr_wfk_dense,ngfftf,grid,comm) 58 59 use defs_basis 60 use defs_datatypes 61 use defs_abitypes 62 use m_bs_defs 63 use m_profiling_abi 64 use m_errors 65 use m_xmpi 66 use m_nctk 67 use m_hdr 68 69 use m_gwdefs, only : GW_Q0_DEFAULT 70 use m_io_tools, only : file_exists 71 use m_crystal, only : crystal_t 72 use m_bz_mesh, only : kmesh_t, kmesh_init, kmesh_print, find_qmesh, make_mesh 73 use m_double_grid, only : double_grid_t, double_grid_init 74 use m_ebands, only : ebands_init, ebands_print, ebands_copy, ebands_free, ebands_update_occ, & 75 apply_scissor, ebands_report_gap 76 use m_vcoul, only : vcoul_t, vcoul_init 77 use m_gsphere, only : gsphere_t, gsph_init, print_gsphere, gsph_extend 78 use m_wfk, only : wfk_read_eigenvalues 79 80 !This section has been created automatically by the script Abilint (TD). 81 !Do not modify the following lines by hand. 82 #undef ABI_FUNC 83 #define ABI_FUNC 'setup_bse_interp' 84 use interfaces_14_hidewrite 85 !End of the abilint section 86 87 implicit none 88 89 !Arguments ------------------------------------ 90 !scalars 91 integer,intent(in) :: comm 92 type(dataset_type),intent(in) :: Dtset 93 type(datafiles_type),intent(inout) :: Dtfil 94 type(excparam),intent(inout) :: Bsp 95 type(hdr_type),intent(out) :: Hdr_wfk_dense 96 type(crystal_t),intent(in) :: Cryst 97 type(kmesh_t),intent(in) :: Kmesh 98 type(kmesh_t),intent(out) :: Kmesh_dense,Qmesh_dense 99 type(ebands_t),intent(out) :: KS_BSt_dense,QP_Bst_dense 100 type(double_grid_t),intent(out) :: grid 101 type(vcoul_t),intent(out) :: Vcp_dense 102 type(gsphere_t),intent(out) :: Gsph_x,Gsph_c 103 !arrays 104 integer,intent(in) :: ngfftf(18) 105 106 !Local variables ------------------------------ 107 !scalars 108 integer,parameter :: pertcase0=0,master=0 109 integer :: bantot_dense,ib,ibtot,ik_ibz,isppol,jj 110 integer :: nbnds_kss_dense 111 integer :: spin,hexc_size 112 integer :: my_rank 113 integer :: it 114 integer :: nprocs 115 integer :: is1,is2,is3,is4 116 real(dp) :: nelect_hdr_dense 117 logical,parameter :: remove_inv=.FALSE. 118 character(len=500) :: msg 119 character(len=fnlen) :: wfk_fname_dense 120 integer :: nqlwl 121 !arrays 122 integer,allocatable :: npwarr(:) 123 real(dp),allocatable :: shiftk(:,:) 124 real(dp),allocatable :: doccde(:),eigen(:),occfact(:) 125 real(dp),pointer :: energies_p_dense(:,:,:) 126 complex(dpc),allocatable :: gw_energy(:,:,:) 127 integer,allocatable :: nbands_temp(:) 128 integer :: kptrlatt_dense(3,3) 129 real(dp),allocatable :: qlwl(:,:) 130 real(dp) :: minmax_tene(2) 131 132 !************************************************************************ 133 134 DBG_ENTER("COLL") 135 136 kptrlatt_dense = zero 137 138 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 139 140 SELECT CASE(BSp%interp_mode) 141 CASE (1,2,3,4) 142 nbnds_kss_dense = -1 143 wfk_fname_dense = Dtfil%fnameabi_wfkfine 144 call wrtout(std_out,"BSE Interpolation: will read energies from: "//trim(wfk_fname_dense),"COLL") 145 146 if (nctk_try_fort_or_ncfile(wfk_fname_dense, msg) /= 0) then 147 MSG_ERROR(msg) 148 end if 149 150 Dtfil%fnameabi_wfkfine = wfk_fname_dense 151 152 call wfk_read_eigenvalues(wfk_fname_dense,energies_p_dense,Hdr_wfk_dense,comm) 153 nbnds_kss_dense = MAXVAL(Hdr_wfk_dense%nband) 154 CASE DEFAULT 155 MSG_ERROR("Not yet implemented") 156 END SELECT 157 158 nelect_hdr_dense = Hdr_wfk_dense%nelect 159 160 if (ABS(Dtset%nelect-nelect_hdr_dense)>tol6) then 161 write(msg,'(2(a,f8.2))')& 162 & "File contains ", nelect_hdr_dense," electrons but nelect initialized from input is ",Dtset%nelect 163 MSG_ERROR(msg) 164 end if 165 166 ! Setup of the k-point list and symmetry tables in the BZ 167 SELECT CASE(BSp%interp_mode) 168 CASE (1,2,3,4) 169 if(Dtset%chksymbreak == 0) then 170 ABI_MALLOC(shiftk,(3,Dtset%nshiftk)) 171 kptrlatt_dense(:,1) = BSp%interp_kmult(1)*Dtset%kptrlatt(:,1) 172 kptrlatt_dense(:,2) = BSp%interp_kmult(2)*Dtset%kptrlatt(:,2) 173 kptrlatt_dense(:,3) = BSp%interp_kmult(3)*Dtset%kptrlatt(:,3) 174 do jj = 1,Dtset%nshiftk 175 shiftk(:,jj) = Bsp%interp_kmult(:)*Dtset%shiftk(:,jj) 176 end do 177 call make_mesh(Kmesh_dense,Cryst,Dtset%kptopt,kptrlatt_dense,Dtset%nshiftk,shiftk,break_symmetry=.TRUE.) 178 ABI_FREE(shiftk) 179 else 180 !Initialize Kmesh with no wrapping inside ]-0.5;0.5] 181 call kmesh_init(Kmesh_dense,Cryst,Hdr_wfk_dense%nkpt,Hdr_wfk_dense%kptns,Dtset%kptopt) 182 end if 183 CASE DEFAULT 184 MSG_ERROR("Not yet implemented") 185 END SELECT 186 187 ! Init Qmesh 188 call find_qmesh(Qmesh_dense,Cryst,Kmesh_dense) 189 190 call gsph_init(Gsph_c,Cryst,0,ecut=Dtset%ecuteps) 191 192 call double_grid_init(Kmesh,Kmesh_dense,Dtset%kptrlatt,BSp%interp_kmult,grid) 193 194 BSp%nkibz_interp = Kmesh_dense%nibz !We might allow for a smaller number of points.... 195 196 call kmesh_print(Kmesh_dense,"Interpolated K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 197 call kmesh_print(Kmesh_dense,"Interpolated K-mesh for the wavefunctions",ab_out, 0, "COLL") 198 199 if (nbnds_kss_dense < Dtset%nband(1)) then 200 write(msg,'(2(a,i0),3a,i0)')& 201 & 'Interpolated WFK file contains only ', nbnds_kss_dense,' levels instead of ',Dtset%nband(1),' required;',ch10,& 202 & 'The calculation will be done with nbands= ',nbnds_kss_dense 203 MSG_WARNING(msg) 204 MSG_ERROR("Not supported yet !") 205 end if 206 207 ABI_MALLOC(nbands_temp,(Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol)) 208 do isppol=1,Hdr_wfk_dense%nsppol 209 do ik_ibz=1,Hdr_wfk_dense%nkpt 210 nbands_temp(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) = Dtset%nband(1) 211 end do 212 end do 213 214 call gsph_extend(Gsph_c,Cryst,Dtset%ecutwfn,Gsph_x) 215 call print_gsphere(Gsph_x,unit=std_out,prtvol=Dtset%prtvol) 216 217 nqlwl=1 218 ABI_MALLOC(qlwl,(3,nqlwl)) 219 qlwl(:,nqlwl)= GW_Q0_DEFAULT 220 221 ! Compute Coulomb term on the largest G-sphere. 222 if (Gsph_x%ng > Gsph_c%ng ) then 223 call vcoul_init(Vcp_dense,Gsph_x,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%icutcoul,& 224 & Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,ngfftf,comm) 225 else 226 call vcoul_init(Vcp_dense,Gsph_c,Cryst,Qmesh_dense,Kmesh_dense,Dtset%rcut,Dtset%icutcoul,& 227 & Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,ngfftf,comm) 228 end if 229 230 ABI_FREE(qlwl) 231 232 bantot_dense=SUM(Hdr_wfk_dense%nband(1:Hdr_wfk_dense%nkpt*Hdr_wfk_dense%nsppol)) 233 ABI_ALLOCATE(doccde,(bantot_dense)) 234 ABI_ALLOCATE(eigen,(bantot_dense)) 235 ABI_ALLOCATE(occfact,(bantot_dense)) 236 doccde=zero; eigen=zero; occfact=zero 237 238 jj=0; ibtot=0 239 do isppol=1,Hdr_wfk_dense%nsppol 240 do ik_ibz=1,Hdr_wfk_dense%nkpt 241 do ib=1,Hdr_wfk_dense%nband(ik_ibz+(isppol-1)*Hdr_wfk_dense%nkpt) 242 ibtot=ibtot+1 243 if (ib<=BSP%nbnds) then 244 jj=jj+1 245 occfact(jj)=Hdr_wfk_dense%occ(ibtot) 246 eigen (jj)=energies_p_dense(ib,ik_ibz,isppol) 247 end if 248 end do 249 end do 250 end do 251 252 ABI_FREE(energies_p_dense) 253 254 ABI_MALLOC(npwarr,(kmesh_dense%nibz)) 255 npwarr=BSP%npwwfn 256 257 call ebands_init(bantot_dense,KS_BSt_dense,Dtset%nelect,doccde,eigen,Hdr_wfk_dense%istwfk,Kmesh_dense%ibz,nbands_temp,& 258 & Kmesh_dense%nibz,npwarr,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,Hdr_wfk_dense%tphysel,Hdr_wfk_dense%tsmear,& 259 & Hdr_wfk_dense%occopt,occfact,Kmesh_dense%wt,& 260 & hdr_wfk_dense%charge, hdr_wfk_dense%kptopt, hdr_wfk_dense%kptrlatt_orig, hdr_wfk_dense%nshiftk_orig, & 261 & hdr_wfk_dense%shiftk_orig, hdr_wfk_dense%kptrlatt, hdr_wfk_dense%nshiftk, hdr_wfk_dense%shiftk) 262 263 ABI_DEALLOCATE(doccde) 264 ABI_DEALLOCATE(eigen) 265 ABI_DEALLOCATE(npwarr) 266 267 ABI_FREE(nbands_temp) 268 269 ABI_FREE(occfact) 270 271 !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when 272 ! the occupation scheme for semiconductors is used. 273 call ebands_update_occ(KS_BSt_dense,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 274 275 call ebands_print(KS_BSt_dense,"Interpolated band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol) 276 277 call ebands_report_gap(KS_BSt_dense,header="Interpolated KS band structure",unit=std_out,mode_paral="COLL") 278 279 BSp%nkbz_interp = Kmesh_dense%nbz 280 281 call ebands_copy(KS_BSt_dense,QP_bst_dense) 282 283 SELECT CASE (Bsp%calc_type) 284 CASE (BSE_HTYPE_RPA_KS) 285 if (ABS(BSp%mbpt_sciss)>tol6) then 286 write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies." 287 call wrtout(std_out,msg,"COLL") 288 call apply_scissor(QP_BSt_dense,BSp%mbpt_sciss) 289 else 290 write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]." 291 call wrtout(std_out,msg,"COLL") 292 end if 293 ! 294 CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw 295 MSG_ERROR("Not yet implemented with interpolation !") 296 CASE (BSE_HTYPE_RPA_QP) 297 MSG_ERROR("Not implemented error!") 298 CASE DEFAULT 299 write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type 300 MSG_ERROR(msg) 301 END SELECT 302 303 call ebands_report_gap(QP_BSt_dense,header=" Interpolated QP band structure",unit=std_out,mode_paral="COLL") 304 305 ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index. 306 ! FIXME: linewidths not coded. 307 ABI_ALLOCATE(gw_energy,(BSp%nbnds,Kmesh_dense%nibz,Dtset%nsppol)) 308 gw_energy = QP_BSt_dense%eig 309 310 ABI_ALLOCATE(Bsp%nreh_interp,(Hdr_wfk_dense%nsppol)) 311 Bsp%nreh_interp=zero 312 313 call init_transitions(BSp%Trans_interp,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz_interp,Bsp%nbnds,& 314 & Bsp%nkibz_interp,Hdr_wfk_dense%nsppol,Hdr_wfk_dense%nspinor,gw_energy,QP_BSt_dense%occ,Kmesh_dense%tab,minmax_tene,& 315 & Bsp%nreh_interp) 316 317 ABI_DEALLOCATE(gw_energy) 318 319 do spin=1,Dtset%nsppol 320 write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh_interp(spin) 321 call wrtout(std_out,msg,"COLL") 322 end do 323 324 if (ANY(Bsp%nreh_interp/=Bsp%nreh_interp(1))) then 325 write(msg,'(a,(i0))')" BSE code does not support different number of transitions for the two spin channels",Bsp%nreh 326 MSG_ERROR(msg) 327 end if 328 ! 329 ! Create transition table vcks2t 330 is1=BSp%lomo_min;is2=BSp%homo_max;is3=BSp%lumo_min;is4=BSp%humo_max 331 ABI_ALLOCATE(Bsp%vcks2t_interp,(is1:is2,is3:is4,BSp%nkbz_interp,Dtset%nsppol)) 332 Bsp%vcks2t_interp = 0 333 334 do spin=1,Dtset%nsppol 335 do it=1,BSp%nreh_interp(spin) 336 BSp%vcks2t_interp(BSp%Trans_interp(it,spin)%v,BSp%Trans_interp(it,spin)%c,& 337 & BSp%Trans_interp(it,spin)%k,spin) = it 338 end do 339 end do 340 341 hexc_size = SUM(Bsp%nreh_interp); if (Bsp%use_coupling>0) hexc_size=2*hexc_size 342 if (Bsp%nstates_interp<=0) then 343 Bsp%nstates_interp=hexc_size 344 else 345 if (Bsp%nstates_interp>hexc_size) then 346 Bsp%nstates_interp=hexc_size 347 write(msg,'(2(a,i0),2a)')& 348 & "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates_interp,ch10,& 349 & "the number of excitonic states nstates has been modified" 350 MSG_WARNING(msg) 351 end if 352 end if 353 354 DBG_EXIT("COLL") 355 356 end subroutine setup_bse_interp