TABLE OF CONTENTS


ABINIT/emispec_paw [ Functions ]

[ Top ] [ Functions ]

NAME

 emispec_paw

FUNCTION

 This program computes the elements of the emission spectra
 from the Kubo-Greenwood formula for PAW formalism
 largely inspired from the conducti_core_paw routine

INPUTS

  (main routine)

OUTPUT

  (main routine)

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 windows for computations of sigma
  wtk(nkpt)=weight assigned to each k point.

PARENTS

      conducti

CHILDREN

      hdr_free,hdr_io,hdr_read_from_fname,metric,wffclose,wffopen

SOURCE

1771  subroutine emispec_paw(filnam,filnam_out,mpi_enreg)
1772 
1773 
1774 !This section has been created automatically by the script Abilint (TD).
1775 !Do not modify the following lines by hand.
1776 #undef ABI_FUNC
1777 #define ABI_FUNC 'emispec_paw'
1778 !End of the abilint section
1779 
1780  implicit none
1781 
1782 !Arguments -----------------------------------
1783 !scalars
1784  character(len=fnlen) :: filnam,filnam_out
1785  type(MPI_type),intent(in) :: mpi_enreg
1786 
1787 !Local variables-------------------------------
1788 !scalars
1789  integer,parameter :: master=0
1790  integer :: iomode,atnbr,bantot,bdtot_index
1791  integer :: fform2,headform,iatom,iband,icor,ierr,ikpt
1792  integer :: iom,isppol,l1,mband,me,mom
1793  integer :: natom,nband_k,nkpt,nphicor,nspinor,nsppol,ntypat
1794  integer :: occopt,rdwr,spaceComm,iunt,ems_unt,opt2_unt
1795  real(dp) :: del,ecut,fermie
1796  real(dp) :: omin,omax,dom,oml
1797  real(dp) :: Tatm,tsmear,ucvol
1798  character(len=fnlen) :: filnam2,filnam_gen
1799  character(len=500) :: msg
1800  type(hdr_type) :: hdr
1801  type(wffile_type) :: wff2
1802 !arrays
1803  integer,allocatable :: nband(:),ncor(:),lcor(:)
1804  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1805  real(dp),allocatable :: dhdk2_g(:,:,:)
1806  real(dp),allocatable :: eig0_k(:),eigen0(:)
1807  real(dp),allocatable :: energy_cor(:)
1808  real(dp),allocatable :: occ(:),occ_k(:),oml1(:)
1809  real(dp),allocatable :: psinablapsi2(:,:,:,:,:)
1810  real(dp),allocatable :: sigx(:,:,:),sigx_av(:,:),wtk(:)
1811 
1812 ! *********************************************************************************
1813  ABI_UNUSED(mpi_enreg%paral_kgb)
1814 
1815 !Read output file name
1816 !write(std_out,'(a)')' emispec_paw : Please, give the name of the output file ...'
1817 !read(std_in, '(a)') filnam_out
1818 !write(std_out,'(a)')' The name of the output file is :',trim(filnam_out)
1819 !Read data file
1820  if (open_file(filnam,msg,newunit=iunt,form='formatted',action="read") /=0) then
1821    MSG_ERROR(msg)
1822  end if
1823  rewind(iunt)
1824  read(iunt,*)
1825  read(iunt,'(a)')filnam_gen       ! generic name for the files
1826  filnam2=trim(filnam_gen)//'_OPT2'
1827 !Read size of the frequency range
1828  read(iunt,*) dom,omin,omax,mom,atnbr
1829  close(iunt)
1830 
1831  write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width      =',mom,omin,omax,dom,' Ha'
1832  write(std_out,'(a)')'--------------------------------------------'
1833  write(std_out,'(a,i4)') 'selected atom for X ray emission',atnbr
1834  write(std_out,'(a)')'--------------------------------------------'
1835 
1836 !Open the Wavefunction and optic files
1837 !These default values are typical of sequential use
1838  iomode=IO_MODE_FORTRAN; spaceComm=xmpi_comm_self; me=0
1839 
1840 ! Read the header of the OPT2 file.
1841  call hdr_read_from_fname(hdr, filnam2, fform2, spaceComm)
1842  call hdr_free(hdr)
1843 
1844  if (fform2 /= 611) then
1845    MSG_ERROR("Abinit8 requires an OPT2 file with fform = 611")
1846  end if
1847 
1848 !Open the optic files
1849  opt2_unt = get_unit()
1850  call WffOpen(iomode,spaceComm,filnam2,ierr,wff2,master,me,opt2_unt)
1851 
1852 !Read the header
1853  rdwr=1
1854  call hdr_io(fform2,hdr,rdwr,wff2)
1855 
1856 !Extract info from the header
1857  headform=hdr%headform
1858  bantot=hdr%bantot
1859  ecut=hdr%ecut_eff
1860  natom=hdr%natom
1861  nkpt=hdr%nkpt
1862  nspinor=hdr%nspinor
1863  nsppol=hdr%nsppol
1864  ntypat=hdr%ntypat
1865  occopt=hdr%occopt
1866  rprimd(:,:)=hdr%rprimd(:,:)
1867  ABI_ALLOCATE(nband,(nkpt*nsppol))
1868  ABI_ALLOCATE(occ,(bantot))
1869  ABI_ALLOCATE(wtk,(nkpt))
1870  fermie=hdr%fermie
1871  tsmear=hdr%tsmear
1872  occ(1:bantot)=hdr%occ(1:bantot)
1873  wtk(1:nkpt)=hdr%wtk(1:nkpt)
1874  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
1875 
1876 !Get mband, as the maximum value of nband(nkpt)
1877  mband=maxval(nband(:))
1878 
1879  write(std_out,*)
1880  write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1,1:3)
1881  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(2,1:3)
1882  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(3,1:3)
1883  write(std_out,'(a,i8)')       ' natom             =',natom
1884  write(std_out,'(a,3i8)')      ' nkpt,mband,nsppol        =',nkpt,mband,nsppol
1885  write(std_out, '(a, f10.5,a)' ) ' ecut              =',ecut,' Ha'
1886  write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie            =',fermie,' Ha',fermie*Ha_eV,' eV'
1887  Tatm=tsmear*Ha_K
1888  write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp              =',tsmear,' Ha ',Tatm,' Kelvin'
1889 
1890  ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol))
1891  read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
1892 
1893  write(std_out,'(a)')'--------------------------------------------'
1894  read(opt2_unt) nphicor
1895  write(std_out,'(a,i4)') 'Number of core orbitals nc=',nphicor
1896  ABI_ALLOCATE(ncor,(nphicor))
1897  ABI_ALLOCATE(lcor,(nphicor))
1898  ABI_ALLOCATE(energy_cor,(nphicor))
1899  do icor=1,nphicor
1900    read(opt2_unt) ncor(icor),lcor(icor),energy_cor(icor)
1901    write(std_out,'(a,2i4,f10.5)') ' n, l, Energy (Ha): ',ncor(icor),lcor(icor),energy_cor(icor)
1902  end do
1903  write(std_out,'(a)')'--------------------------------------------'
1904 
1905 !---------------------------------------------------------------------------------
1906 !gmet inversion to get ucvol
1907  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1908 
1909 !---------------------------------------------------------------------------------
1910 !frequency range fixed by the occupation of the valence states
1911  bdtot_index=0
1912  omax=0
1913  do isppol=1,nsppol
1914    do ikpt=1,nkpt
1915      nband_k=nband(ikpt+(isppol-1)*nkpt)
1916      do iband=1,nband_k
1917        if(eigen0(1+bdtot_index)>=omax) omax=eigen0(1+bdtot_index)
1918        bdtot_index=bdtot_index+1
1919      end do
1920    end do
1921  end do
1922  omin=eigen0(1)
1923  del=(omax-omin)/(mom-1)
1924  ABI_ALLOCATE(oml1,(mom))
1925  do iom=1,mom
1926    oml1(iom)=omin+dble(iom-1)*del
1927  end do
1928  write(std_out,*) 'Valence state orbital energies: omin,omax',omin,omax
1929  ABI_ALLOCATE(sigx,(natom,mom,nphicor))
1930  ABI_ALLOCATE(sigx_av,(mom,nphicor))
1931 !---------------------------------------------------------------------------------
1932 !emission X  -------
1933 !
1934  ABI_ALLOCATE(psinablapsi2,(2,3,mband,nphicor,natom))
1935  sigx=zero
1936  sigx_av=zero
1937  bdtot_index = 0
1938 
1939 !LOOP OVER SPINS
1940  do isppol=1,nsppol
1941 
1942 !  LOOP OVER K-POINTS
1943    do ikpt=1,nkpt
1944      nband_k=nband(ikpt+(isppol-1)*nkpt)
1945      ABI_ALLOCATE(eig0_k,(nband_k))
1946      ABI_ALLOCATE(occ_k,(nband_k))
1947      ABI_ALLOCATE(dhdk2_g,(natom,nband_k,nphicor))
1948 
1949      dhdk2_g   = zero
1950      psinablapsi2=zero
1951 
1952 !    eigenvalue for k-point
1953      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
1954 
1955 !    occupation numbers for k-point
1956      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1957 
1958 !    dipole matrix elements
1959      do iatom=1,natom
1960        read(opt2_unt) ((psinablapsi2(1:2,1,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor)
1961        read(opt2_unt) ((psinablapsi2(1:2,2,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor)
1962        read(opt2_unt) ((psinablapsi2(1:2,3,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor)
1963      end do
1964 
1965 !    loops over atom, bands, core states
1966      do iatom=1,natom
1967        do iband=1,nband_k
1968          do icor=1,nphicor
1969 
1970            do l1=1,3
1971              dhdk2_g(iatom,iband,icor)=dhdk2_g(iatom,iband,icor)+( &
1972 &             psinablapsi2(1,l1,iband,icor,iatom)*psinablapsi2(1,l1,iband,icor,iatom) &
1973 &             +psinablapsi2(2,l1,iband,icor,iatom)*psinablapsi2(2,l1,iband,icor,iatom))
1974            end do
1975 
1976 !          emission for each omega
1977            do iom=1,mom
1978              oml=-energy_cor(icor)-oml1(iom)
1979              sigx(iatom,iom,icor)=sigx(iatom,iom,icor)+ wtk(ikpt)*dhdk2_g(iatom,iband,icor)&
1980 &             *occ_k(iband)/oml*dexp(-((-energy_cor(icor)+eig0_k(iband)-oml)/dom)**2)
1981            end do
1982          end do !icor
1983        end do  !iband
1984      end do !iatom
1985      bdtot_index=bdtot_index+nband_k
1986      ABI_DEALLOCATE(eig0_k)
1987      ABI_DEALLOCATE(occ_k)
1988      ABI_DEALLOCATE(dhdk2_g)
1989 !    end loop over k
1990    end do
1991 !  end loop over spins
1992  end do
1993  ABI_DEALLOCATE(psinablapsi2)
1994 
1995  do iatom=1,natom
1996    do icor=1,nphicor
1997      do iom=1,mom
1998        if(sigx(iatom,iom,icor)<=tol16) sigx(iatom,iom,icor)=zero
1999      end do
2000    end do
2001  end do ! iatom
2002 
2003  sigx=sigx*two_pi*third*dble(natom)/(dom*ucvol)*half/dsqrt(pi)
2004 
2005  do icor=1,nphicor
2006    do iom=1,mom
2007      do iatom=1,natom
2008        sigx_av(iom,icor) =sigx_av(iom,icor)+sigx(iatom,iom,icor)/dble(natom)
2009      end do
2010    end do
2011  end do
2012 
2013  if (open_file(trim(filnam_out)//'_emisX',msg,newunit=ems_unt,form='formatted', action="write") /= 0) then
2014    MSG_ERROR(msg)
2015  end if
2016  do iom=1,mom
2017    write(ems_unt,'(9(1x,e15.8))') &
2018 &   ((-energy_cor(icor)+oml1(iom)),sigx_av(iom,icor),sigx(atnbr,iom,icor),icor=1,nphicor)
2019  end do
2020  close(ems_unt)
2021 
2022  call WffClose(wff2,ierr)
2023 
2024  ABI_DEALLOCATE(sigx)
2025  ABI_DEALLOCATE(sigx_av)
2026  ABI_DEALLOCATE(ncor)
2027  ABI_DEALLOCATE(lcor)
2028  ABI_DEALLOCATE(energy_cor)
2029  ABI_DEALLOCATE(eigen0)
2030  ABI_DEALLOCATE(nband)
2031  ABI_DEALLOCATE(oml1)
2032  ABI_DEALLOCATE(occ)
2033  ABI_DEALLOCATE(wtk)
2034 
2035  call hdr_free(hdr)
2036 
2037 end subroutine emispec_paw

ABINIT/m_conducti [ Modules ]

[ Top ] [ 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-2018 ABINIT group (VRecoules, PGhosh, SMazevet, SM, SVinko)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

23 #if defined HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26 
27 #include "abi_common.h"
28 
29 module m_conducti
30 
31  use defs_basis
32  use m_errors
33  use m_abicore
34  use defs_abitypes
35  use m_xmpi
36  use m_wffile
37  use m_wfk
38  use m_hdr
39  use m_nctk
40 
41  use m_io_tools,     only : open_file, get_unit
42  use m_fstrings,     only : sjoin
43  use m_symtk,        only : matr3inv
44  use m_hide_lapack,  only : jacobi
45  use m_occ,          only : getnel
46  use m_geometry,     only : metric
47  use m_splines,      only : intrpl
48 
49  implicit none
50 
51  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

PARENTS

      conducti

CHILDREN

      getnel,hdr_free,jacobi,matr3inv,metric,msig,nctk_fort_or_ncfile
      wfk_close,wfk_open_read,wfk_read_eigk

SOURCE

 901 subroutine conducti_nc(filnam,filnam_out,mpi_enreg)
 902 
 903 
 904 !This section has been created automatically by the script Abilint (TD).
 905 !Do not modify the following lines by hand.
 906 #undef ABI_FUNC
 907 #define ABI_FUNC 'conducti_nc'
 908 !End of the abilint section
 909 
 910  implicit none
 911 
 912 !Arguments -----------------------------------
 913 !scalars
 914  character(len=fnlen) :: filnam,filnam_out
 915  type(MPI_type),intent(in) :: mpi_enreg
 916 
 917 !Local variables-------------------------------
 918 !scalars
 919  integer,parameter :: formeig0=0,formeig1=1
 920  integer :: bantot,bd2tot_index,bdtot0_index,bdtot_index
 921  integer :: headform,iband,ii,jj,ikpt,iunt
 922  integer :: index_1,iom,isppol,jband,l1,l2,mband,mom,natom,nband1
 923  integer :: nrot,iomode
 924  integer :: nband_k,nkpt,nlign,nrest,nspinor,nsppol,ntypat
 925  integer :: occopt,comm
 926  integer :: tens_unt,lij_unt,sig_unt,kth_unt,ocond_unt
 927  real(dp) :: deltae,dosdeltae,diff_occ,dom,ecut,entropy,fermie,maxocc
 928  real(dp) :: nelect,np_sum,np_sum_k1,np_sum_k2,omin,oml,socc,socc_k,sig
 929  real(dp) :: tphysel,tsmear,ucvol,wind,Tatm
 930  character(len=fnlen) :: filnam0,filnam1,filnam2,filnam3
 931  character(len=500) :: msg
 932  type(hdr_type) :: hdr
 933  type(wfk_t) :: gswfk,ddk1,ddk2,ddk3
 934 !arrays
 935  integer,allocatable :: nband(:)
 936  real(dp) :: gmet(3,3),gmet_inv(3,3),gprimd(3,3),gprimd_inv(3,3),rmet(3,3),rprimd(3,3)
 937  real(dp),allocatable :: cond_kg(:,:,:),cond_kg_cart(:,:,:),cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:)
 938  real(dp),allocatable :: doccde(:),doccde_k(:),cond_kg_xx(:),cond_kg_yy(:),cond_kg_zz(:),trace(:)
 939  real(dp),allocatable :: eig0_k(:),eig0tmp(:),eig1_k(:,:),eigen0(:),eigen11(:)
 940  real(dp),allocatable :: eigen12(:),eigtmp(:)
 941  real(dp),allocatable :: eigen13(:),occ(:),occ_k(:),wtk(:),cond_tot(:),oml1(:)
 942  real(dp),allocatable :: kin11(:),kin12(:),kin21(:),kin22(:)
 943  real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:)
 944  real(dp) :: cond_kg_w(3,3),z(3,3)
 945  real(dp) :: eig_cond(3)
 946 
 947 ! *********************************************************************************
 948 
 949  ABI_UNUSED(mpi_enreg%paral_kgb)
 950 
 951 !Read data file
 952  if (open_file(filnam,msg,newunit=iunt,form='formatted', status="old") /= 0) then
 953    MSG_ERROR(msg)
 954  end if
 955 
 956  rewind(iunt)
 957  read(iunt,*)
 958  read(iunt,'(a)')filnam1       ! first ddk file
 959  read(iunt,'(a)')filnam2       ! second ddk file
 960  read(iunt,'(a)')filnam3       ! third ddk file
 961  read(iunt,'(a)')filnam0       ! ground-state data
 962 
 963 !Open the GS Wavefunction file and the 3 DDK files.
 964 
 965 ! TODO: one should perform basic consistency tests for the GS WFK and the DDK files, e.g.
 966 ! k-points and their order, spins, number of bands could differ in the four files.
 967 ! Note indeed that we are assuming the same numer of bands in all the files.
 968  comm = xmpi_comm_self
 969  call nctk_fort_or_ncfile(filnam0, iomode, msg)
 970  if (len_trim(msg) /= 0) MSG_ERROR(msg)
 971  call wfk_open_read(gswfk,filnam0, formeig0, iomode, get_unit(), comm)
 972 
 973  call nctk_fort_or_ncfile(filnam1, iomode, msg)
 974  if (len_trim(msg) /= 0) MSG_ERROR(msg)
 975  call wfk_open_read(ddk1,filnam1, formeig1, iomode, get_unit(), comm, hdr_out=hdr)
 976 
 977  call nctk_fort_or_ncfile(filnam2, iomode, msg)
 978  if (len_trim(msg) /= 0) MSG_ERROR(msg)
 979  call wfk_open_read(ddk2,filnam2, formeig1, iomode, get_unit(), comm)
 980 
 981  call nctk_fort_or_ncfile(filnam3, iomode, msg)
 982  if (len_trim(msg) /= 0) MSG_ERROR(msg)
 983  call wfk_open_read(ddk3,filnam3, formeig1, iomode, get_unit(), comm)
 984 
 985  if (wfk_compare(ddk1, ddk2) /= 0) then
 986    MSG_ERROR("ddk1 and ddk2 are not consistent. see above messages")
 987  end if
 988  if (wfk_compare(ddk1, ddk3) /= 0) then
 989    MSG_ERROR("ddk1 and ddk3 are not consistent. see above messages")
 990  end if
 991 
 992 !Extract params from the header of the first ddk file (might have been the GS file ?)
 993 
 994 !Extract info from the header
 995  headform=hdr%headform
 996  bantot=hdr%bantot
 997  ecut=hdr%ecut_eff
 998  natom=hdr%natom
 999  nkpt=hdr%nkpt
1000  nspinor=hdr%nspinor
1001  nsppol=hdr%nsppol
1002  ntypat=hdr%ntypat
1003  occopt=hdr%occopt
1004  rprimd(:,:)=hdr%rprimd(:,:)
1005  ABI_ALLOCATE(nband,(nkpt*nsppol))
1006  ABI_ALLOCATE(occ,(bantot))
1007  fermie=hdr%fermie
1008  occ(1:bantot)=hdr%occ(1:bantot)
1009  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
1010 
1011 !Get mband, as the maximum value of nband(nkpt)
1012  mband=maxval(nband(:))
1013 
1014  write(std_out,*)
1015  write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1:3,1)
1016  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,2)
1017  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,3)
1018  write(std_out,'(a,i8)')       ' natom             =',natom
1019  write(std_out,'(a,2i8)')      ' nkpt,mband        =',nkpt,mband
1020  write(std_out,'(a, f10.5,a)' ) ' ecut              =',ecut,' Ha'
1021  write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie            =',fermie,' Ha',fermie*Ha_eV,' eV'
1022 
1023 !Prepare the reading of ddk Wff files
1024  ABI_ALLOCATE(eigtmp,(2*mband*mband))
1025  ABI_ALLOCATE(eig0tmp,(mband))
1026 
1027 !Read the eigenvalues of ground-state and ddk files
1028  ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol))
1029  ABI_ALLOCATE(eigen11,(2*mband*mband*nkpt*nsppol))
1030  ABI_ALLOCATE(eigen12,(2*mband*mband*nkpt*nsppol))
1031  ABI_ALLOCATE(eigen13,(2*mband*mband*nkpt*nsppol))
1032  bdtot0_index=0 ; bdtot_index=0
1033  do isppol=1,nsppol
1034    do ikpt=1,nkpt
1035      nband1=nband(ikpt+(isppol-1)*nkpt)
1036      call wfk_read_eigk(gswfk,ikpt,isppol,xmpio_single,eig0tmp)
1037      eigen0(1+bdtot0_index:nband1+bdtot0_index)=eig0tmp(1:nband1)
1038 
1039      call wfk_read_eigk(ddk1,ikpt,isppol,xmpio_single,eigtmp)
1040      eigen11(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
1041 
1042      call wfk_read_eigk(ddk2,ikpt,isppol,xmpio_single,eigtmp)
1043      eigen12(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
1044 
1045      call wfk_read_eigk(ddk3,ikpt,isppol,xmpio_single,eigtmp)
1046      eigen13(1+bdtot_index:2*nband1**2+bdtot_index)=eigtmp(1:2*nband1**2)
1047 
1048      bdtot0_index=bdtot0_index+nband1
1049      bdtot_index=bdtot_index+2*nband1**2
1050    end do
1051  end do
1052 
1053  ! Close files
1054  call wfk_close(gswfk)
1055  call wfk_close(ddk1)
1056  call wfk_close(ddk2)
1057  call wfk_close(ddk3)
1058 
1059  ABI_DEALLOCATE(eigtmp)
1060  ABI_DEALLOCATE(eig0tmp)
1061 
1062 !---------------------------------------------------------------------------------
1063 !gmet inversion
1064  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1065  call matr3inv(gmet,gmet_inv)
1066  call matr3inv(gprimd,gprimd_inv)
1067 
1068 !---------------------------------------------------------------------------------
1069 !derivative of occupation wrt the energy.
1070  ABI_ALLOCATE(doccde,(mband*nkpt*nsppol))
1071  ABI_ALLOCATE(wtk,(nkpt))
1072 
1073  read(iunt,*)tsmear
1074  Tatm=tsmear*Ha_K
1075  write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp              =',tsmear,' Ha ',Tatm,' Kelvin'
1076 !
1077  nlign=nkpt/6
1078  nrest=nkpt-6*nlign
1079  index_1=0
1080  do ii=1,nlign
1081    read(iunt,*)wtk(1+index_1:6+index_1)
1082    index_1=index_1+6
1083  end do
1084  if (nrest/=0) then
1085    read(iunt,*)wtk(6*nlign+1:nkpt)
1086  end if
1087 !
1088  if (occopt==1) then
1089    write(std_out,'(a,i4)')  ' occopt            =',occopt
1090    doccde=0.0d0
1091  else
1092    tphysel=zero
1093    maxocc=two/(nsppol*nspinor)
1094    dosdeltae=zero
1095    call getnel(doccde,dosdeltae,eigen0,entropy,fermie,maxocc,mband,nband,&
1096 &   nelect,nkpt,nsppol,occ,occopt,1,tphysel,tsmear,11,wtk)
1097 !  DEBUG
1098 !  write(std_out,'(a,f10.5)')' getnel : nelect   =',nelect
1099 !  ENDDEBUG
1100  end if
1101 !---------------------------------------------------------------------------------
1102 !size of the frequency range
1103  read(iunt,*)dom,wind
1104  close(iunt)
1105  mom=dint(wind/dom)
1106  ABI_ALLOCATE(oml1,(mom))
1107  do iom=1,mom
1108    oml1(iom)=tol10*1000.0d0+dble(iom)*dom
1109  end do
1110 
1111  ABI_ALLOCATE(cond_nd,(mom,3,3))
1112  ABI_ALLOCATE(cond_kg,(mom,3,3))
1113  ABI_ALLOCATE(cond_kg_cart,(mom,3,3))
1114  ABI_ALLOCATE(cond_kg_xx,(mom))
1115  ABI_ALLOCATE(cond_kg_yy,(mom))
1116  ABI_ALLOCATE(trace,(mom))
1117  ABI_ALLOCATE(cond_kg_zz,(mom))
1118  ABI_ALLOCATE(cond_tot,(mom))
1119  ABI_ALLOCATE(kin11,(mom))
1120  ABI_ALLOCATE(kin12,(mom))
1121  ABI_ALLOCATE(kin21,(mom))
1122  ABI_ALLOCATE(kin22,(mom))
1123  ABI_ALLOCATE(kin11_k,(mom))
1124  ABI_ALLOCATE(kin12_k,(mom))
1125  ABI_ALLOCATE(kin21_k,(mom))
1126  ABI_ALLOCATE(kin22_k,(mom))
1127  ABI_ALLOCATE(Kth,(mom))
1128  ABI_ALLOCATE(Stp,(mom))
1129  write(std_out,'(a,i8,2f10.5,a)')' mom,wind,dom      =',mom,wind,dom,' Ha'
1130 
1131 !---------------------------------------------------------------------------------
1132 
1133  kin11   = 0.0d0
1134  kin12   = 0.0d0
1135  kin21   = 0.0d0
1136  kin22   = 0.0d0
1137  np_sum  = 0.0d0
1138  socc    = 0.0d0
1139  cond_kg = 0.0d0
1140 
1141 
1142 !LOOP OVER SPINS
1143  do isppol=1,nsppol
1144 !
1145    bdtot_index = 0
1146    bd2tot_index = 0
1147 !
1148    deltae  = 0.0d0
1149 !
1150 !  BIG FAT k POINT LOOP
1151 !
1152    do ikpt=1,nkpt
1153 !
1154      nband_k=nband(ikpt+(isppol-1)*nkpt)
1155 !
1156      ABI_ALLOCATE(eig0_k,(nband_k))
1157      ABI_ALLOCATE(eig1_k,(2*nband_k**2,3))
1158      ABI_ALLOCATE(occ_k,(nband_k))
1159      ABI_ALLOCATE(doccde_k,(nband_k))
1160      ABI_ALLOCATE(dhdk2_r,(nband_k,nband_k,3,3))
1161      ABI_ALLOCATE(dhdk2_g,(nband_k,nband_k))
1162 
1163      cond_nd   = 0.0d0
1164      kin11_k   = 0.0d0
1165      kin12_k   = 0.0d0
1166      kin21_k   = 0.0d0
1167      kin22_k   = 0.0d0
1168      np_sum_k1 = 0.0d0
1169      np_sum_k2 = 0.0d0
1170      socc_k    = 0.0d0
1171      dhdk2_r   = 0.0d0
1172      dhdk2_g   = 0.0d0
1173 !
1174 !    eigenvalue for k-point
1175      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
1176 !    first derivative eigenvalues for k-point
1177      eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_k**2+bd2tot_index)
1178      eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_k**2+bd2tot_index)
1179      eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_k**2+bd2tot_index)
1180 !    occupation numbers for k-point
1181      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1182 !    derivative of occupation number for k-point
1183      doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index)
1184 !
1185 !    DEBUG
1186 !    write(16,*)
1187 !    write(16,*)' conducti : ikpt=',ikpt
1188 !    do iband=1,nband_k
1189 !    write(16, '(i4,4es22.12)' )iband,wtk(ikpt),occ_k(iband),&
1190 !    &                            doccde_k(iband),eig0_k(iband)
1191 !    end do
1192 !    write(16,*)
1193 !    ENDDEBUG
1194 !
1195 !    LOOP OVER BAND
1196      do iband=1,nband_k
1197        do jband=1,nband_k
1198 !
1199          do l1=1,3
1200            do l2=1,3
1201              do ii=1,3
1202                do jj=1,3
1203                  dhdk2_r(iband,jband,l1,l2)=dhdk2_r(iband,jband,l1,l2)+(rprimd(l1,ii)&
1204 &                 *eig1_k(2*iband-1+(jband-1)*2*nband_k,ii)*&
1205 &                 rprimd(l2,jj)*eig1_k(2*iband-1+(jband-1)*2*nband_k,jj)&
1206 &                 +rprimd(l1,ii)*eig1_k(2*iband  +(jband-1)*2*nband_k,ii)*&
1207 &                 rprimd(l2,jj)*eig1_k(2*iband+(jband-1)*2*nband_k,jj))
1208                end do
1209              end do
1210            end do
1211          end do
1212 
1213          do l1=1,3
1214            do l2=1,3
1215              dhdk2_r(iband,jband,l1,l2)=dhdk2_r(iband,jband,l1,l2)/two_pi/two_pi
1216            end do
1217          end do
1218 !
1219          do l1=1,3
1220            do l2=1,3
1221              dhdk2_g(iband,jband)=dhdk2_g(iband,jband)+gmet_inv(l1,l2)*( &
1222 &             eig1_k(2*iband-1+(jband-1)*2*nband_k,l1)*&
1223 &             eig1_k(2*iband-1+(jband-1)*2*nband_k,l2) &
1224 &             +eig1_k(2*iband  +(jband-1)*2*nband_k,l1)*&
1225 &             eig1_k(2*iband  +(jband-1)*2*nband_k,l2))
1226            end do
1227          end do
1228          dhdk2_g(iband,jband)=dhdk2_g(iband,jband)/two_pi/two_pi
1229 !
1230          diff_occ = occ_k(iband)-occ_k(jband)
1231 !        if (dabs(diff_occ)>=tol8) then
1232 !
1233 !        Conductivity for each omega
1234          omin = 0.0d0
1235          do iom=1,mom
1236            oml=oml1(iom)
1237            if (jband>iband) then
1238              sig= dhdk2_g(iband,jband)&
1239 &             *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)&
1240 &             -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2))
1241              kin11_k(iom)=kin11_k(iom)+sig
1242              kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie)
1243              kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie)
1244              kin22_k(iom)=kin22_k(iom) + &
1245 &             sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie)
1246            end if
1247            do l1=1,3
1248              do l2=1,3
1249                cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(iband,jband,l1,l2)&
1250 &               *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)
1251              end do
1252            end do
1253 
1254          end do
1255 !
1256 !        Sumrule start
1257          if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then
1258            np_sum_k1=np_sum_k1 -dhdk2_g(iband,jband)&
1259 &           *(diff_occ)/(eig0_k(iband)-eig0_k(jband))
1260          else
1261            np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(iband,jband)
1262          end if
1263 !
1264 
1265 !        end loop over band
1266 !        end if
1267        end do
1268        socc_k=socc_k+occ_k(iband)
1269      end do
1270 !
1271      do iom=1,mom
1272        kin11(iom)=kin11(iom)+wtk(ikpt)*kin11_k(iom)
1273        kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom)
1274        kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom)
1275        kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom)
1276        do l1=1,3
1277          do l2=1,3
1278            cond_kg(iom,l1,l2)=cond_kg(iom,l1,l2)+wtk(ikpt)*cond_nd(iom,l1,l2)
1279          end do
1280        end do
1281      end do
1282 
1283      np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2)
1284      socc=socc+wtk(ikpt)*socc_k
1285 !
1286 !    validity limit
1287      deltae=deltae+(eig0_k(nband_k)-fermie)
1288 
1289      bd2tot_index=bd2tot_index+2*nband_k**2
1290      bdtot_index=bdtot_index+nband_k
1291      ABI_DEALLOCATE(eig0_k)
1292      ABI_DEALLOCATE(eig1_k)
1293      ABI_DEALLOCATE(occ_k)
1294      ABI_DEALLOCATE(doccde_k)
1295      ABI_DEALLOCATE(dhdk2_r)
1296      ABI_DEALLOCATE(dhdk2_g)
1297 !    End loop over k
1298    end do
1299 
1300    write(std_out,'(a,3f10.5)')' sumrule           =',np_sum/socc/three,socc
1301    write(std_out,'(a,f10.5,a,f10.5,a)')&
1302 &   ' Emax-Efermi       =',deltae/dble(nkpt),' Ha',deltae/dble(nkpt)*Ha_eV,' eV'
1303 
1304 !  End loop over spins
1305  end do
1306 
1307  cond_kg=cond_kg*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
1308 
1309 
1310 !Check that new output file does NOT exist
1311 !Keep this line : prevent silly (compiler ?) bug on HP 8000
1312  write(std_out,*)' conducti : call isfile '
1313 !
1314  if (open_file(trim(filnam_out)//'_tens',msg,newunit=tens_unt,form='formatted',action="write") /= 0) then
1315    MSG_ERROR(msg)
1316  end if
1317  if (open_file(trim(filnam_out)//'_Lij',msg,newunit=lij_unt,form='formatted',action="write") /= 0) then
1318    MSG_ERROR(msg)
1319  end if
1320  write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22'
1321 
1322  if (open_file(trim(filnam_out)//'_sig',msg,newunit=sig_unt,form='formatted',action="write") /= 0) then
1323    MSG_ERROR(msg)
1324  end if
1325  write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV)    cond(ua)             cond(ohm.cm)-1'
1326 
1327  if (open_file(trim(filnam_out)//'_Kth',msg,newunit=kth_unt,form='formatted',action="write") /= 0) then
1328    MSG_ERROR(msg)
1329  end if
1330  write(kth_unt,'(a)')&
1331 & ' #omega(ua) hbar*omega(eV)  thermal cond(ua)   Kth(W/m/K)   thermopower(ua)   Stp(microohm/K)'
1332 
1333  if (open_file(trim(filnam_out)//'.out',msg,newunit=ocond_unt,form='formatted',action="write") /= 0) then
1334    MSG_ERROR(msg)
1335  end if
1336  write(ocond_unt,'(a)' )' Conducti output file:'
1337  write(ocond_unt,'(a)' )' Contains all results produced by conducti utility'
1338  write(ocond_unt,'(a)' )' '
1339  write(ocond_unt,'(a)')' # omega(ua)       cond(ua)             thermal cond(ua)       thermopower(ua)'
1340 !
1341 !call isfile(filnam_out,'new')
1342 
1343 !Keep this line : prevent silly (compiler ?) bug on HP 8000
1344  write(std_out,*)' conducti : after call isfile '
1345 !
1346 !Compute thermal conductivity and thermopower
1347  do iom=1,mom
1348    oml=oml1(iom)
1349    kin11(iom)=kin11(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
1350    kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
1351    kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
1352    kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
1353    if (dabs(kin11(iom))<10.0d-20) kin11(iom)=zero
1354    Kth(iom)=kin22(iom)
1355    Stp(iom)=zero
1356    if(kin11(iom)/=zero)  then
1357      Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/kin11(iom))
1358      Stp(iom)=kin12(iom)/(kin11(iom)*Tatm)
1359    end if
1360    if (dabs(Kth(iom))<10.0d-20) Kth(iom)=0.0d0
1361    if (dabs(Stp(iom))<10.0d-20) Stp(iom)=0.0d0
1362    if (dabs(kin12(iom))<10.0d-20) kin12(iom)=zero
1363    if (dabs(kin21(iom))<10.0d-20) kin21(iom)=zero
1364    if (dabs(kin22(iom))<10.0d-20) kin22(iom)=zero
1365 
1366    write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9
1367    write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,kin11(iom),kin11(iom)*Ohmcm
1368    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
1369    write(ocond_unt,'(1f12.5,3es22.12)') oml,kin11(iom),Kth(iom),Stp(iom)
1370  end do
1371 !
1372 
1373  write(tens_unt,'(a)' )' Conductivity file '
1374  write(tens_unt,'(a)' )' ----------------- '
1375  write(tens_unt,'(a)' )' Contain first the full conductivity tensor, for the desired set of energies,'
1376  write(tens_unt,'(a)' )' then, the three principal values, for the desired set of energies'
1377  write(tens_unt,'(a)' )' (note that eigenvalues are not directly associated with xx,yy,zz)'
1378  write(tens_unt,'(a)' )' '
1379 
1380  write(ocond_unt,'(a)' )' '
1381  write(ocond_unt,'(a)' )' full conductivity tensor, for the desired set of energies'
1382  write(ocond_unt,'(a)' )' then, the three principal values, for the desired set of energies:'
1383 
1384  do iom=1,mom
1385    oml=oml1(iom)*Ha_eV
1386    write(tens_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :'
1387    write(ocond_unt, '(a,es16.6,a)' ) ' energy (in eV) =',oml,', conductivity tensor (in Ohm.cm-1) follows :'
1388    do l1=1,3
1389      write(tens_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3)
1390      write(ocond_unt,"(3f25.15)") (cond_kg(iom,l1,l2)*Ohmcm,l2=1,3)
1391    end do
1392  end do
1393 
1394 !Diagonalizing the conductivity matrix for sigma_xx,sigma_yy,sigma_zz
1395  cond_kg_xx=0d0
1396  cond_kg_yy=0d0
1397  cond_kg_zz=0d0
1398 !trace=0d0    used for checking with the original version of the code
1399  do iom=1,mom
1400    oml=oml1(iom)*Ha_eV
1401    cond_kg_w=0d0
1402    do l1=1,3
1403      do l2=1,3
1404        cond_kg_w(l1,l2)=cond_kg(iom,l1,l2)
1405      end do
1406    end do
1407    call jacobi(cond_kg_w,3,3,eig_cond,z,nrot)
1408 
1409 !  When the value is too small, set it to zero before printing
1410    if(abs(eig_cond(1))<tol10)eig_cond(1)=zero
1411    if(abs(eig_cond(2))<tol10)eig_cond(2)=zero
1412    if(abs(eig_cond(3))<tol10)eig_cond(3)=zero
1413 
1414    cond_kg_xx(iom)=eig_cond(1)
1415    cond_kg_yy(iom)=eig_cond(2)
1416    cond_kg_zz(iom)=eig_cond(3)
1417 !  trace(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom)
1418  end do
1419 
1420 !DEBUG Keep this line : prevent silly (compiler ?) bug on HP 8000
1421 !write(std_out,*)' conducti : after open '
1422 !ENDDEBUG
1423 
1424  write(tens_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.'
1425  write(tens_unt,'(a)')' '
1426  write(tens_unt,'(a)')' #omega(ua)   cond_1(ua)     cond_2(ua) cond_3(ua)  cond_tot(ua)'
1427 
1428  write(ocond_unt,'(a)')' '
1429  write(ocond_unt,'(a,a)')ch10,' Now, print principal values of the conductivity tensor.'
1430  write(ocond_unt,'(a)')' '
1431  write(ocond_unt,'(a)')' #omega(ua)   cond_1(ua)     cond_2(ua) cond_3(ua)  cond_tot(ua)'
1432 
1433 
1434  do iom=1,mom
1435    cond_tot(iom)=cond_kg_xx(iom)+cond_kg_yy(iom)+cond_kg_zz(iom)
1436    write(tens_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
1437    write(ocond_unt,'(f12.5,4es22.12)')oml1(iom),cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
1438  end do
1439 
1440  write(tens_unt,*)
1441  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'
1442  write(ocond_unt,*)
1443  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'
1444 
1445  do iom=1,mom
1446    oml=oml1(iom)*Ha_eV
1447    cond_tot(iom)=cond_tot(iom)*Ohmcm
1448    cond_kg_xx(iom)=cond_kg_xx(iom)*Ohmcm
1449    cond_kg_yy(iom)=cond_kg_yy(iom)*Ohmcm
1450    cond_kg_zz(iom)=cond_kg_zz(iom)*Ohmcm
1451    write(tens_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
1452    write(ocond_unt,'(f12.5,4es22.12)')oml,cond_kg_xx(iom),cond_kg_yy(iom),cond_kg_zz(iom),cond_tot(iom)
1453  end do
1454 !Calculate the imaginary part of the conductivity (principal value)
1455 !+derived optical properties.
1456 
1457  call msig (kin11,mom,oml1,filnam_out)
1458 
1459  close(tens_unt)
1460  close(lij_unt)
1461  close(sig_unt)
1462  close(kth_unt)
1463  close(ocond_unt)
1464 
1465  ABI_DEALLOCATE(nband)
1466  ABI_DEALLOCATE(oml1)
1467  ABI_DEALLOCATE(occ)
1468  ABI_DEALLOCATE(eigen11)
1469  ABI_DEALLOCATE(eigen12)
1470  ABI_DEALLOCATE(eigen13)
1471  ABI_DEALLOCATE(eigen0)
1472  ABI_DEALLOCATE(doccde)
1473  ABI_DEALLOCATE(wtk)
1474  ABI_DEALLOCATE(cond_nd)
1475  ABI_DEALLOCATE(cond_kg)
1476  ABI_DEALLOCATE(cond_kg_cart)
1477  ABI_DEALLOCATE(cond_kg_xx)
1478  ABI_DEALLOCATE(cond_kg_yy)
1479  ABI_DEALLOCATE(trace)
1480  ABI_DEALLOCATE(cond_kg_zz)
1481  ABI_DEALLOCATE(cond_tot)
1482  ABI_DEALLOCATE(kin11)
1483  ABI_DEALLOCATE(kin22)
1484  ABI_DEALLOCATE(kin12)
1485  ABI_DEALLOCATE(kin21)
1486  ABI_DEALLOCATE(kin11_k)
1487  ABI_DEALLOCATE(kin22_k)
1488  ABI_DEALLOCATE(kin12_k)
1489  ABI_DEALLOCATE(kin21_k)
1490  ABI_DEALLOCATE(Stp)
1491  ABI_DEALLOCATE(Kth)
1492 
1493  call hdr_free(hdr)
1494 
1495  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

  (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,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

PARENTS

      conducti

CHILDREN

      hdr_free,hdr_io,hdr_read_from_fname,metric,msig,wffclose,wffopen

SOURCE

128  subroutine conducti_paw(filnam,filnam_out,mpi_enreg)
129 
130 
131 !This section has been created automatically by the script Abilint (TD).
132 !Do not modify the following lines by hand.
133 #undef ABI_FUNC
134 #define ABI_FUNC 'conducti_paw'
135 !End of the abilint section
136 
137  implicit none
138 
139 !Arguments -----------------------------------
140 !scalars
141  character(len=fnlen) :: filnam,filnam_out
142  type(MPI_type),intent(in) :: mpi_enreg
143 
144 !Local variables-------------------------------
145 !scalars
146  integer,parameter :: master=0
147  integer :: iomode,bantot,bdtot_index
148  integer :: fform1,headform,iband,ierr,ikpt
149  integer :: iom,isppol,jband,l1,l2,mband,me,mom
150  integer :: natom,nband_k,nkpt,nspinor,nsppol,ntypat
151  integer :: occopt,rdwr,spaceComm,iunt,opt_unt
152  integer :: lij_unt,sig_unt,kth_unt,ocond_unt
153  real(dp) :: del,deltae,diff_occ,ecut,fermie,maxocc
154  real(dp) :: np_sum,np_sum_k1,np_sum_k2,omin,omax,dom,oml,sig,socc,socc_k
155  real(dp) :: Tatm,tphysel,tsmear,ucvol
156  character(len=fnlen) :: filnam1,filnam_gen
157  character(len=500) :: msg
158  type(hdr_type) :: hdr
159  type(wffile_type) :: wff1
160 !arrays
161  integer,allocatable :: nband(:)
162  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
163  real(dp),allocatable :: cond_nd(:,:,:),dhdk2_r(:,:,:,:),dhdk2_g(:,:,:)
164  real(dp),allocatable :: doccde(:),doccde_k(:),eig0_k(:),eigen0(:)
165  real(dp),allocatable :: occ(:),occ_k(:),wtk(:),oml1(:)
166  real(dp),allocatable :: kin11(:,:),kin12(:),kin21(:),kin22(:)
167  real(dp),allocatable :: kin11_k(:),kin12_k(:),kin21_k(:),kin22_k(:),Kth(:),Stp(:)
168  real(dp),allocatable :: psinablapsi(:,:,:,:),sig_abs(:)
169 
170 ! *********************************************************************************
171  ABI_UNUSED(mpi_enreg%paral_kgb)
172 
173 !write(std_out,'(a)')' The name of the output file is :',trim(filnam_out)
174 !Read data file
175  if (open_file(filnam, msg, newunit=iunt, form='formatted', action="read", status="old") /=0 ) then
176    MSG_ERROR(msg)
177  end if
178  rewind(iunt)
179  read(iunt,*)
180  read(iunt,'(a)')filnam_gen       ! generic name for the files
181  filnam1=trim(filnam_gen)//'_OPT'
182 !Read size of the frequency range
183  read(iunt,*) dom,omin,omax,mom
184  close(iunt)
185  write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width      =',mom,omin,omax,dom,' Ha'
186 
187 !These default values are typical of sequential use
188  iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0
189 
190 ! Read the header of the optic files
191  call hdr_read_from_fname(hdr, filnam1, fform1, spaceComm)
192  call hdr_free(hdr)
193  if (fform1 /= 610) then
194    MSG_ERROR("Abinit8 requires an OPT file with fform = 610")
195  end if
196 
197 !Open the conducti and/or optic files
198  opt_unt = get_unit()
199  call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,opt_unt)
200 
201 !Read the header from Ground state file
202  rdwr=1
203  call hdr_io(fform1,hdr,rdwr,wff1)
204  ABI_CHECK(fform1/=0,"Error opening wff1")
205 
206 !Extract info from the header
207  headform=hdr%headform
208  bantot=hdr%bantot
209  ecut=hdr%ecut_eff
210  natom=hdr%natom
211  nkpt=hdr%nkpt
212  nspinor=hdr%nspinor
213  nsppol=hdr%nsppol
214  ntypat=hdr%ntypat
215  occopt=hdr%occopt
216  rprimd(:,:)=hdr%rprimd(:,:)
217  ABI_ALLOCATE(nband,(nkpt*nsppol))
218  ABI_ALLOCATE(occ,(bantot))
219  ABI_ALLOCATE(wtk,(nkpt))
220  fermie=hdr%fermie
221  tsmear=hdr%tsmear
222  occ(1:bantot)=hdr%occ(1:bantot)
223  wtk(1:nkpt)=hdr%wtk(1:nkpt)
224  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
225 
226 !Get mband, as the maximum value of nband(nkpt)
227  mband=maxval(nband(:))
228 
229  write(std_out,*)
230  write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1:3,1)
231  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,2)
232  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,3)
233  write(std_out,'(a,i8)')       ' natom             =',natom
234  write(std_out,'(a,3i8)')      ' nkpt,mband,nsppol        =',nkpt,mband,nsppol
235  write(std_out, '(a, f10.5,a)' ) ' ecut              =',ecut,' Ha'
236  write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie            =',fermie,' Ha',fermie*Ha_eV,' eV'
237  Tatm=tsmear*Ha_K
238  write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp              =',tsmear,' Ha ',Tatm,' Kelvin'
239 
240  ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol))
241  read(opt_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
242 !
243 !
244 !---------------------------------------------------------------------------------
245 !gmet inversion to get ucvol
246  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
247 
248 !---------------------------------------------------------------------------------
249 !derivative of occupation wrt the energy.
250  ABI_ALLOCATE(doccde,(mband*nkpt*nsppol))
251  if (occopt==1) then
252    write(std_out,'(a,i4)')  ' occopt            =',occopt
253    doccde=zero
254  else
255    tphysel=zero
256    maxocc=two/(nsppol*nspinor)
257  end if
258 !---------------------------------------------------------------------------------
259 !size of the frequency range
260  del=(omax-omin)/(mom-1)
261  ABI_ALLOCATE(oml1,(mom))
262  do iom=1,mom
263    oml1(iom)=omin+dble(iom-1)*del
264  end do
265  ABI_ALLOCATE(cond_nd,(mom,3,3))
266  ABI_ALLOCATE(kin11,(mom,nsppol))
267  ABI_ALLOCATE(kin12,(mom))
268  ABI_ALLOCATE(kin21,(mom))
269  ABI_ALLOCATE(kin22,(mom))
270  ABI_ALLOCATE(sig_abs,(mom))
271  ABI_ALLOCATE(kin11_k,(mom))
272  ABI_ALLOCATE(kin12_k,(mom))
273  ABI_ALLOCATE(kin21_k,(mom))
274  ABI_ALLOCATE(kin22_k,(mom))
275  ABI_ALLOCATE(Kth,(mom))
276  ABI_ALLOCATE(Stp,(mom))
277 
278 !---------------------------------------------------------------------------------
279 !Conductivity -------
280 !
281  ABI_ALLOCATE(psinablapsi,(2,3,mband,mband))
282  kin11   = zero
283  kin12   = zero
284  kin21   = zero
285  kin22   = zero
286  np_sum  = zero
287  socc    = zero
288  sig_abs = zero
289 
290  bdtot_index = 0
291 
292 !LOOP OVER SPINS/K
293  deltae  = zero
294  do isppol=1,nsppol
295    do ikpt=1,nkpt
296      nband_k=nband(ikpt+(isppol-1)*nkpt)
297      ABI_ALLOCATE(eig0_k,(nband_k))
298      ABI_ALLOCATE(occ_k,(nband_k))
299      ABI_ALLOCATE(doccde_k,(nband_k))
300      ABI_ALLOCATE(dhdk2_r,(nband_k,nband_k,3,3))
301      ABI_ALLOCATE(dhdk2_g,(natom,nband_k,nband_k))
302 
303      cond_nd   = zero
304      kin11_k   = zero
305      kin12_k   = zero
306      kin21_k   = zero
307      kin22_k   = zero
308      np_sum_k1 = zero
309      np_sum_k2 = zero
310      socc_k    = zero
311      dhdk2_r   = zero
312      dhdk2_g   = zero
313 
314 !    eigenvalue for k-point
315      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
316 !    first derivative eigenvalues for k-point
317      psinablapsi=zero
318      read(opt_unt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k)
319      read(opt_unt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k)
320      read(opt_unt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k)
321 !    DEBUG
322 !    write(963,*)isppol,ikpt,((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k)
323 !    write(963,*)isppol,ikpt,((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k)
324 !    write(963,*)isppol,ikpt,((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k)
325 !    ENDDEBUG
326 
327 !    occupation numbers for k-point
328      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
329 !    derivative of occupation number for k-point
330      doccde_k(:)=doccde(1+bdtot_index:nband_k+bdtot_index)
331 
332 !    LOOP OVER BANDS
333      do iband=1,nband_k
334        do jband=1,nband_k
335          do l1=1,3
336            do l2=1,3
337              dhdk2_r(iband,jband,l1,l2)=dhdk2_r(iband,jband,l1,l2)+(&
338 &             psinablapsi(1,l1,iband,jband)*psinablapsi(1,l2,iband,jband)&
339 &             +psinablapsi(2,l1,iband,jband)*psinablapsi(2,l2,iband,jband))
340            end do
341          end do
342 
343          do l1=1,3
344            dhdk2_g(1,iband,jband)=dhdk2_g(1,iband,jband)+( &
345 &           psinablapsi(1,l1,iband,jband)*psinablapsi(1,l1,iband,jband) &
346 &           +psinablapsi(2,l1,iband,jband)*psinablapsi(2,l1,iband,jband))
347          end do
348 
349          diff_occ = occ_k(iband)-occ_k(jband)
350          if (dabs(diff_occ)>=tol8) then
351 
352 !          Conductivity for each omega
353 !          omin = zero
354            do iom=1,mom
355              oml=oml1(iom)
356              if (jband>iband) then
357                sig= dhdk2_g(1,iband,jband)&
358 &               *(diff_occ)/oml*(dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)&
359 &               -dexp(-((eig0_k(iband)-eig0_k(jband)-oml)/dom)**2))
360                kin11_k(iom)=kin11_k(iom)+sig
361                kin12_k(iom)=kin12_k(iom)-sig*(eig0_k(jband)-fermie)
362                kin21_k(iom)=kin21_k(iom)-sig*(eig0_k(iband)-fermie)
363                kin22_k(iom)=kin22_k(iom) + &
364 &               sig*(eig0_k(iband)-fermie)*(eig0_k(jband)-fermie)
365              end if
366              do l1=1,3
367                do l2=1,3
368                  cond_nd(iom,l1,l2)=cond_nd(iom,l1,l2) +dhdk2_r(iband,jband,l1,l2)&
369 &                 *(diff_occ)/oml*dexp(-((eig0_k(jband)-eig0_k(iband)-oml)/dom)**2)
370                end do
371              end do
372            end do
373 
374 !          Sumrule start
375            if (dabs(eig0_k(iband)-eig0_k(jband))>=tol10) then
376              np_sum_k1=np_sum_k1 -dhdk2_g(1,iband,jband)&
377 &             *(diff_occ)/(eig0_k(iband)-eig0_k(jband))
378            else
379              np_sum_k2=np_sum_k2 - doccde_k(iband)*dhdk2_g(1,iband,jband)
380            end if
381 
382 !          end loop over band
383          end if
384        end do
385        socc_k=socc_k+occ_k(iband)
386      end do
387      do iom=1,mom
388        kin11(iom,isppol)=kin11(iom,isppol)+wtk(ikpt)*kin11_k(iom)
389        kin12(iom)=kin12(iom)+wtk(ikpt)*kin12_k(iom)
390        kin21(iom)=kin21(iom)+wtk(ikpt)*kin21_k(iom)
391        kin22(iom)=kin22(iom)+wtk(ikpt)*kin22_k(iom)
392      end do
393      np_sum=np_sum + wtk(ikpt)*(np_sum_k1+np_sum_k2)
394      socc=socc+wtk(ikpt)*socc_k
395 
396 !    Validity limit
397      deltae=deltae+(eig0_k(nband_k)-fermie)
398 
399      bdtot_index=bdtot_index+nband_k
400      ABI_DEALLOCATE(eig0_k)
401      ABI_DEALLOCATE(occ_k)
402      ABI_DEALLOCATE(doccde_k)
403      ABI_DEALLOCATE(dhdk2_r)
404      ABI_DEALLOCATE(dhdk2_g)
405 !    End loop over k
406    end do
407 !  End loop over Spin
408  end do
409 
410  write(std_out,'(a,3f10.5)')' sumrule           =',np_sum/socc/three/dble(nsppol),socc
411  write(std_out,'(a,f10.5,a,f10.5,a)')&
412 & ' Emax-Efermi       =',deltae/dble(nkpt*nsppol),' Ha',deltae/dble(nkpt*nsppol)*Ha_eV,' eV'
413 
414  if (open_file(trim(filnam_out)//'_Lij',msg, newunit=lij_unt, form='formatted', action="write") /= 0) then
415    MSG_ERROR(msg)
416  end if
417  write(lij_unt,'(a)')' # omega(ua) L12 L21 L22 L22'
418 
419  if (open_file(trim(filnam_out)//'_sig', msg, newunit=sig_unt, form='formatted', action="write") /= 0) then
420    MSG_ERROR(msg)
421  end if
422  if (nsppol==1) then
423    write(sig_unt,'(a)')' # omega(ua) hbar*omega(eV)    cond(ua)             cond(ohm.cm)-1'
424  else
425    write(sig_unt,'(2a)')' # omega(ua) hbar*omega(eV)      cond(ua)            cond(ohm.cm)-1',&
426 &   '      cond(ohm.cm)-1 UP      cond(ohm.cm)-1 DN'
427  end if
428 
429  if (open_file(trim(filnam_out)//'_Kth', msg, newunit=kth_unt, form='formatted', action="write") /=0) then
430    MSG_ERROR(msg)
431  end if
432  write(kth_unt,'(a)')&
433 & ' #omega(ua) hbar*omega(eV)  thermal cond(ua)   Kth(W/m/K)   thermopower(ua)   Stp(microohm/K)'
434 
435  if (open_file(trim(filnam_out)//'.out', msg, newunit=ocond_unt, form='formatted', action="write") /= 0) then
436    MSG_ERROR(msg)
437  end if
438  write(ocond_unt,'(a)' )' #Conducti output file:'
439  write(ocond_unt,'(a)' )' #Contains all results produced by conducti utility'
440  write(ocond_unt,'(a)' )' '
441  write(ocond_unt,'(a)')' # omega(ua)       cond(ua)             thermal cond(ua)       thermopower(ua)'
442 
443 !call isfile(filnam_out,'new')
444 
445 !Compute thermal conductivity and thermopower
446  do iom=1,mom
447    oml=oml1(iom)
448    do isppol=1,nsppol
449      kin11(iom,isppol)=kin11(iom,isppol)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
450      if (dabs(kin11(iom,isppol))<10.0d-20) kin11(iom,isppol)=zero
451      sig_abs(iom)=sig_abs(iom)+kin11(iom,isppol)
452    end do
453    kin21(iom)=kin21(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
454    kin12(iom)=kin12(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
455    kin22(iom)=kin22(iom)*two_pi*third/(dom*ucvol)*half/dsqrt(pi)
456    Kth(iom)=kin22(iom)
457    Stp(iom)=zero
458    if(sig_abs(iom)/=zero)  then
459      Kth(iom)=Kth(iom)-(kin12(iom)*kin21(iom)/sig_abs(iom))
460      Stp(iom)=kin12(iom)/(sig_abs(iom)*Tatm)
461    end if
462    if (dabs(Kth(iom))<10.0d-20) Kth(iom)=zero
463    if (dabs(Stp(iom))<10.0d-20) Stp(iom)=zero
464    if (abs(kin12(iom))<10.0d-80) kin12=zero
465    if (abs(kin21(iom))<10.0d-80) kin21=zero
466    if (abs(kin22(iom))<10.0d-80) kin22=zero
467 
468    write(lij_unt,'(f12.5,4es22.12)')oml,kin12(iom),kin21(iom),kin22(iom),kin22(iom)/Tatm*3.4057d9
469 
470    if (nsppol==1) then
471      write(sig_unt,'(2f12.5,2es22.12)') oml,oml*Ha_eV,sig_abs(iom),sig_abs(iom)*Ohmcm
472    else
473      write(sig_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,sig_abs(iom),sig_abs(iom)*Ohmcm,&
474 &     kin11(iom,1)*Ohmcm,kin11(iom,2)*Ohmcm
475    end if
476    write(kth_unt,'(2f12.5,4es22.12)') oml,oml*Ha_eV,Kth(iom),Kth(iom)*3.4057d9/Tatm,&
477 &   Stp(iom),Stp(iom)*3.6753d-2
478    write(ocond_unt,'(1f12.5,3es22.12)') oml,sig_abs(iom),Kth(iom),Stp(iom)
479  end do
480 
481 !Calculate the imaginary part of the conductivity (principal value)
482 !+derived optical properties.
483  call msig (sig_abs,mom,oml1,filnam_out)
484 
485  close(lij_unt)
486  close(sig_unt)
487  close(kth_unt)
488  close(ocond_unt)
489 
490  write(std_out,'(2a)')ch10,'OUTPUT'
491  write(std_out,'(a)')trim(filnam_out)//'_Lij : Onsager kinetic coefficients'
492  write(std_out,'(a)')trim(filnam_out)//'_sig : Optical conductivity'
493  write(std_out,'(a)')trim(filnam_out)//'_Kth : Thermal conductivity and thermopower'
494  write(std_out,'(a)')trim(filnam_out)//'_eps : Dielectric fonction'
495  write(std_out,'(a)')trim(filnam_out)//'_abs : n, k, reflectivity, absorption'
496 
497  call WffClose(wff1,ierr)
498 
499  ABI_DEALLOCATE(psinablapsi)
500  ABI_DEALLOCATE(kin11)
501  ABI_DEALLOCATE(kin22)
502  ABI_DEALLOCATE(kin12)
503  ABI_DEALLOCATE(kin21)
504  ABI_DEALLOCATE(kin11_k)
505  ABI_DEALLOCATE(kin22_k)
506  ABI_DEALLOCATE(kin12_k)
507  ABI_DEALLOCATE(kin21_k)
508  ABI_DEALLOCATE(Stp)
509  ABI_DEALLOCATE(Kth)
510  ABI_DEALLOCATE(cond_nd)
511  ABI_DEALLOCATE(sig_abs)
512  ABI_DEALLOCATE(eigen0)
513  ABI_DEALLOCATE(nband)
514  ABI_DEALLOCATE(oml1)
515  ABI_DEALLOCATE(occ)
516  ABI_DEALLOCATE(doccde)
517  ABI_DEALLOCATE(wtk)
518 
519  call hdr_free(hdr)
520 
521 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

  (main routine)

OUTPUT

  (main routine)

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 windows for computations of sigma
  wtk(nkpt)=weight assigned to each k point.

PARENTS

      conducti

CHILDREN

      hdr_free,hdr_io,hdr_read_from_fname,metric,wffclose,wffopen

SOURCE

572  subroutine conducti_paw_core(filnam,filnam_out,mpi_enreg)
573 
574 
575 !This section has been created automatically by the script Abilint (TD).
576 !Do not modify the following lines by hand.
577 #undef ABI_FUNC
578 #define ABI_FUNC 'conducti_paw_core'
579 !End of the abilint section
580 
581  implicit none
582 
583 !Arguments -----------------------------------
584 !scalars
585  character(len=fnlen) :: filnam,filnam_out
586  type(MPI_type),intent(in) :: mpi_enreg
587 
588 !Local variables-------------------------------
589 !scalars
590  integer,parameter :: master=0
591  integer :: iomode,atnbr,bantot,bdtot_index
592  integer :: fform2,headform,iatom,iband,icor,ierr,ikpt
593  integer :: iom,isppol,l1,mband,me,mom
594  integer :: natom,nband_k,nkpt,nphicor,nspinor,nsppol,ntypat
595  integer :: occopt,rdwr,spaceComm,iunt,opt2_unt,sigx_unt
596  real(dp) :: del,diff_occ,ecut,fermie
597  real(dp) :: omin,omax,dom,oml
598  real(dp) :: Tatm,tsmear,ucvol
599  character(len=fnlen) :: filnam2,filnam_gen
600  character(len=500) :: msg
601  type(hdr_type) :: hdr
602  type(wffile_type) :: wff2
603 !arrays
604  integer,allocatable :: nband(:),ncor(:),lcor(:)
605  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
606  real(dp),allocatable :: dhdk2_g(:,:,:)
607  real(dp),allocatable :: eig0_k(:),eigen0(:)
608  real(dp),allocatable :: energy_cor(:)
609  real(dp),allocatable :: occ(:),occ_k(:),oml1(:)
610  real(dp),allocatable :: psinablapsi2(:,:,:,:,:)
611  real(dp),allocatable :: sigx(:,:,:),sigx_av(:,:),wtk(:)
612 
613 ! *********************************************************************************
614  ABI_UNUSED(mpi_enreg%paral_kgb)
615 
616 !Read output file name
617 !write(std_out,'(a)')' conducti_paw_core : Please, give the name of the output file ...'
618 !read(5, '(a)') filnam_out
619 !write(std_out,'(a)')' The name of the output file is :',trim(filnam_out)
620 !Read data file
621  if (open_file(filnam,msg, newunit=iunt, form='formatted', action="read", status="old") /= 0) then
622    MSG_ERROR(msg)
623  end if
624  rewind(iunt)
625  read(iunt,*)
626  read(iunt,'(a)')filnam_gen       ! generic name for the files
627  filnam2=trim(filnam_gen)//'_OPT2'
628 !Read size of the frequency range
629  read(iunt,*) dom,omin,omax,mom,atnbr
630  close(iunt)
631 
632  write(std_out,'(a,i8,3f10.5,a)')' npts,omin,omax,width      =',mom,omin,omax,dom,' Ha'
633  write(std_out,'(a)')'--------------------------------------------'
634  write(std_out,'(a,i4)') 'selected atom for spectro X',atnbr
635  write(std_out,'(a)')'--------------------------------------------'
636 
637 !These default values are typical of sequential use
638  iomode=IO_MODE_FORTRAN; spaceComm=xmpi_comm_self; me=0
639 
640 ! Read the header of the OPT2 file.
641  call hdr_read_from_fname(hdr, filnam2, fform2, spaceComm)
642  call hdr_free(hdr)
643 
644  if (fform2 /= 611) then
645    MSG_ERROR("Abinit8 requires an OPT2 file with fform = 611")
646  end if
647 
648 !Open the optic files
649  opt2_unt = get_unit()
650  call WffOpen(iomode,spaceComm,filnam2,ierr,wff2,master,me,opt2_unt)
651 
652 !Read the header
653  rdwr=1
654  call hdr_io(fform2,hdr,rdwr,wff2)
655 
656 !Extract info from the header
657  headform=hdr%headform
658  bantot=hdr%bantot
659  ecut=hdr%ecut_eff
660  natom=hdr%natom
661  nkpt=hdr%nkpt
662  nspinor=hdr%nspinor
663  nsppol=hdr%nsppol
664  ntypat=hdr%ntypat
665  occopt=hdr%occopt
666  rprimd(:,:)=hdr%rprimd(:,:)
667  ABI_ALLOCATE(nband,(nkpt*nsppol))
668  ABI_ALLOCATE(occ,(bantot))
669  ABI_ALLOCATE(wtk,(nkpt))
670  fermie=hdr%fermie
671  tsmear=hdr%tsmear
672  occ(1:bantot)=hdr%occ(1:bantot)
673  wtk(1:nkpt)=hdr%wtk(1:nkpt)
674  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
675 
676 !Get mband, as the maximum value of nband(nkpt)
677  mband=maxval(nband(:))
678 
679  write(std_out,*)
680  write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1,1:3)
681  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(2,1:3)
682  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(3,1:3)
683  write(std_out,'(a,i8)')       ' natom             =',natom
684  write(std_out,'(a,3i8)')      ' nkpt,mband,nsppol        =',nkpt,mband,nsppol
685  write(std_out, '(a, f10.5,a)' ) ' ecut              =',ecut,' Ha'
686  write(std_out,'(a,f10.5,a,f10.5,a)' )' fermie            =',fermie,' Ha',fermie*Ha_eV,' eV'
687  Tatm=tsmear*Ha_K
688  write(std_out,'(a,f12.5,a,f12.5,a)') ' Temp              =',tsmear,' Ha ',Tatm,' Kelvin'
689 
690  ABI_ALLOCATE(eigen0,(mband*nkpt*nsppol))
691  read(opt2_unt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
692 !
693  write(std_out,'(a)')'--------------------------------------------'
694  read(opt2_unt) nphicor
695  write(std_out,'(a,i4)') 'Number of core orbitals nc=',nphicor
696  ABI_ALLOCATE(ncor,(nphicor))
697  ABI_ALLOCATE(lcor,(nphicor))
698  ABI_ALLOCATE(energy_cor,(nphicor))
699  do icor=1,nphicor
700    read(opt2_unt) ncor(icor),lcor(icor),energy_cor(icor)
701    write(std_out,'(a,2i4,f10.5)') ' n, l, Energy (Ha): ',ncor(icor),lcor(icor),energy_cor(icor)
702  end do
703  write(std_out,'(a)')'--------------------------------------------'
704 
705 !---------------------------------------------------------------------------------
706 !gmet inversion to get ucvol
707  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
708 
709 !---------------------------------------------------------------------------------
710 !size of the frequency range
711  del=(omax-omin)/(mom-1)
712  ABI_ALLOCATE(oml1,(mom))
713  do iom=1,mom
714    oml1(iom)=omin+dble(iom-1)*del
715  end do
716 
717  ABI_ALLOCATE(sigx,(natom,mom,nphicor))
718  ABI_ALLOCATE(sigx_av,(mom,nphicor))
719 !---------------------------------------------------------------------------------
720 !SpectroX  -------
721 !
722  ABI_ALLOCATE(psinablapsi2,(2,3,mband,nphicor,natom))
723  sigx=zero
724  sigx_av=zero
725  bdtot_index = 0
726 
727 !LOOP OVER SPINS
728  do isppol=1,nsppol
729 
730 !  LOOP OVER K-POINTS
731    do ikpt=1,nkpt
732      nband_k=nband(ikpt+(isppol-1)*nkpt)
733      ABI_ALLOCATE(eig0_k,(nband_k))
734      ABI_ALLOCATE(occ_k,(nband_k))
735      ABI_ALLOCATE(dhdk2_g,(natom,nband_k,nphicor))
736 
737      dhdk2_g   = zero
738      psinablapsi2=zero
739 
740 !    eigenvalue for k-point
741      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
742 
743 !    occupation numbers for k-point
744      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
745 
746      do iatom=1,natom
747 !      first derivative eigenvalues for k-point
748        read(opt2_unt) ((psinablapsi2(1:2,1,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor)
749        read(opt2_unt) ((psinablapsi2(1:2,2,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor)
750        read(opt2_unt) ((psinablapsi2(1:2,3,iband,icor,iatom),iband=1,nband_k),icor=1,nphicor)
751      end do
752 
753 !    LOOP OVER ATOMS/BANDS
754      do iatom=1,natom
755        do iband=1,nband_k
756          do icor=1,nphicor
757 
758            do l1=1,3
759              dhdk2_g(iatom,iband,icor)=dhdk2_g(iatom,iband,icor)+( &
760 &             psinablapsi2(1,l1,iband,icor,iatom)*psinablapsi2(1,l1,iband,icor,iatom) &
761 &             +psinablapsi2(2,l1,iband,icor,iatom)*psinablapsi2(2,l1,iband,icor,iatom))
762            end do
763 
764            diff_occ = (two/dble(nsppol))-occ_k(iband)
765 !          Spectro for each omega
766            omin = -1.0
767            do iom=1,mom
768              oml=-energy_cor(icor)+oml1(iom)+omin
769              sigx(iatom,iom,icor)=sigx(iatom,iom,icor)+ wtk(ikpt)*dhdk2_g(iatom,iband,icor)&
770 &             *(diff_occ)/oml*dexp(-((-energy_cor(icor)+eig0_k(iband)-oml)/dom)**2)
771            end do
772          end do !icor
773        end do  !iband
774      end do !iatom
775      bdtot_index=bdtot_index+nband_k
776      ABI_DEALLOCATE(eig0_k)
777      ABI_DEALLOCATE(occ_k)
778      ABI_DEALLOCATE(dhdk2_g)
779 !    end loop over k
780    end do
781 !  end loop over spins
782  end do
783  ABI_DEALLOCATE(psinablapsi2)
784 
785  do iatom=1,natom
786    do icor=1,nphicor
787      do iom=1,mom
788        if(sigx(iatom,iom,icor)<=tol16) sigx(iatom,iom,icor)=zero
789      end do
790    end do
791  end do ! iatom
792 
793  sigx=sigx*two_pi*third*dble(natom)/(dom*ucvol)*half/dsqrt(pi)
794 
795  do icor=1,nphicor
796    do iom=1,mom
797      do iatom=1,natom
798        sigx_av(iom,icor) =sigx_av(iom,icor)+sigx(iatom,iom,icor)/dble(natom)
799      end do
800    end do
801  end do
802 
803  if (open_file(trim(filnam_out)//'_sigX', msg, newunit=sigx_unt, form='formatted', action="write") /= 0) then
804    MSG_ERROR(msg)
805  end if
806  do iom=1,mom
807    write(sigx_unt,'(9(1x,e14.8))') &
808 &   ((-energy_cor(icor)+oml1(iom)+omin),sigx_av(iom,icor),sigx(atnbr,iom,icor),icor=1,nphicor)
809  end do
810  close(sigx_unt)
811 
812  call WffClose(wff2,ierr)
813 
814  ABI_DEALLOCATE(sigx)
815  ABI_DEALLOCATE(sigx_av)
816  ABI_DEALLOCATE(ncor)
817  ABI_DEALLOCATE(lcor)
818  ABI_DEALLOCATE(energy_cor)
819  ABI_DEALLOCATE(eigen0)
820  ABI_DEALLOCATE(nband)
821  ABI_DEALLOCATE(oml1)
822  ABI_DEALLOCATE(occ)
823  ABI_DEALLOCATE(wtk)
824 
825  call hdr_free(hdr)
826 
827 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

PARENTS

      conducti_nc,conducti_paw

CHILDREN

      intrpl

SOURCE

1534 subroutine msig(fcti,npti,xi,filnam_out_sig)
1535 
1536 
1537 !This section has been created automatically by the script Abilint (TD).
1538 !Do not modify the following lines by hand.
1539 #undef ABI_FUNC
1540 #define ABI_FUNC 'msig'
1541 !End of the abilint section
1542 
1543  implicit none
1544 
1545 !Arguments -----------------------------------
1546 !scalars
1547  integer,intent(in) :: npti
1548 !arrays
1549  real(dp),intent(in) :: fcti(npti),xi(npti)
1550  character(len=fnlen),intent(in) :: filnam_out_sig
1551 
1552 !Local variables-------------------------------
1553 !scalars
1554  integer,parameter :: npt=10000
1555  integer :: ii,ip,npt1,npt2,eps_unt,abs_unt
1556  real(dp),parameter :: del=0.001_dp,ohmtosec=9.d11
1557  real(dp) :: dx,dx1,dx2,eps1,eps2,idel,komega,pole,refl,sigma2,xsum
1558  character(len=500) :: msg
1559 !arrays
1560  real(dp),allocatable :: abso(:),fct(:),fct1(:),fct2(:),fct3(:),fct4(:),fct5(:)
1561  real(dp),allocatable :: fctii(:),fp(:),fpp(:),fppp(:),nomega(:),ppsig(:)
1562  real(dp),allocatable :: x1(:),x2(:)
1563 
1564 ! *********************************************************************************
1565 !BEGIN EXECUTABLE SECTION
1566 
1567  write(std_out,'(2a)')ch10,'Calculate the principal value and related optical properties'
1568  write(std_out,'(a)')'following W.J. Thomson computer in physics vol 12 p94 1998 for '
1569  write(std_out,'(a)')'the principal value. S. Mazevet'
1570  write(std_out,'(a)')'OPTIONS'
1571  write(std_out,'(a)')'use default number of integration pts: npt=10000'
1572  write(std_out,'(a)')'Use default value for delta interval: del=1e-3'
1573 
1574  if (open_file(trim(filnam_out_sig)//'_eps',msg, newunit=eps_unt,status='replace',action="write") /= 0) then
1575    MSG_ERROR(msg)
1576  end if
1577  write(eps_unt,'(a)')'#energy (eV),sigma_1(Ohm-1cm-1),sigma_2(Ohm-1cm-1),epsilon_1,epsilon_2'
1578 
1579  if (open_file(trim(filnam_out_sig)//'_abs', msg, newunit=abs_unt, status='replace',action="write") /= 0) then
1580    MSG_ERROR(msg)
1581  end if
1582  write(abs_unt,'(a)')'#energy(eV),nomega,komega,refl.,abso.(cm-1)'
1583 
1584  ABI_ALLOCATE(fct,(npt))
1585  ABI_ALLOCATE(fct2,(npt))
1586  ABI_ALLOCATE(fct3,(npt))
1587  ABI_ALLOCATE(fct4,(npt))
1588  ABI_ALLOCATE(fct5,(npt))
1589  ABI_ALLOCATE(fp,(npt))
1590  ABI_ALLOCATE(fpp,(npt))
1591  ABI_ALLOCATE(fppp,(npt))
1592  ABI_ALLOCATE(x1,(npt))
1593  ABI_ALLOCATE(x2,(npt))
1594  ABI_ALLOCATE(fct1,(npt))
1595  ABI_ALLOCATE(ppsig,(npt))
1596  ABI_ALLOCATE(fctii,(npt))
1597  ABI_ALLOCATE(abso,(npt))
1598  ABI_ALLOCATE(nomega,(npt))
1599 
1600  if (npti > npt) then
1601    write (std_out,*) 'msig: input npti is too large for hard coded npt array size = ', npt
1602    MSG_ERROR("Aborting now")
1603  end if
1604 
1605 !loop on the initial energy grid
1606  do ip=1,npti
1607 
1608 !  adjust the interval before and after the pole to reflect range/npt interval
1609    xsum=zero
1610    dx=(xi(npti)-xi(1))/dble(npt-1)
1611    pole=xi(ip)
1612    npt1=int((pole-del)/dx)
1613    dx1=zero
1614    if(npt1/=1) dx1=(pole-del)/(npt1-1)
1615    npt2=int((xi(npti)-pole-del)/dx)
1616    dx2=(xi(npti)-pole-del)/(npt2-1)
1617 
1618 !  for the moment skip the pp calculation when the pole if too close to the end of the range
1619    if(npt1<=1.or.npt2<=1) then
1620 
1621      xsum=zero
1622      ppsig(ip)=zero
1623 
1624    else
1625 
1626 !    define the fct for which the pp calculation is needed using xi^2-pole^2 factorization
1627      fctii(1:npti) = zero
1628      fctii(1:npti)=fcti(1:npti)*pole/(xi(1:npti)+pole)
1629 
1630 !    define the grid on each side of the pole x1 before x2 after
1631      do ii=1,npt1
1632        x1(ii)=dx1*dble(ii-1)
1633      end do
1634      do ii=1,npt2
1635        x2(ii)=pole+del+dx2*dble(ii-1)
1636      end do
1637 
1638 !    interpolate the initial fct fii on the new grids x1 and x2 (cubic spline)
1639 !    write(std_out,*) npti,npt1
1640 
1641 !    MJV 6/12/2008:
1642 !    for each use of fctii should ensure that npt1 npt2 etc... are less than
1643 !    npt=len(fctii)
1644      call intrpl(npti,xi,fctii,npt1,x1,fct4,fct1,fct5,1)
1645      call intrpl(npti,xi,fctii,npt2,x2,fct3,fct2,fct5,1)
1646 
1647 !    calculate the two integrals from 0-->pole-lamda and pole+lamda--> end range
1648 !    trapezoidal integration
1649      do ii=1,npt1
1650        fct1(ii)=fct4(ii)/(x1(ii)-pole)
1651      end do
1652      do ii=1,npt2
1653        fct2(ii)=fct3(ii)/(x2(ii)-pole)
1654      end do
1655 
1656      do ii=2,npt1-1
1657        xsum=xsum+fct1(ii)*dx1
1658      end do
1659      do ii=2,npt2-1
1660        xsum=xsum+fct2(ii)*dx2
1661      end do
1662      xsum=xsum+half*(fct1(1)+fct1(npt1))*dx1+half*(fct2(1)+fct2(npt2))*dx2
1663 
1664 !    calculate the first and third derivative at the pole and add the taylor expansion
1665      call intrpl(npti,xi,fctii,npti,xi,fct3,fct4,fct5,1)
1666      call intrpl(npti,xi,fct4,1,(/pole/),fp,fpp,fppp,1)
1667 
1668      idel=two*fp(1)*(del)+fppp(1)*(del**3)/nine
1669      xsum=xsum+idel
1670 
1671    end if
1672 
1673 !  calculate the derivated optical quantities and output the value
1674    sigma2=(-two/pi)*xsum
1675    eps1=one-(four_pi*sigma2/(pole))
1676    eps2=four*fcti(ip)*pi/(pole)
1677 
1678 !  A special treatment of the case where eps2 is very small compared to eps1 is needed
1679    if(eps2**2 > eps1**2 * tol12)then
1680      nomega(ip)=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2)))
1681      komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2)))
1682      abso(ip)=four_pi*fcti(ip)*ohmtosec*Ohmcm/nomega(ip)/(Sp_Lt_SI*100._dp)
1683    else if(eps1>zero)then
1684      nomega(ip)=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2)))
1685      komega=half*abs(eps2/sqrt(eps1))
1686      abso(ip)=four_pi*fcti(ip)*ohmtosec*Ohmcm/nomega(ip)/(Sp_Lt_SI*100._dp)
1687    else if(eps1<zero)then
1688      nomega(ip)=half*abs(eps2/sqrt(-eps1))
1689      komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2)))
1690      abso(ip)=two*sqrt(-eps1)*pole*ohmtosec*Ohmcm/(Sp_Lt_SI*100._dp)
1691    end if
1692 
1693    refl=((one-nomega(ip))**2 + komega**2)/ &
1694 &   ((one+nomega(ip))**2 + komega**2)
1695 
1696    write(eps_unt,'(5e18.10)') Ha_eV*pole,fcti(ip)*Ohmcm,sigma2*Ohmcm,eps1,eps2
1697    write(abs_unt,'(5e18.10)') Ha_eV*pole,nomega(ip),komega,refl,abso(ip)
1698 
1699  end do
1700 
1701  close(eps_unt)
1702  close(abs_unt)
1703 
1704  ABI_DEALLOCATE(fct)
1705  ABI_DEALLOCATE(x1)
1706  ABI_DEALLOCATE(x2)
1707  ABI_DEALLOCATE(fct2)
1708  ABI_DEALLOCATE(fct3)
1709  ABI_DEALLOCATE(fct4)
1710  ABI_DEALLOCATE(fct5)
1711  ABI_DEALLOCATE(fp)
1712  ABI_DEALLOCATE(fpp)
1713  ABI_DEALLOCATE(fppp)
1714  ABI_DEALLOCATE(fct1)
1715  ABI_DEALLOCATE(ppsig)
1716  ABI_DEALLOCATE(fctii)
1717  ABI_DEALLOCATE(abso)
1718  ABI_DEALLOCATE(nomega)
1719 
1720 end subroutine msig