Monte Carlo Simulation

1
Monte Carlo Simulation

• Ken Hawick
• 2009

2
Monte Carlo Methods

• Monte Carlo Casino
• Relates to use of Random Numbers
• MC or sometimes MCMC (Markov Chain Monte Carlo)
  methods used in variety of simulations

3
Models of Phase Transitions

• Ising Model example model of a magnet
• Thermal Heatbath and use of Random Noise
• Other systems such as networks such as the
  Internet social networks small-world networks
  distorted lattices

4
Phase Transitions Complex Networks

• A phase transition is a sudden or sometimes
  unexpected change in some measurable or
  observable property
• Usually cannot be analysed with simple
  mathematics
• Complex network - a graph of interacting
  components, such as agents
• A general way to try to quantify complexity

5
Monte Carlo Method

• Impose a Markov dynamics on the Ising Hamiltonian
  (energy function)
• ie construct a path through phase space starting
  from a random point
• Ising phase space has 2N points (too huge, so
  can only sample)
• Transition probabilities based on Boltzman
  probability exp( -?E / kT)

• Metropolis approach updates a single spin at a
  time
• Slow to converge - we say dynamical exponent z is
  large
• Time T-Tc z
• But does give realistic growth or quenching
  behaviour

6
Monte Carlo Timescales

• Time in the simulation is Markov time
• We can (cautiously) ascribe some universal
  dynamical properties to it
• Key idea is that length-scales and correlation
  times are important

• L x L System is capable of supporting
  length-scales up to L if non periodic
• If apply periodicity, only L/2 capable of being
  supported
• This shows up as a finite size effect

7
Lattice Models and Graphs

• Spin Models - the Ising Model
• Monte-Carlo approach to dynamics
• Wormholes and graph shortcuts
• Small-World Graphs and properties
• Probing the small-world regime
• Technicalities - Fast Simulations
• Lattice and Graph Regimes to study

8
The Ising Model

• Model of magnet
• Spins - up or down on each site of a mesh
• Energy(Hamiltonian) relates how they are coupled
• Like-Like or Like-Unlike etc
• Order varies with Temperature

9
Programming Ising Systems

• Array (2d or 3d) of spins (bits or bytes or
  ints)
• Need to know nearest neighbours
• Metropolis Algorithm
• Pick a random site (one out of the N)
• Count many like-like nearest neighbour bonds
  would change if we flip it (change in the
  systems energy)
• Compute exp(-?E/kT) and compare with a random
  number probability (often use ? to denote
  change)
• If the hit is accepted do the flip
• Repeat for all N spins

10
Computing the Changes
11
Computing Probabilities

• Compute exp( -2 / T ) from previous slide
• Throw a random uniform deviate dice
• Compare the two.
• Decreases in energy will always succeed, but
  small increases in energy can still occur
  depending upon the temperature
• Hence idea of thermal noise or heatbath

12
Practicalities

• Need to initialise the system of spins
• Just make them up/down randomly with 50/50
  probability
• This is effectively a very hot (random) system
• Now need to equilibrate for a particular
  temperature
• Keep running the algorithm until (E or M) stops
  changing much
• Often make the lattice periodic so all sites
  have same number of neighbours

13
More Programming Practicalities

• We need a good quality and fast random number
  generator
• As we will see later the bigger our model size
  the more acccurately we can compute its
  properties, so we need an efficient memory model
• We need to compute expensive functions like
  exp() but we can use a look-up table

14
Shape Function - Scattering

• Simulated systems can be compared with eg neutron
  scattering experiments
• Scattering pattern is radially averaged Fourier
  Transform of simulation

15
What to measure?

• Correlation length can be done using Fourier
  transform
• Other properties such as the number of clusters
  or size of the largest cluster

• Simple measurement is the magnetisation M
• Just the sum of all the spins so ranges from 0
  to N if there are N spins in the system
• Can also count total energy of the system, E,
  the total number of like-like bonds

16
Binder Cumulant

• Sometimes easier to measure B the Binder
  cumulant rather than M to get a more precise fix
  on the transition
• B is computed from M but forms S-shaped
  curves that are easier to plot

• It turns out that B is sensitive to the system
  size, so we can run say 3 different system
  sizes and the point where the the three
  s-shaped curves intercept is where the phase
  transition is.

17
Binder Cumulant

• For a system of size N
• B 1 ( lt M4gt / 3 lt M2gt2 )
• Where lt . gt denotes the average, and 2 etc
  means to the power of 2
• B will form s-shaped curves with T that become
  sharper for bigger systems sizes N
• (sometime in the literature it is called U_N)

18
Correlation length

• Ising system known to behave like a critical
  magnet
• Above Tc, no longer a magnet
• Phase transition
• Near Tc all length scales and esp infinite length
  in infinite system are present

19
Geometries

• Crystal lattice in 2d (square, or hexagonal)
• Lattice in 3d (cube or face-centred cubic or
  body centred cubic)
• But no reason why cannot study arbitrary graph
  or network

• Instead of simple array with 4 or 6 nearest
  neighbours, need to store as a graph structure
  with neighbour lists

20
What if we Mess with the geometry?

• Known that Periodic boundary conditions are
  preferable for large systems
• Large is better as no finite size rounding
  effects

• Various Tc values for different dimensionalities
• eg in 2d Tc 2.269185317 (Onsager, exact)
• In 3d Tc 4.511536
• (Baillie et al, Monte Carlo)

21
What if We Crinkle Space?

• Place shortcuts in the graph?
• What do these mean?
• Wormholes?

• Sure enough when we apply some shortcuts, the
  length scales available are shorter and we see
  rounding effects - shifts in Tc

22
Small-Worlds

• Our wormholes idea was serendipitous - its
  analogous to the small-world effects seen in
  networks (graphs)
• (Watts Strogatz)

• Our wormhole is effectively a small world
  shortcut
• We only looked at one or two - but we can
  systematically look at regimes in between random
  graph and regular meshes

23
Small-World Rewiring Models for 2D Lattices
24
Small-World Parameter p

• Probability that given link is an arbitrary link
  rather than a normal mesh neighbour link
• Or probability that given site connects to
  another given site

• Can either add links (shortcuts) to a mesh
  (small-world regime)
• Or can damage a mesh or lattice - removing links
  with probability p
• Large p - mean field regime
• Medium p - interesting
• Small p - weird!

25
Small-World Ising Systems

• Herrero looked at small 2d and 3d systems in
  large to medium p regime
• Scaling controversy in small p regime - how does
  the shift in Tc scale with p?

• Newman looked at scaling and percolation theory
  and makes some hypotheses - but it is not obvious
  which are correct
• Hence need to probe small-p regime on large
  systems

26
Small-World Shift in Tc
27
Measurements

• Binder cumulant method is useful for computing Tc
• It is sensitive to correlations in the data - a
  lot of data is needed
• There are other things we can look for in the
  data so we are logging E and M and Wolff Number

• At large p the shift is easy to estimate - its is
  very sensitive to noise at low p and hence, we
  need to fight the Monte Carlo convergence effect.

28
Intercepts of Binder Cumulant
29
Technicalities - Speed 1

• N runs - errors for MC goes as 1/ vN
• So best way to parallelise is job-wise scripting
• Unix shell or python or ruby
• Script generates scripts for the jobs, spreads
  them onto Helix queues
• Data logs gathered up and averaged
• Using simple statistical methods - blocks assumed
  to be independent samples

30
Speed 2 - Random Numbers

• Parallel Generator harder to reliably seed - we
  use scatter approach
• Various generator algorithms OK, but combinations
  best
• Ising model is a good test of generators!

• My favourite is currently Marsaglias
  Lagged-Fibonacci generator
• Uses lag-table to mix deviates and achieves long
  period
• No known problems with it

31
Speed Depends on T
32
Speed 3

• Never trust your compiler!
• It is worth worrying about
• Optimisation levels -O6, -DNDEBUG (x3)
• Inline functions (x2)
• Pre-computed neighbours lists (x6)

33
Neighbour Lists 1

• Computing neighbours involves arithmetic -
  especially if using periodic boundary conditions
  (modulo)
• Can speed up if pre-compute these, still faster
  despite cache effects

34
Neighbour Lists

• Simple neighbour list helps speed but hard to
  generalise
• What do you call them in 4d?
• In d-dimensions?

35
Graph Structure

• Graph data structure allows arbitrary neighbour
  lists
• Arbitrary small-world effects - not just mesh
  rewiring?
• Allows us to investigate many different effects
  using the same basic simulation code

36
However - Memory Limited

• Neighbour Lists save speed but make storage x7 in
  3d simulation
• Bit encoding saves space - slower since more
  calculations needed to extract a specific bit
  from a word
• High Performance Computing is still an art!

37
Target Architecture?

• Helix nodes typically 1 or 2 GByte - but less
  than this is available. Also the dual processor
  per node effect means two potential jobs are
  sharing the memory
• Problem - eventually operating System can only
  address 232 4GB of memory
• In practice 32 bit PCs limited to circa 3.5GB

38
Limitations on System Size

• We want to explore regimes of 256x256x256 spin
  sites
• Helix nodes capable of this just, but even if
  simulation were fast enough, 512x512x512 regime
  not accessible to us
• Small-p regime requires large system sizes
• So plan a new architecture - 64 bit nodes with
  16GBytes of memory each
• Need 64-bit random number generator

39
Cluster Update Methods

• Swendson and Wang invented a method of updating
  whole clusters at once
• Clusters are constructed according to Boltzmann
  link probabilities
• Good z but expensive to identify the clusters

• Wolff invented a simpler one cluster method which
  has even better z and is lightly easier to
  compute
• However, z in the different regimes is unknown,
  so we are having to experiment to tune the
  simulation
• Our data structure works with either

40
3D Ising System Results
41
Dimensional Dependence
42
Critical Exponent Beta
43
Dijkstra All-Pairs
44
Some Results

• Wolff Cluster update works on small-world
• At small-p scaling behaviour looks consistent
• Can characterise a complex network by its
  departure behaviour properties from known
  (simple) networks

• In 2D shift in Tc goes as p0.5
• In 3D still a power law px
• Looks like also in 4D, and possibly 5D (still
  working on that)
• Beta (and nu) dependence quite subtle - challenge
  universality class ideas

45
Summary

• arXiv.org/abs/cond-mat/0611763
• www.massey.ac.nz/kahawick/cstn/036/abs-036.html
• More work in progress
• Interesting Scaling behaviour
• Network Models interesting aside from Ising
  effect
• All the High Performance Computing skills are
  needed for this one
• Contact
• Ken Hawick, Massey University
• k.a.hawick_at_massey.ac.nz