主题:循环出错?
m为整数,比如4;h2(m,m),tt(m,m),h0(m,m)都是已知矩阵,为什么下面两个操作看似相似,得到的结果却不同呢?E为单位矩阵
(1)gl=-h2.x.tt-h0
gl=E+gl
gl=.i.gl
(2)gl=-h2.x.tt
gl=gl-h0
gl=E+gl
gl=.i.gl
而且,还不能直接这样操作:
gl=.i.(-h2.x.tt-h0+E)
why?
昨晚回去睡觉了,抱歉啊。呵呵。
程序代码:
program main
use IMSL
implicit none
parameter pi=3.1415926
real e,conductance,eps
integer:: i,n,m,j,u
parameter(u=3,n=2**u,m=6)
real h(2*m*n,2*m*n),h0(2*m,2*m),h2(2*m,2*m),h1(2*m,2*m),hl(2*m*n,2*m),hr(2*m*n,2*m),hclc(2*m,2*m*n),hcrc(2*m,2*m*n) !哈密顿矩阵相关
real gr(2*m,2*m),gl(2*m,2*m),gcdu(2*m*n,2*m*n),gcr(2*m*n,2*m*n),gca(2*m*n,2*m*n),gg(2*m,2*m) !格林函数相关矩阵
real l(2*m*n,2*m*n),r(2*m*n,2*m*n),lc(2*m*n,2*m*n),rc(2*m*n,2*m*n) !自能矩阵及其厄米
real t0(2*m,2*m),t1(2*m,2*m),tt(2*m,2*m),t2(2*m,2*m),tt_(2*m,2*m),t3(2*m,2*m),t4(2*m,2*m) ,t(2*m,2*m),t5(2*m,2*m)!转移矩阵相关
real energy1(2*m,2*m),energy2(2*m*n,2*m*n) !与能量相关对角矩阵
real tl(2*m*n,2*m*n),tr(2*m*n,2*m*n) !耦合矩阵
real identity(2*m,2*m),eigenvalue(2*m*n)
eps=1.0e-3
hl=0.0
hr=0.0
identity=0.0
do i=1,2*m
identity(i,i)=1
end do
call hmax(h,h0,h1,h2,n,m)
! open(12,file='tezz.txt')
! eigenvalue=eig(h)
! write(12,*) eigenvalue
! close(12)
!对左右两边相互作用矩阵Hcl,Hcr及其共轭
hl(1:2*m,1:2*m)=h2
hr(2*m*(n-1)+1:2*m*n,1:2*m)=h1
hclc=transpose(hl)
hcrc=transpose(hr)
! do e=-10.0,10.0,0.0043
energy1=0.0
energy2=0.0
conductance=0.0
gcdu=0.0
e=0.03
do i=1,2*m
energy1(i,i)=e
end do
do i=1,2*m*n
energy2(i,i)=e
end do
!求转移矩阵tt,tt_
t0=energy1-h0
t0=.i.t0
t0=t0.x.h2 !求t0,初始值
tt=t0
t1=energy1-h0
t1=.i.t1
t1=t1.x.h1 !求t0_,初始值
t2=t1
do i=2,u
t4=-t1.x.t0
t5=-t0.x.t1
t4=identity+t4+t5
t4=.i.t4
t0=t4.x.(t0**2) !ti
t2=t2.x.t0
tt=tt+t2
t1=t4.x.(t1**2) !t1_
if(all(t1<eps) .and. all(t0<eps)) exit
t2=t2.x.t1
end do
t0=energy1-h0
t0=.i.t0
t0=t0.x.h2 !求t0,初始值
t2=t0
t1=energy1-h0
t1=.i.t1
t1=t1.x.h1 !求t0_,初始值
tt_=t1
do i=2,u
t3=-t1.x.t0
t5=-t0.x.t1
t3=identity+t5+t3
t3=.i.t3
t1=t3.x.(t1**2) !ti_
t2=t2.x.t1
tt_=tt_+t2
t0=t3.x.(t0**2) !ti
if(all(t1<eps) .and. all(t0<eps)) exit
t2=t2.x.t0
end do
!表面格林函数
gr=-h1.x.tt-h0
gr=energy1+gr
gr=.i.gr
! gr=.i.(energy1-h0-h1.x.tt)
! gl=.i.(energy1-h0-h1.x.tt_)
gl=-h2.x.tt_-h0
gl=energy1+gl
gl=.i.gl
!求自能函数及其厄米
l=(hl.x.gl).x.hclc
r=(hr.x.gr).x.hcrc
lc=transpose(l)
rc=transpose(r)
! open(12,file='ttt_.txt')
! write(12,*) lc
! close(12)
!材料自身格林函数Gcr,Gca
gcr=-h-l-r+gcr
! gcr=energy2+gcr
open(12,file='ttt_.txt')
write(12,*) gcr
close(12)
gcr=.i.gcr
! gcr=.i.(energy2-h-l-r)
gca=transpose(gcr)
!求耦合函数TL,TR
TL=l-lc
TR=r-rc
!求电导
open(1,file='econ.txt')
gcdu=-((TL.x.gcr).x.(TR.x.gca))
! eigenvalue=eig(gcdu)
write(1,*) gcdu
! do i=1,2*m*n
! conductance=conductance+gcdu(i,i)
! end do
! write(1,*) e,2.0*pi*conductance
! end do
close(1)
stop
contains
subroutine hmax(h,h0,h1,h2,n,m)
real v
integer n,m,i,j
real::h(2*m*n,2*m*n)
real h00(m,m),h01(m,m),h10(m,m)
real h0(2*m,2*m)
real h1(2*m,2*m),h2(2*m,2*m)
v=-3.0
h00=0.0
h01=0.0
h=0.0
h0=0.0
h1=0.0
h2=0.0
h10=0.0
!y方向有m条横向链,x有n条带,所以一个大单胞总碳原子个数为2mn
do i=1,m
if(mod(i,2)==0) then
h00(i,i-1)=v
else
h00(i,i+1)=v
end if
end do
!右上角h01,左边列对右边列作用
do i=1,m
if(i==1) then
h01(i,m)=v
else
h01(i,i-1)=v
end if
end do
!左下角h10,同一条带右边列对左边列作用矩阵
! h10=transpose(h01)
do i=1,m
if(i==m) then
h10(i,1)=v
else
h10(i,i+1)=v
end if
end do
h0(1:m,1:m)=h00
h0(1:m,m+1:2*m)=h01
h0(m+1:2*m,1:m)=h10
h0(m+1:2*m,m+1:2*m)=h00
h1(m+1:2*m,1:m)=h10 !与右边带相互作用矩阵
h2(1:m,m+1:2*m)=h01 !与左边带相互作用矩阵
do i=1,2*m*(n-1)+1,2*m
h(i:i+2*m-1,i:i+2*m-1)=h0
if(i==1) then
h(1:2*m,2*m+1:4*m)=h1
! h(1:2*m,(n-1)*2*m+1:2*m*n)=h2
else if(i==(n-1)*2*m+1) then
h(i:2*m*n,i-2*m:i-1)=h2
! h(i:2*m*n,1:2*m)=h1
else
h(i:i+2*m-1,i+2*m:i+4*m-1)=h1
h(i:i+2*m-1,i-2*m:i-1)=h2
end if
end do
end subroutine
end program main
(1)gl=-h2.x.tt-h0
gl=E+gl
gl=.i.gl
(2)gl=-h2.x.tt
gl=gl-h0
gl=E+gl
gl=.i.gl
而且,还不能直接这样操作:
gl=.i.(-h2.x.tt-h0+E)
why?
昨晚回去睡觉了,抱歉啊。呵呵。
程序代码:
program main
use IMSL
implicit none
parameter pi=3.1415926
real e,conductance,eps
integer:: i,n,m,j,u
parameter(u=3,n=2**u,m=6)
real h(2*m*n,2*m*n),h0(2*m,2*m),h2(2*m,2*m),h1(2*m,2*m),hl(2*m*n,2*m),hr(2*m*n,2*m),hclc(2*m,2*m*n),hcrc(2*m,2*m*n) !哈密顿矩阵相关
real gr(2*m,2*m),gl(2*m,2*m),gcdu(2*m*n,2*m*n),gcr(2*m*n,2*m*n),gca(2*m*n,2*m*n),gg(2*m,2*m) !格林函数相关矩阵
real l(2*m*n,2*m*n),r(2*m*n,2*m*n),lc(2*m*n,2*m*n),rc(2*m*n,2*m*n) !自能矩阵及其厄米
real t0(2*m,2*m),t1(2*m,2*m),tt(2*m,2*m),t2(2*m,2*m),tt_(2*m,2*m),t3(2*m,2*m),t4(2*m,2*m) ,t(2*m,2*m),t5(2*m,2*m)!转移矩阵相关
real energy1(2*m,2*m),energy2(2*m*n,2*m*n) !与能量相关对角矩阵
real tl(2*m*n,2*m*n),tr(2*m*n,2*m*n) !耦合矩阵
real identity(2*m,2*m),eigenvalue(2*m*n)
eps=1.0e-3
hl=0.0
hr=0.0
identity=0.0
do i=1,2*m
identity(i,i)=1
end do
call hmax(h,h0,h1,h2,n,m)
! open(12,file='tezz.txt')
! eigenvalue=eig(h)
! write(12,*) eigenvalue
! close(12)
!对左右两边相互作用矩阵Hcl,Hcr及其共轭
hl(1:2*m,1:2*m)=h2
hr(2*m*(n-1)+1:2*m*n,1:2*m)=h1
hclc=transpose(hl)
hcrc=transpose(hr)
! do e=-10.0,10.0,0.0043
energy1=0.0
energy2=0.0
conductance=0.0
gcdu=0.0
e=0.03
do i=1,2*m
energy1(i,i)=e
end do
do i=1,2*m*n
energy2(i,i)=e
end do
!求转移矩阵tt,tt_
t0=energy1-h0
t0=.i.t0
t0=t0.x.h2 !求t0,初始值
tt=t0
t1=energy1-h0
t1=.i.t1
t1=t1.x.h1 !求t0_,初始值
t2=t1
do i=2,u
t4=-t1.x.t0
t5=-t0.x.t1
t4=identity+t4+t5
t4=.i.t4
t0=t4.x.(t0**2) !ti
t2=t2.x.t0
tt=tt+t2
t1=t4.x.(t1**2) !t1_
if(all(t1<eps) .and. all(t0<eps)) exit
t2=t2.x.t1
end do
t0=energy1-h0
t0=.i.t0
t0=t0.x.h2 !求t0,初始值
t2=t0
t1=energy1-h0
t1=.i.t1
t1=t1.x.h1 !求t0_,初始值
tt_=t1
do i=2,u
t3=-t1.x.t0
t5=-t0.x.t1
t3=identity+t5+t3
t3=.i.t3
t1=t3.x.(t1**2) !ti_
t2=t2.x.t1
tt_=tt_+t2
t0=t3.x.(t0**2) !ti
if(all(t1<eps) .and. all(t0<eps)) exit
t2=t2.x.t0
end do
!表面格林函数
gr=-h1.x.tt-h0
gr=energy1+gr
gr=.i.gr
! gr=.i.(energy1-h0-h1.x.tt)
! gl=.i.(energy1-h0-h1.x.tt_)
gl=-h2.x.tt_-h0
gl=energy1+gl
gl=.i.gl
!求自能函数及其厄米
l=(hl.x.gl).x.hclc
r=(hr.x.gr).x.hcrc
lc=transpose(l)
rc=transpose(r)
! open(12,file='ttt_.txt')
! write(12,*) lc
! close(12)
!材料自身格林函数Gcr,Gca
gcr=-h-l-r+gcr
! gcr=energy2+gcr
open(12,file='ttt_.txt')
write(12,*) gcr
close(12)
gcr=.i.gcr
! gcr=.i.(energy2-h-l-r)
gca=transpose(gcr)
!求耦合函数TL,TR
TL=l-lc
TR=r-rc
!求电导
open(1,file='econ.txt')
gcdu=-((TL.x.gcr).x.(TR.x.gca))
! eigenvalue=eig(gcdu)
write(1,*) gcdu
! do i=1,2*m*n
! conductance=conductance+gcdu(i,i)
! end do
! write(1,*) e,2.0*pi*conductance
! end do
close(1)
stop
contains
subroutine hmax(h,h0,h1,h2,n,m)
real v
integer n,m,i,j
real::h(2*m*n,2*m*n)
real h00(m,m),h01(m,m),h10(m,m)
real h0(2*m,2*m)
real h1(2*m,2*m),h2(2*m,2*m)
v=-3.0
h00=0.0
h01=0.0
h=0.0
h0=0.0
h1=0.0
h2=0.0
h10=0.0
!y方向有m条横向链,x有n条带,所以一个大单胞总碳原子个数为2mn
do i=1,m
if(mod(i,2)==0) then
h00(i,i-1)=v
else
h00(i,i+1)=v
end if
end do
!右上角h01,左边列对右边列作用
do i=1,m
if(i==1) then
h01(i,m)=v
else
h01(i,i-1)=v
end if
end do
!左下角h10,同一条带右边列对左边列作用矩阵
! h10=transpose(h01)
do i=1,m
if(i==m) then
h10(i,1)=v
else
h10(i,i+1)=v
end if
end do
h0(1:m,1:m)=h00
h0(1:m,m+1:2*m)=h01
h0(m+1:2*m,1:m)=h10
h0(m+1:2*m,m+1:2*m)=h00
h1(m+1:2*m,1:m)=h10 !与右边带相互作用矩阵
h2(1:m,m+1:2*m)=h01 !与左边带相互作用矩阵
do i=1,2*m*(n-1)+1,2*m
h(i:i+2*m-1,i:i+2*m-1)=h0
if(i==1) then
h(1:2*m,2*m+1:4*m)=h1
! h(1:2*m,(n-1)*2*m+1:2*m*n)=h2
else if(i==(n-1)*2*m+1) then
h(i:2*m*n,i-2*m:i-1)=h2
! h(i:2*m*n,1:2*m)=h1
else
h(i:i+2*m-1,i+2*m:i+4*m-1)=h1
h(i:i+2*m-1,i-2*m:i-1)=h2
end if
end do
end subroutine
end program main