2dhysteresisloop.f90

ISING MODEL模拟磁滞回线

 !研究2D ISING MODEL程序模拟中，临界温度下，磁场强度随外场增加再减小的规律，
 !顺便也考察其他量随外加磁场变换的情况。
module global	 !!常量声明
  integer,parameter:: nrows=20,ncols=20	!2Dsize	若改变尺寸，只需要改变这里就可以
  integer,parameter::Tmcs=80000,Emcs=50000	 !Total monte carlo steps
  real,parameter::exHB=10.0,exLB=-5.0,intervalB=0.1 !high temperature
  real,parameter::T=3.0 !取定温度，可以改变温度
end module global

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

 open(unit=11,file='parameter.dat',status='replace',action='write')
   write(11,'(A)')'all the parameters of 2Dhysteresisloop 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,'T=',T
   write(11,'(A,f5.2,2X,A,f5.2,2X,A,f5.2)')'exHB=',exHB,'exLB=',exLB,'intervalB=',intervalB
 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='magnetexB.dat',status='replace',action='write') 
 !单个自旋的磁化强度与温度的函数关系(绝对值也输出到该文件)
 open(unit=16,file='susceptibility.dat',status='replace',action='write')	
 !单个自旋的极化率与温度的关系
 open(unit=17,file='energyexB.dat',status='replace',action='write')	
 !单个自旋的能量与温度的函数关系
 open(unit=18,file='heatcapacity.dat',status='replace',action='write')
 !单个自旋的比热容与温度的函数关系
 
   write(*,*)'2D Ising model in different external B with one temperature  '
   write(*,'(2x,a,i4,2x,a,i4,2x,a,f5.2)')'nrows=',nrows,'ncols=',ncols,'T=',T
   call date_and_time(time=t1)
   do i=1,60  !100  NT=int((exHB-exLB)/intervalB)
	 ! if(i.le.51) exB=exLB+intervalB*(i-1)
	 ! if(i.gt.51) exB=exHB-intervalB*(i-51)
	 !0---Hs,Hs--0,0--(-Hs),(-Hs)--0--Hs
	  exB=-intervalB*(i-1)
	  write(*,'(2x,a,f6.3)')'exB=',exB
	  call initial(exB,s)
	  call isingMC(exB,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  isinghysteresisloop

subroutine initial(exB,s)	 !!初始位型,cold start
 use global
 implicit none
 real::exB
 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(exB.eq.3.0)write(12,'(i3,2x,i3,2x,i3)')i,j,s(i,j)
		 !只输出一次就可以了
	  enddo 
   enddo
endsubroutine initial

subroutine isingMC(exB,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::exB,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(exB*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(exB.eq.-3.0)write(14,'(i5,2x,f8.3)')imcs,magnet
	   if(exB.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((exB.eq.-3.0).and.(s(i,j).eq.1))write(13,'(i3,2x,i3,2x,i2)')i,j,s(i,j)
		 if((exB.eq.3.0).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)')exB,magnet_ave,abs(magnet_ave)
  write(16,'(f5.2,2x,f8.3)')exB, susceptibility
  write(17,'(f5.2,2x,f8.3)')exB, energy_ave
  write(18,'(f5.2,2x,f12.6)')exB, 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