Branch data 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 : : !!***
|
