TABLE OF CONTENTS
ABINIT/get_tau_k [ 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
SOURCE
1653 subroutine get_tau_k(Cryst,ifc,Bst,elph_ds,elph_tr_ds,eigenGS,max_occ) 1654 1655 !Arguments ------------------------------------ 1656 type(crystal_t),intent(in) :: Cryst 1657 type(ifc_type),intent(in) :: ifc 1658 type(ebands_t),intent(inout) :: Bst 1659 type(elph_type),intent(inout) :: elph_ds 1660 type(elph_tr_type), intent(inout) :: elph_tr_ds 1661 real(dp),intent(in) :: max_occ 1662 real(dp),intent(in) :: eigenGS(elph_ds%nband,elph_ds%k_phon%nkpt,elph_ds%nsppol) 1663 1664 !Local variables------------------------------- 1665 !scalars 1666 character(len=500) :: message 1667 character(len=fnlen) :: fname 1668 integer :: ntemper,nsppol,nbranch,nband,natom 1669 integer :: nkpt,nqpt,nkptirr,nqptirr,new_nkptirr 1670 integer :: isppol,iFSkpt,iFSqpt,iqpt,iqpt_fullbz,imqpt_fullbz,ikpt_kpq,ikpt_kmq 1671 integer :: iband,jband,jpband,jbeff,ibranch,jbranch,itemp 1672 integer :: irec,ierr,nrpt,ik_this_proc 1673 integer :: unit_tau,unit_invtau 1674 integer :: nene,nene_all,iene,iene_fine,unit_taue,unit_mfp 1675 integer :: icomp,jcomp,itensor 1676 integer :: ikpt_irr,iomega,unit_cond,unit_therm,unit_sbk 1677 integer :: nskip,nspline 1678 real(dp) :: occ_omega,occ_e 1679 real(dp) :: xx,Temp,therm_factor 1680 real(dp) :: factor,dfermide 1681 real(dp) :: e_k,chu_tau,rate_e,mfp_e 1682 real(dp) :: ene,enemin,enemax,deltaene 1683 real(dp) :: omega,omega_min,omega_max,domega 1684 real(dp) :: diagerr 1685 real(dp) :: chu_mfp,chu_cond,chu_cth,chu_sbk,femto 1686 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 1687 real(dp) :: eigval(elph_ds%nbranch),eigval2(elph_ds%nbranch) 1688 real(dp) :: imeigval(elph_ds%nbranch) 1689 real(dp) :: tmp_wtkpq, tmp_wtkmq, tol_wtk 1690 real(dp) :: yp1,ypn 1691 !arrays 1692 integer,allocatable :: FSfullpktofull(:,:),mqtofull(:) 1693 integer,allocatable :: kpttokpt(:,:,:) 1694 real(dp) :: cond_inv(3,3) 1695 real(dp),allocatable :: fermie(:) 1696 real(dp),allocatable :: tmp_eigenGS(:,:,:) 1697 real(dp),allocatable :: tmp_gkk_qpt(:,:,:),tmp_gkk_rpt(:,:,:),tmp_gkk_kpt(:,:) 1698 real(dp),allocatable :: tmp_gkk_kpt2(:,:,:), gkk_kpt(:,:,:) 1699 real(dp),allocatable :: tau_k(:,:,:,:),inv_tau_k(:,:,:,:),tmp_tau_k(:,:,:,:) 1700 real(dp),allocatable :: phfrq(:,:),pheigvec(:,:) 1701 real(dp),allocatable :: displ(:,:,:,:) 1702 real(dp),allocatable :: a2f_2d(:),a2f_2d2(:) 1703 real(dp),allocatable :: tmp_wtk(:,:,:,:),tmp2_wtk(:),tmp_wtk1(:),tmp_wtk2(:) 1704 real(dp),allocatable :: ene_pt(:),ene_ptfine(:),ff2(:) 1705 real(dp),allocatable :: wtq(:,:,:),tmp_wtq(:,:,:),tmp2_wtq(:,:) 1706 real(dp),allocatable :: dos_e(:,:) 1707 real(dp),allocatable :: coskr1(:,:),sinkr1(:,:) 1708 real(dp),allocatable :: coskr2(:,:),sinkr2(:,:) 1709 real(dp),allocatable :: cond_e(:,:,:,:),cond(:,:,:,:),sbk(:,:,:,:),seebeck(:,:,:,:),cth(:,:,:,:) 1710 1711 ! ************************************************************************* 1712 1713 write(std_out,*) 'get_tau_k : enter ' 1714 1715 nrpt = ifc%nrpt 1716 natom = cryst%natom 1717 1718 nsppol = elph_ds%nsppol 1719 nbranch = elph_ds%nbranch 1720 nband = elph_ds%ngkkband 1721 nkpt = elph_ds%k_phon%nkpt 1722 nqpt = elph_ds%nqpt_full 1723 nkptirr = elph_ds%k_phon%nkptirr 1724 new_nkptirr = elph_ds%k_phon%new_nkptirr 1725 nqptirr = elph_ds%nqptirred 1726 ntemper = elph_ds%ntemper 1727 nene = 2*elph_ds%na2f-1 ! only need e_k +- omega_max range, take deltaene=delta_oemga 1728 1729 chu_tau = 2.4188843265*1.0d-17 1730 chu_mfp = 5.291772*1.0d-11 1731 chu_cond = 4.59988159904764*1.0d6 1732 chu_cth = 1.078637439971599*1.0d4 1733 chu_sbk = 8.617343101*1.0d-5 1734 femto = 1.0d-15 1735 1736 tol_wtk = tol7/nkptirr/nband 1737 1738 ABI_MALLOC(fermie ,(ntemper)) 1739 ABI_MALLOC(tmp_gkk_qpt ,(2,nbranch**2,nqpt)) 1740 ABI_MALLOC(tmp_gkk_rpt ,(2,nbranch**2,nrpt)) 1741 ABI_MALLOC(tmp_gkk_kpt ,(2,nbranch**2)) 1742 ABI_MALLOC(tmp_gkk_kpt2 ,(2,nbranch,nbranch)) 1743 ABI_MALLOC(gkk_kpt ,(2,nbranch,nbranch)) 1744 ABI_MALLOC(a2f_2d, (nene)) 1745 ABI_MALLOC(a2f_2d2, (nene)) 1746 ABI_MALLOC(inv_tau_k, (ntemper,nsppol,nkpt,nband)) 1747 ABI_MALLOC(tau_k, (ntemper,nsppol,nkpt,nband)) 1748 ABI_MALLOC(tmp_tau_k ,(ntemper,nsppol,new_nkptirr,nband)) 1749 1750 if (elph_ds%gkqwrite == 0) then 1751 call wrtout(std_out,' get_tau_k : keeping gkq matrices in memory','COLL') 1752 else if (elph_ds%gkqwrite == 1) then 1753 fname=trim(elph_ds%elph_base_name) // '_GKKQ' 1754 write (message,'(2a)')' get_tau_k : reading gkq matrices from file ',trim(fname) 1755 call wrtout(std_out,message,'COLL') 1756 else 1757 write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite 1758 ABI_BUG(message) 1759 end if 1760 1761 !========================================================= 1762 !Get equivalence between a kpt_phon pair and a qpt in qpt_full 1763 !only works if the qpt grid is complete (identical to 1764 !the kpt one, with a basic shift of (0,0,0) 1765 !========================================================= 1766 1767 !mapping of k + q onto k' for k and k' in full BZ 1768 !for dense k grid 1769 ABI_MALLOC(FSfullpktofull,(nkpt,nkpt)) 1770 ABI_MALLOC(mqtofull,(nkpt)) 1771 1772 !kpttokpt(itim,isym,iqpt) = kpoint index which transforms to ikpt under isym and with time reversal itim. 1773 ABI_MALLOC(kpttokpt,(2,Cryst%nsym,nkpt)) 1774 1775 call wrtout(std_out,'get_tau_k: calling mkqptequiv to set up the FS kpoint set',"COLL") 1776 1777 call mkqptequiv (FSfullpktofull,Cryst,elph_ds%k_phon%kpt,nkpt,nkpt,kpttokpt,elph_ds%k_phon%kpt,mqtofull) 1778 1779 !========================================================= 1780 !========================================================= 1781 1782 omega_max = elph_ds%omega_max 1783 omega_min = elph_ds%omega_min 1784 domega = elph_ds%domega 1785 enemax = maxval(eigenGS(elph_ds%maxFSband,:,:)) 1786 enemin = minval(eigenGS(elph_ds%minFSband,:,:)) 1787 1788 if (enemin < (elph_ds%fermie-0.2)) then 1789 enemin = elph_ds%fermie-0.2 1790 end if 1791 if (enemax > (elph_ds%fermie+0.2)) then 1792 enemax = elph_ds%fermie+0.2 1793 end if 1794 1795 nspline = elph_ds%ep_nspline 1796 nene_all = INT((enemax-enemin+domega)/(nspline*domega)) + 1 1797 deltaene = domega 1798 write(std_out,*) 'E_min= ',enemin, 'E_max= ',enemax 1799 write(std_out,*) 'Number of energy points= ',nene_all 1800 write(std_out,'(a,I8)') 'scale factor for spline interpolation in RTA = ', elph_ds%ep_nspline 1801 write(std_out,*) 'delta_ene before spline interpolation= ',deltaene*nspline 1802 write(std_out,*) 'delta_ene after spline interpolation= ',deltaene 1803 write(std_out,*) 'Omega_min= ',omega_min, 'Omega_max= ',omega_max 1804 write(std_out,*) 'Number of phonon points= ',elph_ds%na2f 1805 write(std_out,*) 'delta_omega= ',domega 1806 write(std_out,*) 'number of bands= ', elph_ds%nband, nband 1807 1808 ABI_MALLOC(tmp_wtk,(nband,nkpt,nsppol,nene_all)) 1809 ABI_MALLOC(tmp2_wtk,(nene_all)) 1810 ABI_MALLOC(ff2,(nene_all)) 1811 ABI_MALLOC(ene_pt,(nene_all)) 1812 ABI_MALLOC(ene_ptfine,(nene_all*nspline)) 1813 ABI_MALLOC(tmp_wtk1,(nene_all*nspline)) 1814 ABI_MALLOC(tmp_wtk2,(nene_all*nspline)) 1815 ABI_MALLOC(dos_e,(nsppol,nene_all)) 1816 1817 !Get energy points for spline interpolation 1818 do iene = 1, nene_all 1819 ene_pt(iene) = enemin + (iene-1)*nspline*deltaene 1820 end do 1821 1822 do iene = 1, nene_all*nspline 1823 ene_ptfine(iene) = enemin + (iene-1)*deltaene 1824 end do 1825 1826 ABI_MALLOC(tmp_wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f+1)) 1827 ABI_MALLOC(wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f)) 1828 ABI_MALLOC(tmp2_wtq,(elph_ds%nbranch, elph_ds%na2f)) 1829 1830 !phonon 1831 ABI_MALLOC(phfrq,(nbranch, nkptirr)) 1832 ABI_MALLOC(displ,(2, nbranch, nbranch, nkptirr)) 1833 ABI_MALLOC(pheigvec,(2*nbranch*nbranch, nkptirr)) 1834 1835 do iFSqpt = 1, nkptirr 1836 call ifc%fourq(cryst,elph_ds%k_phon%kptirr(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 1837 end do 1838 1839 omega_min = omega_min - domega 1840 1841 !bxu, obtain wtq for the q_fine, then condense to q_phon 1842 call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,Cryst%gprimd,elph_ds%kptrlatt, & 1843 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_phon,tmp_wtq) 1844 omega_min = omega_min + domega 1845 1846 do iomega = 1, elph_ds%na2f 1847 wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1) 1848 !write(1005,*) omega_min+(iomega-1)*domega, sum(tmp_wtq(:,:,iomega+1))/nkpt 1849 end do 1850 ABI_FREE(tmp_wtq) 1851 1852 ! electron 1853 tmp_wtk =zero 1854 dos_e = zero 1855 call ep_el_weights(elph_ds%ep_b_min, elph_ds%ep_b_max, eigenGS(elph_ds%minFSband:elph_ds%minFSband+nband-1,:,:), & 1856 & elph_ds%elphsmear, & 1857 & enemin, enemax, nene_all, Cryst%gprimd, elph_ds%k_phon%irredtoGS, elph_ds%kptrlatt, max_occ, & 1858 & 1, nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk) 1859 !& elph_ds%minFSband, elph_ds%nband, elph_ds%nFSband, nsppol, elph_ds%telphint, elph_ds%k_phon, tmp_wtk) 1860 1861 do isppol = 1, nsppol 1862 do iene = 1, nene_all 1863 dos_e(isppol,iene) = sum(tmp_wtk(:,:,isppol,iene))/nkpt 1864 end do 1865 end do 1866 1867 ABI_MALLOC(coskr1, (nqpt,nrpt)) 1868 ABI_MALLOC(sinkr1, (nqpt,nrpt)) 1869 call ftgam_init(ifc%gprim, nqpt, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr1, sinkr1) 1870 ABI_MALLOC(coskr2, (nkptirr,nrpt)) 1871 ABI_MALLOC(sinkr2, (nkptirr,nrpt)) 1872 call ftgam_init(ifc%gprim, nkptirr, nrpt, elph_ds%k_phon%kpt, Ifc%rpt, coskr2, sinkr2) 1873 1874 !get fermie for itemp 1875 fermie = elph_ds%fermie 1876 do itemp=1,ntemper ! runs over termperature in K 1877 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 1878 1879 Bst%occopt = 3 1880 Bst%tsmear = Temp*kb_HaK 1881 call ebands_update_occ(Bst,-99.99_dp) 1882 write(message,'(a,f12.6,a,E20.12)')'At T=',Temp,' Fermi level is:',Bst%fermie 1883 call wrtout(std_out,message,'COLL') 1884 1885 if (abs(elph_ds%fermie) < tol10) then 1886 fermie(itemp) = Bst%fermie 1887 end if 1888 end do 1889 1890 inv_tau_k = zero 1891 !get a2f_2d = \sum_{q,nbranch,jband'} |gkk|^2*\delta(\epsilon_{k'j'}-\epsilon')*\delta(\omega_q-\omega) 1892 do isppol=1,nsppol 1893 write (std_out,*) '##############################################' 1894 write (std_out,*) 'get_tau_k : Treating spin polarization ', isppol 1895 write (std_out,*) '##############################################' 1896 1897 ! do iFSkpt =1,nkpt 1898 do ik_this_proc =1,elph_ds%k_phon%my_nkpt 1899 iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc) 1900 write (std_out,*) 'get_tau_k : working on kpt # ', iFSkpt, '/', nkpt 1901 do jband = 1, nband 1902 ! write(*,*)'i am here 1 ', isppol,iFSkpt,jband 1903 a2f_2d = zero 1904 a2f_2d2 = zero 1905 1906 !sum from here 1907 nskip = 0 1908 do jpband = 1, nband 1909 jbeff = jpband+(jband-1)*nband 1910 1911 if (elph_ds%gkqwrite == 0) then 1912 tmp_gkk_qpt(:,:,:) = elph_ds%gkk_qpt(:,jbeff,:,ik_this_proc,isppol,:) 1913 else if (elph_ds%gkqwrite == 1) then 1914 irec = (ik_this_proc-1)*elph_ds%k_phon%my_nkpt + iqpt 1915 if (iFSkpt == 1) then 1916 write (std_out,*) ' get_tau_k read record ', irec 1917 end if 1918 read (elph_ds%unitgkq,REC=irec) tmp_gkk_qpt(:,:,iqpt_fullbz) 1919 end if 1920 1921 !FT to real space 1922 call ftgam(Ifc%wghatm,tmp_gkk_qpt,tmp_gkk_rpt,natom,nqpt,nrpt,1,coskr1,sinkr1) 1923 1924 !sum over irred q over k_phon, with corresponding weights 1925 do iFSqpt = 1, nkptirr 1926 iqpt_fullbz = elph_ds%k_phon%irredtoGS(iFSqpt) 1927 ikpt_kpq = FSfullpktofull(iFSkpt,iqpt_fullbz) 1928 1929 imqpt_fullbz = mqtofull(iqpt_fullbz) 1930 ikpt_kmq = FSfullpktofull(iFSkpt,imqpt_fullbz) 1931 1932 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr 1933 call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(iqpt_fullbz,:),sinkr2(iqpt_fullbz,:)) 1934 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt) 1935 1936 !if ep_scalprod==0 we have to dot in the displacement vectors here 1937 if (elph_ds%ep_scalprod==0) then 1938 1939 call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red) 1940 1941 tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/)) 1942 call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt) 1943 1944 do jbranch=1,nbranch 1945 eigval(jbranch) = gkk_kpt(1, jbranch, jbranch) 1946 imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch) 1947 1948 if (abs(imeigval(jbranch)) > tol10) then 1949 write (message,'(a,i0,a,es16.8)')" real values branch = ",jbranch,' eigval = ',eigval(jbranch) 1950 ABI_WARNING(message) 1951 write (message,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch) 1952 ABI_WARNING(message) 1953 end if 1954 1955 end do 1956 1957 ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated. 1958 else if (elph_ds%ep_scalprod == 1) then 1959 1960 ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here 1961 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,& 1962 & pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom) 1963 1964 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,& 1965 & tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom) 1966 1967 diagerr = zero 1968 do ibranch=1,nbranch 1969 eigval(ibranch) = gkk_kpt(1,ibranch,ibranch) 1970 do jbranch=1,ibranch-1 1971 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 1972 end do 1973 do jbranch=ibranch+1,nbranch 1974 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 1975 end do 1976 end do 1977 1978 if (diagerr > tol12) then 1979 write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr 1980 ABI_WARNING(message) 1981 end if 1982 1983 else 1984 write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod 1985 ABI_BUG(message) 1986 end if ! end ep_scalprod if 1987 1988 !For k'=k-q 1989 !Do FT from real-space gamma grid to 1 kpt in k_phon%new_kptirr 1990 call ftgam(Ifc%wghatm,tmp_gkk_kpt,tmp_gkk_rpt,natom,1,nrpt,0,coskr2(imqpt_fullbz,:),sinkr2(imqpt_fullbz,:)) 1991 !tmp_gkk_kpt(:,:)=tmp_gkk_qpt(:,:,iFSqpt) 1992 1993 !if ep_scalprod==0 we have to dot in the displacement vectors here 1994 if (elph_ds%ep_scalprod==0) then 1995 1996 call phdispl_cart2red(natom,Cryst%gprimd,displ(:,:,:,iFSqpt),displ_red) 1997 1998 tmp_gkk_kpt2 = reshape (tmp_gkk_kpt(:,:), (/2,nbranch,nbranch/)) 1999 call gam_mult_displ(nbranch, displ_red, tmp_gkk_kpt2, gkk_kpt) 2000 2001 do jbranch=1,nbranch 2002 eigval2(jbranch) = gkk_kpt(1, jbranch, jbranch) 2003 imeigval(jbranch) = gkk_kpt(2, jbranch, jbranch) 2004 2005 if (abs(imeigval(jbranch)) > tol10) then 2006 write (message,'(a,i0,a,es16.8)')" real values branch = ",jbranch,' eigval = ',eigval2(jbranch) 2007 ABI_WARNING(message) 2008 write (message,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch) 2009 ABI_WARNING(message) 2010 end if 2011 2012 end do 2013 2014 ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated. 2015 else if (elph_ds%ep_scalprod == 1) then 2016 2017 ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here 2018 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, tmp_gkk_kpt, 3*natom,& 2019 & pheigvec(:,iFSqpt), 3*natom, czero, tmp_gkk_kpt2, 3*natom) 2020 2021 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec(:,iFSqpt), 3*natom,& 2022 & tmp_gkk_kpt2, 3*natom, czero, gkk_kpt, 3*natom) 2023 2024 diagerr = zero 2025 do ibranch=1,nbranch 2026 eigval2(ibranch) = gkk_kpt(1,ibranch,ibranch) 2027 do jbranch=1,ibranch-1 2028 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 2029 end do 2030 do jbranch=ibranch+1,nbranch 2031 diagerr = diagerr + abs(gkk_kpt(1,jbranch,ibranch)) 2032 end do 2033 end do 2034 2035 if (diagerr > tol12) then 2036 write(message,'(a,es15.8)') 'get_tau_k: residual in diagonalization of gamma with phon eigenvectors: ', diagerr 2037 ABI_WARNING(message) 2038 end if 2039 2040 else 2041 write (message,'(a,i0)')' Wrong value for ep_scalprod = ',elph_ds%ep_scalprod 2042 ABI_BUG(message) 2043 end if ! end ep_scalprod if 2044 2045 tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kpq,isppol,:) 2046 yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene 2047 ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene 2048 call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2) 2049 call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk1) 2050 2051 tmp2_wtk(:) = tmp_wtk(jpband,ikpt_kmq,isppol,:) 2052 yp1 = (tmp2_wtk(2)-tmp2_wtk(1))/nspline/deltaene 2053 ypn = (tmp2_wtk(nene_all)-tmp2_wtk(nene_all-1))/nspline/deltaene 2054 call spline(ene_pt,tmp2_wtk,nene_all,yp1,ypn,ff2) 2055 call splint(nene_all,ene_pt,tmp2_wtk,ff2,nene_all*nspline,ene_ptfine,tmp_wtk2) 2056 2057 tmp2_wtq(:,:) = wtq(:,iFSqpt,:) 2058 do iene=1,nene 2059 e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol) 2060 ene = e_k - omega_max + (iene-1)*deltaene 2061 if (ene<enemin .or. ene>enemax) cycle 2062 iene_fine = NINT((ene-enemin+deltaene)/deltaene) 2063 tmp_wtkpq = tmp_wtk1(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt) 2064 tmp_wtkmq = tmp_wtk2(iene_fine) * elph_ds%k_phon%wtkirr(iFSqpt) 2065 2066 if (tmp_wtkpq+tmp_wtkmq < tol_wtk ) then 2067 nskip = nskip +1 2068 cycle 2069 end if 2070 2071 do ibranch = 1, nbranch 2072 if (abs(phfrq(ibranch,iFSqpt)) < tol7) cycle 2073 2074 if (ene > e_k) then 2075 omega = ene - e_k 2076 if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle 2077 iomega = NINT((omega-omega_min+domega)/domega) 2078 2079 a2f_2d(iene) = a2f_2d(iene) +& 2080 & eigval(ibranch)/phfrq(ibranch,iFSqpt)*& 2081 & tmp_wtkpq * tmp2_wtq(ibranch,iomega) 2082 end if 2083 2084 if (ene < e_k) then 2085 omega = e_k - ene 2086 if (abs(omega) < tol7 .or. abs(omega) > omega_max) cycle 2087 iomega = NINT((omega-omega_min+domega)/domega) 2088 2089 a2f_2d2(iene) = a2f_2d2(iene) +& 2090 & eigval(ibranch)/phfrq(ibranch,iFSqpt)*& 2091 & tmp_wtkmq * tmp2_wtq(ibranch,iomega) 2092 end if 2093 2094 end do ! ibranch 3 2095 end do ! nene 800 2096 end do ! kptirr 216 2097 end do ! j' band 3 2098 ! print *, ' skipped ', nskip, ' energy points out of ', nene*nband*nkptirr 2099 2100 ! get inv_tau_k 2101 do itemp=1,ntemper ! runs over termperature in K 2102 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 2103 do iene=1,nene 2104 e_k = eigenGS(elph_ds%minFSband+jband-1,iFSkpt,isppol) 2105 ene = e_k - omega_max + (iene-1)*deltaene 2106 if (ene<enemin .or. ene>enemax) cycle 2107 2108 xx=(ene-fermie(itemp))/(kb_HaK*Temp) 2109 occ_e=1.0_dp/(exp(xx)+1.0_dp) 2110 if (ene > e_k .and. (ene-e_k) .le. omega_max) then 2111 omega = ene - e_k 2112 if (abs(omega) < tol7) cycle 2113 xx = omega/(kb_HaK*Temp) 2114 occ_omega=1.0_dp/(exp(xx)-1.0_dp) 2115 2116 therm_factor = occ_e + occ_omega 2117 2118 inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +& 2119 a2f_2d(iene)*therm_factor*deltaene 2120 end if 2121 if (ene < e_k .and. (e_k-ene) .le. omega_max) then 2122 omega = e_k - ene 2123 if (abs(omega) < tol7) cycle 2124 xx = omega/(kb_HaK*Temp) 2125 occ_omega=1.0_dp/(exp(xx)-1.0_dp) 2126 2127 therm_factor = 1 - occ_e + occ_omega 2128 2129 inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband) +& 2130 a2f_2d2(iene)*therm_factor*deltaene 2131 end if 2132 2133 end do ! nene 2134 end do ! Temp 2135 ! write(*,*)'i am here 2 ', isppol,iFSkpt,jband 2136 end do ! jband 2137 end do ! kpt 2138 end do ! nsppol 2139 2140 !write (300+mpi_enreg%me,*) inv_tau_k 2141 call xmpi_sum (inv_tau_k, xmpi_world, ierr) 2142 2143 ABI_FREE(phfrq) 2144 ABI_FREE(displ) 2145 ABI_FREE(pheigvec) 2146 ABI_FREE(tmp2_wtk) 2147 ABI_FREE(ff2) 2148 ABI_FREE(ene_pt) 2149 ABI_FREE(ene_ptfine) 2150 ABI_FREE(tmp_wtk1) 2151 ABI_FREE(tmp_wtk2) 2152 ABI_FREE(tmp2_wtq) 2153 ABI_FREE(wtq) 2154 ABI_FREE(coskr1) 2155 ABI_FREE(sinkr1) 2156 ABI_FREE(coskr2) 2157 ABI_FREE(sinkr2) 2158 ABI_FREE(kpttokpt) 2159 ABI_FREE(FSfullpktofull) 2160 ABI_FREE(mqtofull) 2161 ABI_FREE(tmp_gkk_qpt) 2162 ABI_FREE(tmp_gkk_rpt) 2163 ABI_FREE(tmp_gkk_kpt) 2164 ABI_FREE(tmp_gkk_kpt2) 2165 ABI_FREE(gkk_kpt) 2166 ABI_FREE(a2f_2d) 2167 ABI_FREE(a2f_2d2) 2168 2169 !output inv_tau_k and tau_k 2170 fname = trim(elph_ds%elph_base_name) // '_INVTAUK' 2171 if (open_file(fname,message,newunit=unit_invtau,status='unknown') /= 0) then 2172 ABI_ERROR(message) 2173 end if 2174 2175 !print header to relaxation time file 2176 write (unit_invtau,*) '# k-dep inverse of the relaxation time as a function of temperature.' 2177 write (unit_invtau,*) '# ' 2178 write (unit_invtau,*) '# nkptirr= ', nkptirr, 'nband= ', nband 2179 write (unit_invtau,*) '# number of temperatures= ', ntemper 2180 write (unit_invtau,*) '# tau [femtosecond^-1] ' 2181 2182 fname = trim(elph_ds%elph_base_name) // '_TAUK' 2183 if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then 2184 ABI_ERROR(message) 2185 end if 2186 2187 !print header to relaxation time file 2188 write (unit_tau,*) '# k-dep relaxation time as a function of temperature.' 2189 write (unit_tau,*) '# ' 2190 write (unit_tau,*) '# nkptirr= ', nkptirr, 'nband= ', nband 2191 write (unit_tau,*) '# number of temperatures= ', ntemper 2192 write (unit_tau,*) '# tau [femtosecond] ' 2193 2194 tau_k = zero 2195 do itemp=1,ntemper ! runs over termperature in K 2196 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 2197 write(unit_invtau,'(a,f16.8)') '# Temperature = ', Temp 2198 write(unit_tau,'(a,f16.8)') '# Temperature = ', Temp 2199 do isppol=1,nsppol 2200 write(unit_invtau,'(a,i6)') '# For isppol = ', isppol 2201 write(unit_tau,'(a,i6)') '# For isppol = ', isppol 2202 do iFSkpt = 1,nkpt 2203 !FIXME: check when tau_k is too small, whether there should be a phonon 2204 !scattering or not, and should tau_k be zero or not. 2205 do jband = 1,nband 2206 if (abs(inv_tau_k(itemp,isppol,iFSkpt,jband)) < tol9) then 2207 inv_tau_k(itemp,isppol,iFSkpt,jband) = zero 2208 tau_k(itemp,isppol,iFSkpt,jband) = zero 2209 else 2210 !no need to *nkpt due to wtkirr, as we need /nkpt for the sum 2211 !no need to *two_pi due to the missing prefactor in gkk (see mka2f_tr_lova) 2212 inv_tau_k(itemp,isppol,iFSkpt,jband) = inv_tau_k(itemp,isppol,iFSkpt,jband)*elph_ds%occ_factor 2213 tau_k(itemp,isppol,iFSkpt,jband) = one/inv_tau_k(itemp,isppol,iFSkpt,jband) 2214 end if 2215 end do ! nband 2216 write(unit_invtau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, ' kpt=', elph_ds%k_phon%kptirr(:,iFSkpt) 2217 write(unit_invtau,'(100D16.8)') (inv_tau_k(itemp,isppol,iFSkpt,iband)*femto/chu_tau,iband=1,nband) 2218 write(unit_tau,'(a,i8,a,3f12.6)') '# kpt# ', iFSkpt, ' kpt=', elph_ds%k_phon%kptirr(:,iFSkpt) 2219 write(unit_tau,'(100D16.8)') (tau_k(itemp,isppol,iFSkpt,iband)*chu_tau/femto,iband=1,nband) 2220 end do ! nkptirr 2221 write(unit_invtau,*) ' ' 2222 write(unit_tau,*) ' ' 2223 end do ! nsppol 2224 write(unit_invtau,*) ' ' 2225 write(unit_invtau,*) ' ' 2226 write(unit_tau,*) ' ' 2227 write(unit_tau,*) ' ' 2228 end do ! ntemper 2229 2230 ! Only use the irred k for eigenGS and tau_k 2231 ABI_MALLOC(tmp_eigenGS,(elph_ds%nband,elph_ds%k_phon%new_nkptirr,elph_ds%nsppol)) 2232 2233 do ikpt_irr = 1, new_nkptirr 2234 tmp_eigenGS(:,ikpt_irr,:) = eigenGS(:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:) 2235 tmp_tau_k(:,:,ikpt_irr,:) = tau_k(:,:,elph_ds%k_phon%new_irredtoGS(ikpt_irr),:)*chu_tau 2236 end do 2237 2238 !BoltzTraP output files in SIESTA format 2239 if (elph_ds%prtbltztrp == 1) then 2240 !Prevent use in case occopt = 9 2241 if (Bst%occopt==9) then 2242 ABI_ERROR("Boltztrap outputting not possible with occopt = 9 at the moment") 2243 end if 2244 call ebands_prtbltztrp_tau_out (tmp_eigenGS(elph_ds%minFSband:elph_ds%maxFSband,:,:),& 2245 & elph_ds%tempermin,elph_ds%temperinc,ntemper,fermie, & 2246 & elph_ds%elph_base_name,elph_ds%k_phon%new_kptirr,nband,elph_ds%nelect,new_nkptirr, & 2247 & elph_ds%nspinor,nsppol,Cryst%nsym,Cryst%rprimd,Cryst%symrel,tmp_tau_k) 2248 end if !prtbltztrp 2249 ABI_FREE(tmp_eigenGS) 2250 ABI_FREE(tmp_tau_k) 2251 2252 !Get the energy dependence of tau. 2253 !Eq. (6) in Restrepo et al. Appl. Phys. Lett. 94, 212103 (2009) [[cite:Restrepo2009]] 2254 2255 fname = trim(elph_ds%elph_base_name) // '_TAUE' 2256 if (open_file(fname,message,newunit=unit_taue,status='unknown') /= 0) then 2257 ABI_ERROR(message) 2258 end if 2259 2260 !print header to relaxation time file 2261 write (unit_taue,*) '# Energy-dep relaxation time as a function of temperature.' 2262 write (unit_taue,*) '# ' 2263 write (unit_taue,*) '# number of temperatures= ', ntemper 2264 write (unit_taue,*) '# ene[Ha] tau [femtosecond] DOS[au] ' 2265 2266 fname = trim(elph_ds%elph_base_name) // '_MFP' 2267 if (open_file(fname,message,newunit=unit_mfp,status='unknown') /= 0) then 2268 ABI_ERROR(message) 2269 end if 2270 2271 write (unit_mfp,*) '# Energy-dep mean free path as a function of temperature.' 2272 write (unit_mfp,*) '# ' 2273 write (unit_mfp,*) '# number of temperatures= ', ntemper 2274 write (unit_mfp,*) '# ene[Ha] mfp [femtometer] ' 2275 2276 do itemp=1,ntemper ! runs over termperature in K 2277 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 2278 write(unit_taue,'(a,f16.8)') '# Temperature = ', Temp 2279 do isppol = 1, nsppol 2280 write(unit_taue,*) '# Tau_e for isppol = ',isppol 2281 do iene = 1, nene_all 2282 rate_e = zero 2283 do iFSkpt = 1, nkpt 2284 do jband = 1, nband 2285 rate_e = rate_e + inv_tau_k(itemp,isppol,iFSkpt,jband)* & 2286 & tmp_wtk(jband,iFSkpt,isppol,iene) 2287 end do ! jband 2288 end do ! kpt 2289 if (dabs(dos_e(isppol,iene)) < tol7) then 2290 rate_e = zero 2291 else 2292 rate_e = rate_e/nkpt/dos_e(isppol,iene) 2293 end if 2294 write(unit_taue,"(3D16.8)") enemin+(iene-1)*deltaene*nspline, rate_e*femto/chu_tau, dos_e(isppol,iene) 2295 end do ! number of energies 2296 write(unit_taue,*) ' ' 2297 end do ! nsppol 2298 write(unit_taue,*) ' ' 2299 end do ! ntemperature 2300 2301 ! calculate and output mean free path 2302 do itemp=1,ntemper ! runs over termperature in K 2303 Temp=elph_ds%tempermin+elph_ds%temperinc*dble(itemp) 2304 write(unit_mfp,'(a,f16.8)') '# Temperature = ', Temp 2305 do isppol = 1, nsppol 2306 do icomp = 1, 3 2307 write(unit_mfp,*) '# Mean free path for isppol, icomp= ',isppol,icomp 2308 do iene = 1, nene_all 2309 mfp_e = zero 2310 do iFSkpt = 1, nkptirr 2311 do jband = 1, nband 2312 mfp_e = mfp_e + tau_k(itemp,isppol,iFSkpt,jband)* & 2313 & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* & 2314 & tmp_wtk(jband,iFSkpt,isppol,iene) 2315 !& elph_ds%k_phon%new_wtkirr(iFSqpt) 2316 end do ! jband 2317 end do ! kpt 2318 if (dabs(dos_e(isppol,iene)) < tol7) then 2319 mfp_e = zero 2320 else 2321 mfp_e = mfp_e/nkptirr/dos_e(isppol,iene) 2322 end if 2323 write(unit_mfp,"(2D16.8)") enemin+(iene-1)*deltaene*nspline, mfp_e*chu_mfp/femto 2324 end do ! number of energies 2325 write(unit_mfp,*) ' ' 2326 end do ! icomp 2327 write(unit_mfp,*) ' ' 2328 end do ! nsppol 2329 write(unit_mfp,*) ' ' 2330 end do ! ntemperature 2331 2332 ABI_MALLOC(cond_e ,(ntemper,nsppol,nene_all,9)) 2333 2334 !get cond_e 2335 cond_e = zero 2336 do itemp=1,ntemper ! runs over termperature in K 2337 do isppol = 1, nsppol 2338 do iene = 1, nene_all 2339 ! do iFSkpt =1,nkpt 2340 do ik_this_proc =1,elph_ds%k_phon%my_nkpt 2341 iFSkpt = elph_ds%k_phon%my_ikpt(ik_this_proc) 2342 do jband = 1, nband 2343 do icomp = 1, 3 2344 do jcomp = 1, 3 2345 itensor = (icomp-1)*3+jcomp 2346 cond_e(itemp,isppol,iene,itensor) = cond_e(itemp,isppol,iene,itensor) + & 2347 & tau_k(itemp,isppol,iFSkpt,jband)* & 2348 & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,icomp,isppol)* & 2349 & elph_tr_ds%el_veloc(iFSkpt,elph_ds%minFSband+jband-1,jcomp,isppol)* & 2350 & tmp_wtk(jband,iFSkpt,isppol,iene) 2351 end do 2352 end do 2353 end do ! jband 2354 end do ! kpt 2355 end do ! number of energies 2356 end do ! nsppol 2357 end do ! ntemperature 2358 2359 ! MG FIXME: Why xmpi_world, besides only master should perform IO in the section below. 2360 call xmpi_sum (cond_e, xmpi_world, ierr) 2361 2362 cond_e = cond_e/nkpt 2363 2364 !get transport coefficients 2365 2366 fname = trim(elph_ds%elph_base_name) // '_COND' 2367 if (open_file(fname,message,newunit=unit_cond,status='unknown') /= 0) then 2368 ABI_ERROR(message) 2369 end if 2370 2371 !print header to conductivity file 2372 write (unit_cond,*) '# Conductivity as a function of temperature.' 2373 write (unit_cond,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 2374 write (unit_cond,*) '# ' 2375 write (unit_cond,*) '# Columns are: ' 2376 write (unit_cond,*) '# temperature[K] cond[au] cond [SI] ' 2377 write (unit_cond,*) '# ' 2378 2379 fname = trim(elph_ds%elph_base_name) // '_CTH' 2380 if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then 2381 ABI_ERROR(message) 2382 end if 2383 2384 !print header to thermal conductivity file 2385 write (unit_therm,'(a)') '# Thermal conductivity as a function of temperature.' 2386 write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong' 2387 write (unit_therm,'(a)') '# ' 2388 write (unit_therm,'(a)') '# Columns are: ' 2389 write (unit_therm,'(a)') '# temperature[K] thermal cond [au] thermal cond [SI]' 2390 write (unit_therm,'(a)') '# ' 2391 2392 fname = trim(elph_ds%elph_base_name) // '_SBK' 2393 if (open_file(fname,message,newunit=unit_sbk,status='unknown') /=0) then 2394 ABI_ERROR(message) 2395 end if 2396 2397 !print header to relaxation time file 2398 write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.' 2399 write (unit_sbk,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 2400 write (unit_sbk,*) '# ' 2401 write (unit_sbk,*) '# Columns are: ' 2402 write (unit_sbk,*) '# temperature[K] S [au] S [SI] ' 2403 write (unit_sbk,*) '# ' 2404 2405 ABI_MALLOC(cond ,(ntemper,nsppol,3,3)) 2406 ABI_MALLOC(cth ,(ntemper,nsppol,3,3)) 2407 ABI_MALLOC(sbk ,(ntemper,nsppol,3,3)) 2408 ABI_MALLOC(seebeck ,(ntemper,nsppol,3,3)) 2409 2410 cond = zero 2411 cth = zero 2412 sbk = zero 2413 seebeck = zero 2414 do isppol=1,nsppol 2415 do icomp=1, 3 2416 do jcomp=1, 3 2417 itensor=(icomp-1)*3+jcomp 2418 do itemp=1,ntemper 2419 Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp) 2420 do iene = 1, nene_all 2421 factor = (enemin+(iene-1)*deltaene*nspline - fermie(itemp))/(kb_HaK*Temp) 2422 if (factor < -40.0d0) then 2423 dfermide = zero 2424 else if (factor > 40.0d0) then 2425 dfermide = zero 2426 else 2427 dfermide = EXP(factor)/(kb_HaK*Temp*(EXP(factor)+one)**2.0d0) 2428 end if 2429 cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp) + & 2430 & cond_e(itemp,isppol,iene,itensor)*dfermide*deltaene*nspline 2431 cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* & 2432 & (enemin+(iene-1)*deltaene*nspline - fermie(itemp))**2.0d0*dfermide*deltaene*nspline 2433 sbk(itemp,isppol,icomp,jcomp) = sbk(itemp,isppol,icomp,jcomp) + cond_e(itemp,isppol,iene,itensor)* & 2434 & (enemin+(iene-1)*deltaene*nspline - fermie(itemp))*dfermide*deltaene*nspline 2435 end do 2436 end do ! temperature 2437 end do ! jcomp 2438 end do ! icomp 2439 end do !end isppol 2440 2441 do isppol=1,nsppol 2442 do itemp=1,ntemper 2443 cond_inv(:,:)=cond(itemp,isppol,:,:) 2444 call matrginv(cond_inv,3,3) 2445 call DGEMM('N','N',3,3,3,one,sbk(itemp,isppol,:,:),3,cond_inv,& 2446 & 3,zero,seebeck(itemp,isppol,:,:),3) 2447 end do 2448 end do 2449 2450 do isppol=1,nsppol 2451 do icomp=1, 3 2452 do jcomp=1, 3 2453 itensor=(icomp-1)*3+jcomp 2454 write(unit_cond,*) '# Conductivity for isppol, itrten= ',isppol,itensor 2455 write(unit_therm,*) '# Thermal conductivity for isppol, itrten= ',isppol,itensor 2456 write(unit_sbk,*) '# Seebeck coefficient for isppol, itrten= ',isppol,itensor 2457 do itemp=1,ntemper 2458 Temp=elph_ds%tempermin + elph_ds%temperinc*dble(itemp) 2459 2460 seebeck(itemp,isppol,icomp,jcomp) = -1.0d0*seebeck(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp) 2461 cond(itemp,isppol,icomp,jcomp) = cond(itemp,isppol,icomp,jcomp)/cryst%ucvol 2462 cth(itemp,isppol,icomp,jcomp) = cth(itemp,isppol,icomp,jcomp)/(kb_HaK*Temp)/cryst%ucvol 2463 write(unit_cond,'(3D20.10)')Temp,cond(itemp,isppol,icomp,jcomp),cond(itemp,isppol,icomp,jcomp)*chu_cond 2464 write(unit_therm,'(3D20.10)')Temp,cth(itemp,isppol,icomp,jcomp),cth(itemp,isppol,icomp,jcomp)*chu_cth 2465 write(unit_sbk,'(3D20.10)')Temp,seebeck(itemp,isppol,icomp,jcomp),seebeck(itemp,isppol,icomp,jcomp)*chu_sbk 2466 end do ! temperature 2467 write(unit_cond,*) 2468 write(unit_therm,*) 2469 write(unit_sbk,*) 2470 end do ! jcomp 2471 end do ! icomp 2472 end do !end isppol 2473 2474 2475 ABI_FREE(inv_tau_k) 2476 ABI_FREE(tau_k) 2477 ABI_FREE(tmp_wtk) 2478 ABI_FREE(dos_e) 2479 ABI_FREE(cond_e) 2480 ABI_FREE(fermie) 2481 ABI_FREE(cond) 2482 ABI_FREE(sbk) 2483 ABI_FREE(cth) 2484 ABI_FREE(seebeck) 2485 2486 close (unit=unit_tau) 2487 close (unit=unit_taue) 2488 close (unit=unit_mfp) 2489 close (unit=unit_invtau) 2490 close (unit=unit_cond) 2491 close (unit=unit_therm) 2492 close (unit=unit_sbk) 2493 2494 end subroutine get_tau_k
ABINIT/m_a2ftr [ Modules ]
NAME
m_a2ftr
FUNCTION
COPYRIGHT
Copyright (C) 2004-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_a2ftr 23 24 use defs_basis 25 use defs_elphon 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_splines 30 use m_ebands 31 32 use defs_datatypes, only : ebands_t 33 use m_io_tools, only : open_file 34 use m_numeric_tools, only : simpson_int 35 use m_hide_lapack, only : matrginv 36 use m_geometry, only : phdispl_cart2red 37 use m_crystal, only : crystal_t 38 use m_ifc, only : ifc_type 39 use m_dynmat, only : ftgam_init, ftgam 40 use m_epweights, only : d2c_wtq, ep_ph_weights, ep_el_weights, ep_ph_weights 41 42 implicit none 43 44 private
ABINIT/mka2f_tr [ 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
NOTES
copied from ftiaf9.f
SOURCE
95 subroutine mka2f_tr(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,pair2red,elph_tr_ds) 96 97 !Arguments ------------------------------------ 98 !scalars 99 integer,intent(in) :: ntemper 100 real(dp),intent(in) :: tempermin,temperinc 101 type(ifc_type),intent(in) :: ifc 102 type(crystal_t),intent(in) :: crystal 103 type(elph_tr_type),intent(inout) :: elph_tr_ds 104 type(elph_type),intent(inout) :: elph_ds 105 !arrays 106 integer,intent(in) :: pair2red(elph_ds%nenergy,elph_ds%nenergy) 107 108 !Local variables ------------------------- 109 !x =w/(2kbT) 110 !scalars 111 integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr 112 integer :: unit_a2f_tr,natom,ii,jj 113 integer :: idir, iatom, k1, kdir 114 integer :: unit_lor,unit_rho,unit_tau,unit_sbk, unit_therm 115 integer :: itemp, tmp_nenergy 116 integer :: itrtensor, icomp, jcomp!, kcomp 117 integer :: ie, ie_1, ie2, ie_2, ie1, ie_tmp, ssp, s1(4), s2(4) 118 integer :: ie2_left, ie2_right 119 integer :: ik_this_proc, ierr,nrpt 120 logical,parameter :: debug=.False. 121 real(dp) :: Temp,chgu,chtu,chsu,chwu,diagerr,ucvol 122 real(dp) :: a2fprefactor, gtemp 123 real(dp) :: lambda_tr,lor0,lorentz,maxerr,omega 124 real(dp) :: rho,tau,wtherm,xtr 125 real(dp) :: lambda_tr_trace 126 real(dp) :: domega, omega_min, omega_max 127 real(dp) :: gaussval, gaussprefactor, gaussfactor, gaussmaxarg, xx 128 real(dp) :: qnorm2, tmp_fermie 129 real(dp) :: e1, e2, diff, xe 130 real(dp) :: occ_omega, occ_e1, occ_e2 131 real(dp) :: nv1, nv2, sigma1, sigma2 132 real(dp) :: dos_n_e2, veloc_sq_icomp, veloc_sq_jcomp 133 real(dp) :: tointegq00_1, tointegq00_2, tointegq01_1, tointegq01_2,tointegq11_1, tointegq11_2 134 real(dp) :: j00, j01, j11 135 real(dp) :: tointegq00,tointegq01,tointegq11 136 real(dp) :: pref_s, pref_w, tmp_veloc_sq0, tmp_veloc_sq1, tmp_veloc_sq2 137 character(len=500) :: message 138 character(len=fnlen) :: fname 139 !arrays 140 complex(dpc),parameter :: c0=dcmplx(0.d0,0.d0),c1=dcmplx(1.d0,0.d0) 141 real(dp) :: gprimd(3,3) 142 real(dp) :: eigval(elph_ds%nbranch) 143 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 144 real(dp) :: gam_now(2,elph_ds%nbranch*elph_ds%nbranch) 145 real(dp) :: tmpa2f(elph_ds%na2f) 146 real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch) 147 real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch) 148 real(dp) :: q11_inv(3,3) 149 !real(dp) :: fullq(6,6) 150 real(dp),allocatable :: phfrq(:,:) 151 real(dp),allocatable :: tmp_a2f_1d_tr(:,:,:,:,:) 152 real(dp),allocatable :: displ(:,:,:,:) 153 real(dp),allocatable :: pheigvec(:,:) 154 real(dp),allocatable :: tmp_wtq(:,:,:) 155 real(dp),allocatable :: integrho(:), tointegrho(:) 156 real(dp),allocatable :: integrand_q00(:),integrand_q01(:),integrand_q11(:) 157 real(dp),allocatable :: q00(:,:,:,:), q01(:,:,:,:),q11(:,:,:,:) 158 real(dp),allocatable :: seebeck(:,:,:,:)!, rho_nm(:,:,:,:) 159 real(dp),allocatable :: rho_T(:) 160 real(dp),allocatable :: coskr(:,:), sinkr(:,:) 161 162 ! ********************************************************************* 163 !calculate a2f_tr for frequencies between 0 and omega_max 164 165 166 write(std_out,*) 'mka2f_tr : enter ' 167 168 ucvol = crystal%ucvol 169 natom = crystal%natom 170 gprimd = crystal%gprimd 171 172 ! number of real-space points for FT interpolation 173 nrpt = ifc%nrpt 174 ! 175 !MG: the step should be calculated locally using nomega and the extrema of the spectrum. 176 !One should not rely on previous calls for the setup of elph_ds%domega 177 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t 178 domega =elph_ds%domega 179 omega_min = elph_ds%omega_min 180 omega_max = elph_ds%omega_max 181 182 if (elph_ds%ep_lova .eq. 1) then 183 tmp_nenergy = 1 184 else if (elph_ds%ep_lova .eq. 0) then 185 tmp_nenergy = elph_ds%nenergy 186 end if 187 188 !! defaults for number of temperature steps and max T (all in Kelvin...) 189 !ntemper=1000 190 !tempermin=zero 191 !temperinc=one 192 193 ABI_MALLOC(rho_T,(ntemper)) 194 195 196 gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear 197 gaussfactor = one / elph_ds%a2fsmear 198 gaussmaxarg = sqrt(-log(1.d-90)) 199 !lor0=(pi*kb_HaK)**2/3. 200 lor0=pi**2/3.0_dp 201 202 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F) 203 !WARNING! supposes this value has been set in mkelph_linwid. 204 205 !ENDMG 206 207 maxerr=0. 208 nerr=0 209 210 ABI_MALLOC(tmp_wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f+1)) 211 ABI_MALLOC(elph_ds%k_fine%wtq,(elph_ds%nbranch, elph_ds%k_fine%nkpt, elph_ds%na2f)) 212 ABI_MALLOC(elph_ds%k_phon%wtq,(elph_ds%nbranch, elph_ds%k_phon%nkpt, elph_ds%na2f)) 213 214 ABI_MALLOC(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt)) 215 ABI_MALLOC(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt)) 216 ABI_MALLOC(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt)) 217 218 do iFSqpt=1,elph_ds%k_fine%nkpt 219 call ifc%fourq(crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 220 end do 221 omega_min = omega_min - domega 222 223 !bxu, obtain wtq for the q_fine, then condense to q_phon 224 call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, & 225 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_fine,tmp_wtq) 226 !call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,elph_ds%na2f+1,gprimd,elph_ds%kptrlatt_fine, & 227 !& elph_ds%nbranch,1,elph_ds%k_fine,tmp_wtq) 228 omega_min = omega_min + domega 229 230 do iomega = 1, elph_ds%na2f 231 elph_ds%k_fine%wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1) 232 end do 233 ABI_FREE(tmp_wtq) 234 235 if (elph_ds%use_k_fine == 1) then 236 call d2c_wtq(elph_ds) 237 end if 238 239 ABI_FREE(phfrq) 240 ABI_FREE(displ) 241 ABI_FREE(pheigvec) 242 243 !reduce the dimension from fine to phon for phfrq and pheigvec 244 ABI_MALLOC(phfrq,(elph_ds%nbranch, elph_ds%k_phon%nkpt)) 245 ABI_MALLOC(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_phon%nkpt)) 246 ABI_MALLOC(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_phon%nkpt)) 247 248 do iFSqpt=1,elph_ds%k_phon%nkpt 249 call ifc%fourq(crystal,elph_ds%k_phon%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 250 end do 251 252 ABI_MALLOC(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2,ntemper)) 253 ABI_MALLOC(tmp_a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,4,tmp_nenergy**2)) 254 255 ! prepare phase factors 256 ABI_MALLOC(coskr, (elph_ds%k_phon%nkpt, nrpt)) 257 ABI_MALLOC(sinkr, (elph_ds%k_phon%nkpt, nrpt)) 258 call ftgam_init(ifc%gprim, elph_ds%k_phon%nkpt, nrpt, elph_ds%k_phon%kpt, ifc%rpt, coskr, sinkr) 259 260 elph_tr_ds%a2f_1d_tr = zero 261 tmp_a2f_1d_tr = zero 262 263 264 do ie = 1, elph_ds%n_pair 265 do ssp = 1,4 266 do isppol = 1, elph_ds%nsppol 267 268 ! loop over qpoint in full kpt grid (presumably dense) 269 do ik_this_proc =1,elph_ds%k_phon%my_nkpt 270 iFSqpt = elph_ds%k_phon%my_ikpt(ik_this_proc) 271 272 qnorm2 = sum(elph_ds%k_phon%kpt(:,iFSqpt)**2) 273 ! if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero 274 do itrtensor=1,9 275 276 if (elph_ds%ep_int_gkk == 1) then 277 gam_now(:,:) = elph_tr_ds%gamma_qpt_tr(:,itrtensor,:,isppol,iFSqpt) 278 else 279 ! Do FT from real-space gamma grid to 1 qpt. 280 call ftgam(ifc%wghatm,gam_now,elph_tr_ds%gamma_rpt_tr(:,itrtensor,:,isppol,:,ssp,ie),natom,1,nrpt,0,& 281 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 282 end if 283 284 ! Diagonalize gamma matrix at this qpoint (complex matrix). 285 286 ! if ep_scalprod==0 we have to dot in the displacement vectors here 287 if (elph_ds%ep_scalprod==0) then 288 289 displ_red(:,:,:) = zero 290 do jbranch=1,elph_ds%nbranch 291 do iatom=1,natom 292 do idir=1,3 293 ibranch=idir+3*(iatom-1) 294 do kdir=1,3 295 k1 = kdir+3*(iatom-1) 296 displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + & 297 & gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt) 298 displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + & 299 & gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt) 300 end do 301 end do 302 end do 303 end do 304 305 tmpgam2 = reshape (gam_now, (/2,elph_ds%nbranch,elph_ds%nbranch/)) 306 call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1) 307 do jbranch=1,elph_ds%nbranch 308 eigval(jbranch) = tmpgam1(1, jbranch, jbranch) 309 end do 310 311 else if (elph_ds%ep_scalprod == 1) then 312 313 314 ! NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care 315 316 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now, 3*natom,& 317 & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom) 318 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,& 319 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 320 diagerr = zero 321 do ibranch=1,elph_ds%nbranch 322 eigval(ibranch) = tmpgam2(1,ibranch,ibranch) 323 do jbranch=1,ibranch-1 324 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 325 end do 326 do jbranch=ibranch+1,elph_ds%nbranch 327 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 328 end do 329 end do 330 if (diagerr > tol12) then 331 nerr=nerr+1 332 maxerr=max(diagerr, maxerr) 333 end if 334 end if ! end ep_scalprod if 335 336 ! Add all contributions from the phonon modes at this qpoint to a2f and the phonon dos. 337 do ibranch=1,elph_ds%nbranch 338 ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then 339 if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. & 340 & (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then 341 ! note: this should depend on the velocity of sound, to accept acoustic modes! 342 a2fprefactor = zero 343 else 344 ! a2fprefactor = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol)) 345 ! Use the dos_n0 at the lowest input temperature, assuming to be low 346 a2fprefactor = eigval (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))) 347 end if 348 349 omega = omega_min 350 tmpa2f(:) = zero 351 do iomega=1,elph_ds%na2f 352 xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor 353 omega = omega+domega 354 if (abs(xx) > gaussmaxarg) cycle 355 356 gaussval = gaussprefactor*exp(-xx*xx) 357 gtemp = gaussval*a2fprefactor 358 359 if (dabs(gtemp) < 1.0d-50) gtemp = zero 360 tmpa2f(iomega) = tmpa2f(iomega) + gtemp 361 end do 362 363 ! tmpa2f(:) = zero 364 ! do iomega=1,elph_ds%na2f 365 ! gtemp = a2fprefactor*elph_ds%k_phon%wtq(ibranch,iFSqpt,iomega) 366 ! if (dabs(gtemp) < 1.0d-50) gtemp = zero 367 ! tmpa2f(iomega) = tmpa2f(iomega) + gtemp 368 ! end do 369 370 tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) = tmp_a2f_1d_tr (:,itrtensor,isppol,ssp,ie) + tmpa2f(:) 371 372 end do ! end ibranch 373 end do ! end itrtensor 374 end do ! end iFSqpt - loop done in parallel 375 end do ! end isppol 376 end do ! ss' 377 end do ! n_pair 378 379 ! MG: FIXME: Why xmpi_world? besides only one CPU should perform IO (see below) 380 ! Likely this routine is never executed in parallel 381 call xmpi_sum (tmp_a2f_1d_tr, xmpi_world, ierr) 382 383 ABI_FREE(coskr) 384 ABI_FREE(sinkr) 385 386 do itemp=1,ntemper ! runs over termperature in K 387 do isppol=1,elph_ds%nsppol 388 do jj=1,tmp_nenergy**2 389 do ii=1,4 390 elph_tr_ds%a2f_1d_tr(:,:,isppol,ii,jj,itemp) = tmp_a2f_1d_tr(:,:,isppol,ii,jj)/elph_tr_ds%dos_n0(itemp,isppol) 391 end do 392 end do 393 end do 394 end do 395 396 ABI_FREE(tmp_a2f_1d_tr) 397 398 !second 1 / elph_ds%k_phon%nkpt factor for the integration weights 399 elph_tr_ds%a2f_1d_tr = elph_tr_ds%a2f_1d_tr / elph_ds%k_phon%nkpt 400 401 if (elph_ds%ep_scalprod == 1) then 402 write(std_out,*) 'mka2f_tr: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr 403 end if 404 405 !output the elph_tr_ds%a2f_1d_tr 406 fname = trim(elph_ds%elph_base_name) // '_A2F_TR' 407 if (open_file(fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then 408 ABI_ERROR(message) 409 end if 410 411 write (unit_a2f_tr,'(a)') '#' 412 write (unit_a2f_tr,'(a)') '# ABINIT package : a2f_tr file' 413 write (unit_a2f_tr,'(a)') '#' 414 write (unit_a2f_tr,'(a)') '# a2f_tr function integrated over the FS. omega in a.u.' 415 write (unit_a2f_tr,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_phon%nkpt 416 write (unit_a2f_tr,'(a,I10)') '# number of energy points : ',elph_ds%na2f 417 write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 418 write (unit_a2f_tr,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 419 write (unit_a2f_tr,'(a)') '#' 420 421 !done with header 422 do isppol=1,elph_ds%nsppol 423 write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_tr_ds%dos_n0(1,isppol) 424 omega = omega_min 425 do iomega=1,elph_ds%na2f 426 ! bxu, at which eps and eps' should I save it 427 ! better to save them all, but could be too many 428 write (unit_a2f_tr, '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr(iomega,:,isppol,1,INT(elph_ds%n_pair/2)+1,1) 429 omega=omega+domega 430 end do 431 write (unit_a2f_tr,*) 432 end do !isppol 433 434 close (unit=unit_a2f_tr) 435 436 !calculation of transport properties 437 ABI_MALLOC(integrho,(elph_ds%na2f)) 438 ABI_MALLOC(tointegrho,(elph_ds%na2f)) 439 ABI_MALLOC(integrand_q00,(elph_ds%na2f)) 440 ABI_MALLOC(integrand_q01,(elph_ds%na2f)) 441 ABI_MALLOC(integrand_q11,(elph_ds%na2f)) 442 ABI_MALLOC(q00,(ntemper,3,3,elph_ds%nsppol)) 443 ABI_MALLOC(q01,(ntemper,3,3,elph_ds%nsppol)) 444 ABI_MALLOC(q11,(ntemper,3,3,elph_ds%nsppol)) 445 ABI_MALLOC(seebeck,(elph_ds%nsppol,ntemper,3,3)) 446 !ABI_MALLOC(rho_nm,(elph_ds%nsppol,ntemper,3,3)) 447 448 fname = trim(elph_ds%elph_base_name) // '_RHO' 449 if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then 450 ABI_ERROR(message) 451 end if 452 !print header to resistivity file 453 write (unit_rho,*) '# Resistivity as a function of temperature.' 454 write (unit_rho,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 455 write (unit_rho,*) '# ' 456 write (unit_rho,*) '# Columns are: ' 457 write (unit_rho,*) '# temperature[K] rho[au] rho [SI] rho/temp [au]' 458 write (unit_rho,*) '# ' 459 460 fname = trim(elph_ds%elph_base_name) // '_TAU' 461 if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then 462 ABI_ERROR(message) 463 end if 464 !print header to relaxation time file 465 write (unit_tau,*) '# Relaxation time as a function of temperature.' 466 write (unit_tau,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 467 write (unit_tau,*) '# ' 468 write (unit_tau,*) '# Columns are: ' 469 write (unit_tau,*) '# temperature[K] tau[au] tau [SI] ' 470 write (unit_tau,*) '# ' 471 472 fname = trim(elph_ds%elph_base_name) // '_SBK' 473 if (open_file(fname,message,newunit=unit_sbk,status='unknown') /= 0) then 474 ABI_ERROR(message) 475 end if 476 !print header to relaxation time file 477 write (unit_sbk,*) '# Seebeck Coefficint as a function of temperature.' 478 write (unit_sbk,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 479 write (unit_sbk,*) '# ' 480 write (unit_sbk,*) '# Columns are: ' 481 write (unit_sbk,*) '# temperature[K] S [au] S [SI] ' 482 write (unit_sbk,*) '# ' 483 484 fname = trim(elph_ds%elph_base_name) // '_WTH' 485 if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then 486 ABI_ERROR(message) 487 end if 488 489 !print header to thermal conductivity file 490 write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.' 491 write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong' 492 write (unit_therm,'(a)') '# ' 493 write (unit_therm,'(a)') '# Columns are: ' 494 write (unit_therm,'(a)') '# temperature[K] thermal rho[au] thermal cond [au] thermal rho [SI] thermal cond [SI]' 495 write (unit_therm,'(a)') '# ' 496 497 fname = trim(elph_ds%elph_base_name) // '_LOR' 498 if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then 499 ABI_ERROR(message) 500 end if 501 502 !print header to lorentz file 503 write (unit_lor,*) '# Lorentz number as a function of temperature.' 504 write (unit_lor,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 505 write (unit_lor,*) '# ' 506 write (unit_lor,*) '# Columns are: ' 507 write (unit_lor,*) '# temperature[K] Lorentz number[au] Lorentz quantum = (pi*kb_HaK)**2/3' 508 write (unit_lor,*) '# ' 509 510 do isppol=1,elph_ds%nsppol 511 lambda_tr_trace = zero 512 do itrtensor=1,9 513 omega = omega_min 514 tointegrho = zero 515 do iomega=1,elph_ds%na2f 516 if(omega<=0) then 517 omega=omega+domega 518 cycle 519 end if 520 ! bxu, agian, which eps and eps' to use? 521 ! sometimes Ef is in the gap 522 tointegrho(iomega)=two*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega 523 omega=omega+domega 524 end do 525 526 integrho = zero 527 call simpson_int(elph_ds%na2f,domega,tointegrho,integrho) 528 lambda_tr=integrho(elph_ds%na2f) 529 write (message, '(a,2i3,a,es16.6)' )& 530 & ' mka2f_tr: TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' = ', lambda_tr 531 call wrtout(std_out,message,'COLL') 532 if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr 533 end do !end itrtensor do 534 535 lambda_tr_trace = lambda_tr_trace / three 536 write (message, '(a,i3,a,es16.6)' )& 537 & ' mka2f_tr: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' = ', lambda_tr_trace 538 call wrtout(std_out,message,'COLL') 539 call wrtout(ab_out,message,'COLL') 540 end do !end isppol do 541 542 !constant to change units of rho from au to SI 543 chgu=2.173969*1.0d-7 544 chtu=2.4188843265*1.0d-17 545 chsu=8.617343101*1.0d-5 ! au to J/C/K 546 chwu=9.270955772*1.0d-5 ! au to mK/W 547 548 !change the fermi level to zero, as required for q01 to vanish. 549 tmp_fermie = elph_ds%fermie 550 !Get Q00, Q01, Q11, and derive rho, tau 551 q00 = zero 552 q01 = zero 553 q11 = zero 554 ! prepare s1 and s2 arrays 555 s1 = (/1, 1, -1, -1/) 556 s2 = (/1, -1, 1, -1/) 557 558 do isppol=1,elph_ds%nsppol 559 do icomp=1, 3 560 do jcomp=1, 3 561 itrtensor=(icomp-1)*3+jcomp 562 563 write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor 564 write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor 565 566 do itemp=1,ntemper ! runs over termperature in K 567 Temp=tempermin+temperinc*dble(itemp) 568 tmp_veloc_sq0 = sqrt(elph_tr_ds%veloc_sq0(itemp,icomp,isppol)*elph_tr_ds%veloc_sq0(itemp,jcomp,isppol)) 569 570 integrand_q00 = zero 571 integrand_q01 = zero 572 integrand_q11 = zero 573 574 omega = omega_min 575 do iomega=1,elph_ds%na2f 576 if(omega .le. 0) then 577 omega=omega+domega 578 cycle 579 end if 580 xtr=omega/(kb_HaK*Temp) 581 occ_omega=1.0_dp/(exp(xtr)-1.0_dp) 582 583 tmp_veloc_sq1 = zero 584 tmp_veloc_sq2 = zero 585 do ie1=1,elph_ds%nenergy 586 e1 = elph_tr_ds%en_all(isppol,ie1) 587 588 !! BXU, the tolerance here needs to be used with caution 589 !! which depends on the dimensions of the system 590 !! E.g. in 2D system, DOS can be much smaller 591 if (elph_tr_ds%dos_n(ie1,isppol)/natom .lt. 0.1d0) cycle ! energy in the gap forbidden 592 593 xtr=(e1-tmp_fermie)/(kb_HaK*Temp) 594 occ_e1=1.0_dp/(exp(xtr)+1.0_dp) 595 596 e2 = e1 + omega 597 xtr=(e2-tmp_fermie)/(kb_HaK*Temp) 598 occ_e2=1.0_dp/(exp(xtr)+1.0_dp) 599 ! Do we need to change the fermie to the one with T dependence? 600 ! find which ie2 give the closest energy 601 if (e2 .gt. elph_tr_ds%en_all(isppol,elph_ds%nenergy)) then 602 ie2 = 0 603 cycle 604 else 605 ie_tmp = 1 606 diff = dabs(e2-elph_tr_ds%en_all(isppol,1)) 607 do ie2 = 2, elph_ds%nenergy 608 if (dabs(e2-elph_tr_ds%en_all(isppol,ie2)) .lt. diff) then 609 diff = dabs(e2-elph_tr_ds%en_all(isppol,ie2)) 610 ie_tmp = ie2 611 end if 612 end do 613 ie2 = ie_tmp 614 615 if (e2 < elph_tr_ds%en_all(isppol,ie2)) then 616 ie2_right = ie2 617 ie2_left = ie2-1 618 else 619 ie2_right = ie2+1 620 ie2_left = ie2 621 end if 622 623 end if 624 625 if (elph_tr_ds%dos_n(ie2,isppol)/natom .lt. 0.1d0) cycle 626 627 tointegq00 = zero 628 tointegq01 = zero 629 tointegq11 = zero 630 631 ! BXU linear interpolation of volec_sq and dos_n 632 xe=(e2-elph_tr_ds%en_all(isppol,ie2_left))/ & 633 & (elph_tr_ds%en_all(isppol,ie2_right)-elph_tr_ds%en_all(isppol,ie2_left)) 634 veloc_sq_icomp = elph_tr_ds%veloc_sq(icomp,isppol,ie2_left)*(1.0d0-xe) + & 635 & elph_tr_ds%veloc_sq(icomp,isppol,ie2_right)*xe 636 veloc_sq_jcomp = elph_tr_ds%veloc_sq(jcomp,isppol,ie2_left)*(1.0d0-xe) + & 637 & elph_tr_ds%veloc_sq(jcomp,isppol,ie2_right)*xe 638 dos_n_e2 = elph_tr_ds%dos_n(ie2_left,isppol)*(1.0d0-xe) + & 639 & elph_tr_ds%dos_n(ie2_right,isppol)*xe 640 641 tmp_veloc_sq1 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie1)*elph_tr_ds%veloc_sq(jcomp,isppol,ie1)) 642 ! tmp_veloc_sq2 = sqrt(elph_tr_ds%veloc_sq(icomp,isppol,ie2)*elph_tr_ds%veloc_sq(jcomp,isppol,ie2)) 643 tmp_veloc_sq2 = sqrt(veloc_sq_icomp*veloc_sq_jcomp) 644 645 ! ie_1 = (ie1-1)*elph_ds%nenergy + ie2 646 ! ie_2 = (ie2-1)*elph_ds%nenergy + ie1 647 ie_1 = pair2red(ie1,ie2) 648 ie_2 = pair2red(ie2,ie1) 649 if (ie_1 .eq. 0 .or. ie_2 .eq. 0) then 650 ABI_BUG('CHECK pair2red!') 651 end if 652 653 do ssp=1,4 ! (s,s'=+/-1, condense the indices) 654 655 nv1 = 1.0_dp/(elph_tr_ds%dos_n(ie1,isppol)*sqrt(tmp_veloc_sq1)) 656 sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK) 657 !DEBUG 658 if (elph_ds%ep_lova .eq. 1) then 659 nv1 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0)) 660 sigma1 = sqrt(3.0_dp)*(e1-tmp_fermie)/(pi*Temp*kb_HaK) 661 end if 662 !ENDDEBUG 663 664 tointegq00_1 = zero 665 tointegq01_1 = zero 666 tointegq11_1 = zero 667 668 !DEBUG 669 if (elph_ds%ep_lova .eq. 1) then 670 nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0)) 671 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 672 j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp 673 j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 674 j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 675 tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 676 & occ_e1*(1.0_dp-occ_e2)*j00*occ_omega 677 tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 678 & occ_e1*(1.0_dp-occ_e2)*j01*occ_omega 679 tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 680 & occ_e1*(1.0_dp-occ_e2)*j11*occ_omega 681 !END DEBUG 682 else if (elph_ds%ep_lova .eq. 0) then 683 nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2)) 684 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 685 j00 = (nv1 + s1(ssp)*nv2)*(nv1 + s2(ssp)*nv2)/4.0_dp 686 j01 = (nv1 + s1(ssp)*nv2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 687 j11 = (nv1*sigma1 + s1(ssp)*nv2*sigma2)*(nv1*sigma1 + s2(ssp)*nv2*sigma2)/4.0_dp 688 ! bxu TEST 689 if (debug) then 690 if (ssp .eq. 1 .and. itrtensor .eq. 1) then 691 write(21,"(3i5,4E20.12)") iomega, ie1, ie2, & 692 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 693 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 694 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 695 end if 696 if (ssp .eq. 2 .and. itrtensor .eq. 1) then 697 write(22,"(3i5,4E20.12)") iomega, ie1, ie2, & 698 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 699 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 700 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 701 end if 702 if (ssp .eq. 3 .and. itrtensor .eq. 1) then 703 write(23,"(3i5,4E20.12)") iomega, ie1, ie2, & 704 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 705 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 706 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 707 end if 708 if (ssp .eq. 4 .and. itrtensor .eq. 1) then 709 write(24,"(3i5,4E20.12)") iomega, ie1, ie2, & 710 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp), j01, & 711 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01, & 712 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)*j01*occ_e1*(1.0_dp-occ_e2)*occ_omega 713 end if 714 end if 715 tointegq00_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* & 716 & occ_e1*(1.0_dp-occ_e2)*j00*occ_omega 717 tointegq01_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* & 718 & occ_e1*(1.0_dp-occ_e2)*j01*occ_omega 719 tointegq11_1 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_1,itemp)* & 720 & occ_e1*(1.0_dp-occ_e2)*j11*occ_omega 721 end if 722 723 tointegq00_2 = zero 724 tointegq01_2 = zero 725 tointegq11_2 = zero 726 727 !DEBUG 728 if (elph_ds%ep_lova .eq. 1) then 729 nv2 = 1.0_dp/(elph_tr_ds%dos_n0(itemp,isppol)*sqrt(tmp_veloc_sq0)) 730 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 731 j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp 732 j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 733 j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 734 tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 735 & occ_e1*(1.0_dp-occ_e2)*j00*(occ_omega+1) 736 tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 737 & occ_e1*(1.0_dp-occ_e2)*j01*(occ_omega+1) 738 tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,1,itemp)* & 739 & occ_e1*(1.0_dp-occ_e2)*j11*(occ_omega+1) 740 !END DEBUG 741 else if (elph_ds%ep_lova .eq. 0) then 742 nv2 = 1.0_dp/(dos_n_e2*sqrt(tmp_veloc_sq2)) 743 sigma2 = sqrt(3.0_dp)*(e2-tmp_fermie)/(pi*Temp*kb_HaK) 744 j00 = (nv2 + s1(ssp)*nv1)*(nv2 + s2(ssp)*nv1)/4.0_dp 745 j01 = (nv2 + s1(ssp)*nv1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 746 j11 = (nv2*sigma2 + s1(ssp)*nv1*sigma1)*(nv2*sigma2 + s2(ssp)*nv1*sigma1)/4.0_dp 747 !DEBUG bxu TEST 748 if (debug) then 749 if (ssp .eq. 1 .and. itrtensor .eq. 1) then 750 write(31,"(3i5,4E20.12)") iomega, ie2, ie1, & 751 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 752 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 753 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 754 end if 755 if (ssp .eq. 2 .and. itrtensor .eq. 1) then 756 write(32,"(3i5,4E20.12)") iomega, ie2, ie1, & 757 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 758 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 759 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 760 end if 761 if (ssp .eq. 3 .and. itrtensor .eq. 1) then 762 write(33,"(3i5,4E20.12)") iomega, ie2, ie1, & 763 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 764 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 765 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 766 end if 767 if (ssp .eq. 4 .and. itrtensor .eq. 1) then 768 write(34,"(3i5,4E20.12)") iomega, ie2, ie1, & 769 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp), j01, & 770 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01, & 771 & elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)*j01*occ_e2*(1.0_dp-occ_e1)*(occ_omega+1) 772 end if 773 end if 774 !ENDDEBUG 775 tointegq00_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* & 776 & occ_e2*(1.0_dp-occ_e1)*j00*(occ_omega+1) 777 tointegq01_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* & 778 & occ_e2*(1.0_dp-occ_e1)*j01*(occ_omega+1) 779 tointegq11_2 = elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,ssp,ie_2,itemp)* & 780 & occ_e2*(1.0_dp-occ_e1)*j11*(occ_omega+1) 781 end if ! elph_ds%ep_lova 782 783 tointegq00 = tointegq00 + tointegq00_1 + tointegq00_2 784 tointegq01 = tointegq01 + tointegq01_1 + tointegq01_2 785 tointegq11 = tointegq11 + tointegq11_1 + tointegq11_2 786 787 end do ! ss' = 4 788 integrand_q00(iomega) = integrand_q00(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq00 789 integrand_q01(iomega) = integrand_q01(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq01 790 integrand_q11(iomega) = integrand_q11(iomega) + elph_tr_ds%de_all(isppol,ie1)*tointegq11 791 end do ! ie1 ~ 20 792 omega=omega+domega 793 q00(itemp,icomp,jcomp,isppol) = q00(itemp,icomp,jcomp,isppol) +& 794 & domega*integrand_q00(iomega) 795 q01(itemp,icomp,jcomp,isppol) = q01(itemp,icomp,jcomp,isppol) +& 796 & domega*integrand_q01(iomega) 797 q11(itemp,icomp,jcomp,isppol) = q11(itemp,icomp,jcomp,isppol) +& 798 & domega*integrand_q11(iomega) 799 end do ! omega ~ 400 800 801 q00(itemp,icomp,jcomp,isppol)=q00(itemp,icomp,jcomp,isppol)* & 802 & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp) 803 q01(itemp,icomp,jcomp,isppol)=q01(itemp,icomp,jcomp,isppol)* & 804 & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp) 805 q11(itemp,icomp,jcomp,isppol)=q11(itemp,icomp,jcomp,isppol)* & 806 & ucvol*two_pi*elph_tr_ds%dos_n0(itemp,isppol)/(kb_HaK*Temp) 807 808 rho = 0.5_dp*q00(itemp,icomp,jcomp,isppol) 809 ! is tau energy dependent? 810 tau = 2.0d0*ucvol/(q00(itemp,icomp,jcomp,isppol)*elph_tr_ds%dos_n0(itemp,isppol)*tmp_veloc_sq0) 811 write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp 812 write(unit_tau,'(3D20.10)')temp,tau,tau*chtu 813 rho_T(itemp)=rho 814 end do ! temperature = 1? 815 write(unit_rho,*) 816 write(unit_tau,*) 817 818 end do ! jcomp = 3 819 end do ! icomp = 3 820 end do ! isppol = 2 821 822 !----------------------------- 823 824 seebeck = zero 825 !rho_nm = zero 826 827 !do isppol=1,elph_ds%nsppol 828 !do itemp=1,ntemper 829 !q11_inv(:,:)=q11(itemp,:,:,isppol) 830 !call matrginv(q11_inv,3,3) 831 !do icomp=1,3 832 !do jcomp=1,3 833 !do kcomp=1,3 834 !seebeck(isppol,itemp,icomp,jcomp) = seebeck(isppol,itemp,icomp,jcomp) + & 835 !& q01(itemp,icomp,kcomp,isppol)*q11_inv(kcomp,jcomp) 836 !end do 837 !end do 838 !end do 839 !end do 840 !end do 841 842 do isppol=1,elph_ds%nsppol 843 do itemp=1,ntemper 844 q11_inv(:,:)=q11(itemp,:,:,isppol) 845 call matrginv(q11_inv,3,3) 846 call DGEMM('N','N',3,3,3,one,q01(itemp,:,:,isppol),3,q11_inv,& 847 & 3,zero,seebeck(isppol,itemp,:,:),3) 848 ! call DGEMM('N','N',3,3,3,one,seebeck(isppol,itemp,:,:),3,q01(itemp,:,:,isppol),& 849 ! & 3,zero,rho_nm(isppol,itemp,:,:),3) 850 end do 851 end do 852 pref_s = pi/sqrt(3.0_dp) 853 seebeck=pref_s*seebeck 854 855 !fullq = zero 856 !do icomp=1,3 857 !do jcomp=1,3 858 !fullq(icomp,jcomp) = q00(1,icomp,jcomp,1) 859 !end do 860 !end do 861 !do icomp=1,3 862 !do jcomp=4,6 863 !fullq(icomp,jcomp) = q01(1,icomp,jcomp-3,1) 864 !end do 865 !end do 866 !do icomp=4,6 867 !do jcomp=1,3 868 !fullq(icomp,jcomp) = q01(1,icomp-3,jcomp,1) 869 !end do 870 !end do 871 !do icomp=4,6 872 !do jcomp=4,6 873 !fullq(icomp,jcomp) = q11(1,icomp-3,jcomp-3,1) 874 !end do 875 !end do 876 !write(*,*)' fullq' 877 !write(*,"(6E20.12)") (fullq(1,jcomp),jcomp=1,6) 878 !write(*,"(6E20.12)") (fullq(2,jcomp),jcomp=1,6) 879 !write(*,"(6E20.12)") (fullq(3,jcomp),jcomp=1,6) 880 !write(*,"(6E20.12)") (fullq(4,jcomp),jcomp=1,6) 881 !write(*,"(6E20.12)") (fullq(5,jcomp),jcomp=1,6) 882 !write(*,"(6E20.12)") (fullq(6,jcomp),jcomp=1,6) 883 write(message,'(a)') 'q00:' 884 call wrtout(std_out,message,'COLL') 885 write(message,'(3E20.12)') (q00(1,1,jcomp,1),jcomp=1,3) 886 call wrtout(std_out,message,'COLL') 887 write(message,'(3E20.12)') (q00(1,2,jcomp,1),jcomp=1,3) 888 call wrtout(std_out,message,'COLL') 889 write(message,'(3E20.12)') (q00(1,3,jcomp,1),jcomp=1,3) 890 call wrtout(std_out,message,'COLL') 891 write(message,'(a)') 'q01:' 892 call wrtout(std_out,message,'COLL') 893 write(message,'(3E20.12)') (q01(1,1,jcomp,1),jcomp=1,3) 894 call wrtout(std_out,message,'COLL') 895 write(message,'(3E20.12)') (q01(1,2,jcomp,1),jcomp=1,3) 896 call wrtout(std_out,message,'COLL') 897 write(message,'(3E20.12)') (q01(1,3,jcomp,1),jcomp=1,3) 898 call wrtout(std_out,message,'COLL') 899 write(message,'(a)') 'q11, q11_inv:' 900 call wrtout(std_out,message,'COLL') 901 write(message,'(6E20.12)') (q11(1,1,jcomp,1),jcomp=1,3),(q11_inv(1,jcomp),jcomp=1,3) 902 call wrtout(std_out,message,'COLL') 903 write(message,'(6E20.12)') (q11(1,2,jcomp,1),jcomp=1,3),(q11_inv(2,jcomp),jcomp=1,3) 904 call wrtout(std_out,message,'COLL') 905 write(message,'(6E20.12)') (q11(1,3,jcomp,1),jcomp=1,3),(q11_inv(3,jcomp),jcomp=1,3) 906 call wrtout(std_out,message,'COLL') 907 !q11_inv = zero 908 !do icomp = 1, 3 909 !q11_inv(icomp,icomp) = 2.0_dp 910 !end do 911 912 !call matrginv(fullq,6,6) 913 914 !do isppol=1,elph_ds%nsppol 915 !do itemp=1,ntemper 916 !rho_nm(isppol,itemp,:,:) = q00(itemp,:,:,isppol) - rho_nm(isppol,itemp,:,:) 917 !end do 918 !end do 919 !rho_nm = 0.5_dp*rho_nm 920 921 !Output of Seebeck coefficient 922 do isppol=1,elph_ds%nsppol 923 do icomp=1,3 924 do jcomp=1,3 925 itrtensor=(icomp-1)*3+jcomp 926 write(unit_sbk,*) '# Seebeck for isppol, itrten = ', isppol, itrtensor 927 ! write(88,*) '# Rho for isppol, itrten = ', isppol, itrtensor 928 ! write(89,*) '# Rho for isppol, itrten = ', isppol, itrtensor 929 do itemp=1,ntemper 930 Temp=tempermin+temperinc*dble(itemp) 931 write(unit_sbk,'(3D20.10)')temp, seebeck(isppol,itemp,icomp,jcomp), seebeck(isppol,itemp,icomp,jcomp)*chsu 932 ! write(88,'(3D20.10)')temp, rho_nm(isppol,itemp,icomp,jcomp), rho_nm(isppol,itemp,icomp,jcomp)*chgu 933 ! write(89,'(3D20.10)')temp, 0.5_dp/fullq(1,1), 0.5_dp*chgu/fullq(1,1) 934 end do ! temperature 935 write(unit_sbk,*) 936 ! write(88,*) 937 ! write(89,*) 938 end do ! jcomp 939 end do ! icomp 940 end do ! isppol 941 942 !Get thermal resistivity, based on eqn. (52) in Allen's PRB 17, 3725 (1978) [[cite:Allen1978]] 943 !WARNING: before 6.13.1 the thermal resistivity and Lorentz number were not in 944 !atomic units, BUT the SI units are good. 945 pref_w = 3.0_dp/(2.0_dp*pi**2.0d0) 946 do isppol=1,elph_ds%nsppol 947 do icomp=1, 3 948 do jcomp=1, 3 949 itrtensor=(icomp-1)*3+jcomp 950 951 write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol 952 write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol 953 954 do itemp=1,ntemper 955 956 Temp=tempermin + temperinc*dble(itemp) 957 958 wtherm = pref_w*q11(itemp,icomp,jcomp,isppol)/(kb_HaK*Temp) 959 960 ! write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9 961 write(unit_therm,'(5D20.10)')temp,wtherm,1.0_dp/wtherm,wtherm*chwu,1.0_dp/(wtherm*chwu) 962 963 lorentz=rho_T(itemp)/(wtherm*kb_HaK*Temp) 964 write(unit_lor,*)temp,lorentz,lor0 965 966 end do 967 write(unit_therm,*) 968 write(unit_lor,*) 969 end do ! jcomp 970 end do ! icomp 971 end do ! isppol 972 973 974 ABI_FREE(phfrq) 975 ABI_FREE(displ) 976 ABI_FREE(pheigvec) 977 ABI_FREE(integrand_q00) 978 ABI_FREE(integrand_q01) 979 ABI_FREE(integrand_q11) 980 ABI_FREE(q00) 981 ABI_FREE(q01) 982 ABI_FREE(q11) 983 ABI_FREE(seebeck) 984 ABI_FREE(rho_T) 985 ABI_FREE(integrho) 986 ABI_FREE(tointegrho) 987 988 close (unit=unit_lor) 989 close (unit=unit_rho) 990 close (unit=unit_tau) 991 close (unit=unit_sbk) 992 close (unit=unit_therm) 993 994 ABI_FREE(elph_ds%k_fine%wtq) 995 ABI_FREE(elph_ds%k_phon%wtq) 996 997 ABI_FREE(elph_tr_ds%a2f_1d_tr) 998 999 ABI_FREE(elph_tr_ds%gamma_qpt_tr) 1000 ABI_FREE(elph_tr_ds%gamma_rpt_tr) 1001 write(std_out,*) ' mka2f_tr : end ' 1002 1003 end subroutine mka2f_tr
ABINIT/mka2f_tr_lova [ 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
NOTES
copied from ftiaf9.f
SOURCE
1042 subroutine mka2f_tr_lova(crystal,ifc,elph_ds,ntemper,tempermin,temperinc,elph_tr_ds) 1043 1044 !Arguments ------------------------------------ 1045 !scalars 1046 integer,intent(in) :: ntemper 1047 real(dp),intent(in) :: tempermin,temperinc 1048 type(crystal_t),intent(in) :: crystal 1049 type(ifc_type),intent(in) :: ifc 1050 type(elph_tr_type),intent(inout) :: elph_tr_ds 1051 type(elph_type),intent(inout) :: elph_ds 1052 1053 !Local variables ------------------------- 1054 !x =w/(2kbT) 1055 !scalars 1056 integer :: iFSqpt,ibranch,iomega,isppol,jbranch,nerr 1057 integer :: unit_a2f_tr, unit_a2f_trout, unit_a2f_trin, natom 1058 integer :: idir, iatom, k1, kdir,unit_lor,unit_rho,unit_tau,unit_therm 1059 integer :: itemp,nrpt,itrtensor, icomp, jcomp 1060 real(dp) :: Temp,chgu,chtu,femto,diagerr,firh,firhT,gaussfactor,domega 1061 real(dp) :: firh_tau,firhT_tau ! added by BX to get Tau 1062 real(dp) :: a2fprefactor_in, temp_in 1063 real(dp) :: a2fprefactor_out, temp_out 1064 real(dp) :: gaussprefactor,gaussval,lambda_tr,lor0,lorentz,maxerr,maxx,omega 1065 real(dp) :: rho,tau,tolexp,wtherm,xtr,xx 1066 real(dp) :: lambda_tr_trace,omega_min, omega_max,qnorm2,spinfact 1067 character(len=500) :: message 1068 character(len=fnlen) :: fname 1069 !arrays 1070 complex(dpc),parameter :: c0=dcmplx(0.d0,0.d0),c1=dcmplx(1.d0,0.d0) 1071 real(dp) :: gprimd(3,3) 1072 real(dp) :: eigval_in(elph_ds%nbranch) 1073 real(dp) :: eigval_out(elph_ds%nbranch) 1074 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 1075 real(dp) :: gam_now_in (2,elph_ds%nbranch*elph_ds%nbranch) 1076 real(dp) :: gam_now_out(2,elph_ds%nbranch*elph_ds%nbranch) 1077 real(dp) :: tmpa2f_in (elph_ds%na2f) 1078 real(dp) :: tmpa2f_out(elph_ds%na2f) 1079 real(dp) :: tmpgam1(2,elph_ds%nbranch,elph_ds%nbranch) 1080 real(dp) :: tmpgam2(2,elph_ds%nbranch,elph_ds%nbranch) 1081 real(dp),allocatable :: phfrq(:,:) 1082 real(dp),allocatable :: displ(:,:,:,:) 1083 real(dp),allocatable :: pheigvec(:,:) 1084 real(dp),allocatable :: integrho(:),integtau(:),tointegrho(:),tointega2f(:),tointegtau(:) 1085 real(dp),allocatable :: rho_T(:),tau_T(:) 1086 real(dp),allocatable :: coskr(:,:) 1087 real(dp),allocatable :: sinkr(:,:) 1088 1089 ! ********************************************************************* 1090 1091 !calculate a2f_tr for frequencies between 0 and omega_max 1092 write(std_out,*) 'mka2f_tr_lova : enter ' 1093 ! 1094 !MG: the step should be calculated locally using nomega and the extrema of the spectrum. 1095 !One should not rely on previous calls for the setup of elph_ds%domega 1096 !I will remove elph_ds%domega since mka2f.F90 will become a method of gamma_t 1097 domega =elph_ds%domega 1098 1099 ! Number of points for FFT interpolation 1100 nrpt = ifc%nrpt 1101 natom = crystal%natom 1102 gprimd = crystal%gprimd 1103 1104 ABI_MALLOC(elph_tr_ds%a2f_1d_tr,(elph_ds%na2f,9,elph_ds%nsppol,1,1,1)) 1105 ABI_MALLOC(elph_tr_ds%a2f_1d_trin,(elph_ds%na2f,9,elph_ds%nsppol)) 1106 ABI_MALLOC(elph_tr_ds%a2f_1d_trout,(elph_ds%na2f,9,elph_ds%nsppol)) 1107 1108 !! defaults for number of temperature steps and max T (all in Kelvin...) 1109 !ntemper=1000 1110 !tempermin=zero 1111 !temperinc=one 1112 ABI_MALLOC(rho_T,(ntemper)) 1113 ABI_MALLOC(tau_T,(ntemper)) 1114 1115 1116 !tolerance on gaussian being = 0 1117 tolexp = 1.d-100 1118 maxx = sqrt(-log(tolexp)) 1119 lor0=(pi*kb_HaK)**2/3. 1120 1121 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F) 1122 !WARNING! supposes this value has been set in mkelph_linwid. 1123 1124 gaussprefactor = sqrt(piinv) / elph_ds%a2fsmear 1125 gaussfactor = one / elph_ds%a2fsmear 1126 1127 !spinfact should be 1 for a normal non sppol calculation without spinorbit 1128 !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 1129 !for sppol 2 it should be 0.5 as we have 2 spin channels to sum 1130 spinfact = one / elph_ds%nsppol !/ elph_ds%nspinor 1131 1132 !ENDMG 1133 1134 elph_tr_ds%a2f_1d_tr = zero 1135 elph_tr_ds%a2f_1d_trin = zero 1136 elph_tr_ds%a2f_1d_trout = zero 1137 1138 maxerr=0. 1139 nerr=0 1140 1141 ABI_MALLOC(phfrq,(elph_ds%nbranch, elph_ds%k_fine%nkpt)) 1142 ABI_MALLOC(displ,(2, elph_ds%nbranch, elph_ds%nbranch, elph_ds%k_fine%nkpt)) 1143 ABI_MALLOC(pheigvec,(2*elph_ds%nbranch*elph_ds%nbranch, elph_ds%k_fine%nkpt)) 1144 1145 do iFSqpt=1,elph_ds%k_fine%nkpt 1146 call ifc%fourq(crystal,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ(:,:,:,iFSqpt),out_eigvec=pheigvec(:,iFSqpt)) 1147 end do 1148 1149 omega_min = minval(phfrq) 1150 omega_max = maxval(phfrq) 1151 1152 ABI_MALLOC(coskr, (elph_ds%k_fine%nkpt,nrpt)) 1153 ABI_MALLOC(sinkr, (elph_ds%k_fine%nkpt,nrpt)) 1154 call ftgam_init(Ifc%gprim, elph_ds%k_fine%nkpt, nrpt, elph_ds%k_fine%kpt, Ifc%rpt, coskr, sinkr) 1155 1156 do isppol=1,elph_ds%nsppol 1157 1158 ! loop over qpoint in full kpt grid (presumably dense) 1159 do iFSqpt=1,elph_ds%k_fine%nkpt 1160 qnorm2 = sum(elph_ds%k_fine%kpt(:,iFSqpt)**2) 1161 ! if (flag_to_exclude_soft_modes = .false.) qnorm2 = zero 1162 do itrtensor=1,9 1163 ! Do FT from real-space gamma grid to 1 qpt. 1164 1165 if (elph_ds%ep_int_gkk == 1) then 1166 gam_now_in(:,:) = elph_tr_ds%gamma_qpt_trin(:,itrtensor,:,isppol,iFSqpt) 1167 gam_now_out(:,:) = elph_tr_ds%gamma_qpt_trout(:,itrtensor,:,isppol,iFSqpt) 1168 else 1169 call ftgam(Ifc%wghatm,gam_now_in, elph_tr_ds%gamma_rpt_trin(:,itrtensor,:,isppol,:),natom,1,nrpt,0,& 1170 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 1171 call ftgam(Ifc%wghatm,gam_now_out,elph_tr_ds%gamma_rpt_trout(:,itrtensor,:,isppol,:),natom,1,nrpt,0,& 1172 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 1173 end if 1174 1175 ! Diagonalize gamma matrix at this qpoint (complex matrix). 1176 1177 ! if ep_scalprod==0 we have to dot in the displacement vectors here 1178 if (elph_ds%ep_scalprod==0) then 1179 1180 displ_red(:,:,:) = zero 1181 do jbranch=1,elph_ds%nbranch 1182 do iatom=1,natom 1183 do idir=1,3 1184 ibranch=idir+3*(iatom-1) 1185 do kdir=1,3 1186 k1 = kdir+3*(iatom-1) 1187 displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + & 1188 & gprimd(kdir,idir)*displ(1,k1,jbranch,iFSqpt) 1189 displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + & 1190 & gprimd(kdir,idir)*displ(2,k1,jbranch,iFSqpt) 1191 end do 1192 end do 1193 end do 1194 end do 1195 1196 tmpgam2 = reshape (gam_now_in, (/2,elph_ds%nbranch,elph_ds%nbranch/)) 1197 call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1) 1198 do jbranch=1,elph_ds%nbranch 1199 eigval_in(jbranch) = tmpgam1(1, jbranch, jbranch) 1200 end do 1201 1202 tmpgam2 = reshape (gam_now_out, (/2,elph_ds%nbranch,elph_ds%nbranch/)) 1203 call gam_mult_displ(elph_ds%nbranch, displ_red, tmpgam2, tmpgam1) 1204 do jbranch=1,elph_ds%nbranch 1205 eigval_out(jbranch) = tmpgam1(1, jbranch, jbranch) 1206 end do 1207 1208 else if (elph_ds%ep_scalprod == 1) then 1209 1210 ! 1211 ! NOTE: in these calls gam_now and pheigvec do not have the right rank, but blas usually does not care 1212 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_in, 3*natom,& 1213 & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom) 1214 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,& 1215 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 1216 diagerr = zero 1217 1218 do ibranch=1,elph_ds%nbranch 1219 eigval_in(ibranch) = tmpgam2(1,ibranch,ibranch) 1220 do jbranch=1,ibranch-1 1221 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 1222 end do 1223 do jbranch=ibranch+1,elph_ds%nbranch 1224 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 1225 end do 1226 end do 1227 if (diagerr > tol12) then 1228 nerr=nerr+1 1229 maxerr=max(diagerr, maxerr) 1230 end if 1231 1232 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now_out, 3*natom,& 1233 & pheigvec(:,iFSqpt), 3*natom, c0, tmpgam1, 3*natom) 1234 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec(:,iFSqpt), 3*natom,& 1235 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 1236 diagerr = zero 1237 1238 do ibranch=1,elph_ds%nbranch 1239 eigval_out(ibranch) = tmpgam2(1,ibranch,ibranch) 1240 do jbranch=1,ibranch-1 1241 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 1242 end do 1243 do jbranch=ibranch+1,elph_ds%nbranch 1244 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch)) 1245 end do 1246 end do 1247 if (diagerr > tol12) then 1248 nerr=nerr+1 1249 maxerr=max(diagerr, maxerr) 1250 end if 1251 end if 1252 ! end ep_scalprod if 1253 1254 ! Add all contributions from the phonon modes at this qpoint to 1255 ! a2f and the phonon dos. 1256 do ibranch=1,elph_ds%nbranch 1257 ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then 1258 if ( abs(phfrq(ibranch,iFSqpt)) < tol7 .or. & 1259 & (phfrq(ibranch,iFSqpt) < tol4 .and. qnorm2 > 0.03 )) then ! 1260 ! note: this should depend on the velocity of sound, to accept acoustic 1261 ! modes! 1262 a2fprefactor_in = zero 1263 a2fprefactor_out= zero 1264 else 1265 a2fprefactor_in = eigval_in (ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol)) 1266 a2fprefactor_out = eigval_out(ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*elph_ds%n0(isppol)) 1267 end if 1268 1269 omega = omega_min 1270 tmpa2f_in (:) = zero 1271 tmpa2f_out(:) = zero 1272 do iomega=1,elph_ds%na2f 1273 xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor 1274 gaussval = gaussprefactor*exp(-xx*xx) 1275 1276 temp_in = gaussval*a2fprefactor_in 1277 temp_out = gaussval*a2fprefactor_out 1278 1279 if (dabs(temp_in) < 1.0d-50) temp_in = zero 1280 if (dabs(temp_out) < 1.0d-50) temp_out = zero 1281 tmpa2f_in (iomega) = tmpa2f_in (iomega) + temp_in 1282 tmpa2f_out(iomega) = tmpa2f_out(iomega) + temp_out 1283 omega = omega+domega 1284 end do 1285 1286 elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trin (:,itrtensor,isppol) + tmpa2f_in(:) 1287 elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) = elph_tr_ds%a2f_1d_trout(:,itrtensor,isppol) + tmpa2f_out(:) 1288 1289 end do ! end ibranch do 1290 end do ! end itrtensor do 1291 end do ! end iFSqpt do 1292 end do ! end isppol 1293 1294 ABI_FREE(coskr) 1295 ABI_FREE(sinkr) 1296 1297 !second 1 / elph_ds%k_fine%nkpt factor for the integration weights 1298 elph_tr_ds%a2f_1d_trin = elph_tr_ds%a2f_1d_trin / elph_ds%k_fine%nkpt 1299 elph_tr_ds%a2f_1d_trout = elph_tr_ds%a2f_1d_trout / elph_ds%k_fine%nkpt 1300 1301 if (elph_ds%ep_scalprod == 1) then 1302 write(std_out,*) 'mka2f_tr_lova: errors in diagonalization of gamma_tr with phon eigenvectors: ', nerr,maxerr 1303 end if 1304 1305 elph_tr_ds%a2f_1d_tr(:,:,:,1,1,1) = elph_tr_ds%a2f_1d_trout(:,:,:) - elph_tr_ds%a2f_1d_trin(:,:,:) 1306 1307 !output the elph_tr_ds%a2f_1d_tr 1308 fname = trim(elph_ds%elph_base_name) // '_A2F_TR' 1309 if (open_file (fname,message,newunit=unit_a2f_tr,status='unknown') /= 0) then 1310 ABI_ERROR(message) 1311 end if 1312 1313 fname = trim(elph_ds%elph_base_name) // '_A2F_TRIN' 1314 if (open_file(fname,message,newunit=unit_a2f_trin,status='unknown') /= 0) then 1315 ABI_ERROR(message) 1316 end if 1317 1318 fname = trim(elph_ds%elph_base_name) // '_A2F_TROUT' 1319 if (open_file (fname,message,newunit=unit_a2f_trout,status='unknown') /=0) then 1320 ABI_ERROR(message) 1321 end if 1322 1323 write (unit_a2f_tr,'(a)') '#' 1324 write (unit_a2f_tr,'(a)') '# ABINIT package : a2f_tr file' 1325 write (unit_a2f_tr,'(a)') '#' 1326 write (unit_a2f_tr,'(a)') '# a2f_tr function integrated over the FS. omega in a.u.' 1327 write (unit_a2f_tr,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt 1328 write (unit_a2f_tr,'(a,I10)') '# number of energy points : ',elph_ds%na2f 1329 write (unit_a2f_tr,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 1330 write (unit_a2f_tr,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 1331 write (unit_a2f_tr,'(a)') '#' 1332 1333 write (unit_a2f_trin,'(a)') '#' 1334 write (unit_a2f_trin,'(a)') '# ABINIT package : a2f_trin file' 1335 write (unit_a2f_trin,'(a)') '#' 1336 write (unit_a2f_trin,'(a)') '# a2f_trin function integrated over the FS. omega in a.u.' 1337 write (unit_a2f_trin,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt 1338 write (unit_a2f_trin,'(a,I10)') '# number of energy points : ',elph_ds%na2f 1339 write (unit_a2f_trin,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 1340 write (unit_a2f_trin,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 1341 write (unit_a2f_trin,'(a)') '#' 1342 1343 write (unit_a2f_trout,'(a)') '#' 1344 write (unit_a2f_trout,'(a)') '# ABINIT package : a2f_trout file' 1345 write (unit_a2f_trout,'(a)') '#' 1346 write (unit_a2f_trout,'(a)') '# a2f_trout function integrated over the FS. omega in a.u.' 1347 write (unit_a2f_trout,'(a,I10)') '# number of kpoints integrated over : ', elph_ds%k_fine%nkpt 1348 write (unit_a2f_trout,'(a,I10)') '# number of energy points : ',elph_ds%na2f 1349 write (unit_a2f_trout,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', omega_min,' Ha and omega_max = ', omega_max, ' Ha' 1350 write (unit_a2f_trout,'(a,E16.6)') '# and the smearing width for gaussians is ', elph_ds%a2fsmear 1351 write (unit_a2f_trout,'(a)') '#' 1352 1353 !done with header 1354 do isppol=1,elph_ds%nsppol 1355 write (unit_a2f_tr,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol) 1356 write (unit_a2f_trin,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol) 1357 write (unit_a2f_trout,'(a,E16.6)') '# The DOS at Fermi level is ', elph_ds%n0(isppol) 1358 ! omega = zero 1359 omega = omega_min 1360 do iomega=1,elph_ds%na2f 1361 write (unit_a2f_tr, '(10D16.6)') omega, elph_tr_ds%a2f_1d_tr (iomega,:,isppol,1,1,1) 1362 write (unit_a2f_trin, '(10D16.6)') omega, elph_tr_ds%a2f_1d_trin (iomega,:,isppol) 1363 write (unit_a2f_trout,'(10D16.6)') omega, elph_tr_ds%a2f_1d_trout(iomega,:,isppol) 1364 omega=omega+domega 1365 end do 1366 write (unit_a2f_tr,*) 1367 write (unit_a2f_trin,*) 1368 write (unit_a2f_trout,*) 1369 end do !isppol 1370 1371 close (unit=unit_a2f_tr) 1372 close (unit=unit_a2f_trin) 1373 close (unit=unit_a2f_trout) 1374 1375 !calculation of transport properties 1376 ABI_MALLOC(integrho,(elph_ds%na2f)) 1377 ABI_MALLOC(tointegrho,(elph_ds%na2f)) 1378 ABI_MALLOC(tointega2f,(elph_ds%na2f)) 1379 ABI_MALLOC(integtau,(elph_ds%na2f)) 1380 ABI_MALLOC(tointegtau,(elph_ds%na2f)) 1381 1382 fname = trim(elph_ds%elph_base_name) // '_RHO' 1383 if (open_file(fname,message,newunit=unit_rho,status='unknown') /= 0) then 1384 ABI_ERROR(message) 1385 end if 1386 1387 !print header to resistivity file 1388 write (unit_rho,*) '# Resistivity as a function of temperature.' 1389 write (unit_rho,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 1390 write (unit_rho,*) '# ' 1391 write (unit_rho,*) '# Columns are: ' 1392 write (unit_rho,*) '# temperature[K] rho[au] rho [SI] rho/temp [au]' 1393 write (unit_rho,*) '# ' 1394 1395 fname = trim(elph_ds%elph_base_name) // '_TAU' 1396 if (open_file(fname,message,newunit=unit_tau,status='unknown') /= 0) then 1397 ABI_ERROR(message) 1398 end if 1399 1400 !print header to relaxation time file 1401 write (unit_tau,*) '# Relaxation time as a function of temperature.' 1402 write (unit_tau,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 1403 write (unit_tau,*) '# ' 1404 write (unit_tau,*) '# Columns are: ' 1405 write (unit_tau,*) '# temperature[K] tau[au] tau [femtosecond] ' 1406 write (unit_tau,*) '# ' 1407 1408 fname = trim(elph_ds%elph_base_name) // '_WTH' 1409 if (open_file(fname,message,newunit=unit_therm,status='unknown') /= 0) then 1410 ABI_ERROR(message) 1411 end if 1412 1413 !print header to thermal conductivity file 1414 write (unit_therm,'(a)') '# Thermal conductivity/resistivity as a function of temperature.' 1415 write (unit_therm,'(a)') '# the formalism is isotropic, so non-cubic crystals may be wrong' 1416 write (unit_therm,'(a)') '# ' 1417 write (unit_therm,'(a)') '# Columns are: ' 1418 write (unit_therm,'(a)') '# temperature[K] thermal rho[au] thermal cond [au] thermal rho [SI] thermal cond [SI]' 1419 write (unit_therm,'(a)') '# ' 1420 1421 fname = trim(elph_ds%elph_base_name) // '_LOR' 1422 if (open_file(fname,message,newunit=unit_lor,status='unknown') /= 0) then 1423 ABI_ERROR(message) 1424 end if 1425 1426 !print header to lorentz file 1427 write (unit_lor,*) '# Lorentz number as a function of temperature.' 1428 write (unit_lor,*) '# the formalism is isotropic, so non-cubic crystals may be wrong' 1429 write (unit_lor,*) '# ' 1430 write (unit_lor,*) '# Columns are: ' 1431 write (unit_lor,*) '# temperature[K] Lorentz number[au] Lorentz quantum = (pi*kb_HaK)**2/3' 1432 write (unit_lor,*) '# ' 1433 1434 do isppol=1,elph_ds%nsppol 1435 lambda_tr_trace = zero 1436 do itrtensor=1,9 1437 omega = omega_min 1438 tointega2f = zero 1439 do iomega=1,elph_ds%na2f 1440 if(omega<=0) then 1441 omega=omega+domega 1442 cycle 1443 end if 1444 tointega2f(iomega)=elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)/omega 1445 omega=omega+domega 1446 end do 1447 1448 integrho = zero 1449 call simpson_int(elph_ds%na2f,domega,tointega2f,integrho) 1450 lambda_tr = two * spinfact * integrho(elph_ds%na2f) 1451 write (message, '(a,2i3,a,es16.6)' )& 1452 & ' mka2f_tr_lova : TRANSPORT lambda for isppol itrtensor', isppol, itrtensor, ' = ', lambda_tr 1453 call wrtout(std_out,message,'COLL') 1454 if (itrtensor == 1 .or. itrtensor == 5 .or. itrtensor == 9) lambda_tr_trace = lambda_tr_trace + lambda_tr 1455 end do !end itrtensor do 1456 1457 lambda_tr_trace = lambda_tr_trace / three 1458 write (message, '(a,i3,a,es16.6)' )& 1459 & ' mka2f_tr_lova: 1/3 trace of TRANSPORT lambda for isppol ', isppol, ' = ', lambda_tr_trace 1460 call wrtout(std_out,message,'COLL') 1461 call wrtout(ab_out,message,'COLL') 1462 end do !end isppol do 1463 1464 !constant to change units of rho from au to SI 1465 chgu=2.173969d-7 1466 chtu=2.4188843265d-17 1467 femto=1.0d-15 1468 1469 do isppol=1,elph_ds%nsppol 1470 do icomp=1, 3 1471 do jcomp=1, 3 1472 itrtensor=(icomp-1)*3+jcomp 1473 1474 ! prefactor for resistivity integral 1475 ! firh=6.d0*pi*crystal%ucvol*kb_HaK/(elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol)) 1476 ! FIXME: check factor of 2 which is different from Savrasov paper. 6 below for thermal conductivity is correct. 1477 firh=2.d0*pi*crystal%ucvol*kb_HaK/elph_ds%n0(isppol)/& 1478 & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol)) 1479 1480 ! Add by BX to get Tau_elph 1481 firh_tau = 2.0d0*pi*kb_HaK 1482 ! End Adding 1483 1484 write(unit_rho,*) '# Rho for isppol, itrten = ', isppol, itrtensor 1485 write(unit_tau,*) '# Tau for isppol, itrten = ', isppol, itrtensor 1486 1487 ! jmb 1488 tointegtau(:)=0. 1489 tointegrho(:)=0. 1490 do itemp=1,ntemper ! runs over termperature in K 1491 Temp=tempermin+temperinc*dble(itemp) 1492 firhT=firh*Temp 1493 firhT_tau=firh_tau*Temp 1494 omega = omega_min 1495 do iomega=1,elph_ds%na2f 1496 if(omega<=0) then 1497 omega=omega+domega 1498 cycle 1499 end if 1500 xtr=omega/(2*kb_HaK*Temp) 1501 if(xtr < log(huge(zero)*tol16)/2)then 1502 tointegrho(iomega)=spinfact*firhT*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1) & 1503 & /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2) 1504 ! Add by BX to get Tau 1505 tointegtau(iomega)=spinfact*firhT_tau*omega*elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1) & 1506 & /(((2*Temp*kb_HaK)**2)*((exp(xtr)-exp(-xtr))/2)**2) 1507 else 1508 tointegrho(iomega)=zero 1509 tointegtau(iomega)=zero 1510 end if 1511 omega=omega+domega 1512 end do 1513 1514 call simpson_int(elph_ds%na2f,domega,tointegrho,integrho) 1515 call simpson_int(elph_ds%na2f,domega,tointegtau,integtau) 1516 rho=integrho(elph_ds%na2f) 1517 tau=1.0d99 1518 if(dabs(integtau(elph_ds%na2f)) < tol7) then 1519 write(message,'(a)') ' Cannot get a physical relaxation time ' 1520 ABI_WARNING(message) 1521 else 1522 tau=1.0d0/integtau(elph_ds%na2f) 1523 end if 1524 ! if(elph_ds%na2f < 350.0) then 1525 ! tau=1.0d0/integtau(elph_ds%na2f) 1526 ! end if 1527 write(unit_rho,'(4D20.10)')temp,rho,rho*chgu,rho/temp 1528 write(unit_tau,'(3D20.10)')temp,tau,tau*chtu/femto 1529 rho_T(itemp)=rho 1530 tau_T(itemp)=tau 1531 end do ! temperature 1532 write(unit_rho,*) 1533 write(unit_tau,*) 1534 1535 end do ! jcomp 1536 end do ! icomp 1537 end do ! isppol 1538 1539 !----------------------------- 1540 1541 1542 do isppol=1,elph_ds%nsppol 1543 do icomp=1, 3 1544 do jcomp=1, 3 1545 itrtensor=(icomp-1)*3+jcomp 1546 ! prefactor for integral of thermal conductivity 1547 ! firh=(18.*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol)*elph_tr_ds%FSelecveloc_sq(isppol)) 1548 firh=(6.d0*crystal%ucvol)/(pi*kb_HaK*elph_ds%n0(isppol))/ & 1549 & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol)) 1550 1551 1552 write(unit_therm,*) '# Thermal resistivity for isppol, itrten= ', isppol 1553 write(unit_lor,*) '# Lorentz coefficient for isppol, itrten= ', isppol 1554 1555 tointegrho(:)=0. 1556 do itemp=1,ntemper 1557 1558 Temp=tempermin + temperinc*dble(itemp) 1559 omega = omega_min 1560 do iomega=1,elph_ds%na2f 1561 if(omega<=0) then 1562 omega=omega+domega 1563 cycle 1564 end if 1565 xtr=omega/(2*kb_HaK*Temp) 1566 if(xtr < log(huge(zero)*tol16)/2)then 1567 tointegrho(iomega) = spinfact*xtr**2/omega*& 1568 & ( elph_tr_ds%a2f_1d_tr(iomega,itrtensor,isppol,1,1,1)+& 1569 & 4*xtr**2*elph_tr_ds%a2f_1d_trout(iomega,itrtensor,isppol)/pi**2+ & 1570 & 2*xtr**2*elph_tr_ds%a2f_1d_trin(iomega,itrtensor,isppol)/pi**2) & 1571 & /(((exp(xtr)-exp(-xtr))/2)**2) 1572 else 1573 tointegrho(iomega) = zero 1574 end if 1575 omega=omega+domega 1576 end do 1577 1578 call simpson_int(elph_ds%na2f,domega,tointegrho,integrho) 1579 wtherm=integrho(elph_ds%na2f)*firh 1580 1581 if(abs(wtherm) > tol12)then 1582 write(unit_therm,'(5D20.10)')temp,wtherm,1./wtherm,wtherm/3.4057d9,1./(wtherm) *3.4057d9 1583 1584 lorentz=rho_T(itemp)/(wtherm*temp) 1585 write(unit_lor,*)temp,lorentz,lor0 1586 else 1587 write(unit_therm,'(5D20.10)')temp,zero,huge(one),zero,huge(one) 1588 write(unit_lor,*)temp,huge(one),lor0 1589 end if 1590 1591 end do 1592 write(unit_therm,*) 1593 write(unit_lor,*) 1594 end do ! jcomp 1595 end do ! icomp 1596 end do !end isppol do 1597 1598 1599 ABI_FREE(phfrq) 1600 ABI_FREE(displ) 1601 ABI_FREE(pheigvec) 1602 ABI_FREE(rho_T) 1603 ABI_FREE(tau_T) 1604 1605 close (unit=unit_lor) 1606 close (unit=unit_rho) 1607 close (unit=unit_tau) 1608 close (unit=unit_therm) 1609 1610 ABI_FREE(integrho) 1611 ABI_FREE(integtau) 1612 ABI_FREE(tointega2f) 1613 ABI_FREE(tointegrho) 1614 ABI_FREE(tointegtau) 1615 ABI_FREE(elph_tr_ds%a2f_1d_tr) 1616 ABI_FREE(elph_tr_ds%a2f_1d_trin) 1617 ABI_FREE(elph_tr_ds%a2f_1d_trout) 1618 1619 ABI_FREE(elph_tr_ds%gamma_qpt_trin) 1620 ABI_FREE(elph_tr_ds%gamma_qpt_trout) 1621 ABI_FREE(elph_tr_ds%gamma_rpt_trin) 1622 ABI_FREE(elph_tr_ds%gamma_rpt_trout) 1623 1624 !DEBUG 1625 write(std_out,*) ' mka2f_tr_lova : end ' 1626 !ENDDEBUG 1627 1628 end subroutine mka2f_tr_lova