Computes the following quantities: (i) Anomalous Hall conductivity (from Berry curvature) (ii) Complex optical conductivity (Kubo-Greenwood) & JDOS (iii) Orbital magnetization (iv) Nonlinear shift current (v) Spin Hall conductivity
subroutine berry_main
!============================================================!
! !
!! Computes the following quantities:
!! (i) Anomalous Hall conductivity (from Berry curvature)
!! (ii) Complex optical conductivity (Kubo-Greenwood) & JDOS
!! (iii) Orbital magnetization
!! (iv) Nonlinear shift current
!! (v) Spin Hall conductivity
! !
!============================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i, elem_charge_SI, hbar_SI, &
eV_au, bohr, pi, eV_seconds
use w90_comms, only: on_root, num_nodes, my_node_id, comms_reduce
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, berry_kmesh, &
berry_curv_adpt_kmesh, &
berry_curv_adpt_kmesh_thresh, &
wanint_kpoint_file, cell_volume, transl_inv, &
berry_task, berry_curv_unit, spin_decomp, &
kubo_nfreq, kubo_freq_list, nfermi, &
fermi_energy_list, shc_freq_scan, &
kubo_adpt_smr, kubo_adpt_smr_fac, &
kubo_adpt_smr_max, kubo_smr_fixed_en_width, &
scissors_shift, num_valence_bands, &
shc_bandshift, shc_bandshift_firstband, shc_bandshift_energyshift
use w90_get_oper, only: get_HH_R, get_AA_R, get_BB_R, get_CC_R, &
get_SS_R, get_SHC_R
real(kind=dp), allocatable :: adkpt(:, :)
! AHC and orbital magnetization, calculated for a list of Fermi levels
!
! First index labels J0,J1,J2 terms, second labels the Cartesian component
!
real(kind=dp) :: imf_k_list(3, 3, nfermi), imf_list(3, 3, nfermi), imf_list2(3, 3, nfermi)
real(kind=dp) :: img_k_list(3, 3, nfermi), img_list(3, 3, nfermi)
real(kind=dp) :: imh_k_list(3, 3, nfermi), imh_list(3, 3, nfermi)
real(kind=dp) :: ahc_list(3, 3, nfermi)
real(kind=dp) :: LCtil_list(3, 3, nfermi), ICtil_list(3, 3, nfermi), &
Morb_list(3, 3, nfermi)
real(kind=dp) :: imf_k_list_dummy(3, 3, nfermi) ! adaptive refinement of AHC
! shift current
real(kind=dp), allocatable :: sc_k_list(:, :, :)
real(kind=dp), allocatable :: sc_list(:, :, :)
! Complex optical conductivity, dividided into Hermitean and
! anti-Hermitean parts
!
complex(kind=dp), allocatable :: kubo_H_k(:, :, :)
complex(kind=dp), allocatable :: kubo_H(:, :, :)
complex(kind=dp), allocatable :: kubo_AH_k(:, :, :)
complex(kind=dp), allocatable :: kubo_AH(:, :, :)
! decomposition into up-up, down-down and spin-flip transitions
complex(kind=dp), allocatable :: kubo_H_k_spn(:, :, :, :)
complex(kind=dp), allocatable :: kubo_H_spn(:, :, :, :)
complex(kind=dp), allocatable :: kubo_AH_k_spn(:, :, :, :)
complex(kind=dp), allocatable :: kubo_AH_spn(:, :, :, :)
! Joint density of states
!
real(kind=dp), allocatable :: jdos_k(:)
real(kind=dp), allocatable :: jdos(:)
! decomposition into up-up, down-down and spin-flip transitions
real(kind=dp), allocatable :: jdos_k_spn(:, :)
real(kind=dp), allocatable :: jdos_spn(:, :)
! Spin Hall conductivity
real(kind=dp), allocatable :: shc_fermi(:), shc_k_fermi(:)
complex(kind=dp), allocatable :: shc_freq(:), shc_k_freq(:)
! for fermi energy scan, adaptive kmesh
real(kind=dp), allocatable :: shc_k_fermi_dummy(:)
real(kind=dp) :: kweight, kweight_adpt, kpt(3), kpt_ad(3), &
db1, db2, db3, fac, freq, rdum, vdum(3)
integer :: n, i, j, k, jk, ikpt, if, ispn, ierr, loop_x, loop_y, loop_z, &
loop_xyz, loop_adpt, adpt_counter_list(nfermi), ifreq, &
file_unit
character(len=120) :: file_name
logical :: eval_ahc, eval_morb, eval_kubo, not_scannable, eval_sc, eval_shc
logical :: ladpt_kmesh
logical :: ladpt(nfermi)
if (nfermi == 0) call io_error( &
'Must specify one or more Fermi levels when berry=true')
if (timing_level > 1 .and. on_root) call io_stopwatch('berry: prelims', 1)
! Mesh spacing in reduced coordinates
!
db1 = 1.0_dp/real(berry_kmesh(1), dp)
db2 = 1.0_dp/real(berry_kmesh(2), dp)
db3 = 1.0_dp/real(berry_kmesh(3), dp)
eval_ahc = .false.
eval_morb = .false.
eval_kubo = .false.
eval_sc = .false.
eval_shc = .false.
if (index(berry_task, 'ahc') > 0) eval_ahc = .true.
if (index(berry_task, 'morb') > 0) eval_morb = .true.
if (index(berry_task, 'kubo') > 0) eval_kubo = .true.
if (index(berry_task, 'sc') > 0) eval_sc = .true.
if (index(berry_task, 'shc') > 0) eval_shc = .true.
! Wannier matrix elements, allocations and initializations
!
if (eval_ahc) then
call get_HH_R
call get_AA_R
imf_list = 0.0_dp
adpt_counter_list = 0
endif
if (eval_morb) then
call get_HH_R
call get_AA_R
call get_BB_R
call get_CC_R
imf_list2 = 0.0_dp
img_list = 0.0_dp
imh_list = 0.0_dp
endif
! List here berry_tasks that assume nfermi=1
!
not_scannable = eval_kubo .or. (eval_shc .and. shc_freq_scan)
if (not_scannable .and. nfermi .ne. 1) call io_error( &
'The berry_task(s) you chose require that you specify a single ' &
//'Fermi energy: scanning the Fermi energy is not implemented')
if (eval_kubo) then
call get_HH_R
call get_AA_R
allocate (kubo_H_k(3, 3, kubo_nfreq))
allocate (kubo_H(3, 3, kubo_nfreq))
allocate (kubo_AH_k(3, 3, kubo_nfreq))
allocate (kubo_AH(3, 3, kubo_nfreq))
allocate (jdos_k(kubo_nfreq))
allocate (jdos(kubo_nfreq))
kubo_H = cmplx_0
kubo_AH = cmplx_0
jdos = 0.0_dp
if (spin_decomp) then
call get_SS_R
allocate (kubo_H_k_spn(3, 3, 3, kubo_nfreq))
allocate (kubo_H_spn(3, 3, 3, kubo_nfreq))
allocate (kubo_AH_k_spn(3, 3, 3, kubo_nfreq))
allocate (kubo_AH_spn(3, 3, 3, kubo_nfreq))
allocate (jdos_k_spn(3, kubo_nfreq))
allocate (jdos_spn(3, kubo_nfreq))
kubo_H_spn = cmplx_0
kubo_AH_spn = cmplx_0
jdos_spn = 0.0_dp
endif
endif
if (eval_sc) then
call get_HH_R
call get_AA_R
allocate (sc_k_list(3, 6, kubo_nfreq))
allocate (sc_list(3, 6, kubo_nfreq))
sc_k_list = 0.0_dp
sc_list = 0.0_dp
endif
if (eval_shc) then
call get_HH_R
call get_AA_R
call get_SS_R
call get_SHC_R
if (shc_freq_scan) then
allocate (shc_freq(kubo_nfreq))
allocate (shc_k_freq(kubo_nfreq))
shc_freq = 0.0_dp
shc_k_freq = 0.0_dp
else
allocate (shc_fermi(nfermi))
allocate (shc_k_fermi(nfermi))
allocate (shc_k_fermi_dummy(nfermi))
shc_fermi = 0.0_dp
shc_k_fermi = 0.0_dp
!only used for fermiscan & adpt kmesh
shc_k_fermi_dummy = 0.0_dp
adpt_counter_list = 0
endif
endif
if (on_root) then
write (stdout, '(/,/,1x,a)') &
'Properties calculated in module b e r r y'
write (stdout, '(1x,a)') &
'------------------------------------------'
if (eval_ahc) write (stdout, '(/,3x,a)') &
'* Anomalous Hall conductivity'
if (eval_morb) write (stdout, '(/,3x,a)') '* Orbital magnetization'
if (eval_kubo) then
if (spin_decomp) then
write (stdout, '(/,3x,a)') &
'* Complex optical conductivity and its spin-decomposition'
write (stdout, '(/,3x,a)') &
'* Joint density of states and its spin-decomposition'
else
write (stdout, '(/,3x,a)') '* Complex optical conductivity'
write (stdout, '(/,3x,a)') '* Joint density of states'
endif
endif
if (eval_sc) write (stdout, '(/,3x,a)') &
'* Shift current'
if (eval_shc) then
write (stdout, '(/,3x,a)') '* Spin Hall Conductivity'
if (shc_freq_scan) then
write (stdout, '(/,3x,a)') ' Frequency scan'
else
write (stdout, '(/,3x,a)') ' Fermi energy scan'
endif
endif
if (transl_inv) then
if (eval_morb) &
call io_error('transl_inv=T disabled for morb')
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('berry: prelims', 2)
call io_stopwatch('berry: k-interpolation', 1)
endif
end if !on_root
! Set up adaptive refinement mesh
!
allocate (adkpt(3, berry_curv_adpt_kmesh**3), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating adkpt in berry')
ikpt = 0
!
! OLD VERSION (only works correctly for odd grids including original point)
!
! do i=-(berry_curv_adpt_kmesh-1)/2,(berry_curv_adpt_kmesh-1)/2
! do j=-(berry_curv_adpt_kmesh-1)/2,(berry_curv_adpt_kmesh-1)/2
! do k=-(berry_curv_adpt_kmesh-1)/2,(berry_curv_adpt_kmesh-1)/2
! ikpt=ikpt+1
! adkpt(1,ikpt)=i*db1/berry_curv_adpt_kmesh
! adkpt(2,ikpt)=j*db2/berry_curv_adpt_kmesh
! adkpt(3,ikpt)=k*db3/berry_curv_adpt_kmesh
! end do
! end do
! end do
!
! NEW VERSION (both even and odd grids)
!
do i = 0, berry_curv_adpt_kmesh - 1
do j = 0, berry_curv_adpt_kmesh - 1
do k = 0, berry_curv_adpt_kmesh - 1
ikpt = ikpt + 1
adkpt(1, ikpt) = db1*((i + 0.5_dp)/berry_curv_adpt_kmesh - 0.5_dp)
adkpt(2, ikpt) = db2*((j + 0.5_dp)/berry_curv_adpt_kmesh - 0.5_dp)
adkpt(3, ikpt) = db3*((k + 0.5_dp)/berry_curv_adpt_kmesh - 0.5_dp)
end do
end do
end do
! Loop over interpolation k-points
!
if (wanint_kpoint_file) then
! NOTE: still need to specify berry_kmesh in the input file
!
! - Must use the correct nominal value in order to
! correctly set up adaptive smearing in kubo
if (on_root) write (stdout, '(/,1x,a,i10,a)') &
'Reading interpolation grid from file kpoint.dat: ', &
sum(num_int_kpts_on_node), ' points'
! Loop over k-points on the irreducible wedge of the Brillouin
! zone, read from file 'kpoint.dat'
!
do loop_xyz = 1, num_int_kpts_on_node(my_node_id)
kpt(:) = int_kpts(:, loop_xyz)
kweight = weight(loop_xyz)
kweight_adpt = kweight/berry_curv_adpt_kmesh**3
! .
! ***BEGIN COPY OF CODE BLOCK 1***
!
if (eval_ahc) then
call berry_get_imf_klist(kpt, imf_k_list)
ladpt = .false.
do if = 1, nfermi
vdum(1) = sum(imf_k_list(:, 1, if))
vdum(2) = sum(imf_k_list(:, 2, if))
vdum(3) = sum(imf_k_list(:, 3, if))
if (berry_curv_unit == 'bohr2') vdum = vdum/bohr**2
rdum = sqrt(dot_product(vdum, vdum))
if (rdum > berry_curv_adpt_kmesh_thresh) then
adpt_counter_list(if) = adpt_counter_list(if) + 1
ladpt(if) = .true.
else
imf_list(:, :, if) = imf_list(:, :, if) + imf_k_list(:, :, if)*kweight
endif
enddo
if (any(ladpt)) then
do loop_adpt = 1, berry_curv_adpt_kmesh**3
! Using imf_k_list here would corrupt values for other
! frequencies, hence dummy. Only if-th element is used
call berry_get_imf_klist(kpt(:) + adkpt(:, loop_adpt), &
imf_k_list_dummy, ladpt=ladpt)
do if = 1, nfermi
if (ladpt(if)) then
imf_list(:, :, if) = imf_list(:, :, if) &
+ imf_k_list_dummy(:, :, if)*kweight_adpt
endif
enddo
end do
endif
end if
if (eval_morb) then
call berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list)
imf_list2 = imf_list2 + imf_k_list*kweight
img_list = img_list + img_k_list*kweight
imh_list = imh_list + imh_k_List*kweight
endif
if (eval_kubo) then
if (spin_decomp) then
call berry_get_kubo_k(kpt, kubo_H_k, kubo_AH_k, jdos_k, &
kubo_H_k_spn, kubo_AH_k_spn, jdos_k_spn)
else
call berry_get_kubo_k(kpt, kubo_H_k, kubo_AH_k, jdos_k)
endif
kubo_H = kubo_H + kubo_H_k*kweight
kubo_AH = kubo_AH + kubo_AH_k*kweight
jdos = jdos + jdos_k*kweight
if (spin_decomp) then
kubo_H_spn = kubo_H_spn + kubo_H_k_spn*kweight
kubo_AH_spn = kubo_AH_spn + kubo_AH_k_spn*kweight
jdos_spn = jdos_spn + jdos_k_spn*kweight
endif
endif
if (eval_sc) then
call berry_get_sc_klist(kpt, sc_k_list)
sc_list = sc_list + sc_k_list*kweight
end if
!
! ***END COPY OF CODE BLOCK 1***
if (eval_shc) then
! print calculation progress, from 0%, 10%, ... to 100%
! Note the 1st call to berry_get_shc_klist will be much longer
! than later calls due to the time spent on
! berry_get_shc_klist -> wham_get_eig_deleig ->
! pw90common_fourier_R_to_k -> ws_translate_dist
call berry_print_progress(loop_xyz, 1, num_int_kpts_on_node(my_node_id), 1)
if (.not. shc_freq_scan) then
call berry_get_shc_klist(kpt, shc_k_fermi=shc_k_fermi)
!check whether needs to tigger adpt kmesh or not.
!Since the calculated shc_k at one Fermi energy can be reused
!by all the Fermi energies, if we find out that at a specific
!Fermi energy shc_k(if) > thresh, then we will update shc_k at
!all the Fermi energies as well.
!This also avoids repeated calculation if shc_k(if) > thresh
!is satisfied at more than one Fermi energy.
ladpt_kmesh = .false.
!if adpt_kmesh==1, no need to calculate on the same kpt again.
!This happens if adpt_kmesh==1 while adpt_kmesh_thresh is low.
if (berry_curv_adpt_kmesh > 1) then
do if = 1, nfermi
rdum = abs(shc_k_fermi(if))
if (berry_curv_unit == 'bohr2') rdum = rdum/bohr**2
if (rdum > berry_curv_adpt_kmesh_thresh) then
adpt_counter_list(1) = adpt_counter_list(1) + 1
ladpt_kmesh = .true.
exit
endif
enddo
else
ladpt_kmesh = .false.
end if
if (ladpt_kmesh) then
do loop_adpt = 1, berry_curv_adpt_kmesh**3
!Using shc_k here would corrupt values for other
!kpt, hence dummy. Only if-th element is used.
call berry_get_shc_klist(kpt(:) + adkpt(:, loop_adpt), &
shc_k_fermi=shc_k_fermi_dummy)
shc_fermi = shc_fermi + kweight_adpt*shc_k_fermi_dummy
end do
else
shc_fermi = shc_fermi + kweight*shc_k_fermi
end if
else ! freq_scan, no adaptive kmesh
call berry_get_shc_klist(kpt, shc_k_freq=shc_k_freq)
shc_freq = shc_freq + kweight*shc_k_freq
end if
end if
end do !loop_xyz
else ! Do not read 'kpoint.dat'. Loop over a regular grid in the full BZ
kweight = db1*db2*db3
kweight_adpt = kweight/berry_curv_adpt_kmesh**3
do loop_xyz = my_node_id, PRODUCT(berry_kmesh) - 1, num_nodes
loop_x = loop_xyz/(berry_kmesh(2)*berry_kmesh(3))
loop_y = (loop_xyz - loop_x*(berry_kmesh(2) &
*berry_kmesh(3)))/berry_kmesh(3)
loop_z = loop_xyz - loop_x*(berry_kmesh(2)*berry_kmesh(3)) &
- loop_y*berry_kmesh(3)
kpt(1) = loop_x*db1
kpt(2) = loop_y*db2
kpt(3) = loop_z*db3
! ***BEGIN CODE BLOCK 1***
!
if (eval_ahc) then
call berry_get_imf_klist(kpt, imf_k_list)
ladpt = .false.
do if = 1, nfermi
vdum(1) = sum(imf_k_list(:, 1, if))
vdum(2) = sum(imf_k_list(:, 2, if))
vdum(3) = sum(imf_k_list(:, 3, if))
if (berry_curv_unit == 'bohr2') vdum = vdum/bohr**2
rdum = sqrt(dot_product(vdum, vdum))
if (rdum > berry_curv_adpt_kmesh_thresh) then
adpt_counter_list(if) = adpt_counter_list(if) + 1
ladpt(if) = .true.
else
imf_list(:, :, if) = imf_list(:, :, if) + imf_k_list(:, :, if)*kweight
endif
enddo
if (any(ladpt)) then
do loop_adpt = 1, berry_curv_adpt_kmesh**3
! Using imf_k_list here would corrupt values for other
! frequencies, hence dummy. Only if-th element is used
call berry_get_imf_klist(kpt(:) + adkpt(:, loop_adpt), &
imf_k_list_dummy, ladpt=ladpt)
do if = 1, nfermi
if (ladpt(if)) then
imf_list(:, :, if) = imf_list(:, :, if) &
+ imf_k_list_dummy(:, :, if)*kweight_adpt
endif
enddo
end do
endif
end if
if (eval_morb) then
call berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list)
imf_list2 = imf_list2 + imf_k_list*kweight
img_list = img_list + img_k_list*kweight
imh_list = imh_list + imh_k_List*kweight
endif
if (eval_kubo) then
if (spin_decomp) then
call berry_get_kubo_k(kpt, kubo_H_k, kubo_AH_k, jdos_k, &
kubo_H_k_spn, kubo_AH_k_spn, jdos_k_spn)
else
call berry_get_kubo_k(kpt, kubo_H_k, kubo_AH_k, jdos_k)
endif
kubo_H = kubo_H + kubo_H_k*kweight
kubo_AH = kubo_AH + kubo_AH_k*kweight
jdos = jdos + jdos_k*kweight
if (spin_decomp) then
kubo_H_spn = kubo_H_spn + kubo_H_k_spn*kweight
kubo_AH_spn = kubo_AH_spn + kubo_AH_k_spn*kweight
jdos_spn = jdos_spn + jdos_k_spn*kweight
endif
endif
if (eval_sc) then
call berry_get_sc_klist(kpt, sc_k_list)
sc_list = sc_list + sc_k_list*kweight
end if
!
! ***END CODE BLOCK 1***
if (eval_shc) then
! print calculation progress, from 0%, 10%, ... to 100%
! Note the 1st call to berry_get_shc_klist will be much longer
! than later calls due to the time spent on
! berry_get_shc_klist -> wham_get_eig_deleig ->
! pw90common_fourier_R_to_k -> ws_translate_dist
call berry_print_progress(loop_xyz, my_node_id, PRODUCT(berry_kmesh) - 1, num_nodes)
if (.not. shc_freq_scan) then
call berry_get_shc_klist(kpt, shc_k_fermi=shc_k_fermi)
!check whether needs to tigger adpt kmesh or not.
!Since the calculated shc_k at one Fermi energy can be reused
!by all the Fermi energies, if we find out that at a specific
!Fermi energy shc_k(if) > thresh, then we will update shc_k at
!all the Fermi energies as well.
!This also avoids repeated calculation if shc_k(if) > thresh
!is satisfied at more than one Fermi energy.
ladpt_kmesh = .false.
!if adpt_kmesh==1, no need to calculate on the same kpt again.
!This happens if adpt_kmesh==1 while adpt_kmesh_thresh is low.
if (berry_curv_adpt_kmesh > 1) then
do if = 1, nfermi
rdum = abs(shc_k_fermi(if))
if (berry_curv_unit == 'bohr2') rdum = rdum/bohr**2
if (rdum > berry_curv_adpt_kmesh_thresh) then
adpt_counter_list(1) = adpt_counter_list(1) + 1
ladpt_kmesh = .true.
exit
endif
enddo
else
ladpt_kmesh = .false.
end if
if (ladpt_kmesh) then
do loop_adpt = 1, berry_curv_adpt_kmesh**3
!Using shc_k here would corrupt values for other
!kpt, hence dummy. Only if-th element is used.
call berry_get_shc_klist(kpt(:) + adkpt(:, loop_adpt), &
shc_k_fermi=shc_k_fermi_dummy)
shc_fermi = shc_fermi + kweight_adpt*shc_k_fermi_dummy
end do
else
shc_fermi = shc_fermi + kweight*shc_k_fermi
end if
else ! freq_scan, no adaptive kmesh
call berry_get_shc_klist(kpt, shc_k_freq=shc_k_freq)
shc_freq = shc_freq + kweight*shc_k_freq
end if
end if
end do !loop_xyz
end if !wanint_kpoint_file
! Collect contributions from all nodes
!
if (eval_ahc) then
call comms_reduce(imf_list(1, 1, 1), 3*3*nfermi, 'SUM')
call comms_reduce(adpt_counter_list(1), nfermi, 'SUM')
endif
if (eval_morb) then
call comms_reduce(imf_list2(1, 1, 1), 3*3*nfermi, 'SUM')
call comms_reduce(img_list(1, 1, 1), 3*3*nfermi, 'SUM')
call comms_reduce(imh_list(1, 1, 1), 3*3*nfermi, 'SUM')
end if
if (eval_kubo) then
call comms_reduce(kubo_H(1, 1, 1), 3*3*kubo_nfreq, 'SUM')
call comms_reduce(kubo_AH(1, 1, 1), 3*3*kubo_nfreq, 'SUM')
call comms_reduce(jdos(1), kubo_nfreq, 'SUM')
if (spin_decomp) then
call comms_reduce(kubo_H_spn(1, 1, 1, 1), 3*3*3*kubo_nfreq, 'SUM')
call comms_reduce(kubo_AH_spn(1, 1, 1, 1), 3*3*3*kubo_nfreq, 'SUM')
call comms_reduce(jdos_spn(1, 1), 3*kubo_nfreq, 'SUM')
endif
endif
if (eval_sc) then
call comms_reduce(sc_list(1, 1, 1), 3*6*kubo_nfreq, 'SUM')
end if
if (eval_shc) then
if (shc_freq_scan) then
call comms_reduce(shc_freq(1), kubo_nfreq, 'SUM')
else
call comms_reduce(shc_fermi(1), nfermi, 'SUM')
call comms_reduce(adpt_counter_list(1), nfermi, 'SUM')
end if
end if
if (on_root) then
if (timing_level > 1) call io_stopwatch('berry: k-interpolation', 2)
write (stdout, '(1x,a)') ' '
if (eval_ahc .and. berry_curv_adpt_kmesh .ne. 1) then
if (.not. wanint_kpoint_file) write (stdout, '(1x,a28,3(i0,1x))') &
'Regular interpolation grid: ', berry_kmesh
write (stdout, '(1x,a28,3(i0,1x))') 'Adaptive refinement grid: ', &
berry_curv_adpt_kmesh, berry_curv_adpt_kmesh, berry_curv_adpt_kmesh
if (berry_curv_unit == 'ang2') then
write (stdout, '(1x,a28,a17,f6.2,a)') &
'Refinement threshold: ', 'Berry curvature >', &
berry_curv_adpt_kmesh_thresh, ' Ang^2'
elseif (berry_curv_unit == 'bohr2') then
write (stdout, '(1x,a28,a17,f6.2,a)') &
'Refinement threshold: ', 'Berry curvature >', &
berry_curv_adpt_kmesh_thresh, ' bohr^2'
endif
if (nfermi == 1) then
if (wanint_kpoint_file) then
write (stdout, '(1x,a30,i5,a,f5.2,a)') &
' Points triggering refinement: ', &
adpt_counter_list(1), '(', &
100*real(adpt_counter_list(1), dp) &
/sum(num_int_kpts_on_node), '%)'
else
write (stdout, '(1x,a30,i5,a,f5.2,a)') &
' Points triggering refinement: ', &
adpt_counter_list(1), '(', &
100*real(adpt_counter_list(1), dp)/product(berry_kmesh), '%)'
endif
endif
elseif (eval_shc) then
if (berry_curv_adpt_kmesh .ne. 1) then
if (.not. wanint_kpoint_file) write (stdout, '(1x,a28,3(i0,1x))') &
'Regular interpolation grid: ', berry_kmesh
if (.not. shc_freq_scan) then
write (stdout, '(1x,a28,3(i0,1x))') &
'Adaptive refinement grid: ', &
berry_curv_adpt_kmesh, berry_curv_adpt_kmesh, berry_curv_adpt_kmesh
if (berry_curv_unit == 'ang2') then
write (stdout, '(1x,a28,f12.2,a)') &
'Refinement threshold: ', &
berry_curv_adpt_kmesh_thresh, ' Ang^2'
elseif (berry_curv_unit == 'bohr2') then
write (stdout, '(1x,a28,f12.2,a)') &
'Refinement threshold: ', &
berry_curv_adpt_kmesh_thresh, ' bohr^2'
endif
if (wanint_kpoint_file) then
write (stdout, '(1x,a30,i8,a,f6.2,a)') &
' Points triggering refinement: ', adpt_counter_list(1), '(', &
100*real(adpt_counter_list(1), dp)/sum(num_int_kpts_on_node), '%)'
else
write (stdout, '(1x,a30,i8,a,f6.2,a)') &
' Points triggering refinement: ', adpt_counter_list(1), '(', &
100*real(adpt_counter_list(1), dp)/product(berry_kmesh), '%)'
endif
endif
else
if (.not. wanint_kpoint_file) write (stdout, '(1x,a20,3(i0,1x))') &
'Interpolation grid: ', berry_kmesh(1:3)
endif
write (stdout, '(a)') ''
if (kubo_adpt_smr) then
write (stdout, '(1x,a)') 'Using adaptive smearing'
write (stdout, '(7x,a,f8.3)') 'adaptive smearing prefactor ', kubo_adpt_smr_fac
write (stdout, '(7x,a,f8.3,a)') 'adaptive smearing max width ', kubo_adpt_smr_max, ' eV'
else
write (stdout, '(1x,a)') 'Using fixed smearing'
write (stdout, '(7x,a,f8.3,a)') 'fixed smearing width ', &
kubo_smr_fixed_en_width, ' eV'
endif
write (stdout, '(a)') ''
if (abs(scissors_shift) > 1.0e-7_dp) then
write (stdout, '(1X,A,I0,A,G18.10,A)') "Using scissors_shift to shift energy bands with index > ", &
num_valence_bands, " by ", scissors_shift, " eV."
endif
if (shc_bandshift) then
write (stdout, '(1X,A,I0,A,G18.10,A)') "Using shc_bandshift to shift energy bands with index >= ", &
shc_bandshift_firstband, " by ", shc_bandshift_energyshift, " eV."
endif
else
if (.not. wanint_kpoint_file) write (stdout, '(1x,a20,3(i0,1x))') &
'Interpolation grid: ', berry_kmesh(1:3)
endif
if (eval_ahc) then
!
! --------------------------------------------------------------------
! At this point imf contains
!
! (1/N) sum_k Omega_{alpha beta}(k),
!
! an approximation to
!
! V_c.int dk/(2.pi)^3 Omega_{alpha beta}(k) dk
!
! (V_c is the cell volume). We want
!
! sigma_{alpha beta}=-(e^2/hbar) int dk/(2.pi)^3 Omega(k) dk
!
! Hence need to multiply by -(e^2/hbar.V_c).
! To get a conductivity in units of S/cm,
!
! (i) Divide by V_c to obtain (1/N) sum_k omega(k)/V_c, with units
! of [L]^{-1} (Berry curvature Omega(k) has units of [L]^2)
! (ii) [L] = Angstrom. Multiply by 10^8 to convert to (cm)^{-1}
! (iii) Multiply by -e^2/hbar in SI, with has units ofconductance,
! (Ohm)^{-1}, or Siemens (S), to get the final result in S/cm
!
! ===========================
! fac = -e^2/(hbar.V_c*10^-8)
! ===========================
!
! with 'V_c' in Angstroms^3, and 'e', 'hbar' in SI units
! --------------------------------------------------------------------
!
fac = -1.0e8_dp*elem_charge_SI**2/(hbar_SI*cell_volume)
ahc_list(:, :, :) = imf_list(:, :, :)*fac
if (nfermi > 1) then
write (stdout, '(/,1x,a)') &
'---------------------------------'
write (stdout, '(1x,a)') &
'Output data files related to AHC:'
write (stdout, '(1x,a)') &
'---------------------------------'
file_name = trim(seedname)//'-ahc-fermiscan.dat'
write (stdout, '(/,3x,a)') '* '//file_name
file_unit = io_file_unit()
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
endif
do if = 1, nfermi
if (nfermi > 1) write (file_unit, '(4(F12.6,1x))') &
fermi_energy_list(if), sum(ahc_list(:, 1, if)), &
sum(ahc_list(:, 2, if)), sum(ahc_list(:, 3, if))
write (stdout, '(/,1x,a18,F10.4)') 'Fermi energy (ev):', &
fermi_energy_list(if)
if (nfermi > 1) then
if (wanint_kpoint_file) then
write (stdout, '(1x,a30,i5,a,f5.2,a)') &
' Points triggering refinement: ', &
adpt_counter_list(if), '(', &
100*real(adpt_counter_list(if), dp) &
/sum(num_int_kpts_on_node), '%)'
else
write (stdout, '(1x,a30,i5,a,f5.2,a)') &
' Points triggering refinement: ', &
adpt_counter_list(if), '(', &
100*real(adpt_counter_list(if), dp) &
/product(berry_kmesh), '%)'
endif
endif
write (stdout, '(/,1x,a)') &
'AHC (S/cm) x y z'
if (iprint > 1) then
write (stdout, '(1x,a)') &
'=========='
write (stdout, '(1x,a9,2x,3(f10.4,1x))') 'J0 term :', &
ahc_list(1, 1, if), ahc_list(1, 2, if), ahc_list(1, 3, if)
write (stdout, '(1x,a9,2x,3(f10.4,1x))') 'J1 term :', &
ahc_list(2, 1, if), ahc_list(2, 2, if), ahc_list(2, 3, if)
write (stdout, '(1x,a9,2x,3(f10.4,1x))') 'J2 term :', &
ahc_list(3, 1, if), ahc_list(3, 2, if), ahc_list(3, 3, if)
write (stdout, '(1x,a)') &
'-------------------------------------------'
write (stdout, '(1x,a9,2x,3(f10.4,1x),/)') 'Total :', &
sum(ahc_list(:, 1, if)), sum(ahc_list(:, 2, if)), &
sum(ahc_list(:, 3, if))
else
write (stdout, '(1x,a10,1x,3(f10.4,1x),/)') '==========', &
sum(ahc_list(:, 1, if)), sum(ahc_list(:, 2, if)), &
sum(ahc_list(:, 3, if))
endif
enddo
if (nfermi > 1) close (file_unit)
endif
if (eval_morb) then
!
! --------------------------------------------------------------------
! At this point X=img_ab(:)-fermi_energy*imf_ab(:) and
! Y=imh_ab(:)-fermi_energy*imf_ab(:)
! contain, eg,
!
! (1/N) sum_k X(k), where X(k)=-2*Im[g(k)-E_F.f(k)]
!
! This is an approximation to
!
! V_c.int dk/(2.pi)^3 X(k) dk
!
! (V_c is the cell volume). We want a magnetic moment per cell,
! in units of the Bohr magneton. The magnetization-like quantity is
!
! \tilde{M}^LC=-(e/2.hbar) int dk/(2.pi)^3 X(k) dk
!
! So we take X and
!
! (i) The summand is an energy in eV times a Berry curvature in
! Ang^2. To convert to a.u., divide by 27.2 and by 0.529^2
! (ii) Multiply by -(e/2.hbar)=-1/2 in atomic units
! (iii) At this point we have a magnetic moment (per cell) in atomic
! units. 1 Bohr magneton = 1/2 atomic unit, so need to multiply
! by 2 to convert it to Bohr magnetons
! --------------------------------------------------------------------
!
fac = -eV_au/bohr**2
if (nfermi > 1) then
write (stdout, '(/,1x,a)') &
'---------------------------------'
write (stdout, '(1x,a)') &
'Output data files related to the orbital magnetization:'
write (stdout, '(1x,a)') &
'---------------------------------'
file_name = trim(seedname)//'-morb-fermiscan.dat'
write (stdout, '(/,3x,a)') '* '//file_name
file_unit = io_file_unit()
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
endif
do if = 1, nfermi
LCtil_list(:, :, if) = (img_list(:, :, if) &
- fermi_energy_list(if)*imf_list2(:, :, if))*fac
ICtil_list(:, :, if) = (imh_list(:, :, if) &
- fermi_energy_list(if)*imf_list2(:, :, if))*fac
Morb_list(:, :, if) = LCtil_list(:, :, if) + ICtil_list(:, :, if)
if (nfermi > 1) write (file_unit, '(4(F12.6,1x))') &
fermi_energy_list(if), sum(Morb_list(1:3, 1, if)), &
sum(Morb_list(1:3, 2, if)), sum(Morb_list(1:3, 3, if))
write (stdout, '(/,/,1x,a,F12.6)') 'Fermi energy (ev) =', &
fermi_energy_list(if)
write (stdout, '(/,/,1x,a)') &
'M_orb (bohr magn/cell) x y z'
if (iprint > 1) then
write (stdout, '(1x,a)') &
'======================'
write (stdout, '(1x,a22,2x,3(f10.4,1x))') 'Local circulation :', &
sum(LCtil_list(1:3, 1, if)), sum(LCtil_list(1:3, 2, if)), &
sum(LCtil_list(1:3, 3, if))
write (stdout, '(1x,a22,2x,3(f10.4,1x))') &
'Itinerant circulation:', &
sum(ICtil_list(1:3, 1, if)), sum(ICtil_list(1:3, 2, if)), &
sum(ICtil_list(1:3, 3, if))
write (stdout, '(1x,a)') &
'--------------------------------------------------------'
write (stdout, '(1x,a22,2x,3(f10.4,1x),/)') 'Total :', &
sum(Morb_list(1:3, 1, if)), sum(Morb_list(1:3, 2, if)), &
sum(Morb_list(1:3, 3, if))
else
write (stdout, '(1x,a22,2x,3(f10.4,1x),/)') &
'======================', &
sum(Morb_list(1:3, 1, if)), sum(Morb_list(1:3, 2, if)), &
sum(Morb_list(1:3, 3, if))
endif
enddo
if (nfermi > 1) close (file_unit)
endif
! -----------------------------!
! Complex optical conductivity !
! -----------------------------!
!
if (eval_kubo) then
!
! Convert to S/cm
fac = 1.0e8_dp*elem_charge_SI**2/(hbar_SI*cell_volume)
kubo_H = kubo_H*fac
kubo_AH = kubo_AH*fac
if (spin_decomp) then
kubo_H_spn = kubo_H_spn*fac
kubo_AH_spn = kubo_AH_spn*fac
endif
!
write (stdout, '(/,1x,a)') &
'----------------------------------------------------------'
write (stdout, '(1x,a)') &
'Output data files related to complex optical conductivity:'
write (stdout, '(1x,a)') &
'----------------------------------------------------------'
!
! Symmetric: real (imaginary) part is Hermitean (anti-Hermitean)
!
do n = 1, 6
i = alpha_S(n)
j = beta_S(n)
file_name = trim(seedname)//'-kubo_S_'// &
achar(119 + i)//achar(119 + j)//'.dat'
file_name = trim(file_name)
file_unit = io_file_unit()
write (stdout, '(/,3x,a)') '* '//file_name
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
do ifreq = 1, kubo_nfreq
if (spin_decomp) then
write (file_unit, '(9E16.8)') real(kubo_freq_list(ifreq), dp), &
real(0.5_dp*(kubo_H(i, j, ifreq) + kubo_H(j, i, ifreq)), dp), &
aimag(0.5_dp*(kubo_AH(i, j, ifreq) + kubo_AH(j, i, ifreq))), &
real(0.5_dp*(kubo_H_spn(i, j, 1, ifreq) &
+ kubo_H_spn(j, i, 1, ifreq)), dp), &
aimag(0.5_dp*(kubo_AH_spn(i, j, 1, ifreq) &
+ kubo_AH_spn(j, i, 1, ifreq))), &
real(0.5_dp*(kubo_H_spn(i, j, 2, ifreq) &
+ kubo_H_spn(j, i, 2, ifreq)), dp), &
aimag(0.5_dp*(kubo_AH_spn(i, j, 2, ifreq) &
+ kubo_AH_spn(j, i, 2, ifreq))), &
real(0.5_dp*(kubo_H_spn(i, j, 3, ifreq) &
+ kubo_H_spn(j, i, 3, ifreq)), dp), &
aimag(0.5_dp*(kubo_AH_spn(i, j, 3, ifreq) &
+ kubo_AH_spn(j, i, 3, ifreq)))
else
write (file_unit, '(3E16.8)') real(kubo_freq_list(ifreq), dp), &
real(0.5_dp*(kubo_H(i, j, ifreq) + kubo_H(j, i, ifreq)), dp), &
aimag(0.5_dp*(kubo_AH(i, j, ifreq) + kubo_AH(j, i, ifreq)))
endif
enddo
close (file_unit)
enddo
!
! Antisymmetric: real (imaginary) part is anti-Hermitean (Hermitean)
!
do n = 1, 3
i = alpha_A(n)
j = beta_A(n)
file_name = trim(seedname)//'-kubo_A_'// &
achar(119 + i)//achar(119 + j)//'.dat'
file_name = trim(file_name)
file_unit = io_file_unit()
write (stdout, '(/,3x,a)') '* '//file_name
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
do ifreq = 1, kubo_nfreq
if (spin_decomp) then
write (file_unit, '(9E16.8)') real(kubo_freq_list(ifreq), dp), &
real(0.5_dp*(kubo_AH(i, j, ifreq) - kubo_AH(j, i, ifreq)), dp), &
aimag(0.5_dp*(kubo_H(i, j, ifreq) - kubo_H(j, i, ifreq))), &
real(0.5_dp*(kubo_AH_spn(i, j, 1, ifreq) &
- kubo_AH_spn(j, i, 1, ifreq)), dp), &
aimag(0.5_dp*(kubo_H_spn(i, j, 1, ifreq) &
- kubo_H_spn(j, i, 1, ifreq))), &
real(0.5_dp*(kubo_AH_spn(i, j, 2, ifreq) &
- kubo_AH_spn(j, i, 2, ifreq)), dp), &
aimag(0.5_dp*(kubo_H_spn(i, j, 2, ifreq) &
- kubo_H_spn(j, i, 2, ifreq))), &
real(0.5_dp*(kubo_AH_spn(i, j, 3, ifreq) &
- kubo_AH_spn(j, i, 3, ifreq)), dp), &
aimag(0.5_dp*(kubo_H_spn(i, j, 3, ifreq) &
- kubo_H_spn(j, i, 3, ifreq)))
else
write (file_unit, '(3E16.8)') real(kubo_freq_list(ifreq), dp), &
real(0.5_dp*(kubo_AH(i, j, ifreq) - kubo_AH(j, i, ifreq)), dp), &
aimag(0.5_dp*(kubo_H(i, j, ifreq) - kubo_H(j, i, ifreq)))
endif
enddo
close (file_unit)
enddo
!
! Joint density of states
!
file_name = trim(seedname)//'-jdos.dat'
write (stdout, '(/,3x,a)') '* '//file_name
file_unit = io_file_unit()
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
do ifreq = 1, kubo_nfreq
if (spin_decomp) then
write (file_unit, '(5E16.8)') real(kubo_freq_list(ifreq), dp), &
jdos(ifreq), jdos_spn(:, ifreq)
else
write (file_unit, '(2E16.8)') real(kubo_freq_list(ifreq), dp), &
jdos(ifreq)
endif
enddo
close (file_unit)
endif
if (eval_sc) then
! -----------------------------!
! Nonlinear shift current
! -----------------------------!
! --------------------------------------------------------------------
! At this point sc_list contains
!
! (1/N) sum_k (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})(k) delta(w),
!
! an approximation to
!
! V_c.int dk/(2.pi)^3 (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})(k) delta(w) dk
!
! (V_c is the cell volume). We want
!
! sigma_{abc}=( pi.e^3/(4.hbar^2) ) int dk/(2.pi)^3 Im[ (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})(k) delta(w) ] dk
!
! Note factor 1/4 instead of 1/2 as compared to SS PRB 61 5337 (2000) (Eq. 57),
! because we introduce 2 delta functions instead of 1.
! Hence we need to multiply by pi.e^3/(4.hbar^2.V_c).
! To get the nonlinear response in units of A/V^2,
!
! (i) Divide by V_c to obtain (1/N) sum_k (r_^{b}r^{c}_{a}+r_^{c}r^{b}_{a})delta(w)/V_c, with units
! of [T] (integrand terms r_^{b}r^{c}_{a} delta(w) have units of [T].[L]^3)
! (ii) Multiply by eV_seconds to convert the units of [T] from eV to seconds (coming from delta function)
! (iii) Multiply by ( pi.e^3/(4.hbar^2) ) in SI, which multiplied by [T] in seconds from (ii), gives final
! units of A/V^2
!
! ===========================
! fac = eV_seconds.( pi.e^3/(4.hbar^2.V_c) )
! ===========================
!
! with 'V_c' in Angstroms^3, and 'e', 'hbar' in SI units
! --------------------------------------------------------------------
fac = eV_seconds*pi*elem_charge_SI**3/(4*hbar_SI**(2)*cell_volume)
write (stdout, '(/,1x,a)') &
'----------------------------------------------------------'
write (stdout, '(1x,a)') &
'Output data files related to shift current: '
write (stdout, '(1x,a)') &
'----------------------------------------------------------'
do i = 1, 3
do jk = 1, 6
j = alpha_S(jk)
k = beta_S(jk)
file_name = trim(seedname)//'-sc_'// &
achar(119 + i)//achar(119 + j)//achar(119 + k)//'.dat'
file_name = trim(file_name)
file_unit = io_file_unit()
write (stdout, '(/,3x,a)') '* '//file_name
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
do ifreq = 1, kubo_nfreq
write (file_unit, '(2E18.8E3)') real(kubo_freq_list(ifreq), dp), &
fac*sc_list(i, jk, ifreq)
enddo
close (file_unit)
enddo
enddo
endif
! -----------------------!
! Spin Hall conductivity !
! -----------------------!
!
if (eval_shc) then
!
! Convert to the unit: (hbar/e) S/cm
! at this point, we need to
! (i) multiply -e^2/hbar/(V*N_k) as in the QZYZ18 Eq.(5),
! note 1/N_k has already been applied by the kweight
! (ii) convert charge current to spin current:
! divide the result by -e and multiply hbar/2 to
! recover the spin current, so the overall
! effect is -hbar/2/e
! (iii) multiply 1e8 to convert it to the unit S/cm
! So, the overall factor is
! fac = 1.0e8 * e^2 / hbar / V / 2.0
! and the final unit of spin Hall conductivity is (hbar/e)S/cm
!
fac = 1.0e8_dp*elem_charge_SI**2/(hbar_SI*cell_volume)/2.0_dp
if (shc_freq_scan) then
shc_freq = shc_freq*fac
else
shc_fermi = shc_fermi*fac
endif
!
write (stdout, '(/,1x,a)') &
'----------------------------------------------------------'
write (stdout, '(1x,a)') &
'Output data files related to Spin Hall conductivity:'
write (stdout, '(1x,a)') &
'----------------------------------------------------------'
!
if (.not. shc_freq_scan) then
file_name = trim(seedname)//'-shc-fermiscan'//'.dat'
else
file_name = trim(seedname)//'-shc-freqscan'//'.dat'
endif
file_name = trim(file_name)
file_unit = io_file_unit()
write (stdout, '(/,3x,a)') '* '//file_name
open (file_unit, FILE=file_name, STATUS='UNKNOWN', FORM='FORMATTED')
if (.not. shc_freq_scan) then
write (file_unit, '(a,3x,a,3x,a)') &
'#No.', 'Fermi energy(eV)', 'SHC((hbar/e)*S/cm)'
do n = 1, nfermi
write (file_unit, '(I4,1x,F12.6,1x,E17.8)') &
n, fermi_energy_list(n), shc_fermi(n)
enddo
else
write (file_unit, '(a,3x,a,3x,a,3x,a)') '#No.', 'Frequency(eV)', &
'Re(sigma)((hbar/e)*S/cm)', 'Im(sigma)((hbar/e)*S/cm)'
do n = 1, kubo_nfreq
write (file_unit, '(I4,1x,F12.6,1x,1x,2(E17.8,1x))') n, &
real(kubo_freq_list(n), dp), real(shc_freq(n), dp), aimag(shc_freq(n))
enddo
endif
close (file_unit)
endif
end if !on_root
end subroutine berry_main