Ising Model Starter

program isingmodel

! JS: NOTE:  There is a book on reserve in the science library,
! JS:  David Chandler, "Introduction to Modern Statistical Mechanics"
! JS:  that describes both the Ising model and the Monte Carlo procedure
! JS:  in some depth.  It might be a useful resource if you get stuck.

! array, defining variables

implicit none
real :: kT, E, E2, M, C, paccept,ptrial, Estart,Eflip

integer :: col, row, east,west,north,south
integer :: nMC
integer, parameter :: J = -1    ! define a ferromagnetic system
integer, dimension(16,16) :: spin 

! random number generator
call random_seed()

! initialized our spin array
do col=1,16
   do row=1,16
      call random_number(paccept)
      if (paccept .lt. 0.5) then
         spin(col,row) = +1
      else 
         spin(col,row) = -1
      end if
   end do
end do

! program structure:

! defining the loops (and what do we mean by loops)

! loop over temperature (kT in units of the coupling, J)

do kT=4.0, 0.1, -0.05      !JS: Start at kT=4, and go to kT=0.1 in steps of -0.05

! JS: First, run a set of MC steps to equilibrate our system. 
! JS: In other words, so that the configuration of the spins is
! JS: in an equilibrium configuration consistent with Boltzman's laws
! JS: for the particular temperature 

! loop MC trials
   do nMC=1,10000

! spatial loop in the 2D plane (incorporate pbc)

      do col=1,16
         do row=1,16

            east = col+1
            west = col-1
            north = row-1
            south = row+1

            if (col.eq.1) west = 16
            if (col.eq.16) east = 1
            if (row.eq.1) north = 16
            if (row.eq.16) south = 1

! try flip a spin at position (col,row)
! calculate energy (function?)

            Estart = J*(spin(col,row)*spin(east,row) &
                 + spin(col,row)*spin(west,row) &
                 + spin(col,row)*spin(col,north) &
                 + spin(col,row)*spin(col,south) )

            Eflip =  J*((-spin(col,row))*spin(east,row) &
                 + (-spin(col,row))*spin(west,row) &
                 + (-spin(col,row))*spin(col,north) &
                 + (-spin(col,row))*spin(col,south) )

! acceptance/rejectance for MC algorithm (define as subroutine?)

            if (Eflip .lt. Estart) then 
               spin(col,row) = -spin(col,row)
            else  
               call random_number(paccept)
               ptrial = EXP(-(Eflip-Estart)/kT)

               if (paccept .lt. ptrial) then
                  spin(col,row) = -spin(col,row)
               end if
            end if

         end do
      end do

   end do  ! loop over nMC

! JS: Now that the above has finished, and we are in equilibrium,
! JS: run another MC cycle again, this time to collect statistics
! JS:  In other words, another do-loop over nMC, row, col, etc.
! JS:  flipping spins as appropriate, etc.

! JS: BUT now, every time (mod(nMC, 50) .eq. 0), collect E, E^2, M to calculate
! JS: averages.  This E will have to be the total energy of the entire system
! JS: calculated over all the spins.  You might want to create functions to 
! JS: give you the value of E, M, E^2, or you could just do it inline.

! JS: When the datacollection is over (you have completed this second run of
! JS: nMC steps), it will be time to output the averages.

!OUTPUT:  Graphs; (Averages) Heat Capacity (<E^2>-<E>^2, ;derivative of <E(T)>, 
! Temperature, Energy <E> (to get Heat Capacity), Magnetization <M>
! JS: (You will want to report these as an intensive property,  
! JS:  i.e., <E>/N, C/N, <M>/N, where N is the total number of spins (16*16)

end do     !JS: loop over kT values (slowly cooling the system by lowering kT)

end program isingmodel