{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation of the Ising model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from numba import jit\n", "from numpy import *\n", "from numpy import random\n", "\n", "N = 20 # size of the Ising system N x N\n", "\n", "@jit(nopython=True)\n", "def CEnergy(latt):\n", " \"Energy of a 2D Ising lattice at particular configuration\"\n", " Ene = 0\n", " for i in range(len(latt)):\n", " for j in range(len(latt)):\n", " S = latt[i,j]\n", " # right below left above\n", " WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]\n", " Ene += -WF*S # Each neighbor gives energy -J==-1\n", " return Ene/2. # Each pair counted twice\n", "\n", "\n", "def RandomL(N):\n", " \"Radom lattice, corresponding to infinite temerature\"\n", " latt = sign(2*random.random((N,N))-1)\n", " latt = array( latt, dtype=int)\n", " return latt\n", "\n", "def Exponents(T):\n", " PW = zeros(9, dtype=float)\n", " # Precomputed exponents : PW[4+x]=exp(-x*2/T)\n", " PW[4+4] = exp(-4.*2/T)\n", " PW[4+2] = exp(-2.*2/T)\n", " PW[4+0] = 1.0\n", " PW[4-2] = exp( 2.*2/T)\n", " PW[4-4] = exp( 4.*2/T)\n", " return PW" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "latt = RandomL(N)\n", "print(latt)\n", "print('Energy=', CEnergy(latt))\n", "T=2.\n", "PW = Exponents(T)\n", "\n", "print('Exponents=', PW)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "@jit(nopython=True)\n", "def SamplePython(Nitt, latt, PW, warm, measure):\n", " \"Monte Carlo sampling for the Ising model in Pythons\"\n", " Ene = CEnergy(latt) # Starting energy\n", " Mn =sum(latt) # Starting magnetization\n", " N = len(latt)\n", " \n", " N1=0 # Measurements\n", " E1,M1,E2,M2=0.0,0.0,0.0,0.0\n", " \n", " N2 = N*N\n", " for itt in range(Nitt):\n", " t = int(random.rand()*N2)\n", " (i,j) = (t % N, int(t/N))\n", " S = latt[i,j]\n", " WF = latt[(i+1)%N, j] + latt[i,(j+1)%N] + latt[(i-1)%N,j] + latt[i,(j-1)%N]\n", " # new configuration -S, old configuration S => magnetization change -2*S\n", " # energy change = (-J)*WF*(-S) - (-J)*WF*(S) = 2*J**WF*S\n", " # We will prepare : PW[4+x]=exp(-x*2/T)\n", " # P = exp(-2*WF*S/T) = exp(-(WF*S)*2/T) == PW[4+WF*S],\n", " # because PW[4+x]=exp(-x*2/T)\n", " P = PW[4+S*WF]\n", " if P>random.rand(): # flip the spin\n", " latt[i,j] = -S\n", " Ene += 2*S*WF\n", " Mn -= 2*S\n", " \n", " if itt>warm and itt%measure==0:\n", " N1 += 1\n", " E1 += Ene\n", " M1 += Mn\n", " E2 += Ene*Ene\n", " M2 += Mn*Mn\n", " \n", " E,M = E1/N1, M1/N1\n", " cv = (E2/N1-E**2)/T**2 # cv =(-^2)/T^2\n", " chi = (M2/N1-M**2)/T # chi=(-^2)/T\n", " return (M/N2, E/N2, cv/N2, chi/N2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Nitt = 5000000 # total number of Monte Carlo steps\n", "warm = 1000 # Number of warmup steps\n", "measure=5 # How often to take a measurement\n", "\n", "\n", "(M, E, cv, chi) = SamplePython(Nitt, latt, PW, warm, measure)\n", "#(aM, aE, cv, chi) = SampleCPP(Nitt, latt, PW,T, warm, measure)\n", "print('=', M/(N*N), '=', E/(N*N) , 'cv=', cv, 'chi=',chi)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wT = linspace(5,0.5,30)\n", "wMag=[]\n", "wEne=[]\n", "wCv=[]\n", "wChi=[]\n", "N2 = N*N\n", "for T in wT:\n", " # Precomputed exponents : PW[4+x]=exp(-x*2/T)\n", " PW = Exponents(T)\n", " \n", " (M, E, cv, chi) = SamplePython(Nitt, latt, PW, warm, measure)\n", " wMag.append( M )\n", " wEne.append( E )\n", " wCv.append( cv )\n", " wChi.append( chi )\n", " print('T=',T, 'M=', M, 'E=', E, 'cv=', cv, 'chi=', chi)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from pylab import *\n", "%matplotlib inline\n", "\n", "plot(wT, wEne, label='E(T)')\n", "plot(wT, wCv, label='cv(T)')\n", "plot(wT, wMag, label='M(T)')\n", "xlabel('T')\n", "legend(loc='best')\n", "show()\n", "plot(wT, wChi, label='chi(T)')\n", "xlabel('T')\n", "legend(loc='best')\n", "show() " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }