!研究2D ISING MODEL程序模擬中，有外場情況下，并外場不斷改變時(shí)，
 !對各種物理量可能帶來的影響
 !觀察各量在取定外場，隨溫度變化的發(fā)展規(guī)律，提取臨界點(diǎn)，并總結(jié)不同外場下的對比規(guī)律
 !
 !
module global	 !!常量聲明
  integer,parameter:: nrows=20,ncols=20	!2Dsize	若改變尺寸，只需要改變這里就可以
  integer,parameter::Tmcs=80000,Emcs=50000	 !Total monte carlo steps
  real,parameter::HT=15.0,LT=0.4,intervalT=0.2 !high temperature
  !!!!!!!!!!!!!!!!!!!!!!!
  real::exB=-5.0 !from 0.0 to 10,exB=mu*B/J,包含外場信息的比例系數(shù)
  !!!!!!!!!!!!!!!!!!!!!!!
end module global

 program  isingsizeEffect  !!主程序
 use global
 implicit none
 character(10)::t1,t2
 integer::i
 real::T
 integer,dimension(nrows+2,ncols+2)::s

 open(unit=11,file='parameter.dat',status='replace',action='write')
   write(11,'(A)')'all the parameters of 2DexternalHising program'
   write(11,'(A,I3,2X,A,I3,2X)')'nrows=',nrows,'ncols=',ncols
   write(11,'(A,I7,2X,A,I7,2X,a,f5.2)')'Tmcs=',Tmcs,'Emcs=',Emcs,'exB=',exB
   write(11,'(A,f5.2,2X,A,f5.2,2X,A,f5.2)')'HT=',HT,'LT=',LT,'intervalT=',intervalT
 open(unit=12,file='initialconfig.dat',status='replace',action='write')
 !取定溫度的初始位型，溫度在臨界點(diǎn)上或下
 open(unit=13,file='lastconfig1.dat',status='replace',action='write')
 open(unit=133,file='lastconfig2.dat',status='replace',action='write')
 !取定溫度的最終位型，溫度在臨界點(diǎn)上或下
 open(unit=14,file='magnetTime1.dat',status='replace',action='write') 
 open(unit=144,file='magnetTime2.dat',status='replace',action='write') 
 !取定溫度，磁化強(qiáng)度與時(shí)間（MCS）的關(guān)系，溫度在臨界點(diǎn)上或下
 open(unit=15,file='magnetTemp.dat',status='replace',action='write') 
 !單個(gè)自旋的磁化強(qiáng)度與溫度的函數(shù)關(guān)系(絕對值也輸出到該文件)
 open(unit=16,file='susceptibility.dat',status='replace',action='write')	
 !單個(gè)自旋的極化率與溫度的關(guān)系
 open(unit=17,file='energyTemp.dat',status='replace',action='write')	
 !單個(gè)自旋的能量與溫度的函數(shù)關(guān)系
 open(unit=18,file='heatcapacity.dat',status='replace',action='write')
 !單個(gè)自旋的比熱容與溫度的函數(shù)關(guān)系
 
   write(*,*)'2D Ising model in different temperature with external B'
   write(*,'(2x,a,i4,2x,a,i4,2x,a,f5.2)')'nrows=',nrows,'ncols=',ncols,'exB=',exB
   call date_and_time(time=t1)
   do i=1,70  !70  NT=int((HT-LT)/intervalT)
	  T=LT+intervalT*(i-1)
	  write(*,'(2x,a,f6.3)')'T=',T
	  call initial(T,s)
	  call isingMC(T,s)
   enddo
   call date_and_time(time=t2)
   write(11,'(2X,A,A,A,A)')"time=",t2,"-",t1
   close(11)   
   close(12)   
   close(13)
   close(133)    
   close(14)
   close(144)   
   close(15)  
   close(16) 
   close(17)
   close(18)
end program   isingsizeEffect

subroutine initial(T,s)	 !!初始位型,cold start
 use global
 implicit none
 real::T
 integer::i,j 
 integer,dimension(nrows+2,ncols+2)::s
   
   s=0
   write(12,'(a,2x,a,2x,a)')'i','j','s(i,j)'
   do i=1,nrows+2
      do j=1,ncols+2
	     s(i,j)=1
		 if(T.eq.5.5)write(12,'(i3,2x,i3,2x,i3)')i,j,s(i,j)
		 !只輸出一次就可以了
	  enddo 
   enddo
endsubroutine initial

subroutine isingMC(T,s)		!!MC模擬
 use global
 implicit none
 external duni	!調(diào)用外部函數(shù)
 double precision::duni		 !函數(shù)關(guān)鍵詞聲明，數(shù)據(jù)類型要一致
 integer::i,j,k,m,n,stat_mcs,imcs
 real::energy,energy_ave,energy2_ave
 real::magnet,magnet_ave,magnet2_ave
 real::T,heatcapacity,susceptibility,xxduni
 real::deltU,try_spin
 integer,dimension(nrows+2,ncols+2)::s

  stat_mcs=0
  energy=0
  energy_ave=0
  energy2_ave=0
  magnet=0
  magnet_ave=0
  magnet2_ave=0
  heatcapacity=0
  susceptibility=0
  !這些量在不同溫度進(jìn)行清零，相同溫度下，進(jìn)行累加
  do k=1,int(T*100+10)
	 xxduni=duni()
  enddo

  !!Monte Carlo抽樣
  do imcs=1,Tmcs
     if((imcs.ge.Emcs).and.(mod(imcs,10).eq.0))then
	 !平衡后，隔步抽樣，減小關(guān)聯(lián)效應(yīng)
	   stat_mcs=stat_mcs+1
	   magnet=sum(s(2:nrows+1,2:ncols+1))/(nrows*ncols*1.0)	!平均到單個(gè)自旋上
	   if(T.eq.1.0)write(14,'(i5,2x,f8.3)')imcs,magnet
	   if(T.eq.3.0)write(144,'(i5,2x,f8.3)')imcs,magnet
	   !統(tǒng)計(jì)臨界溫度上和臨界溫度下平均磁化強(qiáng)度與時(shí)間的關(guān)系
	   magnet_ave=magnet_ave+magnet
	   magnet2_ave=magnet2_ave+magnet**2
	   do i=2,nrows+1
	      do j=2,ncols+1
             energy=energy-s(i,j)*(s(i-1,j)+s(i+1,j)+s(i,j-1)+s(i,j+1));
		  enddo 
       enddo
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	   energy=energy/(nrows*ncols*2.0)-exB*magnet !考慮外場情況，平均到單個(gè)自旋上
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	   energy_ave=energy_ave+energy	
	   energy2_ave=energy2_ave+energy**2
	 endif

	 do i=1,nrows*ncols !one monte carlo step
        n=nint((nrows-1)*duni())+2
        m=nint((ncols-1)*duni())+2
        try_spin=-s(n,m)!做粒子的嘗試翻轉(zhuǎn)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        deltU=-try_spin*(s(n-1,m)+s(n+1,m)+s(n,m-1)+s(n,m+1))*2-exB*try_spin*2
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!考慮外場情況,再減去自旋翻轉(zhuǎn)后，磁場強(qiáng)度的變化量
		!計(jì)算出因?yàn)樵摿Ｗ拥淖孕D(zhuǎn)而引起的前后能量變化
		if(exp(-deltU/T).gt.duni())then	 !!取約化單位為kT/J
		  s(n,m)=try_spin
          if (n==2) s(nrows+2,m)=try_spin
          if (n==nrows+1) s(1,m)=try_spin
		  if (m==2) s(n,ncols+2)=try_spin
          if (m==ncols+1) s(n,1)=try_spin
		  !周期性邊界條件，因?yàn)槭亲罱泥徬嗷プ饔茫砸ＷC
          !第一行與倒數(shù)第二行相同分布
          !最后一行與第二行相同分布
          !第一列與倒數(shù)第二列相同分布
          !最后一列與第二列相同分布
        endif 
     enddo 
  enddo !The end of total MCS
   do i=1,nrows+2
      do j=1,ncols+2
		 if((T.eq.5.5).and.(s(i,j).eq.1))write(13,'(i3,2x,i3,2x,i2)')i,j,s(i,j)
		 if((T.eq.13.5).and.(s(i,j).eq.1))write(133,'(i3,2x,i3,2x,i2)')i,j,s(i,j)
	  enddo 
   enddo

  !!取定溫度下，各量在monte carlo抽樣統(tǒng)計(jì)平均值的計(jì)算
  !為方便對比不同尺寸，已經(jīng)平均到單個(gè)自旋上
  magnet_ave=(magnet_ave/stat_mcs); !平均磁化強(qiáng)度計(jì)算
  magnet2_ave=magnet2_ave/stat_mcs; !磁化強(qiáng)度的方均值
  susceptibility=(magnet2_ave-(magnet_ave)**2)/T; !磁化率

  energy_ave=energy_ave/stat_mcs; !平均能量計(jì)算
  energy2_ave=energy2_ave/stat_mcs; !能量方均值
  heatcapacity=((energy2_ave)-((energy_ave)**2))/(T**2); !比熱

  write(15,'(f5.2,2x,f8.3,2x,f8.3)')T,magnet_ave,abs(magnet_ave)
  write(16,'(f5.2,2x,f8.3)')T, susceptibility
  write(17,'(f5.2,2x,f8.3)')T, energy_ave
  write(18,'(f5.2,2x,f12.6)')T, heatcapacity

endsubroutine isingMC

DOUBLE PRECISION FUNCTION DUNI()
      DOUBLE PRECISION CSAVE,CD,CM
      PARAMETER(&
     &CSAVE=0.9162596898123D+13/0.140737488355328D+15,&
     &CD=0.76543212345678D+14/0.140737488355328D+15,&
     &CM=0.140737488355213D+15/0.140737488355328D+15)
      DOUBLE PRECISION U(17),S,T,DUSTAR,C,DUNIB
      INTEGER I,J,II,JJ,K,KK,I1,J1,K1,L1,M1,ISEED 
      SAVE U,I,J,K,C
      DATA U/&
     &0.471960981577884755837789724978D+00,&
     &0.930323453205669578433639632431D+00,&
     &0.110161790933730836587127944899D+00,&
     &0.571501996273139518362638757010D-01,&
     &0.402467554779738266237538503137D+00,&
     &0.451181953427459489458279456915D+00,&
     &0.296076152342721102174129954053D+00,&
     &0.128202189325888116466879622359D-01,&
     &0.314274693850973603980853259266D+00,&
     &0.335521366752294932468163594171D-02,&
     &0.488685045200439371607850367840D+00,&
     &0.195470426865656758693860613516D+00,&
     &0.864162706791773556901599326053D+00,&
     &0.335505955815259203596381170316D+00,&
     &0.377190200199058085469526470541D+00,&
     &0.400780392114818314671676525916D+00,&
     &0.374224214182207466262750307281D+00/
      DATA I,J,K,C/17,5,47,CSAVE/
      DUNI = U(I)-U(J)
      IF(DUNI.LT.0.0D0)DUNI = DUNI+1.0D0
      U(I) = DUNI
      I = I-1
      IF(I.EQ.0)I = 17
      J = J-1
      IF(J.EQ.0)J = 17
      C = C-CD
      IF(C.LT.0.0D0) C=C+CM
      DUNI = DUNI-C 
      IF(DUNI.LT.0.0D0)DUNI = DUNI+1.0D0
      RETURN
      ENTRY DUSTAR(ISEED)
        I1 = MOD(ABS(ISEED),177)+1
        J1 = MOD(ABS(ISEED),167)+1
        K1 = MOD(ABS(ISEED),157)+1
        L1 = MOD(ABS(ISEED),147)+1
        DO 2 II = 1,17
          S = 0.0D0 
          T = 0.5D0 
          DO 3 JJ = 1,K
                  M1 = MOD(MOD(I1*J1,179)*K1,179) 
                  I1 = J1
                  J1 = K1
                  K1 = M1
                  L1 = MOD(53*L1+1,169) 
                  IF(MOD(L1*M1,64).GE.32)S=S+T
    3             T = 0.5D0*T 
    2   U(II) = S
        DUSTAR = FLOAT(ISEED) 
        RETURN
      ENTRY DUNIB(KK)
        IF(KK.LE.47)THEN
             K=47
        ELSE
             K=KK
        ENDIF
        DUNIB=FLOAT(K)
      END