diff --git a/bin/wt2.x b/bin/wt2.x new file mode 100755 index 00000000..7c9e4158 Binary files /dev/null and b/bin/wt2.x differ diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 00000000..8e37dd70 --- /dev/null +++ b/src/Makefile @@ -0,0 +1,47 @@ +OBJ = module.o sparse.o wt_aux.o math_lib.o symmetry.o readHmnR.o inverse.o proteus.o \ + eigen.o ham_qlayer2qlayer.o psi.o unfolding.o rand.o \ + ham_slab.o ham_bulk.o ek_slab.o ek_bulk_polar.o ek_bulk.o \ + readinput.o fermisurface.o surfgreen.o surfstat.o \ + mat_mul.o ham_ribbon.o ek_ribbon.o \ + fermiarc.o berrycurvature.o \ + wanniercenter.o dos.o orbital_momenta.o \ + landau_level_sparse.o landau_level.o lanczos_sparse.o \ + berry.o wanniercenter_adaptive.o \ + effective_mass.o findnodes.o \ + sigma_OHE.o sigma.o Boltz_transport_anomalous.o \ + main.o + +# compiler +F90 = mpiifort -fpp -DMPI -DINTELMKL -fpe3 + +INCLUDE = -I${MKLROOT}/include +WFLAG = -nogen-interface +OFLAG = -O3 -static-intel +FFLAG = $(OFLAG) $(WFLAG) +LFLAG = $(OFLAG) + +# ARPACK LIBRARY +ARPACK=/home/yanyi2024/0-tools/0-Basic/ARPACK/libarpack_MAC.a + +# blas and lapack libraries +# static linking +#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 = ${ARPACK} -L/${MKLROOT}/lib/intel64 -lmkl_core -lmkl_sequential -lmkl_intel_lp64 -lpthread + +main : $(OBJ) + $(F90) $(LFLAG) $(OBJ) -o wt2.x $(LIBS) + cp -f wt2.x ../bin + +.SUFFIXES: .o .f90 + +.f90.o : + $(F90) $(FFLAG) $(INCLUDE) -c $*.f90 + +clean : + rm -f *.o *.mod *~ wt2.x + diff --git a/src/ham_qlayer2qlayer.f90 b/src/ham_qlayer2qlayer.f90 index 4f0559f4..41e94b2d 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 @@ -443,7 +442,7 @@ subroutine ham_qlayer2qlayer2(k,Hij) implicit none ! loop index - integer :: iR, inew_ic + integer :: iR, inew_ic, int_ic real(dp) :: ia, ib, ic ! new index used to sign irvec @@ -460,22 +459,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 28597d3e..740d7754 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 70ed2f71..5bd990e9 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,109 +359,142 @@ subroutine get_hmnr_cell(cell) allocate(allirs(nrpts, num_wann, num_wann, 3)) - call cart_direct_real(shift_to_topsurface_cart, shift_vec_direct, cell%lattice) + 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) !> Get new Hrs rpts_array=0 - nrpts_new=0 + nrpts_surfacecell=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 Problem: Just cut the nrpts_surfacecell, but the irvec is not well ordered. The order (iter) can be anywhere. !> write to new_hr.dat - nrpts_new=sum(rpts_array)-iter + ! nrpts_surfacecell=sum(rpts_array)-iter + + nrpts_surfacecell=sum(rpts_array) outfileindex= outfileindex+ 1 open(unit=outfileindex, file='wannier90_hr_newcell.dat') write(outfileindex, '(a,1X,a,1X,a,1X, a, a)') & 'HmnR for new cell generated at ', time_now, date_now, 'UTC', zone_now write(outfileindex, *)Num_wann - 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 + write(outfileindex, *)nrpts_surfacecell + write(outfileindex, '(15I5)')(1 , i=1, nrpts_surfacecell) + do ir=1, nrpts_surfacecell + max_val=maxval(abs(HmnR_surfacecell(:, :, ir)))/eV2Hartree !if (max_val