编辑代码

program main
	integer :: i,a,b,nn,ll,nst,x,y,m
    integer, allocatable :: spin(:),H(:),mm(:) !!自旋&哈密顿量(能量)&磁化强度
    real(8) :: Z,T,k,E
    real(8), allocatable :: W(:)

    k=1 !!玻尔兹曼常数
    Z=0 
    E=0 !!用于计算能量密度
    ll=2 !!长度
    nn=ll*ll
    nst=2**nn !!2的n次方个态，指数用双星号
    T=2.00 !!温度
    m=0 !!磁化强度绝对值的统计平均
    allocate (spin(0:nn-1))
    allocate (W(0:nst-1))
    allocate (H(0:nst-1))
    allocate (mm(0:nst-1))
    do a=0,nst-1
        H(a)=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)+y*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)
            H(a)=(-1)*spin(i)*(s1+s2+s3+s4)+H(a)
        enddo
        mm(a)=abs(sum(spin))/nn
        H(a)=H(a)/2
        W(a)=exp(-H(a)/k/T)
        E=H(a)*W(a)+E
        Z=Z+W(a)
        m=mm(a)*W(a)+m
    enddo
    E=E/Z
    print *, "能量密度为：",E
    print *, "<|m|>为：",m/Z

    do i=0,nst-1
    print *,H(i),W(i)
    enddo

    
end