jinlin2.f90

ISING model的monte carlo 不同温度模拟程序

  module global
  integer::m,n,q,w
  integer,allocatable::A(:,:)
  integer::nrows,ncols
  integer::npass
  integer::ipass
  integer::nequil
  integer::trial_spin
  integer::nscans
  integer::iscan
  integer::update
  real::i,j
  real::temp,beta
  real::deltaU
  real::high_temp
  real::low_temp
  real::temp_interval

  real::log_eta
  real::magnetization
  real::magnetization_ave
  real::magnetization2_ave
  real::energy
  real::energy_ave
  real::energy2_ave
  integer::output_count,sameT_mcs_times
  real::seed1,seed2
  character(10)::t1,t2
end module global
  
  program ising
    use global
	implicit none
	call date_and_time(time=t1)
    call initializeParameter
	!!Attention the difference temperature
	nscans=int((high_temp-low_temp)/temp_interval)+1
 scan_loop:do iscan=1,nscans
      seed1=1
	  seed2=w+1
      output_count=0
      energy_ave=0.0
      energy2_ave=0.0
      magnetization_ave=0.0
      magnetization2_ave=0.0
      temp=high_temp-temp_interval*(iscan-1)
      beta=1.0/temp
	  call initializeConfiguration
      call Evolve_ising
      call output
enddo scan_loop
    close(31)
    close(32)
	close(33)
	close(34)
    close(35)
	close(36)
	close(37)
	close(38)
	close(40)
	close(41)
	call date_and_time(time=t2)
	write(39,*)"time=",t2,"-",t1 
	close(39)
	deallocate(A)
 end program ising

subroutine initializeParameter
  use global
  implicit none
write(*,*)"nrows="
  read(*,*) nrows
write(*,*)"ncols="
  read(*,*) ncols
write(*,*)"npass="
   read(*,*) npass
write(*,*)"nequil="
   read(*,*) nequil
write(*,*)"high_temp="
  read(*,*) high_temp
write(*,*)"low_temp="
  read(*,*) low_temp
write(*,*)"temp_interval="
  read(*,*) temp_interval

   allocate(A(nrows+2,ncols+2))
   update=0
   !one_mcs loop times
        w=nrows*ncols
   !!
open(unit=31,file='lastspin-array2.dat',status='replace',action='write')
  write(31,*) '          nrows      ncols      spin'       
open(unit=32,file='lastspin-array.dat',status='replace',action='write')
  write(32,*) '          nrows      ncols      spin' 
open(unit=33,file='magnetization.dat',status='replace',action='write')
  write(33,*)	"     temp    ave_magnetization  ave_magnetization^2   susceptibility"
open(unit=34,file='energy.dat',status='replace',action='write')
  write(34,*) "       temp    ave_energy         ave_energy^2       C_v"
open(unit=35,file='sameTenergy.dat',status='replace',action='write')
  write(35,*)" sameT_mcs_times  Totalenergy  energy_ave/output_count energy2_ave/output_count"
open(unit=36,file='startspin-array.dat',status='replace',action='write')
  write(36,*)'           i       j      spin'
  open(unit=37,file='startspin-array2.dat',status='replace',action='write')
  write(37,*)'           i       j      spin'

open(unit=38,file='random and seeds.dat',status='replace',action='write')
  write(38,*)"        m         n        seed3      -deltaU      log_eta"
open(unit=39,file='parameter.dat',status='replace',action='write')
  write(39,*)"nrows=",nrows,    "      ncols=",ncols   
  write(39,*)"MCS-npass=",npass,"     nequil=",nequil
  write(39,*)"high_temp=" ,high_temp,"   low_temp=",low_temp,"   temp_interval=",temp_interval
  write(39,*)"  time=t2-t1 "
open(unit=40,file='sameTenergy2.dat',status='replace',action='write')
  write(40,*)" sameT_mcs_times  Totalenergy  energy_ave/output_count energy2_ave/output_count"
open(unit=41,file='sameTenergy3.dat',status='replace',action='write')
  write(41,*)" sameT_mcs_times  Totalenergy  energy_ave/output_count energy2_ave/output_count"

end subroutine initializeParameter

subroutine initializeConfiguration
use global
implicit none
external rand1
real rand1
	  do i=2,nrows+1
	    do j=2,ncols-1
		 A(i,j)=1
		enddo
		do j=ncols,ncols+1
		 A(i,j)=-1    
        enddo
	  enddo
	 do i=2,nrows+1
	    A(i,1)=A(i,nrows+1)
		A(i,nrows+2)=A(i,2)
	 enddo
	 do j=2,ncols+1
	    A(1,j)=A(ncols+1,j)
		A(ncols+2,j)=A(2,j)
	 enddo
	 
	 write(36,*)"temp=",temp
	  do i=2,nrows+1
	    do j=2,ncols+1
		 if(A(i,j)==1) write(36,*)i,j,A(i,j)
		enddo
	  enddo

	  write(37,*)"temp=",temp
	  do i=2,nrows+1
	    do j=2,ncols+1
		 if(A(i,j)==-1) write(37,*)i,j,A(i,j)
		enddo
	  enddo

end subroutine initializeConfiguration

subroutine Evolve_ising
   use global
   implicit none
   external rand1
 real::rand1
MC_passes: do ipass=0,npass
 update=update+1
 if (mod(update,10)==0) then !!消除样本之间可能的关联
 if (ipass>nequil) then
       output_count=output_count+1
       magnetization=sum(A(2:nrows+1,2:ncols+1))/(ncols*nrows*1.0)
       magnetization_ave=magnetization_ave+magnetization
       magnetization2_ave=magnetization2_ave+magnetization**2
       energy=0.0
  do i=2,nrows+1
    do j=2,ncols+1
      energy=energy-A(i,j)*(A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1))
    enddo
  enddo
!output the difference T energy information
!T=high_temp-temp_interval*(iscan-1)
    if (iscan==1)then
	  sameT_mcs_times=ipass
      write(35,*)sameT_mcs_times, energy,   energy_ave/output_count,   energy2_ave/output_count
	endif
    if (iscan==2)then
	  sameT_mcs_times=ipass
      write(40,*)sameT_mcs_times, energy,   energy_ave/output_count,   energy2_ave/output_count
	endif

    if (iscan==7)then
	  sameT_mcs_times=ipass
      write(41,*)sameT_mcs_times, energy,   energy_ave/output_count,   energy2_ave/output_count
	endif
!!
   energy=energy/(nrows*ncols*2.0)
   energy_ave=energy_ave+energy
   energy2_ave=energy2_ave+energy**2
  endif
endif

one_mcs: do m=2,nrows+1
           do n=2,ncols+1
            seed2=seed2+1
         trial_spin=-A(m,n)
         deltaU=-trial_spin*(A(m-1,n)+A(m+1,n)+A(m,n-1)+A(m,n+1))*2
         log_eta=log(rand1(seed2)+1.0d-10)
        if (-beta*deltaU > log_eta) then
           A(m,n)=trial_spin
           if (m==2)       A(nrows+2,n)=trial_spin
           if (m==nrows+1) A(1,n)=trial_spin
           if (n==2)       A(m,ncols+2)=trial_spin
           if (n==ncols+1) A(m,1)=trial_spin
       endif
	enddo
   enddo one_mcs
enddo MC_passes

end subroutine Evolve_ising


subroutine output
  use global
  implicit none
    
 write(32,*)"temp=",temp 
  do i=2,nrows+1
    do j=2,ncols+1
    if (A(i,j)==1) write(32,*)  i,  j,   A(i,j)
    enddo
  enddo

  write(31,*)"temp=",temp 
  do i=2,nrows+1
    do j=2,ncols+1
    if (A(i,j)==-1) write(31,*)  i,  j,   A(i,j)
    enddo
  enddo

  !!output every T's parameter amount
 write(33,*) temp,    abs(magnetization_ave/output_count),  magnetization2_ave/output_count,&
   beta*(magnetization2_ave/output_count-(magnetization_ave/output_count)**2)
 write(34,*) temp,    energy_ave/output_count,   energy2_ave/output_count,&
   (beta**2)*(energy2_ave/output_count-(energy_ave/output_count)**2)
 
end subroutine output

!random number generator
real function rand1(seed)
   implicit none
   integer a,m,q,p,n,ndiv,j,k
   real rm,rmax,seed
   parameter (a = 16807, m = 2147483647, rm = 1.0/m)
   parameter (q = 127773, p = 2836, n = 32, ndiv = 1 + (m-1)/n)
   parameter (rmax = 1.0 - 1.2e-7)
   integer r(n),r0,r1
   save r,r0,r1
   data r/n*0/
      if(seed .ne. 0) then
        r1 = abs(seed)
           do j = n+8,1,-1
             k = r1/q
             r1 = a*(r1-k*q) - p*k
             if (r1 .lt. 0) r1 = r1 + m
             if (j .le. n) r(j) = r1
           enddo
        r0 = r(1)
      end if
      k = r1/q
      r1 = a*(r1 - k*q) - p*k
      if (r1 .lt. 0) r1 = r1 + m
      j = 1 + r0/ndiv
      r0 = r(j)
      r(j) = r1
      rand1 = min(rm*r0,rmax)
  end function rand1