From 2902d2e75f41653a70e5d1365f63f27b0b5436f7 Mon Sep 17 00:00:00 2001 From: QuanSheng Wu Date: Sat, 9 Dec 2017 22:21:42 +0100 Subject: [PATCH] Don't need this module any more, see example of MoS2 --- 2d/Makefile | 29 - 2d/aux.f90 | 57 -- 2d/dos.f90 | 187 ------- 2d/eigen.f90 | 234 -------- 2d/ek_bulk.f90 | 275 ---------- 2d/ek_ribbon.f90 | 177 ------ 2d/ekb_ribbon.f90 | 132 ----- 2d/fermisurface.f90 | 452 ---------------- 2d/ham_bulk.f90 | 155 ------ 2d/ham_qlayer2qlayer.f90 | 181 ------- 2d/ham_ribbon.f90 | 54 -- 2d/inverse.f90 | 113 ---- 2d/main.f90 | 209 -------- 2d/mat_mul.f90 | 133 ----- 2d/module.f90 | 225 -------- 2d/psi.f90 | 145 ----- 2d/readHmnR.f90 | 99 ---- 2d/readinput.f90 | 544 ------------------- 2d/surfgreen.f90 | 320 ----------- 2d/surfstat.f90 | 228 -------- 2d/surfstat2.f90 | 377 ------------- 2d/symmetry.f90 | 198 ------- 2d/wanniercenter.f90 | 412 -------------- 2d/wanniercenter2.f90 | 1099 -------------------------------------- 24 files changed, 6035 deletions(-) delete mode 100644 2d/Makefile delete mode 100644 2d/aux.f90 delete mode 100644 2d/dos.f90 delete mode 100644 2d/eigen.f90 delete mode 100644 2d/ek_bulk.f90 delete mode 100644 2d/ek_ribbon.f90 delete mode 100644 2d/ekb_ribbon.f90 delete mode 100644 2d/fermisurface.f90 delete mode 100644 2d/ham_bulk.f90 delete mode 100644 2d/ham_qlayer2qlayer.f90 delete mode 100644 2d/ham_ribbon.f90 delete mode 100644 2d/inverse.f90 delete mode 100644 2d/main.f90 delete mode 100644 2d/mat_mul.f90 delete mode 100644 2d/module.f90 delete mode 100644 2d/psi.f90 delete mode 100644 2d/readHmnR.f90 delete mode 100644 2d/readinput.f90 delete mode 100644 2d/surfgreen.f90 delete mode 100644 2d/surfstat.f90 delete mode 100644 2d/surfstat2.f90 delete mode 100644 2d/symmetry.f90 delete mode 100644 2d/wanniercenter.f90 delete mode 100644 2d/wanniercenter2.f90 diff --git a/2d/Makefile b/2d/Makefile deleted file mode 100644 index 5847f6c7..00000000 --- a/2d/Makefile +++ /dev/null @@ -1,29 +0,0 @@ -obj = module.o aux.o symmetry.o readHmnR.o inverse.o\ - eigen.o ham_qlayer2qlayer.o \ - readinput.o surfgreen.o surfstat.o \ - mat_mul.o ham_ribbon.o ek_ribbon.o ham_bulk.o ek_bulk.o \ - wanniercenter.o fermisurface.o dos.o \ - main.o - -# compiler -f90 = mpif90 # -check all -pg -traceback - -#FLAGS = -O3 -nogen-interface -warn all -flag = -O3 -nogen-interface #-warn all - -# blas and lapack libraries -libs = -L/opt/intel/mkl/lib/ \ - -lmkl_intel_lp64 -lmkl_sequential \ - -lmkl_core -liomp5 - -main : $(obj) - $(f90) $(obj) -o wann_tools_2d $(libs) - cp -f wann_tools_2d ../bin - -.SUFFIXES: .o .f90 - -.f90.o : - $(f90) -c $(flag) $(includes) $*.f90 - -clean : - rm -f *.o *.mod *~ wann_tools_2d diff --git a/2d/aux.f90 b/2d/aux.f90 deleted file mode 100644 index 52335bc0..00000000 --- a/2d/aux.f90 +++ /dev/null @@ -1,57 +0,0 @@ -!>> some auxilary subroutines for time or others - -!>> Convert the current wall-time into a real number -!> The unit is second. - - subroutine now(time_now) - - use wmpi - use para, only : dp - implicit none - integer :: time_new(8) - real(dp) :: time_now - call Date_and_time(values=time_new) - time_now= time_new(3)*24*3600+time_new(5)*3600+& - time_new(6)*60+time_new(7)+time_new(8)/1000d0 - - return - end subroutine now - - !> The unit is second. - subroutine print_time_cost(time_start, time_end, subname) - - use wmpi - use para, only : dp, stdout, cpuid - implicit none - real(dp), intent(in) :: time_start - real(dp), intent(in) :: time_end - character(*) :: subname - - if (cpuid==0)write(stdout, 101)'Time cost for ', subname, ' is about ', & - time_end- time_start, ' s' - - 101 format(1x, 3a, f16.3, a) - return - end subroutine print_time_cost - - !> print header - subroutine header - use wmpi - use para, only : stdout, cpuid - implicit none - if (cpuid==0) then - write(stdout, '(a)') " W W W TTTTTTTTTTTTTTTTTTTT " - write(stdout, '(a)') " W W W W TT " - write(stdout, '(a)') " W W W W TT " - write(stdout, '(a)') " W W W W TT " - write(stdout, '(a)') " W W W W TT " - write(stdout, '(a)') " W W W W TT " - write(stdout, '(a)') " WW WW TT " - write(stdout, '(a)') " W W TT " - write(stdout, '(a)') " " - write(stdout, '(a)') " Welcome to Wannier_tools. " - write(stdout, '(a)') " Tools for topological novel materials. " - write(stdout, '(a)') " Enjoy it and good luck. " - endif - end subroutine header - diff --git a/2d/dos.f90 b/2d/dos.f90 deleted file mode 100644 index f6b3f463..00000000 --- a/2d/dos.f90 +++ /dev/null @@ -1,187 +0,0 @@ - -!> calculate density of state and joint density of state for 3D bulk system -!> JDOS(\omega)= \sum_k (f_c(k)-f_v(k) \delta(\omega- Ec(k)+ Ev(k)) -subroutine dos_joint_dos - - use wmpi - use para - implicit none - - !> the integration k space - real(dp) :: kxmin - real(dp) :: kxmax - real(dp) :: kymin - real(dp) :: kymax - real(dp) :: emin - real(dp) :: emax - - integer :: ik - integer :: ik_start - integer :: ik_end - integer :: Nk_local - - integer :: ie - integer :: ib, ib1, ib2 - integer :: ikx - integer :: iky - integer :: NE - integer :: ierr - - !> integration for band - integer :: iband_low - integer :: iband_high - integer :: iband_tot - - real(dp) :: x - real(dp) :: dk2 - - real(dp) :: k(2) - - real(dp), allocatable :: W(:) - real(dp), allocatable :: omega_dos(:) - real(dp), allocatable :: omega_jdos(:) - real(dp), allocatable :: dos(:) - real(dp), allocatable :: dos_mpi(:) - real(dp), allocatable :: jdos(:) - real(dp), allocatable :: jdos_mpi(:) - complex(dp), allocatable :: Hk(:, :) - - !> fermi distribution - real(dp), allocatable :: fermi_dis(:) - - !> delta function - real(dp), external :: delta - - NE= OmegaNum - iband_low= Numoccupied- 40 - iband_high= Numoccupied+ 40 - - if (iband_low <1) iband_low = 1 - if (iband_high >Num_wann) iband_high = Num_wann - - iband_tot= iband_high- iband_low+ 1 - - - allocate(dos(NE)) - allocate(dos_mpi(NE)) - allocate(jdos(NE)) - allocate(jdos_mpi(NE)) - allocate(omega_dos(NE)) - allocate(omega_jdos(NE)) - allocate(W(Num_wann)) - allocate(Hk(Num_wann, Num_wann)) - allocate(fermi_dis(Num_wann)) - W= 0d0 - Hk= 0d0 - fermi_dis= 0d0 - jdos= 0d0 - jdos_mpi= 0d0 - dos= 0d0 - dos_mpi= 0d0 - omega_dos= 0d0 - omega_jdos= 0d0 - - - dk2= 1d0/dble(Nk1*Nk2) - - emin= 0d0 - emax= OmegaMax - eta= (emax- emin)/ dble(NE)*5d0 - - !> energy - do ie=1, NE - omega_jdos(ie)= emin+ (emax-emin)* (ie-1d0)/dble(NE-1) - enddo ! ie - - emin= OmegaMin - emax= OmegaMax - - !> energy - do ie=1, NE - omega_dos(ie)= emin+ (emax-emin)* (ie-1d0)/dble(NE-1) - enddo ! ie - - - !> get eigenvalue - dos_mpi= 0d0 - jdos_mpi= 0d0 - do ik=1+cpuid, Nk1*Nk2, num_cpu - if (cpuid.eq.0) write(stdout, *) 'ik, knv2', ik, Nk1*Nk2 - ikx= (ik- 1)/(Nk2)+ 1 - iky= ik- (ikx-1)*Nk2 - - k= K2D_start+ K2D_vec1*(ikx-1)/dble(nk1-1) & - + K2D_vec2*(iky-1)/dble(nk2-1) - call ham_bulk_old(k, Hk) - W= 0d0 - call eigensystem_c('N', 'U', Num_wann ,Hk, W) - - !> calculate fermi-dirac distribution - do ib=iband_low, iband_high - if (W(ib)<0) then - fermi_dis(ib)= 1d0 - else - fermi_dis(ib)= 0d0 - endif - enddo !ib - - !> get density of state - do ie= 1, NE - do ib1= iband_low, iband_high-1 - do ib2= ib1+1, iband_high - x= omega_jdos(ie)- W(ib2) + W(ib1) - jdos_mpi(ie)= jdos_mpi(ie)+ delta(eta, x)* (fermi_dis(ib1)- fermi_dis(ib2)) - enddo ! ib2 - enddo ! ib1 - enddo ! ie - - !> get density of state - do ie= 1, NE - !> intergrate with k - do ib= iband_low, iband_high-1 - x= omega_dos(ie)- W(ib) - dos_mpi(ie) = dos_mpi(ie)+ delta(eta, x) - enddo ! ib - enddo ! ie - - enddo ! ik - - call mpi_allreduce(dos_mpi,dos,size(dos),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - call mpi_allreduce(jdos_mpi,jdos,size(jdos),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - - if (cpuid.eq.0) then - open(unit=100, file='jdos.dat') - do ie=1, NE - write(100, *)omega_jdos(ie), jdos(ie)*dk2 - enddo ! ie - close(100) - endif - - if (cpuid.eq.0) then - open(unit=100, file='dos.dat') - do ie=1, NE - write(100, *)omega_dos(ie), dos(ie)*dk2 - enddo ! ie - close(100) - endif - - return -end subroutine dos_joint_dos - - -function delta(eta, x) - implicit none - integer, parameter :: dp=kind(1d0) - integer, parameter :: pi= 3.1415926535d0 - real(dp), intent(in) :: eta - real(dp), intent(in) :: x - real(dp) :: delta - - delta= 1d0/pi*eta/(eta*eta+x*x) - !delta= exp(-x*x/eta/eta/2d0)/sqrt(2d0*pi)/eta - -end function - - diff --git a/2d/eigen.f90 b/2d/eigen.f90 deleted file mode 100644 index 690e878e..00000000 --- a/2d/eigen.f90 +++ /dev/null @@ -1,234 +0,0 @@ -! complex version -! a subroutine to calculate eigenvector and eigenvalue - subroutine eigensystem_c(JOBZ,UPLO,N,A,W) - - use para, only : Dp, stdout - implicit none - -! JOBZ (input) CHARACTER*1 -! = 'N': Compute eigenvalues only; -! = 'V': Compute eigenvalues and eigenvectors. -! - character*1, intent(in) :: JOBZ - -! UPLO (input) CHARACTER*1 -! = 'U': Upper triangle of A is stored; -! = 'L': Lower triangle of A is stored. - character*1, intent(in) :: UPLO - -! N (input) INTEGER -! The order of the matrix A. N >= 0. - integer, intent(in) :: N - -! A (input/output) COMPLEX*16 array, dimension (LDA, N) -! On entry, the Hermitian matrix A. If UPLO = 'U', the -! leading N-by-N upper triangular part of A contains the -! upper triangular part of the matrix A. If UPLO = 'L', -! the leading N-by-N lower triangular part of A contains -! the lower triangular part of the matrix A. -! On exit, if JOBZ = 'V', then if INFO = 0, A contains the -! orthonormal eigenvectors of the matrix A. -! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') -! or the upper triangle (if UPLO='U') of A, including the -! diagonal, is destroyed. - - complex(Dp),intent(inout) :: A(N,N) - -! W (output) DOUBLE PRECISION array, dimension (N) -! If INFO = 0, the eigenvalues in ascending order. - - real(Dp), intent(inout) :: W(N) - - integer :: info - - real(Dp),allocatable :: rwork(:) - - complex(Dp),allocatable :: work(:) - - allocate(rwork(16*N)) - allocate( work(16*N)) - rwork= 0d0 - work= 0d0 - - info=0 - W=0.0d0 - - if (N==1) then - W=A(1, 1) - A(1, 1)= 1d0 - return - endif - - call zheev( JOBZ, UPLO, N, A, N, & - W, work, 16*N, rwork, info ) - - if (info.ne.0) then - write(stdout, *) 'ERROR : something wrong with zheev' - stop - endif - - deallocate(rwork, work) - return - end subroutine eigensystem_c - -! real version -! a subroutine to calculate eigenvector and eigenvalue - subroutine eigensystem_r (JOBZ,UPLO,N,A,W) - - use para, only : Dp, stdout - implicit none - -! JOBZ (input) CHARACTER*1 -! = 'N': Compute eigenvalues only; -! = 'V': Compute eigenvalues and eigenvectors. -! - character*1, intent(in) :: JOBZ - -! UPLO (input) CHARACTER*1 -! = 'U': Upper triangle of A is stored; -! = 'L': Lower triangle of A is stored. - character*1, intent(in) :: UPLO - -! N (input) INTEGER -! The order of the matrix A. N >= 0. - integer, intent(in) :: N - -! A (input/output) COMPLEX*16 array, dimension (LDA, N) -! On entry, the Hermitian matrix A. If UPLO = 'U', the -! leading N-by-N upper triangular part of A contains the -! upper triangular part of the matrix A. If UPLO = 'L', -! the leading N-by-N lower triangular part of A contains -! the lower triangular part of the matrix A. -! On exit, if JOBZ = 'V', then if INFO = 0, A contains the -! orthonormal eigenvectors of the matrix A. -! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') -! or the upper triangle (if UPLO='U') of A, including the -! diagonal, is destroyed. - - real(Dp),intent(inout) :: A(N,N) - -! W (output) DOUBLE PRECISION array, dimension (N) -! If INFO = 0, the eigenvalues in ascending order. - - real(Dp), intent(inout) :: W(N) - - integer :: info - - real(Dp),allocatable :: work(:) - - allocate( work(16*N)) - work= 0d0 - - if (N==1) then - W=A(1, 1) - A(1, 1)= 1d0 - return - endif - - info=0 - W=0.0d0 - call dsyev( JOBZ, UPLO, N, A, N, & - W, work, 16*N, info ) - - if (info.ne.0) then - write(stdout, *) 'ERROR : something wrong with dsyev' - stop - endif - - deallocate(work) - return - end subroutine eigensystem_r - -! a subroutine to calculate eigenvector and eigenvalue - subroutine zgeev_pack(N, A, W) - - use para, only : Dp - implicit none - -! JOBVL (input) CHARACTER*1 -! JOBVR (input) CHARACTER*1 -! = 'N': Compute eigenvalues only; -! = 'V': Compute eigenvalues and eigenvectors. -! - character*1 :: JOBVL - character*1 :: JOBVR - -! N (input) INTEGER -! The order of the matrix A. N >= 0. - integer, intent(in) :: N - -! A (input/output) COMPLEX*16 array, dimension (LDA, N) -! On entry, the Hermitian matrix A. If UPLO = 'U', the -! leading N-by-N upper triangular part of A contains the -! upper triangular part of the matrix A. If UPLO = 'L', -! the leading N-by-N lower triangular part of A contains -! the lower triangular part of the matrix A. -! On exit, if JOBZ = 'V', then if INFO = 0, A contains the -! orthonormal eigenvectors of the matrix A. -! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') -! or the upper triangle (if UPLO='U') of A, including the -! diagonal, is destroyed. - - complex(Dp),intent(inout) :: A(N, N) - -! W (output) DOUBLE PRECISION array, dimension (N) -! If INFO = 0, the eigenvalues in ascending order. - -! eigenvalues - complex(Dp), intent(out) :: W(N) - -! left eigenvectors - complex(dp), allocatable :: VL(:, :) - -! right eigenvectors - complex(dp), allocatable :: VR(:, :) - - integer :: info - - integer :: lda - integer :: ldvl - integer :: ldvr - - integer :: lwork - - real(dp),allocatable :: rwork(:) - - complex(dp),allocatable :: work(:) - - !> only calculate eigenvalues - JOBVL= 'N' - JOBVR= 'N' - - lda=N - ldvl=N - ldvr=N - lwork= 16*N - allocate(VL(N, N)) - allocate(VR(N, N)) - - allocate(rwork(16*N)) - allocate( work(lwork)) - VL= 0d0 - VR= 0d0 - rwork= 0d0 - work= 0d0 - - info=0 - W=0.0d0 - - if (N==1) then - W=A(1, 1) - VL(1, 1)= 1d0 - VR(1, 1)= 1d0 - return - endif - - call ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, & - WORK, LWORK, RWORK, INFO ) - if (info /= 0) then - stop ">>> Error : something wrong happens in zgeev_pack" - endif - - return - end subroutine zgeev_pack - diff --git a/2d/ek_bulk.f90 b/2d/ek_bulk.f90 deleted file mode 100644 index 602eaa62..00000000 --- a/2d/ek_bulk.f90 +++ /dev/null @@ -1,275 +0,0 @@ -! calculate bulk's energy band using wannier TB method -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - - subroutine ek_bulk - - use wmpi - use para - - implicit none - - integer :: ik, i, j - integer :: Nwann - integer :: ierr - real(dp) :: emin - real(dp) :: emax - real(Dp) :: k(2) - real(Dp) :: W(Num_wann) - - ! Hamiltonian of bulk system - complex(Dp) :: Hamk_bulk(Num_wann,Num_wann) - - ! eigen value of H - real(dp), allocatable :: eigv(:,:) - real(dp), allocatable :: eigv_mpi(:,:) - real(dp), allocatable :: weight(:,:,:) - real(dp), allocatable :: weight_mpi(:,:,:) - - allocate( eigv (Num_wann, knv2)) - allocate( eigv_mpi(Num_wann, knv2)) - allocate( weight (Num_wann,Num_wann, knv2)) - allocate( weight_mpi(Num_wann,Num_wann, knv2)) - eigv = 0d0 - eigv_mpi= 0d0 - weight = 0d0 - weight_mpi = 0d0 - - do ik= 1+cpuid, knv2, num_cpu - if (cpuid==0) write(stdout, *)'BulkBand, ik, knv2 ', ik, knv2 - - k = k2_path(ik, :) - - ! calculation bulk hamiltonian - Hamk_bulk= 0d0 - call ham_bulk_old(k, Hamk_bulk) - !call ham_bulk (k, Hamk_bulk) - - !> diagonalization by call zheev in lapack - W= 0d0 - call eigensystem_c( 'V', 'U', Num_wann ,Hamk_bulk, W) - - eigv(:, ik)= W - do i=1, Num_wann !> band - do j=1, Num_wann/2 !> projector - weight(j, i, ik)= sqrt(abs(Hamk_bulk(j, i))**2+ & - abs(Hamk_bulk(j+ Num_wann/2, i))**2) - enddo ! j - enddo ! i - enddo ! ik - - 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) - - weight= weight_mpi/maxval(weight_mpi)*255d0 - - if (cpuid==0)then - open(unit=14, file='bulkek.dat') - - do i=1, Num_wann - do ik=1, knv2 - write(14, '(2f19.9, 1000i5)')k2len(ik),eigv_mpi(i, ik), & - int(weight(:, i, ik)) - enddo - write(14, *)' ' - enddo - close(14) - endif - - !> minimum and maximum value of energy bands - !emin= minval(eigv_mpi)-0.5d0 - !emax= maxval(eigv_mpi)+0.5d0 - emin= -5d0 - emax= 5d0 - - !> write script for gnuplot - if (cpuid==0) then - open(unit=101, file='bulkek.gnu') - write(101, '(a)') 'set terminal postscript enhanced color font ",30"' - write(101,'(2a)') '#set palette defined ( 0 "green", ', & - '5 "yellow", 10 "red" )' - write(101, '(a)')"set output 'bulkek.eps'" - write(101, '(a)')'set style data linespoints' - write(101, '(a)')'unset ztics' - write(101, '(a)')'unset key' - write(101, '(a)')'set pointsize 0.8' - write(101, '(a)')'set view 0,0' - write(101, '(a)')'set xtics font ",24"' - write(101, '(a)')'set ytics font ",24"' - write(101, '(a)')'set ylabel font ",24"' - write(101, '(a)')'set ylabel "Energy (eV)"' - write(101, '(a)')'set ylabel offset 1.5,0' - write(101, '(a, f10.5, a)')'set xrange [0: ', maxval(k2len), ']' - write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', emin, ':', emax, ']' - write(101, 202, advance="no") (k2line_name(i), k2line_stop(i), i=1, nk2lines) - write(101, 203)k2line_name(nk2lines+1), k2line_stop(nk2lines+1) - - do i=1, nk2lines-1 - write(101, 204)k2line_stop(i+1), emin, k2line_stop(i+1), emax - enddo - write(101, '(2a)')"plot 'bulkek.dat' u 1:2 ", & - "w lp lw 2 pt 7 ps 0.2 lc rgb 'black', 0 w l lw 2" - close(101) - endif - - 202 format('set xtics (',:20('"',A3,'" ',F10.5,',')) - 203 format(A3,'" ',F10.5,')') - 204 format('set arrow from ',F10.5,',',F10.5,' to ',F10.5,',',F10.5, ' nohead') - - return - end subroutine ek_bulk - - - ! calculate bulk's energy band using wannier TB method - !> calculate spin direction for each band and each kpoint - subroutine ek_bulk_spin - - use wmpi - use para - - implicit none - - integer :: ik, i, j, ib - integer :: ierr - integer :: nwann - real(dp) :: sx, sy, sz - real(dp) :: emin - real(dp) :: emax - real(dp) :: k(2) - real(dp) :: W(Num_wann) - - ! Hamiltonian of bulk system - complex(dp), allocatable :: Hamk_bulk(:, :) - complex(dp), allocatable :: sigmax(:, :) - complex(dp), allocatable :: sigmay(:, :) - complex(dp), allocatable :: sigmaz(:, :) - complex(dp), allocatable :: psi(:) - - ! eigen value of H - real(dp), allocatable :: eigv(:,:) - real(dp), allocatable :: eigv_mpi(:,:) - real(dp), allocatable :: spin(:, :, :) - real(dp), allocatable :: spin_mpi(:, :, :) - - allocate( psi( Num_wann)) - allocate( Hamk_bulk( Num_wann, Num_wann)) - allocate( sigmax( Num_wann, Num_wann)) - allocate( sigmay( Num_wann, Num_wann)) - allocate( sigmaz( Num_wann, Num_wann)) - allocate( eigv (Num_wann, knv2)) - allocate( eigv_mpi(Num_wann, knv2)) - allocate( spin (2, Num_wann, knv2)) - allocate( spin_mpi(2, Num_wann, knv2)) - psi = 0d0 - eigv = 0d0 - eigv_mpi= 0d0 - spin= 0d0 - spin_mpi= 0d0 - sigmax= 0d0 - sigmay= 0d0 - sigmaz= 0d0 - - nwann= Num_wann/2 - - if (SOC==0) stop 'you should set soc=0 in the input file' - !> construct spin matrix - do i= 1, nwann - sigmax(i, i+ nwann)= 1d0 - sigmax(i+ nwann, i)= 1d0 - sigmay(i, i+ nwann)= -zi - sigmay(i+ nwann, i)= zi - sigmaz(i, i)= 1d0 - sigmaz(i+ nwann, i+ nwann)= -1d0 - enddo - - do ik= 1+cpuid, knv2, num_cpu - if (cpuid==0) write(stdout, *)'BulkBandSpin, ik, knv2 ', ik, knv2 - - k = k2_path(ik, :) - - ! calculation bulk hamiltonian - Hamk_bulk= 0d0 - call ham_bulk(k, Hamk_bulk) - - !> diagonalization by call zheev in lapack - W= 0d0 - call eigensystem_c( 'V', 'U', Num_wann ,Hamk_bulk, W) - - eigv(:, ik)= W - - do ib= 1, Num_wann - psi(:)= Hamk_bulk(:, ib) - sx= 0d0 - sy= 0d0 - sz= 0d0 - do i= 1, Num_wann - do j= 1, Num_wann - sx= sx+ conjg(psi(i))* sigmax(i, j)* psi(j) - sy= sy+ conjg(psi(i))* sigmay(i, j)* psi(j) - enddo ! j - enddo ! i - spin(1, ib, ik)= sx - spin(2, ib, ik)= sy - enddo ! ib - - enddo ! ik - - call mpi_allreduce(eigv,eigv_mpi,size(eigv),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - call mpi_allreduce(spin,spin_mpi,size(spin),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - - if (cpuid==0)then - open(unit=14, file='bulkek.dat') - - do i=1, Num_wann - do ik=1, knv2 - write(14, '(1000f19.9)')k2len(ik),eigv_mpi(i, ik), spin(:, i, ik) - enddo - write(14, *)' ' - enddo - close(14) - endif - - !> minimum and maximum value of energy bands - !emin= minval(eigv_mpi)-0.5d0 - !emax= maxval(eigv_mpi)+0.5d0 - emin= -5d0 - emax= 5d0 - - !> write script for gnuplot - if (cpuid==0) then - open(unit=101, file='bulkek.gnu') - write(101, '(a)') 'set terminal postscript enhanced color font 24' - write(101,'(2a)') '#set palette defined ( 0 "green", ', & - '5 "yellow", 10 "red" )' - write(101, '(a)')"set output 'bulkek.eps'" - write(101, '(a)')'set style data linespoints' - write(101, '(a)')'unset ztics' - write(101, '(a)')'unset key' - write(101, '(a)')'set pointsize 0.8' - write(101, '(a)')'set view 0,0' - write(101, '(a)')'set ylabel "Energy (eV)"' - write(101, '(a, f10.5, a)')'set xrange [0: ', maxval(k2len), ']' - write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', emin, ':', emax, ']' - write(101, 202, advance="no") (k2line_name(i), k2line_stop(i), i=1, nk2lines) - write(101, 203)k2line_name(nk2lines+1), k2line_stop(nk2lines+1) - - do i=1, nk2lines-1 - write(101, 204)k2line_stop(i+1), emin, k2line_stop(i+1), emax - enddo - write(101, '(2a)')"plot 'bulkek.dat' u 1:2 ", & - "w lp lw 2 pt 7 ps 0.2, \" - write(101, '(2a)')" 'bulkek.dat' u 1:2:($3/6):($4/6) ", & - "w vec" - close(101) - endif - - 202 format('set xtics (',:20('"',A3,'" ',F10.5,',')) - 203 format(A3,'" ',F10.5,')') - 204 format('set arrow from ',F10.5,',',F10.5,' to ',F10.5,',',F10.5, ' nohead') - - return - end subroutine ek_bulk_spin - diff --git a/2d/ek_ribbon.f90 b/2d/ek_ribbon.f90 deleted file mode 100644 index d0d4bde5..00000000 --- a/2d/ek_ribbon.f90 +++ /dev/null @@ -1,177 +0,0 @@ -!> This subroutine is used for calculating energy -!> dispersion with wannier functions -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - subroutine ek_ribbon - - use wmpi - use para - implicit none - -! loop index - integer :: i - integer :: j - integer :: l - - integer :: lwork - -! wave vector - real(Dp) :: k - - real(dp) :: emin, emax - -! parameters for zheev - integer :: ierr - real(Dp), allocatable :: rwork(:) - complex(Dp), allocatable :: work(:) - -! eigenvalue - real(Dp), allocatable :: eigenvalue(:) - -! energy dispersion - real(Dp),allocatable :: ekribbon(:,:) - real(Dp),allocatable :: ekribbon_mpi(:,:) - - - !> color for plot, surface state weight - real(dp), allocatable :: surf_weight(:, :) - real(dp), allocatable :: surf_weight_mpi(:, :) - -! hamiltonian ribbon - complex(Dp),allocatable ::CHamk(:,:) - - lwork= 16*Nslab*Num_wann - ierr = 0 - - - allocate(eigenvalue(nslab*Num_wann)) - allocate( surf_weight (Nslab* Num_wann, knv1)) - allocate( surf_weight_mpi (Nslab* Num_wann, knv1)) - allocate(ekribbon(Nslab*Num_wann,knv1)) - allocate(ekribbon_mpi(Nslab*Num_wann,knv1)) - allocate(CHamk(nslab*Num_wann,nslab*Num_wann)) - allocate(work(lwork)) - allocate(rwork(lwork)) - - surf_weight= 0d0 - surf_weight_mpi= 0d0 - - ! sweep k - ekribbon=0.0d0 - ekribbon_mpi=0.0d0 - do i=1+cpuid,knv1,num_cpu - if (cpuid==0) write(stdout, *), 'ribbonEk, ik ', i, knv1 - k= k1_path(i) - chamk=0.0d0 - - !> no magnetgic field - if (abs(Bx) in-plane magnetic field - elseif (abs(Bx)>eps9 .or. abs(By)>eps9)then - write(stdout, *)'>> magnetic field is larger than zero' - call ham_ribbon(k,Chamk) - !call ham_slab_parallel_B(k,Chamk) - !> vertical magnetic field - else - print *, 'Error: we only support in-plane magnetic field at present' - stop 'please set Bz= 0' - endif - - - eigenvalue=0.0d0 - - ! diagonal Chamk - call eigensystem_c('V', 'U', Num_wann*Nslab, CHamk, eigenvalue) - - ekribbon(:,i)=eigenvalue - - do j=1, Nslab*Num_wann - do l=1, Num_wann - surf_weight(j, i)= surf_weight(j, i) & - + abs(CHamk(l, j))**2 & ! first slab - + abs(CHamk(Num_wann*Nslab- l+ 1, j))**2 !& ! last slab - !+ abs(CHamk(Num_wann+ l, j))**2 & ! the second slab - !+ abs(CHamk(Num_wann*(Nslab-1)- l, j))**2 ! last second slab - enddo ! l - surf_weight(j, i)= sqrt(surf_weight(j, i)) - enddo ! j - enddo ! i - - call mpi_allreduce(ekribbon,ekribbon_mpi,size(ekribbon),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - call mpi_allreduce(surf_weight, surf_weight_mpi,size(surf_weight),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - - - ekribbon=ekribbon_mpi - if (maxval(surf_weight_mpi)<0.00001d0)surf_weight_mpi=1d0 - surf_weight= surf_weight_mpi/ maxval(surf_weight_mpi) - - - if(cpuid==0)then - open(unit=100, file='ribbonek.dat') - do j=1, Num_wann*Nslab - do i=1, knv2 - !write(100,'(3f15.7, i8)')k1len(i), ekslab(j,i), & - ! (surf_weight(j, i)) - write(100,'(2f15.7, i8)')k1len(i), ekribbon(j,i), & - int(255-surf_weight(j, i)*255d0) - enddo - write(100 , *)'' - enddo - close(100) - write(stdout,*) 'calculate energy band done' - endif - - emin= minval(ekribbon)-0.5d0 - emax= maxval(ekribbon)+0.5d0 - !> write script for gnuplot - if (cpuid==0) then - open(unit=113, file='ribbonek.gnu') - write(113, '(a)')"set encoding iso_8859_1" - write(113, '(a)')'#set terminal postscript enhanced color' - write(113, '(a)')"#set output 'ribbonek.eps'" - write(113, '(3a)')'#set terminal pngcairo truecolor enhanced', & - ' font ",60" size 1920, 1680' - write(113, '(3a)')'set terminal png truecolor enhanced', & - ' font ",60" size 1920, 1680' - write(113, '(a)')"set output 'ribbonek.png'" - write(113,'(2a)') 'set palette defined ( 0 "green", ', & - '5 "yellow", 10 "red" )' - write(113, '(a)')'set style data linespoints' - write(113, '(a)')'unset ztics' - write(113, '(a)')'unset key' - write(113, '(a)')'set pointsize 0.8' - write(113, '(a)')'set border lw 3 ' - write(113, '(a)')'set view 0,0' - write(113, '(a)')'#set xtics font ",36"' - write(113, '(a)')'#set ytics font ",36"' - write(113, '(a)')'#set ylabel font ",36"' - write(113, '(a)')'#set xtics offset 0, -1' - write(113, '(a)')'set ylabel offset -1, 0 ' - write(113, '(a)')'set ylabel "Energy (eV)"' - write(113, '(a, f10.5, a)')'set xrange [0: ', maxval(k1len), ']' - write(113, '(a, f10.5, a, f10.5, a)')'set yrange [', emin, ':', emax, ']' - write(113, 202, advance="no") (trim(k1line_name(i)), k1line_stop(i), i=1, nk1lines) - write(113, 203)trim(k1line_name(nk1lines+1)), k1line_stop(nk1lines+1) - - do i=1, nk1lines-1 - write(113, 204)k1line_stop(i+1), emin, k1line_stop(i+1), emax - enddo - write(113, '(a)')'rgb(r,g,b) = int(r)*65536 + int(g)*256 + int(b)' - write(113, '(2a)')"plot 'ribbonek.dat' u 1:2:(rgb(255,$3, 3)) ", & - "w lp lw 2 pt 7 ps 1 lc rgb variable" - - !write(113, '(2a)')"splot 'ribbonek.dat' u 1:2:3 ", & - ! "w lp lw 2 pt 13 palette" - close(113) - endif - - 202 format('set xtics (',:20('"',A3,'" ',F10.5,',')) - 203 format(A3,'" ',F10.5,')') - 204 format('set arrow from ',F10.5,',',F10.5, & - ' to ',F10.5,',',F10.5, ' nohead') - - return - end subroutine ek_ribbon - diff --git a/2d/ekb_ribbon.f90 b/2d/ekb_ribbon.f90 deleted file mode 100644 index 71cf62d6..00000000 --- a/2d/ekb_ribbon.f90 +++ /dev/null @@ -1,132 +0,0 @@ -! This subroutine is used to caculate energy dispersion for -! slab Bi2Se3 -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - - subroutine ekb_ribbon - - use para - implicit none - - integer :: mdim - integer :: ndim1 - - ! loop index - integer :: i - - integer :: lwork - - ! the indices of the smallest and largest - ! eigenvalues to be returned - integer :: il,iu - - - ! wave vector - real(Dp) :: k - - real(Dp) :: kmax=0.5d0 - - !the lower and upper bounds of the interval - !to be searched for eigenvalues - real(Dp) :: vl,vu - - !The absolute error tolerance for the eigenvalues - real(Dp) :: abstol - - ! parameters for zheev - integer :: info,ierr,ifail - - ! number of magnetic fields - integer :: Nb - - ! magnetic field - real(dp) :: Bmag - real(dp), allocatable :: magnetic(:) - - ! energy dispersion - real(Dp),allocatable :: ekbribbon(:,:) - real(Dp),allocatable :: rwork(:) - integer,allocatable :: iwork(:) - complex(Dp),allocatable :: work(:) - - ! eigenvalue - real(Dp),allocatable :: eigenvalue(:) - - ! hamiltonian slab - complex(Dp),allocatable ::z(:,:) - complex(Dp),allocatable ::CHamk(:,:) - - - Nb= 400 - Ndim1=Num_wann*nslab1*nslab2 - lwork=64*Ndim1 - - allocate(ekbribbon(Ndim1,Nb )) - allocate(z(Ndim1,Ndim1)) - allocate(CHamk(Ndim1,Ndim1)) - allocate(rwork(7*Ndim1)) - allocate(work(lwork)) - allocate(iwork(5*Ndim1)) - allocate(eigenvalue(Ndim1)) - allocate(magnetic(Nb)) - - - ierr = 0 - - ! sweep k - ekbribbon=0.0d0 - - kmax=0.5d0 - - abstol=1e-10 - vl=-0.6/27.2114d0 - vu=0.5/27.2114d0 - il=17*Nslab1*Nslab2 - iu=20*Nslab1*Nslab2 - mdim=iu-il+1 - - do i=1, Nb - magnetic(i)= (i-1)*1.0d0/real(Nb) - enddo - - do i=1,Nb - Bmag= magnetic(i) - k=0d0 - chamk=0.0d0 - call ham_ribbon_b(bmag, k,Chamk) - eigenvalue=0.0d0 - - ! diagonal Chamk - call eigensystem_c('N', 'U', Ndim1, CHamk, eigenvalue) - - ! only eigenvalues are computed - ! the eigenvalues with indices il through iu will be found - - !call zheevx('N','I','U',Ndim1,Chamk,Ndim1,vl,vu,il,iu,abstol,& - !mdim,eigenvalue,z,Ndim1,work,lwork,rwork,iwork,ifail,info) - - - ekbribbon(:,i)=eigenvalue - - write(stdout,'(a2,i4,f12.5,f10.2,a2)')'B, Nb', i, Nb, ekbribbon(1,i) - enddo - - !open(unit=100, file='ribbonek.dat',status='unknown') - !do i=1,Nk1 - ! k=-kmax*real(Nk1-i)/(Nk1-1) - ! write(100,'(60000f15.7)')k,ekbribbon(:,Nk1-i+1) - !enddo - - !do i=1,Nk1 - ! k=kmax*real(i-1)/(Nk1-1) - ! write(100,'(60000f15.7)')k,ekbribbon(:,i) - !enddo - !close(100) - - open(unit=100, file='ribbonekb.dat',status='unknown') - do i=1, Nb - write(100, '(50000f15.7)')magnetic(i), ekbribbon(:, i) - enddo - !write(stdout,*) 'calculate energy band done' - - return - end diff --git a/2d/fermisurface.f90 b/2d/fermisurface.f90 deleted file mode 100644 index 80837c2e..00000000 --- a/2d/fermisurface.f90 +++ /dev/null @@ -1,452 +0,0 @@ -! calculate bulk's energy band using wannier TB method - subroutine fermisurface - - use wmpi - use para - - implicit none - - integer :: ik, i, j, ikx, iky - integer :: nkx - integer :: nky - - integer :: ierr - real(dp) :: kz - real(Dp) :: k(2) - - ! Hamiltonian of bulk system - complex(Dp) :: Hamk_bulk(Num_wann,Num_wann) - - real(dp) :: zmin, zmax - real(dp) :: kxmin, kxmax, kymin, kymax - real(dp) :: kxmin_shape, kxmax_shape, kymin_shape, kymax_shape - - real(dp), allocatable :: kxy(:,:) - real(dp), allocatable :: kxy_shape(:,:) - - real(dp), allocatable :: dos(:) - real(dp), allocatable :: dos_mpi(:) - - complex(dp), allocatable :: ones(:,:) - - nkx= Nk1 - nky= Nk2 - allocate( kxy(2, nkx*nky)) - allocate( kxy_shape(2, nkx*nky)) - kxy=0d0 - kxy_shape=0d0 - - ik=0 - do ikx=1, Nk1 - do iky=1, Nk2 - ik= ik+ 1 - kxy(:, ik)= K2D_start+ K2D_vec1*(ikx-1)/dble(nk1-1) & - + K2D_vec2*(iky-1)/dble(nk2-1) - kxy_shape(:, ik)= kxy(1,ik)*kua+ kxy(2, ik)*kub - enddo - enddo - - kxmin_shape= minval(kxy_shape(1, :)) - kxmax_shape= maxval(kxy_shape(1, :)) - kymin_shape= minval(kxy_shape(2, :)) - kymax_shape= maxval(kxy_shape(2, :)) - - allocate( dos (Nk1*Nk2)) - allocate( dos_mpi(Nk1*Nk2)) - dos = 0d0 - dos_mpi= 0d0 - - allocate(ones(Num_wann, Num_wann)) - ones= 0d0 - do i=1, Num_wann - ones(i, i)= 1d0 - enddo - do ik= 1+cpuid, Nk1*Nk2, num_cpu - if (cpuid==0) write(stdout, *),'FS, ik, knv3' , ik, Nk1*Nk2 - - k(1) = kxy(1, ik) - k(2) = kxy(2, ik) - - ! calculation bulk hamiltonian - Hamk_bulk= 0d0 - call ham_bulk(k, Hamk_bulk) - - Hamk_bulk= (E_arc -zi* eta_arc)* ones - Hamk_bulk - call inv(Num_wann, Hamk_bulk) - do i=1, Num_wann - dos(ik)= dos(ik)+ aimag(Hamk_bulk(i, i))/pi - enddo - - enddo - - call mpi_allreduce(dos,dos_mpi,size(dos),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - - if (cpuid==0)then - open(unit=14, file='fs.dat') - - do ik=1, Nk1*Nk2 - write(14, '(3f16.8)')kxy_shape(:, ik), log(dos_mpi(ik)) - if (mod(ik, nky)==0) write(14, *)' ' - enddo - close(14) - endif - zmax= maxval(log(dos_mpi)) - zmin= minval(log(dos_mpi)) - - !> minimum and maximum value of energy bands - - !> write script for gnuplot - if (cpuid==0) then - open(unit=101, file='fs.gnu') - write(101, '(a)')"set encoding iso_8859_1" - write(101, '(a)')'#set terminal postscript enhanced color' - write(101, '(a)')"#set output 'fs.eps'" - write(101, '(3a)')'set terminal pngcairo truecolor enhanced', & - ' size 1920, 1680 font ",36"' - write(101, '(a)')"set output 'fs.png'" - write(101,'(a, f10.4, 2a, f10.4, a)') & - 'set palette defined ( ', zmin, ' "white", ', & - '0 "black", ', zmax,' "red" )' - write(101, '(a)')'#set palette rgbformulae 33,13,10' - write(101, '(a)')'unset ztics' - write(101, '(a)')'unset key' - write(101, '(a)')'set pm3d' - write(101, '(a)')'#set view equal xyz' - write(101, '(a)')'set view map' - write(101, '(a)')'set border lw 3' - write(101, '(a)')'#set xtics font ",24"' - write(101, '(a)')'#set ytics font ",24"' - write(101, '(a)')'set size ratio -1' - write(101, '(a)')'unset xtics' - write(101, '(a)')'unset ytics' - write(101, '(a)')'set colorbox' - !write(101, '(a, f10.5, a, f10.5, a)')'set xrange [', kxmin, ':', kxmax, ']' - !write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', kymin, ':', kymax, ']' - write(101, '(a, f10.5, a, f10.5, a)')'set xrange [', kxmin_shape, ':', kxmax_shape, ']' - write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', kymin_shape, ':', kymax_shape, ']' - write(101, '(a)')'set pm3d interpolate 2,2' - write(101, '(2a)')"splot 'fs.dat' u 1:2:3 w pm3d" - - close(101) - endif - - - return - end subroutine fermisurface - -! calculate bulk's energy band using wannier TB method - subroutine gapshape - - use wmpi - use para - - implicit none - - integer :: ik, i, j, ikx, iky - integer :: nkx - integer :: nky - - integer :: ierr, i1, i2 - real(Dp) :: k(2) - - ! Hamiltonian of bulk system - complex(Dp) :: Hamk_bulk(Num_wann,Num_wann) - - real(dp) :: zmin, zmax - real(dp) :: kxmin_shape, kxmax_shape, kymin_shape, kymax_shape - - real(dp), allocatable :: kxy(:,:) - real(dp), allocatable :: kxy_shape(:,:) - - real(dp), allocatable :: gap(:, :) - real(dp), allocatable :: gap_mpi(:, :) - real(dp), allocatable :: W(:) - - complex(dp), allocatable :: ones(:,:) - - nkx= Nk - nky= Nk - allocate( kxy(2, nkx*nky)) - allocate( kxy_shape(2, nkx*nky)) - kxy=0d0 - kxy_shape=0d0 - - - ik=0 - do ikx=1, Nk1 - do iky=1, Nk2 - ik= ik+ 1 - kxy(:, ik)= K2D_start+ K2D_vec1*(ikx-1)/dble(nk1-1) & - + K2D_vec2*(iky-1)/dble(nk2-1) - kxy_shape(:, ik)= kxy(1,ik)*kua+ kxy(2, ik)*kub - enddo - enddo - - kxmin_shape=minval(kxy_shape(1,:)) - kxmax_shape=maxval(kxy_shape(1,:)) - kymin_shape=minval(kxy_shape(2,:)) - kymax_shape=maxval(kxy_shape(2,:)) - - - allocate( gap (3, Nk1*Nk2)) - allocate( gap_mpi(3, Nk1*Nk2)) - gap = 0d0 - gap_mpi= 0d0 - - allocate(W(Num_wann)) - allocate(ones(Num_wann, Num_wann)) - W= 0d0 - ones= 0d0 - do i=1, Num_wann - ones(i, i)= 1d0 - enddo - - if (Numoccupied> Num_wann) then - stop 'Numoccupied should less than Num_wann' - endif - - do ik= 1+cpuid, Nk1*Nk2, num_cpu - if (cpuid==0)write(stdout, *)'Gap plane', ik, Nk1*Nk2 - - k(1) = kxy(1, ik) - k(2) = kxy(2, ik) - - !> calculation bulk hamiltonian - Hamk_bulk= 0d0 - call ham_bulk(k, Hamk_bulk) - - !> diagonalization - call eigensystem_c( 'N', 'U', Num_wann ,Hamk_bulk, W) - gap(1, ik)= W(Numoccupied+1)- W(Numoccupied) - gap(2, ik)= W(Numoccupied) - gap(3, ik)= W(Numoccupied+1) - - enddo - - call mpi_allreduce(gap,gap_mpi,size(gap),& - mpi_dp,mpi_sum,mpi_cmw,ierr) - - if (cpuid==0)then - open(unit=14, file='GapPlane.dat') - - write(14, '(100a16)')'# kx', 'ky', 'gap', 'Ev', 'Ec', 'k1', 'k2' - do ik=1, Nk1*Nk2 - write(14, '(30f16.8)')kxy_shape(:, ik), (gap_mpi(:, ik)), kxy(:, ik) - if (mod(ik, nky)==0) write(14, *)' ' - enddo - close(14) - - open(unit=15, file='gap2d.dat') - write(15, '(100a16)')'kx', 'ky', 'gap', 'Ev', 'Ec', 'k1', 'k2' - do ik=1, Nk1*Nk2 - if (abs(gap_mpi(1, ik))< 0.10d0) then - write(15, '(8f16.8)')kxy_shape(:, ik), (gap_mpi(:, ik)), kxy(:, ik) - endif - enddo - close(15) - - open(unit=1116, file='GapPlane_matlab.dat') - - write(1116, '(100a16)')'% kx', 'ky', 'gap', 'Ev', 'Ec', 'k1', 'k2' - do ik=1, Nk1*Nk2 - write(1116, '(30f16.8)')kxy_shape(:, ik), (gap_mpi(:, ik)), kxy(:, ik) - enddo - close(1116) - - endif - - !> minimum and maximum value of energy bands - - zmax= maxval(gap_mpi(1, :)) - zmin= minval(gap_mpi(1, :)) - - !> write script for gnuplot - if (cpuid==0) then - open(unit=103, file='GapPlane.gnu') - write(103, '(a)')"set encoding iso_8859_1" - write(103, '(a)')'#set terminal postscript enhanced color' - write(103, '(a)')"#set output 'GapPlane.eps'" - write(103, '(3a)')'#set terminal pngcairo truecolor enhanced', & - ' size 1920, 1680 font ",60"' - write(103, '(3a)')'set terminal png truecolor enhanced', & - ' size 1920, 1680 font ",60"' - write(103, '(a)')"set output 'GapPlane.png'" - write(103,'(a, f10.4, a, f10.4, a, f10.4, a)') & - 'set palette defined ( ', zmin, ' "black", ', & - (zmin+zmax)/20d0,' "orange", ',zmax,' "white" )' - write(103, '(a)')"set origin 0.10, 0.0" - write(103, '(a)')"set size 0.85, 1.0" - write(103, '(a)')'unset ztics' - write(103, '(a)')'unset key' - write(103, '(a)')'set pm3d' - write(103, '(a)')'set view map' - write(103, '(a)')'set border lw 3' - write(103, '(a)')'#set size ratio -1' - write(103, '(a)')'set title "Gap in k plane"' - write(103, '(a)')'set xtics nomirror scale 0.5' - write(103, '(a)')'set ytics nomirror scale 0.5' - write(103, '(a)')"set xlabel 'k (1/{\305})'" - write(103, '(a)')"set ylabel 'k (1/{\305})'" - write(103, '(a)')'set colorbox' - write(103, '(a)')'set xrange [ ] noextend' - write(103, '(a)')'set yrange [ ] noextend' - write(103, '(a)')'set pm3d interpolate 2,2' - write(103, '(2a)')"splot 'GapPlane.dat' u 1:2:3 w pm3d" - - close(103) - endif - - - return - end subroutine gapshape - - - !> get fermilevel for the given hamiltonian - subroutine get_fermilevel - use wmpi - use para - implicit none - - integer :: i1 - integer :: i2 - integer :: i3 - integer :: io - integer :: ik - - !> number of k points - - integer :: ierr - integer :: iter - integer :: itermax - - !> fermi level - real(dp) :: EF - - real(dp) :: k(3) - - real(dp) :: Beta - - real(dp) :: lmin - real(dp) :: lmax - real(dp) :: tot - - - !> fermi-dirac distribution function - real(dp), external :: fermi - - !> kpoint coordinates - real(dp), allocatable :: kpoints(:, :) - - !> eigen value for each kpoint - real(dp), allocatable :: W(:) - real(dp), allocatable :: eigvals(:, :) - real(dp), allocatable :: eigvals_mpi(:, :) - - complex(dp), allocatable :: ham(:, :) - - - allocate(W(Num_wann)) - allocate(eigvals(Num_wann, Nk1*Nk2)) - allocate(eigvals_mpi(Num_wann, Nk1*Nk2)) - allocate(ham(Num_wann, Num_wann)) - allocate(kpoints(2, Nk1*Nk2)) - eigvals= 0d0 - eigvals_mpi= 0d0 - ham= 0d0 - kpoints= 0d0 - Beta= 200d0 - - ik= 0 - do i1=1, Nk - do i2=1, Nk - do i3=1, Nk - ik= ik+ 1 - kpoints(1, ik)= (i1-1d0)/dble(Nk) - kpoints(2, ik)= (i2-1d0)/dble(Nk) - kpoints(3, ik)= (i3-1d0)/dble(Nk) - enddo - enddo - enddo - - do ik=1+ cpuid, Nk1*Nk2, num_cpu - - ham= 0d0 - k= kpoints(:, ik) - call ham_bulk(k, ham) - call eigensystem_c( 'N', 'U', num_wann, ham, W) - eigvals_mpi(:, ik)= W - enddo ! ik - - call mpi_allreduce(eigvals_mpi, eigvals, size(eigvals), & - mpi_dp, mpi_sum, mpi_cmw, ierr) - - ! using bisection algorithm to search the fermi level - iter= 0 - itermax= 100 - tot= 9999d0 - lmin= minval(eigvals) - lmax= maxval(eigvals) - if (cpuid==0) print *, 'Emin= ', lmin - if (cpuid==0) print *, 'Emax= ', lmax - do while( abs(tot- Ntotch).gt. eps6 .and. iter.lt.itermax) - - iter= iter+ 1 - - EF= (lmin+ lmax)* half - - tot= 0d0 - do ik=1, Nk1*Nk2 - do io=1, Num_wann - tot= tot+ fermi(eigvals(io, ik)- EF, Beta) - enddo ! io - enddo ! ik - tot= tot/dble(Nk1*Nk2) - - if (tot > Ntotch)then - lmax= EF - else - lmin= EF - endif - - if (cpuid==0) then - write(stdout, 100)iter, tot-Ntotch, EF, ' Charge: ', tot - endif - 100 format(2x,">iter",i4,2x,"diff:",f12.6,2x,"EF: ",f12.6,a,f12.6) - - enddo ! bisection - - E_fermi= EF - - return - end subroutine get_fermilevel - - !------------+------------+------------+------------+------------+--------+! - ! calculate the Fermi-Dirac distribution - !------------+------------+------------+------------+------------+--------+! - function fermi(omega, Beta) result(value) - - use para - implicit none - - ! >> inout variables - real(dp), intent(in) :: omega - real(dp), intent(in) :: Beta - - ! return value - real(dp) :: value - - ! avoid numerical instability - if (beta*omega .ge. 20d0) then - value = zero - elseif (beta*omega.le. -20d0)then - value = one - else - value= one/(one+exp(beta*omega)) - endif - - return - end function fermi - - - diff --git a/2d/ham_bulk.f90 b/2d/ham_bulk.f90 deleted file mode 100644 index 6823c149..00000000 --- a/2d/ham_bulk.f90 +++ /dev/null @@ -1,155 +0,0 @@ -! This subroutine is used to caculate Hamiltonian for -! bulk system . - -! History -! May/29/2011 by Quansheng Wu -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - subroutine ham_bulk(k,Hamk_bulk) - - use para - implicit none - -! loop index - integer :: i1,i2,ia,ib,iR, ia1, ia2 - integer :: istart1, istart2, iend1, iend2 - integer :: nwann - - real(dp) :: kdotr - -! wave vector in 2d - real(Dp) :: k(2) - - ! coordinates of R vector - real(Dp) :: R(2), R1(2), R2(2), R3(2) - -! Hamiltonian of bulk system - complex(Dp),intent(out) ::Hamk_bulk(Num_wann, Num_wann) - complex(dp), allocatable :: mat1(:, :) - complex(dp), allocatable :: mat2(:, :) - allocate(mat1(Num_wann, Num_wann)) - allocate(mat2(Num_wann, Num_wann)) - - nwann= Num_wann/2 - Hamk_bulk=0d0 - do iR=1, Nrpts - ia=irvec(1,iR) - ib=irvec(2,iR) - - R(1)=dble(ia) - R(2)=dble(ib) - - - do ia1= 1, Num_atoms - istart1= sum(nprojs(1:ia1-1))+1 - iend1= sum(nprojs(1:ia1)) - do ia2= 1, Num_atoms - istart2= sum(nprojs(1:ia2-1))+1 - iend2= sum(nprojs(1:ia2)) - R1= Atom_position_direct(:, ia1) - R2= Atom_position_direct(:, ia2) - R3= R+ R1- R2 - kdotr=k(1)*R3(1) + k(2)*R3(2) - - Hamk_bulk(istart1:iend1,istart2:iend2)=& - Hamk_bulk(istart1:iend1,istart2:iend2) & - +HmnR(istart1:iend1,istart2:iend2,iR)*Exp(2d0*pi*zi*kdotr)/ndegen(iR) - - if (soc>0) then - Hamk_bulk(istart1+nwann:iend1+nwann,istart2:iend2)=& - Hamk_bulk(istart1+nwann:iend1+nwann,istart2:iend2) & - +HmnR(istart1+nwann:iend1+nwann,istart2:iend2,iR)*Exp(2d0*pi*zi*kdotr)/ndegen(iR) - - Hamk_bulk(istart1+nwann:iend1+nwann,istart2+nwann:iend2+nwann)=& - Hamk_bulk(istart1+nwann:iend1+nwann,istart2+nwann:iend2+nwann) & - +HmnR(istart1+nwann:iend1+nwann,istart2+nwann:iend2+nwann,iR)*Exp(2d0*pi*zi*kdotr)/ndegen(iR) - - Hamk_bulk(istart1:iend1,istart2+nwann:iend2+nwann)=& - Hamk_bulk(istart1:iend1,istart2+nwann:iend2+nwann) & - +HmnR(istart1:iend1,istart2+nwann:iend2+nwann,iR)*Exp(2d0*pi*zi*kdotr)/ndegen(iR) - - endif - - enddo ! ia2 - enddo ! ia1 - enddo ! iR - - !call mat_mul(Num_wann, mirror_z, Hamk_bulk, mat1) - !call mat_mul(Num_wann, mat1, mirror_z, mat2) - !Hamk_bulk= (Hamk_bulk+ mat2)/2d0 - - ! 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,*)'Hamk_bulk is not hermitian' - write(stdout,*)i1, i2, Hamk_bulk(i1,i2), Hamk_bulk(i2,i1) - stop - endif - enddo - enddo - - return - end subroutine ham_bulk - - - -! This subroutine is used to caculate Hamiltonian for -! bulk system . - -! History -! May/29/2011 by Quansheng Wu - subroutine ham_bulk_old(k,Hamk_bulk) - - use para - implicit none - -! loop index - integer :: i1,i2,ia,ib,iR - integer :: nwann - - real(dp) :: kdotr - -! wave vector in 2d - real(Dp) :: k(2) - - ! coordinates of R vector - real(Dp) :: R(2), R1(2), R2(2) - -! Hamiltonian of bulk system - complex(Dp),intent(out) ::Hamk_bulk(Num_wann, Num_wann) - complex(dp), allocatable :: mat1(:, :) - complex(dp), allocatable :: mat2(:, :) - allocate(mat1(Num_wann, Num_wann)) - allocate(mat2(Num_wann, Num_wann)) - - Hamk_bulk=0d0 - do iR=1, Nrpts - ia=irvec(1,iR) - ib=irvec(2,iR) - - R(1)=dble(ia) - R(2)=dble(ib) - kdotr=k(1)*R (1) + k(2)*R (2) - - Hamk_bulk(:, :)=& - Hamk_bulk(:, :) & - + HmnR(:, :, iR)*Exp(2d0*pi*zi*kdotr)/ndegen(iR) - enddo ! iR - - !call mat_mul(Num_wann, mirror_z, Hamk_bulk, mat1) - !call mat_mul(Num_wann, mat1, mirror_z, mat2) - !Hamk_bulk= (Hamk_bulk+ mat2)/2d0 - - ! 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,*)'Hamk_bulk is not hermitian' - write(stdout,*)i1, i2, Hamk_bulk(i1,i2), Hamk_bulk(i2,i1) - stop - endif - enddo - enddo - - return - end diff --git a/2d/ham_qlayer2qlayer.f90 b/2d/ham_qlayer2qlayer.f90 deleted file mode 100644 index 6870412b..00000000 --- a/2d/ham_qlayer2qlayer.f90 +++ /dev/null @@ -1,181 +0,0 @@ -! this subroutine is used to caculate Hamiltonian between -! slabs -! 4/23/2010 by QS Wu -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - - subroutine ham_qlayer2qlayer(k,H00new,H01new) - - use para - - implicit none - -! loop index - integer :: i,j,iR - -! index used to sign irvec - integer :: ia,ib - -! - -! new index used to sign irvec - real(dp) :: new_ia,new_ib - integer :: inew_ib - -! wave vector k times lattice vector R - real(Dp) :: kdotr - -! input wave vector k's cooridinates - real(Dp),intent(in) :: k - -! H00 Hamiltonian between nearest neighbour-quintuple-layers -! the factor 2 is induced by spin - - complex(dp) :: ratio - - complex(Dp), allocatable :: Hij(:, :, :) - - -! H00 Hamiltonian between nearest neighbour-quintuple-layers -! the factor 2 is induced by spin - - -! complex(Dp),allocatable,intent(out) :: H00new(:,:) - complex(Dp),intent(out) :: H00new(Ndim,Ndim) - -! H01 Hamiltonian between next-nearest neighbour-quintuple-layers -! complex(Dp),allocatable,intent(out) :: H01new(:,:) - complex(Dp),intent(out) :: H01new(Ndim,Ndim) - - allocate(Hij(-ijmax:ijmax,Num_wann,Num_wann)) - - Hij=0.0d0 - do iR=1,Nrpts - ia=irvec(1,iR) - ib=irvec(2,iR) - - call latticetransform(ia, ib, new_ia, new_ib) - - inew_ib= int(new_ib) - if (abs(new_ib).le.ijmax)then - kdotr=k*new_ia - ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr) - - Hij(inew_ib, 1:Num_wann, 1:Num_wann )& - =Hij(inew_ib, 1:Num_wann, 1:Num_wann )& - +HmnR(:,:,iR)*ratio/ndegen(iR) - endif - - enddo - - H00new=0.0d0 - H01new=0.0d0 - -! nslab's principle layer -! H00new - do i=1,Np - do j=1,Np - if (abs(i-j).le.(ijmax)) then - H00new(Num_wann*(i-1)+1:Num_wann*i,Num_wann*(j-1)+1:Num_wann*j)& - =Hij(j-i,:,:) - endif - enddo - enddo - -! H01new - do i=1,Np - do j=Np+1,Np*2 - if (j-i.le.ijmax) then - H01new(Num_wann*(i-1)+1:Num_wann*i,& - Num_wann*(j-1-Np)+1:Num_wann*(j-Np))=Hij(j-i,:,:) - endif - enddo - enddo - - do i=1,Ndim - do j=1,Ndim - if(abs(H00new(i,j)-conjg(H00new(j,i))).ge.1e-4)then - write(stdout,*)'there are something wrong with ham_qlayer2qlayer' - stop - endif - - enddo - enddo - - return - end subroutine ham_qlayer2qlayer - - ! for slab - subroutine ham_qlayer2qlayer2(k,Hij) - - use para - - implicit none - -! loop index - integer :: iR - -! index used to sign irvec - integer :: ia,ib - -! - -! new index used to sign irvec - real(dp) :: new_ia,new_ib - integer :: inew_ib - -! wave vector k times lattice vector R - real(Dp) :: kdotr - -! input wave vector k's cooridinates - real(Dp),intent(in) :: k - -! H00 Hamiltonian between nearest neighbour-quintuple-layers -! the factor 2 is induced by spin - - complex(dp) :: ratio - - complex(Dp), intent(out) :: Hij(-ijmax:ijmax,Num_wann,Num_wann) - - Hij=0.0d0 - do iR=1,Nrpts - ia=irvec(1,iR) - ib=irvec(2,iR) - - !> new lattice - call latticetransform(ia, ib, new_ia, new_ib) - - inew_ib= int(new_ib) - if (abs(new_ib).le.ijmax)then - kdotr=k*new_ia - ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr) - - Hij(inew_ib, 1:Num_wann, 1:Num_wann )& - =Hij(inew_ib, 1:Num_wann, 1:Num_wann )& - +HmnR(:,:,iR)*ratio/ndegen(iR) - endif - - enddo - - return - end subroutine ham_qlayer2qlayer2 - - !> use Umatrix to get the new representation of a vector in new basis - !> R= a*R1+b*R2+c*R3= x*R1'+y*R2'+z*R3' - subroutine latticetransform(a, b, x, y) - use para - implicit none - - integer, intent(in) :: a, b - real(dp), intent(out) :: x, y - - real(dp) :: Uinv(2, 2) - - Uinv= Umatrix - - call inv_r(2, Uinv) - - x= a*Uinv(1, 1)+ b*Uinv(2, 1) - y= a*Uinv(1, 2)+ b*Uinv(2, 2) - - return - end subroutine latticetransform diff --git a/2d/ham_ribbon.f90 b/2d/ham_ribbon.f90 deleted file mode 100644 index 22150956..00000000 --- a/2d/ham_ribbon.f90 +++ /dev/null @@ -1,54 +0,0 @@ -! This subroutine is used to caculate Hamiltonian for -! ribbon system . - -! History -! 4/18/2010 by Quansheng Wu -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - subroutine ham_ribbon(k,Hamk_ribbon) - - use para - implicit none - - ! loop index - integer :: i1, i2, i, j - integer :: ir - - ! wave vector in 2d - real(Dp), intent(inout) :: k(2) - - ! Hamiltonian of ribbon system - complex(Dp),intent(out) ::Hamk_ribbon(Num_wann*nslab,Num_wann*nslab) - - ! the factor 2 is induced by spin - complex(Dp), allocatable :: Hij(:, :, :) - - allocate( Hij(-ijmax:ijmax,Num_wann,Num_wann)) - call ham_qlayer2qlayer2(k,Hij) - - Hamk_ribbon=0.0d0 - ! i1 column index - do i1=1, nslab - ! i2 row index - do i2=1, nslab - if (abs(i2-i1).le.ijmax)then - Hamk_ribbon((i2-1)*Num_wann+1:(i2-1)*Num_wann+Num_wann,& - (i1-1)*Num_wann+1:(i1-1)*Num_wann+Num_wann )& - = Hij(i2-i1,1:Num_wann,1:Num_wann) - endif - enddo ! i2 - enddo ! i1 - - ! check hermitcity - - do i1=1,nslab*Num_wann - do i2=1,nslab*Num_wann - if(abs(Hamk_ribbon(i1,i2)-conjg(Hamk_ribbon(i2,i1))).ge.1e-6)then - write(stdout,*)'there are something wrong with Hamk_ribbon' - stop - endif - enddo - enddo - - return - end subroutine ham_ribbon - diff --git a/2d/inverse.f90 b/2d/inverse.f90 deleted file mode 100644 index aa0d8426..00000000 --- a/2d/inverse.f90 +++ /dev/null @@ -1,113 +0,0 @@ - subroutine inv(ndim,Amat) - - - implicit none - - integer,parameter :: dp=8 - - integer :: i - integer :: info - - integer,intent(in):: ndim - -! IPIV : INTEGER. Array, DIMENSION at least max(1, n). The pivot indices that define -! the permutation matrix P; row i of the matrix was interchanged with -! row ipiv(i). Corresponds to the single precision factorization (if info= -! 0 and iter ≥ 0) or the double precision factorization (if info= 0 and -! iter < 0). - integer,allocatable :: ipiv(:) - - - complex(dp),parameter :: zone=(1.0d0,0.0d0) - - - -! Amat : -! Overwritten by the factors L and U from the factorization of A = P*L*U; -! the unit diagonal elements of L are not stored. -! If iterative refinement has been successfully used (info= 0 and iter≥ -! 0), then A is unchanged. -! If double precision factorization has been used (info= 0 and iter < -! 0), then the array A contains the factors L and U from the factorization -! A = P*L*U; the unit diagonal elements of L are not stored. - complex(dp),intent(inout):: Amat(ndim,ndim) - -! Bmat : -! Overwritten by the solution matrix X for dgesv, sgesv,zgesv,zgesv. - complex(dp),allocatable :: Bmat(:,:) - - - allocate(ipiv(ndim)) - allocate(Bmat(ndim,ndim)) - - ipiv=0 - - ! unit matrix - Bmat= (0d0, 0d0) - do i=1,ndim - Bmat(i,i)= zone - enddo - - call zgesv(ndim,ndim,Amat,ndim,ipiv,Bmat,ndim,info) - - if(info.ne.0)print *,'something wrong with zgesv' - - Amat=Bmat - - return - end subroutine inv - - !> get the inverse of a real matrix - subroutine inv_r(ndim,Amat) - - implicit none - - integer,parameter :: dp=8 - - integer :: i - integer :: info - - integer,intent(in):: ndim - -! IPIV : INTEGER. Array, DIMENSION at least max(1, n). The pivot indices that define -! the permutation matrix P; row i of the matrix was interchanged with -! row ipiv(i). Corresponds to the single precision factorization (if info= -! 0 and iter ≥ 0) or the double precision factorization (if info= 0 and -! iter < 0). - integer,allocatable :: ipiv(:) - - -! Amat : -! Overwritten by the factors L and U from the factorization of A = P*L*U; -! the unit diagonal elements of L are not stored. -! If iterative refinement has been successfully used (info= 0 and iter≥ -! 0), then A is unchanged. -! If double precision factorization has been used (info= 0 and iter < -! 0), then the array A contains the factors L and U from the factorization -! A = P*L*U; the unit diagonal elements of L are not stored. - real(dp),intent(inout):: Amat(ndim,ndim) - -! Bmat : -! Overwritten by the solution matrix X for dgesv, sgesv,zgesv,zgesv. - real(dp),allocatable :: Bmat(:,:) - - - allocate(ipiv(ndim)) - allocate(Bmat(ndim,ndim)) - - ipiv=0 - - ! unit matrix - Bmat= 0d0 - do i=1,ndim - Bmat(i,i)= 1d0 - enddo - - call dgesv(ndim,ndim,Amat,ndim,ipiv,Bmat,ndim,info) - - if(info.ne.0)print *,'something wrong with dgesv' - - Amat=Bmat - - return - end subroutine inv_r diff --git a/2d/main.f90 b/2d/main.f90 deleted file mode 100644 index cb7624f3..00000000 --- a/2d/main.f90 +++ /dev/null @@ -1,209 +0,0 @@ -!--------+--------+--------+--------+--------+--------+--------+------! -! main program of a set of tools based on Wannier90 TB -! constructed by Q.S.Wu on 4/9/2010 -! change by Q.S.Wu on 4/22/2010 -! changed by Q.S.wu on July/15/2010 -! version HmnR.data contains soc -! mpi-version is not test yet, take you own risk. -! mpi-version is tested , please report bugs to QSWU -! Jan 25 2015 by Q.S.Wu at ETH Zurich -! wuquansheng@gmail.com -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. -!--------+--------+--------+--------+--------+--------+--------+------! - - program main - - use wmpi - use para - implicit none - - !> file existence - logical :: exists - integer :: ierr - character(8) :: cht - - !> time measure - real(dp) :: time_start, time_end, time_init - - !> initial the environment of mpi - call mpi_init(ierr) - call mpi_comm_rank(mpi_cmw,cpuid,ierr) - call mpi_comm_size(mpi_cmw,num_cpu,ierr) - - if (cpuid==0) open(unit=stdout, file='WT.out') - - !> if mpi initial wrong, alarm - if (cpuid==0.and.ierr.ne.0)then - write(stdout,*)'mpi initialize wrong' - stop - endif - - call now(time_init) - call header - - !> print information for mpi - if (cpuid==0) then - write(stdout, '(1x, a, i5, a)')'You are using ', num_cpu, ' CPU cores' - write(stdout, *)' ' - endif - - !> readin the control parameters for this program - call readinput - - !> open file Hmn_R.data to get Num_wann and Nrpts - if (cpuid.eq.0)then - write(stdout,*)'' - inquire (file =Hrfile, EXIST = exists) - if (exists)then - if (index(Hrfile, 'HWR')==0) then - write(stdout,'(2x,a,a,a)')'File ',trim(Hrfile), & - ' exist, We are using HmnR from wannier90' - open(unit=1001,file=Hrfile,status='old') - read(1001,*) - read(1001,'(i)') Num_wann - read(1001,'(i)') Nrpts - write(stdout,*)'>> Num_wann', Num_wann - write(stdout,*)'>> NRPTS', NRPTS - close(1001) - else - write(stdout,'(2x, 3a)')'File ',trim(Hrfile), & - ' exist, We are using HmnR from HWR' - open(unit=1001,file=Hrfile,status='old') - read(1001,*) - read(1001,'(a26, i3)')cht, Num_wann - read(1001,'(a32,i10)')cht,Nrpts - write(stdout,*)'>> Num_wann', Num_wann - write(stdout,*)'>> NRPTS', NRPTS - close(1001) - endif ! hwr or not - else - write(stdout,'(2x,a25)')'>>> Error : no HmnR input' - stop - endif ! exists or not - if ((soc==0 .and. sum(nprojs)/=Num_wann) .or. & - (soc>0 .and. sum(nprojs)/=Num_wann/2))then - print *, sum(nprojs), num_wann, num_wann/2 - stop 'projectors are wrong' - endif - endif ! cpuid - - - !> broadcast and Nrpts to every cpu - call MPI_bcast(Num_wann,1,mpi_in,0,mpi_cmw,ierr) - call MPI_bcast(Nrpts,1,mpi_in,0,mpi_cmw,ierr) - - !> dimension for surface green's function - Ndim= Num_wann* Np - - - !> allocate necessary arrays for tight binding hamiltonians - allocate(irvec(2,nrpts)) - allocate(ndegen(nrpts)) - allocate(HmnR(num_wann,num_wann,nrpts)) - - if (cpuid==0)then - write(stdout,*) ' >> Begin to read Hmn_R.data' - endif - call readHmnR() - if (cpuid==0)then - write(stdout,*) ' << Read Hmn_R.data successfully' - endif - - - !> broadcast data to every cpu - call MPI_bcast(irvec,size(irvec),mpi_in,0,mpi_cmw,ierr) - call MPI_bcast(HmnR,size(HmnR),mpi_dc,0,mpi_cmw,ierr) - call MPI_bcast(ndegen,size(ndegen),mpi_in,0,mpi_cmw,ierr) - - !> bulk band - if (Dos_calc) then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating Dos' - call now(time_start) - call dos_joint_dos - call now(time_end) - call print_time_cost(time_start, time_end, 'Dos') - if(cpuid.eq.0)write(stdout, *)'<< End of calculating DOS' - endif - - !> FS - if (Fs_calc) then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating FS' - call now(time_start) - call fermisurface - call now(time_end) - call print_time_cost(time_start, time_end, 'FS') - if(cpuid.eq.0)write(stdout, *)'<< End of calculating FS' - endif - - !> FS - if (GapPlane_calc) then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating GapPlane' - call now(time_start) - call gapshape - call now(time_end) - call print_time_cost(time_start, time_end, 'GapPlane') - if(cpuid.eq.0)write(stdout, *)'<< End of calculating GapPlane' - endif - - - - - !> bulk band - if (BulkBand_calc) then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating bulk band' - call now(time_start) - call ek_bulk - call now(time_end) - call print_time_cost(time_start, time_end, 'BulkBand') - if(cpuid.eq.0)write(stdout, *)'<< End of calculating bulk band' - endif - - !> Ribbon band - if (RibbonBand_calc)then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the slab band structure' - call now(time_start) - call ek_ribbon - !call psik - call now(time_end) - call print_time_cost(time_start, time_end, 'SlabBand') - if(cpuid.eq.0)write(stdout, *)'<< End of calculating the slab band structure' - endif - - !> wannier center calculate - if (wanniercenter_calc)then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Wilson loop' - call now(time_start) - call wannier_center2D_plane - call now(time_end) - call print_time_cost(time_start, time_end, 'WannierCenter') - if(cpuid.eq.0)write(stdout, *)'<< End of calculating the Wilson loop' - endif - - !> surface state - if (SlabSS_calc) then - if(cpuid.eq.0)write(stdout, *)' ' - if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the surface state' - call now(time_start) - call surfstat - call now(time_end) - call print_time_cost(time_start, time_end, 'surfacestate') - if(cpuid.eq.0)write(stdout, *)'End of calculating the surface state' - endif - - call now(time_end) - - if(cpuid.eq.0)write(stdout, *)' ' - call print_time_cost(time_init, time_end, 'whole program') - if (cpuid.eq.0)write(stdout,*)'Congratulations! you finished the calculation.' - if (cpuid.eq.0)write(stdout,*)'See you next time :)' - - - call mpi_finalize(ierr) - - end !<< end of main program diff --git a/2d/mat_mul.f90 b/2d/mat_mul.f90 deleted file mode 100644 index 4901e168..00000000 --- a/2d/mat_mul.f90 +++ /dev/null @@ -1,133 +0,0 @@ -! performs matrix-matrix multiply -! C=A*B - subroutine mat_mul(nmatdim,A,B,C) - - use para, only : Dp - implicit none - - - integer,intent(in) :: nmatdim - - complex(Dp) :: ALPHA - complex(Dp) :: BETA - - - complex(Dp), intent(in) :: A(nmatdim ,nmatdim) - complex(Dp), intent(in) :: B(nmatdim ,nmatdim) - !complex(Dp) :: mat_mul(nmatdim,nmatdim) - complex(Dp), intent(out) :: C(nmatdim,nmatdim) - - ALPHA=1.0d0 - BETA=0.0D0 - - C(:,:)=(0.0d0,0.0d0) - - call ZGEMM('N','N',nmatdim,nmatdim,nmatdim,ALPHA, & - & A,nmatdim,B,nmatdim,BETA,C,nmatdim) - - return - end subroutine mat_mul - - !> ZGESVD computes the singular value decomposition (SVD) for GE matrices - !> In this pack, we assume the matrix A is a square matrix, the dimension - !> of row and column are the same - !> A = U * SIGMA * conjugate-transpose(V) - !> VT= conjugate-transpose(V) - subroutine zgesvd_pack(M, A, U, S, VT) - - use para, only : Dp - implicit none - - integer, intent(in) :: M - complex(dp), intent(inout) :: A(M, M) - complex(dp), intent(out) :: U(M, M) - real(dp) , intent(out) :: S(M, M) - complex(dp), intent(out) :: VT(M, M) - - character :: JOBU - character :: JOBVT - integer :: N - integer :: LDA - integer :: LDU - integer :: LDVT - integer :: LWORK - real(dp), allocatable :: RWORK(:) - complex(dp), allocatable :: WORK(:) - integer :: INFO - - N= M - LDA= M - LDU= M - LDVT= M - allocate(RWORK(5*M)) - - allocate(work(5*M)) - - JOBU= 'A' - JOBVT= 'A' - - LWORK = -1 - call zgesvd (JOBU, JOBVT, M, N, A, LDA, S, U, LDU, & - VT, LDVT, WORK, LWORK, RWORK, INFO) - if (INFO==0 .and. real(WORK(1))>0 )then - LWORK= WORK(1) - deallocate(work) - allocate(WORK(LWORK)) - else - write(*, *)'something wrong with zgesvd' - endif - - - call zgesvd (JOBU, JOBVT, M, N, A, LDA, S, U, LDU, & - VT, LDVT, WORK, LWORK, RWORK, INFO) - if (INFO /= 0) write(*, *)'something wrong with zgesvd' - - return - end subroutine zgesvd_pack - - subroutine zhpevx_pack(mat,ndim,eig,rot) - ! ! - ! Diagonalize the ndim x ndim hermitian matrix 'mat' and ! - ! return the eigenvalues 'eig' and the unitary rotation 'rot'! - ! ! - !============================================================! - - use para, only : dp, stdout - - integer, intent(in) :: ndim - complex(dp), intent(in) :: mat(ndim,ndim) - real(dp), intent(out) :: eig(ndim) - complex(dp), intent(out) :: rot(ndim,ndim) - - complex(dp), allocatable :: mat_pack(:),cwork(:) - real(dp), allocatable :: rwork(:) - integer :: i,j,info,nfound - integer, allocatable :: iwork(:),ifail(:) - - allocate(mat_pack((ndim*(ndim+1))/2)) - allocate(cwork(2*ndim)) - allocate(rwork(7*ndim)) - allocate(iwork(5*ndim)) - allocate(ifail(ndim)) - do j=1,ndim - do i=1,j - mat_pack(i+((j-1)*j)/2)=mat(i,j) - enddo - enddo - rot=0d0;eig=0.0_dp;cwork=0d0;rwork=0.0_dp;iwork=0 - call ZHPEVX('V','A','U',ndim,mat_pack,0.0_dp,0.0_dp,0,0,-1.0_dp, & - nfound,eig(1),rot,ndim,cwork,rwork,iwork,ifail,info) - if(info < 0) then - write(stdout,'(a,i3,a)') 'THE ',-info,& - ' ARGUMENT OF ZHPEVX HAD AN ILLEGAL VALUE' - stop 'Error in zhpevx_pack' - endif - if(info > 0) then - write(stdout,'(i3,a)') info,' EIGENVECTORS FAILED TO CONVERGE' - stop 'Error in zhpevx_pack' - endif - - return - end subroutine zhpevx_pack - - diff --git a/2d/module.f90 b/2d/module.f90 deleted file mode 100644 index 799c6cd8..00000000 --- a/2d/module.f90 +++ /dev/null @@ -1,225 +0,0 @@ -!> some global parameters -!> Copyright (c) 2010 QuanSheng Wu. All rights reserved. -!> add namelist for convenience June 5th 2016 by QuanSheng Wu - - module wmpi - include 'mpif.h' - end module wmpi - - - module para - - use wmpi - implicit none - - integer,parameter :: stdout= 8 - - character*80 :: Hrfile - namelist / TB_FILE / Hrfile - - !> control parameters - logical :: BulkBand_calc - logical :: RibbonBand_calc - logical :: SlabSS_calc - logical :: Dos_calc - logical :: FS_calc - logical :: GapPlane_calc - logical :: WannierCenter_calc - - namelist / Control / BulkBand_calc, & - RibbonBand_calc, & - SlabSS_calc, & - Dos_calc, & - FS_calc, & - GapPlane_calc, & - WannierCenter_calc - - - ! double precision - integer,parameter :: Dp=kind(1.0d0) - - ! number of slabs of Bi2Se3 - ! slab=1 means there is a quintuple layer of Bi2Se3 system - integer :: Nslab - - !> number of princple layers for surface green's function - integer :: Np - - integer, public, save :: ijmax=6 - - !> leading dimension of surface green's function - integer :: Ndim - - !> number of occupied bands for bulk unit cell - integer :: Numoccupied - - !> number of electrons - integer :: Ntotch - - integer :: Num_wann - - ! number of R points - integer :: Nrpts - - ! number of k points used in ek_slab - integer :: Nk - integer :: Nk1 - integer :: Nk2 - - integer, public, save :: Nr1=5 - integer, public, save :: Nr2=5 - - ! number of k points used for spintexture - integer,parameter :: kmesh(2)=(/200 , 200/) - integer,parameter :: knv=kmesh(1)*kmesh(2) - - - ! a parameter to control soc - ! Soc=0 means no spin-orbit coupling - ! Soc!=0 means no spin-orbit coupling - integer :: Soc - - ! used to calculate dos epsilon+i eta - real(Dp) :: eta - real(Dp) :: Eta_Arc - - ! the number of omega - integer :: OmegaNum - - ! omega interval - real(dp) :: OmegaMin, OmegaMax - - ! Fermi energy for arc calculation - real(Dp) :: E_arc - - ! threshold value for output the gap data for Gap3D - real(Dp) :: Gap_threshold - - !> namelist parameters - namelist /PARAMETERS/ Eta_Arc, OmegaNum, OmegaMin, OmegaMax, & - E_arc, Nk1, Nk2, NP, Gap_threshold - - ! Fermi energy - real(Dp) :: E_fermi - - !> surface onsite energy shift - real(dp) :: surf_onsite - - !> magnetic field (Tesla) - real(dp) :: Bx, By, Bz - - !> system parameters namelist - namelist / SYSTEM / Soc, E_fermi, Bx, By, Bz, surf_onsite, & - Nslab, Numoccupied, Ntotch - - !> e/2/h*a*a a=1d-10m, h is the planck constant - !> then the flux equals alpha*B*s - real(dp),parameter :: alpha= 1.20736d0*1D-6 - - ! circumference ratio pi - real(dp),parameter :: Pi= 3.14159265359d0 - real(dp),parameter :: half= 0.5d0 - real(dp),parameter :: zero= 0.0d0 - real(dp),parameter :: one = 1.0d0 - real(dp),parameter :: eps6= 1e-6 - real(dp),parameter :: eps9= 1e-9 - - real(Dp),parameter :: Ka(2)=(/1.0d0,0.0d0/) - real(Dp),parameter :: Kb(2)=(/0.0d0,1.0d0/) - - real(Dp),public, save :: Ra2(2) - real(Dp),public, save :: Rb2(2) - - real(Dp),public, save :: Ka2(2) - real(Dp),public, save :: Kb2(2) - - ! three primitive vectors in Cartsien coordinatec - real(dp),public, save :: Rua(2) - real(dp),public, save :: Rub(2) - - !> three primitive vectors in new coordinate system, see slab part - real(dp),public, save :: Rua_new(2) - real(dp),public, save :: Rub_new(2) - - ! three reciprocal primitive vectors - real(dp),public, save :: Kua(2) - real(dp),public, save :: Kub(2) - - real(dp),public, save :: Urot(2, 2) - - !> klist for 2D case include all 2D system - integer :: nk2lines - integer :: knv2 - real(dp) :: kp(2, 32) - real(dp) :: ke(2, 32) - real(dp) :: k2line_stop(32) - character(4) :: k2line_name(32) - real(dp),allocatable :: k2len(:) - real(dp),allocatable :: k2_path(:, :) - - !> kpoints plane for 2D system--> arcs - real(dp) :: K2D_start(2) - real(dp) :: K2D_vec1(2) - real(dp) :: K2D_vec2(2) - - !> klist for 1D case include all 1D system - integer :: nk1lines - integer :: knv1 - real(dp) :: kp1(32) - real(dp) :: ke1(32) - real(dp) :: k1line_stop(32) - character(4) :: k1line_name(32) - real(dp),allocatable :: k1len(:) - real(dp),allocatable :: k1_path(:) - - ! R coordinates - integer, allocatable :: irvec(:,:) - - ! Hamiltonian m,n are band indexes - complex(dp), allocatable :: HmnR(:,:,:) - - ! degree of degeneracy of R point - integer, allocatable :: ndegen(:) - - ! complex constant 0+1*i - complex(dp),parameter :: zi=(0.0d0, 1.0d0) - complex(dp),parameter :: pi2zi=(0.0d0, 6.283185307179586d0) - - integer :: cpuid - integer :: num_cpu - integer, parameter :: mpi_in= mpi_integer - integer, parameter :: mpi_dp= mpi_double_precision - integer, parameter :: mpi_dc= mpi_double_complex - integer, parameter :: mpi_cmw= mpi_comm_world - - !> a matrix change old primitive cell to new primitive cell - !> which can define new surface - !> a 2*2 matrix - real(dp), public, save :: Umatrix(2, 2) - - !> number of atoms in one primitive cell - integer :: Num_atoms - character(10) :: AngOrBohr - character(10) :: DirectOrCart - character(10), allocatable :: Atom_name(:) - real(dp), allocatable :: Atom_position(:, :) - real(dp), allocatable :: Atom_position_direct(:, :) - - integer :: max_projs - integer, allocatable :: nprojs(:) - character(10), allocatable :: proj_name(:, :) - - !> symmetry operator apply on function basis - complex(dp), allocatable :: inversion(:, :) - complex(dp), allocatable :: mirror_x(:, :) - complex(dp), allocatable :: mirror_z(:, :) - complex(dp), allocatable :: glide(:, :) - - !> symmetry operator apply on coordinate system - real(dp), allocatable :: inv_op(:, :) - real(dp), allocatable :: mirror_z_op(:, :) - real(dp), allocatable :: mirror_x_op(:, :) - real(dp), allocatable :: mirror_y_op(:, :) - real(dp), allocatable :: glide_y_op(:, :) - - end module para diff --git a/2d/psi.f90 b/2d/psi.f90 deleted file mode 100644 index 99fd8881..00000000 --- a/2d/psi.f90 +++ /dev/null @@ -1,145 +0,0 @@ -! This subroutine is used to caculate energy dispersion for -! slab Bi2Se3 - - subroutine psik() - - use para,only : Dp,Num_wann,Nslab, stdout, cpuid - - implicit none - -! loop index - integer :: i,j - integer :: info - -! wave vector - real(Dp) :: k(2) - -! parameters for zheev - real(Dp) :: rwork(5*Num_wann*nslab) - complex(Dp) :: work (5*Num_wann*nslab) - -! eigenvalue - real(Dp) :: eigenvalue (nslab*Num_wann) - - -! energy dispersion - real(Dp) :: ekslab(Nslab*Num_wann) - real(Dp) :: psi2 (Nslab) - complex(Dp) :: psi (Nslab*Num_wann) - -! hamiltonian slab - complex(Dp),allocatable :: CHamk(:,:) - complex(Dp),allocatable :: eigenvector(:,:) - - - allocate(CHamk(nslab*Num_wann,nslab*Num_wann)) - allocate(eigenvector(nslab*Num_wann,nslab*Num_wann)) - -! sweep k - ekslab=0.0d0 - -! Ka direction - k=0.0d0 - chamk=0.0d0 - -! calculate Hamiltonian - call ham_slab(k,Chamk) - eigenvalue=0.0d0 - eigenvector=Chamk - -! diagonal Chamk - call eigensystem_c('V', 'U', Num_wann*nslab, eigenvector, eigenvalue) - - ekslab=eigenvalue - - !> with given band index - info=18*Nslab+2 - - psi(:)=eigenvector(:, info) - if (cpuid.eq.0) write(stdout,*) 'eigenvalue',info,ekslab(info) - - j=0 - psi2=0.0d0 - do i=1,Nslab - do j=1,Num_wann/2 - psi2(i)=psi2(i)+abs(psi((i-1)*Num_wann+j))**2 - enddo - enddo - - open(unit=100, file='squarepsi.dat') - do i=1,Nslab - write(100,'(2f16.9)')real(i),psi2(i) - enddo - close(100) - - write(stdout,*) 'calculate psi done' - - end - -! This subroutine is used to caculate energy dispersion for -! bulk system - - subroutine psik_bulk() - - use para,only : Dp,Num_wann, stdout, Numoccupied, cpuid - implicit none - -! loop index - integer :: i,ik, ib - integer :: info - -! wave vector - real(Dp) :: k(3) - real(Dp) :: kpoints(3, 8) - -! eigenvalue - real(Dp) :: eigenvalue (Num_wann) - - -! energy dispersion - complex(Dp) :: psi (Num_wann) - -! hamiltonian slab - complex(Dp),allocatable :: CHamk(:,:) - complex(Dp),allocatable :: eigenvector(:,:) - - - allocate(CHamk(Num_wann,Num_wann)) - allocate(eigenvector(Num_wann,Num_wann)) - - - kpoints(:, 1)= (/0.0d0, 0.0d0, 0.0d0/) - kpoints(:, 2)= (/0.5d0, 0.0d0, 0.0d0/) - kpoints(:, 3)= (/0.0d0, 0.5d0, 0.0d0/) - kpoints(:, 4)= (/0.0d0, 0.0d0, 0.5d0/) - kpoints(:, 5)= (/0.5d0, 0.5d0, 0.5d0/) - kpoints(:, 6)= (/0.5d0, 0.5d0, 0.0d0/) - kpoints(:, 7)= (/0.0d0, 0.5d0, 0.5d0/) - kpoints(:, 8)= (/0.5d0, 0.0d0, 0.5d0/) - - ib= Numoccupied - if (cpuid==0)open(unit=100, file='wavefunction.dat') - do ik=1, 8 - - k=kpoints(:, ik) - ! calculate Hamiltonian - call ham_bulk(k,Chamk) - eigenvalue=0.0d0 - eigenvector=Chamk - - ! diagonal Chamk - call eigensystem_c('V', 'U', Num_wann, eigenvector, eigenvalue) - - psi(:)=eigenvector(:, ib) - if (cpuid==0)write(stdout,*) 'eigenvalue',info,eigenvalue(ib) - - if (cpuid==0)write(100, '(a,3f8.4)')'K point ', k - do i=1, Num_wann - if (cpuid==0)write(100,'(i, 30f16.9)')i, eigenvector(i, ib-1), eigenvector(i, ib) - enddo - if (cpuid==0)write(100,*)' ' - - enddo ! ik - - if (cpuid==0)close(100) - end subroutine psik_bulk diff --git a/2d/readHmnR.f90 b/2d/readHmnR.f90 deleted file mode 100644 index d312c419..00000000 --- a/2d/readHmnR.f90 +++ /dev/null @@ -1,99 +0,0 @@ -! read data from HmnR.data constructed by quansheng wu 4/2/2010 -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - - subroutine readHmnR() - - use para - - implicit none - - character*4 :: c_temp - -! file existenct - logical :: exists - - integer :: i, j, ir, ia - integer :: n, m - integer :: i1, i2, i3 - integer :: i4, i5 - integer :: nwan - real(dp) :: r1, r2 - -! real and imag of HmnR - real(dp) :: rh,ih - - if(cpuid.eq.0)write(stdout,*)'' - open(12, file=Hrfile) - if (index(Hrfile, 'HWR')==0) then - read(12, *) - read(12, *)nwan - read(12, *)nrpts - read(12,*)(ndegen(i), i=1, nrpts) - do ir=1, nrpts - do i=1, nwan - do j=1, nwan - read(12,*)i1, i2, i4, i5, r1, r2 - HmnR(i4,i5,ir)= dcmplx(r1, r2) ! in eV - !write(*,'(5i5,2f10.5)')i1, i2, i3, i4, i5, r1, r2 - enddo - enddo - irvec(1, ir)=i1 - irvec(2, ir)=i2 - enddo - - else - - !File *.HWR exist, We are using HmnR from WHM - ! skip 8 lines - do i=1,8 - read(12,*) - enddo - read(12,'(a11,f18.7)')c_temp,E_fermi - do iR=1,Nrpts - read(12,'(a3,2i5,a3,i4)')c_temp,irvec(1:2,iR),c_temp,ndegen(iR) - do i=1, Num_wann*Num_wann - read(12,*)n,m,rh,ih - HmnR(n,m,iR)=rh+ zi*ih ! in Hartree - enddo - if (sum(abs(irvec(:, ir)))==0) then - do i=1, Num_wann - HmnR(i,i,iR)= HmnR(i,i,iR)-E_fermi ! in Hartree - enddo - endif - enddo - HmnR= HmnR*27.2114d0 - - nwan= Num_wann - if (cpuid==0) then - open(unit=105, file='wannier90_hr.dat') - write(105, *)'hr file transformed from HWR' - write(105, *)Nwan - write(105, *)nrpts - write(105, '(15I5)')(ndegen(i), i=1, nrpts) - do ir=1, nrpts - do i=1, Nwan - do j=1, Nwan - write(105, '(4I5, 2f16.8)')irvec(:, ir), i, j, HmnR(i, j, ir) - enddo - enddo - enddo - close(105) - endif - - - - endif ! HWR or not - close(12) - - !call get_fermilevel - - do iR=1,Nrpts - if (Irvec(1,iR).eq.0.and.Irvec(2,iR).eq.0)then - do i=1, Num_wann - HmnR(i,i,iR)=HmnR(i,i,iR)-E_fermi - enddo - endif - enddo - - return - end subroutine readHmnR diff --git a/2d/readinput.f90 b/2d/readinput.f90 deleted file mode 100644 index 598685d4..00000000 --- a/2d/readinput.f90 +++ /dev/null @@ -1,544 +0,0 @@ -! -! this subroutine is used to read some paramters from -! input.dat -! constructed on 4/22/2010 by QS.Wu - - - subroutine readinput - - use wmpi - use para - implicit none - - character*12 :: fname='input.dat' - character*25 :: char_temp - character*80 :: inline - logical :: exists - logical :: lfound - real(dp) :: cell_volume - real(dp) :: cell_volume2 - - integer :: stat - integer :: i, ia, n - integer :: j - integer :: NN - integer :: nwann - real(dp) :: t1, temp - real(dp) :: pos(2) - real(dp) :: k1(2), k2(2) - real(dp) :: kstart(2), kend(2), kstart1, kend1 - real(dp) :: R1(2), R2(2) - real(dp), external :: norm - - - inquire(file=fname,exist=exists) - if (exists)then - if(cpuid==0)write(stdout,*) ' ' - if(cpuid==0)write(stdout,*) '>>>Read some paramters from input.dat' - open(unit=1001,file=fname,status='old') - else - if(cpuid==0)write(stdout,*)'file' ,fname, 'dosnot exist' - stop - endif - - !> inout file - - - - !read(1001,*)Hrfile - read(1001, TB_FILE, iostat= stat) - if (stat/=0) then - Hrfile='wannier90_hr.dat' - inquire(file='wannier90_hr.dat',exist=exists) - if (.not.exists) stop "TB_FIlE namelist should be given or wannier90_hr.dat should exist" - endif - if(cpuid==0)write(stdout,'(1x, a, a25)')"Tight binding Hamiltonian file: ",Hrfile - - - BulkBand_calc = .FALSE. - RibbonBand_calc = .FALSE. - SlabSS_calc = .FALSE. - Dos_calc = .FALSE. - FS_calc = .FALSE. - GapPlane_calc = .FALSE. - wanniercenter_calc = .FALSE. - - read(1001, CONTROL, iostat=stat) - - if (stat/=0) then - write(*, *)"ERROR: namelist CONTROL should be set" - write(*, *)"You should set one of these functions to be T" - write(*, *)"BulkBand_calc,BulkFS_calc,BulkGap_cube_calc,BulkGap_plane_calc" - write(*, *)"RibbonBand_calc,WireBand_calc,SlabSS_calc,SlabArc_calc " - write(*, *)"SlabSpintexture,wanniercenter_calc" - write(*, *)"BerryPhase_calc,BerryCurvature_calc" - write(*, *)"The default Vaule is F" - stop - endif - - !> control parameters - if (cpuid==0) then - write(stdout, *) " " - write(stdout, *) ">>>Control parameters: " - write(stdout, *) "BulkBand_calc : ", BulkBand_calc - write(stdout, *) "RibbonBand_calc : ", RibbonBand_calc - write(stdout, *) "SlabSS_calc : ", SlabSS_calc - write(stdout, *) "wanniercenter_calc : ", wanniercenter_calc - endif - - !> set system parameters by default - Nslab= 10 - Numoccupied = 0 - Ntotch = 0 - SOC = -1 - E_FERMI = 0 - Bx = 0 - By = 0 - Bz = 0 - surf_onsite = 0 - - !> read system parameters from file - read(1001, SYSTEM, iostat=stat) - if (stat/=0) then - write(*, *)"ERROR: namelist SYSTEM is wrong and should be set correctly" - !stop - endif - - if (SOC == -1) then - write(*, *)"ERROR: you should set SOC in namelist SYSTEM correctly" - stop - endif - - if (Numoccupied == 0) then - write(*, *)"ERROR: you should set Numoccupied in namelist SYSTEM correctly" - stop - endif - - if (Ntotch == 0) then - Ntotch = Numoccupied - endif - - if (cpuid==0) then - write(stdout, *) " " - write(stdout, *) ">>>System parameters: " - write(stdout, '(1x, a, i6 )')"NumSlabs :", Nslab - write(stdout, '(1x, a, i6)')"Number of Occupied bands:", NumOccupied - write(stdout, '(1x, a, i6)')"Number of total electrons:", Ntotch - write(stdout, '(1x, a, i6)')"With SOC or not in Hrfile:", SOC - write(stdout, '(1x, a, 3f16.6)')"Fermi energy :", E_FERMI - write(stdout, '(1x, a, 3f16.6)')"Bx, By, Bz :", Bx, By, Bz - write(stdout, '(1x, a, 3f16.6)')"surf_onsite :", surf_onsite - endif - - !> set up parameters for calculation - E_arc = 0.0d0 - Eta_Arc= 0.001d0 - OmegaNum = 100 - OmegaMin = -1d0 - OmegaMax = 1d0 - Nk1 = 50 - Nk2 = 50 - NP = 2 - Gap_threshold= 0.01d0 - - read(1001, PARAMETERS, iostat= stat) - - if (cpuid==0) then - write(stdout, *) " " - write(stdout, *) ">>>calculation parameters : " - write(stdout, '(1x, a, f16.5)')'E_arc : ', E_arc - write(stdout, '(1x, a, f16.5)')'Eta_arc : ', Eta_arc - write(stdout, '(1x, a, f16.5)')'Gap_threshold', Gap_threshold - write(stdout, '(1x, a, f16.5)')'OmegaMin : ', OmegaMin - write(stdout, '(1x, a, f16.5)')'OmegaMax : ', OmegaMax - write(stdout, '(1x, a, i6 )')'OmegaNum : ', OmegaNum - write(stdout, '(1x, a, i6 )')'Nk1 : ', Nk1 - write(stdout, '(1x, a, i6 )')'Nk2 : ', Nk2 - write(stdout, '(1x, a, i6 )')'NP number of principle layers : ', Np - endif - - NK = Nk1 - - !> read lattice information - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 100)inline - if (trim(adjustl(inline))=='LATTICE') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found LATTICE card' - exit - endif - enddo - 100 continue - - if (lfound) then - read(1001, *)inline ! The unit of lattice vector - AngOrBohr=trim(adjustl(inline)) - read(1001, *)Rua - read(1001, *)Rub - else - stop 'ERROR: please set lattice information' - endif - - if (index(AngOrBohr, 'Bohr')>0) then - Rua= Rua*0.529177d0 - Rub= Rub*0.529177d0 - endif - - !> transform lattice from direct space to reciprocal space - - Kua= 0d0 - Kub= 0d0 - cell_volume=Rua(1)*Rub(2)- Rub(1)*Rua(2) - cell_volume= abs(cell_volume) - cell_volume= 2d0*3.1415926535d0/cell_volume - - Kua(1)= cell_volume*Rub(2) - Kua(2)=-cell_volume*Rub(1) - - Kub(1)=-cell_volume*Rua(2) - Kub(2)= cell_volume*Rua(1) - - if(cpuid==0)write(stdout, '(a)') '>> lattice information (Angstrom)' - if(cpuid==0)write(stdout, '(2f12.6)')Rua - if(cpuid==0)write(stdout, '(2f12.6)')Rub - - if(cpuid==0)write(stdout, '(a)') '>> Reciprocal lattice information (1/Angstrom)' - if(cpuid==0)write(stdout, '(2f12.6)')Kua - if(cpuid==0)write(stdout, '(2f12.6)')Kub - - !> Read atom positions information - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 101)inline - if (trim(adjustl(inline))=='ATOM_POSITIONS') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found ATOM_POSITIONS card' - exit - endif - enddo - 101 continue - - if (lfound) then - read(1001, *)Num_atoms ! The unit of lattice vector - if(cpuid==0)write(stdout, '(a, i)')'Num_atoms', Num_atoms - allocate(atom_name(Num_atoms)) - allocate(Atom_position(2, Num_atoms)) - allocate(Atom_position_direct(2, Num_atoms)) - read(1001, *)inline ! The unit of lattice vector - DirectOrCart= trim(adjustl(inline)) - do i=1, Num_atoms - read(1001, *) atom_name(i), Atom_position(:, i) - if(cpuid==0)write(stdout, '(a4,2f12.6)')atom_name(i), Atom_position(:, i) - if (index(DirectOrCart, "D")>0)then - pos= Atom_position(:, i) - Atom_position(:, i)= pos(1)*Rua+ pos(2)*Rub - endif - enddo - if(cpuid==0)write(stdout,'(a)')'Atom position in cartisen coordinate' - do i=1, Num_atoms - if(cpuid==0)write(stdout, '(a4,2f12.6)')atom_name(i), Atom_position(:, i) - enddo - - if(cpuid==0)write(stdout,'(a)')'Atom position in direct coordinate' - do ia=1, Num_atoms - call cart_direct_real(Atom_position(:, ia), Atom_position_direct(:, ia)) - if(cpuid==0)write(stdout, '(a4,2f12.6)')atom_name(ia), Atom_position_direct(:, ia) - enddo - else - stop "ERROR: please set atom's positions information" - endif - - !> Read projectors information - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 102)inline - if (trim(adjustl(inline))=='PROJECTORS'.or.& - trim(adjustl(inline))=='PROJECTOR') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found PROJECTORS card' - exit - endif - enddo - 102 continue - - if (lfound) then - allocate(nprojs(Num_atoms)) - nprojs= 0 - read(1001, *)nprojs - if(cpuid==0)write(stdout, '(a, 40i5)')'nprojs', nprojs - - max_projs= maxval(nprojs) - allocate(proj_name(max_projs, Num_atoms)) - proj_name= ' ' - do i=1, Num_atoms - read(1001, *)char_temp, proj_name(1:nprojs(i), i) - if(cpuid==0)write(stdout, '(40a8)') & - char_temp, proj_name(1:nprojs(i), i) - enddo - else - stop "ERROR: please set projectors for Wannier functions information" - endif - - - !> read edge information - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 103)inline - if (trim(adjustl(inline))=='EDGE') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found EDGE card' - exit - endif - enddo - 103 continue - - if (.not.lfound) then - print *, inline - stop 'ERROR: please set edge information' - endif - - !> read information for new lattice - !> in order to get different edge state - !> R1'=U11*R1+U12*R2 - !> R2'=U21*R1+U22*R2 - read(1001, *)Umatrix(1, :) - read(1001, *)Umatrix(2, :) - - if (cpuid==0) then - write(stdout, '(a)')'>> new vectors to define the edge (in unit of lattice vector)' - write(stdout, '(a, 2f12.6)')' The 1st vector on edge :', Umatrix(1, :) - write(stdout, '(a, 2f12.6)')' The 2nd vector on edge :', Umatrix(2, :) - endif - - !> check whether Umatrix is right - !> the volume of the new cell should be the same as the old ones - R1= Umatrix(1, 1)*Rua+ Umatrix(1, 2)*Rub - R2= Umatrix(2, 1)*Rua+ Umatrix(2, 2)*Rub - - - cell_volume2= R1(1)*R2(2)- R1(2)*R2(1) - cell_volume2= 2d0*3.1415926535d0/cell_volume2 - - if (cell_volume2<0) then - R2=-R2 - Umatrix(2, :)= -Umatrix(2, :) - endif - - if (abs(abs(cell_volume2)-abs(cell_volume))> 0.001d0.and.cpuid==0) then - write(stdout, *)' ERROR The Umatrix is wrong, the new cell', & - 'volume should be the same as the old ones' - write(stdout, '(a,2f10.4)')'cell_volume vs cell_volume-new', cell_volume, cell_volume2 - stop - endif - - !> read kpath_slab information - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 105)inline - if (trim(adjustl(inline))=='KPATH_BULK') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found KPATH_BULK card' - exit - endif - enddo - - !> read in k lines for 2D system - k2line_name= ' ' - if (cpuid==0) write(stdout, *)'k lines for 2D system' - read(1001, *)nk2lines - if (cpuid==0) write(stdout, *)'No. of k lines for 2D system : ', nk2lines - do i=1, nk2lines - read(1001, *) k2line_name(i), kp(:, i), & - char_temp, ke(:, i) - if(cpuid==0) write(stdout, '(a6, 2f9.5, 4x, a6, 2f9.5)')& - k2line_name(i), kp(:, i), & - char_temp, ke(:, i) - enddo - k2line_name(nk2lines+1) = char_temp - - NN= Nk - knv2= NN*nk2lines - allocate( k2_path(knv2, 2)) - allocate( k2len (knv2)) - k2_path= 0d0 - k2len= 0d0 - - t1=0d0 - k2len=0d0 - k2line_stop= 0d0 - do j=1, nk2lines - do i=1, NN - kstart(1:2)= kp(:, j) - kend(1:2) = ke(:, j) - k1(1:2)= kstart(1)*Kua+ kstart(2)*Kub - k2(1:2)= kend(1)*Kua+ kend(2)*Kub - k2_path(i+(j-1)*NN,:)= kstart(1:2)+ (kend(1:2)-kstart(1:2))*(i-1)/dble(NN-1) - - temp= dsqrt((k2(1)- k1(1))**2 & - +(k2(2)- k1(2))**2)/dble(NN-1) - - if (i.gt.1) then - t1=t1+temp - endif - k2len(i+(j-1)*NN)= t1 - enddo - k2line_stop(j+1)= t1 - - enddo - - - 105 continue - - !> read kpath_ribbon information - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 116)inline - if (trim(adjustl(inline))=='KPATH_RIBBON') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found KPATH_RIBBON card' - exit - endif - enddo - - !> read in k lines for 1D system - k1line_name= ' ' - if (cpuid==0) write(stdout, *)'k lines for 1D system' - read(1001, *)nk1lines - if (cpuid==0) write(stdout, *)'No. of k lines for 1D system : ', nk1lines - do i=1, nk1lines - read(1001, *) k1line_name(i), kp1(i), & - char_temp, ke1(i) - if(cpuid==0) write(stdout, '(a6, 1f9.5, 4x, a6, 1f9.5)')& - k1line_name(i), kp1(i), & - char_temp, ke1(i) - enddo - k1line_name(nk1lines+1) = char_temp - - NN= Nk - knv1= NN*nk1lines - allocate( k1_path(knv2)) - allocate( k1len (knv2)) - k1_path= 0d0 - k1len= 0d0 - - t1=0d0 - k1len=0d0 - k1line_stop= 0d0 - do j=1, nk1lines - do i=1, NN - kstart1= kp1(j) - kend1 = ke1(j) - k1_path(i+(j-1)*NN)= kstart1+ (kend1-kstart1)*(i-1)/dble(NN-1) - - temp= abs((kend1- kstart1))/dble(NN-1) - - if (i.gt.1) then - t1=t1+temp - endif - k1len(i+(j-1)*NN)= t1 - enddo - k1line_stop(j+1)= t1 - - enddo - - - 116 continue - if (.not.lfound .and.RibbonBand_calc) then - stop 'ERROR: please set KPATH_SLAB for slab band structure calculation' - endif - - - !> read kplane_slab information - !> default value for KPLANE_SLAB - K2D_start= (/-0.5, -0.5/) - K2D_vec1 = (/ 1.0, 0.0/) - K2D_vec2 = (/ 0.0, 1.0/) - - rewind(1001) - lfound = .false. - do while (.true.) - read(1001, *, end= 106)inline - if (trim(adjustl(inline))=='KPLANE_BULK') then - lfound= .true. - if (cpuid==0) write(stdout, *)' ' - if (cpuid==0) write(stdout, *)'We found KPLANE_BULK card' - exit - endif - enddo - - !> kpoints plane for 2D system--> arcs - read(1001, *)K2D_start - read(1001, *)K2D_vec1 - read(1001, *)K2D_vec2 - 106 continue - - if (cpuid==0) write(stdout, *)'>> Kpoints plane for 2D system--> arcs ' - if (cpuid==0) write(stdout, '((a, 2f8.4))')'K2D_start:', K2D_start - if (cpuid==0) write(stdout, '((a, 2f8.4))')'The first vector: ', K2D_vec1 - if (cpuid==0) write(stdout, '((a, 2f8.4))')'The second vector: ', K2D_vec2 - - - !> close input.dat - close(1001) - - eta=(omegamax- omegamin)/omeganum*2d0 - - if(cpuid==0)write(stdout,*)'<<> local variables - ! iteration number - integer :: iter - - ! maximun iteration - integer ,parameter:: itermax=100 - - ! accuracy control - real(Dp) :: accuracy=1e-16 - - ! a real type temp variable - real(Dp) :: real_temp - - ! omegac=omega(i)+I * eta - complex(Dp) :: omegac - - - ! some variables in Eq.(11) - complex(Dp), allocatable :: alphai(:, :) - complex(Dp), allocatable :: betai(:, :) - complex(Dp), allocatable :: epsiloni(:, :) - complex(Dp), allocatable :: epsilons(:, :) - complex(Dp), allocatable :: epsilons_t(:, :) - - complex(Dp), allocatable :: mat1 (:, :) - complex(Dp), allocatable :: mat2 (:, :) - - ! g0= inv(w-e_i) - complex(Dp), allocatable :: g0 (:, :) - - ! allocate some variables - allocate(alphai(Ndim, Ndim)) - allocate(betai (Ndim, Ndim)) - allocate(epsiloni (Ndim, Ndim)) - allocate(epsilons (Ndim, Ndim)) - allocate(epsilons_t(Ndim, Ndim)) - allocate(mat1(Ndim, Ndim)) - allocate(mat2(Ndim, Ndim)) - allocate(g0(Ndim, Ndim)) - - epsiloni= H00 - epsilons= H00 - epsilons_t= H00 - alphai = H01 - betai = conjg(transpose(H01)) - !print *, sqrt(sum(abs(H00)**2)), 'H00' - - ! w+i*0^+ - omegac= dcmplx(omega, eta) - !print *, omegac - - ! begin iteration - do iter=1, itermax - - g0= omegac*ones- epsiloni - call inv(Ndim, g0, ones) - - ! a_i-1*(w-e_i-1)^-1 - call mat_mul(Ndim, alphai, g0, mat1 ) - - ! b_i-1*(w-e_i-1)^-1 - call mat_mul(Ndim, betai, g0, mat2 ) - - ! a_i-1*(w-e_i-1)^-1*b_i-1 - call mat_mul(Ndim, mat1, betai, g0) - epsiloni= epsiloni+ g0 - !print *, sqrt(sum(abs(epsiloni)**2)), 'ei' - ! es_i= es_i-1 + a_i-1*(w-e_i-1)^-1*b_i-1 - epsilons= epsilons+ g0 - !print *, sqrt(sum(abs(epsilons)**2)), 'es' - !pause - - ! b_i-1*(w-e_i-1)^-1*a_i-1 - call mat_mul(Ndim, mat2, alphai, g0) - epsiloni= epsiloni+ g0 - ! es_i= es_i-1 + a_i-1*(w-e_i-1)^-1*b_i-1 - epsilons_t= epsilons_t+ g0 - - ! a_i= a_i-1*(w-e_i-1)^-1*a_i-1 - call mat_mul(Ndim, mat1, alphai, g0) - alphai= g0 - ! b_i= b_i-1*(w-e_i-1)^-1*b_i-1 - call mat_mul(Ndim, mat2, betai, g0) - betai= g0 - - !real_temp=maxval(abs(alphai)) - real_temp=sum(abs(alphai)) - !if (cpuid.eq.0) print *, iter, real_temp - if (real_temp.le.accuracy) exit - - enddo ! end of iteration - - ! calculate surface green's function - GLL= omegac*ones- epsilons - call inv(Ndim, GLL, ones) - - GRR= omegac*ones- epsilons_t - call inv(Ndim, GRR, ones) - - return - end subroutine surfgreen_1985 - - -!+---------+---------+---------+---------+---------+---------+--------+! -! this subroutine is used to calculate surface state using ! -! green's function method --- J.Phys.F.Met.Phys.14(1984)1205-1215 ! -! Quick iterative scheme for the calculation of transfer matrices: -! application to Mo (100) -! History: -! by Quan Sheng Wu on 4/20/2010 ! -! mpi version 4/21/2010 -! little change to cpu version in Zurich Swiss Jan 25 2015 -!+---------+---------+---------+---------+---------+---------+--------+! - subroutine surfgreen_1984(omega,GLL,GRR,H00,H01,ones) - - use wmpi - use para - implicit none - - - ! general loop index - integer :: i,j - - ! iteration loop index - integer :: it - - ! iteration number - integer :: iter - - ! maximun iteration - integer ,parameter :: itermax=100 - - ! accuracy control - real(Dp),parameter :: accuracy=1e-6 - - ! a real type temp variable - real(Dp) :: real_temp - - ! frequency - real(Dp),intent(in) :: omega - - ! energy energy=omega(i)+I * eta - complex(Dp) :: energy - - ! surface green fuction an output variable - complex(Dp),intent(out) :: GLL(Ndim,Ndim) - complex(Dp),intent(out) :: GRR(Ndim,Ndim) - - ! H00 Hamiltonian between nearest neighbour-quintuple-layers - ! the factor 2 is induced by spin - complex(Dp),intent(in) :: H00(Ndim,Ndim) - - ! H01 Hamiltonian between next-nearest neighbour-quintuple-layers - complex(Dp),intent(in) :: H01(Ndim,Ndim) - - ! unit matrix ones - complex(Dp),intent(in) :: ones(Ndim,Ndim) - - complex(Dp),allocatable :: H01dag(:,:) - - !> surface hamiltonian - complex(Dp),allocatable :: Hs(:,:) - - ! temp hamiltonian - complex(Dp),allocatable :: t0(:,:) - complex(Dp),allocatable :: Tmatrix(:,:) - complex(Dp),allocatable :: Tmatrixt(:,:) - complex(Dp),allocatable :: t0tilde(:,:) - complex(Dp),allocatable :: tnew(:,:) - complex(Dp),allocatable :: tnewtilde(:,:) - complex(Dp) ,allocatable :: told(:,:) - complex(Dp),allocatable :: toldtilde(:,:) - complex(Dp),allocatable :: temp(:,:) - complex(Dp),allocatable :: Tmat_temp(:,:) - complex(Dp),allocatable :: Tmat_tempt(:,:) - - real(Dp),allocatable :: abs_told(:,:) - - ! allocate variables - allocate(Hs(ndim,ndim)) - allocate(t0(ndim,ndim)) - allocate(Tmatrix(ndim,ndim)) - allocate(Tmatrixt(ndim,ndim)) - allocate(t0tilde(ndim,ndim)) - allocate(tnew(ndim,ndim)) - allocate(tnewtilde(ndim,ndim)) - allocate(told(ndim,ndim)) - allocate(toldtilde(ndim,ndim)) - allocate(temp(ndim,ndim)) - allocate(Tmat_temp(ndim,ndim)) - allocate(Tmat_tempt(ndim,ndim)) - allocate(abs_told(ndim,ndim)) - - allocate(H01dag(ndim,ndim)) - - Hs=0.0d0 - t0=0.0d0 - t0tilde=0.0d0 - told=0.0d0 - toldtilde=0.0d0 - tnew=0.0d0 - tnewtilde=0.0d0 - Tmatrix=0.0d0 - Tmatrixt=0.0d0 - temp=0.0d0 - abs_told=0.0d0 - Tmat_temp=0.0d0 - GLL=0d0 - GRR=0d0 - Tmat_tempt=0.0d0 - H01dag=0.0d0 - - ! H01dag=H01^dag - do i=1,ndim - do j=1,ndim - H01dag(i,j)=conjg(H01(j,i)) - enddo - enddo - - energy=omega+zi*eta - temp=energy*ones-H00 - - - call inv(ndim,temp) - t0= matmul(temp,H01dag) - t0tilde= matmul(temp,H01) - told=t0 - toldtilde=t0tilde - - ! begin iteration - iter=0 - Tmatrix=0.0d0 - Tmat_temp=ones - Tmatrixt=0.0d0 - Tmat_tempt=ones - ITER1 : do it=1,itermax - iter=iter+1 - - temp=ones-matmul(told,toldtilde)-matmul(toldtilde,told) - call inv(ndim,temp) - - tnew=matmul(temp,told) - tnew=matmul(tnew,told) - tnewtilde=matmul(temp,toldtilde) - tnewtilde=matmul(tnewtilde,toldtilde) - - Tmat_temp=matmul(Tmat_temp,toldtilde) - Tmatrix=Tmatrix+matmul(Tmat_temp,tnew) - - Tmat_tempt=matmul(Tmat_tempt,told) - Tmatrixt=Tmatrixt+matmul(Tmat_tempt,tnewtilde) - - told=tnew - toldtilde=tnewtilde - - do i=1,ndim - do j=1,ndim - abs_told(i,j)=abs(told(i,j))+abs(toldtilde(i,j)) - enddo - enddo - real_temp=maxval(abs_told) - if (real_temp.le.accuracy) exit ITER1 - - end do ITER1 - - !print *,'iter,acc',iter,real_temp - - !> set up surface hamiltonian - !> usually Hs= H00 - !> but you can add static potential on the surface - Hs= H00 - do i=1, ndim - Hs(i, i)=Hs(i, i)+ surf_onsite - enddo - - Tmatrix=t0+Tmatrix - temp=energy*ones-Hs-matmul(H01,Tmatrix) - call inv(ndim,temp) - - ! g_00=(epsilon-Hs -H01*T)^-1 - GLL(1:ndim,1:ndim)=temp - - Tmatrixt=t0tilde+Tmatrixt - temp=energy*ones-Hs -matmul(conjg(transpose(H01)),Tmatrixt) - call inv(ndim,temp) - - ! g_00=(epsilon-Hs -H01*T)^-1 - GRR(1:ndim,1:ndim)=temp - - return - end subroutine surfgreen_1984 diff --git a/2d/surfstat.f90 b/2d/surfstat.f90 deleted file mode 100644 index 51281e6a..00000000 --- a/2d/surfstat.f90 +++ /dev/null @@ -1,228 +0,0 @@ -!+---------+---------+---------+---------+---------+---------+--------+! -! this subroutine is used to calculate surface state using ! -! green's function method --- J.Phys.F.Met.Phys.15(1985)851-858 ! -! -! History: -! by Quan Sheng Wu on 4/20/2010 ! -! mpi version 4/21/2010 -! 2d version 8/02/2016 -! change Kb to K=(Ka+Kb)/3 direction 4/22/2010 -! Quansheng Wu on Jan 30 2015 at ETH Zurich -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. -!+---------+---------+---------+---------+---------+---------+--------+! - subroutine surfstat - - use wmpi - use para - implicit none - - integer :: ierr - - ! general loop index - integer :: i,j - - ! kpoint loop index - integer :: ikp - - real(dp) :: emin - real(dp) :: emax - real(dp) :: k, w - - real(dp), allocatable :: omega(:) - - real(dp), allocatable :: dos_l(:,:) - real(dp), allocatable :: dos_r(:,:) - real(dp), allocatable :: dos_l_mpi(:,:) - real(dp), allocatable :: dos_r_mpi(:,:) - - complex(dp), allocatable :: GLL(:,:) - complex(dp), allocatable :: GRR(:,:) - complex(dp), allocatable :: H00(:,:) - complex(dp), allocatable :: H01(:,:) - complex(dp), allocatable :: ones(:,:) - - - allocate( omega(omeganum)) - allocate( dos_l(knv1, omeganum)) - allocate( dos_r(knv1, omeganum)) - allocate( dos_l_mpi(knv1, omeganum)) - allocate( dos_r_mpi(knv1, omeganum)) - omega=0d0 - dos_l=0d0 - dos_r=0d0 - dos_l_mpi=0d0 - dos_r_mpi=0d0 - - eta=(omegamax- omegamin)/dble(omeganum)*1.5d0 - - do i= 1, omeganum - omega(i)=omegamin+(i-1)*(omegamax-omegamin)/dble(omeganum) - enddo - - allocate(GLL(Ndim, Ndim)) - allocate(GRR(Ndim, Ndim)) - allocate(H00(Ndim, Ndim)) - allocate(H01(Ndim, Ndim)) - allocate(ones(Ndim, Ndim)) - GLL= 0d0 - GRR= 0d0 - H00= 0d0 - H01= 0d0 - ones= 0d0 - - do i=1,Ndim - ones(i,i)=1.0d0 - enddo - - do ikp= 1+cpuid, knv1, num_cpu - if (cpuid==0) write(stdout, *) 'SurfaceSS, ik', ikp, 'Nk', knv1 - k= k1_path(ikp) - - !> get the hopping matrix between two principle layers - call ham_qlayer2qlayer(k,H00,H01) - - do j = 1, omeganum - w=omega(j) - - !> calculate surface green function - ! 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,H00,H01,ones) - !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)) - enddo ! i - enddo ! j - enddo ! ikp - - !> we don't have to do allreduce operation - call mpi_reduce(dos_l, dos_l_mpi, size(dos_l), mpi_double_precision,& - mpi_sum, 0, mpi_comm_world, ierr) - call mpi_reduce(dos_r, dos_r_mpi, size(dos_r), mpi_double_precision,& - mpi_sum, 0, mpi_comm_world, ierr) - dos_l=log(dos_l_mpi) - dos_r=log(dos_r_mpi) - - if (cpuid.eq.0)then - open (unit=12, file='dos.dat_l') - open (unit=13, file='dos.dat_r') - do ikp=1, knv1 - do j=1, omeganum - write(12, '(3f16.8)')k1len(ikp), omega(j), dos_l(ikp, j) - write(13, '(3f16.8)')k1len(ikp), omega(j), dos_r(ikp, j) - enddo - write(12, *) ' ' - write(13, *) ' ' - enddo - close(12) - close(13) - write(stdout,*)'ndim',ndim - write(stdout,*) 'knv1,omeganum,eta',knv1, omeganum, eta - write(stdout,*)'calculate density of state successfully' - endif - - emin= minval(omega) - emax= maxval(omega) - !> write script for gnuplot - if (cpuid==0) then - open(unit=116, file='surfdos_l.gnu') - write(116, '(a)')"set encoding iso_8859_1" - write(116, '(a)')'#set terminal postscript enhanced color' - write(116, '(a)')"#set output 'surfdos_l.eps'" - write(116, '(3a)')'#set terminal pngcairo truecolor enhanced', & - ' font ", 60" size 1920, 1680' - write(116, '(3a)')'set terminal png truecolor enhanced', & - ' font ", 60" size 1920, 1680' - write(116, '(a)')"set output 'surfdos_l.png'" - write(116,'(2a)') 'set palette defined (-10 "#194eff", ', & - '0 "white", 10 "red" )' - write(116, '(a)')'#set palette rgbformulae 33,13,10' - write(116, '(a)')'set style data linespoints' - write(116, '(a)')'set size 0.8, 1' - write(116, '(a)')'set origin 0.1, 0' - write(116, '(a)')'unset ztics' - write(116, '(a)')'unset key' - write(116, '(a)')'set pointsize 0.8' - write(116, '(a)')'set pm3d' - write(116, '(a)')'#set view equal xyz' - write(116, '(a)')'set view map' - write(116, '(a)')'set border lw 3' - write(116, '(a)')'#set cbtics font ",48"' - write(116, '(a)')'#set xtics font ",48"' - write(116, '(a)')'#set ytics font ",48"' - write(116, '(a)')'#set ylabel font ",48"' - write(116, '(a)')'set ylabel "Energy (eV)"' - write(116, '(a)')'#set xtics offset 0, -1' - write(116, '(a)')'#set ylabel offset -6, 0 ' - write(116, '(a, f8.5, a)')'set xrange [0: ', maxval(k1len), ']' - write(116, '(a, f8.5, a, f8.5, a)')'set yrange [', emin, ':', emax, ']' - write(116, 202, advance="no") (k1line_name(i), k1line_stop(i), i=1, nk1lines) - write(116, 203)k1line_name(nk1lines+1), k1line_stop(nk1lines+1) - - do i=1, nk1lines-1 - write(116, 204)k1line_stop(i+1), emin, k1line_stop(i+1), emax - enddo - write(116, '(a)')'set pm3d interpolate 2,2' - write(116, '(2a)')"splot 'dos.dat_l' u 1:2:3 w pm3d" - close(116) - - endif - - !> write script for gnuplot - if (cpuid==0) then - open(unit=117, file='surfdos_r.gnu') - write(117, '(a)')"set encoding iso_8859_1" - write(117, '(a)')'#set terminal postscript enhanced color' - write(117, '(a)')"#set output 'surfdos_r.eps'" - write(117, '(3a)')'#set terminal pngcairo truecolor enhanced', & - ' font ", 60" size 1920, 1680' - write(117, '(3a)')'set terminal png truecolor enhanced', & - ' font ", 60" size 1920, 1680' - write(117, '(a)')"set output 'surfdos_r.png'" - write(117,'(2a)') 'set palette defined (-10 "#194eff", ', & - '0 "white", 10 "red" )' - write(117, '(a)')'#set palette rgbformulae 33,13,10' - write(117, '(a)')'set style data linespoints' - write(117, '(a)')'unset ztics' - write(117, '(a)')'unset key' - write(117, '(a)')'set pointsize 0.8' - write(117, '(a)')'set pm3d' - write(117, '(a)')'set border lw 3' - write(117, '(a)')'set size 0.8, 1' - write(117, '(a)')'set origin 0.1, 0' - write(117, '(a)')'#set size ratio -1' - write(117, '(a)')'#set view equal xyz' - write(117, '(a)')'set view map' - write(117, '(a)')'#set cbtics font ",48"' - write(117, '(a)')'#set xtics font ",48"' - write(117, '(a)')'#set ytics font ",48"' - write(117, '(a)')'#set ylabel font ",48"' - write(117, '(a)')'set ylabel "Energy (eV)"' - write(117, '(a)')'#set xtics offset 0, -1' - write(117, '(a)')'#set ylabel offset -6, 0 ' - write(117, '(a, f8.5, a)')'set xrange [0: ', maxval(k1len), ']' - write(117, '(a, f8.5, a, f8.5, a)')'set yrange [', emin, ':', emax, ']' - write(117, 202, advance="no") (k1line_name(i), k1line_stop(i), i=1, nk1lines) - write(117, 203)k1line_name(nk1lines+1), k1line_stop(nk1lines+1) - - do i=1, nk1lines-1 - write(117, 204)k1line_stop(i+1), emin, k1line_stop(i+1), emax - enddo - write(117, '(a)')'set pm3d interpolate 2,2' - write(117, '(2a)')"splot 'dos.dat_r' u 1:2:3 w pm3d" - close(117) - - endif - - - 202 format('set xtics (',:20('"',A1,'" ',F8.5,',')) - 203 format(A1,'" ',F8.5,')') - 204 format('set arrow from ',F8.5,',',F10.5,' to ',F8.5,',',F10.5, ' nohead front lw 3') - - - return - end subroutine diff --git a/2d/surfstat2.f90 b/2d/surfstat2.f90 deleted file mode 100644 index d02f6504..00000000 --- a/2d/surfstat2.f90 +++ /dev/null @@ -1,377 +0,0 @@ -!+---------+---------+---------+---------+---------+---------+--------+! -! this subroutine is used to calculate surface state using ! -! green's function method --- J.Phys.F.Met.Phys.15(1985)851-858 ! -! -! History: -! by Quan Sheng Wu on 4/20/2010 ! -! mpi version 4/21/2010 -! change Kb to K=(Ka+Kb)/3 direction 4/22/2010 -! Quansheng Wu on Jan 30 2015 at ETH Zurich -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. -!+---------+---------+---------+---------+---------+---------+--------+! - subroutine surfstat2 - - use wmpi - use para - implicit none - - integer :: ierr - - ! general loop index - integer :: i,j - - integer :: knv2 - - ! kpoint loop index - integer :: ikp - - integer :: NN, nlines - - integer :: Nwann - - real(dp) :: emin - real(dp) :: emax - real(dp) :: t1, temp - real(dp) :: rx0, ry0, x0, y0 - real(dp) :: k(2), w, s - real(dp) :: k1(2) - real(dp) :: k2(2) - - real(dp) :: kp(16, 2) - real(dp) :: ke(16, 2) - real(dp) :: kpath_stop(16) - character(4) :: kpath_name(17) - - real(dp) :: kstart(2), kend(2) - - real(dp), allocatable :: kpoint(:,:) - real(dp), allocatable :: omega(:) - - real(dp), allocatable :: dos_l(:, :,:) - real(dp), allocatable :: dos_r(:, :,:) - real(dp), allocatable :: dos_l_mpi(:, :,:) - real(dp), allocatable :: dos_r_mpi(:, :,:) - - complex(dp), allocatable :: GLL(:,:) - complex(dp), allocatable :: GRR(:,:) - complex(dp), allocatable :: H00(:,:) - complex(dp), allocatable :: H01(:,:) - complex(dp), allocatable :: ones(:,:) - - real(dp), allocatable :: k_len(:) - - kpath_name= ' ' - kp(1,:)=(/0.5d0, 0.00d0/) ; kpath_name(1)= 'X' - ke(1,:)=(/0.0d0, 0.00d0/) - kp(2,:)=(/0.0d0, 0.00d0/) ; kpath_name(2)= 'G' - ke(2,:)=(/0.0d0, 0.50d0/) ! K - kp(3,:)=(/0.0d0, 0.50d0/) ; kpath_name(3)= 'Y' - ke(3,:)=(/0.5d0, 0.5d0/) ! K - kp(4,:)=(/0.5d0, 0.5d0/) ; kpath_name(4)= 'M' - ke(4,:)=(/0.0d0, 0.0d0/) ; kpath_name(5)= 'G' - - !kp(1,:)=(/0.5d0, 0.00d0/) ; kpath_name(1)= 'X' - !ke(1,:)=(/0.0d0, 0.00d0/) - !kp(1,:)=(/0.266666d0, 0.2666660d0/) ; kpath_name(1)= 'G' - !ke(1,:)=(/0.311111d0, 0.3111111d0/) ! K - !kp(2,:)=(/0.333333d0, 0.3333330d0/) ; kpath_name(2)= 'K' - !ke(2,:)=(/0.5d0, 0.0d0/) ! K - !kp(3,:)=(/0.5d0, 0.0d0/) ; kpath_name(3)= 'Y' - !ke(3,:)=(/0.333333d0, 0.3333330d0/) ! K - !kp(5,:)=(/0.333333d0, 0.3333330d0/) ! K - !ke(5,:)=(/0.0d0, 0.0d0/) ; kpath_name(5)= 'G' - - - !kp(1,:)=(/0.2d0, 0.2d0/) ! Gamma - !ke(1,:)=(/0.5d0, 0.5d0/) ! Z - - !kp(1,:)=(/-0.5d0, 0.00d0/) ; kpath_name(1)= 'A' - !ke(1,:)=(/0.0d0, 0.00d0/) - !kp(2,:)=(/0.0d0, 0.00d0/) ; kpath_name(2)= 'G' - !ke(2,:)=(/0.5d0, 0.00d0/) ; kpath_name(3)= 'A' - - !kp(1,:)=(/0.3d0, 0.00d0/) ; kpath_name(1)= 'X' - !ke(1,:)=(/0.0d0, 0.00d0/) - !kp(1, 1)= surf_onsite - !ke(1, 2)= surf_onsite - !kp(2, :)=(/0.0d0, 0.00d0/) ; kpath_name(2)= 'G' - - !kp(1,:)=(/-0.5d0, 0.00d0/) ; kpath_name(1)= 'T' - !ke(1,:)=(/0.0d0, 0.20d0/) - !kp(2,:)=(/0.0d0, 0.20d0/) ; kpath_name(2)= 'P' - !ke(2,:)=(/0.5d0, 0.50d0/) ; kpath_name(3)= 'M' - - !> 0.19 0.195 0.2 0.21 0.23 0.25 0.26 0.29 - ! 0.105 0.1076 0.11 0.116 0.127 0.138 0.1435 0.16 - !kp(1,:)=(/ 0.1380d0,-0.18d0/) ; kpath_name(1)= 'T' - !ke(1,:)=(/ 0.1380d0, 0.00d0/) - !kp(2,:)=(/ 0.1380d0, 0.00d0/) ; kpath_name(2)= 'G' - !ke(2,:)=(/ 0.1380d0, 0.18d0/) - - kp(1,:)=(/ 0.0000d0, 0.00d0/) ; kpath_name(1)= 'X' - ke(1,:)=(/ 0.3000d0, 0.00d0/) - - nlines=1 - NN= Nk - knv2=NN*nlines - allocate( kpoint(knv2, 2)) - allocate( k_len (knv2)) - kpoint= 0d0 - - t1=0d0 - k_len=0d0 - kpath_stop= 0d0 - do j=1, nlines - do i=1, NN - kstart= kp(j,:) - kend = ke(j,:) - k1= kstart(1)*Ka2+ kstart(2)*Kb2 - k2= kend(1)*Ka2+ kend(2)*Kb2 - kpoint(i+(j-1)*NN,:)= kstart+ (kend-kstart)*(i-1)/dble(NN-1) - - temp= dsqrt((k2(1)- k1(1))**2 & - +(k2(2)- k1(2))**2)/dble(NN-1) - - if (i.gt.1) then - t1=t1+temp - endif - k_len(i+(j-1)*NN)= t1 - enddo - kpath_stop(j+1)= t1 - enddo - - !> lines across two weyl points - !k1= -0.5d0*ka2+0.5d0*kb2 - !t1= 0d0 - !do i=1, NN - ! s= (i-1d0)/(NN-1d0)-0.5d0 - ! k2= k1 - ! kpoint(i, 1)= s - ! kpoint(i, 2)= -2.818d0*s**3-0.2955*s - ! k1= kpoint(i, 1)*ka2+ kpoint(i, 2)*kb2 - ! temp= dsqrt((k2(1)- k1(1))**2 & - ! +(k2(2)- k1(2))**2) - - ! if (i.gt.1) then - ! t1=t1+temp - ! endif - ! k_len(i)= t1 - !enddo - - !> k lines around one weyl points - !> for WTe2, W1 (0.121480, 0.0447), W2= 0.121831, 0.0388) - !> W0 (0.1216, 0.042) - - !> for MoTe2 W1 (0.10296, 0.05388) E=0.06eV - !x0= 0.10296d0 - !y0= 0.05388d0 - - !> W2 (0.10516, 0.01378) E=0.005eV - !x0= 0.10516d0 - !y0= 0.01378d0 - - !rx0= 0.001d0 - !ry0= 0.005d0 - !s=0 - !k2(1)= x0+ rx0*cos(s+pi) - !k2(2)= y0 + ry0*sin(s+pi) - !k1= k2(1)*ka2+ k2(2)*kb2 - !t1=0d0 - !do i=1, NN - ! k2= k1 - ! s= (i-1d0)/(NN-1d0)*2d0*pi+pi - ! kpoint(i, 1)= x0 + rx0*cos(s) - ! kpoint(i, 2)= y0 + ry0*sin(s) - ! k1= kpoint(i, 1)*ka2+ kpoint(i, 2)*kb2 - ! temp= dsqrt((k2(1)- k1(1))**2 & - ! + (k2(2)- k1(2))**2) - - ! if (i.gt.1) then - ! t1=t1+temp - ! endif - ! k_len(i)= t1 - !enddo - - - Nwann= Num_wann/2 - if (SOC ==0 ) Nwann= Num_wann - allocate( omega(omeganum)) - allocate( dos_l(Nwann, knv2, omeganum )) - allocate( dos_r(Nwann, knv2, omeganum )) - allocate( dos_l_mpi(Nwann, knv2, omeganum )) - allocate( dos_r_mpi(Nwann, knv2, omeganum )) - omega=0d0 - dos_l=0d0 - dos_r=0d0 - dos_l_mpi=0d0 - dos_r_mpi=0d0 - - eta=(omegamax- omegamin)/dble(omeganum)*3.0d0 - - do i= 1, omeganum - omega(i)=omegamin+(i-1)*(omegamax-omegamin)/dble(omeganum) - enddo - - allocate(GLL(Ndim, Ndim)) - allocate(GRR(Ndim, Ndim)) - allocate(H00(Ndim, Ndim)) - allocate(H01(Ndim, Ndim)) - allocate(ones(Ndim, Ndim)) - GLL= 0d0 - GRR= 0d0 - H00= 0d0 - H01= 0d0 - ones= 0d0 - - do i=1,Ndim - ones(i,i)=1.0d0 - enddo - - do ikp= 1+cpuid, knv2, num_cpu - if (cpuid==0) write(*, *) ikp, 'ik', knv2 - if (cpuid==0) write(stdout, *) ikp, 'ik', knv2 - k= kpoint(ikp,:) - - !> get the hopping matrix between two principle layers - call ham_qlayer2qlayer(k,H00,H01) - - do j = 1, omeganum - w=omega(j) - - !> calculate surface green function - ! 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,H00,H01,ones) - !call surfgreen_1984(w,GLL,GRR,H00,H01,ones) - - ! calculate spectral function - do i= 1, Nwann - dos_l(i, ikp, j)=dos_l(i, ikp, j)- aimag(GLL(i,i)) & - - aimag(GLL(i+Nwann,i+Nwann)) - enddo ! i - do i=Ndim-Num_wann+1, Ndim- Nwann - dos_r(i, ikp, j)=dos_r(i, ikp, j)- aimag(GRR(i,i)) & - - aimag(GRR(i+Nwann,i+Nwann)) - enddo ! i - enddo ! j - enddo ! ikp - - !> we don't have to do allreduce operation - call mpi_reduce(dos_l, dos_l_mpi, size(dos_l), mpi_double_precision,& - mpi_sum, 0, mpi_comm_world, ierr) - call mpi_reduce(dos_r, dos_r_mpi, size(dos_r), mpi_double_precision,& - mpi_sum, 0, mpi_comm_world, ierr) - !dos_l=log(dos_l_mpi) - !dos_r=log(dos_r_mpi) - dos_l= dos_l_mpi/maxval(dos_l_mpi)*255d0 - dos_r= dos_r_mpi/maxval(dos_r_mpi)*255d0 - - if (cpuid.eq.0)then - open (unit=16, file='dos.dat_l') - open (unit=13, file='dos.dat_r') - do ikp=1, knv2 - do j=1, omeganum - write(16, '(2f16.8, 30000i5)')k_len(ikp), omega(j), int8(dos_l(:, ikp, j)) - write(13, '(2f16.8, 30000i5)')k_len(ikp), omega(j), int8(dos_r(:, ikp, j)) - enddo - write(16, *) ' ' - write(13, *) ' ' - enddo - close(16) - close(13) - write(stdout,*)'ndim',ndim - write(stdout,*) 'knv2,omeganum,eta',knv2, omeganum, eta - write(stdout,*)'calculate density of state successfully' - endif - - emin= minval(omega) - emax= maxval(omega) - !> write script for gnuplot - if (cpuid==0) then - open(unit=104, file='surfdos_l.gnu') - write(104, '(a)')'#set terminal postscript enhanced color' - write(104, '(a)')"#set output 'surfdos_l.eps'" - write(104, '(3a)')'set terminal pngcairo truecolor enhanced', & - ' font ", 36" size 1920, 1680' - write(104, '(a)')"set output 'surfdos_l.png'" - write(104,'(2a)') 'set palette defined (-20 "#194eff", ', & - '0 "white", 10 "red" )' - write(104, '(a)')'#set palette rgbformulae 33,13,10' - write(104, '(a)')'set style data linespoints' - write(104, '(a)')'#set size ratio -1' - write(104, '(a)')'unset ztics' - write(104, '(a)')'unset key' - write(104, '(a)')'set pointsize 0.8' - write(104, '(a)')'set pm3d' - write(104, '(a)')'#set view equal xyz' - write(104, '(a)')'set view map' - write(104, '(a)')'set border lw 3' - write(104, '(a)')'set cbtics font ",48"' - write(104, '(a)')'set xtics font ",48"' - write(104, '(a)')'set ytics font ",48"' - write(104, '(a)')'set ylabel font ",48"' - write(104, '(a)')'set ylabel "Energy (eV)"' - write(104, '(a)')'#set xtics offset 0, -1' - write(104, '(a)')'#set ylabel offset -6, 0 ' - write(104, '(a, f8.5, a)')'set xrange [0: ', maxval(k_len), ']' - write(104, '(a, f8.5, a, f8.5, a)')'set yrange [', emin, ':', emax, ']' - write(104, 202, advance="no") (kpath_name(i), kpath_stop(i), i=1, nlines) - write(104, 203)kpath_name(nlines+1), kpath_stop(nlines+1) - - do i=1, nlines-1 - write(104, 204)kpath_stop(i+1), emin, kpath_stop(i+1), emax - enddo - write(104, '(a)')'set pm3d interpolate 2,2' - write(104, '(2a)')"splot 'dos.dat_l' u 1:2:3 w pm3d" - - endif - - !> write script for gnuplot - if (cpuid==0) then - open(unit=105, file='surfdos_r.gnu') - write(105, '(a)')'#set terminal postscript enhanced color' - write(105, '(a)')"#set output 'surfdos_r.eps'" - write(105, '(3a)')'set terminal pngcairo truecolor enhanced', & - ' font ", 36" size 1920, 1680' - write(105, '(a)')"set output 'surfdos_r.png'" - write(105,'(2a)') 'set palette defined (-20 "#194eff", ', & - '0 "white", 10 "red" )' - write(105, '(a)')'#set palette rgbformulae 33,13,10' - write(105, '(a)')'set style data linespoints' - write(105, '(a)')'unset ztics' - write(105, '(a)')'unset key' - write(105, '(a)')'set pointsize 0.8' - write(105, '(a)')'set pm3d' - write(105, '(a)')'set border lw 3' - write(105, '(a)')'#set size ratio -1' - write(105, '(a)')'#set view equal xyz' - write(105, '(a)')'set view map' - write(105, '(a)')'set cbtics font ",48"' - write(105, '(a)')'set xtics font ",48"' - write(105, '(a)')'set ytics font ",48"' - write(105, '(a)')'set ylabel font ",48"' - write(105, '(a)')'set ylabel "Energy (eV)"' - write(105, '(a)')'#set xtics offset 0, -1' - write(105, '(a)')'#set ylabel offset -6, 0 ' - write(105, '(a, f8.5, a)')'set xrange [0: ', maxval(k_len), ']' - write(105, '(a, f8.5, a, f8.5, a)')'set yrange [', emin, ':', emax, ']' - write(105, 202, advance="no") (kpath_name(i), kpath_stop(i), i=1, nlines) - write(105, 203)kpath_name(nlines+1), kpath_stop(nlines+1) - - do i=1, nlines-1 - write(105, 204)kpath_stop(i+1), emin, kpath_stop(i+1), emax - enddo - write(105, '(a)')'set pm3d interpolate 2,2' - write(105, '(2a)')"splot 'dos.dat_r' u 1:2:3 w pm3d" - - endif - - - 202 format('set xtics (',:20('"',A1,'" ',F8.5,',')) - 203 format(A1,'" ',F8.5,')') - 204 format('set arrow from ',F8.5,',',F10.5,' to ',F8.5,',',F10.5, ' nohead front lw 3') - - - return - end subroutine diff --git a/2d/symmetry.f90 b/2d/symmetry.f90 deleted file mode 100644 index 082e9d5d..00000000 --- a/2d/symmetry.f90 +++ /dev/null @@ -1,198 +0,0 @@ -!> set symmetry - subroutine symmetry - use para - implicit none - - integer :: nwan - integer :: ia, i, n - - !> get the atom afterr the mirror_x operation - integer, allocatable :: iatom_mirror_x(:) - - !> set up symmetry operators - !> here we assume that the wannier functions have the symmetry - !> as atomic orbitals - - nwan= Num_wann/2 - - allocate(iatom_mirror_x(Num_atoms)) - allocate(inversion(Num_wann, Num_wann)) - allocate(mirror_z(Num_wann, Num_wann)) - allocate(mirror_x(Num_wann, Num_wann)) - inversion= 0d0 - mirror_z = 0d0 - mirror_x = 0d0 - iatom_mirror_x= 0 - - !> inversion symmetry - !> s-> s; p-> -p; d-> -d - n= 0 - do ia=1, Num_atoms - do i=1, nprojs(ia) - n= n+ 1 - select case (proj_name(i, ia)) - case ('s', 'S') - inversion(n, n)= 1 - inversion(n+ nwan, n+ nwan)= 1 - case ('px', 'Px', 'PX') - inversion(n, n)=-1 - inversion(n+ nwan, n+ nwan)=-1 - case ('py', 'Py', 'PY') - inversion(n, n)=-1 - inversion(n+ nwan, n+ nwan)=-1 - case ('pz', 'Pz', 'PZ') - inversion(n, n)=-1 - inversion(n+ nwan, n+ nwan)=-1 - case ('dxy', 'Dxy', 'DXY') - inversion(n, n)= 1 - inversion(n+ nwan, n+ nwan)= 1 - case ('dyz', 'Dyz', 'DYZ') - inversion(n, n)= 1 - inversion(n+ nwan, n+ nwan)= 1 - case ('dxz', 'Dxz', 'DXZ', 'dzx', 'Dzx', 'DZX') - inversion(n, n)= 1 - inversion(n+ nwan, n+ nwan)= 1 - case ('dx2-y2', 'Dx2-y2', 'DX2-Y2', 'dx2', 'DX2') - inversion(n, n)= 1 - inversion(n+ nwan, n+ nwan)= 1 - case ('dz2', 'Dz2', 'DZ2') - inversion(n, n)= 1 - inversion(n+ nwan, n+ nwan)= 1 - case default - write(*, *) "ERROR: only support s px py pz dxy dyz dxz dx2-y2 dz2 orbitals" - stop - end select - enddo ! i - enddo ! ia - - !> mirror_x symmetry - !> s-> s; px->-px, py->py, pz-> pz - !> dxy-> -dxy, dyz-> dyz, dxz-> -dxz, dx2-> dx2 dz2->dz2 - !> up-> -i dn dn-> -i up - !> here we throw away the phase -i, it is just a constant, leading M*M= -1 - - n= 0 - mirror_x= 0d0 - do ia=1, Num_atoms - do i=1, nprojs(ia) - n= n+ 1 - select case (proj_name(i, ia)) - case ('s', 'S') - mirror_x(n, n+ nwan)= 1d0 - mirror_x(n+ nwan, n)= 1d0 - case ('px', 'Px', 'PX') - mirror_x(n, n+ nwan)=-1d0 - mirror_x(n+ nwan, n)=-1d0 - case ('py', 'Py', 'PY') - mirror_x(n, n+ nwan)= 1d0 - mirror_x(n+ nwan, n)= 1d0 - case ('pz', 'Pz', 'PZ') - mirror_x(n, n+ nwan)= 1d0 - mirror_x(n+ nwan, n)= 1d0 - case ('dxy', 'Dxy', 'DXY') - mirror_x(n, n+ nwan)=-1d0 - mirror_x(n+ nwan, n)=-1d0 - case ('dyz', 'Dyz', 'DYZ') - mirror_x(n, n+ nwan)= 1d0 - mirror_x(n+ nwan, n)= 1d0 - case ('dxz', 'Dxz', 'DXZ', 'dzx', 'Dzx', 'DZX') - mirror_x(n, n+ nwan)=-1d0 - mirror_x(n+ nwan, n)=-1d0 - case ('dx2-', 'dx2-y2', 'Dx2-y2', 'DX2-Y2', 'dx2', 'DX2') - mirror_x(n, n+ nwan)= 1d0 - mirror_x(n+ nwan, n)= 1d0 - case ('dz2', 'Dz2', 'DZ2') - mirror_x(n, n+ nwan)= 1d0 - mirror_x(n+ nwan, n)= 1d0 - case default - write(*, *) "ERROR: only support s px py pz dxy dyz dxz dx2-y2 dz2 orbitals" - stop - end select - enddo ! i - enddo ! ia - - do i=1, Num_wann - !write(*, '(1000i2)')int(real(mirror_x(:, i))) - enddo - - - !> mirror_z symmetry - !> s-> s; px->px, py->py, pz-> -pz - !> dxy-> dxy, dyz-> -dyz, dxz-> -dxz, dx2-> dx2 dz2->dz2 - !> up-> up dn-> -dn - n= 0 - do ia=1, Num_atoms - do i=1, nprojs(ia) - n= n+ 1 - select case (proj_name(i, ia)) - case ('s', 'S') - mirror_z(n, n)= 1 - mirror_z(n+ nwan, n+ nwan)=-1 - case ('px', 'Px', 'PX') - mirror_z(n, n)= 1 - mirror_z(n+ nwan, n+ nwan)=-1 - case ('py', 'Py', 'PY') - mirror_z(n, n)= 1 - mirror_z(n+ nwan, n+ nwan)=-1 - case ('pz', 'Pz', 'PZ') - mirror_z(n, n)=-1 - mirror_z(n+ nwan, n+ nwan)= 1 - case ('dxy', 'Dxy', 'DXY') - mirror_z(n, n)= 1 - mirror_z(n+ nwan, n+ nwan)=-1 - case ('dyz', 'Dyz', 'DYZ') - - mirror_z(n, n)=-1 - mirror_z(n+ nwan, n+ nwan)= 1 - case ('dxz', 'Dxz', 'DXZ', 'dzx', 'Dzx', 'DZX') - mirror_z(n, n)=-1 - mirror_z(n+ nwan, n+ nwan)= 1 - case ('dx2-', 'dx2-y2', 'Dx2-y2', 'DX2-Y2', 'dx2', 'DX2') - mirror_z(n, n)= 1 - mirror_z(n+ nwan, n+ nwan)=-1 - case ('dz2', 'Dz2', 'DZ2') - mirror_z(n, n)= 1 - mirror_z(n+ nwan, n+ nwan)=-1 - case default - write(*, *) "ERROR: only support s px py pz dxy dyz dxz dx2-y2 dz2 orbitals" - stop - end select - enddo ! i - enddo ! ia - - - !> set up symmetry operators - allocate(mirror_x_op(3,3)) - allocate(mirror_y_op(3,3)) - allocate(mirror_z_op(3,3)) - - !> for glide symmetry, (1:3, 1:3) shows the mirror operation, (1:3, 4) - !> gives the shift - allocate(glide_y_op(3,4)) - - mirror_x_op= 0d0 - mirror_y_op= 0d0 - mirror_z_op= 0d0 - glide_y_op= 0d0 - - mirror_x_op(1, 1)=-1d0 - mirror_x_op(2, 2)= 1d0 - mirror_x_op(3, 3)= 1d0 - - mirror_y_op(1, 1)= 1d0 - mirror_y_op(2, 2)=-1d0 - mirror_y_op(3, 3)= 1d0 - - mirror_z_op(1, 1)= 1d0 - mirror_z_op(2, 2)= 1d0 - mirror_z_op(3, 3)=-1d0 - - glide_y_op(1, 1) = 1d0 - glide_y_op(2, 2) =-1d0 - glide_y_op(3, 3) = 1d0 - glide_y_op(:, 4) = (/0.5d0, 0.0d0, 0.5d0/) - - - - return - end subroutine symmetry diff --git a/2d/wanniercenter.f90 b/2d/wanniercenter.f90 deleted file mode 100644 index 03c8d7bd..00000000 --- a/2d/wanniercenter.f90 +++ /dev/null @@ -1,412 +0,0 @@ - !> this suboutine is used for wannier center calculation for 3D system - !> only for one plane - - subroutine wannier_center2D_plane - use para - use wmpi - implicit none - - integer :: i - integer :: j - integer :: l - integer :: m - integer :: ia - integer :: nfill - integer :: imax - - integer :: ik1 - integer :: ik2 - - integer :: ierr - - !> k points in kx-ky plane - real(dp), allocatable :: kpoints(:, :, :) - - !> hamiltonian for each k point - !> and also the eigenvector of hamiltonian after eigensystem_c - complex(dp), allocatable :: Hamk(:, :) - complex(dp), allocatable :: Hamk_dag(:, :) - - !> eigenvector for each kx - complex(dp), allocatable :: Eigenvector(:, :, :) - - !> Mmnkb= - !> |u_n(k)> is the periodic part of wave function - complex(dp), allocatable :: Mmnkb(:, :) - complex(dp), allocatable :: Mmnkb_com(:, :) - complex(dp), allocatable :: Mmnkb_full(:, :) - - !> - complex(dp), allocatable :: Lambda_eig(:) - complex(dp), allocatable :: Lambda(:, :) - complex(dp), allocatable :: Lambda0(:, :) - - !> three matrix for SVD - !> M= U.Sigma.V^\dag - !> VT= V^\dag - complex(dp), allocatable :: U(:, :) - real (dp), allocatable :: Sigma(:, :) - complex(dp), allocatable :: VT(:, :) - - !> wannier centers for each ky, bands - real(dp), allocatable :: WannierCenterKy(:, :) - real(dp), allocatable :: WannierCenterKy_mpi(:, :) - - !> larges gap of each two wannier centers for a given k point - !> dim= Nky - real(dp), allocatable :: largestgap(:) - real(dp), allocatable :: largestgap_mpi(:) - - !> eigenvalue - real(dp), allocatable :: eigenvalue(:) - - !> for each orbital, it correspond to an atom - !> dim= Num_wann - integer, allocatable :: AtomIndex_orbital(:) - - real(dp) :: Umatrix_t(2, 2) - - !> b.r - real(dp) :: br - - !> exp(-i*b.R) - complex(dp) :: ratio - - real(dp) :: k(2) - real(dp) :: b(2) - - real(dp) :: maxgap - real(dp) :: maxgap0 - - !> Z2 calculation for time reversal invariant system - integer :: Z2 - integer :: Delta - - real(dp) :: g - real(dp) :: phi1 - real(dp) :: phi2 - real(dp) :: phi3 - real(dp) :: zm - real(dp) :: zm1 - real(dp) :: xnm1 - real(dp) :: Deltam - real(dp), allocatable :: xnm(:) - real(dp) :: k0(2), k1(2), k2(2) - - - - nfill= Numoccupied - - allocate(kpoints(2, Nk1, Nk2)) - kpoints= 0d0 - - allocate(Lambda_eig(nfill)) - allocate(Lambda(nfill, nfill)) - allocate(Lambda0(nfill, nfill)) - allocate(Mmnkb(nfill, nfill)) - allocate(Mmnkb_com(nfill, nfill)) - allocate(Mmnkb_full(Num_wann, Num_wann)) - allocate(hamk(Num_wann, Num_wann)) - allocate(hamk_dag(Num_wann, Num_wann)) - allocate(Eigenvector(Num_wann, Num_wann, Nk1)) - allocate(eigenvalue(Num_wann)) - allocate(U(nfill, nfill)) - allocate(Sigma(nfill, nfill)) - allocate(VT(nfill, nfill)) - allocate(WannierCenterKy(nfill, Nk2)) - allocate(WannierCenterKy_mpi(nfill, Nk2)) - allocate(AtomIndex_orbital(Num_wann)) - allocate(xnm(nfill)) - allocate(largestgap(Nk2)) - allocate(largestgap_mpi(Nk2)) - largestgap= 0d0 - largestgap_mpi= 0d0 - WannierCenterKy= 0d0 - WannierCenterKy_mpi= 0d0 - hamk=0d0 - eigenvalue=0d0 - Eigenvector=0d0 - Mmnkb_full=0d0 - Mmnkb=0d0 - Mmnkb_com=0d0 - Lambda =0d0 - Lambda0=0d0 - U= 0d0 - Sigma= 0d0 - VT= 0d0 - - !> set k plane - !> the first dimension should be in one primitive cell, [0, 2*pi] - !> the first dimension is the integration direction - !> the WCCs are calculated along the second k line - k0= K2D_start ! - k1= K2D_vec1 ! - k2= K2D_vec2 ! - - do ik2=1, Nk2 - do ik1=1, Nk1 - kpoints(:, ik1, ik2)= k0+ k1*(ik1-1d0)/dble(Nk1)+ k2*(ik2-1d0)/dble(Nk2-1) - enddo - enddo - b= k1/dble(Nk1) - b= b(1)*kua+b(2)*kub - - !> set up atom index for each orbitals in the basis - if (soc>0) then !> with spin orbital coupling - l= 0 - do ia=1, Num_atoms !> spin up - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia - enddo ! l - enddo ! ia - do ia=1, Num_atoms !> spin down - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia - enddo ! l - enddo ! ia - else !> without spin orbital coupling - l= 0 - do ia=1, Num_atoms !> spin down - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia - enddo ! l - enddo ! ia - - endif - - Umatrix_t= transpose(Umatrix) - call inv_r(2, Umatrix_t) - - !>> Get wannier center for ky=0 plane - !> for each ky, we can get wanniercenter - do ik2=1+ cpuid, Nk2, num_cpu - if (cpuid.eq.0) write(stdout, *)' Wilson loop ', 'ik, Nk', ik2, nk2 - Lambda0=0d0 - do i=1, nfill - Lambda0(i, i)= 1d0 - enddo - - !> for each k1, we get the eigenvectors - do ik1=1, Nk1 - !if (cpuid==0) print *, ik1, ik2, Nk1, Nk2 - k= kpoints(:, ik1, ik2) - - call ham_bulk_old(k,hamk) - !call ham_bulk (k,hamk) - - !> diagonal hamk - call eigensystem_c('V', 'U', Num_wann, hamk, eigenvalue) - - Eigenvector(:, :, ik1)= hamk - enddo - - !> sum over k1 to get wanniercenters - do ik1=1, Nk1 - !> - Mmnkb= 0d0 - hamk_dag= Eigenvector(:, :, ik1) - if (ik1==Nk1) then - hamk= Eigenvector(:, :, 1) - else - hamk= Eigenvector(:, :, ik1+ 1) - endif - do m=1, Num_wann - ia= AtomIndex_orbital(m) - br= b(1)*Atom_position(1, ia)+ & - b(2)*Atom_position(2, ia) - ratio= cos(br)- zi* sin(br) - !ratio= 1d0 - - do j=1, nfill - do i=1, nfill - Mmnkb(i, j)= Mmnkb(i, j)+ & - conjg(hamk_dag(m, i))* hamk(m, j)* ratio - enddo ! i - enddo ! j - enddo ! m - - !> perform Singluar Value Decomposed of Mmnkb - call zgesvd_pack(nfill, Mmnkb, U, Sigma, VT) - - !> after the calling of zgesvd_pack, Mmnkb becomes a temporal matrix - U= conjg(transpose(U)) - VT= conjg(transpose(VT)) - call mat_mul(nfill, VT, U, Mmnkb) - - call mat_mul(nfill, Mmnkb, Lambda0, Lambda) - Lambda0 = Lambda - enddo !< ik1 - - !> diagonalize Lambda to get the eigenvalue - call zgeev_pack(nfill, Lambda, Lambda_eig) - do i=1, nfill - WannierCenterKy(i, ik2)= aimag(log(Lambda_eig(i)))/2d0/pi - WannierCenterKy(i, ik2)= mod(WannierCenterKy(i, ik2)+10d0, 1d0) - enddo - - call sortheap(nfill, WannierCenterKy(:, ik2)) - - maxgap0= -99999d0 - imax= nfill - do i=1, nfill - if (i/=nfill) then - maxgap= WannierCenterKy(i+1, ik2)- WannierCenterKy(i, ik2) - else - maxgap=1d0+ WannierCenterKy(1, ik2)- WannierCenterKy(nfill, ik2) - endif - - if (maxgap>maxgap0) then - maxgap0= maxgap - imax= i - endif - - enddo - - if (imax==nfill) then - largestgap(ik2)= (WannierCenterKy(1, ik2)+ & - WannierCenterKy(nfill, ik2) +1d0)/2d0 - largestgap(ik2)= mod(largestgap(ik2), 1d0) - else - largestgap(ik2)= (WannierCenterKy(imax+1, ik2)+ & - WannierCenterKy(imax, ik2))/2d0 - endif - - - enddo !< ik2 - - WannierCenterKy_mpi= 0d0 - largestgap_mpi= 0d0 - call mpi_allreduce(WannierCenterKy, WannierCenterKy_mpi, & - size(WannierCenterKy), mpi_dp, mpi_sum, mpi_cmw, ierr) - call mpi_allreduce(largestgap, largestgap_mpi, & - size(largestgap), mpi_dp, mpi_sum, mpi_cmw, ierr) - - if (cpuid==0) then - open(unit=110, file='wcc.dat') - - do ik2=1, Nk2 - write(110, '(10000f16.8)') dble(ik2-1)/dble(Nk2-1)/2d0, & - largestgap_mpi(ik2), dmod(sum(WannierCenterKy_mpi(:, ik2)), 1d0), & - WannierCenterKy_mpi(:, ik2) - enddo - close(110) - endif - - - - !> Z2 calculation Alexey Soluyanov arXiv:1102.5600 - - Delta= 0 - !> for each iky, we get a Deltam - do ik2=1, nk2-1 - - !> largestgap position - zm= largestgap_mpi(ik2) - !if (ik2==nk2) then - ! zm1= largestgap_mpi(1) - ! xnm= WannierCenterKy_mpi(1:nfill, 1) - !else - zm1= largestgap_mpi(ik2+1) - xnm= WannierCenterKy_mpi(1:nfill, ik2+1) - !endif - Deltam= 1 - do i=1, nfill - xnm1= xnm(i) - phi1= 2d0*pi*zm - phi2= 2d0*pi*zm1 - phi3= 2d0*pi*xnm1 - - g= sin(phi2-phi1)+ sin(phi3-phi2)+ sin(phi1-phi3) - Deltam= Deltam* sign(1d0, g) - enddo !i - if (Deltam<0) then - Delta= Delta+ 1 - endif - enddo !ik2 - - Z2= mod(Delta, 2) - - if (cpuid==0) write(stdout, *)'Z2 for the plane you choose: ', Z2 - - - !> generate gnu script for wannier charge center plots - if (cpuid==0) then - open(unit=301, file='wcc.gnu') - write(301, '(a)')"set encoding iso_8859_1" - write(301, '(a)')'set terminal postscript enhanced color font ",30"' - write(301, '(a)')"set output 'wcc.eps'" - write(301, '(a)')'unset key ' - write(301, '(a)')'set border lw 3 ' - write(301, '(a)')'set xtics offset 0, 0.2' - write(301, '(a)')'set xtics format "%4.1f" nomirror out ' - write(301, '(a)')'set xlabel "k" ' - write(301, '(a)')'set xlabel offset 0, 0.7 ' - write(301, '(a)')'set ytics 0.5 ' - write(301, '(a)')'set ytics format "%4.1f" nomirror out' - write(301, '(a)')'set ylabel "WCC"' - write(301, '(a)')'set ylabel offset 2, 0.0 ' - write(301, '(a)')'set xrange [0: 0.5]' - write(301, '(a)')'set yrange [0:1]' - write(301, '(a)')"plot 'wcc.dat' u 1:2 w l lw 2 lc 'blue', \" - write(301, '(a, i5, a)')" for [i=4: ", nfill+3, "] 'wcc.dat' u 1:i w p pt 7 ps 1.1 lc 'red'" - close(301) - endif - - return - end subroutine wannier_center2D_plane - - - subroutine sortheap(n, arr) - implicit none - integer, intent(in) :: n - real(8), intent(inout) :: arr(n) - - !* local variables - integer :: i - - do i=n/2, 1, -1 - call sift_down(i, n) - enddo - - do i=n, 2, -1 - call swap(arr(1), arr(i)) - call sift_down(1, i-1) - enddo - contains - subroutine sift_down(l, r) - integer, intent(in) :: l, r - integer :: j, jold - real(8) :: a - a= arr(l) - jold= l - j= l+ l - - do while (j<=r) - if (j= arr(j))exit - arr(jold)= arr(j) - jold= j - j= j+ j - enddo - arr(jold)= a - return - end subroutine sift_down - end subroutine sortheap - - !>> swap two real numbers - subroutine swap(a, b) - real(8), intent(inout) :: a - real(8), intent(inout) :: b - real(8) :: c - c=a - a=b - b=c - return - end subroutine swap - - diff --git a/2d/wanniercenter2.f90 b/2d/wanniercenter2.f90 deleted file mode 100644 index 56b5940a..00000000 --- a/2d/wanniercenter2.f90 +++ /dev/null @@ -1,1099 +0,0 @@ -! this suboutine is used for wannier center calculation for slab system -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - - subroutine wannier_center2D - use para - use wmpi - implicit none - - integer :: Nkx - integer :: Nky - integer :: knv2 - - integer :: i - integer :: j - integer :: l - integer :: nfill - - integer :: ikx - integer :: iky - - integer :: ierr - - !> k points in kx-ky plane - real(dp), allocatable :: kpoints(:, :, :) - - !> hamiltonian for each k point - !> and also the eigenvector of hamiltonian after eigensystem_c - complex(dp), allocatable :: Hamk(:, :) - complex(dp), allocatable :: Hamk_dag(:, :) - - !> eigenvector for each kx - complex(dp), allocatable :: Eigenvector(:, :, :) - - !> Mmnkb= - !> |u_n(k)> is the periodic part of wave function - complex(dp), allocatable :: Mmnkb(:, :) - complex(dp), allocatable :: Mmnkb_com(:, :) - complex(dp), allocatable :: Mmnkb_full(:, :) - - !> - complex(dp), allocatable :: Lambda_eig(:) - complex(dp), allocatable :: Lambda(:, :) - complex(dp), allocatable :: Lambda0(:, :) - - !> three matrix for SVD - !> M= U.Sigma.V^\dag - !> VT= V^\dag - complex(dp), allocatable :: U(:, :) - real (dp), allocatable :: Sigma(:, :) - complex(dp), allocatable :: VT(:, :) - - !> wannier centers for each ky, bands - real(dp), allocatable :: WannierCenterKy(:, :) - real(dp), allocatable :: WannierCenterKy_mpi(:, :) - - !> eigenvalue - real(dp), allocatable :: eigenvalue(:) - - real(dp) :: kx - real(dp) :: ky - real(dp) :: k(2) - real(dp) :: k11(3), k12(3) - real(dp) :: k21(3), k22(3) - - Nkx= Nk - Nky= 40 - knv2= Nkx*Nky - - nfill= Numoccupied*Nslab - - allocate(kpoints(2, Nkx, Nky)) - kpoints= 0d0 - - allocate(Lambda_eig(nfill)) - allocate(Lambda(nfill, nfill)) - allocate(Lambda0(nfill, nfill)) - allocate(Mmnkb(nfill, nfill)) - allocate(Mmnkb_com(nfill, nfill)) - allocate(Mmnkb_full(Num_wann*Nslab, Num_wann*Nslab)) - allocate(hamk(Num_wann*Nslab, Num_wann*Nslab)) - allocate(hamk_dag(Num_wann*Nslab, Num_wann*Nslab)) - allocate(Eigenvector(Num_wann*Nslab, Num_wann*Nslab, Nkx)) - allocate(eigenvalue(Num_wann*Nslab)) - allocate(U(nfill, nfill)) - allocate(Sigma(nfill, nfill)) - allocate(VT(nfill, nfill)) - allocate(WannierCenterKy(nfill, Nky)) - allocate(WannierCenterKy_mpi(nfill, Nky)) - WannierCenterKy= 0d0 - WannierCenterKy_mpi= 0d0 - hamk=0d0 - eigenvalue=0d0 - Eigenvector=0d0 - Mmnkb_full=0d0 - Mmnkb=0d0 - Mmnkb_com=0d0 - Lambda =0d0 - Lambda0=0d0 - U= 0d0 - Sigma= 0d0 - VT= 0d0 - - !> set k plane - !> the second dimension should be in one primitive plane - k11=(/ 0.0d0, 0.0d0, 0.0d0/) ! - k12=(/ 1.0d0, 0.0d0, 0.0d0/) ! X - k21=(/ 0.0d0, 0.0d0, 0.0d0/) ! - k22=(/ 0.0d0, 1.0d0, 0.0d0/) ! Z - - - do iky=1, Nky - do ikx=1, Nkx - kxy(:, ikx, iky)= k11+(k12-k11)*(i-1)/dble(nkx)+ k21+ (k22-k21)*(j-1)/dble(nky) - enddo - enddo - - !> for each ky, we can get wanniercenter - do iky=1+ cpuid, nky, num_cpu - Lambda0=0d0 - do i=1, nfill - Lambda0(i, i)= 1d0 - enddo - - if (cpuid==0) print *, iky, nky - !> for each kx, we get the eigenvectors - do ikx=1, nkx - k(1)= kpoints(1, ikx, iky) - k(2)= kpoints(2, ikx, iky) - - call ham_slab(k,hamk) - - !> diagonal hamk - call eigensystem_c('V', 'U', Num_wann*Nslab, hamk, eigenvalue) - - Eigenvector(:, :, ikx)= hamk - enddo - - !> sum over kx to get wanniercenters - do ikx=1, nkx - hamk= Eigenvector(:, :, ikx) - hamk_dag= conjg(transpose(hamk)) - if (ikx==nkx) then - hamk= Eigenvector(:, :, 1) - else - hamk= Eigenvector(:, :, ikx+ 1) - endif - - !> - call mat_mul(Num_wann*Nslab, hamk_dag, hamk, Mmnkb_full) - Mmnkb= Mmnkb_full(1:nfill, 1:nfill) - !Mmnkb_com= 0d0 - !hamk_dag= Eigenvector(:, :, ikx) - !hamk= Eigenvector(:, :, ikx+1) - !do i=1, nfill - ! do j=1, nfill - ! do l= 1, Num_wann*Nslab - ! Mmnkb_com(i, j)= Mmnkb_com(i, j)+ conjg(hamk_dag(l, i))* hamk(l, j) - ! enddo - ! enddo - !enddo - - !print *, maxval(real(Mmnkb-Mmnkb_com)) - !stop - - - !> perform Singluar Value Decomposed of Mmnkb - call zgesvd_pack(nfill, Mmnkb, U, Sigma, VT) - - !> after the calling of zgesvd_pack, Mmnkb becomes a temporal matrix - U= conjg(transpose(U)) - VT= conjg(transpose(VT)) - call mat_mul(nfill, VT, U, Mmnkb) - - !> check hermicity - !do i=1, nfill - ! do j=i, nfill - ! if (abs(Mmnkb(i, j)-conjg(Mmnkb(j, i)))>0.0001d0)then - ! print *, 'Mmnkb is not Hermitian' - ! print*, i, j, Mmnkb(i, j), Mmnkb(j, i) - - ! endif - ! enddo - !enddo - - !stop - - - call mat_mul(nfill, Mmnkb, Lambda0, Lambda) - Lambda0 = Lambda - enddo !< ikx - - !> diagonalize Lambda to get the eigenvalue - call zgeev_pack(nfill, Lambda, Lambda_eig) - do i=1, nfill - WannierCenterKy(i, iky)= aimag(log(Lambda_eig(i)))/2d0/pi - enddo - - enddo !< iky - - call mpi_allreduce(WannierCenterKy, WannierCenterKy_mpi, & - size(WannierCenterKy), mpi_dp, mpi_sum, mpi_cmw, ierr) - - if (cpuid==0) then - open(unit=101, file='wanniercenter.dat') - - do iky=1, Nky - write(101, '(10000f16.8)') kpoints(2, 1, iky), & - dmod(sum(WannierCenterKy_mpi(:, iky)), 1d0), & - WannierCenterKy_mpi(:, iky) - enddo - close(101) - endif - - return - end subroutine wannier_center2D - - -! this suboutine is used for wannier center calculation for slab system -!> calculate z2 - - subroutine wannier_center2D_alt - use para - use wmpi - implicit none - - integer :: Nkx - integer :: Nky - integer :: knv2 - - integer :: i - integer :: j - integer :: l - integer :: m - integer :: ia - integer :: ia1 - integer :: nfill - - integer :: ikx - integer :: iky - - integer :: ierr - - !> k points in kx-ky plane - real(dp), allocatable :: kpoints(:, :, :) - - !> hamiltonian for each k point - !> and also the eigenvector of hamiltonian after eigensystem_c - complex(dp), allocatable :: Hamk(:, :) - complex(dp), allocatable :: Hamk_dag(:, :) - - !> eigenvector for each kx - complex(dp), allocatable :: Eigenvector(:, :, :) - - !> Mmnkb= - !> |u_n(k)> is the periodic part of wave function - complex(dp), allocatable :: Mmnkb(:, :) - complex(dp), allocatable :: Mmnkb_com(:, :) - complex(dp), allocatable :: Mmnkb_full(:, :) - - !> - complex(dp), allocatable :: Lambda_eig(:) - complex(dp), allocatable :: Lambda(:, :) - complex(dp), allocatable :: Lambda0(:, :) - - !> three matrix for SVD - !> M= U.Sigma.V^\dag - !> VT= V^\dag - complex(dp), allocatable :: U(:, :) - real (dp), allocatable :: Sigma(:, :) - complex(dp), allocatable :: VT(:, :) - - !> wannier centers for each ky, bands - real(dp), allocatable :: WannierCenterKy(:, :) - real(dp), allocatable :: WannierCenterKy_mpi(:, :) - - !> eigenvalue - real(dp), allocatable :: eigenvalue(:) - - !> atom position in the unit cell - !> for slab system, dim=Nslab*Num_atoms - real(dp), allocatable :: AtomsPosition_unitcell(:, :) - real(dp), allocatable :: AtomsPosition_supercell(:, :) - - !> for each orbital, it correspond to an atom - !> dim= Num_wann*Nslab - integer, allocatable :: AtomIndex_orbital(:) - - real(dp) :: Umatrix_t(3, 3) - - integer :: imax - real(dp) :: maxgap - real(dp) :: maxgap0 - !> b.R - real(dp) :: br - - real(dp) :: kx - real(dp) :: ky - real(dp) :: k(2) - real(dp) :: b(2) - - !> larges gap of each two wannier centers for a given k point - !> dim= Nky - real(dp), allocatable :: largestgap(:) - real(dp), allocatable :: largestgap_mpi(:) - - !> exp(-i*b.R) - complex(dp) :: ratio - - !> Z2 calculation for time reversal invariant system - integer :: Z2 - - !> Chern number - real(dp) :: Chern - - Nkx= Nk - Nky= 40 - knv2= Nkx*Nky - - nfill= Numoccupied*Nslab - - allocate(kpoints(2, Nkx, Nky)) - kpoints= 0d0 - - allocate(Lambda_eig(nfill)) - allocate(Lambda(nfill, nfill)) - allocate(Lambda0(nfill, nfill)) - allocate(Mmnkb(nfill, nfill)) - allocate(Mmnkb_com(nfill, nfill)) - allocate(Mmnkb_full(Num_wann*Nslab, Num_wann*Nslab)) - allocate(hamk(Num_wann*Nslab, Num_wann*Nslab)) - allocate(hamk_dag(Num_wann*Nslab, Num_wann*Nslab)) - allocate(Eigenvector(Num_wann*Nslab, Num_wann*Nslab, Nkx)) - allocate(eigenvalue(Num_wann*Nslab)) - allocate(U(nfill, nfill)) - allocate(Sigma(nfill, nfill)) - allocate(VT(nfill, nfill)) - allocate(WannierCenterKy(nfill, Nky)) - allocate(WannierCenterKy_mpi(nfill, Nky)) - allocate(AtomsPosition_unitcell(3, Num_atoms)) - allocate(AtomsPosition_supercell(3, Nslab*Num_atoms)) - allocate(AtomIndex_orbital(Num_wann*Nslab)) - allocate(largestgap(Nky)) - allocate(largestgap_mpi(Nky)) - largestgap= 0d0 - largestgap_mpi= 0d0 - WannierCenterKy= 0d0 - WannierCenterKy_mpi= 0d0 - hamk=0d0 - eigenvalue=0d0 - Eigenvector=0d0 - Mmnkb_full=0d0 - Mmnkb=0d0 - Mmnkb_com=0d0 - Lambda =0d0 - Lambda0=0d0 - U= 0d0 - Sigma= 0d0 - VT= 0d0 - - !> setup kpoints - do iky=1, Nky - do ikx=1, Nkx - kx= (ikx-1d0)/real(Nkx) - ky= (iky-1d0)/real(Nky) - kpoints(1, ikx, iky)= kx - kpoints(2, ikx, iky)= ky - b(1)= 1.d0/real(Nkx) - b(2)= 0.d0 - enddo - enddo - - !> set up atom index for each orbitals in the basis - if (soc>0) then !> with spin orbital coupling - l= 0 - do i=1, Nslab - do ia=1, Num_atoms !> spin up - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia+ (i-1)*Num_atoms - enddo ! l - enddo ! ia - do ia=1, Num_atoms !> spin down - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia+ (i-1)*Num_atoms - enddo ! l - enddo ! ia - enddo ! i - else !> without spin orbital coupling - l= 0 - do i=1, Nslab - do ia=1, Num_atoms !> spin down - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia+ (i-1)*Num_atoms - enddo ! l - enddo ! ia - enddo ! i - - endif - - Umatrix_t= transpose(Umatrix) - call inv_r(3, Umatrix_t) - - !> set up atoms' position in the unit cell in the new basis - !> only for 2D slab system - do ia=1, Num_atoms - do i=1, 3 - do j=1, 3 - AtomsPosition_unitcell(i, ia)= AtomsPosition_unitcell(i, ia)+ & - Umatrix_t(i, j)*Atom_position(j, ia) - enddo ! j - enddo ! i - enddo ! ia - - !> set up atoms' position in the supercell - !> actually, we only need the first two coordinates - ia1= 0 - do i=1, Nslab - do ia=1, Num_atoms - ia1= ia1+ 1 - AtomsPosition_supercell(1, ia1)= AtomsPosition_unitcell(1, ia) - AtomsPosition_supercell(2, ia1)= AtomsPosition_unitcell(2, ia) - enddo ! ia - enddo ! i - - !> for each ky, we can get wanniercenter - do iky=1+ cpuid, nky, num_cpu - Lambda0=0d0 - do i=1, nfill - Lambda0(i, i)= 1d0 - enddo - - if (cpuid==0) print *, iky, nky - !> for each kx, we get the eigenvectors - do ikx=1, nkx - k(1)= kpoints(1, ikx, iky) - k(2)= kpoints(2, ikx, iky) - - call ham_slab(k,hamk) - - !> diagonal hamk - call eigensystem_c('V', 'U', Num_wann*Nslab, hamk, eigenvalue) - - Eigenvector(:, :, ikx)= hamk - enddo - - !> sum over kx to get wanniercenters - do ikx=1, nkx - !hamk= Eigenvector(:, :, ikx) - !hamk_dag= conjg(transpose(hamk)) - !if (ikx==nkx) then - ! hamk= Eigenvector(:, :, 1) - !else - ! hamk= Eigenvector(:, :, ikx+ 1) - !endif - - !> - !call mat_mul(Num_wann*Nslab, hamk_dag, hamk, Mmnkb_full) - !Mmnkb= Mmnkb_full(1:nfill, 1:nfill) - Mmnkb= 0d0 - hamk_dag= Eigenvector(:, :, ikx) - if (ikx==nkx) then - hamk= Eigenvector(:, :, 1) - else - hamk= Eigenvector(:, :, ikx+ 1) - endif - do l=1, Nslab - do m=1, Num_wann - ia= AtomIndex_orbital(m+(l-1)*Num_wann) - br= b(1)*AtomsPosition_supercell(1, ia)+ & - b(2)*AtomsPosition_supercell(2, ia) - ratio= cos(br)- zi* sin(br) - - do i=1, nfill - do j=1, nfill - Mmnkb(i, j)= Mmnkb(i, j)+ & - conjg(hamk_dag((l-1)*Num_wann+m, i))* & - hamk((l-1)*Num_wann+m, j)* ratio - enddo ! m - enddo ! l - enddo ! j - enddo ! i - - !> perform Singluar Value Decomposed of Mmnkb - call zgesvd_pack(nfill, Mmnkb, U, Sigma, VT) - - !> after the calling of zgesvd_pack, Mmnkb becomes a temporal matrix - U= conjg(transpose(U)) - VT= conjg(transpose(VT)) - call mat_mul(nfill, VT, U, Mmnkb) - - call mat_mul(nfill, Mmnkb, Lambda0, Lambda) - Lambda0 = Lambda - enddo !< ikx - - !> diagonalize Lambda to get the eigenvalue - call zgeev_pack(nfill, Lambda, Lambda_eig) - do i=1, nfill - WannierCenterKy(i, iky)= aimag(log(Lambda_eig(i)))/2d0/pi - enddo - - call sortheap(nfill, WannierCenterKy(:, iky)) - - maxgap0= -99999d0 - do i=1, nfill - if (i/=nfill) then - maxgap= WannierCenterKy(i+1, iky)- WannierCenterKy(i, iky) - else - maxgap=1d0+ WannierCenterKy(1, iky)- WannierCenterKy(nfill, iky) - endif - - if (maxgap>maxgap0) then - maxgap0= maxgap - imax= i - endif - - enddo - - if (imax==nfill) then - largestgap(iky)= (WannierCenterKy(1, iky)+ & - WannierCenterKy(nfill, iky) -1d0)/2d0 - else - largestgap(iky)= (WannierCenterKy(imax+1, iky)+ & - WannierCenterKy(imax, iky))/2d0 - endif - enddo !< iky - - call mpi_allreduce(WannierCenterKy, WannierCenterKy_mpi, & - size(WannierCenterKy), mpi_dp, mpi_sum, mpi_cmw, ierr) - call mpi_allreduce(largestgap, largestgap_mpi, & - size(largestgap), mpi_dp, mpi_sum, mpi_cmw, ierr) - - - - if (cpuid==0) then - open(unit=101, file='wanniercenter.dat') - open(unit=102, file='largestgap.dat') - - do iky=1, Nky - write(101, '(10000f16.8)') kpoints(2, 1, iky), & - dmod(sum(WannierCenterKy_mpi(:, iky)), 1d0), & - WannierCenterKy_mpi(:, iky) - write(102, '(10000f16.8)') kpoints(2, 1, iky), & - largestgap_mpi(iky) - enddo - close(101) - close(102) - endif - - return - end subroutine wannier_center2D_alt - - - -! this suboutine is used for wannier center calculation for slab system - - subroutine wannier_center3D - use para - use wmpi - implicit none - - integer :: Nk1 - integer :: Nk2 - integer :: knv2 - - integer :: i - integer :: j - integer :: l - integer :: m - integer :: ia - integer :: ia1 - integer :: nfill - integer :: imax - - integer :: ik1 - integer :: ik2 - - integer :: ierr - - !> k points in kx-ky plane - real(dp), allocatable :: kpoints(:, :, :) - - !> hamiltonian for each k point - !> and also the eigenvector of hamiltonian after eigensystem_c - complex(dp), allocatable :: Hamk(:, :) - complex(dp), allocatable :: Hamk_dag(:, :) - - !> eigenvector for each kx - complex(dp), allocatable :: Eigenvector(:, :, :) - - !> Mmnkb= - !> |u_n(k)> is the periodic part of wave function - complex(dp), allocatable :: Mmnkb(:, :) - complex(dp), allocatable :: Mmnkb_com(:, :) - complex(dp), allocatable :: Mmnkb_full(:, :) - - !> - complex(dp), allocatable :: Lambda_eig(:) - complex(dp), allocatable :: Lambda(:, :) - complex(dp), allocatable :: Lambda0(:, :) - - !> three matrix for SVD - !> M= U.Sigma.V^\dag - !> VT= V^\dag - complex(dp), allocatable :: U(:, :) - real (dp), allocatable :: Sigma(:, :) - complex(dp), allocatable :: VT(:, :) - - !> wannier centers for each ky, bands - real(dp), allocatable :: WannierCenterKy(:, :) - real(dp), allocatable :: WannierCenterKy_mpi(:, :) - - !> larges gap of each two wannier centers for a given k point - !> dim= Nky - real(dp), allocatable :: largestgap(:) - real(dp), allocatable :: largestgap_mpi(:) - - !> eigenvalue - real(dp), allocatable :: eigenvalue(:) - - !> atom position in the unit cell - !> for slab system, dim=Num_atoms - real(dp), allocatable :: AtomsPosition_unitcell(:, :) - real(dp), allocatable :: AtomsPosition_supercell(:, :) - - !> for each orbital, it correspond to an atom - !> dim= Num_wann - integer, allocatable :: AtomIndex_orbital(:) - - real(dp) :: Umatrix_t(3, 3) - - !> b.R - real(dp) :: br - - !> exp(-i*b.R) - complex(dp) :: ratio - - real(dp) :: kx - real(dp) :: ky - real(dp) :: k(3) - real(dp) :: b(3) - - real(dp) :: maxgap - real(dp) :: maxgap0 - - !> Z2 calculation for time reversal invariant system - integer :: Z2 - integer :: Delta - - real(dp) :: g - real(dp) :: phi1 - real(dp) :: phi2 - real(dp) :: phi3 - real(dp) :: zm - real(dp) :: zm1 - real(dp) :: xnm1 - real(dp) :: Deltam - real(dp), allocatable :: xnm(:) - - - Nk1= Nk - Nk2= Nk - knv2= Nk1*Nk2 - - nfill= Numoccupied - - allocate(kpoints(2, Nk1, Nk2)) - kpoints= 0d0 - - allocate(Lambda_eig(nfill)) - allocate(Lambda(nfill, nfill)) - allocate(Lambda0(nfill, nfill)) - allocate(Mmnkb(nfill, nfill)) - allocate(Mmnkb_com(nfill, nfill)) - allocate(Mmnkb_full(Num_wann, Num_wann)) - allocate(hamk(Num_wann, Num_wann)) - allocate(hamk_dag(Num_wann, Num_wann)) - allocate(Eigenvector(Num_wann, Num_wann, Nk1)) - allocate(eigenvalue(Num_wann)) - allocate(U(nfill, nfill)) - allocate(Sigma(nfill, nfill)) - allocate(VT(nfill, nfill)) - allocate(WannierCenterKy(nfill, Nk2)) - allocate(WannierCenterKy_mpi(nfill, Nk2)) - allocate(AtomIndex_orbital(Num_wann)) - allocate(xnm(nfill)) - allocate(largestgap(Nk2)) - allocate(largestgap_mpi(Nk2)) - largestgap= 0d0 - largestgap_mpi= 0d0 - WannierCenterKy= 0d0 - WannierCenterKy_mpi= 0d0 - hamk=0d0 - eigenvalue=0d0 - Eigenvector=0d0 - Mmnkb_full=0d0 - Mmnkb=0d0 - Mmnkb_com=0d0 - Lambda =0d0 - Lambda0=0d0 - U= 0d0 - Sigma= 0d0 - VT= 0d0 - - !> setup kpoints - do ik2=1, Nk2 - do ik1=1, Nk1 - kx= (ik1-1d0)/real(Nk1) - ky= (ik2-1d0)/real(Nk2-1)/2d0 - kpoints(1, ik1, ik2)= kx - kpoints(2, ik1, ik2)= ky - b(1)= 1.d0/real(Nk1) - b(2)= 0.d0 - b(3)= 0.d0 - enddo - enddo - - !> set up atom index for each orbitals in the basis - if (soc>0) then !> with spin orbital coupling - l= 0 - do ia=1, Num_atoms !> spin up - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia - enddo ! l - enddo ! ia - do ia=1, Num_atoms !> spin down - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia - enddo ! l - enddo ! ia - else !> without spin orbital coupling - l= 0 - do ia=1, Num_atoms !> spin down - do j=1, nprojs(ia) - l= l+ 1 - AtomIndex_orbital(l)= ia - enddo ! l - enddo ! ia - - endif - - Umatrix_t= transpose(Umatrix) - call inv_r(3, Umatrix_t) - - !>> Get wannier center for ky=0 plane - !> for each ky, we can get wanniercenter - do ik2=1+ cpuid, Nk2, num_cpu - Lambda0=0d0 - do i=1, nfill - Lambda0(i, i)= 1d0 - enddo - - !> for each k1, we get the eigenvectors - do ik1=1, Nk1 - if (cpuid==0) print *, ik1, ik2, Nk1, Nk2 - k(1)= kpoints(1, ik1, ik2) - k(2)= 0d0 - k(3)= kpoints(2, ik1, ik2) - - call ham_bulk(k,hamk) - - !> diagonal hamk - call eigensystem_c('V', 'U', Num_wann, hamk, eigenvalue) - - Eigenvector(:, :, ik1)= hamk - enddo - - !> sum over k1 to get wanniercenters - do ik1=1, Nk1 - !> - Mmnkb= 0d0 - hamk_dag= Eigenvector(:, :, ik1) - if (ik1==Nk1) then - hamk= Eigenvector(:, :, 1) - else - hamk= Eigenvector(:, :, ik1+ 1) - endif - do m=1, Num_wann - ia= AtomIndex_orbital(m) - br= b(1)*Atom_position(1, ia)+ & - b(2)*Atom_position(2, ia)+ & - b(3)*Atom_position(3, ia) - !ratio= cos(br)- zi* sin(br) - ratio= 1d0 - - do j=1, nfill - do i=1, nfill - Mmnkb(i, j)= Mmnkb(i, j)+ & - conjg(hamk_dag(m, i))* hamk(m, j)* ratio - enddo ! i - enddo ! j - enddo ! m - - !> perform Singluar Value Decomposed of Mmnkb - call zgesvd_pack(nfill, Mmnkb, U, Sigma, VT) - - !> after the calling of zgesvd_pack, Mmnkb becomes a temporal matrix - U= conjg(transpose(U)) - VT= conjg(transpose(VT)) - call mat_mul(nfill, VT, U, Mmnkb) - - call mat_mul(nfill, Mmnkb, Lambda0, Lambda) - Lambda0 = Lambda - enddo !< ik1 - - !> diagonalize Lambda to get the eigenvalue - call zgeev_pack(nfill, Lambda, Lambda_eig) - do i=1, nfill - WannierCenterKy(i, ik2)= aimag(log(Lambda_eig(i)))/2d0/pi - WannierCenterKy(i, ik2)= mod(WannierCenterKy(i, ik2)+10d0, 1d0) - enddo - - call sortheap(nfill, WannierCenterKy(:, ik2)) - - maxgap0= -99999d0 - imax= nfill - do i=1, nfill - if (i/=nfill) then - maxgap= WannierCenterKy(i+1, ik2)- WannierCenterKy(i, ik2) - else - maxgap=1d0+ WannierCenterKy(1, ik2)- WannierCenterKy(nfill, ik2) - endif - - if (maxgap>maxgap0) then - maxgap0= maxgap - imax= i - endif - - enddo - - if (imax==nfill) then - largestgap(ik2)= (WannierCenterKy(1, ik2)+ & - WannierCenterKy(nfill, ik2) +1d0)/2d0 - largestgap(ik2)= mod(largestgap(ik2), 1d0) - else - largestgap(ik2)= (WannierCenterKy(imax+1, ik2)+ & - WannierCenterKy(imax, ik2))/2d0 - endif - - - enddo !< ik2 - - WannierCenterKy_mpi= 0d0 - largestgap_mpi= 0d0 - call mpi_allreduce(WannierCenterKy, WannierCenterKy_mpi, & - size(WannierCenterKy), mpi_dp, mpi_sum, mpi_cmw, ierr) - call mpi_allreduce(largestgap, largestgap_mpi, & - size(largestgap), mpi_dp, mpi_sum, mpi_cmw, ierr) - - if (cpuid==0) then - open(unit=101, file='wanniercenterky0.dat') - - do ik2=1, Nk2 - write(101, '(10000f16.8)') kpoints(2, 1, ik2), & - largestgap_mpi(ik2), dmod(sum(WannierCenterKy_mpi(:, ik2)), 1d0), & - WannierCenterKy_mpi(:, ik2) - enddo - close(101) - endif - - - !> Z2 calculation Alexey Soluyanov arXiv:1102.5600 - - Delta= 0 - !> for each iky, we get a Deltam - do ik2=1, nk2-1 - - !> largestgap position - zm= largestgap_mpi(ik2) - !if (ik2==nk2) then - ! zm1= largestgap_mpi(1) - ! xnm= WannierCenterKy_mpi(1:nfill, 1) - !else - zm1= largestgap_mpi(ik2+1) - xnm= WannierCenterKy_mpi(1:nfill, ik2+1) - !endif - Deltam= 1 - do i=1, nfill - xnm1= xnm(i) - phi1= 2d0*pi*zm - phi2= 2d0*pi*zm1 - phi3= 2d0*pi*xnm1 - - g= sin(phi2-phi1)+ sin(phi3-phi2)+ sin(phi1-phi3) - Deltam= Deltam* sign(1d0, g) - enddo !i - if (Deltam<0) then - Delta= Delta+ 1 - endif - enddo !ik2 - - Z2= mod(Delta, 2) - - if (cpuid==0) print*,'Z2 for ky=0: ', Z2, Delta - - - !>> Get wannier center for ky=0.5 plane - !> for each ky, we can get wanniercenter - do ik2=1+ cpuid, Nk2, num_cpu - Lambda0=0d0 - do i=1, nfill - Lambda0(i, i)= 1d0 - enddo - - !> for each k1, we get the eigenvectors - do ik1=1, Nk1 - if (cpuid==0) print *, ik1, ik2, Nk1, Nk2 - k(1)= kpoints(1, ik1, ik2) - k(2)= 0.5d0 - k(3)= kpoints(2, ik1, ik2) - - call ham_bulk(k,hamk) - - !> diagonal hamk - call eigensystem_c('V', 'U', Num_wann, hamk, eigenvalue) - - Eigenvector(:, :, ik1)= hamk - enddo - - !> sum over k1 to get wanniercenters - do ik1=1, Nk1 - !> - Mmnkb= 0d0 - hamk_dag= Eigenvector(:, :, ik1) - if (ik1==Nk1) then - hamk= Eigenvector(:, :, 1) - else - hamk= Eigenvector(:, :, ik1+ 1) - endif - do m=1, Num_wann - ia= AtomIndex_orbital(m) - br= b(1)*Atom_position(1, ia)+ & - b(2)*Atom_position(2, ia)+ & - b(3)*Atom_position(3, ia) - !ratio= cos(br)- zi* sin(br) - ratio= 1d0 - - do j=1, nfill - do i=1, nfill - Mmnkb(i, j)= Mmnkb(i, j)+ & - conjg(hamk_dag(m, i))* hamk(m, j)* ratio - enddo ! i - enddo ! j - enddo ! m - - !> perform Singluar Value Decomposed of Mmnkb - call zgesvd_pack(nfill, Mmnkb, U, Sigma, VT) - - !> after the calling of zgesvd_pack, Mmnkb becomes a temporal matrix - U= conjg(transpose(U)) - VT= conjg(transpose(VT)) - call mat_mul(nfill, VT, U, Mmnkb) - - call mat_mul(nfill, Mmnkb, Lambda0, Lambda) - Lambda0 = Lambda - enddo !< ik1 - - !> diagonalize Lambda to get the eigenvalue - call zgeev_pack(nfill, Lambda, Lambda_eig) - do i=1, nfill - WannierCenterKy(i, ik2)= aimag(log(Lambda_eig(i)))/2d0/pi - WannierCenterKy(i, ik2)= mod(WannierCenterKy(i, ik2)+10d0, 1d0) - enddo - - call sortheap(nfill, WannierCenterKy(:, ik2)) - - maxgap0= -99999d0 - imax= nfill - do i=1, nfill - if (i/=nfill) then - maxgap= WannierCenterKy(i+1, ik2)- WannierCenterKy(i, ik2) - else - maxgap=1d0+ WannierCenterKy(1, ik2)- WannierCenterKy(nfill, ik2) - endif - - if (maxgap>maxgap0) then - maxgap0= maxgap - imax= i - endif - - enddo - - if (imax==nfill) then - largestgap(ik2)= (WannierCenterKy(1, ik2)+ & - WannierCenterKy(nfill, ik2) +1d0)/2d0 - largestgap(ik2)= mod(largestgap(ik2), 1d0) - else - largestgap(ik2)= (WannierCenterKy(imax+1, ik2)+ & - WannierCenterKy(imax, ik2))/2d0 - endif - - enddo !< ik2 - - WannierCenterKy_mpi= 0d0 - largestgap_mpi= 0d0 - call mpi_allreduce(WannierCenterKy, WannierCenterKy_mpi, & - size(WannierCenterKy), mpi_dp, mpi_sum, mpi_cmw, ierr) - call mpi_allreduce(largestgap, largestgap_mpi, & - size(largestgap), mpi_dp, mpi_sum, mpi_cmw, ierr) - - if (cpuid==0) then - open(unit=101, file='wanniercenterky05.dat') - - do ik2=1, Nk2 - write(101, '(10000f16.8)') kpoints(2, 1, ik2), & - largestgap_mpi(ik2), & - dmod(sum(WannierCenterKy_mpi(:, ik2)), 1d0), & - WannierCenterKy_mpi(:, ik2) - enddo - close(101) - endif - - - !> Z2 calculation Alexey Soluyanov arXiv:1102.5600 - - Delta= 0 - !> for each iky, we get a Deltam - do ik2=1, nk2- 1 - - !> largestgap position - zm= largestgap_mpi(ik2) - zm1= largestgap_mpi(ik2+1) - xnm= WannierCenterKy_mpi(1:nfill, ik2+1) - Deltam= 1 - do i=1, nfill - xnm1= xnm(i) - phi1= 2*pi*zm - phi2= 2*pi*zm1 - phi3= 2*pi*xnm1 - - g= sin(phi2-phi1)+ sin(phi3-phi2)+ sin(phi1-phi3) - Deltam= Deltam* sign(1d0, g) - enddo !i - if (Deltam<0) then - Delta= Delta+ 1 - endif - enddo !ik2 - - Z2= mod(Delta, 2) - - if (cpuid==0) print*,'Z2 for ky=0.5: ', Z2 - - return - end subroutine wannier_center3D - - subroutine sortheap(n, arr) - implicit none - integer, intent(in) :: n - real(8), intent(inout) :: arr(n) - - !* local variables - integer :: i - - do i=n/2, 1, -1 - call sift_down(i, n) - enddo - - do i=n, 2, -1 - call swap(arr(1), arr(i)) - call sift_down(1, i-1) - enddo - contains - subroutine sift_down(l, r) - integer, intent(in) :: l, r - integer :: j, jold - real(8) :: a - a= arr(l) - jold= l - j= l+ l - - do while (j<=r) - if (j= arr(j))exit - arr(jold)= arr(j) - jold= j - j= j+ j - enddo - arr(jold)= a - return - end subroutine sift_down - end subroutine sortheap - - !>> swap two real numbers - subroutine swap(a, b) - real(8), intent(inout) :: a - real(8), intent(inout) :: b - real(8) :: c - c=a - a=b - b=c - return - end subroutine swap - -