From bf576c4b9603bce75f4652b9742a015eb8e8ce30 Mon Sep 17 00:00:00 2001 From: QuanSheng Wu Date: Sat, 9 Nov 2024 09:44:28 +0800 Subject: [PATCH] added DOS_SLAB_CALC tags to calculate the density of states of slab systems --- examples/Bi2Se3/wt.in | 6 +- src/Makefile | 2 +- src/dos.f90 | 188 +++++++++++++++++- src/ek_bulk.f90 | 4 +- src/ham_bulk.f90 | 113 +++++++++++ src/landau_level.f90 | 176 +++++++++++++++- src/main.f90 | 10 + src/module.f90 | 15 +- src/readinput.f90 | 1 + src/surfstat.f90 | 2 +- .../Makefile | 2 +- 11 files changed, 506 insertions(+), 13 deletions(-) diff --git a/examples/Bi2Se3/wt.in b/examples/Bi2Se3/wt.in index add9757b..228fa9d6 100644 --- a/examples/Bi2Se3/wt.in +++ b/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 diff --git a/src/Makefile b/src/Makefile index 07eba542..943d4229 100644 --- a/src/Makefile +++ b/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 diff --git a/src/dos.f90 b/src/dos.f90 index e029898e..9d6e4d6d 100644 --- a/src/dos.f90 +++ b/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 ! diff --git a/src/ek_bulk.f90 b/src/ek_bulk.f90 index 5217dbdc..fedf4672 100644 --- a/src/ek_bulk.f90 +++ b/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 diff --git a/src/ham_bulk.f90 b/src/ham_bulk.f90 index 03fd510f..9ae0c57f 100644 --- a/src/ham_bulk.f90 +++ b/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 diff --git a/src/landau_level.f90 b/src/landau_level.f90 index 48c01248..3f036be2 100644 --- a/src/landau_level.f90 +++ b/src/landau_level.f90 @@ -637,13 +637,14 @@ subroutine landau_level_B integer :: Nq, Ndimq integer :: Nmag - integer :: ib, i, j, ie, ierr, iq, ig + integer :: ib, i, j, ie, ierr, iq, ig, ieta - real(dp) :: B0, B0Tesla, B0Tesla_quantumflux_magsupcell + real(dp) :: B0, B0Tesla, B0Tesla_quantumflux_magsupcell, x, eta real(dp) :: time_start, time_end, time_start0 ! wave vector real(dp) :: k3(3) + real(dp), external :: delta !> dim= Ndimq, knv3 real(dp), allocatable :: mag(:), mag_Tesla(:) @@ -659,6 +660,18 @@ subroutine landau_level_B real(dp), allocatable :: dos_selected(:, :, :) real(dp), allocatable :: dos_selected_mpi(:, :, :) + !> energy interval + !> OmegaNum is defined in the module.f90 and read from the input.dat or wt.in + real, allocatable :: omega(:) + + !> spectrum calculated + real(dp), allocatable :: dos_B_omega(:, :, :), dos_B_omega_mpi(:, :, :) + real(dp), allocatable :: n_int(:, :) + + integer :: NumberofEta, ie_Earc + real(dp), allocatable :: eta_array(:), n_Earc(:) + + !> dim= Ndimq*Ndimq complex(dp), allocatable :: ham_landau(:, :) @@ -683,6 +696,22 @@ subroutine landau_level_B eigv = 0d0 ham_landau= 0d0 + NumberofEta=9 + allocate(eta_array(NumberofEta)) + allocate(n_Earc(NumberofEta)) + eta_array=(/0.1d0, 0.2d0, 0.4d0, 0.8d0, 1.0d0, 2d0, 4d0, 8d0, 10d0/) + eta_array= eta_array*Fermi_broadening + + allocate(omega(OmegaNum)) + allocate(n_int(0:Omeganum, NumberofEta)) + allocate(dos_B_omega(Nmag, OmegaNum, NumberofEta), dos_B_omega_mpi(Nmag, OmegaNum, NumberofEta)) + dos_B_omega= 0d0; dos_B_omega_mpi= 0d0 + + !> energy + do ie=1, OmegaNum + omega(ie)= OmegaMin+ (OmegaMax-OmegaMin)* (ie-1d0)/dble(OmegaNum-1) + enddo ! ie + !> deal with the magnetic field !> first transform the Bx By into B*Cos\theta, B*Sin\theta !~ if (abs(By)<1e-8) then @@ -751,7 +780,6 @@ subroutine landau_level_B call ham_3Dlandau(Ndimq, Nq, k3, ham_landau) - ham_landau= ham_landau/eV2Hartree !> diagonalization by call zheev in lapack W= 0d0 @@ -777,6 +805,18 @@ subroutine landau_level_B eigv(:, ib)= W + !> calculate density of states using the eigenvalues + do ieta=1, NumberofEta + eta= eta_array(ieta) + do ie=1, OmegaNum + do i = 1, Ndimq + x= omega(ie)- W(i) + dos_B_omega(ib, ie, ieta)= dos_B_omega(ib, ie, ieta)+ & + delta(eta, x) + enddo + enddo + enddo + call now(time_end) enddo !ib @@ -785,10 +825,14 @@ subroutine landau_level_B mpi_dp,mpi_sum,mpi_cmw,ierr) call mpi_allreduce(dos_selected, dos_selected_mpi,size(dos_selected),& mpi_dp,mpi_sum,mpi_cmw,ierr) + call mpi_allreduce(dos_B_omega, dos_B_omega_mpi, size(dos_B_omega), & + mpi_dp, mpi_sum, mpi_cmw, ierr) #else eigv_mpi= eigv dos_selected_mpi= dos_selected + dos_B_omega_mpi= dos_B_omega #endif + eigv_mpi = eigv_mpi/eV2Hartree outfileindex= outfileindex+ 1 if (cpuid.eq.0) then @@ -852,6 +896,132 @@ subroutine landau_level_B close(outfileindex) endif + outfileindex= outfileindex+ 1 + if (cpuid.eq.0)then + open (unit=outfileindex, file='LandauLevel_B_dos.dat') + open (unit=outfileindex+1,file='wannierdiagram.dat') + write(outfileindex, '("#", a)')'Hofstadter butterfly (Landau level in (E,B) axis) ' + write(outfileindex, '("#", a, i6)')'Magnetic supercell size : ', Magq + write(outfileindex+1, '("#", a)')'Wannier diagram ' + write(outfileindex+1, '("#", a, i6)')'Magnetic supercell size : ', Magq + write(outfileindex, '("#", a, 3f12.6, a)')'k points in fractional coordinates (', k3, ')' + write(outfileindex, '("# Column", I5, 100I16)')(i, i=1, 12) + write(outfileindex, '("#", a15, 5a16)', advance='NO')'Flux', 'B(Tesla)', 'E(eV)' + write(outfileindex+1, '("# Column", I5, 100I16)')(i, i=1, 20) + write(outfileindex+1, '("#", a15, 2a16)', advance='NO')'Flux', 'B(Tesla)' + do ieta=1, NumberofEta-1 + write(outfileindex, '(a16)', advance='NO') 'A(B, E)' + write(outfileindex+1, '(2a16)', advance='NO') 'n ', 'A(n, E)' + enddo + write(outfileindex, '(a16)') 'A(B, E)' + write(outfileindex+1, '(2a16)') 'n ', 'A(n, E)' + + write(outfileindex, '("#", a, 20X, 300f16.2)')'Broadening \eta (meV): ', Eta_array(:)*1000d0/eV2Hartree + write(outfileindex+1, '("#", a, 5X, f26.2,300f32.2)')'Broadening \eta (meV): ', Eta_array(:)*1000d0/eV2Hartree + + !> find n_int at EF + do ie=1, omeganum + if (omega(ie)>iso_energy) then + ie_Earc= ie- 1 + exit + endif + enddo + + do ib=1, Nmag + n_int=0 + do ie=1, omeganum + write(outfileindex, '(300E16.4)')mag(ib)/2d0/pi, mag_Tesla(ib), omega(ie)/eV2Hartree, dos_B_omega_mpi(ib, ie, :) + enddo + + !> get number of electrons between the lowest energy level and iso_energy + n_Earc= 0d0 + do ie=1, ie_Earc + n_Earc(:)= n_Earc(:)+ dos_B_omega_mpi(ib, ie, :) + enddo + + !> get number of electrons between the lowest energy level and omega(ie) + do ie=1, omeganum + n_int(ie, :)=n_int(ie-1, :)+ dos_B_omega_mpi(ib, ie, :) + enddo + + !> set n(E)= n_int- n_Earc= \int_Earc^E \rho(\epsilon)d\epsilon + do ie=1, omeganum + n_int(ie, :)= n_int(ie, :)-n_Earc(:) + enddo + + !do ieta=1, NumberofEta + ! n_int(:, ieta)=n_int(:, ieta)/n_int(omeganum, ieta) + !enddo + + do ie=1, omeganum + write(outfileindex+1,'(300E16.4)') mag(ib)/2d0/pi, mag_Tesla(ib), & + (n_int(ie, ieta), dos_B_omega_mpi(ib, ie, ieta), ieta=1, NumberofEta) + enddo + + write(outfileindex+1, *) ' ' + write(outfileindex, *) ' ' + enddo + close(outfileindex) + write(stdout,*)'calculate Landau level spectrum in B-E mode successfully' + endif + + outfileindex= outfileindex+ 2 + if (cpuid.eq.0)then + open (unit=outfileindex, file='LandauLevel_B_dos.gnu') + write(outfileindex, '(a)')'#set terminal postscript enhanced color font ",24"' + write(outfileindex, '(a)')'set terminal pngcairo enhanced color font ",60" size 1920, 1680' + write(outfileindex, '(a)')"set output 'LandauLevel_B_dos.png'" + write(outfileindex, '(a)')'set pm3d' + write(outfileindex, '(a)')'set palette rgb 21,22,23' + write(outfileindex, '(a)')'#set isosamples 50,50' + write(outfileindex, '(a)')'set size 0.9,1' + write(outfileindex, '(a)')'set origin 0.05,0' + write(outfileindex, '(a)')'set view map' + write(outfileindex, '(a)')'unset ztics' + write(outfileindex, '(a)')'unset surface' + write(outfileindex, '(a)')'unset key' + write(outfileindex, '(a)')'#set xtics font ",24"' + write(outfileindex, '(a)')'#set ytics font ",24"' + write(outfileindex, '(a)')'#set xlabel font ",24"' + write(outfileindex, '(a)')'set ylabel "Energy (eV)"' + write(outfileindex, '(a, i6,a)')'set title "Hofstadter butterfly with Nq=', Nq, '" font ",40"' + write(outfileindex, '(a, f12.6, a, f12.6, a)')'set yrange [',OmegaMin/eV2Hartree,':',OmegaMax/eV2Hartree,']' + write(outfileindex, '(a)')'set xlabel "Phi per unit cell"' + write(outfileindex,*) '#set xlabel "{/Symbol F}/{/Symbol F}_0"' + write(outfileindex, '(a, f12.6, a)') 'set xrange [ 0.0000 :', maxval(mag/2d0/pi) ,']' + write(outfileindex, '(a)') "splot 'LandauLevel_B_dos.dat' u 1:3:(log($4)) w pm3d" + write(outfileindex, '(a)')'set xlabel "B (Tesla)"' + write(outfileindex, '(a, f12.6, a)') '#set xrange [ 0.0000 :', maxval(mag_Tesla) ,']' + write(outfileindex, '(a)') "#splot 'LandauLevel_B_dos.dat' u 2:3:(log($4)) w pm3d" + endif + + outfileindex=outfileindex+1 + if(cpuid == 0) then + open(unit=outfileindex, file='wannierdiagram.gnu') + write(outfileindex,*) 'set terminal pngcairo enhanced color font ",60" size 1920, 1680' + write(outfileindex,*) "set output 'wannierdiagram.png'" + write(outfileindex,*) 'set pm3d' + write(outfileindex, '(a)')'set size 0.9,1' + write(outfileindex, '(a)')'set origin 0.05,0' + write(outfileindex, '(a)')'set palette rgb 21,22,23' + write(outfileindex,*) '#set isosamples 50,50' + write(outfileindex,*) 'set view map' + write(outfileindex,*) 'unset ztics' + write(outfileindex,*) 'unset surface' + write(outfileindex,*) 'unset key' + + write(outfileindex,*) 'set ylabel "n"' + write(outfileindex, '(a, i6, a)') 'set title "Wannier diagram with Nq=', Nq, '"' + write(outfileindex,*) '#set yrange [ ] noextend' + write(outfileindex,*) '#set xlabel "{/Symbol F}/{/Symbol F}_0"' + write(outfileindex,*) 'set xlabel "Phi/Phi_0 per unit cell"' + write(outfileindex,*) '#set xlabel "{/Symbol F}/{/Symbol F}_0"' + write(outfileindex,*) '#set xrange [ ] noextend' + + write(outfileindex,*) "splot 'wannierdiagram.dat' u 1:3:(log($4)) w pm3d #lc palette" + end if + + #if defined (MPI) call mpi_barrier(mpi_cmw, ierr) #endif diff --git a/src/main.f90 b/src/main.f90 index 3a4cce96..d56559d8 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -358,6 +358,16 @@ program main if(cpuid.eq.0)write(stdout, *)'<< End of calculating the bulk FS' endif + if (Dos_slab_calc) then + if(cpuid.eq.0)write(stdout, *)' ' + if(cpuid.eq.0)write(stdout, *)'>> Start of calculating DOS for slab system' + call now(time_start) + call dos_slab + call now(time_end) + call print_time_cost(time_start, time_end, 'Dos_calc and Jdos_calc') + if(cpuid.eq.0)write(stdout, *)'<< End of calculating the DOS for slab system' + endif + !> calculate density of state and joint density of state if (JDos_calc.and.Dos_calc) then if(cpuid.eq.0)write(stdout, *)' ' diff --git a/src/module.f90 b/src/module.f90 index a8793ac0..83692009 100644 --- a/src/module.f90 +++ b/src/module.f90 @@ -393,6 +393,7 @@ module para logical :: WireBand_calc ! Flag for 1D wire energy band calculation logical :: SlabSS_calc ! Flag for surface state ARPES spectrum calculation logical :: Dos_calc ! Flag for density of state calculation + logical :: Dos_slab_calc ! Flag for density of state calculation for slab system logical :: ChargeDensity_selected_bands_calc ! Flag for charge density logical :: ChargeDensity_selected_energies_calc ! Flag for charge density logical :: JDos_calc ! Flag for joint density of state calculation @@ -476,7 +477,7 @@ module para BerryCurvature_plane_selectedbands_calc, & BerryCurvature_slab_calc, MirrorChern_calc, BerryCurvature_Cube_calc, & Z2_3D_calc, Chern_3D_calc, WeylChirality_calc, NLChirality_calc, & - Dos_calc, JDos_calc, EffectiveMass_calc, & + Dos_calc, Dos_slab_calc, JDos_calc, EffectiveMass_calc, & FindNodes_calc, TBtoKP_calc, LOTO_correction, & BulkBand_plane_calc, Hof_Butt_calc, Symmetry_Import_calc, & Boltz_Berry_correction, & @@ -1175,7 +1176,17 @@ module element_table real, parameter :: a_0 = 5.29177210903E-1 ! Bohr radius in the unit of angstrom - character(len=10) :: element_name(lenth_of_table) = ['H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf'] + character(len=10) :: element_name(lenth_of_table) = & + ['H','He','Li','Be','B','C','N','O','F',& + 'Ne','Na','Mg','Al','Si','P','S','Cl','Ar',& + 'K','Ca','Sc','Ti','V','Cr','Mn','Fe','Co','Ni',& + 'Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr',& + 'Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In',& + 'Sn','Sb','Te','I','Xe','Cs','Ba','La','Ce','Pr','Nd',& + 'Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu',& + 'Hf','Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb',& + 'Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa','U','Np',& + 'Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf'] ! Each element's electron configuration character(len=64) :: element_electron_config(lenth_of_table) = [& diff --git a/src/readinput.f90 b/src/readinput.f90 index d0c05784..9c4762ae 100644 --- a/src/readinput.f90 +++ b/src/readinput.f90 @@ -129,6 +129,7 @@ subroutine readinput BerryCurvature_kpath_Occupied_calc = .FALSE. MirrorChern_calc = .FALSE. Dos_calc = .FALSE. + Dos_slab_calc = .FALSE. JDos_calc = .FALSE. EffectiveMass_calc = .FALSE. FindNodes_calc = .FALSE. diff --git a/src/surfstat.f90 b/src/surfstat.f90 index 064c33fc..438ecb7a 100644 --- a/src/surfstat.f90 +++ b/src/surfstat.f90 @@ -143,7 +143,7 @@ subroutine surfstat ! there are two method to calculate surface green's function ! the method in 1985 is better, you can find the ref in the ! subroutine - call surfgreen_1985(w,GLL,GRR,GB,H00,H01,ones, eta_broadening) + call surfgreen_1985(w,GLL,GRR,GB,H00,H01,ones, eta_broadening) ! call surfgreen_1984(w,GLL,GRR,H00,H01,ones, eta_broadening) ! calculate spectral function diff --git a/utility/twisted_graphene_system_tight_binding/Makefile b/utility/twisted_graphene_system_tight_binding/Makefile index c5cab830..e0532e4a 100644 --- a/utility/twisted_graphene_system_tight_binding/Makefile +++ b/utility/twisted_graphene_system_tight_binding/Makefile @@ -6,7 +6,7 @@ obj = module.o readinput.o generate_crystal_structure.o \ #fcheck = -ffree-line-length-0 -g -Wextra -Warray-temporaries -Wconversion -fimplicit-none -fbacktrace -ffree-line-length-0 -fcheck=all -ffpe-trap=zero,overflow,underflow -finit-real=nan f90=ifort -fpp -fcheck = -nogen-interface #-warn all -check all +fcheck = -nogen-interface -warn all -check all flag = -O3 ${fcheck} #-nogen-interface # -warn all # blas and lapack libraries