The Ising Model

This is an example taken from the Fortran 90 standard, pp 296-298; the code has been modified to make it a main program instead of a subroutine.

The code is a 3-D Monte Carlo simulation of state transition. Each gridpoint is a logical variable whose value can be interpreted as spin-up or spin-down. The transition between states is governed by a local probabilistic process where all points change state at the same time. Each spin either flips to the opposite state or not depending on the state of its six nearest neighbours. Gridpoints on the edge of the cube are defined by cubic periodicity -- in other words the grid is taken to be replicated in all dimensions in space.

The code (which is available by clicking here) given below will perform the desired transitions. Add HPF directives to allow the problem to be solved with minimal communications:

      MODULE Funkt
      CONTAINS
       FUNCTION RAND (m)
        INTEGER m
        REAL, DIMENSION(m,m,m) :: RAND
         CALL RANDOM_NUMBER(HARVEST = RAND)
        RETURN
       END FUNCTION RAND
      END MODULE Funkt

      PROGRAM TRANSITION
      USE Funkt
      IMPLICIT NONE
       INTEGER, PARAMETER :: n = 16
       INTEGER            :: iterations, i
       LOGICAL, DIMENSION(n,n,n) :: ising,     flips
       INTEGER, DIMENSION(n,n,n) :: ones,      count
       REAL, DIMENSION(n,n,n)    :: threshold
       REAL, DIMENSION(6)        :: p

       p = (/ 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 /)

       iterations  = 10
       ising = RAND(n) .LE. 0.5

       DO i = 1,iterations
        ones = 0
        WHERE (ising) ones = 1
        count = CSHIFT(ones, -1, 1) + CSHIFT(ones, 1, 1) &
              + CSHIFT(ones, -1, 2) + CSHIFT(ones, 1, 2) &
              + CSHIFT(ones, -1, 3) + CSHIFT(ones, 1, 3)
        WHERE (.NOT.ising) count = 6 - count
        threshold = 1.0
        WHERE (count == 4) threshold = p(4)
        WHERE (count == 5) threshold = p(5)
        WHERE (count == 6) threshold = p(6)
        flips = RAND(n) .LE. threshold
        WHERE (flips) ising = .NOT. ising
       ENDDO
      END PROGRAM TRANSITION

Note that CSHIFT is used in SOR problems and performs a circular shift on an array, for example, if

then

   CSHIFT(A,-1)

is A shifted one place to the left with the left-most number wrapping around to the right,

and is A shifted one place to the right

   CSHIFT(A,1)

is

It is also possible to specify a dimension for 2D and upward arrays. If

then

   CSHIFT(B,1,1)

shifts the array one position in dimension 1 (downwards)

and

   CSHIFT(B,1,2)

and so on.

• Solution