162 lines
18 KiB
HTML
162 lines
18 KiB
HTML
<!doctype html>
|
||
<html lang="en">
|
||
<head>
|
||
<meta charset="utf-8">
|
||
<meta http-equiv="x-ua-compatible" content="ie=edge">
|
||
<meta name="viewport" content="width=device-width, initial-scale=1, user-scalable=yes">
|
||
<meta name="description" content="Dimitri Lozeve's blog: Ising model simulation">
|
||
|
||
<title>Dimitri Lozeve - Ising model simulation</title>
|
||
<link rel="stylesheet" href="../css/tufte.css" />
|
||
<link rel="stylesheet" href="../css/pandoc.css" />
|
||
<link rel="stylesheet" href="../css/default.css" />
|
||
<link rel="stylesheet" href="../css/syntax.css" />
|
||
|
||
<!-- RSS feed -->
|
||
<link rel="alternate" type="application/rss+xml" title="Dimitri Lozeve's blog" href="../rss.xml" />
|
||
|
||
<!-- MathJax -->
|
||
<script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
|
||
<script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
|
||
</head>
|
||
<body>
|
||
<article>
|
||
<header>
|
||
<nav>
|
||
<a href="../">Home</a>
|
||
<a href="../archive.html">Posts</a>
|
||
<a href="../projects.html">Projects</a>
|
||
<a href="../cv.html">CV</a>
|
||
<a href="../contact.html">Contact</a>
|
||
</nav>
|
||
|
||
|
||
<h1 class="title">Ising model simulation</h1>
|
||
|
||
|
||
<p class="byline">February 5, 2018 – Dimitri Lozeve</p>
|
||
|
||
|
||
</header>
|
||
|
||
<section>
|
||
<div id="toc"><h2>Table of Contents</h2><ul>
|
||
<li><a href="#mathematical-definition">Mathematical definition</a></li>
|
||
<li><a href="#simulation">Simulation</a></li>
|
||
<li><a href="#implementation">Implementation</a></li>
|
||
<li><a href="#conclusion">Conclusion</a></li>
|
||
</ul></div>
|
||
<p>The <a href="https://en.wikipedia.org/wiki/Ising_model">Ising model</a> is a model used to represent magnetic dipole moments in statistical physics. Physical details are on the Wikipedia page, but what is interesting is that it follows a complex probability distribution on a lattice, where each site can take the value +1 or -1.</p>
|
||
<p><img src="../images/ising.gif" /></p>
|
||
<h2 id="mathematical-definition">Mathematical definition</h2>
|
||
<p>We have a lattice <span class="math inline">\(\Lambda\)</span> consisting of sites <span class="math inline">\(k\)</span>. For each site, there is a moment <span class="math inline">\(\sigma_k \in \{ -1, +1 \}\)</span>. <span class="math inline">\(\sigma =
|
||
(\sigma_k)_{k\in\Lambda}\)</span> is called the <em>configuration</em> of the lattice.</p>
|
||
<p>The total energy of the configuration is given by the <em>Hamiltonian</em> <span class="math display">\[
|
||
H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j,
|
||
\]</span> where <span class="math inline">\(i\sim j\)</span> denotes <em>neighbours</em>, and <span class="math inline">\(J\)</span> is the <em>interaction matrix</em>.</p>
|
||
<p>The <em>configuration probability</em> is given by: <span class="math display">\[
|
||
\pi_\beta(\sigma) = \frac{e^{-\beta H(\sigma)}}{Z_\beta}
|
||
\]</span> where <span class="math inline">\(\beta = (k_B T)^{-1}\)</span> is the inverse temperature, and <span class="math inline">\(Z_\beta\)</span> the normalisation constant.</p>
|
||
<p>For our simulation, we will use a constant interaction term <span class="math inline">\(J > 0\)</span>. If <span class="math inline">\(\sigma_i = \sigma_j\)</span>, the probability will be proportional to <span class="math inline">\(\exp(\beta J)\)</span>, otherwise it would be <span class="math inline">\(\exp(\beta J)\)</span>. Thus, adjacent spins will try to align themselves.</p>
|
||
<h2 id="simulation">Simulation</h2>
|
||
<p>The Ising model is generally simulated using Markov Chain Monte Carlo (MCMC), with the <a href="https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm">Metropolis-Hastings</a> algorithm.</p>
|
||
<p>The algorithm starts from a random configuration and runs as follows:</p>
|
||
<ol>
|
||
<li>Select a site <span class="math inline">\(i\)</span> at random and reverse its spin: <span class="math inline">\(\sigma'_i = -\sigma_i\)</span></li>
|
||
<li>Compute the variation in energy (hamiltonian) <span class="math inline">\(\Delta E = H(\sigma') - H(\sigma)\)</span></li>
|
||
<li>If the energy is lower, accept the new configuration</li>
|
||
<li>Otherwise, draw a uniform random number <span class="math inline">\(u \in ]0,1[\)</span> and accept the new configuration if <span class="math inline">\(u < \min(1, e^{-\beta \Delta E})\)</span>.</li>
|
||
</ol>
|
||
<h2 id="implementation">Implementation</h2>
|
||
<p>The simulation is in Clojure, using the <a href="http://quil.info/">Quil library</a> (a <a href="https://processing.org/">Processing</a> library for Clojure) to display the state of the system.</p>
|
||
<p>This post is “literate Clojure”, and contains <a href="https://github.com/dlozeve/ising-model/blob/master/src/ising_model/core.clj"><code>core.clj</code></a>. The complete project can be found on <a href="https://github.com/dlozeve/ising-model">GitHub</a>.</p>
|
||
<div class="sourceCode" id="cb1"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true"></a>(<span class="kw">ns</span> ising-model.core</span>
|
||
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true"></a> (<span class="at">:require</span> [quil.core <span class="at">:as</span> q]</span>
|
||
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true"></a> [quil.middleware <span class="at">:as</span> m]))</span></code></pre></div>
|
||
<p>The application works with Quil’s <a href="https://github.com/quil/quil/wiki/Functional-mode-(fun-mode)">functional mode</a>, with each function taking a state and returning an updated state at each time step.</p>
|
||
<p>The <code>setup</code> function generates the initial state, with random initial spins. It also sets the frame rate. The matrix is a single vector in row-major mode. The state also holds relevant parameters for the simulation: <span class="math inline">\(\beta\)</span>, <span class="math inline">\(J\)</span>, and the iteration step.</p>
|
||
<div class="sourceCode" id="cb2"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> setup </span>[size]</span>
|
||
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true"></a> <span class="st">"Setup the display parameters and the initial state"</span></span>
|
||
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true"></a> (q/frame-rate <span class="dv">300</span>)</span>
|
||
<span id="cb2-4"><a href="#cb2-4" aria-hidden="true"></a> (q/color-mode <span class="at">:hsb</span>)</span>
|
||
<span id="cb2-5"><a href="#cb2-5" aria-hidden="true"></a> (<span class="kw">let</span> [matrix (<span class="kw">vec</span> (<span class="kw">repeatedly</span> (<span class="kw">*</span> size size) #(<span class="kw">-</span> (<span class="kw">*</span> <span class="dv">2</span> (<span class="kw">rand-int</span> <span class="dv">2</span>)) <span class="dv">1</span>)))]</span>
|
||
<span id="cb2-6"><a href="#cb2-6" aria-hidden="true"></a> {<span class="at">:grid-size</span> size</span>
|
||
<span id="cb2-7"><a href="#cb2-7" aria-hidden="true"></a> <span class="at">:matrix</span> matrix</span>
|
||
<span id="cb2-8"><a href="#cb2-8" aria-hidden="true"></a> <span class="at">:beta</span> <span class="dv">10</span></span>
|
||
<span id="cb2-9"><a href="#cb2-9" aria-hidden="true"></a> <span class="at">:intensity</span> <span class="dv">10</span></span>
|
||
<span id="cb2-10"><a href="#cb2-10" aria-hidden="true"></a> <span class="at">:iteration</span> <span class="dv">0</span>}))</span></code></pre></div>
|
||
<p>Given a site <span class="math inline">\(i\)</span>, we reverse its spin to generate a new configuration state.</p>
|
||
<div class="sourceCode" id="cb3"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> toggle-state </span>[state i]</span>
|
||
<span id="cb3-2"><a href="#cb3-2" aria-hidden="true"></a> <span class="st">"Compute the new state when we toggle a cell's value"</span></span>
|
||
<span id="cb3-3"><a href="#cb3-3" aria-hidden="true"></a> (<span class="kw">let</span> [matrix (<span class="at">:matrix</span> state)]</span>
|
||
<span id="cb3-4"><a href="#cb3-4" aria-hidden="true"></a> (<span class="kw">assoc</span> state <span class="at">:matrix</span> (<span class="kw">assoc</span> matrix i (<span class="kw">*</span> <span class="dv">-1</span> (matrix i))))))</span></code></pre></div>
|
||
<p>In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site <span class="math inline">\(i\)</span>: <span class="math display">\[ \Delta E =
|
||
J\sigma_i \sum_{j\sim i} \sigma_j. \]</span></p>
|
||
<p>The <code>filter some?</code> is required to eliminate sites outside of the boundaries of the lattice.</p>
|
||
<div class="sourceCode" id="cb4"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> get-neighbours </span>[state idx]</span>
|
||
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true"></a> <span class="st">"Return the values of a cell's neighbours"</span></span>
|
||
<span id="cb4-3"><a href="#cb4-3" aria-hidden="true"></a> [(<span class="kw">get</span> (<span class="at">:matrix</span> state) (<span class="kw">-</span> idx (<span class="at">:grid-size</span> state)))</span>
|
||
<span id="cb4-4"><a href="#cb4-4" aria-hidden="true"></a> (<span class="kw">get</span> (<span class="at">:matrix</span> state) (<span class="kw">dec</span> idx))</span>
|
||
<span id="cb4-5"><a href="#cb4-5" aria-hidden="true"></a> (<span class="kw">get</span> (<span class="at">:matrix</span> state) (<span class="kw">inc</span> idx))</span>
|
||
<span id="cb4-6"><a href="#cb4-6" aria-hidden="true"></a> (<span class="kw">get</span> (<span class="at">:matrix</span> state) (<span class="kw">+</span> (<span class="at">:grid-size</span> state) idx))])</span>
|
||
<span id="cb4-7"><a href="#cb4-7" aria-hidden="true"></a></span>
|
||
<span id="cb4-8"><a href="#cb4-8" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> delta-e </span>[state i]</span>
|
||
<span id="cb4-9"><a href="#cb4-9" aria-hidden="true"></a> <span class="st">"Compute the energy difference introduced by a particular cell"</span></span>
|
||
<span id="cb4-10"><a href="#cb4-10" aria-hidden="true"></a> (<span class="kw">*</span> (<span class="at">:intensity</span> state) ((<span class="at">:matrix</span> state) i)</span>
|
||
<span id="cb4-11"><a href="#cb4-11" aria-hidden="true"></a> (<span class="kw">reduce</span> <span class="kw">+</span> (<span class="kw">filter</span> some? (get-neighbours state i)))))</span></code></pre></div>
|
||
<p>We also add a function to compute directly the hamiltonian for the entire configuration state. We can use it later to log its values across iterations.</p>
|
||
<div class="sourceCode" id="cb5"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> hamiltonian </span>[state]</span>
|
||
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true"></a> <span class="st">"Compute the Hamiltonian of a configuration state"</span></span>
|
||
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true"></a> (<span class="kw">-</span> (<span class="kw">reduce</span> <span class="kw">+</span> (<span class="kw">for</span> [i (<span class="kw">range</span> (<span class="kw">count</span> (<span class="at">:matrix</span> state)))</span>
|
||
<span id="cb5-4"><a href="#cb5-4" aria-hidden="true"></a> j (<span class="kw">filter</span> some? (get-neighbours state i))]</span>
|
||
<span id="cb5-5"><a href="#cb5-5" aria-hidden="true"></a> (<span class="kw">*</span> (<span class="at">:intensity</span> state) ((<span class="at">:matrix</span> state) i) j)))))</span></code></pre></div>
|
||
<p>Finally, we put everything together in the <code>update-state</code> function, which will decide whether to accept or reject the new configuration.</p>
|
||
<div class="sourceCode" id="cb6"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> update-state </span>[state]</span>
|
||
<span id="cb6-2"><a href="#cb6-2" aria-hidden="true"></a> <span class="st">"Accept or reject a new state based on energy</span></span>
|
||
<span id="cb6-3"><a href="#cb6-3" aria-hidden="true"></a><span class="st"> difference (Metropolis-Hastings)"</span></span>
|
||
<span id="cb6-4"><a href="#cb6-4" aria-hidden="true"></a> (<span class="kw">let</span> [i (<span class="kw">rand-int</span> (<span class="kw">count</span> (<span class="at">:matrix</span> state)))</span>
|
||
<span id="cb6-5"><a href="#cb6-5" aria-hidden="true"></a> new-state (toggle-state state i)</span>
|
||
<span id="cb6-6"><a href="#cb6-6" aria-hidden="true"></a> alpha (q/exp (<span class="kw">-</span> (<span class="kw">*</span> (<span class="at">:beta</span> state) (delta-e state i))))]</span>
|
||
<span id="cb6-7"><a href="#cb6-7" aria-hidden="true"></a> <span class="co">;;(println (hamiltonian new-state))</span></span>
|
||
<span id="cb6-8"><a href="#cb6-8" aria-hidden="true"></a> (<span class="kw">update</span> (<span class="kw">if</span> (<span class="kw"><</span> (<span class="kw">rand</span>) alpha) new-state state)</span>
|
||
<span id="cb6-9"><a href="#cb6-9" aria-hidden="true"></a> <span class="at">:iteration</span> <span class="kw">inc</span>)))</span></code></pre></div>
|
||
<p>The last thing to do is to draw the new configuration:</p>
|
||
<div class="sourceCode" id="cb7"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true"></a> (<span class="bu">defn</span><span class="fu"> draw-state </span>[state]</span>
|
||
<span id="cb7-2"><a href="#cb7-2" aria-hidden="true"></a> <span class="st">"Draw a configuration state as a grid"</span></span>
|
||
<span id="cb7-3"><a href="#cb7-3" aria-hidden="true"></a> (q/background <span class="dv">255</span>)</span>
|
||
<span id="cb7-4"><a href="#cb7-4" aria-hidden="true"></a> (<span class="kw">let</span> [cell-size (<span class="kw">quot</span> (q/width) (<span class="at">:grid-size</span> state))]</span>
|
||
<span id="cb7-5"><a href="#cb7-5" aria-hidden="true"></a> (<span class="kw">doseq</span> [[i v] (map-indexed <span class="kw">vector</span> (<span class="at">:matrix</span> state))]</span>
|
||
<span id="cb7-6"><a href="#cb7-6" aria-hidden="true"></a>(<span class="kw">let</span> [x (<span class="kw">*</span> cell-size (<span class="kw">rem</span> i (<span class="at">:grid-size</span> state)))</span>
|
||
<span id="cb7-7"><a href="#cb7-7" aria-hidden="true"></a> y (<span class="kw">*</span> cell-size (<span class="kw">quot</span> i (<span class="at">:grid-size</span> state)))]</span>
|
||
<span id="cb7-8"><a href="#cb7-8" aria-hidden="true"></a> (q/no-stroke)</span>
|
||
<span id="cb7-9"><a href="#cb7-9" aria-hidden="true"></a> (q/fill</span>
|
||
<span id="cb7-10"><a href="#cb7-10" aria-hidden="true"></a> (<span class="kw">if</span> (<span class="kw">=</span> <span class="dv">1</span> v) <span class="dv">0</span> <span class="dv">255</span>))</span>
|
||
<span id="cb7-11"><a href="#cb7-11" aria-hidden="true"></a> (q/rect x y cell-size cell-size))))</span>
|
||
<span id="cb7-12"><a href="#cb7-12" aria-hidden="true"></a> <span class="co">;;(when (zero? (mod (:iteration state) 50)) (q/save-frame "img/ising-######.jpg"))</span></span>
|
||
<span id="cb7-13"><a href="#cb7-13" aria-hidden="true"></a> )</span></code></pre></div>
|
||
<p>And to reset the simulation when the user clicks anywhere on the screen:</p>
|
||
<div class="sourceCode" id="cb8"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true"></a>(<span class="bu">defn</span><span class="fu"> mouse-clicked </span>[state event]</span>
|
||
<span id="cb8-2"><a href="#cb8-2" aria-hidden="true"></a> <span class="st">"When the mouse is clicked, reset the configuration to a random one"</span></span>
|
||
<span id="cb8-3"><a href="#cb8-3" aria-hidden="true"></a> (setup <span class="dv">100</span>))</span></code></pre></div>
|
||
<div class="sourceCode" id="cb9"><pre class="sourceCode clojure"><code class="sourceCode clojure"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true"></a>(q/defsketch ising-model</span>
|
||
<span id="cb9-2"><a href="#cb9-2" aria-hidden="true"></a> <span class="at">:title</span> <span class="st">"Ising model"</span></span>
|
||
<span id="cb9-3"><a href="#cb9-3" aria-hidden="true"></a> <span class="at">:size</span> [<span class="dv">300</span> <span class="dv">300</span>]</span>
|
||
<span id="cb9-4"><a href="#cb9-4" aria-hidden="true"></a> <span class="at">:setup</span> #(setup <span class="dv">100</span>)</span>
|
||
<span id="cb9-5"><a href="#cb9-5" aria-hidden="true"></a> <span class="at">:update</span> update-state</span>
|
||
<span id="cb9-6"><a href="#cb9-6" aria-hidden="true"></a> <span class="at">:draw</span> draw-state</span>
|
||
<span id="cb9-7"><a href="#cb9-7" aria-hidden="true"></a> <span class="at">:mouse-clicked</span> mouse-clicked</span>
|
||
<span id="cb9-8"><a href="#cb9-8" aria-hidden="true"></a> <span class="at">:features</span> [<span class="at">:keep-on-top</span> <span class="at">:no-bind-output</span>]</span>
|
||
<span id="cb9-9"><a href="#cb9-9" aria-hidden="true"></a> <span class="at">:middleware</span> [m/fun-mode])</span></code></pre></div>
|
||
<h2 id="conclusion">Conclusion</h2>
|
||
<p>The Ising model is a really easy (and common) example use of MCMC and Metropolis-Hastings. It allows to easily and intuitively understand how the algorithm works, and to make nice visualizations!</p>
|
||
</section>
|
||
|
||
|
||
</article>
|
||
|
||
<footer>
|
||
Site proudly generated by
|
||
<a href="http://jaspervdj.be/hakyll">Hakyll</a>
|
||
</footer>
|
||
</body>
|
||
</html>
|