Line data Source code
1 : !!****m* ABINIT/m_frohlich
2 : !! NAME
3 : !! m_frohlich
4 : !!
5 : !! FUNCTION
6 : !! Description
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2018-2026 ABINIT group (VV, XG)
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 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_frohlich
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_crystal
28 : use m_ebands
29 : use m_efmas_defs
30 : use m_ifc
31 : use m_dtset
32 :
33 : use m_fstrings, only : sjoin, itoa
34 : use m_gaussian_quadrature, only : cgqf
35 :
36 : implicit none
37 :
38 : private
39 : !!***
40 :
41 : !----------------------------------------------------------------------
42 :
43 : !!****t* m_frohlich/frohlich_t
44 : !! NAME
45 : !! frohlich_t
46 : !!
47 : !! FUNCTION
48 : !! Description
49 : !!
50 : !! SOURCE
51 :
52 : type,public :: frohlich_t
53 :
54 : integer :: kind = 0
55 : ! Type of the Fr\"ohlich model
56 : ! used to access the calculations of various properties
57 : ! 0 -> ndeg != 1 && ndeg != 3, generalized Fr\"ohlich model
58 : ! 1 -> ndeg = 1, standard (possibly anisotropic) Fr\"ohlich model
59 : ! 2 -> ndeg = 3, cubic generalized Fr\"ohlich model with 3-fold degeneracy
60 :
61 : ! Geometry ------------------------------------
62 :
63 : real(dp) :: ucvol
64 : ! Real space unit cell volume
65 :
66 : real(dp) :: gmet(3,3)
67 : ! Reciprocal space metric
68 :
69 :
70 : ! Electroncic subspace -----------------------
71 :
72 : logical :: isinitel = .false.
73 : ! Flag indicating that the electronic subspace has been initialized
74 :
75 : real(dp) :: kpt(3)
76 : ! k-point characterizing the electronic subspace, i.e. the k-point at which
77 : ! the effective mass tensor is obtained (usually, CBM or VBM)
78 :
79 : integer :: ndeg
80 : ! Number of degenerate bands taken into account
81 :
82 : complex(dp), allocatable :: eig2_diag_cart(:,:,:,:)
83 : ! Band curvature double tensor in Cartesian coordinates
84 : ! (3, 3, ndeg, ndeg)
85 :
86 : real(dp) :: band_params(3)
87 : ! Parameters describing electronic bands
88 : ! The meaning of this arrayd depends on the value of the kind variable
89 : ! kind = 0 -> undefined
90 : ! kind = 1 -> inverse effetive masses along the 100, 010 and 001 directions
91 : ! kind = 2 -> Luttinger-Kohn parameters A, B, C
92 :
93 : logical :: lutt_warn(3) = .false.
94 :
95 : real(dp), allocatable :: sqrt_efmas_avg(:)
96 : ! Square root of effective mass averaged over q-sphere for each band
97 : ! (ndeg)
98 :
99 : logical, allocatable :: saddle_warn(:)
100 : ! Signals if the k-point characterizing the electronic subspace is a
101 : ! saddle-point for each of the degenerate bands
102 : ! (ndeg)
103 :
104 : real(dp) :: sqrt_efmas_tot
105 : ! Total square root of effective mass average
106 :
107 : real(dp), allocatable :: invefmas(:,:)
108 : ! Inverse electronic effective masses along
109 : ! kind = 1 -> 100, 010, 001 directions
110 : ! kind = 2 -> 100, 110, 111 directions
111 : ! (ndeg, 3)
112 :
113 :
114 : ! Vibrational (phonon) subspace ---------------
115 :
116 : logical :: isinitph = .false.
117 : ! Flag indicating that the phonon subspace has been initialized
118 :
119 : integer :: natom
120 : ! Number of atoms in the cell
121 :
122 : integer :: nqdir
123 : ! Number of points for spherical integration used to compute ZPR and other
124 : ! quantities
125 :
126 : real(dp), allocatable :: unit_qdir(:,:)
127 : ! Unit q-vectors representing reciprocal space directions used to compute
128 : ! ZPR and other quantities
129 : ! (3, nqdir)
130 :
131 : real(dp), allocatable :: weights_qdir(:)
132 : ! Gaussian quadrature weights used for spherical intragraion over q-vectors
133 : ! (nqdir)
134 :
135 : real(dp), allocatable :: dielt_qdir(:)
136 : ! High-frequency dielectric constant for each q-vector direction
137 : ! (nqdir)
138 :
139 : real(dp), allocatable :: phfreq_qdir(:,:)
140 : ! Phonon frequencies for each mode and q-vector direction
141 : ! (3*natom, nqdir)
142 :
143 : real(dp), allocatable :: polarity_qdir(:,:,:)
144 : ! Mode polarity vectors for each mode and q-vector direction
145 : ! (3, 3*natom, nqdir)
146 :
147 : real(dp), allocatable :: proj_polarity_qdir(:,:)
148 : ! Projections of mode polarity vectors for each mode and q-vector direction
149 : ! (3*natom, nqdir)
150 :
151 : real(dp), allocatable :: investar(:,:)
152 : ! Inverse effective dielectric constant for each mode and q-vector
153 : ! (3*natom, nqdir)
154 :
155 : real(dp), allocatable :: dielavg(:)
156 : ! Dielectric average over q-vectors for each mode (Eq. (26) of [deMelo2023])
157 : ! (3*natom)
158 :
159 : logical, allocatable :: isiractive(:)
160 : ! Flags to detect the infrared-active phonon modes
161 : ! (3*natom)
162 :
163 : real(dp) :: dielt_eff
164 : ! Effective dielectric constant in the strong-coupling regime
165 :
166 : real(dp) :: phfreq_eff
167 : ! Effective LO phonon frequency in the strong-coupling regime
168 :
169 :
170 : ! Weak-coupling parameters --------------------
171 :
172 : real(dp), allocatable :: zpr_band(:)
173 : ! Zero-point renormalization energy for each band (Eq. (17) of [deMelo2023])
174 : ! (ndeg)
175 :
176 : logical :: sign_warn = .false.
177 : ! Sginals an error if a saddle-point is encountered or bands contribute to
178 : ! the ZPR with different signs
179 :
180 : real(dp) :: zpr_gamma
181 : ! Correction to the ZPR taking into account the infrared divergence of the
182 : ! electron-phonon coupling in the Fr\"ohlich model
183 :
184 : real(dp) :: zpr
185 : ! Total ZPR for these bands and k-point
186 :
187 : real(dp) :: enpol_wc
188 : ! Fr\"ohlich polaron formation energy in the weak-coupling regime
189 : ! (Eq. (26) of [deMelo2023])
190 :
191 : real(dp), allocatable :: zpr_k(:,:,:)
192 : ! Direction dependent ZPR (Eq. (86) of [Guster2021])
193 : ! (3, ndeg, 3)
194 :
195 : real(dp), allocatable :: invpolmas(:,:)
196 : ! Inverse Fr\"ohlich polaron effective masses in the weak-coupling regime
197 : ! (Sec. III A, B of [Guster2021])
198 : ! kind = 1 -> 100, 010, 001 directions
199 : ! kind = 2 -> 100, 110, 111 directions
200 : ! (ndeg, 3)
201 :
202 :
203 : contains
204 :
205 : procedure :: init_ph => frohlich_init_ph
206 : ! Initialization of a vibrational (phonon) subspace parameters
207 :
208 : procedure :: init_el => frohlich_init_el
209 : ! Initialization of an electronic subspace parameters
210 :
211 : procedure :: free_ph => frohlich_free_ph
212 : ! Free memory allocated to the vibrational subspace
213 :
214 : procedure :: free_el => frohlich_free_el
215 : ! Free memory allocated to the electronic subspace
216 : ! and other related quantities
217 :
218 : procedure :: calc_zpr => frohlich_calc_zpr
219 : ! Calculate the zero-point renormalization energy corresponding to the
220 : ! weak-coupling treatment of the Fr\"ohlich model
221 :
222 : procedure :: calc_polaronmass => frohlich_calc_polaronmass
223 : ! Calculate the polaron effective mass corresponding to the weak-coupling
224 : ! treatment of the Fr\"ohlich model; available for kind = 1 or 2
225 :
226 : ! procedure :: ncwrite => frohlich_ncwrite
227 : ! Write main dimensions and header on a netcdf file
228 :
229 :
230 : end type frohlich_t
231 : !!***
232 :
233 : public :: frohlichmodel_zpr ! Main routine to compute ZPR
234 : public :: frohlichmodel_polaronmass ! Main routine to compute polaron effective
235 : ! masses
236 :
237 : contains !=====================================================================
238 : !!***
239 :
240 : !!****f* m_frohlich/frohlichmodel_polaronmass
241 : !! NAME
242 : !! frohlichmodel_polaronmass
243 : !!
244 : !! FUNCTION
245 : !! Main routine to compute the polaron effective masses of the generalized
246 : !! Fr\"ohlich model and other related quantities
247 : !!
248 : !! INPUTS
249 : !! cryst<crystal_t>=Structure defining the unit cell
250 : !! dtset<dataset_type>=All input variables for this dataset.
251 : !! efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band
252 : !! degeneracy at each k point
253 : !! efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure
254 : !! efmasval(:,:)%eig2_diag band curvature double tensor
255 : !! ifc<ifc_type>=contains the dynamical matrix and the IFCs.
256 : !!
257 : !! NOTES
258 : !! This routine has to be merged with the frohlichmodel_zpr routine, and their
259 : !! text output needs to be refined. For now, they are being kept to be
260 : !! compatible with the legacy unit tests.
261 : !!
262 : !! SOURCE
263 :
264 1 : subroutine frohlichmodel_polaronmass(frohlich, cryst, dtset, efmasdeg, efmasval, ifc)
265 :
266 : !Arguments ------------------------------------
267 : !scalars
268 : class(frohlich_t),intent(inout) :: frohlich
269 : type(crystal_t),intent(in) :: cryst
270 : type(dataset_type),intent(in) :: dtset
271 : type(ifc_type),intent(in) :: ifc
272 : !arrays
273 : type(efmasdeg_type), intent(in) :: efmasdeg(:)
274 : type(efmasval_type), intent(in) :: efmasval(:,:)
275 :
276 : !Local variables-------------------------------
277 : !scalar
278 : integer :: nu, ikpt, ideg, ndeg
279 : integer :: iqdir
280 : !arrays
281 : real(dp) :: kpt(3)
282 :
283 : ! *************************************************************************
284 :
285 : ! Initalize phonon and dielectric subspace
286 1 : call frohlich%init_ph(cryst, dtset%efmas_ntheta, ifc)
287 :
288 : ! For each k-point (with possible degeneracy), initialize an inistance of the
289 : ! generalize Fr\"ohlich model and calculate related quantities
290 2 : do ikpt=1,dtset%nkpt
291 4 : kpt(:) = dtset%kptns(:, ikpt)
292 :
293 4 : do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
294 : ndeg = efmasdeg(ikpt)%degs_bounds(2, ideg) - &
295 2 : efmasdeg(ikpt)%degs_bounds(1, ideg) + 1
296 :
297 2 : call frohlich%init_el(cryst, kpt, ndeg, efmasval(ideg, ikpt)%eig2_diag)
298 2 : call frohlich%calc_zpr()
299 2 : call frohlich%calc_polaronmass()
300 :
301 : ! Luttinger parameters if applicable
302 2 : if (frohlich%ndeg == 3) then
303 :
304 4 : if (.not. (any(frohlich%saddle_warn))) then
305 4 : if (any(frohlich%lutt_warn)) then
306 :
307 0 : write(ab_out, '(2a)') ch10, &
308 0 : ' Luttinger parameters could not be determined:'
309 :
310 0 : if (frohlich%lutt_warn(1)) then
311 : write(ab_out, '(a)') ' Predicted degeneracies for &
312 0 : °_dim = 3 are not met for (100) direction.'
313 : endif
314 :
315 0 : if (frohlich%lutt_warn(2)) then
316 : write(ab_out, '(a)') ' Predicted degeneracies for &
317 0 : °_dim = 3 are not met for (111) direction.'
318 : endif
319 :
320 0 : if (frohlich%lutt_warn(3)) then
321 : write(ab_out, '(a)') ' Predicted degeneracies for &
322 0 : °_dim = 3 are not met for (110) direction.'
323 : endif
324 :
325 0 : write(ab_out, '(a)') ch10
326 : else
327 : write(ab_out, '(a,3f14.6)') &
328 1 : ' Luttinger parameters (A, B, C) (a.u.): ', &
329 2 : frohlich%band_params(:)
330 : endif
331 : endif
332 :
333 : ! Print inverse electronic effective masses in the output
334 1 : write(ab_out, '(a)') repeat('-', 80)
335 : write(ab_out, '(a)') ' Polaron properties from the generalized &
336 1 : &Froehlich model'
337 1 : write(ab_out, '(a)') repeat('-', 80)
338 1 : write(ab_out, '(a)') ' Polar modes'
339 1 : write(ab_out, '(a)') ' ## Frequency(meV) Epsilon*'
340 :
341 : ! For a cubic material, dielectric tensor and phonon frequencies at Gamma
342 : ! do not depend on the q-vector direction
343 1 : iqdir = 1
344 4 : do nu = 4,3*cryst%natom
345 4 : if (frohlich%isiractive(nu)) then
346 1 : write(ab_out,'(2x,i3,5x,f15.6,5x,f15.6)') nu, &
347 1 : frohlich%phfreq_qdir(nu, iqdir)*Ha_eV*1000.0_dp, &
348 2 : one/frohlich%investar(nu, iqdir)
349 : endif
350 : enddo
351 :
352 1 : write(ab_out, '(a)') ' '
353 1 : write(ab_out, '(a,f10.2)') ' ZPR (meV): ', &
354 2 : frohlich%zpr_k(1, 1, 1)*Ha_eV*1000.0_dp
355 1 : write(ab_out, '(a)') ' '
356 : write(ab_out, '(a)') ' Electronic effective mass (a.u.) &
357 1 : &along 3 directions'
358 1 : write(ab_out, '(a, 3f15.6)')' Direction 100: ', &
359 5 : one/frohlich%invefmas(:, 1)
360 1 : write(ab_out, '(a, 3f15.6)')' Direction 110: ', &
361 5 : one/frohlich%invefmas(:, 2)
362 1 : write(ab_out, '(a, 3f15.6)')' Direction 111: ', &
363 5 : one/frohlich%invefmas(:, 3)
364 :
365 : ! Print inverse polaron effective masses in the output
366 1 : write(ab_out, '(a)') ' '
367 1 : write(ab_out, '(a)') ' Polaron effective mass (a.u.) along 3 directions'
368 1 : write(ab_out, '(a, 3f15.6)') ' Direction 100: ', &
369 5 : one/frohlich%invpolmas(:, 1)
370 1 : write(ab_out, '(a, 3f15.6)') ' Direction 110: ', &
371 5 : one/frohlich%invpolmas(:, 2)
372 1 : write(ab_out, '(a, 3f15.6)') ' Direction 111: ', &
373 5 : one/frohlich%invpolmas(:, 3)
374 1 : write(ab_out, '(a)')' '
375 : write(ab_out, '(a)')' Sum rule of inverse polaron masses check-up &
376 1 : &(for convergence purposes):'
377 1 : write(ab_out,'(a, 3f15.6)')' Direction 100: ', &
378 5 : sum(frohlich%invpolmas(:, 1))
379 1 : write(ab_out,'(a, 3f15.6)')' Direction 110: ', &
380 5 : sum(frohlich%invpolmas(:, 2))
381 1 : write(ab_out,'(a, 3f15.6)')' Direction 111: ', &
382 5 : sum(frohlich%invpolmas(:, 3))
383 : endif
384 :
385 3 : call frohlich%free_el()
386 : enddo
387 : enddo
388 :
389 1 : call frohlich%free_ph()
390 :
391 1 : end subroutine frohlichmodel_polaronmass
392 : !!***
393 :
394 :
395 : !!****f* m_frohlich/frohlichmodel_zpr
396 : !! NAME
397 : !! frohlichmodel_zpr
398 : !!
399 : !! FUNCTION
400 : !! Main routine to compute the ZPR of the generalized Fr\"ohlich model and
401 : !! other related quantities
402 : !!
403 : !! INPUTS
404 : !! cryst<crystal_t>=Structure defining the unit cell
405 : !! dtset<dataset_type>=All input variables for this dataset.
406 : !! efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band
407 : !! degeneracy at each k point
408 : !! efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure
409 : !! efmasval(:,:)%eig2_diag band curvature double tensor
410 : !! ifc<ifc_type>=contains the dynamical matrix and the IFCs.
411 : !!
412 : !! SOURCE
413 :
414 2 : subroutine frohlichmodel_zpr(frohlich, cryst, dtset, efmasdeg, efmasval, ifc)
415 :
416 : !Arguments ------------------------------------
417 : !scalars
418 : class(frohlich_t),intent(inout) :: frohlich
419 : type(crystal_t),intent(in) :: cryst
420 : type(dataset_type),intent(in) :: dtset
421 : type(ifc_type),intent(in) :: ifc
422 : !arrays
423 : type(efmasdeg_type), intent(in) :: efmasdeg(:)
424 : type(efmasval_type), intent(in) :: efmasval(:,:)
425 :
426 : !Local variables-------------------------------
427 : !scalar
428 : integer :: nu, ikpt, ideg, ndeg, iband
429 : !arrays
430 : real(dp) :: kpt(3)
431 :
432 : ! *************************************************************************
433 :
434 : ! Initalize phonon and dielectric subspace
435 2 : call frohlich%init_ph(cryst, dtset%efmas_ntheta, ifc)
436 :
437 : ! Dielectric average
438 2 : write(ab_out, '(a)') repeat('-', 80)
439 2 : write(ab_out, '(a)') ' Dielectric average (EQ. 25 Melo2022)'
440 2 : write(ab_out, '(a)') repeat('-', 80)
441 : write(ab_out, '(a)') &
442 2 : ' Mode <1/epsilon*SQRT(w_LO/2)> Cumulative sum'
443 :
444 11 : do nu=4,3*cryst%natom
445 9 : write(ab_out, '(i5,f28.12,f28.12)') nu, &
446 74 : frohlich%dielavg(nu), sum(frohlich%dielavg(1:nu))
447 : enddo
448 2 : write(ab_out, '(a)') repeat('-', 80)
449 :
450 : ! Infrared ZPR correction (does not depend on the electronic subspace)
451 2 : write(ab_out, '(6a,f14.6,a,f14.6,a)') ch10, &
452 : ' Rough correction to the ZPR, to take into account the missing q=0 piece &
453 2 : &using Frohlich model:', ch10, &
454 : ' (+ for occupied states, - for unoccupied states) * zpr_q0_fact / &
455 2 : &(Nqpt_full_bz)**(1/3) ', ch10, &
456 2 : ' where Nqpt_full_bz=number of q wavevectors in full BZ, and zpr_q0_fact=',&
457 4 : frohlich%zpr_gamma, ' Ha=', frohlich%zpr_gamma*Ha_eV, ' eV'
458 :
459 : ! For each k-point (with possible degeneracy), initialize an inistance of the
460 : ! generalize Fr\"ohlich model and calculate ZPR and other related quantities
461 6 : do ikpt=1,dtset%nkpt
462 16 : kpt(:) = dtset%kptns(:, ikpt)
463 :
464 17 : do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
465 : ndeg = efmasdeg(ikpt)%degs_bounds(2, ideg) - &
466 11 : efmasdeg(ikpt)%degs_bounds(1, ideg) + 1
467 :
468 11 : call frohlich%init_el(cryst, kpt, ndeg, efmasval(ideg, ikpt)%eig2_diag)
469 11 : call frohlich%calc_zpr()
470 :
471 : ! Print ZPR results for each k-point
472 11 : if (ndeg == 1) then
473 9 : write(ab_out, '(2a,3(f6.3,a),i5)') ch10, &
474 9 : ' - At k-point (', kpt(1), ',', kpt(2), ',', kpt(3), '), band ', &
475 18 : efmasdeg(ikpt)%degs_bounds(1,ideg)
476 : else
477 2 : write(ab_out, '(2a,3(f6.3,a),i5,a,i5)') ch10, &
478 2 : ' - At k-point (', kpt(1), ',', kpt(2), ',', kpt(3), '), bands ', &
479 2 : efmasdeg(ikpt)%degs_bounds(1,ideg), ' through ', &
480 4 : efmasdeg(ikpt)%degs_bounds(2,ideg)
481 : endif
482 :
483 : ! Luttinger parameters if applicable
484 11 : if (frohlich%kind == 2) then
485 :
486 4 : if (.not. (any(frohlich%saddle_warn))) then
487 4 : if (any(frohlich%lutt_warn)) then
488 :
489 0 : write(ab_out, '(2a)') ch10, &
490 0 : ' Luttinger parameters could not be determined:'
491 :
492 0 : if (frohlich%lutt_warn(1)) then
493 : write(ab_out, '(a)') ' Predicted degeneracies for &
494 0 : °_dim = 3 are not met for (100) direction.'
495 : endif
496 :
497 0 : if (frohlich%lutt_warn(2)) then
498 : write(ab_out, '(a)') ' Predicted degeneracies for &
499 0 : °_dim = 3 are not met for (111) direction.'
500 : endif
501 :
502 0 : if (frohlich%lutt_warn(3)) then
503 : write(ab_out, '(a)') ' Predicted degeneracies for &
504 0 : °_dim = 3 are not met for (110) direction.'
505 : endif
506 :
507 0 : write(ab_out, '(a)') ch10
508 : else
509 : write(ab_out, '(a,3f14.6)') &
510 1 : ' Luttinger parameters (A, B, C) [at. units]: ', &
511 2 : frohlich%band_params(:)
512 : endif
513 : endif
514 : endif
515 :
516 : ! Effective mass average and ZPR
517 25 : do iband=1,ndeg
518 25 : if (frohlich%saddle_warn(iband)) then
519 5 : write(ab_out, '(a,i5,a)') ' Band ', &
520 5 : efmasdeg(ikpt)%degs_bounds(1, ideg) + iband - 1, ' SADDLE POINT - &
521 10 : &Frohlich effective mass and ZPR cannot be defined. '
522 : else
523 9 : write(ab_out, '(a,i5,a,f14.10)') ' Band ', &
524 9 : efmasdeg(ikpt)%degs_bounds(1, ideg) + iband - 1, ' Angular average &
525 9 : &effective mass for Frohlich model (<m**0.5>)**2= ', &
526 : sign(frohlich%sqrt_efmas_avg(iband)**2, &
527 18 : frohlich%sqrt_efmas_avg(iband))
528 : endif
529 : enddo
530 :
531 11 : if (.not. frohlich%sign_warn) then
532 : write(ab_out, '(2a)') &
533 7 : ' Angular and band average effective mass and ZPR for Frohlich model.'
534 :
535 7 : write(ab_out, '(a,es16.6)') ' Value of (<<m**0.5>>)**2 = ', &
536 23 : (sum(abs(frohlich%sqrt_efmas_avg(:))) / ndeg)**2
537 :
538 7 : write(ab_out, '(a,es16.6)') ' Absolute Value of <<m**0.5>> = ', &
539 23 : (sum(abs(frohlich%sqrt_efmas_avg(:))) / ndeg)
540 :
541 : write(ab_out, '(a,es16.6,a,es16.6,a)') &
542 7 : ' ZPR from Frohlich model = ', frohlich%zpr, ' Ha=', &
543 14 : frohlich%zpr*Ha_eV,' eV'
544 : else
545 : write(ab_out, '(a)') ' Angular and band average effective mass for &
546 4 : &Frohlich model cannot be defined because of a sign problem.'
547 : endif
548 :
549 15 : call frohlich%free_el()
550 : enddo
551 : enddo
552 :
553 2 : call frohlich%free_ph()
554 :
555 2 : end subroutine frohlichmodel_zpr
556 : !!***
557 :
558 :
559 : !!****f* m_frohlich/frohlich_free_el
560 : !! NAME
561 : !! frohlich_free_el
562 : !!
563 : !! FUNCTION
564 : !! Deallocate dynamic memory related to the electronic subspace and other
565 : !! related quantities
566 : !!
567 : !! INPUTS
568 : !!
569 : !! SOURCE
570 :
571 13 : subroutine frohlich_free_el(self)
572 :
573 : !Arguments ------------------------------------
574 : class(frohlich_t),intent(inout) :: self
575 : ! *************************************************************************
576 :
577 13 : self%isinitel = .false.
578 :
579 : ! real
580 13 : ABI_SFREE(self%sqrt_efmas_avg)
581 13 : ABI_SFREE(self%zpr_band)
582 13 : ABI_SFREE(self%invpolmas)
583 13 : ABI_SFREE(self%zpr_k)
584 13 : ABI_SFREE(self%invefmas)
585 :
586 : ! complex
587 13 : ABI_SFREE(self%eig2_diag_cart)
588 :
589 : ! logical
590 13 : ABI_SFREE(self%saddle_warn)
591 :
592 13 : end subroutine frohlich_free_el
593 : !!***
594 :
595 :
596 : !!****f* m_frohlich/frohlich_free_ph
597 : !! NAME
598 : !! frohlich_free_ph
599 : !!
600 : !! FUNCTION
601 : !! Deallocate dynamic memory related to the vibrational subspace
602 : !!
603 : !! INPUTS
604 : !!
605 : !! SOURCE
606 :
607 3 : subroutine frohlich_free_ph(self)
608 :
609 : !Arguments ------------------------------------
610 : class(frohlich_t),intent(inout) :: self
611 : ! *************************************************************************
612 :
613 3 : self%isinitph = .false.
614 :
615 : ! real
616 3 : ABI_SFREE(self%unit_qdir)
617 3 : ABI_SFREE(self%weights_qdir)
618 3 : ABI_SFREE(self%dielt_qdir)
619 3 : ABI_SFREE(self%phfreq_qdir)
620 3 : ABI_SFREE(self%polarity_qdir)
621 3 : ABI_SFREE(self%proj_polarity_qdir)
622 3 : ABI_SFREE(self%investar)
623 3 : ABI_SFREE(self%dielavg)
624 :
625 : ! logical
626 3 : ABI_SFREE(self%isiractive)
627 :
628 3 : end subroutine frohlich_free_ph
629 : !!***
630 :
631 :
632 : !!****f* m_frohlich/frohlich_calc_polaronmass
633 : !! NAME
634 : !! frohlich_calc_polaronmass
635 : !!
636 : !! FUNCTION
637 : !! Description
638 : !!
639 : !! INPUTS
640 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
641 : !!
642 : !! OUTPUT
643 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
644 : !!
645 : !! SOURCE
646 :
647 2 : subroutine frohlich_calc_polaronmass(self)
648 :
649 : !Arguments ------------------------------------
650 : !scalars
651 : class(frohlich_t), intent(inout) :: self
652 :
653 : !Local variables-------------------------------
654 : !scalars
655 : integer :: ii
656 : integer :: ikdir, nkdir
657 : integer :: ik, nkgrid
658 : integer :: lwork, info
659 : integer :: nu, iband
660 : integer :: ixi, nxi
661 : integer :: iqdir
662 : integer :: sigma
663 : real(dp) :: deltak, minefmas
664 : real(dp) :: xi, qlen
665 : real(dp) :: tmp
666 : !arrays
667 : real(dp) :: id33(3, 3)
668 : real(dp) :: kpt(3), unit_kdir(3, 3)
669 : real(dp) :: k_plus_q(3)
670 2 : real(dp), allocatable :: eigenvec(:,:), eigenval(:)
671 2 : real(dp), allocatable :: invefmas(:,:)
672 2 : real(dp), allocatable :: lutt_eigenvec(:,:,:,:), lutt_eigenval(:,:,:)
673 2 : real(dp), allocatable :: lk_h(:,:)
674 2 : real(dp), allocatable :: phfreq(:), investar(:)
675 2 : real(dp), allocatable :: intsum(:,:,:,:)
676 2 : real(dp), allocatable :: zpr(:,:,:), zpr_ddk(:,:)
677 2 : real(dp), allocatable :: intsuminv(:,:)
678 2 : complex(dp), allocatable :: work(:)
679 : ! *************************************************************************
680 :
681 : ! TODO: check if the electronic and phonon parts are initialized
682 6 : ABI_MALLOC(self%invpolmas, (self%ndeg, 3))
683 :
684 2 : if (self%kind == 2) then
685 : ! Polaron effective masses in the triply-degenerate case
686 : ! (Eq. (86) of [Guster2021])
687 :
688 : ! Initialization of the diagonalization routine
689 : ! TODO: this part is common for many routines -> initialize once and store
690 : ! necessary variables as local state constants in the datatype
691 4 : ABI_MALLOC(eigenvec, (self%ndeg, self%ndeg))
692 3 : ABI_MALLOC(eigenval, (self%ndeg))
693 3 : ABI_MALLOC(work, (3*self%ndeg - 2))
694 1 : lwork = -1
695 : call dsyev('V', 'U', self%ndeg, eigenvec(:,:), self%ndeg, eigenval(:), &
696 1 : work, lwork, info)
697 1 : lwork = int(work(1))
698 1 : ABI_FREE(work)
699 3 : ABI_MALLOC(work, (lwork))
700 :
701 : ! Inverse effectiv mass tensor for 100, 110 and 111 directions
702 : ! used to obtain the material-dependent characteristic wavevector length
703 : nkdir = 3
704 4 : unit_kdir(:, 1) = (/1, 0, 0/)
705 4 : unit_kdir(:, 2) = (/1, 1, 0/)/sqrt(2.0)
706 4 : unit_kdir(:, 3) = (/1, 1, 1/)/sqrt(3.0)
707 3 : ABI_MALLOC(invefmas, (self%ndeg, nkdir))
708 2 : ABI_MALLOC(self%invefmas, (self%ndeg, nkdir))
709 :
710 : ! 100: 2A, 2B (two-fold)
711 1 : invefmas(1, 1) = two*self%band_params(1)
712 1 : invefmas(2, 1) = two*self%band_params(2)
713 1 : invefmas(3, 1) = two*self%band_params(2)
714 : ! 110: A + B + C, 2B, A + B - C
715 : invefmas(1, 2) = self%band_params(1) + self%band_params(2) + &
716 1 : self%band_params(3)
717 : invefmas(2, 2) = min(two*self%band_params(2), self%band_params(1) + &
718 1 : self%band_params(2) - self%band_params(3))
719 : invefmas(3, 2) = max(two*self%band_params(2), self%band_params(1) + &
720 1 : self%band_params(2) - self%band_params(3))
721 : ! 111: 2/3*(A + 2B + 2C), 2/3*(A + 2B - C) (two-fold)
722 : invefmas(1, 3) = two*(self%band_params(1) + two*self%band_params(2) + &
723 1 : two*self%band_params(3)) / three
724 : invefmas(2, 3) = two*(self%band_params(1) + two*self%band_params(2) - &
725 1 : self%band_params(3)) / three
726 1 : invefmas(3, 3) = invefmas(2, 3)
727 :
728 13 : self%invefmas(:,:) = invefmas(:,:)
729 :
730 : ! Setting up the parameter to obtain the second derivative of the ZPR with
731 : ! the finite differences (FD) method (Eqs. (80), (86) of [Guster2021])
732 : ! FD: number of points
733 : nkgrid = 3
734 : ! FD: material-dependent length scale <--> lowest optical frequency
735 4 : iqdir = 1
736 13 : minefmas = one / maxval(abs(invefmas(:,:)))
737 1 : deltak = sqrt(two*minefmas*self%phfreq_qdir(4, iqdir) / 1000.0)
738 :
739 : ! FD: Luttinger eigenvalues and eigenvectors in the symmetry inequivalent
740 : ! cubic directions to be used in the ZPR calculations
741 6 : ABI_MALLOC(lutt_eigenvec, (self%ndeg, self%ndeg, nkgrid, nkdir))
742 4 : ABI_MALLOC(lutt_eigenval, (self%ndeg, nkgrid, nkdir))
743 :
744 40 : lutt_eigenval(:,:,:) = zero
745 4 : do ikdir=1,nkdir
746 13 : do ik=1,nkgrid
747 36 : kpt(:) = (ik - one)*deltak*unit_kdir(:, ikdir)
748 : call lk_hamiltonian(self%band_params(:), kpt(:), &
749 9 : lutt_eigenvec(:,:, ik, ikdir))
750 : call dsyev('V', 'U', self%ndeg, lutt_eigenvec(:,:, ik, ikdir), &
751 12 : self%ndeg, lutt_eigenval(:, ik, ikdir), work(:), lwork, info)
752 : enddo
753 : enddo
754 :
755 : ! Cubic system is assumed: phonon frequencies and permitivity at Gamma
756 : ! do not depend on the q-vector direction
757 1 : iqdir = 1
758 3 : ABI_MALLOC(phfreq, (3*self%natom))
759 2 : ABI_MALLOC(investar, (3*self%natom))
760 7 : phfreq(:) = self%phfreq_qdir(:, iqdir)
761 7 : investar(:) = self%investar(:, iqdir)
762 :
763 : ! Effective mass sign
764 1 : sigma = 1
765 1 : if (self%band_params(1) < 0) sigma = -1
766 :
767 : ! 3x3 Identity matrix
768 1 : id33(:,:) = zero
769 4 : do ii=1,3
770 4 : id33(ii, ii) = one
771 : enddo
772 :
773 : ! FD: main loop
774 : ! nqidr = 2*ntheta**2 (ntheta = nxi)
775 1 : nxi = sqrt(one*self%nqdir/2)
776 6 : ABI_MALLOC(intsum, (self%ndeg, self%ndeg, nkgrid, nkdir))
777 4 : ABI_MALLOC(zpr, (nkgrid, self%ndeg, nkdir))
778 2 : ABI_MALLOC(self%zpr_k, (nkgrid, self%ndeg, nkdir))
779 4 : ABI_MALLOC(lk_h, (self%ndeg, self%ndeg))
780 3 : ABI_MALLOC(intsuminv, (self%ndeg, self%ndeg))
781 :
782 : ! Values of ZPR (Eq. (86) of [Guster2021]) to be used in FD
783 40 : zpr(:,:,:) = zero
784 4 : do nu=4,3*self%natom
785 :
786 : ! Summation over the infrared-active phonon modes
787 4 : if (self%isiractive(nu)) then
788 :
789 4 : do ikdir=1,nkdir
790 13 : do iband=1,self%ndeg
791 39 : do ik=1,nkgrid
792 108 : kpt(:) = (ik - one)*deltak*unit_kdir(:, ikdir)
793 :
794 2763 : do ixi=0,nxi
795 2727 : xi = ixi*pi / (two*nxi)
796 2727 : if (ixi == nxi) xi = xi - tol8
797 :
798 : ! q-vector length = (omega/A)^{1/2}*tan(xi)
799 2727 : qlen = sqrt(phfreq(nu) / abs(self%band_params(1)))*tan(xi)
800 54542754 : do iqdir=1,self%nqdir
801 218160000 : k_plus_q(:) = kpt(:) + qlen*self%unit_qdir(:, iqdir)
802 : call lk_hamiltonian(self%band_params(:), k_plus_q(:), &
803 54540000 : lk_h(:,:))
804 :
805 : intsum(:,:, ik, ikdir) = &
806 : abs(lutt_eigenval(iband, ik, ikdir)*id33(:,:)) - (sigma * &
807 709020000 : lk_h(:,:) + phfreq(nu)*id33(:,:))
808 :
809 54540000 : call mat3inv(intsum(:,:, ik, ikdir), intsuminv(:,:))
810 :
811 : tmp = dot_product(lutt_eigenvec(:, iband, 2, ikdir), &
812 1145340000 : matmul(intsuminv(:,:), lutt_eigenvec(:, iband, 2, ikdir)))
813 :
814 : zpr(ik, iband, ikdir) = zpr(ik, iband, ikdir) + &
815 : investar(nu)*phfreq(nu)*tmp*sqrt(phfreq(nu) / &
816 : abs(self%band_params(1))) / cos(xi)**2 * &
817 54542727 : self%weights_qdir(iqdir)
818 : enddo
819 : enddo
820 : enddo
821 : enddo
822 : enddo
823 :
824 : endif
825 : enddo
826 40 : zpr(:,:,:) = half*quarter*piinv*sigma/nxi * zpr(:,:,:)
827 40 : self%zpr_k(:,:,:) = zpr(:,:,:)
828 :
829 : ! FD: actual finite differences to obtain the second derivative of ZPR
830 : ! (Eq. (80) of [Guster2021])
831 3 : ABI_MALLOC(zpr_ddk, (self%ndeg, nkdir))
832 :
833 13 : zpr_ddk(:,:) = zero
834 4 : do ikdir=1,nkdir
835 13 : do iband=1,self%ndeg
836 : ! Copied from the previous implementation
837 : ! Does not look like a usual FD expression...
838 : zpr_ddk(iband, ikdir) = four/three *(zpr(2, iband, ikdir) - &
839 : zpr(1, iband, ikdir) - (zpr(3, iband, ikdir) - &
840 12 : zpr(1, iband, ikdir)) / 16.0_dp) * two/deltak**2
841 : enddo
842 : enddo
843 :
844 4 : do ikdir=1,nkdir
845 13 : do iband=1,self%ndeg
846 : self%invpolmas(iband, ikdir) = &
847 12 : two/deltak**2 * lutt_eigenval(iband, 2, ikdir) + zpr_ddk(iband, ikdir)
848 : enddo
849 : enddo
850 :
851 1 : ABI_FREE(eigenvec)
852 1 : ABI_FREE(eigenval)
853 1 : ABI_FREE(work)
854 1 : ABI_FREE(invefmas)
855 1 : ABI_FREE(lutt_eigenvec)
856 1 : ABI_FREE(lutt_eigenval)
857 1 : ABI_FREE(phfreq)
858 1 : ABI_FREE(investar)
859 1 : ABI_FREE(intsum)
860 1 : ABI_FREE(zpr)
861 1 : ABI_FREE(lk_h)
862 1 : ABI_FREE(zpr_ddk)
863 1 : ABI_FREE(intsuminv)
864 : endif
865 :
866 2 : end subroutine frohlich_calc_polaronmass
867 : !!***
868 :
869 :
870 : !!****f* m_frohlich/frohlich_calc_zpr
871 : !! NAME
872 : !! frohlich_calc_zpr
873 : !!
874 : !! FUNCTION
875 : !! Description
876 : !!
877 : !! INPUTS
878 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
879 : !!
880 : !! OUTPUT
881 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
882 : !!
883 : !! SOURCE
884 :
885 13 : subroutine frohlich_calc_zpr(self)
886 :
887 : !Arguments ------------------------------------
888 : !scalars
889 : class(frohlich_t), intent(inout) :: self
890 :
891 : !Local variables-------------------------------
892 : !scalars
893 : integer :: iqdir, nu
894 : integer :: iband, jband
895 : integer :: lwork, info
896 : !arrays
897 13 : real(dp), allocatable :: efmas_qdir(:,:)
898 13 : real(dp), allocatable :: invefmas_avg(:)
899 13 : real(dp), allocatable :: eigenval(:), rwork(:)
900 13 : complex(dp), allocatable :: eigenvec(:,:), work(:)
901 13 : complex(dp), allocatable :: f3d(:,:)
902 13 : logical, allocatable :: efmas_pos(:)
903 : ! *************************************************************************
904 :
905 : ! TODO: check if the electronic and phonon parts are initialized
906 :
907 : ! Initialization of the diagonalization routine for degenerate case
908 39 : ABI_MALLOC(eigenval, (self%ndeg))
909 13 : if (self%ndeg > 1) then
910 12 : ABI_MALLOC(eigenvec, (self%ndeg, self%ndeg))
911 3 : lwork = -1
912 9 : ABI_MALLOC(rwork, (3*self%ndeg - 2))
913 3 : ABI_MALLOC(work, (1))
914 : call zheev('V', 'U', self%ndeg, eigenvec(:,:), self%ndeg, eigenval(:), &
915 3 : work(:), lwork, rwork(:), info)
916 3 : lwork = int(work(1))
917 3 : ABI_FREE(work)
918 9 : ABI_MALLOC(work, (lwork))
919 : endif
920 :
921 52 : ABI_MALLOC(f3d, (self%ndeg, self%ndeg))
922 39 : ABI_MALLOC(invefmas_avg, (self%ndeg))
923 39 : ABI_MALLOC(efmas_pos, (self%ndeg))
924 52 : ABI_MALLOC(efmas_qdir, (self%ndeg, self%nqdir))
925 26 : ABI_MALLOC(self%sqrt_efmas_avg, (self%ndeg))
926 26 : ABI_MALLOC(self%saddle_warn, (self%ndeg))
927 26 : ABI_MALLOC(self%zpr_band, (self%ndeg))
928 :
929 : ! Effective mass tensor in q-vector directions
930 260013 : do iqdir=1,self%nqdir
931 :
932 : ! Band curvature tensor (inverse effective mass)
933 620000 : do iband=1,self%ndeg
934 1260000 : do jband=1,self%ndeg
935 1920000 : f3d(iband, jband) = dot_product(self%unit_qdir(:, iqdir), matmul( &
936 20840000 : self%eig2_diag_cart(:,:, iband, jband), self%unit_qdir(:, iqdir)))
937 : enddo
938 : enddo
939 :
940 : ! Band curvature tensor diagonalization
941 260000 : if (self%ndeg == 1) then
942 200000 : eigenval(1) = f3d(1, 1)
943 : else
944 660000 : eigenvec(:,:) = f3d(:,:)
945 220000 : eigenval(:) = zero
946 320000 : work(:) = zero
947 420000 : rwork(:) = zero
948 : call zheev('V', 'U', self%ndeg, eigenvec(:,:), self%ndeg, eigenval(:), &
949 60000 : work(:), lwork, rwork(:), info)
950 60000 : ABI_CHECK(info == 0, sjoin("zheev returned info: ", itoa(info)))
951 : endif
952 :
953 620013 : efmas_qdir(:, iqdir) = one/eigenval(:)
954 : enddo
955 :
956 : ! Integration over the q-sphere: square root of the effective mass average and
957 : ! ZPR for each band (Eq. (17) of [deMelo2023])
958 31 : self%saddle_warn(:) = .false.
959 31 : self%sqrt_efmas_avg(:) = zero
960 31 : self%zpr_band(:) = zero
961 31 : invefmas_avg(:) = zero
962 260013 : do iqdir=1,self%nqdir
963 : ! square root of the effective mass average
964 : self%sqrt_efmas_avg(:) = self%sqrt_efmas_avg(:) + &
965 620000 : self%weights_qdir(iqdir)*sqrt(abs(efmas_qdir(:, iqdir)))
966 :
967 : ! ZPR
968 1460000 : do nu=4,3*self%natom
969 : self%zpr_band(:) = self%zpr_band(:) + (self%weights_qdir(iqdir) * &
970 : self%investar(nu, iqdir)*sqrt(self%phfreq_qdir(nu, iqdir) * &
971 2960000 : abs(efmas_qdir(:, iqdir))))
972 : enddo
973 :
974 : ! Inverse effective mass average: used to obtain the sign of the ZPR
975 : invefmas_avg(:) = invefmas_avg(:) + &
976 620000 : self%weights_qdir(iqdir)/efmas_qdir(:, iqdir)
977 :
978 : ! Check for saddle-points
979 260018 : if (iqdir == 1) efmas_pos(:) = (efmas_qdir(:, iqdir) > 0)
980 620013 : do iband=1,self%ndeg
981 620000 : if (efmas_pos(iband) .neqv. (efmas_qdir(iband, iqdir) > 0)) then
982 26700 : self%saddle_warn(iband) = .true.
983 : endif
984 : enddo
985 : enddo
986 31 : self%sqrt_efmas_avg(:) = quarter*piinv*self%sqrt_efmas_avg(:)
987 31 : self%zpr_band(:) = sqrthalf*quarter*piinv*self%zpr_band(:)
988 31 : invefmas_avg(:) = quarter*piinv*invefmas_avg(:)
989 :
990 : ! Check for the sign problem and caclulate total ZPR and polaron formation
991 : ! energy in the weak-coupling regime
992 13 : self%sign_warn = .false.
993 31 : do iband=1,self%ndeg
994 18 : if (self%saddle_warn(iband) .or. &
995 13 : (efmas_pos(iband) .neqv. efmas_pos(1))) then
996 5 : self%sign_warn = .true.
997 : else
998 : self%sqrt_efmas_avg(iband) = &
999 13 : sign(self%sqrt_efmas_avg(iband), invefmas_avg(iband))
1000 :
1001 13 : self%zpr_band(iband) = sign(self%zpr_band(iband), invefmas_avg(iband))
1002 : endif
1003 : enddo
1004 :
1005 13 : if (.not. self%sign_warn) then
1006 22 : self%zpr = -sum(self%zpr_band(:)) / self%ndeg
1007 22 : self%sqrt_efmas_tot = sum(self%sqrt_efmas_avg(:)) / self%ndeg
1008 75 : self%enpol_wc = -self%sqrt_efmas_tot*sum(self%dielavg(:))
1009 : endif
1010 :
1011 13 : ABI_FREE(eigenval)
1012 13 : ABI_FREE(f3d)
1013 13 : ABI_FREE(invefmas_avg)
1014 13 : ABI_FREE(efmas_pos)
1015 13 : ABI_FREE(efmas_qdir)
1016 13 : if (self%ndeg > 1) then
1017 3 : ABI_FREE(eigenvec)
1018 3 : ABI_FREE(rwork)
1019 3 : ABI_FREE(work)
1020 : endif
1021 :
1022 13 : end subroutine frohlich_calc_zpr
1023 : !!***
1024 :
1025 :
1026 : !!****f* m_frohlich/frohlich_init_el
1027 : !! NAME
1028 : !! frohlich_init_el
1029 : !!
1030 : !! FUNCTION
1031 : !! Description
1032 : !!
1033 : !! INPUTS
1034 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
1035 : !! cryst<crystal_t> = Structure defining the unit cell
1036 : !! %rprimd(3, 3) = real space primitvie vectors
1037 : !! kpt(3) = k-point characterizing the electronic subspace
1038 : !! ndeg = Number of degenerate bands
1039 : !! eig2_diag(3, 3, ndeg, ndeg) = Band curvature double tensor
1040 : !!
1041 : !! OUTPUT
1042 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
1043 : !!
1044 : !! SOURCE
1045 :
1046 13 : subroutine frohlich_init_el(self, cryst, kpt, ndeg, eig2_diag)
1047 :
1048 : !Arguments ------------------------------------
1049 : !scalars
1050 : class(frohlich_t), intent(inout) :: self
1051 : type(crystal_t), intent(in) :: cryst
1052 : integer,intent(in) :: ndeg
1053 : !arrays
1054 : real(dp), intent(in) :: kpt(3)
1055 : complex(dp), intent(in) :: eig2_diag(3, 3, ndeg, ndeg)
1056 :
1057 : !Local variables-------------------------------
1058 : !scalars
1059 : integer :: lwork, info
1060 : integer :: iband, jband
1061 : integer :: idir, ipar
1062 : !arrays
1063 : real(dp) :: unit_kdir(3, 3)
1064 26 : real(dp) :: eigenval(ndeg), lutt_eigenval(ndeg, ndeg)
1065 13 : real(dp), allocatable :: rwork(:)
1066 26 : complex(dp) :: eigenvec(ndeg, ndeg), lutt_eigenvec(ndeg, ndeg)
1067 13 : complex(dp), allocatable :: work(:)
1068 : logical :: lutt_found(3)
1069 : ! *************************************************************************
1070 :
1071 52 : self%kpt(:) = kpt(:)
1072 13 : self%ndeg = ndeg
1073 :
1074 : ! Initialize the band curvatutre double tensor in Cartesian coordiantes
1075 52 : ABI_MALLOC(self%eig2_diag_cart, (3, 3, ndeg, ndeg))
1076 :
1077 31 : do iband=1,ndeg
1078 63 : do jband=1,ndeg
1079 : self%eig2_diag_cart(:,:, iband, jband) = one/two_pi**2 * &
1080 : matmul(matmul(cryst%rprimd(:,:),eig2_diag(:,:, iband, jband)), &
1081 4114 : transpose(cryst%rprimd(:,:)))
1082 : enddo
1083 : enddo
1084 :
1085 : ! Determine the type of the Fr\"ohlich model
1086 : ! TODO: check if a material actually has cubic symmetry for kind = 3?
1087 13 : if (ndeg == 1) then
1088 : ! Standard or Anisotropic Fr\"ohlich model
1089 10 : self%kind = 1
1090 :
1091 : ! Assuming the single band is parabolic in the three Cartesian directions
1092 40 : unit_kdir(:, 1) = (/1, 0, 0/)
1093 40 : unit_kdir(:, 2) = (/0, 1, 0/)
1094 40 : unit_kdir(:, 3) = (/0, 0, 1/)
1095 :
1096 40 : do idir=1,3
1097 : self%band_params(idir) = dot_product(unit_kdir(:, idir), &
1098 670 : matmul(self%eig2_diag_cart(:,:, 1, 1), unit_kdir(:, idir)))
1099 : enddo
1100 :
1101 3 : else if (ndeg == 3) then
1102 : ! Cubic generalized Fr\"ohlich model with triply degenerate bands
1103 2 : self%kind = 2
1104 :
1105 : ! Symmety inequivalent cubic directions to obtain the Luttinger parameters
1106 8 : unit_kdir(:, 1) = (/1, 0, 0/)
1107 8 : unit_kdir(:, 2) = (/1, 1, 0/)/sqrt(2.0)
1108 8 : unit_kdir(:, 3) = (/1, 1, 1/)/sqrt(3.0)
1109 :
1110 : ! Initialize the diagonalization routine
1111 2 : lwork = -1
1112 6 : ABI_MALLOC(rwork, (3*ndeg - 2))
1113 2 : ABI_MALLOC(work, (1))
1114 : call zheev('V', 'U', ndeg, eigenvec(:,:), ndeg, eigenval, work(:), lwork, &
1115 2 : rwork(:), info)
1116 2 : lwork=int(work(1))
1117 2 : ABI_FREE(work)
1118 6 : ABI_MALLOC(work, (lwork))
1119 :
1120 26 : lutt_eigenval(:, :) = zero
1121 : ! Inverse effective mass tensor in the symmetry inequivalent directions
1122 8 : do idir=1,3
1123 :
1124 24 : do iband=1,ndeg
1125 78 : do jband=1,ndeg
1126 : lutt_eigenvec(iband, jband) = dot_product(unit_kdir(:, idir), &
1127 1800 : matmul(self%eig2_diag_cart(:, :, iband, jband), unit_kdir(:, idir)))
1128 : enddo
1129 : enddo
1130 :
1131 36 : work(:) = zero
1132 48 : rwork(:) = zero
1133 : call zheev('V', 'U', ndeg, lutt_eigenvec(:,:), ndeg, &
1134 42 : lutt_eigenval(idir, :), work(:), lwork, rwork(:), info)
1135 8 : ABI_CHECK(info == 0, sjoin("zheev returned info: ", itoa(info)))
1136 : enddo
1137 :
1138 2 : ABI_FREE(work)
1139 2 : ABI_FREE(rwork)
1140 :
1141 : ! TODO: this looks ugly
1142 : ! is there any better wa to analyze and set the Luttinger parameters?
1143 :
1144 : ! Check degeneracies in the (100) direction, get A and B luttinger params
1145 : ! Inverse effective masses are: 2*A and 2*B (twofold)
1146 2 : if (abs(lutt_eigenval(1, 2) - lutt_eigenval(1, 3)) < tol5) then
1147 2 : self%band_params(1)=half*lutt_eigenval(1, 1)
1148 2 : self%band_params(2)=half*half*(lutt_eigenval(1, 2) + lutt_eigenval(1, 3))
1149 0 : else if (abs(lutt_eigenval(1, 2) - lutt_eigenval(1, 1)) < tol5) then
1150 0 : self%band_params(1) = half*half*(lutt_eigenval(1, 2) + lutt_eigenval(1, 1))
1151 0 : self%band_params(2) = half*lutt_eigenval(1, 3)
1152 : else
1153 0 : self%lutt_warn(1) = .true.
1154 : endif
1155 :
1156 : ! Check degeneracies in the (111) direction, get C luttinger parameter
1157 : ! Inverse effective masses are 2/3*(A + 2B - C) (twofold) and
1158 : ! 2/3*(A + 2B + 2C)
1159 2 : if (abs(lutt_eigenval(3, 2) - lutt_eigenval(3, 3)) < tol5) then
1160 : self%band_params(3) = self%band_params(1) + 2*self%band_params(2) - &
1161 2 : onehalf*half*(lutt_eigenval(3, 2) + lutt_eigenval(3, 3))
1162 0 : else if (abs(lutt_eigenval(3, 2) - lutt_eigenval(3, 1)) < tol5) then
1163 : self%band_params(3) = self%band_params(1) + 2*self%band_params(2) - &
1164 0 : onehalf*half*(lutt_eigenval(3, 2) + lutt_eigenval(3, 1))
1165 : else
1166 0 : self%lutt_warn(2) = .true.
1167 : endif
1168 :
1169 : ! Verifty that the (110) direction inverse effective masses are coherent
1170 : ! with the Luttinger parameters: 2*B, A + B + C, A + B - C
1171 2 : lutt_found(:) = .false.
1172 8 : do ipar=1,ndeg
1173 8 : if (abs(lutt_eigenval(2, ipar) - 2*self%band_params(2)) < tol4) then
1174 2 : lutt_found(1) = .true.
1175 4 : else if (abs(lutt_eigenval(2, ipar) - (self%band_params(1) + &
1176 : self%band_params(2) - self%band_params(3))) < tol4) then
1177 2 : lutt_found(2) = .true.
1178 2 : else if (abs(lutt_eigenval(2, ipar) - (self%band_params(1) + &
1179 : self%band_params(2) + self%band_params(3))) < tol4) then
1180 2 : lutt_found(3) = .true.
1181 : endif
1182 : enddo
1183 :
1184 8 : if (.not. (all(lutt_found))) self%lutt_warn(3) = .true.
1185 :
1186 : else
1187 : ! Generalized Fr\"ohlich model
1188 1 : self%kind = 0
1189 : endif
1190 :
1191 : ! The electronic subspace has been initialized
1192 13 : self%isinitel = .true.
1193 :
1194 13 : end subroutine frohlich_init_el
1195 : !!***
1196 :
1197 :
1198 : !!****f* m_frohlich/frohlich_init_ph
1199 : !! NAME
1200 : !! frohlich_init_ph
1201 : !!
1202 : !! FUNCTION
1203 : !! Description
1204 : !!
1205 : !! INPUTS
1206 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
1207 : !! cryst<crystal_t> = Structure defining the unit cell
1208 : !! %natom = number of atoms in the unit cell
1209 : !! %ucvol = real space unit cell volume
1210 : !! efmas_ntheta = number of points required for spherical integration used
1211 : !! to obtain the effective mass tensor [Laflamme2016]
1212 : !! ifc<ifc_type> = Dynamical matrix and interatomic force constants
1213 : !!
1214 : !! OUTPUT
1215 : !! self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
1216 : !!
1217 : !! SOURCE
1218 :
1219 3 : subroutine frohlich_init_ph(self, cryst, efmas_ntheta, ifc)
1220 :
1221 : !Arguments ------------------------------------
1222 : !scalars
1223 : class(frohlich_t), intent(inout) :: self
1224 : type(crystal_t), intent(in) :: cryst
1225 : type(ifc_type), intent(in) :: ifc
1226 : integer, intent(in) :: efmas_ntheta
1227 : !arrays
1228 :
1229 : !Local variables-------------------------------
1230 : !scalars
1231 : integer :: ntheta, nphi, nqdir
1232 : integer :: iphi, itheta, iqdir, nu
1233 : ! integer :: ikpt, ideg, iband, jband
1234 : ! integer :: ndeg
1235 : real(dp) :: weight
1236 : real(dp) :: weight_phi, phi_radians
1237 : real(dp) :: costheta, sintheta, cosphi, sinphi
1238 : !arrays
1239 : ! real(dp) :: kpt(3)
1240 3 : real(dp), allocatable :: gq_points_theta(:), gq_weights_theta(:)
1241 3 : real(dp), allocatable :: gq_points_cosphi(:), gq_points_sinphi(:)
1242 :
1243 : ! *************************************************************************
1244 :
1245 3 : self%natom = cryst%natom
1246 3 : self%ucvol = cryst%ucvol
1247 39 : self%gmet(:,:) = cryst%gmet(:,:)
1248 :
1249 : ! Initialization of integrals
1250 : ! This part allocates and initializes arrays unit_qdir(3,nqdir) and
1251 : ! weights_qdir(nqdir) used for spherical integration over q-points directions
1252 3 : ntheta = efmas_ntheta
1253 3 : nphi = 2*ntheta
1254 3 : nqdir = nphi*ntheta
1255 :
1256 9 : ABI_MALLOC(gq_points_theta, (ntheta))
1257 6 : ABI_MALLOC(gq_weights_theta, (ntheta))
1258 9 : ABI_MALLOC(gq_points_cosphi, (nphi))
1259 6 : ABI_MALLOC(gq_points_sinphi, (nphi))
1260 :
1261 3 : self%nqdir = nqdir
1262 9 : ABI_MALLOC(self%unit_qdir, (3, nqdir))
1263 9 : ABI_MALLOC(self%weights_qdir, (nqdir))
1264 :
1265 3 : call cgqf(ntheta, 1, zero, zero, zero, pi, gq_points_theta, gq_weights_theta)
1266 3 : weight_phi = two*pi/real(nphi, dp)
1267 :
1268 603 : do iphi=1,nphi
1269 600 : phi_radians = weight_phi * (iphi-1)
1270 600 : gq_points_cosphi(iphi) = cos(phi_radians)
1271 603 : gq_points_sinphi(iphi) = sin(phi_radians)
1272 : enddo
1273 :
1274 3 : nqdir = 0
1275 303 : do itheta=1,ntheta
1276 300 : costheta = cos(gq_points_theta(itheta))
1277 300 : sintheta = sin(gq_points_theta(itheta))
1278 300 : weight = gq_weights_theta(itheta)*weight_phi*sintheta
1279 :
1280 60303 : do iphi=1,nphi
1281 60000 : cosphi = gq_points_cosphi(iphi)
1282 60000 : sinphi = gq_points_sinphi(iphi)
1283 60000 : nqdir = nqdir + 1
1284 :
1285 60000 : self%unit_qdir(1, nqdir) = sintheta*cosphi
1286 60000 : self%unit_qdir(2, nqdir) = sintheta*sinphi
1287 60000 : self%unit_qdir(3, nqdir) = costheta
1288 60300 : self%weights_qdir(nqdir) = weight
1289 : enddo
1290 : enddo
1291 :
1292 3 : ABI_FREE(gq_points_theta)
1293 3 : ABI_FREE(gq_weights_theta)
1294 3 : ABI_FREE(gq_points_cosphi)
1295 3 : ABI_FREE(gq_points_sinphi)
1296 :
1297 : ! Initialization of the generalized Fr\"ohlich parameters
1298 :
1299 9 : ABI_MALLOC(self%dielt_qdir, (nqdir))
1300 12 : ABI_MALLOC(self%phfreq_qdir, (3*self%natom, nqdir))
1301 12 : ABI_MALLOC(self%polarity_qdir, (3, 3*self%natom, nqdir))
1302 9 : ABI_MALLOC(self%proj_polarity_qdir, (3*self%natom, nqdir))
1303 9 : ABI_MALLOC(self%investar, (3*self%natom, nqdir))
1304 9 : ABI_MALLOC(self%dielavg, (3*self%natom))
1305 9 : ABI_MALLOC(self%isiractive, (3*self%natom))
1306 :
1307 : ! Phonon frequencies and mode polarity vectors for each q-vector direction
1308 : call ifc%calcnwrite_nana_terms(cryst, nqdir, self%unit_qdir(:,:), &
1309 3 : phfrq2l=self%phfreq_qdir, polarity2l=self%polarity_qdir)
1310 :
1311 : ! High-frequency dielectric constant for each q-vector direction
1312 60003 : do iqdir=1,nqdir
1313 : self%dielt_qdir(iqdir) = dot_product(self%unit_qdir(:, iqdir), &
1314 960003 : matmul(ifc%dielt(:,:), self%unit_qdir(:, iqdir)))
1315 : enddo
1316 :
1317 480003 : self%investar(:,:) = zero
1318 24 : self%dielavg(:) = zero
1319 : ! Prjections of mode polarities on unit q-vectors,
1320 : ! inverse effective dielectric constants for each mode and q-vector directions
1321 : ! and dielectric average over q-vectors (Eqs. (22), (26) of [deMelo2023])
1322 60003 : do iqdir=1,nqdir
1323 300003 : do nu=4,3*self%natom
1324 : self%proj_polarity_qdir(nu, iqdir) = &
1325 960000 : dot_product(self%unit_qdir(:, iqdir), self%polarity_qdir(:, nu, iqdir))
1326 :
1327 : self%investar(nu, iqdir) = (self%proj_polarity_qdir(nu, iqdir) / &
1328 : self%dielt_qdir(iqdir) / self%phfreq_qdir(nu, iqdir))**2 * &
1329 240000 : four*pi/self%ucvol
1330 :
1331 : self%dielavg(nu) = self%dielavg(nu) + (self%weights_qdir(iqdir) * &
1332 300000 : self%investar(nu, iqdir)*sqrt(self%phfreq_qdir(nu, iqdir)))
1333 : enddo
1334 : enddo
1335 24 : self%dielavg(:) = sqrthalf*quarter*piinv * self%dielavg(:)
1336 :
1337 24 : self%isiractive(:) = .false.
1338 : ! Check for the infrared-active phonon modes: dielectric average is non-zero
1339 15 : do nu=4,3*self%natom
1340 15 : if (self%dielavg(nu) > tol10) self%isiractive(nu) = .true.
1341 : enddo
1342 :
1343 : ! Rough correction for the ZPR in the infrared limit of the Fr\"ohlich model
1344 3 : self%zpr_gamma = zero
1345 15 : do nu=4,3*self%natom
1346 240015 : do iqdir=1,self%nqdir
1347 : self%zpr_gamma = self%zpr_gamma + &
1348 240012 : self%investar(nu, iqdir)*self%weights_qdir(iqdir)
1349 : enddo
1350 : enddo
1351 : !Do not remove this useless line: a bug in the ifx2025 compiler is avoided thanks to it ...
1352 3 : write(std_out,'(a,es16.6)')' frohlich_init_ph : after the loop, self%zpr_gamma=',self%zpr_gamma
1353 3 : self%zpr_gamma = quarter*piinv*self%zpr_gamma
1354 3 : self%zpr_gamma = two*(three*quarter*piinv/self%ucvol)**third * self%zpr_gamma
1355 :
1356 : ! The phonon subspace has been initialized
1357 3 : self%isinitph = .true.
1358 :
1359 3 : end subroutine frohlich_init_ph
1360 : !!***
1361 :
1362 :
1363 : !!****f* m_frohlich/mat3inv
1364 : !! NAME
1365 : !! mat3inv
1366 : !!
1367 : !! FUNCTION
1368 : !! Inverts a 3x3 matrix of real elements
1369 : !!
1370 : !! INPUTS
1371 : !! mm(3, 3) = real matrix to be inverted
1372 : !!
1373 : !! OUTPUT
1374 : !! mit(3, 3) = inverse of the input matrix
1375 : !!
1376 : !! SOURCE
1377 :
1378 54540000 : subroutine mat3inv(mm, mit)
1379 :
1380 : !Arguments ------------------------------------
1381 : !arrays
1382 : real(dp), intent(in) :: mm(3, 3)
1383 : real(dp), intent(out) :: mit(3, 3)
1384 :
1385 : !Local variables-------------------------------
1386 : !scalars
1387 : real(dp) :: dd
1388 : character(len=500) :: msg
1389 : !arrays
1390 : real(dp) :: tt(3, 3)
1391 :
1392 : ! *************************************************************************
1393 :
1394 : ! Minors and determinant
1395 54540000 : tt(1, 1) = mm(2, 2)*mm(3, 3) - mm(3, 2)*mm(2, 3)
1396 54540000 : tt(1, 2) = mm(1, 3)*mm(3, 2) - mm(1, 2)*mm(3, 3)
1397 54540000 : tt(1, 3) = mm(1, 2)*mm(2, 3) - mm(1, 3)*mm(2, 2)
1398 54540000 : tt(2, 1) = mm(2, 3)*mm(3, 1) - mm(2, 1)*mm(3, 3)
1399 54540000 : tt(2, 2) = mm(1, 1)*mm(3, 3) - mm(3, 1)*mm(1, 3)
1400 54540000 : tt(2, 3) = mm(1, 3)*mm(2, 1) - mm(1, 1)*mm(2, 3)
1401 54540000 : tt(3, 1) = mm(2, 1)*mm(3, 2) - mm(2, 2)*mm(3, 1)
1402 54540000 : tt(3, 2) = mm(1, 2)*mm(3, 1) - mm(1, 1)*mm(3, 2)
1403 54540000 : tt(3, 3) = mm(1, 1)*mm(2, 2) - mm(2, 1)*mm(1, 2)
1404 :
1405 54540000 : dd = mm(1, 1)*tt(1, 1) + mm(2, 1)*tt(2, 1) + mm(3, 1)*tt(3, 1)
1406 :
1407 : ! Make sure the matrix is not singular
1408 54540000 : if (dd /=0) then
1409 709020000 : mit(:,:) = tt(:,:) / dd
1410 : else
1411 0 : write(msg, '(2a,2x,9(i0,1x),a)') 'Attempting to invert real array',ch10, &
1412 0 : mm,' ==> determinant is zero.'
1413 0 : ABI_ERROR(msg)
1414 : endif
1415 :
1416 54540000 : end subroutine mat3inv
1417 : !!***
1418 :
1419 :
1420 : !!****f* m_frohlich/lk_hamiltonian
1421 : !! NAME
1422 : !! lk_hamiltonian
1423 : !!
1424 : !! FUNCTION
1425 : !! Construct a Luttinger-Kohn Hamiltonian at a given k-point for a set of the
1426 : !! Luttinger-Kohn parameters A, B, C
1427 : !!
1428 : !! INPUTS
1429 : !! band_params(3) = Luttinger-Kohn parameters A, B, C
1430 : !! kpt(3) = wavevector at which the Hamiltonian is obtained (dimensional)
1431 : !!
1432 : !! OUTPUT
1433 : !! h_lk(3, 3) = Luttingher-Kohn Hamiltonian, H(k)
1434 : !!
1435 : !! SOURCE
1436 :
1437 54540009 : subroutine lk_hamiltonian(band_params, kpt, h_lk)
1438 :
1439 : !Arguments ------------------------------------
1440 : !arrays
1441 : real(dp), intent(in) :: band_params(3)
1442 : real(dp), intent(in) :: kpt(3)
1443 : real(dp), intent(out) :: h_lk(3, 3)
1444 :
1445 : ! *************************************************************************
1446 :
1447 54540009 : h_lk(1, 1) = band_params(1)*kpt(1)**2 + band_params(2)*(kpt(2)**2 + kpt(3)**2)
1448 54540009 : h_lk(2, 2) = band_params(1)*kpt(2)**2 + band_params(2)*(kpt(1)**2 + kpt(3)**2)
1449 54540009 : h_lk(3, 3) = band_params(1)*kpt(3)**2 + band_params(2)*(kpt(1)**2 + kpt(2)**2)
1450 54540009 : h_lk(1, 2) = band_params(3)*kpt(1)*kpt(2)
1451 54540009 : h_lk(1, 3) = band_params(3)*kpt(1)*kpt(3)
1452 54540009 : h_lk(2, 3) = band_params(3)*kpt(2)*kpt(3)
1453 : ! Symmetric matrix
1454 54540009 : h_lk(2, 1) = h_lk(1, 2)
1455 54540009 : h_lk(3, 1) = h_lk(1, 3)
1456 54540009 : h_lk(3, 2) = h_lk(2, 3)
1457 :
1458 54540009 : end subroutine lk_hamiltonian
1459 : !!***
1460 :
1461 55180084 : end module m_frohlich
1462 : !!***
|
