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) [[cite:Restrepo2009]])

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

1699 subroutine get_tau_k(Cryst,ifc,Bst,elph_ds,elph_tr_ds,eigenGS,max_occ)
1700 
1701 
1702 !This section has been created automatically by the script Abilint (TD).
1703 !Do not modify the following lines by hand.
1704 #undef ABI_FUNC
1705 #define ABI_FUNC 'get_tau_k'
1706 !End of the abilint section
1707 
1708  implicit none
1709 
1710 !Arguments ------------------------------------
1711  type(crystal_t),intent(in) :: Cryst
1712  type(ifc_type),intent(in) :: ifc
1713  type(ebands_t),intent(inout)   :: Bst
1714  type(elph_type),intent(inout) :: elph_ds
1715  type(elph_tr_type), intent(inout) :: elph_tr_ds
1716  real(dp),intent(in) :: max_occ
1717  real(dp),intent(in) :: eigenGS(elph_ds%nband,elph_ds%k_phon%nkpt,elph_ds%nsppol)
1718 
1719 !Local variables-------------------------------
1720 !scalars
1721  character(len=500) :: message
1722  character(len=fnlen) :: fname
1723  integer :: ntemper,nsppol,nbranch,nband,natom
1724  integer :: nkpt,nqpt,nkptirr,nqptirr,new_nkptirr
1725  integer :: isppol,iFSkpt,iFSqpt,iqpt,iqpt_fullbz,imqpt_fullbz,ikpt_kpq,ikpt_kmq
1726  integer :: iband,jband,jpband,jbeff,ibranch,jbranch,itemp
1727  integer :: irec,ierr,nrpt,ik_this_proc
1728  integer :: unit_tau,unit_invtau
1729  integer :: nene,nene_all,iene,iene_fine,unit_taue,unit_mfp
1730  integer :: icomp,jcomp,itensor
1731  integer :: ikpt_irr,iomega,unit_cond,unit_therm,unit_sbk
1732  integer :: nskip,nspline
1733  real(dp) :: occ_omega,occ_e
1734  real(dp) :: xx,Temp,therm_factor
1735  real(dp) :: factor,dfermide
1736  real(dp) :: e_k,chu_tau,rate_e,mfp_e
1737  real(dp) :: ene,enemin,enemax,deltaene
1738  real(dp) :: omega,omega_min,omega_max,domega
1739  real(dp) :: diagerr
1740  real(dp) :: chu_mfp,chu_cond,chu_cth,chu_sbk,femto
1741  real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1742  real(dp) :: eigval(elph_ds%nbranch),eigval2(elph_ds%nbranch)
1743  real(dp) :: imeigval(elph_ds%nbranch)
1744  real(dp) :: tmp_wtkpq, tmp_wtkmq, tol_wtk
1745  real(dp) :: yp1,ypn
1746 !arrays
1747  integer,allocatable :: FSfullpktofull(:,:),mqtofull(:)
1748  integer,allocatable :: kpttokpt(:,:,:)
1749  real(dp) :: cond_inv(3,3)
1750  real(dp),allocatable :: fermie(:)
1751  real(dp),allocatable :: tmp_eigenGS(:,:,:)
1752  real(dp),allocatable :: tmp_gkk_qpt(:,:,:),tmp_gkk_rpt(:,:,:),tmp_gkk_kpt(:,:)
1753  real(dp),allocatable :: tmp_gkk_kpt2(:,:,:), gkk_kpt(:,:,:)
1754  real(dp),allocatable :: tau_k(:,:,:,:),inv_tau_k(:,:,:,:),tmp_tau_k(:,:,:,:)
1755  real(dp),allocatable :: phfrq(:,:),pheigvec(:,:)
1756  real(dp),allocatable :: displ(:,:,:,:)
1757  real(dp),allocatable :: a2f_2d(:),a2f_2d2(:)
1758  real(dp),allocatable :: tmp_wtk(:,:,:,:),tmp2_wtk(:),tmp_wtk1(:),tmp_wtk2(:)
1759  real(dp),allocatable :: ene_pt(:),ene_ptfine(:),ff2(:)
1760  real(dp),allocatable :: wtq(:,:,:),tmp_wtq(:,:,:),tmp2_wtq(:,:)
1761  real(dp),allocatable :: dos_e(:,:)
1762  real(dp),allocatable :: coskr1(:,:),sinkr1(:,:)
1763  real(dp),allocatable :: coskr2(:,:),sinkr2(:,:)
1764  real(dp),allocatable :: cond_e(:,:,:,:),cond(:,:,:,:),sbk(:,:,:,:),seebeck(:,:,:,:),cth(:,:,:,:)
1765 
1766 ! *************************************************************************
1767 
1768  write(std_out,*) 'get_tau_k : enter '
1769 
1770  nrpt = ifc%nrpt
1771  natom = cryst%natom
1772 
1773  nsppol   = elph_ds%nsppol
1774  nbranch  = elph_ds%nbranch
1775  nband    = elph_ds%ngkkband
1776  nkpt     = elph_ds%k_phon%nkpt
1777  nqpt     = elph_ds%nqpt_full
1778  nkptirr  = elph_ds%k_phon%nkptirr
1779  new_nkptirr  = elph_ds%k_phon%new_nkptirr
1780  nqptirr  = elph_ds%nqptirred
1781  ntemper  = elph_ds%ntemper
1782  nene = 2*elph_ds%na2f-1 ! only need e_k +- omega_max range, take deltaene=delta_oemga
1783 
1784  chu_tau  = 2.4188843265*1.0d-17
1785  chu_mfp  = 5.291772*1.0d-11
1786  chu_cond = 4.59988159904764*1.0d6
1787  chu_cth  = 1.078637439971599*1.0d4
1788  chu_sbk  = 8.617343101*1.0d-5
1789  femto    = 1.0d-15
1790 
1791  tol_wtk = tol7/nkptirr/nband
1792 
1793  ABI_ALLOCATE(fermie ,(ntemper))
1794  ABI_ALLOCATE(tmp_gkk_qpt ,(2,nbranch**2,nqpt))
1795  ABI_ALLOCATE(tmp_gkk_rpt ,(2,nbranch**2,nrpt))
1796  ABI_ALLOCATE(tmp_gkk_kpt ,(2,nbranch**2))
1797  ABI_ALLOCATE(tmp_gkk_kpt2 ,(2,nbranch,nbranch))
1798  ABI_ALLOCATE(gkk_kpt ,(2,nbranch,nbranch))
1799  ABI_ALLOCATE(a2f_2d, (nene))
1800  ABI_ALLOCATE(a2f_2d2, (nene))
1801  ABI_ALLOCATE(inv_tau_k, (ntemper,nsppol,nkpt,nband))
1802  ABI_ALLOCATE(tau_k, (ntemper,nsppol,nkpt,nband))
1803  ABI_ALLOCATE(tmp_tau_k ,(ntemper,nsppol,new_nkptirr,nband))
1804 
1805  if (elph_ds%gkqwrite == 0) then
1806    call wrtout(std_out,' get_tau_k : keeping gkq matrices in memory','COLL')
1807  else if (elph_ds%gkqwrite == 1) then
1808    fname=trim(elph_ds%elph_base_name) // '_GKKQ'
1809    write (message,'(2a)')' get_tau_k : reading gkq matrices from file ',trim(fname)
1810    call wrtout(std_out,message,'COLL')
1811  else
1812    write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite
1813    MSG_BUG(message)
1814  end if
1815 
1816 !=========================================================
1817 !Get equivalence between a kpt_phon pair and a qpt in qpt_full
1818 !only works if the qpt grid is complete (identical to
1819 !the kpt one, with a basic shift of (0,0,0)
1820 !=========================================================
1821 
1822 !mapping of k + q onto k' for k and k' in full BZ
1823 !for dense k grid
1824  ABI_ALLOCATE(FSfullpktofull,(nkpt,nkpt))
1825  ABI_ALLOCATE(mqtofull,(nkpt))
1826 
1827 !kpttokpt(itim,isym,iqpt) = kpoint index which transforms to ikpt under isym and with time reversal itim.
1828  ABI_ALLOCATE(kpttokpt,(2,Cryst%nsym,nkpt))
1829 
1830  call wrtout(std_out,'get_tau_k: calling mkqptequiv to set up the FS kpoint set',"COLL")
1831 
1832  call mkqptequiv (FSfullpktofull,Cryst,elph_ds%k_phon%kpt,nkpt,nkpt,kpttokpt,elph_ds%k_phon%kpt,mqtofull)
1833 
1834 !=========================================================
1835 !=========================================================
1836 
1837  omega_max       = elph_ds%omega_max
1838  omega_min       = elph_ds%omega_min
1839  domega          = elph_ds%domega
1840  enemax = maxval(eigenGS(elph_ds%maxFSband,:,:))
1841  enemin = minval(eigenGS(elph_ds%minFSband,:,:))
1842 
1843  if (enemin < (elph_ds%fermie-0.2)) then
1844    enemin = elph_ds%fermie-0.2
1845  end if
1846  if (enemax > (elph_ds%fermie+0.2)) then
1847    enemax = elph_ds%fermie+0.2
1848  end if
1849 
1850  nspline = elph_ds%ep_nspline
1851  nene_all = INT((enemax-enemin+domega)/(nspline*domega)) + 1
1852  deltaene = domega
1853  write(std_out,*) 'E_min= ',enemin, 'E_max= ',enemax
1854  write(std_out,*) 'Number of energy points= ',nene_all
1855  write(std_out,'(a,I8)') 'scale factor for spline interpolation in RTA = ', elph_ds%ep_nspline
1856  write(std_out,*) 'delta_ene before spline interpolation= ',deltaene*nspline
1857  write(std_out,*) 'delta_ene after spline interpolation= ',deltaene
1858  write(std_out,*) 'Omega_min= ',omega_min, 'Omega_max= ',omega_max
1859  write(std_out,*) 'Number of phonon points= ',elph_ds%na2f
1860  write(std_out,*) 'delta_omega= ',domega
1861  write(std_out,*) 'number of bands= ', elph_ds%nband, nband
1862 
1863  ABI_ALLOCATE(tmp_wtk,(nband,nkpt,nsppol,nene_all))
1864  ABI_ALLOCATE(tmp2_wtk,(nene_all))
1865  ABI_ALLOCATE(ff2,(nene_all))
1866  ABI_ALLOCATE(ene_pt,(nene_all))
1867  ABI_ALLOCATE(ene_ptfine,(nene_all*nspline))
1868  ABI_ALLOCATE(tmp_wtk1,(nene_all*nspline))
1869  ABI_ALLOCATE(tmp_wtk2,(nene_all*nspline))
1870  ABI_ALLOCATE(dos_e,(nsppol,nene_all))
1871 
1872 !Get energy points for spline interpolation
1873  do iene = 1, nene_all
1874    ene_pt(iene) = enemin + (iene-1)*nspline*deltaene
1875  end do
1876 
1877  do iene = 1, nene_all*nspline
1878    ene_ptfine(iene) = enemin + (iene-1)*deltaene
1879  end do
1880 
1881  ABI_ALLOCATE(tmp_wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f+1))
1882  ABI_ALLOCATE(wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f))
1883  ABI_ALLOCATE(tmp2_wtq,(elph_ds%nbranch, elph_ds%na2f))
1884 
1885 !phonon
1886  ABI_ALLOCATE(phfrq,(nbranch, nkptirr))
1887  ABI_ALLOCATE(displ,(2, nbranch, nbranch, nkptirr))
1888  ABI_ALLOCATE(pheigvec,(2*nbranch*nbranch, nkptirr))
1889 
1890  do iFSqpt = 1, nkptirr
1891    call ifc_fourq(ifc,cryst,elph_ds%k_phon%kptirr(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
1892  end do
1893 
1894  omega_min = omega_min - domega
1895 
1896 !bxu, obtain wtq for the q_fine, then condense to q_phon
1897  call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,Cryst%gprimd,elph_ds%kptrlatt, &
1898 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_phon,tmp_wtq)
1899  omega_min = omega_min + domega
1900 
1901  do iomega = 1, elph_ds%na2f
1902    wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1)
1903    !write(1005,*) omega_min+(iomega-1)*domega, sum(tmp_wtq(:,:,iomega+1))/nkpt
1904  end do
1905  ABI_DEALLOCATE(tmp_wtq)
1906 
1907 ! electron
1908  tmp_wtk =zero
1909  dos_e = zero
1910  call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS(elph_ds%minFSband:elph_ds%minFSband+nband-1,:,:), &
1911 & elph_ds%elphsmear, &
1912 & enemin, enemax, nene_all, Cryst%gprimd, elph_ds%k_phon%irredtoGS, elph_ds%kptrlatt, max_occ, &
1913 & 1, nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk)
1914 !& elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk)
1915 
1916  do isppol = 1, nsppol
1917    do iene = 1, nene_all
1918      dos_e(isppol,iene) = sum(tmp_wtk(:,:,isppol,iene))/nkpt
1919    end do
1920  end do
1921 
1922  ABI_ALLOCATE(coskr1, (nqpt,nrpt))
1923  ABI_ALLOCATE(sinkr1, (nqpt,nrpt))
1924  call ftgam_init(ifc%gprim, nqpt, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr1, sinkr1)
1925  ABI_ALLOCATE(coskr2, (nkptirr,nrpt))
1926  ABI_ALLOCATE(sinkr2, (nkptirr,nrpt))
1927  call ftgam_init(ifc%gprim, nkptirr, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr2, sinkr2)
1928 
1929 !get fermie for itemp
1930  fermie = elph_ds%fermie
1931  do itemp=1,ntemper  ! runs over termperature in K
1932    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
1933 
1934    Bst%occopt = 3
1935    Bst%tsmear = Temp*kb_HaK
1936    call ebands_update_occ(Bst,-99.99_dp)
1937    write(message,'(a,f12.6,a,E20.12)')'At T=',Temp,' Fermi level is:',Bst%fermie
1938    call wrtout(std_out,message,'COLL')
1939 
1940    if (abs(elph_ds%fermie) < tol10) then
1941      fermie(itemp) = Bst%fermie
1942    end if
1943  end do
1944 
1945  inv_tau_k = zero
1946 !get a2f_2d = \sum_{q,nbranch,jband'} |gkk|^2*\delta(\epsilon_{k'j'}-\epsilon')*\delta(\omega_q-\omega)
1947  do isppol=1,nsppol
1948    write (std_out,*) '##############################################'
1949    write (std_out,*) 'get_tau_k : Treating spin polarization ', isppol
1950    write (std_out,*) '##############################################'
1951 
1952 !   do iFSkpt =1,nkpt
1953    do ik_this_proc =1,elph_ds%k_phon%my_nkpt
1954      iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
1955      write (std_out,*) 'get_tau_k : working on kpt # ', iFSkpt, '/', nkpt
1956      do jband = 1, nband
1957 !          write(*,*)'i am here 1 ', isppol,iFSkpt,jband
1958        a2f_2d = zero
1959        a2f_2d2 = zero
1960 
1961 !sum from here
1962        nskip = 0
1963        do jpband = 1, nband
1964          jbeff = jpband+(jband-1)*nband
1965 
1966          if (elph_ds%gkqwrite == 0) then
1967            tmp_gkk_qpt(:,:,:) = elph_ds%gkk_qpt(:,jbeff,:,ik_this_proc,isppol,:)
1968          else if (elph_ds%gkqwrite == 1) then
1969            irec = (ik_this_proc-1)*elph_ds%k_phon%my_nkpt + iqpt
1970            if (iFSkpt == 1) then
1971              write (std_out,*) ' get_tau_k  read record ', irec
1972            end if
1973            read (elph_ds%unitgkq,REC=irec) tmp_gkk_qpt(:,:,iqpt_fullbz)
1974          end if
1975 
1976 !FT to real space
1977          call ftgam(Ifc%wghatm,tmp_gkk_qpt,tmp_gkk_rpt,natom,nqpt,nrpt,1,coskr1,sinkr1)
1978 
1979 !sum over irred q over k_phon, with corresponding weights
1980          do iFSqpt = 1, nkptirr
1981            iqpt_fullbz = elph_ds%k_phon%irredtoGS(iFSqpt)
1982            ikpt_kpq = FSfullpktofull(iFSkpt,iqpt_fullbz)
1983 
1984            imqpt_fullbz = mqtofull(iqpt_fullbz)
1985            ikpt_kmq = FSfullpktofull(iFSkpt,imqpt_fullbz)
1986 
1987 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr
1988            call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(iqpt_fullbz,:),sinkr2(iqpt_fullbz,:))
1989 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt)
1990 
1991 !if ep_scalprod==0 we have to dot in the displacement vectors here
1992            if (elph_ds%ep_scalprod==0) then
1993 
1994              call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red)
1995 
1996              tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/))
1997              call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt)
1998 
1999              do jbranch=1,nbranch
2000                eigval(jbranch) = gkk_kpt(1, jbranch, jbranch)
2001                imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch)
2002 
2003                if (abs(imeigval(jbranch)) > tol10) then
2004                  write (message,'(a,i0,a,es16.8)')" real values  branch = ",jbranch,' eigval = ',eigval(jbranch)
2005                  MSG_WARNING(message)
2006                  write (message,'(a,i0,a,es16.8)')" imaginary values  branch = ",jbranch,' imeigval = ',imeigval(jbranch)
2007                  MSG_WARNING(message)
2008                end if
2009 
2010              end do
2011 
2012 !            if ep_scalprod==1 we have to diagonalize the matrix we interpolated.
2013            else if (elph_ds%ep_scalprod == 1) then
2014 
2015 !            MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here
2016              call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,&
2017 &             pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom)
2018 
2019              call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,&
2020 &             tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom)
2021 
2022              diagerr = zero
2023              do ibranch=1,nbranch
2024                eigval(ibranch) = gkk_kpt(1,ibranch,ibranch)
2025                do jbranch=1,ibranch-1
2026                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
2027                end do
2028                do jbranch=ibranch+1,nbranch
2029                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
2030                end do
2031              end do
2032 
2033              if (diagerr > tol12) then
2034                write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr
2035                MSG_WARNING(message)
2036              end if
2037 
2038            else
2039              write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod
2040              MSG_BUG(message)
2041            end if ! end ep_scalprod if
2042 
2043 !For k'=k-q
2044 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr
2045            call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(imqpt_fullbz,:),sinkr2(imqpt_fullbz,:))
2046 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt)
2047 
2048 !if ep_scalprod==0 we have to dot in the displacement vectors here
2049            if (elph_ds%ep_scalprod==0) then
2050 
2051              call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red)
2052 
2053              tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/))
2054              call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt)
2055 
2056              do jbranch=1,nbranch
2057                eigval2(jbranch) = gkk_kpt(1, jbranch, jbranch)
2058                imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch)
2059 
2060                if (abs(imeigval(jbranch)) > tol10) then
2061                  write (message,'(a,i0,a,es16.8)')" real values  branch = ",jbranch,' eigval = ',eigval2(jbranch)
2062                  MSG_WARNING(message)
2063                  write (message,'(a,i0,a,es16.8)')" imaginary values  branch = ",jbranch,' imeigval = ',imeigval(jbranch)
2064                  MSG_WARNING(message)
2065                end if
2066 
2067              end do
2068 
2069 !            if ep_scalprod==1 we have to diagonalize the matrix we interpolated.
2070            else if (elph_ds%ep_scalprod == 1) then
2071 
2072 !            MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here
2073              call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,&
2074 &             pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom)
2075 
2076              call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,&
2077 &             tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom)
2078 
2079              diagerr = zero
2080              do ibranch=1,nbranch
2081                eigval2(ibranch) = gkk_kpt(1,ibranch,ibranch)
2082                do jbranch=1,ibranch-1
2083                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
2084                end do
2085                do jbranch=ibranch+1,nbranch
2086                  diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch))
2087                end do
2088              end do
2089 
2090              if (diagerr > tol12) then
2091                write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr
2092                MSG_WARNING(message)
2093              end if
2094 
2095            else
2096              write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod
2097              MSG_BUG(message)
2098            end if ! end ep_scalprod if
2099 
2100            tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kpq,isppol,:)
2101            yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene
2102            ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene
2103            call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2)
2104            call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk1)
2105 
2106            tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kmq,isppol,:)
2107            yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene
2108            ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene
2109            call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2)
2110            call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk2)
2111 
2112            tmp2_wtq(:,:) = wtq(:,iFSqpt,:)
2113            do iene=1,nene
2114              e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol)
2115              ene = e_k - omega_max + (iene-1)*deltaene
2116              if (ene<enemin .or. ene>enemax) cycle
2117              iene_fine = NINT((ene-enemin+deltaene)/deltaene)
2118              tmp_wtkpq = tmp_wtk1(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt)
2119              tmp_wtkmq = tmp_wtk2(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt)
2120 
2121              if (tmp_wtkpq+tmp_wtkmq < tol_wtk ) then
2122                nskip = nskip +1
2123                cycle
2124              end if
2125 
2126              do ibranch = 1, nbranch
2127                if (abs(phfrq(ibranch,iFSqpt)) < tol7) cycle
2128 
2129                if (ene > e_k) then
2130                  omega = ene - e_k
2131                  if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle
2132                  iomega = NINT((omega-omega_min+domega)/domega)
2133 
2134                  a2f_2d(iene) = a2f_2d(iene) +&
2135 &                 eigval(ibranch)/phfrq(ibranch,iFSqpt)*&
2136 &                 tmp_wtkpq * tmp2_wtq(ibranch,iomega)
2137                end if
2138 
2139                if (ene < e_k) then
2140                  omega = e_k - ene
2141                  if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle
2142                  iomega = NINT((omega-omega_min+domega)/domega)
2143 
2144                  a2f_2d2(iene) = a2f_2d2(iene) +&
2145 &                 eigval(ibranch)/phfrq(ibranch,iFSqpt)*&
2146 &                 tmp_wtkmq * tmp2_wtq(ibranch,iomega)
2147                end if
2148 
2149              end do ! ibranch 3
2150            end do ! nene  800
2151          end do ! kptirr 216
2152        end do ! j' band 3
2153 !      print *, ' skipped ',  nskip, ' energy points out of ', nene*nband*nkptirr
2154 
2155 ! get inv_tau_k
2156        do itemp=1,ntemper  ! runs over termperature in K
2157          Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2158          do iene=1,nene
2159            e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol)
2160            ene = e_k - omega_max + (iene-1)*deltaene
2161            if (ene<enemin .or. ene>enemax) cycle
2162 
2163            xx=(ene-fermie(itemp))/(kb_HaK*Temp)
2164            occ_e=1.0_dp/(exp(xx)+1.0_dp)
2165            if (ene > e_k .and. (ene-e_k) .le. omega_max) then
2166              omega = ene - e_k
2167              if (abs(omega) < tol7) cycle
2168              xx = omega/(kb_HaK*Temp)
2169              occ_omega=1.0_dp/(exp(xx)-1.0_dp)
2170 
2171              therm_factor = occ_e + occ_omega
2172 
2173              inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +&
2174              a2f_2d(iene)*therm_factor*deltaene
2175            end if
2176            if (ene < e_k .and. (e_k-ene) .le. omega_max) then
2177              omega = e_k - ene
2178              if (abs(omega) < tol7) cycle
2179              xx = omega/(kb_HaK*Temp)
2180              occ_omega=1.0_dp/(exp(xx)-1.0_dp)
2181 
2182              therm_factor = 1 - occ_e + occ_omega
2183 
2184              inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +&
2185              a2f_2d2(iene)*therm_factor*deltaene
2186            end if
2187 
2188          end do ! nene
2189        end do ! Temp
2190 !          write(*,*)'i am here 2 ', isppol,iFSkpt,jband
2191      end do ! jband
2192    end do ! kpt
2193  end do ! nsppol
2194 
2195 !write (300+mpi_enreg%me,*) inv_tau_k
2196  call xmpi_sum (inv_tau_k, xmpi_world, ierr)
2197 
2198  ABI_DEALLOCATE(phfrq)
2199  ABI_DEALLOCATE(displ)
2200  ABI_DEALLOCATE(pheigvec)
2201  ABI_DEALLOCATE(tmp2_wtk)
2202  ABI_DEALLOCATE(ff2)
2203  ABI_DEALLOCATE(ene_pt)
2204  ABI_DEALLOCATE(ene_ptfine)
2205  ABI_DEALLOCATE(tmp_wtk1)
2206  ABI_DEALLOCATE(tmp_wtk2)
2207  ABI_DEALLOCATE(tmp2_wtq)
2208  ABI_DEALLOCATE(wtq)
2209  ABI_DEALLOCATE(coskr1)
2210  ABI_DEALLOCATE(sinkr1)
2211  ABI_DEALLOCATE(coskr2)
2212  ABI_DEALLOCATE(sinkr2)
2213  ABI_DEALLOCATE(kpttokpt)
2214  ABI_DEALLOCATE(FSfullpktofull)
2215  ABI_DEALLOCATE(mqtofull)
2216  ABI_DEALLOCATE(tmp_gkk_qpt)
2217  ABI_DEALLOCATE(tmp_gkk_rpt)
2218  ABI_DEALLOCATE(tmp_gkk_kpt)
2219  ABI_DEALLOCATE(tmp_gkk_kpt2)
2220  ABI_DEALLOCATE(gkk_kpt)
2221  ABI_DEALLOCATE(a2f_2d)
2222  ABI_DEALLOCATE(a2f_2d2)
2223 
2224 !output inv_tau_k and tau_k
2225  fname = trim(elph_ds%elph_base_name) // '_INVTAUK'
2226  if (open_file(fname,message,newunit=unit_invtau,status='unknown') /= 0) then
2227    MSG_ERROR(message)
2228  end if
2229 
2230 !print header to relaxation time file
2231  write (unit_invtau,*) '# k-dep inverse of the relaxation time as a function of temperature.'
2232  write (unit_invtau,*) '# '
2233  write (unit_invtau,*) '# nkptirr= ', nkptirr, 'nband= ', nband
2234  write (unit_invtau,*) '# number of temperatures=  ', ntemper
2235  write (unit_invtau,*) '# tau [femtosecond^-1]     '
2236 
2237  fname = trim(elph_ds%elph_base_name) // '_TAUK'
2238  if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
2239    MSG_ERROR(message)
2240  end if
2241 
2242 !print header to relaxation time file
2243  write (unit_tau,*) '# k-dep relaxation time as a function of temperature.'
2244  write (unit_tau,*) '# '
2245  write (unit_tau,*) '# nkptirr= ', nkptirr, 'nband= ', nband
2246  write (unit_tau,*) '# number of temperatures=  ', ntemper
2247  write (unit_tau,*) '# tau [femtosecond]     '
2248 
2249  tau_k = zero
2250  do itemp=1,ntemper  ! runs over termperature in K
2251    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2252    write(unit_invtau,'(a,f16.8)') '# Temperature = ', Temp
2253    write(unit_tau,'(a,f16.8)') '# Temperature = ', Temp
2254    do isppol=1,nsppol
2255      write(unit_invtau,'(a,i6)') '# For isppol = ', isppol
2256      write(unit_tau,'(a,i6)') '# For isppol = ', isppol
2257      do iFSkpt = 1,nkpt
2258 !FIXME: check when tau_k is too small, whether there should be a phonon
2259 !scattering or not, and should tau_k be zero or not.
2260        do jband = 1,nband
2261          if (abs(inv_tau_k(itemp,isppol,iFSkpt,jband)) < tol9) then
2262            inv_tau_k(itemp,isppol,iFSkpt,jband) = zero
2263            tau_k(itemp,isppol,iFSkpt,jband) = zero
2264          else
2265 !no need to *nkpt due to wtkirr, as we need /nkpt for the sum
2266 !no need to *two_pi due to the missing prefactor in gkk (see mka2f_tr_lova)
2267            inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband)*elph_ds%occ_factor
2268            tau_k(itemp,isppol,iFSkpt,jband) = one/inv_tau_k(itemp,isppol,iFSkpt,jband)
2269          end if
2270        end do ! nband
2271        write(unit_invtau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, '   kpt=', elph_ds%k_phon%kptirr(:,iFSkpt)
2272        write(unit_invtau,'(100D16.8)') (inv_tau_k(itemp,isppol,iFSkpt,iband)*femto/chu_tau,iband=1,nband)
2273        write(unit_tau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, '   kpt=', elph_ds%k_phon%kptirr(:,iFSkpt)
2274        write(unit_tau,'(100D16.8)') (tau_k(itemp,isppol,iFSkpt,iband)*chu_tau/femto,iband=1,nband)
2275      end do ! nkptirr
2276      write(unit_invtau,*) ' '
2277      write(unit_tau,*) ' '
2278    end do ! nsppol
2279    write(unit_invtau,*) ' '
2280    write(unit_invtau,*) ' '
2281    write(unit_tau,*) ' '
2282    write(unit_tau,*) ' '
2283  end do ! ntemper
2284 
2285 ! Only use the irred k for eigenGS and tau_k
2286  ABI_ALLOCATE(tmp_eigenGS,(elph_ds%nband,elph_ds%k_phon%new_nkptirr,elph_ds%nsppol))
2287 
2288  do ikpt_irr = 1, new_nkptirr
2289    tmp_eigenGS(:,ikpt_irr,:) = eigenGS(:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)
2290    tmp_tau_k(:,:,ikpt_irr,:) = tau_k(:,:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)*chu_tau
2291  end do
2292 
2293 !BoltzTraP output files in SIESTA format
2294  if (elph_ds%prtbltztrp == 1) then
2295    call ebands_prtbltztrp_tau_out (tmp_eigenGS(elph_ds%minFSband:elph_ds%maxFSband,:,:),&
2296 &   elph_ds%tempermin,elph_ds%temperinc,ntemper,fermie, &
2297 &   elph_ds%elph_base_name,elph_ds%k_phon%new_kptirr,natom,nband,elph_ds%nelect,new_nkptirr, &
2298 &   elph_ds%nspinor,nsppol,Cryst%nsym,Cryst%rprimd,Cryst%symrel,tmp_tau_k)
2299  end if !prtbltztrp
2300  ABI_DEALLOCATE(tmp_eigenGS)
2301  ABI_DEALLOCATE(tmp_tau_k)
2302 
2303 !Get the energy dependence of tau.
2304 !Eq. (6) in  Restrepo et al. Appl. Phys. Lett. 94, 212103 (2009) [[cite:Restrepo2009]]
2305 
2306  fname = trim(elph_ds%elph_base_name) // '_TAUE'
2307  if (open_file(fname,message,newunit=unit_taue,status='unknown') /= 0) then
2308    MSG_ERROR(message)
2309  end if
2310 
2311 !print header to relaxation time file
2312  write (unit_taue,*) '# Energy-dep relaxation time as a function of temperature.'
2313  write (unit_taue,*) '# '
2314  write (unit_taue,*) '# number of temperatures=  ', ntemper
2315  write (unit_taue,*) '# ene[Ha] tau [femtosecond] DOS[au]    '
2316 
2317  fname = trim(elph_ds%elph_base_name) // '_MFP'
2318  if (open_file(fname,message,newunit=unit_mfp,status='unknown') /= 0) then
2319    MSG_ERROR(message)
2320  end if
2321 
2322  write (unit_mfp,*) '# Energy-dep mean free path as a function of temperature.'
2323  write (unit_mfp,*) '# '
2324  write (unit_mfp,*) '# number of temperatures=  ', ntemper
2325  write (unit_mfp,*) '# ene[Ha] mfp [femtometer]   '
2326 
2327  do itemp=1,ntemper  ! runs over termperature in K
2328    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2329    write(unit_taue,'(a,f16.8)') '# Temperature = ', Temp
2330    do isppol = 1, nsppol
2331      write(unit_taue,*) '# Tau_e for isppol = ',isppol
2332      do iene = 1, nene_all
2333        rate_e = zero
2334        do iFSkpt = 1, nkpt
2335          do jband = 1, nband
2336            rate_e = rate_e + inv_tau_k(itemp,isppol,iFSkpt,jband)* &
2337 &           tmp_wtk(jband,iFSkpt,isppol,iene)
2338          end do ! jband
2339        end do ! kpt
2340        if (dabs(dos_e(isppol,iene)) < tol7) then
2341          rate_e = zero
2342        else
2343          rate_e = rate_e/nkpt/dos_e(isppol,iene)
2344        end if
2345        write(unit_taue,"(3D16.8)") enemin+(iene-1)*deltaene*nspline, rate_e*femto/chu_tau, dos_e(isppol,iene)
2346      end do ! number of energies
2347      write(unit_taue,*) ' '
2348    end do ! nsppol
2349    write(unit_taue,*) ' '
2350  end do ! ntemperature
2351 
2352 ! calculate and output mean free path
2353  do itemp=1,ntemper  ! runs over termperature in K
2354    Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp)
2355    write(unit_mfp,'(a,f16.8)') '# Temperature = ', Temp
2356    do isppol = 1, nsppol
2357      do icomp = 1, 3
2358        write(unit_mfp,*) '# Mean free path for isppol, icomp= ',isppol,icomp
2359        do iene = 1, nene_all
2360          mfp_e = zero
2361          do iFSkpt = 1, nkptirr
2362            do jband = 1, nband
2363              mfp_e = mfp_e + tau_k(itemp,isppol,iFSkpt,jband)* &
2364 &             elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* &
2365 &             tmp_wtk(jband,iFSkpt,isppol,iene)
2366 !&                          elph_ds%k_phon%new_wtkirr(iFSqpt)
2367            end do ! jband
2368          end do ! kpt
2369          if (dabs(dos_e(isppol,iene)) < tol7) then
2370            mfp_e = zero
2371          else
2372            mfp_e = mfp_e/nkptirr/dos_e(isppol,iene)
2373          end if
2374          write(unit_mfp,"(2D16.8)") enemin+(iene-1)*deltaene*nspline, mfp_e*chu_mfp/femto
2375        end do ! number of energies
2376        write(unit_mfp,*) ' '
2377      end do ! icomp
2378      write(unit_mfp,*) ' '
2379    end do ! nsppol
2380    write(unit_mfp,*) ' '
2381  end do ! ntemperature
2382 
2383  ABI_ALLOCATE(cond_e ,(ntemper,nsppol,nene_all,9))
2384 
2385 !get cond_e
2386  cond_e = zero
2387  do itemp=1,ntemper  ! runs over termperature in K
2388    do isppol = 1, nsppol
2389      do iene = 1, nene_all
2390 !       do iFSkpt =1,nkpt
2391        do ik_this_proc =1,elph_ds%k_phon%my_nkpt
2392          iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
2393          do jband = 1, nband
2394            do icomp = 1, 3
2395              do jcomp = 1, 3
2396                itensor = (icomp-1)*3+jcomp
2397                cond_e(itemp,isppol,iene,itensor) = cond_e(itemp,isppol,iene,itensor) + &
2398 &               tau_k(itemp,isppol,iFSkpt,jband)* &
2399 &               elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* &
2400 &               elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,jcomp,isppol)* &
2401 &               tmp_wtk(jband,iFSkpt,isppol,iene)
2402              end do
2403            end do
2404          end do ! jband
2405        end do ! kpt
2406      end do ! number of energies
2407    end do ! nsppol
2408  end do ! ntemperature
2409 
2410  ! MG FIXME: Why xmpi_world, besides only master should perform IO in the section below.
2411  call xmpi_sum (cond_e, xmpi_world, ierr)
2412 
2413  cond_e = cond_e/nkpt
2414 
2415 !get transport coefficients
2416 
2417  fname = trim(elph_ds%elph_base_name) // '_COND'
2418  if (open_file(fname,message,newunit=unit_cond,status='unknown') /= 0) then
2419    MSG_ERROR(message)
2420  end if
2421 
2422 !print header to conductivity file
2423  write (unit_cond,*) '#  Conductivity as a function of temperature.'
2424  write (unit_cond,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
2425  write (unit_cond,*) '#  '
2426  write (unit_cond,*) '#  Columns are: '
2427  write (unit_cond,*) '#  temperature[K]   cond[au]   cond [SI]    '
2428  write (unit_cond,*) '#  '
2429 
2430  fname = trim(elph_ds%elph_base_name) // '_CTH'
2431  if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
2432    MSG_ERROR(message)
2433  end if
2434 
2435 !print header to thermal conductivity file
2436  write (unit_therm,'(a)') '# Thermal conductivity as a function of temperature.'
2437  write (unit_therm,'(a)') '#  the formalism is isotropic, so non-cubic crystals may be wrong'
2438  write (unit_therm,'(a)') '#  '
2439  write (unit_therm,'(a)') '#  Columns are: '
2440  write (unit_therm,'(a)') '#  temperature[K]   thermal cond [au]   thermal cond [SI]'
2441  write (unit_therm,'(a)') '#  '
2442 
2443  fname = trim(elph_ds%elph_base_name) // '_SBK'
2444  if (open_file(fname,message,newunit=unit_sbk,status='unknown') /=0) then
2445    MSG_ERROR(message)
2446  end if
2447 
2448 !print header to relaxation time file
2449  write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.'
2450  write (unit_sbk,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
2451  write (unit_sbk,*) '#  '
2452  write (unit_sbk,*) '#  Columns are: '
2453  write (unit_sbk,*) '#  temperature[K]   S [au]   S [SI]     '
2454  write (unit_sbk,*) '#  '
2455 
2456  ABI_ALLOCATE(cond ,(ntemper,nsppol,3,3))
2457  ABI_ALLOCATE(cth ,(ntemper,nsppol,3,3))
2458  ABI_ALLOCATE(sbk ,(ntemper,nsppol,3,3))
2459  ABI_ALLOCATE(seebeck ,(ntemper,nsppol,3,3))
2460 
2461  cond = zero
2462  cth = zero
2463  sbk = zero
2464  seebeck = zero
2465  do isppol=1,nsppol
2466    do icomp=1, 3
2467      do jcomp=1, 3
2468        itensor=(icomp-1)*3+jcomp
2469        do itemp=1,ntemper
2470          Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp)
2471          do iene = 1, nene_all
2472            factor = (enemin+(iene-1)*deltaene*nspline - fermie(itemp))/(kb_HaK*Temp)
2473            if (factor < -40.0d0) then
2474              dfermide = zero
2475            else if (factor > 40.0d0) then
2476              dfermide = zero
2477            else
2478              dfermide = EXP(factor)/(kb_HaK*Temp*(EXP(factor)+one)**2.0d0)
2479            end if
2480            cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp) + &
2481 &           cond_e(itemp,isppol,iene,itensor)*dfermide*deltaene*nspline
2482            cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* &
2483 &           (enemin+(iene-1)*deltaene*nspline - fermie(itemp))**2.0d0*dfermide*deltaene*nspline
2484            sbk(itemp,isppol,icomp,jcomp) = sbk(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* &
2485 &           (enemin+(iene-1)*deltaene*nspline - fermie(itemp))*dfermide*deltaene*nspline
2486          end do
2487        end do ! temperature
2488      end do ! jcomp
2489    end do ! icomp
2490  end do !end isppol
2491 
2492  do isppol=1,nsppol
2493    do itemp=1,ntemper
2494      cond_inv(:,:)=cond(itemp,isppol,:,:)
2495      call matrginv(cond_inv,3,3)
2496      call DGEMM('N','N',3,3,3,one,sbk(itemp,isppol,:,:),3,cond_inv,&
2497 &     3,zero,seebeck(itemp,isppol,:,:),3)
2498    end do
2499  end do
2500 
2501  do isppol=1,nsppol
2502    do icomp=1, 3
2503      do jcomp=1, 3
2504        itensor=(icomp-1)*3+jcomp
2505        write(unit_cond,*) '# Conductivity for isppol, itrten= ',isppol,itensor
2506        write(unit_therm,*) '# Thermal conductivity for isppol, itrten= ',isppol,itensor
2507        write(unit_sbk,*) '# Seebeck coefficient for isppol, itrten= ',isppol,itensor
2508        do itemp=1,ntemper
2509          Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp)
2510 
2511          seebeck(itemp,isppol,icomp,jcomp) = -1.0d0*seebeck(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)
2512          cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp)/cryst%ucvol
2513          cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)/cryst%ucvol
2514          write(unit_cond,'(3D20.10)')Temp,cond(itemp,isppol,icomp,jcomp),cond(itemp,isppol,icomp,jcomp)*chu_cond
2515          write(unit_therm,'(3D20.10)')Temp,cth(itemp,isppol,icomp,jcomp),cth(itemp,isppol,icomp,jcomp)*chu_cth
2516          write(unit_sbk,'(3D20.10)')Temp,seebeck(itemp,isppol,icomp,jcomp),seebeck(itemp,isppol,icomp,jcomp)*chu_sbk
2517        end do ! temperature
2518        write(unit_cond,*)
2519        write(unit_therm,*)
2520        write(unit_sbk,*)
2521      end do ! jcomp
2522    end do ! icomp
2523  end do !end isppol
2524 
2525 
2526  ABI_DEALLOCATE(inv_tau_k)
2527  ABI_DEALLOCATE(tau_k)
2528  ABI_DEALLOCATE(tmp_wtk)
2529  ABI_DEALLOCATE(dos_e)
2530  ABI_DEALLOCATE(cond_e)
2531  ABI_DEALLOCATE(fermie)
2532  ABI_DEALLOCATE(cond)
2533  ABI_DEALLOCATE(sbk)
2534  ABI_DEALLOCATE(cth)
2535  ABI_DEALLOCATE(seebeck)
2536 
2537  close (unit=unit_tau)
2538  close (unit=unit_taue)
2539  close (unit=unit_mfp)
2540  close (unit=unit_invtau)
2541  close (unit=unit_cond)
2542  close (unit=unit_therm)
2543  close (unit=unit_sbk)
2544 
2545 end subroutine get_tau_k

ABINIT/m_a2ftr [ Modules ]

[ Top ] [ Modules ]

NAME

 m_a2ftr

FUNCTION

COPYRIGHT

   Copyright (C) 2004-2018 ABINIT group (JPC, MJV, 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_a2ftr
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_elphon
32  use m_errors
33  use m_abicore
34  use m_xmpi
35  use m_kptrank
36  use m_splines
37  use m_ebands
38 
39  use m_io_tools,        only : open_file
40  use m_numeric_tools,   only : simpson_int
41  use m_hide_lapack,     only : matrginv
42  use m_geometry,        only : phdispl_cart2red
43  use m_crystal,         only : crystal_t
44  use m_ifc,             only : ifc_type, ifc_fourq
45  use m_dynmat,          only : ftgam_init, ftgam
46  use m_epweights,       only : d2c_wtq, ep_ph_weights, ep_el_weights, ep_ph_weights
47  use m_fstab,           only : mkqptequiv
48 
49  implicit none
50 
51  private

ABINIT/mka2f_tr [ Functions ]

[ Top ] [ Functions ]

NAME

 mka2f_tr

FUNCTION

  calculates the FS averaged Transport alpha^2F_tr function
  calculates and outputs the associated electrical conductivity, relaxation time, and Seebeck coefficient
  and thermal conductivities
  for the first task : copied from mka2F

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%gkk2 = gkk2 matrix elements on full FS grid for each phonon mode
    elph_ds%nbranch = number of phonon branches = 3*natom
    elph_ds%nFSband = number of bands included in the FS integration
    elph_ds%k_phon%nkpt = number of kpts included in the FS integration
    elph_ds%k_phon%wtk = integration weights on the FS
    delph_ds%n0 = DOS at the Fermi level calculated from the k_phon integration weights
    elph_ds%k_phon%kpt = coordinates of all FS kpoints
  mustar = coulomb pseudopotential parameter eventually for 2 spin channels
  natom = number of atoms
  ntemper = number of temperature points to calculate, from tempermin to tempermin+ntemper*temperinc
  tempermin = minimum temperature at which resistivity etc are calculated (in K)
  temperinc = interval for temperature grid on which resistivity etc are calculated (in K)
  elph_tr_ds%dos_n0 = DOS at the Fermi level calculated from the k_phon integration
           weights, but has a temperature dependence
  elph_tr_ds%dos_n  = DOS at varied energy levels around Fermi level
  elph_tr_ds%veloc_sq0 = Fermi velocity square with T dependence

OUTPUT

  elph_ds

PARENTS

      elphon

CHILDREN

      d2c_wtq,dgemm,ep_ph_weights,ftgam,ftgam_init,gam_mult_displ,ifc_fourq
      matrginv,simpson_int,wrtout,xmpi_sum,zgemm

NOTES

   copied from ftiaf9.f

SOURCE

 109 subroutine mka2f_tr(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,pair2red,elph_tr_ds)
 110 
 111 
 112 !This section has been created automatically by the script Abilint (TD).
 113 !Do not modify the following lines by hand.
 114 #undef ABI_FUNC
 115 #define ABI_FUNC 'mka2f_tr'
 116 !End of the abilint section
 117 
 118  implicit none
 119 
 120 !Arguments ------------------------------------
 121 !scalars
 122  integer,intent(in) :: ntemper
 123  real(dp),intent(in) :: tempermin,temperinc
 124  type(ifc_type),intent(in) :: ifc
 125  type(crystal_t),intent(in) :: crystal
 126  type(elph_tr_type),intent(inout) :: elph_tr_ds
 127  type(elph_type),intent(inout) :: elph_ds
 128 !arrays
 129  integer,intent(in) :: pair2red(elph_ds%nenergy,elph_ds%nenergy)
 130 
 131 !Local variables -------------------------
 132 !x =w/(2kbT)
 133 !scalars
 134  integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr
 135  integer :: unit_a2f_tr,natom,ii,jj
 136  integer :: idir, iatom, k1, kdir
 137  integer :: unit_lor,unit_rho,unit_tau,unit_sbk, unit_therm
 138  integer :: itemp, tmp_nenergy
 139  integer :: itrtensor, icomp, jcomp!, kcomp
 140  integer :: ie, ie_1, ie2, ie_2, ie1, ie_tmp, ssp, s1(4), s2(4)
 141  integer :: ie2_left, ie2_right
 142  integer :: ik_this_proc, ierr,nrpt
 143  logical,parameter :: debug=.False.
 144  real(dp) :: Temp,chgu,chtu,chsu,chwu,diagerr,ucvol
 145  real(dp) :: a2fprefactor, gtemp
 146  real(dp) :: lambda_tr,lor0,lorentz,maxerr,omega
 147  real(dp) :: rho,tau,wtherm,xtr
 148  real(dp) :: lambda_tr_trace
 149  real(dp) :: domega, omega_min, omega_max
 150  real(dp) :: gaussval, gaussprefactor, gaussfactor, gaussmaxarg, xx
 151  real(dp) :: qnorm2, tmp_fermie
 152  real(dp) :: e1, e2, diff, xe
 153  real(dp) :: occ_omega, occ_e1, occ_e2
 154  real(dp) :: nv1, nv2, sigma1, sigma2
 155  real(dp) :: dos_n_e2, veloc_sq_icomp, veloc_sq_jcomp
 156  real(dp) :: tointegq00_1, tointegq00_2, tointegq01_1, tointegq01_2,tointegq11_1, tointegq11_2
 157  real(dp) :: j00, j01, j11
 158  real(dp) :: tointegq00,tointegq01,tointegq11
 159  real(dp) :: pref_s, pref_w, tmp_veloc_sq0, tmp_veloc_sq1, tmp_veloc_sq2
 160  character(len=500) :: message
 161  character(len=fnlen) :: fname
 162 !arrays
 163  real(dp),parameter :: c0(2)=(/0.d0,0.d0/),c1(2)=(/1.d0,0.d0/)
 164  real(dp) ::  gprimd(3,3)
 165  real(dp) :: eigval(elph_ds%nbranch)
 166  real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
 167  real(dp) :: gam_now(2,elph_ds%nbranch*elph_ds%nbranch)
 168  real(dp) :: tmpa2f(elph_ds%na2f)
 169  real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch)
 170  real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch)
 171  real(dp) :: q11_inv(3,3)
 172 !real(dp) ::  fullq(6,6)
 173  real(dp),allocatable :: phfrq(:,:)
 174  real(dp),allocatable :: tmp_a2f_1d_tr(:,:,:,:,:)
 175  real(dp),allocatable :: displ(:,:,:,:)
 176  real(dp),allocatable :: pheigvec(:,:)
 177  real(dp),allocatable :: tmp_wtq(:,:,:)
 178  real(dp),allocatable :: integrho(:), tointegrho(:)
 179  real(dp),allocatable :: integrand_q00(:),integrand_q01(:),integrand_q11(:)
 180  real(dp),allocatable :: q00(:,:,:,:), q01(:,:,:,:),q11(:,:,:,:)
 181  real(dp),allocatable :: seebeck(:,:,:,:)!, rho_nm(:,:,:,:)
 182  real(dp),allocatable :: rho_T(:)
 183  real(dp),allocatable :: coskr(:,:), sinkr(:,:)
 184 
 185 ! *********************************************************************
 186 !calculate a2f_tr for frequencies between 0 and omega_max
 187 
 188 
 189  write(std_out,*) 'mka2f_tr : enter '
 190 
 191  ucvol = crystal%ucvol
 192  natom = crystal%natom
 193  gprimd = crystal%gprimd
 194 
 195  ! number of real-space points for FT interpolation
 196  nrpt = ifc%nrpt
 197 !
 198 !MG: the step should be calculated locally using nomega and the extrema of the spectrum.
 199 !One should not rely on previous calls for the setup of elph_ds%domega
 200 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t
 201  domega =elph_ds%domega
 202  omega_min       = elph_ds%omega_min
 203  omega_max       = elph_ds%omega_max
 204 
 205  if (elph_ds%ep_lova .eq. 1) then
 206    tmp_nenergy = 1
 207  else if (elph_ds%ep_lova .eq. 0) then
 208    tmp_nenergy = elph_ds%nenergy
 209  end if
 210 
 211 !! defaults for number of temperature steps and max T (all in Kelvin...)
 212 !ntemper=1000
 213 !tempermin=zero
 214 !temperinc=one
 215 
 216  ABI_ALLOCATE(rho_T,(ntemper))
 217 
 218 
 219  gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear
 220  gaussfactor = one / elph_ds%a2fsmear
 221  gaussmaxarg = sqrt(-log(1.d-90))
 222 !lor0=(pi*kb_HaK)**2/3.
 223  lor0=pi**2/3.0_dp
 224 
 225 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F)
 226 !WARNING! supposes this value has been set in mkelph_linwid.
 227 
 228 !ENDMG
 229 
 230  maxerr=0.
 231  nerr=0
 232 
 233  ABI_ALLOCATE(tmp_wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f+1))
 234  ABI_ALLOCATE(elph_ds%k_fine%wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f))
 235  ABI_ALLOCATE(elph_ds%k_phon%wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f))
 236 
 237  ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt))
 238  ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt))
 239  ABI_ALLOCATE(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt))
 240 
 241  do iFSqpt=1,elph_ds%k_fine%nkpt
 242    call ifc_fourq(ifc,crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
 243  end do
 244  omega_min = omega_min - domega
 245 
 246 !bxu, obtain wtq for the q_fine, then condense to q_phon
 247  call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, &
 248 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_fine,tmp_wtq)
 249 !call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, &
 250 !& elph_ds%nbranch,1,elph_ds%k_fine,tmp_wtq)
 251  omega_min = omega_min + domega
 252 
 253  do iomega = 1, elph_ds%na2f
 254    elph_ds%k_fine%wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1)
 255  end do
 256  ABI_DEALLOCATE(tmp_wtq)
 257 
 258  if (elph_ds%use_k_fine == 1) then
 259    call d2c_wtq(elph_ds)
 260  end if
 261 
 262  ABI_DEALLOCATE(phfrq)
 263  ABI_DEALLOCATE(displ)
 264  ABI_DEALLOCATE(pheigvec)
 265 
 266 !reduce the dimension from fine to phon for phfrq and pheigvec
 267  ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_phon%nkpt))
 268  ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_phon%nkpt))
 269  ABI_ALLOCATE(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_phon%nkpt))
 270 
 271  do iFSqpt=1,elph_ds%k_phon%nkpt
 272    call ifc_fourq(ifc,crystal,elph_ds%k_phon%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
 273  end do
 274 
 275  ABI_ALLOCATE(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2,ntemper))
 276  ABI_ALLOCATE(tmp_a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2))
 277 
 278 ! prepare phase factors
 279  ABI_ALLOCATE(coskr, (elph_ds%k_phon%nkpt, nrpt))
 280  ABI_ALLOCATE(sinkr, (elph_ds%k_phon%nkpt, nrpt))
 281  call ftgam_init(ifc%gprim, elph_ds%k_phon%nkpt, nrpt, elph_ds%k_phon%kpt, ifc%rpt, coskr, sinkr)
 282 
 283  elph_tr_ds%a2f_1d_tr = zero
 284  tmp_a2f_1d_tr = zero
 285 
 286 
 287  do ie = 1, elph_ds%n_pair
 288    do ssp = 1,4
 289      do isppol = 1, elph_ds%nsppol
 290 
 291 !      loop over qpoint in full kpt grid (presumably dense)
 292        do ik_this_proc =1,elph_ds%k_phon%my_nkpt
 293          iFSqpt = elph_ds%k_phon%my_ikpt(ik_this_proc)
 294 
 295          qnorm2 = sum(elph_ds%k_phon%kpt(:,iFSqpt)**2)
 296 !        if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero
 297          do itrtensor=1,9
 298 
 299            if (elph_ds%ep_int_gkk == 1) then
 300              gam_now(:,:) = elph_tr_ds%gamma_qpt_tr(:,itrtensor,:,isppol,iFSqpt)
 301            else
 302 !            Do FT from real-space gamma grid to 1 qpt.
 303              call ftgam(ifc%wghatm,gam_now,elph_tr_ds%gamma_rpt_tr(:,itrtensor,:,isppol,:,ssp,ie),natom,1,nrpt,0,&
 304 &             coskr(iFSqpt,:), sinkr(iFSqpt,:))
 305            end if
 306 
 307 !          Diagonalize gamma matrix at this qpoint (complex matrix).
 308 
 309 !          if ep_scalprod==0 we have to dot in the displacement vectors here
 310            if (elph_ds%ep_scalprod==0) then
 311 
 312              displ_red(:,:,:) = zero
 313              do jbranch=1,elph_ds%nbranch
 314                do iatom=1,natom
 315                  do idir=1,3
 316                    ibranch=idir+3*(iatom-1)
 317                    do kdir=1,3
 318                      k1 = kdir+3*(iatom-1)
 319                      displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + &
 320 &                     gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt)
 321                      displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + &
 322 &                     gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt)
 323                    end do
 324                  end do
 325                end do
 326              end do
 327 
 328              tmpgam2 = reshape (gam_now, (/2,elph_ds%nbranch,elph_ds%nbranch/))
 329              call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1)
 330              do jbranch=1,elph_ds%nbranch
 331                eigval(jbranch) = tmpgam1(1, jbranch, jbranch)
 332              end do
 333 
 334            else if (elph_ds%ep_scalprod == 1) then
 335 
 336 
 337 !            NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care
 338 
 339              call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now, 3*natom,&
 340 &             pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom)
 341              call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,&
 342 &             tmpgam1, 3*natom, c0, tmpgam2, 3*natom)
 343              diagerr = zero
 344              do ibranch=1,elph_ds%nbranch
 345                eigval(ibranch) = tmpgam2(1,ibranch,ibranch)
 346                do jbranch=1,ibranch-1
 347                  diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
 348                end do
 349                do jbranch=ibranch+1,elph_ds%nbranch
 350                  diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
 351                end do
 352              end do
 353              if (diagerr > tol12) then
 354                nerr=nerr+1
 355                maxerr=max(diagerr, maxerr)
 356              end if
 357            end if  ! end ep_scalprod if
 358 
 359 !          Add all contributions from the phonon modes at this qpoint to a2f and the phonon dos.
 360            do ibranch=1,elph_ds%nbranch
 361 !            if (abs(phfrq(ibranch,iFSqpt)) < tol10) then
 362              if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. &
 363 &             (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then
 364 !              note: this should depend on the velocity of sound, to accept acoustic modes!
 365                a2fprefactor = zero
 366              else
 367 !              a2fprefactor  = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol))
 368 !              Use the dos_n0 at the lowest input temperature, assuming to be low
 369                a2fprefactor  = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt)))
 370              end if
 371 
 372              omega = omega_min
 373              tmpa2f(:) = zero
 374              do iomega=1,elph_ds%na2f
 375                xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor
 376                omega = omega+domega
 377                if (abs(xx) > gaussmaxarg) cycle
 378 
 379                gaussval = gaussprefactor*exp(-xx*xx)
 380                gtemp = gaussval*a2fprefactor
 381 
 382                if (dabs(gtemp) < 1.0d-50) gtemp = zero
 383                tmpa2f(iomega) = tmpa2f(iomega) + gtemp
 384              end do
 385 
 386 !            tmpa2f(:) = zero
 387 !            do iomega=1,elph_ds%na2f
 388 !            gtemp = a2fprefactor*elph_ds%k_phon%wtq(ibranch,iFSqpt,iomega)
 389 !            if (dabs(gtemp) < 1.0d-50) gtemp = zero
 390 !            tmpa2f(iomega) = tmpa2f(iomega) + gtemp
 391 !            end do
 392 
 393              tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) = tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) + tmpa2f(:)
 394 
 395            end do ! end ibranch
 396          end do ! end itrtensor
 397        end do ! end iFSqpt  - loop done in parallel
 398      end do ! end isppol
 399    end do ! ss'
 400  end do ! n_pair
 401 
 402  ! MG: FIXME: Why xmpi_world? besides only one CPU should perform IO (see below)
 403  ! Likely this routine is never executed in parallel
 404  call xmpi_sum (tmp_a2f_1d_tr, xmpi_world, ierr)
 405 
 406  ABI_DEALLOCATE(coskr)
 407  ABI_DEALLOCATE(sinkr)
 408 
 409  do itemp=1,ntemper  ! runs over termperature in K
 410    do isppol=1,elph_ds%nsppol
 411      do jj=1,tmp_nenergy**2
 412        do ii=1,4
 413          elph_tr_ds%a2f_1d_tr(:,:,isppol,ii,jj,itemp) = tmp_a2f_1d_tr(:,:,isppol,ii,jj)/elph_tr_ds%dos_n0(itemp,isppol)
 414        end do
 415      end do
 416    end do
 417  end do
 418 
 419  ABI_DEALLOCATE(tmp_a2f_1d_tr)
 420 
 421 !second 1 / elph_ds%k_phon%nkpt factor for the integration weights
 422  elph_tr_ds%a2f_1d_tr  = elph_tr_ds%a2f_1d_tr  / elph_ds%k_phon%nkpt
 423 
 424  if (elph_ds%ep_scalprod == 1) then
 425    write(std_out,*) 'mka2f_tr: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr
 426  end if
 427 
 428 !output the elph_tr_ds%a2f_1d_tr
 429  fname = trim(elph_ds%elph_base_name) // '_A2F_TR'
 430  if (open_file(fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then
 431    MSG_ERROR(message)
 432  end if
 433 
 434  write (unit_a2f_tr,'(a)')       '#'
 435  write (unit_a2f_tr,'(a)')       '# ABINIT package : a2f_tr file'
 436  write (unit_a2f_tr,'(a)')       '#'
 437  write (unit_a2f_tr,'(a)')       '# a2f_tr function integrated over the FS. omega in a.u.'
 438  write (unit_a2f_tr,'(a,I10)')   '#     number of kpoints integrated over : ', elph_ds%k_phon%nkpt
 439  write (unit_a2f_tr,'(a,I10)')   '#     number of energy points : ',elph_ds%na2f
 440  write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '#       between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
 441  write (unit_a2f_tr,'(a,E16.6)') '#   and the smearing width for gaussians is ', elph_ds%a2fsmear
 442  write (unit_a2f_tr,'(a)')       '#'
 443 
 444 !done with header
 445  do isppol=1,elph_ds%nsppol
 446    write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_tr_ds%dos_n0(1,isppol)
 447    omega = omega_min
 448    do iomega=1,elph_ds%na2f
 449 !    bxu, at which eps and eps' should I save it
 450 !    better to save them all, but could be too many
 451      write (unit_a2f_tr,   '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr(iomega,:,isppol,1,INT(elph_ds%n_pair/2)+1,1)
 452      omega=omega+domega
 453    end do
 454    write (unit_a2f_tr,*)
 455  end do !isppol
 456 
 457  close (unit=unit_a2f_tr)
 458 
 459 !calculation of transport properties
 460  ABI_ALLOCATE(integrho,(elph_ds%na2f))
 461  ABI_ALLOCATE(tointegrho,(elph_ds%na2f))
 462  ABI_ALLOCATE(integrand_q00,(elph_ds%na2f))
 463  ABI_ALLOCATE(integrand_q01,(elph_ds%na2f))
 464  ABI_ALLOCATE(integrand_q11,(elph_ds%na2f))
 465  ABI_ALLOCATE(q00,(ntemper,3,3,elph_ds%nsppol))
 466  ABI_ALLOCATE(q01,(ntemper,3,3,elph_ds%nsppol))
 467  ABI_ALLOCATE(q11,(ntemper,3,3,elph_ds%nsppol))
 468  ABI_ALLOCATE(seebeck,(elph_ds%nsppol,ntemper,3,3))
 469 !ABI_ALLOCATE(rho_nm,(elph_ds%nsppol,ntemper,3,3))
 470 
 471  fname = trim(elph_ds%elph_base_name) // '_RHO'
 472  if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then
 473    MSG_ERROR(message)
 474  end if
 475 !print header to resistivity file
 476  write (unit_rho,*) '# Resistivity as a function of temperature.'
 477  write (unit_rho,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
 478  write (unit_rho,*) '#  '
 479  write (unit_rho,*) '#  Columns are: '
 480  write (unit_rho,*) '#  temperature[K]   rho[au]   rho [SI]        rho/temp [au]'
 481  write (unit_rho,*) '#  '
 482 
 483  fname = trim(elph_ds%elph_base_name) // '_TAU'
 484  if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
 485    MSG_ERROR(message)
 486  end if
 487 !print header to relaxation time file
 488  write (unit_tau,*) '# Relaxation time as a function of temperature.'
 489  write (unit_tau,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
 490  write (unit_tau,*) '#  '
 491  write (unit_tau,*) '#  Columns are: '
 492  write (unit_tau,*) '#  temperature[K]   tau[au]   tau [SI]     '
 493  write (unit_tau,*) '#  '
 494 
 495  fname = trim(elph_ds%elph_base_name) // '_SBK'
 496  if (open_file(fname,message,newunit=unit_sbk,status='unknown') /= 0) then
 497    MSG_ERROR(message)
 498  end if
 499 !print header to relaxation time file
 500  write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.'
 501  write (unit_sbk,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
 502  write (unit_sbk,*) '#  '
 503  write (unit_sbk,*) '#  Columns are: '
 504  write (unit_sbk,*) '#  temperature[K]   S [au]   S [SI]     '
 505  write (unit_sbk,*) '#  '
 506 
 507  fname = trim(elph_ds%elph_base_name) // '_WTH'
 508  if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
 509    MSG_ERROR(message)
 510  end if
 511 
 512 !print header to thermal conductivity file
 513  write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.'
 514  write (unit_therm,'(a)') '#  the formalism is isotropic, so non-cubic crystals may be wrong'
 515  write (unit_therm,'(a)') '#  '
 516  write (unit_therm,'(a)') '#  Columns are: '
 517  write (unit_therm,'(a)') '#  temperature[K]   thermal rho[au]   thermal cond [au]   thermal rho [SI]   thermal cond [SI]'
 518  write (unit_therm,'(a)') '#  '
 519 
 520  fname = trim(elph_ds%elph_base_name) // '_LOR'
 521  if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then
 522    MSG_ERROR(message)
 523  end if
 524 
 525 !print header to lorentz file
 526  write (unit_lor,*) '# Lorentz number as a function of temperature.'
 527  write (unit_lor,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
 528  write (unit_lor,*) '#  '
 529  write (unit_lor,*) '#  Columns are: '
 530  write (unit_lor,*) '#  temperature[K]   Lorentz number[au]   Lorentz quantum = (pi*kb_HaK)**2/3'
 531  write (unit_lor,*) '#  '
 532 
 533  do isppol=1,elph_ds%nsppol
 534    lambda_tr_trace = zero
 535    do itrtensor=1,9
 536      omega = omega_min
 537      tointegrho = zero
 538      do iomega=1,elph_ds%na2f
 539        if(omega<=0) then
 540          omega=omega+domega
 541          cycle
 542        end if
 543 !      bxu, agian, which eps and eps' to use?
 544 !      sometimes Ef is in the gap
 545        tointegrho(iomega)=two*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega
 546        omega=omega+domega
 547      end do
 548 
 549      integrho = zero
 550      call simpson_int(elph_ds%na2f,domega,tointegrho,integrho)
 551      lambda_tr=integrho(elph_ds%na2f)
 552      write (message, '(a,2i3,a,es16.6)' )&
 553 &     ' mka2f_tr: TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' =  ', lambda_tr
 554      call wrtout(std_out,message,'COLL')
 555      if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr
 556    end do !end itrtensor do
 557 
 558    lambda_tr_trace = lambda_tr_trace / three
 559    write (message, '(a,i3,a,es16.6)' )&
 560 &   ' mka2f_tr: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' =  ', lambda_tr_trace
 561    call wrtout(std_out,message,'COLL')
 562    call wrtout(ab_out,message,'COLL')
 563  end do !end isppol do
 564 
 565 !constant to change units of rho from au to SI
 566  chgu=2.173969*1.0d-7
 567  chtu=2.4188843265*1.0d-17
 568  chsu=8.617343101*1.0d-5 ! au to J/C/K
 569  chwu=9.270955772*1.0d-5 ! au to mK/W
 570 
 571 !change the fermi level to zero, as required for q01 to vanish.
 572  tmp_fermie = elph_ds%fermie
 573 !Get Q00, Q01, Q11, and derive rho, tau
 574  q00 = zero
 575  q01 = zero
 576  q11 = zero
 577 ! prepare s1 and s2 arrays
 578  s1 = (/1, 1, -1, -1/)
 579  s2 = (/1, -1, 1, -1/)
 580 
 581  do isppol=1,elph_ds%nsppol
 582    do icomp=1, 3
 583      do jcomp=1, 3
 584        itrtensor=(icomp-1)*3+jcomp
 585 
 586        write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor
 587        write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor
 588 
 589        do itemp=1,ntemper  ! runs over termperature in K
 590          Temp=tempermin+temperinc*dble(itemp)
 591          tmp_veloc_sq0 = sqrt(elph_tr_ds%veloc_sq0(itemp,icomp,isppol)*elph_tr_ds%veloc_sq0(itemp,jcomp,isppol))
 592 
 593          integrand_q00 = zero
 594          integrand_q01 = zero
 595          integrand_q11 = zero
 596 
 597          omega = omega_min
 598          do iomega=1,elph_ds%na2f
 599            if(omega .le. 0) then
 600              omega=omega+domega
 601              cycle
 602            end if
 603            xtr=omega/(kb_HaK*Temp)
 604            occ_omega=1.0_dp/(exp(xtr)-1.0_dp)
 605 
 606            tmp_veloc_sq1 = zero
 607            tmp_veloc_sq2 = zero
 608            do ie1=1,elph_ds%nenergy
 609              e1 = elph_tr_ds%en_all(isppol,ie1)
 610 
 611 !! BXU, the tolerance here needs to be used with caution
 612 !! which depends on the dimensions of the system
 613 !! E.g. in 2D system, DOS can be much smaller
 614              if (elph_tr_ds%dos_n(ie1,isppol)/natom .lt. 0.1d0) cycle ! energy in the gap forbidden
 615 
 616              xtr=(e1-tmp_fermie)/(kb_HaK*Temp)
 617              occ_e1=1.0_dp/(exp(xtr)+1.0_dp)
 618 
 619              e2 = e1 + omega
 620              xtr=(e2-tmp_fermie)/(kb_HaK*Temp)
 621              occ_e2=1.0_dp/(exp(xtr)+1.0_dp)
 622 !            Do we need to change the fermie to the one with T dependence?
 623 !            find which ie2 give the closest energy
 624              if (e2 .gt. elph_tr_ds%en_all(isppol,elph_ds%nenergy)) then
 625                ie2 = 0
 626                cycle
 627              else
 628                ie_tmp = 1
 629                diff = dabs(e2-elph_tr_ds%en_all(isppol,1))
 630                do ie2 = 2, elph_ds%nenergy
 631                  if (dabs(e2-elph_tr_ds%en_all(isppol,ie2)) .lt. diff) then
 632                    diff = dabs(e2-elph_tr_ds%en_all(isppol,ie2))
 633                    ie_tmp = ie2
 634                  end if
 635                end do
 636                ie2 = ie_tmp
 637 
 638                if (e2 < elph_tr_ds%en_all(isppol,ie2)) then
 639                  ie2_right = ie2
 640                  ie2_left  = ie2-1
 641                else
 642                  ie2_right = ie2+1
 643                  ie2_left  = ie2
 644                end if
 645 
 646              end if
 647 
 648              if (elph_tr_ds%dos_n(ie2,isppol)/natom .lt. 0.1d0) cycle
 649 
 650              tointegq00 = zero
 651              tointegq01 = zero
 652              tointegq11 = zero
 653 
 654 ! BXU linear interpolation of volec_sq and dos_n
 655              xe=(e2-elph_tr_ds%en_all(isppol,ie2_left))/ &
 656 &             (elph_tr_ds%en_all(isppol,ie2_right)-elph_tr_ds%en_all(isppol,ie2_left))
 657              veloc_sq_icomp = elph_tr_ds%veloc_sq(icomp,isppol,ie2_left)*(1.0d0-xe) + &
 658 &             elph_tr_ds%veloc_sq(icomp,isppol,ie2_right)*xe
 659              veloc_sq_jcomp = elph_tr_ds%veloc_sq(jcomp,isppol,ie2_left)*(1.0d0-xe) + &
 660 &             elph_tr_ds%veloc_sq(jcomp,isppol,ie2_right)*xe
 661              dos_n_e2 = elph_tr_ds%dos_n(ie2_left,isppol)*(1.0d0-xe) + &
 662 &             elph_tr_ds%dos_n(ie2_right,isppol)*xe
 663 
 664              tmp_veloc_sq1 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie1)*elph_tr_ds%veloc_sq(jcomp,isppol,ie1))
 665 !             tmp_veloc_sq2 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie2)*elph_tr_ds%veloc_sq(jcomp,isppol,ie2))
 666              tmp_veloc_sq2 = sqrt(veloc_sq_icomp*veloc_sq_jcomp)
 667 
 668 !            ie_1 = (ie1-1)*elph_ds%nenergy + ie2
 669 !            ie_2 = (ie2-1)*elph_ds%nenergy + ie1
 670              ie_1 = pair2red(ie1,ie2)
 671              ie_2 = pair2red(ie2,ie1)
 672              if (ie_1 .eq. 0 .or. ie_2 .eq. 0) then
 673                MSG_BUG('CHECK pair2red!')
 674              end if
 675 
 676              do ssp=1,4  ! (s,s'=+/-1, condense the indices)
 677 
 678                nv1 = 1.0_dp/(elph_tr_ds%dos_n(ie1,isppol)*sqrt(tmp_veloc_sq1))
 679                sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK)
 680 !DEBUG
 681                if (elph_ds%ep_lova .eq. 1) then
 682                  nv1 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0))
 683                  sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK)
 684                end if
 685 !ENDDEBUG
 686 
 687                tointegq00_1 = zero
 688                tointegq01_1 = zero
 689                tointegq11_1 = zero
 690 
 691 !DEBUG
 692                if (elph_ds%ep_lova .eq. 1) then
 693                  nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0))
 694                  sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
 695                  j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp
 696                  j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
 697                  j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
 698                  tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
 699 &                 occ_e1*(1.0_dp-occ_e2)*j00*occ_omega
 700                  tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
 701 &                 occ_e1*(1.0_dp-occ_e2)*j01*occ_omega
 702                  tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
 703 &                 occ_e1*(1.0_dp-occ_e2)*j11*occ_omega
 704 !END DEBUG
 705                else if (elph_ds%ep_lova .eq. 0) then
 706                  nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2))
 707                  sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
 708                  j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp
 709                  j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
 710                  j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp
 711 !                bxu TEST
 712                  if (debug) then
 713                    if (ssp .eq. 1 .and. itrtensor .eq. 1) then
 714                      write(21,"(3i5,4E20.12)") iomega, ie1, ie2, &
 715 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
 716 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
 717 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
 718                    end if
 719                    if (ssp .eq. 2 .and. itrtensor .eq. 1) then
 720                      write(22,"(3i5,4E20.12)") iomega, ie1, ie2, &
 721 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
 722 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
 723 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
 724                    end if
 725                    if (ssp .eq. 3 .and. itrtensor .eq. 1) then
 726                      write(23,"(3i5,4E20.12)") iomega, ie1, ie2, &
 727 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
 728 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
 729 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
 730                    end if
 731                    if (ssp .eq. 4 .and. itrtensor .eq. 1) then
 732                      write(24,"(3i5,4E20.12)") iomega, ie1, ie2, &
 733 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, &
 734 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, &
 735 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega
 736                    end if
 737                  end if
 738                  tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* &
 739 &                 occ_e1*(1.0_dp-occ_e2)*j00*occ_omega
 740                  tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* &
 741 &                 occ_e1*(1.0_dp-occ_e2)*j01*occ_omega
 742                  tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* &
 743 &                 occ_e1*(1.0_dp-occ_e2)*j11*occ_omega
 744                end if
 745 
 746                tointegq00_2 = zero
 747                tointegq01_2 = zero
 748                tointegq11_2 = zero
 749 
 750 !DEBUG
 751                if (elph_ds%ep_lova .eq. 1) then
 752                  nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0))
 753                  sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
 754                  j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp
 755                  j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
 756                  j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
 757                  tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
 758 &                 occ_e1*(1.0_dp-occ_e2)*j00*(occ_omega+1)
 759                  tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
 760 &                 occ_e1*(1.0_dp-occ_e2)*j01*(occ_omega+1)
 761                  tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* &
 762 &                 occ_e1*(1.0_dp-occ_e2)*j11*(occ_omega+1)
 763 !END DEBUG
 764                else if (elph_ds%ep_lova .eq. 0) then
 765                  nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2))
 766                  sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK)
 767                  j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp
 768                  j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
 769                  j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp
 770 !DEBUG           bxu TEST
 771                  if (debug) then
 772                    if (ssp .eq. 1 .and. itrtensor .eq. 1) then
 773                      write(31,"(3i5,4E20.12)") iomega, ie2, ie1, &
 774 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
 775 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
 776 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
 777                    end if
 778                    if (ssp .eq. 2 .and. itrtensor .eq. 1) then
 779                      write(32,"(3i5,4E20.12)") iomega, ie2, ie1, &
 780 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
 781 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
 782 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
 783                    end if
 784                    if (ssp .eq. 3 .and. itrtensor .eq. 1) then
 785                      write(33,"(3i5,4E20.12)") iomega, ie2, ie1, &
 786 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
 787 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
 788 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
 789                    end if
 790                    if (ssp .eq. 4 .and. itrtensor .eq. 1) then
 791                      write(34,"(3i5,4E20.12)") iomega, ie2, ie1, &
 792 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, &
 793 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, &
 794 &                     elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1)
 795                    end if
 796                  end if
 797 !ENDDEBUG
 798                  tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* &
 799 &                 occ_e2*(1.0_dp-occ_e1)*j00*(occ_omega+1)
 800                  tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* &
 801 &                 occ_e2*(1.0_dp-occ_e1)*j01*(occ_omega+1)
 802                  tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* &
 803 &                 occ_e2*(1.0_dp-occ_e1)*j11*(occ_omega+1)
 804                end if ! elph_ds%ep_lova
 805 
 806                tointegq00 = tointegq00 + tointegq00_1 + tointegq00_2
 807                tointegq01 = tointegq01 + tointegq01_1 + tointegq01_2
 808                tointegq11 = tointegq11 + tointegq11_1 + tointegq11_2
 809 
 810              end do ! ss' = 4
 811              integrand_q00(iomega) = integrand_q00(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq00
 812              integrand_q01(iomega) = integrand_q01(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq01
 813              integrand_q11(iomega) = integrand_q11(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq11
 814            end do ! ie1 ~ 20
 815            omega=omega+domega
 816            q00(itemp,icomp,jcomp,isppol) = q00(itemp,icomp,jcomp,isppol) +&
 817 &           domega*integrand_q00(iomega)
 818            q01(itemp,icomp,jcomp,isppol) = q01(itemp,icomp,jcomp,isppol) +&
 819 &           domega*integrand_q01(iomega)
 820            q11(itemp,icomp,jcomp,isppol) = q11(itemp,icomp,jcomp,isppol) +&
 821 &           domega*integrand_q11(iomega)
 822          end do ! omega ~ 400
 823 
 824          q00(itemp,icomp,jcomp,isppol)=q00(itemp,icomp,jcomp,isppol)* &
 825 &         ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp)
 826          q01(itemp,icomp,jcomp,isppol)=q01(itemp,icomp,jcomp,isppol)* &
 827 &         ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp)
 828          q11(itemp,icomp,jcomp,isppol)=q11(itemp,icomp,jcomp,isppol)* &
 829 &         ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp)
 830 
 831          rho = 0.5_dp*q00(itemp,icomp,jcomp,isppol)
 832 !        is tau energy dependent?
 833          tau = 2.0d0*ucvol/(q00(itemp,icomp,jcomp,isppol)*elph_tr_ds%dos_n0(itemp,isppol)*tmp_veloc_sq0)
 834          write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp
 835          write(unit_tau,'(3D20.10)')temp,tau,tau*chtu
 836          rho_T(itemp)=rho
 837        end do ! temperature = 1?
 838        write(unit_rho,*)
 839        write(unit_tau,*)
 840 
 841      end do ! jcomp = 3
 842    end do ! icomp = 3
 843  end do ! isppol = 2
 844 
 845 !-----------------------------
 846 
 847  seebeck = zero
 848 !rho_nm  = zero
 849 
 850 !do isppol=1,elph_ds%nsppol
 851 !do itemp=1,ntemper
 852 !q11_inv(:,:)=q11(itemp,:,:,isppol)
 853 !call matrginv(q11_inv,3,3)
 854 !do icomp=1,3
 855 !do jcomp=1,3
 856 !do kcomp=1,3
 857 !seebeck(isppol,itemp,icomp,jcomp) = seebeck(isppol,itemp,icomp,jcomp) + &
 858 !&                             q01(itemp,icomp,kcomp,isppol)*q11_inv(kcomp,jcomp)
 859 !end do
 860 !end do
 861 !end do
 862 !end do
 863 !end do
 864 
 865  do isppol=1,elph_ds%nsppol
 866    do itemp=1,ntemper
 867      q11_inv(:,:)=q11(itemp,:,:,isppol)
 868      call matrginv(q11_inv,3,3)
 869      call DGEMM('N','N',3,3,3,one,q01(itemp,:,:,isppol),3,q11_inv,&
 870 &     3,zero,seebeck(isppol,itemp,:,:),3)
 871 !    call DGEMM('N','N',3,3,3,one,seebeck(isppol,itemp,:,:),3,q01(itemp,:,:,isppol),&
 872 !    &     3,zero,rho_nm(isppol,itemp,:,:),3)
 873    end do
 874  end do
 875  pref_s = pi/sqrt(3.0_dp)
 876  seebeck=pref_s*seebeck
 877 
 878 !fullq = zero
 879 !do icomp=1,3
 880 !do jcomp=1,3
 881 !fullq(icomp,jcomp) = q00(1,icomp,jcomp,1)
 882 !end do
 883 !end do
 884 !do icomp=1,3
 885 !do jcomp=4,6
 886 !fullq(icomp,jcomp) = q01(1,icomp,jcomp-3,1)
 887 !end do
 888 !end do
 889 !do icomp=4,6
 890 !do jcomp=1,3
 891 !fullq(icomp,jcomp) = q01(1,icomp-3,jcomp,1)
 892 !end do
 893 !end do
 894 !do icomp=4,6
 895 !do jcomp=4,6
 896 !fullq(icomp,jcomp) = q11(1,icomp-3,jcomp-3,1)
 897 !end do
 898 !end do
 899 !write(*,*)' fullq'
 900 !write(*,"(6E20.12)") (fullq(1,jcomp),jcomp=1,6)
 901 !write(*,"(6E20.12)") (fullq(2,jcomp),jcomp=1,6)
 902 !write(*,"(6E20.12)") (fullq(3,jcomp),jcomp=1,6)
 903 !write(*,"(6E20.12)") (fullq(4,jcomp),jcomp=1,6)
 904 !write(*,"(6E20.12)") (fullq(5,jcomp),jcomp=1,6)
 905 !write(*,"(6E20.12)") (fullq(6,jcomp),jcomp=1,6)
 906  write(message,'(a)') 'q00:'
 907  call wrtout(std_out,message,'COLL')
 908  write(message,'(3E20.12)') (q00(1,1,jcomp,1),jcomp=1,3)
 909  call wrtout(std_out,message,'COLL')
 910  write(message,'(3E20.12)') (q00(1,2,jcomp,1),jcomp=1,3)
 911  call wrtout(std_out,message,'COLL')
 912  write(message,'(3E20.12)') (q00(1,3,jcomp,1),jcomp=1,3)
 913  call wrtout(std_out,message,'COLL')
 914  write(message,'(a)') 'q01:'
 915  call wrtout(std_out,message,'COLL')
 916  write(message,'(3E20.12)') (q01(1,1,jcomp,1),jcomp=1,3)
 917  call wrtout(std_out,message,'COLL')
 918  write(message,'(3E20.12)') (q01(1,2,jcomp,1),jcomp=1,3)
 919  call wrtout(std_out,message,'COLL')
 920  write(message,'(3E20.12)') (q01(1,3,jcomp,1),jcomp=1,3)
 921  call wrtout(std_out,message,'COLL')
 922  write(message,'(a)') 'q11, q11_inv:'
 923  call wrtout(std_out,message,'COLL')
 924  write(message,'(6E20.12)') (q11(1,1,jcomp,1),jcomp=1,3),(q11_inv(1,jcomp),jcomp=1,3)
 925  call wrtout(std_out,message,'COLL')
 926  write(message,'(6E20.12)') (q11(1,2,jcomp,1),jcomp=1,3),(q11_inv(2,jcomp),jcomp=1,3)
 927  call wrtout(std_out,message,'COLL')
 928  write(message,'(6E20.12)') (q11(1,3,jcomp,1),jcomp=1,3),(q11_inv(3,jcomp),jcomp=1,3)
 929  call wrtout(std_out,message,'COLL')
 930 !q11_inv = zero
 931 !do icomp = 1, 3
 932 !q11_inv(icomp,icomp) = 2.0_dp
 933 !end do
 934 
 935 !call matrginv(fullq,6,6)
 936 
 937 !do isppol=1,elph_ds%nsppol
 938 !do itemp=1,ntemper
 939 !rho_nm(isppol,itemp,:,:) = q00(itemp,:,:,isppol) - rho_nm(isppol,itemp,:,:)
 940 !end do
 941 !end do
 942 !rho_nm = 0.5_dp*rho_nm
 943 
 944 !Output of Seebeck coefficient
 945  do isppol=1,elph_ds%nsppol
 946    do icomp=1,3
 947      do jcomp=1,3
 948        itrtensor=(icomp-1)*3+jcomp
 949        write(unit_sbk,*) '# Seebeck for isppol, itrten = ', isppol, itrtensor
 950 !      write(88,*) '# Rho for isppol, itrten = ', isppol, itrtensor
 951 !      write(89,*) '# Rho for isppol, itrten = ', isppol, itrtensor
 952        do itemp=1,ntemper
 953          Temp=tempermin+temperinc*dble(itemp)
 954          write(unit_sbk,'(3D20.10)')temp, seebeck(isppol,itemp,icomp,jcomp), seebeck(isppol,itemp,icomp,jcomp)*chsu
 955 !        write(88,'(3D20.10)')temp, rho_nm(isppol,itemp,icomp,jcomp), rho_nm(isppol,itemp,icomp,jcomp)*chgu
 956 !        write(89,'(3D20.10)')temp, 0.5_dp/fullq(1,1), 0.5_dp*chgu/fullq(1,1)
 957        end do ! temperature
 958        write(unit_sbk,*)
 959 !      write(88,*)
 960 !      write(89,*)
 961      end do ! jcomp
 962    end do ! icomp
 963  end do ! isppol
 964 
 965 !Get thermal resistivity, based on eqn. (52) in Allen's PRB 17, 3725 (1978) [[cite:Allen1978]]
 966 !WARNING: before 6.13.1 the thermal resistivity and Lorentz number were not in
 967 !atomic units, BUT the SI units are good.
 968  pref_w = 3.0_dp/(2.0_dp*pi**2.0d0)
 969  do isppol=1,elph_ds%nsppol
 970    do icomp=1, 3
 971      do jcomp=1, 3
 972        itrtensor=(icomp-1)*3+jcomp
 973 
 974        write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol
 975        write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol
 976 
 977        do itemp=1,ntemper
 978 
 979          Temp=tempermin + temperinc*dble(itemp)
 980 
 981          wtherm = pref_w*q11(itemp,icomp,jcomp,isppol)/(kb_HaK*Temp)
 982 
 983 !        write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9
 984          write(unit_therm,'(5D20.10)')temp,wtherm,1.0_dp/wtherm,wtherm*chwu,1.0_dp/(wtherm*chwu)
 985 
 986          lorentz=rho_T(itemp)/(wtherm*kb_HaK*Temp)
 987          write(unit_lor,*)temp,lorentz,lor0
 988 
 989        end do
 990        write(unit_therm,*)
 991        write(unit_lor,*)
 992      end do ! jcomp
 993    end do ! icomp
 994  end do ! isppol
 995 
 996 
 997  ABI_DEALLOCATE(phfrq)
 998  ABI_DEALLOCATE(displ)
 999  ABI_DEALLOCATE(pheigvec)
1000  ABI_DEALLOCATE(integrand_q00)
1001  ABI_DEALLOCATE(integrand_q01)
1002  ABI_DEALLOCATE(integrand_q11)
1003  ABI_DEALLOCATE(q00)
1004  ABI_DEALLOCATE(q01)
1005  ABI_DEALLOCATE(q11)
1006  ABI_DEALLOCATE(seebeck)
1007  ABI_DEALLOCATE(rho_T)
1008  ABI_DEALLOCATE(integrho)
1009  ABI_DEALLOCATE(tointegrho)
1010 
1011  close (unit=unit_lor)
1012  close (unit=unit_rho)
1013  close (unit=unit_tau)
1014  close (unit=unit_sbk)
1015  close (unit=unit_therm)
1016 
1017  ABI_DEALLOCATE(elph_ds%k_fine%wtq)
1018  ABI_DEALLOCATE(elph_ds%k_phon%wtq)
1019 
1020  ABI_DEALLOCATE(elph_tr_ds%a2f_1d_tr)
1021 
1022  ABI_DEALLOCATE(elph_tr_ds%gamma_qpt_tr)
1023  ABI_DEALLOCATE(elph_tr_ds%gamma_rpt_tr)
1024  write(std_out,*) ' mka2f_tr : end '
1025 
1026 end subroutine mka2f_tr

ABINIT/mka2f_tr_lova [ Functions ]

[ Top ] [ Functions ]

NAME

 mka2f_tr_lova

FUNCTION

  calculates the FS averaged Transport alpha^2F_tr alpha^2F_trout alpha^2F_trin functions
  calculates and outputs the associated electrical and thermal conductivities
  for the first task: copied from mka2F

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%gkk2 = gkk2 matrix elements on full FS grid for each phonon mode
    elph_ds%nbranch = number of phonon branches = 3*natom
    elph_ds%nFSband = number of bands included in the FS integration
    elph_ds%k_fine%nkpt = number of kpts included in the FS integration
    elph_ds%k_fine%wtk = integration weights on the FS
    delph_ds%n0 = DOS at the Fermi level calculated from the k_fine integration weights
    elph_ds%k_fine%kpt = coordinates of all FS kpoints
  mustar = coulomb pseudopotential parameter
       eventually for 2 spin channels
  ntemper = number of temperature points to calculate, from tempermin to tempermin+ntemper*temperinc
  tempermin = minimum temperature at which resistivity etc are calculated (in K)
  temperinc = interval for temperature grid on which resistivity etc are calculated (in K)

OUTPUT

  elph_ds

PARENTS

      elphon

CHILDREN

      ftgam,ftgam_init,gam_mult_displ,ifc_fourq,simpson_int,wrtout,zgemm

NOTES

   copied from ftiaf9.f

SOURCE

1071 subroutine mka2f_tr_lova(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,elph_tr_ds)
1072 
1073 
1074 !This section has been created automatically by the script Abilint (TD).
1075 !Do not modify the following lines by hand.
1076 #undef ABI_FUNC
1077 #define ABI_FUNC 'mka2f_tr_lova'
1078 !End of the abilint section
1079 
1080  implicit none
1081 
1082 !Arguments ------------------------------------
1083 !scalars
1084  integer,intent(in) :: ntemper
1085  real(dp),intent(in) :: tempermin,temperinc
1086  type(crystal_t),intent(in) :: crystal
1087  type(ifc_type),intent(in) :: ifc
1088  type(elph_tr_type),intent(inout) :: elph_tr_ds
1089  type(elph_type),intent(inout) :: elph_ds
1090 
1091 !Local variables -------------------------
1092 !x =w/(2kbT)
1093 !scalars
1094  integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr
1095  integer :: unit_a2f_tr, unit_a2f_trout, unit_a2f_trin, natom
1096  integer :: idir, iatom, k1, kdir,unit_lor,unit_rho,unit_tau,unit_therm
1097  integer :: itemp,nrpt,itrtensor, icomp, jcomp
1098  real(dp) :: Temp,chgu,chtu,femto,diagerr,firh,firhT,gaussfactor,domega
1099  real(dp) :: firh_tau,firhT_tau ! added by BX to get Tau
1100  real(dp) :: a2fprefactor_in, temp_in
1101  real(dp) :: a2fprefactor_out, temp_out
1102  real(dp) :: gaussprefactor,gaussval,lambda_tr,lor0,lorentz,maxerr,maxx,omega
1103  real(dp) :: rho,tau,tolexp,wtherm,xtr,xx
1104  real(dp) :: lambda_tr_trace,omega_min, omega_max,qnorm2,spinfact
1105  character(len=500) :: message
1106  character(len=fnlen) :: fname
1107 !arrays
1108  real(dp),parameter :: c0(2)=(/0.d0,0.d0/),c1(2)=(/1.d0,0.d0/)
1109  real(dp) :: gprimd(3,3)
1110  real(dp) :: eigval_in(elph_ds%nbranch)
1111  real(dp) :: eigval_out(elph_ds%nbranch)
1112  real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1113  real(dp) :: gam_now_in (2,elph_ds%nbranch*elph_ds%nbranch)
1114  real(dp) :: gam_now_out(2,elph_ds%nbranch*elph_ds%nbranch)
1115  real(dp) :: tmpa2f_in (elph_ds%na2f)
1116  real(dp) :: tmpa2f_out(elph_ds%na2f)
1117  real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch)
1118  real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch)
1119  real(dp),allocatable :: phfrq(:,:)
1120  real(dp),allocatable :: displ(:,:,:,:)
1121  real(dp),allocatable :: pheigvec(:,:)
1122  real(dp),allocatable :: integrho(:),integtau(:),tointegrho(:),tointega2f(:),tointegtau(:)
1123  real(dp),allocatable :: rho_T(:),tau_T(:)
1124  real(dp),allocatable :: coskr(:,:)
1125  real(dp),allocatable :: sinkr(:,:)
1126 
1127 ! *********************************************************************
1128 
1129 !calculate a2f_tr for frequencies between 0 and omega_max
1130  write(std_out,*) 'mka2f_tr_lova : enter '
1131 !
1132 !MG: the step should be calculated locally using nomega and the extrema of the spectrum.
1133 !One should not rely on previous calls for the setup of elph_ds%domega
1134 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t
1135  domega =elph_ds%domega
1136 
1137  ! Number of points for FFT interpolation
1138  nrpt = ifc%nrpt
1139  natom = crystal%natom
1140  gprimd = crystal%gprimd
1141 
1142  ABI_ALLOCATE(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,1,1,1))
1143  ABI_ALLOCATE(elph_tr_ds%a2f_1d_trin,(elph_ds%na2f,9,elph_ds%nsppol))
1144  ABI_ALLOCATE(elph_tr_ds%a2f_1d_trout,(elph_ds%na2f,9,elph_ds%nsppol))
1145 
1146 !! defaults for number of temperature steps and max T (all in Kelvin...)
1147 !ntemper=1000
1148 !tempermin=zero
1149 !temperinc=one
1150  ABI_ALLOCATE(rho_T,(ntemper))
1151  ABI_ALLOCATE(tau_T,(ntemper))
1152 
1153 
1154 !tolerance on gaussian being = 0
1155  tolexp = 1.d-100
1156  maxx = sqrt(-log(tolexp))
1157  lor0=(pi*kb_HaK)**2/3.
1158 
1159 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F)
1160 !WARNING! supposes this value has been set in mkelph_linwid.
1161 
1162  gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear
1163  gaussfactor = one / elph_ds%a2fsmear
1164 
1165 !spinfact should be 1 for a normal non sppol calculation without spinorbit
1166 !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1167 !for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1168  spinfact = one / elph_ds%nsppol !/ elph_ds%nspinor
1169 
1170 !ENDMG
1171 
1172  elph_tr_ds%a2f_1d_tr = zero
1173  elph_tr_ds%a2f_1d_trin = zero
1174  elph_tr_ds%a2f_1d_trout = zero
1175 
1176  maxerr=0.
1177  nerr=0
1178 
1179  ABI_ALLOCATE(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt))
1180  ABI_ALLOCATE(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt))
1181  ABI_ALLOCATE(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt))
1182 
1183  do iFSqpt=1,elph_ds%k_fine%nkpt
1184    call ifc_fourq(ifc,crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt))
1185  end do
1186 
1187  omega_min = minval(phfrq)
1188  omega_max = maxval(phfrq)
1189 
1190  ABI_ALLOCATE(coskr, (elph_ds%k_fine%nkpt,nrpt))
1191  ABI_ALLOCATE(sinkr, (elph_ds%k_fine%nkpt,nrpt))
1192  call ftgam_init(Ifc%gprim, elph_ds%k_fine%nkpt, nrpt, elph_ds%k_fine%kpt, Ifc%rpt, coskr, sinkr)
1193 
1194  do isppol=1,elph_ds%nsppol
1195 
1196 !  loop over qpoint in full kpt grid (presumably dense)
1197    do iFSqpt=1,elph_ds%k_fine%nkpt
1198      qnorm2 = sum(elph_ds%k_fine%kpt(:,iFSqpt)**2)
1199 !    if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero
1200      do itrtensor=1,9
1201 !      Do FT from real-space gamma grid to 1 qpt.
1202 
1203        if (elph_ds%ep_int_gkk == 1) then
1204          gam_now_in(:,:) = elph_tr_ds%gamma_qpt_trin(:,itrtensor,:,isppol,iFSqpt)
1205          gam_now_out(:,:) = elph_tr_ds%gamma_qpt_trout(:,itrtensor,:,isppol,iFSqpt)
1206        else
1207          call ftgam(Ifc%wghatm,gam_now_in, elph_tr_ds%gamma_rpt_trin(:,itrtensor,:,isppol,:),natom,1,nrpt,0,&
1208 &         coskr(iFSqpt,:), sinkr(iFSqpt,:))
1209          call ftgam(Ifc%wghatm,gam_now_out,elph_tr_ds%gamma_rpt_trout(:,itrtensor,:,isppol,:),natom,1,nrpt,0,&
1210 &         coskr(iFSqpt,:), sinkr(iFSqpt,:))
1211        end if
1212 
1213 !      Diagonalize gamma matrix at this qpoint (complex matrix).
1214 
1215 !      if ep_scalprod==0 we have to dot in the displacement vectors here
1216        if (elph_ds%ep_scalprod==0) then
1217 
1218          displ_red(:,:,:) = zero
1219          do jbranch=1,elph_ds%nbranch
1220            do iatom=1,natom
1221              do idir=1,3
1222                ibranch=idir+3*(iatom-1)
1223                do kdir=1,3
1224                  k1 = kdir+3*(iatom-1)
1225                  displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + &
1226 &                 gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt)
1227                  displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + &
1228 &                 gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt)
1229                end do
1230              end do
1231            end do
1232          end do
1233 
1234          tmpgam2 = reshape (gam_now_in, (/2,elph_ds%nbranch,elph_ds%nbranch/))
1235          call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1)
1236          do jbranch=1,elph_ds%nbranch
1237            eigval_in(jbranch)   = tmpgam1(1, jbranch, jbranch)
1238          end do
1239 
1240          tmpgam2 = reshape (gam_now_out, (/2,elph_ds%nbranch,elph_ds%nbranch/))
1241          call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1)
1242          do jbranch=1,elph_ds%nbranch
1243            eigval_out(jbranch)   = tmpgam1(1, jbranch, jbranch)
1244          end do
1245 
1246        else if (elph_ds%ep_scalprod == 1) then
1247 
1248 !
1249 !        NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care
1250          call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_in, 3*natom,&
1251 &         pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom)
1252          call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,&
1253 &         tmpgam1, 3*natom, c0, tmpgam2, 3*natom)
1254          diagerr = zero
1255 
1256          do ibranch=1,elph_ds%nbranch
1257            eigval_in(ibranch) = tmpgam2(1,ibranch,ibranch)
1258            do jbranch=1,ibranch-1
1259              diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1260            end do
1261            do jbranch=ibranch+1,elph_ds%nbranch
1262              diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1263            end do
1264          end do
1265          if (diagerr > tol12) then
1266            nerr=nerr+1
1267            maxerr=max(diagerr, maxerr)
1268          end if
1269 
1270          call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_out, 3*natom,&
1271 &         pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom)
1272          call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,&
1273 &         tmpgam1, 3*natom, c0, tmpgam2, 3*natom)
1274          diagerr = zero
1275 
1276          do ibranch=1,elph_ds%nbranch
1277            eigval_out(ibranch) = tmpgam2(1,ibranch,ibranch)
1278            do jbranch=1,ibranch-1
1279              diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1280            end do
1281            do jbranch=ibranch+1,elph_ds%nbranch
1282              diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))
1283            end do
1284          end do
1285          if (diagerr > tol12) then
1286            nerr=nerr+1
1287            maxerr=max(diagerr, maxerr)
1288          end if
1289        end if
1290 !      end ep_scalprod if
1291 
1292 !      Add all contributions from the phonon modes at this qpoint to
1293 !      a2f and the phonon dos.
1294        do ibranch=1,elph_ds%nbranch
1295 !        if (abs(phfrq(ibranch,iFSqpt)) < tol10) then
1296          if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. &
1297 &         (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then !
1298 !          note: this should depend on the velocity of sound, to accept acoustic
1299 !          modes!
1300            a2fprefactor_in = zero
1301            a2fprefactor_out= zero
1302          else
1303            a2fprefactor_in  = eigval_in (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol))
1304            a2fprefactor_out = eigval_out(ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol))
1305          end if
1306 
1307          omega = omega_min
1308          tmpa2f_in (:) = zero
1309          tmpa2f_out(:) = zero
1310          do iomega=1,elph_ds%na2f
1311            xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor
1312            gaussval = gaussprefactor*exp(-xx*xx)
1313 
1314            temp_in = gaussval*a2fprefactor_in
1315            temp_out = gaussval*a2fprefactor_out
1316 
1317            if (dabs(temp_in) < 1.0d-50) temp_in = zero
1318            if (dabs(temp_out) < 1.0d-50) temp_out = zero
1319            tmpa2f_in (iomega) = tmpa2f_in (iomega) + temp_in
1320            tmpa2f_out(iomega) = tmpa2f_out(iomega) + temp_out
1321            omega = omega+domega
1322          end do
1323 
1324          elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) + tmpa2f_in(:)
1325          elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) + tmpa2f_out(:)
1326 
1327        end do ! end ibranch do
1328      end do ! end itrtensor do
1329    end do ! end iFSqpt do
1330  end do ! end isppol
1331 
1332  ABI_DEALLOCATE(coskr)
1333  ABI_DEALLOCATE(sinkr)
1334 
1335 !second 1 / elph_ds%k_fine%nkpt factor for the integration weights
1336  elph_tr_ds%a2f_1d_trin  = elph_tr_ds%a2f_1d_trin  / elph_ds%k_fine%nkpt
1337  elph_tr_ds%a2f_1d_trout = elph_tr_ds%a2f_1d_trout / elph_ds%k_fine%nkpt
1338 
1339  if (elph_ds%ep_scalprod == 1) then
1340    write(std_out,*) 'mka2f_tr_lova: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr
1341  end if
1342 
1343  elph_tr_ds%a2f_1d_tr(:,:,:,1,1,1) = elph_tr_ds%a2f_1d_trout(:,:,:) - elph_tr_ds%a2f_1d_trin(:,:,:)
1344 
1345 !output the elph_tr_ds%a2f_1d_tr
1346  fname = trim(elph_ds%elph_base_name) // '_A2F_TR'
1347  if (open_file (fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then
1348    MSG_ERROR(message)
1349  end if
1350 
1351  fname = trim(elph_ds%elph_base_name) // '_A2F_TRIN'
1352  if (open_file(fname,message,newunit=unit_a2f_trin,status='unknown') /= 0) then
1353    MSG_ERROR(message)
1354  end if
1355 
1356  fname = trim(elph_ds%elph_base_name) // '_A2F_TROUT'
1357  if (open_file (fname,message,newunit=unit_a2f_trout,status='unknown') /=0) then
1358    MSG_ERROR(message)
1359  end if
1360 
1361  write (unit_a2f_tr,'(a)')       '#'
1362  write (unit_a2f_tr,'(a)')       '# ABINIT package : a2f_tr file'
1363  write (unit_a2f_tr,'(a)')       '#'
1364  write (unit_a2f_tr,'(a)')       '# a2f_tr function integrated over the FS. omega in a.u.'
1365  write (unit_a2f_tr,'(a,I10)')   '#     number of kpoints integrated over : ', elph_ds%k_fine%nkpt
1366  write (unit_a2f_tr,'(a,I10)')   '#     number of energy points : ',elph_ds%na2f
1367  write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '#       between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
1368  write (unit_a2f_tr,'(a,E16.6)') '#   and the smearing width for gaussians is ', elph_ds%a2fsmear
1369  write (unit_a2f_tr,'(a)')       '#'
1370 
1371  write (unit_a2f_trin,'(a)')       '#'
1372  write (unit_a2f_trin,'(a)')       '# ABINIT package : a2f_trin file'
1373  write (unit_a2f_trin,'(a)')       '#'
1374  write (unit_a2f_trin,'(a)')       '# a2f_trin function integrated over the FS. omega in a.u.'
1375  write (unit_a2f_trin,'(a,I10)')   '#     number of kpoints integrated over : ', elph_ds%k_fine%nkpt
1376  write (unit_a2f_trin,'(a,I10)')   '#     number of energy points : ',elph_ds%na2f
1377  write (unit_a2f_trin,'(a,E16.6,a,E16.6,a)') '#       between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
1378  write (unit_a2f_trin,'(a,E16.6)') '#   and the smearing width for gaussians is ', elph_ds%a2fsmear
1379  write (unit_a2f_trin,'(a)')       '#'
1380 
1381  write (unit_a2f_trout,'(a)')       '#'
1382  write (unit_a2f_trout,'(a)')       '# ABINIT package : a2f_trout file'
1383  write (unit_a2f_trout,'(a)')       '#'
1384  write (unit_a2f_trout,'(a)')       '# a2f_trout function integrated over the FS. omega in a.u.'
1385  write (unit_a2f_trout,'(a,I10)')   '#     number of kpoints integrated over : ', elph_ds%k_fine%nkpt
1386  write (unit_a2f_trout,'(a,I10)')   '#     number of energy points : ',elph_ds%na2f
1387  write (unit_a2f_trout,'(a,E16.6,a,E16.6,a)') '#       between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha'
1388  write (unit_a2f_trout,'(a,E16.6)') '#   and the smearing width for gaussians is ', elph_ds%a2fsmear
1389  write (unit_a2f_trout,'(a)')       '#'
1390 
1391 !done with header
1392  do isppol=1,elph_ds%nsppol
1393    write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol)
1394    write (unit_a2f_trin,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol)
1395    write (unit_a2f_trout,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol)
1396 !  omega = zero
1397    omega = omega_min
1398    do iomega=1,elph_ds%na2f
1399      write (unit_a2f_tr,   '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr   (iomega,:,isppol,1,1,1)
1400      write (unit_a2f_trin, '(10D16.6)') omega, elph_tr_ds%a2f_1d_trin (iomega,:,isppol)
1401      write (unit_a2f_trout,'(10D16.6)') omega, elph_tr_ds%a2f_1d_trout(iomega,:,isppol)
1402      omega=omega+domega
1403    end do
1404    write (unit_a2f_tr,*)
1405    write (unit_a2f_trin,*)
1406    write (unit_a2f_trout,*)
1407  end do !isppol
1408 
1409  close (unit=unit_a2f_tr)
1410  close (unit=unit_a2f_trin)
1411  close (unit=unit_a2f_trout)
1412 
1413 !calculation of transport properties
1414  ABI_ALLOCATE(integrho,(elph_ds%na2f))
1415  ABI_ALLOCATE(tointegrho,(elph_ds%na2f))
1416  ABI_ALLOCATE(tointega2f,(elph_ds%na2f))
1417  ABI_ALLOCATE(integtau,(elph_ds%na2f))
1418  ABI_ALLOCATE(tointegtau,(elph_ds%na2f))
1419 
1420  fname = trim(elph_ds%elph_base_name) // '_RHO'
1421  if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then
1422    MSG_ERROR(message)
1423  end if
1424 
1425 !print header to resistivity file
1426  write (unit_rho,*) '# Resistivity as a function of temperature.'
1427  write (unit_rho,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
1428  write (unit_rho,*) '#  '
1429  write (unit_rho,*) '#  Columns are: '
1430  write (unit_rho,*) '#  temperature[K]   rho[au]   rho [SI]        rho/temp [au]'
1431  write (unit_rho,*) '#  '
1432 
1433  fname = trim(elph_ds%elph_base_name) // '_TAU'
1434  if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then
1435    MSG_ERROR(message)
1436  end if
1437 
1438 !print header to relaxation time file
1439  write (unit_tau,*) '# Relaxation time as a function of temperature.'
1440  write (unit_tau,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
1441  write (unit_tau,*) '#  '
1442  write (unit_tau,*) '#  Columns are: '
1443  write (unit_tau,*) '#  temperature[K]   tau[au]   tau [femtosecond]     '
1444  write (unit_tau,*) '#  '
1445 
1446  fname = trim(elph_ds%elph_base_name) // '_WTH'
1447  if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then
1448    MSG_ERROR(message)
1449  end if
1450 
1451 !print header to thermal conductivity file
1452  write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.'
1453  write (unit_therm,'(a)') '#  the formalism is isotropic, so non-cubic crystals may be wrong'
1454  write (unit_therm,'(a)') '#  '
1455  write (unit_therm,'(a)') '#  Columns are: '
1456  write (unit_therm,'(a)') '#  temperature[K]   thermal rho[au]   thermal cond [au]   thermal rho [SI]   thermal cond [SI]'
1457  write (unit_therm,'(a)') '#  '
1458 
1459  fname = trim(elph_ds%elph_base_name) // '_LOR'
1460  if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then
1461    MSG_ERROR(message)
1462  end if
1463 
1464 !print header to lorentz file
1465  write (unit_lor,*) '# Lorentz number as a function of temperature.'
1466  write (unit_lor,*) '#  the formalism is isotropic, so non-cubic crystals may be wrong'
1467  write (unit_lor,*) '#  '
1468  write (unit_lor,*) '#  Columns are: '
1469  write (unit_lor,*) '#  temperature[K]   Lorentz number[au]   Lorentz quantum = (pi*kb_HaK)**2/3'
1470  write (unit_lor,*) '#  '
1471 
1472  do isppol=1,elph_ds%nsppol
1473    lambda_tr_trace = zero
1474    do itrtensor=1,9
1475      omega = omega_min
1476      tointega2f = zero
1477      do iomega=1,elph_ds%na2f
1478        if(omega<=0) then
1479          omega=omega+domega
1480          cycle
1481        end if
1482        tointega2f(iomega)=elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega
1483        omega=omega+domega
1484      end do
1485 
1486      integrho = zero
1487      call simpson_int(elph_ds%na2f,domega,tointega2f,integrho)
1488      lambda_tr = two * spinfact * integrho(elph_ds%na2f)
1489      write (message, '(a,2i3,a,es16.6)' )&
1490 &     ' mka2f_tr_lova : TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' =  ', lambda_tr
1491      call wrtout(std_out,message,'COLL')
1492      if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr
1493    end do !end itrtensor do
1494 
1495    lambda_tr_trace = lambda_tr_trace / three
1496    write (message, '(a,i3,a,es16.6)' )&
1497 &   ' mka2f_tr_lova: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' =  ', lambda_tr_trace
1498    call wrtout(std_out,message,'COLL')
1499    call wrtout(ab_out,message,'COLL')
1500  end do !end isppol do
1501 
1502 !constant to change units of rho from au to SI
1503  chgu=2.173969d-7
1504  chtu=2.4188843265d-17
1505  femto=1.0d-15
1506 
1507  do isppol=1,elph_ds%nsppol
1508    do icomp=1, 3
1509      do jcomp=1, 3
1510        itrtensor=(icomp-1)*3+jcomp
1511 
1512 !      prefactor for resistivity integral
1513 !      firh=6.d0*pi*crystal%ucvol*kb_HaK/(elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol))
1514 !      FIXME: check factor of 2 which is different from Savrasov paper. 6 below for thermal conductivity is correct.
1515        firh=2.d0*pi*crystal%ucvol*kb_HaK/elph_ds%n0(isppol)/&
1516 &       sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol))
1517 
1518 !      Add by BX to get Tau_elph
1519        firh_tau = 2.0d0*pi*kb_HaK
1520 !      End Adding
1521 
1522        write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor
1523        write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor
1524 
1525 ! jmb
1526        tointegtau(:)=0.
1527        tointegrho(:)=0.
1528        do itemp=1,ntemper  ! runs over termperature in K
1529          Temp=tempermin+temperinc*dble(itemp)
1530          firhT=firh*Temp
1531          firhT_tau=firh_tau*Temp
1532          omega = omega_min
1533          do iomega=1,elph_ds%na2f
1534            if(omega<=0) then
1535              omega=omega+domega
1536              cycle
1537            end if
1538            xtr=omega/(2*kb_HaK*Temp)
1539            if(xtr < log(huge(zero)*tol16)/2)then
1540              tointegrho(iomega)=spinfact*firhT*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)  &
1541 &             /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2)
1542 !            Add by BX to get Tau
1543              tointegtau(iomega)=spinfact*firhT_tau*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)  &
1544 &             /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2)
1545            else
1546              tointegrho(iomega)=zero
1547              tointegtau(iomega)=zero
1548            end if
1549            omega=omega+domega
1550          end do
1551 
1552          call simpson_int(elph_ds%na2f,domega,tointegrho,integrho)
1553          call simpson_int(elph_ds%na2f,domega,tointegtau,integtau)
1554          rho=integrho(elph_ds%na2f)
1555          tau=1.0d99
1556          if(dabs(integtau(elph_ds%na2f)) < tol7) then
1557            write(message,'(a)') ' Cannot get a physical relaxation time '
1558            MSG_WARNING(message)
1559          else
1560            tau=1.0d0/integtau(elph_ds%na2f)
1561          end if
1562 !         if(elph_ds%na2f < 350.0) then
1563 !           tau=1.0d0/integtau(elph_ds%na2f)
1564 !         end if
1565          write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp
1566          write(unit_tau,'(3D20.10)')temp,tau,tau*chtu/femto
1567          rho_T(itemp)=rho
1568          tau_T(itemp)=tau
1569        end do ! temperature
1570        write(unit_rho,*)
1571        write(unit_tau,*)
1572 
1573      end do ! jcomp
1574    end do ! icomp
1575  end do ! isppol
1576 
1577 !-----------------------------
1578 
1579 
1580  do isppol=1,elph_ds%nsppol
1581    do icomp=1, 3
1582      do jcomp=1, 3
1583        itrtensor=(icomp-1)*3+jcomp
1584 !      prefactor for integral of thermal conductivity
1585 !      firh=(18.*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol))
1586        firh=(6.d0*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol))/ &
1587 &       sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol))
1588 
1589 
1590        write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol
1591        write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol
1592 
1593        tointegrho(:)=0.
1594        do itemp=1,ntemper
1595 
1596          Temp=tempermin + temperinc*dble(itemp)
1597          omega = omega_min
1598          do iomega=1,elph_ds%na2f
1599            if(omega<=0) then
1600              omega=omega+domega
1601              cycle
1602            end if
1603            xtr=omega/(2*kb_HaK*Temp)
1604            if(xtr < log(huge(zero)*tol16)/2)then
1605              tointegrho(iomega) = spinfact*xtr**2/omega*&
1606 &             ( elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)+&
1607 &             4*xtr**2*elph_tr_ds%a2f_1d_trout(iomega,itrtensor,isppol)/pi**2+   &
1608 &             2*xtr**2*elph_tr_ds%a2f_1d_trin(iomega,itrtensor,isppol)/pi**2)  &
1609 &             /(((exp(xtr)-exp(-xtr))/2)**2)
1610            else
1611              tointegrho(iomega) = zero
1612            end if
1613            omega=omega+domega
1614          end do
1615 
1616          call simpson_int(elph_ds%na2f,domega,tointegrho,integrho)
1617          wtherm=integrho(elph_ds%na2f)*firh
1618 
1619          if(abs(wtherm) > tol12)then
1620            write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9
1621 
1622            lorentz=rho_T(itemp)/(wtherm*temp)
1623            write(unit_lor,*)temp,lorentz,lor0
1624          else
1625            write(unit_therm,'(5D20.10)')temp,zero,huge(one),zero,huge(one)
1626            write(unit_lor,*)temp,huge(one),lor0
1627          end if
1628 
1629        end do
1630        write(unit_therm,*)
1631        write(unit_lor,*)
1632      end do ! jcomp
1633    end do ! icomp
1634  end do !end isppol do
1635 
1636 
1637  ABI_DEALLOCATE(phfrq)
1638  ABI_DEALLOCATE(displ)
1639  ABI_DEALLOCATE(pheigvec)
1640  ABI_DEALLOCATE(rho_T)
1641  ABI_DEALLOCATE(tau_T)
1642 
1643  close (unit=unit_lor)
1644  close (unit=unit_rho)
1645  close (unit=unit_tau)
1646  close (unit=unit_therm)
1647 
1648  ABI_DEALLOCATE(integrho)
1649  ABI_DEALLOCATE(integtau)
1650  ABI_DEALLOCATE(tointega2f)
1651  ABI_DEALLOCATE(tointegrho)
1652  ABI_DEALLOCATE(tointegtau)
1653  ABI_DEALLOCATE(elph_tr_ds%a2f_1d_tr)
1654  ABI_DEALLOCATE(elph_tr_ds%a2f_1d_trin)
1655  ABI_DEALLOCATE(elph_tr_ds%a2f_1d_trout)
1656 
1657  ABI_DEALLOCATE(elph_tr_ds%gamma_qpt_trin)
1658  ABI_DEALLOCATE(elph_tr_ds%gamma_qpt_trout)
1659  ABI_DEALLOCATE(elph_tr_ds%gamma_rpt_trin)
1660  ABI_DEALLOCATE(elph_tr_ds%gamma_rpt_trout)
1661 
1662 !DEBUG
1663  write(std_out,*) ' mka2f_tr_lova : end '
1664 !ENDDEBUG
1665 
1666 end subroutine mka2f_tr_lova