2d Ising Model

hw4prob1.r

# hw4prob1.r
# 2d ising model

hamiltonian <- function (E, J, sig, N) {
  nearest_neighbor_sum = sum(sig[1:N-1,]*sig[2:N,], sig[,1:N-1]*sig[,2:N])
  return(-E * nearest_neighbor_sum - J * sum(sig))
}

ising <- function (E, J, N, beta, iterations) {
  number = N*N
  sig = matrix(0, N, N)

  # Set our baseline.
  for (i in 1:N) {
    for (j in 1:N) {
      u = runif(1)

      if (u < 0.5) {
        sig[i,j] = -1
      } else {
        sig[i,j] = 1
      }
    }
  }

  H = hamiltonian(E, J, sig, N)
  for (i in 1:iterations) {
    print(i)
    x = ceiling(runif(1)*N)
    y = ceiling(runif(1)*N)

    sigt = sig
    sigt[x, y] = -sig[x, y]
    new_neighbor_score = sum(sigt[x,y]*sigt[x-1,y], sigt[x,y]*sigt[x,y-1])

    if (x != N) {
      new_neighbor_score = new_neighbor_score + sigt[x,y]*sigt[x+1,y]
    }

    if (y != N) {
      new_neighbor_score = new_neighbor_score + sigt[x,y]*sigt[x,y+1]
    }
    tmpH = H[i] - 2*E*new_neighbor_score

    if (tmpH - H[i] < 0) {
      sig = sigt
      H[i+1] = tmpH
    } else if (runif(1) < exp(-beta*(tmpH-H[i]))) {#exp(-beta * (H2 - H1))) {
      sig = sigt
      H[i+1] = tmpH
    } else {
      H[i+1] = H[i]
    }
  }
  # image(sig)
  plot((1:(iterations+1)), H)

}