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