{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "ebb8ceea", "metadata": {}, "outputs": [], "source": [ "using Distributed" ] }, { "cell_type": "code", "execution_count": null, "id": "3ffaf32e", "metadata": {}, "outputs": [], "source": [ "if nprocs() == 1\n", " addprocs(4)\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "459365e4", "metadata": {}, "outputs": [], "source": [ "@everywhere begin\n", " using Pkg\n", " Pkg.activate(\".\")\n", " Pkg.instantiate()\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "db328410", "metadata": {}, "outputs": [], "source": [ "using PyPlot\n", "using ProgressMeter\n", "@everywhere using Random" ] }, { "cell_type": "code", "execution_count": null, "id": "8e5c9108", "metadata": {}, "outputs": [], "source": [ "@everywhere struct IsingSystem\n", " Lx::Int64 # number of spins in x direction \n", " Ly::Int64 # number of spins in y direction \n", " N::Int64 # total number of spins, N = Lx*Ly\n", " S::Matrix{Int64} # the array of spins\n", " seed::Int64 # rng seed\n", " r::AbstractRNG # random number generator, r = MersenneTwister(seed)\n", " function IsingSystem(Lx, seed)\n", " Ly = Lx\n", " N = Lx*Ly\n", " r = MersenneTwister(seed)\n", " S = rand(r, [-1, 1], Ly, Lx) # hot start\n", " new(Lx, Ly, N, S, seed, r)\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "13d2d2ee", "metadata": {}, "outputs": [], "source": [ "@everywhere mutable struct Simulation\n", " T::Float64 # temperature, in dimensionless units kT/J\n", " nsweeps::Int64 # number of sweeps\n", " therm_sweeps::Int64 # number of thermalization sweeps \n", " function Simulation(T, nsweeps)\n", " therm_sweeps = div(nsweeps, 4) \n", " new(T, nsweeps, therm_sweeps)\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "c274c380", "metadata": {}, "outputs": [], "source": [ "@everywhere function nearest_neighbours(sys)\n", " lattice = reshape(1:sys.N, (sys.Ly, sys.Lx))\n", " ups = circshift(lattice, (-1,0))\n", " rights = circshift(lattice, (0,-1))\n", " downs = circshift(lattice,(1,0))\n", " lefts = circshift(lattice,(0,1))\n", " #\n", " leftups = circshift(lattice,(-1,1))\n", " rightdowns = circshift(lattice,(1,-1))\n", " #\n", " neighs = vcat(ups[:]', rights[:]', downs[:]', \n", " lefts[:]', leftups[:]', rightdowns[:]')\n", " return neighs\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "0700e0ba", "metadata": {}, "outputs": [], "source": [ "@everywhere function boltzmann_factors!(sim)\n", " t = sim.T\n", " w = zeros(13)\n", " for i = -6:2:6\n", " w[i + 7] = exp(-2*i/t)\n", " end\n", " return w\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "8d212435", "metadata": {}, "outputs": [], "source": [ "@everywhere function energy_per_spin(sys, neighs)\n", " \n", " conf = sys.S\n", " energ = 0.0\n", "\n", " # Convenience functions to obtain neighbors\n", " up(i) = neighs[1, i]\n", " right(i) = neighs[2, i]\n", " leftup(i) = neighs[5, i]\n", "\n", " for i in eachindex(conf)\n", " energ += conf[i] * \n", " (conf[up(i)] + conf[right(i)] + conf[leftup(i)])\n", " end\n", " return -energ / sys.N\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "414812f0", "metadata": {}, "outputs": [], "source": [ "@everywhere function metropolis_one_sweep!(sys, neighs, w)\n", "\n", " conf = sys.S\n", " rng = sys.r\n", "\n", " for i in eachindex(conf)\n", " # Sum of neighbor spins\n", " sumNeigh = conf[neighs[1,i]] + conf[neighs[2,i]] + \n", " conf[neighs[3,i]] + conf[neighs[4,i]] + \n", " conf[neighs[5,i]] + conf[neighs[6,i]]\n", " ss = conf[i]*sumNeigh # possible values -6, -4, -2, 0, 2, 4, 6\n", " # Ratio of Boltzmann factors\n", " ratio = w[ss+7]\n", " # Metropolis criterium\n", " if ratio > 1 || rand(rng) < ratio\n", " conf[i] *= -1 # flip the spin\n", " end\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "4a584395", "metadata": {}, "outputs": [], "source": [ "@everywhere function metropolis!(sys, sim, neighs)\n", "\n", " # Set parameters & initialize\n", " conf = sys.S\n", " \n", " # Precompute Boltzmann factors\n", " w = boltzmann_factors!(sim)\n", " \n", " # thermalize configuration\n", " for j in 1:sim.therm_sweeps\n", " # Move over the lattice and propose to flip each spin\n", " metropolis_one_sweep!(sys, neighs, w)\n", " end\n", "\n", " eav = 0.0\n", " e2av = 0.0\n", " for j in 1:sim.nsweeps\n", " # Move over the lattice and try to flip each spin\n", " metropolis_one_sweep!(sys, neighs, w)\n", "\n", " e = energy_per_spin(sys, neighs)\n", "\n", " eav += e;\n", " e2av += e * e;\n", " end\n", " \n", " # Average value of energy per spin\n", " eav /= sim.nsweeps\n", " e2av /= sim.nsweeps\n", "\n", " # Specific heat dE/dt\n", " cv = sys.N*(e2av - eav*eav)/(sim.T)^2\n", "\n", " return cv\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "83829751", "metadata": {}, "outputs": [], "source": [ "Tmin = 3.2\n", "Tmax = 4.0\n", "np = 40\n", "Ts = range(Tmin, Tmax, length=np)" ] }, { "cell_type": "code", "execution_count": null, "id": "c1139db5", "metadata": {}, "outputs": [], "source": [ "@everywhere function monte_carlo2(t)\n", " sys = IsingSystem(32, myid())\n", " sim = Simulation(t, 1024*256)\n", "\n", " # Build nearest neighbor lookup table\n", " neighs = nearest_neighbours(sys)\n", "\n", " return metropolis!(sys, sim, neighs)\n", "end\n", "\n", "Cvs2 = @showprogress pmap(monte_carlo2, Ts);" ] }, { "cell_type": "code", "execution_count": null, "id": "317efc4d", "metadata": {}, "outputs": [], "source": [ "@everywhere function monte_carlo3(t)\n", " sys = IsingSystem(16, (4+1)*myid())\n", " sim = Simulation(t, 1024*128)\n", "\n", " # Build nearest neighbor lookup table\n", " neighs = nearest_neighbours(sys)\n", "\n", " return metropolis!(sys, sim, neighs)\n", "end\n", "\n", "Cvs3 = @showprogress pmap(monte_carlo3, Ts);" ] }, { "cell_type": "code", "execution_count": null, "id": "b70f0fbb", "metadata": {}, "outputs": [], "source": [ "figure(2)\n", "plot(Ts, Cvs2, \".-\", label=L\"$32 \\times 32$\")\n", "plot(Ts, Cvs3, \".-\", label=L\"$16 \\times 16$\")\n", "grid(true)\n", "legend()\n", "xlabel(\"temperature\")\n", "ylabel(\"heat capacity per spin\")" ] }, { "cell_type": "code", "execution_count": null, "id": "1d694cb1", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.6.2", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" } }, "nbformat": 4, "nbformat_minor": 5 }