Line data Source code
1 : !!****m* ABINIT/m_efmas
2 : !! NAME
3 : !! m_efmas
4 : !!
5 : !! FUNCTION
6 : !! This module contains datatypes for efmas functionalities.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2001-2026 ABINIT group (JLJ)
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_efmas
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use netcdf
28 : use m_efmas_defs
29 : use m_nctk
30 : use m_cgtools
31 : use m_dtset
32 :
33 : use defs_abitypes, only : MPI_type
34 : use m_gaussian_quadrature, only : cgqf
35 : use m_io_tools, only : get_unit
36 :
37 : implicit none
38 :
39 : private
40 :
41 : !public procedures.
42 : public :: efmasval_free
43 : public :: efmasval_free_array
44 : public :: efmasdeg_free
45 : public :: efmasdeg_free_array
46 : public :: efmas_ncread
47 : public :: check_degeneracies
48 : public :: print_tr_efmas
49 : public :: print_efmas
50 : public :: efmas_main
51 : public :: efmas_analysis
52 :
53 : !private procedures.
54 : private :: MATMUL_ ! Workaround to make tests pass on ubu/buda slaves
55 : interface MATMUL_
56 : module procedure MATMUL_DP
57 : module procedure MATMUL_DPC
58 : end interface MATMUL_
59 :
60 : !!***
61 :
62 : CONTAINS
63 :
64 : !===========================================================
65 :
66 : !!****f* m_efmas/efmasval_free
67 : !! NAME
68 : !! efmasval_free
69 : !!
70 : !! FUNCTION
71 : !! This routine deallocates an efmasval_type.
72 : !!
73 : !! INPUTS
74 : !!
75 : !! OUTPUT
76 : !!
77 : !! SOURCE
78 :
79 398 : subroutine efmasval_free(efmasval)
80 :
81 : !Arguments ------------------------------------
82 : type(efmasval_type),intent(inout) :: efmasval
83 : ! *********************************************************************
84 :
85 398 : ABI_SFREE(efmasval%ch2c)
86 398 : ABI_SFREE(efmasval%eig2_diag)
87 :
88 398 : end subroutine efmasval_free
89 : !!***
90 :
91 : !----------------------------------------------------------------------
92 :
93 : !!****f* ABINIT/efmasval_free_array
94 : !! NAME
95 : !! efmasval_free_array
96 : !!
97 : !! FUNCTION
98 : !! This routine deallocates an efmasval_type or, optionally, an array of efmasval_type.
99 : !!
100 : !! INPUTS
101 : !!
102 : !! OUTPUT
103 : !!
104 : !! SOURCE
105 :
106 719 : subroutine efmasval_free_array(efmasval)
107 :
108 : !Arguments ------------------------------------
109 : type(efmasval_type),allocatable,intent(inout) :: efmasval(:,:)
110 :
111 : !!!Local variables-------------------------------
112 : integer :: i,j,n(2)
113 :
114 : ! *********************************************************************
115 :
116 : !XG20180810: please do not remove. Otherwise, I get an error on my Mac.
117 : !write(std_out,*)' efmasval_free_array : enter '
118 :
119 719 : if(allocated(efmasval)) then
120 60 : n=shape(efmasval)
121 292 : do i=1,n(1)
122 690 : do j=1,n(2)
123 670 : call efmasval_free(efmasval(i,j))
124 : end do
125 : end do
126 418 : ABI_FREE(efmasval)
127 : end if
128 :
129 719 : end subroutine efmasval_free_array
130 : !!***
131 :
132 : !----------------------------------------------------------------------
133 :
134 : !!****f* m_efmas/efmasdeg_free
135 : !! NAME
136 : !! efmasdeg_free
137 : !!
138 : !! FUNCTION
139 : !! This routine deallocates an efmasdeg_type.
140 : !!
141 : !! INPUTS
142 : !!
143 : !! OUTPUT
144 : !!
145 : !! SOURCE
146 :
147 29 : subroutine efmasdeg_free(efmasdeg)
148 :
149 : !Arguments ------------------------------------
150 : type(efmasdeg_type),intent(inout) :: efmasdeg
151 :
152 : ! *********************************************************************
153 :
154 29 : ABI_SFREE(efmasdeg%degs_bounds)
155 29 : ABI_SFREE(efmasdeg%ideg)
156 :
157 29 : end subroutine efmasdeg_free
158 : !!***
159 :
160 : !----------------------------------------------------------------------
161 :
162 : !!****f* m_efmas/efmasdeg_free_array
163 : !! NAME
164 : !! efmasdeg_free_array
165 : !!
166 : !! FUNCTION
167 : !! This routine deallocates an efmasdeg_type or, optionally, an array of efmasdeg_type.
168 : !!
169 : !! INPUTS
170 : !!
171 : !! OUTPUT
172 : !!
173 : !! SOURCE
174 :
175 719 : subroutine efmasdeg_free_array(efmasdeg)
176 :
177 : !Arguments ------------------------------------
178 : type(efmasdeg_type),allocatable,intent(inout) :: efmasdeg(:)
179 :
180 : !!!Local variables-------------------------------
181 : integer :: i,n
182 :
183 : ! *********************************************************************
184 719 : if(allocated(efmasdeg)) then
185 20 : n=size(efmasdeg)
186 49 : do i=1,n
187 49 : call efmasdeg_free(efmasdeg(i))
188 : end do
189 49 : ABI_FREE(efmasdeg)
190 : end if
191 :
192 719 : end subroutine efmasdeg_free_array
193 : !!***
194 :
195 : !----------------------------------------------------------------------
196 :
197 : !!****f* m_efmas/check_degeneracies
198 : !! NAME
199 : !! check_degeneracies
200 : !!
201 : !! FUNCTION
202 : !! This routine check for 0th order band degeneracies at given k-point.
203 : !!
204 : !! INPUTS
205 : !!
206 : !! OUTPUT
207 : !!
208 : !! SOURCE
209 :
210 24 : subroutine check_degeneracies(efmasdeg,bands,nband,eigen,deg_tol)
211 :
212 : !Arguments ------------------------------------
213 : type(efmasdeg_type),intent(out) :: efmasdeg
214 : integer,intent(in) :: bands(2),nband
215 : real(dp),intent(in) :: eigen(nband)
216 : real(dp),intent(in),optional :: deg_tol
217 :
218 : !!!Local variables-------------------------------
219 : integer :: deg_dim,iband, ideg
220 24 : integer, allocatable :: degs_bounds(:,:)
221 : real(dp) :: tol
222 24 : real(dp) :: eigen_tmp(nband)
223 : logical :: treated
224 :
225 : ! *********************************************************************
226 :
227 24 : tol=tol5; if(present(deg_tol)) tol=deg_tol
228 :
229 : !!! Determine sets of degenerate states in eigen0, i.e., at 0th order.
230 24 : efmasdeg%ndegs=1
231 24 : efmasdeg%nband=nband
232 72 : ABI_MALLOC(degs_bounds,(2,nband))
233 72 : ABI_MALLOC(efmasdeg%ideg, (nband))
234 1038 : degs_bounds=0; degs_bounds(1,1)=1
235 362 : efmasdeg%ideg=0; efmasdeg%ideg(1)=1
236 :
237 362 : eigen_tmp(:) = eigen(:)
238 :
239 338 : do iband=2,nband
240 314 : if (ABS(eigen_tmp(iband)-eigen_tmp(iband-1))>tol) then
241 160 : degs_bounds(2,efmasdeg%ndegs) = iband-1
242 160 : efmasdeg%ndegs=efmasdeg%ndegs+1
243 160 : degs_bounds(1,efmasdeg%ndegs) = iband
244 : end if
245 338 : efmasdeg%ideg(iband) = efmasdeg%ndegs
246 : end do
247 24 : degs_bounds(2,efmasdeg%ndegs)=nband
248 72 : ABI_MALLOC(efmasdeg%degs_bounds,(2,efmasdeg%ndegs))
249 576 : efmasdeg%degs_bounds(1:2,1:efmasdeg%ndegs) = degs_bounds(1:2,1:efmasdeg%ndegs)
250 24 : ABI_FREE(degs_bounds)
251 :
252 : !!! Determine if treated bands are part of a degeneracy at 0th order.
253 72 : efmasdeg%deg_range=0
254 24 : deg_dim=0
255 24 : treated=.false.
256 24 : write(std_out,'(a,i6)') 'Number of sets of bands for this k-point:',efmasdeg%ndegs
257 : write(std_out,'(a)') 'Set index; range of bands included in the set; is the set degenerate?(T/F); &
258 24 : & is the set treated by EFMAS?(T/F):'
259 208 : do ideg=1,efmasdeg%ndegs
260 184 : deg_dim = efmasdeg%degs_bounds(2,ideg) - efmasdeg%degs_bounds(1,ideg) + 1
261 : !If there is some level in the set that is inside the interval defined by bands(1:2), treat such set
262 : !The band range might be larger than the nband interval: it includes it, and also include degenerate states
263 184 : if(efmasdeg%degs_bounds(1,ideg)<=bands(2) .and. efmasdeg%degs_bounds(2,ideg)>=bands(1)) then
264 42 : treated = .true.
265 42 : if(efmasdeg%degs_bounds(1,ideg)<=bands(1)) then
266 24 : efmasdeg%deg_range(1) = ideg
267 : end if
268 42 : if(efmasdeg%degs_bounds(2,ideg)>=bands(2)) then
269 24 : efmasdeg%deg_range(2) = ideg
270 : end if
271 : end if
272 184 : write(std_out,'(2i6,a,i6,2l4)') ideg, efmasdeg%degs_bounds(1,ideg), ' -', efmasdeg%degs_bounds(2,ideg), &
273 392 : & (deg_dim>1), treated
274 : end do
275 :
276 : ! write(std_out,*)'ndegs=', efmasdeg%ndegs
277 : ! write(std_out,*)'degs_bounds=', efmasdeg%degs_bounds
278 : ! write(std_out,*)'ideg=', efmasdeg%ideg
279 : ! write(std_out,*)'deg_range=', efmasdeg%deg_range
280 :
281 : !!This first attempt WORKS, but only if the symmetries are enabled, see line 1578 of dfpt_looppert.F90.
282 : !use m_crystal, only : crystal_t, crystal_init, crystal_free, crystal_print
283 : !use m_esymm
284 : !integer :: timrev
285 : !character(len=132),allocatable :: title(:)
286 : !type(crystal_t) :: Cryst
287 : !type(esymm_t) :: Bsym
288 :
289 : !timrev = 1
290 : !if(dtset%istwfk(1)/=1) timrev=2
291 : !ABI_MALLOC(title,(dtset%ntypat))
292 : !title(:) = "Bloup"
293 : !call crystal_init(Cryst,dtset%spgroup,dtset%natom,dtset%npsp,dtset%ntypat,dtset%nsym,dtset%rprimd_orig(:,:,1),&
294 : !& dtset%typat,dtset%xred_orig(:,:,1),dtset%ziontypat,dtset%znucl,timrev,.false.,.false.,title,&
295 : !& dtset%symrel,dtset%tnons,dtset%symafm)
296 : !call crystal_print(Cryst)
297 : !ABI_FREE(title)
298 : !call esymm_init(Bsym,kpt_rbz(:,ikpt),Cryst,.false.,nspinor,1,mband,tol5,eigen0,dtset%tolsym)
299 : !write(std_out,*) 'DEBUG : Bsym. ndegs=',Bsym%ndegs
300 : !do iband=1,Bsym%ndegs
301 : ! write(std_out,*) Bsym%degs_bounds(:,iband)
302 : !end do
303 :
304 : !call crystal_free(Cryst)
305 : !call esymm_free(Bsym)
306 :
307 24 : end subroutine check_degeneracies
308 : !!***
309 :
310 : !----------------------------------------------------------------------
311 :
312 : !!****f* m_efmas/print_efmas
313 : !! NAME
314 : !! print_efmas
315 : !!
316 : !! FUNCTION
317 : !! This routine prints the information needed to compute rapidly the band effective masses,
318 : !! namely, the generalized second-order k-derivatives of the eigenenergies,
319 : !! see Eq.(66) of Laflamme2016.
320 : !!
321 : !! INPUTS
322 : !!
323 : !! OUTPUT
324 : !!
325 : !! SOURCE
326 :
327 17 : subroutine print_efmas(efmasdeg,efmasval,kpt,ncid)
328 :
329 : !Arguments ------------------------------------
330 : !scalars
331 : integer, intent(in) :: ncid
332 : !arrays
333 : real(dp), intent(in) :: kpt(:,:)
334 : type(efmasdeg_type), intent(in) :: efmasdeg(:)
335 : type(efmasval_type), intent(in) :: efmasval(:,:)
336 :
337 : !Local variables-------------------------------
338 : integer :: deg_dim,eig2_diag_arr_dim, ncerr
339 : integer :: iband,ideg,ideg_tot,ieig,ikpt
340 : integer :: jband,mband,ndegs_tot,nkpt,nkptdeg,nkptval
341 17 : integer, allocatable :: nband_arr(:), ndegs_arr(:), degs_range_arr(:,:)
342 17 : integer, allocatable :: ideg_arr(:,:), degs_bounds_arr(:,:)
343 17 : real(dp), allocatable :: ch2c_arr(:,:,:,:), eig2_diag_arr(:,:,:,:), max_abs_eigen1(:)
344 : character(len=500) :: msg
345 : !----------------------------------------------------------------------
346 :
347 : !XG20180519 Here, suppose that dtset%nkpt=nkpt_rbz (as done by Jonathan).
348 : !To be reexamined/corrected at the time of parallelization.
349 :
350 17 : nkptdeg=size(efmasdeg,1)
351 17 : nkptval=size(efmasval,2)
352 17 : if(nkptdeg/=nkptval) then
353 0 : write(msg,'(a,i8,a,i8,a)') ' nkptdeg and nkptval =',nkptdeg,' and ',nkptval,' differ, which is inconsistent.'
354 0 : ABI_ERROR(msg)
355 : end if
356 17 : nkpt=nkptdeg
357 17 : if(nkpt/=size(kpt,2)) then
358 0 : write(msg,'(a,i8,a,i8,a)') ' nkptdeg and nkpt =',nkptdeg,' and ',nkpt,' differ, which is inconsistent.'
359 0 : ABI_ERROR(msg)
360 : end if
361 :
362 17 : mband=size(efmasval,1)
363 :
364 : !Total number of (degenerate) sets over all k points
365 41 : ndegs_tot=sum(efmasdeg%ndegs)
366 : !Total number of generalized second-order k-derivatives
367 : eig2_diag_arr_dim=zero
368 41 : do ikpt=1,nkpt
369 83 : do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
370 42 : deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
371 66 : eig2_diag_arr_dim = eig2_diag_arr_dim + deg_dim**2
372 : enddo
373 : enddo
374 :
375 : !Allocate the arrays to be nc-written
376 51 : ABI_MALLOC(nband_arr, (nkpt) )
377 34 : ABI_MALLOC(ndegs_arr, (nkpt) )
378 51 : ABI_MALLOC(degs_range_arr, (2,nkpt) )
379 68 : ABI_MALLOC(ideg_arr, (mband,nkpt) )
380 51 : ABI_MALLOC(degs_bounds_arr, (2,ndegs_tot) )
381 51 : ABI_MALLOC(ch2c_arr, (2,3,3,eig2_diag_arr_dim) )
382 34 : ABI_MALLOC(eig2_diag_arr, (2,3,3,eig2_diag_arr_dim) )
383 51 : ABI_MALLOC(max_abs_eigen1, (nkpt))
384 :
385 : !Prepare the arrays to be nc-written
386 17 : ideg_tot=1
387 17 : ieig=1
388 41 : do ikpt=1,nkpt
389 24 : max_abs_eigen1(ikpt) = efmasdeg(ikpt)%max_abs_eigen1
390 24 : nband_arr(ikpt)=efmasdeg(ikpt)%nband
391 24 : ndegs_arr(ikpt)=efmasdeg(ikpt)%ndegs
392 72 : degs_range_arr(:,ikpt)=efmasdeg(ikpt)%deg_range(:)
393 362 : ideg_arr(:,ikpt)=0
394 362 : ideg_arr(1:efmasdeg(ikpt)%nband,ikpt)=efmasdeg(ikpt)%ideg(:)
395 208 : do ideg=1,efmasdeg(ikpt)%ndegs
396 552 : degs_bounds_arr(:,ideg_tot)=efmasdeg(ikpt)%degs_bounds(:,ideg)
397 208 : ideg_tot=ideg_tot+1
398 : enddo
399 83 : do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
400 42 : deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
401 145 : do jband=1,deg_dim
402 280 : do iband=1,deg_dim
403 2613 : ch2c_arr(1,:,:,ieig+iband-1)=real(efmasval(ideg,ikpt)%ch2c(:,:,iband,jband))
404 2613 : ch2c_arr(2,:,:,ieig+iband-1)=aimag(efmasval(ideg,ikpt)%ch2c(:,:,iband,jband))
405 2613 : eig2_diag_arr(1,:,:,ieig+iband-1)=real(efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband))
406 2692 : eig2_diag_arr(2,:,:,ieig+iband-1)=aimag(efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband))
407 : enddo
408 121 : ieig=ieig+deg_dim
409 : enddo
410 : enddo
411 : enddo
412 :
413 : !Define dimensions
414 : ncerr=nctk_def_dims(ncid, [ &
415 : & nctkdim_t("number_of_reduced_dimensions", 3), &
416 : & nctkdim_t("real_or_complex", 2), &
417 : & nctkdim_t("number_of_kpoints", nkpt), &
418 : & nctkdim_t("max_number_of_states", mband), &
419 : & nctkdim_t("total_number_of_degenerate_sets", ndegs_tot), &
420 : & nctkdim_t("eig2_diag_arr_dim", eig2_diag_arr_dim)&
421 119 : & ], defmode=.True.)
422 17 : NCF_CHECK(ncerr)
423 :
424 : ncerr = nctk_def_arrays(ncid, [ &
425 : & nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), &
426 : & nctkarr_t("number_of_states", "int", "number_of_kpoints"), &
427 : & nctkarr_t("number_of_degenerate_sets", "int", "number_of_kpoints"), &
428 : & nctkarr_t("degs_range_arr", "int", "two, number_of_kpoints"), &
429 : & nctkarr_t("ideg_arr", "int", "max_number_of_states, number_of_kpoints"), &
430 : & nctkarr_t("degs_bounds_arr", "int", "two, total_number_of_degenerate_sets"), &
431 : & nctkarr_t("max_abs_eigen1", "dp", "number_of_kpoints"), &
432 : & nctkarr_t("ch2c_arr", "dp", "real_or_complex, number_of_reduced_dimensions, number_of_reduced_dimensions, eig2_diag_arr_dim"), &
433 : & nctkarr_t("eig2_diag_arr","dp","real_or_complex, number_of_reduced_dimensions, number_of_reduced_dimensions, eig2_diag_arr_dim")&
434 170 : ])
435 17 : NCF_CHECK(ncerr)
436 :
437 : ! Write data.
438 17 : NCF_CHECK(nctk_set_datamode(ncid))
439 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "reduced_coordinates_of_kpoints"), kpt))
440 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "number_of_states"), nband_arr))
441 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "number_of_degenerate_sets"), ndegs_arr))
442 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "degs_range_arr"), degs_range_arr))
443 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ideg_arr"), ideg_arr))
444 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "degs_bounds_arr"), degs_bounds_arr))
445 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ch2c_arr"), ch2c_arr))
446 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eig2_diag_arr"), eig2_diag_arr))
447 17 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "max_abs_eigen1"), max_abs_eigen1))
448 :
449 : !Deallocate the arrays
450 17 : ABI_FREE(nband_arr)
451 17 : ABI_FREE(ndegs_arr)
452 17 : ABI_FREE(degs_range_arr)
453 17 : ABI_FREE(ideg_arr)
454 17 : ABI_FREE(degs_bounds_arr)
455 17 : ABI_FREE(ch2c_arr)
456 17 : ABI_FREE(eig2_diag_arr)
457 17 : ABI_FREE(max_abs_eigen1)
458 :
459 17 : end subroutine print_efmas
460 : !!***
461 :
462 : !----------------------------------------------------------------------
463 :
464 : !!****f* m_efmas/efmas_ncread
465 : !! NAME
466 : !! efmas_ncread
467 : !!
468 : !! FUNCTION
469 : !! This routine reads from an EFMAS NetCDF file the information needed
470 : !! to compute rapidly the band effective masses,
471 : !! namely, the generalized second-order k-derivatives of the eigenenergies,
472 : !! see Eq.(66) of Laflamme2016.
473 : !!
474 : !! INPUTS
475 : !!
476 : !! OUTPUT
477 : !!
478 : !! SOURCE
479 :
480 3 : subroutine efmas_ncread(efmasdeg,efmasval,kpt,ncid)
481 :
482 : !Arguments ------------------------------------
483 : !scalars
484 : integer,intent(in) :: ncid
485 : !arrays
486 : real(dp), allocatable,intent(out) :: kpt(:,:)
487 : type(efmasdeg_type), allocatable, intent(out) :: efmasdeg(:)
488 : type(efmasval_type), allocatable, intent(out) :: efmasval(:,:)
489 :
490 : !Local variables-------------------------------
491 : integer :: deg_dim,eig2_diag_arr_dim
492 : integer :: iband,ideg,ideg_tot,ieig,ikpt
493 : integer :: jband,mband,nband,ndegs,ndegs_tot,nkpt
494 3 : integer, allocatable :: nband_arr(:), ndegs_arr(:), degs_range_arr(:,:)
495 3 : integer, allocatable :: ideg_arr(:,:), degs_bounds_arr(:,:)
496 3 : real(dp), allocatable :: ch2c_arr(:,:,:,:), eig2_diag_arr(:,:,:,:), max_abs_eigen1(:)
497 : !----------------------------------------------------------------------
498 :
499 3 : NCF_CHECK(nctk_set_datamode(ncid))
500 3 : NCF_CHECK(nctk_get_dim(ncid, "number_of_kpoints", nkpt))
501 3 : NCF_CHECK(nctk_get_dim(ncid, "max_number_of_states", mband))
502 3 : NCF_CHECK(nctk_get_dim(ncid, "total_number_of_degenerate_sets", ndegs_tot))
503 3 : NCF_CHECK(nctk_get_dim(ncid, "eig2_diag_arr_dim", eig2_diag_arr_dim))
504 :
505 : !Allocate the arrays to be read from NetCDF file
506 9 : ABI_MALLOC(kpt, (3,nkpt) )
507 9 : ABI_MALLOC(nband_arr, (nkpt) )
508 6 : ABI_MALLOC(ndegs_arr, (nkpt) )
509 9 : ABI_MALLOC(degs_range_arr, (2,nkpt) )
510 12 : ABI_MALLOC(ideg_arr, (mband,nkpt) )
511 9 : ABI_MALLOC(degs_bounds_arr, (2,ndegs_tot) )
512 9 : ABI_MALLOC(ch2c_arr, (2,3,3,eig2_diag_arr_dim) )
513 6 : ABI_MALLOC(eig2_diag_arr, (2,3,3,eig2_diag_arr_dim) )
514 9 : ABI_MALLOC(max_abs_eigen1, (nkpt))
515 :
516 : !Read from NetCDF file
517 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "reduced_coordinates_of_kpoints"), kpt))
518 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "number_of_states"), nband_arr))
519 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "number_of_degenerate_sets"), ndegs_arr))
520 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "degs_range_arr"), degs_range_arr))
521 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "ideg_arr"), ideg_arr))
522 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "degs_bounds_arr"), degs_bounds_arr))
523 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "ch2c_arr"), ch2c_arr))
524 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "eig2_diag_arr"), eig2_diag_arr))
525 3 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "max_abs_eigen1"), max_abs_eigen1))
526 :
527 : !Prepare the efmas* datastructures
528 14 : ABI_MALLOC(efmasdeg,(nkpt))
529 77 : ABI_MALLOC(efmasval,(mband,nkpt))
530 :
531 3 : ideg_tot=1
532 3 : ieig=1
533 8 : do ikpt=1,nkpt
534 15 : efmasdeg(ikpt)%deg_range(:)=degs_range_arr(:,ikpt)
535 5 : nband=nband_arr(ikpt)
536 5 : efmasdeg(ikpt)%nband=nband
537 5 : efmasdeg(ikpt)%max_abs_eigen1 = max_abs_eigen1(ikpt)
538 15 : ABI_MALLOC(efmasdeg(ikpt)%ideg, (nband))
539 70 : efmasdeg(ikpt)%ideg=ideg_arr(1:nband,ikpt)
540 5 : ndegs=ndegs_arr(ikpt)
541 5 : efmasdeg(ikpt)%ndegs=ndegs
542 15 : ABI_MALLOC(efmasdeg(ikpt)%degs_bounds,(2,nband))
543 53 : do ideg=1,ndegs
544 135 : efmasdeg(ikpt)%degs_bounds(:,ideg)=degs_bounds_arr(:,ideg_tot)
545 45 : ideg_tot=ideg_tot+1
546 50 : if( efmasdeg(ikpt)%deg_range(1) <= ideg .and. ideg <= efmasdeg(ikpt)%deg_range(2) ) then
547 13 : deg_dim=efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
548 52 : ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(3,3,deg_dim,deg_dim))
549 39 : ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(3,3,deg_dim,deg_dim))
550 447 : efmasval(ideg,ikpt)%ch2c=zero
551 447 : efmasval(ideg,ikpt)%eig2_diag=zero
552 31 : do jband=1,deg_dim
553 50 : do iband=1,deg_dim
554 : efmasval(ideg,ikpt)%ch2c(:,:,iband,jband)=&
555 416 : & dcmplx(ch2c_arr(1,:,:,ieig+iband-1),ch2c_arr(2,:,:,ieig+iband-1))
556 : efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)=&
557 434 : & dcmplx(eig2_diag_arr(1,:,:,ieig+iband-1),eig2_diag_arr(2,:,:,ieig+iband-1))
558 : enddo
559 31 : ieig=ieig+deg_dim
560 : enddo
561 : else
562 32 : ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(0,0,0,0))
563 32 : ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(0,0,0,0))
564 : end if
565 : end do
566 : enddo
567 :
568 : !Deallocate the arrays
569 3 : ABI_FREE(nband_arr)
570 3 : ABI_FREE(ndegs_arr)
571 3 : ABI_FREE(degs_range_arr)
572 3 : ABI_FREE(ideg_arr)
573 3 : ABI_FREE(degs_bounds_arr)
574 3 : ABI_FREE(ch2c_arr)
575 3 : ABI_FREE(eig2_diag_arr)
576 3 : ABI_FREE(max_abs_eigen1)
577 :
578 3 : end subroutine efmas_ncread
579 : !!***
580 :
581 : !----------------------------------------------------------------------
582 :
583 : !!****f* m_efmas/print_tr_efmas
584 : !! NAME
585 : !! print_tr_efmas
586 : !!
587 : !! FUNCTION
588 : !! This routine prints the transport equivalent effective mass and others info
589 : !! for a degenerate set of bands
590 : !!
591 : !! INPUTS
592 : !!
593 : !! OUTPUT
594 : !!
595 : !! SOURCE
596 :
597 188 : subroutine print_tr_efmas(io_unit,kpt,band,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,efmas_tensor,ntheta, &
598 184 : & m_avg,m_avg_frohlich,saddle_warn,efmas_eigval,efmas_eigvec,transport_tensor_scale)
599 :
600 : !Arguments ------------------------------------
601 : integer, intent(in) :: io_unit, band, deg_dim, mdim, ndirs
602 : real(dp), intent(in) :: m_cart(ndirs,deg_dim), kpt(3), dirs(3,ndirs), rprimd(3,3), efmas_tensor(mdim,mdim,deg_dim)
603 : integer, intent(in) :: ntheta
604 : real(dp), intent(in) :: m_avg(deg_dim),m_avg_frohlich(deg_dim)
605 : logical, intent(in) :: saddle_warn(deg_dim)
606 : real(dp), intent(in), optional :: efmas_eigval(mdim,deg_dim)
607 : real(dp), intent(in), optional :: efmas_eigvec(mdim,mdim,deg_dim)
608 : real(dp), intent(in), optional :: transport_tensor_scale(deg_dim)
609 :
610 : !Local variables ------------------------------
611 : logical :: extras
612 : integer :: iband, adir
613 : character(len=22) :: format_eigvec
614 : character(len=500) :: msg, tmpstr
615 : real(dp) :: vec(3),mat(3,3)
616 :
617 94 : if(deg_dim>1) then
618 36 : extras = present(efmas_eigval) .and. present(efmas_eigvec)
619 36 : if(mdim==3 .and. .not. extras) then
620 0 : write(msg,'(a,l1,a,i1,a)') 'Subroutine print_tr_efmas called with degenerate=',deg_dim>1,&
621 0 : & ' and mdim=',mdim,', but missing required arguments for this case.'
622 0 : ABI_ERROR(msg)
623 : end if
624 36 : if(mdim==2 .and. .not. (extras .or. present(transport_tensor_scale))) then
625 0 : write(msg,'(a,l1,a,i1,a)') 'Subroutine print_tr_efmas called with degenerate=',deg_dim>1,&
626 0 : & ' and mdim=',mdim,', but missing required arguments for this case.'
627 0 : ABI_ERROR(msg)
628 : end if
629 : else
630 58 : extras = present(efmas_eigval) .and. present(efmas_eigvec)
631 58 : if(mdim>1 .and. .not. extras) then
632 0 : write(msg,'(a,l1,a,i1,a)') 'Subroutine print_tr_efmas called with degenerate=',deg_dim>1,&
633 0 : & ' and mdim=',mdim,', but missing required arguments for this case.'
634 0 : ABI_ERROR(msg)
635 : end if
636 : end if
637 :
638 94 : if(deg_dim>1) then
639 36 : write(io_unit,'(2a)') ch10,' COMMENTS: '
640 36 : write(io_unit,'(a,3(f6.3,a),i5,a,i5)') ' - At k-point (',kpt(1),',',kpt(2),',',kpt(3),'), bands ',band,' through ',&
641 72 : & band+deg_dim-1
642 36 : if(mdim>1) then
643 34 : write(io_unit,'(a)') ' are DEGENERATE (effective mass tensor is therefore not defined).'
644 34 : if(mdim==3) then
645 32 : write(io_unit,'(a)') ' See Section IIIB Eqs. (67)-(70) and Appendix E of PRB 93 205147 (2016).' ! [[cite:Laflamme2016]]
646 : write(io_unit,'(a,i7,a)') &
647 32 : & ' - Angular average effective mass for Frohlich model is to be averaged over degenerate bands. See later.'
648 2 : elseif(mdim==2) then
649 2 : write(io_unit,'(a)') ' - Also, 2D requested (perpendicular to Z axis).'
650 2 : write(io_unit,'(a)') ' See Section IIIB and Appendix F, Eqs. (F12)-(F14) of PRB 93 205147 (2016).' ! [[cite:Laflamme2016]]
651 : end if
652 34 : write(io_unit,'(a,i7,a)') ' - Associated theta integrals calculated with ntheta=',ntheta,' points.'
653 : else
654 2 : write(io_unit,'(a)') ' are DEGENERATE.'
655 2 : write(io_unit,'(a)') ' - Also, 1D requested (parallel to X axis).'
656 : end if
657 : end if
658 :
659 238 : if(ANY(saddle_warn)) then
660 10 : write(msg,'(2a)') ch10,'Band(s)'
661 24 : do iband=1,deg_dim
662 24 : if(saddle_warn(iband)) then
663 14 : write(tmpstr,'(i5)') band+iband-1
664 14 : msg = TRIM(msg)//' '//TRIM(tmpstr)//','
665 : end if
666 : end do
667 10 : write(tmpstr,'(6a)') ch10,'are not band extrema, but saddle points;',ch10, &
668 10 : & 'the transport equivalent formalism breaks down in these conditions.',ch10, &
669 20 : & 'The associated tensor(s) will therefore not be printed.'
670 10 : ABI_WARNING_UNIT(TRIM(msg)//TRIM(tmpstr), io_unit)
671 : end if
672 :
673 94 : if(deg_dim>1 .and. mdim>1) then
674 34 : write(msg,'(a)') ' Transport equivalent effective mass tensor'
675 : else
676 60 : write(msg,'(a)') ' Effective mass tensor'
677 : end if
678 :
679 94 : if(mdim>1) then
680 92 : write(format_eigvec,'(a,i1,a)') '(i3,',mdim,'f14.10,a,3f14.10)'
681 : end if
682 :
683 252 : do iband=1,deg_dim
684 158 : write(io_unit,'(2a,3(f6.3,a),i5)') ch10,' K-point (',kpt(1),',',kpt(2),',',kpt(3),') | band = ',band+iband-1
685 158 : write(io_unit,'(a)') trim(msg)//':'
686 158 : if(.not. saddle_warn(iband)) then
687 556 : do adir=1,mdim
688 556 : write(io_unit,'(3f26.10)') efmas_tensor(adir,:,iband)
689 : end do
690 144 : if(present(efmas_eigval))then
691 138 : write(io_unit,'(a)') trim(msg)//' eigenvalues:'
692 138 : write(io_unit,'(3f26.10)') efmas_eigval(:,iband)
693 : endif
694 : else
695 14 : write(io_unit,'(a)') ' *** SADDLE POINT: TRANSPORT EQV. EFF. MASS NOT DEFINED (see WARNING above) ***'
696 : end if
697 :
698 158 : if(mdim>1) then
699 152 : if(mdim==2 .and. deg_dim>1) then
700 6 : write(io_unit,'(a,f26.10)') 'Scaling of transport tensor (Eq. (FXX)) = ',transport_tensor_scale(iband)
701 : end if
702 152 : if(.not. saddle_warn(iband)) then
703 138 : if(io_unit == std_out) then
704 69 : write(io_unit,'(a)') trim(msg)//' eigenvectors in cartesian / reduced coord.:'
705 272 : do adir=1,mdim
706 873 : if( count( abs(efmas_eigval(adir,iband)-efmas_eigval(:,iband))<tol4 ) > 1 ) then
707 132 : write(io_unit,'(i3,a)') adir, ' Eigenvalue degenerate => eigenvector undefined'
708 : else
709 284 : vec=zero; vec(1:mdim)=efmas_eigvec(adir,:,iband)
710 923 : mat = transpose(rprimd)/two_pi
711 1349 : vec=matmul(mat,vec); vec=vec/sqrt(sum(vec**2))
712 71 : write(io_unit,format_eigvec) adir, efmas_eigvec(adir,:,iband), ' / ', vec
713 : end if
714 : end do
715 : end if
716 : end if
717 : end if
718 :
719 158 : if(mdim==3)then
720 : !An exactly zero average effective masse is artificial (or from a saddle point with symmetry). Does not print.
721 144 : if(abs(m_avg(iband))>tol8)then
722 : write(io_unit,'(a,f14.10)') &
723 144 : & ' Angular average effective mass 1/(<1/m>)= ',m_avg(iband)
724 : endif
725 144 : if(abs(m_avg_frohlich(iband))>tol8)then
726 : write(io_unit,'(a,f14.10)') &
727 144 : & ' Angular average effective mass for Frohlich model (<m**0.5>)**2= ',m_avg_frohlich(iband)
728 : endif
729 : endif
730 :
731 158 : write(io_unit,'(a)') ' Effective masses along directions: (cart. coord. / red. coord. -> eff. mass)'
732 940 : do adir=1,ndirs
733 2752 : vec=dirs(:,adir)
734 8944 : mat = transpose(rprimd)/two_pi
735 13072 : vec=matmul(mat,vec); vec=vec/sqrt(sum(vec**2))
736 846 : write(io_unit,'(i5,a,3f10.6,a,3f10.6,a,f14.10)') adir,': ', dirs(:,adir), ' / ', vec, ' -> ', m_cart(adir,iband)
737 : end do
738 : end do
739 :
740 94 : if(deg_dim>1 .and. mdim==3) then
741 32 : write(io_unit,'(2a)') ch10,&
742 64 : & ' Angular average effective mass for Frohlich model, averaged over degenerate bands.'
743 : write(io_unit,'(a,es16.6)') &
744 120 : & ' Value of (<<m**0.5>>)**2 = ',(sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim)**2
745 : write(io_unit,'(a,es16.6,a)') &
746 120 : & ' Absolute Value of <<m**0.5>> = ', sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim,ch10
747 : endif
748 :
749 186 : end subroutine print_tr_efmas
750 : !!***
751 :
752 : !----------------------------------------------------------------------
753 :
754 : !!****f* m_efmas/efmas_main
755 : !! NAME
756 : !! efmas_main
757 : !!
758 : !! FUNCTION
759 : !! This routine calculates the generalized second-order k-derivative, Eq.66 of Laflamme2016,
760 : !! in reduced coordinates.
761 : !!
762 : !! INPUTS
763 : !! cg(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz)=pw coefficients of GS wavefunctions at k.
764 : !! cg1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppo*nkpt_rbz,3,mpert) = first-order wf in G
765 : !! space for each perturbation. The wavefunction is orthogonal to the
766 : !! active space.
767 : !! dim_eig2rf = 1 if cg1_pert, gh0c1_pert and gh1c_pert are allocated.
768 : !! 0 otherwise.
769 : !! dtset = dataset structure containing the input variable of the calculation.
770 : !! efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point
771 : !! eigen0(nkpt_rbz*dtset%mband*dtset%nsppol) = 0-order eigenvalues at all K-points:
772 : !! <k,n'|H(0)|k,n'> (hartree).
773 : !! eigen1(nkpt_rbz*2*dtset%nsppol*dtset%mband**2,3,mpert) = matrix of first-order:
774 : !! <k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
775 : !! gh0c1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz,3,mpert) = matrix containing the
776 : !! vector: <G|H(0)|psi(1)>, for each perturbation.
777 : !! gh1c_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz,3,mpert)) = matrix containing the
778 : !! vector: <G|H(1)|n,k>, for each perturbation. The wavefunction is
779 : !! orthogonal to the active space.
780 : !! istwfk_pert(nkpt_rbz,3,mpert) = integer for choice of storage of wavefunction at
781 : !! each k point for each perturbation.
782 : !! mpert = maximum number of perturbations.
783 : !! mpi_enreg = information about MPI parallelization.
784 : !! nkpt_rbz = number of k-points for each perturbation.
785 : !! npwarr(nkpt_rbz,mpert) = array of numbers of plane waves for each k-point
786 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
787 : !!
788 : !! OUTPUT
789 : !!
790 : !! SIDE EFFECTS
791 : !! efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= generalized 2nd-order k-derivatives of eigenvalues
792 : !! efmasval(:,:)%ch2c INPUT : frozen wavefunction H2 contribution to generalized 2nd order k-derivatives of eigenenergy
793 : !! efmasval(:,:)%eig2_diag OUTPUT : generalized 2nd order k-derivatives of eigenenergy
794 : !!
795 : !! SOURCE
796 :
797 17 : subroutine efmas_main(cg,cg1_pert,dim_eig2rf,dtset,efmasdeg,efmasval,eigen0,&
798 17 : & eigen1,gh0c1_pert,gh1c_pert,istwfk_pert,mpert,mpi_enreg,nkpt_rbz,npwarr,rprimd)
799 :
800 : !Arguments ------------------------------------
801 : !scalars
802 : integer, intent(in) :: dim_eig2rf,mpert,nkpt_rbz
803 : type(dataset_type), intent(in) :: dtset
804 : type(MPI_type), intent(in) :: mpi_enreg
805 : !arrays
806 : integer, intent(in) :: istwfk_pert(nkpt_rbz,3,mpert)
807 : integer, intent(in) :: npwarr(nkpt_rbz,mpert)
808 : real(dp), intent(in) :: cg1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz*dim_eig2rf,3,mpert)
809 : real(dp), intent(in) :: gh0c1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz*dim_eig2rf,3,mpert)
810 : real(dp), intent(in) :: gh1c_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz*dim_eig2rf,3,mpert)
811 : real(dp), intent(in) :: eigen0(nkpt_rbz*dtset%mband*dtset%nsppol)
812 : real(dp), intent(in) :: eigen1(nkpt_rbz*2*dtset%nsppol*dtset%mband**2,3,mpert)
813 : real(dp), intent(in) :: rprimd(3,3)
814 : real(dp), intent(in) :: cg(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz)
815 : type(efmasdeg_type), allocatable,intent(inout) :: efmasdeg(:)
816 : type(efmasval_type), allocatable,intent(inout) :: efmasval(:,:)
817 :
818 : !Local variables-------------------------------
819 : logical :: degenerate, debug
820 : integer :: ipert, isppol
821 : integer :: icg2 !TODOM : Reactivate the sections for icg2 / allow choice of k-point other than the first in the list.
822 : integer :: npw_k, nband_k, nspinor, ideg, ikpt
823 : integer :: istwf_k, master,me,spaceworld
824 : integer :: band2tot_index, bandtot_index, iband, jband, kband
825 : integer :: adir,bdir, deg_dim, degl
826 : character(len=500) :: msg
827 : real(dp) :: deltae, dot2i,dot2r,dot3i,dot3r,doti,dotr
828 17 : real(dp), allocatable :: cg0(:,:), cg1_pert2(:,:),cg1_pert1(:,:)
829 17 : real(dp), allocatable :: gh1c_pert2(:,:),gh1c_pert1(:,:),gh0c1_pert1(:,:)
830 : complex(dp) :: eig2_part(3,3), eig2_ch2c(3,3), eig2_paral(3,3), eig2_gauge_change(3,3)
831 : complex(dp) :: eig1a, eig1b, g_ch
832 17 : complex(dp), allocatable :: eigen1_deg(:,:), eig2_diag(:,:,:,:), eig2_diag_cart(:,:,:,:)
833 : ! *********************************************************************
834 :
835 17 : debug = .false. ! Prints additional info to std_out
836 :
837 : ! Init parallelism
838 17 : master =0
839 17 : spaceworld=mpi_enreg%comm_cell
840 17 : me=mpi_enreg%me_kpt
841 :
842 17 : write(msg,'(4a)') ch10,&
843 17 : & ' CALCULATION OF EFFECTIVE MASSES',ch10,&
844 34 : & ' NOTE : Additional infos (eff. mass eigenvalues, eigenvectors and, if degenerate, average mass) are available in stdout.'
845 17 : call wrtout(std_out,msg,'COLL')
846 17 : call wrtout(ab_out,msg,'COLL')
847 :
848 17 : if(dtset%nsppol/=1)then
849 0 : write(msg,'(a,i3,a)') 'nsppol=',dtset%nsppol,' is not yet treated in m_efmas.'
850 0 : ABI_ERROR(msg)
851 : end if
852 17 : if(dtset%nspden/=1)then
853 0 : write(msg,'(a,i3,a)') 'nspden=',dtset%nspden,' is not yet treated in m_efmas.'
854 0 : ABI_ERROR(msg)
855 : end if
856 17 : if(dtset%efmas_deg==0) then
857 1 : write(msg,'(a)') 'efmas_deg==0 is for debugging; the results for degenerate bands will be garbage.'
858 1 : ABI_WARNING(msg)
859 1 : ABI_WARNING_UNIT(msg, ab_out)
860 : end if
861 :
862 17 : ipert = dtset%natom+1
863 17 : isppol = 1
864 :
865 17 : icg2 = 0
866 17 : band2tot_index=0
867 17 : bandtot_index=0
868 :
869 :
870 : !XG20180519 : in the original coding by Jonathan, there is a lack of care about using dtset%nkpt or nkpt_rbz ...
871 : !Not important in the sequential case (?!) but likely problematic in the parallel case.
872 41 : do ikpt=1,dtset%nkpt
873 24 : npw_k = npwarr(ikpt,ipert)
874 24 : nband_k = dtset%nband(ikpt)
875 24 : nspinor = dtset%nspinor
876 24 : efmasdeg(ikpt)%max_abs_eigen1 = zero
877 :
878 72 : ABI_MALLOC(cg1_pert2,(2,npw_k*nspinor))
879 48 : ABI_MALLOC(cg1_pert1,(2,npw_k*nspinor))
880 48 : ABI_MALLOC(gh1c_pert2,(2,npw_k*nspinor))
881 48 : ABI_MALLOC(gh1c_pert1,(2,npw_k*nspinor))
882 48 : ABI_MALLOC(gh0c1_pert1,(2,npw_k*nspinor))
883 48 : ABI_MALLOC(cg0,(2,npw_k*nspinor))
884 :
885 66 : do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
886 42 : deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
887 42 : degenerate = (deg_dim>1) .and. (dtset%efmas_deg/=0)
888 42 : degl = efmasdeg(ikpt)%degs_bounds(1,ideg)-1
889 :
890 168 : ABI_MALLOC(eigen1_deg,(deg_dim,deg_dim))
891 : !!! If treated band degenerate at 0th order, check that we are at extrema.
892 42 : if(degenerate) then
893 72 : do adir=1,3
894 204 : do iband=1,deg_dim
895 672 : do jband=1,deg_dim
896 : eigen1_deg(iband,jband) = cmplx(eigen1(2*(jband+degl)-1+(iband+degl-1)*2*nband_k,adir,ipert),&
897 618 : & eigen1(2*(jband+degl) +(iband+degl-1)*2*nband_k,adir,ipert),dp)
898 : end do
899 : end do
900 :
901 672 : efmasdeg(ikpt)%max_abs_eigen1 = max(efmasdeg(ikpt)%max_abs_eigen1, maxval(abs(eigen1_deg)))
902 690 : if (.not.(ALL(ABS(eigen1_deg)<tol5))) then
903 0 : write(msg,'(a,a)') ' Effective masses calculations require given k-point(s) to be band extrema for given bands, ',&
904 0 : & 'but max abs gradient of band(s) was found to be greater than 1e-5. Abinit will continue anyway.'
905 0 : ABI_WARNING(TRIM(msg))
906 0 : ABI_WARNING(msg)
907 0 : ABI_WARNING_UNIT(msg, ab_out)
908 : end if
909 : end do !adir=1,3
910 : end if !degenerate(1)
911 42 : ABI_FREE(eigen1_deg)
912 :
913 168 : ABI_MALLOC(eig2_diag,(3,3,deg_dim,deg_dim))
914 126 : ABI_MALLOC(eig2_diag_cart,(3,3,deg_dim,deg_dim))
915 2734 : eig2_diag = zero
916 :
917 121 : do iband=1,deg_dim
918 79 : write(std_out,*)" In the set (possibly degenerate) compute band ",iband ! This line here to avoid weird
919 136180 : cg0(:,:) = cg(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2)
920 322 : do jband=1,deg_dim
921 201 : eig2_part = zero
922 201 : eig2_ch2c = zero
923 201 : eig2_paral = zero
924 201 : eig2_gauge_change = zero
925 804 : do adir=1,3
926 603 : istwf_k = istwfk_pert(ikpt,adir,ipert)
927 2613 : do bdir=1,3
928 :
929 : ! Calculate the gauge change (to be subtracted to go from parallel transport to diagonal gauge).
930 32130 : do kband=1,nband_k
931 : !!! Equivalent to the gauge change in eig2stern.F90, but works also for other choices than the parallel gauge.
932 : eig1a = cmplx( eigen1(2*kband-1+(degl+iband-1)*2*nband_k+band2tot_index,adir,ipert), &
933 30321 : & -eigen1(2*kband+(degl+iband-1)*2*nband_k+band2tot_index,adir,ipert), kind=dp )
934 : eig1b = cmplx( eigen1(2*kband-1+(degl+jband-1)*2*nband_k+band2tot_index,bdir,ipert), &
935 30321 : & eigen1(2*kband+(degl+jband-1)*2*nband_k+band2tot_index,bdir,ipert), kind=dp )
936 30321 : g_ch = eig1a*eig1b
937 : eig1a = cmplx( eigen1(2*kband-1+(degl+iband-1)*2*nband_k+band2tot_index,bdir,ipert), &
938 30321 : & -eigen1(2*kband+(degl+iband-1)*2*nband_k+band2tot_index,bdir,ipert), kind=dp )
939 : eig1b = cmplx( eigen1(2*kband-1+(degl+jband-1)*2*nband_k+band2tot_index,adir,ipert), &
940 30321 : & eigen1(2*kband+(degl+jband-1)*2*nband_k+band2tot_index,adir,ipert), kind=dp )
941 30321 : g_ch = g_ch + eig1a*eig1b
942 :
943 30321 : deltae = eigen0(kband+bandtot_index) - eigen0((degl+iband)+bandtot_index)
944 30321 : if( kband<=degl.or.kband>degl+deg_dim) then
945 24372 : g_ch = g_ch/deltae
946 : else
947 : g_ch = zero
948 : end if
949 32130 : eig2_gauge_change(adir,bdir) = eig2_gauge_change(adir,bdir) + g_ch
950 : end do !kband
951 :
952 2881710 : cg1_pert2(:,:) = cg1_pert(:,1+(degl+jband-1)*npw_k*nspinor+icg2:(degl+jband)*npw_k*nspinor+icg2,bdir,ipert)
953 2881710 : cg1_pert1(:,:) = cg1_pert(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2,adir,ipert)
954 2881710 : gh1c_pert1(:,:) = gh1c_pert(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2,adir,ipert)
955 2881710 : gh1c_pert2(:,:) = gh1c_pert(:,1+(degl+jband-1)*npw_k*nspinor+icg2:(degl+jband)*npw_k*nspinor+icg2,bdir,ipert)
956 2881710 : gh0c1_pert1(:,:) = gh0c1_pert(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2,adir,ipert)
957 :
958 : ! The first two dotprod corresponds to: <Psi(1)|H(1)|Psi(0)> + cc.
959 : ! They are calculated using wavefunctions <Psi(1)| that are orthogonal to the active space.
960 : dotr=zero ; doti=zero
961 : call dotprod_g(dotr,doti,istwf_k,npw_k*nspinor,2,cg1_pert1,gh1c_pert2,mpi_enreg%me_g0,&
962 1809 : & mpi_enreg%comm_spinorfft)
963 : dot2r=zero ; dot2i=zero
964 : call dotprod_g(dot2r,dot2i,istwf_k,npw_k*nspinor,2,gh1c_pert1,cg1_pert2,mpi_enreg%me_g0,&
965 1809 : & mpi_enreg%comm_spinorfft)
966 :
967 : ! This dotprod corresponds to : <Psi(1)|H(0)- E(0)|Psi(1)>
968 : ! It is calculated using wavefunctions that are orthogonal to the active space.
969 : dot3r=zero ; dot3i=zero
970 : call dotprod_g(dot3r,dot3i,istwf_k,npw_k*nspinor,2,gh0c1_pert1,cg1_pert2,mpi_enreg%me_g0,&
971 1809 : & mpi_enreg%comm_spinorfft)
972 :
973 1809 : eig2_part(adir,bdir) = cmplx(dotr+dot2r+dot3r,doti+dot2i+dot3i,kind=dp)
974 : !eig2_part(adir,bdir) = cmplx(dotr+dot2r,doti+dot2i,kind=dp) !DEBUG
975 : !eig2_part(adir,bdir) = cmplx(dotr,doti,kind=dp) !DEBUG
976 :
977 2412 : eig2_ch2c(adir,bdir) = efmasval(ideg,ikpt)%ch2c(adir,bdir,iband,jband)
978 :
979 : end do !bdir
980 : end do !adir
981 :
982 804 : do adir=1,3
983 2613 : do bdir=1,3
984 2412 : eig2_paral(adir,bdir) = eig2_part(adir,bdir) + eig2_part(bdir,adir) + eig2_ch2c(adir,bdir)
985 : end do
986 : end do
987 :
988 2613 : eig2_diag(:,:,iband,jband) = eig2_paral - eig2_gauge_change
989 2613 : efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)=eig2_diag(:,:,iband,jband)
990 :
991 : !!! Decomposition of the hessian in into its different contributions.
992 : if(debug) then
993 :
994 : eig2_diag_cart(:,:,iband,jband) = matmul(matmul(rprimd,eig2_diag(:,:,iband,jband)),transpose(rprimd))/two_pi**2
995 : eig2_paral = matmul(matmul(rprimd,eig2_paral), transpose(rprimd))/two_pi**2
996 : eig2_gauge_change = matmul(matmul(rprimd,eig2_gauge_change), transpose(rprimd))/two_pi**2
997 : eig2_ch2c = matmul(matmul(rprimd,eig2_ch2c), transpose(rprimd))/two_pi**2
998 : eig2_part = matmul(matmul(rprimd,eig2_part), transpose(rprimd))/two_pi**2
999 :
1000 : write(std_out,'(a)') 'Hessian of eigenvalues = H. in parallel gauge - Gauge transformation'
1001 : do adir=1,3
1002 : write(std_out,'(3f12.8,2(a,3f12.8))')&
1003 : & real(eig2_diag_cart(adir,:,iband,jband),dp),' |',real(eig2_paral(adir,:),dp),' |',real(eig2_gauge_change(adir,:),dp)
1004 : end do
1005 : write(std_out,'(a)') 'H. in parallel gauge = Second der. of H + First derivatives + First derivatives^T'
1006 : do adir=1,3
1007 : write(std_out,'(3f12.8,2(a,3f12.8))') real(eig2_paral(adir,:),dp),' |',real(eig2_ch2c(adir,:),dp),' |', &
1008 : & real(eig2_part(adir,:),dp)
1009 : end do
1010 : end if !debug
1011 :
1012 280 : if(.not. degenerate .and. iband==jband) then
1013 29 : write(std_out,'(a,3f20.16)') 'Gradient of eigenvalues = ',&
1014 493 : & matmul(rprimd,eigen1(2*(degl+iband)-1+(degl+iband-1)*2*nband_k+band2tot_index,:,ipert))/two_pi
1015 : end if !.not.degenerate
1016 :
1017 : end do !jband
1018 : end do !iband
1019 :
1020 42 : ABI_FREE(eig2_diag)
1021 66 : ABI_FREE(eig2_diag_cart)
1022 : end do !ideg
1023 :
1024 24 : ABI_FREE(cg1_pert2)
1025 24 : ABI_FREE(cg1_pert1)
1026 24 : ABI_FREE(gh1c_pert2)
1027 24 : ABI_FREE(gh1c_pert1)
1028 24 : ABI_FREE(gh0c1_pert1)
1029 24 : ABI_FREE(cg0)
1030 :
1031 24 : icg2=icg2+npw_k*dtset%nspinor*nband_k
1032 24 : bandtot_index=bandtot_index+nband_k
1033 41 : band2tot_index=band2tot_index+2*nband_k**2
1034 : end do ! ikpt
1035 :
1036 17 : end subroutine efmas_main
1037 : !!***
1038 :
1039 : !----------------------------------------------------------------------
1040 :
1041 : !!****f* m_efmas/efmas_analysis
1042 : !! NAME
1043 : !! efmas_analysis
1044 : !!
1045 : !! FUNCTION
1046 : !! This routine analyzes the generalized second-order k-derivatives of eigenenergies,
1047 : !! and compute the effective mass tensor
1048 : !! (inverse of hessian of eigenvalues with respect to the wavevector)
1049 : !! in cartesian coordinates along different directions in k-space, or also the transport equivalent effective mass.
1050 : !!
1051 : !! INPUTS
1052 : !! dtset = dataset structure containing the input variable of the calculation.
1053 : !! efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point
1054 : !! efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure
1055 : !! efmasval(:,:)%eig2_diag band curvature double tensor
1056 : !! kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
1057 : !! mpert = maximum number of perturbations.
1058 : !! mpi_enreg = information about MPI parallelization.
1059 : !! nkpt_rbz = number of k-points for each perturbation.
1060 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1061 : !!
1062 : !! OUTPUT
1063 : !!
1064 : !! SOURCE
1065 :
1066 17 : subroutine efmas_analysis(dtset,efmasdeg,efmasval,kpt_rbz,mpi_enreg,nkpt_rbz,rprimd)
1067 :
1068 : !Arguments ------------------------------------
1069 : !scalars
1070 : integer, intent(in) :: nkpt_rbz
1071 : type(dataset_type), intent(in) :: dtset
1072 : type(MPI_type), intent(in) :: mpi_enreg
1073 : !arrays
1074 : real(dp), intent(in) :: rprimd(3,3)
1075 : real(dp), intent(in) :: kpt_rbz(3,nkpt_rbz)
1076 : type(efmasdeg_type), intent(in) :: efmasdeg(:)
1077 : type(efmasval_type), intent(in) :: efmasval(:,:)
1078 :
1079 : !Local variables-------------------------------
1080 : logical :: degenerate
1081 : logical :: debug
1082 : logical :: print_fsph
1083 17 : logical, allocatable :: saddle_warn(:), start_eigf3d_pos(:)
1084 : integer :: info, isppol, ideg,jdeg, ikpt, master,me,spaceworld
1085 : integer :: iband, jband, adir,bdir, deg_dim, degl, lwork
1086 : integer :: itheta, iphi, ntheta, nphi
1087 : integer :: mdim, cdirs, ndirs, io_unit
1088 : integer :: ipiv(3)
1089 : character(len=500) :: msg, filename
1090 : real(dp) :: cosph,costh,sinph,sinth,f3d_scal,weight
1091 : real(dp) :: gprimd(3,3)
1092 17 : real(dp), allocatable :: unit_r(:), dr_dth(:), dr_dph(:)
1093 17 : real(dp), allocatable :: eigenval(:), rwork(:)
1094 17 : real(dp), allocatable :: eigf3d(:)
1095 17 : real(dp), allocatable :: m_avg(:), m_avg_frohlich(:),m_cart(:,:)
1096 17 : real(dp), allocatable :: deigf3d_dth(:), deigf3d_dph(:)
1097 17 : real(dp), allocatable :: unit_speed(:,:), transport_tensor(:,:,:)
1098 17 : real(dp), allocatable :: cart_rotation(:,:), transport_tensor_eig(:)
1099 17 : real(dp), allocatable :: transport_eqv_m(:,:,:), transport_eqv_eigval(:,:), transport_eqv_eigvec(:,:,:)
1100 17 : real(dp), allocatable :: transport_tensor_scale(:)
1101 17 : real(dp), allocatable :: gq_points_th(:),gq_points_costh(:),gq_points_sinth(:),gq_weights_th(:)
1102 17 : real(dp), allocatable :: gq_points_ph(:),gq_points_cosph(:),gq_points_sinph(:),gq_weights_ph(:)
1103 17 : real(dp), allocatable :: dirs(:,:)
1104 17 : real(dp),allocatable :: prodr(:,:)
1105 : !real(dp), allocatable :: f3dfd(:,:,:)
1106 : complex(dp) :: matr2d(2,2)
1107 17 : complex(dp), allocatable :: eigenvec(:,:), work(:)
1108 17 : complex(dp), allocatable :: eig2_diag_cart(:,:,:,:)
1109 17 : complex(dp), allocatable :: f3d(:,:), df3d_dth(:,:), df3d_dph(:,:)
1110 17 : complex(dp), allocatable :: unitary_tr(:,:), eff_mass(:,:)
1111 17 : complex(dp),allocatable :: prodc(:,:)
1112 :
1113 : ! *********************************************************************
1114 :
1115 17 : debug = .false. ! Prints additional info to std_out
1116 17 : print_fsph = .false. ! Open a file and print the angle dependent curvature f(\theta,\phi)
1117 : ! for each band & kpts treated; 1 file per degenerate ensemble of bands.
1118 : ! Angles are those used in the numerical integration.
1119 :
1120 : ! Init parallelism
1121 17 : master =0
1122 17 : spaceworld=mpi_enreg%comm_cell
1123 17 : me=mpi_enreg%me_kpt
1124 :
1125 17 : isppol = 1
1126 :
1127 17 : mdim = dtset%efmas_dim
1128 :
1129 : !HERE ALLOCATE
1130 :
1131 17 : gprimd = rprimd
1132 17 : call dgetrf(mdim,mdim,gprimd,mdim,ipiv,info)
1133 17 : ABI_MALLOC(rwork,(3))
1134 17 : call dgetri(mdim,gprimd,mdim,ipiv,rwork,3,info)
1135 17 : ABI_FREE(rwork)
1136 425 : gprimd = two_pi*transpose(gprimd)
1137 :
1138 17 : cdirs = dtset%efmas_calc_dirs
1139 17 : ndirs = mdim
1140 17 : if(cdirs/=0) ndirs = dtset%efmas_n_dirs
1141 51 : ABI_MALLOC(dirs,(3,ndirs))
1142 17 : if(cdirs==0) then
1143 49 : dirs = zero
1144 16 : do adir=1,ndirs
1145 16 : dirs(adir,adir)=1.0_dp
1146 : end do
1147 12 : elseif(cdirs==1) then
1148 210 : dirs(:,:) = dtset%efmas_dirs(:,1:ndirs)
1149 60 : do adir=1,ndirs
1150 360 : dirs(:,adir) = dirs(:,adir)/sqrt(sum(dirs(:,adir)**2))
1151 : end do
1152 2 : elseif(cdirs==2) then
1153 18 : dirs(:,:) = matmul(gprimd,dtset%efmas_dirs(:,1:ndirs))
1154 2 : do adir=1,ndirs
1155 8 : dirs(:,adir) = dirs(:,adir)/sqrt(sum(dirs(:,adir)**2))
1156 : end do
1157 1 : elseif(cdirs==3) then
1158 3 : dirs(1,:) = sin(dtset%efmas_dirs(1,1:ndirs)*pi/180)*cos(dtset%efmas_dirs(2,1:ndirs)*pi/180)
1159 3 : dirs(2,:) = sin(dtset%efmas_dirs(1,1:ndirs)*pi/180)*sin(dtset%efmas_dirs(2,1:ndirs)*pi/180)
1160 3 : dirs(3,:) = cos(dtset%efmas_dirs(1,1:ndirs)*pi/180)
1161 : end if
1162 :
1163 : !!! Initialization of integrals for the degenerate case.
1164 17 : ntheta = dtset%efmas_ntheta
1165 17 : nphi = 2*ntheta
1166 51 : ABI_MALLOC(gq_points_th,(ntheta))
1167 34 : ABI_MALLOC(gq_points_costh,(ntheta))
1168 34 : ABI_MALLOC(gq_points_sinth,(ntheta))
1169 34 : ABI_MALLOC(gq_weights_th,(ntheta))
1170 51 : ABI_MALLOC(gq_points_ph,(nphi))
1171 34 : ABI_MALLOC(gq_points_cosph,(nphi))
1172 34 : ABI_MALLOC(gq_points_sinph,(nphi))
1173 34 : ABI_MALLOC(gq_weights_ph,(nphi))
1174 17 : call cgqf(ntheta,1,zero,zero,zero,pi,gq_points_th,gq_weights_th)
1175 : !XG180501 : TODO : There is no need to make a Gauss-Legendre integral for the phi variable,
1176 : !since the function to be integrated is periodic...
1177 17 : call cgqf(nphi,1,zero,zero,zero,2*pi,gq_points_ph,gq_weights_ph)
1178 1717 : do itheta=1,ntheta
1179 1700 : gq_points_costh(itheta)=cos(gq_points_th(itheta))
1180 1717 : gq_points_sinth(itheta)=sin(gq_points_th(itheta))
1181 : enddo
1182 3417 : do iphi=1,nphi
1183 3400 : gq_points_cosph(iphi)=cos(gq_points_ph(iphi))
1184 3417 : gq_points_sinph(iphi)=sin(gq_points_ph(iphi))
1185 : enddo
1186 :
1187 68 : ABI_MALLOC(eff_mass,(mdim,mdim))
1188 :
1189 : !XG20180519 : incoherent, efmasdeg is dimensioned at nkpt_rbz, and not at dtset%nkpt ...
1190 41 : do ikpt=1,dtset%nkpt
1191 83 : do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
1192 :
1193 42 : deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
1194 42 : degenerate = (deg_dim>1) .and. (dtset%efmas_deg/=0)
1195 42 : degl = efmasdeg(ikpt)%degs_bounds(1,ideg)-1
1196 :
1197 : !!! Allocations
1198 168 : ABI_MALLOC(eigenvec,(deg_dim,deg_dim))
1199 126 : ABI_MALLOC(eigenval,(deg_dim))
1200 :
1201 168 : ABI_MALLOC(eig2_diag_cart,(3,3,deg_dim,deg_dim))
1202 :
1203 121 : do iband=1,deg_dim
1204 79 : write(std_out,*)" Compute band ",iband ! This line here to avoid weird
1205 322 : do jband=1,deg_dim
1206 :
1207 2463 : eff_mass=zero
1208 :
1209 2613 : eig2_diag_cart(:,:,iband,jband)=efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)
1210 28140 : eig2_diag_cart(:,:,iband,jband) = matmul(matmul(rprimd,eig2_diag_cart(:,:,iband,jband)),transpose(rprimd))/two_pi**2
1211 :
1212 :
1213 280 : if(.not. degenerate .and. iband==jband) then
1214 :
1215 : !Compute effective mass tensor from second derivative matrix. Simple inversion.
1216 371 : eff_mass(:,:) = eig2_diag_cart(1:mdim,1:mdim,iband,jband)
1217 29 : call zgetrf(mdim,mdim,eff_mass(1:mdim,1:mdim),mdim,ipiv,info)
1218 29 : ABI_MALLOC(work,(3))
1219 29 : call zgetri(mdim,eff_mass(1:mdim,1:mdim),mdim,ipiv,work,3,info)
1220 29 : ABI_FREE(work)
1221 :
1222 : !DIAGONALIZATION
1223 145 : ABI_MALLOC(transport_eqv_eigvec,(mdim,mdim,deg_dim))
1224 608 : transport_eqv_eigvec=zero
1225 116 : ABI_MALLOC(transport_eqv_eigval,(mdim,deg_dim))
1226 208 : transport_eqv_eigval=zero
1227 371 : transport_eqv_eigvec(:,:,iband) = real(eff_mass(1:mdim,1:mdim),dp)
1228 29 : lwork=-1
1229 29 : ABI_MALLOC(rwork,(1))
1230 29 : call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_eqv_eigval(:,iband),rwork,lwork,info)
1231 29 : lwork = max(1, 3*mdim-1) ! lwork >= max(1, 3*mdim-1)
1232 29 : ABI_FREE(rwork)
1233 :
1234 87 : ABI_MALLOC(rwork,(lwork))
1235 258 : rwork=zero
1236 29 : call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_eqv_eigval(:,iband),rwork,lwork,info)
1237 29 : ABI_FREE(rwork)
1238 713 : transport_eqv_eigvec(:,:,iband) = transpose(transport_eqv_eigvec(:,:,iband)) !So that lines contain eigenvectors.
1239 :
1240 : !Frohlich average effective mass
1241 29 : ABI_MALLOC(m_avg,(1))
1242 29 : ABI_MALLOC(m_avg_frohlich,(1))
1243 29 : ABI_MALLOC(saddle_warn,(1))
1244 87 : ABI_MALLOC(unit_r,(mdim))
1245 29 : ABI_MALLOC(start_eigf3d_pos,(1))
1246 :
1247 58 : m_avg=zero
1248 58 : m_avg_frohlich=zero
1249 58 : saddle_warn=.false.
1250 :
1251 29 : if(mdim==3)then
1252 : !One has to perform the integral over the sphere
1253 2828 : do itheta=1,ntheta
1254 2800 : costh=gq_points_costh(itheta) ; sinth=gq_points_sinth(itheta)
1255 562828 : do iphi=1,nphi
1256 560000 : cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
1257 560000 : weight=gq_weights_th(itheta)*gq_weights_ph(iphi)
1258 :
1259 560000 : unit_r(1)=sinth*cosph
1260 560000 : unit_r(2)=sinth*sinph
1261 560000 : unit_r(3)=costh
1262 :
1263 17920000 : f3d_scal=dot_product(unit_r(:),matmul(real(eig2_diag_cart(:,:,iband,jband),dp),unit_r(:)))
1264 1120000 : m_avg = m_avg + weight*sinth*f3d_scal
1265 1120000 : m_avg_frohlich = m_avg_frohlich + weight*sinth/(abs(f3d_scal)**half)
1266 :
1267 560028 : if(itheta==1 .and. iphi==1) start_eigf3d_pos = f3d_scal > 0
1268 562800 : if(start_eigf3d_pos(1) .neqv. (f3d_scal>0)) then
1269 26755 : saddle_warn(1)=.true.
1270 : end if
1271 : enddo
1272 : enddo
1273 56 : m_avg = quarter/pi*m_avg
1274 56 : m_avg = one/m_avg
1275 56 : m_avg_frohlich = quarter/pi*m_avg_frohlich
1276 56 : m_avg_frohlich = m_avg_frohlich**2
1277 28 : m_avg_frohlich(1) = DSIGN(m_avg_frohlich(1),m_avg(1))
1278 :
1279 : endif ! mdim==3
1280 :
1281 : !EFMAS_DIRS
1282 116 : ABI_MALLOC(m_cart,(ndirs,deg_dim))
1283 261 : m_cart=zero
1284 168 : do adir=1,ndirs
1285 4338 : m_cart(adir,1)=1.0_dp/dot_product(dirs(:,adir),matmul(real(eig2_diag_cart(:,:,iband,jband),dp),dirs(:,adir)))
1286 : end do
1287 :
1288 : !PRINTING RESULTS
1289 : call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+iband,1,mdim,ndirs,dirs,m_cart,rprimd,real(eff_mass,dp), &
1290 : & ntheta,m_avg,m_avg_frohlich,saddle_warn,&
1291 371 : & transport_eqv_eigval(:,iband:iband),transport_eqv_eigvec(:,:,iband:iband))
1292 : call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+iband,1,mdim,ndirs,dirs,m_cart,rprimd,real(eff_mass,dp), &
1293 : & ntheta,m_avg,m_avg_frohlich,saddle_warn,&
1294 371 : & transport_eqv_eigval(:,iband:iband),transport_eqv_eigvec(:,:,iband:iband))
1295 29 : ABI_FREE(m_cart)
1296 29 : ABI_FREE(transport_eqv_eigvec)
1297 29 : ABI_FREE(transport_eqv_eigval)
1298 29 : ABI_FREE(m_avg)
1299 29 : ABI_FREE(m_avg_frohlich)
1300 29 : ABI_FREE(unit_r)
1301 29 : ABI_FREE(saddle_warn)
1302 29 : ABI_FREE(start_eigf3d_pos)
1303 :
1304 : end if !.not.degenerate
1305 : end do !jband
1306 : end do !iband
1307 :
1308 : !!! EQV_MASS
1309 42 : if(degenerate .and. mdim==3) then
1310 96 : ABI_CALLOC(unit_r,(mdim))
1311 80 : ABI_CALLOC(dr_dth,(mdim))
1312 80 : ABI_CALLOC(dr_dph,(mdim))
1313 230 : ABI_CALLOC(f3d,(deg_dim,deg_dim))
1314 230 : ABI_CALLOC(df3d_dth,(deg_dim,deg_dim))
1315 230 : ABI_CALLOC(df3d_dph,(deg_dim,deg_dim))
1316 230 : ABI_CALLOC(unitary_tr,(deg_dim,deg_dim))
1317 76 : ABI_CALLOC(eigf3d,(deg_dim))
1318 48 : ABI_MALLOC(saddle_warn,(deg_dim))
1319 32 : ABI_MALLOC(start_eigf3d_pos,(deg_dim))
1320 76 : ABI_CALLOC(m_avg,(deg_dim))
1321 76 : ABI_CALLOC(m_avg_frohlich,(deg_dim))
1322 304 : ABI_CALLOC(m_cart,(ndirs,deg_dim))
1323 76 : ABI_CALLOC(deigf3d_dth,(deg_dim))
1324 76 : ABI_CALLOC(deigf3d_dph,(deg_dim))
1325 240 : ABI_CALLOC(unit_speed,(mdim,deg_dim))
1326 652 : ABI_CALLOC(transport_tensor,(mdim,mdim,deg_dim))
1327 80 : ABI_CALLOC(transport_tensor_eig,(mdim))
1328 636 : ABI_CALLOC(transport_eqv_m,(mdim,mdim,deg_dim))
1329 224 : ABI_CALLOC(transport_eqv_eigval,(mdim,deg_dim))
1330 636 : ABI_CALLOC(transport_eqv_eigvec,(mdim,mdim,deg_dim))
1331 48 : ABI_MALLOC(prodc,(deg_dim,deg_dim))
1332 64 : ABI_MALLOC(prodr,(mdim,mdim))
1333 : !ABI_MALLOC(f3dfd,(2,nphi,deg_dim))
1334 60 : saddle_warn=.false.
1335 60 : start_eigf3d_pos=.true.
1336 :
1337 : !Hack to print f(theta,phi) & weights to a file
1338 : if(print_fsph) then
1339 : write(msg,*) degl+1
1340 : filename='f_band_'//TRIM(ADJUSTL(msg))//'-'
1341 : write(msg,*) degl+deg_dim
1342 : filename=TRIM(filename)//TRIM(ADJUSTL(msg))//'.dat'
1343 : io_unit = get_unit()
1344 : open(io_unit,file=TRIM(filename),status='replace')
1345 : write(io_unit,*) 'ntheta=',ntheta,', nphi=',nphi
1346 : write(io_unit,*) 'itheta, iphi, weight, f_n(theta,phi)'
1347 : end if
1348 :
1349 1616 : do itheta=1,ntheta
1350 1600 : costh=gq_points_costh(itheta) ; sinth=gq_points_sinth(itheta)
1351 321616 : do iphi=1,nphi
1352 320000 : cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
1353 320000 : weight=gq_weights_th(itheta)*gq_weights_ph(iphi)
1354 :
1355 320000 : unit_r(1)=sinth*cosph
1356 320000 : unit_r(2)=sinth*sinph
1357 320000 : unit_r(3)=costh
1358 :
1359 320000 : dr_dth(1)=costh*cosph
1360 320000 : dr_dth(2)=costh*sinph
1361 320000 : dr_dth(3)=-sinth
1362 :
1363 320000 : dr_dph(1)=-sinph !sin(theta)*
1364 320000 : dr_dph(2)=cosph !cos(theta)*
1365 320000 : dr_dph(3)=zero
1366 :
1367 1200000 : do iband=1,deg_dim
1368 3960000 : do jband=1,deg_dim
1369 63480000 : f3d(iband,jband)=DOT_PRODUCT(unit_r,MATMUL(eig2_diag_cart(:,:,iband,jband),unit_r))
1370 2760000 : df3d_dth(iband,jband)=DOT_PRODUCT(dr_dth,MATMUL(eig2_diag_cart(:,:,iband,jband),unit_r))+&
1371 126960000 : & DOT_PRODUCT(unit_r,MATMUL(eig2_diag_cart(:,:,iband,jband),dr_dth))
1372 2760000 : df3d_dph(iband,jband)=DOT_PRODUCT(dr_dph,MATMUL(eig2_diag_cart(:,:,iband,jband),unit_r))+&
1373 127840000 : & DOT_PRODUCT(unit_r,MATMUL(eig2_diag_cart(:,:,iband,jband),dr_dph))
1374 : end do
1375 : end do
1376 : !DIAGONALIZATION
1377 4280000 : eigenvec = f3d !IN
1378 320000 : lwork=-1
1379 320000 : ABI_MALLOC(work,(1))
1380 960000 : ABI_MALLOC(rwork,(3*deg_dim-2))
1381 320000 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1382 320000 : lwork=int(work(1))
1383 320000 : ABI_FREE(work)
1384 1200000 : eigenval = zero
1385 960000 : ABI_MALLOC(work,(lwork))
1386 3760000 : work=zero; rwork=zero
1387 320000 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1388 320000 : ABI_FREE(rwork)
1389 320000 : ABI_FREE(work)
1390 4280000 : unitary_tr = eigenvec !OUT
1391 1520000 : eigf3d = eigenval !OUT
1392 320060 : if(itheta==1 .and. iphi==1) start_eigf3d_pos = eigf3d > 0
1393 1200000 : do iband=1,deg_dim
1394 1200000 : if(start_eigf3d_pos(iband) .neqv. (eigf3d(iband)>0)) then
1395 188 : saddle_warn(iband)=.true.
1396 : end if
1397 : end do
1398 :
1399 : !Hack to print f(theta,phi)
1400 : if(print_fsph) write(io_unit,*) gq_points_th(itheta), gq_points_ph(iphi), weight, eigf3d(:)
1401 :
1402 : !!DEBUG-Mech.
1403 : !!A=-4.20449; B=0.378191; C=5.309 !Mech's fit
1404 : !A=-4.62503023; B=0.68699088; C=5.20516873 !My fit
1405 : !R = sqrt(B**2 + C**2*sin(theta)**2*(cos(theta)**2 + sin(theta)**2*sin(phi)**2*cos(phi)**2))
1406 : !eigf3d(1) = A - R
1407 : !eigf3d(2) = A + R
1408 :
1409 : !!angular FD
1410 : !f3dfd(2,iphi,:)=eigf3d(:)
1411 :
1412 1520000 : m_avg = m_avg + weight*sinth*eigf3d
1413 1520000 : m_avg_frohlich = m_avg_frohlich + weight*sinth/(abs(eigf3d))**half
1414 :
1415 320000 : prodc=MATMUL_(f3d,unitary_tr,deg_dim,deg_dim) ; f3d=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1416 : !f3d = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(f3d,unitary_tr))
1417 1200000 : do iband=1,deg_dim
1418 1200000 : eigf3d(iband) = real(f3d(iband,iband),dp)
1419 : end do
1420 320000 : prodc=MATMUL_(df3d_dth,unitary_tr,deg_dim,deg_dim) ; df3d_dth=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1421 : !df3d_dth=MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(df3d_dth,unitary_tr))
1422 1200000 : do iband=1,deg_dim
1423 1200000 : deigf3d_dth(iband) = real(df3d_dth(iband,iband),dp)
1424 : end do
1425 320000 : prodc=MATMUL_(df3d_dph,unitary_tr,deg_dim,deg_dim) ; df3d_dph=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1426 : !df3d_dph = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(df3d_dph,unitary_tr))
1427 1200000 : do iband=1,deg_dim
1428 1200000 : deigf3d_dph(iband) = real(df3d_dph(iband,iband),dp)
1429 : end do
1430 :
1431 : !!DEBUG-Mech.
1432 : !eigf3d(1) = A - R
1433 : !eigf3d(2) = A + R
1434 : !deigf3d_dth(1) = -1./2./R*C**2*(2.*sin(theta)*cos(theta)*(cos(theta)**2 + sin(theta)**2*sin(phi)**2*cos(phi)**2) + 2.*sin(theta)**3*cos(theta)*(sin(phi)**2*cos(phi)**2 - 1))
1435 : !deigf3d_dth(2) = 1./2./R*C**2*(2.*sin(theta)*cos(theta)*(cos(theta)**2 + sin(theta)**2*sin(phi)**2*cos(phi)**2) + 2.*sin(theta)**3*cos(theta)*(sin(phi)**2*cos(phi)**2 - 1))
1436 : !deigf3d_dph(1) = -1./2./R*C**2*sin(theta)**3*(2.*sin(phi)*cos(phi)**3 - 2.*sin(phi)**3*cos(phi))
1437 : !deigf3d_dph(2) = 1./2./R*C**2*sin(theta)**3*(2.*sin(phi)*cos(phi)**3 - 2.*sin(phi)**3*cos(phi))
1438 :
1439 : !!angular FD
1440 : !if(iphi/=1 .and. itheta/=1) then
1441 : ! deigf3d_dph(:) = (f3dfd(2,iphi,:)-f3dfd(2,iphi-1,:))/two_pi*nphi/sin(theta)
1442 : !else
1443 : ! deigf3d_dph(:) = zero
1444 : !end if
1445 : !if(itheta/=1) then
1446 : ! deigf3d_dth(:) = (f3dfd(2,iphi,:)-f3dfd(1,iphi,:))/pi*ntheta
1447 : !else
1448 : ! deigf3d_dth(:) = zero
1449 : !end if
1450 :
1451 1200000 : unit_speed(1,:) = 2._dp*sinth*cosph*eigf3d + costh*cosph*deigf3d_dth - sinph*deigf3d_dph !/sin(theta)
1452 1200000 : unit_speed(2,:) = 2._dp*sinth*sinph*eigf3d + costh*sinph*deigf3d_dth + cosph*deigf3d_dph !/sin(theta)
1453 1200000 : unit_speed(3,:) = 2._dp*costh*eigf3d - sinth*deigf3d_dth
1454 :
1455 1201600 : do jdeg=1,deg_dim
1456 3840000 : do bdir=1,mdim
1457 11440000 : do adir=1,mdim
1458 : transport_tensor(adir,bdir,jdeg) = transport_tensor(adir,bdir,jdeg) + &
1459 10560000 : & weight*sinth*unit_speed(adir,jdeg)*unit_speed(bdir,jdeg)/(ABS(eigf3d(jdeg))**2.5_dp)
1460 : end do
1461 : end do
1462 : end do
1463 : end do !iphi
1464 : !!angular FD
1465 : !f3dfd(1,:,:) = f3dfd(2,:,:)
1466 : end do !itheta
1467 :
1468 : !Hack to print f(theta,phi)
1469 : if(print_fsph) close(io_unit)
1470 :
1471 60 : m_avg = quarter/pi*m_avg
1472 60 : m_avg = one/m_avg
1473 :
1474 60 : m_avg_frohlich = quarter/pi*m_avg_frohlich
1475 60 : m_avg_frohlich = m_avg_frohlich**2
1476 :
1477 588 : transport_tensor = 1.0_dp/2.0_dp*transport_tensor
1478 :
1479 : !Effective masses along directions.
1480 88 : do adir=1,ndirs
1481 268 : do iband=1,deg_dim
1482 858 : do jband=1,deg_dim
1483 13766 : f3d(iband,jband) = dot_product(dirs(:,adir),matmul(eig2_diag_cart(:,:,iband,jband),dirs(:,adir)))
1484 : end do
1485 : end do
1486 : !f3d(:,:) = eig2_diag_cart(adir,adir,:,:)
1487 930 : eigenvec = f3d !IN
1488 72 : lwork=-1
1489 72 : ABI_MALLOC(work,(1))
1490 216 : ABI_MALLOC(rwork,(3*deg_dim-2))
1491 72 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1492 72 : lwork=int(work(1))
1493 72 : ABI_FREE(work)
1494 268 : eigenval = zero
1495 216 : ABI_MALLOC(work,(lwork))
1496 836 : work=zero; rwork=zero
1497 72 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1498 72 : ABI_FREE(rwork)
1499 72 : ABI_FREE(work)
1500 930 : unitary_tr = eigenvec !OUT
1501 340 : eigf3d = eigenval !OUT
1502 284 : m_cart(adir,:)=1._dp/eigf3d(:)
1503 : end do
1504 :
1505 60 : do iband=1,deg_dim
1506 : !DIAGONALIZATION
1507 572 : transport_eqv_eigvec(:,:,iband) = transport_tensor(:,:,iband)
1508 44 : lwork=-1
1509 44 : ABI_MALLOC(rwork,(1))
1510 44 : call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_tensor_eig,rwork,lwork,info)
1511 44 : lwork=int(rwork(1))
1512 44 : ABI_FREE(rwork)
1513 176 : transport_tensor_eig = zero
1514 132 : ABI_MALLOC(rwork,(lwork))
1515 440 : rwork=zero
1516 44 : call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_tensor_eig,rwork,lwork,info)
1517 44 : ABI_FREE(rwork)
1518 1100 : transport_eqv_eigvec(:,:,iband) = transpose(transport_eqv_eigvec(:,:,iband)) !So that lines contain eigenvectors.
1519 :
1520 44 : prodr=MATMUL_(transport_tensor(:,:,iband),transport_eqv_eigvec(:,:,iband),mdim,mdim,transb='t')
1521 44 : transport_tensor(:,:,iband)=MATMUL_(transport_eqv_eigvec(:,:,iband),prodr,mdim,mdim)
1522 : !transport_tensor(:,:,iband) = MATMUL(transport_eqv_eigvec(:,:,iband), &
1523 : ! MATMUL(transport_tensor(:,:,iband),TRANSPOSE(transport_eqv_eigvec(:,:,iband))))
1524 :
1525 44 : transport_eqv_eigval(1,iband) = transport_tensor_eig(2)*transport_tensor_eig(3)*(3._dp/8._dp/pi)**2
1526 44 : transport_eqv_eigval(2,iband) = transport_tensor_eig(3)*transport_tensor_eig(1)*(3._dp/8._dp/pi)**2
1527 44 : transport_eqv_eigval(3,iband) = transport_tensor_eig(1)*transport_tensor_eig(2)*(3._dp/8._dp/pi)**2
1528 : !The transport tensor loses the sign of the effective mass, this restores it.
1529 176 : transport_eqv_eigval(:,iband) = DSIGN(transport_eqv_eigval(:,iband),m_avg(iband))
1530 44 : transport_eqv_m(1,1,iband) = transport_eqv_eigval(1,iband)
1531 44 : transport_eqv_m(2,2,iband) = transport_eqv_eigval(2,iband)
1532 44 : transport_eqv_m(3,3,iband) = transport_eqv_eigval(3,iband)
1533 :
1534 44 : m_avg_frohlich(iband) = DSIGN(m_avg_frohlich(iband),m_avg(iband))
1535 :
1536 44 : prodr=MATMUL_(transport_eqv_m(:,:,iband),transport_eqv_eigvec(:,:,iband),mdim,mdim)
1537 60 : transport_eqv_m(:,:,iband)=MATMUL_(transport_eqv_eigvec(:,:,iband),prodr,mdim,mdim,transa='t')
1538 : !transport_eqv_m(:,:,iband) = MATMUL(TRANSPOSE(transport_eqv_eigvec(:,:,iband)), &
1539 : ! MATMUL(transport_eqv_m(:,:,iband),transport_eqv_eigvec(:,:,iband)))
1540 :
1541 : end do
1542 :
1543 : call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1544 16 : & ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec)
1545 : call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1546 16 : & ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec)
1547 :
1548 16 : ABI_FREE(unit_r)
1549 16 : ABI_FREE(dr_dth)
1550 16 : ABI_FREE(dr_dph)
1551 16 : ABI_FREE(f3d)
1552 16 : ABI_FREE(df3d_dth)
1553 16 : ABI_FREE(df3d_dph)
1554 16 : ABI_FREE(unitary_tr)
1555 16 : ABI_FREE(eigf3d)
1556 16 : ABI_FREE(saddle_warn)
1557 16 : ABI_FREE(start_eigf3d_pos)
1558 16 : ABI_FREE(m_avg)
1559 16 : ABI_FREE(m_avg_frohlich)
1560 16 : ABI_FREE(m_cart)
1561 16 : ABI_FREE(deigf3d_dth)
1562 16 : ABI_FREE(deigf3d_dph)
1563 16 : ABI_FREE(unit_speed)
1564 16 : ABI_FREE(transport_tensor)
1565 16 : ABI_FREE(transport_tensor_eig)
1566 16 : ABI_FREE(transport_eqv_m)
1567 16 : ABI_FREE(transport_eqv_eigval)
1568 16 : ABI_FREE(transport_eqv_eigvec)
1569 16 : ABI_FREE(prodc)
1570 16 : ABI_FREE(prodr)
1571 : !ABI_FREE(f3dfd)
1572 :
1573 2 : elseif (degenerate .and. mdim==2) then
1574 :
1575 5 : ABI_CALLOC(unit_r,(mdim))
1576 4 : ABI_CALLOC(dr_dph,(mdim))
1577 15 : ABI_CALLOC(f3d,(deg_dim,deg_dim))
1578 15 : ABI_CALLOC(df3d_dph,(deg_dim,deg_dim))
1579 15 : ABI_CALLOC(unitary_tr,(deg_dim,deg_dim))
1580 5 : ABI_CALLOC(eigf3d,(deg_dim))
1581 3 : ABI_MALLOC(saddle_warn,(deg_dim))
1582 2 : ABI_MALLOC(start_eigf3d_pos,(deg_dim))
1583 5 : ABI_CALLOC(m_avg,(deg_dim))
1584 5 : ABI_CALLOC(m_avg_frohlich,(deg_dim))
1585 13 : ABI_CALLOC(m_cart,(ndirs,deg_dim))
1586 5 : ABI_CALLOC(deigf3d_dph,(deg_dim))
1587 13 : ABI_CALLOC(unit_speed,(mdim,deg_dim))
1588 26 : ABI_CALLOC(transport_tensor,(mdim,mdim,deg_dim))
1589 10 : ABI_CALLOC(cart_rotation,(mdim,mdim))
1590 4 : ABI_CALLOC(transport_tensor_eig,(mdim))
1591 25 : ABI_CALLOC(transport_eqv_m,(mdim,mdim,deg_dim))
1592 12 : ABI_CALLOC(transport_eqv_eigval,(mdim,deg_dim))
1593 25 : ABI_CALLOC(transport_eqv_eigvec,(mdim,mdim,deg_dim))
1594 5 : ABI_CALLOC(transport_tensor_scale,(deg_dim))
1595 3 : ABI_MALLOC(prodc,(deg_dim,deg_dim))
1596 3 : ABI_MALLOC(prodr,(mdim,mdim))
1597 4 : saddle_warn=.false.
1598 4 : start_eigf3d_pos=.true.
1599 :
1600 201 : do iphi=1,nphi
1601 200 : cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
1602 200 : weight=gq_weights_ph(iphi)
1603 :
1604 200 : unit_r(1)=cosph
1605 200 : unit_r(2)=sinph
1606 :
1607 200 : dr_dph(1)=-sinph
1608 200 : dr_dph(2)=cosph
1609 :
1610 800 : do iband=1,deg_dim
1611 2600 : do jband=1,deg_dim
1612 12600 : matr2d = eig2_diag_cart(1:mdim,1:mdim,iband,jband)
1613 21600 : f3d(iband,jband)=DOT_PRODUCT(unit_r,MATMUL(matr2d,unit_r))
1614 : df3d_dph(iband,jband)=DOT_PRODUCT(dr_dph,MATMUL(matr2d,unit_r))+&
1615 42000 : & DOT_PRODUCT(unit_r,MATMUL(matr2d,dr_dph))
1616 : end do
1617 : end do
1618 :
1619 : !DIAGONALIZATION
1620 2800 : eigenvec = f3d !IN
1621 200 : lwork=-1
1622 200 : ABI_MALLOC(work,(1))
1623 600 : ABI_MALLOC(rwork,(3*deg_dim-2))
1624 200 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1625 200 : lwork=int(work(1))
1626 200 : ABI_FREE(work)
1627 800 : eigenval = zero
1628 600 : ABI_MALLOC(work,(lwork))
1629 2600 : work=zero; rwork=zero
1630 200 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1631 200 : ABI_FREE(rwork)
1632 200 : ABI_FREE(work)
1633 2800 : unitary_tr = eigenvec !OUT
1634 1000 : eigf3d = eigenval !OUT
1635 204 : if(iphi==1) start_eigf3d_pos = eigf3d > 0
1636 800 : do iband=1,deg_dim
1637 800 : if(start_eigf3d_pos(iband) .neqv. (eigf3d(iband)>0)) then
1638 0 : saddle_warn(iband)=.true.
1639 : end if
1640 : end do
1641 :
1642 1000 : m_avg = m_avg + weight*eigf3d
1643 1000 : m_avg_frohlich = m_avg_frohlich + weight/(abs(eigf3d))**half
1644 :
1645 200 : prodc=MATMUL_(f3d,unitary_tr,deg_dim,deg_dim) ; f3d=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1646 : !f3d = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(f3d,unitary_tr))
1647 800 : do iband=1,deg_dim
1648 800 : eigf3d(iband) = real(f3d(iband,iband),dp)
1649 : end do
1650 200 : prodc=MATMUL_(df3d_dph,unitary_tr,deg_dim,deg_dim) ; df3d_dph=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1651 : !df3d_dph = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(df3d_dph,unitary_tr))
1652 800 : do iband=1,deg_dim
1653 800 : deigf3d_dph(iband) = real(df3d_dph(iband,iband),dp)
1654 : end do
1655 :
1656 800 : unit_speed(1,:) = 2._dp*cosph*eigf3d - sinph*deigf3d_dph
1657 800 : unit_speed(2,:) = 2._dp*sinph*eigf3d + cosph*deigf3d_dph
1658 :
1659 801 : do jdeg=1,deg_dim
1660 2000 : do bdir=1,mdim
1661 4200 : do adir=1,mdim
1662 : transport_tensor(adir,bdir,jdeg) = transport_tensor(adir,bdir,jdeg) + &
1663 3600 : & weight*unit_speed(adir,jdeg)*unit_speed(bdir,jdeg)/(ABS(eigf3d(jdeg))**2)
1664 : end do
1665 : end do
1666 : end do
1667 :
1668 : end do !iphi
1669 :
1670 : !!!DEBUG
1671 :
1672 4 : m_avg = half/pi*m_avg
1673 4 : m_avg = one/m_avg
1674 :
1675 4 : m_avg_frohlich = half/pi*m_avg_frohlich
1676 4 : m_avg_frohlich = m_avg_frohlich**2
1677 :
1678 22 : transport_tensor = 1.0_dp/2.0_dp*transport_tensor
1679 :
1680 : !Effective masses along directions.
1681 3 : do adir=1,ndirs
1682 8 : do iband=1,deg_dim
1683 26 : do jband=1,deg_dim
1684 420 : f3d(iband,jband) = dot_product(dirs(:,adir),matmul(eig2_diag_cart(:,:,iband,jband),dirs(:,adir)))
1685 : end do
1686 : end do
1687 : !f3d(:,:) = eig2_diag_cart(adir,adir,:,:)
1688 28 : eigenvec = f3d !IN
1689 2 : lwork=-1
1690 2 : ABI_MALLOC(work,(1))
1691 6 : ABI_MALLOC(rwork,(3*deg_dim-2))
1692 2 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1693 2 : lwork=int(work(1))
1694 2 : ABI_FREE(work)
1695 8 : eigenval = zero
1696 6 : ABI_MALLOC(work,(lwork))
1697 26 : work=zero; rwork=zero
1698 2 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1699 2 : ABI_FREE(rwork)
1700 2 : ABI_FREE(work)
1701 28 : unitary_tr = eigenvec !OUT
1702 10 : eigf3d = eigenval !OUT
1703 9 : m_cart(adir,:)=1._dp/eigf3d(:)
1704 : end do
1705 :
1706 4 : do iband=1,deg_dim
1707 : !DIAGONALIZATION
1708 24 : cart_rotation = transport_tensor(:,:,iband)
1709 3 : lwork=-1
1710 3 : ABI_MALLOC(rwork,(1))
1711 3 : call dsyev('V','U',mdim,cart_rotation,mdim,transport_tensor_eig,rwork,lwork,info)
1712 3 : lwork=int(rwork(1))
1713 3 : ABI_FREE(rwork)
1714 9 : transport_tensor_eig = zero
1715 9 : ABI_MALLOC(rwork,(lwork))
1716 21 : rwork=zero
1717 3 : call dsyev('V','U',mdim,cart_rotation,mdim,transport_tensor_eig,rwork,lwork,info)
1718 3 : ABI_FREE(rwork)
1719 21 : transport_eqv_eigvec(:,:,iband) = transpose(cart_rotation(:,:)) !So that lines contain eigenvectors, not columns.
1720 :
1721 3 : prodr=MATMUL_(transport_tensor(:,:,iband),cart_rotation,mdim,mdim)
1722 3 : transport_tensor(:,:,iband)=MATMUL_(cart_rotation,prodr,mdim,mdim,transa='t')
1723 : !transport_tensor(:,:,iband) = MATMUL(TRANSPOSE(cart_rotation),MATMUL(transport_tensor(:,:,iband),cart_rotation))
1724 :
1725 3 : transport_eqv_eigval(1,iband) = 0.5*m_avg(iband)*(1.0 + transport_tensor_eig(2)/transport_tensor_eig(1))
1726 3 : transport_eqv_eigval(2,iband) = transport_eqv_eigval(1,iband)*transport_tensor_eig(1)/transport_tensor_eig(2)
1727 : !The transport tensor loses the sign of the effective mass, this restores it.
1728 9 : transport_eqv_eigval(:,iband) = SIGN(transport_eqv_eigval(:,iband),m_avg(iband))
1729 3 : transport_eqv_m(1,1,iband) = transport_eqv_eigval(1,iband)
1730 3 : transport_eqv_m(2,2,iband) = transport_eqv_eigval(2,iband)
1731 3 : transport_tensor_scale(iband) = sqrt(transport_tensor_eig(1)*transport_tensor_eig(2))/two_pi
1732 :
1733 3 : m_avg_frohlich(iband) = SIGN(m_avg_frohlich(iband),m_avg(iband))
1734 :
1735 3 : prodr=MATMUL_(transport_eqv_m(:,:,iband),cart_rotation,mdim,mdim,transb='t')
1736 4 : transport_eqv_m(:,:,iband)=MATMUL_(cart_rotation,prodr,mdim,mdim)
1737 : !transport_eqv_m(:,:,iband) = MATMUL(cart_rotation,MATMUL(transport_eqv_m(:,:,iband),TRANSPOSE(cart_rotation)))
1738 :
1739 : end do
1740 :
1741 : call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1742 1 : & ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec,transport_tensor_scale)
1743 : call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1744 1 : & ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec,transport_tensor_scale)
1745 :
1746 1 : ABI_FREE(unit_r)
1747 1 : ABI_FREE(dr_dph)
1748 1 : ABI_FREE(f3d)
1749 1 : ABI_FREE(df3d_dph)
1750 1 : ABI_FREE(unitary_tr)
1751 1 : ABI_FREE(eigf3d)
1752 1 : ABI_FREE(saddle_warn)
1753 1 : ABI_FREE(start_eigf3d_pos)
1754 1 : ABI_FREE(m_avg)
1755 1 : ABI_FREE(m_avg_frohlich)
1756 1 : ABI_FREE(m_cart)
1757 1 : ABI_FREE(deigf3d_dph)
1758 1 : ABI_FREE(unit_speed)
1759 1 : ABI_FREE(transport_tensor)
1760 1 : ABI_FREE(cart_rotation)
1761 1 : ABI_FREE(transport_tensor_eig)
1762 1 : ABI_FREE(transport_eqv_m)
1763 1 : ABI_FREE(transport_eqv_eigval)
1764 1 : ABI_FREE(transport_eqv_eigvec)
1765 1 : ABI_FREE(transport_tensor_scale)
1766 1 : ABI_FREE(prodc)
1767 1 : ABI_FREE(prodr)
1768 :
1769 25 : elseif (degenerate .and. mdim==1) then
1770 :
1771 15 : ABI_CALLOC(f3d,(deg_dim,deg_dim))
1772 15 : ABI_CALLOC(unitary_tr,(deg_dim,deg_dim))
1773 5 : ABI_CALLOC(eigf3d,(deg_dim))
1774 10 : ABI_CALLOC(m_cart,(ndirs,deg_dim))
1775 14 : ABI_CALLOC(transport_eqv_m,(mdim,mdim,deg_dim))
1776 5 : ABI_CALLOC(m_avg,(deg_dim))
1777 5 : ABI_CALLOC(m_avg_frohlich,(deg_dim))
1778 3 : ABI_MALLOC(saddle_warn,(deg_dim))
1779 :
1780 4 : saddle_warn=.false.
1781 :
1782 13 : f3d(:,:) = eig2_diag_cart(1,1,:,:)
1783 :
1784 : !DIAGONALIZATION
1785 14 : eigenvec = f3d !IN
1786 1 : lwork=-1
1787 1 : ABI_MALLOC(work,(1))
1788 3 : ABI_MALLOC(rwork,(3*deg_dim-2))
1789 1 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1790 1 : lwork=int(work(1))
1791 1 : ABI_FREE(work)
1792 4 : eigenval = zero
1793 3 : ABI_MALLOC(work,(lwork))
1794 13 : work=zero; rwork=zero
1795 1 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1796 1 : ABI_FREE(rwork)
1797 1 : ABI_FREE(work)
1798 14 : unitary_tr = eigenvec !OUT
1799 5 : eigf3d = eigenval !OUT
1800 :
1801 4 : transport_eqv_m(1,1,:)=1._dp/eigf3d(:)
1802 :
1803 : !Effective masses along directions.
1804 2 : do adir=1,ndirs
1805 4 : do iband=1,deg_dim
1806 13 : do jband=1,deg_dim
1807 210 : f3d(iband,jband) = dot_product(dirs(:,adir),matmul(eig2_diag_cart(:,:,iband,jband),dirs(:,adir)))
1808 : end do
1809 : end do
1810 14 : eigenvec = f3d !IN
1811 1 : lwork=-1
1812 1 : ABI_MALLOC(work,(1))
1813 3 : ABI_MALLOC(rwork,(3*deg_dim-2))
1814 1 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1815 1 : lwork=int(work(1))
1816 1 : ABI_FREE(work)
1817 4 : eigenval = zero
1818 3 : ABI_MALLOC(work,(lwork))
1819 13 : work=zero; rwork=zero
1820 1 : call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1821 1 : ABI_FREE(rwork)
1822 1 : ABI_FREE(work)
1823 5 : eigf3d = eigenval !OUT
1824 5 : m_cart(adir,:)=1._dp/eigf3d(:)
1825 : end do
1826 :
1827 : call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m,&
1828 1 : & ntheta,m_avg,m_avg_frohlich,saddle_warn)
1829 : call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m,&
1830 1 : & ntheta,m_avg,m_avg_frohlich,saddle_warn)
1831 :
1832 1 : ABI_FREE(f3d)
1833 1 : ABI_FREE(unitary_tr)
1834 1 : ABI_FREE(eigf3d)
1835 1 : ABI_FREE(m_cart)
1836 1 : ABI_FREE(transport_eqv_m)
1837 1 : ABI_FREE(m_avg)
1838 1 : ABI_FREE(m_avg_frohlich)
1839 1 : ABI_FREE(saddle_warn)
1840 : end if !(degenerate)
1841 :
1842 42 : ABI_FREE(eig2_diag_cart)
1843 42 : ABI_FREE(eigenval)
1844 66 : ABI_FREE(eigenvec)
1845 : end do !ideg
1846 :
1847 : end do !ikpt
1848 :
1849 17 : ABI_FREE(eff_mass)
1850 17 : ABI_FREE(dirs)
1851 17 : ABI_FREE(gq_points_th)
1852 17 : ABI_FREE(gq_points_costh)
1853 17 : ABI_FREE(gq_points_sinth)
1854 17 : ABI_FREE(gq_weights_th)
1855 17 : ABI_FREE(gq_points_ph)
1856 17 : ABI_FREE(gq_points_cosph)
1857 17 : ABI_FREE(gq_points_sinph)
1858 17 : ABI_FREE(gq_weights_ph)
1859 :
1860 17 : write(std_out,'(3a)') ch10,' END OF EFFECTIVE MASSES SECTION',ch10
1861 17 : write(ab_out, '(3a)') ch10,' END OF EFFECTIVE MASSES SECTION',ch10
1862 :
1863 17 : end subroutine efmas_analysis
1864 : !!***
1865 :
1866 : !----------------------------------------------------------------------
1867 :
1868 : !!****f* m_efmas/MATMUL_DP
1869 : !! NAME
1870 : !! MATMUL_DP
1871 : !!
1872 : !! FUNCTION
1873 : !! Mimic MATMUL Fortran intrinsic function with BLAS3: C=A.B
1874 : !! This is a temporary workaround to make tests pass on intel/mkl architectures
1875 : !! Real version
1876 : !!
1877 : !! INPUTS
1878 : !! aa(:,:),bb(:,:)= input matrices
1879 : !! mm,nn= sizes of output matrix
1880 : !! [transa,transb]= equivalent to transa, transb args of gemm ('n','t','c')
1881 : !! if not present, default is 'n'.
1882 : !!
1883 : !! OUTPUT
1884 : !! MATMUL_DP(mm,nn)= output matrix A.B
1885 : !!
1886 : !! SOURCE
1887 :
1888 188 : function MATMUL_DP(aa,bb,mm,nn,transa,transb)
1889 :
1890 : !Arguments ------------------------------------
1891 : !scalars
1892 : integer,intent(in) :: mm,nn
1893 : character(len=1),optional,intent(in) :: transa,transb
1894 : !arrays
1895 : real(dp),intent(in) :: aa(:,:),bb(:,:)
1896 : real(dp) :: MATMUL_DP(mm,nn)
1897 :
1898 : !Local variables-------------------------------
1899 : integer :: kk,lda,ldb
1900 : character(len=1) :: transa_,transb_
1901 :
1902 : ! *************************************************************************
1903 :
1904 188 : transa_='n';if (present(transa)) transa_=transa
1905 188 : transb_='n';if (present(transb)) transb_=transb
1906 :
1907 188 : lda=size(aa,1) ; ldb=size(bb,1)
1908 :
1909 188 : if (transa_=='n') then
1910 141 : kk=size(aa,2)
1911 141 : if (size(aa,1)/=mm) then
1912 0 : ABI_BUG('Error in sizes!')
1913 : end if
1914 : else
1915 47 : kk=size(aa,1)
1916 47 : if (size(aa,2)/=mm) then
1917 0 : ABI_BUG('Error in sizes!')
1918 : end if
1919 : end if
1920 :
1921 188 : if (transb_=='n') then
1922 141 : if (size(bb,1)/=kk.or.size(bb,2)/=nn) then
1923 0 : ABI_BUG('Error in sizes!')
1924 : end if
1925 : else
1926 47 : if (size(bb,1)/=nn.or.size(bb,2)/=kk) then
1927 0 : ABI_BUG('Error in sizes!')
1928 : end if
1929 : end if
1930 :
1931 188 : call DGEMM(transa_,transb_,mm,nn,kk,one,aa,lda,bb,ldb,zero,MATMUL_DP,mm)
1932 :
1933 : end function MATMUL_DP
1934 : !!***
1935 :
1936 : !----------------------------------------------------------------------
1937 :
1938 : !!****f* m_efmas/MATMUL_DPC
1939 : !! NAME
1940 : !! MATMUL_DPC
1941 : !!
1942 : !! FUNCTION
1943 : !! Mimic MATMUL Fortran intrinsic function with BLAS3
1944 : !! This is a temporary workaround to make tests pass on intel/mkl architectures
1945 : !! Complex version
1946 : !!
1947 : !! INPUTS
1948 : !! aa(:,:),bb(:,:)= input matrices
1949 : !! mm,nn= sizes of output matrix
1950 : !! [transa,transb]= equivalent to transa, transb args of gemm ('n','t','c')
1951 : !! if not present, default is 'n'.
1952 : !!
1953 : !! OUTPUT
1954 : !! MATMUL_DPC(:,:)= output matrix A.B
1955 : !!
1956 : !! SOURCE
1957 :
1958 1920800 : function MATMUL_DPC(aa,bb,mm,nn,transa,transb)
1959 :
1960 : !Arguments ------------------------------------
1961 : !scalars
1962 : integer,intent(in) :: mm,nn
1963 : character(len=1),optional,intent(in) :: transa,transb
1964 : !arrays
1965 : complex(dp),intent(in) :: aa(:,:),bb(:,:)
1966 : complex(dp) :: MATMUL_DPC(mm,nn)
1967 :
1968 : !Local variables-------------------------------
1969 : integer :: kk,lda,ldb
1970 : character(len=1) :: transa_,transb_
1971 : ! *************************************************************************
1972 :
1973 1920800 : transa_='n';if (present(transa)) transa_=transa
1974 1920800 : transb_='n';if (present(transb)) transb_=transb
1975 :
1976 1920800 : lda=size(aa,1) ; ldb=size(bb,1)
1977 :
1978 1920800 : if (transa_=='n') then
1979 960400 : kk=size(aa,2)
1980 960400 : if (size(aa,1)/=mm) then
1981 0 : ABI_BUG('Error in sizes!')
1982 : end if
1983 : else
1984 960400 : kk=size(aa,1)
1985 960400 : if (size(aa,2)/=mm) then
1986 0 : ABI_BUG('Error in sizes!')
1987 : end if
1988 : end if
1989 :
1990 1920800 : if (transb_=='n') then
1991 1920800 : if (size(bb,1)/=kk.or.size(bb,2)/=nn) then
1992 0 : ABI_BUG('Error in sizes!')
1993 : end if
1994 : else
1995 0 : if (size(bb,1)/=nn.or.size(bb,2)/=kk) then
1996 0 : ABI_BUG('Error in sizes!')
1997 : end if
1998 : end if
1999 :
2000 1920800 : call ZGEMM(transa_,transb_,mm,nn,kk,cone,aa,lda,bb,ldb,czero,MATMUL_DPC,mm)
2001 :
2002 : end function MATMUL_DPC
2003 : !!***
2004 :
2005 14366357 : end module m_efmas
2006 : !!***
|
