scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

This is a simple transport code that makes use of the a2Fsave files. 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@11804 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
burak 2015-10-22 16:38:56 +00:00
parent 7f0c8479a5
commit 812466dcf9
11 changed files with 1608 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
PP/simple_transport/cryst_to_car.f90 Normal file
View File

@ -0,0 +1,57 @@
subroutine cryst_to_cart (nvec, vec, trmat, iflag)
!-----------------------------------------------------------------------
!
! This routine transforms the atomic positions or the k-point
! components from crystallographic to cartesian coordinates
! ( iflag=1 ) and viceversa ( iflag=-1 ).
! Output cartesian coordinates are stored in the input ('vec') array
!
!
implicit none
!
integer, intent(in) :: nvec, iflag
! nvec: number of vectors (atomic positions or k-points)
! to be transformed from crystal to cartesian and vice versa
! iflag: gives the direction of the transformation
double precision, intent(in) :: trmat (3, 3)
! trmat: transformation matrix
! if iflag=1:
! trmat = at , basis of the real-space lattice, for atoms or
! = bg , basis of the reciprocal-space lattice, for k-points
! if iflag=-1: the opposite
double precision, intent(inout) :: vec (3, nvec)
! coordinates of the vector (atomic positions or k-points) to be
! transformed - overwritten on output
!
! local variables
!
integer :: nv, kpol
! counter on vectors
! counter on polarizations
double precision :: vau (3)
! workspace
!
! Compute the cartesian coordinates of each vectors
! (atomic positions or k-points components)
!
do nv = 1, nvec
if (iflag.eq.1) then
do kpol = 1, 3
vau (kpol) = trmat (kpol, 1) * vec (1, nv) + trmat (kpol, 2) &
* vec (2, nv) + trmat (kpol, 3) * vec (3, nv)
end do
else
do kpol = 1, 3
vau (kpol) = trmat (1, kpol) * vec (1, nv) + trmat (2, kpol) &
* vec (2, nv) + trmat (3, kpol) * vec (3, nv)
end do
endif
do kpol = 1, 3
vec (kpol, nv) = vau (kpol)
end do
end do
!
return
end subroutine cryst_to_cart

239
PP/simple_transport/dos.f90 Normal file
View File

@ -0,0 +1,239 @@
program dos
!
! Written by Burak Himmetoglu (burakhmmtgl@gmail.com)
! Uses some parts of the PW distribution
!
! Computation of DOS using the variable smearing technique.
! See: J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, PRB 75, 195121 (2007)
!
! Description of the input card:
!
! fil_a2F : prefix.a2Fsave file
! fil_info : file containing direct and reciprocal lattice vectors
! T : Temperature in K
! Nen : Number of energy values.
! emin, emax : Min and Max energies.
! lsoc : if .true. then the band structure is non-collinear
! lfix : if .true. DOS with constant smearing is also computed.
! degauss : if lsoc is .true. this is the constant smearing value in eV.
! cbm, vbm : Conduction band min. and Valence band max. (in eV)
! aa : Parameter determining the k-dependent smearing deg_nk : deg_nk = aa * delta_k * (de_nk/delta_k).
! Larger aa will lead to smoother DOS. Needs to be tested for various grids.
! nthreads : Number of threads for OpenMP parallelization.
!$ USE omp_lib
USE smearing_mod
implicit none
!
!
integer :: i,j,iq,ik,ie,nk1fit,nk2fit,nk3fit,nkfit, &
& nbnd, nksfit, npk, nsym, &
& s(3,3,48),ns,nrot,ibnd,io,Nen, ind_k
!
! OpenMP
integer :: TID, nthreads
double precision :: t0
!
double precision :: at(3,3), bg(3,3),deg,degauss,wk,T, &
& cbm, vbm, ef_mid, doping, nelec, &
& en(1000), de(1000), ne(1000),emin, emax, aa
!
!
double precision, allocatable :: xkfit(:,:), etfit(:,:), wkfit(:),&
& vk(:,:,:), dfk(:,:,:)
!
integer, allocatable :: eqkfit(:), sfit(:)
!
logical :: lsoc, lfix
!
character*20 :: fil_kp, fil_a2F, fil_info
!
double precision, PARAMETER :: Rytocm1 = 109737.57990d0, &
& RytoGHz = 3.289828D6, &
& RytoTHz = RytoGHz/1000.d0, &
& RytoeV = 13.60569253d0, &
& tpi = 6.283185307d0, &
& convsig = 2.89d5, &
& KtoRy = 1.d0/38.681648/300/RytoeV, &
& autocm = 5.2917721092d-9, &
degmin = 1e-3
! convsig: conversion factor for
! conductivity. From au to (Ohm cm)-1
namelist /input/ fil_info, fil_a2F, Nen, aa, T, degauss, &
& cbm, vbm, emin, emax, nthreads, lsoc, lfix
read(5,input)
!Set number of threads
!$ call omp_set_num_threads(nthreads)
t0 = 0.0
npk = 40000
!
! Energy grid
do i=1,Nen
en(i) = (emax-emin)/(Nen-1) * (i-1) + emin
end do
!
! Convert to Ryd
en = en / RytoeV
T = T * KtoRy
!
! Read a2Fsave
open(11,file=fil_a2F,status='unknown')
!
read(11,*) nbnd, nksfit
!
allocate(etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit))
!
! Read band structure and k-point data
read(11,*) etfit
read(11,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit)
read(11,*) wkfit
read(11,*) nk1fit, nk2fit, nk3fit
read(11,* ) nsym
do ns=1,nsym
read(11,*) ((s(i,j,ns),j=1,3),i=1,3)
enddo
!
close(11)
!
! Regular grid size
nkfit=nk1fit*nk2fit*nk3fit
!
! Read info file on k-points (and lattice)
open(11,file=fil_info,status='unknown')
!
read(11,*)
read(11,*) ((at(i,j),i=1,3),j=1,3)
!
read(11,*)
read(11,*)
!
read(11,*) ((bg(i,j),i=1,3),j=1,3)
!
close(11)
!
allocate(eqkfit(nkfit),sfit(nkfit))
allocate(dfk(nbnd,nkfit,3),vk(nbnd,nkfit,3))
! eqkfit : pointers to band energies in uniform grid
! sfit : pointers to symmetries
!
call lint ( nsym, s, .true., at, bg, npk, 0,0,0, &
& nk1fit,nk2fit,nk3fit, nksfit, xkfit, 1, nkfit, eqkfit, sfit)
!
! eqkfit(nk) : maps IBZ to full grid. The full grid is in crystal coords
!
! Deallocate unnecessary variables
deallocate(wkfit,sfit,xkfit)
!
! Weights for regular grid (include lsoc)
if ( lsoc .eqv. .true.) then
wk = 1.0/nkfit
else
wk = 2.0/nkfit
end if
!
!Start timing here
!$ t0 = omp_get_wtime()
!
! Call band velocities and forward derivatives
call vband_ibz( nk1fit,nk2fit,nk3fit,nbnd,nksfit,etfit,eqkfit,at, vk, dfk)
!
!Check number of electrons
nelec = 0.d0
ef_mid = (cbm+vbm)/2.d0/RytoeV ! Fermi level set to midgap
!
!$omp parallel do default(shared) &
!$omp collapse(2) &
!$omp private(ik,ibnd,ind_k,deg) &
!$omp reduction(+: nelec)
do ibnd=1,nbnd
!
do ik=1,nkfit
!
ind_k = eqkfit(ik)
deg = sig0(nk1fit,nk2fit,nk3fit,dfk(ibnd,ik,:),aa)
nelec = nelec + wk * w1gauss((ef_mid-etfit(ibnd,ind_k))/deg,0)
end do
end do
!$omp end parallel do
!
!Compute density of states and number of electrons at a given energy (en(i))
de = 0.d0
ne = 0.d0
do ie=1,Nen
!
!$omp parallel do default(shared) &
!$omp collapse(2) &
!$omp private(ik,ibnd,deg,ind_k) &
!$omp reduction(+: de, ne)
do ik=1,nkfit
!
do ibnd=1,nbnd
!
ind_k = eqkfit(ik)
deg = sig0(nk1fit,nk2fit,nk3fit,dfk(ibnd,ik,:),aa)
!
de(ie) = de(ie) + wk * w0gauss((etfit(ibnd,ind_k)-en(ie))/deg)/deg
ne(ie) = ne(ie) + wk * w1gauss((etfit(ibnd,ind_k)-en(ie))/deg,0)
!
end do ! ibnd
end do ! ik
!$omp end parallel do
end do ! ie
!
! End timing
!$ t0 = omp_get_wtime()
!
open(11,file='dos.out',status='unknown')
if (t0 >0) then
write(11,"(A,I5)") "#Number of threads = ", nthreads
write(11,"(A,e14.6)") "#Integration time(s) = ", t0
end if
!
do ie=1,Nen
write(11,"(3f14.6)") en(ie)*RytoeV, de(ie)/RytoeV, ne(ie)
end do
close(11)
!
if (lfix .eqv. .true.) then
!Compute density of states and number of electrons at a given energy (en(i)) for
!fixed smearing of degauss
de = 0.d0
ne = 0.d0
deg = degauss / RytoeV
do ie=1,Nen
!
!$omp parallel do default(shared) &
!$omp collapse(2) &
!$omp private(ik,ibnd,ind_k) &
!$omp reduction(+: de,ne)
do ik=1,nkfit
!
do ibnd=1,nbnd
!
ind_k = eqkfit(ik)
de(ie) = de(ie) + wk * w0gauss((etfit(ibnd,ind_k)-en(ie))/deg)/deg
ne(ie) = ne(ie) + wk * w1gauss((etfit(ibnd,ind_k)-en(ie))/deg,0)
!
end do ! ibnd
end do ! ik
!$omp end parallel do
end do ! ie
!
open(11,file='dos.fix.out',status='unknown')
do ie=1,Nen
write(11,"(3f14.6)") en(ie)*RytoeV, de(ie)/RytoeV, ne(ie)
end do
close(11)
end if ! lfix
!
! Free memory
deallocate(eqkfit,etfit,vk,dfk)
!
end program dos

281
PP/simple_transport/ef.f90 Normal file
View File

@ -0,0 +1,281 @@
program efermi
!---------------------------------------------------------------------------------
! Written by Burak Himmetoglu (burakhmmtgl@gmail.com)
! Uses some parts of the PW distribution
!
! This is a simple code that reads the electron concentration and band
! structure (from a2Fsave file) and then computes the Fermi level.
! The system must be insulating.
!
! Description of the input card:
!
! fil_a2F : prefix.a2Fsave file
! fil_info : file containing direct and reciprocal lattice vectors
! T : Temperature in K
! vol : Volume in au^3
! cbm, vbm : Conduction band min. and Valence band max. (in eV)
! doping : Electron concentration (in cm^-3). For hole doping, use negative value
! ndop : Number of doping levels considered
! nthreads : Number of threads for OpenMP parallelization
!$ use omp_lib
implicit none
!
integer :: i,j,iq,ik,imu,itemp,nu,nqtot,nsig,nat,nk1fit,nk2fit, &
& nk3fit, nkfit, nksfit_real, nbnd, nksfit, npk, nsym, &
& s(3,3,48),ns,nrot,ibnd,ndop
!
! OMP
integer :: TID, nthreads
double precision :: t0
!
double precision :: at(3,3), bg(3,3),vol,T, fd, nelec1, &
& cbm, vbm, ef_mid, doping(10), nelec, &
& Ef_IBZ(10)
!
double precision, allocatable :: xkfit(:,:), etfit(:,:), wkfit(:),&
& wk(:)
!
integer, allocatable :: eqkfit(:), sfit(:)
!
!
character*20 :: fil_kp, fil_a2F, fil_info
!
double precision, PARAMETER :: Rytocm1 = 109737.57990d0, &
& RytoGHz = 3.289828D6, &
& RytoTHz = RytoGHz/1000.d0, &
& RytoeV = 13.60569253d0, &
& tpi = 6.283185307d0, &
& convsig = 2.89d5, &
& KtoRy = 1.d0/38.681648/300/RytoeV, &
& autocm = 5.2917721092d-9, &
& convRH = 9.2522431d-13
! convsig: conversion factor for
! conductivity. From au to (Ohm cm)-1
namelist /input/ fil_info, fil_a2F, vol, cbm, vbm, T, &
& doping, ndop, nthreads
read(5,input)
t0 = 0.0
!$ call omp_set_num_threads(nthreads)
!$ t0 = omp_get_wtime()
npk = 40000
!
! Convert temperature from K to Ryd
T = T * KtoRy
!
! Read a2F
!
open(11,file=fil_a2F,status='unknown')
!
read(11,*) nbnd, nksfit
!
allocate(etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit))
!
! Read band structure and k-point data
read(11,*) etfit
read(11,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit)
read(11,*) wkfit
read(11,*) nk1fit, nk2fit, nk3fit
read(11,* ) nsym
do ns=1,nsym
read(11,*) ((s(i,j,ns),j=1,3),i=1,3)
enddo
!
close(11)
!
! Regular grid size
nkfit=nk1fit*nk2fit*nk3fit
! Read info file on k-points (and lattice)
!
open(11,file=fil_info,status='unknown')
!
read(11,*)
read(11,*) ((at(i,j),i=1,3),j=1,3)
!
read(11,*)
read(11,*)
!
read(11,*) ((bg(i,j),i=1,3),j=1,3)
!
close(11)
!
allocate (eqkfit(nkfit), sfit(nkfit),wk(nkfit))
! eqkfit : pointers to band energies in uniform grid
! sfit : pointers to symmetries
!
call lint ( nsym, s, .true., at, bg, npk, 0,0,0, &
& nk1fit,nk2fit,nk3fit, nksfit, xkfit, 1, nkfit, eqkfit, sfit)
!
! eqkfit(nk) : maps IBZ to full grid. The full grid is in crystal coords
!
! Determine the number of electrons (neutral, insulating system)
ef_mid = (cbm + vbm)/2.0/RytoeV ! Fermi level at the middle of the gap
!
nelec1 = 0.d0
!
! Irreducable BZ
nelec1 = sumk (etfit,nbnd,nksfit,wkfit,T,1,ef_mid)
!
do i=1,ndop
! Determine number of electrons from doping levels (given in cm-3)
nelec = doping(i) * vol * autocm**3 + nelec1
!
! Compute efermi using the bisection method
Ef_IBZ(i) = fermi_en (etfit,wkfit,nbnd,nksfit,nelec,T,1)
!
end do
!
! Memory clean
!
deallocate (xkfit,etfit,wkfit,wk,eqkfit,sfit)
!
!$t0 = omp_get_wtime() - t0
!
if (t0 >0) write(6,"(A,e14.6)") 'Walltime= ', t0
!
do i=1,ndop
write(6,"(A,e14.6,f14.6)") 'doping, Efermi', doping(i), Ef_IBZ(i) * RytoeV
end do
!
!
contains
! -----------------------------------------------
! FUNCTIONS
! -----------------------------------------------
double precision function sumk (et,nbnd,nks,wk,degauss,ngauss,ee)
!$ use omp_lib
implicit none
integer, intent(in) :: nbnd, nks, ngauss
double precision, intent(in) :: ee,degauss
! wk: weight of kpt
! ee: E
! degauss: temperature/degauss
double precision, intent(in) :: et(nbnd,nks), wk(nks)
double precision :: fd, x, arg, a, hp, hd, sqrtpm1
integer :: i, j, ik, n, ni
sqrtpm1 = 1.0d0/1.77245385090551602729d0
sumk = 0.d0
n = 1 ! First order Gauss-Hermite
if (ngauss .eq. 1) then ! FD
!$omp parallel do &
!$omp default(shared) &
!$omp private(i,ik,fd) &
!$omp reduction(+ : sumk)
do i=1,nbnd
do ik=1,nks
fd = 1.d0 / ( 1.d0+exp( (et(i,ik)-ee)/degauss ) )
sumk = sumk + wk(ik) * fd
end do
end do
!$omp end parallel do
else if (ngauss .eq. 0) then ! Gaussian (for zero temperature)
!$omp parallel do &
!$omp default(shared) &
!$omp private(i,ik,fd,x) &
!$omp reduction(+ : sumk)
do i=1,nbnd
do ik=1,nks
x = (ee-et(i,ik))/degauss/dsqrt(2.d0)
fd = 0.5d0 * (1.d0 + derf( x ) )
sumk = sumk + wk(ik) * fd
end do
end do
!$omp end parallel do
else if ( ngauss .eq. 2 ) then ! MP
!$omp parallel do &
!$omp default(shared) &
!$omp private(i,j,ik,fd,x) &
!$omp reduction(+ : sumk)
do i=1,nbnd
do ik=1,nks
!
x = (ee-et(i,ik))/degauss/dsqrt(2.d0)
fd = 0.5d0 * (1.d0 + derf(x) )
hd = 0.d0
arg = min (200.d0, x**2)
hp = exp ( - arg)
ni = 0
a = sqrtpm1
!
do j = 1, n
hd = 2.0d0 * x * hp - 2.0d0 * DBLE (ni) * hd
ni = ni + 1
a = - a / (DBLE (j) * 4.0d0)
fd = fd - a * hd
hp = 2.0d0 * x * hd-2.0d0 * DBLE (ni) * hp
ni = ni + 1
end do
!
sumk = sumk + wk(ik) * fd
end do
end do
!$omp end parallel do
end if
end function sumk
!
double precision function fermi_en (eig,wk,nbnd,nktot,nelec,temperature,ngauss)
implicit none
integer, intent(in) :: nbnd, nktot, ngauss
double precision, intent(in) :: eig(nbnd,nktot),wk(nktot), nelec, &
& temperature
! Local variables
integer :: i, ik
double precision :: Elw, Eup, Ef, sumkup, sumklw, sumkmid
integer, PARAMETER :: maxiter = 300
double precision, PARAMETER :: eps = 1d-10
Elw = eig(1,1)
Eup = eig(nbnd,1)
do ik = 2, nktot
Elw = min ( Elw, eig (1, ik) )
Eup = max ( Eup, eig (nbnd, ik) )
end do
Eup = Eup + 2*temperature
Elw = Elw - 2*temperature
sumkup = sumk(eig, nbnd, nktot, wk, temperature, ngauss, Eup)
sumklw = sumk(eig, nbnd, nktot, wk, temperature, ngauss, Elw)
!
if ( (sumkup - nelec) < -eps .or. (sumklw - nelec) > eps ) then
fermi_en = -99
return
end if
!
do 100 i=1, maxiter
Ef = 0.5d0 * (Eup + Elw)
sumkmid = sumk(eig, nbnd, nktot, wk, temperature, ngauss, Ef)
if ( abs (sumkmid-nelec) < eps ) then
fermi_en = Ef
exit
else if ( (sumkmid-nelec) < -eps) then
Elw = Ef
else
Eup = Ef
end if
100 continue
end function fermi_en
!
end program efermi

294
PP/simple_transport/fermi_int_0.f90 Normal file
View File

@ -0,0 +1,294 @@
program fermi_int_0
!
! Written by Burak Himmetoglu (burakhmmtgl@gmail.com)
! Uses some parts of the PW distribution
!
! This program computes the transport integrals using
! constant scattering rate, and results in conductivity and Seebeck coefficients.
! Integrals that are computed:
!
! I0 = \sum_{n,k} (-df_nk/de)
! I1(i,j) = \sum_{n,k} tau (-df_nk/de) v_nk,i v_nk,j
! I2(i,j) = \sum_{n,k} tau (-df_nk/de) (e_nk - e_F) v_nk,i v_nk,j
!
! Description of the input card:
!
! fil_a2F : prefix.a2Fsave file
! fil_info : file containing direct and reciprocal lattice vectors
! T : Temperature in K
! alat : lattice parameter in a.u. (Bohr)
! vol : volume of lattice in a.u.^3
! invtau : inverse tau in eV
! efmin : minimum eF
! efmax : maximum eF
! Nmu : number of eFs
! phband_i, phband_f : starting and ending indices of bands of interest. The bands must lie between efmin & efmax
! lsoc : if .true. then the band structure is non-collinear
! nthreads : Number of threads for OpenMP parallelization
!$ USE omp_lib
implicit none
integer :: i,j,k,nu,ik,ikk,nk1fit,nk2fit,nk3fit,nkfit, &
& nbnd, nksfit, npk, nsym, Nmu, imu, &
& s(3,3,48),ns,nrot,ibnd,io,phband_i, phband_f, &
& nphband, n, nn, jbnd, ibnd_ph, ind_k
!
double precision :: wk, at(3,3), bg(3,3), efermi, alat, &
& T, wo(3), al(3), invtau,aa,cut,deg, &
& invT, fd, dfd, fac, vol, efmin, efmax, &
& inv_I1(3,3)
!
logical :: lsoc
!
double precision, allocatable :: xkfit(:,:),etfit(:,:),wkfit(:), &
& dfk(:,:,:),vk(:,:,:), En(:)
!
double precision, allocatable :: I0(:) ,I1(:,:,:),I2(:,:,:), &
& sig(:,:,:),Se(:,:,:)
!
integer, allocatable :: eqkfit(:), sfit(:)
!
! OMP variables
double precision :: t0
!
integer :: nthreads
!
character*20 :: fil_a2F, fil_info
!
double precision, PARAMETER :: Rytocm1 = 109737.57990d0, &
& RytoGHz = 3.289828D6, &
& RytoTHz = RytoGHz/1000.d0, &
& RytoeV = 13.60569253d0, &
& tpi = 6.283185307d0, &
& convsig = 2.89d5, &
& KtoRy = 1.d0/38.681648/300/RytoeV
! convsig: conversion factor for
! conductivity. From au to (Ohm cm)-1
!
namelist /input/ fil_info, fil_a2F, T, efmin, efmax, Nmu, &
& phband_i, phband_f, &
& invtau, alat, vol, nthreads, lsoc
read(5,input)
t0 = 0.0
!Set number of threads
!$ call omp_set_num_threads(nthreads)
npk = 40000
!
! Convert from eV to Ryd
T = T * KtoRy
efmin = efmin / RytoeV
efmax = efmax / RytoeV
invtau = invtau / RytoeV
!
! Setup the Efs and transport integral dimensions
allocate(En(Nmu),I0(Nmu))
allocate(I1(Nmu,3,3),I2(Nmu,3,3),sig(Nmu,3,3),Se(Nmu,3,3))
!
do i=1,Nmu
En(i) = (i-1)*(efmax-efmin)/(Nmu-1) + efmin
end do
!
! Total number of bands of interest (usually the number of relevant conduction/valence bands)
nphband = phband_f - phband_i + 1
!
! Read the a2Fsave file
open(11,file=fil_a2F,status='unknown')
!
read(11,*) nbnd, nksfit
!
allocate(etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit))
!
! Read band structure and k-point data
read(11,*) etfit
read(11,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit)
read(11,*) wkfit
read(11,*) nk1fit, nk2fit, nk3fit
read(11,* ) nsym
do ns=1,nsym
read(11,*) ((s(i,j,ns),j=1,3),i=1,3)
enddo
!
close(11)
!
! Regular grid size
nkfit=nk1fit*nk2fit*nk3fit
!
! Read info file on k-points (and lattice)
open(11,file=fil_info,status='unknown')
!
read(11,*)
read(11,*) ((at(i,j),i=1,3),j=1,3)
!
read(11,*)
read(11,*)
!
read(11,*) ((bg(i,j),i=1,3),j=1,3)
!
close(11)
!
allocate (eqkfit(nkfit),sfit(nkfit))
allocate (dfk(nphband,nkfit,3),vk(nphband,nkfit,3))
!
! Find the map between IBZ and the full-grid
call lint ( nsym, s, .true., at, bg, npk, 0,0,0, &
& nk1fit,nk2fit,nk3fit, nksfit, xkfit, 1, nkfit, eqkfit, sfit)
!
! Deallocate unnecessary variables
deallocate(wkfit,sfit,xkfit)
!
! Weights for regular grid
if ( lsoc .eqv. .true.) then
wk = 1.0/nkfit
else
wk = 2.0/nkfit
end if
!
! Call band velocities and forward derivatives
call vband_ibz( nk1fit,nk2fit,nk3fit,nphband,nksfit,etfit(phband_i:phband_f,:),eqkfit,at, vk, dfk)
!
! Include the 2pi/a factor
vk = vk / tpi * alat
!
! Initiate Fermi integrals
I0 = 0.0
I1 = 0.0
I2 = 0.0
invT = 1.0/T
!
!$ t0 = omp_get_wtime()
! Loop over Efs, bands and kpoints
do imu=1,Nmu
!
do ibnd=phband_i,phband_f
!
ibnd_ph = ibnd - phband_i + 1
!
!$omp parallel do default(shared) &
!$omp private(ik,i,j,fd,dfd,fac) &
!$omp reduction(+: I0, I1, I2)
do ik=1,nkfit
!
! Fermi factors
fac = etfit(ibnd,eqkfit(ik)) - En(imu)
fd = 1.0/( exp(fac*invT) + 1.0 )
dfd = invT * fd * (1.0 - fd)
!
! Compute I0 (related to Thomas-Fermi screening)
I0(imu) = I0(imu) + wk * dfd
!
! Compute I1, I2
do i=1,3
do j=1,3
I1(imu,i,j) = I1(imu,i,j) + wk * dfd / invtau * vk(ibnd_ph,ik,i) &
& * vk(ibnd_ph,ik,j)
!
I2(imu,i,j) = I2(imu,i,j) + wk * dfd * fac / invtau * vk(ibnd_ph,ik,i) &
& * vk(ibnd_ph,ik,j)
!
end do ! j
end do ! i
!
end do ! ik
!$omp end parallel do
!
end do ! ibnd
!
end do ! imu
!Total integration time
!$ t0 = omp_get_wtime() - t0
!
! Conductivity (convert to units of Ohm^-1 cm^-1)
sig = I1 * convsig / vol
!
!Compute Seebeck
Se = 0.0
do imu=1,Nmu
call inv33(I1(imu,:,:), inv_I1)
do i=1,3
do j=1,3
do k=1,3
Se(imu,i,j) = Se(imu,i,j) + I2(imu,i,k)*inv_I1(k,j)
end do
end do
end do
end do
! Units (convert to units of V/K)
Se = Se * (-RytoeV * KtoRy / T)
!
! Write D(Ef), conductivity and Seebeck into file
open(11,file='sig_0.out',status='unknown')
open(12,file='Se_0.out',status='unknown')
open(13,file='Def_0.out',status='unknown')
do imu=1,Nmu
write(11,"(10e14.6)") En(imu) * RytoeV, ((sig(imu,i,j),i=1,3),j=1,3)
write(12,"(10e14.6)") En(imu) * RytoeV, ((Se(imu,i,j),i=1,3),j=1,3)
write(13,"(2e14.6)") En(imu) * RytoeV, I0(imu) / RytoeV
end do
close(11)
close(12)
close(13)
!
!
open(11,file='report_0',status='unknown')
if ( t0 > 0 ) then
write(11,"(A,I2)") "Number of threads = ", nthreads
write(11,"(A,e14.6)") "Integration time(s) =", t0
end if
!
write(11,"(A,I5)") "Number of kpoints (regular grid) = ", nkfit
close(11)
!
! Free memory
deallocate(etfit,eqkfit,dfk,vk,En,I0,I1,I2,sig,Se)
contains
!
! Determinant of a 3x3 matrix
double precision function det33 (a)
implicit none
double precision, intent(in) :: a(3,3)
! Determinant of a
det33 = a(1,1)*( a(2,2)*a(3,3)-a(2,3)*a(3,2) ) - &
& a(1,2)*( a(1,1)*a(3,3)-a(1,3)*a(3,1) ) + &
& a(1,3)*( a(2,1)*a(3,2)-a(2,2)*a(3,1) )
end function det33
!
!
! Inverse of a 3x3 matrix (analytical)
subroutine inv33(a, inv_a)
implicit none
double precision, intent(in) :: a(3,3)
double precision, intent(out) :: inv_a(3,3)
double precision :: deta_inv
deta_inv = 1.0 / det33(a)
inv_a(1,1) = deta_inv * ( a(2,2)*a(3,3)-a(2,3)*a(3,2) )
inv_a(1,2) = deta_inv * ( a(1,3)*a(3,2)-a(1,2)*a(3,3) )
inv_a(1,3) = deta_inv * ( a(1,2)*a(2,3)-a(1,3)*a(2,2) )
inv_a(2,1) = deta_inv * ( a(2,3)*a(3,1)-a(2,1)*a(3,3) )
inv_a(2,2) = deta_inv * ( a(1,1)*a(3,3)-a(1,3)*a(3,1) )
inv_a(2,3) = deta_inv * ( a(1,3)*a(2,1)-a(1,1)*a(2,3) )
inv_a(3,1) = deta_inv * ( a(2,1)*a(3,2)-a(2,2)*a(3,1) )
inv_a(3,2) = deta_inv * ( a(1,2)*a(3,1)-a(1,1)*a(3,2) )
inv_a(3,3) = deta_inv * ( a(1,1)*a(2,2)-a(1,2)*a(2,1) )
end subroutine inv33
end program fermi_int_0

288
PP/simple_transport/fermi_int_1.f90 Normal file
View File

@ -0,0 +1,288 @@
program fermi_int_1
!
! Written by Burak Himmetoglu (burakhmmtgl@gmail.com)
! Uses some parts of the PW distribution
!
! This program computes the transport integrals at a given Fermi level using
! constant scattering rate, and results in conductivity and Seebeck coefficients.
!
! Integrals that are computed:
!
! I0 = \sum_{n,k} (-df_nk/de)
! I1(i,j) = \sum_{n,k} tau (-df_nk/de) v_nk,i v_nk,j
! I2(i,j) = \sum_{n,k} tau (-df_nk/de) (e_nk - e_F) v_nk,i v_nk,j
!
! The difference between fermi_int_0 is that the calculation proceeds in a reduced
! grid such that |e_nk-Ef| <= cut*kT. This allows faster integration times for a
! given Fermi level (Ef).
!
! Description of the input card:
!
! fil_a2F : prefix.a2Fsave file
! fil_info : file containing direct and reciprocal lattice vectors
! T : Temperature in K
! alat : lattice parameter in a.u. (Bohr)
! vol : volume of lattice in a.u.^3
! cut : Reduction parameter for integrals (i.e. Integration reduced to the region satisfying |e_nk-Ef| <= cut*kT)
! invtau : inverse tau in eV
! phband_i, phband_f : starting and ending indices of bands of interest. The bands must lie between efmin & efmax
! lsoc : if .true. then the band structure is non-collinear
! nthreads : Number of threads for OpenMP parallelization
!$ USE omp_lib
implicit none
integer :: i,j,k,nu,ik,ikk,nk1fit,nk2fit,nk3fit,nkfit, &
& nbnd, nksfit, npk, nsym, Nmu, imu, &
& s(3,3,48),ns,nrot,ibnd,io,phband_i, phband_f, &
& nphband, n, nn, jbnd, ibnd_ph, ind_k
!
double precision :: wk, at(3,3), bg(3,3), efermi, alat, &
& T, wo(3), al(3), invtau,aa,cut,deg, &
& invT, fd, dfd, fac, vol
!
logical :: lsoc
!
double precision, allocatable :: xkfit(:,:),etfit(:,:),wkfit(:), &
& dfk(:,:,:),vk(:,:,:)
!
double precision :: I0,I1(3,3),I2(3,3),sig(3,3),Se(3,3),inv_I1(3,3)
!
integer, allocatable :: eqkfit(:), sfit(:), iflag(:,:),nkeff(:)
!
! OMP variables
double precision :: t0
!
integer :: nthreads
!
character*20 :: fil_a2F, fil_info
!
double precision, PARAMETER :: Rytocm1 = 109737.57990d0, &
& RytoGHz = 3.289828D6, &
& RytoTHz = RytoGHz/1000.d0, &
& RytoeV = 13.60569253d0, &
& tpi = 6.283185307d0, &
& convsig = 2.89d5, &
& KtoRy = 1.d0/38.681648/300/RytoeV
! convsig: conversion factor for
! conductivity. From au to (Ohm cm)-1
!
namelist /input/ fil_info, fil_a2F, T, phband_i, phband_f,cut, &
& efermi, invtau, alat, vol, nthreads, lsoc
read(5,input)
!Set number of threads
!$ call omp_set_num_threads(nthreads)
t0 = 0.0
npk = 40000
!
! Convert from eV to Ryd
T = T * KtoRy
invtau = invtau / RytoeV
efermi = efermi / RytoeV
!
! Total number of bands of interest (usually the number of relevant conduction/valence bands)
nphband = phband_f - phband_i + 1
!
! Read the a2Fsave file
open(11,file=fil_a2F,status='unknown')
!
read(11,*) nbnd, nksfit
!
allocate(etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit))
!
! Read band structure and k-point data
read(11,*) etfit
read(11,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit)
read(11,*) wkfit
read(11,*) nk1fit, nk2fit, nk3fit
read(11,* ) nsym
do ns=1,nsym
read(11,*) ((s(i,j,ns),j=1,3),i=1,3)
enddo
!
close(11)
!
! Regular grid size
nkfit=nk1fit*nk2fit*nk3fit
!
! Read info file on k-points (and lattice)
open(11,file=fil_info,status='unknown')
!
read(11,*)
read(11,*) ((at(i,j),i=1,3),j=1,3)
!
read(11,*)
read(11,*)
!
read(11,*) ((bg(i,j),i=1,3),j=1,3)
!
close(11)
!
allocate (eqkfit(nkfit),sfit(nkfit))
allocate (iflag(nphband,nkfit),nkeff(nphband))
allocate (dfk(nphband,nkfit,3),vk(nphband,nkfit,3))
!
! Find the map between IBZ and the full-grid
call lint ( nsym, s, .true., at, bg, npk, 0,0,0, &
& nk1fit,nk2fit,nk3fit, nksfit, xkfit, 1, nkfit, eqkfit, sfit)
!
! Deallocate unnecessary variables
deallocate(wkfit,sfit,xkfit)
!
! Weights for regular grid
if ( lsoc .eqv. .true.) then
wk = 1.0/nkfit
else
wk = 2.0/nkfit
end if
!
! Call band velocities and forward derivatives
call vband_ibz( nk1fit,nk2fit,nk3fit,nphband,nksfit,etfit(phband_i:phband_f,:),eqkfit,at, vk, dfk)
!
! Include the 2pi/a factor
vk = vk / tpi * alat
!
! Determine the reduced grid for each band
call reducegrid(nkfit,nksfit,nphband,etfit(phband_i:phband_f,:),eqkfit,efermi, &
& cut,T, nkeff,iflag)
!
! Initiate Fermi integrals
I0 = 0.0
I1 = 0.0
I2 = 0.0
invT = 1.0/T
!
!$ t0 = omp_get_wtime()
! Loop over bands and kpoints
do ibnd=phband_i,phband_f
!
ibnd_ph = ibnd - phband_i + 1
!
!$omp parallel do default(shared) &
!$omp private(ik,ikk,ind_k,i,j,fd,dfd,fac) &
!$omp reduction(+: I0, I1, I2)
do ik=1,nkeff(ibnd_ph)
!
ikk = iflag(ibnd_ph,ik) ! ikk is in full-grid (just reduced)
ind_k = eqkfit(ikk) ! ind_k is in IBZ
!
! Fermi factors
fac = etfit(ibnd,ind_k) - efermi
fd = 1.0/( exp(fac*invT) + 1.0 )
dfd = invT * fd * (1.0 - fd)
!
! Compute I0 (related to Thomas-Fermi screening)
I0 = I0 + wk * dfd
!
! Compute I1, I2
do i=1,3
do j=1,3
I1(i,j) = I1(i,j) + wk * dfd / invtau * vk(ibnd_ph,ikk,i) &
& * vk(ibnd_ph,ikk,j)
!
I2(i,j) = I2(i,j) + wk * dfd * fac / invtau * vk(ibnd_ph,ikk,i) &
& * vk(ibnd_ph,ikk,j)
!
end do ! j
end do ! i
!
end do ! ik
!$omp end parallel do
!
end do ! ibnd
!
!Total integration time
!$ t0 = omp_get_wtime() - t0
!
! Conductivity (convert to units of Ohm^-1 cm^-1)
sig = I1 * convsig / vol
!
!Compute Seebeck
Se = 0.0
call inv33(I1, inv_I1)
do i=1,3
do j=1,3
do k=1,3
Se(i,j) = Se(i,j) + I2(i,k)*inv_I1(k,j)
end do
end do
end do
! Units (convert to units of V/K)
Se = Se * (-RytoeV * KtoRy / T)
!
! Write D(Ef), Conductivity and Seebeck into file
open(11,file='sig_1.out',status='unknown')
open(12,file='Se_1.out',status='unknown')
open(13,file='Def_1.out',status='unknown')
write(11,"(10e14.6)") efermi * RytoeV, ((sig(i,j),i=1,3),j=1,3)
write(12,"(10e14.6)") efermi * RytoeV, ((Se(i,j),i=1,3),j=1,3)
write(12,"(2e14.6)") efermi * RytoeV, I0 / RytoeV
close(11)
close(12)
close(13)
!
!
open(11,file='report_1',status='unknown')
if ( t0 > 0 ) then
write(11,"(A,I2)") "Number of threads = ", nthreads
write(11,"(A,e14.6)") "Integration time(s) =", t0
end if
!
write(11,"(A,I5)") "Number of kpoints (regular grid) = ", nkfit
do ibnd=1,nphband
write(11,"(A,2I5)") "band, reduced grid size", ibnd, nkeff(ibnd)
end do
close(11)
! Free memory
deallocate(etfit,eqkfit,dfk,vk,nkeff,iflag)
contains
!
! Determinant of a 3x3 matrix
double precision function det33 (a)
implicit none
double precision, intent(in) :: a(3,3)
! Determinant of a
det33 = a(1,1)*( a(2,2)*a(3,3)-a(2,3)*a(3,2) ) - &
& a(1,2)*( a(1,1)*a(3,3)-a(1,3)*a(3,1) ) + &
& a(1,3)*( a(2,1)*a(3,2)-a(2,2)*a(3,1) )
end function det33
!
!
! Inverse of a 3x3 matrix (analytical)
subroutine inv33(a, inv_a)
implicit none
double precision, intent(in) :: a(3,3)
double precision, intent(out) :: inv_a(3,3)
double precision :: deta_inv
deta_inv = 1.0 / det33(a)
inv_a(1,1) = deta_inv * ( a(2,2)*a(3,3)-a(2,3)*a(3,2) )
inv_a(1,2) = deta_inv * ( a(1,3)*a(3,2)-a(1,2)*a(3,3) )
inv_a(1,3) = deta_inv * ( a(1,2)*a(2,3)-a(1,3)*a(2,2) )
inv_a(2,1) = deta_inv * ( a(2,3)*a(3,1)-a(2,1)*a(3,3) )
inv_a(2,2) = deta_inv * ( a(1,1)*a(3,3)-a(1,3)*a(3,1) )
inv_a(2,3) = deta_inv * ( a(1,3)*a(2,1)-a(1,1)*a(2,3) )
inv_a(3,1) = deta_inv * ( a(2,1)*a(3,2)-a(2,2)*a(3,1) )
inv_a(3,2) = deta_inv * ( a(1,2)*a(3,1)-a(1,1)*a(3,2) )
inv_a(3,3) = deta_inv * ( a(1,1)*a(2,2)-a(1,2)*a(2,1) )
end subroutine inv33
end program fermi_int_1

111
PP/simple_transport/lint.f90 Normal file
View File

@ -0,0 +1,111 @@
subroutine lint ( nsym, s, minus_q, at, bg, npk, k1,k2,k3, &
& nk1,nk2,nk3, nks, xk, kunit, nkBZ, eqBZ, sBZ)
implicit none
integer, intent (IN) :: nks, nsym, s(3,3,48), npk, k1, k2, k3, &
nk1, nk2, nk3, kunit, nkBZ
logical, intent (IN) :: minus_q ! use .true.
!
double precision, intent(IN):: at(3,3), bg(3,3), xk(3,npk)
integer, intent(OUT) :: eqBZ(nkBZ), sBZ(nkBZ)
!
double precision :: xkr(3), deltap(3), deltam(3)
double precision, parameter:: eps=1.0d-4
!
double precision, allocatable :: xkg(:,:), xp(:,:)
integer :: i,j,k, ns, n, nk
integer :: nkh
!
! Re-generate a uniform grid of k-points xkg
!
allocate (xkg( 3,nkBZ))
!
! if(kunit < 1 .or. kunit > 2) call errore('lint','bad kunit value',kunit)
!
! kunit=2: get only "true" k points, not k+q points, from the list
!
!
nkh = nks/kunit
!
allocate (xp(3,nkh))
!
if (kunit == 1) then
xp(:,1:nkh) = xk(:,1:nkh)
else
do j=1,nkh
xp(:,j) = xk(:,2*j-1)
enddo
end if
!write(6,*) 'nkh ', 'nkBZ ', 'nsym'
!write(6,*) nkh, nkBZ, nsym
!
! Generate a uniform mesh
!
do i=1,nk1
do j=1,nk2
do k=1,nk3
n = (k-1) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
xkg(1,n) = dble(i-1)/nk1 + dble(k1)/2/nk1
xkg(2,n) = dble(j-1)/nk2 + dble(k2)/2/nk2
xkg(3,n) = dble(k-1)/nk3 + dble(k3)/2/nk3
end do
end do
end do
!
!
call cryst_to_cart (nkh,xp,at,-1)
!
!Maybe add here some OpenMP clauses
do nk=1,nkBZ
do n=1,nkh
do ns=1,nsym
do i=1,3
!
xkr(i) = s(i,1,ns) * xp(1,n) + s(i,2,ns) * xp(2,n) + &
& s(i,3,ns) * xp(3,n)
!
end do
!
do i=1,3
deltap(i) = xkr(i)-xkg(i,nk) - nint (xkr(i)-xkg(i,nk))
deltam(i) = xkr(i)+xkg(i,nk) - nint (xkr(i)+xkg(i,nk))
end do
!
if ( sqrt ( deltap(1)**2 + deltap(2)**2 + &
& deltap(3)**2 ) < eps .or. ( minus_q .and. &
& sqrt ( deltam(1)**2 + deltam(2)**2 + &
& deltam(3)**2 ) < eps ) ) then
!
eqBZ(nk) = n
sBZ(nk) = ns
!
go to 15
end if
!
end do
end do
!
!call errore('lint','cannot locate k point xk',nk)
write(6,*) 'somethings wrong', nk
!
15 continue
end do
!
do n=1,nkh
do nk=1,nkBZ
if (eqBZ(nk) == n) go to 20
end do
! this failure of the algorithm may indicate that the displaced grid
! (with k1,k2,k3.ne.0) does not have the full symmetry of the lattice
!
!call errore('lint','cannot remap grid on k-point list',n)
write(6,*) 'somethings wrong-2', nk, n
20 continue
end do
deallocate(xkg)
deallocate(xp)
return
end subroutine lint

20
PP/simple_transport/makefile Normal file
View File

@ -0,0 +1,20 @@
#
#
F90 = gfortran
FFLAGS = -fopenmp
FFLAGS2 = -fbounds-check
ef: ef.f90
${F90} ${FFLAGS} cryst_to_car.f90 lint.f90 ef.f90 -o ef.x
dos: dos.f90
${F90} ${FFLAGS} smearing_mod.f90 cryst_to_car.f90 lint.f90 vband_ibz.f90 dos.f90 -o dos.x
tauef: tauef.f90
${F90} ${FFLAGS2} smearing_mod.f90 cryst_to_car.f90 lint.f90 vband_ibz.f90 kplusq.f90 reducegrid.f90 invtau_nk.f90 tauef.f90 -o tauef.x
fermi_int_0: fermi_int_0.f90
${F90} ${FFLAGS} cryst_to_car.f90 lint.f90 vband_ibz.f90 fermi_int_0.f90 -o fermi_int_0.x
fermi_int_1: fermi_int_1.f90
${F90} ${FFLAGS} cryst_to_car.f90 reducegrid.f90 lint.f90 vband_ibz.f90 fermi_int_1.f90 -o fermi_int_1.x

51
PP/simple_transport/reducegrid.f90 Normal file
View File

@ -0,0 +1,51 @@
! This subroutine generates a reduced grid around the Fermi level
! for efficient calculation of transport integrals
!
subroutine reducegrid (nktot,nks,nbnd,etk,ekq,ef,cut,deg,nkeff,iflag)
implicit none
integer, intent(in) :: nktot, nks, nbnd
! nktot : total number of k-points (full grid)
! nks : total number of k-points (IBZ)
! nbnd : number of bands of interest
double precision, intent(in) :: etk(nbnd,nks),ef,deg,cut
! etk : energy eigenvalues
! ef, deg : Fermi level and smearing (or temperature)
! cut : scanning parameter (recommended value is ~ 10)
integer, intent(in) :: ekq(nktot)
! Map of IBZ to full-grid
integer, intent(out) :: nkeff(nbnd), iflag(nbnd,nktot)
! nkeff : effective reduced grid size for each band of interest
! iflag : index of the reduced grid point
! LOCAL parameters
integer :: ik,ikk,ibnd, ind_k
! Initiate nkeff
nkeff = 0
iflag = 0
do ibnd = 1, nbnd
!
ikk = 1
do ik = 1, nktot
!
ind_k = ekq(ik) ! map between IBZ and full grid
!
if ( abs(etk(ibnd,ind_k)-ef) < cut*deg) then
!
iflag(ibnd,ikk) = ik
nkeff(ibnd) = nkeff(ibnd) + 1
ikk = ikk + 1
!
end if
end do ! ik
end do ! ibnd
end subroutine reducegrid

86
PP/simple_transport/smearing_mod.f90 Normal file
View File

@ -0,0 +1,86 @@
! Author: Burak Himmetoglu (burakhmmtgl@gmail.com)
! This module contains smearing methods.
!
! Methods:
! w0gauss : Standard Gaussian smearing
! w1gauss : FD smearing for n=0, Integral of Gaussian if n=1
! sig_nk : Variable smearing for Dirac-delta involving e_{nk} - e_{m,k+q}
! sig0 : Variable smearing for Dirac-delta involving e_{nk}
module smearing_mod
implicit none
double precision, parameter :: sqrtpm1 = 1.0d0/1.7724538509055160d0
contains
double precision function w0gauss(x)
implicit none
double precision, intent(in) :: x
double precision :: arg
arg = min(200.d0,x**2.d0)
w0gauss = sqrtpm1 * exp( -arg )
end function w0gauss
!
!
double precision function w1gauss(x,ngauss)
implicit none
double precision, intent(in) :: x
integer, intent(in) :: ngauss
! ngauss=0 (gaussian smearing)
! ngauss=1 FD
if ( ngauss .eq. 0 ) then
w1gauss = 0.5d0 * ( 1.d0 + derf( x * sqrtpm1 ) )
else if ( ngauss .eq. 1 ) then
w1gauss = 1.d0 / ( exp(-x) + 1.d0)
end if
end function w1gauss
!
!
double precision function sig_nk(nk1,nk2,nk3,vk,vkq,aa)
implicit none
integer, intent(in) :: nk1, nk2, nk3
double precision, intent(in) :: vk(3), vkq(3), aa
double precision :: dF1, dF2, dF3!, dFarr(4)
dF1 = ( vk(1) - vkq(1) ) * 1.0/nk1
dF2 = ( vk(2) - vkq(2) ) * 1.0/nk2
dF3 = ( vk(3) - vkq(3) ) * 1.0/nk3
!
!dFarr(1) = abs(dF1+dF2+dF3)
!dFarr(2) = abs(-dF1+dF2+dF3)
!dFarr(3) = abs(dF1+dF2-dF3)
!dFarr(4) = abs(dF1-dF2+dF3)
!
!sig_nk = aa * maxval(dFarr)
sig_nk = aa * sqrt(abs(dF1**2+dF2**2+dF3**2))
end function sig_nk
!
!
double precision function sig0(nk1,nk2,nk3,vk,aa)
implicit none
integer, intent(in) :: nk1, nk2, nk3
double precision, intent(in) :: vk(3), aa
double precision :: dF1, dF2, dF3!, dFarr(4)
dF1 = vk(1) * 1.0/nk1
dF2 = vk(2) * 1.0/nk2
dF3 = vk(3) * 1.0/nk3
!
!dFarr(1) = abs(dF1+dF2+dF3)
!dFarr(2) = abs(-dF1+dF2+dF3)
!dFarr(3) = abs(dF1+dF2-dF3)
!dFarr(4) = abs(dF1-dF2+dF3)
!
!sig0 = aa * maxval(dFarr)
sig0 = aa * sqrt(abs(dF1**2+dF2**2+dF3**2))
end function sig0
end module smearing_mod

72
PP/simple_transport/smearing_mod.mod Normal file
View File

@ -0,0 +1,72 @@
GFORTRAN module version '10' created from smearing_mod.f90
MD5:109d4839832f6dee33625924ab198f78 -- If you edit this, you'll get what you deserve.
(() () () () () () () () () () () () () () () () () () () () () () () ()
() () ())
()
()
()
()
()
(2 'sig0' 'smearing_mod' '' 1 ((PROCEDURE UNKNOWN-INTENT MODULE-PROC
DECL UNKNOWN 0 0 FUNCTION IMPLICIT_PURE) (REAL 8 0 0 0 REAL ()) 3 0 (4 5
6 7 8) () 2 () () () 0 0)
9 'sig_nk' 'smearing_mod' '' 1 ((PROCEDURE UNKNOWN-INTENT MODULE-PROC
DECL UNKNOWN 0 0 FUNCTION IMPLICIT_PURE) (REAL 8 0 0 0 REAL ()) 10 0 (
11 12 13 14 15 16) () 9 () () () 0 0)
17 'smearing_mod' 'smearing_mod' '' 1 ((MODULE UNKNOWN-INTENT
UNKNOWN-PROC UNKNOWN UNKNOWN 0 0) (UNKNOWN 0 0 0 0 UNKNOWN ()) 0 0 () ()
0 () () () 0 0)
18 'sqrtpm1' 'smearing_mod' '' 1 ((PARAMETER UNKNOWN-INTENT UNKNOWN-PROC
UNKNOWN IMPLICIT-SAVE 0 0) (REAL 8 0 0 0 REAL ()) 0 0 () (CONSTANT (
REAL 8 0 0 0 REAL ()) 0 '0.906eba8214db68@0') () 0 () () () 0 0)
19 'w0gauss' 'smearing_mod' '' 1 ((PROCEDURE UNKNOWN-INTENT MODULE-PROC
DECL UNKNOWN 0 0 FUNCTION IMPLICIT_PURE) (REAL 8 0 0 0 REAL ()) 20 0 (
21) () 19 () () () 0 0)
22 'w1gauss' 'smearing_mod' '' 1 ((PROCEDURE UNKNOWN-INTENT MODULE-PROC
DECL UNKNOWN 0 0 FUNCTION IMPLICIT_PURE) (REAL 8 0 0 0 REAL ()) 23 0 (
24 25) () 22 () () () 0 0)
4 'nk1' '' '' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
5 'nk2' '' '' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
6 'nk3' '' '' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
7 'vk' '' '' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DIMENSION
DUMMY) (REAL 8 0 0 0 REAL ()) 0 0 () (1 0 EXPLICIT (CONSTANT (INTEGER 4
0 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 0 INTEGER ()) 0 '3'))
0 () () () 0 0)
8 'aa' '' '' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
REAL 8 0 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
11 'nk1' '' '' 10 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
(INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
12 'nk2' '' '' 10 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
(INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
13 'nk3' '' '' 10 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
(INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
14 'vk' '' '' 10 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
DIMENSION DUMMY) (REAL 8 0 0 0 REAL ()) 0 0 () (1 0 EXPLICIT (CONSTANT (
INTEGER 4 0 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 0 INTEGER ())
0 '3')) 0 () () () 0 0)
15 'vkq' '' '' 10 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
DIMENSION DUMMY) (REAL 8 0 0 0 REAL ()) 0 0 () (1 0 EXPLICIT (CONSTANT (
INTEGER 4 0 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 0 INTEGER ())
0 '3')) 0 () () () 0 0)
16 'aa' '' '' 10 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
REAL 8 0 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
21 'x' '' '' 20 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
REAL 8 0 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
24 'x' '' '' 23 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
REAL 8 0 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
25 'ngauss' '' '' 23 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
DUMMY) (INTEGER 4 0 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
)
('sig0' 0 2 'sig_nk' 0 9 'smearing_mod' 0 17 'sqrtpm1' 0 18 'w0gauss' 0
19 'w1gauss' 0 22)

109
PP/simple_transport/vband_ibz.f90 Normal file
View File

@ -0,0 +1,109 @@
subroutine vband_ibz ( nk1,nk2,nk3, nbnd, nksfit, etk, eqk, at, vk, dfk )
!
! Written by Burak Himmetoglu (burakhmmtgl@gmail.com)
! Uses some parts of the PW distribution
!
! Computes band velocities and forward derivatives using finite
! differences.
!
!$ use omp_lib
implicit none
integer, intent(in) :: nk1,nk2,nk3, nbnd, nksfit, eqk(nk1*nk2*nk3)
double precision, intent(in) :: etk(nbnd,nksfit), at(3,3)
double precision, intent(out) :: vk(nbnd,nk1*nk2*nk3,3), dfk(nbnd,nk1*nk2*nk3,3)
! vk is band velocity
! dfk is forward derivative (used in adaptive smearing)
integer :: i,j,k,n,ik,ibnd,nktot,nm,np
double precision :: temp, temp2, vaux(3,nk1*nk2*nk3)
! Create the k-points, then compute velocities
nktot = nk1*nk2*nk3
!$omp parallel do default(none) &
!$omp collapse(4) &
!$omp private(ibnd,i,j,k,n,np,nm) &
!$omp shared(nbnd,nk1,nk2,nk3,vk,dfk,etk,eqk)
do ibnd=1,nbnd
do i=1,nk1
do j=1,nk2
do k=1,nk3
n = (k-1) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
!
! velocity along k1
!
if ( i .eq. 1 ) then
np = (k-1) + (j-1)*nk3 + (i)*nk2*nk3 + 1
nm = (k-1) + (j-1)*nk3 + (i-2+nk1)*nk2*nk3 + 1
else if (i .eq. nk1) then
np = (k-1) + (j-1)*nk3 + (i-nk1)*nk2*nk3 + 1
nm = (k-1) + (j-1)*nk3 + (i-2)*nk2*nk3 + 1
else
np = (k-1) + (j-1)*nk3 + (i)*nk2*nk3 + 1
nm = (k-1) + (j-1)*nk3 + (i-2)*nk2*nk3 + 1
end if
vk(ibnd,n,1) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(nm)) )/(2.0/nk1)
dfk(ibnd,n,1) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(n)) )/(2.0/nk1)
!
! velocity along k2
!
if ( j .eq. 1 ) then
np = (k-1) + (j)*nk3 + (i-1)*nk2*nk3 + 1
nm = (k-1) + (j-2+nk2)*nk3 + (i-1)*nk2*nk3 + 1
else if (j .eq. nk2) then
np = (k-1) + (j-nk2)*nk3 + (i-1)*nk2*nk3 + 1
nm = (k-1) + (j-2)*nk3 + (i-1)*nk2*nk3 + 1
else
np = (k-1) + (j)*nk3 + (i-1)*nk2*nk3 + 1
nm = (k-1) + (j-2)*nk3 + (i-1)*nk2*nk3 + 1
end if
vk(ibnd,n,2) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(nm)) )/(2.0/nk2)
dfk(ibnd,n,2) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(n)) )/(2.0/nk2)
!
! velocity along k3
!
if ( k .eq. 1 ) then
np = (k) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
nm = (k-2+nk3) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
else if (k .eq. nk3) then
np = (k-nk3) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
nm = (k-2) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
else
np = (k) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
nm = (k-2) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1
end if
vk(ibnd,n,3) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(nm)) )/(2.0/nk3)
dfk(ibnd,n,3) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(n)) )/(2.0/nk3)
!
end do ! k
end do ! j
end do ! i
end do ! ibnd
!$omp end parallel do
!
! Convert the derivatives in crystal coordinates to cartesian coordinates
!
do ibnd=1,nbnd
!
do j=1,3
vaux(j,:) = vk(ibnd,:,j) ! temporary vector
end do
!
! vaux is updated after this call
!
call cryst_to_cart (nktot,vaux,at,1)
!
do j=1,3
vk(ibnd,:,j) = vaux(j,:)
end do
!
end do
!
end subroutine vband_ibz