dimitri/blog
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add MathJax rendering

Browse source
This commit is contained in:
Dimitri Lozeve 2018-11-13 22:48:36 +01:00
parent 15111fd8be
commit 528b427c43
11 changed files with 55 additions and 37 deletions
Download patch file Download diff file Expand all files Collapse all files

30
_site/posts/ising-model.html
View file

@ -6,6 +6,8 @@
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>Dimitri Lozeve - Ising model simulation</title>
<link rel="stylesheet" href="../css/default.css" />
<link rel="stylesheet" href="../css/syntax.css" />
<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_CHTML" async></script>
</head>
<body>
<header>
@ -32,23 +34,23 @@
<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>
<h1 id="mathematical-definition">Mathematical definition</h1>
<p>We have a lattice <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>Λ</mi><annotation encoding="application/x-tex">\Lambda</annotation></semantics></math> consisting of sites <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>k</mi><annotation encoding="application/x-tex">k</annotation></semantics></math>. For each site, there is a moment <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msub><mi>σ</mi><mi>k</mi></msub><mo>∈</mo><mo stretchy="false" form="prefix">{</mo><mo>−</mo><mn>1</mn><mo>,</mo><mo>+</mo><mn>1</mn><mo stretchy="false" form="postfix">}</mo></mrow><annotation encoding="application/x-tex">\sigma_k \in \{ -1, +1 \}</annotation></semantics></math>. <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>σ</mi><mo>=</mo><mo stretchy="false" form="prefix">(</mo><msub><mi>σ</mi><mi>k</mi></msub><msub><mo stretchy="false" form="postfix">)</mo><mrow><mi>k</mi><mo>∈</mo><mi>Λ</mi></mrow></msub></mrow><annotation encoding="application/x-tex">\sigma =
(\sigma_k)_{k\in\Lambda}</annotation></semantics></math> is called the <em>configuration</em> of the lattice.</p>
<p>The total energy of the configuration is given by the <em>Hamiltonian</em> <math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>H</mi><mo stretchy="false" form="prefix">(</mo><mi>σ</mi><mo stretchy="false" form="postfix">)</mo><mo>=</mo><mo>−</mo><munder><mo>∑</mo><mrow><mi>i</mi><mo>∼</mo><mi>j</mi></mrow></munder><msub><mi>J</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mspace width="0.167em"></mspace><msub><mi>σ</mi><mi>i</mi></msub><mspace width="0.167em"></mspace><msub><mi>σ</mi><mi>j</mi></msub><mo>,</mo></mrow><annotation encoding="application/x-tex">
<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,
</annotation></semantics></math> where <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>i</mi><mo>∼</mo><mi>j</mi></mrow><annotation encoding="application/x-tex">i\sim j</annotation></semantics></math> denotes <em>neighbours</em>, and <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>J</mi><annotation encoding="application/x-tex">J</annotation></semantics></math> is the <em>interaction matrix</em>.</p>
<p>The <em>configuration probability</em> is given by: <math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msub><mi>π</mi><mi>β</mi></msub><mo stretchy="false" form="prefix">(</mo><mi>σ</mi><mo stretchy="false" form="postfix">)</mo><mo>=</mo><mfrac><msup><mi>e</mi><mrow><mo>−</mo><mi>β</mi><mi>H</mi><mo stretchy="false" form="prefix">(</mo><mi>σ</mi><mo stretchy="false" form="postfix">)</mo></mrow></msup><msub><mi>Z</mi><mi>β</mi></msub></mfrac></mrow><annotation encoding="application/x-tex">
\]</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}
</annotation></semantics></math> where <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>β</mi><mo>=</mo><mo stretchy="false" form="prefix">(</mo><msub><mi>k</mi><mi>B</mi></msub><mi>T</mi><msup><mo stretchy="false" form="postfix">)</mo><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow><annotation encoding="application/x-tex">\beta = (k_B T)^{-1}</annotation></semantics></math> is the inverse temperature, and <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><msub><mi>Z</mi><mi>β</mi></msub><annotation encoding="application/x-tex">Z_\beta</annotation></semantics></math> the normalisation constant.</p>
<p>For our simulation, we will use a constant interaction term <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>J</mi><mo>&gt;</mo><mn>0</mn></mrow><annotation encoding="application/x-tex">J &gt; 0</annotation></semantics></math>. If <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><msub><mi>σ</mi><mi>i</mi></msub><mo>=</mo><msub><mi>σ</mi><mi>j</mi></msub></mrow><annotation encoding="application/x-tex">\sigma_i = \sigma_j</annotation></semantics></math>, the probability will be proportional to <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mo>exp</mo><mo stretchy="false" form="prefix">(</mo><mi>β</mi><mi>J</mi><mo stretchy="false" form="postfix">)</mo></mrow><annotation encoding="application/x-tex">\exp(\beta J)</annotation></semantics></math>, otherwise it would be <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mo>exp</mo><mo stretchy="false" form="prefix">(</mo><mi>β</mi><mi>J</mi><mo stretchy="false" form="postfix">)</mo></mrow><annotation encoding="application/x-tex">\exp(\beta J)</annotation></semantics></math>. Thus, adjacent spins will try to align themselves.</p>
\]</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 &gt; 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>
<h1 id="simulation">Simulation</h1>
<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 <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math> at random and reverse its spin: <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>σ</mi><msub><mi>′</mi><mi>i</mi></msub><mo>=</mo><mo>−</mo><msub><mi>σ</mi><mi>i</mi></msub></mrow><annotation encoding="application/x-tex">\sigma'_i = -\sigma_i</annotation></semantics></math></li>
<li>Compute the variation in energy (hamiltonian) <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>Δ</mi><mi>E</mi><mo>=</mo><mi>H</mi><mo stretchy="false" form="prefix">(</mo><mi>σ</mi><mi>′</mi><mo stretchy="false" form="postfix">)</mo><mo>−</mo><mi>H</mi><mo stretchy="false" form="prefix">(</mo><mi>σ</mi><mo stretchy="false" form="postfix">)</mo></mrow><annotation encoding="application/x-tex">\Delta E = H(\sigma') - H(\sigma)</annotation></semantics></math></li>
<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 <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>u</mi><mo>∈</mo><mo stretchy="false" form="postfix">]</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo stretchy="false" form="prefix">[</mo></mrow><annotation encoding="application/x-tex">u \in ]0,1[</annotation></semantics></math> and accept the new configuration if <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>u</mi><mo>&lt;</mo><mo>min</mo><mo stretchy="false" form="prefix">(</mo><mn>1</mn><mo>,</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>β</mi><mi>Δ</mi><mi>E</mi></mrow></msup><mo stretchy="false" form="postfix">)</mo></mrow><annotation encoding="application/x-tex">u &lt; \min(1, e^{-\beta \Delta E})</annotation></semantics></math>.</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 &lt; \min(1, e^{-\beta \Delta E})\)</span>.</li>
</ol>
<h1 id="implementation">Implementation</h1>
<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>
@ -57,7 +59,7 @@ H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j,
<a class="sourceLine" id="cb1-2" title="2"> (<span class="at">:require</span> [quil.core <span class="at">:as</span> q]</a>
<a class="sourceLine" id="cb1-3" title="3"> [quil.middleware <span class="at">:as</span> m]))</a></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: <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>β</mi><annotation encoding="application/x-tex">\beta</annotation></semantics></math>, <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>J</mi><annotation encoding="application/x-tex">J</annotation></semantics></math>, and the iteration 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"><a class="sourceLine" id="cb2-1" title="1">(<span class="bu">defn</span><span class="fu"> setup </span>[size]</a>
<a class="sourceLine" id="cb2-2" title="2"> <span class="st">&quot;Setup the display parameters and the initial state&quot;</span></a>
<a class="sourceLine" id="cb2-3" title="3"> (q/frame-rate <span class="dv">300</span>)</a>
@ -68,13 +70,13 @@ H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j,
<a class="sourceLine" id="cb2-8" title="8"> <span class="at">:beta</span> <span class="dv">10</span></a>
<a class="sourceLine" id="cb2-9" title="9"> <span class="at">:intensity</span> <span class="dv">10</span></a>
<a class="sourceLine" id="cb2-10" title="10"> <span class="at">:iteration</span> <span class="dv">0</span>}))</a></code></pre></div>
<p>Given a site <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math>, we reverse its spin to generate a new configuration state.</p>
<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"><a class="sourceLine" id="cb3-1" title="1">(<span class="bu">defn</span><span class="fu"> toggle-state </span>[state i]</a>
<a class="sourceLine" id="cb3-2" title="2"> <span class="st">&quot;Compute the new state when we toggle a cell's value&quot;</span></a>
<a class="sourceLine" id="cb3-3" title="3"> (<span class="kw">let</span> [matrix (<span class="at">:matrix</span> state)]</a>
<a class="sourceLine" id="cb3-4" title="4"> (<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))))))</a></code></pre></div>
<p>In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site <math display="inline" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mi>i</mi><annotation encoding="application/x-tex">i</annotation></semantics></math>: <math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>Δ</mi><mi>E</mi><mo>=</mo><mi>J</mi><msub><mi>σ</mi><mi>i</mi></msub><munder><mo>∑</mo><mrow><mi>j</mi><mo>∼</mo><mi>i</mi></mrow></munder><msub><mi>σ</mi><mi>j</mi></msub><mi>.</mi></mrow><annotation encoding="application/x-tex"> \Delta E =
J\sigma_i \sum_{j\sim i} \sigma_j. </annotation></semantics></math></p>
<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"><a class="sourceLine" id="cb4-1" title="1">(<span class="bu">defn</span><span class="fu"> get-neighbours </span>[state idx]</a>
<a class="sourceLine" id="cb4-2" title="2"> <span class="st">&quot;Return the values of a cell's neighbours&quot;</span></a>