diff --git a/bin/wt2.x b/bin/wt2.x index 4e8ce0b..179745a 100755 Binary files a/bin/wt2.x and b/bin/wt2.x differ diff --git a/src/ham_qlayer2qlayer.f90 b/src/ham_qlayer2qlayer.f90 index 4f0559f..f1d990d 100644 --- a/src/ham_qlayer2qlayer.f90 +++ b/src/ham_qlayer2qlayer.f90 @@ -17,6 +17,7 @@ subroutine ham_qlayer2qlayer(k,H00new,H01new) ! new index used to sign irvec real(dp) :: new_ia,new_ib,new_ic integer :: inew_ic + integer :: int_ic ! wave vector k times lattice vector R real(Dp) :: kdotr @@ -44,21 +45,19 @@ subroutine ham_qlayer2qlayer(k,H00new,H01new) allocate(Hij(Num_wann,Num_wann,-ijmax:ijmax)) Hij=0.0d0 - do iR=1,Nrpts - ia=irvec(1,iR) - ib=irvec(2,iR) - ic=irvec(3,iR) - - call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic) - - inew_ic= int(new_ic) - if (abs(new_ic).le.ijmax)then - kdotr=k(1)*new_ia+ k(2)*new_ib + do iR=1,nrpts_surfacecell + ia=irvec_surfacecell(1,iR) + ib=irvec_surfacecell(2,iR) + ic=irvec_surfacecell(3,iR) + + int_ic= int(ic) + if (abs(ic).le.ijmax)then + kdotr=k(1)*ia+ k(2)*ib ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr) - Hij(1:Num_wann, 1:Num_wann, inew_ic )& - =Hij(1:Num_wann, 1:Num_wann, inew_ic)& - +HmnR(:,:,iR)*ratio/ndegen(iR) + Hij(1:Num_wann, 1:Num_wann, int_ic )& + =Hij(1:Num_wann, 1:Num_wann, int_ic)& + +HmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(iR) endif enddo @@ -121,6 +120,7 @@ subroutine ham_qlayer2qlayer_bak(k,H00new,H01new) ! new index used to sign irvec real(dp) :: new_ia,new_ib,new_ic integer :: inew_ic + integer :: int_ic ! wave vector k times lattice vector R real(Dp) :: kdotr @@ -148,21 +148,19 @@ subroutine ham_qlayer2qlayer_bak(k,H00new,H01new) allocate(Hij(-ijmax:ijmax,Num_wann,Num_wann)) Hij=0.0d0 - do iR=1,Nrpts - ia=irvec(1,iR) - ib=irvec(2,iR) - ic=irvec(3,iR) - - call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic) - - inew_ic= int(new_ic) - if (abs(new_ic).le.ijmax)then - kdotr=k(1)*new_ia+ k(2)*new_ib + do iR=1,nrpts_surfacecell + ia=irvec_surfacecell(1,iR) + ib=irvec_surfacecell(2,iR) + ic=irvec_surfacecell(3,iR) + + int_ic= int(ic) + if (abs(ic).le.ijmax)then + kdotr=k(1)*ia+ k(2)*ib ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr) - Hij(inew_ic, 1:Num_wann, 1:Num_wann )& - =Hij(inew_ic, 1:Num_wann, 1:Num_wann )& - +HmnR(:,:,iR)*ratio/ndegen(iR) + Hij(1:Num_wann, 1:Num_wann, int_ic )& + =Hij(1:Num_wann, 1:Num_wann, int_ic)& + +HmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(iR) endif enddo @@ -441,13 +439,13 @@ subroutine ham_qlayer2qlayer2(k,Hij) use para implicit none - + ! loop index - integer :: iR, inew_ic + integer :: iR, int_ic real(dp) :: ia, ib, ic - ! new index used to sign irvec - real(dp) :: new_ia,new_ib,new_ic + ! ! new index used to sign irvec + ! real(dp) :: new_ia,new_ib,new_ic ! wave vector k times lattice vector R real(Dp) :: kdotr @@ -460,22 +458,21 @@ subroutine ham_qlayer2qlayer2(k,Hij) 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) - ic=irvec(3,iR) - !> new lattice - call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic) + do iR=1,nrpts_surfacecell + ia=irvec_surfacecell(1,iR) + ib=irvec_surfacecell(2,iR) + ic=irvec_surfacecell(3,iR) - inew_ic= int(new_ic) - if (abs(new_ic).le.ijmax)then - kdotr=k(1)*new_ia+k(2)*new_ib + int_ic= int(ic) + if (abs(ic).le.ijmax)then + kdotr=k(1)*ia+k(2)*ib ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr) - Hij(inew_ic, 1:Num_wann, 1:Num_wann )& - =Hij(inew_ic, 1:Num_wann, 1:Num_wann )& - +HmnR(:,:,iR)*ratio/ndegen(iR) + Hij(int_ic, 1:Num_wann, 1:Num_wann )& + =Hij(int_ic, 1:Num_wann, 1:Num_wann )& + +HmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(iR) + endif enddo diff --git a/src/module.f90 b/src/module.f90 index 28597d3..342f5a3 100644 --- a/src/module.f90 +++ b/src/module.f90 @@ -897,11 +897,18 @@ module para integer, allocatable :: ndegen(:) ! degree of degeneracy of R point - complex(dp), allocatable :: HmnR_newcell(:,:,:) ! Hamiltonian m,n are band indexes + complex(dp), allocatable :: HmnR_surfacecell(:,:,:) ! Hamiltonian m,n are band indexes real(dp), allocatable :: Atom_position_cart_newcell(:,:) ! Hamiltonian m,n are band indexes real(dp), allocatable :: Atom_position_direct_newcell(:,:) ! Hamiltonian m,n are band indexes - integer, allocatable :: irvec_newcell(:,:) ! R coordinates - integer, allocatable :: ndegen_newcell(:) ! degree of degeneracy of R point + integer, allocatable :: irvec_surfacecell(:,:) ! R coordinates + integer, allocatable :: ndegen_surfacecell(:) ! degree of degeneracy of R point + + real(dp), allocatable :: Rmn_old(:) + real(dp), allocatable :: Rmn_new(:) + real(dp), allocatable :: irvec_new(:) + integer, allocatable :: irvec_new_int(:) + integer, allocatable :: nrpts_surfacecell + real(dp),public, save :: Rua_newcell(3) !> three rotated primitive vectors in old coordinate system real(dp),public, save :: Rub_newcell(3) !> three rotated primitive vectors in old coordinate system real(dp),public, save :: Ruc_newcell(3) !> three rotated primitive vectors in old coordinate system diff --git a/src/readHmnR.f90 b/src/readHmnR.f90 index 70ed2f7..ff6a53c 100644 --- a/src/readHmnR.f90 +++ b/src/readHmnR.f90 @@ -337,7 +337,7 @@ subroutine get_hmnr_cell(cell) real(dp) :: shift_vec_direct(3) !> for newcell - real(dp) :: apos1d(3),apos2d(3) + real(dp) :: apos1d(3),apos2d(3),apos1dprime(3),apos2dprime(3) !>count newcell nrpts integer :: max_ir integer :: nir1,nir2,nir3, ir_cell @@ -349,6 +349,9 @@ subroutine get_hmnr_cell(cell) integer, allocatable :: allirs(:, :, :, :) integer :: irn1(3),irn2(3) + !> The Allocate Process + integer :: ia1,ia2,ia1prime,ia2prime + max_ir=8 nrpts_max=(2*max_ir+1)**3 allocate( rpts_array(-max_ir:max_ir,-max_ir:max_ir,-max_ir:max_ir)) @@ -356,6 +359,11 @@ subroutine get_hmnr_cell(cell) allocate(allirs(nrpts, num_wann, num_wann, 3)) + + allocate(Rmn_old(3)) + allocate(Rmn_new(3)) + allocate(irvec_new_int(3)) + call cart_direct_real(shift_to_topsurface_cart, shift_vec_direct, cell%lattice) call date_and_time(DATE=date_now,ZONE=zone_now, TIME=time_now) @@ -364,83 +372,110 @@ subroutine get_hmnr_cell(cell) nrpts_new=0 rpts_map= 0 - !> get number of R points for the new cell first + !> 1. Get number of R points for the new cell do ir=1, nrpts - do i=1, Num_wann - do j=1, Num_wann - apos1d= Origin_cell%wannier_centers_direct(:, i)- shift_vec_direct - apos2d= Origin_cell%wannier_centers_direct(:, j)- shift_vec_direct+irvec(:, ir) - call latticetransform(apos1d(1),apos1d(2),apos1d(3),new_ia,new_ib,new_ic) - irn1=floor([new_ia,new_ib,new_ic]) - - call latticetransform(apos2d(1),apos2d(2),apos2d(3),new_ia,new_ib,new_ic) - irn2=floor([new_ia,new_ib,new_ic]) - - nir1=irn2(1)-irn1(1) - nir2=irn2(2)-irn1(2) - nir3=irn2(3)-irn1(3) - if (abs(nir1)>max_ir .or. abs(nir2)>max_ir .or. abs(nir3)>max_ir) cycle - rpts_array(nir1,nir2,nir3)=1 + do j=1, Num_wann + do i=1, Num_wann + ia1 = Origin_cell%spinorbital_to_atom_index(i) + ia2 = Origin_cell%spinorbital_to_atom_index(j) + apos1d = Origin_cell%Atom_position_direct(:, ia1) + apos2d = Origin_cell%Atom_position_direct(:, ia2) + + ia1prime = Cell_defined_by_surface%spinorbital_to_atom_index(i) + ia2prime = Cell_defined_by_surface%spinorbital_to_atom_index(j) + apos1dprime = Cell_defined_by_surface%Atom_position_direct(:,ia1prime) + apos2dprime = Cell_defined_by_surface%Atom_position_direct(:,ia2prime) + + Rmn_old = irvec(:, ir) + apos2d - apos1d + call latticetransform(Rmn_old(1),Rmn_old(2),Rmn_old(3),Rmn_new(1),Rmn_new(2),Rmn_new(3)) + irvec_new = Rmn_new - (apos2dprime - apos1dprime) + + !> Due to the accuracy of computing, need to rounding (but always tiny) + irvec_new_int(1) = ANINT(irvec_new(1)) + irvec_new_int(2) = ANINT(irvec_new(2)) + irvec_new_int(3) = ANINT(irvec_new(3)) + + if (abs(irvec_new_int(1))>max_ir .or. abs(irvec_new_int(2))>max_ir .or. abs(irvec_new_int(3))>max_ir) cycle + rpts_array(irvec_new_int(1),irvec_new_int(2),irvec_new_int(3))=1 + enddo enddo enddo - !> find all irvec - nrpts_new= sum(rpts_array) - allocate(irvec_newcell(3, Nrpts_new)) - iter= 0 + !> The total number of lattice points searched above + nrpts_surfacecell= sum(rpts_array) + + allocate(irvec_surfacecell(3, nrpts_surfacecell)) + irvec_surfacecell = 0 + + !> 2. Create an order map to sign the R points generated in step1 + + iter = 0 do nir3=-max_ir,max_ir do nir2=-max_ir,max_ir do nir1=-max_ir,max_ir if (rpts_array(nir1, nir2, nir3)==1) then iter=iter+1 - irvec_newcell(:, iter)=[nir1, nir2, nir3] + irvec_surfacecell(:, iter)=[nir1, nir2, nir3] rpts_map(nir1, nir2, nir3)=iter endif enddo enddo enddo - !>> hamiltonian for the new cell - allocate(HmnR_newcell(Num_wann, Num_wann, Nrpts_new)) - allocate(ndegen_newcell(Nrpts_new)) - HmnR_newcell= 0d0 - ndegen_newcell= 1 - !> get number of R points for the new cell first + allocate(HmnR_surfacecell(Num_wann, Num_wann, nrpts_surfacecell)) + allocate(ndegen_surfacecell(nrpts_surfacecell)) + HmnR_surfacecell= 0d0 + ndegen_surfacecell= 1 + + !> 3. Allocate the old HmnR to the new HmnR. + !> Note: The cell can't be treated as a mass point. We need to consider the relative coordinates of atoms. do ir=1, nrpts do j=1, Num_wann do i=1, Num_wann - !ia1=Origin_cell%spinorbital_to_atom_index(i) - !ia2=Origin_cell%spinorbital_to_atom_index(j) - !apos1d= Origin_cell%Atom_position_direct(:, ia1)- shift_vec_direct - !apos2d= Origin_cell%Atom_position_direct(:, ia2)- shift_vec_direct+irvec(:, ir) - apos1d= Origin_cell%wannier_centers_direct(:, i)- shift_vec_direct - apos2d= Origin_cell%wannier_centers_direct(:, j)- shift_vec_direct+irvec(:, ir) - call latticetransform(apos1d(1),apos1d(2),apos1d(3),new_ia,new_ib,new_ic) - irn1=floor([new_ia,new_ib,new_ic]) - - call latticetransform(apos2d(1),apos2d(2),apos2d(3),new_ia,new_ib,new_ic) - irn2=floor([new_ia,new_ib,new_ic]) - - nir1=irn2(1)-irn1(1) - nir2=irn2(2)-irn1(2) - nir3=irn2(3)-irn1(3) - if (abs(nir1)>max_ir .or. abs(nir2)>max_ir .or. abs(nir3)>max_ir) cycle - ir_cell= rpts_map(nir1, nir2, nir3) - HmnR_newcell(i,j,ir_cell)=HmnR(i,j,ir)/ndegen(ir) + + ia1 = Origin_cell%spinorbital_to_atom_index(i) + ia2 = Origin_cell%spinorbital_to_atom_index(j) + apos1d = Origin_cell%Atom_position_direct(:, ia1) + apos2d = Origin_cell%Atom_position_direct(:, ia2) + + ia1prime = Cell_defined_by_surface%spinorbital_to_atom_index(i) + ia2prime = Cell_defined_by_surface%spinorbital_to_atom_index(j) + apos1dprime = Cell_defined_by_surface%Atom_position_direct(:,ia1prime) + apos2dprime = Cell_defined_by_surface%Atom_position_direct(:,ia2prime) + + !> R'_mn = R' + tau'_2 - tau'_1 + !> R_mn = R + tau_2 - tau_1 + !> R'_mn = Pinv * R_mn + !> R' = (Pinv * R_mn) - (tau'_2 - tau'_1) + + Rmn_old = irvec(:, ir) + apos2d - apos1d + call latticetransform(Rmn_old(1),Rmn_old(2),Rmn_old(3),Rmn_new(1),Rmn_new(2),Rmn_new(3)) + irvec_new = Rmn_new - (apos2dprime - apos1dprime) + + !> For safety, we perform a rounding operation due to the accuracy of computing + irvec_new_int(1) = ANINT(irvec_new(1)) + irvec_new_int(2) = ANINT(irvec_new(2)) + irvec_new_int(3) = ANINT(irvec_new(3)) + + if (abs(irvec_new_int(1))>max_ir .or. abs(irvec_new_int(2))>max_ir .or. abs(irvec_new_int(3))>max_ir) cycle + ir_cell = rpts_map(irvec_new_int(1), irvec_new_int(2), irvec_new_int(3)) + HmnR_surfacecell(i,j,ir_cell) = HmnR(i,j,ir)/ndegen(ir) + enddo enddo enddo - !> do cut-off according to the hopping value - iter= 0 - do ir=1, nrpts_new - max_val=maxval(abs(HmnR_newcell(:, :, ir)))/eV2Hartree - if (max_val do cut-off according to the hopping value + ! iter= 0 + ! do ir=1, nrpts_new + ! max_val=maxval(abs(HmnR_newcell(:, :, ir)))/eV2Hartree + ! if (max_val write to new_hr.dat @@ -453,12 +488,12 @@ subroutine get_hmnr_cell(cell) write(outfileindex, *)nrpts_new write(outfileindex, '(15I5)')(1 , i=1, nrpts_new) do ir=1, nrpts_new - max_val=maxval(abs(HmnR_newcell(:, :, ir)))/eV2Hartree + max_val=maxval(abs(HmnR_surfacecell(:, :, ir)))/eV2Hartree !if (max_val