program main
integer :: i,a,nn,ll,nst,x,y,ii
integer :: N,N1,N2,N3,N4 !!三角晶格数
real(8) :: jj !!哈密顿量系数虚部
integer, allocatable :: spin(:) !!自旋
real(8) :: Z,T,k,op1,J,sign
real(8), allocatable :: W(:),OP(:),Hr(:,:),Hi(:,:) !!权重，序参量，哈密顿量(能量)
integer :: s1,s2,s3,s4,s5,s6 !!邻居
real(8) :: r1,r2,r3,r4 !!位形占比
real(8) :: ee !!期望值
Z=0
J=0.4 !!哈密顿量系数实部
jj=0.5d0*acos(sinh(2.d0*J)) !! fi
ll=3 !!长度
nn=ll*ll
N=nn*2 !!三角晶格总数
nst=2**nn !!2的n次方个态，指数用双星号
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))
Hr(:,0)=0
Hr(:,1)=0
Hi(:,0)=0
Hi(:,1)=0
ee=0
do a=0,nst-1
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)=dble(spin(i)*s5*s6)+Hr(a,0)
Hr(a,1)=dble(spin(i)*s1*s3)+Hr(a,1)
Hi(a,0)=dble(spin(i)*s5*s6)+Hi(a,0)
Hi(a,1)=dble(spin(i)*s1*s3)+Hi(a,1)
ii=spin(i)
if(ii*s1*s3==1.or.ii*s5*s6==1) then
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
if(ii*s1==1.and.ii*s3==1) N1=N1+1
if(ii*s5==1.and.ii*s6==1) N1=N1+1
if(ii*s1==-1.and.ii*s3==-1) N2=N2+1
if(ii*s5==-1.and.ii*s6==-1) N2=N2+1
if(ii*s1==-1.and.ii*s3==1) N3=N3+1
if(ii*s5==-1.and.ii*s6==1) N3=N3+1
if(ii*s1==1.and.ii*s3==-1) N4=N4+1
if(ii*s5==1.and.ii*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
if(ii*s1==1.and.ii*s3==1) N1=N1+1
if(ii*s5==1.and.ii*s6==1) N1=N1+1
if(ii*s1==1.and.ii*s3==-1) N2=N2+1
if(ii*s5==1.and.ii*s6==-1) N2=N2+1
if(ii*s1==-1.and.ii*s3==-1) N3=N3+1
if(ii*s5==-1.and.ii*s6==-1) N3=N3+1
if(ii*s1==-1.and.ii*s3==1) N4=N4+1
if(ii*s5==-1.and.ii*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
if(ii*s1==1.and.ii*s3==1) N1=N1+1
if(ii*s5==1.and.ii*s6==1) N1=N1+1
if(ii*s1==-1.and.ii*s3==1) N2=N2+1
if(ii*s5==-1.and.ii*s6==1) N2=N2+1
if(ii*s1==1.and.ii*s3==-1) N3=N3+1
if(ii*s5==1.and.ii*s6==-1) N3=N3+1
if(ii*s1==-1.and.ii*s3==-1) N4=N4+1
if(ii*s5==-1.and.ii*s6==-1) N4=N4+1
end if
end if
enddo
Hr(a,2)=-Hr(a,0)-Hr(a,1)
Hi(a,2)=-Hi(a,0)+Hi(a,1)
W(a)=exp(J*Hr(a,2))*cos(jj*Hi(a,2))
sign=cos(jj*Hi(a,2))/abs(cos(jj*Hi(a,2)))
r1=dble(N1)/dble(N)
r2=dble(N2)/dble(N)
r3=dble(N3)/dble(N)
r4=dble(N4)/dble(N)
OP(a)=(r1-r2)**2+(r1-r3)**2+(r1-r4)**2
OP(a)=OP(a)+(r2-r3)**2+(r2-r4)**2+(r3-r4)**2
OP(a)=OP(a)/3.d0*W(a)
ee=ee+W(a)*sign
enddo
Z=sum(W(:))
ee=ee/Z
op1=sum(OP(:))/Z
print *,op1
end