TABLE OF CONTENTS
ABINIT/m_frohlichmodel [ Modules ]
NAME
m_frohlichmodel
FUNCTION
Compute ZPR, temperature-dependent electronic structure, and other properties using the Frohlich model
COPYRIGHT
Copyright (C) 2018-2022 ABINIT group (XG) 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_frohlichmodel 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_crystal 29 use m_ebands 30 use m_efmas_defs 31 use m_ifc 32 use m_dtset 33 34 use m_fstrings, only : sjoin, itoa 35 use m_gaussian_quadrature, only : cgqf 36 37 implicit none 38 39 private 40 41 public :: hamiltonian, frohlichmodel, polaronmass 42 43 contains
m_frohlichmodel/frohlichmodel [ Functions ]
[ Top ] [ m_frohlichmodel ] [ Functions ]
NAME
frohlichmodel
FUNCTION
Main routine to compute properties based on the Frohlich model
INPUTS
cryst<crystal_t>=Structure defining the unit cell dtset<dataset_type>=All input variables for this dataset. 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 ifc<ifc_type>=contains the dynamical matrix and the IFCs.
SOURCE
63 subroutine frohlichmodel(cryst, dtset, efmasdeg, efmasval, ifc) 64 65 !Arguments ------------------------------------ 66 !scalars 67 type(crystal_t),intent(in) :: cryst 68 type(dataset_type),intent(in) :: dtset 69 type(ifc_type),intent(in) :: ifc 70 !arrays 71 type(efmasdeg_type), intent(in) :: efmasdeg(:) 72 type(efmasval_type), intent(in) :: efmasval(:,:) 73 74 !Local variables ------------------------------ 75 !scalars 76 logical :: sign_warn 77 integer :: deg_dim,iband,ideg,idir,ikpt,imode,info,ipar,iphi,iqdir,itheta 78 integer :: jband,lwork,nphi,nqdir,ntheta 79 real(dp) :: angle_phi,cosph,costh,sinph,sinth,weight,weight_phi 80 real(dp) :: zpr_frohlich,zpr_q0_avg,zpr_q0_fact 81 !character(len=500) :: msg 82 !arrays 83 logical, allocatable :: saddle_warn(:), start_eigf3d_pos(:) 84 logical :: lutt_found(3), lutt_warn(3) 85 real(dp) :: kpt(3), lutt_params(3), lutt_unit_kdir(3,3) 86 real(dp), allocatable :: eigenval(:), rwork(:), unit_qdir(:,:) 87 real(dp), allocatable :: lutt_dij(:,:), lutt_eigenval(:,:) 88 real(dp), allocatable :: m_avg(:), m_avg_frohlich(:) 89 real(dp), allocatable :: gq_points_th(:),gq_weights_th(:) 90 real(dp), allocatable :: gq_points_cosph(:),gq_points_sinph(:) 91 real(dp), allocatable :: weight_qdir(:) 92 real(dp), allocatable :: polarity_qdir(:,:,:) 93 real(dp), allocatable :: proj_polarity_qdir(:,:) 94 real(dp), allocatable :: zpr_q0_phononfactor_qdir(:) 95 real(dp), allocatable :: frohlich_phononfactor_qdir(:) 96 real(dp), allocatable :: phfrq_qdir(:,:) 97 real(dp), allocatable :: dielt_qdir(:) 98 real(dp), allocatable :: dielt_avg(:) 99 real(dp), allocatable :: zpr_frohlich_avg(:) 100 complex(dpc), allocatable :: eigenvec(:,:), work(:) 101 complex(dpc), allocatable :: eig2_diag_cart(:,:,:,:) 102 complex(dpc), allocatable :: f3d(:,:) 103 104 !************************************************************************ 105 106 !!! Initialization of integrals 107 ntheta = dtset%efmas_ntheta 108 nphi = 2*ntheta 109 nqdir = nphi*ntheta 110 111 ABI_MALLOC(gq_points_th,(ntheta)) 112 ABI_MALLOC(gq_weights_th,(ntheta)) 113 ABI_MALLOC(gq_points_cosph,(nphi)) 114 ABI_MALLOC(gq_points_sinph,(nphi)) 115 116 ABI_MALLOC(unit_qdir,(3,nqdir)) 117 ABI_MALLOC(weight_qdir,(nqdir)) 118 119 call cgqf(ntheta,1,zero,zero,zero,pi,gq_points_th,gq_weights_th) 120 weight_phi=two*pi/real(nphi,dp) 121 do iphi=1,nphi 122 angle_phi=weight_phi*(iphi-1) 123 gq_points_cosph(iphi)=cos(angle_phi) 124 gq_points_sinph(iphi)=sin(angle_phi) 125 enddo 126 nqdir=0 127 do itheta=1,ntheta 128 costh=cos(gq_points_th(itheta)) 129 sinth=sin(gq_points_th(itheta)) 130 weight=gq_weights_th(itheta)*weight_phi*sinth 131 do iphi=1,nphi 132 cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi) 133 nqdir=nqdir+1 134 135 unit_qdir(1,nqdir)=sinth*cosph 136 unit_qdir(2,nqdir)=sinth*sinph 137 unit_qdir(3,nqdir)=costh 138 weight_qdir(nqdir)=weight 139 140 enddo 141 enddo 142 143 ABI_FREE(gq_points_th) 144 ABI_FREE(gq_weights_th) 145 ABI_FREE(gq_points_cosph) 146 ABI_FREE(gq_points_sinph) 147 148 ABI_MALLOC(polarity_qdir,(3,3*cryst%natom,nqdir)) 149 ABI_MALLOC(proj_polarity_qdir,(3*cryst%natom,nqdir)) 150 ABI_MALLOC(zpr_q0_phononfactor_qdir,(nqdir)) 151 ABI_MALLOC(frohlich_phononfactor_qdir,(nqdir)) 152 ABI_MALLOC(phfrq_qdir,(3*cryst%natom,nqdir)) 153 ABI_MALLOC(dielt_qdir,(nqdir)) 154 ABI_MALLOC(dielt_avg,(3*cryst%natom)) 155 156 !Compute phonon frequencies and mode-polarity for each qdir 157 call ifc%calcnwrite_nana_terms(cryst, nqdir, unit_qdir, phfrq2l=phfrq_qdir, polarity2l=polarity_qdir) 158 159 !Compute dielectric tensor for each qdir 160 do iqdir=1,nqdir 161 dielt_qdir(iqdir)=DOT_PRODUCT(unit_qdir(:,iqdir),MATMUL(ifc%dielt(:,:),unit_qdir(:,iqdir))) 162 enddo 163 164 !Compute projection of mode-polarity on qdir, and other derived quantities summed over phonon branches for each iqdir. 165 !Note that acoustic modes are discarded (imode sum starts only from 4) 166 zpr_q0_phononfactor_qdir=zero 167 zpr_q0_avg=zero 168 frohlich_phononfactor_qdir=zero 169 dielt_avg=zero 170 171 do iqdir=1,nqdir 172 do imode=4,3*cryst%natom 173 proj_polarity_qdir(imode,iqdir)=DOT_PRODUCT(unit_qdir(:,iqdir),polarity_qdir(:,imode,iqdir)) 174 zpr_q0_phononfactor_qdir(iqdir)=zpr_q0_phononfactor_qdir(iqdir)+& 175 & proj_polarity_qdir(imode,iqdir)**2 / phfrq_qdir(imode,iqdir) **2 176 frohlich_phononfactor_qdir(iqdir)=frohlich_phononfactor_qdir(iqdir)+& 177 & proj_polarity_qdir(imode,iqdir)**2 / phfrq_qdir(imode,iqdir) **(three*half) 178 dielt_avg(imode)=dielt_avg(imode)+& 179 & weight_qdir(iqdir)*frohlich_phononfactor_qdir(iqdir)/dielt_qdir(iqdir)**2 180 enddo 181 zpr_q0_avg=zpr_q0_avg+& 182 & weight_qdir(iqdir)*zpr_q0_phononfactor_qdir(iqdir)/dielt_qdir(iqdir)**2 183 enddo 184 dielt_avg=dielt_avg*two**(-half)*cryst%ucvol**(-one) 185 zpr_q0_avg=zpr_q0_avg*quarter*piinv 186 zpr_q0_fact=zpr_q0_avg*eight*pi*(three*quarter*piinv)**third*cryst%ucvol**(-four*third) 187 188 write(ab_out,'(a)')'--------------------------------------------------------------------------------' 189 write(ab_out,'(a)')' Dielectric average (EQ. 25 Melo2022)' 190 write(ab_out,'(a)')'--------------------------------------------------------------------------------' 191 write(ab_out,'(a)')' Mode <1/epsilon*SQRT(w_LO/2)> Cumulative sum' 192 do imode=4,3*cryst%natom 193 write(ab_out,'(i5,f28.12,f28.12)') & 194 !Reversed cummulative sum to avoid rewriting the spherical intergration from above 195 & imode,dielt_avg(imode)-dielt_avg(abs(mod(imode-1,3*cryst%natom))),dielt_avg(imode) 196 enddo 197 write(ab_out,'(a)')'--------------------------------------------------------------------------------' 198 199 !DEBUG 200 ! do iqdir=1,nqdir,513 201 ! write(std_out,'(a,3f8.4,3es12.4)')' unit_qdir,dielt_qdir,zpr_q0_phononfactor_qdir,frohlich_phononfactor=',& 202 !& unit_qdir(:,iqdir),dielt_qdir(iqdir),zpr_q0_phononfactor_qdir(iqdir),frohlich_phononfactor_qdir(iqdir) 203 ! do imode=1,3*cryst%natom 204 ! write(std_out,'(a,i5,6es12.4)')' imode,phfrq_qdir,phfrq(cmm1),polarity_qdir=',& 205 !& imode,phfrq_qdir(imode,iqdir),phfrq_qdir(imode,iqdir)*Ha_cmm1,polarity_qdir(:,imode,iqdir),proj_polarity_qdir(imode,iqdir) 206 ! enddo 207 ! enddo 208 ! write(std_out,'(2a,3es12.4)')ch10,& 209 !& ' zpr_q0_avg, zpr_q0_fact, zpr_q0_fact (eV) =',zpr_q0_avg, zpr_q0_fact, zpr_q0_fact*Ha_eV 210 !ENDDEBUG 211 212 write(ab_out,'(6a,f14.6,a,f14.6,a)') ch10,& 213 & ' Rough correction to the ZPR, to take into account the missing q=0 piece using Frohlich model:',ch10,& 214 & ' (+ for occupied states, - for unoccupied states) * zpr_q0_fact / (Nqpt_full_bz)**(1/3) ',ch10,& 215 & ' where Nqpt_full_bz=number of q wavevectors in full BZ, and zpr_q0_fact=',zpr_q0_fact,' Ha=',zpr_q0_fact*Ha_eV,' eV' 216 217 !Compute effective masses, and integrate the Frohlich model 218 do ikpt=1,dtset%nkpt 219 220 kpt(:)=dtset%kptns(:,ikpt) 221 do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2) 222 223 deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1 224 225 ABI_MALLOC(eig2_diag_cart,(3,3,deg_dim,deg_dim)) 226 227 !Convert eig2_diag to cartesian coordinates 228 do iband=1,deg_dim 229 do jband=1,deg_dim 230 eig2_diag_cart(:,:,iband,jband)=efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband) 231 eig2_diag_cart(:,:,iband,jband)=& 232 & matmul(matmul(cryst%rprimd,eig2_diag_cart(:,:,iband,jband)),transpose(cryst%rprimd))/two_pi**2 233 enddo 234 enddo 235 236 ABI_MALLOC(f3d,(deg_dim,deg_dim)) 237 ABI_MALLOC(m_avg,(deg_dim)) 238 ABI_MALLOC(m_avg_frohlich,(deg_dim)) 239 ABI_MALLOC(zpr_frohlich_avg,(deg_dim)) 240 ABI_MALLOC(eigenval,(deg_dim)) 241 ABI_MALLOC(saddle_warn,(deg_dim)) 242 ABI_MALLOC(start_eigf3d_pos,(deg_dim)) 243 244 m_avg=zero 245 m_avg_frohlich=zero 246 saddle_warn=.false. 247 248 !Initializations for the diagonalization routine 249 if(deg_dim>1)then 250 251 ABI_MALLOC(eigenvec,(deg_dim,deg_dim)) 252 lwork=-1 253 ABI_MALLOC(rwork,(3*deg_dim-2)) 254 ABI_MALLOC(work,(1)) 255 call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info) 256 lwork=int(work(1)) 257 ABI_FREE(work) 258 ABI_MALLOC(work,(lwork)) 259 260 endif 261 262 !Compute the Luttinger parameters for the cubic case (deg_dim=3) 263 if(deg_dim==3) then 264 265 ABI_MALLOC(lutt_eigenval, (3,deg_dim)) 266 ABI_MALLOC(lutt_dij, (deg_dim, deg_dim)) 267 268 !Define unit_kdir for Luttinger parameters 269 lutt_unit_kdir(:,1) = (/1,0,0/) 270 lutt_unit_kdir(:,2) = 1/sqrt(2.0)*(/1,1,0/) 271 lutt_unit_kdir(:,3) = 1/sqrt(3.0)*(/1,1,1/) 272 273 !Degeneracy problems warning 274 lutt_warn=(/.false.,.false.,.false./) 275 276 !Inverse effective mass tensor eigenvalues in lutt_unit_kdir directions 277 do idir=1,3 278 do iband=1,deg_dim 279 do jband=1,deg_dim 280 lutt_dij(iband,jband)=& 281 & DOT_PRODUCT(lutt_unit_kdir(:,idir),MATMUL(eig2_diag_cart(:,:,iband,jband),lutt_unit_kdir(:,idir))) 282 enddo 283 enddo 284 285 eigenvec=lutt_dij ; lutt_eigenval(idir,:)=zero 286 work=zero ; rwork=zero 287 call zheev('V','U',deg_dim,eigenvec,deg_dim,lutt_eigenval(idir,:),work,lwork,rwork,info) 288 ABI_CHECK(info == 0, sjoin("zheev returned info:", itoa(info))) 289 enddo 290 291 !Check degeneracies in (100) direction, and evaluate A and B. 292 !Eigenvalues are 2*A (d=1), 2*B (d=2) 293 if(abs(lutt_eigenval(1,2)-lutt_eigenval(1,3))<tol5) then 294 lutt_params(2)=0.5*((lutt_eigenval(1,2)+lutt_eigenval(1,3))/2) 295 lutt_params(1)=0.5*lutt_eigenval(1,1) 296 else if(abs(lutt_eigenval(1,2)-lutt_eigenval(1,1))<tol5) then 297 lutt_params(2)=0.5*((lutt_eigenval(1,2)+lutt_eigenval(1,1))/2) 298 lutt_params(1)=0.5*lutt_eigenval(1,3) 299 else 300 lutt_warn(1)=.true. 301 endif 302 303 !Check degeneracies in (111) direction and evaluate C 304 !Eigenvalues are 2/3*(A+2B-C) (d=2), 2/3*(A+2B+2C) (d=1) 305 if(abs(lutt_eigenval(3,2)-lutt_eigenval(3,3))<tol5) then 306 lutt_params(3)=lutt_params(1)+2*lutt_params(2)-1.5*(0.5*(lutt_eigenval(3,2)+lutt_eigenval(3,3))) 307 else if(abs(lutt_eigenval(3,2)-lutt_eigenval(3,1))<tol5) then 308 lutt_params(3)=lutt_params(1)+2*lutt_params(2)-1.5*(0.5*(lutt_eigenval(3,2)+lutt_eigenval(3,1))) 309 else 310 lutt_warn(2)=.true. 311 endif 312 313 !Verify that the (110) direction eigenvalues are coherent with Luttinger parameters 314 !Eigenvalues are 2B, A+B-C, A+B+C 315 lutt_found=(/.false.,.false.,.false./) 316 do ipar=1,deg_dim 317 if(abs(lutt_eigenval(2,ipar)-2*lutt_params(2))<tol4) then 318 lutt_found(1)=.true. 319 else if(abs(lutt_eigenval(2,ipar)-(lutt_params(1)+lutt_params(2)-lutt_params(3)))<tol4) then 320 lutt_found(2)=.true. 321 else if(abs(lutt_eigenval(2,ipar)-(lutt_params(1)+lutt_params(2)+lutt_params(3)))<tol4) then 322 lutt_found(3)=.true. 323 endif 324 enddo 325 326 if(.not. (all(lutt_found))) then 327 lutt_warn(3)=.true. 328 endif 329 330 ABI_FREE(lutt_eigenval) 331 ABI_FREE(lutt_dij) 332 333 endif !Luttinger parameters 334 335 !Perform the integral over the sphere 336 zpr_frohlich_avg=zero 337 do iqdir=1,nqdir 338 do iband=1,deg_dim 339 do jband=1,deg_dim 340 f3d(iband,jband)=DOT_PRODUCT(unit_qdir(:,iqdir),MATMUL(eig2_diag_cart(:,:,iband,jband),unit_qdir(:,iqdir))) 341 enddo 342 enddo 343 344 if(deg_dim==1)then 345 eigenval(1)=f3d(1,1) 346 else 347 eigenvec = f3d ; eigenval = zero 348 work=zero ; rwork=zero 349 call zheev('V','U',deg_dim,eigenvec,deg_dim,eigenval,work,lwork,rwork,info) 350 ABI_CHECK(info == 0, sjoin("zheev returned info:", itoa(info))) 351 endif 352 353 m_avg = m_avg + weight_qdir(iqdir)*eigenval 354 m_avg_frohlich = m_avg_frohlich + weight_qdir(iqdir)/(abs(eigenval))**half 355 zpr_frohlich_avg = zpr_frohlich_avg + & 356 & weight_qdir(iqdir) * frohlich_phononfactor_qdir(iqdir)/((abs(eigenval))**half *dielt_qdir(iqdir)**2) 357 358 if(iqdir==1) start_eigf3d_pos = (eigenval > 0) 359 360 do iband=1,deg_dim 361 if(start_eigf3d_pos(iband) .neqv. (eigenval(iband)>0)) then 362 saddle_warn(iband)=.true. 363 end if 364 end do 365 366 enddo 367 368 if(deg_dim>1)then 369 ABI_FREE(eigenvec) 370 ABI_FREE(rwork) 371 ABI_FREE(work) 372 endif 373 374 m_avg = quarter*piinv*m_avg 375 m_avg = one/m_avg 376 377 m_avg_frohlich = quarter*piinv * m_avg_frohlich 378 m_avg_frohlich = m_avg_frohlich**2 379 380 zpr_frohlich_avg = quarter*piinv * zpr_frohlich_avg 381 382 if(deg_dim==1)then 383 write(ab_out,'(2a,3(f6.3,a),i5)')ch10,& 384 & ' - At k-point (',kpt(1),',',kpt(2),',',kpt(3),'), band ',& 385 & efmasdeg(ikpt)%degs_bounds(1,ideg) 386 else 387 write(ab_out,'(2a,3(f6.3,a),i5,a,i5)')ch10,& 388 & ' - At k-point (',kpt(1),',',kpt(2),',',kpt(3),'), bands ',& 389 & efmasdeg(ikpt)%degs_bounds(1,ideg),' through ',efmasdeg(ikpt)%degs_bounds(2,ideg) 390 endif 391 392 !Print the Luttinger for the cubic case (deg_dim=3) 393 if(deg_dim==3) then 394 if (.not. (any(saddle_warn))) then 395 if(any(lutt_warn)) then 396 ! Warn for degeneracy breaking in inverse effective mass tensor eigenvalues 397 write(ab_out, '(2a)') ch10, ' Luttinger parameters could not be determined:' 398 if (lutt_warn(1)) then 399 write(ab_out, '(a)') ' Predicted degeneracies for deg_dim = 3 are not met for (100) direction.' 400 endif 401 if (lutt_warn(2)) then 402 write(ab_out, '(a)') ' Predicted degeneracies for deg_dim = 3 are not met for (111) direction.' 403 endif 404 if (lutt_warn(3)) then 405 write(ab_out, '(a)') ' Predicted inverse effective mass tensor eigenvalues for direction (110) are not met.' 406 endif 407 write(ab_out, '(a)') ch10 408 else 409 write(ab_out, '(a,3f14.6)') ' Luttinger parameters (A, B, C) [at. units]: ',lutt_params(:) 410 endif 411 endif 412 endif 413 414 sign_warn=.false. 415 do iband=1,deg_dim 416 if(saddle_warn(iband)) then 417 write(ab_out,'(a,i5,a)') ' Band ',efmasdeg(ikpt)%degs_bounds(1,ideg)+iband-1,& 418 & ' SADDLE POINT - Frohlich effective mass and ZPR cannot be defined. ' 419 sign_warn=.true. 420 else 421 m_avg_frohlich(iband) = DSIGN(m_avg_frohlich(iband),m_avg(iband)) 422 zpr_frohlich_avg(iband) = -DSIGN(zpr_frohlich_avg(iband),m_avg(iband)) 423 write(ab_out,'(a,i5,a,f14.10)') & 424 & ' Band ',efmasdeg(ikpt)%degs_bounds(1,ideg)+iband-1,& 425 & ' Angular average effective mass for Frohlich model (<m**0.5>)**2= ',m_avg_frohlich(iband) 426 endif 427 if(start_eigf3d_pos(iband) .neqv. start_eigf3d_pos(1))then 428 sign_warn=.true. 429 endif 430 enddo 431 432 if(sign_warn .eqv. .false.)then 433 zpr_frohlich = four*pi* two**(-half) * (sum(zpr_frohlich_avg(1:deg_dim))/deg_dim) / cryst%ucvol 434 write(ab_out,'(2a)')& 435 & ' Angular and band average effective mass and ZPR for Frohlich model.' 436 write(ab_out,'(a,es16.6)') & 437 & ' Value of (<<m**0.5>>)**2 = ',(sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim)**2 438 write(ab_out,'(a,es16.6)') & 439 & ' Absolute Value of <<m**0.5>> = ', sum(abs(m_avg_frohlich(1:deg_dim))**0.5)/deg_dim 440 write(ab_out,'(a,es16.6,a,es16.6,a)') & 441 & ' ZPR from Frohlich model = ',zpr_frohlich,' Ha=',zpr_frohlich*Ha_eV,' eV' 442 else 443 write(ab_out,'(a)')& 444 & ' Angular and band average effective mass for Frohlich model cannot be defined because of a sign problem.' 445 endif 446 447 ABI_FREE(eig2_diag_cart) 448 ABI_FREE(f3d) 449 ABI_FREE(m_avg) 450 ABI_FREE(m_avg_frohlich) 451 ABI_FREE(zpr_frohlich_avg) 452 ABI_FREE(eigenval) 453 ABI_FREE(saddle_warn) 454 ABI_FREE(start_eigf3d_pos) 455 456 enddo ! ideg 457 enddo ! ikpt 458 459 ABI_FREE(unit_qdir) 460 ABI_FREE(weight_qdir) 461 ABI_FREE(polarity_qdir) 462 ABI_FREE(proj_polarity_qdir) 463 ABI_FREE(phfrq_qdir) 464 ABI_FREE(dielt_qdir) 465 ABI_FREE(dielt_avg) 466 ABI_FREE(zpr_q0_phononfactor_qdir) 467 ABI_FREE(frohlich_phononfactor_qdir) 468 469 end subroutine frohlichmodel
m_frohlichmodel/polaronmass [ Functions ]
[ Top ] [ m_frohlichmodel ] [ Functions ]
NAME
polaronmass
FUNCTION
Improved routine to compute properties based on the Frohlich model, including effective masses in the cubic case.
INPUTS
cryst<crystal_t>=Structure defining the unit cell dtset<dataset_type>=All input variables for this dataset. 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 ifc<ifc_type>=contains the dynamical matrix and the IFCs.
SOURCE
490 subroutine polaronmass(cryst, dtset, efmasdeg, efmasval, ifc) 491 492 !Arguments ------------------------------------ 493 !scalars 494 type(crystal_t),intent(in) :: cryst 495 type(dataset_type),intent(in) :: dtset 496 type(ifc_type),intent(in) :: ifc 497 !arrays 498 type(efmasdeg_type), intent(in) :: efmasdeg(:) 499 type(efmasval_type), intent(in) :: efmasval(:,:) 500 501 !Local variables ------------------------------ 502 !scalars 503 integer :: deg_dim, signpm 504 integer :: i, iband, jband, ideg, idir, iqdir, ieig 505 integer :: ikpt, ixi, ipar, iphi, iphon, imode, itheta, ik 506 integer :: nkdir, nqdir, ntheta, nphi, nxi, nkgrid 507 integer :: info, lwork 508 real(dp) :: angle_phi,cosph,costh,sinph,sinth,weight,weight_phi 509 real(dp) :: qpt, krange, nq_factor 510 !character(len=500) :: msg 511 !arrays 512 !Electronic 513 real(dp), allocatable :: eigenvec(:,:,:,:), eigenval(:,:,:) 514 real(dp), allocatable :: rwork(:), unit_qdir(:,:) 515 complex(dpc), allocatable :: leigenvec(:,:), work(:) 516 complex(dpc), allocatable :: eig2_diag_cart(:,:,:,:) 517 !Luttinger 518 logical :: lutt_found(3), lutt_warn(3) 519 real(dp) :: lutt_params(3) 520 real(dp) :: kpoint(3) 521 real(dp), allocatable :: lutt_dij(:,:), lutt_eigenval(:,:), leigenval(:) 522 !Dielectric 523 real(dp), allocatable :: dielt_qdir(:) 524 !Phonons 525 real(dp), allocatable :: gq_points_th(:),gq_weights_th(:) 526 real(dp), allocatable :: gq_points_cosph(:),gq_points_sinph(:) 527 real(dp), allocatable :: weight_qdir(:) 528 real(dp), allocatable :: polarity_qdir(:,:,:) 529 real(dp), allocatable :: phfrq_qdir(:,:) 530 !Self-energy and polaron mass 531 real(dp) :: xi 532 real(dp) :: temporary1(3,3), temporary2(3), temporary3 533 real(dp) :: unitary_33(3,3) 534 real(dp) :: minelecmass,eham(3,3) 535 real(dp) :: kpt(3), k_plus_q(3) 536 real(dp), allocatable :: lutt_unit_kdir(:,:) 537 real(dp), allocatable :: omega_zero(:) 538 real(dp), allocatable :: intsum(:,:,:,:) 539 real(dp), allocatable :: sigma(:,:,:), d2sigmadk2(:,:) 540 real(dp), allocatable :: invepsilonstar(:) 541 real(dp), allocatable :: invemass(:,:), invemass_ieig(:,:),invpolmass(:,:) 542 543 !************************************************************************ 544 545 !Define Luttinger and Phonon integration parameters 546 !Based solely solely one value - ntheta 547 ntheta = dtset%efmas_ntheta 548 nphi = 2*ntheta 549 nqdir = nphi*ntheta 550 551 !Define constants 552 ! 3x3 Unitary matrix 553 unitary_33 = 0.0_dp 554 do i=1,3 555 unitary_33(i,i) = 1.0_dp 556 enddo 557 558 !Define unit_kdir for Luttinger parameters 559 nkdir=3 560 ABI_MALLOC(lutt_unit_kdir,(3,nkdir)) 561 lutt_unit_kdir(:,1) = (/1,0,0/) 562 lutt_unit_kdir(:,2) = 1/sqrt(2.0)*(/1,1,0/) 563 lutt_unit_kdir(:,3) = 1/sqrt(3.0)*(/1,1,1/) 564 !These are for testing purpose only 565 ! lutt_unit_kdir(:,4) = (/0,1,0/) 566 ! lutt_unit_kdir(:,5) = (/0,0,1/) 567 568 !Compute effective masses, and integrate the Frohlich model 569 do ikpt=1,dtset%nkpt 570 571 kpt(:)=dtset%kptns(:,ikpt) 572 do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2) 573 574 deg_dim = efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1 575 576 ABI_MALLOC(eig2_diag_cart,(3,3,deg_dim,deg_dim)) 577 578 !Convert eig2_diag to cartesian coordinates 579 do iband=1,deg_dim 580 do jband=1,deg_dim 581 eig2_diag_cart(:,:,iband,jband)=efmasval(ideg,ikpt)%eig2_diag(:,:,iband,jband) 582 eig2_diag_cart(:,:,iband,jband)=& 583 & matmul(matmul(cryst%rprimd,eig2_diag_cart(:,:,iband,jband)),transpose(cryst%rprimd))/two_pi**2 584 enddo ! jband 585 enddo ! iband 586 587 ABI_MALLOC(leigenval,(deg_dim)) 588 589 !Initializations for the diagonalization routine 590 if(deg_dim>1)then 591 592 ABI_MALLOC(leigenvec,(deg_dim,deg_dim)) 593 lwork=-1 594 ABI_MALLOC(rwork,(3*deg_dim-2)) 595 ABI_MALLOC(work,(1)) 596 call zheev('V','U',deg_dim,leigenvec,deg_dim,leigenval,work,lwork,rwork,info) 597 lwork=int(work(1)) 598 ABI_FREE(work) 599 ABI_MALLOC(work,(lwork)) 600 601 endif 602 603 !Compute the Luttinger parameters for the cubic case (deg_dim=3) 604 if(deg_dim==3) then 605 606 ABI_MALLOC(lutt_eigenval, (3,deg_dim)) 607 ABI_MALLOC(lutt_dij, (deg_dim, deg_dim)) 608 609 !Degeneracy problems warning 610 lutt_warn=(/.false.,.false.,.false./) 611 612 !Inverse effective mass tensor eigenvalues in lutt_unit_kdir directions 613 do idir=1,3 614 do iband=1,deg_dim 615 do jband=1,deg_dim 616 lutt_dij(iband,jband)=& 617 & DOT_PRODUCT(lutt_unit_kdir(:,idir),MATMUL(eig2_diag_cart(:,:,iband,jband),lutt_unit_kdir(:,idir))) 618 enddo 619 enddo 620 621 leigenvec=lutt_dij ; lutt_eigenval(idir,:)=zero 622 work=zero ; rwork=zero 623 call zheev('V','U',deg_dim,leigenvec,deg_dim,lutt_eigenval(idir,:),work,lwork,rwork,info) 624 ABI_CHECK(info == 0, sjoin("zheev returned info:", itoa(info))) 625 enddo 626 627 !Check degeneracies in (100) direction, and evaluate A and B. 628 !Eigenvalues are 2*A (d=1), 2*B (d=2) 629 if(abs(lutt_eigenval(1,2)-lutt_eigenval(1,3))<tol3) then 630 lutt_params(2)=0.5*((lutt_eigenval(1,2)+lutt_eigenval(1,3))/2) 631 lutt_params(1)=0.5*lutt_eigenval(1,1) 632 else if(abs(lutt_eigenval(1,2)-lutt_eigenval(1,1))<tol3) then 633 lutt_params(2)=0.5*((lutt_eigenval(1,2)+lutt_eigenval(1,1))/2) 634 lutt_params(1)=0.5*lutt_eigenval(1,3) 635 else 636 lutt_warn(1)=.true. 637 endif 638 639 !Check degeneracies in (111) direction and evaluate C 640 !Eigenvalues are 2/3*(A+2B-C) (d=2), 2/3*(A+2B+2C) (d=1) 641 if(abs(lutt_eigenval(3,2)-lutt_eigenval(3,3))<tol3) then 642 lutt_params(3)=lutt_params(1)+2*lutt_params(2)-1.5*(0.5*(lutt_eigenval(3,2)+lutt_eigenval(3,3))) 643 else if(abs(lutt_eigenval(3,2)-lutt_eigenval(3,1))<tol3) then 644 lutt_params(3)=lutt_params(1)+2*lutt_params(2)-1.5*(0.5*(lutt_eigenval(3,2)+lutt_eigenval(3,1))) 645 else 646 lutt_warn(2)=.true. 647 endif 648 649 !Verify that the (110) direction eigenvalues are coherent with Luttinger parameters 650 !Eigenvalues are 2B, A+B-C, A+B+C 651 lutt_found=(/.false.,.false.,.false./) 652 do ipar=1,deg_dim 653 if(abs(lutt_eigenval(2,ipar)-2*lutt_params(2))<tol3) then 654 lutt_found(1)=.true. 655 else if(abs(lutt_eigenval(2,ipar)-(lutt_params(1)+lutt_params(2)-lutt_params(3)))<tol3) then 656 lutt_found(2)=.true. 657 else if(abs(lutt_eigenval(2,ipar)-(lutt_params(1)+lutt_params(2)+lutt_params(3)))<tol3) then 658 lutt_found(3)=.true. 659 endif 660 enddo 661 662 if(.not. (all(lutt_found))) then 663 lutt_warn(3)=.true. 664 endif 665 666 ABI_FREE(lutt_eigenval) 667 ABI_FREE(lutt_dij) 668 669 endif !Luttinger parameters 670 671 if(deg_dim>1)then 672 ABI_FREE(leigenvec) 673 ABI_FREE(rwork) 674 ABI_FREE(work) 675 endif 676 677 !Print the Luttinger for the cubic case (deg_dim=3) 678 if(deg_dim==3) then 679 if(any(lutt_warn)) then 680 ! Warn for degeneracy breaking in inverse effective mass tensor eigenvalues 681 write(ab_out, '(2a)') ch10, ' Luttinger parameters could not be determined:' 682 if (lutt_warn(1)) then 683 write(ab_out, '(a)') ' Predicted degeneracies for deg_dim = 3 are not met for (100) direction.' 684 endif 685 if (lutt_warn(2)) then 686 write(ab_out, '(a)') ' Predicted degeneracies for deg_dim = 3 are not met for (111) direction.' 687 endif 688 if (lutt_warn(3)) then 689 write(ab_out, '(a)') ' Predicted inverse effective mass tensor eigenvalues for direction (110) are not met.' 690 endif 691 write(ab_out, '(a)') ch10 692 else 693 write(ab_out, '(a,3f14.6)') ' Luttinger parameters (A, B, C) (a.u.): ',lutt_params(:) 694 endif 695 endif 696 697 ABI_FREE(eig2_diag_cart) 698 ABI_FREE(leigenval) 699 700 enddo ! ideg 701 enddo ! ikpt 702 703 !Build inverse electronic effective mass for different directions from Luttinger params 704 deg_dim = 3 705 ABI_MALLOC(invemass,(deg_dim,nkdir)) 706 ABI_MALLOC(invemass_ieig,(deg_dim,nkdir)) 707 708 ! 100 -direction 709 invemass(1,1) = two*lutt_params(1) ! 2A 710 invemass(2,1) = two*lutt_params(2) ! 2B 711 invemass(3,1) = two*lutt_params(2) ! 2B 712 ! 110 -direction 713 invemass(1,2) = lutt_params(1) + lutt_params(2) + lutt_params(3) ! A + B + C 714 invemass(2,2) = MIN(two*lutt_params(2), lutt_params(1) + lutt_params(2) - lutt_params(3) ) ! 2B 715 invemass(3,2) = MAX(two*lutt_params(2), lutt_params(1) + lutt_params(2) - lutt_params(3) ) ! A + B - C 716 ! 111 -direction 717 invemass(1,3) = two*( lutt_params(1) + two*lutt_params(2) + two*lutt_params(3) ) / three ! 2(A + 2B + 2C)/3 718 invemass(2,3) = two*( lutt_params(1) + two*lutt_params(2) - lutt_params(3) ) / three ! 2(A + 2B - 2C)/3 719 invemass(3,3) = two*( lutt_params(1) + two*lutt_params(2) - lutt_params(3) ) / three ! 2(A + 2B - 2C)/3 720 721 !END Diagonalize 3x3 Luttinger-Kohn Hamiltonian 722 723 ABI_MALLOC(unit_qdir,(3,nqdir)) 724 ABI_MALLOC(polarity_qdir,(3,3*cryst%natom,nqdir)) 725 ABI_MALLOC(phfrq_qdir,(3*cryst%natom,nqdir)) 726 ABI_MALLOC(dielt_qdir,(nqdir)) 727 728 ABI_MALLOC(gq_points_th,(ntheta)) 729 ABI_MALLOC(gq_weights_th,(ntheta)) 730 ABI_MALLOC(gq_points_cosph,(nphi)) 731 ABI_MALLOC(gq_points_sinph,(nphi)) 732 ABI_MALLOC(weight_qdir,(nqdir)) 733 734 call cgqf(ntheta,1,zero,zero,-one,one,gq_points_th,gq_weights_th) 735 weight_phi=two*pi/real(nphi,dp) 736 do iphi=1,nphi 737 angle_phi=weight_phi*(iphi-1) 738 gq_points_cosph(iphi)=cos(angle_phi) 739 gq_points_sinph(iphi)=sin(angle_phi) 740 enddo 741 nqdir=0 742 do itheta=1,ntheta 743 costh=gq_points_th(itheta) 744 sinth=sqrt(one-costh**2) 745 weight=gq_weights_th(itheta)*weight_phi 746 do iphi=1,nphi 747 cosph=gq_points_cosph(iphi) ; sinph=gq_points_sinph(iphi) 748 nqdir=nqdir+1 749 750 unit_qdir(1,nqdir)=sinth*cosph 751 unit_qdir(2,nqdir)=sinth*sinph 752 unit_qdir(3,nqdir)=costh 753 weight_qdir(nqdir)=weight 754 755 enddo 756 enddo 757 758 ABI_FREE(gq_points_th) 759 ABI_FREE(gq_weights_th) 760 ABI_FREE(gq_points_cosph) 761 ABI_FREE(gq_points_sinph) 762 763 !In the following, compute invepsilonstar(imode) benefiting from the cubic symmetry. 764 !Retrieve IR active phonon frequencies 765 !Compute phonon frequencies and mode-polarity for each qdir (actually, only the first would suffice) 766 call ifc%calcnwrite_nana_terms(cryst, nqdir, unit_qdir, phfrq2l=phfrq_qdir, polarity2l=polarity_qdir) 767 768 !Calculate inverse epsilon* for each optical phonon mode 769 ! (epsilon*)**(-1) = 4Pi/Omega_0 (p_j0/(dielt_qdir*omeja_j0))**2 770 ABI_MALLOC(invepsilonstar,(3*cryst%natom)) 771 ABI_MALLOC(omega_zero,(3*cryst%natom)) 772 invepsilonstar = zero 773 774 do imode = 1,3*cryst%natom 775 !For ease of treatment loop over all phonon branches but... 776 !Avoid the acoustic branches 777 if(imode > 3) then 778 invepsilonstar(imode) = & 779 four*pi/cryst%ucvol*(dot_product(unit_qdir(:,1),polarity_qdir(:,imode,1)) & 780 /( ifc%dielt(1,1)*phfrq_qdir(imode,1)) )**two 781 endif 782 enddo ! imode 783 784 !Set nkgrid and krange 785 !For the time being we need 3 points for the finite difference 786 nkgrid = 3 787 !Set material dependent length scale for the finite difference 788 !Lowest optical phonon frequency to be used 789 minelecmass=1.0_dp/maxval(abs(invemass)) 790 krange = sqrt(two*minelecmass*phfrq_qdir(4,1)/1000.0) 791 792 !Diagonalize 3x3 Luttinger-Kohn Hamiltonian 793 794 !Initializations for the diagonalization routine 795 ABI_MALLOC(eigenval,(deg_dim,nkgrid,nkdir)) 796 ABI_MALLOC(eigenvec,(deg_dim,deg_dim,nkgrid,nkdir)) 797 lwork=-1 798 ABI_MALLOC(rwork,(3*deg_dim-2)) 799 ABI_MALLOC(work,(1)) 800 801 !Initialize eigenval 802 eigenval = zero 803 do idir = 1,nkdir 804 do ik = 1, nkgrid 805 kpoint(:) = (ik -1.0_dp)*krange*lutt_unit_kdir(:,idir) 806 eham = hamiltonian(lutt_params, kpoint) 807 call dsyev('V','U',3,eham,3,eigenval(1:3,ik,idir),work,lwork,info) 808 eigenvec(:,:,ik,idir) = eham 809 lwork=int(work(1)) 810 ABI_FREE(work) 811 ABI_MALLOC(work,(lwork)) 812 enddo 813 enddo 814 815 ABI_FREE(rwork) 816 ABI_FREE(work) 817 818 ABI_MALLOC(intsum,(deg_dim,deg_dim,nkgrid,nkdir)) 819 ABI_MALLOC(sigma,(nkgrid,deg_dim,nkdir)) 820 821 !main loop 822 !Dummy value for omega_zero 823 omega_zero = phfrq_qdir(:,1) 824 !Establish VB or CB sign 825 signpm = 1 826 if(lutt_params(1) .LT. zero) then 827 signpm = -1 828 endif 829 830 !Define Luttinger and Phonon integration parameters 831 !One value input - effmass_ntheta input variable 832 nxi = ntheta 833 834 !Normalization factor for the integral 835 nq_factor = nxi*two/pi*(four*pi) 836 837 !Initialize self-energy 838 sigma = zero 839 840 !Summation over IR active phonon modes 841 do iphon = 1, 3*cryst%natom 842 if( invepsilonstar(iphon) > 1E-10 ) then 843 !Summation over relevant directions (100,110,111) 844 do idir=1,nkdir 845 !Summation over electronic eigenvalues 846 do ieig=1,3 847 !Summation over k-range 848 !Here consider a small BZ region around Gamma 849 !A convergene study might be performed around this value 850 do ik = 1,nkgrid 851 !Define k-point vector around which one integrates 852 !Needs at least 2 considering TRS and finite difference employed later 853 !A three point formula is implemented later. 854 kpoint(:) = ( ik - 1.0_dp) * krange*lutt_unit_kdir(:,idir) 855 !Perform hyperbolic tangent integration for the semi-infinite q domain 856 !Use a mapping to a tangent function - faster convergence wrt qpt sampling 857 do ixi = 0,nxi 858 xi = ixi*pi/( two * nxi) 859 if ( ixi .EQ. nxi ) xi = xi - tol8 860 !Wave-vector length 861 qpt = ( omega_zero(iphon)/abs(lutt_params(1)) )**half*tan(xi) 862 ! XG 20211106 : the best integration scheme for a finite interval is based on Gauss-Legendre approach, 863 ! which was set up previously with nqdir point. So replaced the itheta and iphi loop by the loop over nqdir 864 do iqdir=1,nqdir 865 k_plus_q = kpoint + qpt*unit_qdir(:,iqdir) 866 intsum(:,:,ik, idir) = & 867 abs(eigenval(ieig,ik,idir))*unitary_33 - & 868 ( signpm*hamiltonian(lutt_params, k_plus_q ) + omega_zero(iphon)*unitary_33 ) 869 temporary1 = invmat3( intsum(:,:,ik, idir) ) 870 !Use the eigenvector from the second k point along the direction. 871 temporary2 = matmul( temporary1(1:3,1:3), eigenvec(1:3,ieig,2,idir) ) 872 temporary3 = dot_product( eigenvec(1:3,ieig,2,idir), temporary2(1:3) ) 873 sigma(ik,ieig,idir) = sigma(ik,ieig,idir) & 874 + signpm*piinv*invepsilonstar(iphon)*omega_zero(iphon) & 875 *temporary3*(omega_zero(iphon)/abs(lutt_params(1)))**half/(cos(xi))**two * weight_qdir(iqdir) 876 enddo ! iqdir 877 enddo ! ixi 878 enddo ! ik 879 enddo ! ieig 880 enddo ! idir 881 endif 882 enddo ! iphon 883 884 !Normalize self-energy integral 885 sigma = sigma/nq_factor 886 887 ABI_MALLOC(d2sigmadk2,(deg_dim,nkdir)) 888 889 d2sigmadk2 = zero 890 891 do idir = 1,nkdir 892 do ieig=1,3 893 !Compute effective mass for the correct eigenvector 894 invemass_ieig(ieig,idir)=2.0_dp*eigenval(ieig,2,idir)/krange**2.0_dp 895 !Summation over k-range 896 !Here consider a small BZ region around Gamma 897 !A convergene study might be performed around this value 898 !Finite difference for the 2nd derivative of the self energy 899 d2sigmadk2(ieig,idir) = 2.0_dp/krange**2.0_dp* (four/three)* & 900 & ( sigma(2,ieig,idir) - sigma(1,ieig,idir) - (sigma(3,ieig,idir) - sigma(1,ieig,idir))/16.0_dp ) 901 enddo 902 enddo 903 904 ABI_MALLOC(invpolmass,(deg_dim,nkdir)) 905 906 do idir = 1,nkdir 907 do ieig=1,3 908 !Calculate the inverse polaron mass 909 invpolmass(ieig,idir) = invemass_ieig(ieig,idir) + d2sigmadk2( ieig, idir) 910 !DEBUG 911 ! write(std_out,*)' idir, ieig, invemass(ieig,idir),invemass_ieig(ieig,idir), d2sigmadk2( ieig, idir), invpolmass(ieig,idir)=',& 912 !& idir, ieig, invemass(ieig,idir), invemass_ieig(ieig,idir), d2sigmadk2( ieig, idir), invpolmass(ieig,idir) 913 !ENDDEBUG 914 enddo 915 enddo 916 917 !Print inverse electronic effective masses in the output 918 write(ab_out,'(a)')'--------------------------------------------------------------------------------' 919 write(ab_out,'(a)')' Polaron properties from the generalized Froehlich model' 920 write(ab_out,'(a)')'--------------------------------------------------------------------------------' 921 write(ab_out,'(a)')' Polar modes' 922 write(ab_out,'(a)')' ## Frequency(meV) Epsilon*' 923 do imode = 1,3*cryst%natom 924 if(invepsilonstar(imode) > tol10) then 925 write(ab_out,'(2x,i3,5x,f15.6,5x,f15.6)')imode,omega_zero(imode)*Ha_eV*1000.0_dp,1.0_dp/invepsilonstar(imode) 926 endif 927 enddo 928 write(ab_out,'(a)')' ' 929 write(ab_out,'(a,f10.2)')' ZPR (meV): ',sigma(1,1,1)*Ha_eV*1000.0_dp 930 write(ab_out,'(a)')' ' 931 write(ab_out,'(a)')' Electronic effective mass (a.u.) along 3 directions' 932 write(ab_out,'(a, 3f15.6)')' Direction 100: ',one/invemass_ieig(:,1) 933 write(ab_out,'(a, 3f15.6)')' Direction 110: ',one/invemass_ieig(:,2) 934 write(ab_out,'(a, 3f15.6)')' Direction 111: ',one/invemass_ieig(:,3) 935 936 !Print inverse polaron effective masses in the output 937 write(ab_out,'(a)')' ' 938 write(ab_out,'(a)')' Polaron effective mass (a.u.) along 3 directions' 939 write(ab_out,'(a, 3f15.6)')' Direction 100: ',one/invpolmass(:,1) 940 write(ab_out,'(a, 3f15.6)')' Direction 110: ',one/invpolmass(:,2) 941 write(ab_out,'(a, 3f15.6)')' Direction 111: ',one/invpolmass(:,3) 942 write(ab_out,'(a)')' ' 943 write(ab_out,'(a)')' Sum rule of inverse polaron masses check-up (for convergence purposes):' 944 write(ab_out,'(a, 3f15.6)')' Direction 100: ',SUM(invpolmass(:,1)) 945 write(ab_out,'(a, 3f15.6)')' Direction 110: ',SUM(invpolmass(:,2)) 946 write(ab_out,'(a, 3f15.6)')' Direction 111: ',SUM(invpolmass(:,3)) 947 948 949 ABI_FREE(lutt_unit_kdir) 950 ABI_FREE(eigenvec) 951 ABI_FREE(eigenval) 952 953 ABI_FREE(weight_qdir) 954 ABI_FREE(unit_qdir) 955 ABI_FREE(polarity_qdir) 956 ABI_FREE(phfrq_qdir) 957 ABI_FREE(dielt_qdir) 958 ABI_FREE(invepsilonstar) 959 ABI_FREE(omega_zero) 960 961 ABI_FREE(intsum) 962 ABI_FREE(sigma) 963 ABI_FREE(d2sigmadk2) 964 ABI_FREE(invpolmass) 965 ABI_FREE(invemass) 966 ABI_FREE(invemass_ieig) 967 968 end subroutine polaronmass