fixed several bugs in sigma_OHE.f90

examples/Graphene/wt.in-landaulevel

@@ -3,10 +3,10 @@ Hrfile = 'Graphene_hr.dat'
 /
 
 &CONTROL
-LandauLevel_B_dos_calc = T   ! Hofstadter butterfly by Lanczos method
-LandauLevel_B_calc = F  ! Hofstadter butterfly by diagonalization
-LandauLevel_k_calc = F  ! energy bands at a fixed magnetic flux magp/Nslab*\Phi_0 using diagonalization method
-LandauLevel_k_dos_calc = F  ! energy bands at a fixed magnetic flux magp/Nslab*\Phi_0 using Laczos method
+LandauLevel_B_dos_calc = F   ! Hofstadter butterfly by Lanczos method
+LandauLevel_B_calc = T       ! Hofstadter butterfly by diagonalization
+LandauLevel_k_calc = F       ! energy bands at a fixed magnetic flux magp/Nslab*\Phi_0 using diagonalization method
+LandauLevel_k_dos_calc = F   ! energy bands at a fixed magnetic flux magp/Nslab*\Phi_0 using Laczos method
 /
 
 &SYSTEM

src/Makefile.gfortran

@@ -12,17 +12,18 @@ obj =  module.o sparse.o wt_aux.o math_lib.o symmetry.o readHmnR.o inverse.o pro
        main.o
 
 # compiler
-f90  = gfortran -cpp                                              
+f90  = gfortran -cpp -DINTELMKL
 
 flag = -O3 -ffree-line-length-512  #-g -Wall -Wextra -Warray-temporaries -Wconversion -fimplicit-none -fbacktrace -ffree-line-length-0 -fcheck=all -ffpe-trap=zero,overflow,underflow -finit-real=nan
 
+ARPACK=/Users/quan/work/workplace/ARPACK/libarpack_gfortran.a
 
 # blas and lapack libraries
-libs = -L/usr/lib/ -llapack -lblas
+#libs = -L/usr/lib/ -llapack -lblas
 
 # Intel MKL also supports with gfortran comipler
 # dynamic linking
-# libs = -L/${MKLROOT}/lib/intel64 -lmkl_core -lmkl_sequential -lmkl_intel_lp64 -lpthread
+libs = ${ARPACK} -L/${MKLROOT}/lib -lmkl_core -lmkl_sequential -lmkl_intel_lp64 -lpthread
 
 
 main :  $(obj)

src/Makefile.intel-mpi

@@ -20,17 +20,18 @@ OFLAG = -O3 -static-intel
 FFLAG = $(OFLAG) $(WFLAG)
 LFLAG = $(OFLAG)
 
+# ARPACK LIBRARY
 ARPACK=/Users/quan/work/workplace/ARPACK/libarpack_MAC.a
 
 # blas and lapack libraries
 # static linking
-LIBS = -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \
+#LIBS = -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \
         ${MKLROOT}/lib/intel64/libmkl_sequential.a \
         ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl \
 		  ${ARPACK}
 
 # dynamic linking
-# LIBS = -L/${MKLROOT}/lib/intel64 -lmkl_core -lmkl_sequential -lmkl_intel_lp64 -lpthread
+LIBS = ${ARPACK} -L/${MKLROOT}/lib/intel64 -lmkl_core -lmkl_sequential -lmkl_intel_lp64 -lpthread
 
 main : $(OBJ)
 	$(F90) $(LFLAG) $(OBJ) -o wt.x $(LIBS)

src/eigen.f90

@@ -280,6 +280,7 @@ subroutine zheevx_pack(JOBZ, UPLO, N, il, iu, A, eigval, eigvec)
 
      mdim= iu-il+1
      lwork= 64*N
+     vl=0d0; vu=0d0 
 
      allocate(ifail(N))
      allocate(iwork(5*N))
@@ -295,9 +296,14 @@ subroutine zheevx_pack(JOBZ, UPLO, N, il, iu, A, eigval, eigvec)
      eigvec= 0d0
      vl=0d0; vu=0d0
 
+     call zheevx(JOBZ,'I',UPLO,N,A,N,vl,vu,il,iu,abstol,&
+         mdim,eigenvalues,eigvec,N,work,-1,rwork,iwork,ifail,info)
+     lwork=int(work(1))
+
      call zheevx(JOBZ,'I',UPLO,N,A,N,vl,vu,il,iu,abstol,&
          mdim,eigenvalues,eigvec,N,work,lwork,rwork,iwork,ifail,info)
 
+
      if (info/=0) then
         print *, ' Error info in zheevx: ', info
         stop ' Error happens in zheev_pack'
@@ -310,5 +316,91 @@ subroutine zheevx_pack(JOBZ, UPLO, N, il, iu, A, eigval, eigvec)
      return
   end subroutine  zheevx_pack
 
+  subroutine zgeev_sys(N, A, W,JOBVL,VL,JOBVR,VR )
+     ! a pack of Lapack subroutine zgeev, which is a subroutine to
+     ! calculate eigenvector and eigenvalue of a complex hermite matrix
+
+     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 is COMPLEX*16 array, dimension (LDA,N)
+!~           On entry, the N-by-N matrix A.
+!~           On exit, A has been overwritten.
+
+     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) :: VL(N,N)
+
+!    right eigenvectors
+     complex(dp) :: VR(N,N)
+
+     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(rwork(16*N))
+     allocate( work(lwork))
+     VL= (0d0, 0d0)
+     VR= (0d0, 0d0)
+     rwork= 0d0
+     work= (0d0, 0d0)
+
+     info=0
+     W=(0d0, 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
 
 
+     deallocate(rwork)
+     deallocate( work)
+
+     return
+  end subroutine zgeev_sys
+

src/ek_bulk.f90

@@ -2757,7 +2757,7 @@ subroutine generate_ek_kpath_gnu(datafilename, gnufilename, gnuoutfilename, &
    outfileindex= outfileindex+ 1
    if (cpuid==0) then
       open(unit=outfileindex, file=gnufilename)
-      write(outfileindex, '(a)') 'set terminal pdf enhanced color font ",30"'
+      write(outfileindex, '(a)') 'set terminal pdf enhanced color font ",24"'
       write(outfileindex,'(2a)') 'set palette defined ( 0  "green", ', &
          '5 "yellow", 10 "red" )'
       write(outfileindex, '(3a)')"set output '", trim(adjustl(gnuoutfilename)), "' "
@@ -2766,8 +2766,8 @@ subroutine generate_ek_kpath_gnu(datafilename, gnufilename, gnuoutfilename, &
       write(outfileindex, '(a)')'set pointsize 0.8'
       write(outfileindex, '(a)')'#set xtics font ",24"'
       write(outfileindex, '(a)')'#set ytics font ",24"'
-      write(outfileindex, '(a)')'set ylabel font ",24"'
-      write(outfileindex, '(a)')'set ylabel offset 1.5,0'
+      write(outfileindex, '(a)')'#set ylabel font ",24"'
+      write(outfileindex, '(a)')'set ylabel offset 0.5,0'
       write(outfileindex, '(a, f10.5, a)')'set xrange [0: ', maxval(klen*Angstrom2atomic), ']'
       write(outfileindex, '(a,f12.6)')'emin=', emin
       write(outfileindex, '(a,f12.6)')'emax=', emax

src/inverse.f90

@@ -105,7 +105,7 @@ subroutine inv_r(ndim,Amat)
 
      call dgesv(ndim,ndim,Amat,ndim,ipiv,Bmat,ndim,info)
 
-     if(info.ne.0)print *,'something wrong with dgesv'
+     if(info.ne.0)stop 'something wrong with dgesv in inverse.f90'
 
      Amat=Bmat
 

src/main.f90

@@ -21,6 +21,7 @@
 ! version     2.5.0  At EPFL, Switzerland, Dec. 9. 2019, magnetoresistance, band unfolding
 ! version     2.5.1  At EPFL, Switzerland, Mar. 6. 2020, For WannierTools tutorial 2020
 ! version     2.6.0  At EPFL, Switzerland, Feb.15. 2021, Landau level, sparse Hamiltonian, TBG
+! version     2.6.1  At Beijing, China, April 10. 2022 clean for the Wannier90 tutorial 2022
 !
 ! Corresponding to : [email protected]
 !
@@ -43,7 +44,7 @@ program main
 
 
      !> version of WannierTools
-     version='2.6.0'
+     version='2.6.1'
 
      ierr = 0
      cpuid= 0

src/readinput.f90

@@ -221,7 +221,7 @@ subroutine readinput
    BulkBand_line_calc= BulkBand_line_calc.or.BulkBand_calc
    BulkBand_unfold_line_calc= BulkBand_unfold_line_calc.or.Landaulevel_unfold_line_calc
 
-   if (MirrorChern_calc.or.Boltz_OHE_calc) Symmetry_Import_calc = .true.
+   if (MirrorChern_calc) Symmetry_Import_calc = .true.
 
    !> control parameters
    if (cpuid==0) then
@@ -278,6 +278,7 @@ subroutine readinput
       write(stdout, *) "LandauLevel_k_dos_calc            : ", LandauLevel_k_dos_calc
       write(stdout, *) "LandauLevel_B_dos_calc            : ", LandauLevel_B_dos_calc
       write(stdout, *) "FermiLevel_calc                   : ", FermiLevel_calc
+      write(stdout, *) "Symmetry_Import_calc              : ", Symmetry_Import_calc
       write(stdout, *) "ChargeDensity_selected_bands_calc : ", ChargeDensity_selected_bands_calc
       write(stdout, *) "ChargeDensity_selected_energies_calc : ", ChargeDensity_selected_energies_calc
    endif
@@ -3828,7 +3829,7 @@ subroutine cart_direct_rec(k1, k2)
    mata(3, :)= Origin_cell%Kuc
 
    call inv_r(3, mata)
-   K2= k1(1)*mata(1, :)+ k1(2)*mata(2, :)+ k1(3)*mata(3, :)
+   k2= k1(1)*mata(1, :)+ k1(2)*mata(2, :)+ k1(3)*mata(3, :)
 
    deallocate(mata)
 

src/sigma.f90

@@ -159,23 +159,8 @@ subroutine sigma_resistivity
       time_start= 0d0
       time_end= 0d0
 
-      if (Symmetry_Import_calc) then
-         call sigma_ohe_calc_symm(mu_array, KBT_array, BTau_array, Nband_Fermi_Level, bands_fermi_level, sigma_ohe_tensor)
-      else
-         call sigma_ohe_calc     (mu_array, KBT_array, BTau_array, Nband_Fermi_Level, bands_fermi_level, sigma_ohe_tensor)
-      endif
-     !if (cpuid.eq.0) then
-     !   if (NumT==1.and.OmegaNum==1)then
-     !      !> write out the conductivity/tau into file
-     !      do i=1, NBTau
-     !      do ib=1, Nband_Fermi_Level
-     !         write(myfileindex(ib), '(20E16.6)')BTau_array(i), BTau_array(i)*0.175874356d0, sigma_ohe_tensor(:, i, 1, 1, ib)
-     !      enddo ! ib, band
-     !      enddo ! i, NBTau
-     !   endif !> only one T and one mu
-     !endif ! cpuid=0
-
-
+      call sigma_ohe_calc_symm(mu_array, KBT_array, BTau_array, Nband_Fermi_Level, bands_fermi_level, sigma_ohe_tensor)
+
      !if (cpuid.eq.0) then
      !   do ib=1, Nband_Fermi_Level
      !      close(myfileindex(ib))