Nothing
#
# 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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.