diff --git a/soc/fermiarc.f90 b/soc/fermiarc.f90 index 12749957..63706dc9 100644 --- a/soc/fermiarc.f90 +++ b/soc/fermiarc.f90 @@ -19,7 +19,7 @@ subroutine fermiarc integer :: ierr ! general loop index - integer :: i,j + integer :: i,j,io ! kpoint loop index integer :: ikp @@ -103,12 +103,16 @@ subroutine fermiarc call surfgreen_1985(omega,GLL,GRR,H00,H01,ones) ! call surfgreen_1984(omega,GLL,GRR,H00,H01,ones) - ! calculate spectral function - do i= 1, ndim - dos_l(ikp)=dos_l(ikp)- aimag(GLL(i,i)) - dos_r(ikp)=dos_r(ikp)- aimag(GRR(i,i)) - enddo + do i= 1, NtopOrbitals + io= TopOrbitals(i) + dos_l(ikp)=dos_l(ikp)- aimag(GLL(io,io)) + enddo ! i + do i= 1, NBottomOrbitals + io= Ndim- Num_wann+ BottomOrbitals(i) + dos_r(ikp)=dos_r(ikp)- aimag(GRR(io,io)) + enddo ! i + enddo !> we don't have to do allreduce operation @@ -129,7 +133,7 @@ subroutine fermiarc close(12) close(13) write(stdout,*)'ndim',ndim - write(stdout,*) 'Nk1,Nk2,eta',Nk1, Nk2, eta + write(stdout,*)'Nk1,Nk2,eta', Nk1, Nk2, eta write(stdout,*)'calculate density of state successfully' endif @@ -220,7 +224,7 @@ subroutine fermiarc_jdos integer :: ierr ! general loop index - integer :: i,j + integer :: i,j,io integer :: Nwann ! kpoint loop index @@ -382,8 +386,6 @@ subroutine fermiarc_jdos enddo enddo - - omega = E_arc eta= eta_arc @@ -405,10 +407,15 @@ subroutine fermiarc_jdos ik2= ik12(2, ikp) ! calculate spectral function - do i= 1, ndim - dos_l(ik1, ik2)=dos_l(ik1, ik2)- aimag(GLL(i,i)) - dos_r(ik1, ik2)=dos_r(ik1, ik2)- aimag(GRR(i,i)) - enddo + do i= 1, NtopOrbitals + io= TopOrbitals(i) + dos_l(ik1, ik2)=dos_l(ik1, ik2)- aimag(GLL(io,io)) + enddo ! i + do i= 1, NBottomOrbitals + io= Ndim- Num_wann+ BottomOrbitals(i) + dos_r(ik1, ik2)=dos_r(ik1, ik2)- aimag(GRR(io,io)) + enddo ! i + !! calculate spectral function !do i= 1, Nwann ! dos_l(ik1, ik2)=dos_l(ik1, ik2)- aimag(GLL(i,i)) & @@ -422,39 +429,68 @@ subroutine fermiarc_jdos !ctemp=matmul(surfgreen,sigma_x) call mat_mul(ndim,GLL,sigma_x,ctemp) do j=1,ndim - sx_l(ik1, ik2)=sx_l(ik1, ik2)-Aimag(ctemp(j,j)) + ! sx_l(ik1, ik2)=sx_l(ik1, ik2)-Aimag(ctemp(j,j)) enddo + do i= 1, NtopOrbitals + io= TopOrbitals(i) + sx_l(ik1, ik2)=sx_l(ik1, ik2)- aimag(ctemp(io,io)) + enddo ! i !ctemp=matmul(surfgreen,sigma_y) call mat_mul(ndim,GLL,sigma_y,ctemp) do j=1,ndim - sy_l(ik1, ik2)=sy_l(ik1, ik2)-Aimag(ctemp(j,j)) + ! sy_l(ik1, ik2)=sy_l(ik1, ik2)-Aimag(ctemp(j,j)) enddo + do i= 1, NtopOrbitals + io= TopOrbitals(i) + sy_l(ik1, ik2)=sy_l(ik1, ik2)- aimag(ctemp(io,io)) + enddo ! i + !ctemp=matmul(surfgreen,sigma_z) call mat_mul(ndim,GLL,sigma_z,ctemp) do j=1,ndim - sz_l(ik1, ik2)=sz_l(ik1, ik2)-Aimag(ctemp(j,j)) + ! sz_l(ik1, ik2)=sz_l(ik1, ik2)-Aimag(ctemp(j,j)) enddo + do i= 1, NtopOrbitals + io= TopOrbitals(i) + sz_l(ik1, ik2)=sz_l(ik1, ik2)- aimag(ctemp(io,io)) + enddo ! i + !ctemp=matmul(surfgreen,sigma_x) call mat_mul(ndim,GRR,sigma_x,ctemp) do j=1,ndim - sx_r(ik1, ik2)=sx_r(ik1, ik2)-Aimag(ctemp(j,j)) + ! sx_r(ik1, ik2)=sx_r(ik1, ik2)-Aimag(ctemp(j,j)) enddo + do i= 1, NBottomOrbitals + io= Ndim- Num_wann+ BottomOrbitals(i) + sx_r(ik1, ik2)=sx_r(ik1, ik2)- aimag(ctemp(io,io)) + enddo ! i + !ctemp=matmul(surfgreen,sigma_y) call mat_mul(ndim,GRR,sigma_y,ctemp) do j=1,ndim - sy_r(ik1, ik2)=sy_r(ik1, ik2)-Aimag(ctemp(j,j)) + ! sy_r(ik1, ik2)=sy_r(ik1, ik2)-Aimag(ctemp(j,j)) enddo + do i= 1, NBottomOrbitals + io= Ndim- Num_wann+ BottomOrbitals(i) + sy_r(ik1, ik2)=sy_r(ik1, ik2)- aimag(ctemp(io,io)) + enddo ! i + !ctemp=matmul(surfgreen,sigma_z) call mat_mul(ndim,GRR,sigma_z,ctemp) do j=1,ndim - sz_r(ik1, ik2)=sz_r(ik1, ik2)-Aimag(ctemp(j,j)) + ! sz_r(ik1, ik2)=sz_r(ik1, ik2)-Aimag(ctemp(j,j)) enddo + do i= 1, NBottomOrbitals + io= Ndim- Num_wann+ BottomOrbitals(i) + sz_r(ik1, ik2)=sz_r(ik1, ik2)- aimag(ctemp(io,io)) + enddo ! i + enddo !> we don't have to do allreduce operation diff --git a/soc/module.f90 b/soc/module.f90 index 2af38e5c..3638536d 100644 --- a/soc/module.f90 +++ b/soc/module.f90 @@ -180,6 +180,17 @@ module para namelist / KPATH_BULK / nk3lines, k3line_name, k3line_start + + !>> top surface atoms + integer :: NtopAtoms, NtopOrbitals + integer, allocatable :: TopAtoms(:) + integer, allocatable :: TopOrbitals(:) + + !>> bottom surface atoms + integer :: NBottomAtoms, NBottomOrbitals + integer, allocatable :: BottomAtoms(:) + integer, allocatable :: BottomOrbitals(:) + !>> effective mass !> k step for effective mass calculation @@ -256,6 +267,9 @@ module para integer, allocatable :: nprojs(:) character(10), allocatable :: proj_name(:, :) + !> the start index for each atoms, only consider the spinless component + integer, allocatable :: orbitals_start(:) + !> symmetry operator apply on function basis complex(dp), allocatable :: inversion(:, :) complex(dp), allocatable :: mirror_x(:, :) diff --git a/soc/readinput.f90 b/soc/readinput.f90 index 692a788c..037abbc4 100644 --- a/soc/readinput.f90 +++ b/soc/readinput.f90 @@ -20,7 +20,7 @@ subroutine readinput integer :: stat integer :: i, ia, n - integer :: j + integer :: j, io integer :: NN integer :: nwann real(dp) :: t1, temp @@ -41,10 +41,6 @@ subroutine readinput stop endif - !> inout file - - - !read(1001,*)Hrfile read(1001, TB_FILE, iostat= stat) if (stat/=0) then @@ -322,6 +318,14 @@ subroutine readinput stop "ERROR: please set projectors for Wannier functions information" endif + !> set up orbitals_start + allocate(orbitals_start(Num_atoms)) + orbitals_start= 1 + do i=1, Num_atoms-1 + orbitals_start(i+1)= orbitals_start(i)+ nprojs(i) + enddo + + !> read Wannier centres nwann= sum(nprojs) if (SOC>0) nwann= 2*nwann @@ -794,7 +798,84 @@ subroutine readinput k1=k_mass call direct_cart_rec(k1, k_mass) if (cpuid==0) write(stdout, '(a, 3f7.4, a)')'>> k points ', k_mass, " in unit of 1/Angstrom" - + + !> setup the atoms on top and bottom surface that used for output the + !> surface-state spectrum + !> by default we output all the atoms' weight + NtopAtoms = Num_atoms + NbottomAtoms= Num_atoms + allocate(TopAtoms(NtopAtoms)) + allocate(BottomAtoms(NbottomAtoms)) + do i=1, NTopAtoms + TopAtoms(i)= i + enddo + do i=1, NBottomAtoms + BottomAtoms(i)= i + enddo + + rewind(1001) + lfound = .false. + do while (.true.) + read(1001, *, end= 112)inline + if (trim(adjustl(inline))=='TOPBOTTOMATOMS') then + lfound= .true. + if (cpuid==0) write(stdout, *)' ' + if (cpuid==0) write(stdout, *)'We found TOPBOTTOMATOMS card' + exit + endif + enddo + read(1001, *)NtopAtoms + deallocate(TopAtoms) + allocate(TopAtoms(NtopAtoms)) + read(1001, *)TopAtoms + read(1001, *)NbottomAtoms + deallocate(BottomAtoms) + allocate(BottomAtoms(NbottomAtoms)) + read(1001, *)BottomAtoms + + 112 continue + if (.not.lfound.and.cpuid==0)write(stdout, *)'>> Output all atoms weight for surface state spectrum' + if (cpuid==0) write(stdout, *)'> NtopAtoms ', NtopAtoms + if (cpuid==0) write(stdout, '(a,100i3)')'> TopAtoms ', TopAtoms + if (cpuid==0) write(stdout, '(a,100i3)')'> NbottomAtoms ', NbottomAtoms + if (cpuid==0) write(stdout, '(a,100i3)')'> BottomAtoms ',BottomAtoms + + NtopOrbitals=0 + do i=1, NTopAtoms + NtopOrbitals= NtopOrbitals+ nprojs(TopAtoms(i)) + enddo + if (SOC>0) NtopOrbitals= NtopOrbitals*2 + allocate(TopOrbitals(NtopOrbitals)) + TopOrbitals= 1 + + !> set up top surface orbitals for output the surface spectrum + io=0 + do i=1, NTopAtoms + do j=1, nprojs(TopAtoms(i)) + io =io+ 1 + TopOrbitals(io)= orbitals_start(i)+ j- 1 + if (SOC>0)TopOrbitals(io+ NtopOrbitals/2 )= orbitals_start(i)+ j- 1+ Num_wann/2 + enddo ! j + enddo ! i + + NBottomOrbitals=0 + do i=1, NBottomAtoms + NBottomOrbitals= NBottomOrbitals+ nprojs(BottomAtoms(i)) + enddo + if (SOC>0) NBottomOrbitals= NBottomOrbitals*2 + allocate(BottomOrbitals(NBottomOrbitals)) + BottomOrbitals= 1 + + !> set up Bottom surface orbitals for output the surface spectrum + io=0 + do i=1, NBottomAtoms + do j=1, nprojs(BottomAtoms(i)) + io =io+ 1 + BottomOrbitals(io)= orbitals_start(i)+ j- 1 + if (SOC>0)BottomOrbitals(io+ NBottomOrbitals/2)= orbitals_start(i)+ j- 1+ Num_wann/2 + enddo ! j + enddo ! i + !> close input.dat close(1001) diff --git a/soc/surfstat.f90 b/soc/surfstat.f90 index 64fb410c..c0f647be 100644 --- a/soc/surfstat.f90 +++ b/soc/surfstat.f90 @@ -18,7 +18,7 @@ subroutine surfstat integer :: ierr ! general loop index - integer :: i,j + integer :: i, j, io ! kpoint loop index integer :: ikp @@ -115,9 +115,13 @@ subroutine surfstat !call surfgreen_1984(w,GLL,GRR,H00,H01,ones) ! calculate spectral function - do i= 1, ndim - dos_l(ikp, j)=dos_l(ikp,j)- aimag(GLL(i,i)) - dos_r(ikp, j)=dos_r(ikp,j)- aimag(GRR(i,i)) + do i= 1, NtopOrbitals + io= TopOrbitals(i) + dos_l(ikp, j)=dos_l(ikp,j)- aimag(GLL(io,io)) + enddo ! i + do i= 1, NBottomOrbitals + io= Ndim- Num_wann+ BottomOrbitals(i) + dos_r(ikp, j)=dos_r(ikp,j)- aimag(GRR(io,io)) enddo ! i enddo ! j enddo ! ikp