TABLE OF CONTENTS


ABINIT/get_nv_fs_en [ Functions ]

[ Top ] [ 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