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