From 27b1a3de62cdfd6bcbe1f220cc8fca3038f9196d Mon Sep 17 00:00:00 2001 From: QuanSheng Wu Date: Sat, 4 Nov 2023 23:55:52 +0800 Subject: [PATCH] try to add valley operator, not tested yet --- examples/TBG-7.34degree/bulkek.gnu0 | 8 +- src/ek_bulk.f90 | 308 ++++++++++++++++++++++++++++ src/ham_bulk.f90 | 85 ++++++++ src/main.f90 | 15 +- src/module.f90 | 14 +- src/readHmnR.f90 | 158 ++++++++++++++ src/readinput.f90 | 3 + 7 files changed, 583 insertions(+), 8 deletions(-) diff --git a/examples/TBG-7.34degree/bulkek.gnu0 b/examples/TBG-7.34degree/bulkek.gnu0 index f7071e56..f38df767 100644 --- a/examples/TBG-7.34degree/bulkek.gnu0 +++ b/examples/TBG-7.34degree/bulkek.gnu0 @@ -9,18 +9,18 @@ set pointsize 0.8 #set ylabel font ",24" set ylabel offset 0.5,0 set xrange [0: 0.51666] -emin= -1.0 -emax= 1.0 +emin= -1.882175 +emax= 1.922040 set ylabel "Energy (eV)" set yrange [ emin : emax ] set xtics ("G " 0.00000,"K " 0.21836,"M " 0.32755,"G " 0.51666) set arrow from 0.21836, emin to 0.21836, emax nohead set arrow from 0.32755, emin to 0.32755, emax nohead # please comment the following lines to plot the fatband -plot 'bulkek.dat' u 1:2 w p pt 7 ps 0.2 lc rgb 'black', 0 w l lw 2 +#plot 'bulkek.dat' u 1:2 w p pt 7 ps 0.2 lc rgb 'black', 0 w l lw 2 # uncomment the following lines to plot the fatband -#plot 'bulkek.dat' u 1:2:3 w lp lw 2 pt 7 ps 0.2 lc palette, 0 w l lw 2 +plot 'bulkek.dat' u 1:2:3 w lp lw 2 pt 7 ps 0.2 lc palette, 0 w l lw 2 # uncomment the following lines to plot the spin if necessary #plot 'bulkek.dat' u 1:2 w lp lw 2 pt 7 ps 0.2, \ 'bulkek.dat' u 1:2:($3/6):($4/6) w vec diff --git a/src/ek_bulk.f90 b/src/ek_bulk.f90 index 074ca2df..c7382864 100644 --- a/src/ek_bulk.f90 +++ b/src/ek_bulk.f90 @@ -175,6 +175,115 @@ subroutine ek_bulk_line return end subroutine ek_bulk_line + +subroutine ek_bulk_line_valley + ! Calculate bulk's energy bands using wannier TB method + ! Line mode + ! Copyright (c) 2010 QuanSheng Wu. All rights reserved. + + use wmpi + use para + + implicit none + + integer :: ik, il, ig, io, i, j, knv3, ierr + real(dp) :: emin, emax, k(3) + character*40 :: filename + complex(dp), external :: zdotc + + !> eigenvalues of H + real(Dp), allocatable :: W(:) + + ! Hamiltonian of bulk system + complex(Dp), allocatable :: psi(:), vpsi(:) + complex(Dp), allocatable :: Hamk_bulk(:, :), valley_k(:, :) + + ! eigenectors of H + real(dp), allocatable :: eigv(:,:), eigv_mpi(:,:) + real(dp), allocatable :: weight(:,:), weight_mpi(:,:) + + knv3= nk3_band + allocate(W(Num_wann)) + allocate(Hamk_bulk(Num_wann, Num_wann)) + allocate(valley_k(Num_wann, Num_wann)) + allocate( eigv (Num_wann, knv3)) + allocate( eigv_mpi(Num_wann, knv3)) + allocate( weight (Num_wann, knv3)) + allocate( weight_mpi(Num_wann, knv3)) + allocate( psi(Num_wann), vpsi(Num_wann)) + W = 0d0; Hamk_bulk = 0d0 + eigv = 0d0; eigv_mpi= 0d0 + weight = 0d0; weight_mpi = 0d0 + psi = 0d0; vpsi= 0d0 + + + do ik= 1+cpuid, knv3, num_cpu + + k = k3points(:, ik) + + ! calculation bulk hamiltonian + Hamk_bulk= 0d0 + + call ham_bulk_atomicgauge(k, Hamk_bulk) + + !> diagonalization by call zheev in lapack + W= 0d0 + call eigensystem_c('V', 'U', Num_wann ,Hamk_bulk, W) + + !> get the fourier transform of a valley operator + call valley_k_atomicgauge(k, valley_k) + + eigv(:, ik)= W + do j=1, Num_wann !> band + psi= Hamk_bulk(:, j) + call zgemv('N', Num_wann, Num_wann, One_complex, psi, Num_wann, valley_k, 1, zzero, vpsi, 1) + weight(j, ik)= real(zdotc(Num_wann, psi, 1, vpsi, 1)) + enddo ! i + enddo ! ik + +#if defined (MPI) + call mpi_allreduce(eigv,eigv_mpi,size(eigv),& + mpi_dp,mpi_sum,mpi_cmw,ierr) + call mpi_allreduce(weight, weight_mpi,size(weight),& + mpi_dp,mpi_sum,mpi_cmw,ierr) +#else + eigv_mpi= eigv + weight_mpi= weight +#endif + eigv_mpi= eigv_mpi/eV2Hartree + + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek.dat') + + do i=1, Num_wann + do ik=1, knv3 + write(outfileindex, '(200E16.5)')k3len(ik)*Angstrom2atomic,eigv_mpi(i, ik), & + weight_mpi(i, ik) + enddo + write(outfileindex, *)' ' + enddo + close(outfileindex) + endif + + !> minimum and maximum value of energy bands + emin= minval(eigv_mpi)-0.5d0 + emax= maxval(eigv_mpi)+0.5d0 + + call generate_ek_kpath_gnu('bulkek.dat', 'bulkek.gnu', 'bulkek.pdf', & + emin, emax, knv3, Nk3lines, & + k3line_name, k3line_stop, k3len) + + deallocate(W) + deallocate(Hamk_bulk) + deallocate( eigv ) + deallocate( eigv_mpi) + deallocate( weight ) + deallocate( weight_mpi) + + return +end subroutine ek_bulk_line_valley + !>> calculate energy band levels at given kpoints subroutine ek_bulk_point_mode @@ -1107,6 +1216,205 @@ subroutine sparse_ekbulk end subroutine sparse_ekbulk +!> only test for valley projection +subroutine sparse_ekbulk_valley + use sparse + use para + implicit none + + + !> some temporary integers + integer :: ik, i, j, ierr, ib, ig + + ! wave vector + real(dp) :: k3(3) + + !> dim= Num_wann, knv3 + real(dp), allocatable :: W(:) + real(dp), allocatable :: eigv(:, :) + real(dp), allocatable :: eigv_mpi(:, :) + + real(dp) :: emin, emax + + complex(dp) :: alpha_mv, beta_mv + + !> dim= Num_wann*Num_wann + integer :: nnzmax, nnz, nnzmax_valley, nnz_valley + complex(dp), allocatable :: acoo(:), acoo_valley(:) + integer, allocatable :: jcoo(:), jcoo_valley(:) + integer, allocatable :: icoo(:), icoo_valley(:) + + !> eigenvector of the sparse matrix acoo. Dim=(Num_wann, neval) + complex(dp), allocatable :: psi(:), vpsi(:) + complex(dp), allocatable :: psi_project(:) + complex(dp), allocatable :: zeigv(:, :) + + !> print the weight for the Selected_WannierOrbitals + real(dp), allocatable :: weight_valley(:, :), weight_valley_mpi(:, :) + + !number of ARPACK eigenvalues + integer :: neval + + ! number of Arnoldi vectors + integer :: nvecs + + !> calculate eigenvector or not + logical :: ritzvec + + !shift-invert sigma + complex(dp) :: sigma=(0d0,0d0) + + !> time measurement + real(dp) :: time1, time2, time3 + + !> + character :: matdescra(6) + + ! vector multiply vector v1*v2 + complex(dp), external :: zdotc + + + !if (OmegaNum==0) OmegaNum= Num_wann + !if (NumSelectedEigenVals==0) NumSelectedEigenVals=OmegaNum + + !> first use NumSelectedEigenVals, if NumSelectedEigenVals is not set, + !> then use OmegaNum; if OmegaNum is also not set, + !> then use Num_wann + if (NumSelectedEigenVals>0) then + neval= NumSelectedEigenVals + else if (OmegaNum>0) then + neval= OmegaNum + else + neval = Num_wann + endif + + if (neval>Num_wann-2) neval= Num_wann- 2 + + !> ncv + nvecs=int(2*neval) + + if (nvecs<20) nvecs= 20 + if (nvecs>Num_wann) nvecs= Num_wann + + sigma=(1d0,0d0)*E_arc + nnzmax=splen+Num_wann + nnz=splen + allocate( acoo(nnzmax)) + allocate( jcoo(nnzmax)) + allocate( icoo(nnzmax)) + + nnzmax_valley = splen_valley_input + allocate( acoo_valley(nnzmax_valley)) + allocate( jcoo_valley(nnzmax_valley)) + allocate( icoo_valley(nnzmax_valley)) + + allocate( W( neval)) + allocate( eigv( neval, nk3_band)) + allocate( eigv_mpi( neval, nk3_band)) + allocate( psi(Num_wann), vpsi(Num_wann)) + allocate( psi_project(Num_wann)) + allocate( zeigv(Num_wann,nvecs)) + allocate( weight_valley(neval, nk3_band), weight_valley_mpi(neval, nk3_band)) + psi=0d0; vpsi=0d0; psi_project= 0d0; zeigv= 0d0 + weight_valley_mpi=0d0; weight_valley=0d0 + + eigv_mpi= 0d0; eigv = 0d0 + acoo= 0d0; icoo=0; jcoo=0 + + ritzvec= .True. + + matdescra(1)= 'H' + matdescra(2)= 'U' + matdescra(3)= 'U' + matdescra(4)= 'F' + + !> calculate the energy bands along special k line + k3= 0 + do ik=1+ cpuid, nk3_band, num_cpu + if (cpuid==0) write(stdout, '(a, 2i10)') 'BulkBand_calc in sparse mode:', ik,nk3_band + k3 = K3points(:, ik) + call now(time1) + call ham_bulk_coo_sparsehr(k3,acoo,icoo,jcoo) + acoo= acoo/eV2Hartree + nnz= splen + call now(time2) + + !> diagonalization by call zheev in lapack + W= 0d0 + !> after arpack_sparse_coo_eigs, nnz will be updated. + call arpack_sparse_coo_eigs(Num_wann,nnzmax,nnz,acoo,jcoo,icoo,neval,nvecs,W,sigma, zeigv, ritzvec) + call now(time3) + eigv(1:neval, ik)= W(1:neval) + + !> get valley operator in coo format + call valley_k_coo_sparsehr(nnzmax_valley, k3, acoo_valley, icoo_valley, jcoo_valley) + + !> get the valley projection + do ib= 1, neval + psi(:)= zeigv(:, ib) !> the eigenvector of ib'th band + + !> weight_valley= + call mkl_zcoomv('N', Num_wann, Num_wann, One_complex, matdescra, acoo_valley, & + icoo_valley, jcoo_valley, nnzmax_valley, psi, zzero, vpsi) + + weight_valley(ib, ik)= real(zdotc(Num_wann, psi, 1, vpsi, 1)) + + enddo + + if (cpuid==0)write(stdout, '(a, f20.2, a)')' >> Time cost for constructing H: ', time2-time1, ' s' + if (cpuid==0)write(stdout, '(a, f20.2, a)')' >> Time cost for diagonalize H: ', time3-time2, ' s' + enddo !ik + +#if defined (MPI) + call mpi_allreduce(eigv,eigv_mpi,size(eigv),& + mpi_dp,mpi_sum,mpi_cmw,ierr) + call mpi_allreduce(weight_valley, weight_valley_mpi,size(weight_valley),& + mpi_dp,mpi_sum,mpi_cmw,ierr) +#else + eigv_mpi= eigv + weight_valley_mpi= weight_valley +#endif + + !> minimum and maximum value of energy bands + emin= minval(eigv_mpi)+0.05*(maxval(eigv_mpi)-minval(eigv_mpi)) + emax= maxval(eigv_mpi)-0.05*(maxval(eigv_mpi)-minval(eigv_mpi)) + + + + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek.dat') + do i=1, neval + do ik=1, nk3_band + write(outfileindex, '(300f16.9)')k3len(ik)*Angstrom2atomic,eigv_mpi(i, ik), & + weight_valley_mpi(i, ik) + enddo + write(outfileindex, *)' ' + enddo + close(outfileindex) + endif + + call generate_ek_kpath_gnu('bulkek.dat', 'bulkek.gnu', 'bulkek.pdf', & + emin, emax, nk3_band, Nk3lines, & + k3line_name, k3line_stop, k3len) + +#if defined (MPI) + call mpi_barrier(mpi_cmw, ierr) +#endif + + deallocate( acoo) + deallocate( jcoo) + deallocate( icoo) + deallocate( W) + deallocate( eigv) + deallocate( eigv_mpi) + deallocate( zeigv) + deallocate( weight_valley, weight_valley_mpi) + + return +end subroutine sparse_ekbulk_valley + + subroutine sparse_ekbulk_plane use sparse use para diff --git a/src/ham_bulk.f90 b/src/ham_bulk.f90 index 1d79e98e..8c553c99 100644 --- a/src/ham_bulk.f90 +++ b/src/ham_bulk.f90 @@ -92,6 +92,63 @@ subroutine ham_bulk_atomicgauge(k,Hamk_bulk) return end subroutine ham_bulk_atomicgauge + +subroutine valley_k_atomicgauge(k,valley_k) + ! This subroutine caculates Hamiltonian for + ! bulk system with the consideration of the atom's position + ! + ! History + ! + ! May/29/2011 by Quansheng Wu + ! Atomic gauge Guan Yifei 2019 + ! Lattice gauge Hl + ! Atomic gauge Ha= U* Hl U + ! where U = e^ik.wc(i) on diagonal + + use para + implicit none + + integer :: i1,i2,iR + + ! wave vector in 3d + real(Dp) :: k(3), kdotr, pos0(3), dis + + complex(dp) :: factor + + real(dp) :: pos(3), pos1(3), pos2(3), pos_cart(3), pos_direct(3) + ! Hamiltonian of bulk system + complex(Dp),intent(out) :: valley_k(Num_wann, Num_wann) + real(dp), external :: norm + + + valley_k= 0d0 + !> the first atom in home unit cell + do iR=1, Nrpts_valley + do i2=1, Num_wann + pos2= Origin_cell%wannier_centers_direct(:, i2) + !> the second atom in unit cell R + do i1=1, Num_wann + pos1= Origin_cell%wannier_centers_direct(:, i1) + pos_direct= irvec_valley(:, iR) + pos_direct= pos_direct+ pos2- pos1 + + call direct_cart_real(pos_direct, pos_cart, Origin_cell%lattice) + + dis= norm(pos_cart) + if (dis> Rcut) cycle + + kdotr=k(1)*pos_direct(1) + k(2)*pos_direct(2) + k(3)*pos_direct(3) + factor= (cos(twopi*kdotr)+zi*sin(twopi*kdotr)) + + valley_k(i1, i2)= valley_k(i1, i2) & + + valley_operator_R(i1, i2, iR)*factor + enddo ! i1 + enddo ! i2 + enddo ! iR + + return +end subroutine valley_k_atomicgauge + subroutine d2Hdk2_atomicgauge(k, DHDk2_wann) !> second derivatve of H(k) use para, only : Nrpts, irvec, crvec, Origin_cell, HmnR, ndegen, & @@ -746,6 +803,34 @@ subroutine ham_bulk_coo_sparsehr(k,acoo,icoo,jcoo) return end subroutine ham_bulk_coo_sparsehr +subroutine valley_k_coo_sparsehr(nnz, k,acoo,icoo,jcoo) + !> This subroutine use sparse hr format + !> Here we use atomic gauge which means the atomic position is taken into account + !> in the Fourier transformation + use para + implicit none + + real(dp) :: k(3), posij(3) + real(dp) :: kdotr + integer, intent(in) :: nnz + integer,intent(inout) :: icoo(nnz),jcoo(nnz) + complex(dp),intent(inout) :: acoo(nnz) + complex(dp) :: ratio + integer :: i,j,ir + + do i=1,nnz + ir=valley_operator_irv(i) + icoo(i)=valley_operator_icoo(i) + jcoo(i)=valley_operator_jcoo(i) + posij=irvec_valley(:, ir)+ Origin_cell%wannier_centers_direct(:, jcoo(i))- Origin_cell%wannier_centers_direct(:, icoo(i)) + kdotr=posij(1)*k(1)+posij(2)*k(2)+posij(3)*k(3) + ratio= (cos(twopi*kdotr)+zi*sin(twopi*kdotr)) + acoo(i)=ratio*valley_operator_acoo(i) + end do + + return +end subroutine valley_k_coo_sparsehr + subroutine rotation_to_Ham_basis(UU, mat_wann, mat_ham) !> this subroutine rotate the matrix from Wannier basis to Hamiltonian basis diff --git a/src/main.f90 b/src/main.f90 index 82633ec5..423b45fc 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -163,9 +163,12 @@ program main allocate(HmnR(num_wann,num_wann,nrpts)) HmnR= 0d0 call readNormalHmnR() + if (valley_projection_calc) call read_valley_operator !> sparse hmnr input else call readSparseHmnR() + if (valley_projection_calc) call readsparse_valley_operator + end if else stop "We only support Is_HrFile=.true. for this version" @@ -216,12 +219,20 @@ program main if(cpuid.eq.0)write(stdout, *)' ' if(cpuid.eq.0)write(stdout, *)'>> Start of calculating bulk band' call now(time_start) - if(Is_Sparse_Hr) then + if (Is_Sparse_Hr) then #if defined (INTELMKL) + if (valley_projection_calc) then + call sparse_ekbulk_valley + else call sparse_ekbulk + endif #endif else - call ek_bulk_line + if (valley_projection_calc) then + call ek_bulk_line_valley + else + call ek_bulk_line + endif !call ek_bulk_spin !call ek_bulk_mirror_z end if diff --git a/src/module.f90 b/src/module.f90 index bb813c04..3fcee347 100644 --- a/src/module.f90 +++ b/src/module.f90 @@ -445,6 +445,7 @@ module para logical :: landau_chern_calc = .false. logical :: export_newhr = .false. logical :: export_maghr =.false. + logical :: valley_projection_calc ! for valley projection, only for sparse hamiltonina, you need to provide the valley operator logical :: w3d_nested_calc = .false. @@ -478,7 +479,8 @@ module para LanczosSeqDOS_calc, Translate_to_WS_calc, LandauLevel_k_dos_calc, & LandauLevel_B_dos_calc,LanczosBand_calc,LanczosDos_calc, & LandauLevel_B_calc, LandauLevel_kplane_calc,landau_chern_calc, & - FermiLevel_calc,ANE_calc, export_newhr,export_maghr,w3d_nested_calc + FermiLevel_calc,ANE_calc, export_newhr,export_maghr,w3d_nested_calc, & + valley_projection_calc integer :: Nslab ! Number of slabs for 2d Slab system integer :: Nslab1 ! Number of slabs for 1D wire system @@ -498,6 +500,7 @@ module para integer :: Num_wann ! Number of Wannier functions integer :: Nrpts ! Number of R points + integer :: Nrpts_valley ! Number of R points integer :: Nk ! number of k points for different use @@ -651,6 +654,7 @@ module para real(dp),parameter :: eps9= 1e-9 ! 0.000000001 real(dp),parameter :: eps12= 1e-12 ! 0.000000000001 complex(dp),parameter :: zzero= (0.0d0, 0d0) ! 0 + complex(dp),parameter :: One_complex= (1.0d0, 0d0) ! 0 real(Dp),parameter :: Ka(2)=(/1.0d0,0.0d0/) real(Dp),parameter :: Kb(2)=(/0.0d0,1.0d0/) @@ -851,16 +855,22 @@ module para real(dp) :: K3D_vec3_cube(3) ! the 3rd k vector for the k cube integer, allocatable :: irvec(:,:) ! R coordinates in fractional units + integer, allocatable :: irvec_valley(:,:) ! R coordinates in fractional units real(dp), allocatable :: crvec(:,:) ! R coordinates in Cartesian coordinates in units of Angstrom complex(dp), allocatable :: HmnR(:,:,:) ! Hamiltonian m,n are band indexes + complex(dp), allocatable :: valley_operator_R(:,:,:) ! Hamiltonian m,n are band indexes !sparse HmnR arraies integer,allocatable :: hicoo(:),hjcoo(:),hirv(:) complex(dp),allocatable :: hacoo(:) + !> valley operator + integer,allocatable :: valley_operator_icoo(:), valley_operator_jcoo(:), valley_operator_irv(:) + complex(dp),allocatable :: valley_operator_acoo(:) + !> sparse hr length - integer :: splen, splen_input + integer :: splen, splen_input, splen_valley_input integer, allocatable :: ndegen(:) ! degree of degeneracy of R point diff --git a/src/readHmnR.f90 b/src/readHmnR.f90 index 19f1fdd2..240b118a 100644 --- a/src/readHmnR.f90 +++ b/src/readHmnR.f90 @@ -236,6 +236,77 @@ subroutine readNormalHmnR() return end subroutine readNormalHmnR +subroutine read_valley_operator() + !>> Read in the tight-binding model from wannier90_hr.dat + !> The format is defined by the wannier90 software + ! Constructed by quansheng wu 4/2/2010 + ! + ! Yifei Guan added the sparse hr file parsing June/2018 + ! License: GPL V3 + + use para + !> in: N of wann + !> out : nth atom + + implicit none + + character*4 :: c_temp + + integer :: i, j, ir, ia, io + integer :: i1, i2, i3, i4, i5 + integer :: n, m, ir0 + integer :: add_electric_field + integer :: nwann, nwann_nsoc + + real(dp) :: static_potential + real(dp) :: tot, rh, ih + real(dp) :: pos(Origin_cell%Num_atoms) + + + if(cpuid.eq.0)write(stdout,*)' ' + open(12, file="valley_operator.dat", status='OLD') + + !> skip a comment line + read(12, *) + + !> number of Wannier orbitals in the hr file + nwann=0 + read(12, *)nwann + if (nwann==0) then + stop "ERROR : num_wann is zero in hr file" + endif + nwann_nsoc=nwann + if (SOC>0) nwann_nsoc= nwann/2 + + !> number of lattice vectors taken into account + read(12, *)Nrpts_valley + + !> The degeneracy of each R point, this is caused by the Fourier transform in the Wigner-Seitz cell + read(12, *)(j, i=1, Nrpts_valley) + allocate(irvec_valley(3, Nrpts_valley)) + allocate(valley_operator_R(num_wann, num_wann, Nrpts_valley)) + ir=0 + do ir=1,Nrpts_valley + do n=1,nwann + do m=1,nwann + read(12,*,end=1001)i1, i2, i3, i4, i5, rh, ih + irvec_valley(1,ir)=i1 + irvec_valley(2,ir)=i2 + irvec_valley(3,ir)=i3 + valley_operator_R(i4,i5,ir)=dcmplx(rh,ih) + end do + enddo + enddo + stop + + 1001 continue + close(12) + + if (cpuid.eq.0) write(stdout, *) ">>> finished reading of valley operator" + +end subroutine read_valley_operator + + subroutine get_hmnr_cell(cell) !> Get new hmnr for a new cell with the same size as the previous one @@ -650,6 +721,93 @@ subroutine readSparseHmnR end subroutine readSparseHmnR +subroutine readsparse_valley_operator + !> This subroutine not just read the sparse valley operator + !> 2023 Nov. 4 + use para + implicit none + integer:: i,j,nwann,nwann_nsoc,i1,i2,i3,i4,i5,ir,n,m + real(dp) :: r1, r2 + + !> the direction which adding electric field which is also the stacking direction + integer :: add_electric_field + real(dp) :: Bx_in_au, By_in_au, Bz_in_au + complex(dp) :: h_value + + open(13, file='valley_operator.dat') + + !> skip a comment line + read(13, *) + + !> comparing with the standard hr file, we add another line to show howmany + !> lines that Hmn(R) is not zero. + if(Is_Sparse_Hr) then + read(13,*) splen_valley_input !> number of non-zeros lines + end if + + !> number of Wannier orbitals in the hr file + nwann=0 + read(13, *)nwann + if (nwann==0.or.nwann.ne.Num_wann) then + print *, 'nwann, num_wann', nwann, num_wann + stop "ERROR : num_wann is zero or num_wann is not consistent with hr.dat" + endif + nwann_nsoc=nwann + if (SOC>0) nwann_nsoc= nwann/2 + + !> number of lattice vectors taken into account + read(13, *)Nrpts_valley + allocate(irvec_valley(3, Nrpts_valley)) + + !> The degeneracy of each R point + read(13, *)(j, i=1, nrpts) + ir=0 + + allocate(valley_operator_acoo(splen_valley_input)) + allocate(valley_operator_icoo(splen_valley_input)) + allocate(valley_operator_jcoo(splen_valley_input)) + allocate(valley_operator_irv(splen_valley_input)) + valley_operator_acoo=(0d0, 0d0) + valley_operator_icoo=0 + valley_operator_jcoo=0 + valley_operator_irv=0 + + j=0 + ir=1 + irvec_valley=-9999 + read(13,*,end=1002)i1, i2, i3, i4, i5, r1, r2 + irvec_valley(1,ir)=i1 + irvec_valley(2,ir)=i2 + irvec_valley(3,ir)=i3 + !> will reread the above line + backspace(13) + do i=1,Nrpts_valley + do n=1,nwann + do m=1,nwann + read(13,*,end=1002)i1, i2, i3, i4, i5, r1, r2 + if (sum(abs(irvec_valley(:,ir)-[i1,i2,i3]))/=0) ir=ir+1 + j=j+1 + valley_operator_icoo(j)=i4 + valley_operator_jcoo(j)=i5 + valley_operator_irv (j)=ir + valley_operator_acoo(j)=dcmplx(r1,r2) + irvec_valley(1, ir)= i1 + irvec_valley(2, ir)= i2 + irvec_valley(3, ir)= i3 + end do + enddo + enddo +1002 continue + !> correct nrpts + Nrpts_valley=ir + + if (cpuid.eq.0) write(stdout, '(a, i6)')' >> Nrpts_valley is ', Nrpts_valley + if (cpuid.eq.0) write(stdout, '(a)')' >> Valley operator reading finished ' + + return +end subroutine readsparse_valley_operator + + !> return the maximum distance between two neighbour elements of array x !> dx_max= mod(dx_max, 1) !> x is in [0, 1) diff --git a/src/readinput.f90 b/src/readinput.f90 index 17ebb5bb..c6affafb 100644 --- a/src/readinput.f90 +++ b/src/readinput.f90 @@ -153,6 +153,7 @@ subroutine readinput FermiLevel_calc = .FALSE. ANE_calc = .FALSE. w3d_nested_calc =.false. + valley_projection_calc =.false. ChargeDensity_selected_bands_calc= .FALSE. ChargeDensity_selected_energies_calc= .FALSE. @@ -208,6 +209,7 @@ subroutine readinput write(*, *)"Translate_to_WS_calc" write(*, *)"FermiLevel_calc" write(*, *)"ANE_calc" + write(*, *)"valley_projection_calc" write(*, *)"ChargeDensity_selected_energies_calc" write(*, *)"ChargeDensity_selected_bands_calc" write(*, *)"The default Vaule is F" @@ -287,6 +289,7 @@ subroutine readinput write(stdout, *) "Symmetry_Import_calc : ", Symmetry_Import_calc write(stdout, *) "ChargeDensity_selected_bands_calc : ", ChargeDensity_selected_bands_calc write(stdout, *) "ChargeDensity_selected_energies_calc : ", ChargeDensity_selected_energies_calc + write(stdout, *) "valley_projection_calc : " , valley_projection_calc endif !===============================================================================================================!