TABLE OF CONTENTS
- ABINIT/m_conducti
- m_conducti/conducti_nc
- m_conducti/conducti_paw
- m_conducti/conducti_paw_core
- m_conducti/msig
ABINIT/m_conducti [ Modules ]
NAME
m_conducti
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula.
COPYRIGHT
Copyright (C) 2002-2024 ABINIT group (VRecoules, PGhosh, SMazevet, SM, SVinko, NBrouwer) 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
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_conducti 25 26 use defs_basis 27 use m_errors 28 use m_abicore 29 use m_xmpi 30 use m_wfk 31 use m_hdr 32 use m_nctk 33 #ifdef HAVE_NETCDF 34 use netcdf 35 #endif 36 37 use defs_abitypes, only : MPI_type 38 use m_io_tools, only : open_file, close_unit, get_unit 39 use m_fstrings, only : sjoin 40 use m_symtk, only : matr3inv 41 use m_hide_lapack, only : jacobi 42 use m_occ, only : getnel 43 use m_geometry, only : metric 44 use m_splines, only : intrpl,splint,spline 45 use m_mpinfo, only : distrb2,init_mpi_enreg,destroy_mpi_enreg,proc_distrb_cycle,initmpi_band 46 47 implicit none 48 49 private
m_conducti/conducti_nc [ Functions ]
[ Top ] [ m_conducti ] [ Functions ]
NAME
conducti_nc
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula.
INPUTS
(main routine)
OUTPUT
(main routine)
NOTES
bantot doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy. dom=frequency range eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree). eigen11(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 100 eigen12(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 010 eigen13(2*mband*mband*nkpt_rbz*nsppol)=first-order eigenvalues (hartree) in reciprocal direction 001 ecut=kinetic energy planewave cutoff (hartree). entropy= entropy associated with the smearing (adimensional) fermie= fermi energy (Hartree) gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$). gmet_inv(3,3)=inverse of reciprocal space metric ($\textrm{bohr}^{2}$). gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1). kin11= Onsager kinetic coeficient=optical conductivity kin12= Onsager kinetic coeficient kin21= Onsager kinetic coeficient kin22= Onsager kinetic coeficient Kth=thermal conductivity mom=number of frequency for conductivity computation mband=maximum number of bands. natom = number of atoms in the unit cell. nband(nkpt*nsppol)=number of bands at each RF k point for each spin. nelect=number of electrons per unit cell nkpt=number of k points in the IBZ for this perturbation ngfft(3)=integer fft box dimensions. nspinor=number of spinorial components of the wavefunctions. nsppol=1 for unpolarized, 2 for spin-polarized. ntypat = number of atom types. occ(mband*nkpt*nsppol)=occupation number for each band and k. occopt==option for occupancies rmet(3,3)=real space metric ($\textrm{bohr}^{2}$). rprimd(3,3)=real space primitive translations. of primitive translations. Sth=thermopower tsmear=smearing width (or temperature) in Hartree ucvol=unit cell volume in ($\textrm{bohr}^{3}$). wind=frequency windows for computations of sigma wtk(nkpt)=weight assigned to each k point. znucl(natom)=atomic number of atoms np_sum=noziere-pines sumrule cond_kg(mom)=kubo-greenwood conductivity
SOURCE
1769 subroutine conducti_nc(filnam,filnam_out) 1770 1771 !Arguments ----------------------------------- 1772 !scalars 1773 character(len=fnlen) :: filnam,filnam_out 1774 1775 !Local variables------------------------------- 1776 !scalars 1777 integer,parameter :: formeig0=0,formeig1=1 1778 integer :: bantot,bd2tot_index,bdtot0_index,bdtot_index 1779 integer :: headform,iband,ii,jj,ikpt,iunt 1780 integer :: index_1,iom,isppol,jband,l1,l2,mband,mom,natom,nband1 1781 integer :: nrot,iomode 1782 integer :: nband_k,nkpt,nlign,nrest,nspinor,nsppol,ntypat 1783 integer :: occopt,comm 1784 integer :: tens_unt,lij_unt,sig_unt,kth_unt,ocond_unt 1785 real(dp) :: deltae,dosdeltae,diff_occ,dom,ecut,entropy,fermie,maxocc 1786 real(dp) :: nelect,np_sum,np_sum_k1,np_sum_k2,omin,oml,socc,socc_k,sig 1787 real(dp) :: tphysel,tsmear,ucvol,wind,Tatm 1788 character(len=fnlen) :: filnam0,filnam1,filnam2,filnam3 1789 character(len=500) :: msg 1790 type(hdr_type) :: hdr 1791 type(wfk_t) :: gswfk,ddk1,ddk2,ddk3 1792 !arrays 1793 integer,allocatable :: nband(:) 1794 real(dp) :: gmet(3,3),gmet_inv(3,3),gprimd(3,3),gprimd_inv(3,3),rmet(3,3),rprimd(3,3) 1795 real(dp),allocatable :: cond_kg(:,:,:),cond_kg_cart(:,:,:),cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:) 1796 real(dp),allocatable :: doccde(:),doccde_k(:),cond_kg_xx(:),cond_kg_yy(:),cond_kg_zz(:),trace(:) 1797 real(dp),allocatable :: eig0_k(:),eig0tmp(:),eig1_k(:,:),eigen0(:),eigen11(:) 1798 real(dp),allocatable :: eigen12(:),eigtmp(:) 1799 real(dp),allocatable :: eigen13(:),occ(:),occ_k(:),wtk(:),cond_tot(:),oml1(:) 1800 real(dp),allocatable :: kin11(:),kin12(:),kin21(:),kin22(:) 1801 real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:) 1802 real(dp) :: cond_kg_w(3,3),z(3,3) 1803 real(dp) :: eig_cond(3) 1804 1805 ! ********************************************************************************* 1806 1807 !Read data file 1808 if (open_file(filnam,msg,newunit=iunt,form='formatted',status="old")/=0) then 1809 ABI_ERROR(msg) 1810 end if 1811 1812 rewind(iunt) 1813 read(iunt,*) 1814 read(iunt,'(a)')filnam1 ! first ddk file 1815 read(iunt,'(a)')filnam2 ! second ddk file 1816 read(iunt,'(a)')filnam3 ! third ddk file 1817 read(iunt,'(a)')filnam0 ! ground-state data 1818 1819 !Open the GS Wavefunction file and the 3 DDK files. 1820 ! TODO: one should perform basic consistency tests for the GS WFK and the DDK files, e.g. 1821 ! k-points and their order, spins, number of bands could differ in the four files. 1822 ! Note indeed that we are assuming the same numer of bands in all the files. 1823 comm = xmpi_comm_self 1824 call nctk_fort_or_ncfile(filnam0, iomode, msg) 1825 if (len_trim(msg) /= 0) ABI_ERROR(msg) 1826 call wfk_open_read(gswfk,filnam0, formeig0, iomode, get_unit(), comm) 1827 1828 call nctk_fort_or_ncfile(filnam1, iomode, msg) 1829 if (len_trim(msg) /= 0) ABI_ERROR(msg) 1830 call wfk_open_read(ddk1,filnam1, formeig1, iomode, get_unit(), comm, hdr_out=hdr) 1831 1832 call nctk_fort_or_ncfile(filnam2, iomode, msg) 1833 if (len_trim(msg) /= 0) ABI_ERROR(msg) 1834 call wfk_open_read(ddk2,filnam2, formeig1, iomode, get_unit(), comm) 1835 1836 call nctk_fort_or_ncfile(filnam3, iomode, msg) 1837 if (len_trim(msg) /= 0) ABI_ERROR(msg) 1838 call wfk_open_read(ddk3,filnam3, formeig1, iomode, get_unit(), comm) 1839 1840 if (ddk1%compare(ddk2) /= 0) then 1841 ABI_ERROR("ddk1 and ddk2 are not consistent. see above messages") 1842 end if 1843 if (ddk1%compare(ddk3) /= 0) then 1844 ABI_ERROR("ddk1 and ddk3 are not consistent. see above messages") 1845 end if 1846 1847 !Extract params from the header of the first ddk file (might have been the GS file ?) 1848 1849 !Extract info from the header 1850 headform=hdr%headform 1851 bantot=hdr%bantot 1852 ecut=hdr%ecut_eff 1853 natom=hdr%natom 1854 nkpt=hdr%nkpt 1855 nspinor=hdr%nspinor 1856 nsppol=hdr%nsppol 1857 ntypat=hdr%ntypat 1858 occopt=hdr%occopt 1859 rprimd(:,:)=hdr%rprimd(:,:) 1860 ABI_MALLOC(nband,(nkpt*nsppol)) 1861 ABI_MALLOC(occ,(bantot)) 1862 fermie=hdr%fermie 1863 occ(1:bantot)=hdr%occ(1:bantot) 1864 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 1865 1866 !Get mband, as the maximum value of nband(nkpt) 1867 mband=maxval(nband(:)) 1868 1869 write(std_out,*) 1870 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1) 1871 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2) 1872 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3) 1873 write(std_out,'(a,i8)') ' natom =',natom 1874 write(std_out,'(a,2i8)') ' nkpt,mband =',nkpt,mband 1875 write(std_out,'(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 1876 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 1877 1878 !Prepare the reading of ddk Wff files 1879 ABI_MALLOC(eigtmp,(2*mband*mband)) 1880 ABI_MALLOC(eig0tmp,(mband)) 1881 1882 !Read the eigenvalues of ground-state and ddk files 1883 ABI_MALLOC(eigen0,(mband*nkpt*nsppol)) 1884 ABI_MALLOC(eigen11,(2*mband*mband*nkpt*nsppol)) 1885 ABI_MALLOC(eigen12,(2*mband*mband*nkpt*nsppol)) 1886 ABI_MALLOC(eigen13,(2*mband*mband*nkpt*nsppol)) 1887 bdtot0_index=0 ; bdtot_index=0 1888 do isppol=1,nsppol 1889 do ikpt=1,nkpt 1890 nband1=nband(ikpt+(isppol-1)*nkpt) 1891 call gswfk%read_eigk(ikpt,isppol,xmpio_single,eig0tmp) 1892 eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1) 1893 1894 call ddk1%read_eigk(ikpt,isppol,xmpio_single,eigtmp) 1895 eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 1896 1897 call ddk2%read_eigk(ikpt,isppol,xmpio_single,eigtmp) 1898 eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 1899 1900 call ddk3%read_eigk(ikpt,isppol,xmpio_single,eigtmp) 1901 eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2) 1902 1903 bdtot0_index=bdtot0_index+nband1 1904 bdtot_index=bdtot_index+2*nband1**2 1905 end do 1906 end do 1907 1908 !Close files 1909 call gswfk%close() 1910 call ddk1%close() 1911 call ddk2%close() 1912 call ddk3%close() 1913 1914 ABI_FREE(eigtmp) 1915 ABI_FREE(eig0tmp) 1916 1917 !--------------------------------------------------------------------------------- 1918 !Gmet inversion 1919 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1920 call matr3inv(gmet,gmet_inv) 1921 call matr3inv(gprimd,gprimd_inv) 1922 1923 !--------------------------------------------------------------------------------- 1924 !Derivative of occupation wrt the energy. 1925 1926 ABI_MALLOC(doccde,(mband*nkpt*nsppol)) 1927 ABI_MALLOC(wtk,(nkpt)) 1928 1929 read(iunt,*)tsmear 1930 Tatm=tsmear*Ha_K 1931 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 1932 ! 1933 nlign=nkpt/6 1934 nrest=nkpt-6*nlign 1935 index_1=0 1936 do ii=1,nlign 1937 read(iunt,*)wtk(1+index_1:6+index_1) 1938 index_1=index_1+6 1939 end do 1940 if (nrest/=0) then 1941 read(iunt,*)wtk(6*nlign+1:nkpt) 1942 end if 1943 ! 1944 if (occopt==1) then 1945 write(std_out,'(a,i4)') ' occopt =',occopt 1946 doccde=zero 1947 else 1948 tphysel=zero 1949 maxocc=two/(nsppol*nspinor) 1950 dosdeltae=zero 1951 ! CP: using 1 and nband(0) as dummy value, because function 1952 ! not implemented for occopt==9; adding fermih=fermie in the list of arguments as well 1953 call getnel(doccde,dosdeltae,eigen0,entropy,fermie,fermie,maxocc,mband,nband,& 1954 & nelect,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,11,wtk,1,nband(1)) 1955 end if 1956 1957 !--------------------------------------------------------------------------------- 1958 !Size of the frequency range 1959 1960 read(iunt,*)dom,wind 1961 close(iunt) 1962 mom=int(wind/dom) 1963 ABI_MALLOC(oml1,(mom)) 1964 do iom=1,mom 1965 oml1(iom)=tol10*1000._dp+dble(iom)*dom 1966 end do 1967 1968 ABI_MALLOC(cond_nd,(mom,3,3)) 1969 ABI_MALLOC(cond_kg,(mom,3,3)) 1970 ABI_MALLOC(cond_kg_cart,(mom,3,3)) 1971 ABI_MALLOC(cond_kg_xx,(mom)) 1972 ABI_MALLOC(cond_kg_yy,(mom)) 1973 ABI_MALLOC(trace,(mom)) 1974 ABI_MALLOC(cond_kg_zz,(mom)) 1975 ABI_MALLOC(cond_tot,(mom)) 1976 ABI_MALLOC(kin11,(mom)) 1977 ABI_MALLOC(kin12,(mom)) 1978 ABI_MALLOC(kin21,(mom)) 1979 ABI_MALLOC(kin22,(mom)) 1980 ABI_MALLOC(kin11_k,(mom)) 1981 ABI_MALLOC(kin12_k,(mom)) 1982 ABI_MALLOC(kin21_k,(mom)) 1983 ABI_MALLOC(kin22_k,(mom)) 1984 ABI_MALLOC(Kth,(mom)) 1985 ABI_MALLOC(Stp,(mom)) 1986 write(std_out,'(a,i8,2f10.5,a)')' mom,wind,dom =',mom,wind,dom,' Ha' 1987 1988 !--------------------------------------------------------------------------------- 1989 1990 kin11 = zero 1991 kin12 = zero 1992 kin21 = zero 1993 kin22 = zero 1994 np_sum = zero 1995 socc = zero 1996 cond_kg = zero 1997 1998 !LOOP OVER SPINS 1999 do isppol=1,nsppol 2000 ! 2001 bdtot_index = 0 2002 bd2tot_index = 0 2003 2004 deltae = zero 2005 ! 2006 ! BIG FAT k POINT LOOP 2007 ! 2008 do ikpt=1,nkpt 2009 2010 nband_k=nband(ikpt+(isppol-1)*nkpt) 2011 2012 ABI_MALLOC(eig0_k,(nband_k)) 2013 ABI_MALLOC(eig1_k,(2*nband_k**2,3)) 2014 ABI_MALLOC(occ_k,(nband_k)) 2015 ABI_MALLOC(doccde_k,(nband_k)) 2016 ABI_MALLOC(dhdk2_r,(3,3,nband_k,nband_k)) 2017 ABI_MALLOC(dhdk2_g,(nband_k,nband_k)) 2018 2019 cond_nd = zero 2020 kin11_k = zero 2021 kin12_k = zero 2022 kin21_k = zero 2023 kin22_k = zero 2024 np_sum_k1 = zero 2025 np_sum_k2 = zero 2026 socc_k = zero 2027 dhdk2_r = zero 2028 dhdk2_g = zero 2029 2030 ! eigenvalue for k-point 2031 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 2032 ! first derivative eigenvalues for k-point 2033 eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_k**2+bd2tot_index) 2034 eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_k**2+bd2tot_index) 2035 eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_k**2+bd2tot_index) 2036 ! occupation numbers for k-point 2037 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 2038 ! derivative of occupation number for k-point 2039 doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index) 2040 2041 ! LOOP OVER BAND 2042 do iband=1,nband_k 2043 do jband=1,nband_k 2044 ! 2045 ! TODO : replace with BLAS calls 2046 do l1=1,3 2047 do l2=1,3 2048 do ii=1,3 2049 do jj=1,3 2050 dhdk2_r(l1,l2,iband,jband)=dhdk2_r(l1,l2,iband,jband)+(rprimd(l1,ii)& 2051 & *eig1_k(2*iband-1+(jband-1)*2*nband_k,ii)*& 2052 & rprimd(l2,jj)*eig1_k(2*iband-1+(jband-1)*2*nband_k,jj)& 2053 & +rprimd(l1,ii)*eig1_k(2*iband +(jband-1)*2*nband_k,ii)*& 2054 & rprimd(l2,jj)*eig1_k(2*iband+(jband-1)*2*nband_k,jj)) 2055 end do 2056 end do 2057 end do 2058 end do 2059 2060 do l1=1,3 2061 do l2=1,3 2062 dhdk2_r(l1,l2,iband,jband)=dhdk2_r(l1,l2,iband,jband)/two_pi/two_pi 2063 end do 2064 end do 2065 2066 ! TODO: replace with BLAS calls 2067 do l1=1,3 2068 do l2=1,3 2069 dhdk2_g(iband,jband)=dhdk2_g(iband,jband)+gmet_inv(l1,l2)*( & 2070 & eig1_k(2*iband-1+(jband-1)*2*nband_k,l1)*& 2071 & eig1_k(2*iband-1+(jband-1)*2*nband_k,l2) & 2072 & +eig1_k(2*iband +(jband-1)*2*nband_k,l1)*& 2073 & eig1_k(2*iband +(jband-1)*2*nband_k,l2)) 2074 end do 2075 end do 2076 dhdk2_g(iband,jband)=dhdk2_g(iband,jband)/two_pi/two_pi 2077 2078 diff_occ = occ_k(iband)-occ_k(jband) 2079 ! if (dabs(diff_occ)>=tol8) then 2080 2081 ! Conductivity for each omega 2082 omin = zero 2083 do iom=1,mom 2084 oml=oml1(iom) 2085 if (jband>iband) then 2086 sig= dhdk2_g(iband,jband)& 2087 & *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)& 2088 & -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2)) 2089 kin11_k(iom)=kin11_k(iom)+sig 2090 kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie) 2091 kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie) 2092 kin22_k(iom)=kin22_k(iom) + & 2093 & sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie) 2094 end if 2095 do l1=1,3 2096 do l2=1,3 2097 cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(l1,l2,iband,jband)& 2098 & *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2) 2099 end do 2100 end do 2101 2102 end do 2103 2104 ! Sumrule start 2105 if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then 2106 np_sum_k1=np_sum_k1 -dhdk2_g(iband,jband)& 2107 & *(diff_occ)/(eig0_k(iband)-eig0_k(jband)) 2108 else 2109 np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(iband,jband) 2110 end if 2111 2112 2113 ! end loop over band 2114 end do 2115 socc_k=socc_k+occ_k(iband) 2116 end do 2117 2118 do iom=1,mom 2119 kin11(iom)=kin11(iom)+wtk(ikpt)*kin11_k(iom) 2120 kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom) 2121 kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom) 2122 kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom) 2123 do l1=1,3 2124 do l2=1,3 2125 cond_kg(iom,l1,l2)=cond_kg(iom,l1,l2)+wtk(ikpt)*cond_nd(iom,l1,l2) 2126 end do 2127 end do 2128 end do 2129 2130 np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2) 2131 socc=socc+wtk(ikpt)*socc_k 2132 2133 ! Validity limit 2134 deltae=deltae+(eig0_k(nband_k)-fermie) 2135 2136 bd2tot_index=bd2tot_index+2*nband_k**2 2137 bdtot_index=bdtot_index+nband_k 2138 ABI_FREE(eig0_k) 2139 ABI_FREE(eig1_k) 2140 ABI_FREE(occ_k) 2141 ABI_FREE(doccde_k) 2142 ABI_FREE(dhdk2_r) 2143 ABI_FREE(dhdk2_g) 2144 ! End loop over k 2145 end do 2146 2147 write(std_out,'(a,3f10.5)')' sumrule =',np_sum/socc/three,socc 2148 write(std_out,'(a,f10.5,a,f10.5,a)')& 2149 & ' Emax-Efermi =',deltae/dble(nkpt),' Ha',deltae/dble(nkpt)*Ha_eV,' eV' 2150 2151 !End loop over spins 2152 end do 2153 2154 cond_kg=cond_kg*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 2155 2156 2157 !Check that new output file does NOT exist 2158 !Keep this line : prevent silly (compiler ?) bug on HP 8000 2159 write(std_out,*)' conducti : call isfile ' 2160 ! 2161 if (open_file(trim(filnam_out)//'_tens',msg,newunit=tens_unt,form='formatted',action="write")/=0) then 2162 ABI_ERROR(msg) 2163 end if 2164 if (open_file(trim(filnam_out)//'_Lij',msg,newunit=lij_unt,form='formatted',action="write")/=0) then 2165 ABI_ERROR(msg) 2166 end if 2167 write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22' 2168 2169 if (open_file(trim(filnam_out)//'_sig',msg,newunit=sig_unt,form='formatted',action="write")/=0) then 2170 ABI_ERROR(msg) 2171 end if 2172 write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1' 2173 2174 if (open_file(trim(filnam_out)//'_Kth',msg,newunit=kth_unt,form='formatted',action="write")/=0) then 2175 ABI_ERROR(msg) 2176 end if 2177 write(kth_unt,'(a)')& 2178 & ' #omega(ua) hbar*omega(eV) thermal cond(ua) Kth(W/m/K) thermopower(ua) Stp(microohm/K)' 2179 2180 if (open_file(trim(filnam_out)//'.out',msg,newunit=ocond_unt,form='formatted',action="write")/=0) then 2181 ABI_ERROR(msg) 2182 end if 2183 write(ocond_unt,'(a)' )' Conducti output file:' 2184 write(ocond_unt,'(a)' )' Contains all results produced by conducti utility' 2185 write(ocond_unt,'(a)' )' ' 2186 write(ocond_unt,'(a)')' # omega(ua) cond(ua) thermal cond(ua) thermopower(ua)' 2187 2188 !Keep this line : prevent silly (compiler ?) bug on HP 8000 2189 write(std_out,*)' conducti : after call isfile ' 2190 2191 !Compute thermal conductivity and thermopower 2192 do iom=1,mom 2193 oml=oml1(iom) 2194 kin11(iom)=kin11(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 2195 kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 2196 kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 2197 kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 2198 if (dabs(kin11(iom))<10.0d-20) kin11(iom)=zero 2199 Kth(iom)=kin22(iom) 2200 Stp(iom)=zero 2201 if(kin11(iom)/=zero) then 2202 Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/kin11(iom)) 2203 Stp(iom)=kin12(iom)/(kin11(iom)*Tatm) 2204 end if 2205 if (dabs(Kth(iom))<10.0d-20) Kth(iom)=zero 2206 if (dabs(Stp(iom))<10.0d-20) Stp(iom)=zero 2207 if (dabs(kin12(iom))<10.0d-20) kin12(iom)=zero 2208 if (dabs(kin21(iom))<10.0d-20) kin21(iom)=zero 2209 if (dabs(kin22(iom))<10.0d-20) kin22(iom)=zero 2210 2211 write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9 2212 write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,kin11(iom),kin11(iom)*Ohmcm 2213 write(kth_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,Kth(iom),Kth(iom)*3.4057d9/Tatm,Stp(iom),Stp(iom)*3.6753d-2 2214 write(ocond_unt,'(1f12.5,3es22.12)') oml,kin11(iom),Kth(iom),Stp(iom) 2215 end do 2216 2217 write(tens_unt,'(a)' )' Conductivity file ' 2218 write(tens_unt,'(a)' )' ----------------- ' 2219 write(tens_unt,'(a)' )' Contain first the full conductivity tensor, for the desired set of energies,' 2220 write(tens_unt,'(a)' )' then, the three principal values, for the desired set of energies' 2221 write(tens_unt,'(a)' )' (note that eigenvalues are not directly associated with xx,yy,zz)' 2222 write(tens_unt,'(a)' )' ' 2223 2224 write(ocond_unt,'(a)' )' ' 2225 write(ocond_unt,'(a)' )' full conductivity tensor, for the desired set of energies' 2226 write(ocond_unt,'(a)' )' then, the three principal values, for the desired set of energies:' 2227 2228 do iom=1,mom 2229 oml=oml1(iom)*Ha_eV 2230 write(tens_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :' 2231 write(ocond_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :' 2232 do l1=1,3 2233 write(tens_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3) 2234 write(ocond_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3) 2235 end do 2236 end do 2237 2238 !Diagonalizing the conductivity matrix for sigma_xx,sigma_yy,sigma_zz 2239 cond_kg_xx=0d0 2240 cond_kg_yy=0d0 2241 cond_kg_zz=0d0 2242 !trace=0d0 ! Used for checking with the original version of the code 2243 do iom=1,mom 2244 oml=oml1(iom)*Ha_eV 2245 cond_kg_w=0d0 2246 do l1=1,3 2247 do l2=1,3 2248 cond_kg_w(l1,l2)=cond_kg(iom,l1,l2) 2249 end do 2250 end do 2251 call jacobi(cond_kg_w,3,3,eig_cond,z,nrot) 2252 2253 ! When the value is too small, set it to zero before printing 2254 if(abs(eig_cond(1))<tol10)eig_cond(1)=zero 2255 if(abs(eig_cond(2))<tol10)eig_cond(2)=zero 2256 if(abs(eig_cond(3))<tol10)eig_cond(3)=zero 2257 2258 cond_kg_xx(iom)=eig_cond(1) 2259 cond_kg_yy(iom)=eig_cond(2) 2260 cond_kg_zz(iom)=eig_cond(3) 2261 ! trace(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom) 2262 end do 2263 2264 !DEBUG Keep this line : prevent silly (compiler ?) bug on HP 8000 2265 !write(std_out,*)' conducti : after open ' 2266 !ENDDEBUG 2267 2268 write(tens_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.' 2269 write(tens_unt,'(a)')' ' 2270 write(tens_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)' 2271 2272 write(ocond_unt,'(a)')' ' 2273 write(ocond_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.' 2274 write(ocond_unt,'(a)')' ' 2275 write(ocond_unt,'(a)')' #omega(ua) cond_1(ua) cond_2(ua) cond_3(ua) cond_tot(ua)' 2276 2277 2278 do iom=1,mom 2279 cond_tot(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom) 2280 write(tens_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 2281 write(ocond_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 2282 end do 2283 2284 write(tens_unt,*) 2285 write(tens_unt,'(a)')' #hbar*omega(eV) cond_1(ohm.cm)-1 cond_2(ohm.cm)-1 cond_3(ohm.cm)-1 cond_t(ohm.cm)-1' 2286 write(ocond_unt,*) 2287 write(ocond_unt,'(a)')' #hbar*omega(eV) cond_1(ohm.cm)-1 cond_2(ohm.cm)-1 cond_3(ohm.cm)-1 cond_t(ohm.cm)-1' 2288 2289 do iom=1,mom 2290 oml=oml1(iom)*Ha_eV 2291 cond_tot(iom)=cond_tot(iom)*Ohmcm 2292 cond_kg_xx(iom)=cond_kg_xx(iom)*Ohmcm 2293 cond_kg_yy(iom)=cond_kg_yy(iom)*Ohmcm 2294 cond_kg_zz(iom)=cond_kg_zz(iom)*Ohmcm 2295 write(tens_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 2296 write(ocond_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom) 2297 end do 2298 2299 !Calculate the imaginary part of the conductivity (principal value) 2300 !+derived optical properties. 2301 call msig(kin11,mom,oml1,filnam_out) 2302 2303 close(tens_unt) 2304 close(lij_unt) 2305 close(sig_unt) 2306 close(kth_unt) 2307 close(ocond_unt) 2308 2309 ABI_FREE(nband) 2310 ABI_FREE(oml1) 2311 ABI_FREE(occ) 2312 ABI_FREE(eigen11) 2313 ABI_FREE(eigen12) 2314 ABI_FREE(eigen13) 2315 ABI_FREE(eigen0) 2316 ABI_FREE(doccde) 2317 ABI_FREE(wtk) 2318 ABI_FREE(cond_nd) 2319 ABI_FREE(cond_kg) 2320 ABI_FREE(cond_kg_cart) 2321 ABI_FREE(cond_kg_xx) 2322 ABI_FREE(cond_kg_yy) 2323 ABI_FREE(trace) 2324 ABI_FREE(cond_kg_zz) 2325 ABI_FREE(cond_tot) 2326 ABI_FREE(kin11) 2327 ABI_FREE(kin22) 2328 ABI_FREE(kin12) 2329 ABI_FREE(kin21) 2330 ABI_FREE(kin11_k) 2331 ABI_FREE(kin22_k) 2332 ABI_FREE(kin12_k) 2333 ABI_FREE(kin21_k) 2334 ABI_FREE(Stp) 2335 ABI_FREE(Kth) 2336 call hdr%free() 2337 2338 end subroutine conducti_nc
m_conducti/conducti_paw [ Functions ]
[ Top ] [ m_conducti ] [ Functions ]
NAME
conducti_paw
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula for PAW formalism
INPUTS
filnam=generic name for input data filnam_out=generic name for output data [varocc]=if true, read arbitrary occupations from a file
OUTPUT
Only printing
NOTES
bantot doccde(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy. dom=frequency range eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree). eigen11(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree) in direction x eigen12(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree) in direction y eigen13(2,nkpt,mband,mband,nsppol)=first-order eigenvalues (hartree) in direction z ecut=kinetic energy planewave cutoff (hartree). fermie= fermi energy (Hartree) gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{2}$). gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1). kin11= Onsager kinetic coeficient=optical conductivity kin12= Onsager kinetic coeficient kin21= Onsager kinetic coeficient kin22= Onsager kinetic coeficient Kth=thermal conductivity mom=number of frequency for conductivity computation mband=maximum number of bands. natom = number of atoms in the unit cell. nband(nkpt*nsppol)=number of bands at each RF k point for each spin. nkpt=number of k points in the IBZ for this perturbation ngfft(3)=integer fft box dimensions. nspinor=number of spinorial components of the wavefunctions. nsppol=1 for unpolarized, 2 for spin-polarized. ntypat = number of atom types. occ(mband*nkpt*nsppol)=occupation number for each band and k. occopt==option for occupancies rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor)) rprimd(3,3)=real space primitive translations. of primitive translations. Sth=thermopower tsmear=smearing width (or temperature) in Hartree ucvol=unit cell volume in ($\textrm{bohr}^{3}$). wind=frequency windows for computations of sigma wtk(nkpt)=weight assigned to each k point. znucl(natom)=atomic number of atoms np_sum=noziere-pines sumrule
SOURCE
126 subroutine conducti_paw(filnam,filnam_out,varocc) 127 128 !Arguments ----------------------------------- 129 !scalars 130 character(len=fnlen) :: filnam,filnam_out 131 integer, optional :: varocc 132 133 !Local variables------------------------------- 134 !scalars 135 integer,parameter :: master=0 136 integer :: bsize,bd_stride,dimid,iomode,bantot,bdtot_index,ncid,varid,nb_per_proc,etiq 137 integer :: comm,fform1,headform,iband,ijband,ierr,ikpt,master_band 138 integer :: iom,isppol,jband,l1,l2,mband,me,mpierr,mom 139 integer :: natom,nband_k,nkpt,nproc,nspinor,nsppol,ntypat 140 integer :: occopt,iunt,opt_unt,occ_unt,iocc,my_iband,pnp_size 141 integer :: lij_unt,sig_unt,kth_unt,ocond_unt,occunit,occnpt 142 logical :: nc_unlimited,mykpt,myband,iomode_estf_mpiio,read_half_dipoles 143 real(dp) :: dirac,del,deltae,deltae_min,deltae_min_tmp,dhdk2_g,diff_eig,diff_occ 144 real(dp) :: dosdeltae,ecut,entropy,fact_diag,fermie,fermih,kin_fact,maxocc 145 real(dp) :: np_sum,np_sum_k1,np_sum_k2,omin,omax,dom,oml,sig,socc,socc_k 146 real(dp) :: Tatm,tphysel,tsmear,ucvol,eig_in_max,eig_in_min 147 character(len=fnlen) :: filnam1,filnam_gen,occfile 148 character(len=500) :: msg 149 type(hdr_type) :: hdr 150 type(MPI_type) :: mpi_enreg 151 !arrays 152 integer :: nc_count_5(5),nc_count_6(6),nc_start_5(5),nc_start_6(6),nc_stride_5(5),nc_stride_6(6) 153 integer,allocatable :: nband(:) 154 real(dp) :: dhdk2_r(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 155 real(dp),allocatable :: cond_nd(:,:,:),cond_nd_k(:,:,:) 156 real(dp),allocatable :: doccde(:),doccde_k(:),eig0_k(:),eigen0(:),eig0nc(:,:,:) 157 real(dp),allocatable :: occ(:),occ_k(:),wtk(:),oml1(:),occ_in(:),eig_in(:),occ_tmp(:),ypp(:) 158 real(dp),allocatable :: kin11(:,:),kin12(:),kin21(:),kin22(:) 159 real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:) 160 real(dp),allocatable :: psinablapsi(:,:,:),sig_abs(:) 161 162 ! ********************************************************************************* 163 164 ! --------------------------------------------------------------------------------- 165 ! Read input data 166 167 !Global MPI communicator 168 comm = xmpi_world 169 nproc = xmpi_comm_size(comm) 170 me = xmpi_comm_rank(comm) 171 172 !Read input parameters file 173 if (me==master) then 174 if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read",status="old")/=0) then 175 ABI_ERROR(msg) 176 end if 177 rewind(iunt) 178 read(iunt,*) 179 read(iunt,'(a)') filnam_gen ! Generic name for the files 180 filnam1=trim(filnam_gen)//'_OPT' 181 ! Read frequency range 182 read(iunt,*) dom,omin,omax,mom 183 ! In case of varocc read filename of occupation datafile 184 if (present(varocc)) then 185 if (me==master) then 186 write(std_out,*) 'Warning, this undocumented feature is highly experimental' 187 write(std_out,*) 'and of limited physical validity, proceed with extreme caution!!!' 188 end if 189 read(iunt,*) occfile !filename of the non-eq distribution function 190 end if 191 close(iunt) 192 end if 193 194 !Send data to all procs 195 call xmpi_bcast(dom,master,comm,mpierr) 196 call xmpi_bcast(omin,master,comm,mpierr) 197 call xmpi_bcast(omax,master,comm,mpierr) 198 call xmpi_bcast(mom,master,comm,mpierr) 199 200 ! --------------------------------------------------------------------------------- 201 ! Read OPT file 202 203 !Check for FORTRAN/.nc OPT file and set iomode to IO_MODE_FORTRAN_MASTER/IO_MODE_ETSF 204 if (me==master) then 205 call nctk_fort_or_ncfile(filnam1,iomode,msg) 206 if (iomode/=IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER 207 end if 208 call xmpi_bcast(filnam1,master,comm,mpierr) 209 call xmpi_bcast(iomode,master,comm,mpierr) 210 211 !Open OPT file and read HEADER 212 if (me==master) then 213 if (iomode==IO_MODE_ETSF) then 214 NCF_CHECK(nctk_open_read(ncid,filnam1,xmpi_comm_self)) 215 call hdr_ncread(hdr,ncid,fform1) 216 else 217 if (open_file(filnam1,msg,newunit=opt_unt,form="unformatted",status="old")/=0) then 218 ABI_ERROR(msg) 219 end if 220 call hdr_fort_read(hdr,opt_unt,fform1,rewind=.true.) 221 end if 222 ABI_CHECK(fform1/=0,sjoin("Error while reading ",filnam1)) 223 ABI_CHECK(fform1==610.or.fform1==620,"Conducti requires an OPT file with fform=610 or 620!") 224 end if 225 call hdr%bcast(master,me,comm) 226 call xmpi_bcast(fform1,master,comm,mpierr) 227 228 !Extract info from the header 229 headform=hdr%headform 230 bantot=hdr%bantot 231 ecut=hdr%ecut_eff 232 natom=hdr%natom 233 nkpt=hdr%nkpt 234 nspinor=hdr%nspinor 235 nsppol=hdr%nsppol 236 ntypat=hdr%ntypat 237 occopt=hdr%occopt 238 rprimd(:,:)=hdr%rprimd(:,:) 239 fermie=hdr%fermie 240 tsmear=hdr%tsmear 241 ABI_MALLOC(nband,(nkpt*nsppol)) 242 ABI_MALLOC(occ,(bantot)) 243 ABI_MALLOC(wtk,(nkpt)) 244 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 245 occ(1:bantot)=hdr%occ(1:bantot) 246 wtk(1:nkpt)=hdr%wtk(1:nkpt) 247 mband=maxval(nband(:)) ! Get mband, as the maximum value of nband(nkpt) 248 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) ! Get metrics of simulation cell 249 250 !Read eigenvalues 251 ABI_MALLOC(eigen0,(mband*nkpt*nsppol)) 252 if (me==master) then 253 if (iomode==IO_MODE_ETSF) then 254 varid=nctk_idname(ncid,"eigenvalues") 255 ABI_MALLOC(eig0nc,(mband,nkpt,nsppol)) 256 #ifdef HAVE_NETCDF 257 NCF_CHECK(nf90_get_var(ncid,varid,eig0nc)) 258 #endif 259 eigen0 = reshape(eig0nc,[mband*nkpt*nsppol]) 260 ABI_FREE(eig0nc) 261 !Close file here because the rest will possibly be read with collective I/O 262 #ifdef HAVE_NETCDF 263 NCF_CHECK(nf90_close(ncid)) 264 #endif 265 else 266 read(opt_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol) 267 end if 268 end if 269 call xmpi_bcast(eigen0,master,comm,mpierr) 270 271 !In case of varocc, use arbitrary occupations 272 !Read occfile and overwrite occ by interpolating the data in the input file 273 if (present(varocc)) then 274 if (me==master) then 275 ABI_MALLOC(occ_tmp,(bantot)) 276 if (open_file(occfile,msg,newunit=occ_unt,form='formatted')/=0) then 277 ABI_ERROR(msg) 278 end if 279 ! Get units used in occfile (1=Ha, 2= eV) 280 read(occ_unt,*) occunit 281 read(occ_unt,*) occnpt 282 ABI_MALLOC(occ_in,(occnpt)) 283 ABI_MALLOC(eig_in,(occnpt)) 284 do iocc=1,occnpt 285 read(occ_unt,*) eig_in(iocc),occ_in(iocc) 286 if (occunit==2) then !Convert units from eV to Hartree 287 eig_in(iocc)=eig_in(iocc)/Ha_eV 288 occ_in(iocc)=occ_in(iocc) 289 end if 290 end do 291 close(occ_unt) 292 ! Interpolation 293 ABI_MALLOC(ypp,(occnpt)) 294 call spline(eig_in,occ_in,occnpt,zero,zero,ypp) 295 ! Interpolate neccessary values 296 eig_in_max=maxval(eig_in) 297 eig_in_min=minval(eig_in) 298 do iocc=1,bantot 299 ! Check for Extrapolation and set them to physically sound values 300 if (eigen0(iocc)<eig_in_min) then 301 occ_tmp(iocc)=one 302 else if (eigen0(iocc)>eig_in_max) then 303 occ_tmp(iocc)=zero 304 else 305 call splint(occnpt,eig_in,occ_in,ypp,1,eigen0(iocc),occ_tmp(iocc),ierr) 306 end if 307 end do 308 occ=occ_tmp 309 ! Clean up 310 ABI_FREE(ypp) 311 ABI_FREE(occ_tmp) 312 ABI_FREE(occ_in) 313 ABI_FREE(eig_in) 314 end if 315 call xmpi_bcast(occ,master,comm,mpierr) 316 occopt=2 317 end if ! varocc? 318 319 !--------------------------------------------------------------------------------- 320 ! Prepare kpt/band parallelization 321 322 call init_mpi_enreg(mpi_enreg) 323 mpi_enreg%comm_kpt=comm 324 mpi_enreg%me_kpt=me 325 mpi_enreg%nproc_spkpt=nproc 326 mpi_enreg%paralbd=1 327 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt,mband,nsppol)) 328 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt)) 329 call distrb2(mband,nb_per_proc,nband,nkpt,nproc,nsppol,mpi_enreg) 330 call initmpi_band(nkpt,mpi_enreg,nband,nkpt,nsppol) 331 332 !--------------------------------------------------------------------------------- 333 !Print some data 334 335 Tatm=tsmear*Ha_K 336 if (me==master) then 337 write(std_out,*) 338 write(std_out,'(a)' )' Input data:' 339 write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 340 write(std_out,*) 341 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1:3,1) 342 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,2) 343 write(std_out,'(a,3f10.5,a)' )' ',rprimd(1:3,3) 344 write(std_out,'(a,i8)') ' natom =',natom 345 write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol 346 write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 347 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 348 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 349 end if 350 351 ! --------------------------------------------------------------------------------- 352 ! Compute derivative of occupations wrt the energy 353 354 ABI_MALLOC(doccde,(mband*nkpt*nsppol)) 355 if (occopt==1) then 356 if (me==master) then 357 write(std_out,'(a,i4)') ' occopt =',occopt 358 end if 359 doccde=zero 360 else 361 tphysel=zero 362 maxocc=two/(nsppol*nspinor) 363 end if 364 if (occopt==2) then 365 tsmear=one 366 entropy=one 367 if (me==master) then 368 write(std_out,'(a)') ' occopt=2 fixed occupation case' 369 end if 370 else 371 call getnel(doccde,dosdeltae,eigen0,entropy,fermie,fermih,maxocc,mband,nband,& 372 & socc,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,12,wtk) 373 entropy=tsmear*entropy 374 if (me==master) then 375 write(std_out, '(a,es22.12)') ' tsmear*entropy (Ha) =',entropy 376 end if 377 entropy=entropy/socc 378 if (me==master) then 379 write(std_out, '(a,es22.12)') ' tsmear*entropy (Ha/e) =',entropy 380 end if 381 endif 382 383 !--------------------------------------------------------------------------------- 384 ! Determine the frequency range and allocate frequency-dependent arrays 385 386 del=(omax-omin)/(mom-1) 387 ABI_MALLOC(oml1,(mom)) 388 do iom=1,mom 389 oml1(iom)=omin+dble(iom-1)*del 390 end do 391 392 ABI_MALLOC(kin11,(mom,nsppol)) 393 ABI_MALLOC(kin12,(mom)) 394 ABI_MALLOC(kin21,(mom)) 395 ABI_MALLOC(kin22,(mom)) 396 ABI_MALLOC(cond_nd,(3,3,mom)) 397 ABI_MALLOC(sig_abs,(mom)) 398 ABI_MALLOC(Kth,(mom)) 399 ABI_MALLOC(Stp,(mom)) 400 kin11 = zero 401 kin12 = zero 402 kin21 = zero 403 kin22 = zero 404 cond_nd = zero 405 sig_abs = zero 406 Kth = zero 407 Stp = zero 408 409 !--------------------------------------------------------------------------------- 410 !Prepare valence-valence dipoles reading 411 412 iomode_estf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio) 413 414 !In case of netCDF access to OPT file, prepare collective I/O 415 if (iomode == IO_MODE_ETSF) then 416 if (iomode_estf_mpiio) then 417 NCF_CHECK(nctk_open_read(ncid,filnam1,comm)) 418 varid=nctk_idname(ncid,"dipole_valence_valence") 419 if (nproc>1) then 420 NCF_CHECK(nctk_set_collective(ncid,varid)) 421 end if 422 #ifdef HAVE_NETCDF 423 nc_unlimited=(nf90_inq_dimid(ncid,"unlimited_bands",dimid)==NF90_NOERR) 424 read_half_dipoles=(nf90_inq_dimid(ncid,"max_number_of_state_pairs",dimid)==NF90_NOERR) 425 #endif 426 else 427 if (me==master) then 428 NCF_CHECK(nctk_open_read(ncid,filnam1,xmpi_comm_self)) 429 varid=nctk_idname(ncid,"dipole_valence_valence") 430 !if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then 431 ! NCF_CHECK(nctk_set_collective(ncid,varid)) 432 !end if 433 #ifdef HAVE_NETCDF 434 nc_unlimited=(nf90_inq_dimid(ncid,"unlimited_bands",dimid)==NF90_NOERR) 435 read_half_dipoles=(nf90_inq_dimid(ncid,"max_number_of_state_pairs",dimid)==NF90_NOERR) 436 #endif 437 end if 438 call xmpi_bcast(nc_unlimited,master,comm,ierr) 439 call xmpi_bcast(read_half_dipoles,master,comm,ierr) 440 end if 441 if (nc_unlimited.and.read_half_dipoles) then 442 msg="The OPT file has a wrong format!" 443 ABI_BUG(msg) 444 end if 445 else 446 read_half_dipoles=(fform1==620) 447 end if 448 449 if (iomode_estf_mpiio) then 450 !If MPI-IO, store only ib elements for each jb 451 ABI_MALLOC(psinablapsi,(2,3,mband)) 452 else 453 !If not, store all pairs (or half) 454 if (read_half_dipoles) then ! only ib>=jb 455 ABI_MALLOC(psinablapsi,(2,3,(mband*(mband+1))/2)) 456 else 457 ABI_MALLOC(psinablapsi,(2,3,mband*mband)) 458 end if 459 end if 460 pnp_size=size(psinablapsi) 461 462 !--------------------------------------------------------------------------------- 463 ! Compute conductivity 464 465 ABI_MALLOC(kin11_k,(mom)) 466 ABI_MALLOC(kin12_k,(mom)) 467 ABI_MALLOC(kin21_k,(mom)) 468 ABI_MALLOC(kin22_k,(mom)) 469 ABI_MALLOC(cond_nd_k,(3,3,mom)) 470 471 np_sum = zero 472 socc = zero 473 deltae = zero 474 deltae_min = 1.d99 475 476 !LOOP OVER SPINS/K 477 bdtot_index = 0 478 do isppol=1,nsppol 479 do ikpt=1,nkpt 480 etiq=ikpt+(isppol-1)*nkpt 481 nband_k=nband(ikpt+(isppol-1)*nkpt) 482 mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) 483 master_band=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)) 484 485 ! In case of non MPI-IO, has to read all (n,m) dipoles for this k-point 486 ! Master node reads and send to relevant processor 487 if (.not.iomode_estf_mpiio.and.me==master) then 488 if (iomode==IO_MODE_ETSF) then 489 psinablapsi=zero 490 #ifdef HAVE_NETCDF 491 if (nc_unlimited) then 492 nc_start_6=[1,1,1,ikpt,isppol,1] ; nc_count_6=[2,3,mband,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1] 493 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6)) 494 else if (.not.read_half_dipoles) then 495 nc_start_6=[1,1,1,1,ikpt,isppol] ; nc_count_6=[2,3,mband,mband,1,1] ; nc_stride_6=[1,1,1,1,1,1] 496 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6)) 497 else 498 nc_start_5=[1,1,1,ikpt,isppol] ; nc_count_5=[2,3,(mband*(mband+1))/2,1,1] ; nc_stride_5=[1,1,1,1,1] 499 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5)) 500 end if 501 #endif 502 else 503 psinablapsi=zero 504 bsize=nband_k**2;if (read_half_dipoles) bsize=(nband_k*(nband_k+1))/2 505 read(opt_unt)(psinablapsi(1:2,1,ijband),ijband=1,bsize) 506 read(opt_unt)(psinablapsi(1:2,2,ijband),ijband=1,bsize) 507 read(opt_unt)(psinablapsi(1:2,3,ijband),ijband=1,bsize) 508 end if 509 if (.not.mykpt) then 510 call xmpi_exch(psinablapsi,pnp_size,master,psinablapsi,master_band,comm,etiq,ierr) 511 end if 512 end if 513 514 ! Select k-points for current proc 515 if (mykpt) then 516 517 ABI_MALLOC(eig0_k,(nband_k)) 518 ABI_MALLOC(occ_k,(nband_k)) 519 ABI_MALLOC(doccde_k,(nband_k)) 520 521 cond_nd_k = zero 522 kin11_k = zero 523 kin12_k = zero 524 kin21_k = zero 525 kin22_k = zero 526 np_sum_k1 = zero 527 np_sum_k2 = zero 528 socc_k = zero 529 530 ! k-dependent data 531 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 532 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 533 doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index) 534 535 ! In case of non MPI-IO, receive all (n,m) dipoles from master proc 536 ! Then broadcast them to all band processors 537 if (.not.iomode_estf_mpiio) then 538 if (me/=master.and.me==master_band) then 539 call xmpi_exch(psinablapsi,pnp_size,master,psinablapsi,me,comm,etiq,ierr) 540 end if 541 call xmpi_bcast(psinablapsi,master,mpi_enreg%comm_band,mpierr) 542 end if 543 544 ! LOOP OVER BANDS n 545 do iband=1,nband_k 546 547 !If MPI-IO, store only ib elements for each jb 548 !If not, store all (ib,jb) pairs 549 my_iband=merge(1,iband,iomode_estf_mpiio) 550 551 ! Select bands for current proc 552 myband=(mpi_enreg%proc_distrb(ikpt,iband,isppol)==me) 553 if (myband) then 554 555 ! In case of MPI-IO, read valence-valence dipoles for band n 556 if (iomode_estf_mpiio) then 557 #ifdef HAVE_NETCDF 558 if (nc_unlimited) then 559 nc_start_6=[1,1,iband,ikpt,isppol,1] ; nc_count_6=[2,3,1,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1] 560 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6)) 561 else if (.not.read_half_dipoles) then 562 nc_start_6=[1,1,1,iband,ikpt,isppol] ; nc_count_6=[2,3,mband,1,1,1] ; nc_stride_6=[1,1,1,1,1,1] 563 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6)) 564 else 565 nc_start_5=[1,1,(iband*(iband-1))/2+1,ikpt,isppol] ; nc_count_5=[2,3,iband,1,1] ; nc_stride_5=[1,1,1,1,1] 566 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5)) 567 end if 568 #endif 569 end if 570 571 ! LOOP OVER BANDS m 572 do jband=1,iband-1 573 diff_occ = occ_k(iband)-occ_k(jband) 574 diff_eig = eig0_k(iband)-eig0_k(jband) 575 if (dabs(diff_occ)>=tol8) then 576 577 dhdk2_r = zero 578 dhdk2_g = zero 579 580 if (read_half_dipoles) then 581 ijband=(my_iband*(my_iband-1))/2+jband 582 else 583 !psinablapsi size is mband for netCDF I/O, nband_k for Fortran I/O 584 bd_stride=merge(mband,nband_k,iomode==IO_MODE_ETSF) 585 ijband=(my_iband-1)*bd_stride+jband 586 end if 587 588 do l2=1,3 589 do l1=1,3 590 dhdk2_r(l1,l2)=dhdk2_r(l1,l2)+(& 591 & psinablapsi(1,l1,ijband)*psinablapsi(1,l2,ijband)& 592 & +psinablapsi(2,l1,ijband)*psinablapsi(2,l2,ijband)) 593 end do 594 end do 595 do l1=1,3 596 dhdk2_g=dhdk2_g & 597 & +(psinablapsi(1,l1,ijband)*psinablapsi(1,l1,ijband) & 598 & +psinablapsi(2,l1,ijband)*psinablapsi(2,l1,ijband)) 599 end do 600 601 602 !Minimal validity limit 603 deltae_min_tmp=dabs(diff_eig) 604 if ((deltae_min_tmp>=tol5).and.(deltae_min_tmp<=deltae_min)) deltae_min=deltae_min_tmp 605 606 !Conductivity for each omega - Apply KG formula 607 kin_fact=(eig0_k(iband)+eig0_k(jband))*half-(fermie+entropy) 608 do iom=1,mom 609 oml=oml1(iom) 610 dirac=(dexp(-((diff_eig+oml)/dom)**2)-dexp(-((diff_eig-oml)/dom)**2))/oml ! Take into account (n,m) and (m,n) 611 sig=dhdk2_g*diff_occ*dirac 612 kin11_k(iom)=kin11_k(iom)+sig 613 kin12_k(iom)=kin12_k(iom)-sig*kin_fact 614 kin21_k(iom)=kin21_k(iom)-sig*kin_fact 615 kin22_k(iom)=kin22_k(iom)+sig*kin_fact**2 616 do l2=1,3 617 do l1=1,3 618 cond_nd_k(l1,l2,iom)=cond_nd_k(l1,l2,iom)+dhdk2_r(l1,l2)*diff_occ*dirac 619 end do 620 end do 621 end do 622 623 !Evaluate sumrule 624 fact_diag=merge(one,two,iband==jband) 625 if (dabs(diff_eig)>=tol10) then 626 np_sum_k1=np_sum_k1 - fact_diag*dhdk2_g*diff_occ/diff_eig 627 else 628 np_sum_k2=np_sum_k2 - fact_diag*doccde_k(iband)*dhdk2_g 629 end if 630 631 end if ! diff_occ>tol8 632 end do !jband 633 634 socc_k=socc_k+occ_k(iband) 635 636 end if ! my band? 637 end do ! iband 638 639 ! Accumulate k-point contribution 640 do iom=1,mom 641 kin11(iom,isppol)=kin11(iom,isppol)+wtk(ikpt)*kin11_k(iom) 642 kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom) 643 kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom) 644 kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom) 645 cond_nd(:,:,iom)=cond_nd(:,:,iom)+wtk(ikpt)*cond_nd_k(:,:,iom) 646 end do 647 np_sum=np_sum+wtk(ikpt)*(np_sum_k1+np_sum_k2) 648 socc=socc+wtk(ikpt)*socc_k 649 650 ! Validity limit 651 deltae=deltae+(eig0_k(nband_k)-fermie) 652 653 ABI_FREE(eig0_k) 654 ABI_FREE(occ_k) 655 ABI_FREE(doccde_k) 656 657 ! End loop over kpt/spin 658 end if ! My kpt? 659 bdtot_index=bdtot_index+nband_k 660 end do ! ikpt 661 end do ! isppol 662 663 !Accumulate kpt/band contributions over processors 664 call xmpi_sum(kin11,comm,mpierr) 665 call xmpi_sum(kin12,comm,mpierr) 666 call xmpi_sum(kin21,comm,mpierr) 667 call xmpi_sum(kin22,comm,mpierr) 668 call xmpi_sum(np_sum,comm,mpierr) 669 call xmpi_sum(socc,comm,mpierr) 670 call xmpi_sum(deltae,comm,mpierr) 671 call xmpi_min(deltae_min,comm,mpierr) 672 deltae=deltae/mpi_enreg%nproc_band 673 674 ABI_FREE(psinablapsi) 675 676 !--------------------------------------------------------------------------------- 677 ! Output results 678 679 !Print file headers (only master node) 680 if (me==master) then 681 682 ! Standard output 683 write(std_out,'(a,3f10.5)')' sumrule =',np_sum/socc/three/dble(nsppol),socc 684 write(std_out,'(a,f10.5,a,f10.5,a)')& 685 & ' Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha', & 686 & deltae/dble(nkpt*nsppol)*Ha_eV,' eV' 687 write(std_out,'(a,f10.5,a,f10.5,a)')& 688 & ' DeltaE min =',deltae_min,' Ha',deltae_min*Ha_eV,' eV' 689 690 ! _Lij file 691 if (open_file(trim(filnam_out)//'_Lij',msg, newunit=lij_unt, form='formatted', action="write") /= 0) then 692 ABI_ERROR(msg) 693 end if 694 write(lij_unt,'(a)')' # omega(ua) L12 L22 L22' 695 696 ! _sig file 697 if (open_file(trim(filnam_out)//'_sig', msg, newunit=sig_unt, form='formatted', action="write") /= 0) then 698 ABI_ERROR(msg) 699 end if 700 if (nsppol==1) then 701 write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1' 702 else 703 write(sig_unt,'(2a)')' # omega(ua) hbar*omega(eV) cond(ua) cond(ohm.cm)-1',& 704 & ' cond(ohm.cm)-1 UP cond(ohm.cm)-1 DN' 705 end if 706 707 ! _Kth file 708 if (open_file(trim(filnam_out)//'_Kth', msg, newunit=kth_unt, form='formatted', action="write") /=0) then 709 ABI_ERROR(msg) 710 end if 711 write(kth_unt,'(a)')& 712 & " #Thermal conductivity following B. Holst et al Phys. Rev. B 83 (2011) 235120" 713 write(kth_unt,'(a)')& 714 & ' #omega(ua) hbar*omega(eV) thermal cond(ua) Kth(W/m/K) thermopower(ua) Stp(microohm/K)' 715 716 ! Output file 717 if (open_file(trim(filnam_out)//'.out', msg, newunit=ocond_unt, form='formatted', action="write") /= 0) then 718 ABI_ERROR(msg) 719 end if 720 write(ocond_unt,'(a)' )'#Conducti output file:' 721 write(ocond_unt,'(a)' )'#Contains all results produced by conducti utility' 722 write(ocond_unt,'(a)' )'# ' 723 write(ocond_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =' ,mom,omin,omax,dom,' Ha' 724 write(ocond_unt,'(a,3f10.5,a)' )'# rprimd(bohr) =',rprimd(1:3,1) 725 write(ocond_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,2) 726 write(ocond_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,3) 727 write(ocond_unt,'(a,i8)' ) '# natom =',natom 728 write(ocond_unt,'(a,3i8)' ) '# nkpt,mband,nsppol =',nkpt,mband,nsppol 729 write(ocond_unt,'(a, f10.5,a)' )'# ecut =',ecut,' Ha' 730 write(ocond_unt,'(a,f10.5,a,f10.5,a)' )'# fermie =',fermie,' Ha ',fermie*Ha_eV,' eV' 731 732 write(ocond_unt,'(a,f12.5,a,f12.5,a)') '# Temp =',tsmear,' Ha ',Tatm,' Kelvin' 733 write(ocond_unt,'(a,f15.5,2x,f15.5)' )'# sumrule =',np_sum/socc/three/dble(nsppol),socc 734 write(ocond_unt,'(a,f10.5,a,f10.5,a)' )& 735 & '# Emax-Efermi =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV' 736 write(ocond_unt,'(a)' )'# ' 737 write(ocond_unt,'(a)')'# omega(ua) cond(ua) thermal cond(ua) thermopower(ua)' 738 739 end if ! me==master? 740 741 !Compute (and print) thermal conductivity and thermopower 742 do iom=1,mom 743 oml=oml1(iom) 744 745 do isppol=1,nsppol 746 kin11(iom,isppol)=kin11(iom,isppol)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 747 if (dabs(kin11(iom,isppol))<tol19) kin11(iom,isppol)=zero 748 sig_abs(iom)=sig_abs(iom)+kin11(iom,isppol) 749 end do 750 751 kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 752 kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 753 kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi) 754 755 Kth(iom)=kin22(iom) 756 Stp(iom)=zero 757 if (sig_abs(iom)/=zero) then 758 Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/sig_abs(iom)) 759 Stp(iom)=kin12(iom)/(sig_abs(iom)*tsmear) 760 end if 761 762 if (dabs(Kth(iom))<tol19) Kth(iom)=zero 763 if (dabs(Stp(iom))<tol19) Stp(iom)=zero 764 if (abs(kin12(iom))<1.d-80) kin12=zero 765 if (abs(kin21(iom))<1.d-80) kin21=zero 766 if (abs(kin22(iom))<1.d-80) kin22=zero 767 768 if (me==master) then 769 write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9 770 if (nsppol==1) then 771 write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,sig_abs(iom),sig_abs(iom)*Ohmcm 772 else 773 write(sig_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,sig_abs(iom),sig_abs(iom)*Ohmcm, & 774 & kin11(iom,1)*Ohmcm,kin11(iom,2)*Ohmcm 775 end if 776 write(kth_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,Kth(iom),Kth(iom)*3.4057d9/Tatm, & 777 & Stp(iom),Stp(iom)*3.6753d-2 778 write(ocond_unt,'(1f12.5,3es22.12)') oml,sig_abs(iom),Kth(iom),Stp(iom) 779 end if 780 781 end do 782 783 !Compute the imaginary part of the conductivity (principal value) 784 ! +derived optical properties. 785 if (me==master) then 786 call msig(sig_abs,mom,oml1,filnam_out) 787 end if 788 789 !--------------------------------------------------------------------------------- 790 ! End 791 792 !Close all files 793 if (me==master) then 794 write(std_out,'(2a)')ch10,'OUTPUT' 795 write(std_out,'(a)')trim(filnam_out)//'_Lij : Onsager kinetic coefficients' 796 write(std_out,'(a)')trim(filnam_out)//'_sig : Optical conductivity' 797 write(std_out,'(a)')trim(filnam_out)//'_Kth : Thermal conductivity and thermopower' 798 write(std_out,'(a)')trim(filnam_out)//'_eps : Dielectric function' 799 write(std_out,'(a)')trim(filnam_out)//'_abs : n, k, reflectivity, absorption' 800 close(lij_unt) 801 close(sig_unt) 802 close(kth_unt) 803 close(ocond_unt) 804 end if 805 if (iomode == IO_MODE_ETSF) then 806 if (iomode_estf_mpiio.or.me==master) then 807 #ifdef HAVE_NETCDF 808 NCF_CHECK(nf90_close(ncid)) 809 #endif 810 end if 811 else if (me==master) then 812 ierr=close_unit(opt_unt,msg) 813 ABI_CHECK(ierr==0,sjoin("Error while closing ",filnam1)) 814 end if 815 816 !Release memory space 817 ABI_FREE(kin11) 818 ABI_FREE(kin22) 819 ABI_FREE(kin12) 820 ABI_FREE(kin21) 821 ABI_FREE(kin11_k) 822 ABI_FREE(kin22_k) 823 ABI_FREE(kin12_k) 824 ABI_FREE(kin21_k) 825 ABI_FREE(Stp) 826 ABI_FREE(Kth) 827 ABI_FREE(cond_nd) 828 ABI_FREE(cond_nd_k) 829 ABI_FREE(sig_abs) 830 ABI_FREE(eigen0) 831 ABI_FREE(nband) 832 ABI_FREE(oml1) 833 ABI_FREE(occ) 834 ABI_FREE(doccde) 835 ABI_FREE(wtk) 836 call hdr%free() 837 call destroy_mpi_enreg(mpi_enreg) 838 839 end subroutine conducti_paw
m_conducti/conducti_paw_core [ Functions ]
[ Top ] [ m_conducti ] [ Functions ]
NAME
conducti_paw_core
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula for PAW formalism
INPUTS
filnam=generic name for input data filnam_out=generic name for output data [with_absorption]=optiona lflag to activate the computation of absorption (_sigX file) (default=TRUE) [with_emissivity]=optiona lflag to activate the computation of emissivity (_emisX file) (default=FALSE)
OUTPUT
Only printing
NOTES
bantot dom=frequency range eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree). ecut=kinetic energy planewave cutoff (hartree). fermie= fermi energy (Hartree) mom=number of frequency for conductivity computation mband=maximum number of bands. natom = number of atoms in the unit cell. nband(nkpt*nsppol)=number of bands at each RF k point for each spin. nkpt=number of k points in the IBZ for this perturbation ngfft(3)=integer fft box dimensions. nspinor=number of spinorial components of the wavefunctions. nsppol=1 for unpolarized, 2 for spin-polarized. ntypat = number of atom types. occ(mband*nkpt*nsppol)=occupation number for each band and k. occopt==option for occupancies psinablapsi2(2,3,mband,nphicor,natom)Matrix elements = <Phi_core|Nabla|Phi_i> rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).sigx(mom,nphicor)) rprimd(3,3)=real space primitive translations. of primitive translations. ucvol=unit cell volume in ($\textrm{bohr}^{3}$). wind=frequency window for computations of sigma wtk(nkpt)=weight assigned to each k point.
SOURCE
889 subroutine conducti_paw_core(filnam,filnam_out,with_absorption,with_emissivity) 890 891 !Arguments ----------------------------------- 892 !scalars 893 character(len=fnlen) :: filnam,filnam_out 894 logical,intent(in),optional :: with_absorption,with_emissivity 895 896 !Local variables------------------------------- 897 !scalars 898 integer,parameter :: master=0 899 integer :: iomode,atnbr,bantot,bdtot_index,comm,my_iband 900 integer :: fform2,headform,iatom,iband,icor,ierr,ikpt 901 integer :: iom,isppol,l1,mband,me,mom,mpierr,j2,etiq,pnp_size 902 integer :: natom,nband_k,nkpt,nphicor,nproc,nspinor,nsppol,ntypat 903 integer :: occopt,iunt,opt2_unt,ncid,varid,master_band,nb_per_proc 904 integer :: sigx1_unt,sigx1_s_unt,sigx1_tot_unt,ems_unt,ems_s_unt,ems_tot_unt 905 logical :: iomode_estf_mpiio,myband,mykpt,need_absorption,need_emissivity 906 real(dp) :: del_sig,del_emis,deltae,diff_occ,ecut,fermie 907 real(dp) :: omin,omax,omin_sig,omax_sig,omin_emis,omax_emis 908 real(dp) :: oml,dom,dom_ctr,dom_max,dom_tan1,dom_tan2 909 real(dp) :: Tatm,tsmear,ucvol 910 character(len=fnlen) :: filnam2,filnam_gen 911 character(len=500) :: msg 912 type(hdr_type) :: hdr 913 type(MPI_type) :: mpi_enreg 914 !arrays 915 integer :: nc_count(7),nc_start(7),nc_stride(7) 916 integer,allocatable :: nband(:),ncor(:),lcor(:),kappacor(:) 917 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 918 real(dp),allocatable :: dom_var1(:,:),dhdk2_g(:) 919 real(dp),allocatable :: eig0_k(:),eigen0(:),eig0nc(:,:,:) 920 real(dp),allocatable :: energy_cor(:),edge(:) 921 real(dp),allocatable :: occ(:),occ_k(:),wtk(:) 922 real(dp),allocatable :: oml_edge(:,:),oml_emis(:) 923 real(dp),allocatable :: psinablapsi2(:,:,:,:,:) 924 real(dp),allocatable :: sigx1(:,:,:,:),sigx1_av(:,:,:),sigx1_k(:,:,:) 925 real(dp),allocatable :: sum_spin_sigx1(:,:,:),sum_spin_sigx1_av(:,:) 926 real(dp),allocatable :: emisx(:,:,:,:),emisx_av(:,:,:),emisx_k(:,:,:) 927 real(dp),allocatable :: sum_spin_emisx(:,:,:),sum_spin_emisx_av(:,:) 928 929 ! ********************************************************************************* 930 931 !optional flags 932 need_absorption=.true. ;if (present(with_absorption)) need_absorption=with_absorption 933 need_emissivity=.false.;if (present(with_emissivity)) need_emissivity=with_emissivity 934 if ((.not.need_absorption).and.(.not.need_emissivity)) return 935 936 ! --------------------------------------------------------------------------------- 937 ! Read input data 938 939 !Global MPI communicators 940 comm = xmpi_world 941 nproc = xmpi_comm_size(comm) 942 me = xmpi_comm_rank(comm) 943 944 !Read input parameters file 945 if (me==master) then 946 if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read",status="old")/=0) then 947 ABI_ERROR(msg) 948 end if 949 rewind(iunt) 950 read(iunt,*) 951 read(iunt,'(a)')filnam_gen ! Generic name for the files 952 filnam2=trim(filnam_gen)//'_OPT2' 953 ! Read frequency range 954 if (need_absorption) then 955 read(iunt,err=11,end=11,fmt=*) dom,omin,omax,mom,atnbr,dom_max,dom_ctr 956 goto 12 957 11 backspace(iunt) ; read(iunt,*) dom,omin,omax,mom,atnbr ; dom_max=zero ; dom_ctr=zero 958 12 continue 959 else if (need_emissivity) then 960 read(iunt,*) dom,omin,omax,mom,atnbr 961 dom_max=zero;dom_ctr=zero 962 end if 963 close(iunt) 964 if (abs(dom_max)>tol10.and.dom_max<dom) then 965 msg = 'dom_max must be higher than dom!' 966 ABI_ERROR(msg) 967 end if 968 end if 969 970 !Send data to all procs 971 call xmpi_bcast(dom,master,comm,mpierr) 972 call xmpi_bcast(omin,master,comm,mpierr) 973 call xmpi_bcast(omax,master,comm,mpierr) 974 call xmpi_bcast(mom,master,comm,mpierr) 975 call xmpi_bcast(atnbr,master,comm,mpierr) 976 call xmpi_bcast(dom_max,master,comm,mpierr) 977 call xmpi_bcast(dom_ctr,master,comm,mpierr) 978 979 ! --------------------------------------------------------------------------------- 980 ! Read OPT2 file 981 982 !Check for FORTRAN/.nc OPT2 file and set iomode to IO_MODE_FORTRAN_MASTER/IO_MODE_ETSF 983 if (me==master) then 984 call nctk_fort_or_ncfile(filnam2,iomode,msg) 985 if (iomode/=IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER 986 end if 987 call xmpi_bcast(filnam2,master,comm,mpierr) 988 call xmpi_bcast(iomode,master,comm,mpierr) 989 990 !Open OPT2 file and read HEADER 991 if (me==master) then 992 if (iomode==IO_MODE_ETSF) then 993 NCF_CHECK(nctk_open_read(ncid,filnam2,xmpi_comm_self)) 994 call hdr_ncread(hdr,ncid,fform2) 995 else 996 if (open_file(filnam2,msg,newunit=opt2_unt,form="unformatted",status="old")/=0) then 997 ABI_ERROR(msg) 998 end if 999 call hdr_fort_read(hdr,opt2_unt,fform2,rewind=.true.) 1000 end if 1001 ABI_CHECK(fform2/=0,sjoin("Error while reading ",filnam2)) 1002 ABI_CHECK(fform2==611.or.fform2==612.or.fform2==613,"OPT2 file format should be fform=611/612/613!") 1003 end if 1004 call hdr%bcast(master,me,comm) 1005 call xmpi_bcast(fform2,master,comm,mpierr) 1006 1007 !Extract info from the header 1008 headform=hdr%headform 1009 bantot=hdr%bantot 1010 ecut=hdr%ecut_eff 1011 natom=hdr%natom 1012 nkpt=hdr%nkpt 1013 nspinor=hdr%nspinor 1014 nsppol=hdr%nsppol 1015 ntypat=hdr%ntypat 1016 occopt=hdr%occopt 1017 rprimd(:,:)=hdr%rprimd(:,:) 1018 fermie=hdr%fermie 1019 tsmear=hdr%tsmear 1020 ABI_MALLOC(nband,(nkpt*nsppol)) 1021 ABI_MALLOC(occ,(bantot)) 1022 ABI_MALLOC(wtk,(nkpt)) 1023 occ(1:bantot)=hdr%occ(1:bantot) 1024 wtk(1:nkpt)=hdr%wtk(1:nkpt) 1025 nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol) 1026 mband=maxval(nband(:)) ! Get mband, as the maximum value of nband(nkpt) 1027 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) ! Get metrics of simulation cell 1028 1029 !Read eigenvalues 1030 ABI_MALLOC(eigen0,(mband*nkpt*nsppol)) 1031 if (me==master) then 1032 if (iomode==IO_MODE_ETSF) then 1033 varid=nctk_idname(ncid,"eigenvalues") 1034 ABI_MALLOC(eig0nc,(mband,nkpt,nsppol)) 1035 #ifdef HAVE_NETCDF 1036 NCF_CHECK(nf90_get_var(ncid,varid,eig0nc)) 1037 #endif 1038 eigen0 = reshape(eig0nc,[mband*nkpt*nsppol]) 1039 ABI_FREE(eig0nc) 1040 else 1041 read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol) 1042 end if 1043 end if 1044 call xmpi_bcast(eigen0,master,comm,mpierr) 1045 1046 !Read core states 1047 if (me==master) then 1048 if (iomode==IO_MODE_ETSF) then 1049 #ifdef HAVE_NETCDF 1050 NCF_CHECK(nctk_get_dim(ncid,"number_of_core_states",nphicor)) 1051 #endif 1052 else 1053 read(opt2_unt) nphicor 1054 end if 1055 end if 1056 call xmpi_bcast(nphicor,master,comm,mpierr) 1057 1058 ABI_MALLOC(ncor,(nphicor)) 1059 ABI_MALLOC(lcor,(nphicor)) 1060 ABI_MALLOC(kappacor,(nphicor)) 1061 ABI_MALLOC(energy_cor,(nphicor)) 1062 ABI_MALLOC(edge,(nphicor)) 1063 1064 if (me==master) then 1065 if (iomode==IO_MODE_ETSF) then 1066 #ifdef HAVE_NETCDF 1067 varid=nctk_idname(ncid,"n_quantum_number_core") 1068 NCF_CHECK(nf90_get_var(ncid,varid,ncor)) 1069 varid=nctk_idname(ncid,"l_quantum_number_core") 1070 NCF_CHECK(nf90_get_var(ncid,varid,lcor)) 1071 varid=nctk_idname(ncid,"kappa_core") 1072 NCF_CHECK(nf90_get_var(ncid,varid,kappacor)) 1073 varid=nctk_idname(ncid,"eigenvalues_core") 1074 NCF_CHECK(nf90_get_var(ncid,varid,energy_cor)) 1075 !Close here netcdf file here because the rest has to be read with collective I/O 1076 NCF_CHECK(nf90_close(ncid)) 1077 #endif 1078 else 1079 do icor=1,nphicor 1080 read(unit=opt2_unt,err=23,end=23) ncor(icor),lcor(icor),kappacor(icor),energy_cor(icor) 1081 goto 24 1082 23 backspace(opt2_unt) ; read(unit=opt2_unt) ncor(icor),lcor(icor),energy_cor(icor) 1083 kappacor(icor)=0 1084 24 continue 1085 end do 1086 end if 1087 end if ! master 1088 call xmpi_bcast(nphicor,master,comm,mpierr) 1089 call xmpi_bcast(ncor,master,comm,mpierr) 1090 call xmpi_bcast(lcor,master,comm,mpierr) 1091 call xmpi_bcast(kappacor,master,comm,mpierr) 1092 call xmpi_bcast(energy_cor,master,comm,mpierr) 1093 edge(1:nphicor)=fermie-energy_cor(1:nphicor) 1094 1095 !--------------------------------------------------------------------------------- 1096 ! Prepare kpt/band parallelization 1097 1098 call init_mpi_enreg(mpi_enreg) 1099 mpi_enreg%comm_kpt=comm 1100 mpi_enreg%me_kpt=me 1101 mpi_enreg%nproc_spkpt=nproc 1102 mpi_enreg%paralbd=1 1103 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt,mband,nsppol)) 1104 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt)) 1105 call distrb2(mband,nb_per_proc,nband,nkpt,nproc,nsppol,mpi_enreg) 1106 call initmpi_band(nkpt,mpi_enreg,nband,nkpt,nsppol) 1107 1108 !--------------------------------------------------------------------------------- 1109 !Print some data 1110 1111 Tatm=tsmear*Ha_K 1112 if (me==master) then 1113 write(std_out,*) 1114 write(std_out,'(a)')'--------------------------------------------' 1115 write(std_out,'(a,i4)') 'selected atom for X ray emission',atnbr 1116 write(std_out,'(a)')'--------------------------------------------' 1117 if (need_absorption) then 1118 if (abs(dom_max)>tol10) then 1119 write(std_out,'(a)')'************************ ARCTAN SMEARING' 1120 write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 1121 write(std_out,'(a,2f10.5,a)')' dom_max,center =',dom_max,dom_ctr,' Ha' 1122 else 1123 write(std_out,'(a)')'************************ FIXED SMEARING' 1124 write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 1125 endif 1126 end if 1127 write(std_out,*) 1128 write(std_out,'(a,3f10.5,a)' )' rprimd(bohr) =',rprimd(1,1:3) 1129 write(std_out,'(a,3f10.5,a)' )' ',rprimd(2,1:3) 1130 write(std_out,'(a,3f10.5,a)' )' ',rprimd(3,1:3) 1131 write(std_out,'(a,i8)') ' natom =',natom 1132 write(std_out,'(a,3i8)') ' nkpt,mband,nsppol =',nkpt,mband,nsppol 1133 write(std_out, '(a, f10.5,a)' ) ' ecut =',ecut,' Ha' 1134 write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie =',fermie,' Ha',fermie*Ha_eV,' eV' 1135 write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp =',tsmear,' Ha ',Tatm,' Kelvin' 1136 write(std_out,*) 1137 write(std_out,*) 1138 write(std_out,'(a)')'--------------------------------------------' 1139 write(std_out,'(a,i4)') ' Number of core orbitals nc=',nphicor 1140 do icor=1,nphicor 1141 if (kappacor(icor)==0) then 1142 write(std_out,'(a,2i4,2f10.5)') ' n, l, Energy(Ha), Edge(Ha): ', & 1143 ncor(icor),lcor(icor),energy_cor(icor),edge(icor) 1144 else 1145 if (kappacor(icor)>0) then 1146 j2=2*lcor(icor)-1 1147 write(std_out,'(a,i4,i4,a,i4,2f10.5)') ' n, j, l, Energy(Ha), Edge(Ha): ',& 1148 & ncor(icor),j2,' / 2',lcor(icor),energy_cor(icor),edge(icor) 1149 else 1150 if (kappacor(icor)<-1) then 1151 j2=2*lcor(icor)+1 1152 write(std_out,'(a,i4,i4,a,i4,2f10.5)') ' n, j, l, Energy(Ha), Edge(Ha): ',& 1153 & ncor(icor),j2,'/2',lcor(icor),energy_cor(icor),edge(icor) 1154 else 1155 write(std_out,'(a,i4,a,i4,2f10.5)') ' n, j, l, Energy(Ha), Edge(Ha): ',& 1156 & ncor(icor),' 1/2',lcor(icor),energy_cor(icor),edge(icor) 1157 end if 1158 end if 1159 end if 1160 end do 1161 write(std_out,'(a)')'--------------------------------------------' 1162 end if ! master 1163 1164 !--------------------------------------------------------------------------------- 1165 ! Determine the frequency range and allocate frequency-dependent arrays 1166 1167 if (need_absorption) then 1168 ABI_MALLOC(oml_edge,(nphicor,mom)) 1169 ABI_MALLOC(dom_var1,(nphicor,mom)) 1170 omax_sig=omax ; omin_sig=omin 1171 del_sig=(omax_sig-omin_sig)/(mom-1) 1172 do iom=1,mom 1173 do icor=1,nphicor 1174 oml_edge(icor,iom)=-energy_cor(icor)+ dble(iom-1)*del_sig + omin_sig - one 1175 if ((oml_edge(icor,iom)<=edge(icor)).or.(abs(dom_max)<=tol10)) then 1176 dom_var1(icor,iom)= dom 1177 else 1178 dom_tan1= (oml_edge(icor,iom)-edge(icor))/dom_ctr 1179 dom_tan2=dom_tan1-one/dom_tan1**2 1180 dom_var1(icor,iom)=dom+dom_max*(half+piinv*datan(dom_tan2)) 1181 endif 1182 enddo 1183 enddo 1184 ABI_MALLOC(sigx1,(nphicor,mom,natom,nsppol)) 1185 sigx1=zero 1186 end if 1187 1188 if (need_emissivity) then 1189 ABI_MALLOC(oml_emis,(mom)) 1190 omin_emis=minval(eigen0) 1191 omax_emis=maxval(eigen0) 1192 del_emis=(omax_emis-omin_emis)/(mom-1) 1193 do iom=1,mom 1194 oml_emis(iom)=omin_emis+dble(iom-1)*del_emis 1195 end do 1196 ABI_MALLOC(emisx,(nphicor,mom,natom,nsppol)) 1197 emisx=zero 1198 end if 1199 1200 !--------------------------------------------------------------------------------- 1201 ! Prepare core-valence dipoles reading 1202 1203 iomode_estf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio) 1204 1205 !In case of netCDF access to OPT file, prepare collective I/O 1206 if (iomode == IO_MODE_ETSF) then 1207 if (iomode_estf_mpiio) then 1208 NCF_CHECK(nctk_open_read(ncid,filnam2,comm)) 1209 varid=nctk_idname(ncid,"dipole_core_valence") 1210 if (nproc>1) then 1211 NCF_CHECK(nctk_set_collective(ncid,varid)) 1212 end if 1213 else if (me==master) then 1214 NCF_CHECK(nctk_open_read(ncid,filnam2,xmpi_comm_self)) 1215 varid=nctk_idname(ncid,"dipole_core_valence") 1216 !if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then 1217 ! NCF_CHECK(nctk_set_collective(ncid,varid)) 1218 !end if 1219 end if 1220 end if 1221 1222 if (iomode_estf_mpiio) then 1223 !If MPI-IO, store only elements for one band 1224 ABI_MALLOC(psinablapsi2,(2,3,nphicor,natom,1)) 1225 else 1226 !If not, store the elements for all bands 1227 ABI_MALLOC(psinablapsi2,(2,3,nphicor,natom,mband)) 1228 end if 1229 pnp_size=size(psinablapsi2) 1230 1231 !--------------------------------------------------------------------------------- 1232 ! Compute X absorption coefficient and/or X emissivity 1233 1234 if (need_absorption) then 1235 ABI_MALLOC(sigx1_k,(nphicor,mom,natom)) 1236 end if 1237 if (need_emissivity) then 1238 ABI_MALLOC(emisx_k,(nphicor,mom,natom)) 1239 end if 1240 ABI_MALLOC(dhdk2_g,(nphicor)) 1241 1242 deltae = zero 1243 1244 !LOOP OVER SPINS/K 1245 bdtot_index = 0 1246 do isppol=1,nsppol 1247 do ikpt=1,nkpt 1248 etiq=ikpt+(isppol-1)*nkpt 1249 nband_k=nband(ikpt+(isppol-1)*nkpt) 1250 mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) 1251 master_band=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)) 1252 1253 ! In case of non MPI-IO, has to read all (n,m) dipoles for this k-point 1254 ! Master node reads and send to relevant processor 1255 if (.not.iomode_estf_mpiio.and.me==master) then 1256 if (iomode==IO_MODE_ETSF) then 1257 nc_start=[1,1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1,1] 1258 nc_count=[2,3,nphicor,natom,mband,1,1] 1259 #ifdef HAVE_NETCDF 1260 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count)) 1261 #endif 1262 else 1263 psinablapsi2=zero 1264 if (fform2==612) then ! New OPT2 file format 1265 read(opt2_unt) (((psinablapsi2(1:2,1,icor,iatom,iband),icor=1,nphicor),iatom=1,natom),iband=1,nband_k) 1266 read(opt2_unt) (((psinablapsi2(1:2,2,icor,iatom,iband),icor=1,nphicor),iatom=1,natom),iband=1,nband_k) 1267 read(opt2_unt) (((psinablapsi2(1:2,3,icor,iatom,iband),icor=1,nphicor),iatom=1,natom),iband=1,nband_k) 1268 else if (fform2==613) then ! Large OPT2 file format 1269 do iband=1,nband_k 1270 read(opt2_unt) ((psinablapsi2(1:2,1,icor,iatom,iband),icor=1,nphicor),iatom=1,natom) 1271 read(opt2_unt) ((psinablapsi2(1:2,2,icor,iatom,iband),icor=1,nphicor),iatom=1,natom) 1272 read(opt2_unt) ((psinablapsi2(1:2,3,icor,iatom,iband),icor=1,nphicor),iatom=1,natom) 1273 end do 1274 else 1275 !The old writing was not efficient (indexes order is bad) 1276 do iatom=1,natom 1277 read(opt2_unt) ((psinablapsi2(1:2,1,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor) 1278 read(opt2_unt) ((psinablapsi2(1:2,2,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor) 1279 read(opt2_unt) ((psinablapsi2(1:2,3,icor,iatom,iband),iband=1,nband_k),icor=1,nphicor) 1280 end do 1281 end if 1282 end if 1283 if (.not.mykpt) then 1284 call xmpi_exch(psinablapsi2,pnp_size,master,psinablapsi2,master_band,comm,etiq,ierr) 1285 end if 1286 end if 1287 1288 ! Select k-points for current proc 1289 if (mykpt) then 1290 1291 ABI_MALLOC(eig0_k,(nband_k)) 1292 ABI_MALLOC(occ_k,(nband_k)) 1293 1294 if (need_absorption) sigx1_k=zero 1295 if (need_emissivity) emisx_k=zero 1296 1297 ! k-dependent data 1298 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 1299 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 1300 1301 ! In case of non MPI-IO, receive all (n,m) dipoles from master proc 1302 ! Then broadcast them to all band processors 1303 if (.not.iomode_estf_mpiio) then 1304 if (me/=master.and.me==master_band) then 1305 call xmpi_exch(psinablapsi2,pnp_size,master,psinablapsi2,me,comm,etiq,ierr) 1306 end if 1307 call xmpi_bcast(psinablapsi2,master,mpi_enreg%comm_band,mpierr) 1308 end if 1309 1310 ! LOOP OVER BANDS 1311 do iband=1,nband_k 1312 1313 !If MPI-IO, store only ib elements for each iband 1314 !If not, store all elements 1315 my_iband=merge(1,iband,iomode_estf_mpiio) 1316 dhdk2_g = zero 1317 1318 ! Select bands for current proc 1319 myband=(mpi_enreg%proc_distrb(ikpt,iband,isppol)==me) 1320 if (myband) then 1321 1322 diff_occ = (two/dble(nsppol*nspinor))-occ_k(iband) 1323 1324 ! In case of MPI-IO, read core-valence dipoles for band n 1325 if (iomode_estf_mpiio) then 1326 nc_start=[1,1,1,1,iband,ikpt,isppol];nc_stride=[1,1,1,1,1,1,1] 1327 nc_count=[2,3,nphicor,natom,1,1,1] 1328 #ifdef HAVE_NETCDF 1329 NCF_CHECK(nf90_get_var(ncid,varid,psinablapsi2,start=nc_start,stride=nc_stride,count=nc_count)) 1330 #endif 1331 end if 1332 1333 ! LOOP OVER ATOMS 1334 do iatom=1,natom 1335 1336 dhdk2_g = zero 1337 do icor=1,nphicor 1338 do l1=1,3 1339 dhdk2_g(icor)=dhdk2_g(icor) & 1340 & +(psinablapsi2(1,l1,icor,iatom,my_iband)*psinablapsi2(1,l1,icor,iatom,my_iband) & 1341 & +psinablapsi2(2,l1,icor,iatom,my_iband)*psinablapsi2(2,l1,icor,iatom,my_iband)) 1342 end do 1343 end do 1344 1345 ! Absoption for each omega 1346 if (need_absorption) then 1347 do iom=1,mom 1348 do icor=1,nphicor 1349 !Lorentzian var arctan 1350 sigx1_k(icor,iom,iatom)=sigx1_k(icor,iom,iatom) & 1351 & +dhdk2_g(icor)*diff_occ*oml_edge(icor,iom) & 1352 & *(dom_var1(icor,iom)/((-energy_cor(icor)+eig0_k(iband) & 1353 & -oml_edge(icor,iom))**2+(dom_var1(icor,iom)/two)**2)) 1354 end do 1355 end do 1356 end if 1357 1358 ! Emissivity for each omega 1359 if (need_emissivity) then 1360 do iom=1,mom 1361 do icor=1,nphicor 1362 oml=-energy_cor(icor)-oml_emis(iom) 1363 emisx_k(icor,iom,iatom)=emisx_k(icor,iom,iatom) & 1364 & +dhdk2_g(icor)*occ_k(iband)/oml & 1365 & *dexp(-((-energy_cor(icor)+eig0_k(iband)-oml)/dom)**2) 1366 end do 1367 end do 1368 end if 1369 1370 end do ! iatom 1371 1372 end if ! my band? 1373 end do ! iband 1374 1375 ! Accumulate k-point contribution 1376 if (need_absorption) then 1377 sigx1(1:nphicor,1:mom,1:natom,isppol)=sigx1(1:nphicor,1:mom,1:natom,isppol) & 1378 & +wtk(ikpt)*sigx1_k(1:nphicor,1:mom,1:natom) 1379 end if 1380 if (need_emissivity) then 1381 emisx(1:nphicor,1:mom,1:natom,isppol)=emisx(1:nphicor,1:mom,1:natom,isppol) & 1382 & +wtk(ikpt)*emisx_k(1:nphicor,1:mom,1:natom) 1383 end if 1384 1385 ! Validity limit 1386 deltae=deltae+eig0_k(nband_k) 1387 1388 ABI_FREE(eig0_k) 1389 ABI_FREE(occ_k) 1390 1391 ! End loop over kpt/spin 1392 end if ! My kpt? 1393 bdtot_index=bdtot_index+nband_k 1394 end do ! ikpt 1395 end do ! isppol 1396 1397 !Accumulate kpt/band contributions over processors 1398 if (need_absorption) then 1399 call xmpi_sum(sigx1,comm,mpierr) 1400 end if 1401 if (need_emissivity) then 1402 call xmpi_sum(emisx,comm,mpierr) 1403 end if 1404 call xmpi_sum(deltae,comm,mpierr) 1405 deltae=deltae/mpi_enreg%nproc_band 1406 1407 !Release some memory 1408 ABI_FREE(dhdk2_g) 1409 ABI_FREE(psinablapsi2) 1410 if (need_absorption) then 1411 ABI_FREE(sigx1_k) 1412 end if 1413 if (need_emissivity) then 1414 ABI_FREE(emisx_k) 1415 end if 1416 1417 !Close core-valence dipoles file 1418 if (iomode == IO_MODE_ETSF) then 1419 if (iomode_estf_mpiio.or.me==master) then 1420 #ifdef HAVE_NETCDF 1421 NCF_CHECK(nf90_close(ncid)) 1422 #endif 1423 end if 1424 else if (me==master) then 1425 ierr=close_unit(opt2_unt,msg) 1426 ABI_CHECK(ierr==0,sjoin("Error while closing ",filnam2)) 1427 end if 1428 1429 !--------------------------------------------------------------------------------- 1430 ! Post-processing 1431 1432 !Absorption: post-processing 1433 if (need_absorption) then 1434 1435 ! Absorption: apply scaling factor 1436 do isppol=1,nsppol 1437 do iatom=1,natom 1438 do iom=1,mom 1439 do icor=1,nphicor 1440 sigx1(icor,iom,iatom,isppol)=sigx1(icor,iom,iatom,isppol)*two_pi*dble(natom)/ucvol/two/(dabs(energy_cor(icor))) 1441 end do 1442 end do 1443 end do 1444 end do 1445 1446 ! Absorption: average over atoms 1447 ABI_MALLOC(sigx1_av,(nphicor,mom,nsppol)) 1448 sigx1_av=zero 1449 do isppol=1,nsppol 1450 do iatom=1,natom 1451 do iom=1,mom 1452 do icor=1,nphicor 1453 sigx1_av(icor,iom,isppol)=sigx1_av(icor,iom,isppol)+sigx1(icor,iom,iatom,isppol) 1454 end do 1455 end do 1456 end do 1457 end do 1458 1459 ! Absorption: spin treatment 1460 if(nsppol==2) then 1461 ABI_MALLOC(sum_spin_sigx1,(nphicor,mom,natom)) 1462 ABI_MALLOC(sum_spin_sigx1_av,(nphicor,mom)) 1463 sum_spin_sigx1=zero ; sum_spin_sigx1_av=zero 1464 do isppol=1,nsppol 1465 do iatom=1,natom 1466 do iom=1,mom 1467 do icor=1,nphicor 1468 sum_spin_sigx1(icor,iom,iatom)=sum_spin_sigx1(icor,iom,iatom) & 1469 & +sigx1(icor,iom,iatom,isppol) 1470 end do 1471 end do 1472 end do 1473 end do 1474 do isppol=1,nsppol 1475 do icor=1,nphicor 1476 do iom=1,mom 1477 sum_spin_sigx1_av(icor,iom)=sum_spin_sigx1_av(icor,iom)+sigx1_av(icor,iom,isppol) 1478 end do 1479 end do 1480 end do 1481 endif 1482 end if 1483 1484 !Emissivity: post-processing 1485 if (need_emissivity) then 1486 1487 ! Emissivity: filter low values 1488 do isppol=1,nsppol 1489 do iatom=1,natom 1490 do iom=1,mom 1491 do icor=1,nphicor 1492 if (emisx(icor,iom,iatom,isppol)<=tol16) emisx(icor,iom,iatom,isppol)=zero 1493 end do 1494 end do 1495 end do 1496 end do 1497 1498 ! Emissivity: apply scaling factor 1499 emisx=emisx*two_pi*third*dble(natom)/(dom*ucvol)*half/dsqrt(pi) 1500 1501 ! Emissivity: average over atoms 1502 ABI_MALLOC(emisx_av,(nphicor,mom,nsppol)) 1503 emisx_av=zero 1504 do isppol=1,nsppol 1505 do iatom=1,natom 1506 do iom=1,mom 1507 do icor=1,nphicor 1508 emisx_av(icor,iom,isppol)=emisx_av(icor,iom,isppol)+emisx(icor,iom,iatom,isppol) 1509 end do 1510 end do 1511 end do 1512 end do 1513 emisx_av=emisx_av/dble(natom) 1514 1515 ! Emissivity: spin treatment 1516 if(nsppol==2) then 1517 ABI_MALLOC(sum_spin_emisx,(nphicor,mom,natom)) 1518 ABI_MALLOC(sum_spin_emisx_av,(nphicor,mom)) 1519 sum_spin_emisx=zero ; sum_spin_emisx_av=zero 1520 do isppol=1,nsppol 1521 do iatom=1,natom 1522 do iom=1,mom 1523 do icor=1,nphicor 1524 sum_spin_emisx(icor,iom,iatom)=sum_spin_emisx(icor,iom,iatom) & 1525 & +emisx(icor,iom,iatom,isppol) 1526 end do 1527 end do 1528 end do 1529 end do 1530 do isppol=1,nsppol 1531 do icor=1,nphicor 1532 do iom=1,mom 1533 sum_spin_emisx_av(icor,iom)=sum_spin_emisx_av(icor,iom)+emisx_av(icor,iom,isppol) 1534 end do 1535 end do 1536 end do 1537 endif 1538 end if 1539 1540 !--------------------------------------------------------------------------------- 1541 ! Output results 1542 1543 !Only master node outputs results in files (only master node) 1544 if (me==master) then 1545 1546 ! Standard output 1547 if (need_absorption) then 1548 write(std_out,*) 'Absorption: valence state orbital energies: omin,omax',omin_sig,omax_sig 1549 end if 1550 if (need_emissivity) then 1551 write(std_out,*) 'Emissivity: valence state orbital energies: omin,omax',omin_emis,omax_emis 1552 end if 1553 if (need_absorption) then 1554 write(std_out,'(a,f10.5,a,f10.5,a)')& 1555 & ' Emax =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV' 1556 end if 1557 1558 ! _sigX file 1559 if (need_absorption) then 1560 if (open_file(trim(filnam_out)//'_sigX',msg,newunit=sigx1_unt,form='formatted',action="write")/=0) then 1561 ABI_ERROR(msg) 1562 end if 1563 if (abs(dom_max)>tol10) then 1564 write(sigx1_unt,'(a)')'#***************************************************** ARCTAN SMEARING ********' 1565 write(sigx1_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 1566 write(sigx1_unt,'(a,2f10.5,a)')'# dom_max,center =',dom_max,dom_ctr,' Ha' 1567 else 1568 write(sigx1_unt,'(a)')'#***************************************************** FIXED SMEARING ********' 1569 write(sigx1_unt,'(a,i8,3f10.5,a)')'# npts,omin,omax,width =',mom,omin,omax,dom,' Ha' 1570 endif 1571 write(sigx1_unt,'(a,3f10.5,a)' )'# rprimd(bohr) =',rprimd(1:3,1) 1572 write(sigx1_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,2) 1573 write(sigx1_unt,'(a,3f10.5,a)' )'# ',rprimd(1:3,3) 1574 write(sigx1_unt,'(a,i8)' ) '# natom =',natom 1575 write(sigx1_unt,'(a,3i8)' ) '# nkpt,mband,nsppol =',nkpt,mband,nsppol 1576 write(sigx1_unt,'(a, f10.5,a)' )'# ecut =',ecut,' Ha' 1577 write(sigx1_unt,'(a,f10.5,a,f10.5,a)' )'# fermie =',fermie,' Ha ',fermie*Ha_eV,' eV' 1578 write(sigx1_unt,'(a,f12.5,a,f12.5,a)') '# Temp =',tsmear,' Ha ',Tatm,' Kelvin' 1579 write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------' 1580 write(sigx1_unt,'(a,i4)') '# Number of core orbitals nc=',nphicor 1581 do icor=1,nphicor 1582 if (kappacor(icor)==0) then 1583 write(sigx1_unt,'(a,2i4,2f10.5)') '# n, l, Energy(Ha), Edge(Ha): ',ncor(icor),lcor(icor),energy_cor(icor),edge(icor) 1584 else 1585 if (kappacor(icor)>0) then 1586 j2=2*lcor(icor)-1 1587 write(sigx1_unt,'(a,i4,i4,a,i4,2f10.5)') '# n, j, l, Energy(Ha), Edge(Ha): ', & 1588 & ncor(icor),j2,' / 2',lcor(icor),energy_cor(icor),edge(icor) 1589 else 1590 if (kappacor(icor)<-1) then 1591 j2=2*lcor(icor)+1 1592 write(sigx1_unt,'(a,i4,i4,a,i4,2f10.5)') '# n, j, l, Energy(Ha), Edge(Ha): ', & 1593 & ncor(icor),j2,' / 2',lcor(icor),energy_cor(icor),edge(icor) 1594 else 1595 write(sigx1_unt,'(a,i4,a,i4,2f10.5)') '# n, j, l, Energy(Ha), Edge(Ha): ', & 1596 & ncor(icor),' 1 / 2',lcor(icor),energy_cor(icor),edge(icor) 1597 end if 1598 end if 1599 end if 1600 end do 1601 write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------' 1602 write(sigx1_unt,'(a,f10.5,a,f10.5,a)')& 1603 & '# Emax =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV' 1604 write(sigx1_unt,'(a)')'#----------------------------------------------------------------------------' 1605 do iom=1,mom 1606 write(sigx1_unt,'(100(1x,e15.8))') & 1607 & (oml_edge(icor,iom),sigx1_av(icor,iom,1)/dble(natom),sigx1(icor,iom,atnbr,1),icor=1,nphicor) 1608 end do 1609 close(sigx1_unt) 1610 end if 1611 1612 ! _s_sigX and tot_sigX files 1613 if (need_absorption.and.nsppol==2) then 1614 if (open_file(trim(filnam_out)//'_s_sigX',msg,newunit=sigx1_s_unt,form='formatted',action="write")/=0) then 1615 ABI_ERROR(msg) 1616 end if 1617 if (open_file(trim(filnam_out)//'_tot_sigX',msg,newunit=sigx1_tot_unt,form='formatted',action="write")/=0) then 1618 ABI_ERROR(msg) 1619 end if 1620 do iom=1,mom 1621 write(sigx1_s_unt,'(100(1x,e15.8))') & 1622 & (oml_edge(icor,iom),sigx1_av(icor,iom,nsppol)/dble(natom),sigx1(icor,iom,atnbr,nsppol),icor=1,nphicor) 1623 write(sigx1_tot_unt,'(100(1x,e15.8))') & 1624 & (oml_edge(icor,iom),sum_spin_sigx1_av(icor,iom)/dble(natom),sum_spin_sigx1(icor,iom,atnbr),icor=1,nphicor) 1625 end do 1626 close(sigx1_s_unt) 1627 close(sigx1_tot_unt) 1628 end if 1629 1630 ! _emisX file 1631 if (need_emissivity) then 1632 if (open_file(trim(filnam_out)//'_emisX',msg,newunit=ems_unt,form='formatted',action="write")/=0) then 1633 ABI_ERROR(msg) 1634 end if 1635 write(ems_unt,*) '# conducti: Xray emission spectrum, all in atomic units by default ' 1636 write(ems_unt,*) '# One block of 3 columns per core wavefunction' 1637 write(ems_unt,*) '# energy, sigx_av, sigx, etc... ' 1638 do iom=1,mom 1639 write(ems_unt,'(3(3(1x,e15.8),2x))') & 1640 & ((-energy_cor(icor)+oml_emis(iom)),emisx_av(icor,iom,1),emisx(icor,iom,atnbr,1),icor=1,nphicor) 1641 end do 1642 close(ems_unt) 1643 end if 1644 1645 ! _s_emisX and tot_emisX files 1646 if (need_emissivity.and.nsppol==2) then 1647 if (open_file(trim(filnam_out)//'_s_emisX',msg,newunit=ems_s_unt,form='formatted',action="write")/=0) then 1648 ABI_ERROR(msg) 1649 end if 1650 if (open_file(trim(filnam_out)//'_tot_emisX',msg,newunit=ems_tot_unt,form='formatted',action="write")/=0) then 1651 ABI_ERROR(msg) 1652 end if 1653 do iom=1,mom 1654 write(ems_s_unt,'(3(3(1x,e15.8),2x))') & 1655 & ((-energy_cor(icor)+oml_emis(iom)),emisx_av(icor,iom,nsppol)/dble(natom),emisx(icor,iom,atnbr,nsppol),icor=1,nphicor) 1656 write(ems_tot_unt,'(3(3(1x,e15.8),2x))') & 1657 & ((-energy_cor(icor)+oml_emis(iom)),sum_spin_emisx_av(icor,iom)/dble(natom),sum_spin_emisx(icor,iom,atnbr),icor=1,nphicor) 1658 end do 1659 close(ems_s_unt) 1660 close(ems_tot_unt) 1661 end if 1662 1663 endif ! master node 1664 1665 !--------------------------------------------------------------------------------- 1666 ! End 1667 1668 !Release memory space 1669 if (need_absorption) then 1670 ABI_FREE(sigx1) 1671 ABI_FREE(sigx1_av) 1672 if (nsppol==2) then 1673 ABI_FREE(sum_spin_sigx1) 1674 ABI_FREE(sum_spin_sigx1_av) 1675 end if 1676 ABI_FREE(dom_var1) 1677 ABI_FREE(oml_edge) 1678 end if 1679 if (need_emissivity) then 1680 ABI_FREE(emisx) 1681 ABI_FREE(emisx_av) 1682 if (nsppol==2) then 1683 ABI_FREE(sum_spin_emisx) 1684 ABI_FREE(sum_spin_emisx_av) 1685 end if 1686 ABI_FREE(oml_emis) 1687 end if 1688 ABI_FREE(ncor) 1689 ABI_FREE(lcor) 1690 ABI_FREE(kappacor) 1691 ABI_FREE(energy_cor) 1692 ABI_FREE(edge) 1693 ABI_FREE(eigen0) 1694 ABI_FREE(nband) 1695 ABI_FREE(occ) 1696 ABI_FREE(wtk) 1697 call hdr%free() 1698 call destroy_mpi_enreg(mpi_enreg) 1699 1700 end subroutine conducti_paw_core
m_conducti/msig [ Functions ]
[ Top ] [ m_conducti ] [ Functions ]
NAME
msig
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula for PAW formalism
INPUTS
fcti(npti)= conductivity, as calculated in conducti npti= number of points to calculate conductivity xi(npti)= energies where the conductivity is calculated
OUTPUT
no output, only files
NOTES
this program calculates the imaginary part of the conductivity (principal value)
+derived optical properties.
the calculation is performed on the same grid as the initial input
to calculate the principal value, a trapezoidale integration +taylor expansion to
third order is used (W.J. Thomson computer in physics vol 12 p94 1998)
two input files are needed inppv.dat (parameters) and sigma.dat (energy,sigma_1)
two output files ppsigma.dat (energy,sigma_1,sigma_2,epsilon_1,epsilon_2)
abs.dat (energy,nomega,komega,romega,absomega)
march 2002 s.mazevet
SOURCE
2373 subroutine msig(fcti,npti,xi,filnam_out_sig) 2374 2375 !Arguments ----------------------------------- 2376 !scalars 2377 integer,intent(in) :: npti 2378 !arrays 2379 real(dp),intent(in) :: fcti(npti),xi(npti) 2380 character(len=fnlen),intent(in) :: filnam_out_sig 2381 2382 !Local variables------------------------------- 2383 !scalars 2384 integer :: npt = 10000 2385 integer :: ii,ip,npt1,npt2,eps_unt,abs_unt 2386 real(dp),parameter :: del=0.001_dp,ohmtosec=9.d11 2387 real(dp) :: dx,dx1,dx2,eps1,eps2,idel,komega,pole,refl,sigma2,xsum 2388 character(len=500) :: msg 2389 !arrays 2390 real(dp),allocatable :: abso(:),fct(:),fct1(:),fct2(:),fct3(:),fct4(:),fct5(:) 2391 real(dp),allocatable :: fctii(:),fp(:),fpp(:),fppp(:),nomega(:),ppsig(:) 2392 real(dp),allocatable :: x1(:),x2(:) 2393 2394 ! ********************************************************************************* 2395 2396 if (npti > 12000) then 2397 msg = "Sorry - the interpolator INTRPL is hard coded for maximum 12000 points." // & 2398 & ch10 // " Reduce the conducti input npti, or implement a better interpolator!" 2399 ABI_ERROR(msg) 2400 end if 2401 2402 write(std_out,'(2a)')ch10,'Calculate the principal value and related optical properties' 2403 write(std_out,'(a)')'following W.J. Thomson computer in physics vol 12 p94 1998 for ' 2404 write(std_out,'(a)')'the principal value.' 2405 write(std_out,'(a)')'OPTIONS' 2406 write(std_out,'(a)')'use default number of integration pts: npt=10000' 2407 write(std_out,'(a)')'Use default value for delta interval: del=1e-3' 2408 2409 if (open_file(trim(filnam_out_sig)//'_eps',msg,newunit=eps_unt,status='replace',action="write")/=0) then 2410 ABI_ERROR(msg) 2411 end if 2412 write(eps_unt,'(a)')'#energy (eV),sigma_1(Ohm-1cm-1),sigma_2(Ohm-1cm-1),epsilon_1,epsilon_2' 2413 2414 if (open_file(trim(filnam_out_sig)//'_abs',msg,newunit=abs_unt,status='replace',action="write")/=0) then 2415 ABI_ERROR(msg) 2416 end if 2417 write(abs_unt,'(a)')'#energy(eV),nomega,komega,refl.,abso.(cm-1)' 2418 2419 ABI_MALLOC(fct,(npt)) 2420 ABI_MALLOC(fct2,(npt)) 2421 ABI_MALLOC(fct3,(npt)) 2422 ABI_MALLOC(fct4,(npt)) 2423 ABI_MALLOC(fct5,(npt)) 2424 ABI_MALLOC(fp,(npt)) 2425 ABI_MALLOC(fpp,(npt)) 2426 ABI_MALLOC(fppp,(npt)) 2427 ABI_MALLOC(x1,(npt)) 2428 ABI_MALLOC(x2,(npt)) 2429 ABI_MALLOC(fct1,(npt)) 2430 ABI_MALLOC(ppsig,(npt)) 2431 ABI_MALLOC(fctii,(npt)) 2432 ABI_MALLOC(abso,(npt)) 2433 ABI_MALLOC(nomega,(npt)) 2434 2435 !loop on the initial energy grid 2436 do ip=1,npti 2437 2438 ! adjust the interval before and after the pole to reflect range/npt interval 2439 xsum=zero 2440 dx=(xi(npti)-xi(1))/dble(npt-1) 2441 pole=xi(ip) 2442 npt1=int((pole-del)/dx) 2443 dx1=zero 2444 if(npt1/=1) dx1=(pole-del)/(npt1-1) 2445 npt2=int((xi(npti)-pole-del)/dx) 2446 dx2=(xi(npti)-pole-del)/(npt2-1) 2447 2448 ! for the moment skip the pp calculation when the pole if too close to the end of the range 2449 if (npt1<=1.or.npt2<=1) then 2450 2451 xsum=zero 2452 ppsig(ip)=zero 2453 2454 else 2455 2456 ! define the fct for which the pp calculation is needed using xi^2-pole^2 factorization 2457 fctii(1:npti) = zero 2458 fctii(1:npti)=fcti(1:npti)*pole/(xi(1:npti)+pole) 2459 2460 ! define the grid on each side of the pole x1 before x2 after 2461 do ii=1,npt1 2462 x1(ii)=dx1*dble(ii-1) 2463 end do 2464 do ii=1,npt2 2465 x2(ii)=pole+del+dx2*dble(ii-1) 2466 end do 2467 2468 ! interpolate the initial fct fii on the new grids x1 and x2 (cubic spline) 2469 ! write(std_out,*) npti,npt1 2470 2471 ! MJV 6/12/2008: 2472 ! For each use of fctii should ensure that npt1 npt2 etc... are less than 2473 ! npt=len(fctii) 2474 ! TODO: move to spline/splint routines with no memory limitation 2475 call intrpl(npti,xi,fctii,npt1,x1,fct4,fct1,fct5,1) 2476 call intrpl(npti,xi,fctii,npt2,x2,fct3,fct2,fct5,1) 2477 2478 ! calculate the two integrals from 0-->pole-lamda and pole+lamda--> end range 2479 ! trapezoidal integration 2480 do ii=1,npt1 2481 fct1(ii)=fct4(ii)/(x1(ii)-pole) 2482 end do 2483 do ii=1,npt2 2484 fct2(ii)=fct3(ii)/(x2(ii)-pole) 2485 end do 2486 2487 do ii=2,npt1-1 2488 xsum=xsum+fct1(ii)*dx1 2489 end do 2490 do ii=2,npt2-1 2491 xsum=xsum+fct2(ii)*dx2 2492 end do 2493 xsum=xsum+half*(fct1(1)+fct1(npt1))*dx1+half*(fct2(1)+fct2(npt2))*dx2 2494 2495 ! Calculate the first and third derivative at the pole and add the taylor expansion 2496 ! TODO: move to spline/splint routines with no memory limitation 2497 call intrpl(npti,xi,fctii,npti,xi,fct3,fct4,fct5,1) 2498 call intrpl(npti,xi,fct4,1,(/pole/),fp,fpp,fppp,1) 2499 2500 idel=two*fp(1)*(del)+fppp(1)*(del**3)/nine 2501 xsum=xsum+idel 2502 2503 end if 2504 2505 ! Calculate the derivated optical quantities and output the value 2506 sigma2=(-two/pi)*xsum 2507 eps1=one-(four_pi*sigma2/(pole)) 2508 eps2=four*fcti(ip)*pi/(pole) 2509 2510 ! A special treatment of the case where eps2 is very small compared to eps1 is needed 2511 if(eps2**2 > eps1**2 * tol12)then 2512 nomega(ip)=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2))) 2513 komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2))) 2514 abso(ip)=four_pi*fcti(ip)*ohmtosec*Ohmcm/nomega(ip)/(Sp_Lt_SI*100._dp) 2515 else if(eps1>zero)then 2516 nomega(ip)=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2))) 2517 komega=half*abs(eps2/sqrt(eps1)) 2518 abso(ip)=four_pi*fcti(ip)*ohmtosec*Ohmcm/nomega(ip)/(Sp_Lt_SI*100._dp) 2519 else if(eps1<zero)then 2520 nomega(ip)=half*abs(eps2/sqrt(-eps1)) 2521 komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2))) 2522 abso(ip)=two*sqrt(-eps1)*pole*ohmtosec*Ohmcm/(Sp_Lt_SI*100._dp) 2523 end if 2524 2525 refl=((one-nomega(ip))**2 + komega**2)/ & 2526 & ((one+nomega(ip))**2 + komega**2) 2527 2528 write(eps_unt,'(5e18.10)') Ha_eV*pole,fcti(ip)*Ohmcm,sigma2*Ohmcm,eps1,eps2 2529 write(abs_unt,'(5e18.10)') Ha_eV*pole,nomega(ip),komega,refl,abso(ip) 2530 2531 end do 2532 2533 close(eps_unt) 2534 close(abs_unt) 2535 2536 ABI_FREE(fct) 2537 ABI_FREE(x1) 2538 ABI_FREE(x2) 2539 ABI_FREE(fct2) 2540 ABI_FREE(fct3) 2541 ABI_FREE(fct4) 2542 ABI_FREE(fct5) 2543 ABI_FREE(fp) 2544 ABI_FREE(fpp) 2545 ABI_FREE(fppp) 2546 ABI_FREE(fct1) 2547 ABI_FREE(ppsig) 2548 ABI_FREE(fctii) 2549 ABI_FREE(abso) 2550 ABI_FREE(nomega) 2551 2552 end subroutine msig