编辑代码

program main
    integer :: i,a,J,nn,ll,nst,x,y
    integer :: N,N1,N2,N3,N4             !!三角晶格数
    integer :: jj                        !!哈密顿量系数虚部
    integer, allocatable :: spin(:),Hr(:,:),Hi(:,:) !!自旋&哈密顿量(能量)
    real(8) :: Z,T,k
    real(8), allocatable :: W(:),OP(:)    !!权重，序参量
    integer :: s1,s2,s3,s4,s5,s6          !!邻居
    real(8) :: r1,r2,r3,r4                !!位形占比

    !k=1 !!玻尔兹曼常数
    Z=0 
    J=1  !!哈密顿量系数实部
    jj=1
    ll=2 !!长度
    nn=ll*ll
    N=(ll-1)**2*2    !!三角晶格总数
    nst=2**nn !!2的n次方个态，指数用双星号
    T=2.00 !!温度
    allocate (spin(0:nn-1))
    allocate (W(0:nst-1))
    allocate (Hr(0:nst-1,0:2))        !!哈密顿量实部，0上三角，1下三角,2和
    allocate (Hi(0:nst-1,0:2))        !!哈密顿量虚部，0上三角，1下三角,2和
    allocate (OP(0:nst-1))
    do a=0,nst-1
        Hr(a,0)=0
        Hr(a,1)=0
        Hi(a,0)=0
        Hi(a,1)=0
        N1=0
        N2=0
        N3=0
        N4=0
        do i=0,nn-1
            spin(i)=btest(a,i)
            spin(i)=2*spin(i)-1
        enddo
        do i=0,nn-1
            x=mod(i,ll)
            y=i/ll
            s1=spin(mod(x-1+ll,ll)+mod(y-1+ll,ll)*ll)   !!左 上
            s2=spin(x+mod(y+1,ll)*ll)                   !!右
            s3=spin(mod(x-1+ll,ll)+y*ll)                !!上
            s4=spin(x+mod(y-1+ll,ll)*ll)                !!左
            s5=spin(mod(x+1,ll)+y*ll)                   !!下
            s6=spin(mod(x+1,ll)+mod(y+1,ll)*ll)         !!右 下
            !!这种记录方法已经考虑边界周期性
            Hr(a,0)=(-1)*J*spin(i)*s5*s6+Hr(a,0)
            Hr(a,1)=(-1)*J*spin(i)*s1*s3+Hr(a,1)
            Hi(a,0)=(-1)*jj*spin(i)*s5*s6+Hi(a,0)
            Hi(a,1)=(-1)*jj*spin(i)*s1*s3+Hi(a,1)
            if((mod(x,3)==0.and.mod(y,3)==0).or.(mod(x,3)==1.and.mod(y,3)==2).or.(mod(x,3)==2.and.mod(y,3)==1) then
                !判断格点为A
                if((i*s1==1.and.i*s3==1).or.(i*s5==1.and.i*s6==1)) N1=N1+1
                if((i*s1==-1.and.i*s3==-1).or.(i*s5==-1.and.i*s6==-1)) N2=N2+1
                if((i*s1==-1.and.i*s3==1).or.(i*s5==-1.and.i*s6==1)) N3=N3+1
                if((i*s1==1.and.i*s3==-1).or.(i*s5==1.and.i*s6==-1)) N4=N4+1
            elseif((mod(x,3)==0.and.mod(y,3)==1).or.(mod(x,3)==1.and.mod(y,3)==0).or.(mod(x,3)==2.and.mod(y,3)==2) then
                !判断格点为B
                if((i*s1==1.and.i*s3==1).or.(i*s5==1.and.i*s6==1)) N1=N1+1
                if((i*s1==-1.and.i*s3==-1).or.(i*s5==-1.and.i*s6==-1)) N2=N2+1
                if((i*s1==-1.and.i*s3==1).or.(i*s5==-1.and.i*s6==1)) N3=N3+1
                if((i*s1==1.and.i*s3==-1).or.(i*s5==1.and.i*s6==-1)) N4=N4+1
            elseif((mod(x,3)==0.and.mod(y,3)==2).or.(mod(x,3)==1.and.mod(y,3)==1).or.(mod(x,3)==2.and.mod(y,3)==0) then
                !判断格点为C
                if((i*s1==1.and.i*s3==1).or.(i*s5==1.and.i*s6==1)) N1=N1+1
                if((i*s1==-1.and.i*s3==-1).or.(i*s5==-1.and.i*s6==-1)) N2=N2+1
                if((i*s1==-1.and.i*s3==1).or.(i*s5==-1.and.i*s6==1)) N3=N3+1
                if((i*s1==1.and.i*s3==-1).or.(i*s5==1.and.i*s6==-1)) N4=N4+1
            end if
        enddo
        Hr(a,2)=Hr(a,0)+Hr(a,1)
        Hi(a,2)=Hi(a,0)-Hi(a,1)
        W(a)=exp(-2.d0*Hr(a,2)/T)*cos(dble(Hi(a,2)))
        !r1=dble(N1)/dble(N)
        !r2=dble(N2)/dble(N)
        !r3=dble(N3)/dble(N)
        !r4=dble(N4)/dble(N)
    enddo
    Z=sum(W(:))
    do i=0,nst-1
    print *,Hr(i,2),Hi(i,2),W(i)
    enddo
end