TABLE OF CONTENTS


ABINIT/efmasval_free_array [ Functions ]

[ Top ] [ Functions ]

NAME

 efmasval_free_array

FUNCTION

 This routine deallocates an efmasval_type or, optionally, an array of efmasval_type.

INPUTS

OUTPUT

PARENTS

      respfn

CHILDREN

SOURCE

126  subroutine efmasval_free_array(efmasval)
127 
128 
129 !This section has been created automatically by the script Abilint (TD).
130 !Do not modify the following lines by hand.
131 #undef ABI_FUNC
132 #define ABI_FUNC 'efmasval_free_array'
133 !End of the abilint section
134 
135    implicit none
136 
137    !Arguments ------------------------------------
138    type(efmasval_type),allocatable,intent(inout) :: efmasval(:,:)
139 
140    !!!Local variables-------------------------------
141    integer :: i,j,n(2)
142 
143    ! *********************************************************************
144 
145    !XG20180810: please do not remove. Otherwise, I get an error on my Mac.
146    write(std_out,*)' efmasval_free_array : enter '
147 
148    if(allocated(efmasval)) then
149      n=shape(efmasval)
150      do i=1,n(1)
151        do j=1,n(2)
152          call efmasval_free(efmasval(i,j))
153        end do
154      end do
155      ABI_DATATYPE_DEALLOCATE(efmasval)
156    end if
157 
158  end subroutine efmasval_free_array

ABINIT/m_efmas [ Modules ]

[ Top ] [ Modules ]

NAME

 m_efmas

FUNCTION

 This module contains datatypes for efmas functionalities.

COPYRIGHT

 Copyright (C) 2001-2018 ABINIT group (JLJ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_efmas
24 
25  use defs_basis
26  use defs_abitypes
27  use m_errors
28 #ifdef HAVE_NETCDF
29  use netcdf, only : nf90_noerr, nf90_put_var, nf90_get_var
30 #endif
31  use m_efmas_defs
32  use m_nctk
33 
34  implicit none
35 
36  private
37 
38 !public procedures.
39  public :: efmasval_free
40  public :: efmasval_free_array
41  public :: efmasdeg_free
42  public :: efmasdeg_free_array
43  public :: efmas_ncread
44  public :: check_degeneracies
45  public :: print_tr_efmas
46  public :: print_efmas
47  public :: efmas_main
48  public :: efmas_analysis 
49 
50 !private procedures.
51  private :: MATMUL_ ! Workaround to make tests pass on ubu/buda slaves
52  interface MATMUL_
53   module procedure MATMUL_DP
54   module procedure MATMUL_DPC
55  end interface MATMUL_

m_efmas/check_degeneracies [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 check_degeneracies

FUNCTION

 This routine check for 0th order band degeneracies at given k-point.

INPUTS

OUTPUT

PARENTS

      d2frnl

CHILDREN

SOURCE

273  subroutine check_degeneracies(efmasdeg,bands,nband,eigen,deg_tol)
274 
275 
276 !This section has been created automatically by the script Abilint (TD).
277 !Do not modify the following lines by hand.
278 #undef ABI_FUNC
279 #define ABI_FUNC 'check_degeneracies'
280 !End of the abilint section
281 
282    implicit none
283 
284    !Arguments ------------------------------------
285    type(efmasdeg_type),intent(out) :: efmasdeg
286    integer,intent(in) :: bands(2),nband
287    real(dp),intent(in) :: eigen(nband)
288    real(dp),intent(in),optional :: deg_tol
289 
290    !!!Local variables-------------------------------
291    integer :: deg_dim,iband, ideg
292    integer, allocatable :: degs_bounds(:,:)
293    real(dp) :: tol
294    real(dp) :: eigen_tmp(nband)
295    logical :: treated
296 
297    ! *********************************************************************
298 
299    tol=tol5; if(present(deg_tol)) tol=deg_tol
300 
301    !!! Determine sets of degenerate states in eigen0, i.e., at 0th order.
302    efmasdeg%ndegs=1
303    efmasdeg%nband=nband
304    ABI_MALLOC(degs_bounds,(2,nband))
305    ABI_MALLOC(efmasdeg%ideg, (nband))
306    degs_bounds=0; degs_bounds(1,1)=1
307    efmasdeg%ideg=0; efmasdeg%ideg(1)=1
308 
309    eigen_tmp(:) = eigen(:)
310 
311    do iband=2,nband
312      if (ABS(eigen_tmp(iband)-eigen_tmp(iband-1))>tol) then
313        degs_bounds(2,efmasdeg%ndegs) = iband-1
314        efmasdeg%ndegs=efmasdeg%ndegs+1
315        degs_bounds(1,efmasdeg%ndegs) = iband
316      end if
317      efmasdeg%ideg(iband) = efmasdeg%ndegs
318    end do
319    degs_bounds(2,efmasdeg%ndegs)=nband
320    ABI_MALLOC(efmasdeg%degs_bounds,(2,efmasdeg%ndegs))
321    efmasdeg%degs_bounds(1:2,1:efmasdeg%ndegs) = degs_bounds(1:2,1:efmasdeg%ndegs)
322    ABI_FREE(degs_bounds)
323 
324    !!! Determine if treated bands are part of a degeneracy at 0th order.
325    efmasdeg%deg_range=0
326    deg_dim=0
327    treated=.false.
328    write(std_out,'(a,i6)') 'Number of sets of bands for this k-point:',efmasdeg%ndegs
329    write(std_out,'(a)') 'Set index; range of bands included in the set; is the set degenerate?(T/F); &
330 &                        is the set treated by EFMAS?(T/F):'
331    do ideg=1,efmasdeg%ndegs
332      deg_dim = efmasdeg%degs_bounds(2,ideg) - efmasdeg%degs_bounds(1,ideg) + 1
333      !If there is some level in the set that is inside the interval defined by bands(1:2), treat such set
334      !The band range might be larger than the nband interval: it includes it, and also include degenerate states
335      if(efmasdeg%degs_bounds(1,ideg)<=bands(2) .and. efmasdeg%degs_bounds(2,ideg)>=bands(1)) then
336        treated = .true.
337        if(efmasdeg%degs_bounds(1,ideg)<=bands(1)) then
338          efmasdeg%deg_range(1) = ideg
339        end if
340        if(efmasdeg%degs_bounds(2,ideg)>=bands(2)) then
341          efmasdeg%deg_range(2) = ideg
342        end if
343      end if
344      write(std_out,'(2i6,a,i6,2l4)') ideg, efmasdeg%degs_bounds(1,ideg), ' -', efmasdeg%degs_bounds(2,ideg), &
345 &                                    (deg_dim>1), treated
346    end do
347 
348 !   write(std_out,*)'ndegs=',          efmasdeg%ndegs
349 !   write(std_out,*)'degs_bounds=',    efmasdeg%degs_bounds
350 !   write(std_out,*)'ideg=',           efmasdeg%ideg
351 !   write(std_out,*)'deg_range=',      efmasdeg%deg_range
352 
353   !!This first attempt WORKS, but only if the symmetries are enabled, see line 1578 of dfpt_looppert.F90.
354   !use m_crystal,          only : crystal_t, crystal_init, crystal_free, crystal_print
355   !use m_esymm
356   !integer :: timrev
357   !character(len=132),allocatable :: title(:)
358   !type(crystal_t) :: Cryst
359   !type(esymm_t) :: Bsym
360 
361   !timrev = 1
362   !if(dtset%istwfk(1)/=1) timrev=2
363   !ABI_ALLOCATE(title,(dtset%ntypat))
364   !title(:) = "Bloup"
365   !call crystal_init(Cryst,dtset%spgroup,dtset%natom,dtset%npsp,dtset%ntypat,dtset%nsym,dtset%rprimd_orig(:,:,1),&
366   !&                 dtset%typat,dtset%xred_orig(:,:,1),dtset%ziontypat,dtset%znucl,timrev,.false.,.false.,title,&
367   !&                 dtset%symrel,dtset%tnons,dtset%symafm)
368   !call crystal_print(Cryst)
369   !ABI_DEALLOCATE(title)
370   !call esymm_init(Bsym,kpt_rbz(:,ikpt),Cryst,.false.,nspinor,1,mband,tol5,eigen0,dtset%tolsym)
371   !write(std_out,*) 'DEBUG : Bsym. ndegs=',Bsym%ndegs
372   !do iband=1,Bsym%ndegs
373   !  write(std_out,*) Bsym%degs_bounds(:,iband)
374   !end do
375 
376   !call crystal_free(Cryst)
377   !call esymm_free(Bsym)
378 
379  end subroutine check_degeneracies

m_efmas/efmas_analysis [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 efmas_analysis

FUNCTION

 This routine analyzes the generalized second-order k-derivatives of eigenenergies, 
 and compute the effective mass tensor
 (inverse of hessian of eigenvalues with respect to the wavevector)
 in cartesian coordinates along different directions in k-space, or also the transport equivalent effective mass.

INPUTS

  dtset = dataset structure containing the input variable of the calculation.
  efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point
  efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure
    efmasval(:,:)%eig2_diag band curvature double tensor
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
  mpert = maximum number of perturbations.
  mpi_enreg = informations about MPI parallelization.
  nkpt_rbz = number of k-points for each perturbation.
  rprimd(3,3)=dimensional primitive translations for real space (bohr)

OUTPUT

PARENTS

      dfpt_looppert

CHILDREN

      cgqf,dgemm,dgetrf,dgetri,dotprod_g,dsyev,print_tr_efmas,zgemm,zgetrf
      zgetri,zheev

SOURCE

1232  subroutine efmas_analysis(dtset,efmasdeg,efmasval,kpt_rbz,mpi_enreg,nkpt_rbz,rprimd)
1233 
1234   use m_cgtools
1235   use m_gaussian_quadrature, only : cgqf
1236   use m_io_tools, only : get_unit
1237 
1238 !This section has been created automatically by the script Abilint (TD).
1239 !Do not modify the following lines by hand.
1240 #undef ABI_FUNC
1241 #define ABI_FUNC 'efmas_analysis'
1242 !End of the abilint section
1243 
1244   implicit none
1245 
1246  !Arguments ------------------------------------
1247  !scalars
1248   integer,            intent(in)    :: nkpt_rbz
1249   type(dataset_type), intent(in)    :: dtset
1250   type(MPI_type),     intent(in) :: mpi_enreg
1251  !arrays
1252   real(dp), intent(in) :: rprimd(3,3)
1253   real(dp), intent(in) :: kpt_rbz(3,nkpt_rbz)
1254   type(efmasdeg_type), intent(in) :: efmasdeg(:)
1255   type(efmasval_type), intent(in) :: efmasval(:,:)
1256 
1257  !Local variables-------------------------------
1258   logical :: degenerate
1259   logical :: debug
1260   logical :: print_fsph
1261   logical, allocatable :: saddle_warn(:), start_eigf3d_pos(:)
1262   integer :: info
1263   integer :: ipert
1264   integer :: isppol
1265   integer :: nband_k
1266   integer :: ideg,jdeg
1267   integer :: ikpt
1268   integer :: master,me,spaceworld
1269   integer :: iband, jband
1270   integer :: adir,bdir
1271   integer :: deg_dim
1272   integer :: degl
1273   integer :: lwork
1274   integer :: itheta, iphi
1275   integer :: ntheta, nphi
1276   integer :: mdim
1277   integer :: cdirs, ndirs
1278   integer :: ipiv(3)
1279   integer :: io_unit
1280   character(len=500) :: msg, filename
1281   real(dp) :: cosph,costh,sinph,sinth
1282   real(dp) :: f3d_scal
1283   real(dp) :: weight
1284   real(dp) :: gprimd(3,3)
1285   real(dp), allocatable :: unit_r(:), dr_dth(:), dr_dph(:)
1286   real(dp), allocatable :: eigenval(:), rwork(:)
1287   real(dp), allocatable :: eigf3d(:)
1288   real(dp), allocatable :: m_avg(:), m_avg_frohlich(:),m_cart(:,:)
1289   real(dp), allocatable :: deigf3d_dth(:), deigf3d_dph(:)
1290   real(dp), allocatable :: unit_speed(:,:), transport_tensor(:,:,:)
1291   real(dp), allocatable :: cart_rotation(:,:), transport_tensor_eig(:)
1292   real(dp), allocatable :: transport_eqv_m(:,:,:), transport_eqv_eigval(:,:), transport_eqv_eigvec(:,:,:)
1293   real(dp), allocatable :: transport_tensor_scale(:)
1294   real(dp), allocatable :: gq_points_th(:),gq_points_costh(:),gq_points_sinth(:),gq_weights_th(:)
1295   real(dp), allocatable :: gq_points_ph(:),gq_points_cosph(:),gq_points_sinph(:),gq_weights_ph(:)
1296   real(dp), allocatable :: dirs(:,:)
1297   real(dp),allocatable :: prodr(:,:)
1298   !real(dp), allocatable :: f3dfd(:,:,:)
1299   complex(dpc) :: matr2d(2,2)
1300   complex(dpc), allocatable :: eigen1_deg(:,:), eigenvec(:,:), work(:)
1301   complex(dpc), allocatable :: eig2_diag_cart(:,:,:,:)
1302   complex(dpc), allocatable :: f3d(:,:), df3d_dth(:,:), df3d_dph(:,:)
1303   complex(dpc), allocatable :: unitary_tr(:,:), eff_mass(:,:)
1304   complex(dpc),allocatable :: prodc(:,:)
1305 
1306 ! *********************************************************************
1307 
1308   debug = .false. ! Prints additional info to std_out
1309   print_fsph = .false. ! Open a file and print the angle dependent curvature f(\theta,\phi)
1310                        ! for each band & kpts treated; 1 file per degenerate ensemble of bands.
1311                        ! Angles are those used in the numerical integration.
1312 
1313 ! Init parallelism
1314   master =0
1315   spaceworld=mpi_enreg%comm_cell
1316   me=mpi_enreg%me_kpt
1317 
1318   isppol = 1
1319 
1320   mdim = dtset%efmas_dim
1321 
1322 !HERE ALLOCATE
1323 
1324   gprimd = rprimd
1325   call dgetrf(mdim,mdim,gprimd,mdim,ipiv,info)
1326   ABI_ALLOCATE(rwork,(3))
1327   call dgetri(mdim,gprimd,mdim,ipiv,rwork,3,info)
1328   ABI_DEALLOCATE(rwork)
1329   gprimd = two_pi*transpose(gprimd)
1330 
1331   cdirs = dtset%efmas_calc_dirs
1332   ndirs = mdim
1333   if(cdirs/=0) ndirs = dtset%efmas_n_dirs
1334   ABI_ALLOCATE(dirs,(3,ndirs))
1335   if(cdirs==0) then
1336     dirs = zero
1337     do adir=1,ndirs
1338       dirs(adir,adir)=1.0_dp
1339     end do
1340   elseif(cdirs==1) then
1341     dirs(:,:) = dtset%efmas_dirs(:,1:ndirs)
1342     do adir=1,ndirs
1343       dirs(:,adir) = dirs(:,adir)/sqrt(sum(dirs(:,adir)**2))
1344     end do
1345   elseif(cdirs==2) then
1346     dirs(:,:) = matmul(gprimd,dtset%efmas_dirs(:,1:ndirs))
1347     do adir=1,ndirs
1348       dirs(:,adir) = dirs(:,adir)/sqrt(sum(dirs(:,adir)**2))
1349     end do
1350   elseif(cdirs==3) then
1351     dirs(1,:) = sin(dtset%efmas_dirs(1,1:ndirs)*pi/180)*cos(dtset%efmas_dirs(2,1:ndirs)*pi/180)
1352     dirs(2,:) = sin(dtset%efmas_dirs(1,1:ndirs)*pi/180)*sin(dtset%efmas_dirs(2,1:ndirs)*pi/180)
1353     dirs(3,:) = cos(dtset%efmas_dirs(1,1:ndirs)*pi/180)
1354   end if
1355 
1356   !!! Initialization of integrals for the degenerate case.
1357   ntheta   = dtset%efmas_ntheta
1358   nphi     = 2*ntheta
1359   ABI_ALLOCATE(gq_points_th,(ntheta))
1360   ABI_ALLOCATE(gq_points_costh,(ntheta))
1361   ABI_ALLOCATE(gq_points_sinth,(ntheta))
1362   ABI_ALLOCATE(gq_weights_th,(ntheta))
1363   ABI_ALLOCATE(gq_points_ph,(nphi))
1364   ABI_ALLOCATE(gq_points_cosph,(nphi))
1365   ABI_ALLOCATE(gq_points_sinph,(nphi))
1366   ABI_ALLOCATE(gq_weights_ph,(nphi))
1367   call cgqf(ntheta,1,zero,zero,zero,pi,gq_points_th,gq_weights_th)
1368   !XG180501 : TODO : There is no need to make a Gauss-Legendre integral for the phi variable,
1369   !since the function to be integrated is periodic...
1370   call cgqf(nphi,1,zero,zero,zero,2*pi,gq_points_ph,gq_weights_ph)
1371   do itheta=1,ntheta
1372     gq_points_costh(itheta)=cos(gq_points_th(itheta))
1373     gq_points_sinth(itheta)=sin(gq_points_th(itheta))
1374   enddo
1375   do iphi=1,nphi
1376     gq_points_cosph(iphi)=cos(gq_points_ph(iphi))
1377     gq_points_sinph(iphi)=sin(gq_points_ph(iphi))
1378   enddo
1379 
1380   ABI_ALLOCATE(eff_mass,(mdim,mdim))
1381 
1382 !XG20180519 : incoherent, efmasdeg is dimensioned at nkpt_rbz, and not at dtset%nkpt ...
1383   do ikpt=1,dtset%nkpt
1384     do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
1385 
1386      deg_dim    = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
1387      degenerate = (deg_dim>1) .and. (dtset%efmas_deg/=0)
1388      degl       = efmasdeg(ikpt)%degs_bounds(1,ideg)-1
1389  
1390      !!! Allocations
1391      ABI_ALLOCATE(eigenvec,(deg_dim,deg_dim))
1392      ABI_ALLOCATE(eigenval,(deg_dim))
1393 
1394      ABI_ALLOCATE(eig2_diag_cart,(3,3,deg_dim,deg_dim))
1395 
1396      do iband=1,deg_dim
1397         write(std_out,*)"  Compute band ",iband  ! This line here to avoid weird
1398         do jband=1,deg_dim
1399 
1400           eff_mass=zero
1401 
1402           eig2_diag_cart(:,:,iband,jband)=efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)          
1403           eig2_diag_cart(:,:,iband,jband) = matmul(matmul(rprimd,eig2_diag_cart(:,:,iband,jband)),transpose(rprimd))/two_pi**2
1404 
1405 
1406           if(.not. degenerate .and. iband==jband) then
1407 
1408             !Compute effective mass tensor from second derivative matrix. Simple inversion.
1409             eff_mass(:,:) = eig2_diag_cart(1:mdim,1:mdim,iband,jband)
1410             call zgetrf(mdim,mdim,eff_mass(1:mdim,1:mdim),mdim,ipiv,info)
1411             ABI_ALLOCATE(work,(3))
1412             call zgetri(mdim,eff_mass(1:mdim,1:mdim),mdim,ipiv,work,3,info)
1413             ABI_DEALLOCATE(work)
1414 
1415             !DIAGONALIZATION
1416             ABI_ALLOCATE(transport_eqv_eigvec,(mdim,mdim,deg_dim))
1417             transport_eqv_eigvec=zero
1418             ABI_ALLOCATE(transport_eqv_eigval,(mdim,deg_dim))
1419             transport_eqv_eigval=zero
1420             transport_eqv_eigvec(:,:,iband) = real(eff_mass(1:mdim,1:mdim),dp)
1421             lwork=-1
1422             ABI_ALLOCATE(rwork,(1))
1423             call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_eqv_eigval(:,iband),rwork,lwork,info)
1424             lwork = max(1, 3*mdim-1) ! lwork >= max(1, 3*mdim-1)
1425             ABI_DEALLOCATE(rwork)
1426 
1427             ABI_ALLOCATE(rwork,(lwork))
1428             rwork=zero
1429             call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_eqv_eigval(:,iband),rwork,lwork,info)
1430             ABI_DEALLOCATE(rwork)
1431             transport_eqv_eigvec(:,:,iband) = transpose(transport_eqv_eigvec(:,:,iband)) !So that lines contain eigenvectors.
1432 
1433             !Frohlich average effective mass
1434             ABI_ALLOCATE(m_avg,(1))
1435             ABI_ALLOCATE(m_avg_frohlich,(1))
1436             ABI_ALLOCATE(saddle_warn,(1))
1437             ABI_ALLOCATE(unit_r,(mdim))
1438             ABI_ALLOCATE(start_eigf3d_pos,(1))
1439 
1440             m_avg=zero
1441             m_avg_frohlich=zero
1442             saddle_warn=.false.
1443 
1444             if(mdim==3)then
1445               !One has to perform the integral over the sphere
1446               do itheta=1,ntheta
1447                 costh=gq_points_costh(itheta) ; sinth=gq_points_sinth(itheta)
1448                 do iphi=1,nphi
1449                   cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
1450                   weight=gq_weights_th(itheta)*gq_weights_ph(iphi)
1451 
1452                   unit_r(1)=sinth*cosph
1453                   unit_r(2)=sinth*sinph
1454                   unit_r(3)=costh
1455 
1456                   f3d_scal=dot_product(unit_r(:),matmul(real(eig2_diag_cart(:,:,iband,jband),dp),unit_r(:))) 
1457                   m_avg = m_avg + weight*sinth*f3d_scal
1458                   m_avg_frohlich = m_avg_frohlich + weight*sinth/(abs(f3d_scal)**half)
1459 
1460                   if(itheta==1 .and. iphi==1) start_eigf3d_pos = f3d_scal > 0
1461                   if(start_eigf3d_pos(1) .neqv. (f3d_scal>0)) then
1462                     saddle_warn(1)=.true.
1463                   end if
1464                 enddo
1465               enddo
1466               m_avg = quarter/pi*m_avg
1467               m_avg = one/m_avg
1468               m_avg_frohlich = quarter/pi*m_avg_frohlich
1469               m_avg_frohlich = m_avg_frohlich**2
1470               m_avg_frohlich(1) = DSIGN(m_avg_frohlich(1),m_avg(1))
1471 
1472             endif ! mdim==3
1473 
1474             !EFMAS_DIRS
1475             ABI_ALLOCATE(m_cart,(ndirs,deg_dim))
1476             m_cart=zero
1477             do adir=1,ndirs
1478               m_cart(adir,1)=1.0_dp/dot_product(dirs(:,adir),matmul(real(eig2_diag_cart(:,:,iband,jband),dp),dirs(:,adir)))
1479             end do
1480 
1481             !PRINTING RESULTS
1482             call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+iband,1,mdim,ndirs,dirs,m_cart,rprimd,real(eff_mass,dp), &
1483 &             ntheta,m_avg,m_avg_frohlich,saddle_warn,&
1484 &             transport_eqv_eigval(:,iband:iband),transport_eqv_eigvec(:,:,iband:iband))
1485             call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+iband,1,mdim,ndirs,dirs,m_cart,rprimd,real(eff_mass,dp), &
1486 &             ntheta,m_avg,m_avg_frohlich,saddle_warn,&
1487 &             transport_eqv_eigval(:,iband:iband),transport_eqv_eigvec(:,:,iband:iband))
1488             ABI_DEALLOCATE(m_cart)
1489             ABI_DEALLOCATE(transport_eqv_eigvec)
1490             ABI_DEALLOCATE(transport_eqv_eigval)
1491             ABI_DEALLOCATE(m_avg)
1492             ABI_DEALLOCATE(m_avg_frohlich)
1493             ABI_DEALLOCATE(unit_r)
1494             ABI_DEALLOCATE(saddle_warn)
1495             ABI_DEALLOCATE(start_eigf3d_pos)
1496 
1497           end if !.not.degenerate
1498         end do !jband
1499       end do !iband
1500 
1501       !!! EQV_MASS
1502       if(degenerate .and. mdim==3) then
1503         ABI_ALLOCATE(unit_r,(mdim))
1504         ABI_ALLOCATE(dr_dth,(mdim))
1505         ABI_ALLOCATE(dr_dph,(mdim))
1506         ABI_ALLOCATE(f3d,(deg_dim,deg_dim))
1507         ABI_ALLOCATE(df3d_dth,(deg_dim,deg_dim))
1508         ABI_ALLOCATE(df3d_dph,(deg_dim,deg_dim))
1509         ABI_ALLOCATE(unitary_tr,(deg_dim,deg_dim))
1510         ABI_ALLOCATE(eigf3d,(deg_dim))
1511         ABI_ALLOCATE(saddle_warn,(deg_dim))
1512         ABI_ALLOCATE(start_eigf3d_pos,(deg_dim))
1513         ABI_ALLOCATE(m_avg,(deg_dim))
1514         ABI_ALLOCATE(m_avg_frohlich,(deg_dim))
1515         ABI_ALLOCATE(m_cart,(ndirs,deg_dim))
1516         ABI_ALLOCATE(deigf3d_dth,(deg_dim))
1517         ABI_ALLOCATE(deigf3d_dph,(deg_dim))
1518         ABI_ALLOCATE(unit_speed,(mdim,deg_dim))
1519         ABI_ALLOCATE(transport_tensor,(mdim,mdim,deg_dim))
1520         ABI_ALLOCATE(transport_tensor_eig,(mdim))
1521         ABI_ALLOCATE(transport_eqv_m,(mdim,mdim,deg_dim))
1522         ABI_ALLOCATE(transport_eqv_eigval,(mdim,deg_dim))
1523         ABI_ALLOCATE(transport_eqv_eigvec,(mdim,mdim,deg_dim))
1524         ABI_ALLOCATE(prodc,(deg_dim,deg_dim))
1525         ABI_ALLOCATE(prodr,(mdim,mdim))
1526         !ABI_ALLOCATE(f3dfd,(2,nphi,deg_dim))
1527         unit_r=zero
1528         dr_dth=zero
1529         dr_dph=zero
1530         f3d=zero
1531         df3d_dth=zero
1532         df3d_dph=zero
1533         unitary_tr=zero
1534         eigf3d=zero
1535         saddle_warn=.false.
1536         start_eigf3d_pos=.true.
1537         m_avg=zero
1538         m_avg_frohlich=zero
1539         m_cart=zero
1540         deigf3d_dth=zero
1541         deigf3d_dph=zero
1542         unit_speed=zero
1543         transport_tensor=zero
1544         transport_tensor_eig=zero
1545         transport_eqv_m=zero
1546         transport_eqv_eigval=zero
1547         transport_eqv_eigvec=zero
1548 
1549         !Hack to print f(theta,phi) & weights to a file
1550         if(print_fsph) then
1551           write(msg,*) degl+1
1552           filename='f_band_'//TRIM(ADJUSTL(msg))//'-'
1553           write(msg,*) degl+deg_dim
1554           filename=TRIM(filename)//TRIM(ADJUSTL(msg))//'.dat'
1555           io_unit = get_unit()
1556           open(io_unit,file=TRIM(filename),status='replace')
1557           write(io_unit,*) 'ntheta=',ntheta,', nphi=',nphi
1558           write(io_unit,*) 'itheta, iphi, weight, f_n(theta,phi)'
1559         end if
1560 
1561         do itheta=1,ntheta
1562           costh=gq_points_costh(itheta) ; sinth=gq_points_sinth(itheta)
1563           do iphi=1,nphi
1564             cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
1565             weight=gq_weights_th(itheta)*gq_weights_ph(iphi)
1566 
1567             unit_r(1)=sinth*cosph 
1568             unit_r(2)=sinth*sinph 
1569             unit_r(3)=costh
1570 
1571             dr_dth(1)=costh*cosph
1572             dr_dth(2)=costh*sinph 
1573             dr_dth(3)=-sinth
1574 
1575             dr_dph(1)=-sinph  !sin(theta)*
1576             dr_dph(2)=cosph   !cos(theta)*
1577             dr_dph(3)=zero
1578 
1579             do iband=1,deg_dim
1580               do jband=1,deg_dim
1581                 f3d(iband,jband)=DOT_PRODUCT(unit_r,MATMUL(eig2_diag_cart(:,:,iband,jband),unit_r))
1582                 df3d_dth(iband,jband)=DOT_PRODUCT(dr_dth,MATMUL(eig2_diag_cart(:,:,iband,jband),unit_r))+&
1583 &                DOT_PRODUCT(unit_r,MATMUL(eig2_diag_cart(:,:,iband,jband),dr_dth))
1584                 df3d_dph(iband,jband)=DOT_PRODUCT(dr_dph,MATMUL(eig2_diag_cart(:,:,iband,jband),unit_r))+&
1585 &                DOT_PRODUCT(unit_r,MATMUL(eig2_diag_cart(:,:,iband,jband),dr_dph))
1586               end do
1587             end do
1588             !DIAGONALIZATION
1589             eigenvec = f3d        !IN
1590             lwork=-1
1591             ABI_ALLOCATE(work,(1))
1592             ABI_ALLOCATE(rwork,(3*deg_dim-2))
1593             call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1594             lwork=int(work(1))
1595             ABI_DEALLOCATE(work)
1596             eigenval = zero
1597             ABI_ALLOCATE(work,(lwork))
1598             work=zero; rwork=zero
1599             call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1600             ABI_DEALLOCATE(rwork)
1601             ABI_DEALLOCATE(work)
1602             unitary_tr = eigenvec !OUT
1603             eigf3d = eigenval     !OUT
1604             if(itheta==1 .and. iphi==1) start_eigf3d_pos = eigf3d > 0
1605             do iband=1,deg_dim
1606               if(start_eigf3d_pos(iband) .neqv. (eigf3d(iband)>0)) then
1607                 saddle_warn(iband)=.true.
1608               end if
1609             end do
1610 
1611             !Hack to print f(theta,phi)
1612             if(print_fsph) write(io_unit,*) gq_points_th(itheta), gq_points_ph(iphi), weight, eigf3d(:)
1613 
1614             !!DEBUG-Mech.
1615             !!A=-4.20449; B=0.378191; C=5.309  !Mech's fit
1616             !A=-4.62503023; B=0.68699088; C=5.20516873 !My fit
1617             !R = sqrt(B**2 + C**2*sin(theta)**2*(cos(theta)**2 + sin(theta)**2*sin(phi)**2*cos(phi)**2))
1618             !eigf3d(1) = A - R
1619             !eigf3d(2) = A + R
1620 
1621             !!angular FD
1622             !f3dfd(2,iphi,:)=eigf3d(:)
1623 
1624             m_avg = m_avg + weight*sinth*eigf3d
1625             m_avg_frohlich = m_avg_frohlich + weight*sinth/(abs(eigf3d))**half
1626 
1627             prodc=MATMUL_(f3d,unitary_tr,deg_dim,deg_dim) ; f3d=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1628             !f3d = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(f3d,unitary_tr))
1629             do iband=1,deg_dim
1630               eigf3d(iband) = real(f3d(iband,iband),dp)
1631             end do
1632             prodc=MATMUL_(df3d_dth,unitary_tr,deg_dim,deg_dim) ; df3d_dth=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1633             !df3d_dth=MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(df3d_dth,unitary_tr))
1634             do iband=1,deg_dim
1635               deigf3d_dth(iband) = real(df3d_dth(iband,iband),dp)
1636             end do
1637             prodc=MATMUL_(df3d_dph,unitary_tr,deg_dim,deg_dim) ; df3d_dph=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1638             !df3d_dph = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(df3d_dph,unitary_tr))
1639             do iband=1,deg_dim
1640               deigf3d_dph(iband) = real(df3d_dph(iband,iband),dp)
1641             end do
1642 
1643             !!DEBUG-Mech.
1644             !eigf3d(1) = A - R
1645             !eigf3d(2) = A + R
1646             !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))
1647             !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))
1648             !deigf3d_dph(1) = -1./2./R*C**2*sin(theta)**3*(2.*sin(phi)*cos(phi)**3 - 2.*sin(phi)**3*cos(phi))
1649             !deigf3d_dph(2) =  1./2./R*C**2*sin(theta)**3*(2.*sin(phi)*cos(phi)**3 - 2.*sin(phi)**3*cos(phi))
1650 
1651             !!angular FD
1652             !if(iphi/=1 .and. itheta/=1) then
1653             !  deigf3d_dph(:) = (f3dfd(2,iphi,:)-f3dfd(2,iphi-1,:))/two_pi*nphi/sin(theta)
1654             !else
1655             !  deigf3d_dph(:) = zero
1656             !end if
1657             !if(itheta/=1) then
1658             !  deigf3d_dth(:) = (f3dfd(2,iphi,:)-f3dfd(1,iphi,:))/pi*ntheta
1659             !else
1660             !  deigf3d_dth(:) = zero
1661             !end if
1662 
1663             unit_speed(1,:) = 2._dp*sinth*cosph*eigf3d + costh*cosph*deigf3d_dth - sinph*deigf3d_dph      !/sin(theta)
1664             unit_speed(2,:) = 2._dp*sinth*sinph*eigf3d + costh*sinph*deigf3d_dth + cosph*deigf3d_dph      !/sin(theta)
1665             unit_speed(3,:) = 2._dp*costh*eigf3d - sinth*deigf3d_dth
1666 
1667             do jdeg=1,deg_dim
1668               do bdir=1,mdim
1669                 do adir=1,mdim
1670                   transport_tensor(adir,bdir,jdeg) = transport_tensor(adir,bdir,jdeg) + &
1671 &                weight*sinth*unit_speed(adir,jdeg)*unit_speed(bdir,jdeg)/(ABS(eigf3d(jdeg))**2.5_dp)
1672                 end do
1673               end do
1674             end do
1675           end do !iphi
1676           !!angular FD
1677           !f3dfd(1,:,:) = f3dfd(2,:,:)
1678         end do !itheta
1679 
1680         !Hack to print f(theta,phi)
1681         if(print_fsph) close(io_unit)
1682 
1683         m_avg = quarter/pi*m_avg
1684         m_avg = one/m_avg
1685 
1686         m_avg_frohlich = quarter/pi*m_avg_frohlich
1687         m_avg_frohlich = m_avg_frohlich**2
1688 
1689         transport_tensor = 1.0_dp/2.0_dp*transport_tensor
1690 
1691         !Effective masses along directions.
1692         do adir=1,ndirs
1693           do iband=1,deg_dim
1694             do jband=1,deg_dim
1695               f3d(iband,jband) = dot_product(dirs(:,adir),matmul(eig2_diag_cart(:,:,iband,jband),dirs(:,adir)))
1696             end do
1697           end do
1698           !f3d(:,:) = eig2_diag_cart(adir,adir,:,:)
1699           eigenvec = f3d        !IN
1700           lwork=-1
1701           ABI_ALLOCATE(work,(1))
1702           ABI_ALLOCATE(rwork,(3*deg_dim-2))
1703           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1704           lwork=int(work(1))
1705           ABI_DEALLOCATE(work)
1706           eigenval = zero
1707           ABI_ALLOCATE(work,(lwork))
1708           work=zero; rwork=zero
1709           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1710           ABI_DEALLOCATE(rwork)
1711           ABI_DEALLOCATE(work)
1712           unitary_tr = eigenvec !OUT
1713           eigf3d = eigenval     !OUT
1714           m_cart(adir,:)=1._dp/eigf3d(:)
1715         end do
1716 
1717         do iband=1,deg_dim
1718           !DIAGONALIZATION
1719           transport_eqv_eigvec(:,:,iband) = transport_tensor(:,:,iband)
1720           lwork=-1
1721           ABI_ALLOCATE(rwork,(1))
1722           call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_tensor_eig,rwork,lwork,info)
1723           lwork=int(rwork(1))
1724           ABI_DEALLOCATE(rwork)
1725           transport_tensor_eig = zero
1726           ABI_ALLOCATE(rwork,(lwork))
1727           rwork=zero
1728           call dsyev('V','U',mdim,transport_eqv_eigvec(:,:,iband),mdim,transport_tensor_eig,rwork,lwork,info)
1729           ABI_DEALLOCATE(rwork)
1730           transport_eqv_eigvec(:,:,iband) = transpose(transport_eqv_eigvec(:,:,iband)) !So that lines contain eigenvectors.
1731 
1732           prodr=MATMUL_(transport_tensor(:,:,iband),transport_eqv_eigvec(:,:,iband),mdim,mdim,transb='t')
1733           transport_tensor(:,:,iband)=MATMUL_(transport_eqv_eigvec(:,:,iband),prodr,mdim,mdim)
1734           !transport_tensor(:,:,iband) = MATMUL(transport_eqv_eigvec(:,:,iband), &
1735           !                              MATMUL(transport_tensor(:,:,iband),TRANSPOSE(transport_eqv_eigvec(:,:,iband))))
1736 
1737           transport_eqv_eigval(1,iband) = transport_tensor_eig(2)*transport_tensor_eig(3)*(3._dp/8._dp/pi)**2
1738           transport_eqv_eigval(2,iband) = transport_tensor_eig(3)*transport_tensor_eig(1)*(3._dp/8._dp/pi)**2
1739           transport_eqv_eigval(3,iband) = transport_tensor_eig(1)*transport_tensor_eig(2)*(3._dp/8._dp/pi)**2
1740           !The transport tensor loses the sign of the effective mass, this restores it.
1741           transport_eqv_eigval(:,iband) = DSIGN(transport_eqv_eigval(:,iband),m_avg(iband))
1742           transport_eqv_m(1,1,iband) = transport_eqv_eigval(1,iband)
1743           transport_eqv_m(2,2,iband) = transport_eqv_eigval(2,iband)
1744           transport_eqv_m(3,3,iband) = transport_eqv_eigval(3,iband)
1745 
1746           m_avg_frohlich(iband) = DSIGN(m_avg_frohlich(iband),m_avg(iband))
1747 
1748           prodr=MATMUL_(transport_eqv_m(:,:,iband),transport_eqv_eigvec(:,:,iband),mdim,mdim)
1749           transport_eqv_m(:,:,iband)=MATMUL_(transport_eqv_eigvec(:,:,iband),prodr,mdim,mdim,transa='t')
1750           !transport_eqv_m(:,:,iband) = MATMUL(TRANSPOSE(transport_eqv_eigvec(:,:,iband)), &
1751           !                             MATMUL(transport_eqv_m(:,:,iband),transport_eqv_eigvec(:,:,iband)))
1752 
1753         end do
1754 
1755         call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1756 &                        ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec)
1757         call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1758 &                        ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec)
1759 
1760         ABI_DEALLOCATE(unit_r)
1761         ABI_DEALLOCATE(dr_dth)
1762         ABI_DEALLOCATE(dr_dph)
1763         ABI_DEALLOCATE(f3d)
1764         ABI_DEALLOCATE(df3d_dth)
1765         ABI_DEALLOCATE(df3d_dph)
1766         ABI_DEALLOCATE(unitary_tr)
1767         ABI_DEALLOCATE(eigf3d)
1768         ABI_DEALLOCATE(saddle_warn)
1769         ABI_DEALLOCATE(start_eigf3d_pos)
1770         ABI_DEALLOCATE(m_avg)
1771         ABI_DEALLOCATE(m_avg_frohlich)
1772         ABI_DEALLOCATE(m_cart)
1773         ABI_DEALLOCATE(deigf3d_dth)
1774         ABI_DEALLOCATE(deigf3d_dph)
1775         ABI_DEALLOCATE(unit_speed)
1776         ABI_DEALLOCATE(transport_tensor)
1777         ABI_DEALLOCATE(transport_tensor_eig)
1778         ABI_DEALLOCATE(transport_eqv_m)
1779         ABI_DEALLOCATE(transport_eqv_eigval)
1780         ABI_DEALLOCATE(transport_eqv_eigvec)
1781         ABI_DEALLOCATE(prodc)
1782         ABI_DEALLOCATE(prodr)
1783         !ABI_DEALLOCATE(f3dfd)
1784 
1785       elseif (degenerate .and. mdim==2) then
1786 
1787         ABI_ALLOCATE(unit_r,(mdim))
1788         ABI_ALLOCATE(dr_dph,(mdim))
1789         ABI_ALLOCATE(f3d,(deg_dim,deg_dim))
1790         ABI_ALLOCATE(df3d_dph,(deg_dim,deg_dim))
1791         ABI_ALLOCATE(unitary_tr,(deg_dim,deg_dim))
1792         ABI_ALLOCATE(eigf3d,(deg_dim))
1793         ABI_ALLOCATE(saddle_warn,(deg_dim))
1794         ABI_ALLOCATE(start_eigf3d_pos,(deg_dim))
1795         ABI_ALLOCATE(m_avg,(deg_dim))
1796         ABI_ALLOCATE(m_avg_frohlich,(deg_dim))
1797         ABI_ALLOCATE(m_cart,(ndirs,deg_dim))
1798         ABI_ALLOCATE(deigf3d_dph,(deg_dim))
1799         ABI_ALLOCATE(unit_speed,(mdim,deg_dim))
1800         ABI_ALLOCATE(transport_tensor,(mdim,mdim,deg_dim))
1801         ABI_ALLOCATE(cart_rotation,(mdim,mdim))
1802         ABI_ALLOCATE(transport_tensor_eig,(mdim))
1803         ABI_ALLOCATE(transport_eqv_m,(mdim,mdim,deg_dim))
1804         ABI_ALLOCATE(transport_eqv_eigval,(mdim,deg_dim))
1805         ABI_ALLOCATE(transport_eqv_eigvec,(mdim,mdim,deg_dim))
1806         ABI_ALLOCATE(transport_tensor_scale,(deg_dim))
1807         ABI_ALLOCATE(prodc,(deg_dim,deg_dim))
1808         ABI_ALLOCATE(prodr,(mdim,mdim))
1809         unit_r=zero
1810         dr_dph=zero
1811         f3d=zero
1812         df3d_dph=zero
1813         unitary_tr=zero
1814         eigf3d=zero
1815         saddle_warn=.false.
1816         start_eigf3d_pos=.true.
1817         m_avg=zero
1818         m_avg_frohlich=zero
1819         m_cart=zero
1820         deigf3d_dph=zero
1821         unit_speed=zero
1822         transport_tensor=zero
1823         cart_rotation=zero
1824         transport_tensor_eig=zero
1825         transport_eqv_m=zero
1826         transport_eqv_eigval=zero
1827         transport_eqv_eigvec=zero
1828         transport_tensor_scale=zero
1829 
1830         do iphi=1,nphi
1831           cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi)
1832           weight=gq_weights_ph(iphi)
1833 
1834           unit_r(1)=cosph 
1835           unit_r(2)=sinph 
1836 
1837           dr_dph(1)=-sinph 
1838           dr_dph(2)=cosph 
1839 
1840           do iband=1,deg_dim
1841             do jband=1,deg_dim
1842               matr2d = eig2_diag_cart(1:mdim,1:mdim,iband,jband)
1843               f3d(iband,jband)=DOT_PRODUCT(unit_r,MATMUL(matr2d,unit_r))
1844               df3d_dph(iband,jband)=DOT_PRODUCT(dr_dph,MATMUL(matr2d,unit_r))+&
1845 &              DOT_PRODUCT(unit_r,MATMUL(matr2d,dr_dph))
1846             end do
1847           end do
1848 
1849           !DIAGONALIZATION
1850           eigenvec = f3d        !IN
1851           lwork=-1
1852           ABI_ALLOCATE(work,(1))
1853           ABI_ALLOCATE(rwork,(3*deg_dim-2))
1854           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1855           lwork=int(work(1))
1856           ABI_DEALLOCATE(work)
1857           eigenval = zero
1858           ABI_ALLOCATE(work,(lwork))
1859           work=zero; rwork=zero
1860           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1861           ABI_DEALLOCATE(rwork)
1862           ABI_DEALLOCATE(work)
1863           unitary_tr = eigenvec !OUT
1864           eigf3d = eigenval     !OUT
1865           if(iphi==1) start_eigf3d_pos = eigf3d > 0
1866           do iband=1,deg_dim
1867             if(start_eigf3d_pos(iband) .neqv. (eigf3d(iband)>0)) then
1868               saddle_warn(iband)=.true.
1869             end if
1870           end do
1871 
1872           m_avg = m_avg + weight*eigf3d
1873           m_avg_frohlich = m_avg_frohlich + weight/(abs(eigf3d))**half
1874 
1875           prodc=MATMUL_(f3d,unitary_tr,deg_dim,deg_dim) ; f3d=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1876           !f3d = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(f3d,unitary_tr))
1877           do iband=1,deg_dim
1878             eigf3d(iband) = real(f3d(iband,iband),dp)
1879           end do
1880           prodc=MATMUL_(df3d_dph,unitary_tr,deg_dim,deg_dim) ; df3d_dph=MATMUL_(unitary_tr,prodc,deg_dim,deg_dim,transa='c')
1881           !df3d_dph = MATMUL(CONJG(TRANSPOSE(unitary_tr)),MATMUL(df3d_dph,unitary_tr))
1882           do iband=1,deg_dim
1883             deigf3d_dph(iband) = real(df3d_dph(iband,iband),dp)
1884           end do
1885 
1886           unit_speed(1,:) = 2._dp*cosph*eigf3d - sinph*deigf3d_dph
1887           unit_speed(2,:) = 2._dp*sinph*eigf3d + cosph*deigf3d_dph
1888 
1889           do jdeg=1,deg_dim
1890             do bdir=1,mdim
1891               do adir=1,mdim
1892                 transport_tensor(adir,bdir,jdeg) = transport_tensor(adir,bdir,jdeg) + &
1893 &                weight*unit_speed(adir,jdeg)*unit_speed(bdir,jdeg)/(ABS(eigf3d(jdeg))**2)
1894               end do
1895             end do
1896           end do
1897 
1898         end do !iphi
1899 
1900         !!!DEBUG
1901 
1902         m_avg = half/pi*m_avg
1903         m_avg = one/m_avg
1904 
1905         m_avg_frohlich = half/pi*m_avg_frohlich
1906         m_avg_frohlich = m_avg_frohlich**2
1907 
1908         transport_tensor = 1.0_dp/2.0_dp*transport_tensor
1909 
1910         !Effective masses along directions.
1911         do adir=1,ndirs
1912           do iband=1,deg_dim
1913             do jband=1,deg_dim
1914               f3d(iband,jband) = dot_product(dirs(:,adir),matmul(eig2_diag_cart(:,:,iband,jband),dirs(:,adir)))
1915             end do
1916           end do
1917           !f3d(:,:) = eig2_diag_cart(adir,adir,:,:)
1918           eigenvec = f3d        !IN
1919           lwork=-1
1920           ABI_ALLOCATE(work,(1))
1921           ABI_ALLOCATE(rwork,(3*deg_dim-2))
1922           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1923           lwork=int(work(1))
1924           ABI_DEALLOCATE(work)
1925           eigenval = zero
1926           ABI_ALLOCATE(work,(lwork))
1927           work=zero; rwork=zero
1928           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
1929           ABI_DEALLOCATE(rwork)
1930           ABI_DEALLOCATE(work)
1931           unitary_tr = eigenvec !OUT
1932           eigf3d = eigenval     !OUT
1933           m_cart(adir,:)=1._dp/eigf3d(:)
1934         end do
1935 
1936         do iband=1,deg_dim
1937             !DIAGONALIZATION
1938           cart_rotation = transport_tensor(:,:,iband)
1939           lwork=-1
1940           ABI_ALLOCATE(rwork,(1))
1941           call dsyev('V','U',mdim,cart_rotation,mdim,transport_tensor_eig,rwork,lwork,info)
1942           lwork=int(rwork(1))
1943           ABI_DEALLOCATE(rwork)
1944           transport_tensor_eig = zero
1945           ABI_ALLOCATE(rwork,(lwork))
1946           rwork=zero
1947           call dsyev('V','U',mdim,cart_rotation,mdim,transport_tensor_eig,rwork,lwork,info)
1948           ABI_DEALLOCATE(rwork)
1949           transport_eqv_eigvec(:,:,iband) = transpose(cart_rotation(:,:)) !So that lines contain eigenvectors, not columns.
1950 
1951           prodr=MATMUL_(transport_tensor(:,:,iband),cart_rotation,mdim,mdim)
1952           transport_tensor(:,:,iband)=MATMUL_(cart_rotation,prodr,mdim,mdim,transa='t')
1953           !transport_tensor(:,:,iband) = MATMUL(TRANSPOSE(cart_rotation),MATMUL(transport_tensor(:,:,iband),cart_rotation))
1954 
1955           transport_eqv_eigval(1,iband) = 0.5*m_avg(iband)*(1.0 + transport_tensor_eig(2)/transport_tensor_eig(1))
1956           transport_eqv_eigval(2,iband) = transport_eqv_eigval(1,iband)*transport_tensor_eig(1)/transport_tensor_eig(2)
1957           !The transport tensor loses the sign of the effective mass, this restores it.
1958           transport_eqv_eigval(:,iband) = SIGN(transport_eqv_eigval(:,iband),m_avg(iband))
1959           transport_eqv_m(1,1,iband) = transport_eqv_eigval(1,iband)
1960           transport_eqv_m(2,2,iband) = transport_eqv_eigval(2,iband)
1961           transport_tensor_scale(iband) = sqrt(transport_tensor_eig(1)*transport_tensor_eig(2))/two_pi
1962 
1963           m_avg_frohlich(iband) = SIGN(m_avg_frohlich(iband),m_avg(iband))
1964 
1965           prodr=MATMUL_(transport_eqv_m(:,:,iband),cart_rotation,mdim,mdim,transb='t')
1966           transport_eqv_m(:,:,iband)=MATMUL_(cart_rotation,prodr,mdim,mdim)
1967           !transport_eqv_m(:,:,iband) = MATMUL(cart_rotation,MATMUL(transport_eqv_m(:,:,iband),TRANSPOSE(cart_rotation)))
1968 
1969         end do
1970 
1971         call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1972 &                        ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec,transport_tensor_scale)
1973         call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m, &
1974 &                        ntheta,m_avg,m_avg_frohlich,saddle_warn,transport_eqv_eigval,transport_eqv_eigvec,transport_tensor_scale)
1975 
1976         ABI_DEALLOCATE(unit_r)
1977         ABI_DEALLOCATE(dr_dph)
1978         ABI_DEALLOCATE(f3d)
1979         ABI_DEALLOCATE(df3d_dph)
1980         ABI_DEALLOCATE(unitary_tr)
1981         ABI_DEALLOCATE(eigf3d)
1982         ABI_DEALLOCATE(saddle_warn)
1983         ABI_DEALLOCATE(start_eigf3d_pos)
1984         ABI_DEALLOCATE(m_avg)
1985         ABI_DEALLOCATE(m_avg_frohlich)
1986         ABI_DEALLOCATE(m_cart)
1987         ABI_DEALLOCATE(deigf3d_dph)
1988         ABI_DEALLOCATE(unit_speed)
1989         ABI_DEALLOCATE(transport_tensor)
1990         ABI_DEALLOCATE(cart_rotation)
1991         ABI_DEALLOCATE(transport_tensor_eig)
1992         ABI_DEALLOCATE(transport_eqv_m)
1993         ABI_DEALLOCATE(transport_eqv_eigval)
1994         ABI_DEALLOCATE(transport_eqv_eigvec)
1995         ABI_DEALLOCATE(transport_tensor_scale)
1996         ABI_DEALLOCATE(prodc)
1997         ABI_DEALLOCATE(prodr)
1998 
1999       elseif (degenerate .and. mdim==1) then
2000 
2001         ABI_ALLOCATE(f3d,(deg_dim,deg_dim))
2002         ABI_ALLOCATE(unitary_tr,(deg_dim,deg_dim))
2003         ABI_ALLOCATE(eigf3d,(deg_dim))
2004         ABI_ALLOCATE(m_cart,(ndirs,deg_dim))
2005         ABI_ALLOCATE(transport_eqv_m,(mdim,mdim,deg_dim))
2006         ABI_ALLOCATE(m_avg,(deg_dim))
2007         ABI_ALLOCATE(m_avg_frohlich,(deg_dim))
2008         ABI_ALLOCATE(saddle_warn,(deg_dim))
2009 
2010 
2011         f3d=zero
2012         unitary_tr=zero
2013         eigf3d=zero
2014         m_cart=zero
2015         transport_eqv_m=zero
2016         m_avg=zero
2017         m_avg_frohlich=zero
2018         saddle_warn=.false.
2019 
2020         f3d(:,:) = eig2_diag_cart(1,1,:,:)
2021 
2022         !DIAGONALIZATION
2023         eigenvec = f3d        !IN
2024         lwork=-1
2025         ABI_ALLOCATE(work,(1))
2026         ABI_ALLOCATE(rwork,(3*deg_dim-2))
2027         call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
2028         lwork=int(work(1))
2029         ABI_DEALLOCATE(work)
2030         eigenval = zero
2031         ABI_ALLOCATE(work,(lwork))
2032         work=zero; rwork=zero
2033         call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
2034         ABI_DEALLOCATE(rwork)
2035         ABI_DEALLOCATE(work)
2036         unitary_tr = eigenvec !OUT
2037         eigf3d = eigenval     !OUT
2038 
2039         transport_eqv_m(1,1,:)=1._dp/eigf3d(:)
2040 
2041         !Effective masses along directions.
2042         do adir=1,ndirs
2043           do iband=1,deg_dim
2044             do jband=1,deg_dim
2045               f3d(iband,jband) = dot_product(dirs(:,adir),matmul(eig2_diag_cart(:,:,iband,jband),dirs(:,adir)))
2046             end do
2047           end do
2048           eigenvec = f3d        !IN
2049           lwork=-1
2050           ABI_ALLOCATE(work,(1))
2051           ABI_ALLOCATE(rwork,(3*deg_dim-2))
2052           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
2053           lwork=int(work(1))
2054           ABI_DEALLOCATE(work)
2055           eigenval = zero
2056           ABI_ALLOCATE(work,(lwork))
2057           work=zero; rwork=zero
2058           call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info)
2059           ABI_DEALLOCATE(rwork)
2060           ABI_DEALLOCATE(work)
2061           eigf3d = eigenval     !OUT
2062           m_cart(adir,:)=1._dp/eigf3d(:)
2063         end do
2064 
2065         call print_tr_efmas(std_out,kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m,&
2066 &          ntheta,m_avg,m_avg_frohlich,saddle_warn)
2067         call print_tr_efmas(ab_out, kpt_rbz(:,ikpt),degl+1,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,transport_eqv_m,&
2068 &          ntheta,m_avg,m_avg_frohlich,saddle_warn)
2069 
2070         ABI_DEALLOCATE(f3d)
2071         ABI_DEALLOCATE(unitary_tr)
2072         ABI_DEALLOCATE(eigf3d)
2073         ABI_DEALLOCATE(m_cart)
2074         ABI_DEALLOCATE(transport_eqv_m)
2075         ABI_DEALLOCATE(m_avg)
2076         ABI_DEALLOCATE(m_avg_frohlich)
2077         ABI_DEALLOCATE(saddle_warn)
2078 
2079 
2080       end if !(degenerate)
2081 
2082       ABI_DEALLOCATE(eig2_diag_cart)
2083 
2084       ABI_DEALLOCATE(eigenval)
2085 
2086       ABI_DEALLOCATE(eigenvec)
2087     end do !ideg
2088 
2089   end do !ikpt
2090 
2091   ABI_DEALLOCATE(eff_mass)
2092   ABI_DEALLOCATE(dirs)
2093   ABI_DEALLOCATE(gq_points_th)
2094   ABI_DEALLOCATE(gq_points_costh)
2095   ABI_DEALLOCATE(gq_points_sinth)
2096   ABI_DEALLOCATE(gq_weights_th)
2097   ABI_DEALLOCATE(gq_points_ph)
2098   ABI_DEALLOCATE(gq_points_cosph)
2099   ABI_DEALLOCATE(gq_points_sinph)
2100   ABI_DEALLOCATE(gq_weights_ph)
2101 
2102   write(std_out,'(3a)') ch10,' END OF EFFECTIVE MASSES SECTION',ch10
2103   write(ab_out, '(3a)') ch10,' END OF EFFECTIVE MASSES SECTION',ch10
2104 
2105  end subroutine efmas_analysis

m_efmas/efmas_main [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 efmas_main

FUNCTION

 This routine calculates the generalized second-order k-derivative, Eq.66 of Laflamme2016,
 in reduced coordinates.

INPUTS

  cg(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz)=pw coefficients of GS wavefunctions at k.
  cg1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppo*nkpt_rbz,3,mpert) = first-order wf in G
            space for each perturbation. The wavefunction is orthogonal to the
            active space.
  dim_eig2rf = 1 if cg1_pert, gh0c1_pert and gh1c_pert are allocated.
               0 otherwise.
  dtset = dataset structure containing the input variable of the calculation.
  efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point
  eigen0(nkpt_rbz*dtset%mband*dtset%nsppol) = 0-order eigenvalues at all K-points:
            <k,n'|H(0)|k,n'> (hartree).
  eigen1(nkpt_rbz*2*dtset%nsppol*dtset%mband**2,3,mpert) = matrix of first-order:
            <k+Q,n'|H(1)|k,n> (hartree) (calculated in dfpt_cgwf).
  gh0c1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz,3,mpert) = matrix containing the
            vector:  <G|H(0)|psi(1)>, for each perturbation.
  gh1c_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz,3,mpert)) = matrix containing the
            vector:  <G|H(1)|n,k>, for each perturbation. The wavefunction is
            orthogonal to the active space.
  istwfk_pert(nkpt_rbz,3,mpert) = integer for choice of storage of wavefunction at
            each k point for each perturbation.
  mpert = maximum number of perturbations.
  mpi_enreg = informations about MPI parallelization.
  nkpt_rbz = number of k-points for each perturbation.
  npwarr(nkpt_rbz,mpert) = array of numbers of plane waves for each k-point
  rprimd(3,3)=dimensional primitive translations for real space (bohr)

OUTPUT

SIDE EFFECTS

  efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= generalized 2nd-order k-derivatives of eigenvalues
    efmasval(:,:)%ch2c INPUT : frozen wavefunction H2 contribution to generalized 2nd order k-derivatives of eigenenergy
    efmasval(:,:)%eig2_diag OUTPUT : generalized 2nd order k-derivatives of eigenenergy

PARENTS

      dfpt_looppert

CHILDREN

      cgqf,dgemm,dgetrf,dgetri,dotprod_g,dsyev,print_tr_efmas,zgemm,zgetrf
      zgetri,zheev

SOURCE

 929  subroutine efmas_main(cg,cg1_pert,dim_eig2rf,dtset,efmasdeg,efmasval,eigen0,&
 930 &   eigen1,gh0c1_pert,gh1c_pert,istwfk_pert,mpert,mpi_enreg,nkpt_rbz,npwarr,rprimd)
 931 
 932   use m_cgtools
 933   use m_gaussian_quadrature, only : cgqf
 934   use m_io_tools, only : get_unit
 935 
 936 !This section has been created automatically by the script Abilint (TD).
 937 !Do not modify the following lines by hand.
 938 #undef ABI_FUNC
 939 #define ABI_FUNC 'efmas_main'
 940 !End of the abilint section
 941 
 942   implicit none
 943 
 944  !Arguments ------------------------------------
 945  !scalars
 946   integer,            intent(in)    :: dim_eig2rf,mpert,nkpt_rbz
 947   type(dataset_type), intent(in)    :: dtset
 948   type(MPI_type),     intent(in) :: mpi_enreg
 949  !arrays
 950   integer,  intent(in) :: istwfk_pert(nkpt_rbz,3,mpert)
 951   integer,  intent(in) :: npwarr(nkpt_rbz,mpert)
 952   real(dp), intent(in) :: cg1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz*dim_eig2rf,3,mpert)
 953   real(dp), intent(in) :: gh0c1_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz*dim_eig2rf,3,mpert)
 954   real(dp), intent(in) :: gh1c_pert(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz*dim_eig2rf,3,mpert)
 955   real(dp), intent(in) :: eigen0(nkpt_rbz*dtset%mband*dtset%nsppol)
 956   real(dp), intent(in) :: eigen1(nkpt_rbz*2*dtset%nsppol*dtset%mband**2,3,mpert)
 957   real(dp), intent(in) :: rprimd(3,3)
 958   real(dp), intent(in) :: cg(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%nsppol*nkpt_rbz)
 959   type(efmasdeg_type), allocatable,intent(in) :: efmasdeg(:)
 960   type(efmasval_type),  allocatable,intent(inout) :: efmasval(:,:)
 961 
 962  !Local variables-------------------------------
 963   logical :: degenerate
 964   logical :: debug
 965   integer :: ipert
 966   integer :: isppol
 967   integer :: icg2   !TODOM : Reactivate the sections for icg2 / allow choice of k-point other than the first in the list.
 968   integer :: npw_k
 969   integer :: nband_k
 970   integer :: nspinor
 971   integer :: ideg
 972   integer :: ikpt
 973   integer :: istwf_k
 974   integer :: master,me,spaceworld
 975   integer :: band2tot_index
 976   integer :: bandtot_index
 977   integer :: iband, jband, kband
 978   integer :: adir,bdir
 979   integer :: deg_dim
 980   integer :: degl
 981   integer :: lwork
 982   character(len=500) :: msg
 983   real(dp) :: deltae
 984   real(dp) :: dot2i,dot2r,dot3i,dot3r,doti,dotr
 985   real(dp), allocatable :: cg0(:,:)
 986   real(dp), allocatable :: cg1_pert2(:,:),cg1_pert1(:,:)
 987   real(dp), allocatable :: gh1c_pert2(:,:),gh1c_pert1(:,:),gh0c1_pert1(:,:)
 988   complex(dpc) :: eig2_part(3,3)
 989   complex(dpc) :: eig2_ch2c(3,3)
 990   complex(dpc) :: eig2_paral(3,3)
 991   complex(dpc) :: eig2_gauge_change(3,3)
 992   complex(dpc) :: eig1a, eig1b, g_ch
 993   complex(dpc), allocatable :: eigen1_deg(:,:)
 994   complex(dpc), allocatable :: eig2_diag(:,:,:,:)
 995   complex(dpc), allocatable :: eig2_diag_cart(:,:,:,:)
 996 
 997 ! *********************************************************************
 998 
 999   debug = .false. ! Prints additional info to std_out
1000 
1001 ! Init parallelism
1002   master =0
1003   spaceworld=mpi_enreg%comm_cell
1004   me=mpi_enreg%me_kpt
1005 
1006   write(msg,'(4a)') ch10,&
1007 &  ' CALCULATION OF EFFECTIVE MASSES',ch10,&
1008 &  ' NOTE : Additional infos (eff. mass eigenvalues, eigenvectors and, if degenerate, average mass) are available in stdout.'
1009   call wrtout(std_out,msg,'COLL')
1010   call wrtout(ab_out,msg,'COLL')
1011 
1012   if(dtset%nsppol/=1)then
1013     write(msg,'(a,i3,a)') 'nsppol=',dtset%nsppol,' is not yet treated in m_efmas.'
1014     MSG_ERROR(msg)
1015   end if
1016   if(dtset%nspden/=1)then
1017     write(msg,'(a,i3,a)') 'nspden=',dtset%nspden,' is not yet treated in m_efmas.'
1018     MSG_ERROR(msg)
1019   end if
1020   if(dtset%efmas_deg==0) then
1021     write(msg,'(a)') 'efmas_deg==0 is for debugging; the results for degenerate bands will be garbage.'
1022     MSG_WARNING(msg)
1023     write(ab_out,'(6a)') ch10,'--- !WARNING',ch10,TRIM(msg),ch10,'---'
1024   end if
1025 
1026   ipert = dtset%natom+1
1027   isppol = 1
1028 
1029   icg2 = 0
1030   band2tot_index=0
1031   bandtot_index=0
1032 
1033   !XG20180519 : in the original coding by Jonathan, there is a lack of care about using dtset%nkpt or nkpt_rbz ...
1034   !Not important in the sequential case (?!) but likely problematic in the parallel case.
1035   do ikpt=1,dtset%nkpt
1036     npw_k = npwarr(ikpt,ipert)
1037     nband_k = dtset%nband(ikpt)
1038     nspinor = dtset%nspinor
1039 
1040     ABI_ALLOCATE(cg1_pert2,(2,npw_k*nspinor))
1041     ABI_ALLOCATE(cg1_pert1,(2,npw_k*nspinor))
1042     ABI_ALLOCATE(gh1c_pert2,(2,npw_k*nspinor))
1043     ABI_ALLOCATE(gh1c_pert1,(2,npw_k*nspinor))
1044     ABI_ALLOCATE(gh0c1_pert1,(2,npw_k*nspinor))
1045     ABI_ALLOCATE(cg0,(2,npw_k*nspinor))
1046 
1047     do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
1048       deg_dim    = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
1049       degenerate = (deg_dim>1) .and. (dtset%efmas_deg/=0)
1050       degl       = efmasdeg(ikpt)%degs_bounds(1,ideg)-1
1051 
1052       ABI_ALLOCATE(eigen1_deg,(deg_dim,deg_dim))
1053       !!! If treated band degenerate at 0th order, check that we are at extrema.
1054       if(degenerate) then
1055         do adir=1,3
1056           do iband=1,deg_dim
1057             do jband=1,deg_dim
1058               eigen1_deg(iband,jband) = cmplx(eigen1(2*(jband+degl)-1+(iband+degl-1)*2*nband_k,adir,ipert),&
1059 &              eigen1(2*(jband+degl)+(iband+degl-1)*2*nband_k,adir,ipert),dpc)
1060             end do
1061           end do
1062           if (.not.(ALL(ABS(eigen1_deg)<tol5))) then
1063             write(msg,'(a,a)') ' Effective masses calculations require given k-point(s) to be band extrema for given bands, ',&
1064 &            'but gradient of band(s) was found to be nonzero.'
1065             MSG_ERROR(msg)
1066           end if
1067         end do !adir=1,3
1068       end if !degenerate(1)
1069       ABI_DEALLOCATE(eigen1_deg)
1070 
1071       ABI_ALLOCATE(eig2_diag,(3,3,deg_dim,deg_dim))
1072       ABI_ALLOCATE(eig2_diag_cart,(3,3,deg_dim,deg_dim))
1073       eig2_diag = zero
1074 
1075       do iband=1,deg_dim
1076         write(std_out,*)"  In the set (possibly degenerate) compute band ",iband  ! This line here to avoid weird
1077         cg0(:,:) = cg(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2)
1078         do jband=1,deg_dim
1079           eig2_part = zero
1080           eig2_ch2c = zero
1081           eig2_paral = zero
1082           eig2_gauge_change = zero
1083           do adir=1,3
1084             istwf_k = istwfk_pert(ikpt,adir,ipert)
1085             do bdir=1,3
1086 
1087               ! Calculate the gauge change (to be subtracted to go from parallel transport to diagonal gauge).
1088               do kband=1,nband_k
1089                 !!! Equivalent to the gauge change in eig2stern.F90, but works also for other choices than the parallel gauge.
1090                 eig1a = cmplx( eigen1(2*kband-1+(degl+iband-1)*2*nband_k+band2tot_index,adir,ipert), &
1091 &                -eigen1(2*kband+(degl+iband-1)*2*nband_k+band2tot_index,adir,ipert), kind=dpc )
1092                 eig1b = cmplx( eigen1(2*kband-1+(degl+jband-1)*2*nband_k+band2tot_index,bdir,ipert), &
1093 &                eigen1(2*kband+(degl+jband-1)*2*nband_k+band2tot_index,bdir,ipert), kind=dpc )
1094                 g_ch = eig1a*eig1b
1095                 eig1a = cmplx( eigen1(2*kband-1+(degl+iband-1)*2*nband_k+band2tot_index,bdir,ipert), &
1096 &                -eigen1(2*kband+(degl+iband-1)*2*nband_k+band2tot_index,bdir,ipert), kind=dpc )
1097                 eig1b = cmplx( eigen1(2*kband-1+(degl+jband-1)*2*nband_k+band2tot_index,adir,ipert), &
1098 &                eigen1(2*kband+(degl+jband-1)*2*nband_k+band2tot_index,adir,ipert), kind=dpc )
1099                 g_ch = g_ch + eig1a*eig1b
1100 
1101                 deltae = eigen0(kband+bandtot_index) - eigen0((degl+iband)+bandtot_index)
1102                 if( kband<=degl.or.kband>degl+deg_dim) then
1103                   g_ch = g_ch/deltae
1104                 else
1105                   g_ch = zero
1106                 end if
1107                 eig2_gauge_change(adir,bdir) = eig2_gauge_change(adir,bdir) + g_ch
1108               end do !kband
1109 
1110               cg1_pert2(:,:)   = cg1_pert(:,1+(degl+jband-1)*npw_k*nspinor+icg2:(degl+jband)*npw_k*nspinor+icg2,bdir,ipert)
1111               cg1_pert1(:,:)   = cg1_pert(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2,adir,ipert)
1112               gh1c_pert1(:,:)  = gh1c_pert(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2,adir,ipert)
1113               gh1c_pert2(:,:)  = gh1c_pert(:,1+(degl+jband-1)*npw_k*nspinor+icg2:(degl+jband)*npw_k*nspinor+icg2,bdir,ipert)
1114               gh0c1_pert1(:,:) = gh0c1_pert(:,1+(degl+iband-1)*npw_k*nspinor+icg2:(degl+iband)*npw_k*nspinor+icg2,adir,ipert)
1115 
1116               ! The first two dotprod corresponds to:  <Psi(1)|H(1)|Psi(0)> + cc.
1117               ! They are calculated using wavefunctions <Psi(1)| that are orthogonal to the active space.
1118               dotr=zero ; doti=zero
1119               call dotprod_g(dotr,doti,istwf_k,npw_k*nspinor,2,cg1_pert1,gh1c_pert2,mpi_enreg%me_g0,&
1120 &              mpi_enreg%comm_spinorfft)
1121               dot2r=zero ; dot2i=zero
1122               call dotprod_g(dot2r,dot2i,istwf_k,npw_k*nspinor,2,gh1c_pert1,cg1_pert2,mpi_enreg%me_g0,&
1123 &              mpi_enreg%comm_spinorfft)
1124 
1125               ! This dotprod corresponds to : <Psi(1)|H(0)- E(0)|Psi(1)>
1126               ! It is calculated using wavefunctions that are orthogonal to the active space.
1127               dot3r=zero ; dot3i=zero
1128               call dotprod_g(dot3r,dot3i,istwf_k,npw_k*nspinor,2,gh0c1_pert1,cg1_pert2,mpi_enreg%me_g0,&
1129 &              mpi_enreg%comm_spinorfft)
1130 
1131               eig2_part(adir,bdir) = cmplx(dotr+dot2r+dot3r,doti+dot2i+dot3i,kind=dpc)
1132               !eig2_part(adir,bdir) = cmplx(dotr+dot2r,doti+dot2i,kind=dpc)  !DEBUG
1133               !eig2_part(adir,bdir) = cmplx(dotr,doti,kind=dpc)              !DEBUG
1134 
1135               eig2_ch2c(adir,bdir) = efmasval(ideg,ikpt)%ch2c(adir,bdir,iband,jband)
1136 
1137             end do !bdir
1138           end do  !adir
1139 
1140           do adir=1,3
1141             do bdir=1,3
1142               eig2_paral(adir,bdir) = eig2_part(adir,bdir) + eig2_part(bdir,adir) + eig2_ch2c(adir,bdir)
1143             end do
1144           end do
1145 
1146           eig2_diag(:,:,iband,jband) = eig2_paral - eig2_gauge_change
1147           efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)=eig2_diag(:,:,iband,jband)
1148 
1149           !!! Decomposition of the hessian in into its different contributions.
1150           if(debug) then
1151 
1152             eig2_diag_cart(:,:,iband,jband) = matmul(matmul(rprimd,eig2_diag(:,:,iband,jband)),transpose(rprimd))/two_pi**2
1153             eig2_paral                 = matmul(matmul(rprimd,eig2_paral),                transpose(rprimd))/two_pi**2
1154             eig2_gauge_change          = matmul(matmul(rprimd,eig2_gauge_change),         transpose(rprimd))/two_pi**2
1155             eig2_ch2c                  = matmul(matmul(rprimd,eig2_ch2c),                 transpose(rprimd))/two_pi**2
1156             eig2_part                  = matmul(matmul(rprimd,eig2_part),                 transpose(rprimd))/two_pi**2
1157 
1158             write(std_out,'(a)') 'Hessian of eigenvalues               = H. in parallel gauge - Gauge transformation'
1159             do adir=1,3
1160               write(std_out,'(3f12.8,2(a,3f12.8))')&
1161 &               real(eig2_diag_cart(adir,:,iband,jband),dp),' |',real(eig2_paral(adir,:),dp),' |',real(eig2_gauge_change(adir,:),dp)
1162             end do
1163             write(std_out,'(a)') 'H. in parallel gauge  = Second der. of H     + First derivatives    + First derivatives^T'
1164             do adir=1,3
1165               write(std_out,'(3f12.8,2(a,3f12.8))') real(eig2_paral(adir,:),dp),' |',real(eig2_ch2c(adir,:),dp),' |', &
1166 &                                                   real(eig2_part(adir,:),dp)
1167             end do
1168           end if !debug
1169 
1170           if(.not. degenerate .and. iband==jband) then
1171             write(std_out,'(a,3f20.16)') 'Gradient of eigenvalues = ',&
1172 &            matmul(rprimd,eigen1(2*(degl+iband)-1+(degl+iband-1)*2*nband_k+band2tot_index,:,ipert))/two_pi
1173           end if !.not.degenerate
1174 
1175         end do !jband
1176       end do !iband
1177 
1178       ABI_DEALLOCATE(eig2_diag)
1179       ABI_DEALLOCATE(eig2_diag_cart)
1180 
1181     end do !ideg
1182 
1183     ABI_DEALLOCATE(cg1_pert2)
1184     ABI_DEALLOCATE(cg1_pert1)
1185     ABI_DEALLOCATE(gh1c_pert2)
1186     ABI_DEALLOCATE(gh1c_pert1)
1187     ABI_DEALLOCATE(gh0c1_pert1)
1188     ABI_DEALLOCATE(cg0)
1189 
1190     icg2=icg2+npw_k*dtset%nspinor*nband_k
1191     bandtot_index=bandtot_index+nband_k
1192     band2tot_index=band2tot_index+2*nband_k**2
1193   end do ! ikpt
1194 
1195  end subroutine efmas_main

m_efmas/efmas_ncread [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 efmas_ncread

FUNCTION

 This routine reads from an EFMAS NetCDF file the information needed 
 to compute rapidly the band effective masses,
 namely, the generalized second-order k-derivatives of the eigenenergies,
 see Eq.(66) of Laflamme2016.

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

573  subroutine efmas_ncread(efmasdeg,efmasval,kpt,ncid)
574 
575 
576 !This section has been created automatically by the script Abilint (TD).
577 !Do not modify the following lines by hand.
578 #undef ABI_FUNC
579 #define ABI_FUNC 'efmas_ncread'
580 !End of the abilint section
581 
582  implicit none
583 
584 !Arguments ------------------------------------
585 !scalars
586  integer,            intent(in) :: ncid
587 !arrays
588  real(dp), allocatable,            intent(out) :: kpt(:,:)
589  type(efmasdeg_type), allocatable, intent(out) :: efmasdeg(:)
590  type(efmasval_type), allocatable, intent(out) :: efmasval(:,:)
591 
592 !Local variables-------------------------------
593  integer :: deg_dim,eig2_diag_arr_dim
594  integer :: iband,ideg,ideg_tot,ieig,ikpt
595  integer :: jband,mband,nband,ndegs,ndegs_tot,nkpt
596  integer, allocatable :: nband_arr(:)
597  integer, allocatable :: ndegs_arr(:)
598  integer, allocatable :: degs_range_arr(:,:)
599  integer, allocatable :: ideg_arr(:,:)
600  integer, allocatable :: degs_bounds_arr(:,:)
601  real(dp), allocatable :: ch2c_arr(:,:,:,:)
602  real(dp), allocatable :: eig2_diag_arr(:,:,:,:)
603 #ifdef HAVE_NETCDF
604  integer :: ncerr
605 #endif
606 !----------------------------------------------------------------------
607 
608 #ifdef HAVE_NETCDF
609  NCF_CHECK(nctk_set_datamode(ncid))
610  NCF_CHECK(nctk_get_dim(ncid, "number_of_kpoints", nkpt))
611  NCF_CHECK(nctk_get_dim(ncid, "max_number_of_states", mband))
612  NCF_CHECK(nctk_get_dim(ncid, "total_number_of_degenerate_sets", ndegs_tot))
613  NCF_CHECK(nctk_get_dim(ncid, "eig2_diag_arr_dim", eig2_diag_arr_dim))
614 
615 !Allocate the arrays to be read from NetCDF file
616  ABI_MALLOC(kpt, (3,nkpt) )
617  ABI_MALLOC(nband_arr, (nkpt) )
618  ABI_MALLOC(ndegs_arr, (nkpt) )
619  ABI_MALLOC(degs_range_arr, (2,nkpt) )
620  ABI_MALLOC(ideg_arr, (mband,nkpt) )
621  ABI_MALLOC(degs_bounds_arr, (2,ndegs_tot) )
622  ABI_MALLOC(ch2c_arr, (2,3,3,eig2_diag_arr_dim) )
623  ABI_MALLOC(eig2_diag_arr, (2,3,3,eig2_diag_arr_dim) )
624 
625 !Read from NetCDF file
626  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "reduced_coordinates_of_kpoints"), kpt))
627  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "number_of_states"), nband_arr))
628  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "number_of_degenerate_sets"), ndegs_arr))
629  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "degs_range_arr"),            degs_range_arr))
630  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "ideg_arr"),                  ideg_arr))
631  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "degs_bounds_arr"),           degs_bounds_arr))
632  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "ch2c_arr"),                  ch2c_arr))
633  NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "eig2_diag_arr"),             eig2_diag_arr))
634 
635 !Prepare the efmas* datastructures
636  ABI_DT_MALLOC(efmasdeg,(nkpt))
637  ABI_DT_MALLOC(efmasval,(mband,nkpt))
638 
639  ideg_tot=1
640  ieig=1
641  do ikpt=1,nkpt
642    efmasdeg(ikpt)%deg_range(:)=degs_range_arr(:,ikpt)
643    nband=nband_arr(ikpt)   
644    efmasdeg(ikpt)%nband=nband
645    ABI_MALLOC(efmasdeg(ikpt)%ideg, (nband))
646    efmasdeg(ikpt)%ideg=ideg_arr(1:nband,ikpt)
647    ndegs=ndegs_arr(ikpt) 
648    efmasdeg(ikpt)%ndegs=ndegs
649    ABI_MALLOC(efmasdeg(ikpt)%degs_bounds,(2,nband))
650    do ideg=1,ndegs
651      efmasdeg(ikpt)%degs_bounds(:,ideg)=degs_bounds_arr(:,ideg_tot)
652      ideg_tot=ideg_tot+1
653      if( efmasdeg(ikpt)%deg_range(1) <= ideg .and. ideg <= efmasdeg(ikpt)%deg_range(2) ) then
654        deg_dim=efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
655        ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(3,3,deg_dim,deg_dim))
656        ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(3,3,deg_dim,deg_dim))
657        efmasval(ideg,ikpt)%ch2c=zero
658        efmasval(ideg,ikpt)%eig2_diag=zero
659        do jband=1,deg_dim
660          do iband=1,deg_dim
661            efmasval(ideg,ikpt)%ch2c(:,:,iband,jband)=&
662 &           dcmplx(ch2c_arr(1,:,:,ieig+iband-1),ch2c_arr(2,:,:,ieig+iband-1))
663            efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband)=& 
664 &           dcmplx(eig2_diag_arr(1,:,:,ieig+iband-1),eig2_diag_arr(2,:,:,ieig+iband-1))
665          enddo
666          ieig=ieig+deg_dim
667        enddo
668      else
669        ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(0,0,0,0))
670        ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(0,0,0,0))
671      end if
672    end do
673  enddo
674 
675 !Deallocate the arrays
676  ABI_FREE(nband_arr)
677  ABI_FREE(ndegs_arr)
678  ABI_FREE(degs_range_arr)
679  ABI_FREE(ideg_arr)
680  ABI_FREE(degs_bounds_arr)
681  ABI_FREE(ch2c_arr)
682  ABI_FREE(eig2_diag_arr)
683 #endif
684 
685  end subroutine efmas_ncread

m_efmas/efmasdeg_free [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 efmasdeg_free

FUNCTION

 This routine deallocates an efmasdeg_type.

INPUTS

OUTPUT

PARENTS

      m_efmas

CHILDREN

SOURCE

181  subroutine efmasdeg_free(efmasdeg)
182 
183 
184 !This section has been created automatically by the script Abilint (TD).
185 !Do not modify the following lines by hand.
186 #undef ABI_FUNC
187 #define ABI_FUNC 'efmasdeg_free'
188 !End of the abilint section
189 
190    implicit none
191 
192    !Arguments ------------------------------------
193    type(efmasdeg_type),intent(inout) :: efmasdeg
194 
195    ! *********************************************************************
196    if(allocated(efmasdeg%degs_bounds)) then
197      ABI_FREE(efmasdeg%degs_bounds)
198    end if
199    if(allocated(efmasdeg%ideg)) then
200      ABI_FREE(efmasdeg%ideg)
201    end if
202  end subroutine efmasdeg_free

m_efmas/efmasdeg_free_array [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 efmasdeg_free_array

FUNCTION

 This routine deallocates an efmasdeg_type or, optionally, an array of efmasdeg_type.

INPUTS

OUTPUT

PARENTS

      respfn

CHILDREN

SOURCE

225  subroutine efmasdeg_free_array(efmasdeg)
226 
227 
228 !This section has been created automatically by the script Abilint (TD).
229 !Do not modify the following lines by hand.
230 #undef ABI_FUNC
231 #define ABI_FUNC 'efmasdeg_free_array'
232 !End of the abilint section
233 
234    implicit none
235 
236    !Arguments ------------------------------------
237    type(efmasdeg_type),allocatable,intent(inout) :: efmasdeg(:)
238 
239    !!!Local variables-------------------------------
240    integer :: i,n
241 
242    ! *********************************************************************
243    if(allocated(efmasdeg)) then
244      n=size(efmasdeg)
245      do i=1,n
246        call efmasdeg_free(efmasdeg(i))
247      end do
248      ABI_DATATYPE_DEALLOCATE(efmasdeg)
249    end if
250  end subroutine efmasdeg_free_array

m_efmas/efmasval_free [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 efmasval_free

FUNCTION

 This routine deallocates an efmasval_type.

INPUTS

OUTPUT

PARENTS

      m_efmas

CHILDREN

SOURCE

 80  subroutine efmasval_free(efmasval)
 81 
 82 
 83 !This section has been created automatically by the script Abilint (TD).
 84 !Do not modify the following lines by hand.
 85 #undef ABI_FUNC
 86 #define ABI_FUNC 'efmasval_free'
 87 !End of the abilint section
 88 
 89    implicit none
 90 
 91    !Arguments ------------------------------------
 92    type(efmasval_type),intent(inout) :: efmasval
 93 
 94    ! *********************************************************************
 95 
 96    if(allocated(efmasval%ch2c)) then
 97      ABI_FREE(efmasval%ch2c)
 98    end if
 99    if(allocated(efmasval%eig2_diag)) then
100      ABI_FREE(efmasval%eig2_diag)
101    end if
102 
103  end subroutine efmasval_free

m_efmas/MATMUL_DP [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 MATMUL_DP

FUNCTION

 Mimic MATMUL Fortran intrinsic function with BLAS3: C=A.B
 This is a temporary workaround to make tests pass on intel/mkl architectures
 Real version

INPUTS

  aa(:,:),bb(:,:)= input matrices
  mm,nn= sizes of output matrix
  [transa,transb]= equivalent to transa, transb args of gemm ('n','t','c')
                   if not present, default is 'n'.

OUTPUT

  MATMUL_DP(mm,nn)= output matrix A.B

SOURCE

2130 function MATMUL_DP(aa,bb,mm,nn,transa,transb)
2131 
2132 
2133 !This section has been created automatically by the script Abilint (TD).
2134 !Do not modify the following lines by hand.
2135 #undef ABI_FUNC
2136 #define ABI_FUNC 'MATMUL_DP'
2137 !End of the abilint section
2138 
2139  implicit none
2140 
2141 !Arguments ------------------------------------
2142 !scalars
2143  integer,intent(in) :: mm,nn
2144  character(len=1),optional,intent(in) :: transa,transb
2145 !arrays
2146  real(dp),intent(in) :: aa(:,:),bb(:,:)
2147  real(dp) :: MATMUL_DP(mm,nn)
2148 
2149 !Local variables-------------------------------
2150  integer :: kk,lda,ldb
2151  character(len=1) :: transa_,transb_
2152 
2153 ! *************************************************************************
2154 
2155  transa_='n';if (present(transa)) transa_=transa
2156  transb_='n';if (present(transb)) transb_=transb
2157 
2158  lda=size(aa,1) ; ldb=size(bb,1)
2159 
2160  if (transa_=='n') then
2161    kk=size(aa,2)
2162    if (size(aa,1)/=mm) then
2163      MSG_BUG('Error in sizes!')
2164    end if
2165  else
2166    kk=size(aa,1)
2167    if (size(aa,2)/=mm) then
2168      MSG_BUG('Error in sizes!')
2169    end if
2170  end if
2171 
2172  if (transb_=='n') then
2173    if (size(bb,1)/=kk.or.size(bb,2)/=nn) then
2174      MSG_BUG('Error in sizes!')
2175    end if
2176  else
2177    if (size(bb,1)/=nn.or.size(bb,2)/=kk) then
2178      MSG_BUG('Error in sizes!')
2179    end if
2180  end if
2181 
2182  call DGEMM(transa_,transb_,mm,nn,kk,one,aa,lda,bb,ldb,zero,MATMUL_DP,mm)
2183 
2184 end function MATMUL_DP

m_efmas/MATMUL_DPC [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 MATMUL_DPC

FUNCTION

 Mimic MATMUL Fortran intrinsic function with BLAS3
 This is a temporary workaround to make tests pass on intel/mkl architectures
 Complex version

INPUTS

  aa(:,:),bb(:,:)= input matrices
  mm,nn= sizes of output matrix
  [transa,transb]= equivalent to transa, transb args of gemm ('n','t','c')
                   if not present, default is 'n'.

OUTPUT

  MATMUL_DPC(:,:)= output matrix A.B

SOURCE

2209 function MATMUL_DPC(aa,bb,mm,nn,transa,transb)
2210 
2211 
2212 !This section has been created automatically by the script Abilint (TD).
2213 !Do not modify the following lines by hand.
2214 #undef ABI_FUNC
2215 #define ABI_FUNC 'MATMUL_DPC'
2216 !End of the abilint section
2217 
2218  implicit none
2219 
2220 !Arguments ------------------------------------
2221 !scalars
2222  integer,intent(in) :: mm,nn
2223  character(len=1),optional,intent(in) :: transa,transb
2224 !arrays
2225  complex(dpc),intent(in) :: aa(:,:),bb(:,:)
2226  complex(dpc) :: MATMUL_DPC(mm,nn)
2227 
2228 !Local variables-------------------------------
2229  integer :: kk,lda,ldb
2230  character(len=1) :: transa_,transb_
2231 
2232 ! *************************************************************************
2233 
2234  transa_='n';if (present(transa)) transa_=transa
2235  transb_='n';if (present(transb)) transb_=transb
2236 
2237  lda=size(aa,1) ; ldb=size(bb,1)
2238 
2239  if (transa_=='n') then
2240    kk=size(aa,2)
2241    if (size(aa,1)/=mm) then
2242      MSG_BUG('Error in sizes!')
2243    end if
2244  else
2245    kk=size(aa,1)
2246    if (size(aa,2)/=mm) then
2247      MSG_BUG('Error in sizes!')
2248    end if
2249  end if
2250 
2251  if (transb_=='n') then
2252    if (size(bb,1)/=kk.or.size(bb,2)/=nn) then
2253      MSG_BUG('Error in sizes!')
2254    end if
2255  else
2256    if (size(bb,1)/=nn.or.size(bb,2)/=kk) then
2257      MSG_BUG('Error in sizes!')
2258    end if
2259  end if
2260 
2261  call ZGEMM(transa_,transb_,mm,nn,kk,cone,aa,lda,bb,ldb,czero,MATMUL_DPC,mm)
2262 
2263 end function MATMUL_DPC

m_efmas/print_efmas [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 print_efmas

FUNCTION

 This routine prints the information needed to compute rapidly the band effective masses, 
 namely, the generalized second-order k-derivatives of the eigenenergies,
 see Eq.(66) of Laflamme2016. 

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

403  subroutine print_efmas(efmasdeg,efmasval,kpt,ncid)
404 
405 
406 !This section has been created automatically by the script Abilint (TD).
407 !Do not modify the following lines by hand.
408 #undef ABI_FUNC
409 #define ABI_FUNC 'print_efmas'
410 !End of the abilint section
411 
412  implicit none
413 
414 !Arguments ------------------------------------
415 !scalars
416  integer,            intent(in) :: ncid
417 !arrays
418  real(dp),            intent(in) :: kpt(:,:)
419  type(efmasdeg_type), intent(in) :: efmasdeg(:)
420  type(efmasval_type), intent(in) :: efmasval(:,:)
421 
422 !Local variables-------------------------------
423  integer :: deg_dim,eig2_diag_arr_dim
424  integer :: iband,ideg,ideg_tot,ieig,ikpt 
425  integer :: jband,mband,ndegs_tot,nkpt,nkptdeg,nkptval 
426  integer, allocatable :: nband_arr(:)
427 integer, allocatable :: ndegs_arr(:)
428  integer, allocatable :: degs_range_arr(:,:)
429  integer, allocatable :: ideg_arr(:,:)
430  integer, allocatable :: degs_bounds_arr(:,:)
431  real(dp), allocatable :: ch2c_arr(:,:,:,:)
432  real(dp), allocatable :: eig2_diag_arr(:,:,:,:)
433  character(len=500) :: msg
434 #ifdef HAVE_NETCDF
435  integer :: ncerr
436 #endif
437 !----------------------------------------------------------------------
438 
439 !XG20180519 Here, suppose that dtset%nkpt=nkpt_rbz (as done by Jonathan). 
440 !To be reexamined/corrected at the time of parallelization.
441 
442 #ifdef HAVE_NETCDF
443  nkptdeg=size(efmasdeg,1)
444  nkptval=size(efmasval,2)
445  if(nkptdeg/=nkptval) then
446    write(msg,'(a,i8,a,i8,a)') ' nkptdeg and nkptval =',nkptdeg,' and ',nkptval,' differ, which is inconsistent.'
447    MSG_ERROR(msg)
448  end if
449  nkpt=nkptdeg
450  if(nkpt/=size(kpt,2)) then
451    write(msg,'(a,i8,a,i8,a)') ' nkptdeg and nkpt =',nkptdeg,' and ',nkpt,' differ, which is inconsistent.'
452    MSG_ERROR(msg)
453  end if
454 
455  mband=size(efmasval,1)
456 
457 !Total number of (degenerate) sets over all k points
458  ndegs_tot=sum(efmasdeg%ndegs)
459 !Total number of generalized second-order k-derivatives
460  eig2_diag_arr_dim=zero
461  do ikpt=1,nkpt
462    do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
463      deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
464      eig2_diag_arr_dim = eig2_diag_arr_dim + deg_dim**2 
465    enddo
466  enddo
467 
468 !Allocate the arrays to be nc-written
469  ABI_MALLOC(nband_arr, (nkpt) )
470  ABI_MALLOC(ndegs_arr, (nkpt) )
471  ABI_MALLOC(degs_range_arr, (2,nkpt) )
472  ABI_MALLOC(ideg_arr, (mband,nkpt) )
473  ABI_MALLOC(degs_bounds_arr, (2,ndegs_tot) )
474  ABI_MALLOC(ch2c_arr, (2,3,3,eig2_diag_arr_dim) )
475  ABI_MALLOC(eig2_diag_arr, (2,3,3,eig2_diag_arr_dim) )
476 
477 !Prepare the arrays to be nc-written
478  ideg_tot=1
479  ieig=1
480  do ikpt=1,nkpt
481    nband_arr(ikpt)=efmasdeg(ikpt)%nband
482    ndegs_arr(ikpt)=efmasdeg(ikpt)%ndegs
483    degs_range_arr(:,ikpt)=efmasdeg(ikpt)%deg_range(:)
484    ideg_arr(:,ikpt)=0
485    ideg_arr(1:efmasdeg(ikpt)%nband,ikpt)=efmasdeg(ikpt)%ideg(:)
486    do ideg=1,efmasdeg(ikpt)%ndegs
487      degs_bounds_arr(:,ideg_tot)=efmasdeg(ikpt)%degs_bounds(:,ideg)
488      ideg_tot=ideg_tot+1
489    enddo
490    do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2)
491      deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
492      do jband=1,deg_dim
493        do iband=1,deg_dim
494          ch2c_arr(1,:,:,ieig+iband-1)=real(efmasval(ideg,ikpt)%ch2c(:,:,iband,jband))
495          ch2c_arr(2,:,:,ieig+iband-1)=aimag(efmasval(ideg,ikpt)%ch2c(:,:,iband,jband))
496          eig2_diag_arr(1,:,:,ieig+iband-1)=real(efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband))
497          eig2_diag_arr(2,:,:,ieig+iband-1)=aimag(efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband))
498        enddo
499        ieig=ieig+deg_dim
500      enddo
501    enddo
502  enddo
503 
504 !Define dimensions
505  ncerr=nctk_def_dims(ncid, [ &
506 &  nctkdim_t("number_of_reduced_dimensions", 3), &
507 &  nctkdim_t("real_or_complex", 2), &
508 &  nctkdim_t("number_of_kpoints", nkpt), &
509 &  nctkdim_t("max_number_of_states", mband), &
510 &  nctkdim_t("total_number_of_degenerate_sets", ndegs_tot), &
511 &  nctkdim_t("eig2_diag_arr_dim", eig2_diag_arr_dim)&
512 &  ], defmode=.True.)
513  NCF_CHECK(ncerr)
514 
515  ncerr = nctk_def_arrays(ncid, [ &
516 & nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), &
517 & nctkarr_t("number_of_states", "int", "number_of_kpoints"), &
518 & nctkarr_t("number_of_degenerate_sets", "int", "number_of_kpoints"), &
519 & nctkarr_t("degs_range_arr", "int", "two, number_of_kpoints"), &
520 & nctkarr_t("ideg_arr", "int", "max_number_of_states, number_of_kpoints"), &
521 & nctkarr_t("degs_bounds_arr", "int", "two, total_number_of_degenerate_sets"), &
522 & nctkarr_t("ch2c_arr", "dp", "real_or_complex, number_of_reduced_dimensions, number_of_reduced_dimensions, eig2_diag_arr_dim"),  &
523 & nctkarr_t("eig2_diag_arr","dp","real_or_complex, number_of_reduced_dimensions, number_of_reduced_dimensions, eig2_diag_arr_dim")&
524   ]) 
525  NCF_CHECK(ncerr)
526 
527  ! Write data.
528  NCF_CHECK(nctk_set_datamode(ncid))
529  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "reduced_coordinates_of_kpoints"), kpt))
530  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "number_of_states"), nband_arr))
531  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "number_of_degenerate_sets"), ndegs_arr))
532  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "degs_range_arr"),            degs_range_arr))
533  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ideg_arr"),                  ideg_arr))
534  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "degs_bounds_arr"),           degs_bounds_arr))
535  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ch2c_arr"),                  ch2c_arr))
536  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eig2_diag_arr"),             eig2_diag_arr))
537 
538 !Deallocate the arrays 
539  ABI_FREE(nband_arr)
540  ABI_FREE(ndegs_arr)
541  ABI_FREE(degs_range_arr)
542  ABI_FREE(ideg_arr)
543  ABI_FREE(degs_bounds_arr)
544  ABI_FREE(ch2c_arr)
545  ABI_FREE(eig2_diag_arr)
546 #endif
547 
548 end subroutine print_efmas

m_efmas/print_tr_efmas [ Functions ]

[ Top ] [ m_efmas ] [ Functions ]

NAME

 print_tr_efmas

FUNCTION

 This routine prints the transport equivalent effective mass and others info
 for a degenerate set of bands

INPUTS

OUTPUT

PARENTS

      m_efmas

CHILDREN

SOURCE

710  subroutine print_tr_efmas(io_unit,kpt,band,deg_dim,mdim,ndirs,dirs,m_cart,rprimd,efmas_tensor,ntheta, &
711 &                       m_avg,m_avg_frohlich,saddle_warn,efmas_eigval,efmas_eigvec,transport_tensor_scale)
712 
713 
714 !This section has been created automatically by the script Abilint (TD).
715 !Do not modify the following lines by hand.
716 #undef ABI_FUNC
717 #define ABI_FUNC 'print_tr_efmas'
718 !End of the abilint section
719 
720    implicit none
721 
722    !Arguments ------------------------------------
723    integer, intent(in) :: io_unit, band, deg_dim, mdim, ndirs
724    real(dp), intent(in) :: m_cart(ndirs,deg_dim), kpt(3), dirs(3,ndirs), rprimd(3,3), efmas_tensor(mdim,mdim,deg_dim)
725    integer, intent(in) :: ntheta
726    real(dp), intent(in) :: m_avg(deg_dim),m_avg_frohlich(deg_dim)
727    logical, intent(in) :: saddle_warn(deg_dim)
728    real(dp), intent(in), optional :: efmas_eigval(mdim,deg_dim)
729    real(dp), intent(in), optional :: efmas_eigvec(mdim,mdim,deg_dim)
730    real(dp), intent(in), optional :: transport_tensor_scale(deg_dim)
731 
732    !Local variables ------------------------------
733    logical :: extras
734    integer :: iband, adir
735    character(len=22) :: format_eigvec
736    character(len=500) :: msg, tmpstr
737    real(dp) :: vec(3)
738 
739    if(deg_dim>1) then
740      extras = present(efmas_eigval) .and. present(efmas_eigvec) 
741      if(mdim==3 .and. .not. extras) then
742        write(msg,'(a,l1,a,i1,a)') 'Subroutine print_tr_efmas called with degenerate=',deg_dim>1,&
743 &            ' and mdim=',mdim,', but missing required arguments for this case.'
744        MSG_ERROR(msg)
745      end if
746      if(mdim==2 .and. .not. (extras .or. present(transport_tensor_scale))) then
747        write(msg,'(a,l1,a,i1,a)') 'Subroutine print_tr_efmas called with degenerate=',deg_dim>1,&
748 &            ' and mdim=',mdim,', but missing required arguments for this case.'
749        MSG_ERROR(msg)
750      end if
751    else
752      extras = present(efmas_eigval) .and. present(efmas_eigvec)
753      if(mdim>1 .and. .not. extras) then
754        write(msg,'(a,l1,a,i1,a)') 'Subroutine print_tr_efmas called with degenerate=',deg_dim>1,&
755 &            ' and mdim=',mdim,', but missing required arguments for this case.'
756        MSG_ERROR(msg)
757      end if
758    end if
759 
760    if(deg_dim>1) then
761      write(io_unit,'(2a)') ch10,' COMMENTS: '
762      write(io_unit,'(a,3(f6.3,a),i5,a,i5)') ' - At k-point (',kpt(1),',',kpt(2),',',kpt(3),'), bands ',band,' through ',&
763 &          band+deg_dim-1
764      if(mdim>1) then
765        write(io_unit,'(a)') '   are DEGENERATE (effective mass tensor is therefore not defined).'
766        if(mdim==3) then
767          write(io_unit,'(a)') '   See Section IIIB Eqs. (67)-(70) and Appendix E of PRB 93 205147 (2016).' ! [[cite:Laflamme2016]]
768          write(io_unit,'(a,i7,a)') &
769 &          ' - Angular average effective mass for Frohlich model is to be averaged over degenerate bands. See later.'
770        elseif(mdim==2) then
771          write(io_unit,'(a)') ' - Also, 2D requested (perpendicular to Z axis).'
772          write(io_unit,'(a)') '   See Section IIIB and Appendix F, Eqs. (F12)-(F14) of PRB 93 205147 (2016).' ! [[cite:Laflamme2016]]
773        end if
774        write(io_unit,'(a,i7,a)') ' - Associated theta integrals calculated with ntheta=',ntheta,' points.'
775      else
776        write(io_unit,'(a)') '   are DEGENERATE.'
777        write(io_unit,'(a)') ' - Also, 1D requested (parallel to X axis).'
778      end if
779    end if
780 
781    if(ANY(saddle_warn)) then
782      write(msg,'(2a)') ch10,'Band(s)'
783      do iband=1,deg_dim
784        if(saddle_warn(iband)) then
785          write(tmpstr,'(i5)') band+iband-1
786          msg = TRIM(msg)//' '//TRIM(tmpstr)//','
787        end if
788      end do
789      write(tmpstr,'(6a)') ch10,'are not band extrema, but saddle points;',ch10, &
790 &                       'the transport equivalent formalism breaks down in these conditions.',ch10, &
791 &                       'The associated tensor(s) will therefore not be printed.'
792      msg = TRIM(msg)//TRIM(tmpstr)
793      if(io_unit==std_out) then
794        MSG_WARNING(msg)
795      else
796        write(io_unit,'(7a)') ch10,'--- !WARNING',ch10,TRIM(msg),ch10,'---'
797      end if
798    end if
799 
800    if(deg_dim>1 .and. mdim>1) then
801      write(msg,'(a)') ' Transport equivalent effective mass tensor'
802    else
803      write(msg,'(a)') ' Effective mass tensor'
804    end if
805 
806    if(mdim>1) then
807      write(format_eigvec,'(a,i1,a)') '(i3,',mdim,'f14.10,a,3f14.10)'
808    end if
809 
810    do iband=1,deg_dim
811      write(io_unit,'(2a,3(f6.3,a),i5)') ch10,' K-point (',kpt(1),',',kpt(2),',',kpt(3),') | band = ',band+iband-1
812      write(io_unit,'(a)') trim(msg)//':'
813      if(.not. saddle_warn(iband)) then
814        do adir=1,mdim
815          write(io_unit,'(3f26.10)') efmas_tensor(adir,:,iband)
816        end do
817        if(present(efmas_eigval))then
818          write(io_unit,'(a)') trim(msg)//' eigenvalues:'
819          write(io_unit,'(3f26.10)') efmas_eigval(:,iband)
820        endif
821      else
822        write(io_unit,'(a)') '     *** SADDLE POINT: TRANSPORT EQV. EFF. MASS NOT DEFINED (see WARNING above) ***'
823      end if
824 
825      if(mdim>1) then
826        if(mdim==2 .and. deg_dim>1) then
827          write(io_unit,'(a,f26.10)') 'Scaling of transport tensor (Eq. (FXX)) = ',transport_tensor_scale(iband)
828        end if
829        if(.not. saddle_warn(iband)) then
830          if(io_unit == std_out) then
831            write(io_unit,'(a)') trim(msg)//' eigenvectors in cartesian / reduced coord.:'
832            do adir=1,mdim
833              if( count( abs(efmas_eigval(adir,iband)-efmas_eigval(:,iband))<tol4 ) > 1 ) then
834                write(io_unit,'(i3,a)') adir, ' Eigenvalue degenerate => eigenvector undefined'
835              else
836                vec=zero; vec(1:mdim)=efmas_eigvec(adir,:,iband)
837                vec=matmul(transpose(rprimd)/two_pi,vec); vec=vec/sqrt(sum(vec**2))
838                write(io_unit,format_eigvec) adir, efmas_eigvec(adir,:,iband), ' / ', vec
839              end if
840            end do
841          end if
842        end if
843      end if
844 
845      if(mdim==3)then
846        !An exactly zero average effective masse is artificial (or from a saddle point with symmetry). Does not print.
847        if(abs(m_avg(iband))>tol8)then
848          write(io_unit,'(a,f14.10)') &
849 &          ' Angular average effective mass 1/(<1/m>)= ',m_avg(iband)
850        endif
851        if(abs(m_avg_frohlich(iband))>tol8)then
852          write(io_unit,'(a,f14.10)') &
853 &          ' Angular average effective mass for Frohlich model (<m**0.5>)**2= ',m_avg_frohlich(iband)
854        endif
855      endif
856 
857      write(io_unit,'(a)') ' Effective masses along directions: (cart. coord. / red. coord. -> eff. mass)'
858      do adir=1,ndirs
859        vec=dirs(:,adir)
860        vec=matmul(transpose(rprimd)/two_pi,vec); vec=vec/sqrt(sum(vec**2))
861        write(io_unit,'(i5,a,3f10.6,a,3f10.6,a,f14.10)') adir,': ', dirs(:,adir), ' / ', vec, ' -> ', m_cart(adir,iband)
862      end do
863    end do
864 
865    if(deg_dim>1 .and. mdim==3) then
866      write(io_unit,'(2a)') ch10,&
867 &     ' Angular average effective mass for Frohlich model, averaged over degenerate bands.'
868      write(io_unit,'(a,es16.6)') &
869 &     ' Value of     (<<m**0.5>>)**2 = ',(sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim)**2
870      write(io_unit,'(a,es16.6,a)') &
871 &     ' Absolute Value of <<m**0.5>> = ', sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim,ch10
872    endif
873 
874  end subroutine print_tr_efmas