you can specify the atoms on the surface to output the surface spectrum

quanshengwu · Sep 26, 2016 · fff25c8 · fff25c8
1 parent c357a1f
commit fff25c8
Show file tree

Hide file tree

Showing 4 changed files with 165 additions and 30 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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