四方点阵沿对角线方向开边界
这里研究一下一个square lattice,如何沿对角线方向取开边界条件,研究这种情况下的边界态是怎样的,并介绍一下如何在一个四方点阵的基础上,变成可以沿对角线开边界的模型.
{:.info}
前言
在通常的研究中,我经常遇到的是一个四方点阵上的紧束缚模型,这个时候想要看边界态,只需要将哈密顿量在一个方向取周期边界条件,另外一个方向取开边界条件即可.关于这种形式的问题,可以参考Chern Insulator边界态及Chern数计算这篇博客,这里主要是研究怎么对一个正方点阵上的紧束缚模型,沿对角线方向开边界.
坐标系旋转
如图,黑色坐标系表示确定四方点阵的直角坐标$(k_x,k_y)$,红色虚线坐标系来确定一个旋转$45^o$后的坐标系$(k_x^{‘},k_y^{‘})$,从坐标系可以清楚的看到,对于$(k_x,k_y)$的直角坐标,其对角线方向正好就是$(k_x^{‘},k_y^{‘})$的$k_x^{‘}$方向.
所以这里最核心的思想就是将原来的直角坐标$(k_x,k_y)$旋转$45^o$变成对应的$(k_x^{‘},k_y^{‘})$坐标,在$(k_x^{‘},k_y^{‘})$的表示下沿着$k_x^{‘}$依照原来四方点阵的方法取开边界条件即可.
{:.success}
接下来以Chern Insulator边界态及Chern数计算这篇博客中的Chern Insulator模型来作为实例来计算,在$(k_x,k_y)$的表示下
在$(k_x,k_y)$旋转$45^o$变成对应的$(k_x^{‘},k_y^{‘})$坐标的时候,它们之间的变化关系为
将这个关系代入之后,即可以将哈密顿量(\ref{eq1})变为
坐标旋转之后,哈密顿量又变成了关于$(k_x^{‘},k_y^{‘})$两个坐标变量的形式,这时候如果想沿$(k_x,k_y)$的对角线方向开边界,则只需要对哈密顿量(\ref{eq2})沿$k_x^{‘}$方向取开边界即可,及第一张示意图所示,剩下的问题就可Chern Insulator边界态及Chern数计算这篇博客中开边界算边界态的过程一样了.
代码
这里我线用fortran写了一下哈密顿量(\ref{eq2})的内容,然后计算了对应的边界态1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165! Author: YuXuanLi
! Email:yxli406@gmail.com
module pub
implicit none
integer yn,kn,hnn
parameter(yn = 50,kn = 30,hnn = 2)
integer,parameter::N = yn*hnn
real,parameter::pi = 3.1415926535
complex,parameter::im = (0.,1.0)
complex::Ham(N,N) = 0
complex g1(hnn,hnn),g2(hnn,hnn),g3(hnn,hnn)
!=================================
real m0,tx,ty,lamx,lamy
!================cheevd===============
integer::lda = N
integer,parameter::lwmax=2*N+N**2
real,allocatable::w(:)
complex,allocatable::work(:)
real,allocatable::rwork(:)
integer,allocatable::iwork(:)
integer lwork
integer lrwork
integer liwork
integer info
end module pub
!============================================================
program sol
use pub
allocate(w(N))
allocate(work(lwmax))
allocate(rwork(1+5*N+2*N**2))
allocate(iwork(3+5*N))
!------------------------------------
m0 = 0.5
tx = 1.0
ty = 1.0
lamx = 1.0
lamy = 1.0
call main1()
stop
end program sol
!============================================================
subroutine main1()
use pub
integer m1
real k
open(3,file="openx.dat")
! open(4,file="openy-m1.dat")
do m1 = -kn,kn
k = pi*m1/kn*sqrt(2.0)
call openx(k)
write(3,999)k/pi/sqrt(2.0),(w(i),i = 1,N)
! call openy(k)
! write(4,999)k/pi,(w(i),i = 1,N)
end do
close(3)
! close(4)
999 format(201f11.6)
end subroutine main1
!============================================================
subroutine openx(ky)
use pub
real ky
call pauli()
Ham = 0
!========== Positive energy ========
do k = 0,yn-1
if (k == 0) then ! Only right block in first line
do m = 1,hnn
do l = 1,hnn
Ham(m,l) = m0*g3(m,l)
Ham(m,l + hnn) = (tx*(cos(sqrt(2.0)/2.0*ky) - 1/im*sin(sqrt(2.0)/2.0*ky)) +&
ty*(cos(sqrt(2.0)/2.0*ky) + 1/im*sin(sqrt(2.0)/2.0*ky)))*g3(m,l)+&
lamx*(1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g1(m,l)+&
lamy*(-1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g2(m,l)
end do
end do
elseif ( k==yn-1 ) then ! Only left block in last line
do m = 1,hnn
do l = 1,hnn
Ham(k*hnn + m,k*hnn + l) = m0*g3(m,l)
Ham(k*hnn + m,k*hnn + l - hnn) = (tx*(cos(sqrt(2.0)/2.0*ky) + 1/im*sin(sqrt(2.0)/2.0*ky)) +&
ty*(cos(sqrt(2.0)/2.0*ky) - 1/im*sin(sqrt(2.0)/2.0*ky)))*g3(m,l)+&
lamx*(-1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g1(m,l)+&
lamy*(1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g2(m,l)
end do
end do
else
do m = 1,hnn
do l = 1,hnn ! k start from 1,matrix block from 2th row
Ham(k*hnn + m,k*hnn + l) = m0*g3(m,l)
Ham(k*hnn + m,k*hnn + l + hnn) = (tx*(cos(sqrt(2.0)/2.0*ky) - 1/im*sin(sqrt(2.0)/2.0*ky)) +&
ty*(cos(sqrt(2.0)/2.0*ky) + 1/im*sin(sqrt(2.0)/2.0*ky)))*g3(m,l)+&
lamx*(1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g1(m,l)+&
lamy*(-1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g2(m,l)
Ham(k*hnn + m,k*hnn + l - hnn) = (tx*(cos(sqrt(2.0)/2.0*ky) + 1/im*sin(sqrt(2.0)/2.0*ky)) +&
ty*(cos(sqrt(2.0)/2.0*ky) - 1/im*sin(sqrt(2.0)/2.0*ky)))*g3(m,l)+&
lamx*(-1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g1(m,l)+&
lamy*(1/im*cos(sqrt(2.0)/2.0*ky) + sin(sqrt(2.0)/2.0*ky))*g2(m,l)
end do
end do
end if
end do
!------------------------
call isHermitian()
call eigsol()
return
end subroutine openx
!============================================================
subroutine pauli()
use pub
g1(1,2) = 1
g1(2,1) = 1
!-----------------
g2(1,2) = -im
g2(2,1) = im
!---------------
g3(1,1) = 1
g3(2,2) = -1
end subroutine pauli
!============================================================
subroutine isHermitian()
use pub
integer i,j
do i = 1,N
do j = 1,N
if (Ham(i,j) .ne. conjg(Ham(j,i)))then
open(16,file = 'hermitian.dat')
write(16,*)i,j
write(16,*)Ham(i,j)
write(16,*)Ham(j,i)
write(16,*)"===================="
write(*,*)"Ham isn't Hermitian"
stop
end if
end do
end do
close(16)
return
end subroutine isHermitian
!================= 矩阵本征值求解 ==============
subroutine eigSol()
use pub
integer m
lwork = -1
liwork = -1
lrwork = -1
call cheevd('V','Upper',N,Ham,lda,w,work,lwork &
,rwork,lrwork,iwork,liwork,info)
lwork = min(2*N+N**2, int( work( 1 ) ) )
lrwork = min(1+5*N+2*N**2, int( rwork( 1 ) ) )
liwork = min(3+5*N, iwork( 1 ) )
call cheevd('V','Upper',N,Ham,lda,w,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 eigSol
计算结果如下图所示
为了检验这个方法的正确性,我利用WannierTools计算Chern绝缘体性质中使用过了WannierTools,在控制参数中将开边界方向取在了对角线方向上,得到了相同的能带图.
{:.info}
1 | &TB_FILE |
至于计算所用的紧束缚模型的数据ChernInsulator_hr.dat,其构造方法可以查阅WannierTools计算Chern绝缘体性质这篇博客,具体数据内容如下
1 | Chern Insulator |
将ChernInsulator_hr.dat与wt.in放置到相同的文件夹中,运行WannierTools即可1
mpirun -np 4 wt.x wt.in &
最后会得到slabek.gnu,slabek.dat这两个文件,利用gnuplot绘图1
gnuplot slabek.gnu
计算结束之后的结果如下图
所以这里的方法是完全正确的.
公众号
相关内容均会在公众号进行同步，若对该Blog感兴趣，欢迎关注微信公众号。
{:.info}
|
yxliphy@gmail.com |
wechat
alipay
