Nothing
#
# Ising Model Dynamics MC
# R wrapper functions to C library functions
# (c) 2013 by Dr.Mehmet Suzen
# GPLv3 or higher
#
# Generate Random Configuration in 1D
# Arguments
# n : number of sites
genConfig1D <- function(n) {
x <- rep(0.0,n)
out <- .C("genConfig1D", n=as.integer(n), x=as.double(x));
out$x
}
# Random single flip configuration in 1D
# Arguments
# x : sites vector
# Return random flip single site
flipConfig1D <- function(x) {
out <- .C("flipConfig1D", n=as.integer(length(x)), x=as.double(x));
out$x
}
# Random many flip configuration in 1D
# Arguments
# x : sites
# upperF : Number of repeats for random flip
flipConfig1Dmany <- function(x, upperF) {
out <- .C("flipConfig1Dmany", n=as.integer(length(x)), x=as.double(x), upperF=as.integer(upperF));
out$x
}
# Nearest-Neighbour energy in periodic boundary conditions in 1D
# Arguments
# ll : number of sites
# vec : site spin configuration
# energy : energy to be returned
lattice1DenergyNN <- function(x) {
energy <-0
if(!is.numeric(x))
stop("argument x must be numeric");
out <- .C("lattice1DenergyNN", n=as.integer(length(x)), vec=as.double(x), energy=as.double(energy));
out$energy
}
#
# Sum given vector
# Arguments
# x : site spin configuration
sumVec <- function(x) {
sum <-0
if(!is.numeric(x))
stop("argument x must be numeric");
out <- .C("sumVec", n=as.integer(length(x)), vec=as.double(x), sum=as.double(sum));
out$sum
}
# Total Energy in 1D
# x : site spin configuration
# J : interaction strength
# H : external field
totalEnergy1D <- function(x, J, H) {
energy <- 0
if(!is.numeric(x))
stop("argument x must be numeric");
out <- .C("totalEnergy1D", n=as.integer(length(x)), vec=as.double(x), J=as.double(J), H=as.double(H), totalEnergy=as.double(energy));
out$totalEnergy
}
# Compute transition probability
# Arguments
# ikBT : 1/kB*T (boltzmann factor)
# x : 1D Spin sites on the lattice
# xflip: 1D Spin sites on the lattice: after a flip
# J : interaction strength
# H : external field
# probSel : which transition probability to use 1 Metropolis 2 Glauber
transitionProbability1D <- function(ikBT, x, xflip, J, H, probSel) {
prob <- 0.0;
if(!is.numeric(x))
stop("argument x must be numeric");
if(!is.numeric(xflip))
stop("argument xflip must be numeric");
out <- .C("transitionProbability1D", ikBT=as.double(ikBT), n=as.integer(length(x)), vec=as.double(x), vecFlip=as.double(xflip),
J=as.double(J), H=as.double(H), prob=as.double(prob), probSel=as.integer(probSel));
out$prob
}
# Take One MC Step : importance sampling!
# Author : msuzen
# Arguments
# ikBT : 1 over Temperature times Boltzmann constant
# x : 1D Spin sites on the lattice
# J : interaction strength
# H : external field
# rout : Pair list, flip states (vec) and if step is accepted (accept)
# probSel : which transition probability to use 1 Metropolis 2 Glauber
isStep1D <- function(ikBT, x, J, H, probSel) {
prob <- 0.0;
accept <- 0.0;
if(!is.numeric(x))
stop("argument x must be numeric");
out <- .C("isStep1D", ikBT=as.double(ikBT), n=as.integer(length(x)), vec=as.double(x),
J=as.double(J), H=as.double(H), prob=as.double(prob), accept=as.integer(accept),
probSel=as.integer(probSel));
rout <- list(vec=out$vec, accept=out$accept)
rout
}
# Perform importance sampling MC on 1D
#
# Author : msuzen
# Arguments
# ikBT : 1 over Temperature times Boltzmann constant
# x : 1D Spin sites on the lattice
# J : interaction strength
# H : external field
# nstep : number of MC steps requested
# ensembleM : Value of the theoretical magnetization (could be thermodynamic limit value)
# Output pair list:
# omegaM : Fluctuating metric vector for Magnetisation (length of naccept)
# naccept : number of MC steps accepted
# nreject : number of MC steps rejected
# probSel : which transition probability to use 1 Metropolis 2 Glauber
#
isPerform1D <- function(ikBT, x, J, H, nstep, ensembleM, probSel) {
if(!is.numeric(x))
stop("argument x must be numeric");
omegaM <- rep(0.0, nstep); # Let R allocate memory; assume all steps will be accepted
# index 1 to naccept will be reported in the vector
times <- rep(1.0, nstep); # Let R allocate memory: this is the so called time to accepted flip!
naccept <- 0;
nreject <- 0;
out <- .C("isPerform1D", ikBT=as.double(ikBT), n=as.integer(length(x)), vec=as.double(x),
J=as.double(J), H=as.double(H), ensembleM=as.double(ensembleM),
omegaM=as.double(omegaM), nstep=as.integer(nstep), naccept=as.integer(naccept),
nreject=as.integer(nreject), times=as.integer(times), probSel=as.integer(probSel));
rout <- list(omegaM=out$omegaM, nreject=out$nreject, naccept=out$naccept, times=out$times);
rout
}
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.