TABLE OF CONTENTS


ABINIT/get_tau_k [ Functions ]

[ Top ] [ Functions ]

NAME

  get_tau_k

FUNCTION

  Calculate the k-dependent relaxation time due to EPC. Impelementation based
  on derivation from Grmvall's book or OD Restrepo's paper (PRB 94 212103 (2009))

COPYRIGHT

  Copyright (C) 2013-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

  Cryst<crystal_t>=Info on the unit cell and on its symmetries.
  Ifc<ifc_type>=Object containing the interatomic force constants.
  elph_ds = elphon datastructure with data and dimensions
  eigenGS = Ground State eigenvalues
  max_occ = maximal occupancy for a band

OUTPUT

  tau_k(nsppol,nkptirr,nband)=mode relaxation time due to electron phonono coupling
  rate_e(nene)= scattering rate due to electron phonono coupling vs. energy

PARENTS

      elphon

CHILDREN

      dgemm,ebands_prtbltztrp_tau_out,ebands_update_occ,ep_el_weights
      ep_ph_weights,ftgam,ftgam_init,gam_mult_displ,ifc_fourq,matrginv
      mkqptequiv,phdispl_cart2red,spline,splint,wrtout,xmpi_sum,zgemm

SOURCE

 37 #if defined HAVE_CONFIG_H
 38 #include "config.h"
 39 #endif
 40 
 41 #include "abi_common.h"
 42 
 43 subroutine get_tau_k(Cryst,ifc,Bst,elph_ds,elph_tr_ds,eigenGS,max_occ)
 44 
 45  use defs_basis
 46  use defs_elphon
 47  use defs_datatypes
 48  use m_kptrank
 49  use m_errors
 50  use m_profiling_abi
 51  use m_xmpi
 52  use m_splines
 53  use m_ifc
 54  use m_ebands
 55 
 56  use m_io_tools,   only : open_file
 57  use m_abilasi,    only : matrginv
 58  use m_geometry,   only : phdispl_cart2red
 59  use m_dynmat,     only : ftgam_init, ftgam
 60  use m_crystal,    only : crystal_t
 61 
 62 !This section has been created automatically by the script Abilint (TD).
 63 !Do not modify the following lines by hand.
 64 #undef ABI_FUNC
 65 #define ABI_FUNC 'get_tau_k'
 66  use interfaces_14_hidewrite
 67  use interfaces_77_ddb, except_this_one => get_tau_k
 68 !End of the abilint section
 69 
 70  implicit none
 71 
 72 !Arguments ------------------------------------
 73  type(crystal_t),intent(in) :: Cryst
 74  type(ifc_type),intent(in) :: ifc
 75  type(ebands_t),intent(inout)   :: Bst
 76  type(elph_type),intent(inout) :: elph_ds
 77  type(elph_tr_type), intent(inout) :: elph_tr_ds
 78  real(dp),intent(in) :: max_occ
 79  real(dp),intent(in) :: eigenGS(elph_ds%nband,elph_ds%k_phon%nkpt,elph_ds%nsppol)
 80 
 81 !Local variables-------------------------------
 82 !scalars
 83  character(len=500) :: message
 84  character(len=fnlen) :: fname
 85  integer :: ntemper,nsppol,nbranch,nband,natom
 86  integer :: nkpt,nqpt,nkptirr,nqptirr,new_nkptirr
 87  integer :: isppol,iFSkpt,iFSqpt,iqpt,iqpt_fullbz,imqpt_fullbz,ikpt_kpq,ikpt_kmq
 88  integer :: iband,jband,jpband,jbeff,ibranch,jbranch,itemp
 89  integer :: irec,ierr,nrpt,ik_this_proc
 90  integer :: unit_tau,unit_invtau
 91  integer :: nene,nene_all,iene,iene_fine,unit_taue,unit_mfp
 92  integer :: icomp,jcomp,itensor
 93  integer :: ikpt_irr,iomega,unit_cond,unit_therm,unit_sbk
 94  integer :: nskip,nspline
 95  real(dp) :: occ_omega,occ_e
 96  real(dp) :: xx,Temp,therm_factor
 97  real(dp) :: factor,dfermide
 98  real(dp) :: e_k,chu_tau,rate_e,mfp_e
 99  real(dp) :: ene,enemin,enemax,deltaene
100  real(dp) :: omega,omega_min,omega_max,domega
101  real(dp) :: diagerr
102  real(dp) :: chu_mfp,chu_cond,chu_cth,chu_sbk,femto
103  real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
104  real(dp) :: eigval(elph_ds%nbranch),eigval2(elph_ds%nbranch)
105  real(dp) :: imeigval(elph_ds%nbranch)
106  real(dp) :: tmp_wtkpq, tmp_wtkmq, tol_wtk
107  real(dp) :: yp1,ypn
108 !arrays
109  integer,allocatable :: FSfullpktofull(:,:),mqtofull(:)
110  integer,allocatable :: kpttokpt(:,:,:)
111  real(dp) :: cond_inv(3,3)
112  real(dp),allocatable :: fermie(:)
113  real(dp),allocatable :: tmp_eigenGS(:,:,:)
114  real(dp),allocatable :: tmp_gkk_qpt(:,:,:),tmp_gkk_rpt(:,:,:),tmp_gkk_kpt(:,:)
115  real(dp),allocatable :: tmp_gkk_kpt2(:,:,:), gkk_kpt(:,:,:)
116  real(dp),allocatable :: tau_k(:,:,:,:),inv_tau_k(:,:,:,:),tmp_tau_k(:,:,:,:)
117  real(dp),allocatable :: phfrq(:,:),pheigvec(:,:)
118  real(dp),allocatable :: displ(:,:,:,:)
119  real(dp),allocatable :: a2f_2d(:),a2f_2d2(:)
120  real(dp),allocatable :: tmp_wtk(:,:,:,:),tmp2_wtk(:),tmp_wtk1(:),tmp_wtk2(:)
121  real(dp),allocatable :: ene_pt(:),ene_ptfine(:),ff2(:)
122  real(dp),allocatable :: wtq(:,:,:),tmp_wtq(:,:,:),tmp2_wtq(:,:)
123  real(dp),allocatable :: dos_e(:,:)
124  real(dp),allocatable :: coskr1(:,:),sinkr1(:,:)
125  real(dp),allocatable :: coskr2(:,:),sinkr2(:,:)
126  real(dp),allocatable :: cond_e(:,:,:,:),cond(:,:,:,:),sbk(:,:,:,:),seebeck(:,:,:,:),cth(:,:,:,:)
127 
128 ! *************************************************************************
129 
130  write(std_out,*) 'get_tau_k : enter '
131 
132  nrpt = ifc%nrpt
133  natom = cryst%natom
134 
135  nsppol   = elph_ds%nsppol
136  nbranch  = elph_ds%nbranch
137  nband    = elph_ds%ngkkband
138  nkpt     = elph_ds%k_phon%nkpt
139  nqpt     = elph_ds%nqpt_full
140  nkptirr  = elph_ds%k_phon%nkptirr
141  new_nkptirr  = elph_ds%k_phon%new_nkptirr
142  nqptirr  = elph_ds%nqptirred
143  ntemper  = elph_ds%ntemper
144  nene = 2*elph_ds%na2f-1 ! only need e_k +- omega_max range, take deltaene=delta_oemga
145 
146  chu_tau  = 2.4188843265*1.0d-17
147  chu_mfp  = 5.291772*1.0d-11
148  chu_cond = 4.59988159904764*1.0d6
149  chu_cth  = 1.078637439971599*1.0d4
150  chu_sbk  = 8.617343101*1.0d-5
151  femto    = 1.0d-15
152 
153  tol_wtk = tol7/nkptirr/nband
154 
155  ABI_ALLOCATE(fermie ,(ntemper))
156  ABI_ALLOCATE(tmp_gkk_qpt ,(2,nbranch**2,nqpt))
157  ABI_ALLOCATE(tmp_gkk_rpt ,(2,nbranch**2,nrpt))
158  ABI_ALLOCATE(tmp_gkk_kpt ,(2,nbranch**2))
159  ABI_ALLOCATE(tmp_gkk_kpt2 ,(2,nbranch,nbranch))
160  ABI_ALLOCATE(gkk_kpt ,(2,nbranch,nbranch))
161  ABI_ALLOCATE(a2f_2d, (nene))
162  ABI_ALLOCATE(a2f_2d2, (nene))
163  ABI_ALLOCATE(inv_tau_k, (ntemper,nsppol,nkpt,nband))
164  ABI_ALLOCATE(tau_k, (ntemper,nsppol,nkpt,nband))
165  ABI_ALLOCATE(tmp_tau_k ,(ntemper,nsppol,new_nkptirr,nband))
166 
167  if (elph_ds%gkqwrite == 0) then
168    call wrtout(std_out,' get_tau_k : keeping gkq matrices in memory','COLL')
169  else if (elph_ds%gkqwrite == 1) then
170    fname=trim(elph_ds%elph_base_name) // '_GKKQ'
171    write (message,'(2a)')' get_tau_k : reading gkq matrices from file ',trim(fname)
172    call wrtout(std_out,message,'COLL')
173  else
174    write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite
175    MSG_BUG(message)
176  end if
177 
178 !=========================================================
179 !Get equivalence between a kpt_phon pair and a qpt in qpt_full
180 !only works if the qpt grid is complete (identical to
181 !the kpt one, with a basic shift of (0,0,0)
182 !=========================================================
183 
184 !mapping of k + q onto k' for k and k' in full BZ
185 !for dense k grid
186  ABI_ALLOCATE(FSfullpktofull,(nkpt,nkpt))
187  ABI_ALLOCATE(mqtofull,(nkpt))
188 
189 !kpttokpt(itim,isym,iqpt) = kpoint index which transforms to ikpt under isym and with time reversal itim.
190  ABI_ALLOCATE(kpttokpt,(2,Cryst%nsym,nkpt))
191 
192  call wrtout(std_out,'get_tau_k: calling mkqptequiv to set up the FS kpoint set',"COLL")
193 
194  call mkqptequiv (FSfullpktofull,Cryst,elph_ds%k_phon%kpt,nkpt,nkpt,kpttokpt,elph_ds%k_phon%kpt,mqtofull)
195 
196 !=========================================================
197 !=========================================================
198 
199  omega_max       = elph_ds%omega_max
200  omega_min       = elph_ds%omega_min
201  domega          = elph_ds%domega
202  enemax = maxval(eigenGS(elph_ds%maxFSband,:,:))
203  enemin = minval(eigenGS(elph_ds%minFSband,:,:))
204 
205  if (enemin < (elph_ds%fermie-0.2)) then
206    enemin = elph_ds%fermie-0.2
207  end if
208  if (enemax > (elph_ds%fermie+0.2)) then
209    enemax = elph_ds%fermie+0.2
210  end if
211 
212  nspline = elph_ds%ep_nspline
213  nene_all = INT((enemax-enemin+domega)/(nspline*domega)) + 1
214  deltaene = domega
215  write(std_out,*) 'E_min= ',enemin, 'E_max= ',enemax
216  write(std_out,*) 'Number of energy points= ',nene_all
217  write(std_out,'(a,I8)') 'scale factor for spline interpolation in RTA = ', elph_ds%ep_nspline
218  write(std_out,*) 'delta_ene before spline interpolation= ',deltaene*nspline
219  write(std_out,*) 'delta_ene after spline interpolation= ',deltaene
220  write(std_out,*) 'Omega_min= ',omega_min, 'Omega_max= ',omega_max
221  write(std_out,*) 'Number of phonon points= ',elph_ds%na2f
222  write(std_out,*) 'delta_omega= ',domega
223  write(std_out,*) 'number of bands= ', elph_ds%nband, nband
224 
225  ABI_ALLOCATE(tmp_wtk,(nband,nkpt,nsppol,nene_all))
226  ABI_ALLOCATE(tmp2_wtk,(nene_all))
227  ABI_ALLOCATE(ff2,(nene_all))
228  ABI_ALLOCATE(ene_pt,(nene_all))
229  ABI_ALLOCATE(ene_ptfine,(nene_all*nspline))
230  ABI_ALLOCATE(tmp_wtk1,(nene_all*nspline))
231  ABI_ALLOCATE(tmp_wtk2,(nene_all*nspline))
232  ABI_ALLOCATE(dos_e,(nsppol,nene_all))
233 
234 !Get energy points for spline interpolation
235  do iene = 1, nene_all
236    ene_pt(iene) = enemin + (iene-1)*nspline*deltaene
237  end do
238 
239  do iene = 1, nene_all*nspline
240    ene_ptfine(iene) = enemin + (iene-1)*deltaene
241  end do
242 
243  ABI_ALLOCATE(tmp_wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f+1))
244  ABI_ALLOCATE(wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f))
245  ABI_ALLOCATE(tmp2_wtq,(elph_ds%nbranch, elph_ds%na2f))
246 
247 !phonon
248  ABI_ALLOCATE(phfrq,(nbranch, nkptirr))
249  ABI_ALLOCATE(displ,(2, nbranch, nbranch, nkptirr))
250  ABI_ALLOCATE(pheigvec,(2*nbranch*nbranch, nkptirr))
251 
252  do iFSqpt = 1, nkptirr
253    call ifc_fourq(ifc,cryst,elph_ds%k_phon%kptirr(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
254  end do
255 
256  omega_min = omega_min - domega
257 
258 !bxu, obtain wtq for the q_fine, then condense to q_phon
259  call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,Cryst%gprimd,elph_ds%kptrlatt, &
260 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_phon,tmp_wtq)
261  omega_min = omega_min + domega
262 
263  do iomega = 1, elph_ds%na2f
264    wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1)
265    !write(1005,*) omega_min+(iomega-1)*domega, sum(tmp_wtq(:,:,iomega+1))/nkpt
266  end do
267  ABI_DEALLOCATE(tmp_wtq)
268 
269 ! electron
270  tmp_wtk =zero
271  dos_e = zero
272  call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS(elph_ds%minFSband:elph_ds%minFSband+nband-1,:,:), &
273 & elph_ds%elphsmear, &
274 & enemin, enemax, nene_all, Cryst%gprimd, elph_ds%k_phon%irredtoGS, elph_ds%kptrlatt, max_occ, &
275 & 1, nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk)
276 !& elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk)
277 
278  do isppol = 1, nsppol
279    do iene = 1, nene_all
280      dos_e(isppol,iene) = sum(tmp_wtk(:,:,isppol,iene))/nkpt
281    end do
282  end do
283 
284  ABI_ALLOCATE(coskr1, (nqpt,nrpt))
285  ABI_ALLOCATE(sinkr1, (nqpt,nrpt))
286  call ftgam_init(ifc%gprim, nqpt, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr1, sinkr1)
287  ABI_ALLOCATE(coskr2, (nkptirr,nrpt))
288  ABI_ALLOCATE(sinkr2, (nkptirr,nrpt))
289  call ftgam_init(ifc%gprim, nkptirr, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr2, sinkr2)
290 
291 !get fermie for itemp
292  fermie = elph_ds%fermie
293  do itemp=1,ntemper  ! runs over termperature in K
294    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
295 
296    Bst%occopt = 3
297    Bst%tsmear = Temp*kb_HaK
298    call ebands_update_occ(Bst,-99.99_dp)
299    write(message,'(a,f12.6,a,E20.12)')'At T=',Temp,' Fermi level is:',Bst%fermie
300    call wrtout(std_out,message,'COLL')
301 
302    if (abs(elph_ds%fermie) < tol10) then
303      fermie(itemp) = Bst%fermie
304    end if
305  end do
306 
307  inv_tau_k = zero
308 !get a2f_2d = \sum_{q,nbranch,jband'} |gkk|^2*\delta(\epsilon_{k'j'}-\epsilon')*\delta(\omega_q-\omega)
309  do isppol=1,nsppol
310    write (std_out,*) '##############################################'
311    write (std_out,*) 'get_tau_k : Treating spin polarization ', isppol
312    write (std_out,*) '##############################################'
313 
314 !   do iFSkpt =1,nkpt
315    do ik_this_proc =1,elph_ds%k_phon%my_nkpt
316      iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
317      write (std_out,*) 'get_tau_k : working on kpt # ', iFSkpt, '/', nkpt
318      do jband = 1, nband
319 !          write(*,*)'i am here 1 ', isppol,iFSkpt,jband
320        a2f_2d = zero
321        a2f_2d2 = zero
322 
323 !sum from here
324        nskip = 0
325        do jpband = 1, nband
326          jbeff = jpband+(jband-1)*nband
327 
328          if (elph_ds%gkqwrite == 0) then
329            tmp_gkk_qpt(:,:,:) = elph_ds%gkk_qpt(:,jbeff,:,ik_this_proc,isppol,:)
330          else if (elph_ds%gkqwrite == 1) then
331            irec = (ik_this_proc-1)*elph_ds%k_phon%my_nkpt + iqpt
332            if (iFSkpt == 1) then
333              write (std_out,*) ' get_tau_k  read record ', irec
334            end if
335            read (elph_ds%unitgkq,REC=irec) tmp_gkk_qpt(:,:,iqpt_fullbz)
336          end if
337 
338 !FT to real space
339          call ftgam(Ifc%wghatm,tmp_gkk_qpt,tmp_gkk_rpt,natom,nqpt,nrpt,1,coskr1,sinkr1)
340 
341 !sum over irred q over k_phon, with corresponding weights
342          do iFSqpt = 1, nkptirr
343            iqpt_fullbz = elph_ds%k_phon%irredtoGS(iFSqpt)
344            ikpt_kpq = FSfullpktofull(iFSkpt,iqpt_fullbz)
345 
346            imqpt_fullbz = mqtofull(iqpt_fullbz)
347            ikpt_kmq = FSfullpktofull(iFSkpt,imqpt_fullbz)
348 
349 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr
350            call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(iqpt_fullbz,:),sinkr2(iqpt_fullbz,:))
351 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt)
352 
353 !if ep_scalprod==0 we have to dot in the displacement vectors here
354            if (elph_ds%ep_scalprod==0) then
355 
356              call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red)
357 
358              tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/))
359              call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt)
360 
361              do jbranch=1,nbranch
362                eigval(jbranch) = gkk_kpt(1, jbranch, jbranch)
363                imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch)
364 
365                if (abs(imeigval(jbranch)) > tol10) then
366                  write (message,'(a,i0,a,es16.8)')" real values  branch = ",jbranch,' eigval = ',eigval(jbranch)
367                  MSG_WARNING(message)
368                  write (message,'(a,i0,a,es16.8)')" imaginary values  branch = ",jbranch,' imeigval = ',imeigval(jbranch)
369                  MSG_WARNING(message)
370                end if
371 
372              end do
373 
374 !            if ep_scalprod==1 we have to diagonalize the matrix we interpolated.
375            else if (elph_ds%ep_scalprod == 1) then
376 
377 !            MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here
378              call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,&
379 &             pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom)
380 
381              call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,&
382 &             tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom)
383 
384              diagerr = zero
385              do ibranch=1,nbranch
386                eigval(ibranch) = gkk_kpt(1,ibranch,ibranch)
387                do jbranch=1,ibranch-1
388                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
389                end do
390                do jbranch=ibranch+1,nbranch
391                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
392                end do
393              end do
394 
395              if (diagerr > tol12) then
396                write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr
397                MSG_WARNING(message)
398              end if
399 
400            else
401              write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod
402              MSG_BUG(message)
403            end if ! end ep_scalprod if
404 
405 !For k'=k-q
406 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr
407            call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(imqpt_fullbz,:),sinkr2(imqpt_fullbz,:))
408 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt)
409 
410 !if ep_scalprod==0 we have to dot in the displacement vectors here
411            if (elph_ds%ep_scalprod==0) then
412 
413              call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red)
414 
415              tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/))
416              call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt)
417 
418              do jbranch=1,nbranch
419                eigval2(jbranch) = gkk_kpt(1, jbranch, jbranch)
420                imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch)
421 
422                if (abs(imeigval(jbranch)) > tol10) then
423                  write (message,'(a,i0,a,es16.8)')" real values  branch = ",jbranch,' eigval = ',eigval2(jbranch)
424                  MSG_WARNING(message)
425                  write (message,'(a,i0,a,es16.8)')" imaginary values  branch = ",jbranch,' imeigval = ',imeigval(jbranch)
426                  MSG_WARNING(message)
427                end if
428 
429              end do
430 
431 !            if ep_scalprod==1 we have to diagonalize the matrix we interpolated.
432            else if (elph_ds%ep_scalprod == 1) then
433 
434 !            MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here
435              call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,&
436 &             pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom)
437 
438              call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,&
439 &             tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom)
440 
441              diagerr = zero
442              do ibranch=1,nbranch
443                eigval2(ibranch) = gkk_kpt(1,ibranch,ibranch)
444                do jbranch=1,ibranch-1
445                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
446                end do
447                do jbranch=ibranch+1,nbranch
448                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
449                end do
450              end do
451 
452              if (diagerr > tol12) then
453                write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr
454                MSG_WARNING(message)
455              end if
456 
457            else
458              write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod
459              MSG_BUG(message)
460            end if ! end ep_scalprod if
461 
462            tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kpq,isppol,:)
463            yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene
464            ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene
465            call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2)
466            call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk1)
467 
468            tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kmq,isppol,:)
469            yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene
470            ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene
471            call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2)
472            call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk2)
473 
474            tmp2_wtq(:,:) = wtq(:,iFSqpt,:)
475            do iene=1,nene
476              e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol)
477              ene = e_k - omega_max + (iene-1)*deltaene
478              if (ene<enemin .or. ene>enemax) cycle
479              iene_fine = NINT((ene-enemin+deltaene)/deltaene)
480              tmp_wtkpq = tmp_wtk1(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt)
481              tmp_wtkmq = tmp_wtk2(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt)
482 
483              if (tmp_wtkpq+tmp_wtkmq < tol_wtk ) then
484                nskip = nskip +1
485                cycle
486              end if
487 
488              do ibranch = 1, nbranch
489                if (abs(phfrq(ibranch,iFSqpt)) < tol7) cycle
490 
491                if (ene > e_k) then
492                  omega = ene - e_k
493                  if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle
494                  iomega = NINT((omega-omega_min+domega)/domega)
495 
496                  a2f_2d(iene) = a2f_2d(iene) +&
497 &                 eigval(ibranch)/phfrq(ibranch,iFSqpt)*&
498 &                 tmp_wtkpq * tmp2_wtq(ibranch,iomega)
499                end if
500 
501                if (ene < e_k) then
502                  omega = e_k - ene
503                  if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle
504                  iomega = NINT((omega-omega_min+domega)/domega)
505 
506                  a2f_2d2(iene) = a2f_2d2(iene) +&
507 &                 eigval(ibranch)/phfrq(ibranch,iFSqpt)*&
508 &                 tmp_wtkmq * tmp2_wtq(ibranch,iomega)
509                end if
510 
511              end do ! ibranch 3
512            end do ! nene  800
513          end do ! kptirr 216
514        end do ! j' band 3
515 !      print *, ' skipped ',  nskip, ' energy points out of ', nene*nband*nkptirr
516 
517 ! get inv_tau_k
518        do itemp=1,ntemper  ! runs over termperature in K
519          Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
520          do iene=1,nene
521            e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol)
522            ene = e_k - omega_max + (iene-1)*deltaene
523            if (ene<enemin .or. ene>enemax) cycle
524 
525            xx=(ene-fermie(itemp))/(kb_HaK*Temp)
526            occ_e=1.0_dp/(exp(xx)+1.0_dp)
527            if (ene > e_k .and. (ene-e_k) .le. omega_max) then
528              omega = ene - e_k
529              if (abs(omega) < tol7) cycle
530              xx = omega/(kb_HaK*Temp)
531              occ_omega=1.0_dp/(exp(xx)-1.0_dp)
532 
533              therm_factor = occ_e + occ_omega
534 
535              inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +&
536              a2f_2d(iene)*therm_factor*deltaene
537            end if
538            if (ene < e_k .and. (e_k-ene) .le. omega_max) then
539              omega = e_k - ene
540              if (abs(omega) < tol7) cycle
541              xx = omega/(kb_HaK*Temp)
542              occ_omega=1.0_dp/(exp(xx)-1.0_dp)
543 
544              therm_factor = 1 - occ_e + occ_omega
545 
546              inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +&
547              a2f_2d2(iene)*therm_factor*deltaene
548            end if
549 
550          end do ! nene
551        end do ! Temp
552 !          write(*,*)'i am here 2 ', isppol,iFSkpt,jband
553      end do ! jband
554    end do ! kpt
555  end do ! nsppol
556 
557 !write (300+mpi_enreg%me,*) inv_tau_k
558  call xmpi_sum (inv_tau_k, xmpi_world, ierr)
559 
560  ABI_DEALLOCATE(phfrq)
561  ABI_DEALLOCATE(displ)
562  ABI_DEALLOCATE(pheigvec)
563  ABI_DEALLOCATE(tmp2_wtk)
564  ABI_DEALLOCATE(ff2)
565  ABI_DEALLOCATE(ene_pt)
566  ABI_DEALLOCATE(ene_ptfine)
567  ABI_DEALLOCATE(tmp_wtk1)
568  ABI_DEALLOCATE(tmp_wtk2)
569  ABI_DEALLOCATE(tmp2_wtq)
570  ABI_DEALLOCATE(wtq)
571  ABI_DEALLOCATE(coskr1)
572  ABI_DEALLOCATE(sinkr1)
573  ABI_DEALLOCATE(coskr2)
574  ABI_DEALLOCATE(sinkr2)
575  ABI_DEALLOCATE(kpttokpt)
576  ABI_DEALLOCATE(FSfullpktofull)
577  ABI_DEALLOCATE(mqtofull)
578  ABI_DEALLOCATE(tmp_gkk_qpt)
579  ABI_DEALLOCATE(tmp_gkk_rpt)
580  ABI_DEALLOCATE(tmp_gkk_kpt)
581  ABI_DEALLOCATE(tmp_gkk_kpt2)
582  ABI_DEALLOCATE(gkk_kpt)
583  ABI_DEALLOCATE(a2f_2d)
584  ABI_DEALLOCATE(a2f_2d2)
585 
586 !output inv_tau_k and tau_k
587  fname = trim(elph_ds%elph_base_name) // '_INVTAUK'
588  if (open_file(fname,message,newunit=unit_invtau,status='unknown') /= 0) then
589    MSG_ERROR(message)
590  end if
591 
592 !print header to relaxation time file
593  write (unit_invtau,*) '# k-dep inverse of the relaxation time as a function of temperature.'
594  write (unit_invtau,*) '# '
595  write (unit_invtau,*) '# nkptirr= ', nkptirr, 'nband= ', nband
596  write (unit_invtau,*) '# number of temperatures=  ', ntemper
597  write (unit_invtau,*) '# tau [femtosecond^-1]     '
598 
599  fname = trim(elph_ds%elph_base_name) // '_TAUK'
600  if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
601    MSG_ERROR(message)
602  end if
603 
604 !print header to relaxation time file
605  write (unit_tau,*) '# k-dep relaxation time as a function of temperature.'
606  write (unit_tau,*) '# '
607  write (unit_tau,*) '# nkptirr= ', nkptirr, 'nband= ', nband
608  write (unit_tau,*) '# number of temperatures=  ', ntemper
609  write (unit_tau,*) '# tau [femtosecond]     '
610 
611  tau_k = zero
612  do itemp=1,ntemper  ! runs over termperature in K
613    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
614    write(unit_invtau,'(a,f16.8)') '# Temperature = ', Temp
615    write(unit_tau,'(a,f16.8)') '# Temperature = ', Temp
616    do isppol=1,nsppol
617      write(unit_invtau,'(a,i6)') '# For isppol = ', isppol
618      write(unit_tau,'(a,i6)') '# For isppol = ', isppol
619      do iFSkpt = 1,nkpt
620 !FIXME: check when tau_k is too small, whether there should be a phonon
621 !scattering or not, and should tau_k be zero or not.
622        do jband = 1,nband
623          if (abs(inv_tau_k(itemp,isppol,iFSkpt,jband)) < tol9) then
624            inv_tau_k(itemp,isppol,iFSkpt,jband) = zero
625            tau_k(itemp,isppol,iFSkpt,jband) = zero
626          else
627 !no need to *nkpt due to wtkirr, as we need /nkpt for the sum
628 !no need to *two_pi due to the missing prefactor in gkk (see mka2f_tr_lova)
629            inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband)*elph_ds%occ_factor
630            tau_k(itemp,isppol,iFSkpt,jband) = one/inv_tau_k(itemp,isppol,iFSkpt,jband)
631          end if
632        end do ! nband
633        write(unit_invtau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, '   kpt=', elph_ds%k_phon%kptirr(:,iFSkpt)
634        write(unit_invtau,'(100D16.8)') (inv_tau_k(itemp,isppol,iFSkpt,iband)*femto/chu_tau,iband=1,nband)
635        write(unit_tau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, '   kpt=', elph_ds%k_phon%kptirr(:,iFSkpt)
636        write(unit_tau,'(100D16.8)') (tau_k(itemp,isppol,iFSkpt,iband)*chu_tau/femto,iband=1,nband)
637      end do ! nkptirr
638      write(unit_invtau,*) ' '
639      write(unit_tau,*) ' '
640    end do ! nsppol
641    write(unit_invtau,*) ' '
642    write(unit_invtau,*) ' '
643    write(unit_tau,*) ' '
644    write(unit_tau,*) ' '
645  end do ! ntemper
646 
647 ! Only use the irred k for eigenGS and tau_k
648  ABI_ALLOCATE(tmp_eigenGS,(elph_ds%nband,elph_ds%k_phon%new_nkptirr,elph_ds%nsppol))
649 
650  do ikpt_irr = 1, new_nkptirr
651    tmp_eigenGS(:,ikpt_irr,:) = eigenGS(:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)
652    tmp_tau_k(:,:,ikpt_irr,:) = tau_k(:,:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)*chu_tau
653  end do
654 
655 !BoltzTraP output files in SIESTA format
656  if (elph_ds%prtbltztrp == 1) then
657    call ebands_prtbltztrp_tau_out (tmp_eigenGS(elph_ds%minFSband:elph_ds%maxFSband,:,:),&
658 &   elph_ds%tempermin,elph_ds%temperinc,ntemper,fermie, &
659 &   elph_ds%elph_base_name,elph_ds%k_phon%new_kptirr,nband,elph_ds%nelect,new_nkptirr, &
660 &   elph_ds%nspinor,nsppol,Cryst%nsym,Cryst%rprimd,Cryst%symrel,tmp_tau_k)
661  end if !prtbltztrp
662  ABI_DEALLOCATE(tmp_eigenGS)
663  ABI_DEALLOCATE(tmp_tau_k)
664 
665 !Get the energy dependence of tau.
666 !Eq. (6) in  Restrepo et al. Appl. Phys. Lett. 94, 212103 (2009)
667 
668  fname = trim(elph_ds%elph_base_name) // '_TAUE'
669  if (open_file(fname,message,newunit=unit_taue,status='unknown') /= 0) then
670    MSG_ERROR(message)
671  end if
672 
673 !print header to relaxation time file
674  write (unit_taue,*) '# Energy-dep relaxation time as a function of temperature.'
675  write (unit_taue,*) '# '
676  write (unit_taue,*) '# number of temperatures=  ', ntemper
677  write (unit_taue,*) '# ene[Ha] tau [femtosecond] DOS[au]    '
678 
679  fname = trim(elph_ds%elph_base_name) // '_MFP'
680  if (open_file(fname,message,newunit=unit_mfp,status='unknown') /= 0) then
681    MSG_ERROR(message)
682  end if
683 
684  write (unit_mfp,*) '# Energy-dep mean free path as a function of temperature.'
685  write (unit_mfp,*) '# '
686  write (unit_mfp,*) '# number of temperatures=  ', ntemper
687  write (unit_mfp,*) '# ene[Ha] mfp [femtometer]   '
688 
689  do itemp=1,ntemper  ! runs over termperature in K
690    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
691    write(unit_taue,'(a,f16.8)') '# Temperature = ', Temp
692    do isppol = 1, nsppol
693      write(unit_taue,*) '# Tau_e for isppol = ',isppol
694      do iene = 1, nene_all
695        rate_e = zero
696        do iFSkpt = 1, nkpt
697          do jband = 1, nband
698            rate_e = rate_e + inv_tau_k(itemp,isppol,iFSkpt,jband)* &
699 &           tmp_wtk(jband,iFSkpt,isppol,iene)
700          end do ! jband
701        end do ! kpt
702        if (dabs(dos_e(isppol,iene)) < tol7) then
703          rate_e = zero
704        else
705          rate_e = rate_e/nkpt/dos_e(isppol,iene)
706        end if
707        write(unit_taue,"(3D16.8)") enemin+(iene-1)*deltaene*nspline, rate_e*femto/chu_tau, dos_e(isppol,iene)
708      end do ! number of energies
709      write(unit_taue,*) ' '
710    end do ! nsppol
711    write(unit_taue,*) ' '
712  end do ! ntemperature
713 
714 ! calculate and output mean free path
715  do itemp=1,ntemper  ! runs over termperature in K
716    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
717    write(unit_mfp,'(a,f16.8)') '# Temperature = ', Temp
718    do isppol = 1, nsppol
719      do icomp = 1, 3
720        write(unit_mfp,*) '# Mean free path for isppol, icomp= ',isppol,icomp
721        do iene = 1, nene_all
722          mfp_e = zero
723          do iFSkpt = 1, nkptirr
724            do jband = 1, nband
725              mfp_e = mfp_e + tau_k(itemp,isppol,iFSkpt,jband)* &
726 &             elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* &
727 &             tmp_wtk(jband,iFSkpt,isppol,iene)
728 !&                          elph_ds%k_phon%new_wtkirr(iFSqpt)
729            end do ! jband
730          end do ! kpt
731          if (dabs(dos_e(isppol,iene)) < tol7) then
732            mfp_e = zero
733          else
734            mfp_e = mfp_e/nkptirr/dos_e(isppol,iene)
735          end if
736          write(unit_mfp,"(2D16.8)") enemin+(iene-1)*deltaene*nspline, mfp_e*chu_mfp/femto
737        end do ! number of energies
738        write(unit_mfp,*) ' '
739      end do ! icomp
740      write(unit_mfp,*) ' '
741    end do ! nsppol
742    write(unit_mfp,*) ' '
743  end do ! ntemperature
744 
745  ABI_ALLOCATE(cond_e ,(ntemper,nsppol,nene_all,9))
746 
747 !get cond_e
748  cond_e = zero
749  do itemp=1,ntemper  ! runs over termperature in K
750    do isppol = 1, nsppol
751      do iene = 1, nene_all
752 !       do iFSkpt =1,nkpt
753        do ik_this_proc =1,elph_ds%k_phon%my_nkpt
754          iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
755          do jband = 1, nband
756            do icomp = 1, 3
757              do jcomp = 1, 3
758                itensor = (icomp-1)*3+jcomp
759                cond_e(itemp,isppol,iene,itensor) = cond_e(itemp,isppol,iene,itensor) + &
760 &               tau_k(itemp,isppol,iFSkpt,jband)* &
761 &               elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* &
762 &               elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,jcomp,isppol)* &
763 &               tmp_wtk(jband,iFSkpt,isppol,iene)
764              end do
765            end do
766          end do ! jband
767        end do ! kpt
768      end do ! number of energies
769    end do ! nsppol
770  end do ! ntemperature
771 
772  ! MG FIXME: Why xmpi_world, besides only master should perform IO in the section below.
773  call xmpi_sum (cond_e, xmpi_world, ierr)
774 
775  cond_e = cond_e/nkpt
776 
777 !get transport coefficients
778 
779  fname = trim(elph_ds%elph_base_name) // '_COND'
780  if (open_file(fname,message,newunit=unit_cond,status='unknown') /= 0) then
781    MSG_ERROR(message)
782  end if
783 
784 !print header to conductivity file
785  write (unit_cond,*) '#  Conductivity as a function of temperature.'
786  write (unit_cond,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
787  write (unit_cond,*) '#  '
788  write (unit_cond,*) '#  Columns are: '
789  write (unit_cond,*) '#  temperature[K]   cond[au]   cond [SI]    '
790  write (unit_cond,*) '#  '
791 
792  fname = trim(elph_ds%elph_base_name) // '_CTH'
793  if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
794    MSG_ERROR(message)
795  end if
796 
797 !print header to thermal conductivity file
798  write (unit_therm,'(a)') '# Thermal conductivity as a function of temperature.'
799  write (unit_therm,'(a)') '#  the formalism is isotropic, so non-cubic crystals may be wrong'
800  write (unit_therm,'(a)') '#  '
801  write (unit_therm,'(a)') '#  Columns are: '
802  write (unit_therm,'(a)') '#  temperature[K]   thermal cond [au]   thermal cond [SI]'
803  write (unit_therm,'(a)') '#  '
804 
805  fname = trim(elph_ds%elph_base_name) // '_SBK'
806  if (open_file(fname,message,newunit=unit_sbk,status='unknown') /=0) then
807    MSG_ERROR(message)
808  end if
809 
810 !print header to relaxation time file
811  write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.'
812  write (unit_sbk,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
813  write (unit_sbk,*) '#  '
814  write (unit_sbk,*) '#  Columns are: '
815  write (unit_sbk,*) '#  temperature[K]   S [au]   S [SI]     '
816  write (unit_sbk,*) '#  '
817 
818  ABI_ALLOCATE(cond ,(ntemper,nsppol,3,3))
819  ABI_ALLOCATE(cth ,(ntemper,nsppol,3,3))
820  ABI_ALLOCATE(sbk ,(ntemper,nsppol,3,3))
821  ABI_ALLOCATE(seebeck ,(ntemper,nsppol,3,3))
822 
823  cond = zero
824  cth = zero
825  sbk = zero
826  seebeck = zero
827  do isppol=1,nsppol
828    do icomp=1, 3
829      do jcomp=1, 3
830        itensor=(icomp-1)*3+jcomp
831        do itemp=1,ntemper
832          Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp)
833          do iene = 1, nene_all
834            factor = (enemin+(iene-1)*deltaene*nspline - fermie(itemp))/(kb_HaK*Temp)
835            if (factor < -40.0d0) then
836              dfermide = zero
837            else if (factor > 40.0d0) then
838              dfermide = zero
839            else
840              dfermide = EXP(factor)/(kb_HaK*Temp*(EXP(factor)+one)**2.0d0)
841            end if
842            cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp) + &
843 &           cond_e(itemp,isppol,iene,itensor)*dfermide*deltaene*nspline
844            cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* &
845 &           (enemin+(iene-1)*deltaene*nspline - fermie(itemp))**2.0d0*dfermide*deltaene*nspline
846            sbk(itemp,isppol,icomp,jcomp) = sbk(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* &
847 &           (enemin+(iene-1)*deltaene*nspline - fermie(itemp))*dfermide*deltaene*nspline
848          end do
849        end do ! temperature
850      end do ! jcomp
851    end do ! icomp
852  end do !end isppol
853 
854  do isppol=1,nsppol
855    do itemp=1,ntemper
856      cond_inv(:,:)=cond(itemp,isppol,:,:)
857      call matrginv(cond_inv,3,3)
858      call DGEMM('N','N',3,3,3,one,sbk(itemp,isppol,:,:),3,cond_inv,&
859 &     3,zero,seebeck(itemp,isppol,:,:),3)
860    end do
861  end do
862 
863  do isppol=1,nsppol
864    do icomp=1, 3
865      do jcomp=1, 3
866        itensor=(icomp-1)*3+jcomp
867        write(unit_cond,*) '# Conductivity for isppol, itrten= ',isppol,itensor
868        write(unit_therm,*) '# Thermal conductivity for isppol, itrten= ',isppol,itensor
869        write(unit_sbk,*) '# Seebeck coefficient for isppol, itrten= ',isppol,itensor
870        do itemp=1,ntemper
871          Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp)
872 
873          seebeck(itemp,isppol,icomp,jcomp) = -1.0d0*seebeck(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)
874          cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp)/cryst%ucvol
875          cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)/cryst%ucvol
876          write(unit_cond,'(3D20.10)')Temp,cond(itemp,isppol,icomp,jcomp),cond(itemp,isppol,icomp,jcomp)*chu_cond
877          write(unit_therm,'(3D20.10)')Temp,cth(itemp,isppol,icomp,jcomp),cth(itemp,isppol,icomp,jcomp)*chu_cth
878          write(unit_sbk,'(3D20.10)')Temp,seebeck(itemp,isppol,icomp,jcomp),seebeck(itemp,isppol,icomp,jcomp)*chu_sbk
879        end do ! temperature
880        write(unit_cond,*)
881        write(unit_therm,*)
882        write(unit_sbk,*)
883      end do ! jcomp
884    end do ! icomp
885  end do !end isppol
886 
887 
888  ABI_DEALLOCATE(inv_tau_k)
889  ABI_DEALLOCATE(tau_k)
890  ABI_DEALLOCATE(tmp_wtk)
891  ABI_DEALLOCATE(dos_e)
892  ABI_DEALLOCATE(cond_e)
893  ABI_DEALLOCATE(fermie)
894  ABI_DEALLOCATE(cond)
895  ABI_DEALLOCATE(sbk)
896  ABI_DEALLOCATE(cth)
897  ABI_DEALLOCATE(seebeck)
898 
899  close (unit=unit_tau)
900  close (unit=unit_taue)
901  close (unit=unit_mfp)
902  close (unit=unit_invtau)
903  close (unit=unit_cond)
904  close (unit=unit_therm)
905  close (unit=unit_sbk)
906 
907 end subroutine get_tau_k