Branch data 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 : 723 : 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 [ + + ]: 723 : 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 : 723 : 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 : 723 : 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 [ + + ]: 723 : 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 : 723 : 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 : : !!***
|
