fixed a bug in unfolding.f90 about the length of kpath

quanshengwu · Dec 9, 2019 · 3fd95cd · 3fd95cd
1 parent 278bf13
commit 3fd95cd
Show file tree

Hide file tree

Showing 4 changed files with 39 additions and 10 deletions.
diff --git a/examples/TBG-13.2degree/wt.in b/examples/TBG-13.2degree/wt.in
@@ -5,7 +5,7 @@
 
  !> bulk band structure calculation flag
  &CONTROL
- BulkBand_calc         = F
+ BulkBand_calc         = T
  BulkBand_plane_calc   = F
  BulkBand_Unfold_line_calc          = T
  BulkBand_Unfold_plane_calc         = F
@@ -21,8 +21,8 @@
  Eta_Arc = 0.050     ! infinite small value, like brodening
  E_arc = 1.00         ! energy contour plot
  OmegaNum_unfold =600       ! omega number
- OmegaMin = -10.0     ! energy interval
- OmegaMax =  8.0     ! energy interval
+ OmegaMin = -13.0     ! energy interval
+ OmegaMax =  7.0     ! energy interval
  Nk1 =200           ! number k points
  Nk2 =200          ! number k points
  /

diff --git a/soc/module.f90 b/soc/module.f90
@@ -714,6 +714,9 @@ module para
      real(dp), allocatable :: k3points_unfold_pointmode_cart(:, :)
      real(dp), allocatable :: k3points_unfold_pointmode_direct(:, :)
 
+     !> kpath for unfoled lattice
+     real(dp), allocatable :: k3len_folded(:)
+     real(dp), allocatable :: k3line_folded_stop(:)
 
      !> kpath for magnetic supercell
      real(dp),allocatable :: K3len_mag(:)  ! put all k points in a line in order to plot the bands 

diff --git a/soc/readinput.f90 b/soc/readinput.f90
@@ -1679,7 +1679,9 @@ subroutine readinput
    allocate(k3line_end(3, nk3lines))
    allocate(k3line_name(nk3lines+1))
    allocate(k3line_stop(nk3lines+1))
+   allocate(k3line_folded_stop(nk3lines+1))
    allocate(k3line_mag_stop(nk3lines+1))
+   k3line_folded_stop= 0d0
    k3line_mag_stop= 0d0
    k3line_stop= 0d0
    k3line_start= 0d0
@@ -1711,9 +1713,11 @@ subroutine readinput
    nk3_band= NN*nk3lines
    allocate(k3len(nk3_band))
    allocate(k3len_mag(nk3_band))
+   allocate(k3len_folded(nk3_band))
    allocate(k3points(3, nk3_band))
    k3len=0d0
    k3len_mag=0d0
+   k3len_folded=0d0
    k3points= 0d0
    t1= 0d0
    do j=1, nk3lines
@@ -1760,6 +1764,28 @@ subroutine readinput
       k3line_mag_stop(j+1)= t1
    enddo
 
+
+   !> for folded cell
+   t1=0
+   do j=1, nk3lines
+      do i=1, NN
+         kstart= k3line_start(:, j)
+         kend  = k3line_end(:, j)
+         k1= kstart(1)*Folded_cell%Kua+ kstart(2)*Folded_cell%Kub+ kstart(3)*Folded_cell%Kuc
+         k2= kend(1)*Folded_cell%Kua+ kend(2)*Folded_cell%Kub+ kend(3)*Folded_cell%Kuc
+
+         temp= dsqrt((k2(1)- k1(1))**2 &
+            +(k2(2)- k1(2))**2  &
+            +(k2(3)- k1(3))**2)/dble(NN-1)
+
+         if (i.gt.1) then
+            t1=t1+temp
+         endif
+         k3len_folded(i+(j-1)*NN)= t1
+      enddo
+      k3line_folded_stop(j+1)= t1
+   enddo
+
 104 continue
    if (.not.lfound .and. (BulkBand_line_calc.or.LandauLevel_k_calc)) then
       stop 'ERROR: please set KPATH_BULK for bulk band stOrigin_cell%Ructure calculation'

diff --git a/soc/unfolding.f90 b/soc/unfolding.f90
@@ -150,7 +150,7 @@ subroutine unfolding_kpath
       write(outfileindex, "('#column', i5, 3000i12)")(i, i=1, 2+NumberofSelectedOrbitals_groups*NumberofEta)
       do ik=1, nk3_band
          do ie=1, omeganum_unfold
-            write(outfileindex, '(300f12.6)')k3len(ik), omega(ie), &
+            write(outfileindex, '(300f12.6)')k3len_folded(ik), omega(ie), &
                ((spectrum_unfold_mpi(ie, ieta, ig, ik), ieta=1, NumberofEta), ig=1, NumberofSelectedOrbitals_groups)
          enddo
          write(outfileindex, *) ' '
@@ -183,22 +183,22 @@ subroutine unfolding_kpath
       write(outfileindex, '(a)')'#set ylabel offset 1.5,0'
       write(outfileindex, '(a,f12.6)')'emin=', omegamin
       write(outfileindex, '(a,f12.6)')'emax=', omegamax
-      write(outfileindex, '(a, f10.5, a)')'set xrange [0: ', maxval(k3len), ']'
+      write(outfileindex, '(a, f10.5, a)')'set xrange [0: ', maxval(k3len_folded), ']'
       if (index(Particle,'phonon')/=0) then
          write(outfileindex, '(a, f10.5, a)')'set yrange [0: emax ]'
          write(outfileindex, '(a)')'set ylabel "Frequency (THz)"'
       else
          write(outfileindex, '(a)')'set ylabel "Energy (eV)"'
          write(outfileindex, '(a)')'set yrange [ emin : emax ]'
       endif
-      write(outfileindex, 202, advance="no") (k3line_name(i), k3line_stop(i), i=1, nk3lines)
-      write(outfileindex, 203)k3line_name(nk3lines+1), k3line_stop(nk3lines+1)
+      write(outfileindex, 202, advance="no") (k3line_name(i), k3line_folded_stop(i), i=1, nk3lines)
+      write(outfileindex, 203)k3line_name(nk3lines+1), k3line_folded_stop(nk3lines+1)
 
       do i=1, nk3lines-1
          if (index(Particle,'phonon')/=0) then
-            write(outfileindex, 204)k3line_stop(i+1), '0.0', k3line_stop(i+1), 'emax'
+            write(outfileindex, 204)k3line_folded_stop(i+1), '0.0', k3line_folded_stop(i+1), 'emax'
          else
-            write(outfileindex, 204)k3line_stop(i+1), 'emin', k3line_stop(i+1), 'emax'
+            write(outfileindex, 204)k3line_folded_stop(i+1), 'emin', k3line_folded_stop(i+1), 'emax'
          endif
       enddo
 
@@ -289,7 +289,7 @@ subroutine unfolding_kplane
          ik =ik +1
          kxy(:, ik)= K3D_start+ K3D_vec1*(i-1)/dble(Nk1-1)+ K3D_vec2*(j-1)/dble(Nk2-1) &
             -(K3D_vec1+K3D_vec2)/2d0
-         kxy_shape(:, ik)= kxy(1, ik)* Origin_cell%Kua+ kxy(2, ik)* Origin_cell%Kub+ kxy(3, ik)* Origin_cell%Kuc 
+         kxy_shape(:, ik)= kxy(1, ik)* Folded_cell%Kua+ kxy(2, ik)* Folded_cell%Kub+ kxy(3, ik)* Folded_cell%Kuc 
          call rotate_k3_to_kplane(kxy_shape(:, ik), kxy_plane(:, ik))
       enddo
    enddo