Skip to content

Commit

Permalink
Implemented .umklapp. to add fractional coordinate wave vectors.
Browse files Browse the repository at this point in the history
  • Loading branch information
nakib committed Nov 30, 2023
1 parent 3541022 commit 0311513
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 10 deletions.
4 changes: 2 additions & 2 deletions src/eliashberg.f90
Original file line number Diff line number Diff line change
Expand Up @@ -203,7 +203,7 @@ subroutine calculate_a2F(wann, el, ph, num, omegas, iso_lambda0, omegalog, exter
!Find interacting phonon wave vector
! Note that q, k, and k' are all on the same mesh
q_indvec = modulo(kp_indvec - k_indvec, el%wvmesh)

! Muxed index of q
iq = mux_vector(q_indvec, el%wvmesh, 0_i64)

Expand All @@ -222,7 +222,7 @@ subroutine calculate_a2F(wann, el, ph, num, omegas, iso_lambda0, omegalog, exter

!Note that the phonon branch index iterates last for a2F_istate
a2F_istate(count, :) = g2_istate(count)*ph_deltas(s, iq, :)

!Sum contribuion to the isotropic a2F
iso_a2F_branches(:, s) = iso_a2F_branches(:, s) + &
a2F_istate(count, :)*WWp
Expand Down
42 changes: 39 additions & 3 deletions src/misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,9 @@ module misc

implicit none

public
public :: operator(.umklapp.)
private :: sort_int, sort_real, Pade_coeffs, twonorm_real_rank1, twonorm_real_rank2, &
invert_complex_square
invert_complex_square, add_and_fold, add_and_fold_array

type timer
!! Container for timing related data and procedures.
Expand Down Expand Up @@ -57,7 +57,11 @@ module misc
interface create_set
module procedure :: create_set_int, create_set_char
end interface create_set


interface operator(.umklapp.)
module procedure add_and_fold
module procedure add_and_fold_array
end interface operator(.umklapp.)
contains

subroutine start_timer(self, event)
Expand Down Expand Up @@ -576,6 +580,38 @@ pure complex(r64) function expi(x)
expi = cmplx(cos(x), sin(x), r64)
end function expi

pure function add_and_fold(q1, q2) result(q3)
!! Adds two fractional coordinate wave vectors
!! and folds back into the 1st Brillouin zone.

real(r64), intent(in) :: q1(3), q2(3)
real(r64) :: q3(3)

integer :: icart

q3 = q1 + q2
do icart = 1, 3
if(q3(icart) >= 1.0_r64) q3(icart) = q3(icart) - 1.0_r64
if(q3(icart) < 0.0_r64) q3(icart) = q3(icart) + 1.0_r64
end do
end function add_and_fold

pure function add_and_fold_array(q1, q2) result(q3)
!! Elementwise adds two fractional coordinate wave vector arrays
!! and folds back into the 1st Brillouin zone.

real(r64), intent(in) :: q1(:, :), q2(:, :)
real(r64) :: q3(3, size(q1, 2))

integer :: iq, nq

nq = size(q3, 2)

do iq = 1, nq
q3(:, iq) = add_and_fold(q1(:, iq), q2(:, iq))
end do
end function add_and_fold_array

pure real(r64) function twonorm_real_rank1(v)
!! 2-norm of a rank-1 real vector

Expand Down
22 changes: 20 additions & 2 deletions src/wannier.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1221,6 +1221,9 @@ subroutine plot_along_path(self, crys, num, scissor)
character(len = 1024) :: filename
character(len=8) :: saux

!DBG
real(r64), allocatable :: el_ens_kp_all(:,:)

call print_message("Plotting bands, dispersions, and e-ph vertex along path...")

call self%reshape_gwann_for_gkRp
Expand Down Expand Up @@ -1270,6 +1273,9 @@ subroutine plot_along_path(self, crys, num, scissor)
el_evecs_kp(1, self%numwannbands, self%numwannbands))
allocate(g2_qpath(nqpath, self%numbranches, self%numwannbands, self%numwannbands))

!DBG
allocate(el_ens_kp_all(nqpath, self%numwannbands))

!Read wave vector of initial electron
open(1, file = trim('initialk.txt'), status = 'old')
read(1,*) k(1, :)
Expand All @@ -1287,7 +1293,7 @@ subroutine plot_along_path(self, crys, num, scissor)
close(1)
!Change back to working directory
call chdir(num%cwd)

call el_wann(self, crys, 1_i64, k, el_ens_k, el_vels_k, el_evecs_k, &
scissor = scissor)

Expand All @@ -1296,13 +1302,16 @@ subroutine plot_along_path(self, crys, num, scissor)
do icart = 1, 3
if(kp(1,icart) >= 1.0_r64) kp(1, icart) = kp(1, icart) - 1.0_r64
end do

!TODO Would be great to have a progress bar here.

!Calculate electrons at this final wave vector
call el_wann(self, crys, 1_i64, kp, el_ens_kp, el_vels_kp, el_evecs_kp, &
scissor = scissor)

!DBG
el_ens_kp_all(i, :) = el_ens_kp(1, :)

do n = 1, self%numwannbands
do m = 1, self%numwannbands
do s = 1, self%numbranches
Expand Down Expand Up @@ -1377,6 +1386,15 @@ subroutine plot_along_path(self, crys, num, scissor)
end do
end do

!DBG
!Output electron dispersions
write(saux,"(I0)") self%numwannbands
open(1, file="el.ens_k+qpath",status="replace")
do i = 1, nqpath
write(1,"("//trim(adjustl(saux))//"E20.10)") el_ens_kp_all(i,:)
end do
close(1)

!Print out |gk(m,n,s,qpath)|
open(1, file = 'gk_qpath',status="replace")
write(1,*) ' m n s |gk|[eV]'
Expand Down
36 changes: 33 additions & 3 deletions test/test_misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -7,20 +7,20 @@ program test_misc
twonorm, binsearch, mux_vector, demux_vector, interpolate, coarse_grained, &
unique, linspace, compsimps, mux_state, demux_state, demux_mesh, expm1, &
Fermi, Bose, Pade_continued, precompute_interpolation_corners_and_weights, &
interpolate_using_precomputed
interpolate_using_precomputed, operator(.umklapp.)

implicit none

integer :: itest
integer, parameter :: num_tests = 24
integer, parameter :: num_tests = 26
type(testify) :: test_array(num_tests), tests_all
integer(i64) :: index, quotient, remainder, int_array(5), v1(3), v2(3), &
v1_muxed, v2_muxed, ik, ik1, ik2, ik3, ib1, ib2, ib3, wvmesh(3), &
mesh_ref_array(3), nk_coarse, ninterp
integer(i64), allocatable :: index_mesh_0(:, :), index_mesh_1(:, :), &
ksint(:, :), idc(:, :), ik_interp(:)
real(r64) :: pauli1(2, 2), ipauli2(2, 2), pauli3(2, 2), &
real_array(5), result
real_array(5), result, q1(3, 4), q2(3, 4), q3(3, 4)
real(r64), allocatable :: integrand(:), domain(:), im_axis(:), real_func(:), &
widc(:, :), f_coarse(:), f_interp(:)

Expand Down Expand Up @@ -281,6 +281,36 @@ program test_misc
Pade_continued(oneI*im_axis, real_func, [0.0_r64, 0.5_r64, 1.0_r64]), &
[-1.0_r64 + 0.0_r64*oneI, -0.8_r64 + 0.4_r64*oneI, -0.5_r64 + 0.5_r64*oneI], &
tol = 1.0e-10_r64)

!.umklapp.
q1(:, 1) = [0.5_r64, 0.0_r64, 0.0_r64]
q1(:, 2) = [0.5_r64, 0.5_r64, 0.0_r64]
q1(:, 3) = [0.5_r64, 0.5_r64, 0.5_r64]
q1(:, 4) = [0.1_r64, 0.2_r64, 0.0_r64]
q2(:, 1) = [0.5_r64, 0.0_r64, 0.0_r64]
q2(:, 2) = [0.5_r64, 0.5_r64, 0.0_r64]
q2(:, 3) = [0.5_r64, 0.5_r64, 0.5_r64]
q2(:, 4) = -[0.5_r64, 0.5_r64, 0.1_r64]
q3(:, 1) = [0.0_r64, 0.0_r64, 0.0_r64]
q3(:, 2) = [0.0_r64, 0.0_r64, 0.0_r64]
q3(:, 3) = [0.0_r64, 0.0_r64, 0.0_r64]
q3(:, 4) = [0.6_r64, 0.7_r64, 0.9_r64]

itest = itest + 1
test_array(itest) = testify(".umklapp.")

call test_array(itest)%assert( &
[ &
q1(:, 1) .umklapp. q2(:, 1), &
q1(:, 2) .umklapp. q2(:, 2), &
q1(:, 3) .umklapp. q2(:, 3), &
q1(:, 4) .umklapp. q2(:, 4)], &
[q3(:, 1), q3(:, 2), q3(:, 3), q3(:, 4)], tol = 1.0e-12_r64)

!.umklapp. elemental
itest = itest + 1
test_array(itest) = testify("elemental .umklapp.")
call test_array(itest)%assert(pack(q1 .umklapp. q2, .true.), reshape(q3, [size(q3)]))

tests_all = testify(test_array)
call tests_all%report
Expand Down

0 comments on commit 0311513

Please sign in to comment.