TABLE OF CONTENTS
ABINIT/get_nv_fs_en [ Functions ]
NAME
get_nv_fs_en
FUNCTION
This routine finds the energy grids for the integration on epsilon and epsilon prime. It then calculates the DOS and FS averaged velocity_sq at these energies. Metals and semiconductors are treated differently, to deal correctly with the gap.
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (BXu) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
crystal<crystal_t>=data type gathering info on the crystalline structure. Ifc<ifc_type>=Object containing the interatomic force constants. elph_ds elph_ds%nband = number of bands in ABINIT elph_ds%k_fine%nkptirr = Number of irreducible points for which there exist at least one band that crosses the Fermi level. elph_ds%nbranch = number of phonon branches = 3*natom elph_ds%k_phon%nkpt = number of k points elph_ds%k_fine%irredtoGS = mapping of elph k-points to ground state grid elph_ds%minFSband = lowest band included in the FS integration elph_ds%nFSband = number of bands included in the FS integration elph_ds%fermie = fermi energy elph_ds%tempermin = minimum temperature at which resistivity etc are calculated (in K) elph_ds%temperinc = interval temperature grid on which resistivity etc are calculated (in K) elph_ds%ep_b_min= first band taken into account in FS integration (if telphint==2) elph_ds%ep_b_max= last band taken into account in FS integration (if telphint==2) elph_ds%telphint = flag for integration over the FS with 0=tetrahedra 1=gaussians elph_ds%elphsmear = smearing width for gaussian integration or buffer in energy for calculations with tetrahedra (telphint=0) elph_tr_ds elph_tr_ds%el_veloc = electronic velocities from the fine k-grid eigenGS = Ground State eigenvalues kptrlatt_fine = k-point grid vectors (if divided by determinant of present matrix) max_occ = maximal occupancy for a band
OUTPUT
elph_ds%nenergy = number of energy points for integration on epsilon elph_tr_ds%en_all = energy points elph_tr_ds%de_all = differences between energy points elph_tr_ds%dos_n = DOS at selected energy points elph_tr_ds%veloc_sq = FS averaged velocity square at selected energy points elph_tr_ds%tmp_gkk_intweight = integration weights at coarse k grid elph_tr_ds%tmp_velocwtk = velocity times integration weights at coarse k grid elph_tr_ds%tmp_vvelocwtk = velocity square times integration weights at coarse k grid
PARENTS
elphon
CHILDREN
d2c_weights,ep_el_weights,ep_fs_weights,get_veloc_tr,ifc_fourq,wrtout
SOURCE
63 #if defined HAVE_CONFIG_H 64 #include "config.h" 65 #endif 66 67 #include "abi_common.h" 68 69 subroutine get_nv_fs_en(crystal,ifc,elph_ds,eigenGS,max_occ,elph_tr_ds,omega_max) 70 71 use defs_basis 72 use defs_elphon 73 use m_io_tools 74 use m_errors 75 use m_profiling_abi 76 use m_ifc 77 78 use m_crystal, only : crystal_t 79 80 !This section has been created automatically by the script Abilint (TD). 81 !Do not modify the following lines by hand. 82 #undef ABI_FUNC 83 #define ABI_FUNC 'get_nv_fs_en' 84 use interfaces_14_hidewrite 85 use interfaces_77_ddb, except_this_one => get_nv_fs_en 86 !End of the abilint section 87 88 implicit none 89 90 !Arguments ------------------------------------ 91 !Scalars 92 real(dp), intent(in) :: max_occ 93 real(dp), intent(out) :: omega_max 94 type(ifc_type),intent(in) :: ifc 95 type(crystal_t),intent(in) :: crystal 96 type(elph_type),intent(inout) :: elph_ds 97 type(elph_tr_type),intent(inout) :: elph_tr_ds 98 !Arrays 99 100 real(dp), intent(in) :: eigenGS(elph_ds%nband,elph_ds%k_fine%nkptirr,elph_ds%nsppol) 101 102 !Local variables------------------------------- 103 !scalars 104 integer :: iFSqpt,isppol,ie1,ierr 105 integer :: i_metal,low_T 106 integer :: in_nenergy, out_nenergy 107 integer :: n_edge1, n_edge2, edge 108 integer :: ie_all, ne_all 109 integer :: sz1, sz2, sz3, sz4 110 real(dp) :: e_vb_max, e_cb_min,ucvol 111 real(dp) :: e1,max_e,fine_range 112 real(dp) :: enemin,enemax 113 real(dp) :: Temp,e_tiny,de0 114 real(dp) :: eff_mass1, eff_mass2, tmp_dos 115 character(len=500) :: message 116 !arrays 117 real(dp) :: gprimd(3,3) 118 real(dp) :: kpt_2nd(3), e_cb_2nd(2), en1(2) 119 real(dp),allocatable :: dos_e1(:,:),tmp_wtk(:,:,:,:) 120 real(dp),allocatable :: phfrq(:,:) 121 real(dp),allocatable :: displ(:,:,:,:) 122 123 ! ************************************************************************* 124 125 gprimd = crystal%gprimd 126 ucvol = crystal%ucvol 127 128 Temp = elph_ds%tempermin+elph_ds%temperinc 129 elph_ds%delta_e = kb_HaK*Temp ! about 1000 cm^-1/100, no need to be omega_max 130 max_e = elph_ds%nenergy*kb_HaK*Temp 131 e_tiny = kb_HaK*0.00001_dp ! this is the min. delta_e 132 de0 = kb_HaK*Temp ! Kb*T 133 134 in_nenergy = elph_ds%nenergy 135 136 ABI_ALLOCATE(tmp_wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%nsppol,4)) 137 ABI_ALLOCATE(dos_e1,(elph_ds%nsppol,3)) 138 139 ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_phon%nkpt)) 140 ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_phon%nkpt)) 141 142 do iFSqpt=1,elph_ds%k_phon%nkpt 143 call ifc_fourq(ifc,crystal,elph_ds%k_phon%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt)) 144 end do 145 146 omega_max = maxval(phfrq)*1.1_dp 147 ABI_DEALLOCATE(phfrq) 148 ABI_DEALLOCATE(displ) 149 150 write(message,'(a,E20.12)')' The max phonon energy is ', omega_max 151 call wrtout(std_out,message,'COLL') 152 153 enemin = elph_ds%fermie - max_e*2 154 enemax = elph_ds%fermie + max_e 155 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 156 & enemin, enemax, 4, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 157 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 158 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine, tmp_wtk) 159 160 do isppol=1,elph_ds%nsppol 161 dos_e1(isppol,1) = sum(tmp_wtk(:,:,isppol,2))/elph_ds%k_fine%nkpt 162 dos_e1(isppol,2) = sum(tmp_wtk(:,:,isppol,3))/elph_ds%k_fine%nkpt 163 dos_e1(isppol,3) = sum(tmp_wtk(:,:,isppol,4))/elph_ds%k_fine%nkpt 164 165 ! ! BXU, only treat metallic case at this moment, as variational method may not 166 ! ! apply to insulators 167 ! i_metal = -1 168 i_metal = 1 169 ! if (dos_e1(isppol,1) .gt. 0.1_dp .and. dos_e1(isppol,2) .gt. 0.1_dp .and. & 170 ! & dos_e1(isppol,3) .gt. 0.1_dp) then ! metal 171 ! i_metal = 1 172 if (i_metal == 1) then 173 write(message,'(a)')' This is a metal.' 174 call wrtout(std_out,message,'COLL') 175 176 fine_range = 1.5_dp 177 e1 = elph_ds%fermie + omega_max*fine_range 178 out_nenergy = 0 179 low_T = 1 180 if (omega_max*fine_range .lt. max_e) then 181 low_T = 0 182 de0 = omega_max*fine_range/in_nenergy ! energy spacing within Ef +/- omega_max 183 do while ((e1-elph_ds%fermie) .lt. max_e) 184 e1 = e1 + elph_ds%delta_e 185 out_nenergy = out_nenergy + 1 186 end do 187 end if 188 189 if (low_T == 0) max_e = e1 - elph_ds%fermie 190 elph_ds%nenergy = in_nenergy*2 + 1 + out_nenergy*2 191 192 else ! semiconductor/insulator, need careful consideration later 193 i_metal = 0 194 ! between CB min and the next k point, use free electron to replace 195 ! The weights will be proportional to the DOS, relative to the weights 196 ! calculated with ep_fs_weights, tetrahedron method prefered 197 198 ! output VB and CB edges for semiconductor/insulator 199 e_vb_max = maxval(eigenGS(elph_ds%minFSband+elph_ds%nFSband/2-1,:,isppol)) 200 e_cb_min = minval(eigenGS(elph_ds%minFSband+elph_ds%nFSband/2,:,isppol)) 201 e_cb_2nd(1) = eigenGS(elph_ds%minFSband+elph_ds%nFSband/2,2,isppol) 202 e_cb_2nd(2) = eigenGS(elph_ds%minFSband+elph_ds%nFSband/2+1,2,isppol) 203 write(message,'(a,E20.12,2x,E20.12)')' elphon : top of VB, bottom of CB = ',& 204 & e_vb_max, e_cb_min 205 call wrtout(std_out,message,'COLL') 206 write(message,'(a,E20.12)')' elphon : energy at the neighbor kpt = ',e_cb_2nd(1) 207 call wrtout(std_out,message,'COLL') 208 209 n_edge1 = 4 ! at the very edge 210 n_edge2 = 8 ! sparse to the end of free-electron part 211 212 kpt_2nd(:) = gprimd(:,1)*elph_ds%k_fine%kptirr(1,2) + & 213 & gprimd(:,2)*elph_ds%k_fine%kptirr(2,2) + & 214 & gprimd(:,3)*elph_ds%k_fine%kptirr(3,2) 215 write(message,'(a,3E20.12)')' The neighbor k point is: ', elph_ds%k_fine%kptirr(:,2) 216 call wrtout(std_out,message,'COLL') 217 218 if (dabs(elph_ds%fermie-e_cb_min) .lt. dabs(elph_ds%fermie-e_vb_max)) then 219 e1 = e_cb_2nd(1) 220 else 221 e1 = e_vb_max 222 end if 223 call ep_fs_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 224 & e1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 225 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 226 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine) 227 228 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 229 230 eff_mass1 = (kpt_2nd(1)*kpt_2nd(1) + kpt_2nd(2)*kpt_2nd(2) + kpt_2nd(3)*kpt_2nd(3)) / & 231 & (2.0_dp*(e_cb_2nd(1)-e_cb_min)) 232 write(message,'(a,E20.12)')' The eff. mass from band1 is: ', eff_mass1 233 call wrtout(std_out,message,'COLL') 234 eff_mass2 = (kpt_2nd(1)*kpt_2nd(1) + kpt_2nd(2)*kpt_2nd(2) + kpt_2nd(3)*kpt_2nd(3)) / & 235 & (2.0_dp*(e_cb_2nd(2)-e_cb_min)) 236 write(message,'(a,E20.12)')' The eff. mass from band2 is: ', eff_mass2 237 call wrtout(std_out,message,'COLL') 238 239 ! bxu, but the eff. mass estimated in this way is too small 240 ! The following is obtained by roughly fitting to the DOS of 48x48x48 241 eff_mass1 = 0.91036 242 write(message,'(a,E20.12)')' The eff. mass we are using is: ', eff_mass1 243 call wrtout(std_out,message,'COLL') 244 245 tmp_dos = (ucvol/2.0_dp/pi**2.0_dp)*(2.0_dp*eff_mass1)**1.5_dp*(e1-e_cb_min)**0.5_dp + & 246 & 2.0_dp*(ucvol/2.0_dp/pi**2.0_dp)*(2.0_dp*eff_mass2)**1.5_dp*(e1-e_cb_min)**0.5_dp 247 write(message,'(a,E20.12)')' The fake DOS at kpt1 = ', tmp_dos 248 call wrtout(std_out,message,'COLL') 249 write(message,'(a,E20.12)')' The calculated DOS at kpt1 = ', elph_ds%n0(isppol) 250 call wrtout(std_out,message,'COLL') 251 252 253 e1 = elph_ds%fermie - max_e 254 ie_all = 1 255 ne_all = 0 256 edge = 0 257 258 call ep_fs_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 259 & e1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 260 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 261 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine) 262 263 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 264 do while ((e1-elph_ds%fermie) .lt. max_e) 265 if (e1 .lt. e_cb_min .and. elph_ds%n0(isppol) .lt. tol9) then 266 e1 = e_cb_2nd(1) 267 edge = 1 268 e1 = e1 + de0 269 end if 270 271 if (e1 .lt. e_cb_2nd(1)) then 272 e1 = e_cb_2nd(1) 273 edge = 1 274 e1 = e1 + de0 275 end if 276 277 if (e1 .gt. e_cb_2nd(1)) then 278 call ep_fs_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 279 & e1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 280 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 281 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine) 282 283 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 284 285 e1 = e1 + de0 286 ie_all = ie_all + 1 287 end if 288 end do ! e_all 289 ne_all = ie_all - 1 + (n_edge1 + n_edge2 - 1)*edge ! energy levels in the free-electron range 290 write(message,'(a,i3,a,i3,a)')' For spin', isppol, ' there are ', & 291 & ne_all, ' energy levels considered ' 292 call wrtout(std_out,message,'COLL') 293 294 elph_ds%nenergy = ne_all 295 end if ! metal or insulator 296 end do ! isppol 297 298 ABI_DEALLOCATE(tmp_wtk) 299 300 if (elph_ds%nenergy .lt. 2) then 301 MSG_ERROR('There are too few energy levels for non-LOVA') 302 end if 303 304 sz1=elph_ds%ngkkband;sz2=elph_ds%k_phon%nkpt 305 sz3=elph_ds%nsppol;sz4=elph_ds%nenergy+1 306 ABI_ALLOCATE(elph_tr_ds%dos_n,(sz4,sz3)) 307 ABI_ALLOCATE(elph_tr_ds%veloc_sq,(3,sz3,sz4)) 308 ABI_ALLOCATE(elph_tr_ds%en_all,(sz3,sz4)) 309 ABI_ALLOCATE(elph_tr_ds%de_all,(sz3,sz4+1)) 310 ABI_ALLOCATE(elph_tr_ds%tmp_gkk_intweight,(sz1,sz2,sz3,sz4)) 311 ABI_ALLOCATE(elph_tr_ds%tmp_velocwtk,(sz1,sz2,3,sz3,sz4)) 312 ABI_ALLOCATE(elph_tr_ds%tmp_vvelocwtk,(sz1,sz2,3,3,sz3,sz4)) 313 314 elph_tr_ds%dos_n = zero 315 elph_tr_ds%veloc_sq = zero 316 elph_tr_ds%tmp_gkk_intweight = zero 317 elph_tr_ds%tmp_velocwtk = zero 318 elph_tr_ds%tmp_vvelocwtk = zero 319 320 ABI_STAT_ALLOCATE(elph_ds%k_phon%velocwtk,(elph_ds%nFSband,elph_ds%k_phon%nkpt,3,elph_ds%nsppol), ierr) 321 ABI_CHECK(ierr==0, 'allocating elph_ds%k_phon%velocwtk') 322 323 ABI_STAT_ALLOCATE(elph_ds%k_phon%vvelocwtk,(elph_ds%nFSband,elph_ds%k_phon%nkpt,3,3,elph_ds%nsppol), ierr) 324 ABI_CHECK(ierr==0, 'allocating elph_ds%k_phon%vvelocwtk') 325 326 elph_ds%k_phon%velocwtk = zero 327 elph_ds%k_phon%vvelocwtk = zero 328 329 !metal 330 if (i_metal .eq. 1) then 331 e1 = elph_ds%fermie - max_e 332 en1(:) = elph_ds%fermie - max_e 333 if (low_T .eq. 1) then 334 enemin = elph_ds%fermie - max_e - elph_ds%delta_e 335 enemax = elph_ds%fermie + max_e 336 337 ABI_ALLOCATE(tmp_wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%nsppol,elph_ds%nenergy+1)) 338 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 339 & enemin, enemax, elph_ds%nenergy+1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 340 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 341 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine, tmp_wtk) 342 343 do isppol=1,elph_ds%nsppol 344 do ie1 = 1, elph_ds%nenergy 345 elph_tr_ds%en_all(isppol,ie1) = en1(isppol) 346 elph_tr_ds%de_all(isppol,ie1) = elph_ds%delta_e 347 348 elph_ds%k_fine%wtk(:,:,isppol) = tmp_wtk(:,:,isppol,ie1+1) 349 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt ! for get_veloc_tr 350 elph_tr_ds%dos_n(ie1,isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 351 352 call get_veloc_tr(elph_ds,elph_tr_ds) 353 elph_tr_ds%veloc_sq(:,isppol,ie1)=elph_tr_ds%FSelecveloc_sq(:,isppol) 354 355 call d2c_weights(elph_ds,elph_tr_ds) 356 357 elph_tr_ds%tmp_gkk_intweight(:,:,isppol,ie1) = elph_ds%k_phon%wtk(:,:,isppol) 358 elph_tr_ds%tmp_velocwtk(:,:,:,isppol,ie1) = elph_ds%k_phon%velocwtk(:,:,:,isppol) 359 elph_tr_ds%tmp_vvelocwtk(:,:,:,:,isppol,ie1) = elph_ds%k_phon%vvelocwtk(:,:,:,:,isppol) 360 en1(isppol) = en1(isppol) + elph_ds%delta_e 361 end do 362 end do 363 ABI_DEALLOCATE(tmp_wtk) 364 365 else ! low_T = 0 366 enemin = e1 - elph_ds%delta_e 367 enemax = e1 + (out_nenergy-1)*elph_ds%delta_e 368 369 ABI_ALLOCATE(tmp_wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%nsppol,out_nenergy+1)) 370 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 371 & enemin, enemax, out_nenergy+1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 372 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 373 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine, tmp_wtk) 374 do isppol=1,elph_ds%nsppol 375 do ie1 = 1, out_nenergy 376 elph_tr_ds%en_all(isppol,ie1) = en1(isppol) 377 elph_tr_ds%de_all(isppol,ie1) = elph_ds%delta_e 378 379 elph_ds%k_fine%wtk(:,:,isppol) = tmp_wtk(:,:,isppol,ie1+1) 380 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt ! for get_veloc_tr 381 elph_tr_ds%dos_n(ie1,isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 382 383 call get_veloc_tr(elph_ds,elph_tr_ds) 384 elph_tr_ds%veloc_sq(:,isppol,ie1)=elph_tr_ds%FSelecveloc_sq(:,isppol) 385 386 call d2c_weights(elph_ds,elph_tr_ds) 387 388 elph_tr_ds%tmp_gkk_intweight(:,:,isppol,ie1) = elph_ds%k_phon%wtk(:,:,isppol) 389 elph_tr_ds%tmp_velocwtk(:,:,:,isppol,ie1) = elph_ds%k_phon%velocwtk(:,:,:,isppol) 390 elph_tr_ds%tmp_vvelocwtk(:,:,:,:,isppol,ie1) = elph_ds%k_phon%vvelocwtk(:,:,:,:,isppol) 391 392 en1(isppol) = en1(isppol) + elph_ds%delta_e 393 end do 394 end do 395 ABI_DEALLOCATE(tmp_wtk) 396 397 e1 = en1(1) 398 enemin = e1 - de0 399 enemax = e1 + in_nenergy*2*de0 400 401 ABI_ALLOCATE(tmp_wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%nsppol,in_nenergy*2+2)) 402 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 403 & enemin, enemax, in_nenergy*2+2, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 404 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 405 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine, tmp_wtk) 406 407 do isppol=1,elph_ds%nsppol 408 do ie1 = out_nenergy+1, out_nenergy+in_nenergy*2+1 409 elph_tr_ds%en_all(isppol,ie1) = en1(isppol) 410 elph_tr_ds%de_all(isppol,ie1) = de0 411 412 elph_ds%k_fine%wtk(:,:,isppol) = tmp_wtk(:,:,isppol,ie1-out_nenergy+1) 413 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt ! for get_veloc_tr 414 elph_tr_ds%dos_n(ie1,isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 415 416 call get_veloc_tr(elph_ds,elph_tr_ds) 417 elph_tr_ds%veloc_sq(:,isppol,ie1)=elph_tr_ds%FSelecveloc_sq(:,isppol) 418 419 call d2c_weights(elph_ds,elph_tr_ds) 420 421 elph_tr_ds%tmp_gkk_intweight(:,:,isppol,ie1) = elph_ds%k_phon%wtk(:,:,isppol) 422 elph_tr_ds%tmp_velocwtk(:,:,:,isppol,ie1) = elph_ds%k_phon%velocwtk(:,:,:,isppol) 423 elph_tr_ds%tmp_vvelocwtk(:,:,:,:,isppol,ie1) = elph_ds%k_phon%vvelocwtk(:,:,:,:,isppol) 424 425 en1(isppol) = en1(isppol) + de0 426 end do 427 end do 428 ABI_DEALLOCATE(tmp_wtk) 429 430 e1 = en1(1) 431 enemin = e1 - elph_ds%delta_e 432 enemax = e1 + (out_nenergy-1)*elph_ds%delta_e 433 434 ABI_ALLOCATE(tmp_wtk,(elph_ds%nFSband,elph_ds%k_fine%nkpt,elph_ds%nsppol,out_nenergy+1)) 435 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 436 & enemin, enemax, out_nenergy+1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 437 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 438 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine, tmp_wtk) 439 440 en1(:) = en1(:) - de0 + elph_ds%delta_e ! adjust to make the points symmetric around Ef 441 do isppol=1,elph_ds%nsppol 442 do ie1 = out_nenergy+in_nenergy*2+2, in_nenergy*2+1+out_nenergy*2 443 elph_tr_ds%en_all(isppol,ie1) = en1(isppol) 444 elph_tr_ds%de_all(isppol,ie1) = elph_ds%delta_e 445 446 elph_ds%k_fine%wtk(:,:,isppol) = tmp_wtk(:,:,isppol,ie1-out_nenergy-in_nenergy*2) 447 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt ! for get_veloc_tr 448 elph_tr_ds%dos_n(ie1,isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 449 450 call get_veloc_tr(elph_ds,elph_tr_ds) 451 elph_tr_ds%veloc_sq(:,isppol,ie1)=elph_tr_ds%FSelecveloc_sq(:,isppol) 452 453 call d2c_weights(elph_ds,elph_tr_ds) 454 455 elph_tr_ds%tmp_gkk_intweight(:,:,isppol,ie1) = elph_ds%k_phon%wtk(:,:,isppol) 456 elph_tr_ds%tmp_velocwtk(:,:,:,isppol,ie1) = elph_ds%k_phon%velocwtk(:,:,:,isppol) 457 elph_tr_ds%tmp_vvelocwtk(:,:,:,:,isppol,ie1) = elph_ds%k_phon%vvelocwtk(:,:,:,:,isppol) 458 459 en1(isppol) = en1(isppol) + elph_ds%delta_e 460 end do 461 end do 462 ABI_DEALLOCATE(tmp_wtk) 463 end if 464 465 !semiconductor 466 else if (i_metal .eq. 0) then 467 e1 = elph_ds%fermie - max_e 468 ie_all = 1 469 470 call ep_fs_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 471 & e1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 472 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 473 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine) 474 475 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 476 do while ((e1-elph_ds%fermie) .lt. max_e) 477 if (e1 .lt. e_cb_min .and. elph_ds%n0(isppol) .lt. tol9) then 478 e1 = e_cb_min 479 end if 480 481 if (ie_all .ge. n_edge1+n_edge2) then 482 if (ie_all .eq. n_edge1+n_edge2) e1 = e1 + de0 483 call ep_fs_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 484 & e1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 485 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 486 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine) 487 488 elph_tr_ds%tmp_gkk_intweight(:,:,isppol,ie_all) = elph_ds%k_fine%wtk(:,:,isppol) 489 elph_tr_ds%dos_n(ie_all,isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt 490 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt ! for get_veloc_tr 491 492 elph_tr_ds%en_all(isppol,ie_all) = e1 493 call get_veloc_tr(elph_ds,elph_tr_ds) 494 elph_tr_ds%veloc_sq(:,isppol,ie_all)=elph_tr_ds%FSelecveloc_sq(:,isppol) 495 ! bxu 496 ! veloc_sq(1,isppol,ie_all) is "1" good and general?? 497 498 elph_tr_ds%de_all(isppol,ie_all) = de0 499 e1 = e1 + elph_tr_ds%de_all(isppol,ie_all) 500 ie_all = ie_all + 1 501 else ! divided according to the 1/DOS (evenly) 502 if (ie_all .lt. n_edge1) then 503 elph_tr_ds%en_all(isppol,ie_all) = e_cb_min + & 504 & (e_tiny**(-0.5_dp) - ie_all*(e_tiny**(-0.5_dp)-(e_cb_2nd(1)-e_cb_min)**(-0.5_dp))/ & 505 & dble(n_edge1))**(-2.0_dp) 506 if (ie_all .gt. 1) then 507 elph_tr_ds%de_all(isppol,ie_all) = elph_tr_ds%en_all(isppol,ie_all) - elph_tr_ds%en_all(isppol,ie_all-1) 508 else 509 elph_tr_ds%de_all(isppol,ie_all) = elph_tr_ds%en_all(isppol,ie_all) - e_cb_min - e_tiny 510 end if 511 e1 = elph_tr_ds%en_all(isppol,ie_all) 512 else 513 elph_tr_ds%en_all(isppol,ie_all) = e_cb_min + & 514 & ((ie_all-n_edge1+1)/dble(n_edge2))**2.0_dp*(e_cb_2nd(1)-e_cb_min) 515 if (ie_all .gt. 1) then 516 elph_tr_ds%de_all(isppol,ie_all) = elph_tr_ds%en_all(isppol,ie_all) - elph_tr_ds%en_all(isppol,ie_all-1) 517 else 518 elph_tr_ds%de_all(isppol,ie_all) = (e_cb_2nd(1)-e_cb_min)/(dble(n_edge2)**2.0_dp) 519 end if 520 e1 = elph_tr_ds%en_all(isppol,ie_all) 521 end if 522 523 call ep_fs_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS, elph_ds%elphsmear, & 524 & e1, gprimd, elph_ds%k_fine%irredtoGS, elph_ds%kptrlatt_fine, max_occ, & 525 & elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, & 526 & elph_ds%nsppol, elph_ds%telphint, elph_ds%k_fine) 527 528 elph_ds%n0(isppol) = sum(elph_ds%k_fine%wtk(:,:,isppol))/elph_ds%k_fine%nkpt ! for get_veloc_tr 529 530 tmp_dos = (ucvol/2.0_dp/pi**2.0_dp)*(2.0_dp*eff_mass1)**1.5_dp*(e1-e_cb_min)**0.5_dp + & 531 & 2.0_dp*(ucvol/2.0_dp/pi**2.0_dp)*(2.0_dp*eff_mass2)**1.5_dp*(e1-e_cb_min)**0.5_dp 532 elph_tr_ds%dos_n(ie_all,isppol) = tmp_dos 533 elph_tr_ds%tmp_gkk_intweight(:,:,isppol,ie_all) = elph_ds%k_fine%wtk(:,:,isppol)*tmp_dos/elph_ds%n0(isppol) 534 535 call get_veloc_tr(elph_ds,elph_tr_ds) 536 elph_tr_ds%veloc_sq(:,isppol,ie_all)=elph_tr_ds%FSelecveloc_sq(:,isppol) 537 538 if (ie_all .eq. (n_edge1+n_edge2)) e1 = e_cb_2nd(1) + de0 539 ie_all = ie_all + 1 540 end if 541 end do ! ie_all 542 else 543 MSG_BUG('check i_metal!') 544 end if ! metal or insulator 545 546 ABI_DEALLOCATE(dos_e1) 547 548 end subroutine get_nv_fs_en