Computes the following quantities: (i) D tensor (ii) K tensor (iii) C tensor (iv) current-induced optical activity (v) natural optical activity
subroutine gyrotropic_main
!============================================================!
! !
!! Computes the following quantities:
!! (i) D tensor
!! (ii) K tensor
!! (iii) C tensor
!! (iv) current-induced optical activity
!! (v) natural optical activity
! !
!============================================================!
use w90_constants, only: dp, cmplx_0, elem_charge_SI, hbar_SI, &
eV_au, bohr, elec_mass_SI, twopi, eps0_SI
use w90_comms, only: on_root, num_nodes, my_node_id, comms_reduce
use w90_utility, only: utility_det3
use w90_io, only: io_error, stdout, io_file_unit, seedname, &
io_stopwatch
use w90_postw90_common, only: nrpts, irvec, num_int_kpts_on_node, int_kpts, &
weight
use w90_parameters, only: timing_level, iprint, num_wann, gyrotropic_kmesh, &
cell_volume, transl_inv, gyrotropic_task, &
gyrotropic_nfreq, gyrotropic_freq_list, nfermi, &
fermi_energy_list, gyrotropic_box, gyrotropic_box_corner, spinors
use w90_get_oper, only: get_HH_R, get_AA_R, get_BB_R, get_CC_R, &
get_SS_R
real(kind=dp), allocatable :: gyro_K_spn(:, :, :)
real(kind=dp), allocatable :: gyro_DOS(:)
real(kind=dp), allocatable :: gyro_K_orb(:, :, :)
real(kind=dp), allocatable :: gyro_C(:, :, :)
real(kind=dp), allocatable :: gyro_D(:, :, :)
real(kind=dp), allocatable :: gyro_Dw(:, :, :, :)
real(kind=dp), allocatable :: gyro_NOA_spn(:, :, :, :)
real(kind=dp), allocatable :: gyro_NOA_orb(:, :, :, :)
character(len=30) :: f_out_name_tmp
character(len=30) :: units_tmp
character(len=120) :: comment_tmp
real(kind=dp) :: kweight, kpt(3), &
db1, db2, db3, fac, freq
integer :: n, i, j, k, ikpt, if, ierr, loop_x, loop_y, loop_z, &
loop_xyz, ifreq, &
file_unit
logical :: eval_K, eval_C, eval_D, eval_Dw, eval_NOA, eval_spn, eval_DOS
if (nfermi == 0) call io_error( &
'Must specify one or more Fermi levels when gyrotropic=true')
if (timing_level > 1 .and. on_root) call io_stopwatch('gyrotropic: prelims', 1)
! Mesh spacing in reduced coordinates
!
db1 = 1.0_dp/real(gyrotropic_kmesh(1), dp)
db2 = 1.0_dp/real(gyrotropic_kmesh(2), dp)
db3 = 1.0_dp/real(gyrotropic_kmesh(3), dp)
eval_K = .false.
eval_C = .false.
eval_D = .false.
eval_Dw = .false.
eval_spn = .false.
eval_NOA = .false.
eval_DOS = .false.
if (index(gyrotropic_task, '-k') > 0) eval_K = .true.
if (index(gyrotropic_task, '-c') > 0) eval_C = .true.
if (index(gyrotropic_task, '-d0') > 0) eval_D = .true.
if (index(gyrotropic_task, '-dw') > 0) eval_Dw = .true.
if (index(gyrotropic_task, '-spin') > 0) eval_spn = .true.
if (index(gyrotropic_task, '-noa') > 0) eval_NOA = .true.
if (index(gyrotropic_task, '-dos') > 0) eval_DOS = .true.
if (index(gyrotropic_task, 'all') > 0) then
eval_K = .true.
eval_C = .true.
eval_D = .true.
eval_Dw = .true.
if (spinors) eval_spn = .true.
eval_NOA = .true.
eval_DOS = .true.
endif
if (.not. (eval_K .or. eval_noa)) eval_spn = .false.
if ((.not. spinors) .and. eval_spn) call io_error( &
"spin contribution requested for gyrotropic, but the wavefunctions are not spinors")
! Wannier matrix elements, allocations and initializations
call get_HH_R
if (eval_D .or. eval_Dw .or. eval_K .or. eval_NOA) then
call get_AA_R
endif
if (eval_spn) then
call get_SS_R
endif
if (eval_K) then
call get_BB_R
call get_CC_R
allocate (gyro_K_orb(3, 3, nfermi))
gyro_K_orb = 0.0_dp
if (eval_spn) then
allocate (gyro_K_spn(3, 3, nfermi))
gyro_K_spn = 0.0_dp
endif
endif
if (eval_D) then
allocate (gyro_D(3, 3, nfermi))
gyro_D = 0.0_dp
endif
if (eval_DOS) then
allocate (gyro_DOS(nfermi))
gyro_DOS = 0.0_dp
endif
if (eval_C) then
allocate (gyro_C(3, 3, nfermi))
gyro_C = 0.0_dp
endif
if (eval_Dw) then
allocate (gyro_Dw(3, 3, nfermi, gyrotropic_nfreq))
gyro_Dw = 0.0_dp
endif
if (eval_NOA) then
allocate (gyro_NOA_orb(3, 3, nfermi, gyrotropic_nfreq))
gyro_NOA_orb = 0.0_dp
if (eval_spn) then
allocate (gyro_NOA_spn(3, 3, nfermi, gyrotropic_nfreq))
gyro_NOA_spn = 0.0_dp
endif
endif
if (on_root) then
flush(stdout)
write (stdout, '(/,/,1x,a)') 'Properties calculated in module g y r o t r o p i c'
write (stdout, '(1x,a)') '------------------------------------------'
if (eval_D) write (stdout, '(/,3x,a)') '* D-tensor --- Eq.2 of TAS17 '
if (eval_dos) write (stdout, '(/,3x,a)') '* density of states '
if (eval_K) then
write (stdout, '(/,3x,a)') '* K-tensor --- Eq.3 of TAS17 '
if (eval_spn) then
write (stdout, '(3x,a)') ' * including spin component '
else
write (stdout, '(3x,a)') ' * excluding spin component '
endif
endif
if (eval_Dw) write (stdout, '(/,3x,a)') '* Dw-tensor --- Eq.12 of TAS17 '
if (eval_C) write (stdout, '(/,3x,a)') '* C-tensor --- Eq.B6 of TAS17 '
if (eval_NOA) then
write (stdout, '(/,3x,a)') '* gamma-tensor of NOA --- Eq.C12 of TAS17 '
if (eval_spn) then
write (stdout, '(3x,a)') ' * including spin component '
else
write (stdout, '(3x,a)') ' * excluding spin component '
endif
endif
if (transl_inv) then
if (eval_K) &
call io_error('transl_inv=T disabled for K-tensor')
write (stdout, '(/,1x,a)') &
'Using a translationally-invariant discretization for the'
write (stdout, '(1x,a)') &
'band-diagonal Wannier matrix elements of r, etc.'
endif
if (timing_level > 1) then
call io_stopwatch('gyrotropic: prelims', 2)
call io_stopwatch('gyrotropic: k-interpolation', 1)
endif
write (stdout, '(1x,a20,3(i0,1x))') 'Interpolation grid: ', gyrotropic_kmesh(1:3)
flush(stdout)
end if !on_root
! Do not read 'kpoint.dat'. Loop over a regular grid in the full BZ
kweight = db1*db2*db3*utility_det3(gyrotropic_box)
do loop_xyz = my_node_id, PRODUCT(gyrotropic_kmesh) - 1, num_nodes
loop_x = loop_xyz/(gyrotropic_kmesh(2)*gyrotropic_kmesh(3))
loop_y = (loop_xyz - loop_x*(gyrotropic_kmesh(2) &
*gyrotropic_kmesh(3)))/gyrotropic_kmesh(3)
loop_z = loop_xyz - loop_x*(gyrotropic_kmesh(2)*gyrotropic_kmesh(3)) &
- loop_y*gyrotropic_kmesh(3)
kpt(1) = loop_x*db1
kpt(2) = loop_y*db2
kpt(3) = loop_z*db3
kpt(:) = gyrotropic_box_corner(:) + matmul(kpt, gyrotropic_box)
call gyrotropic_get_k_list(kpt, kweight, &
gyro_K_spn, gyro_K_orb, gyro_D, gyro_Dw, gyro_C, &
gyro_DOS, gyro_NOA_orb, gyro_NOA_spn, &
eval_K, eval_D, eval_Dw, eval_NOA, eval_spn, eval_C, eval_dos)
end do !loop_xyz
! Collect contributions from all nodes
!
if (eval_K) then
call comms_reduce(gyro_K_orb(1, 1, 1), 3*3*nfermi, 'SUM')
if (eval_spn) call comms_reduce(gyro_K_spn(1, 1, 1), 3*3*nfermi, 'SUM')
endif
if (eval_D) &
call comms_reduce(gyro_D(1, 1, 1), 3*3*nfermi, 'SUM')
if (eval_C) &
call comms_reduce(gyro_C(1, 1, 1), 3*3*nfermi, 'SUM')
if (eval_Dw) &
call comms_reduce(gyro_Dw(1, 1, 1, 1), 3*3*nfermi*gyrotropic_nfreq, 'SUM')
if (eval_dos) &
call comms_reduce(gyro_DOS(1), nfermi, 'SUM')
if (eval_NOA) then
call comms_reduce(gyro_NOA_orb(1, 1, 1, 1), 3*3*nfermi*gyrotropic_nfreq, 'SUM')
if (eval_spn) call comms_reduce(gyro_NOA_spn(1, 1, 1, 1), 3*3*nfermi*gyrotropic_nfreq, 'SUM')
endif
if (on_root) then
if (timing_level > 1) call io_stopwatch('gyrotropic: k-interpolation', 2)
write (stdout, '(1x,a)') ' '
write (stdout, *) 'Calculation finished, writing results'
flush(stdout)
if (eval_K) then
if (eval_spn) then
! At this point gme_spn_list contains
! (1/N) sum_k delta(E_kn-E_f).(d E_{kn}/d k_i).sigma_{kn,j}
! (units of length) in Angstroms.
! ====================================
! To get K in units of Ampere do the following:
! ====================================
! * Divide by V_c in Ang^3 to get a quantity with units of [L]^{-2}
! * Multiply by 10^20 to convert to SI
! * Multiply by -g_s.e.hbar/(4m_e) \simeq e.hbar/(2.m_e) in SI units
! ==============================
! fac = 10^20*e*hbar/(2.m_e.V_c)
! ==============================
fac = -1.0e20_dp*elem_charge_SI*hbar_SI/(2.*elec_mass_SI*cell_volume)
gyro_K_spn(:, :, :) = gyro_K_spn(:, :, :)*fac
f_out_name_tmp = 'K_spin'
units_tmp = "Ampere"
comment_tmp = "spin part of the K tensor -- Eq. 3 of TAS17"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_K_spn, units=units_tmp, &
comment=comment_tmp)
endif ! eval_K && eval_spin
! At this point gme_orb_list contains
! (1/N)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i)
! .Im[<del_k u_kn| x (H_k-E_kn)|del_k u_kn>]
! (units of energy times length^3) in eV.Ang^3.
! ====================================
! To get K in units of Ampere do the following:
! ====================================
! * Divide by V_c in Ang^3 to get a quantity with units of eV
! * Multiply by 'e' in SI to convert to SI (Joules)
! * Multiply by e/(2.hbar) to get K in Ampere
! ====================================
! fac = e^2/(2.hbar.V_c)
! ====================================
fac = elem_charge_SI**2/(2.*hbar_SI*cell_volume)
gyro_K_orb(:, :, :) = gyro_K_orb(:, :, :)*fac
f_out_name_tmp = 'K_orb'
units_tmp = "Ampere"
comment_tmp = "orbital part of the K tensor -- Eq. 3 of TAS17"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_K_orb, units=units_tmp, &
comment=comment_tmp)
endif ! eval_K
if (eval_D) then
fac = 1./cell_volume
gyro_D(:, :, :) = gyro_D(:, :, :)*fac
f_out_name_tmp = 'D'
units_tmp = "dimensionless"
comment_tmp = "the D tensor -- Eq. 2 of TAS17"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_D, units=units_tmp, &
comment=comment_tmp)
endif
if (eval_Dw) then
fac = 1./cell_volume
gyro_Dw(:, :, :, :) = gyro_Dw(:, :, :, :)*fac
f_out_name_tmp = 'tildeD'
units_tmp = "dimensionless"
comment_tmp = "the tildeD tensor -- Eq. 12 of TAS17"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEfW=gyro_Dw, units=units_tmp, &
comment=comment_tmp)
endif
if (eval_C) then
! At this point gyro_C contains
! (1/N)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j)
! (units of energy*length^2) in eV*Ang^2
!
! To get it in Cab = e/h * (1/N*V_cell)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j)
! in units Ampere/cm
!
! divide by V_c in Ang^3 to get eV/Ang
! multiply by 10^8*e in SI to get J/cm
! multiply by e/h in SI
!
fac = 1.0e+8_dp*elem_charge_SI**2/(twopi*hbar_SI*cell_volume)
gyro_C(:, :, :) = gyro_C(:, :, :)*fac
f_out_name_tmp = 'C'
units_tmp = "Ampere/cm"
comment_tmp = "the C tensor -- Eq. B6 of TAS17"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_C, units=units_tmp, &
comment=comment_tmp)
endif
if (eval_noa) then
! at this point gyro_NOA_orb is in eV^-1.Ang^3 !
! We want the result in angstrems !
! * Divide by V_c in Ang^3 to make it eV^{-1}
! * Divide by e in SI to get J^{-1}
! * multiply by e^2/eps_0 to get meters
! *multiply dy 1e10 to get Ang
fac = 1e+10_dp*elem_charge_SI/(cell_volume*eps0_SI)
gyro_NOA_orb = gyro_NOA_orb*fac
f_out_name_tmp = 'NOA_orb'
units_tmp = "Ang"
comment_tmp = "the tensor $gamma_{abc}^{orb}$ (Eq. C12,C14 of TAS17)"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEfW=gyro_NOA_orb, units=units_tmp, &
comment=comment_tmp, symmetrize=.false.)
if (eval_spn) then
! at this point gyro_NOA_spn is in eV^-2.Ang !
! We want the result in angstrems !
! * Divide by V_c in Ang^3 to make it (eV.Ang)^{-2}
! * multiply by 1e20/e^2 in SI to get (J.m)^{-2}
! * multiply by e^2/eps_0 to get (J.m)^{-1}
! *multiply dy hbar^2/m_e to get m
! *multiply by 1e10 to get Ang
fac = 1e+30_dp*hbar_SI**2/(cell_volume*eps0_SI*elec_mass_SI)
gyro_NOA_spn = gyro_NOA_spn*fac
f_out_name_tmp = 'NOA_spin'
units_tmp = "Ang"
comment_tmp = "the tensor $gamma_{abc}^{spin}$ (Eq. C12,C15 of TAS17)"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEfW=gyro_NOA_spn, units=units_tmp, &
comment=comment_tmp, symmetrize=.false.)
endif
endif !eval_NOA
if (eval_DOS) then
! At this point gyro_C contains
! (1/N)sum_{k,n} delta(E_kn-E_f)
! in units of eV^{-1}
! divide by V_c in Ang^3 to get units 1./(eV*Ang^3)
gyro_DOS(:) = gyro_DOS(:)/cell_volume
f_out_name_tmp = 'DOS'
units_tmp = "eV^{-1}.Ang^{-3}"
comment_tmp = "density of states"
call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf1d=gyro_DOS, units=units_tmp, &
comment=comment_tmp)
endif
end if !on_root
end subroutine gyrotropic_main