整理用Fortran读取Wannier90拟合的hr数据，后续设计到极化率以及响应之类的计算还是用Fortran的计算速度是最快的。
前言
在之前博客利用Wannier90的hr数据构建动量空间能带并计算高对称点宇称中用Python实现过读取Wannier90_hr.dat这个数据，当时只是涉及到计算一下能带以及拓扑不变量，此时不会遇到一些计算瓶颈。在涉及到一些响应系数或者RPA极化率计算 的时候Python用起来就有点慢了，所以干脆也就用Fortran写一个读取Wannier90_hr.dat的程序，方面以后自己使用。
代码
module param
implicit none
integer,parameter::dp = kind(1.0)
integer numk_bz,num_wann,nrpts,numk_FS,kn
real(dp) mu
parameter(kn = 100,mu = 0.0)
real(dp),parameter::pi = 3.1415926535897
real(dp),parameter::deltaE = 0.01 ! 费米面展宽
complex(dp),parameter::im = (0.,1.) !Imagine unit
real(dp),allocatable::ones_ham(:,:) ! 为了后续考虑化学势
complex(dp),allocatable::HmnR(:,:,:) ! hr中的hopping
integer,allocatable::irvec(:,:) ! hr中的位置矢量
integer,allocatable::indorbit(:,:,:,:)
real(dp),allocatable::BZklist(:,:)
real(dp),allocatable::FS_points(:,:)
integer FS_index(-kn:kn - 1,-kn:kn - 1)
!------------------------------------------
integer,allocatable::FSindex(:,:) ! 费米点索引指标
end module
!==============================================================================================================================================================
program main
use param
use mpi
implicit none
integer numcore,indcore,ierr
!-------------------------------------------------
real(dp) code_time_start,code_time_end,code_time
integer nki,nkf
character(len = 20)::filename,char1,char2,char3,char4
!-------------------------------------------------------------------
call MPI_INIT(ierr) ! 初始化进程
call MPI_COMM_RANK(MPI_COMM_WORLD, indcore, ierr) ! 得到本进程在通信空间中的rank值,即在组中的逻辑编号(该indcore为0到numcore-1间的整数,相当于进程的ID。)
call MPI_COMM_SIZE(MPI_COMM_WORLD, numcore, ierr) !获得进程数量,用numcoer保存
! 循环分拆
nki = floor(indcore * (2.0 * kn)/numcore) - kn
nkf = floor((indcore + 1) * (2.0 * kn)/numcore) - kn - 1
if(indcore.eq.0)then
call CPU_TIME(code_time_start)
write(*,*)"Start calcation: ",code_time_start
write(*,*)"Number of CPU = ",numcore
write(*,*)"Number of nk = ",2 * kn
end if
! 输出程序的基本参数信息
if(indcore.eq.0)then
write(*,*)"Number of Fermi points = ",numk_FS
write(*,*)"Number of BZ points = ",numk_bz
call squareBZ()
call read_hr()
call fermisurface() ! 确定体系的费米面,其中会给出所有费米点的位置
call CPU_TIME(code_time_end)
write(*,*)"Calcation finished",code_time_end
write(*,*)"Code execution time:",abs(code_time_end - code_time_start),"second"
end if
call MPI_Finalize(ierr)
stop
end program main
!================================================================================================================================================================================================
subroutine fermisurface()
use param
integer ibz,ie,ikx,iky
real(dp) kx,ky,mateigval(num_wann)
complex(dp) Ham_temp(num_wann,num_wann),matvec(num_wann,num_wann)
! call honeycomyBZ()
open(13,file = "FS.dat")
! numk_FS = 0 ! 全局变量,共计费米点的数量
do ibz = 1,numk_bz
kx = BZklist(1,ibz)
ky = BZklist(2,ibz)
call HamR_to_K(kx,ky,0.0,Ham_temp)
call diagonalize_Hermitian_matrix(num_wann,Ham_temp,matvec,mateigval)
do ie = 1,num_wann
if (abs(mateigval(ie)) < deltaE)then
write(13,"(3F12.6)")kx,ky,1.0
numk_FS = numk_FS + 1 ! 统计费米点的数量
end if
end do
end do
close(13)
!--------------------------------------------------------------
!确定费米面点的数量后就可以确定相互作用矩阵的维度
if(numk_FS.eq.0)then
write(*,*)"The Number of Fermi surface points = ",numk_FS
write(*,*)"----------------------------------------------"
write(*,*)" Warnning "
write(*,*)"----------------------------------------------"
stop
else
allocate(FS_points(2,numk_FS))
FS_points = 0.0
end if
!-----------------------------------------------
! 根据动量坐标记录这是第几个费米点
numk_FS = 0
do iky = -kn,kn - 1
ky = 2.0 * pi * iky/kn
do ikx = -kn,kn - 1
kx = 2.0 * pi * ikx/kn
call HamR_to_K(kx,ky,0.0,Ham_temp)
call diagonalize_Hermitian_matrix(num_wann,Ham_temp,matvec,mateigval)
do ie = 1,num_wann
if (abs(mateigval(ie)) < deltaE)then
numk_FS = numk_FS + 1
FS_index(ikx,iky) = numk_FS
FS_points(1,numk_FS) = kx
FS_points(2,numk_FS) = ky
end if
end do
end do
end do
return
end subroutine
!================================================================================================================================================================================================
subroutine squareBZ()
use param
integer ikx,iky,i0
! 对于四方点阵,BZ的点数可以直接确定
numk_bz = (2 * kn)**2
allocate(BZklist(2,numk_bz))
i0 = 0
do ikx = -kn,kn - 1
do iky = -kn,kn - 1
i0 = i0 + 1
BZklist(1,i0) = 2.0 * pi * ikx/kn
BZklist(2,i0) = 2.0 * pi * iky/kn
end do
end do
return
end subroutine
!==============================================================================================================================================================
subroutine HamR_to_K(kx,ky,kz,Ham_temp)
! 将hr变成动量空间
use param
implicit none
real(dp) kx,ky,kz
integer i1,i2,i3,w1,w2,r1,r2,r3
complex(dp) Ham_temp(num_wann,num_wann)
ones_ham = 0.0 ! 对于具有allocatable属性的变量,一定要先赋值为零
! 设置单位矩阵
do i1 = 1,num_wann
ones_ham(i1,i1) = 1.0
end do
Ham_temp = 0.0
do i1 = 1,nrpts
do i2 = 1,num_wann
do i3 = 1,num_wann
r1 = irvec(1,i1)
r2 = irvec(2,i1)
r3 = irvec(3,i1)
w1 = indorbit(1,i1,i2,i3)
w2 = indorbit(2,i1,i2,i3)
Ham_temp(w1,w2) = Ham_temp(w1,w2) + (cos(kx * r1 + ky * r2 + kz * r3) + im * sin(kx * r1 + ky * r2 + kz * r3)) * HmnR(w1,w2,i1)
end do
end do
end do
Ham_temp = Ham_temp - mu * ones_ham ! 在Wannier90_hr的基础上考虑体系的化学势
end subroutine
!==============================================================================================================================================================
subroutine read_hr()
! 读取DFT的hr数据,并返回Hr以及
use param
implicit none
integer :: i,j,k
integer :: r1,r2,r3,w1,w2
real(dp) :: hr,hi
real(dp),allocatable::ndegen(:)
logical stat
inquire(file = 'wannier90_hr.dat', exist = stat)
if(stat)then
open(12,file = 'wannier90_hr.dat',action = 'read')
read(12,*)
read(12,*) num_wann ! 轨道数量,也决定了动量空间哈密顿量维度,这是一个全局变量
! read(12,*)r1
read(12,*) nrpts
allocate(ndegen(nrpts))
! 全局变量
allocate(HmnR(num_wann,num_wann,nrpts),irvec(3,nrpts))
allocate(indorbit(2,nrpts,num_wann,num_wann))
read(12,*) (ndegen(i),i = 1,nrpts)
do i = 1,nrpts
do j = 1,num_wann
do k = 1,num_wann
read(12,*) r1,r2,r3,w1,w2,hr,hi
irvec(1,i) = r1
irvec(2,i) = r2
irvec(3,i) = r3
indorbit(1,i,j,k) = w1 ! 哈密顿量中的轨道索引,后续构建动量空间哈密顿量使用
indorbit(2,i,j,k) = w2
HmnR(w1,w2,i) = hr + im * hi
enddo
end do
end do
close(12)
allocate(ones_ham(num_wann,num_wann))
else
write(*,*)"**************** Error *********************"
write(*,*)"Tight-binding dat is not exist"
write(*,*)"**************** Error *********************"
end if
return
end subroutine read_hr
!==============================================================================================================================================================
function delta(x)
!> Lorentz or gaussian expansion of the Delta function
use param, only : dp, pi
implicit none
! real(dp), intent(in) :: eta
real(dp), intent(in) :: x
real(dp) :: delta, y
real(dp) :: eta = 0.001
!> Lorentz expansion
!delta= 1d0/pi*eta/(eta*eta+x*x)
y= x*x/eta/eta/2d0
!> Gaussian broadening
!> exp(-60) = 8.75651076269652e-27
if (y>60d0) then
delta = 0d0
else
delta= exp(-y)/sqrt(2d0*pi)/eta
endif
return
end function delta
!==============================================================================================================================================================
subroutine diagonalize_complex_matrix(matdim,matin,matout,mateigval)
! 对角化一般复数矩阵,这里的本征值是个复数
! matin 输入矩阵 matout 本征矢量 mateigval 本征值
implicit none
integer,parameter::dp = kind(1.0)
integer matdim,LDA,LDVL,LDVR,LWMAX,INFO,LWORK
complex(dp),intent(in)::matin(matdim,matdim)
complex(dp),intent(out)::matout(matdim,matdim)
complex(dp),intent(out)::mateigval(matdim)
real(dp),allocatable::RWORK(:)
complex(dp),allocatable::WORK(:)
complex(dp),allocatable::VL(:,:)
complex(dp),allocatable::VR(:,:)
LDA = matdim
LDVL = matdim
LDVR = matdim
LWMAX = 2 * matdim + matdim**2
! write(*,*)matin
allocate(RWORK(2 * matdim))
allocate(VL(LDVL,matdim))
allocate(VR(LDVR, matdim))
allocate(WORK(LWMAX))
matout = matin
LWORK = -1
call cgeev( 'V', 'N', matdim, matout, LDA, mateigval, VL, LDVL,VR, LDVR, WORK, LWORK, RWORK, INFO)
LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
call cgeev( 'V', 'N', matdim, matout, LDA, mateigval, VL, LDVL,VR, LDVR, WORK, LWORK, RWORK, INFO )
IF( INFO.GT.0 ) THEN
WRITE(*,*)'The algorithm failed to compute eigenvalues.'
! STOP
END IF
end subroutine diagonalize_complex_matrix
!================================================================================================================================================================================================
subroutine Matrix_Inv(matdim,matin,matout)
! 矩阵求逆
implicit none
integer,parameter::dp = kind(1.0)
integer matdim,info
complex(dp),intent(in) :: matin(matdim,matdim)
complex(dp):: matout(size(matin,1),size(matin,2))
real(dp):: work2(size(matin,1)) ! work2 array for LAPACK
integer::ipiv(size(matin,1)) ! pivot indices
! Store matin in matout to prevent it from being overwritten by LAPACK
matout = matin
! SGETRF computes an LU factorization of a general M - by - N matrix A
! using partial pivoting with row interchanges .
call CGETRF(matdim,matdim,matout,matdim,ipiv,info)
! if (info.ne.0) stop 'Matrix is numerically singular!'
if (info.ne.0) write(*,*)'Matrix is numerically singular!'
! SGETRI computes the inverse of a matrix using the LU factorization
! computed by SGETRF.
call CGETRI(matdim,matout,matdim,ipiv,work2,matdim,info)
! if (info.ne.0) stop 'Matrix inversion failed!'
if (info.ne.0) write(*,*)'Matrix inversion failed!'
return
end subroutine Matrix_Inv
!================================================================================================================================================================================================
subroutine diagonalize_Hermitian_matrix(matdim,matin,matout,mateigval)
! 厄米矩阵对角化
! matin 输入矩阵 matout 本征矢量 mateigval 本征值
implicit none
integer,parameter::dp = kind(1.0)
integer matdim
integer lda0,lwmax0,lwork,lrwork,liwork,info
complex(dp) matin(matdim,matdim),matout(matdim,matdim)
real(dp) mateigval(matdim)
complex(dp),allocatable::work(:)
real(dp),allocatable::rwork(:)
integer,allocatable::iwork(:)
!-----------------
lda0 = matdim
lwmax0 = 2 * matdim + matdim**2
allocate(work(lwmax0))
allocate(rwork(1 + 5 * matdim + 2 * matdim**2))
allocate(iwork(3 + 5 * matdim))
matout = matin
lwork = -1
liwork = -1
lrwork = -1
call cheevd('V','U',matdim,matout,lda0,mateigval,work,lwork ,rwork,lrwork,iwork,liwork,info)
lwork = min(2 * matdim + matdim**2, int( work( 1 ) ) )
lrwork = min(1 + 5 * matdim + 2 * matdim**2, int( rwork( 1 ) ) )
liwork = min(3 + 5 * matdim, iwork( 1 ) )
call cheevd('V','U',matdim,matout,lda0,mateigval,work,lwork,rwork,lrwork,iwork,liwork,info)
if( info .GT. 0 ) then
open(11,file = "mes.dat",status = "unknown")
write(11,*)'The algorithm failed to compute eigenvalues.'
close(11)
end if
return
end subroutine diagonalize_Hermitian_matrix
绘图程序
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.cm as cm
import os
import re
plt.rc('font', family='Times New Roman')
config = {
"font.size": 30,
"mathtext.fontset":'stix',
"font.serif": ['SimSun'],
}
rcParams.update(config) # Latex 字体设置
#----------------------------------------------------------------------------
def WannierDETband():
path = os.path.abspath('.')
band_file = os.path.join(path,'BAND.dat')
knode_file = os.path.join(path,'KLABELS')
fermi_file = os.path.join(path,'FERMI_ENERGY')
plt.figure(figsize = (10,8))
#******************** DFT band ************************************************************
band_data = np.loadtxt(band_file)
plt.plot(band_data[:,0],band_data[:,1:],color = 'b',linewidth = 2.0,label = "DFT")
xmin = np.min(band_data[:,0])
xmax = np.max(band_data[:,0])
ymin = np.min(band_data[:,1:])
ymax = np.max(band_data[:,1:])
plt.xlim(xmin,xmax)
# plt.ylim(ymin,ymax)
plt.ylim(-6,6)
with open(knode_file,"r",encoding = "utf-8") as f:
lines = f.readlines()
lines = lines[1:-3]
knodes = []
for i in range(len(lines)):
knodes.append(str.split(lines[i]))
knodes[i][1] = float(knodes[i][1])
plt.axvline(x = knodes[i][1],linewidth = 1.5,color = 'silver')
plt.axhline(y = 0,linewidth = 2.0,color = 'b',ls = "-.")
knodes = list(map(list, zip(*knodes)))
plt.xticks(knodes[1],list(knodes[0]),fontproperties='Times New Roman')
#******************** Wannier band ********************
wannier_bandfile = os.path.join(path,'wannier90_band.dat')
wannier_kptfile = os.path.join(path,'wannier90_band.kpt')
wannier_labkptfile = os.path.join(path,'wannier90_band.labelinfo.dat')
fermi_file = os.path.join(path,'FERMI_ENERGY')
wannier_print = True
if os.path.exists(wannier_bandfile) is True and wannier_print is True:
wann_k = open(wannier_kptfile,"r",encoding = "utf-8")
wann_l = open(wannier_labkptfile,"r",encoding = "utf-8")
n_k = np.array(int(wann_k.readlines()[0]))
lab = np.array(re.findall('[A-Z]*[A-Z]',wann_l.read()))
fermi_energy = np.loadtxt(fermi_file)
wann_k.close()
wann_l.close()
wb = np.loadtxt(wannier_bandfile)
wj = int(len(wb)/n_k)
for j in range(wj):
if j == wj - 2:
plt.plot(wb[:,0][j*n_k:j*n_k+n_k],wb[:,1][j*n_k:j*n_k+n_k] - fermi_energy,color = 'r',ls = '--',linewidth = 2.0,label = "Wannier")
else:
plt.plot(wb[:,0][j*n_k:j*n_k+n_k],wb[:,1][j*n_k:j*n_k+n_k] - fermi_energy,color = 'r',ls = '--',linewidth = 2.0)
plt.legend(loc='upper right', shadow = True, fancybox = True)
#********************************************************************************
plt.tick_params(axis='x',width = 0,length = 10)
plt.tick_params(axis='y',width = 0,length = 10)
ax = plt.gca()
ax.spines["bottom"].set_linewidth(1.5)
ax.spines["left"].set_linewidth(1.5)
ax.spines["right"].set_linewidth(1.5)
ax.spines["top"].set_linewidth(1.5)
# plt.show()
plt.savefig("bandstructure.png",dpi = 300,bbox_inches = 'tight',transparent=True)
#-----------------------------------------------------------------------------------------------------
def plotfs(mu,numk):
# 费米面以及极化率绘制
dataname = "fermi-mu-" + str(mu) + "-numk-" + str(numk) + ".dat"
dataname = "FS.dat"
# dataname = "Honeycomb-fermi-mu-" + str(mu) + ".dat"
# dataname = "Square-fermi-mu-" + str(mu) + ".dat"
# dataname = "fermi-mu-" + str(mu) + ".dat"
# dataname = "Honeycomb-fermivs-" + str(mu) + ".dat"
# dataname = "maxvals-chi0-" + str(numk) + ".dat"
# picname = os.path.splitext(dataname)[0] + "-" + str(num) + ".png"
picname = os.path.splitext(dataname)[0] + ".png"
da = np.loadtxt(dataname)
plt.figure(figsize = (10,8))
sc = plt.scatter(da[:,0],da[:,1], s = 1, c = da[:,2],cmap = "jet")
#-------------------------------------------------
r0 = 5
hex = [np.sqrt(3)/2 * r0,np.sqrt(3)/4 * r0,-np.sqrt(3)/4 * r0,-np.sqrt(3)/2 * r0,-np.sqrt(3)/4 * r0,np.sqrt(3)/4 * r0,np.sqrt(3)/2 * r0]
hey = [0 * r0 ,3/4 * r0 ,3/4 * r0 ,0 * r0 ,-3/4 * r0 ,-3/4 * r0 ,0 * r0]
plt.plot(hex,hey,c = "black",lw = 2,ls = "--")
#-------------------------------------------------
cb = plt.colorbar(sc,fraction = 0.1,ticks = [np.min(da[:,2]),np.max(da[:,2])]) # 调整colorbar的大小和图之间的间距
# cb.ax.set_yticklabels([r'$d_{3z^2-r^2}$', r'$d_{x^2-y^2}$'])
xtic = [-2 * np.pi,0,2 * np.pi]
xticlab = ["$-2\pi$","$0$","$2\pi$"]
plt.xticks(xtic,list(xticlab),fontproperties='Times New Roman', size = 40)
plt.yticks(xtic,list(xticlab),fontproperties='Times New Roman', size = 40)
xmin = np.min(da[:,0])
xmax = np.max(da[:,0])
ymin = np.min(da[:,1])
ymax = np.max(da[:,1])
# plt.xlim(xmin,xmax)
# plt.ylim(ymin,ymax)
plt.xlim(-2 * np.pi,2 * np.pi)
plt.ylim(-2 * np.pi,2 * np.pi)
plt.tick_params(axis='x',width = 0,length = 10)
plt.tick_params(axis='y',width = 0,length = 10)
ax = plt.gca()
ax.spines["bottom"].set_linewidth(1.5)
ax.spines["left"].set_linewidth(1.5)
ax.spines["right"].set_linewidth(1.5)
ax.spines["top"].set_linewidth(1.5)
# plt.show()
plt.savefig(picname, dpi = 300,bbox_inches = 'tight',transparent=True)
plt.close()
#-----------------------------------------------------------------------------------------------------
# WannierDETband()
numk = 100
mu = -2.15
plotfs(format(mu,".2f"),format(numk,".2f"))
上面的所有文件可以点击这里下载
公众号
相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
|
yxli406@gmail.com |