added DOS_SLAB_CALC tags to calculate the density of states of slab s…

…ystems

examples/Bi2Se3/wt.in

@@ -28,16 +28,20 @@ Se pz px py
 
 !> bulk band structure calculation flag
 &CONTROL
-BulkBand_calc         = T
+BulkBand_calc         = F
 Dos_calc              = F
 BulkBand_points_calc  = F
 SlabSS_calc           = T
+
 /
 
 &SYSTEM
 NumOccupied = 18        ! NumOccupied
 SOC = 1                 ! soc
 E_FERMI = 4.4195        ! e-fermi
+!Add_Zeeman_Field = T
+!Zeeman_energy_in_eV = 0.5
+!Nslab = 10
 /
 
 ! get projected bands onto different orbitals, here we only consider orbital and omit the spin freedom

src/Makefile

@@ -14,7 +14,7 @@ OBJ =  module.o sparse.o wt_aux.o math_lib.o symmetry.o \
 		main.o
 
 # compiler
-F90  = mpif90 -fpp -DMPI -fpe3 -O3 -DARPACK -DINTELMKL # -check all -traceback -g
+F90  = mpif90 -fpp -DMPI -fpe3 -O3 -DARPACK -DINTELMKL #-check all -traceback -g
 #F90  = ifort -fpp  -DINTELMKL -fpe3 -check all -traceback -g
 CC = mpicc -cpp -O3 -DINTELMKL 
 

src/dos.f90

@@ -352,7 +352,6 @@ subroutine charge_density_sparse
 end subroutine charge_density_sparse
 
 
-
 subroutine dos_sub
 !> calculate density of state for 3D bulk system
 !
@@ -446,7 +445,12 @@ subroutine dos_sub
          + K3D_vec2_cube*(iky-1)/dble(nk2)  &
          + K3D_vec3_cube*(ikz-1)/dble(nk3)
 
-      call ham_bulk_atomicgauge(k, Hk)
+      if (index(KPorTB, 'KP')/=0)then
+         call ham_bulk_kp_abcb_graphene(k, Hk)
+      else
+         call ham_bulk_atomicgauge(k, Hk)
+      endif
+
       W= 0d0
       call eigensystem_c( 'N', 'U', Num_wann ,Hk, W)
       eigval(:)= W(iband_low:iband_high)
@@ -531,6 +535,186 @@ subroutine dos_sub
    return
 end subroutine dos_sub
 
+
+
+subroutine dos_slab
+!> calculate density of state for slab system
+!
+!> DOS(\omega)= \sum_k \delta(\omega- E(k))
+
+   use wmpi
+   use para
+   implicit none
+
+   !> the integration k space
+   real(dp) :: emin, emax
+
+   integer :: ik,ie,ib,ikx,iky,ikz,ieta
+   integer :: knv2_slab,NE,ierr, ndim_slab
+   integer :: NumberofEta
+
+   !> integration for band
+   integer :: iband_low,iband_high,iband_tot
+
+   real(dp) :: x, dk2, eta0
+
+   real(dp) :: k(2)
+   real(dp) :: time_start, time_end
+
+   real(dp), allocatable :: eigval(:)
+   real(dp), allocatable :: W(:)
+   real(dp), allocatable :: omega(:)
+   real(dp), allocatable :: dos(:, :)
+   real(dp), allocatable :: dos_mpi(:, :)
+   real(dp), allocatable :: eta_array(:)
+   complex(dp), allocatable :: Hk(:, :)
+
+   !> delta function
+   real(dp), external :: delta
+
+   knv2_slab= Nk1*Nk2
+   ndim_slab= Num_wann*Nslab
+
+   if (OmegaNum<2) OmegaNum=2
+   NE= OmegaNum
+
+   iband_low= Numoccupied- 10000
+   iband_high= Numoccupied+ 10000
+
+   if (iband_low <1) iband_low = 1
+   if (iband_high >ndim_slab) iband_high = ndim_slab
+
+   iband_tot= iband_high- iband_low+ 1
+
+   NumberofEta = 9
+
+   allocate(W(ndim_slab))
+   allocate(Hk(ndim_slab, ndim_slab))
+   allocate(eigval(iband_tot))
+   allocate(eta_array(NumberofEta))
+   allocate(dos(NE, NumberofEta))
+   allocate(dos_mpi(NE, NumberofEta))
+   allocate(omega(NE))
+   dos=0d0
+   dos_mpi=0d0
+   eigval= 0d0
+
+
+   emin= OmegaMin
+   emax= OmegaMax
+   eta_array=(/0.1d0, 0.2d0, 0.4d0, 0.8d0, 1.0d0, 2d0, 4d0, 8d0, 10d0/)
+   eta_array= eta_array*Fermi_broadening
+
+
+   !> energy
+   do ie=1, NE
+      omega(ie)= emin+ (emax-emin)* (ie-1d0)/dble(NE-1)
+   enddo ! ie
+
+   !dk2= kCubeVolume/dble(knv2_slab)
+   dk2= 1d0/dble(knv2_slab)
+
+   !> get eigenvalue
+   time_start= 0d0
+   time_end= 0d0
+   do ik=1+cpuid, knv2_slab, num_cpu
+
+      if (cpuid.eq.0.and. mod(ik/num_cpu, 100).eq.0) &
+         write(stdout, '(a, i18, "/", i18, a, f10.3, "s")') 'ik/knv2_slab', &
+         ik, knv2_slab, ' time left', (knv2_slab-ik)*(time_end-time_start)/num_cpu
+
+      call now(time_start)
+      ikx= (ik-1)/(Nk2)+1
+      iky= ((ik-1-(ikx-1)*Nk2))+1
+      k= K2D_start+ K2D_vec1*(ikx-1)/dble(Nk1)  &
+         + K2D_vec2*(iky-1)/dble(Nk2) 
+
+      !> get the Hamlitonian
+      call ham_slab(k, Hk)
+
+      W= 0d0
+      call eigensystem_c('N', 'U', ndim_slab ,Hk, W)
+      eigval(:)= W(iband_low:iband_high)
+
+      !> get density of state
+      do ie= 1, NE
+         do ib= 1, iband_tot
+            x= omega(ie)- eigval(ib)
+            do ieta= 1, NumberofEta
+               eta0= eta_array(ieta)
+               dos_mpi(ie, ieta) = dos_mpi(ie, ieta)+ delta(eta0, x)
+            enddo
+         enddo ! ib
+      enddo ! ie
+      call now(time_end)
+
+      call now(time_end)
+
+   enddo  ! ik
+
+#if defined (MPI)
+   call mpi_allreduce(dos_mpi,dos,size(dos),&
+      mpi_dp,mpi_sum,mpi_cmw,ierr)
+#else
+   dos= dos_mpi
+#endif
+   dos= dos*dk2
+
+   !> include the spin degeneracy if there is no SOC in the tight binding Hamiltonian.
+   if (SOC<=0) dos=dos*2d0
+
+   outfileindex= outfileindex+ 1
+   if (cpuid.eq.0) then
+      open(unit=outfileindex, file='dos_slab.dat')
+      write(outfileindex, "(a, i10)")'# Density of state of slab system, Nslab= ', Nslab
+      write(outfileindex, '(2a16)')'# E(eV)', 'DOS(E) (1/eV)'
+      write(outfileindex, '("#", a, f6.2, 300f16.2)')'Broadening \eta (meV): ', Eta_array(:)*1000d0/eV2Hartree
+      do ie=1, NE
+         write(outfileindex, '(90f16.6)')omega(ie)/eV2Hartree, dos(ie, :)*eV2Hartree
+      enddo ! ie
+      close(outfileindex)
+   endif
+
+   outfileindex= outfileindex+ 1
+   !> write script for gnuplot
+   if (cpuid==0) then
+      open(unit=outfileindex, file='dos_slab.gnu')
+      write(outfileindex, '(a)')"set encoding iso_8859_1"
+      write(outfileindex, '(a)')'set terminal pdf enhanced color font ",16" size 5,4 '
+      write(outfileindex, '(a)')"set output 'dos_slab.pdf'"
+      write(outfileindex, '(a)')'set border lw 2'
+      write(outfileindex, '(a)')'set autoscale fix'
+      write(outfileindex, '(a, f16.6,a)')'set yrange [0:', maxval(dos)*eV2Hartree+0.5, '1]'
+      write(outfileindex, '(a)')'set key samplen 0.8 spacing 1 font ",12"'
+      write(outfileindex, '(a)')'set xlabel "Energy (eV)"'
+      write(outfileindex, '(a)')'set key title  "Broadening"'
+      write(outfileindex, '(a)')'set title "DOS with different broadenings"'
+      write(outfileindex, '(a)')'set ylabel "DOS (states/eV/unit cell)"'
+      write(outfileindex, '(a, f6.1, a)')"plot 'dos_slab.dat' u 1:2 w l lw 2 title '",&
+         Eta_array(1)*1000/eV2Hartree, "meV', \"
+      do ieta= 2, NumberofEta-1
+         write(outfileindex, 202)" '' u 1:", ieta, " w l lw 2 title '", &
+            Eta_array(ieta)*1000/eV2Hartree, "meV', \"
+      enddo
+      write(outfileindex, '(a, f6.1, a)')" '' u 1:10 w l lw 2 title '",&
+         Eta_array(NumberofEta)*1000/eV2Hartree, "meV'"
+      close(outfileindex)
+   endif
+202 format(a, i3, a, f6.1, a)
+
+#if defined (MPI)
+   call mpi_barrier(mpi_cmw, ierr)
+#endif
+   deallocate(W)
+   deallocate(Hk)
+   deallocate(eigval)
+   deallocate(dos)
+   deallocate(dos_mpi)
+   deallocate(omega)
+
+   return
+end subroutine dos_slab
+
 subroutine joint_dos
 ! calculate joint density of state for 3D bulk system
 !

src/ek_bulk.f90

@@ -49,7 +49,7 @@ subroutine ek_bulk_line
 
       ! generate bulk Hamiltonian
       if (index(KPorTB, 'KP')/=0)then
-         call ham_bulk_kp(k, Hamk_bulk)
+         call ham_bulk_kp_abcb_graphene(k, Hamk_bulk)
       else
          !> deal with phonon system
          if (index(Particle,'phonon')/=0.and.LOTO_correction) then
@@ -870,7 +870,7 @@ subroutine ek_bulk_plane
       ! generate bulk Hamiltonian
       Hamk_bulk= 0d0
       if (index(KPorTB, 'KP')/=0)then
-         call ham_bulk_kp(k, Hamk_bulk)
+         call ham_bulk_kp_abcb_graphene(k, Hamk_bulk)
       else
          !> deal with phonon system
          if (index(Particle,'phonon')/=0.and.LOTO_correction) then

src/ham_bulk.f90

@@ -580,6 +580,119 @@ subroutine ham_bulk_LOTO(k,Hamk_bulk)
    return
 end subroutine ham_bulk_LOTO
 
+subroutine ham_bulk_kp_abcb_graphene(k,Hamk_bulk)
+   ! > construct the kp model at K point of ABCB tetralayer graphene
+
+   use para, only : Num_wann, dp, stdout, zi, zzero, One_complex, &
+      Symmetrical_Electric_field_in_eVpA, eV2Hartree, Angstrom2atomic, &
+      Electric_field_in_eVpA, Origin_cell, cpuid
+   implicit none
+
+   integer :: i1,i2,ia,ib,ic,iR, nwann, io, i
+   integer :: add_electric_field
+   real(dp) :: pos(Origin_cell%Num_atoms)
+   real(dp) :: static_potential
+
+   ! coordinates of R vector
+   real(Dp) :: R(3), R1(3), R2(3), kdotr, kx, ky, kz
+   real(dp) :: v,v3,v4,v6,v7,d1,d2,d3,d4,d5,d6,d7,d8,g1,g2,g3,g4,g5,g6
+
+   complex(dp) :: factor, km, kp
+
+   real(dp), intent(in) :: k(3)
+
+   !> the k point that the kp model is constructed
+   !> in fractional coordinate
+   real(dp) :: k_center_direct(3), l(19)
+
+   ! Hamiltonian of bulk system
+   complex(Dp),intent(out) ::Hamk_bulk(Num_wann, Num_wann)
+
+   if (Num_wann/=8) then
+      print *, "Error : in this kp model, num_wann should be 8"
+      print *, 'Num_wann', Num_wann
+      stop
+   endif
+  !l=(/ 1.26416117e+01, -5.57664351e+00, -5.87940524e+00,  1.60839693e+01, &
+  ! 8.13551651e+00,  2.82753440e-01,  9.16947940e-03, -2.61922210e-03, &
+  !-2.77215942e-03,  7.01879554e-03, -3.38271168e-01, -6.64427749e-02, &
+  !-3.08356540e-02, -8.82150889e-03, -2.57865658e-03,  1.83192652e-01, &
+  ! 4.35695252e-03, -2.28829893e-01, -1.82821963e-02/)
+
+   l=(/1.57937270e+01, -2.49047067e+00, -1.82376120e+00,  1.21854647e-02, &
+    -5.75783391e+00, -2.09787893e-01,  7.48348077e-03,  2.03490675e-01, &
+    -1.82294202e-04,  3.11506594e-02, -1.32325756e-01,  2.72568864e-01, &
+     2.27263042e-02, -2.57806046e-01, -2.01920757e-01, -2.46168084e-01, &
+    -1.45799538e-01, -1.59225389e-01, -2.20335989e-01/)
+
+   v=l(1); v3=l(2); v4=l(3); v6=l(4); v7=l(5);
+   d1=l(6); d2= l(7); d3= l(8); d4=l(9); d5=l(10); d6= l(11); d7=l(12); d8= l(19)
+   g1= l(13); g2= l(14); g3= l(15); g4= l(16); g5= l(17); g6= l(18);
+
+   k_center_direct= (/1d0/3d0, 1d0/3d0, 0d0/)
+   kx=k(1) -k_center_direct(1)
+   ky=k(2) -k_center_direct(2)
+   kz=k(3) -k_center_direct(3)
+
+   km=-kx+ zi*ky
+   kp=-kx- zi*ky
+   !> write out the parameters
+   if (cpuid.eq.0) then
+   !  write(stdout, '(a, 8f8.4)') '> Onsite energies: ', d1, d2, d3, d4, d5, d6, d7, d8
+   endif
+
+   Hamk_bulk= 0d0
+   !> kp
+   Hamk_bulk(:, 1)= (/d1*One_complex, v*kp, -v4*kp, v3*km, v6*km, g6*One_complex, zzero, zzero/)
+   Hamk_bulk(:, 2)= (/v*km, d2*One_complex, g1*One_complex, -v4*kp, v7*kp, v6*km, zzero, zzero/)
+   Hamk_bulk(:, 3)= (/-v4*km, g1*One_complex, d3*One_complex, v*kp, -v4*kp, v3*km, g4*One_complex, zzero/)
+   Hamk_bulk(:, 4)= (/v3*kp, -v4*km, v*km, d4*One_complex, g2*One_complex, -v4*kp, zzero, g5*One_complex/)
+   Hamk_bulk(:, 5)= (/v6*kp, v7*km, -v4*km, g2*One_complex, d5*One_complex, v*kp, -v4*km, g3*One_complex/)
+   Hamk_bulk(:, 6)= (/g6*One_complex, v6*kp, v3*kp, -v4*km, v*km, d6*One_complex, v3*kp, -v4*km/)
+   Hamk_bulk(:, 7)= (/zzero, zzero, g4*One_complex, zzero, -v4*kp, v3*km, d7*One_complex, v*kp/)
+   Hamk_bulk(:, 8)= (/zzero, zzero, zzero, g5*One_complex, g3*One_complex, -v4*kp, v*km, d8*One_complex/)
+
+   !> add electrical field
+
+   call get_stacking_direction_and_pos(add_electric_field, pos)
+
+   if (add_electric_field>0) then
+      io=0
+      do ia=1, Origin_cell%Num_atoms
+         !static_potential= pos(ia)*Origin_cell%cell_parameters(add_electric_field)*Electric_field_in_eVpA
+         !if (Inner_symmetrical_Electric_Field) then
+         !   static_potential= abs(pos(ia)*Origin_cell%cell_parameters(add_electric_field))*Electric_field_in_eVpA
+         !endif
+         static_potential= abs(pos(ia)*Origin_cell%cell_parameters(add_electric_field))&
+            *Symmetrical_Electric_field_in_eVpA*eV2Hartree/Angstrom2atomic+ &
+            (pos(ia)*Origin_cell%cell_parameters(add_electric_field))*&
+            Electric_field_in_eVpA*eV2Hartree/Angstrom2atomic
+
+         do i=1, Origin_cell%nprojs(ia)
+            io=io+1
+            Hamk_bulk(io, io)= Hamk_bulk(io, io)+ static_potential
+         enddo ! nproj
+      enddo ! ia
+   endif  ! add electric field or not
+
+   ! check hermitcity
+   do i1=1, Num_wann
+      do i2=1, Num_wann
+         if(abs(Hamk_bulk(i1,i2)-conjg(Hamk_bulk(i2,i1))).ge.1e-6)then
+            write(stdout,*)'there is something wrong with Hamk_bulk'
+            write(stdout,*)'i1, i2', i1, i2
+            write(stdout,*)'value at (i1, i2)', Hamk_bulk(i1, i2)
+            write(stdout,*)'value at (i2, i1)', Hamk_bulk(i2, i1)
+            !stop
+         endif
+      enddo
+   enddo
+   Hamk_bulk= Hamk_bulk*eV2Hartree
+
+   return
+end subroutine ham_bulk_kp_abcb_graphene
+
+
 subroutine ham_bulk_kp(k,Hamk_bulk)
    ! > construct the kp model at K point of WC system
    ! > space group is 187. The little group is C3h