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