PROGRAM ISING

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!  ISING MONTE CARLO PROGRAM   !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
  
! Monte Carlo for the D=2 Ising model, 
! with beta * H = - JJ * \sum_{<ij>} s_i * s_j ,           
! where s_i = +1 or -1                                    
! [critical coupling is JJ_c = .4406867935...]

  use isinglib

  implicit none
  
  integer, parameter :: prec = selected_real_kind(15,32)

! specify size of the periodic array (may want to change this!):
!
  integer, parameter :: i_max = 30, j_max = 30
  
  integer :: i, j, k, m, ip, jp, im, jm, temp_spin, N_init, N_sweep, &
     spin(i_max,j_max), temp_mag(j_max), temp_ss(j_max)
  
  real (kind=prec):: JJ, expdeltaS, r, ss(j_max), ss_conn(j_max) 
  
  write(6,*) "enter random seed, JJ, N_init, N_sweep"
  read(5,*) ran_seed 
  read(5,*) JJ 
  read(5,*) N_init 
  read(5,*) N_sweep 
  write(6,*) "working"
!
! "hot start"
! 
  do i=1, i_max
    do j=1, j_max
      if(random() < 0.5) then 
        spin(i,j) = 1
      else 
        spin(i,j) = -1
      endif
    enddo
  enddo         
!
! initialization and measurement sweeps:
!
! first initialize the magnetization and spin-spin correlation 
! variables:
!
do m=1, j_max
  temp_mag(m) = 0
  temp_ss(m) = 0
enddo  
!
do k = 1, N_init + N_sweep
  do i=1, i_max
    do j=1, j_max
! generate a temporary spin[i,j], temp_spin: 
      if(random() < 0.5) then 
        temp_spin = 1
      else 
        temp_spin = -1
      endif
! nearest-neighbor coordinates:
      if(i==1) then 
        im=i_max
        ip=i+1
      elseif(i==i_max) then
        im=i-1
        ip=1
      else
        im=i-1
        ip=i+1
      endif        
!
      if(j==1) then 
        jm=j_max
        jp=j+1
      elseif(j==j_max) then
        jm=j-1
        jp=1
      else
        jm=j-1
        jp=j+1
      endif        
!
! calculate exp(-[change in Action]):
!  
      expdeltaS = exp( JJ * (temp_spin-spin(i,j))  &
         * ( spin(i,jm) + spin(i,jp) + spin(im,j) + spin(ip,j) ) )
!
! compare to new random number r; if r > expdeltaS we keep the 
! old configuration; otherwise we keep the new one:
      r = random()
      if(r > expdeltaS) then 
      else 
        spin(i,j) = temp_spin
      endif
    enddo     
  enddo
! if k > N_init, make the measurements:
  if(k > N_init) then
    do m=1, j_max 
      temp_mag(m) = temp_mag(m) + spin(1,m)   
      temp_ss(m) = temp_ss(m) + spin(1,1)*spin(1,m)  
    enddo
  endif
enddo
!
! final results:
!
do m=1, j_max 
  ss(m) = real(temp_ss(m))/N_sweep 
  ss_conn(m) = ss(m) &
     - real(temp_mag(1))/N_sweep * real(temp_mag(m))/N_sweep 
enddo
!
  write(6,*) "spin-spin correlation function along a row"
do m=1, j_max 
  print*,ss(m)  
enddo
  write(6,*) "connected spin-spin correlation function along the row"
do m=1, j_max 
  print*,ss_conn(m)  
enddo

END PROGRAM ISING

