From ee3f1679c923cfb045fee3b088cbac590e700508 Mon Sep 17 00:00:00 2001 From: visitor Date: Thu, 31 Jul 2025 18:35:36 +0200 Subject: [PATCH 1/6] updated .eigval parser, increased number of maximum bands --- src/parser_optics_xatu_dim.f90 | 43 +++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 9 deletions(-) diff --git a/src/parser_optics_xatu_dim.f90 b/src/parser_optics_xatu_dim.f90 index c115833..4238cdc 100644 --- a/src/parser_optics_xatu_dim.f90 +++ b/src/parser_optics_xatu_dim.f90 @@ -185,7 +185,7 @@ subroutine get_exciton_dim() file2open=trim(xatu_states_filepath_in) open(10,file=file2open) read(10,*) - do i=1,50 !we will never use more than 50 bands... + do i=1,5000 !we will never use more than 50 bands... read(10,*) aux1,aux1,aux1,nband_index_aux1(i) if (i.gt.1) then do j=1,i-1 @@ -202,7 +202,7 @@ subroutine get_exciton_dim() !save number of conduction bands open(10,file=file2open) read(10,*) - do i=1,50 !we will never use more than 50 bands... + do i=1,5000 !we will never use more than 50 bands... read(10,*) aux1,aux1,aux1,naux,nband_index_aux2(i) !skip similar values if (nband_ex_aux1.gt.1) then @@ -276,19 +276,44 @@ subroutine get_reciprocal_vectors() !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_exciton_data() implicit none - integer j,nkaka + integer j,nkaka, max_ex, ios integer ib,ibz,ibz_sum,jind real(8) auxr1 dimension auxr1(2*norb_ex) character(len=:), allocatable :: file2open !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !get energies - file2open=trim(xatu_eigval_filepath_in) - open(10,file=file2open) - read(10,*) - read(10,*) nkaka,(e_ex(j), j=1,norb_ex_cut) - close(10) + ! Read eigenvalues from file + file2open = trim(xatu_eigval_filepath_in) + + open(unit=10, file=file2open, status='old', action='read', iostat=ios) + if (ios /= 0) then + print *, 'Error opening eigenvalue file: ', trim(file2open) + stop + end if + + ! skipping headers + read(10,*) ! npointstotal_sq already read + read(10,*) ! naux already read + read(10,*) max_ex + + if (norb_ex_cut > max_ex) then + norb_ex_cut = max_ex + print *, "!!! Exciton cutoff is bigger than available eigvalues in .eigval !!!" + print *, " --> Setting cutoff to maximum number of available excitons: ", max_ex + end if + + ! read norb_ex_cut values, one per line + do j = 1, norb_ex_cut + read(10, *, iostat=ios) e_ex(j) + if (ios /= 0) then + print *, 'Error reading eigenvalue at line', j+1 + stop + end if + end do + + close(10) + file2open=trim(xatu_states_filepath_in) open(10,file=file2open) From 62df4e0d6da0062da7d911f15462122433b82117 Mon Sep 17 00:00:00 2001 From: visitor Date: Wed, 3 Dec 2025 16:00:21 +0100 Subject: [PATCH 2/6] z-direction for Sp --- src/ome_sp.f90 | 7 ++++++- src/sigma_first_ex.f90 | 3 ++- src/sigma_first_sp.f90 | 3 ++- src/sigma_second_ex.f90 | 13 +++++++++---- src/sigma_second_sp.f90 | 20 +++++++++++++------- 5 files changed, 32 insertions(+), 14 deletions(-) diff --git a/src/ome_sp.f90 b/src/ome_sp.f90 index 6ec0a91..2a82475 100644 --- a/src/ome_sp.f90 +++ b/src/ome_sp.f90 @@ -333,7 +333,12 @@ subroutine get_berry_eigen_fourpoint(rkx,rky,rkz,norb,vme_der, & berry_eigen1(2,nn,nnp)=berry_eigen1(2,nn,nnp)+ & complex(0.0d0,1.0d0)*conjg(hk_ev_neigh(5,ialpha,nn))*skernel(ialpha,ialphap)*aux1 berry_eigen2(2,nn,nnp)=berry_eigen2(2,nn,nnp)+ & - conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(2,ialpha,ialphap) + conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(2,ialpha,ialphap) + + !z-dir + berry_eigen1(3,nn,nnp)=0.0d0 + berry_eigen2(3,nn,nnp)=berry_eigen2(3,nn,nnp)+ & + conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(3,ialpha,ialphap) end do end do do nj=1,3 diff --git a/src/sigma_first_ex.f90 b/src/sigma_first_ex.f90 index 0ac72f6..8c18f55 100644 --- a/src/sigma_first_ex.f90 +++ b/src/sigma_first_ex.f90 @@ -136,7 +136,8 @@ subroutine print_sigma_first_ex(nw,wp,sigma_w_ex) do iw=1,nw feps=1.0d0 !use atomic units - write(50,*) wp(iw)*27.211385d0,realpart(feps*sigma_w_ex(1,1,iw)), & + write(50,*) wp(iw)*27.211385d0, & + realpart(feps*sigma_w_ex(1,1,iw)), & realpart(feps*sigma_w_ex(1,2,iw)), & realpart(feps*sigma_w_ex(1,3,iw)), & realpart(feps*sigma_w_ex(2,1,iw)), & diff --git a/src/sigma_first_sp.f90 b/src/sigma_first_sp.f90 index fc9fd19..63a271b 100644 --- a/src/sigma_first_sp.f90 +++ b/src/sigma_first_sp.f90 @@ -107,7 +107,8 @@ subroutine print_sigma_first_sp(nw,wp,sigma_w_sp) do iw=1,nw feps=1.0d0 !use atomic units - write(50,*) wp(iw)*27.211385d0,realpart(feps*sigma_w_sp(1,1,iw)), & + write(50,*) wp(iw)*27.211385d0, & + realpart(feps*sigma_w_sp(1,1,iw)), & realpart(feps*sigma_w_sp(1,2,iw)), & realpart(feps*sigma_w_sp(1,3,iw)), & realpart(feps*sigma_w_sp(2,1,iw)), & diff --git a/src/sigma_second_ex.f90 b/src/sigma_second_ex.f90 index 6bfd791..f8f0aa7 100644 --- a/src/sigma_second_ex.f90 +++ b/src/sigma_second_ex.f90 @@ -279,10 +279,15 @@ subroutine print_sigma_second_ex(nw,wp,sigma_w_ex) !d=3.28d0 !thickness in angstrongs for h-BN !feps=feps/(d/0.52917721067121d0) feps=(6.623618d-03)*(1.0d+06)*(27.211386**(-2))*(5.291772d-11)*(1.0d+09) !%go from au to (\mu A /V^2)*nm - write(90,*) wp(iw)*27.211385d0,realpart(feps*sigma_w_ex(1,1,1,iw)), & - realpart(feps*sigma_w_ex(1,2,2,iw)), & - realpart(feps*sigma_w_ex(2,1,1,iw)), & - realpart(feps*sigma_w_ex(2,2,2,iw)) + write(90,*) wp(iw)*27.211385d0,& + realpart(feps*sigma_w_ex(1,1,1,iw)), & + realpart(feps*sigma_w_ex(1,2,2,iw)), & + realpart(feps*sigma_w_ex(2,1,1,iw)), & + realpart(feps*sigma_w_ex(2,2,2,iw)), & + realpart(feps*sigma_w_ex(1,1,3,iw)), & + realpart(feps*sigma_w_ex(2,2,3,iw)), & + realpart(feps*sigma_w_ex(3,1,1,iw)), & + realpart(feps*sigma_w_ex(3,2,2,iw)) end do close(90) diff --git a/src/sigma_second_sp.f90 b/src/sigma_second_sp.f90 index f37ac1d..6b577ae 100644 --- a/src/sigma_second_sp.f90 +++ b/src/sigma_second_sp.f90 @@ -297,7 +297,7 @@ subroutine print_sigma_second_sp(nw,wp,sigma_w_sp,shift_vector_w) real*8 :: feps !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !write frequency dependent conductivity - open(90,file='shift_sp_shiftvector_'//trim(material_name)//'.dat') + open(90,file='shift_sp_lengthgauge_'//trim(material_name)//'.dat') open(100,file='shift_vector.dat') do iw=1,nw !feps=-6.623618d-03/(27.21138**2)*1.0d06 !%go from au to (\mu A /V^2)*Angstrongs @@ -305,12 +305,18 @@ subroutine print_sigma_second_sp(nw,wp,sigma_w_sp,shift_vector_w) !d=3.28d0 !thickness in angstrongs for h-BN !feps=feps/(d/0.52917721067121d0) feps=(6.623618d-03)*(1.0d+06)*(27.211386**(-2))*(5.291772d-11)*(1.0d+09) !%go from au to (\mu A /V^2)*nm - write(90,*) wp(iw)*27.211385d0,realpart(feps*sigma_w_sp(1,1,1,iw)), & - realpart(feps*sigma_w_sp(1,2,2,iw)), & - realpart(feps*sigma_w_sp(2,1,1,iw)), & - realpart(feps*sigma_w_sp(2,2,2,iw)) - write(100,*) wp(iw)*27.211385d0,shift_vector_w(1,1,iw),shift_vector_w(1,2,iw), & - shift_vector_w(2,1,iw),shift_vector_w(2,2,iw) + write(90,*) wp(iw)*27.211385d0,& + realpart(feps*sigma_w_sp(1,1,1,iw)), & + realpart(feps*sigma_w_sp(1,2,2,iw)), & + realpart(feps*sigma_w_sp(2,1,1,iw)), & + realpart(feps*sigma_w_sp(2,2,2,iw)), & + realpart(feps*sigma_w_sp(1,1,3,iw)), & + realpart(feps*sigma_w_sp(2,2,3,iw)), & + realpart(feps*sigma_w_sp(3,1,1,iw)), & + realpart(feps*sigma_w_sp(3,2,2,iw)) + write(100,*) wp(iw)*27.211385d0,& + shift_vector_w(1,1,iw),shift_vector_w(1,2,iw), & + shift_vector_w(2,1,iw),shift_vector_w(2,2,iw) end do close(90) close(95) From 7ebb16ea61f41d5bb671a902f8a1a844519d1ff4 Mon Sep 17 00:00:00 2001 From: visitor Date: Tue, 10 Feb 2026 11:42:34 +0100 Subject: [PATCH 3/6] printing z components --- src/ome_ex.f90 | 362 +++++++++++++++++++++++++++++++------------------ 1 file changed, 229 insertions(+), 133 deletions(-) diff --git a/src/ome_ex.f90 b/src/ome_ex.f90 index 111fe11..ef364e8 100644 --- a/src/ome_ex.f90 +++ b/src/ome_ex.f90 @@ -37,150 +37,167 @@ module ome_ex contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_ome_ex(iflag_norder) - implicit none - - !here +subroutine get_ome_ex(iflag_norder) + implicit none - integer iflag_norder - integer :: ibz - integer :: nn,nnp - integer :: nj + integer, intent(in) :: iflag_norder + integer :: ibz + integer :: nn, nnp + integer :: nj + + ! For k-resolved excitonic linear output + integer :: u_exk + logical :: do_write_exk + complex*16, allocatable :: vme_ex_k(:,:) ! (3, norb_ex_cut) + + ! auxiliary arrays used to evaluate ex-ome + dimension ek(npointstotal,nband_ex) + dimension xme_ex_band(npointstotal,3,nband_ex,nband_ex) ! only here! provisional + dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) + dimension berry_eigen_ex_band(npointstotal,3,nband_ex,nband_ex) + dimension gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex) + dimension shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex) + + real*8 :: ek + real*8 :: shift_vector_ex_band + complex*16 :: xme_ex_band + complex*16 :: vme_ex_band + complex*16 :: berry_eigen_ex_band + complex*16 :: gen_der_ex_band - !auxiliary arrays used to evaluate ex-ome - dimension ek(npointstotal,nband_ex) - dimension xme_ex_band(npointstotal,3,nband_ex,nband_ex) !only here! provisional - dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) - dimension berry_eigen_ex_band(npointstotal,3,nband_ex,nband_ex) - dimension gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex) - dimension shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex) - - real*8 :: ek - real*8 :: shift_vector_ex_band - complex*16 :: xme_ex_band - complex*16 vme_ex_band - complex*16 berry_eigen_ex_band - complex*16 gen_der_ex_band !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - write(*,*) '6. Entering ome_ex' - - !read SP optical matrix elements from file - write(*,*) ' Reading optical matrix elements (sp)...' - if (iflag_norder.eq.1) then - vme_ex_band=0.0d0 - ek=0.0d0 - call read_ome_sp_linear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek) + write(*,*) '6. Entering ome_ex' + + ! Decide whether to write the k-resolved file: + ! Only makes sense for linear (iflag_norder == 1) + do_write_exk = (iflag_norder .eq. 1) + + ! read SP optical matrix elements from file + write(*,*) ' Reading optical matrix elements (sp)...' + if (iflag_norder .eq. 1) then + vme_ex_band = 0.0d0 + ek = 0.0d0 + call read_ome_sp_linear(iflag_norder, npointstotal, nband_ex, vme_ex_band, ek) + end if + + if (iflag_norder .eq. 2) then + ek = 0.0d0 + vme_ex_band = 0.0d0 + berry_eigen_ex_band = 0.0d0 + gen_der_ex_band = 0.0d0 + shift_vector_ex_band = 0.0d0 + call read_ome_sp_nonlinear(iflag_norder, npointstotal, nband_ex, berry_eigen_ex_band, & + gen_der_ex_band, shift_vector_ex_band, vme_ex_band, ek) + + ! Provisional (19/07/2025): we evaluate here xme_ex_band at a given k-point + xme_ex_band = 0.0d0 + call get_ome_sp_xme_ex_band(ek, vme_ex_band, xme_ex_band) + end if + + ! allocate arrays for ex-ome + ! linear conductivity + if (iflag_norder .eq. 1 .or. iflag_norder .eq. 2) then + allocate(vme_ex(3, norb_ex_cut)) + allocate(xme_ex(3, norb_ex_cut)) + vme_ex = 0.0d0 + xme_ex = 0.0d0 + end if + + ! allocate k-resolved buffer + open file (linear only) + if (do_write_exk) then + allocate(vme_ex_k(3, norb_ex_cut)) + vme_ex_k = 0.0d0 + u_exk = 77 + call write_ome_ex_linear_kresolved_init(u_exk, material_name, norb_ex_cut) + end if + + ! second order ones + if (iflag_norder .eq. 2) then + allocate(qme_ex_inter1(3, norb_ex_cut, norb_ex_cut)) + allocate(qme_ex_inter2(3, norb_ex_cut, norb_ex_cut)) + allocate(qme_ex_inter(3, norb_ex_cut, norb_ex_cut)) + allocate(yme_ex_inter1(3, norb_ex_cut, norb_ex_cut)) + allocate(yme_ex_inter2(3, norb_ex_cut, norb_ex_cut)) + allocate(yme_ex_inter(3, norb_ex_cut, norb_ex_cut)) + allocate(xme_ex_inter(3, norb_ex_cut, norb_ex_cut)) + allocate(vme_ex_inter1(3, norb_ex_cut, norb_ex_cut)) + allocate(vme_ex_inter2(3, norb_ex_cut, norb_ex_cut)) + allocate(vme_ex_inter(3, norb_ex_cut, norb_ex_cut)) + + qme_ex_inter1 = 0.0d0 + qme_ex_inter2 = 0.0d0 + qme_ex_inter = 0.0d0 + yme_ex_inter1 = 0.0d0 + yme_ex_inter2 = 0.0d0 + yme_ex_inter = 0.0d0 + xme_ex_inter = 0.0d0 + vme_ex_inter1 = 0.0d0 + vme_ex_inter2 = 0.0d0 + vme_ex_inter = 0.0d0 + + ! allocate and get derivative of exciton envelope function with respect to k + allocate(fk_ex_der(3, norb_ex, norb_ex_cut)) + call get_fk_ex_der_k() + end if + + if (iflag_norder .eq. 1) then + write(*,*) ' Evaluating excitonic optical matrix elements for linear conductivity...' + end if + if (iflag_norder .eq. 2) then + write(*,*) ' Evaluating excitonic optical matrix elements for nonlinear conductivity...' + end if + + ! k-space integration of excitonic optical matrix elements + do ibz = 1, npointstotal + write(*,*) ' Optical matrix elements (ex): k-point', ibz, '/', npointstotal + + ! Fill V_{0N} and X_{0N} (summed over k) for linear conductivity + if (iflag_norder .eq. 1 .or. iflag_norder .eq. 2) then + call get_ome_gs_ex_sum_k(ibz, ek, xme_ex_band, vme_ex_band) end if - if (iflag_norder.eq.2) then - ek=0.0d0 - vme_ex_band=0.0d0 - berry_eigen_ex_band=0.0d0 - gen_der_ex_band=0.0d0 - shift_vector_ex_band=0.0d0 - call read_ome_sp_nonlinear(iflag_norder,npointstotal,nband_ex,berry_eigen_ex_band, & - gen_der_ex_band,shift_vector_ex_band,vme_ex_band,ek) - - - !Provisional (19/07/2025): we evaluate hete xme_ex_band at a given k-point - !Better to implement in ome_sp later - xme_ex_band=0.0d0 - call get_ome_sp_xme_ex_band(ek,vme_ex_band,xme_ex_band) - end if - - !allocate arrays for ex-ome - !linear conductivity - if (iflag_norder.eq.1 .or. iflag_norder.eq.2) then - allocate(vme_ex(3,norb_ex_cut)) - allocate(xme_ex(3,norb_ex_cut)) - vme_ex=0.0d0 - xme_ex=0.0d0 - end if - !second order ones - if (iflag_norder.eq.2) then - allocate(qme_ex_inter1(3,norb_ex_cut,norb_ex_cut)) - allocate(qme_ex_inter2(3,norb_ex_cut,norb_ex_cut)) - allocate(qme_ex_inter(3,norb_ex_cut,norb_ex_cut)) - allocate(yme_ex_inter1(3,norb_ex_cut,norb_ex_cut)) - allocate(yme_ex_inter2(3,norb_ex_cut,norb_ex_cut)) - allocate(yme_ex_inter(3,norb_ex_cut,norb_ex_cut)) - allocate(xme_ex_inter(3,norb_ex_cut,norb_ex_cut)) - allocate(vme_ex_inter1(3,norb_ex_cut,norb_ex_cut)) - allocate(vme_ex_inter2(3,norb_ex_cut,norb_ex_cut)) - allocate(vme_ex_inter(3,norb_ex_cut,norb_ex_cut)) - qme_ex_inter1=0.0d0 - qme_ex_inter2=0.0d0 - qme_ex_inter=0.0d0 - yme_ex_inter1=0.0d0 - yme_ex_inter2=0.0d0 - yme_ex_inter=0.0d0 - xme_ex_inter=0.0d0 - vme_ex_inter1=0.0d0 - vme_ex_inter2=0.0d0 - vme_ex_inter=0.0d0 - !allocate and getderivative of exciton envelope function with respect to k - allocate (fk_ex_der(3,norb_ex,norb_ex_cut)) - !get exciton envelope function derivative with respect to k - call get_fk_ex_der_k() !get derivative of exciton envelope function with respect to k - end if - - if (iflag_norder.eq.1) then - write(*,*) ' Evaluating excitonic optical matrix elements for linear conductivity...' + ! Also write k-resolved contribution for linear case: + if (do_write_exk) then + call get_ome_gs_ex_kresolved(ibz, ek, vme_ex_band, vme_ex_k) + call write_ome_ex_linear_kresolved_point(u_exk, rkxvector(ibz), rkyvector(ibz), rkzvector(ibz), & + norb_ex_cut, vme_ex_k) end if - if (iflag_norder.eq.2) then - write(*,*) ' Evaluating excitonic optical matrix elements for nonlinear conductivity...' - end if - !k-space integration of excitonic optical marix elements - do ibz=1,npointstotal - write(*,*) ' Optical matrix elements (ex): k-point',ibz,'/',npointstotal - !fill V_{0N} and X_{0N} for linear conductivity - if (iflag_norder.eq.1 .or. iflag_norder.eq.2) then - call get_ome_gs_ex_sum_k(ibz,ek,xme_ex_band,vme_ex_band) !sum over k points - end if - if (iflag_norder.eq.2) then - !fill Q_{NN'} (1,2), Y_{NN'} (1,2) and V_{NN'} (1,2) - call get_ome_inter_ex_sum_k(ibz,ek,xme_ex_band,vme_ex_band,berry_eigen_ex_band) !sum over k points - end if - end do - - if (iflag_norder.eq.2) then - !Sum all terms together for matrix elements - !fill Q_{NN'} (total), Y_{NN'} (total), X_{NN'} (total) and V_{NN'} (total) - do nn=1,norb_ex_cut - do nnp=1,norb_ex_cut - do nj=1,3 - qme_ex_inter(nj,nn,nnp)=qme_ex_inter1(nj,nn,nnp)+qme_ex_inter2(nj,nn,nnp) - yme_ex_inter(nj,nn,nnp)=yme_ex_inter1(nj,nn,nnp)+yme_ex_inter2(nj,nn,nnp) - xme_ex_inter(nj,nn,nnp)=yme_ex_inter(nj,nn,nnp)+qme_ex_inter(nj,nn,nnp) - vme_ex_inter(nj,nn,nnp)=vme_ex_inter1(nj,nn,nnp)+vme_ex_inter2(nj,nn,nnp) - end do - end do - end do - - !do nn=1,10 - !write(*,*) nn,e_ex(nn)*27.211386,qme_ex_inter(1,nn,nn),qme_ex_inter1(1,nn,nn),qme_ex_inter2(1,nn,nn) - !end do - !pause + if (iflag_norder .eq. 2) then + call get_ome_inter_ex_sum_k(ibz, ek, xme_ex_band, vme_ex_band, berry_eigen_ex_band) end if + end do + + if (do_write_exk) then + call write_ome_ex_linear_kresolved_close(u_exk) + deallocate(vme_ex_k) + write(*,*) ' k-resolved linear excitonic matrix elements written (omeexk)' + end if + + if (iflag_norder .eq. 2) then + ! Sum all terms together for matrix elements (N,N') + do nn = 1, norb_ex_cut + do nnp = 1, norb_ex_cut + do nj = 1, 3 + qme_ex_inter(nj, nn, nnp) = qme_ex_inter1(nj, nn, nnp) + qme_ex_inter2(nj, nn, nnp) + yme_ex_inter(nj, nn, nnp) = yme_ex_inter1(nj, nn, nnp) + yme_ex_inter2(nj, nn, nnp) + xme_ex_inter(nj, nn, nnp) = yme_ex_inter(nj, nn, nnp) + qme_ex_inter(nj, nn, nnp) + vme_ex_inter(nj, nn, nnp) = vme_ex_inter1(nj, nn, nnp) + vme_ex_inter2(nj, nn, nnp) + end do + end do + end do + end if - !do ibz=1,npointstotal - !write(*,*) ibz,fk_ex(ibz,1),fk_ex_der(1,ibz,1) - !end do + write(*,*) ' Optical matrix elements (ex) have been evaluated' - !do nn=1,norb_ex_cut - !write(*,*) nn,abs(qme_ex_inter1(1,1,nn)),abs(qme_ex_inter2(1,1,nn)) - !end do - !pause + ! writing excitonic optical matrix elements (summed over k) + if (iflag_norder .eq. 1) then + call write_ome_ex_linear(vme_ex) + end if + write(*,*) ' Optical matrix elements (ex, N -> GS) have been written in file' + write(*,*) ' Optical matrix elements (ex, N -> N^prime) will not be printed in this version' - write(*,*) ' Optical matrix elements (ex) have been evaluated' - !writing excitonic optical matrix elements - if (iflag_norder.eq.1) then - call write_ome_ex_linear(vme_ex) - end if - write(*,*) ' Optical matrix elements (ex, N -> GS) have been written in file' - write(*,*) ' Optical matrix elements (ex, N -> N^prime) will not be printed in this version' - end subroutine get_ome_ex +end subroutine get_ome_ex !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_ome_sp_xme_ex_band(ek,vme_ex_band,xme_ex_band) @@ -251,6 +268,41 @@ subroutine get_ome_gs_ex_sum_k(ibz,ek,xme_ex_band,vme_ex_band) end do end subroutine get_ome_gs_ex_sum_k + !----------------------------------------------------------------- + ! Compute k-resolved contribution to V_{0N} for a single k-point + ! This fills vme_ex_k(:,nn) with the contributions coming from + ! the given k-point `ibz`, using the same index mapping as + ! get_ome_gs_ex_sum_k but without summing over k-points. + !----------------------------------------------------------------- + subroutine get_ome_gs_ex_kresolved(ibz,ek,vme_ex_band,vme_ex_k) + implicit none + integer :: ibz + + integer :: nn,ic,iv,nj + integer :: i_ex_nn + + dimension ek(npointstotal,nband_ex) + dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) + complex*16 :: vme_ex_band + complex*16 :: vme_ex_k(3,norb_ex_cut) + real*8 :: ek + + ! initialize + vme_ex_k = 0.0d0 + + do nn = 1, norb_ex_cut + do ic = 1, nc_ex + do iv = 1, nv_ex + call get_ex_index_first(nf, nv_ex, nc_ex, 0, ibz, i_ex_nn, ic, iv) + do nj = 1, 3 + vme_ex_k(nj, nn) = vme_ex_k(nj, nn) + fk_ex(i_ex_nn, nn) * vme_ex_band(ibz, nj, iv, nv_ex+ic) + end do + end do + end do + end do + + end subroutine get_ome_gs_ex_kresolved + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine get_ome_inter_ex_sum_k(ibz,ek,xme_ex_band,vme_ex_band,berry_eigen_ex_band) implicit none @@ -361,6 +413,47 @@ subroutine write_ome_ex_linear(vme_ex) end subroutine write_ome_ex_linear !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_ome_ex_linear_kresolved_init(unitno, material_name, norb_ex_cut) + implicit none + integer, intent(in) :: unitno, norb_ex_cut + character(len=*), intent(in) :: material_name + + open(unitno, file='ome_linear_ex_k_'//trim(material_name)//'.omeexk', status='replace') + ! Simple header: + ! line 1: tag + ! line 2: norb_ex_cut + write(unitno,*) 1 + write(unitno,*) norb_ex_cut + end subroutine write_ome_ex_linear_kresolved_init + + + subroutine write_ome_ex_linear_kresolved_point(unitno, kx, ky, kz, norb_ex_cut, vme_ex_k) + implicit none + integer, intent(in) :: unitno, norb_ex_cut + real*8, intent(in) :: kx, ky, kz + complex*16, intent(in) :: vme_ex_k(3, norb_ex_cut) + + integer :: nn + ! For each k-point: + ! line: kx ky kz + ! then norb_ex_cut lines: + ! nn Re(Vx) Im(Vx) Re(Vy) Im(Vy) Re(Vz) Im(Vz) + write(unitno,*) kx, ky, kz + do nn = 1, norb_ex_cut + write(unitno,*) nn, dble(vme_ex_k(1,nn)), dimag(vme_ex_k(1,nn)), & + dble(vme_ex_k(2,nn)), dimag(vme_ex_k(2,nn)), & + dble(vme_ex_k(3,nn)), dimag(vme_ex_k(3,nn)) + end do + end subroutine write_ome_ex_linear_kresolved_point + + + subroutine write_ome_ex_linear_kresolved_close(unitno) + implicit none + integer, intent(in) :: unitno + close(unitno) + end subroutine write_ome_ex_linear_kresolved_close + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_ome_sp_linear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek) implicit none @@ -450,4 +543,7 @@ subroutine read_ome_sp_nonlinear(iflag_norder,npointstotal,nband_ex,berry_eigen_ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end subroutine read_ome_sp_nonlinear + + + end module ome_ex \ No newline at end of file From c24728666017d07563610085cb2135c366650f18 Mon Sep 17 00:00:00 2001 From: visitor Date: Thu, 19 Feb 2026 18:26:02 +0100 Subject: [PATCH 4/6] refactored xatu parser and input parser; added lorentzian broadening --- src/ome_sp.f90 | 1450 ++++++++++++++++---------------- src/parser_input_file.f90 | 195 +++-- src/parser_optics_xatu_dim.f90 | 810 +++++++++--------- src/sigma_first_ex.f90 | 18 +- src/sigma_first_sp.f90 | 15 +- src/sigma_second_ex.f90 | 20 +- src/sigma_second_sp.f90 | 16 +- 7 files changed, 1334 insertions(+), 1190 deletions(-) diff --git a/src/ome_sp.f90 b/src/ome_sp.f90 index 2a82475..a747fe2 100644 --- a/src/ome_sp.f90 +++ b/src/ome_sp.f90 @@ -1,744 +1,744 @@ module ome_sp - use constants_math - use parser_wannier90_tb, & - only:material_name,nR,nRvec,norb,R,shop,hhop,rhop_c - use parser_optics_xatu_dim, & - only:npointstotal,rkxvector,rkyvector,rkzvector, & - nband_ex,nband_index,nv_ex,nc_ex - - implicit none - - allocatable vme_ex_band(:,:,:,:) - allocatable ek(:,:) - allocatable gen_der_ex_band(:,:,:,:,:) - allocatable shift_vector_ex_band(:,:,:,:,:) - allocatable berry_eigen_ex_band(:,:,:,:) - - real*8 ek - real*8 shift_vector_ex_band - complex*16 vme_ex_band - complex*16 gen_der_ex_band - complex*16 berry_eigen_ex_band - contains -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_ome_sp(iflag_norder) - implicit none - - integer iflag_norder - integer ibz,ibz_sum - integer i,j,ii,jj,nj - integer naux1 - - allocatable :: skernel(:,:) - allocatable :: hkernel(:,:) - allocatable :: sderkernel(:,:,:) - allocatable :: hderkernel(:,:,:) - allocatable :: akernel(:,:,:) - - allocatable :: e(:) - allocatable :: ecomplex(:) - allocatable :: hk_ev(:,:) - allocatable :: vme(:,:,:) - - allocatable :: abc(:,:,:) - allocatable :: gd1(:,:,:,:) - allocatable :: gd2(:,:,:,:) - allocatable :: gd3(:,:,:,:) - allocatable :: gen_der(:,:,:,:) - allocatable :: vme_der(:,:,:,:) - allocatable :: shift_vector(:,:,:,:) - allocatable :: berry_eigen1(:,:,:) - allocatable :: berry_eigen2(:,:,:) - allocatable :: berry_eigen(:,:,:) - - real(8) rkx,rky,rkz - real(8) :: e - real(8) :: shift_vector - - complex*16 :: skernel,hkernel,sderkernel,hderkernel,akernel - complex*16 :: hk_ev - complex*16 :: ecomplex - complex*16 :: vjseudoa,vjseudob - complex*16 :: vme - complex*16 :: vme_der - complex*16 :: berry_eigen1,berry_eigen2,berry_eigen - - complex*16 abc,gd1,gd2,gd3,gen_der -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - write(*,*) '5. Entering ome_sp' - !allocate to save k-dependent vme and energies between involved band - allocate(vme_ex_band(npointstotal,3,nband_ex,nband_ex)) - allocate(ek(npointstotal,nband_ex)) - allocate(berry_eigen_ex_band(npointstotal,3,nband_ex,nband_ex)) - allocate(gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex)) - allocate(shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex)) - gen_der_ex_band=0.0d0 - shift_vector_ex_band=0.0d0 - berry_eigen_ex_band=0.0d0 - vme_ex_band=0.0d0 - ek=0.0d0 - - - allocate(skernel(norb,norb)) - allocate(hkernel(norb,norb)) - allocate(sderkernel(3,norb,norb)) - allocate(hderkernel(3,norb,norb)) - allocate(akernel(3,norb,norb)) - - allocate(e(norb)) - allocate(ecomplex(norb)) - allocate(hk_ev(norb,norb)) - allocate(vme(3,norb,norb)) - - allocate(abc(3,norb,norb)) - allocate(gd1(3,3,norb,norb)) - allocate(gd2(3,3,norb,norb)) - allocate(gd3(3,3,norb,norb)) - allocate(gen_der(3,3,norb,norb)) - allocate(vme_der(3,3,norb,norb)) - allocate(shift_vector(3,3,norb,norb)) - allocate(berry_eigen1(3,norb,norb)) - allocate(berry_eigen2(3,norb,norb)) - allocate(berry_eigen(3,norb,norb)) - - !if (iflag_norder.eq.2) then + use constants_math + use parser_wannier90_tb, & + only:material_name,nR,nRvec,norb,R,shop,hhop,rhop_c + use parser_optics_xatu_dim, & + only:npointstotal,rkxvector,rkyvector,rkzvector, & + nband_ex,nband_index,nv_ex,nc_ex + + implicit none + + allocatable vme_ex_band(:,:,:,:) + allocatable ek(:,:) + allocatable gen_der_ex_band(:,:,:,:,:) + allocatable shift_vector_ex_band(:,:,:,:,:) + allocatable berry_eigen_ex_band(:,:,:,:) + + real*8 ek + real*8 shift_vector_ex_band + complex*16 vme_ex_band + complex*16 gen_der_ex_band + complex*16 berry_eigen_ex_band +contains + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_ome_sp(iflag_norder) + implicit none + + integer iflag_norder + integer ibz,ibz_sum + integer i,j,ii,jj,nj + integer naux1 + + allocatable :: skernel(:,:) + allocatable :: hkernel(:,:) + allocatable :: sderkernel(:,:,:) + allocatable :: hderkernel(:,:,:) + allocatable :: akernel(:,:,:) + + allocatable :: e(:) + allocatable :: ecomplex(:) + allocatable :: hk_ev(:,:) + allocatable :: vme(:,:,:) + + allocatable :: abc(:,:,:) + allocatable :: gd1(:,:,:,:) + allocatable :: gd2(:,:,:,:) + allocatable :: gd3(:,:,:,:) + allocatable :: gen_der(:,:,:,:) + allocatable :: vme_der(:,:,:,:) + allocatable :: shift_vector(:,:,:,:) + allocatable :: berry_eigen1(:,:,:) + allocatable :: berry_eigen2(:,:,:) + allocatable :: berry_eigen(:,:,:) + + real(8) rkx,rky,rkz + real(8) :: e + real(8) :: shift_vector + + complex*16 :: skernel,hkernel,sderkernel,hderkernel,akernel + complex*16 :: hk_ev + complex*16 :: ecomplex + complex*16 :: vjseudoa,vjseudob + complex*16 :: vme + complex*16 :: vme_der + complex*16 :: berry_eigen1,berry_eigen2,berry_eigen + + complex*16 abc,gd1,gd2,gd3,gen_der + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + write(*,*) '5. Entering ome_sp' + !allocate to save k-dependent vme and energies between involved band + allocate(vme_ex_band(npointstotal,3,nband_ex,nband_ex)) + allocate(ek(npointstotal,nband_ex)) + allocate(berry_eigen_ex_band(npointstotal,3,nband_ex,nband_ex)) + allocate(gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex)) + allocate(shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex)) + gen_der_ex_band=0.0d0 + shift_vector_ex_band=0.0d0 + berry_eigen_ex_band=0.0d0 + vme_ex_band=0.0d0 + ek=0.0d0 + + + allocate(skernel(norb,norb)) + allocate(hkernel(norb,norb)) + allocate(sderkernel(3,norb,norb)) + allocate(hderkernel(3,norb,norb)) + allocate(akernel(3,norb,norb)) + + allocate(e(norb)) + allocate(ecomplex(norb)) + allocate(hk_ev(norb,norb)) + allocate(vme(3,norb,norb)) + + allocate(abc(3,norb,norb)) + allocate(gd1(3,3,norb,norb)) + allocate(gd2(3,3,norb,norb)) + allocate(gd3(3,3,norb,norb)) + allocate(gen_der(3,3,norb,norb)) + allocate(vme_der(3,3,norb,norb)) + allocate(shift_vector(3,3,norb,norb)) + allocate(berry_eigen1(3,norb,norb)) + allocate(berry_eigen2(3,norb,norb)) + allocate(berry_eigen(3,norb,norb)) + + !if (iflag_norder.eq.2) then !allocate (gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex)) !allocate (shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex)) !gen_der_ex_band=0.0d0 !shift_vector_ex_band=0.0d0 - !end if - !open(50,file='coefs_new.dat') - !Brillouin zone sampling - parallelization - - !write - !only:material_name,nR,nRvec,norb,R,shop,hhop,rhop_c - !use parser_optics_xatu_dim, & - !only:npointstotal,rkxvector,rkyvector,rkzvector, & - !nband_ex,nband_index,nv_ex,nc_ex - !write(*,*) material_name - !write(*,*) nR - !do i=1,nR + !end if + !open(50,file='coefs_new.dat') + !Brillouin zone sampling - parallelization + + !write + !only:material_name,nR,nRvec,norb,R,shop,hhop,rhop_c + !use parser_optics_xatu_dim, & + !only:npointstotal,rkxvector,rkyvector,rkzvector, & + !nband_ex,nband_index,nv_ex,nc_ex + !write(*,*) material_name + !write(*,*) nR + !do i=1,nR !write(*,*) nRvec(i,1),nRvec(i,2),nRvec(i,3) - !end do - !do i=1,npointstotal + !end do + !do i=1,npointstotal !write(*,*) rkxvector(i),rkyvector(i),rkzvector(i) - !end do - !pause - - write(*,*) ' Calculating optical matrix elements (sp): sampling BZ...' - ibz_sum=0 !counter for the number of k points in the BZ - - !initializing variables - !$OMP PARALLEL DO PRIVATE(rkx,rky,rkz), & - !$OMP PRIVATE(hkernel,skernel,sderkernel,hderkernel,akernel), & - !$OMP PRIVATE(hk_ev,e,vme), & - !$OMP PRIVATE(abc,gen_der,gd1,gd2,gd3), & - !$OMP PRIVATE(vme_der,shift_vector,berry_eigen1,berry_eigen2,berry_eigen) - do ibz=1,npointstotal - !!$OMP CRITICAL - !ibz_sum=ibz_sum+1 - !call percentage_index(ibz_sum,npointstotal,naux1) - !!$OMP END CRITICAL - write(*,*) ' Optical matrix elements (sp): k-point',ibz,'/',npointstotal + !end do !pause - rkx=rkxvector(ibz) - rky=rkyvector(ibz) - rkz=rkzvector(ibz) - !pause - - !get matrices in the \alpha, \alpha' basis (orbitals,k) - call get_vme_kernels_ome(rkx,rky,rkz,norb,skernel,sderkernel, & - hkernel,hderkernel,akernel) - call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev,e,vme) + write(*,*) ' Calculating optical matrix elements (sp): sampling BZ...' + ibz_sum=0 !counter for the number of k points in the BZ + + !initializing variables + !$OMP PARALLEL DO PRIVATE(rkx,rky,rkz), & + !$OMP PRIVATE(hkernel,skernel,sderkernel,hderkernel,akernel), & + !$OMP PRIVATE(hk_ev,e,vme), & + !$OMP PRIVATE(abc,gen_der,gd1,gd2,gd3), & + !$OMP PRIVATE(vme_der,shift_vector,berry_eigen1,berry_eigen2,berry_eigen) + do ibz=1,npointstotal + !!$OMP CRITICAL + !ibz_sum=ibz_sum+1 + !call percentage_index(ibz_sum,npointstotal,naux1) + !!$OMP END CRITICAL + write(*,*) ' Optical matrix elements (sp): k-point',ibz,'/',npointstotal + !pause + rkx=rkxvector(ibz) + rky=rkyvector(ibz) + rkz=rkzvector(ibz) + + !pause + + !get matrices in the \alpha, \alpha' basis (orbitals,k) + call get_vme_kernels_ome(rkx,rky,rkz,norb,skernel,sderkernel, & + hkernel,hderkernel,akernel) + call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev,e,vme) + + !matrix elements for second-order calculation + if (iflag_norder.eq.2) then + call get_gen_der_sumrule(norb,vme,e,abc,gen_der,gd1,gd2,gd3) + call get_berry_eigen_fourpoint(rkx,rky,rkz,norb,vme_der, & + shift_vector,berry_eigen1,berry_eigen2,berry_eigen) + end if + + !write(*,*) ibz,rkx,rky,shift_vector(1,1,2,3),berry_eigen1(1,2,3),berry_eigen2(1,2,3) + + !saving eigenvalues and optical matrix elements at this k point for the bandlist + do i=1,nband_ex + ii=nband_index(i) + ek(ibz,i)=e(ii) + do nj=1,3 + do j=1,nband_ex + jj=nband_index(j) + vme_ex_band(ibz,nj,i,j)=vme(nj,ii,jj) + + if (iflag_norder.eq.2) then + shift_vector_ex_band(ibz,nj,1,i,j)=shift_vector(nj,1,ii,jj) + shift_vector_ex_band(ibz,nj,2,i,j)=shift_vector(nj,2,ii,jj) + shift_vector_ex_band(ibz,nj,3,i,j)=shift_vector(nj,3,ii,jj) + + gen_der_ex_band(ibz,nj,1,i,j)=gen_der(nj,1,ii,jj) + gen_der_ex_band(ibz,nj,2,i,j)=gen_der(nj,2,ii,jj) + gen_der_ex_band(ibz,nj,3,i,j)=gen_der(nj,3,ii,jj) + + berry_eigen_ex_band(ibz,nj,i,j)=berry_eigen(nj,ii,jj) + berry_eigen_ex_band(ibz,nj,i,j)=berry_eigen(nj,ii,jj) + berry_eigen_ex_band(ibz,nj,i,j)=berry_eigen(nj,ii,jj) + + end if + + end do + end do + end do + + !pause + end do + !$OMP END PARALLEL DO + !close(50) - !matrix elements for second-order calculation + !write matrix elements into file + write(*,*) ' Writing optical matrix elements (sp) into file' + if (iflag_norder.eq.1) then + call write_ome_sp_linear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek) + end if if (iflag_norder.eq.2) then - call get_gen_der_sumrule(norb,vme,e,abc,gen_der,gd1,gd2,gd3) - call get_berry_eigen_fourpoint(rkx,rky,rkz,norb,vme_der, & - shift_vector,berry_eigen1,berry_eigen2,berry_eigen) + call write_ome_sp_nonlinear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek, & + gen_der_ex_band,shift_vector_ex_band,berry_eigen_ex_band) end if - - !write(*,*) ibz,rkx,rky,shift_vector(1,1,2,3),berry_eigen1(1,2,3),berry_eigen2(1,2,3) - - !saving eigenvalues and optical matrix elements at this k point for the bandlist - do i=1,nband_ex - ii=nband_index(i) - ek(ibz,i)=e(ii) - do nj=1,3 - do j=1,nband_ex - jj=nband_index(j) - vme_ex_band(ibz,nj,i,j)=vme(nj,ii,jj) - - if (iflag_norder.eq.2) then - shift_vector_ex_band(ibz,nj,1,i,j)=shift_vector(nj,1,ii,jj) - shift_vector_ex_band(ibz,nj,2,i,j)=shift_vector(nj,2,ii,jj) - shift_vector_ex_band(ibz,nj,3,i,j)=shift_vector(nj,3,ii,jj) - - gen_der_ex_band(ibz,nj,1,i,j)=gen_der(nj,1,ii,jj) - gen_der_ex_band(ibz,nj,2,i,j)=gen_der(nj,2,ii,jj) - gen_der_ex_band(ibz,nj,3,i,j)=gen_der(nj,3,ii,jj) - - berry_eigen_ex_band(ibz,nj,i,j)=berry_eigen(nj,ii,jj) - berry_eigen_ex_band(ibz,nj,i,j)=berry_eigen(nj,ii,jj) - berry_eigen_ex_band(ibz,nj,i,j)=berry_eigen(nj,ii,jj) - - end if - - end do - end do + write(*,*) ' Optical matrix elements (sp) have been written in file' + !pause + !deallocate(vme_ex_band) + !deallocate(ek) + end subroutine get_ome_sp + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !you have to addapt this to 3D + subroutine get_berry_eigen_fourpoint(rkx,rky,rkz,norb,vme_der, & + shift_vector,berry_eigen1,berry_eigen2,berry_eigen) + implicit none + + integer nR,norb + integer nn,nnp + integer ialpha,ialphap + integer nj,njp + + dimension hkernel(norb,norb),skernel(norb,norb) + dimension sderkernel(3,norb,norb),hderkernel(3,norb,norb) + + dimension hk_alpha(norb,norb),hk_ev(norb,norb),e(norb) + dimension akernel(3,norb,norb) + + dimension vjseudoa(3,norb,norb),vjseudob(3,norb,norb),vme(3,norb,norb) + dimension berry_eigen1(3,norb,norb),berry_eigen2(3,norb,norb) + dimension berry_eigen(3,norb,norb) + dimension vme_der(3,3,norb,norb) + dimension vme_der_phase(3,3,norb,norb) + dimension hk_ev_neigh(5,norb,norb) + dimension vme_neigh(5,3,norb,norb) + + dimension shift_vector(3,3,norb,norb) + + real*8 rkx,rky,rkz,rkx_neigh,rky_neigh,rkz_neigh + real*8 e + real*8 vme_der_phase + real*8 shift_vector + real*8 ph1,ph2,ph3,ph4 + + complex*16 hkernel,akernel,skernel,sderkernel,hderkernel + complex*16 hk_alpha,hk_ev,vjseudoa,vjseudob,vme + complex*16 amu,amup,aux1,aux2,aux3,aux4,factor + complex*16 vme_der + complex*16 berry_eigen1,berry_eigen2,berry_eigen + complex*16 hk_ev_neigh,vme_neigh + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + vme=0.0d0 + + hk_ev_neigh=0.0d0 + vme_neigh=0.0d0 + + berry_eigen1=0.0d0 + berry_eigen2=0.0d0 + berry_eigen=0.0d0 + vme_der=0.0d0 + shift_vector=0.0d0 + + rkx_neigh=rkx+dk + rky_neigh=rky + rkz_neigh=rkz + call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & + sderkernel,hkernel,hderkernel,akernel) + !velocity matrix elements + call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev_neigh(1,:,:),e,vme_neigh(1,:,:,:)) + + !lets try to do only a n,n' loop + rkx_neigh=rkx + rky_neigh=rky-dk + rkz_neigh=rkz + call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & + sderkernel,hkernel,hderkernel,akernel) + !velocity matrix elements + call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev_neigh(2,:,:),e,vme_neigh(2,:,:,:)) + + !lets try to do only a n,n' loop + rkx_neigh=rkx-dk + rky_neigh=rky + rkz_neigh=rkz + call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & + sderkernel,hkernel,hderkernel,akernel) + !velocity matrix elements + call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev_neigh(3,:,:),e,vme_neigh(3,:,:,:)) + + !lets try to do only a n,n' loop + rkx_neigh=rkx + rky_neigh=rky+dk + rkz_neigh=rkz + call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & + sderkernel,hkernel,hderkernel,akernel) + !velocity matrix elements + call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev_neigh(4,:,:),e,vme_neigh(4,:,:,:)) + + !lets try to do only a n,n' loop + rkx_neigh=rkx + rky_neigh=rky + rkz_neigh=rkz + call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & + sderkernel,hkernel,hderkernel,akernel) + !velocity matrix elements + call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev_neigh(5,:,:),e,vme_neigh(5,:,:,:)) + + + + do nn=1,norb + do nnp=1,norb + !Direct evaluation of Berry connection. See Esteve-Paredes and Palacios, Scipost 2023 + do ialpha=1,norb + do ialphap=1,norb + !x-dir + aux1=(hk_ev_neigh(1,ialphap,nnp)-hk_ev_neigh(3,ialphap,nnp))/(2.0d0*dk) + berry_eigen1(1,nn,nnp)=berry_eigen1(1,nn,nnp)+ & + complex(0.0d0,1.0d0)*conjg(hk_ev_neigh(5,ialpha,nn))*skernel(ialpha,ialphap)*aux1 + berry_eigen2(1,nn,nnp)=berry_eigen2(1,nn,nnp)+ & + conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(1,ialpha,ialphap) + + !if (nn.eq.1 .and. nnp.eq.1 .and. ialpha.eq.1 .and. ialphap.eq.1) then + !write(*,*) ialpha,ialphap,berry_eigen1(1,nn,nnp),berry_eigen2(1,nn,nnp),akernel(1,ialpha,ialphap) + !pause + !end if + + !y-dir + aux1=(hk_ev_neigh(4,ialphap,nnp)-hk_ev_neigh(2,ialphap,nnp))/(2.0d0*dk) + berry_eigen1(2,nn,nnp)=berry_eigen1(2,nn,nnp)+ & + complex(0.0d0,1.0d0)*conjg(hk_ev_neigh(5,ialpha,nn))*skernel(ialpha,ialphap)*aux1 + berry_eigen2(2,nn,nnp)=berry_eigen2(2,nn,nnp)+ & + conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(2,ialpha,ialphap) + + !z-dir + berry_eigen1(3,nn,nnp)=0.0d0 + berry_eigen2(3,nn,nnp)=berry_eigen2(3,nn,nnp)+ & + conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(3,ialpha,ialphap) + end do + end do + do nj=1,3 + !if BC's is diverging at degenrate points, set it to zero + berry_eigen(nj,nn,nnp)=berry_eigen1(nj,nn,nnp)+berry_eigen2(nj,nn,nnp) + if (abs(berry_eigen(nj,nn,nnp)).gt.50.0d0) then + berry_eigen(nj,nn,nnp)=0.0d0 + end if + + aux1=vme_neigh(1,nj,nn,nnp) + aux3=vme_neigh(3,nj,nn,nnp) + aux4=vme_neigh(4,nj,nn,nnp) + aux2=vme_neigh(2,nj,nn,nnp) + vme_der(1,nj,nn,nnp)=(aux1-aux3)/(2.0d0*dk) + vme_der(2,nj,nn,nnp)=(aux4-aux2)/(2.0d0*dk) + + call get_phase(vme_neigh(1,nj,nn,nnp),ph1) + call get_phase(vme_neigh(3,nj,nn,nnp),ph3) + call get_phase(vme_neigh(4,nj,nn,nnp),ph4) + call get_phase(vme_neigh(2,nj,nn,nnp),ph2) + vme_der_phase(1,nj,nn,nnp)=(ph1-ph3)/(2.0d0*dk) + vme_der_phase(2,nj,nn,nnp)=(ph4-ph2)/(2.0d0*dk) + end do + end do end do - !pause - end do - !$OMP END PARALLEL DO - !close(50) - - !write matrix elements into file - write(*,*) ' Writing optical matrix elements (sp) into file' - if (iflag_norder.eq.1) then - call write_ome_sp_linear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek) - end if - if (iflag_norder.eq.2) then - call write_ome_sp_nonlinear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek, & - gen_der_ex_band,shift_vector_ex_band,berry_eigen_ex_band) - end if - write(*,*) ' Optical matrix elements (sp) have been written in file' - !pause - !deallocate(vme_ex_band) - !deallocate(ek) - end subroutine get_ome_sp -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !you have to addapt this to 3D - subroutine get_berry_eigen_fourpoint(rkx,rky,rkz,norb,vme_der, & - shift_vector,berry_eigen1,berry_eigen2,berry_eigen) - implicit none - - integer nR,norb - integer nn,nnp - integer ialpha,ialphap - integer nj,njp - - dimension hkernel(norb,norb),skernel(norb,norb) - dimension sderkernel(3,norb,norb),hderkernel(3,norb,norb) - - dimension hk_alpha(norb,norb),hk_ev(norb,norb),e(norb) - dimension akernel(3,norb,norb) - - dimension vjseudoa(3,norb,norb),vjseudob(3,norb,norb),vme(3,norb,norb) - dimension berry_eigen1(3,norb,norb),berry_eigen2(3,norb,norb) - dimension berry_eigen(3,norb,norb) - dimension vme_der(3,3,norb,norb) - dimension vme_der_phase(3,3,norb,norb) - dimension hk_ev_neigh(5,norb,norb) - dimension vme_neigh(5,3,norb,norb) - - dimension shift_vector(3,3,norb,norb) - - real*8 rkx,rky,rkz,rkx_neigh,rky_neigh,rkz_neigh - real*8 e - real*8 vme_der_phase - real*8 shift_vector - real*8 ph1,ph2,ph3,ph4 - - complex*16 hkernel,akernel,skernel,sderkernel,hderkernel - complex*16 hk_alpha,hk_ev,vjseudoa,vjseudob,vme - complex*16 amu,amup,aux1,aux2,aux3,aux4,factor - complex*16 vme_der - complex*16 berry_eigen1,berry_eigen2,berry_eigen - complex*16 hk_ev_neigh,vme_neigh -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - vme=0.0d0 - - hk_ev_neigh=0.0d0 - vme_neigh=0.0d0 - - berry_eigen1=0.0d0 - berry_eigen2=0.0d0 - berry_eigen=0.0d0 - vme_der=0.0d0 - shift_vector=0.0d0 - - rkx_neigh=rkx+dk - rky_neigh=rky - rkz_neigh=rkz - call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & - sderkernel,hkernel,hderkernel,akernel) - !velocity matrix elements - call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev_neigh(1,:,:),e,vme_neigh(1,:,:,:)) - - !lets try to do only a n,n' loop - rkx_neigh=rkx - rky_neigh=rky-dk - rkz_neigh=rkz - call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & - sderkernel,hkernel,hderkernel,akernel) - !velocity matrix elements - call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev_neigh(2,:,:),e,vme_neigh(2,:,:,:)) - - !lets try to do only a n,n' loop - rkx_neigh=rkx-dk - rky_neigh=rky - rkz_neigh=rkz - call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & - sderkernel,hkernel,hderkernel,akernel) - !velocity matrix elements - call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev_neigh(3,:,:),e,vme_neigh(3,:,:,:)) - - !lets try to do only a n,n' loop - rkx_neigh=rkx - rky_neigh=rky+dk - rkz_neigh=rkz - call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & - sderkernel,hkernel,hderkernel,akernel) - !velocity matrix elements - call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev_neigh(4,:,:),e,vme_neigh(4,:,:,:)) - - !lets try to do only a n,n' loop - rkx_neigh=rkx - rky_neigh=rky - rkz_neigh=rkz - call get_vme_kernels_ome(rkx_neigh,rky_neigh,rkz_neigh,norb,skernel, & - sderkernel,hkernel,hderkernel,akernel) - !velocity matrix elements - call get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev_neigh(5,:,:),e,vme_neigh(5,:,:,:)) - - - - do nn=1,norb - do nnp=1,norb - !Direct evaluation of Berry connection. See Esteve-Paredes and Palacios, Scipost 2023 - do ialpha=1,norb - do ialphap=1,norb - !x-dir - aux1=(hk_ev_neigh(1,ialphap,nnp)-hk_ev_neigh(3,ialphap,nnp))/(2.0d0*dk) - berry_eigen1(1,nn,nnp)=berry_eigen1(1,nn,nnp)+ & - complex(0.0d0,1.0d0)*conjg(hk_ev_neigh(5,ialpha,nn))*skernel(ialpha,ialphap)*aux1 - berry_eigen2(1,nn,nnp)=berry_eigen2(1,nn,nnp)+ & - conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(1,ialpha,ialphap) - - !if (nn.eq.1 .and. nnp.eq.1 .and. ialpha.eq.1 .and. ialphap.eq.1) then - !write(*,*) ialpha,ialphap,berry_eigen1(1,nn,nnp),berry_eigen2(1,nn,nnp),akernel(1,ialpha,ialphap) - !pause - !end if - - !y-dir - aux1=(hk_ev_neigh(4,ialphap,nnp)-hk_ev_neigh(2,ialphap,nnp))/(2.0d0*dk) - berry_eigen1(2,nn,nnp)=berry_eigen1(2,nn,nnp)+ & - complex(0.0d0,1.0d0)*conjg(hk_ev_neigh(5,ialpha,nn))*skernel(ialpha,ialphap)*aux1 - berry_eigen2(2,nn,nnp)=berry_eigen2(2,nn,nnp)+ & - conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(2,ialpha,ialphap) - - !z-dir - berry_eigen1(3,nn,nnp)=0.0d0 - berry_eigen2(3,nn,nnp)=berry_eigen2(3,nn,nnp)+ & - conjg(hk_ev_neigh(5,ialpha,nn))*hk_ev_neigh(5,ialphap,nnp)*akernel(3,ialpha,ialphap) - end do - end do - do nj=1,3 - !if BC's is diverging at degenrate points, set it to zero - berry_eigen(nj,nn,nnp)=berry_eigen1(nj,nn,nnp)+berry_eigen2(nj,nn,nnp) - if (abs(berry_eigen(nj,nn,nnp)).gt.50.0d0) then - berry_eigen(nj,nn,nnp)=0.0d0 - end if - - aux1=vme_neigh(1,nj,nn,nnp) - aux3=vme_neigh(3,nj,nn,nnp) - aux4=vme_neigh(4,nj,nn,nnp) - aux2=vme_neigh(2,nj,nn,nnp) - vme_der(1,nj,nn,nnp)=(aux1-aux3)/(2.0d0*dk) - vme_der(2,nj,nn,nnp)=(aux4-aux2)/(2.0d0*dk) - - call get_phase(vme_neigh(1,nj,nn,nnp),ph1) - call get_phase(vme_neigh(3,nj,nn,nnp),ph3) - call get_phase(vme_neigh(4,nj,nn,nnp),ph4) - call get_phase(vme_neigh(2,nj,nn,nnp),ph2) - vme_der_phase(1,nj,nn,nnp)=(ph1-ph3)/(2.0d0*dk) - vme_der_phase(2,nj,nn,nnp)=(ph4-ph2)/(2.0d0*dk) - end do - end do - end do - - !shift vector - do nn=1,norb - do nnp=1,norb - do nj=1,3 - do njp=1,3 - shift_vector(nj,njp,nn,nnp)=-vme_der_phase(nj,njp,nn,nnp) & - +(realpart(berry_eigen(nj,nn,nn))-realpart(berry_eigen(nj,nnp,nnp))) - if (abs(shift_vector(nj,njp,nn,nnp)).gt.50.0d0) then - shift_vector(nj,njp,nn,nnp)=0.0d0 - end if - vme_der(nj,njp,nn,nnp)=vme_der(nj,njp,nn,nnp) & - -complex(0.0d0,1.0d0)*vme_neigh(5,njp,nn,nnp) & - *(realpart(berry_eigen(nj,nn,nn))-realpart(berry_eigen(nj,nnp,nnp))) - end do - end do - end do - end do - - end subroutine get_berry_eigen_fourpoint -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_phase(aux1,ph) - real*8 ph - complex*16 aux1 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - if (abs(aux1).lt.10d-6) then - ph=0.0d0 - else - ph=aimag(log(aux1)) - end if - end subroutine get_phase -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_gen_der_sumrule(norb,vme,e,abc,gen_der,gd1,gd2,gd3) - implicit real*8 (a-h,o-z) - - integer norb,norb_inter_cut - integer nn,nnp,nnpp - integer nj,njp - - dimension e(norb) - dimension vme(3,norb,norb) - dimension gen_der(3,3,norb,norb) - dimension gd1(3,3,norb,norb) - dimension gd2(3,3,norb,norb) - dimension gd3(3,3,norb,norb) - dimension abc(3,norb,norb) - - real*8 e - complex*16 vme,gen_der,abc,gd1,gd2,gd3,aux1,aux2,aux3 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - abc=0.0d0 - gd1=0.0d0 - gd2=0.0d0 - gd3=0.0d0 - gen_der=0.0d0 - norb_inter_cut=norb - !write(*,*) 'computing berry connection, generalized derivative' - do nn=1,norb - do nnp=1,norb - do nj=1,3 - !if (nn.eq.nnp) then - if (abs(e(nn)-e(nnp)).lt.1.0d-05) then - abc(nj,nn,nnp)=0.0d0 - else - abc(nj,nn,nnp)=-complex(0.0d0,1.0d0)*vme(nj,nn,nnp)/(e(nn)-e(nnp)) - !abc(nj,nnp,nn)=conjg(abc(nj,nn,nnp)) - end if - do njp=1,3 - !if (nn.eq.nnp) then - if (abs(e(nn)-e(nnp)).lt.1.0d-05) then - gd1(nj,njp,nn,nnp)=0.0d0 - gd2(nj,njp,nn,nnp)=0.0d0 - gd3(nj,njp,nn,nnp)=0.0d0 - else - gd1(nj,njp,nn,nnp)=complex(0.0d0,1.0d0)/((e(nn)-e(nnp))**2)* & - (vme(nj,nn,nnp)*(vme(njp,nn,nn)-vme(njp,nnp,nnp))+ & - vme(njp,nn,nnp)*(vme(nj,nn,nn)-vme(nj,nnp,nnp))) - - gd2(nj,njp,nn,nnp)=0.0d0 - gd3(nj,njp,nn,nnp)=0.0d0 - do nnpp=1,norb_inter_cut - !if (nnpp.eq.nn .or. nnpp.eq.nnp) then - if (abs(e(nnpp)-e(nn)).lt.1.0d-05 .or. & - abs(e(nnpp)-e(nnp)).lt.1.0d-05) then - gd2(nj,njp,nn,nnp)=gd2(nj,njp,nn,nnp)+0.0d0 - gd3(nj,njp,nn,nnp)=gd3(nj,njp,nn,nnp)+0.0d0 - else - gd2(nj,njp,nn,nnp)=gd2(nj,njp,nn,nnp)+ & - complex(0.0d0,1.0d0)/(e(nn)-e(nnp))* & - (vme(nj,nn,nnpp)*vme(njp,nnpp,nnp)/(e(nnpp)-e(nnp))) - - gd3(nj,njp,nn,nnp)=gd3(nj,njp,nn,nnp)+ & - complex(0.0d0,1.0d0)/(e(nn)-e(nnp))* & - (-vme(njp,nn,nnpp)*vme(nj,nnpp,nnp)/(e(nn)-e(nnpp))) - end if - end do - !momentums and A and B term - end if - gen_der(nj,njp,nn,nnp)=gd1(nj,njp,nn,nnp)+gd2(nj,njp,nn,nnp) & - +gd3(nj,njp,nn,nnp) - !gen_der(nj,njp,nnp,nn)=conjg(gen_der(nj,njp,nn,nnp)) - end do - end do - end do - end do - - end subroutine get_gen_der_sumrule - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine write_ome_sp_linear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek) - implicit none - integer iflag_norder - integer npointstotal,nband_ex - integer ibz - integer nj,i,j - - dimension ek(npointstotal,nband_ex) - dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) - - real*8 ek - complex*16 vme_ex_band -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - open(10,file='ome_linear_sp_'//trim(material_name)//'.omesp') - write(10,*) iflag_norder - do ibz=1,npointstotal - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(ek(ibz,j),j=1,nband_ex) - do i=1,nband_ex - do j=1,nband_ex - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz), & - (realpart(vme_ex_band(ibz,nj,i,j)),aimag(vme_ex_band(ibz,nj,i,j)), nj=1,3) - end do + !shift vector + do nn=1,norb + do nnp=1,norb + do nj=1,3 + do njp=1,3 + shift_vector(nj,njp,nn,nnp)=-vme_der_phase(nj,njp,nn,nnp) & + +(realpart(berry_eigen(nj,nn,nn))-realpart(berry_eigen(nj,nnp,nnp))) + if (abs(shift_vector(nj,njp,nn,nnp)).gt.50.0d0) then + shift_vector(nj,njp,nn,nnp)=0.0d0 + end if + vme_der(nj,njp,nn,nnp)=vme_der(nj,njp,nn,nnp) & + -complex(0.0d0,1.0d0)*vme_neigh(5,njp,nn,nnp) & + *(realpart(berry_eigen(nj,nn,nn))-realpart(berry_eigen(nj,nnp,nnp))) + end do + end do + end do end do - end do - close(10) - end subroutine write_ome_sp_linear - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine write_ome_sp_nonlinear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek,gen_der_ex_band, & - shift_vector_ex_band,berry_eigen_ex_band) - implicit none - integer iflag_norder - integer npointstotal,nband_ex - integer ibz - integer nj,i,j - - dimension ek(npointstotal,nband_ex) - dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) - dimension berry_eigen_ex_band(npointstotal,3,nband_ex,nband_ex) - dimension gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex) - dimension shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex) - - real*8 ek - real*8 shift_vector_ex_band - complex*16 vme_ex_band - complex*16 berry_eigen_ex_band - complex*16 gen_der_ex_band -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - open(10,file='ome_nonlinear_sp_'//trim(material_name)//'.omesp') - write(10,*) iflag_norder - do ibz=1,npointstotal - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(ek(ibz,j),j=1,nband_ex) - do i=1,nband_ex - do j=1,nband_ex - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz), & - (realpart(vme_ex_band(ibz,nj,i,j)),aimag(vme_ex_band(ibz,nj,i,j)), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(berry_eigen_ex_band(ibz,nj,i,j)), & - aimag(berry_eigen_ex_band(ibz,nj,i,j)), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(shift_vector_ex_band(ibz,1,nj,i,j), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(shift_vector_ex_band(ibz,2,nj,i,j), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(shift_vector_ex_band(ibz,3,nj,i,j), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(gen_der_ex_band(ibz,1,nj,i,j)), & - aimag(gen_der_ex_band(ibz,1,nj,i,j)), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(gen_der_ex_band(ibz,2,nj,i,j)), & - aimag(gen_der_ex_band(ibz,2,nj,i,j)), nj=1,3) - write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(gen_der_ex_band(ibz,3,nj,i,j)), & - aimag(gen_der_ex_band(ibz,3,nj,i,j)), nj=1,3) - end do + + end subroutine get_berry_eigen_fourpoint + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_phase(aux1,ph) + real*8 ph + complex*16 aux1 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (abs(aux1).lt.10d-6) then + ph=0.0d0 + else + ph=aimag(log(aux1)) + end if + end subroutine get_phase + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_gen_der_sumrule(norb,vme,e,abc,gen_der,gd1,gd2,gd3) + implicit real*8 (a-h,o-z) + + integer norb,norb_inter_cut + integer nn,nnp,nnpp + integer nj,njp + + dimension e(norb) + dimension vme(3,norb,norb) + dimension gen_der(3,3,norb,norb) + dimension gd1(3,3,norb,norb) + dimension gd2(3,3,norb,norb) + dimension gd3(3,3,norb,norb) + dimension abc(3,norb,norb) + + real*8 e + complex*16 vme,gen_der,abc,gd1,gd2,gd3,aux1,aux2,aux3 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + abc=0.0d0 + gd1=0.0d0 + gd2=0.0d0 + gd3=0.0d0 + gen_der=0.0d0 + norb_inter_cut=norb + !write(*,*) 'computing berry connection, generalized derivative' + do nn=1,norb + do nnp=1,norb + do nj=1,3 + !if (nn.eq.nnp) then + if (abs(e(nn)-e(nnp)).lt.1.0d-05) then + abc(nj,nn,nnp)=0.0d0 + else + abc(nj,nn,nnp)=-complex(0.0d0,1.0d0)*vme(nj,nn,nnp)/(e(nn)-e(nnp)) + !abc(nj,nnp,nn)=conjg(abc(nj,nn,nnp)) + end if + do njp=1,3 + !if (nn.eq.nnp) then + if (abs(e(nn)-e(nnp)).lt.1.0d-05) then + gd1(nj,njp,nn,nnp)=0.0d0 + gd2(nj,njp,nn,nnp)=0.0d0 + gd3(nj,njp,nn,nnp)=0.0d0 + else + gd1(nj,njp,nn,nnp)=complex(0.0d0,1.0d0)/((e(nn)-e(nnp))**2)* & + (vme(nj,nn,nnp)*(vme(njp,nn,nn)-vme(njp,nnp,nnp))+ & + vme(njp,nn,nnp)*(vme(nj,nn,nn)-vme(nj,nnp,nnp))) + + gd2(nj,njp,nn,nnp)=0.0d0 + gd3(nj,njp,nn,nnp)=0.0d0 + do nnpp=1,norb_inter_cut + !if (nnpp.eq.nn .or. nnpp.eq.nnp) then + if (abs(e(nnpp)-e(nn)).lt.1.0d-05 .or. & + abs(e(nnpp)-e(nnp)).lt.1.0d-05) then + gd2(nj,njp,nn,nnp)=gd2(nj,njp,nn,nnp)+0.0d0 + gd3(nj,njp,nn,nnp)=gd3(nj,njp,nn,nnp)+0.0d0 + else + gd2(nj,njp,nn,nnp)=gd2(nj,njp,nn,nnp)+ & + complex(0.0d0,1.0d0)/(e(nn)-e(nnp))* & + (vme(nj,nn,nnpp)*vme(njp,nnpp,nnp)/(e(nnpp)-e(nnp))) + + gd3(nj,njp,nn,nnp)=gd3(nj,njp,nn,nnp)+ & + complex(0.0d0,1.0d0)/(e(nn)-e(nnp))* & + (-vme(njp,nn,nnpp)*vme(nj,nnpp,nnp)/(e(nn)-e(nnpp))) + end if + end do + !momentums and A and B term + end if + gen_der(nj,njp,nn,nnp)=gd1(nj,njp,nn,nnp)+gd2(nj,njp,nn,nnp) & + +gd3(nj,njp,nn,nnp) + !gen_der(nj,njp,nnp,nn)=conjg(gen_der(nj,njp,nn,nnp)) + end do + end do + end do end do - end do - close(10) - end subroutine write_ome_sp_nonlinear - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! This routine evaluates the alpha,alpha' matrices - subroutine get_vme_kernels_ome(rkx,rky,rkz,norb,skernel,sderkernel, & - hkernel,hderkernel,akernel) - implicit none - - integer norb - integer ialpha - integer ialphap - integer iRp - integer nj - - dimension skernel(norb,norb) - dimension hkernel(norb,norb) - dimension sderkernel(3,norb,norb) - dimension hderkernel(3,norb,norb) - dimension akernel(3,norb,norb) - - dimension hderhop(3,nR,norb,norb) - dimension sderhop(3,nR,norb,norb) - - real(8) Rx,Ry,Rz - real(8) rkx,rky,rkz - - - complex*16 skernel,sderkernel,hkernel,hderkernel,akernel - complex*16 phase,factor - complex*16 hderhop,sderhop,rhop -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - hkernel=0.0d0 - hderkernel=0.0d0 - skernel=0.0d0 - sderkernel=0.0d0 - akernel=0.0d0 - - hderhop=0.0d0 - sderhop=0.0d0 - - do ialpha=1,norb - do ialphap=1,ialpha - !write(*,*) akernel(1,1,1) - !pause - do iRp=1,nR - Rx=dble(nRvec(iRp,1))*R(1,1)+dble(nRvec(iRp,2))*R(2,1) - Ry=dble(nRvec(iRp,1))*R(1,2)+dble(nRvec(iRp,2))*R(2,2) - Rz=0.0d0 - phase=complex(0.0d0,rkx*Rx+rky*Ry+rkz*Rz) - factor=exp(phase) - - hkernel(ialpha,ialphap)=hkernel(ialpha,ialphap)+ & - factor*hhop(iRp,ialpha,ialphap) - skernel(ialpha,ialphap)=skernel(ialpha,ialphap)+ & - factor*shop(iRp,ialpha,ialphap) - - hderhop(1,iRp,ialpha,ialphap)=complex(0.0d0,Rx)*hhop(iRp,ialpha,ialphap) - hderhop(2,iRp,ialpha,ialphap)=complex(0.0d0,Ry)*hhop(iRp,ialpha,ialphap) - hderhop(3,iRp,ialpha,ialphap)=complex(0.0d0,Rz)*hhop(iRp,ialpha,ialphap) - sderhop(1,iRp,ialpha,ialphap)=complex(0.0d0,Rx)*shop(iRp,ialpha,ialphap) - sderhop(2,iRp,ialpha,ialphap)=complex(0.0d0,Ry)*shop(iRp,ialpha,ialphap) - sderhop(3,iRp,ialpha,ialphap)=complex(0.0d0,Rz)*shop(iRp,ialpha,ialphap) - do nj=1,3 - sderkernel(nj,ialpha,ialphap)=sderkernel(nj,ialpha,ialphap)+ & - factor*sderhop(nj,iRp,ialpha,ialphap) - hderkernel(nj,ialpha,ialphap)=hderkernel(nj,ialpha,ialphap)+ & - factor*hderhop(nj,iRp,ialpha,ialphap) - akernel(nj,ialpha,ialphap)=akernel(nj,ialpha,ialphap)+ & - factor*(rhop_c(nj,iRp,ialpha,ialphap)+ & - complex(0.0d0,1.0d0)*sderhop(nj,iRp,ialpha,ialphap)) - end do - !if (ialpha.eq.1 .and. ialphap.eq.1) then - !write(*,*) iRp,akernel(1,1,1),factor,rhop_c(nj,iRp,ialpha,ialphap), & - !sderhop(nj,iRp,ialpha,ialphap) - !end if - end do - !write(*,*) akernel(1,1,1) - !pause - !complex conjugate - do nj=1,3 - hkernel(ialphap,ialpha)=conjg(hkernel(ialpha,ialphap)) - skernel(ialphap,ialpha)=conjg(skernel(ialpha,ialphap)) - sderkernel(nj,ialphap,ialpha)=conjg(sderkernel(nj,ialpha,ialphap)) - hderkernel(nj,ialphap,ialpha)=conjg(hderkernel(nj,ialpha,ialphap)) - akernel(nj,ialphap,ialpha)=conjg(akernel(nj,ialpha,ialphap))+ & - complex(0.0d0,1.0d0)*conjg(sderkernel(nj,ialpha,ialphap)) - end do + end subroutine get_gen_der_sumrule + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_ome_sp_linear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek) + implicit none + integer iflag_norder + integer npointstotal,nband_ex + integer ibz + integer nj,i,j + + dimension ek(npointstotal,nband_ex) + dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) + + real*8 ek + complex*16 vme_ex_band + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + open(10,file='ome_linear_sp_'//trim(material_name)//'.omesp') + write(10,*) iflag_norder + do ibz=1,npointstotal + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(ek(ibz,j),j=1,nband_ex) + do i=1,nband_ex + do j=1,nband_ex + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz), & + (realpart(vme_ex_band(ibz,nj,i,j)),aimag(vme_ex_band(ibz,nj,i,j)), nj=1,3) + end do + end do + end do + close(10) + end subroutine write_ome_sp_linear + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine write_ome_sp_nonlinear(iflag_norder,npointstotal,nband_ex,vme_ex_band,ek,gen_der_ex_band, & + shift_vector_ex_band,berry_eigen_ex_band) + implicit none + integer iflag_norder + integer npointstotal,nband_ex + integer ibz + integer nj,i,j + + dimension ek(npointstotal,nband_ex) + dimension vme_ex_band(npointstotal,3,nband_ex,nband_ex) + dimension berry_eigen_ex_band(npointstotal,3,nband_ex,nband_ex) + dimension gen_der_ex_band(npointstotal,3,3,nband_ex,nband_ex) + dimension shift_vector_ex_band(npointstotal,3,3,nband_ex,nband_ex) + + real*8 ek + real*8 shift_vector_ex_band + complex*16 vme_ex_band + complex*16 berry_eigen_ex_band + complex*16 gen_der_ex_band + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + open(10,file='ome_nonlinear_sp_'//trim(material_name)//'.omesp') + write(10,*) iflag_norder + do ibz=1,npointstotal + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(ek(ibz,j),j=1,nband_ex) + do i=1,nband_ex + do j=1,nband_ex + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz), & + (realpart(vme_ex_band(ibz,nj,i,j)),aimag(vme_ex_band(ibz,nj,i,j)), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(berry_eigen_ex_band(ibz,nj,i,j)), & + aimag(berry_eigen_ex_band(ibz,nj,i,j)), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(shift_vector_ex_band(ibz,1,nj,i,j), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(shift_vector_ex_band(ibz,2,nj,i,j), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(shift_vector_ex_band(ibz,3,nj,i,j), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(gen_der_ex_band(ibz,1,nj,i,j)), & + aimag(gen_der_ex_band(ibz,1,nj,i,j)), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(gen_der_ex_band(ibz,2,nj,i,j)), & + aimag(gen_der_ex_band(ibz,2,nj,i,j)), nj=1,3) + write(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz),(realpart(gen_der_ex_band(ibz,3,nj,i,j)), & + aimag(gen_der_ex_band(ibz,3,nj,i,j)), nj=1,3) + end do + end do end do - end do - end subroutine get_vme_kernels_ome - - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! This routine evaluates the alpha,alpha' matrices - subroutine get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & - hk_ev,e,vme) - implicit none - - integer :: norb - integer :: ialpha - integer :: ialphap - integer :: iRp - integer :: nj - integer :: i,j,ii,jj,nn,nnp - - dimension skernel(norb,norb) - dimension hkernel(norb,norb) - dimension sderkernel(3,norb,norb) - dimension hderkernel(3,norb,norb) - dimension akernel(3,norb,norb) - - dimension vjseudoa(3,norb,norb) - dimension vjseudob(3,norb,norb) - dimension e(norb) - dimension hk_ev(norb,norb) - dimension vme(3,norb,norb) - - real*8 e - complex*16 skernel,sderkernel,hkernel,hderkernel,akernel - complex*16 hk_ev,vjseudoa,vjseudob,vme - complex*16 amu,amup,aux1,factor - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !diagonalization - e=0.0d0 - call diagoz(norb,e,hkernel) - hk_ev(:,:)=hkernel(:,:) - !Multiply the eigenvectors (C_{\alpha=1 n}(k_0), C_{\alpha=1 n}(k_0), ...) - !by a phase phi_n(k_0), this is new eigenvectors are - !exp(-i*phi_n(k_0))*(C_{\alpha=1 n}(k_0), C_{\alpha=1 n}(k_0), ...) - ! - !Note: right now this phase is used to give a locally smooth phase in the BZ - call phase_eigvec_nk(norb,hk_ev) - vme=0.0d0 - vjseudoa=0.0d0 - vjseudob=0.0d0 - - do nn=1,norb - do nnp=1,nn - !momentums and A and B term - do ialpha=1,norb - do ialphap=1,norb - amu=hk_ev(ialpha,nn) - amup=hk_ev(ialphap,nnp) - do nj=1,3 - vjseudoa(nj,nn,nnp)=vjseudoa(nj,nn,nnp)+ & - conjg(amu)*amup*hderkernel(nj,ialpha,ialphap) - vjseudob(nj,nn,nnp)=vjseudob(nj,nn,nnp)+conjg(amu)*amup* & - (e(nn)*akernel(nj,ialpha,ialphap)-e(nnp)*conjg(akernel(nj,ialphap,ialpha)))* & - complex(0.0d0,1.0d0) - end do - end do - end do - !call cpu_time(time2) - !write(*,*) 'k-sampling time',norb,time2-time1,'s' - !pause - do nj=1,3 - vme(nj,nn,nnp)=vjseudoa(nj,nn,nnp)+vjseudob(nj,nn,nnp) - vme(nj,nnp,nn)=conjg(vme(nj,nn,nnp)) - vjseudoa(nj,nnp,nn)=conjg(vjseudoa(nj,nn,nnp)) - vjseudob(nj,nnp,nn)=conjg(vjseudob(nj,nn,nnp)) - end do + close(10) + end subroutine write_ome_sp_nonlinear + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! This routine evaluates the alpha,alpha' matrices + subroutine get_vme_kernels_ome(rkx,rky,rkz,norb,skernel,sderkernel, & + hkernel,hderkernel,akernel) + implicit none + + integer norb + integer ialpha + integer ialphap + integer iRp + integer nj + + dimension skernel(norb,norb) + dimension hkernel(norb,norb) + dimension sderkernel(3,norb,norb) + dimension hderkernel(3,norb,norb) + dimension akernel(3,norb,norb) + + dimension hderhop(3,nR,norb,norb) + dimension sderhop(3,nR,norb,norb) + + real(8) Rx,Ry,Rz + real(8) rkx,rky,rkz + + + complex*16 skernel,sderkernel,hkernel,hderkernel,akernel + complex*16 phase,factor + complex*16 hderhop,sderhop,rhop + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + hkernel=0.0d0 + hderkernel=0.0d0 + skernel=0.0d0 + sderkernel=0.0d0 + akernel=0.0d0 + + hderhop=0.0d0 + sderhop=0.0d0 + + do ialpha=1,norb + do ialphap=1,ialpha + !write(*,*) akernel(1,1,1) + !pause + do iRp=1,nR + Rx=dble(nRvec(iRp,1))*R(1,1)+dble(nRvec(iRp,2))*R(2,1) + Ry=dble(nRvec(iRp,1))*R(1,2)+dble(nRvec(iRp,2))*R(2,2) + Rz=0.0d0 + phase=complex(0.0d0,rkx*Rx+rky*Ry+rkz*Rz) + factor=exp(phase) + + hkernel(ialpha,ialphap)=hkernel(ialpha,ialphap)+ & + factor*hhop(iRp,ialpha,ialphap) + skernel(ialpha,ialphap)=skernel(ialpha,ialphap)+ & + factor*shop(iRp,ialpha,ialphap) + + hderhop(1,iRp,ialpha,ialphap)=complex(0.0d0,Rx)*hhop(iRp,ialpha,ialphap) + hderhop(2,iRp,ialpha,ialphap)=complex(0.0d0,Ry)*hhop(iRp,ialpha,ialphap) + hderhop(3,iRp,ialpha,ialphap)=complex(0.0d0,Rz)*hhop(iRp,ialpha,ialphap) + sderhop(1,iRp,ialpha,ialphap)=complex(0.0d0,Rx)*shop(iRp,ialpha,ialphap) + sderhop(2,iRp,ialpha,ialphap)=complex(0.0d0,Ry)*shop(iRp,ialpha,ialphap) + sderhop(3,iRp,ialpha,ialphap)=complex(0.0d0,Rz)*shop(iRp,ialpha,ialphap) + do nj=1,3 + sderkernel(nj,ialpha,ialphap)=sderkernel(nj,ialpha,ialphap)+ & + factor*sderhop(nj,iRp,ialpha,ialphap) + hderkernel(nj,ialpha,ialphap)=hderkernel(nj,ialpha,ialphap)+ & + factor*hderhop(nj,iRp,ialpha,ialphap) + akernel(nj,ialpha,ialphap)=akernel(nj,ialpha,ialphap)+ & + factor*(rhop_c(nj,iRp,ialpha,ialphap)+ & + complex(0.0d0,1.0d0)*sderhop(nj,iRp,ialpha,ialphap)) + end do + !if (ialpha.eq.1 .and. ialphap.eq.1) then + !write(*,*) iRp,akernel(1,1,1),factor,rhop_c(nj,iRp,ialpha,ialphap), & + !sderhop(nj,iRp,ialpha,ialphap) + !end if + end do + !write(*,*) akernel(1,1,1) + !pause + !complex conjugate + do nj=1,3 + hkernel(ialphap,ialpha)=conjg(hkernel(ialpha,ialphap)) + skernel(ialphap,ialpha)=conjg(skernel(ialpha,ialphap)) + sderkernel(nj,ialphap,ialpha)=conjg(sderkernel(nj,ialpha,ialphap)) + hderkernel(nj,ialphap,ialpha)=conjg(hderkernel(nj,ialpha,ialphap)) + akernel(nj,ialphap,ialpha)=conjg(akernel(nj,ialpha,ialphap))+ & + complex(0.0d0,1.0d0)*conjg(sderkernel(nj,ialpha,ialphap)) + end do + + end do end do - end do - - end subroutine get_vme_eigen_ome - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine phase_eigvec_nk(norb,hk_ev) - implicit none - - integer norb - integer i,j,ii - - dimension hk_ev(norb,norb) - - real*8 :: arg - complex*16 :: aux1,hk_ev,factor -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !phase election: this is done to smooth the gauge - do j=1,norb - aux1=0.0d0 - do i=1,norb - aux1=aux1+hk_ev(i,j) - end do - !argument of the sym - arg=atan2(aimag(aux1),realpart(aux1)) - factor=exp(complex(0.0d0,-arg)) - !write(*,*) 'sum is now:',aux1*factor - do ii=1,norb - hk_ev(ii,j)=hk_ev(ii,j)*factor - !hk_ev(ii,j)=hk_ev(ii,j)*1.0d0 - end do - end do - - end subroutine -end module ome_sp \ No newline at end of file + end subroutine get_vme_kernels_ome + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! This routine evaluates the alpha,alpha' matrices + subroutine get_vme_eigen_ome(norb,skernel,sderkernel,hkernel,hderkernel,akernel, & + hk_ev,e,vme) + implicit none + + integer :: norb + integer :: ialpha + integer :: ialphap + integer :: iRp + integer :: nj + integer :: i,j,ii,jj,nn,nnp + + dimension skernel(norb,norb) + dimension hkernel(norb,norb) + dimension sderkernel(3,norb,norb) + dimension hderkernel(3,norb,norb) + dimension akernel(3,norb,norb) + + dimension vjseudoa(3,norb,norb) + dimension vjseudob(3,norb,norb) + dimension e(norb) + dimension hk_ev(norb,norb) + dimension vme(3,norb,norb) + + real*8 e + complex*16 skernel,sderkernel,hkernel,hderkernel,akernel + complex*16 hk_ev,vjseudoa,vjseudob,vme + complex*16 amu,amup,aux1,factor + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !diagonalization + e=0.0d0 + call diagoz(norb,e,hkernel) + hk_ev(:,:)=hkernel(:,:) + !Multiply the eigenvectors (C_{\alpha=1 n}(k_0), C_{\alpha=1 n}(k_0), ...) + !by a phase phi_n(k_0), this is new eigenvectors are + !exp(-i*phi_n(k_0))*(C_{\alpha=1 n}(k_0), C_{\alpha=1 n}(k_0), ...) + ! + !Note: right now this phase is used to give a locally smooth phase in the BZ + call phase_eigvec_nk(norb,hk_ev) + vme=0.0d0 + vjseudoa=0.0d0 + vjseudob=0.0d0 + + do nn=1,norb + do nnp=1,nn + !momentums and A and B term + do ialpha=1,norb + do ialphap=1,norb + amu=hk_ev(ialpha,nn) + amup=hk_ev(ialphap,nnp) + do nj=1,3 + vjseudoa(nj,nn,nnp)=vjseudoa(nj,nn,nnp)+ & + conjg(amu)*amup*hderkernel(nj,ialpha,ialphap) + vjseudob(nj,nn,nnp)=vjseudob(nj,nn,nnp)+conjg(amu)*amup* & + (e(nn)*akernel(nj,ialpha,ialphap)-e(nnp)*conjg(akernel(nj,ialphap,ialpha)))* & + complex(0.0d0,1.0d0) + end do + end do + end do + !call cpu_time(time2) + !write(*,*) 'k-sampling time',norb,time2-time1,'s' + !pause + do nj=1,3 + vme(nj,nn,nnp)=vjseudoa(nj,nn,nnp)+vjseudob(nj,nn,nnp) + vme(nj,nnp,nn)=conjg(vme(nj,nn,nnp)) + vjseudoa(nj,nnp,nn)=conjg(vjseudoa(nj,nn,nnp)) + vjseudob(nj,nnp,nn)=conjg(vjseudob(nj,nn,nnp)) + end do + end do + end do + + end subroutine get_vme_eigen_ome + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine phase_eigvec_nk(norb,hk_ev) + implicit none + + integer norb + integer i,j,ii + + dimension hk_ev(norb,norb) + + real*8 :: arg + complex*16 :: aux1,hk_ev,factor + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !phase election: this is done to smooth the gauge + do j=1,norb + aux1=0.0d0 + do i=1,norb + aux1=aux1+hk_ev(i,j) + end do + !argument of the sym + arg=atan2(aimag(aux1),realpart(aux1)) + factor=exp(complex(0.0d0,-arg)) + !write(*,*) 'sum is now:',aux1*factor + do ii=1,norb + hk_ev(ii,j)=hk_ev(ii,j)*factor + !hk_ev(ii,j)=hk_ev(ii,j)*1.0d0 + end do + end do + + end subroutine +end module ome_sp diff --git a/src/parser_input_file.f90 b/src/parser_input_file.f90 index 37eecab..9d3906f 100644 --- a/src/parser_input_file.f90 +++ b/src/parser_input_file.f90 @@ -17,6 +17,7 @@ module parser_input_file public :: get_input_file public :: nband_index public :: norb_ex_cut + public :: broadening_type_text public :: read_line_numbers_int !subroutine character(len=1000) :: material_name_in @@ -24,6 +25,7 @@ module parser_input_file character(len=100) :: iflag_xatu_text character(len=100) :: iflag_ome_sp_text character(len=100) :: iflag_ome_ex_text + character(len=100) :: broadening_type_text character(len=1000) :: xatu_eigval_filepath_in character(len=1000) :: xatu_states_filepath_in character(len=100) :: response_text @@ -46,82 +48,155 @@ module parser_input_file dimension :: nband_index_aux(100) !auxiliary array to save band indeces allocatable :: nband_index(:) - contains + contains + function to_lower(str) result(lower_str) + implicit none + character(len=*), intent(in) :: str + character(len=len(str)) :: lower_str + integer :: i, ic + + do i = 1, len(str) + ic = iachar(str(i:i)) + if (ic >= iachar('A') .and. ic <= iachar('Z')) then + lower_str(i:i) = achar(ic + 32) + else + lower_str(i:i) = str(i:i) + end if + end do + end function to_lower + subroutine get_input_file() implicit none - allocatable :: narray(:) - integer :: i,narray,num_values + integer, allocatable :: narray(:) + integer :: i, num_values, ios + character(len=1000) :: line + character(len=100) :: param_name + logical :: ndim_found, material_found, xatu_found, bandlist_found + logical :: ncells_found, nfermi_found, ome_sp_found, ome_ex_found + logical :: response_found, energy_found, exciton_found !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! write(*,*) '1. Entering parser_input_file' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Initialize flags + ndim_found = .false. + material_found = .false. + xatu_found = .false. + bandlist_found = .false. + ncells_found = .false. + nfermi_found = .false. + ome_sp_found = .false. + ome_ex_found = .false. + response_found = .false. + energy_found = .false. + exciton_found = .false. + ! default broadening + broadening_type_text = 'gaussian' + call get_command_argument(1,filename_input) - !write(*,*) trim(filename_input) - !pause - !open(10,file='./bin/'//trim(filename_input)//'') open(10,file=adjustl(filename_input)) - read(10,*) !# Periodic dimensions - read(10,*) ndim - read(10,*) !# Wannier90_filename - read(10,'(A)') material_name_in - read(10,*) !# Xatu_interface - read(10,*) iflag_xatu_text - if (iflag_xatu_text == 'true') then - iflag_xatu= .true. - !get paths for .eigval and .states files - read(10,'(A)') xatu_eigval_filepath_in - read(10,'(A)') xatu_states_filepath_in - npointstotal_sq=0 - read(10,*) !# Exciton_cutoff - read(10,*) norb_ex_cut - read(10,*) !# Nfermi - read(10,*) nf - read(10,*) !# OME_sp - read(10,*) iflag_ome_sp_text - read(10,*) !# OME_ex - read(10,*) iflag_ome_ex_text - read(10,*) !# Response - read(10,*) response_text - read(10,*) !# Energy_variables - read(10,*) e1,e2,eta,nw - close(10) - else if (iflag_xatu_text == 'false') then - iflag_xatu= .false. - read(10,*) !# Bandlist - !reads a line of numbers and stores them in an allocatable array 'narray' - call read_line_numbers_int(narray,num_values) - !write(*,*) (narray(j),j=1,num_values) - read(10,*) !# Ncells - read(10,*) npointstotal_sq - read(10,*) !# Nfermi - read(10,*) nf - read(10,*) !# OME_sp - read(10,*) iflag_ome_sp_text - read(10,*) !# Response - read(10,*) response_text - read(10,*) !# Energy_variables - read(10,*) e1,e2,eta,nw - close(10) - !fill the nband_index array after knowing the value of num_values - !it will be used in parser_optics_xatu_dim.f90 - allocate (nband_index(num_values)) - nband_index(:)=narray(:) - else - write(*,*) 'Error: Invalid value in the input file. Expected "true" or "false".' - stop + + ! Read file sequentially and process parameters based on their labels + do + read(10,'(A)',iostat=ios) line + if (ios /= 0) exit ! End of file + + line = adjustl(line) + + ! Check if this is a comment line (parameter label) + if (line(1:1) == '#') then + ! Extract parameter name and read corresponding value + param_name = adjustl(line(3:)) ! Remove "# " prefix + + if (index(param_name, 'Periodic dimensions') > 0) then + read(10,*) ndim + ndim_found = .true. + + else if (index(param_name, 'Wannier90_filename') > 0) then + read(10,'(A)') material_name_in + material_found = .true. + + else if (index(param_name, 'Xatu_interface') > 0) then + read(10,*) iflag_xatu_text + xatu_found = .true. + + if (iflag_xatu_text == 'true') then + iflag_xatu = .true. + ! Read the eigval and states file paths that follow + read(10,'(A)') xatu_eigval_filepath_in + read(10,'(A)') xatu_states_filepath_in + else if (iflag_xatu_text == 'false') then + iflag_xatu = .false. + else + write(*,*) 'Error: Invalid value in Xatu_interface. Expected "true" or "false".' + stop + end if + + else if (index(param_name, 'Exciton_cutoff') > 0) then + read(10,*) norb_ex_cut + exciton_found = .true. + + else if (index(param_name, 'Bandlist') > 0) then + call read_line_numbers_int(narray, num_values) + bandlist_found = .true. + + else if (index(param_name, 'Ncells') > 0) then + read(10,*) npointstotal_sq + ncells_found = .true. + + else if (index(param_name, 'Nfermi') > 0) then + read(10,*) nf + nfermi_found = .true. + + else if (index(param_name, 'OME_sp') > 0 .or. index(param_name, 'OME_SP') > 0) then + read(10,*) iflag_ome_sp_text + ome_sp_found = .true. + + else if (index(param_name, 'OME_ex') > 0 .or. index(param_name, 'OME_EX') > 0) then + read(10,*) iflag_ome_ex_text + ome_ex_found = .true. + + else if (index(param_name, 'Response') > 0) then + read(10,*) response_text + response_found = .true. + + else if (index(param_name, 'Energy_variables') > 0) then + read(10,*) e1, e2, eta, nw + energy_found = .true. + else if (index(param_name, 'Broadening_type') > 0 .or. index(param_name,'Broadening')>0) then + read(10,'(A)') broadening_type_text + broadening_type_text = adjustl(broadening_type_text) + broadening_type_text = to_lower(broadening_type_text) + + end if + end if + end do + + close(10) + + ! Handle bandlist case: allocate nband_index if bandlist was found + if (bandlist_found) then + allocate(nband_index(num_values)) + nband_index(:) = narray(:) end if - !declare flags from text strings + ! Set npointstotal_sq to 0 if using xatu interface + if (iflag_xatu) then + npointstotal_sq = 0 + end if + + ! Declare flags from text strings if (iflag_ome_sp_text == 'true') then - iflag_ome_sp= .true. + iflag_ome_sp = .true. else - iflag_ome_sp= .false. + iflag_ome_sp = .false. end if + if (iflag_ome_ex_text == 'true') then - iflag_ome_ex= .true. + iflag_ome_ex = .true. else - iflag_ome_ex= .false. + iflag_ome_ex = .false. end if write(*,*) ' Input file has been read' diff --git a/src/parser_optics_xatu_dim.f90 b/src/parser_optics_xatu_dim.f90 index 4238cdc..80ada8a 100644 --- a/src/parser_optics_xatu_dim.f90 +++ b/src/parser_optics_xatu_dim.f90 @@ -1,388 +1,422 @@ -module parser_optics_xatu_dim - use constants_math - use parser_wannier90_tb, & - only:material_name,R !variables - use parser_input_file, & - only:xatu_eigval_filepath_in,xatu_states_filepath_in, & !filepaths - ndim,npointstotal_sq, & !variables - iflag_xatu,nf,nband_index,norb_ex_cut, & - read_line_numbers_int !subroutine - implicit none - - integer :: nv_ex,nc_ex - integer :: npointstotal - integer :: norb_ex - integer :: norb_ex_band - integer :: nband_ex - integer :: naux - integer :: j - - real(8) G,vcell - real(8) rkxvector,rkyvector,rkzvector - real(8) auxr1 - real(8) e_ex - complex*16 fk_ex - - dimension G(3,3) - - allocatable rkxvector(:) - allocatable rkyvector(:) - allocatable rkzvector(:) - allocatable fk_ex(:,:) - allocatable e_ex(:) - allocatable auxr1(:) - - contains - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !TO BE ADAPTED FOR 3D - ! Here we define some BZ variables, either by reading the output - ! of Xatu or not - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_optics_xatu_dim() - implicit none - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - write(*,*) '3. Entering parser_optics_xatu_dim' - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !get reciprocal lattice vectors (to be adapted for 3D!) - call get_reciprocal_vectors() - !write(*,*) (nband_index(j), j=1,size(nband_index,dim=1)) - !pause - - !Reminder:norb_ex_cut is given as input - !get band and grid dimensions: from XATU-output or opticx-input - if (iflag_xatu .eqv. .true.) then - call get_exciton_dim() - !nband_ex=nc_ex+nv_ex - else - !number of total k-points - if (ndim==1) npointstotal=npointstotal_sq - if (ndim==2) npointstotal=npointstotal_sq**2 - if (ndim==3) npointstotal=npointstotal_sq**3 - norb_ex_band=nv_ex*nc_ex - norb_ex=norb_ex_band*npointstotal - end if - !calculate nv_ex and nc_ex from the array of bands - nband_ex=size(nband_index,dim=1) - nv_ex=0 - do j=1,nband_ex - if (nband_index(j).le.0) nv_ex=nv_ex+1 - end do - nc_ex=nband_ex-nv_ex - - !change syntax for band counting - !XATU: ...-1 0 1 2... to explicit band count - !opticx: ...nf-1,nf,nf+1... - nband_index(:)=nband_index(:)+nf - !write(*,*) (nband_index(j), j=1,nband_ex) - !write(*,*) - !pause - - !allocate grid and exciton arrays - allocate (rkxvector(npointstotal)) - allocate (rkyvector(npointstotal)) - allocate (rkzvector(npointstotal)) - allocate (e_ex(norb_ex_cut)) - allocate (fk_ex(norb_ex,norb_ex_cut)) - - !fill exciton and other arrays: from XATU-output or opticx-input - if (iflag_xatu .eqv. .true.) then - call get_exciton_data() !get grid and exciton wavefunctions - write(*,*) ' Exciton data has been read from XATU output' - else - call get_grid() - !get frid andexciton variables set to zero if XATU interface is not requested - fk_ex=0.0d0 - e_ex=0.0d0 - end if - write(*,*) " Grid and band parameters have been set" - - !write(*,*) rkxvector(1),rkyvector(1),rkzvector(1) - !write(*,*) G(1,1),G(1,2),G(1,3) - !write(*,*) G(2,1),G(2,2),G(2,3) - !pause - end subroutine get_optics_xatu_dim -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! This subroutine prints a part of the exciton wavefunction or -! the total exciton probability density - subroutine print_exciton_wf(isum,iv,ic,nn) - implicit none - integer isum,iv,ic,nn - integer iv_s,ic_s - integer iright,i_ex_nn - integer ibz - real*8 prob_k -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - iright=0 - open(10,file='exciton_wf.dat') - do ibz=1,npointstotal - !write(*,*) ibz,npointstotal,fk_ex(ibz,1),fk_ex(ibz,1),fk_ex(ibz,3) - !pause - !isum=0: print the total probability density - !isum=1: print the probability density of the exciton with index (iv,ic) - if (isum.eq.1) then - prob_k=0.0d0 - do iv_s=1,nv_ex - do ic_s=1,nc_ex - call get_ex_index_first(nf,nv_ex,nc_ex,iright,ibz,i_ex_nn,ic_s,iv_s) - prob_k=prob_k+abs(fk_ex(i_ex_nn,nn))**2 - end do - end do - write(10,*) rkxvector(ibz),rkyvector(ibz),prob_k - else - call get_ex_index_first(nf,nv_ex,nc_ex,iright,ibz,i_ex_nn,ic,iv) - write(10,*) rkxvector(ibz),rkyvector(ibz),abs(fk_ex(i_ex_nn,nn)) !,nf,nv_ex,nc_ex,iright,ibz,i_ex_nn,ic,iv - end if - end do - - close(10) - end subroutine print_exciton_wf -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !this subroutine has not been updated in 2025 - subroutine get_ex_index_first(nf,nv_ex,nc_ex,iright,ibz,i_ex,ic,iv) - implicit none - integer nf,nv_ex,nc_ex,ibz,i_ex,ic,iv - integer iright -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !get band indeces (respect the fermi level) from A_cv index - if (iright.eq.1) then - iv=(nf-nv_ex)+i_ex-int((i_ex-1)/nv_ex)*nv_ex-nf - ic=(nf+1)+int((i_ex-1)/nv_ex)-nf - end if - !get A_cv index from band indeces - if (iright.ne.1) then - i_ex=nc_ex*nv_ex*(ibz-1)+nv_ex*(ic-1)+iv - end if - end subroutine get_ex_index_first -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_exciton_dim() - implicit none - !allocatable :: file2open(:) - !character :: file2open - dimension nband_index_aux1(100) - dimension nband_index_aux2(100) - - - integer :: nband_ex_aux !number of valence bands - integer :: nband_index_aux1,nband_index_aux2 !auxiliary - integer :: nband_ex_aux1,nband_ex_aux2 !auxiliary - integer :: iexit - integer :: i,j,naux,npointstotal_sq - - real(8) aux1 - character(len=:), allocatable :: file2open -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !UPDATE: 2025-07-09. It is no longer necessary to modify .eigval file to read bandlist - !This part gets 'nband_ex' and 'nband_index(nband_ex)' - !Note: nband_ex_aux is used as nband_ex, which will be defined later (historic reasons) - ! but we should have nband_ex=nband_ex_aux - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !save number of valence bands - file2open=trim(xatu_states_filepath_in) - open(10,file=file2open) - read(10,*) - do i=1,5000 !we will never use more than 50 bands... - read(10,*) aux1,aux1,aux1,nband_index_aux1(i) - if (i.gt.1) then - do j=1,i-1 - if (nband_index_aux1(i).eq.nband_index_aux1(j)) then - nband_ex_aux1=i-1 !save number of valence bands - goto 128 - end if - end do - end if - end do -128 continue - close(10) - !write(*,*) ' Number of valence bands:',nband_ex1 - !save number of conduction bands - open(10,file=file2open) - read(10,*) - do i=1,5000 !we will never use more than 50 bands... - read(10,*) aux1,aux1,aux1,naux,nband_index_aux2(i) - !skip similar values - if (nband_ex_aux1.gt.1) then - do j=1,nband_ex_aux1-1 - read(10,*) - end do - end if - if (i.gt.1) then - do j=1,i-1 - if (nband_index_aux2(i).eq.nband_index_aux2(j)) then - nband_ex_aux2=(i-1) !save number of valence bands - goto 129 - end if - end do - end if - end do -129 continue - close(10) - !write(*,*) ' Number of total bands:',nband_ex_aux1,nband_ex_aux2 - nband_ex_aux=nband_ex_aux1+nband_ex_aux2 - allocate(nband_index(nband_ex_aux)) - do i=1,nband_ex_aux1 - nband_index(i)=nband_index_aux1(i)-nf+1 - end do - do i=nband_ex_aux1+1,nband_ex_aux - nband_index(i)=nband_index_aux2(i-nband_ex_aux1)-nf+1 - end do - !old, previous lines replace this call - !reads a line of numbers and stores them in an allocatable array 'narray' - !call read_line_numbers_int(nband_index,nband_ex) - !write(*,*) nband_ex_aux - !write(*,*) (nband_index(j), j=1,nband_ex_aux) - !pause - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !get nk - file2open=trim(xatu_eigval_filepath_in) - open(10,file=file2open) - read(10,*) npointstotal_sq - read(10,*) naux - close(10) - - !get N_BSE=nv_ex*nc_ex*nk**2 variables - if (ndim==1) npointstotal=npointstotal_sq - if (ndim==2) npointstotal=npointstotal_sq**2 - if (ndim==3) npointstotal=npointstotal_sq**3 - norb_ex_band=int(naux/npointstotal) - norb_ex=norb_ex_band*npointstotal - - end subroutine get_exciton_dim -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !This routine needs to be adapted for 3D - subroutine get_reciprocal_vectors() - implicit none - real(8) cx,cy,cz -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - G=0.0d0 - G(1,1)=2.0d0*pi*(-R(1,1)*R(2,2)+R(1,2)*R(2,1))**(-1.0d0) & - *(-R(2,2)) - G(1,2)=2.0d0*pi*(-R(1,1)*R(2,2)+R(1,2)*R(2,1))**(-1.0d0) & - *R(2,1) - G(2,1)=2.0d0*pi*(-R(2,1)*R(1,2)+R(2,2)*R(1,1))**(-1.0d0) & - *(-R(1,2)) - G(2,2)=2.0d0*pi*(-R(2,1)*R(1,2)+R(2,2)*R(1,1))**(-1.0d0) & - *R(1,1) - - call crossproduct(R(1,1),R(1,2),0.0d0,R(2,1), & - R(2,2),0.0d0,cx,cy,cz) - vcell=sqrt(cx**2+cy**2+cz**2) - - end subroutine -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_exciton_data() - implicit none - integer j,nkaka, max_ex, ios - integer ib,ibz,ibz_sum,jind - real(8) auxr1 - - dimension auxr1(2*norb_ex) - character(len=:), allocatable :: file2open -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Read eigenvalues from file - file2open = trim(xatu_eigval_filepath_in) - - open(unit=10, file=file2open, status='old', action='read', iostat=ios) - if (ios /= 0) then - print *, 'Error opening eigenvalue file: ', trim(file2open) - stop - end if - - ! skipping headers - read(10,*) ! npointstotal_sq already read - read(10,*) ! naux already read - read(10,*) max_ex - - if (norb_ex_cut > max_ex) then - norb_ex_cut = max_ex - print *, "!!! Exciton cutoff is bigger than available eigvalues in .eigval !!!" - print *, " --> Setting cutoff to maximum number of available excitons: ", max_ex - end if - - ! read norb_ex_cut values, one per line - do j = 1, norb_ex_cut - read(10, *, iostat=ios) e_ex(j) - if (ios /= 0) then - print *, 'Error reading eigenvalue at line', j+1 - stop - end if - end do - - close(10) - - - file2open=trim(xatu_states_filepath_in) - open(10,file=file2open) - read(10,*) - - !reading k-mesh - do ibz=1,npointstotal - read(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz) - do ib=1,norb_ex_band-1 - read(10,*) - end do - end do - - !reading exciton-wf - ibz_sum=0 !counter to display percentage - write(*,*) ' Reading exciton wavefunctions...' - do ibz=1,norb_ex_cut - ibz_sum=ibz_sum+1 - !display percentaje - if (abs(dble(ibz)/dble(norb_ex_cut))*100.0d0-100.0d0 .lt. 5.0d0) then - call percentage_index(ibz_sum,norb_ex_cut,nkaka) - !write(*,*) ' Reading exciton wf:',int(abs(dble(ibz)/dble(norb_ex_cut))*100.0d0),'%' - !write(*,*) ' Reading exciton wf:',abs(dble(ibz)/dble(norb_ex_cut))*100.0d0,'%' - end if - read(10,*) (auxr1(j),j=1,2*norb_ex) - !write(*,*) (auxr1(j),j=1,2*norb_ex) - !pause - do j=1,norb_ex - jind=2*j-1 - fk_ex(j,ibz)=complex(auxr1(jind),auxr1(jind+1)) - end do - !pause - end do - - close(10) - - !Please I like to work in atomic units! - e_ex=e_ex/27.211385d0 - rkxvector=rkxvector*0.52917721067121d0 - rkyvector=rkyvector*0.52917721067121d0 - rkzvector=rkzvector*0.52917721067121d0 - - end subroutine get_exciton_data - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !TO BE ADAPTED FOR 3D - subroutine get_grid() - implicit none - - integer :: k,i1,i2 - real(8) :: step,r1,r2 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - step=1.0d0/dble(npointstotal_sq-1) - - !atomic units provided that G are in atomic units - k=1 - do i1=1,npointstotal_sq - r1=-0.5d0+dble(i1-1)*step - do i2=1,npointstotal_sq - r2=-0.5d0+dble(i2-1)*step - rkxvector(k)=r1*G(1,1)+r2*G(2,1) - rkyvector(k)=r1*G(1,2)+r2*G(2,2) - rkzvector(k)=0.0d0 - k=k+1 - end do - end do - - end subroutine get_grid - - -end module parser_optics_xatu_dim - +module parser_optics_xatu_dim + use constants_math + use parser_wannier90_tb, & + only:material_name,R !variables + use parser_input_file, & + only:xatu_eigval_filepath_in,xatu_states_filepath_in, & !filepaths + ndim,npointstotal_sq, & !variables + iflag_xatu,nf,nband_index,norb_ex_cut, & + read_line_numbers_int !subroutine + implicit none + + integer :: nv_ex,nc_ex + integer :: npointstotal + integer :: norb_ex + integer :: norb_ex_band + integer :: nband_ex + integer :: naux + integer :: j + + real(8) G,vcell + real(8) rkxvector,rkyvector,rkzvector + real(8) auxr1 + real(8) e_ex + complex*16 fk_ex + + dimension G(3,3) + + allocatable rkxvector(:) + allocatable rkyvector(:) + allocatable rkzvector(:) + allocatable fk_ex(:,:) + allocatable e_ex(:) + allocatable auxr1(:) + + contains + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !TO BE ADAPTED FOR 3D + ! Here we define some BZ variables, either by reading the output + ! of Xatu or not + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_optics_xatu_dim() + implicit none + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + write(*,*) '3. Entering parser_optics_xatu_dim' + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !get reciprocal lattice vectors (to be adapted for 3D!) + call get_reciprocal_vectors() + !write(*,*) (nband_index(j), j=1,size(nband_index,dim=1)) + !pause + + !Reminder:norb_ex_cut is given as input + !get band and grid dimensions: from XATU-output or opticx-input + if (iflag_xatu .eqv. .true.) then + call get_exciton_dim() + !nband_ex=nc_ex+nv_ex + else + !number of total k-points + if (ndim==1) npointstotal=npointstotal_sq + if (ndim==2) npointstotal=npointstotal_sq**2 + if (ndim==3) npointstotal=npointstotal_sq**3 + norb_ex_band=nv_ex*nc_ex + norb_ex=norb_ex_band*npointstotal + end if + !calculate nv_ex and nc_ex from the array of bands + nband_ex=size(nband_index,dim=1) + nv_ex=0 + do j=1,nband_ex + if (nband_index(j).le.0) nv_ex=nv_ex+1 + end do + nc_ex=nband_ex-nv_ex + + !change syntax for band counting + !XATU: ...-1 0 1 2... to explicit band count + !opticx: ...nf-1,nf,nf+1... + nband_index(:)=nband_index(:)+nf + !write(*,*) (nband_index(j), j=1,nband_ex) + !write(*,*) + !pause + + !allocate grid and exciton arrays + allocate (rkxvector(npointstotal)) + allocate (rkyvector(npointstotal)) + allocate (rkzvector(npointstotal)) + allocate (e_ex(norb_ex_cut)) + allocate (fk_ex(norb_ex,norb_ex_cut)) + + !fill exciton and other arrays: from XATU-output or opticx-input + if (iflag_xatu .eqv. .true.) then + call get_exciton_data() !get grid and exciton wavefunctions + write(*,*) ' Exciton data has been read from XATU output' + else + call get_grid() + !get frid andexciton variables set to zero if XATU interface is not requested + fk_ex=0.0d0 + e_ex=0.0d0 + end if + write(*,*) " Grid and band parameters have been set" + + !write(*,*) rkxvector(1),rkyvector(1),rkzvector(1) + !write(*,*) G(1,1),G(1,2),G(1,3) + !write(*,*) G(2,1),G(2,2),G(2,3) + !pause + end subroutine get_optics_xatu_dim +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! This subroutine prints a part of the exciton wavefunction or +! the total exciton probability density + subroutine print_exciton_wf(isum,iv,ic,nn) + implicit none + integer isum,iv,ic,nn + integer iv_s,ic_s + integer iright,i_ex_nn + integer ibz + real*8 prob_k +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + iright=0 + open(10,file='exciton_wf.dat') + do ibz=1,npointstotal + !write(*,*) ibz,npointstotal,fk_ex(ibz,1),fk_ex(ibz,1),fk_ex(ibz,3) + !pause + !isum=0: print the total probability density + !isum=1: print the probability density of the exciton with index (iv,ic) + if (isum.eq.1) then + prob_k=0.0d0 + do iv_s=1,nv_ex + do ic_s=1,nc_ex + call get_ex_index_first(nf,nv_ex,nc_ex,iright,ibz,i_ex_nn,ic_s,iv_s) + prob_k=prob_k+abs(fk_ex(i_ex_nn,nn))**2 + end do + end do + write(10,*) rkxvector(ibz),rkyvector(ibz),prob_k + else + call get_ex_index_first(nf,nv_ex,nc_ex,iright,ibz,i_ex_nn,ic,iv) + write(10,*) rkxvector(ibz),rkyvector(ibz),abs(fk_ex(i_ex_nn,nn)) !,nf,nv_ex,nc_ex,iright,ibz,i_ex_nn,ic,iv + end if + end do + + close(10) + end subroutine print_exciton_wf +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !this subroutine has not been updated in 2025 + subroutine get_ex_index_first(nf,nv_ex,nc_ex,iright,ibz,i_ex,ic,iv) + implicit none + integer nf,nv_ex,nc_ex,ibz,i_ex,ic,iv + integer iright +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !get band indeces (respect the fermi level) from A_cv index + if (iright.eq.1) then + iv=(nf-nv_ex)+i_ex-int((i_ex-1)/nv_ex)*nv_ex-nf + ic=(nf+1)+int((i_ex-1)/nv_ex)-nf + end if + !get A_cv index from band indeces + if (iright.ne.1) then + i_ex=nc_ex*nv_ex*(ibz-1)+nv_ex*(ic-1)+iv + end if + end subroutine get_ex_index_first +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_exciton_dim() + implicit none + !allocatable :: file2open(:) + !character :: file2open + dimension nband_index_aux1(100) + dimension nband_index_aux2(100) + + + integer :: nband_ex_aux !number of valence bands + integer :: nband_index_aux1,nband_index_aux2 !auxiliary + integer :: nband_ex_aux1,nband_ex_aux2 !auxiliary + integer :: iexit + integer :: i,j,naux,npointstotal_sq + integer :: hdr1, hdr2, ios + + real(8) aux1 + character(len=:), allocatable :: file2open +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !UPDATE: 2025-07-09. It is no longer necessary to modify .eigval file to read bandlist + !This part gets 'nband_ex' and 'nband_index(nband_ex)' + !Note: nband_ex_aux is used as nband_ex, which will be defined later (historic reasons) + ! but we should have nband_ex=nband_ex_aux + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !save number of valence bands + file2open=trim(xatu_states_filepath_in) + open(10,file=file2open) + read(10,*) + do i=1,50 !we will never use more than 50 bands... + read(10,*) aux1,aux1,aux1,nband_index_aux1(i) + if (i.gt.1) then + do j=1,i-1 + if (nband_index_aux1(i).eq.nband_index_aux1(j)) then + nband_ex_aux1=i-1 !save number of valence bands + goto 128 + end if + end do + end if + end do +128 continue + close(10) + !write(*,*) ' Number of valence bands:',nband_ex1 + !save number of conduction bands + open(10,file=file2open) + read(10,*) + do i=1,50 !we will never use more than 50 bands... + read(10,*) aux1,aux1,aux1,naux,nband_index_aux2(i) + !skip similar values + if (nband_ex_aux1.gt.1) then + do j=1,nband_ex_aux1-1 + read(10,*) + end do + end if + if (i.gt.1) then + do j=1,i-1 + if (nband_index_aux2(i).eq.nband_index_aux2(j)) then + nband_ex_aux2=(i-1) !save number of valence bands + goto 129 + end if + end do + end if + end do +129 continue + close(10) + !write(*,*) ' Number of total bands:',nband_ex_aux1,nband_ex_aux2 + nband_ex_aux=nband_ex_aux1+nband_ex_aux2 + allocate(nband_index(nband_ex_aux)) + do i=1,nband_ex_aux1 + nband_index(i)=nband_index_aux1(i)-nf+1 + end do + do i=nband_ex_aux1+1,nband_ex_aux + nband_index(i)=nband_index_aux2(i-nband_ex_aux1)-nf+1 + end do + !old, previous lines replace this call + !reads a line of numbers and stores them in an allocatable array 'narray' + !call read_line_numbers_int(nband_index,nband_ex) + !write(*,*) nband_ex_aux + !write(*,*) (nband_index(j), j=1,nband_ex_aux) + !pause + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !get nk + file2open=trim(xatu_eigval_filepath_in) + open(10,file=file2open) + ! Try old format: first line npointstotal_sq, second line naux + read(10,*,iostat=ios) hdr1 + if (ios == 0) then + read(10,*,iostat=ios) hdr2 + if (ios == 0) then + npointstotal_sq = hdr1 + naux = hdr2 + else + ! New xatu Result::writeEigenvalues format: first line = exciton basis size (naux) + naux = hdr1 + ! compute npointstotal from nv_ex*nc_ex (number of pair combinations per k) + if (nv_ex*nc_ex > 0) then + npointstotal = naux / (nv_ex*nc_ex) + else + npointstotal = 0 + end if + ! derive npointstotal_sq based on ndim + if (ndim == 1) then + npointstotal_sq = npointstotal + else if (ndim == 2) then + npointstotal_sq = int(sqrt(dble(npointstotal)) + 0.5d0) + else if (ndim == 3) then + npointstotal_sq = int((dble(npointstotal))**(1.0d0/3.0d0) + 0.5d0) + end if + end if + else + ! Fallback: try older format read without iostat + rewind(10) + read(10,*) npointstotal_sq + read(10,*) naux + end if + close(10) + + !get N_BSE=nv_ex*nc_ex*nk**2 variables + if (npointstotal == 0) then + if (ndim==1) npointstotal=npointstotal_sq + if (ndim==2) npointstotal=npointstotal_sq**2 + if (ndim==3) npointstotal=npointstotal_sq**3 + end if + if (npointstotal > 0) then + norb_ex_band = int(naux / npointstotal) + else + norb_ex_band = 0 + end if + norb_ex = norb_ex_band * npointstotal + + end subroutine get_exciton_dim +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !This routine needs to be adapted for 3D + subroutine get_reciprocal_vectors() + implicit none + real(8) cx,cy,cz +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + G=0.0d0 + G(1,1)=2.0d0*pi*(-R(1,1)*R(2,2)+R(1,2)*R(2,1))**(-1.0d0) & + *(-R(2,2)) + G(1,2)=2.0d0*pi*(-R(1,1)*R(2,2)+R(1,2)*R(2,1))**(-1.0d0) & + *R(2,1) + G(2,1)=2.0d0*pi*(-R(2,1)*R(1,2)+R(2,2)*R(1,1))**(-1.0d0) & + *(-R(1,2)) + G(2,2)=2.0d0*pi*(-R(2,1)*R(1,2)+R(2,2)*R(1,1))**(-1.0d0) & + *R(1,1) + + call crossproduct(R(1,1),R(1,2),0.0d0,R(2,1), & + R(2,2),0.0d0,cx,cy,cz) + vcell=sqrt(cx**2+cy**2+cz**2) + + end subroutine +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine get_exciton_data() + implicit none + integer j,nkaka + integer ib,ibz,ibz_sum,jind + real(8) auxr1 + + dimension auxr1(2*norb_ex) + character(len=:), allocatable :: file2open + integer :: header1, ios +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !get energies + file2open=trim(xatu_eigval_filepath_in) + open(10,file=file2open) + ! Try to read the format produced by xatu Result::writeEigenvalues: + ! Line 1: exciton basis size (ignored here) + ! Line 2: number of eigenvalues (nkaka) + ! Next lines: one eigenvalue per line + read(10,*,iostat=ios) header1 + read(10,*,iostat=ios) + if (ios == 0) then + read(10,*,iostat=ios) nkaka + if (ios == 0) then + do j=1,norb_ex_cut + read(10,*,iostat=ios) e_ex(j) + if (ios /= 0) exit + end do + else + ! fallback to older single-line format + rewind(10) + read(10,*) + read(10,*) nkaka,(e_ex(j), j=1,norb_ex_cut) + end if + else + ! fallback to older single-line format + rewind(10) + read(10,*) + read(10,*) nkaka,(e_ex(j), j=1,norb_ex_cut) + end if + close(10) + + file2open=trim(xatu_states_filepath_in) + open(10,file=file2open) + read(10,*) + + !reading k-mesh + do ibz=1,npointstotal + read(10,*) rkxvector(ibz),rkyvector(ibz),rkzvector(ibz) + do ib=1,norb_ex_band-1 + read(10,*) + end do + end do + + !reading exciton-wf + ibz_sum=0 !counter to display percentage + write(*,*) ' Reading exciton wavefunctions...' + do ibz=1,norb_ex_cut + ibz_sum=ibz_sum+1 + !display percentaje + if (abs(dble(ibz)/dble(norb_ex_cut))*100.0d0-100.0d0 .lt. 5.0d0) then + call percentage_index(ibz_sum,norb_ex_cut,nkaka) + !write(*,*) ' Reading exciton wf:',int(abs(dble(ibz)/dble(norb_ex_cut))*100.0d0),'%' + !write(*,*) ' Reading exciton wf:',abs(dble(ibz)/dble(norb_ex_cut))*100.0d0,'%' + end if + read(10,*) (auxr1(j),j=1,2*norb_ex) + !write(*,*) (auxr1(j),j=1,2*norb_ex) + !pause + do j=1,norb_ex + jind=2*j-1 + fk_ex(j,ibz)=complex(auxr1(jind),auxr1(jind+1)) + end do + !pause + end do + + close(10) + + !Please I like to work in atomic units! + e_ex=e_ex/27.211385d0 + rkxvector=rkxvector*0.52917721067121d0 + rkyvector=rkyvector*0.52917721067121d0 + rkzvector=rkzvector*0.52917721067121d0 + + end subroutine get_exciton_data + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !TO BE ADAPTED FOR 3D + subroutine get_grid() + implicit none + + integer :: k,i1,i2 + real(8) :: step,r1,r2 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + step=1.0d0/dble(npointstotal_sq-1) + + !atomic units provided that G are in atomic units + k=1 + do i1=1,npointstotal_sq + r1=-0.5d0+dble(i1-1)*step + do i2=1,npointstotal_sq + r2=-0.5d0+dble(i2-1)*step + rkxvector(k)=r1*G(1,1)+r2*G(2,1) + rkyvector(k)=r1*G(1,2)+r2*G(2,2) + rkzvector(k)=0.0d0 + k=k+1 + end do + end do + + end subroutine get_grid + + +end module parser_optics_xatu_dim \ No newline at end of file diff --git a/src/sigma_first_ex.f90 b/src/sigma_first_ex.f90 index 8c18f55..bd4d8cf 100644 --- a/src/sigma_first_ex.f90 +++ b/src/sigma_first_ex.f90 @@ -1,7 +1,7 @@ module sigma_first_ex use constants_math use parser_input_file, & - only:iflag_xatu,nf,e1,e2,eta,nw + only:iflag_xatu,nf,e1,e2,eta,nw,broadening_type_text use parser_wannier90_tb, & only:material_name,norb use parser_optics_xatu_dim, & @@ -103,12 +103,21 @@ subroutine get_kubo_intens_ex(vme_ex,nw,wp,eta1,sigma_w_ex) do njp=1,3 !N integrand - skubo_ex_int(nj,njp,nn)=1.0d0/(dble(npointstotal)*vcell) & + skubo_ex_int(nj,njp,nn)=pi/(dble(npointstotal)*vcell) & *conjg(vme_ex(nj,nn))*vme_ex(njp,nn)/e_ex(nn) !pick the correct order of operators !at a given frequency !delta function - delta_n_ex=pi*1.0d0/eta1*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta1**2)*(wp(iw)-e_ex(nn))**2) + if (trim(broadening_type_text) == 'gaussian') then + delta_n_ex = pi*1.0d0/eta1*1.0d0/sqrt(2.0d0*pi)*& + exp(-0.5d0/(eta1**2)*(wp(iw)-e_ex(nn))**2) + else if (trim(broadening_type_text) == 'lorentzian') then + delta_n_ex = 1.0d0/pi*aimag(1.0d0/(wp(iw)-e_ex(nn)-& + complex(0.0d0,eta1))) + else + delta_n_ex = pi*1.0d0/eta1*1.0d0/sqrt(2.0d0*pi)*& + exp(-0.5d0/(eta1**2)*(wp(iw)-e_ex(nn))**2) + end if !sigma_w sigma_w_ex(nj,njp,iw)=sigma_w_ex(nj,njp,iw)+skubo_ex_int(nj,njp,nn)*delta_n_ex @@ -147,7 +156,8 @@ subroutine print_sigma_first_ex(nw,wp,sigma_w_ex) realpart(feps*sigma_w_ex(3,2,iw)), & realpart(feps*sigma_w_ex(3,3,iw)) - write(55,*) wp(iw)*27.211385d0,aimag(feps*sigma_w_ex(1,1,iw)), & + write(55,*) wp(iw)*27.211385d0, & + aimag(feps*sigma_w_ex(1,1,iw)), & aimag(feps*sigma_w_ex(1,2,iw)), & aimag(feps*sigma_w_ex(1,3,iw)), & aimag(feps*sigma_w_ex(2,1,iw)), & diff --git a/src/sigma_first_sp.f90 b/src/sigma_first_sp.f90 index 63a271b..84f8bf5 100644 --- a/src/sigma_first_sp.f90 +++ b/src/sigma_first_sp.f90 @@ -1,7 +1,7 @@ module sigma_first_sp use constants_math use parser_input_file, & - only:iflag_xatu,nf,e1,e2,eta,nw + only:iflag_xatu,nf,e1,e2,eta,nw,broadening_type_text use parser_wannier90_tb, & only:material_name,norb use parser_optics_xatu_dim, & @@ -176,8 +176,17 @@ subroutine get_kubo_intens_sp(nband_ex,npointstotal,vcell,e,vme,nw,wp,eta1,sigma else factor1=(fnn-fnnp)/(e(nn)-e(nnp)) end if - !gaussian broadening. Only real part of delta function is implemented! - delta_nnp=-1.0d0/eta1*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta1**2)*(wp(iw)-e(nn)+e(nnp))**2) + ! Broadening: gaussian or lorentzian based on parser input + if (trim(broadening_type_text) == 'gaussian') then + delta_nnp = -1.0d0/eta1*1.0d0/sqrt(2.0d0*pi)*& + exp(-0.5d0/(eta1**2)*(wp(iw)-e(nn)+e(nnp))**2) + else if (trim(broadening_type_text) == 'lorentzian') then + delta_nnp = 1.0d0/pi*aimag(1.0d0/(-wp(iw)+e(nn)-e(nnp)+& + complex(0.0d0,eta1))) + else + delta_nnp = -1.0d0/eta1*1.0d0/sqrt(2.0d0*pi)*& + exp(-0.5d0/(eta1**2)*(wp(iw)-e(nn)+e(nnp))**2) + end if !save oscillator stregths do nj=1,3 diff --git a/src/sigma_second_ex.f90 b/src/sigma_second_ex.f90 index f8f0aa7..5c540c9 100644 --- a/src/sigma_second_ex.f90 +++ b/src/sigma_second_ex.f90 @@ -1,7 +1,7 @@ module sigma_second_ex use constants_math - use parser_input_file, & - only:e1,e2,eta,nw,response_text + use parser_input_file, & + only:e1,e2,eta,nw,response_text,broadening_type_text use parser_wannier90_tb, & only:material_name use parser_optics_xatu_dim, & @@ -136,12 +136,20 @@ subroutine get_shift_kernel_ex(eta2,nj,njp,njpp,nn,nnp,omegap,omegaq,omega2,shif nj2=njp nj3=njpp - ilorentzian=0 - ihuang=0 - imine=1 !working gaussian-broadedning shift current +if (trim(broadening_type_text) == 'gaussian') then + imine = 1 + ilorentzian = 0 + else if (trim(broadening_type_text) == 'lorentzian') then + imine = 0 + ilorentzian = 1 + else + imine = 1 + ilorentzian = 0 + end if + ihuang=0 - if (imine.eq.1) then + if (imine.eq.1) then !with deltas. Working on 07/11/2023 d1=1.0d0/eta2*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta2**2)*(omegaq-e_ex(nnp))**2) d2=1.0d0/eta2*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta2**2)*(omegaq+e_ex(nnp))**2) diff --git a/src/sigma_second_sp.f90 b/src/sigma_second_sp.f90 index 6b577ae..d3fd6a2 100644 --- a/src/sigma_second_sp.f90 +++ b/src/sigma_second_sp.f90 @@ -1,7 +1,7 @@ module sigma_second_sp use constants_math use parser_input_file, & - only:nf,e1,e2,eta,nw,response_text + only:nf,e1,e2,eta,nw,response_text,broadening_type_text use parser_wannier90_tb, & only:material_name use parser_optics_xatu_dim, & @@ -198,9 +198,17 @@ subroutine get_shift_intens_sp(nband_ex,nw,e_nband,vme_nband,shift_vector_nband, else fnnp=0.0d0 end if - factor1=fnn-fnnp - !delta_nnp=1.0d0/pi*aimag(1.0d0/(wp(iw)-e(nn)+e(nnp)-complex(0.0d0,eta))) - delta_nnp=1.0d0/eta2*1.0d0/sqrt(2.0d0*pi)*exp(-0.5d0/(eta2**2)*(wp(iw)-e_nband(nn)+e_nband(nnp))**2) + factor1=fnn-fnnp + if (trim(broadening_type_text) == 'gaussian') then + delta_nnp = 1.0d0/eta2*1.0d0/sqrt(2.0d0*pi)*& + exp(-0.5d0/(eta2**2)*(wp(iw)-e_nband(nn)+e_nband(nnp))**2) + else if (trim(broadening_type_text) == 'lorentzian') then + delta_nnp = 1.0d0/pi*aimag(1.0d0/(wp(iw)-e_nband(nn)+& + e_nband(nnp)-complex(0.0d0,eta2))) + else + delta_nnp = 1.0d0/eta2*1.0d0/sqrt(2.0d0*pi)*& + exp(-0.5d0/(eta2**2)*(wp(iw)-e_nband(nn)+e_nband(nnp))**2) + end if do nj=1,3 do njp=1,3 shift_vector_w(nj,njp,iw)=shift_vector_w(nj,njp,iw)+ & From 8e3dd96ae32f255b79291b05e7a07c5a670e536a Mon Sep 17 00:00:00 2001 From: visitor Date: Fri, 20 Feb 2026 10:45:45 +0100 Subject: [PATCH 5/6] removed unused ints and arrays --- src/parser_input_file.f90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/parser_input_file.f90 b/src/parser_input_file.f90 index 9d3906f..f136cae 100644 --- a/src/parser_input_file.f90 +++ b/src/parser_input_file.f90 @@ -34,18 +34,14 @@ module parser_input_file logical :: iflag_ome_sp logical :: iflag_ome_ex - integer :: nk1 integer :: ndim integer :: nf integer :: npointstotal_sq - integer :: nband_index_aux integer :: norb_ex_cut - integer :: j !to read stuff integer :: nband_index integer :: nw real(8) :: e1,e2,eta - dimension :: nband_index_aux(100) !auxiliary array to save band indeces allocatable :: nband_index(:) contains @@ -68,7 +64,7 @@ end function to_lower subroutine get_input_file() implicit none integer, allocatable :: narray(:) - integer :: i, num_values, ios + integer :: num_values, ios character(len=1000) :: line character(len=100) :: param_name logical :: ndim_found, material_found, xatu_found, bandlist_found From e89de2d23b1d03274271def96252356fd9b12be7 Mon Sep 17 00:00:00 2001 From: visitor Date: Fri, 20 Feb 2026 19:06:35 +0100 Subject: [PATCH 6/6] fixed indexing of wannier parser --- src/parser_wannier90_tb.f90 | 304 ++++++++++++++++++------------------ 1 file changed, 152 insertions(+), 152 deletions(-) diff --git a/src/parser_wannier90_tb.f90 b/src/parser_wannier90_tb.f90 index 6c241f2..915d67b 100644 --- a/src/parser_wannier90_tb.f90 +++ b/src/parser_wannier90_tb.f90 @@ -1,161 +1,161 @@ module parser_wannier90_tb - use parser_input_file, only:read_line_numbers_int - implicit none - private - public :: norb,nR,n1,n2,n3,nRvec - public :: R,shop,hhop,rhop_c - public :: material_name - public :: wannier90_get - - integer norb - integer nR - integer n1,n2,n3 - integer nRvec - - real(8) R,shop - complex*16 hhop,rhop_c - character(1000) material_name - - dimension R(3,3) - integer, allocatable :: Degen(:) - allocatable nRvec(:,:) - allocatable hhop(:,:,:) - allocatable rhop_c(:,:,:,:) - allocatable shop(:,:,:) - - contains - -!-------------------------------------------------------------------- -pure function to_lower(str) result(out) + use parser_input_file, only:read_line_numbers_int + implicit none + private + public :: norb,nR,n1,n2,n3,nRvec + public :: R,shop,hhop,rhop_c + public :: material_name + public :: wannier90_get + + integer norb + integer nR + integer n1,n2,n3 + integer nRvec + + real(8) R,shop + complex*16 hhop,rhop_c + character(1000) material_name + + dimension R(3,3) + integer, allocatable :: Degen(:) + allocatable nRvec(:,:) + allocatable hhop(:,:,:) + allocatable rhop_c(:,:,:,:) + allocatable shop(:,:,:) + +contains + +!-------------------------------------------------------------------- + pure function to_lower(str) result(out) !> Convert a string to lower‑case (portable; no compiler extension) - character(len=*), intent(in) :: str - character(len=len(str)) :: out - integer :: i - out = str - do i = 1, len(str) - if (out(i:i) >= 'A' .and. out(i:i) <= 'Z') & - out(i:i) = char(iachar(out(i:i)) + 32) - end do -end function to_lower + character(len=*), intent(in) :: str + character(len=len(str)) :: out + integer :: i + out = str + do i = 1, len(str) + if (out(i:i) >= 'A' .and. out(i:i) <= 'Z') & + out(i:i) = char(iachar(out(i:i)) + 32) + end do + end function to_lower !-------------------------------------------------------------------- -subroutine wannier90_get(material_name_in) - implicit none - ! ----------------------------------------------------------------- - character(len=*), intent(in) :: material_name_in ! full path read from input.txt - ! ----------------------------------------------------------------- - integer :: fp, iR, i, ialpha, ialphap - integer :: nkk1, nkk2, nRzero - real(8) :: a1,a2,a3,a4,a5,a6 - character(len=:), allocatable :: file2open, basename - integer :: p, ext_pos - ! ----------------------------------------------------------------- - write(*,*) '2. Entering parser_wannier90_tb' - - ! === 1. Use the path exactly as supplied ======================== - !write(*,*) "MATERIAL NAME PARSED:", material_name_in - file2open = trim(material_name_in) - - ! === 2. Derive clean material name (no dir, no _tb.dat) ========= - p = max( index(file2open,'/',back=.true.), & - index(file2open,'\',back=.true.) ) ! works on Win/Linux - if (p == 0) then - basename = file2open - else - basename = file2open(p+1:) - end if - - ext_pos = len_trim(basename) - len('_tb.dat') + 1 - if ( ext_pos > 0 .and. to_lower(basename(ext_pos:)) == '_tb.dat' ) then - basename = basename(:ext_pos-1) - end if - material_name = adjustl(basename) - - ! === 3. Open Wannier90 TB file ================================== - open(unit=fp, file=file2open, action='read', status='old') - read(fp,*) - read(fp,*) R(1,1),R(1,2),R(1,3) - read(fp,*) R(2,1),R(2,2),R(2,3) - read(fp,*) R(3,1),R(3,2),R(3,3) - read(fp,*) norb - read(fp,*) nR - - !allocate nR, h-,r-, and s-hoppings - allocate (nRvec(nR,3)) - allocate (hhop(nR,norb,norb)) - allocate (shop(nR,norb,norb)) - allocate (rhop_c(3,nR,norb,norb)) - allocate (Degen(nR)) - - ! Read degeneracies - if ((nR / 15) > 1) then - do i = 1, (nR / 15) - read(fp, *) Degen((i - 1) * 15 + 1:(i - 1) * 15 + 15) - end do - end if - read(fp, *) Degen((i - 1) * 15 + 1:(i - 1) * 15 + MOD(nR, 15)) - - !get the hopping matrices - do iR=1,nR - read(fp,*) nRvec(iR,:) - do ialphap=1,norb - do ialpha=1,norb - read(fp,*) nkk1,nkk2,a1,a2 - hhop(iR,ialpha,ialphap)=complex(a1,a2) - end do - end do + subroutine wannier90_get(material_name_in) + implicit none + ! ----------------------------------------------------------------- + character(len=*), intent(in) :: material_name_in ! full path read from input.txt + ! ----------------------------------------------------------------- + integer :: fp, iR, i, ialpha, ialphap + integer :: nkk1, nkk2, nRzero + real(8) :: a1,a2,a3,a4,a5,a6 + character(len=:), allocatable :: file2open, basename + integer :: p, ext_pos + ! ----------------------------------------------------------------- + write(*,*) '2. Entering parser_wannier90_tb' + + ! === 1. Use the path exactly as supplied ======================== + !write(*,*) "MATERIAL NAME PARSED:", material_name_in + file2open = trim(material_name_in) + + ! === 2. Derive clean material name (no dir, no _tb.dat) ========= + p = max( index(file2open,'/',back=.true.), & + index(file2open,'\',back=.true.) ) ! works on Win/Linux + if (p == 0) then + basename = file2open + else + basename = file2open(p+1:) + end if + + ext_pos = len_trim(basename) - len('_tb.dat') + 1 + if ( ext_pos > 0 .and. to_lower(basename(ext_pos:)) == '_tb.dat' ) then + basename = basename(:ext_pos-1) + end if + material_name = adjustl(basename) + + ! === 3. Open Wannier90 TB file ================================== + open(unit=fp, file=file2open, action='read', status='old') read(fp,*) - end do - - !get rhoppings + read(fp,*) R(1,1),R(1,2),R(1,3) + read(fp,*) R(2,1),R(2,2),R(2,3) + read(fp,*) R(3,1),R(3,2),R(3,3) + read(fp,*) norb + read(fp,*) nR + + !allocate nR, h-,r-, and s-hoppings + allocate (nRvec(nR,3)) + allocate (hhop(nR,norb,norb)) + allocate (shop(nR,norb,norb)) + allocate (rhop_c(3,nR,norb,norb)) + allocate (Degen(nR)) + + ! Read degeneracies + if ((nR / 15) > 1) then + do i = 1, (nR / 15) + read(fp, *) Degen((i - 1) * 15 + 1:(i - 1) * 15 + 15) + end do + end if + read(fp, *) Degen((i - 1) * 15 + 1:(i - 1) * 15 + MOD(nR, 15)) + + !get the hopping matrices do iR=1,nR - read(fp,*) !nRvec is already strored - do ialphap=1,norb - do ialpha=1,norb - read(fp,*) nkk1,nkk2,a1,a2,a3,a4,a5,a6 - rhop_c(1,iR,ialpha,ialphap)=complex(a1,a2) - rhop_c(2,iR,ialpha,ialphap)=complex(a3,a4) - rhop_c(3,iR,ialpha,ialphap)=complex(a5,a6) - end do + read(fp,*) nRvec(iR,:) + do ialphap=1,norb + do ialpha=1,norb + read(fp,*) nkk1,nkk2,a1,a2 + hhop(iR,nkk1,nkk2)=complex(a1,a2) + end do + end do + read(fp,*) end do - if (iR /= nR) read(fp,*) !blank line - end do - close(fp) - !get orthogonal overlap: this variable is a reminiscent - !of the interface with the original crystal interface. - !I maintain the overlap matrix though - - !locate the (0,0,0) element of nRvec - do iR=1,nR - if (nRvec(iR,1)==0 .and. nRvec(iR,2)==0 .and. nRvec(iR,3)==0) then - nRzero=iR - exit - end if - end do - !wannier functions are orthonormal - shop=0.0d0 - do ialpha=1,norb - shop(nRzero,ialpha,ialpha)=1.0d0 - end do - - - !APPLY BIAS BY HAND - !do iR=1,nR + + !get rhoppings + do iR=1,nR + read(fp,*) !nRvec is already strored + do ialphap=1,norb + do ialpha=1,norb + read(fp,*) nkk1,nkk2,a1,a2,a3,a4,a5,a6 + rhop_c(1,iR,nkk1,nkk2)=complex(a1,a2) + rhop_c(2,iR,nkk1,nkk2)=complex(a3,a4) + rhop_c(3,iR,nkk1,nkk2)=complex(a5,a6) + end do + end do + if (iR /= nR) read(fp,*) !blank line + end do + close(fp) + !get orthogonal overlap: this variable is a reminiscent + !of the interface with the original crystal interface. + !I maintain the overlap matrix though + + !locate the (0,0,0) element of nRvec + do iR=1,nR + if (nRvec(iR,1)==0 .and. nRvec(iR,2)==0 .and. nRvec(iR,3)==0) then + nRzero=iR + exit + end if + end do + !wannier functions are orthonormal + shop=0.0d0 + do ialpha=1,norb + shop(nRzero,ialpha,ialpha)=1.0d0 + end do + + + !APPLY BIAS BY HAND + !do iR=1,nR !do ialpha=1,norb - !do ialphap=1,norb - !hhop(iR,ialpha,ialphap)=hhop(iR,ialpha,ialphap)-0.02d0*rhop_c(3,iR,ialpha,ialphap) - !end do + !do ialphap=1,norb + !hhop(iR,ialpha,ialphap)=hhop(iR,ialpha,ialphap)-0.02d0*rhop_c(3,iR,ialpha,ialphap) + !end do !end do - !end do - !do ialpha=1,norb + !end do + !do ialpha=1,norb !hhop(nRzero,ialpha,ialpha)=hhop(nRzero,ialpha,ialpha)-0.1d0*rhop_c(3,nRzero,ialpha,ialpha) - !end do - - !convert units: to Hartree and bohrs - hhop=hhop/27.211385d0 - rhop_c=rhop_c/0.52917721067121d0 - R=R/0.52917721067121d0 - write(*,*) ' Wannier hamiltonian has been read' - end subroutine wannier90_get - - -end module parser_wannier90_tb \ No newline at end of file + !end do + + !convert units: to Hartree and bohrs + hhop=hhop/27.211385d0 + rhop_c=rhop_c/0.52917721067121d0 + R=R/0.52917721067121d0 + write(*,*) ' Wannier hamiltonian has been read' + end subroutine wannier90_get + + +end module parser_wannier90_tb