Line data Source code
1 : !!****m* ABINIT/m_bs_defs
2 : !! NAME
3 : !! m_bs_defs
4 : !!
5 : !! FUNCTION
6 : !! This module defines basic structures used for Bethe-Salpeter calculations.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1992-2026 ABINIT and EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_bs_defs
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 :
29 : implicit none
30 :
31 : private
32 : !!***
33 :
34 : ! Algorithm used to solve the BS problem.
35 : integer,public,parameter :: BSE_ALGO_NONE =0 ! Build the BSE Hamiltonian but skip the solution of the BSE equation.
36 : integer,public,parameter :: BSE_ALGO_DDIAGO =1 ! Direct diagonalization.
37 : integer,public,parameter :: BSE_ALGO_HAYDOCK =2 ! Haydock recursion method.
38 : integer,public,parameter :: BSE_ALGO_CG =3 ! Iterative diagonalization via CG method.
39 :
40 : ! Flags defining the content of the files used to restart the Haydock method.
41 : integer,public,parameter :: BSE_HAYD_IMEPS =1
42 : integer,public,parameter :: BSE_HAYD_DOS =2
43 :
44 : ! Approximations for the excitonic Hamiltonian.
45 : integer,public,parameter :: BSE_HTYPE_RPA_KS =1 ! Use KS results to construct the RPA polarizability.
46 : integer,public,parameter :: BSE_HTYPE_RPA_QPENE =2 ! Use KS orbitals and QP energies to construct the RPA polarizability.
47 : integer,public,parameter :: BSE_HTYPE_RPA_QP =3 ! Use QP orbitals and energies to construct the RPA polarizability.
48 :
49 : ! Flags for the treatment of W(G1,G2)
50 : integer,public,parameter :: BSE_WTYPE_NONE =0 ! Coulomb term not included
51 : integer,public,parameter :: BSE_WTYPE_FROM_SCR =1 ! W is read from a SCR file
52 : integer,public,parameter :: BSE_WTYPE_FROM_MDL =2 ! W is approximated with a model dielectric function.
53 :
54 : !$! Treatment of W(omega)
55 : !! integer,public,parameter :: BSE_WFREQ_STATIC=1
56 : !! integer,public,parameter :: BSE_WFREQ_PPM =2
57 : !! integer,public,parameter :: BSE_WFREQ_FULL =3
58 :
59 : ! Flags for the interpolation
60 : integer,public,parameter :: BSE_INTERP_YG =0 ! Interpolation with 8 neighbours
61 : integer,public,parameter :: BSE_INTERP_RL =1 ! Interpolation with 1 neighbour (Rohlfing & Louie 2000)
62 : integer,public,parameter :: BSE_INTERP_RL2 =2 ! Hybrid between RL & YG with 2 neighbours (for debug)
63 :
64 : character(len=fnlen),public,parameter :: BSE_NOFILE="None"
65 :
66 : !----------------------------------------------------------------------
67 :
68 : !!****t* m_bs_defs/transition
69 : !! NAME
70 : !! transition
71 : !!
72 : !! FUNCTION
73 : !! The transition derived data type is used to store the correspondence
74 : !! between the transition index and the set of quantum numbers (ik_bz,v,c)
75 : !! The energy of the transition is stored as well.
76 : !!
77 : !! SOURCE
78 :
79 : type,public :: transition
80 : integer :: k = 0 ! Index of the k-point in the BZ
81 : integer :: v = 0 ! Valence band index.
82 : integer :: c = 0 ! Conduction band index.
83 : complex(dp) :: en=huge(one) ! Transition energy
84 : end type transition
85 :
86 : public :: init_transitions ! Main creation method.
87 : public :: repr_trans ! Returns a string representing the transition or a couple of transitions.
88 : !!***
89 :
90 : interface repr_trans
91 : module procedure repr_1trans
92 : module procedure repr_2trans
93 : end interface repr_trans
94 :
95 : !----------------------------------------------------------------------
96 :
97 : !!****t* m_bs_defs/excparam
98 : !! NAME
99 : !! excparam
100 : !!
101 : !! FUNCTION
102 : !! The excparam derived data type contains the parameters controlling the BS calculation.
103 : !!
104 : !! SOURCE
105 :
106 : type,public :: excparam
107 :
108 : !scalars
109 : integer :: algorithm ! Algorithm used for computing the BS dielectric function.
110 : integer :: calc_type ! Calculation type (see Dtset%bs_calc_type).
111 : integer :: hayd_term ! Option for the terminator used in the Haydock solver.
112 : integer :: use_coupling ! Include off-diagonal block coupling resonant and anti-resonant transitions.
113 : integer :: exchange_term ! Include the exchange term in the BS Hamiltonian.
114 : integer :: inclvkb ! Option for the inclusion of the commutator [Vnl, r] for the optical limit.
115 : integer :: mdlf_type ! Model dielectric function type.
116 : integer :: nline ! Number of line minimization used for CG minimization.
117 : integer :: nbdbuf ! Number of states in the buffer that will be excluded from the convergence check (CG only)
118 : integer :: nstates ! Number of states that will be considered in the CG minimization.
119 : integer :: npweps ! No. of G in the Screening.
120 : integer :: npwwfn ! No. of G for wave functions.
121 : !$integer :: npwx ! No. of G for the exchange part.
122 : integer :: npwvec ! MAX between npwwfn and npweps
123 : integer :: nbnds ! Total number of bands considered.
124 :
125 : !integer :: nbndv ! No. of valence states treated (homo-lomo+1)
126 : !integer :: nbndc ! No. of conduction states (humo-lumo+1)
127 : !integer :: lomo,homo ! Lowest and highest occupied orbital considered.
128 : !integer :: lumo,humo ! Lowest and highest unoccupied orbital considered.
129 :
130 : ! new for spin
131 : ! DO I need these?
132 : integer :: lomo_min,homo_max ! Lowest and highest occupied orbital considered.
133 : integer :: lumo_min,humo_max ! Lowest and highest unoccupied orbital considered.
134 : integer :: maxnbndv, maxnbndc
135 :
136 : integer,allocatable :: lomo_spin(:) ! Lowest occupied orbital considered for the different spins.
137 : integer,allocatable :: homo_spin(:) ! Highest occupied orbital considered for the different spins.
138 : integer,allocatable :: lumo_spin(:) ! Lowest unoccupied orbital considered for the different spins
139 : integer,allocatable :: humo_spin(:) ! Highest unoccupied orbital considered for the different spins
140 : integer,allocatable :: nbndv_spin(:) ! No. of valence states treated (homo-lomo+1)
141 : integer,allocatable :: nbndc_spin(:) ! No. of conduction states (humo-lumo+1)
142 : ! end new
143 :
144 : integer :: niter ! No. of iterations for (Haydock|CG).
145 : integer :: nkibz, nkbz ! No. of k-points in the IBZ and BZ (resp.)
146 : integer :: nomega ! No. of frequencies for epsilon.
147 : integer :: nq ! Number of "small" q for optical limit.
148 : integer :: nsppol ! Number of independent spin polarizations.
149 : integer :: wtype ! Option used for dealing with W (see BSE_WTYPE_) flags
150 :
151 : !Interp@BSE
152 : integer :: interp_mode ! Mode of interpolation
153 : integer :: interp_method ! Method of interpolation
154 : integer :: nkibz_interp,nkbz_interp ! Number of points in the interpolated kmesh
155 : integer :: nstates_interp ! Number of states of interpolation
156 : integer :: rl_nb ! Index of the nb in Rohlfing and Louie technique
157 :
158 : real(dp) :: ecutwfn ! Cutoff energy for wavefunctions.
159 : real(dp) :: ecuteps ! Cutoff energy for W.
160 : real(dp) :: eps_inf ! Electronic dielectric constant used for the model dielectric function.
161 : real(dp) :: mbpt_sciss ! Scissors energy (used if it absolute value is > tol6)
162 : real(dp) :: omegai ! First omega for epsilon.
163 : real(dp) :: omegae ! Last omega for epsilon (defaults to 10eV)
164 : real(dp) :: domega ! Step of the frequency mesh.
165 : real(dp) :: broad ! Lorentzian Broadening.
166 : real(dp) :: ircut ! Infrared cutoff for transitions
167 : real(dp) :: uvcut ! Ultraviolet cutoff for transitions.
168 : real(dp) :: haydock_tol(2) ! Tolerance for stopping the Haydock algorithm.
169 : real(dp) :: cg_tolwfr ! Tolerance for stopping the CG algorithm
170 :
171 : !Interp@BSE
172 : real(dp) :: interp_m3_width ! Width of the interpolated M3 method along the diagonal
173 :
174 : logical :: use_diagonal_Wgg ! Use diagonal approximation for Wgg.
175 : logical :: use_coulomb_term ! Include W term in the BS Hamiltonian.
176 : logical :: have_complex_ene ! .TRUE. if energies have a non-zero imaginary part.
177 :
178 : !Interp@BSE
179 : logical :: use_interp ! .TRUE. if we use interpolation technique
180 : logical :: prep_interp ! .TRUE. if we prepare interpolation with ABC
181 : logical :: sum_overlaps ! .TRUE. if making the sum of the overlaps to 1
182 : logical :: prt_ncham ! .TRUE. if we dump the hamiltonian in NetCDF
183 :
184 : logical :: do_ep_renorm ! .TRUE. for electron-phonon renormalization of the spectrum
185 : logical :: do_lifetime ! .TRUE. if using elphon lifetime (not yet implemented)
186 :
187 : !arrays
188 : integer :: mg0(3) ! For each reduced direction gives the max G0 component
189 : ! to account for umklapp processes
190 :
191 : integer,allocatable :: nreh(:)
192 : ! nreh(nsppol)
193 : ! Number of resonant electron-hole transitions for each spin.
194 :
195 : integer,allocatable :: vcks2t(:,:,:,:)
196 : ! vcks2t(v,c,ik_bz,spin) gives the transition index associated to (v,c,kbz,spin)
197 :
198 : !Interp@BSE
199 : integer :: interp_kmult(3) ! Factor to subdivide kmesh sampling
200 :
201 : integer,allocatable :: nreh_interp(:)
202 : ! nreh_interp(nsppol)
203 : ! Number of transitions for the interpolated mesh
204 :
205 : integer,allocatable :: vcks2t_interp(:,:,:,:)
206 : ! vcks2t(v,c,ik_bz_dense,spin) : Transition index for the dense kmesh
207 :
208 : real(dp),allocatable :: q(:,:) ! Q-points for optical limit (reduced coordinates).
209 :
210 : complex(dp),allocatable :: omega(:)
211 : ! omega(nomega)
212 : ! Frequency mesh for epsilon (including the complex imaginary shift)
213 :
214 : type(transition),allocatable :: Trans(:,:)
215 : ! Trans(max_nreh,nsppol)
216 :
217 : type(transition),allocatable :: Trans_interp(:,:)
218 : ! Transitions for interpolated mesh
219 :
220 : contains
221 : procedure :: free => bs_parameters_free
222 : procedure :: print => print_bs_parameters
223 : procedure :: calctype2str => bsp_calctype2str
224 : end type excparam
225 : !!***
226 :
227 : !!****t* m_bs_defs/excfiles
228 : !! NAME
229 : !! excfiles
230 : !!
231 : !! FUNCTION
232 : !! The excfiles derived data type contains file names and unit numbers used to store
233 : !! temporary or final results of the Bethe-Salpeter calculation.
234 : !!
235 : !! SOURCE
236 :
237 : type,public :: excfiles
238 :
239 : character(len=fnlen) :: in_hreso = BSE_NOFILE
240 : ! Name of the input file with the resonant part of the Hamiltonian (Hermitian).
241 :
242 : character(len=fnlen) :: out_hreso = BSE_NOFILE
243 : ! Name of the output file with the resonant part of the Hamiltonian (Hermitian).
244 :
245 : character(len=fnlen) :: in_hcoup = BSE_NOFILE
246 : ! Name of the input file with the coupling part of the Hamiltonian (Symmetric).
247 :
248 : character(len=fnlen) :: out_hcoup = BSE_NOFILE
249 : ! Name of the output file with the coupling part of the Hamiltonian (Symmetric).
250 :
251 : character(len=fnlen) :: in_eig = BSE_NOFILE
252 : ! Name of the input file with the eigenvalues and the eigenvectors of the Hamiltonian.
253 :
254 : character(len=fnlen) :: out_eig = BSE_NOFILE
255 : ! Name of the output file with the eigenvalues and the eigenvectors of the Hamiltonian.
256 :
257 : character(len=fnlen) :: in_haydock_basename = BSE_NOFILE
258 : ! Name of the input file used to restart Haydock algorithm.
259 :
260 : character(len=fnlen) :: out_basename = BSE_NOFILE
261 : ! Prefix to be used for other output files.
262 :
263 : contains
264 : procedure :: print => print_bs_files
265 : ! Printout of the excfiles data type.
266 : end type excfiles
267 : !!***
268 :
269 : contains
270 : !!***
271 :
272 : !----------------------------------------------------------------------
273 :
274 : !!****f* m_bs_defs/bs_parameters_free
275 : !! NAME
276 : !! bs_parameters_free
277 : !!
278 : !! FUNCTION
279 : !! Free all memory allocated in a structure of type excparam
280 : !!
281 : !! SOURCE
282 :
283 29 : subroutine bs_parameters_free(BSp)
284 :
285 : !Arguments ------------------------------------
286 : class(excparam),intent(inout) :: BSp
287 : !************************************************************************
288 :
289 29 : ABI_SFREE(BSp%q)
290 29 : ABI_SFREE(Bsp%nreh)
291 29 : ABI_SFREE(Bsp%vcks2t)
292 29 : ABI_SFREE(Bsp%omega)
293 29 : ABI_SFREE(Bsp%lomo_spin)
294 29 : ABI_SFREE(Bsp%homo_spin)
295 29 : ABI_SFREE(Bsp%lumo_spin)
296 29 : ABI_SFREE(Bsp%humo_spin)
297 29 : ABI_SFREE(Bsp%nbndv_spin)
298 29 : ABI_SFREE(Bsp%nbndc_spin)
299 29 : ABI_SFREE(Bsp%Trans)
300 29 : ABI_SFREE(Bsp%nreh_interp)
301 29 : ABI_SFREE(Bsp%vcks2t_interp)
302 29 : ABI_SFREE(Bsp%Trans_interp)
303 :
304 29 : end subroutine bs_parameters_free
305 : !!***
306 :
307 : !----------------------------------------------------------------------
308 :
309 : !!****f* m_bs_defs/print_bs_parameters
310 : !! NAME
311 : !! print_bs_parameters
312 : !!
313 : !! FUNCTION
314 : !! Printout of the parameters used for the BS calculation.
315 : !!
316 : !! SOURCE
317 :
318 58 : subroutine print_bs_parameters(BSp, header, unit, mode_paral, prtvol)
319 :
320 : !Arguments ------------------------------------
321 : !scalars
322 : class(excparam),intent(in) :: BSp
323 : integer,optional,intent(in) :: unit,prtvol
324 : character(len=4),optional,intent(in) :: mode_paral
325 : character(len=*),optional,intent(in) :: header
326 :
327 : !Local variables ------------------------------
328 : !scalars
329 : integer :: my_unt,my_prtvol,iq,ii,spin
330 : character(len=4) :: my_mode
331 : character(len=500) :: msg
332 : ! *********************************************************************
333 :
334 58 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
335 58 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
336 58 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
337 :
338 58 : msg=' ==== Parameters of the Bethe-Salpeter run ==== '
339 58 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
340 58 : call wrtout(my_unt,msg,my_mode)
341 :
342 58 : select case (Bsp%algorithm)
343 : case (BSE_ALGO_NONE)
344 0 : msg = " Algorithm: Build Hamiltonian but skip the calculation of the spectrum."
345 : case (BSE_ALGO_DDIAGO)
346 10 : msg = " Algorithm: Direct diagonalization."
347 : case (BSE_ALGO_HAYDOCK)
348 46 : msg = " Algorithm: Haydock technique."
349 : case (BSE_ALGO_CG)
350 2 : msg = " Algorithm: Conjugate gradient."
351 : case default
352 58 : msg = " Algorithm: Unknown!."
353 : end select
354 58 : call wrtout(my_unt,msg,my_mode)
355 :
356 : write(msg,'(4(a,i0,a))')&
357 58 : ' Dimension of the v, W matrices, npweps = ',BSp%npweps,ch10,&
358 58 : ' Cutoff for the wavefunctions, npwwfn = ',BSp%npwwfn,ch10,&
359 58 : ' Number of k-points in the IBZ, nkibz = ',BSp%nkibz,ch10,&
360 116 : ' Highest empty band included, nband = ',BSp%nbnds,""
361 58 : call wrtout(my_unt,msg,my_mode)
362 :
363 118 : do spin=1,Bsp%nsppol
364 60 : msg = " === Spin UP ==="; if (spin == 2) msg = " === Spin DOWN ==="
365 60 : call wrtout(my_unt,msg,my_mode)
366 : write(msg,'(5(a,i0,a))')&
367 60 : ' Number of resonant transitions ',BSp%nreh(spin),ch10,&
368 60 : ' Lowest occupied state ',BSp%lomo_spin(spin),ch10,&
369 60 : ' Highest occupied state ',BSp%homo_spin(spin),ch10,&
370 60 : ' Lowest unoccupied state ',BSp%lumo_spin(spin),ch10,&
371 120 : ' Highest unoccupied state ',BSp%nbnds,""
372 : !' Number of valence bands ',BSp%nbndv,ch10,&
373 : !' Number of conduction bands ',BSp%nbndc,""
374 118 : call wrtout(my_unt,msg,my_mode)
375 : end do
376 :
377 : write(msg,'(3(a,f6.2,a),a,f6.2)')&
378 58 : ' Minimum frequency [eV] Emin = ',BSp%omegai*Ha_eV,ch10,&
379 58 : ' Maximum frequency [eV] Emax = ',BSp%omegae*Ha_eV,ch10,&
380 58 : ' Frequency step [eV] dE = ',BSp%domega*Ha_eV,ch10,&
381 116 : ' Lorentzian broadening [eV] eta = ',BSp%broad*Ha_eV
382 58 : call wrtout(my_unt,msg,my_mode)
383 :
384 : ! Calculation type
385 58 : call bsp_calctype2str(Bsp, msg)
386 58 : call wrtout(my_unt,msg,my_mode)
387 :
388 58 : if (ABS(Bsp%mbpt_sciss)>tol6) then
389 58 : write(msg,'(a,f5.2)')" Scissors operator energy [eV] = ",Bsp%mbpt_sciss*Ha_eV
390 58 : call wrtout(my_unt,msg,my_mode)
391 : end if
392 :
393 58 : msg=' Local fields effects (v term) excluded'
394 58 : if (BSp%exchange_term>0) msg=' Local fields effects (v term) included'
395 58 : call wrtout(my_unt,msg,my_mode)
396 :
397 58 : msg=' Excitonic effects (W term) excluded'
398 58 : if (BSp%use_coulomb_term) msg=' Excitonic effects (W term) included'
399 58 : call wrtout(my_unt,msg,my_mode)
400 :
401 58 : if (BSp%use_coulomb_term) then
402 58 : msg=" Full W_GG' included"
403 58 : if (BSp%use_diagonal_Wgg) msg=' Only diagonal term W_GG included'
404 58 : call wrtout(my_unt,msg,my_mode)
405 58 : if (BSp%wtype==BSE_WTYPE_FROM_SCR) then
406 22 : call wrtout(my_unt," W is read from an external SCR file",my_mode)
407 : end if
408 58 : if (BSp%wtype==BSE_WTYPE_FROM_MDL) then
409 36 : call wrtout(my_unt," W is approximated with the model dielectric function",my_mode)
410 : end if
411 : end if
412 :
413 58 : msg=' Resonant-only calculation (Hermitian case)'
414 58 : if (BSp%use_coupling>0) msg=' Resonant + Coupling calculation'
415 58 : call wrtout(my_unt,msg,my_mode)
416 :
417 58 : if(Bsp%use_interp) then
418 8 : call wrtout(my_unt,' Interpolation technique used',my_mode)
419 : end if
420 :
421 58 : if(Bsp%use_interp) then
422 12 : select case (Bsp%interp_mode)
423 : case (1)
424 4 : msg = ' Interpolation using WFK on the dense mesh'
425 : case (2)
426 2 : msg = ' Interpolation using WFK on the dense mesh + ABC divergence'
427 : case (3)
428 2 : msg = ' Interpolation using WFK on the dense mesh + ABC divergence along diagonal'
429 : case (4)
430 0 : msg = ' Interpolation using WFK on the dense mesh'
431 : case default
432 8 : msg = ' Unknown interpolation technique'
433 : end select
434 8 : call wrtout(my_unt,msg,my_mode)
435 :
436 8 : if(BSp%prep_interp) then
437 0 : call wrtout(my_unt,' Prepare interpolation technique with ABC',my_mode)
438 : end if
439 :
440 8 : if(BSp%interp_method == BSE_INTERP_YG) then
441 6 : write(msg,'(a)') " Use Y. Gillet interpolation with 8 neighbours"
442 6 : call wrtout(my_unt, msg, my_mode)
443 2 : else if(BSP%interp_method == BSE_INTERP_RL2) then
444 0 : write(msg,'(a)') " Use only 2 neighbours to interpolate linearly (Debug mode)"
445 0 : call wrtout(my_unt, msg, my_mode)
446 2 : else if(BSp%interp_method == BSE_INTERP_RL) then
447 2 : write(msg,'(a,i0)') " Use Rohlfing and Louie with nb = ",BSp%rl_nb
448 2 : call wrtout(my_unt, msg, my_mode)
449 : end if
450 :
451 8 : if(BSp%sum_overlaps) then
452 8 : call wrtout(my_unt, " Summing overlaps to 1 in the interpolation",my_mode)
453 : end if
454 : end if
455 :
456 58 : write(msg,'(a)')ch10
457 58 : call wrtout(my_unt,msg,my_mode)
458 :
459 58 : call wrtout(my_unt,' Calculating epsilon_Macro(q-->0,w), along the following directions:',my_mode)
460 406 : do iq=1,BSp%nq
461 1392 : write(msg,'(a,3f10.6,2a)')' q = (',(BSp%q(ii,iq),ii=1,3), ') [r.l.u.]'
462 406 : call wrtout(my_unt,msg,my_mode)
463 : end do
464 :
465 : !TODO
466 : !Add file sizes and size of the buffer used for the matrix.
467 :
468 58 : end subroutine print_bs_parameters
469 : !!***
470 :
471 : !----------------------------------------------------------------------
472 :
473 : !!****f* m_bs_defs/bsp_calctype2str
474 : !! NAME
475 : !! bsp_calctype2str
476 : !!
477 : !! FUNCTION
478 : !! Returns a string with the calculation type.
479 : !!
480 : !! SOURCE
481 :
482 187 : subroutine bsp_calctype2str(BSp, str)
483 :
484 : !Arguments ------------------------------------
485 : !scalars
486 : class(excparam),intent(in) :: BSp
487 : character(len=500),intent(out) :: str
488 :
489 : !************************************************************************
490 :
491 374 : SELECT CASE (Bsp%calc_type)
492 : CASE (BSE_HTYPE_RPA_KS)
493 187 : str = " RPA L0 with KS energies and KS wavefunctions"
494 : CASE (BSE_HTYPE_RPA_QPENE)
495 0 : str = " RPA L0 with QP energies and KS wavefunctions"
496 : CASE (BSE_HTYPE_RPA_QP)
497 0 : str = " RPA L0 with QP energies and QP wavefunctions"
498 : CASE DEFAULT
499 187 : str = " Unknown"
500 : END SELECT
501 :
502 187 : end subroutine bsp_calctype2str
503 : !!***
504 :
505 : !----------------------------------------------------------------------
506 :
507 : !!****f* m_bs_defs/init_transitions
508 : !! NAME
509 : !! init_transitions
510 : !!
511 : !! FUNCTION
512 : !! Main creation method for the transition structured datatype.
513 : !!
514 : !! INPUTS
515 : !! lomo_spin(nsppol)
516 : !! humo_spin(nsppol)
517 : !! ir_cut,uv_cut
518 : !! nkbz=Number of k-points in the BZ.
519 : !! nbnds=Maximum number of bands.
520 : !! nkibz=Number of k-points in the IBZ.
521 : !! nsppol=Number of spins.
522 : !! nspinor=Number of spinor components.
523 : !! gw_energy
524 : !! occ=Occupation factors.
525 : !! ktab
526 : !!
527 : !! OUTPUT
528 : !! max_tene=Maximum transition energy.
529 : !! nreh(nsppol)=Number of resonant transitions for each spin.
530 : !!
531 : !! SIDE EFFECTS
532 : !! Trans(:,:)
533 : !! input: allocatable array
534 : !! output: Trans(max_nreh,nsppol) stores the correspondence t -> (band,kbz,spin) and the transition energy.
535 : !!
536 : !! SOURCE
537 :
538 33 : subroutine init_transitions(Trans,lomo_spin,humo_spin,ir_cut,uv_cut,nkbz,nbnds,nkibz,nsppol,nspinor,gw_energy,occ,ktab,&
539 33 : minmax_tene,nreh)
540 :
541 : !Arguments ------------------------------------
542 : !scalars
543 : integer,intent(in) :: nkbz,nbnds,nkibz,nsppol,nspinor
544 : real(dp),intent(in) :: ir_cut,uv_cut
545 : real(dp),intent(out) :: minmax_tene(2)
546 : type(transition),allocatable,intent(out) :: Trans(:,:)
547 : !arrays
548 : integer,intent(in) :: lomo_spin(nsppol),humo_spin(nsppol)
549 : integer,intent(in) :: ktab(nkbz)
550 : integer,intent(out) :: nreh(nsppol)
551 : real(dp),intent(in) :: occ(nbnds,nkibz,nsppol)
552 : complex(dp),intent(in) :: gw_energy(nbnds,nkibz,nsppol)
553 :
554 : !Local variables ------------------------------
555 : !scalars
556 : integer :: spin,it,ik_bz,ik_ibz,iv,ic,max_occ,sweep,max_nreh,lomo,humo
557 : real(dp) :: tene, delta_f,min_tene,max_tene
558 : complex(dp) :: cplx_enet
559 : logical :: add_transition
560 : !************************************************************************
561 :
562 : ! Find transitions
563 33 : max_occ=2/(nsppol*nspinor)
564 67 : nreh=0
565 : min_tene = -one; max_tene = zero
566 : !
567 : ! sweep=1 calculats the number of resonants transitions taking into
568 : ! account a possible energy cutoff.
569 : ! sweep=2 initializes the tables describing the e-h transition.
570 : !
571 99 : do sweep=1,2
572 : !
573 66 : if (sweep==2) then
574 : ! Allocate Trans structure.
575 67 : max_nreh = MAXVAL(nreh)
576 13638 : ABI_MALLOC(Trans, (max_nreh,nsppol))
577 : end if
578 : !
579 167 : do spin=1,nsppol
580 68 : it=0
581 68 : lomo = lomo_spin(spin)
582 68 : humo = humo_spin(spin)
583 2948 : do ik_bz=1,nkbz
584 2880 : ik_ibz=ktab(ik_bz)
585 : !
586 20532 : do iv=lomo,humo
587 130656 : do ic=lomo,humo
588 110192 : delta_f = ( occ(ic,ik_ibz,spin)-occ(iv,ik_ibz,spin) ) / max_occ
589 110192 : cplx_enet = gw_energy(ic,ik_ibz,spin)-gw_energy(iv,ik_ibz,spin)
590 110192 : tene = DBLE(cplx_enet)
591 :
592 : add_transition = &
593 : (tene > tol12) .and. & ! Resonant transition.
594 : ( ABS(delta_f) > tol12) .and. & ! c-v transition.
595 110192 : (tene < uv_cut .and. tene > ir_cut) ! Energy cutoff.
596 :
597 : if (add_transition) then
598 26944 : it = it + 1
599 26944 : max_tene = MAX(max_tene, tene)
600 26944 : min_tene = MAX(min_tene, tene)
601 : end if
602 127776 : if (add_transition.and.sweep==2) then
603 13472 : Trans(it,spin)%k = ik_bz
604 13472 : Trans(it,spin)%v = iv
605 13472 : Trans(it,spin)%c = ic
606 13472 : Trans(it,spin)%en = cplx_enet
607 : end if
608 :
609 : end do
610 : end do
611 : end do ! ik_bz
612 : ! Save number of transitions for this spin.
613 134 : if (sweep==1) nreh(spin) = it
614 : end do ! spin
615 : !
616 : end do ! sweep
617 :
618 99 : minmax_tene = [min_tene, max_tene]
619 :
620 33 : end subroutine init_transitions
621 : !!***
622 :
623 : !----------------------------------------------------------------------
624 :
625 : !!****f* m_bs_defs/repr_1trans
626 : !! NAME
627 : !! repr_1trans
628 : !!
629 : !! FUNCTION
630 : !! Returns a string with info on the (k,v,c) transition.
631 : !!
632 : !! INPUTS
633 : !! Trans<transition>=structure datatype containing indececes and info on the optical transition.
634 : !! [prtvol]=Verbosity level. Defaults to 0.
635 : !!
636 : !! OUTPUT
637 : !! str(len=500)=The string representing the transition.
638 : !!
639 : !! SOURCE
640 :
641 0 : pure function repr_1trans(Trans, prtvol) result(str)
642 :
643 : !Arguments ------------------------------------
644 : class(transition),intent(in) :: Trans
645 : integer,optional,intent(in) :: prtvol
646 : character(len=500) :: str
647 :
648 : !Local variables ------------------------------
649 : !scalars
650 : integer :: my_prtvol
651 : !************************************************************************
652 :
653 0 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
654 :
655 0 : if (my_prtvol==0) then
656 0 : write(str,'(3(a,i3))')" k= ",Trans%k," v= ",Trans%v," c= ",Trans%c
657 : else
658 0 : write(str,'(3(a,i3),a,2f6.2)')" k= ",Trans%k," v= ",Trans%v," c= ",Trans%c," ene= ",Trans%en*Ha_eV
659 : end if
660 :
661 0 : end function repr_1trans
662 : !!***
663 :
664 : !----------------------------------------------------------------------
665 :
666 : !!****f* m_bs_defs/repr_2trans
667 : !! NAME
668 : !! repr_2trans
669 : !!
670 : !! FUNCTION
671 : !! Returns a string with info on two transitions
672 : !!
673 : !! INPUTS
674 : !! Trans1, Trans2<transition>=structure datatypes containing indececes and info on the optical transitions
675 : !! [prtvol]=Verbosity level. Defaults to 0.
676 : !!
677 : !! OUTPUT
678 : !! string(len=500)=The string representing the transition.
679 : !!
680 : !! SOURCE
681 :
682 0 : pure function repr_2trans(Trans1,Trans2,prtvol) result(string)
683 :
684 : !Arguments ------------------------------------
685 : !scalars
686 : integer,optional,intent(in) :: prtvol
687 : character(len=500) :: string
688 : type(transition),intent(in) :: Trans1,Trans2
689 :
690 : !Local variables ------------------------------
691 : integer :: my_prtvol
692 : !************************************************************************
693 :
694 0 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
695 0 : string = repr_1trans(Trans1,my_prtvol)//" | "//trim(repr_1trans(Trans2,my_prtvol))
696 :
697 0 : end function repr_2trans
698 : !!***
699 :
700 : !----------------------------------------------------------------------
701 :
702 : !!****f* m_bs_defs/print_bs_files
703 : !! NAME
704 : !! print_bs_files
705 : !!
706 : !! FUNCTION
707 : !! Printout of the content of the excfiles structure.
708 : !!
709 : !! INPUTS
710 : !! BS_files<excfiles>=An object of type excfile storing the filenames used in the Bethe-Salpeter code.
711 : !! [unit]=Unit number for output
712 : !! [prtvol]=Verbosity level
713 : !! [mode_paral]=Either "COLL" or "PERS"
714 : !! [header]=String to be printed as header for additional info.
715 : !!
716 : !! OUTPUT
717 : !! Only printing.
718 : !!
719 : !! SOURCE
720 :
721 29 : subroutine print_bs_files(BS_files,header,unit,mode_paral,prtvol)
722 :
723 : !Arguments ------------------------------------
724 : class(excfiles),intent(in) :: BS_files
725 : integer,optional,intent(in) :: unit,prtvol
726 : character(len=4),optional,intent(in) :: mode_paral
727 : character(len=*),optional,intent(in) :: header
728 :
729 : !Local variables ------------------------------
730 : !scalars
731 : integer :: my_unt,my_prtvol
732 : character(len=4) :: my_mode
733 : character(len=500) :: msg
734 : ! *********************************************************************
735 :
736 : !@excfiles
737 29 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
738 29 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
739 29 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
740 :
741 29 : msg=' ==== Files used for the Bethe-Salpeter calculation ==== '
742 29 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
743 29 : call wrtout(my_unt,msg,my_mode)
744 :
745 29 : if (BS_files%in_hreso /= BSE_NOFILE) then
746 7 : call wrtout(my_unt," Resonant block will be read from: "//trim(BS_files%in_hreso),my_mode)
747 : end if
748 :
749 29 : if (BS_files%in_hcoup /= BSE_NOFILE) then
750 0 : call wrtout(my_unt," Coupling block will be read from: "//trim(BS_files%in_hcoup),my_mode)
751 : end if
752 :
753 29 : if (BS_files%in_eig /= BSE_NOFILE) then
754 0 : call wrtout(my_unt," BS eigenstates will be read from: "//trim(BS_files%in_eig),my_mode)
755 : end if
756 :
757 29 : if (BS_files%in_haydock_basename /= BSE_NOFILE) then
758 0 : call wrtout(my_unt," Haydock restart files have basename: "//trim(BS_files%in_haydock_basename),my_mode)
759 : end if
760 :
761 29 : end subroutine print_bs_files
762 : !!***
763 :
764 : !----------------------------------------------------------------------
765 :
766 0 : END MODULE m_bs_defs
767 : !!***
|
