inst/examples/powerLawErgodicity/functions.R

In msuzen/isingLenzMC: Monte Carlo for Classical Ising Model

#
#  Ergodic Dynamics of Ising Model 
#  R simulation functions, data generation, analysis and other utilities
#  (c) 2013, 2014, 2015, 2016 by Dr.Mehmet Suzen
#  GPLv3 or higher
#

require("isingLenzMC"); # 0.2.3

#
# Using Analytic solution for ensemble magnetisation 
# using transfer matrix
#
#  Two state ising model: general case
#  http://micro.stanford.edu/~caiwei/me334/
#  Ch 7 page 35
# 
#Z(h, B, J, N) := exp(N*B*J) * ((cosh(B*h) + sqrt(sinh(B*h)^2+exp(-4*B*J)))^N + (cosh(B*h) - sqrt(sinh(B*h)^2+exp(-4*B*J)))^N);

# now the magnetisation (mulitply with 1/B) compute with maxima
# very long!
#diff(log(Z(h, B, J, N)), h);
#
#(%i4) diff(log(Z(h, B, J, N)), h);
#           B cosh(h B) sinh(h B)
#(%o4) ((---------------------------- + B sinh(h B))
#               - 4 B J       2
#        sqrt(%e        + sinh (h B))
#         - 4 B J       2                   N - 1
# (sqrt(%e        + sinh (h B)) + cosh(h B))      N
#                     B cosh(h B) sinh(h B)
# + (B sinh(h B) - ----------------------------)
#                         - 4 B J       2
#                  sqrt(%e        + sinh (h B))
#                     - 4 B J       2       N - 1
# (cosh(h B) - sqrt(%e        + sinh (h B)))      N)
#          - 4 B J       2                   N
#/((sqrt(%e        + sinh (h B)) + cosh(h B))
#                       - 4 B J       2       N
# + (cosh(h B) - sqrt(%e        + sinh (h B))) )

getEnsembleMagnetisation <- function(n, ikBT, J, H) {
# Derivative of Equation 35 of Wu notes at Standford
# Maxima output above
 if(n < 700) { # This is ad-hoc assumption is that we recover n -> inf solution above this
   cHB   <- cosh(H*ikBT)
   sHB   <- sinh(H*ikBT)
   e4B   <- exp(-4*ikBT*J)
   sq4B  <- sqrt(e4B+sHB^2)
   isq4B <- ikBT*cHB*sHB/sq4B
   # Per site ensemble averaged magnetisation
   eMag  <- isq4B + ikBT*sHB
   eMag  <- n * eMag * (sq4B+cHB)^(n-1)
   eMag  <- eMag + n*(ikBT*sHB - isq4B) * (cHB - sq4B)^(n-1)
   eMag  <- eMag * ((sq4B+cHB)^(n) + (cHB-sq4B)^(n))^(-1)
   if(is.nan(eMag) || is.infinite(eMag)) return(0.9934344) # 0.9934344 infinite solution
   return(eMag/ikBT/n)
  } else {
   return(0.9934344)
 }
}

#
# This is to perform 1D Ising simulation and Getting TM metric
# Thirumalai - Mountain Fluctuation Metric
#
magnetisationMetricTM1D <- function(ikBT, J, H, n, nstep, transP) {
  ensembleM <- getEnsembleMagnetisation(n, ikBT, J, H) # ensemble average magnetisation per-site
  x         <- genConfig1D(n)
  mySim     <- isPerform1D(ikBT, x, J, H, nstep, ensembleM, transP) 
  time      <- 1:mySim$naccept
  # Note that, Metropolis variable names! could be glauber as well depending transP value (1 or 2)
  magnetisationMetric <- mySim$omegaM[time]
  metropolisTime      <- mySim$times[time]
  out                 <- list(magnetisationMetric=magnetisationMetric, metropolisTime=metropolisTime);
  out
}


#
# Running a single sim (time vs. magnetisation metric and 
# generate data file with the parametrisation as fileName
#
runSim <- function(ikBT, J, H, n, nstep, transP) {
    if(transP == 1) dynamicsName <- 'metropolis';
    if(transP == 2) dynamicsName <- 'glauber';
    out          <- magnetisationMetricTM1D(ikBT, J, H, n, nstep, transP)
    # Note that, Metropolis variable names! could be glauber as well depending transP value (1 or 2)
    ll           <- 2*length(out$metropolisTime) 
    timeMag      <- matrix(1:ll, ncol= 2, nrow=length(out$metropolisTime))
    timeMag[,1]  <- out$metropolisTime
    timeMag[,2]  <- out$magnetisationMetric
    timeMag
}