回 帖 发 新 帖 刷新版面

主题:循环出错?

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

回复列表 (共8个回复)

沙发

energy1也是单位矩阵吗?

板凳

energy1也是单位矩阵吗?

3 楼

如果贴代码不方便的话, 发问的时候可否顺手贴出变量的定义?

4 楼

[quote]energy1也是单位矩阵吗?[/quote]
是的,也是单位矩阵[em2][em2][em2][em2][em2][em2]

5 楼

你原来的e跟energy1不是一样的. e一个浮点数而已, e跟矩阵相加与对角矩阵energy1跟矩阵相加肯定不同了.


我明白了~

 gl=-h2.x.tt_-h0  是等价于  gl=-h2.x.(tt_-h0)
我设了一个gtemp, 维度跟gl一样然后
 gl=-h2.x.tt_-h0
 gtemp = -h2.x.(tt_-h0)
write(*, *) gl-gtemp
这样的话结构全部是0, 但如果没有那个括号得到的结构有3, 有9.

BTW: 楼主, 我觉得你大量使用这些自定义运算符的时候测试一下吧, 直接上马程序后期调试死你. 而且对于不清楚优先级的运算或者不记得甚至是对优先级模糊, 顺手加上括号, 括号在编译的时候起作用不会减慢你的计算的. 既然老板要你写fortran希望你踏踏实实来.

6 楼

注意, 你的e跟energy1 都不是单位矩阵, 关键是计算中它们不能替换.

7 楼


初一看,问题就出在E和energy1,不出所料;其他地方,不明白其物理内涵,也就不知道有无问题了。

8 楼

水落石出,看来还是要好好搞的

我来回复

您尚未登录,请登录后再回复。点此登录或注册