I normally write my Monte Carlo codes in Fortran for speed, but I was doing some quick and dirty work and wrote one in Mathematica for the XY model on a square lattice (see Kosterlitz-Thouless Transition).

Monte Carlo is normally all about speed, since we are usually interested in very large systems (i.e. the "thermodynamic limit"). I usually use Mathematica for quick and dirty calculations and for its graphical capabilities, and rarely have considered too carefully the performance of my code. I'm curious if anyone can point out any major shortfalls in this code.

L = 20;
T = 1.0;
\[Alpha] = 1.;
pi = N@\[Pi];

(*acceptance rate array*)
accepts = ConstantArray[1, L^2];

(*mod for periodic boundaries*)
mod[n_] := 1 + Mod[n - 1, L]

(*lists of all lattice indexes *)
pos = Flatten[Table[{i, j}, {i, L}, {j, L}], 1];

(*list of neighbors of site (i,j)*)
neighbors = Table[{{i, mod[j + 1]}, {mod[i + 1], j}, {i, mod[j - 1]}, {mod[i - 1], j}}, {i, 1, L}, {j, 1, L}];

(*initialized lattice of 2d spins*)
XY = Map[Normalize, RandomReal[{-1., 1.}, {L, L, 2}], {-2}];

(*metropolis update for site (i,j)*)
MCUpdate[{i_, j_}] := (
  v = XY[[i, j]]; (*current vector*) 
  u = RotationMatrix[\[Alpha] RandomReal[{-pi,pi}]].v; (*new proposed vector*)
  dE = (Total@Extract[XY, neighbors[[i, j]]]).(u - v); (*change in energy by changing to new vector*)
  n = j + (i - 1) L;
  If[RandomReal[] < Min[1, Quiet@Exp[-dE/T]],
   (*update vector if Metropolis condition is satisfied*)
   XY[[i, j]] = u;
   accepts[[n]] = 1
   (*reject update*)
   accepts[[n]] = 0

(*performs n sweeps of the lattice*)
Sweep[n_] := Do[
  MCUpdate /@ pos;
  (*make sure everything is normalized*)
  Map[Normalize, XY, {-2}];
  If[Mean[accepts] < 0.45, \[Alpha] *= .95, 
   If[Mean[accepts] > 0.55, \[Alpha] *= 1.05]];
  , n]

Some basic explanations:

  • accepts tracks the last L^2 attempts to update one of the variables, so that the magnitude of the proposed changes, \[Alpha], can be adjusted such that the acceptance rate stays around 0.5.
  • MCUpdate proposes a rotation to the vector at site (i,j), computes the change in energy dE due to this change, and then applies the Metropolis condition to determine whether to accept or reject with probability min(1,exp(-dE/T)) where T is the temperature

If you'd like to play with this, the following function will find all the vortices and color them:

(* list of indexes for each plaquette*)
plaq = Flatten[Table[{{i, j}, {i, mod[j + 1]}, {mod[i + 1], mod[j + 1]}, {mod[i + 1], j}}, {i, L}, {j, L}], 1]; (*find vortices and color them*)
vortexcolors = {Blue, Transparent, Red};
vortex := Table[
  p = plaq[[i]];
  spins = Extract[XY, plaq[[i]]];
  \[CapitalDelta] = Round@Sum[
     s1 = spins[[j]]~Append~0; s2 = spins[[1 + Mod[j, 4]]]~Append~0; 
     Sign[Cross[s1, s2].{0, 0, 1}] ArcCos[s1.s2]/(2 \[Pi])
     {j, 4}];
  Graphics[{vortexcolors[[\[CapitalDelta] + 2]], Disk[p[[1]] - {.5, .5}, 0.25]}], {i, L^2}]

and this will allow you to plot all the spins, starting from a random configuration and "sweeping" the lattice some times to come to thermal equilibrium:

(*random starting configuration*)
XY = Map[Normalize, RandomReal[{-1., 1.}, {L, L, 2}], {-2}];
(*arrows to represent vectors*)
plotXY[col_] := (
  Map[(p = First@Position[XY, #] - {1, 1}; Graphics[{col, Arrow[{p, p + #/1.5}]}]) &, XY, {-2}]
(*100 monte carlo sweeps*)
(*plot the configuration with vortices and square lattice grid*)
p1 = {plotXY[Black], vortex};
grid = Graphics[{Lighter@Red, Line[#]}] & /@ (Table[{{i, 0}, {i, L}}, {i, 0, L - 1}]~Join~Table[{{0, i}, {L, i}}, {i, 0, L - 1}]);
Show[{grid, p1}, PlotRange -> {{-1, L}, {-1, L}}]

enter image description here

Edit: As a bonus, is there a better way than Quiet@Exp[-dE/T] to deal with the error that Exp[-large number] cannot be represented by machine precision numbers? I would expect it to simply output zero.

  • $\begingroup$ Have you ever read about Compile? $\endgroup$ Aug 24, 2019 at 1:53
  • $\begingroup$ @HenrikSchumacher I have, I've used it once or twice to great effect but I don't have enough experience with it to really understand when it's useful and how to use it effectively. It seems to me it would not work with global variables, is that true? $\endgroup$
    – Kai
    Aug 24, 2019 at 2:23
  • 2
    $\begingroup$ @Kai Nope. I've used it to simulate Ising models that can perform very fast but remain quite flexible, to do many body simulations that can include things like friction, and to do a very general form of MCMC. No global variables for any of that. $\endgroup$
    – b3m2a1
    Aug 24, 2019 at 2:38
  • $\begingroup$ Re. the bonus question, @HenrikSchumacher showed one way to handle this here. He used Exp[Clip[x, {-700., ∞}]]. The approximation should have a negligible effect on the simulation and at the same time, it avoids underflow, according to his answer. $\endgroup$
    – C. E.
    Aug 24, 2019 at 11:11
  • $\begingroup$ Yes. Good point. Here I used Exp[Min[Max[-700., -dE/T], 0.]] in the compiled code since Clip is not compilable. However, Clip is usually nicer to use because it has attribute Listable. $\endgroup$ Aug 24, 2019 at 14:02

2 Answers 2


You can use the RuntimeTools`Profile function, shown here, to figure out which parts need optimization the most. It's better to compile everything, if possible, but it is also possible to compile just the functions that need it the most.

Applying this function to 100 iterations of your code, we get

enter image description here

RotationMatrix really stands out here, so we can replace it by

rotationMatrix = Compile[{th}, {
    {Cos[th], -Sin[th]},
    {Sin[th], Cos[th]}
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"

This simple change brings the total time down from 9.3 seconds to 5.8 seconds.

Predefining those matrices outside the loop using

rotmats = Map[rotationMatrix, α RandomReal[{-pi, pi}, {L, L}], {2}];

and inside the loop u = rotmats[[i, j]].v brings the time down to 5.0.

Proceeding down the list of things that take time, we find that the condition in the If statement takes some time.

accept = Compile[{dE, T},
   RandomReal[] < Min[1, Exp[-dE/T]],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"

brings the time down to 4.3.

Keep in mind that not all functions are compilable, a list of functions that are compilable is available here. When in doubt whether the code is entirely compilable, use CompiledFunctionTools`CompilePrint and check for calls to the main evaluator. Sometimes it's worth it to rewrite some code just to make it compilable.

Others will show you more optimizations, I'm sure, I just wanted to share this way of finding opportunities for optimizations.

  • $\begingroup$ Very cool! Thanks a bunch, I need to spend more time learning to use Compile, in the few instances I have used it in the past it made a huge difference. I've never seen any of these tricks you've shown here so thank you for sharing. $\endgroup$
    – Kai
    Aug 24, 2019 at 4:27

Here is a compiled version of Sweep that employs complex numbers instead of 2-vectors and rotation.

cSweep = Compile[{{XY0, _Complex, 1}, {\[Alpha]0, _Real}, {T, _Real}, {n, _Integer},
    {top, _Integer, 1}, {bottom, _Integer, 1}, {left, _Integer, 1}, {right, _Integer, 1}
   Block[{XY, accepts, rand, \[Phi], u, v, z, dE, m, lo, hi, \[Alpha]},
    m = Length[XY0];
    lo = 0.45 m;
    hi = 0.55 m;
    XY = XY0;
    \[Alpha] = \[Alpha]0;
     accepts = 0.;
     rand = RandomReal[{0., 1.}, m];
     \[Phi] = RandomReal[\[Alpha] {-Pi, Pi}, m];
      v = Compile`GetElement[XY, i];
      u = v Exp[I Compile`GetElement[\[Phi], i]];
      z = Plus[
        Compile`GetElement[XY, Compile`GetElement[top, i]],
        Compile`GetElement[XY, Compile`GetElement[bottom, i]],
        Compile`GetElement[XY, Compile`GetElement[left, i]],
        Compile`GetElement[XY, Compile`GetElement[right, i]]
      dE = Re[Conjugate[z] (u - v)];
      If[Compile`GetElement[rand, i] < Exp[Min[Max[-700., -dE/T], 0.]],
       XY[[i]] = u;
      , {i, 1, m}];
     If[accepts < lo,
      \[Alpha] *= .95;
      If[accepts > hi, \[Alpha] *= 1.05];
     , {n}];
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"

Usage example:

L = 20;
T = 1.0;
\[Alpha] = 1.;
With[{idx = Partition[Range[L^2], L]},
  bottom = Flatten[RotateLeft[idx]];
  top = Flatten[RotateRight[idx]];
  right = Flatten[idx[[All, RotateLeft[Range[L]]]]];
  left = Flatten[idx[[All, RotateRight[Range[L]]]]];
XY0 = Exp[I RandomReal[{-Pi, Pi}, L^2]];

XY = cSweep[XY0, \[Alpha], T, 100000, top, bottom, left, right]; // AbsoluteTiming // First


Please double-check whether the implementation is correct.

Probably it would be more efficient to adjust the algorithm such that all entries of XY are modified at once. But of course, the rejection step has to be adjusted accordingly in order to guarantee the same equilibrium probability distribution for the stochastic process.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.