2dexternalhising.f90

考虑了外加磁场的ISING MODEL 模拟 要用fortran软件包打开

 !研究2D ISING MODEL程序模拟中，有外场情况下，并外场不断改变时，
 !对各种物理量可能带来的影响
 !观察各量在取定外场，随温度变化的发展规律，提取临界点，并总结不同外场下的对比规律
 !
 !
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,包含外场信息的比例系数
  !!!!!!!!!!!!!!!!!!!!!!!
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')
 !取定温度的初始位型，温度在临界点上或下
 open(unit=13,file='lastconfig1.dat',status='replace',action='write')
 open(unit=133,file='lastconfig2.dat',status='replace',action='write')
 !取定温度的最终位型，温度在临界点上或下
 open(unit=14,file='magnetTime1.dat',status='replace',action='write') 
 open(unit=144,file='magnetTime2.dat',status='replace',action='write') 
 !取定温度，磁化强度与时间（MCS）的关系，温度在临界点上或下
 open(unit=15,file='magnetTemp.dat',status='replace',action='write') 
 !单个自旋的磁化强度与温度的函数关系(绝对值也输出到该文件)
 open(unit=16,file='susceptibility.dat',status='replace',action='write')	
 !单个自旋的极化率与温度的关系
 open(unit=17,file='energyTemp.dat',status='replace',action='write')	
 !单个自旋的能量与温度的函数关系
 open(unit=18,file='heatcapacity.dat',status='replace',action='write')
 !单个自旋的比热容与温度的函数关系
 
   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	!调用外部函数
 double precision::duni		 !函数关键词声明，数据类型要一致
 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
  !这些量在不同温度进行清零，相同温度下，进行累加
  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
	 !平衡后，隔步抽样，减小关联效应
	   stat_mcs=stat_mcs+1
	   magnet=sum(s(2:nrows+1,2:ncols+1))/(nrows*ncols*1.0)	!平均到单个自旋上
	   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
	   !统计临界温度上和临界温度下平均磁化强度与时间的关系
	   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 !考虑外场情况，平均到单个自旋上
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	   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)!做粒子的尝试翻转
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        deltU=-try_spin*(s(n-1,m)+s(n+1,m)+s(n,m-1)+s(n,m+1))*2-exB*try_spin*2
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!考虑外场情况,再减去自旋翻转后，磁场强度的变化量
		!计算出因为该粒子的自旋翻转而引起的前后能量变化
		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
		  !周期性边界条件，因为是最近四邻相互作用，所以要保证
          !第一行与倒数第二行相同分布
          !最后一行与第二行相同分布
          !第一列与倒数第二列相同分布
          !最后一列与第二列相同分布
        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抽样统计平均值的计算
  !为方便对比不同尺寸，已经平均到单个自旋上
  magnet_ave=(magnet_ave/stat_mcs); !平均磁化强度计算
  magnet2_ave=magnet2_ave/stat_mcs; !磁化强度的方均值
  susceptibility=(magnet2_ave-(magnet_ave)**2)/T; !磁化率

  energy_ave=energy_ave/stat_mcs; !平均能量计算
  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