start to add supports for non-orthogonal basis such as Hamiltonian fr…

…om OpenMX

examples/Bi2Se3/wt.in

@@ -29,7 +29,8 @@ Se pz px py
 !> bulk band structure calculation flag
 &CONTROL
 BulkBand_calc         = T
-BulkBand_points_calc  = T
+!BulkBand_points_calc  = T
+Symmetry_Import_calc = T
 /
 
 &SYSTEM

examples/MoTe2/bulk/readme.txt

@@ -0,0 +1,9 @@
+An example to obtain band structure of MoTe2 with the Hamiltonian obtained from the
+OpenMX with non-orthogonal basis. 
+
+The Hamiltonain is stored as two parts in sparse format, H.dat and S.dat.
+
+H.dat contains the Hamiltonian Hmn(R)=<\phi_m(r)|H|\phi_n(r-R)>, only the non-zero entries are stored.
+S.dat contains the overlap matrix Smn(R)=<\phi_m(r)|\phi_n(r-R)>, only the non-zero entries are stored
+
+tar xzvf HS.tar.gz 

examples/MoTe2/bulk/wt.in

@@ -0,0 +1,78 @@
+&TB_FILE
+Hrfile = 'H.dat'
+Overlapfile = 'S.dat'
+Package = 'OPENMX'
+Is_Sparse_Hr = T
+Orthogonal_Basis = F
+/
+
+
+!> bulk band structure calculation flag
+&CONTROL
+BulkBand_calc         = T
+DOS_calc              = F
+SlabBand_calc         = F
+SlabSS_calc           = F
+wanniercenter_calc    = F
+/
+
+&SYSTEM
+NumOccupied = 1         ! NumOccupied
+SOC = 1                 ! soc
+/
+
+&PARAMETERS
+E_arc = -5          ! specify the Fermi energy
+OmegaNum = 32       ! omega number       
+OmegaMin = -0.6     ! energy interval
+OmegaMax =  0.2     ! energy interval
+Nk1 = 10           ! number k points 
+NP = 1              ! number of principle layers
+/
+
+LATTICE
+Angstrom
+     3.5228714356697486    0.0001126967654733   -0.0000045697261942
+    -1.7613380944900303    3.0509520831236592    0.0000043347633975
+    -0.0000361780629472    0.0000168923630182   27.7486008981945815
+
+ATOM_POSITIONS
+6                               ! number of atoms for projectors
+Direct                          ! Direct or Cartisen coordinate
+Mo   0.6666761284131021  0.3333234129710836  0.3966008567881620 
+Mo   0.3333225911561447  0.6666760972210506  0.1441484616741683 
+Te   0.3333296102529972  0.6666699314530090  0.3312669719246579 
+Te   0.6666646131697589  0.3333343705209685  0.0788668529700386 
+Te   0.6666682236074797  0.3333310117300143  0.2094820551190944 
+Te   0.3333330967010965  0.6666663334757505  0.4618839571489489 
+
+PROJECTORS
+2*14 4*19                ! number of projectors
+Mo s s s px py pz px py pz dz2  dx2-y2 dxy dxz dyz
+Mo s s s px py pz px py pz dz2  dx2-y2 dxy dxz dyz
+Te s s s px py pz px py pz dz2  dx2-y2 dxy dxz dyz dz2  dx2-y2 dxy dxz dyz
+Te s s s px py pz px py pz dz2  dx2-y2 dxy dxz dyz dz2  dx2-y2 dxy dxz dyz
+Te s s s px py pz px py pz dz2  dx2-y2 dxy dxz dyz dz2  dx2-y2 dxy dxz dyz
+Te s s s px py pz px py pz dz2  dx2-y2 dxy dxz dyz dz2  dx2-y2 dxy dxz dyz
+
+SURFACE            ! MoS2 conventional (010) surface
+ 0  1  0
+ 0  0  1
+ 1  0  0
+
+KPATH_BULK            ! k point path
+3              ! number of k line only for bulk band
+  G 0.00000  0.00000  0.00000   M   0.50000  0.00000  0.00000   
+  M 0.50000  0.00000  0.00000   K   0.33333  0.33333  0.00000   
+  K 0.33333  0.33333  0.00000   G   0.00000  0.00000  0.00000   
+
+KPATH_SLAB
+2        ! numker of k line for 2D case
+X -0.5  0.0 G  0.0  0.0  ! k path for 2D case
+G  0.0  0.0 X  0.5  0.0 
+
+KPLANE_BULK
+ 0.00  0.00  0.00   ! Original point for 3D k plane 
+ 1.00  0.00  0.00   ! The first vector to define 3d k space plane
+ 0.00  0.50  0.00   ! The second vector to define 3d k space plane
+

examples/TBG-7.34degree/bulkek.gnu0

@@ -20,7 +20,8 @@ set arrow from    0.32755, emin to    0.32755, emax nohead
 #plot 'bulkek.dat' u 1:2  w p pt 7  ps 0.2 lc rgb 'black', 0 w l lw 2
 
 # uncomment the following lines to plot the fatband 
-plot 'bulkek.dat' u 1:2:3  w lp lw 2 pt 7  ps 0.2 lc palette, 0 w l lw 2
+plot 'bulkek_valley_plus.dat' u 1:2:3  w p  pt 7  ps 0.2 lc palette, \
+     'bulkek_valley_minus.dat' u 1:2:3  w p  pt 7  ps 0.2 lc palette, 0 w l lw 2
 # uncomment the following lines to plot the spin if necessary
 #plot 'bulkek.dat' u 1:2 w lp lw 2 pt 7  ps 0.2, \
      'bulkek.dat' u 1:2:($3/6):($4/6) w vec

src/ek_bulk.f90

@@ -64,15 +64,16 @@ subroutine ek_bulk_line
       W= 0d0
       call eigensystem_c('V', 'U', Num_wann ,Hamk_bulk, W)
       eigv(:, ik)= W
-      do j=1, Num_wann  !> band
-           do ig=1, NumberofSelectedOrbitals_groups
-              do i=1, NumberofSelectedOrbitals(ig)
-                 io=Selected_WannierOrbitals(ig)%iarray(i)
+
+      do j= 1, Num_wann  !> band
+           do ig= 1, NumberofSelectedOrbitals_groups
+              do i= 1, NumberofSelectedOrbitals(ig)
+                 io= Selected_WannierOrbitals(ig)%iarray(i)
                  weight(ig, j, ik)= weight(ig, j, ik)+ abs(Hamk_bulk(io, j))**2 
-              enddo
-           enddo
-      enddo ! i
-   enddo ! ik
+              enddo !i
+           enddo !ig
+      enddo !j
+   enddo !ik
 
 #if defined (MPI)
    call mpi_allreduce(eigv,eigv_mpi,size(eigv),&
@@ -1133,6 +1134,13 @@ subroutine sparse_ekbulk
    integer, allocatable :: jcoo(:)
    integer, allocatable :: icoo(:)
 
+   !> storage for overlap matrix
+   integer :: snnzmax, snnz
+   complex(dp), allocatable :: sacoo_k(:)
+   integer, allocatable :: sjcoo_k(:)
+   integer, allocatable :: sicoo_k(:)
+
+
    !> eigenvector of the sparse matrix acoo. Dim=(Num_wann, neval)
    complex(dp), allocatable :: psi(:)
    complex(dp), allocatable :: psi_project(:)
@@ -1187,12 +1195,22 @@ subroutine sparse_ekbulk
       numk=1
    endif
 
-   sigma=(1d0,0d0)*E_arc
    nnzmax=splen+Num_wann
+   if (.not.Orthogonal_Basis) then
+      nnzmax=splen+Num_wann+splen_overlap_input
+   else
+      nnzmax=splen+Num_wann
+   endif
    nnz=splen
+   snnzmax=splen_overlap_input
    allocate( acoo(nnzmax))
    allocate( jcoo(nnzmax))
    allocate( icoo(nnzmax))
+   if (.not.Orthogonal_Basis) then
+      allocate( sacoo_k(snnzmax))
+      allocate( sjcoo_k(snnzmax))
+      allocate( sicoo_k(snnzmax))
+   endif
    allocate( W( neval))
    allocate( eigv( neval, nk3_band))
    allocate( eigv_mpi( neval, nk3_band))
@@ -1207,24 +1225,41 @@ subroutine sparse_ekbulk
 
    eigv_mpi= 0d0;  eigv    = 0d0
    acoo= 0d0; icoo=0; jcoo=0
+   sacoo_k= 0d0; sicoo_k=0; sjcoo_k=0
 
    ritzvec= BulkFatBand_calc
 
    !> calculate the energy bands along special k line
    k3= 0
+
+   !> change the energy unit from Hatree to eV
+   sigma=(1d0,0d0)*E_arc/eV2Hartree
+
    do ik=1+ cpuid, nk3_band, num_cpu
       if (cpuid==0) write(stdout, '(a, 2i10)') 'BulkBand_calc in sparse mode:', ik,nk3_band
       k3 = K3points(:, ik)
       call now(time1)
+      !> use atomic gauge here
       call ham_bulk_coo_sparsehr(k3,acoo,icoo,jcoo)
+      if (.not.Orthogonal_Basis) then
+         call overlap_bulk_coo_sparse(k3, sacoo_k, sicoo_k, sjcoo_k)
+      endif
+      !> change the energy unit from Hatree to eV
       acoo= acoo/eV2Hartree
       nnz= splen
+      snnz=splen_overlap_input
       call now(time2)
 
       !> diagonalization by call zheev in lapack
       W= 0d0
       !> after arpack_sparse_coo_eigs, nnz will be updated.
-      call arpack_sparse_coo_eigs(Num_wann,nnzmax,nnz,acoo,jcoo,icoo,neval,nvecs,W,sigma, zeigv, ritzvec)
+      if (.not.Orthogonal_Basis) then
+         !> non-orthogonal basis like, openmx
+         call arpack_sparse_coo_eigs_nonorth(Num_wann, nnzmax, nnz, acoo, jcoo, icoo, &
+             snnzmax, snnz, sacoo_k, sjcoo_k, sicoo_k, neval,nvecs,W,sigma,zeigv, ritzvec)
+      else
+         call arpack_sparse_coo_eigs(Num_wann,nnzmax,nnz,acoo,jcoo,icoo,neval,nvecs,W,sigma, zeigv, ritzvec)
+      endif
       call now(time3)
       eigv(1:neval, ik)= W(1:neval)
 
@@ -1343,6 +1378,13 @@ subroutine sparse_ekbulk_valley
    integer, allocatable :: jcoo(:), jcoo_valley(:)
    integer, allocatable :: icoo(:), icoo_valley(:)
 
+   !> storage for overlap matrix
+   integer :: snnzmax, snnz
+   complex(dp), allocatable :: sacoo_k(:)
+   integer, allocatable :: sjcoo_k(:)
+   integer, allocatable :: sicoo_k(:)
+
+
    !> eigenvector of the sparse matrix acoo. Dim=(Num_wann, neval)
    complex(dp), allocatable :: psi_project(:)
    complex(dp), allocatable :: zeigv(:, :)
@@ -1401,14 +1443,25 @@ subroutine sparse_ekbulk_valley
    if (nvecs<20) nvecs= 20
    if (nvecs>Num_wann) nvecs= Num_wann
 
-   sigma=(1d0,0d0)*E_arc
    nnzmax=splen+Num_wann
+   if (.not.Orthogonal_Basis) then
+      nnzmax=splen+Num_wann+splen_overlap_input
+   else
+      nnzmax=splen+Num_wann
+   endif
    nnz=splen
+   snnzmax=splen_overlap_input
+
    NDMAX= 12
    allocate( acoo(nnzmax))
    allocate( jcoo(nnzmax))
    allocate( icoo(nnzmax))
-
+   if (.not.Orthogonal_Basis) then
+      allocate( sacoo_k(snnzmax))
+      allocate( sjcoo_k(snnzmax))
+      allocate( sicoo_k(snnzmax))
+   endif
+
    nnzmax_valley = splen_valley_input
    allocate( acoo_valley(nnzmax_valley))
    allocate( jcoo_valley(nnzmax_valley))
@@ -1437,24 +1490,38 @@ subroutine sparse_ekbulk_valley
 
    ritzvec= .True.
 
+   !> change the energy unit from Hatree to eV
+   sigma=(1d0,0d0)*E_arc/eV2Hartree
+
    !> calculate the energy bands along special k line
    k3= 0
    do ik=1+ cpuid, nk3_band, num_cpu
       if (cpuid==0) write(stdout, '(a, 2i10)') 'BulkBand_calc in sparse mode:', ik,nk3_band
       k3 = K3points(:, ik)
       call now(time1)
       call ham_bulk_coo_sparsehr(k3,acoo,icoo,jcoo)
+      !> change the energy unit from Hatree to eV
       acoo= acoo/eV2Hartree
       nnz= splen
+      if (.not.Orthogonal_Basis) then
+         call overlap_bulk_coo_sparse(k3, sacoo_k, sicoo_k, sjcoo_k)
+      endif
       call now(time2)
-      
+
       !> diagonalization by call zheev in lapack
       W= 0d0
       !> after arpack_sparse_coo_eigs, nnz will be updated.
-      call arpack_sparse_coo_eigs(Num_wann,nnzmax,nnz,acoo,jcoo,icoo,neval,nvecs,W,sigma, zeigv, ritzvec)
+      if (.not.Orthogonal_Basis) then
+         !> non-orthogonal basis like, openmx
+         call arpack_sparse_coo_eigs_nonorth(Num_wann, nnzmax, nnz, acoo, jcoo, icoo, &
+             snnzmax, snnz, sacoo_k, sjcoo_k, sicoo_k, neval,nvecs,W,sigma,zeigv, ritzvec)
+      else
+         call arpack_sparse_coo_eigs(Num_wann,nnzmax,nnz,acoo,jcoo,icoo,neval,nvecs,W,sigma, zeigv, ritzvec)
+      endif
       call now(time3)
       eigv(1:neval, ik)= W(1:neval)
 
+
       !> get valley operator in coo format
       call valley_k_coo_sparsehr(nnzmax_valley, k3, acoo_valley, icoo_valley, jcoo_valley)
 
@@ -2246,7 +2313,7 @@ subroutine ek_bulk2D_spin
 
      nwann= Num_wann
      if (soc>0) nwann= Num_wann/2
-     print *, 'nwann', nwann
+    !print *, 'nwann', nwann
 
      if (SOC==0) stop 'you should set soc=0 in the input file'
      !> construct spin matrix

src/ham_bulk.f90

@@ -770,10 +770,9 @@ subroutine ham_bulk_coo_sparsehr_latticegauge(k,acoo,icoo,jcoo)
    integer :: i,j,ir
 
    do i=1,splen
-      ir=hirv(i)
       icoo(i)=hicoo(i)
       jcoo(i)=hjcoo(i)
-      posij=irvec(:, ir)
+      posij= hirv(1:3, i)
       kdotr=posij(1)*k(1)+posij(2)*k(2)+posij(3)*k(3)
       ratio= (cos(twopi*kdotr)+zi*sin(twopi*kdotr))
       acoo(i)=ratio*hacoo(i)
@@ -782,7 +781,6 @@ subroutine ham_bulk_coo_sparsehr_latticegauge(k,acoo,icoo,jcoo)
    return
 end subroutine ham_bulk_coo_sparsehr_latticegauge
 
-
 subroutine ham_bulk_coo_sparsehr(k,acoo,icoo,jcoo)
    !> This subroutine use sparse hr format
    !> Here we use atomic gauge which means the atomic position is taken into account
@@ -797,11 +795,10 @@ subroutine ham_bulk_coo_sparsehr(k,acoo,icoo,jcoo)
    complex(dp) ::  ratio
    integer :: i,j,ir
 
-   do i=1,splen
-      ir=hirv(i)
+   do i=1, splen
       icoo(i)=hicoo(i)
       jcoo(i)=hjcoo(i)
-      posij=irvec(:, ir)+ Origin_cell%wannier_centers_direct(:, jcoo(i))- Origin_cell%wannier_centers_direct(:, icoo(i))
+      posij= hirv(1:3, i)+ Origin_cell%wannier_centers_direct(:, jcoo(i))- Origin_cell%wannier_centers_direct(:, icoo(i))
       kdotr=posij(1)*k(1)+posij(2)*k(2)+posij(3)*k(3)
       ratio= (cos(twopi*kdotr)+zi*sin(twopi*kdotr))
       acoo(i)=ratio*hacoo(i)
@@ -810,6 +807,33 @@ subroutine ham_bulk_coo_sparsehr(k,acoo,icoo,jcoo)
    return
 end subroutine ham_bulk_coo_sparsehr
 
+
+subroutine overlap_bulk_coo_sparse(k, acoo, icoo, jcoo)
+   !> This subroutine use sparse hr format
+   !> Here we use atomic gauge which means the atomic position is taken into account
+   !> in the Fourier transformation
+   use para
+   implicit none
+
+   real(dp) :: k(3), posij(3)
+   real(dp) :: kdotr
+   integer,intent(inout) :: icoo(splen_overlap_input),jcoo(splen_overlap_input)
+   complex(dp),intent(inout) :: acoo(splen_overlap_input)
+   complex(dp) ::  ratio
+   integer :: i,j,ir
+
+   do i=1, splen_overlap_input
+      icoo(i)=sicoo(i)
+      jcoo(i)=sjcoo(i)
+      posij= sirv(1:3, i)+ Origin_cell%wannier_centers_direct(:, sjcoo(i))- Origin_cell%wannier_centers_direct(:, sicoo(i))
+      kdotr=posij(1)*k(1)+posij(2)*k(2)+posij(3)*k(3)
+      ratio= (cos(twopi*kdotr)+zi*sin(twopi*kdotr))
+      acoo(i)=ratio*sacoo(i)
+   end do
+
+   return
+end subroutine overlap_bulk_coo_sparse
+
 subroutine valley_k_coo_sparsehr(nnz, k,acoo,icoo,jcoo)
    !> This subroutine use sparse hr format
    !> Here we use atomic gauge which means the atomic position is taken into account

src/ham_slab.f90

@@ -105,10 +105,9 @@ subroutine ham_slab_sparseHR(nnzmax, k, acoo,jcoo,icoo)
 
          !> sum over R points to get H(k1, k2)
          do ims=1,splen
-            ir=hirv(ims)
-            ia=irvec(1,iR)
-            ib=irvec(2,iR)
-            ic=irvec(3,iR)
+            ia= hirv(1, ims)
+            ib= hirv(2, ims)
+            ic= hirv(3, ims)
 
             !> new lattice
             call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic)