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```
```