{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Computational Physics Lectures: Statistical physics and the Ising Model\n", "\n", " \n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n", "\n", "Date: **Oct 28, 2020**\n", "\n", "Copyright 1999-2020, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n", "\n", "\n", "\n", "\n", "## Ensembles\n", "\n", "In statistical physics the concept of an ensemble is one of the\n", "cornerstones in the definition of thermodynamical quantities. An\n", "ensemble is a collection of microphysics systems from which we derive\n", "expectations values and thermodynamical properties related to\n", "experiment. As an example, the specific heat (which is a measurable\n", "quantity in the laboratory) of a system of infinitely many particles,\n", "can be derived from the basic interactions between the microscopic\n", "constituents. The latter can span from electrons to atoms and\n", "molecules or a system of classical spins. All these microscopic\n", "constituents interact via a well-defined interaction. We say\n", "therefore that statistical physics bridges the gap between the\n", "microscopic world and the macroscopic world. Thermodynamical\n", "quantities such as the specific heat or net magnetization of a system\n", "can all be derived from a microscopic theory.\n", "\n", "\n", "## Famous Ensembles\n", "\n", "The table lists the most used ensembles in statistical physics\n", "together with frequently arising extensive (depend on the size of the\n", "systems such as the number of particles) and intensive variables\n", "(apply to all components of a system), in addition to associated\n", "potentials.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Microcanonical Canonical Grand Canonical Pressure canonical
Exchange of heat no yes yes yes
with the environment
Exchange of particles no no yes no
with the environemt
Thermodynamical $V, \\cal M, \\cal D$ $V, \\cal M, \\cal D$ $V, \\cal M, \\cal D$ $P, \\cal H, \\cal E$
parameters $E$ $T$ $T$ $T$
$N$ $N$ $\\mu$ $N$
Potential Entropy Helmholtz $PV$ Gibbs
$N$ $N$ $\\mu$ $N$
Energy Internal Internal Internal Enthalpy
$N$ $N$ $\\mu$ $N$
\n", "\n", "\n", "## Canonical Ensemble\n", "\n", "One of the most used ensembles is the canonical one, which is related to the microcanonical ensemble\n", "via a Legendre transformation. The temperature is an intensive variable in this ensemble whereas the energy\n", "follows as an expectation value. \n", "In order to calculate expectation values such as the mean energy $\\langle E \\rangle $\n", "at a given temperature, we need a probability distribution.\n", "It is given by the Boltzmann distribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "P_i(\\beta) = \\frac{e^{-\\beta E_i}}{Z}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\beta=1/k_BT$ being the inverse temperature, $k_B$ is the \n", "Boltzmann constant, $E_i$ is the energy of a microstate $i$ while \n", "$Z$ is the partition function for the canonical ensemble\n", "defined as\n", "\n", "\n", "\n", "## The partition function is a normalization constant\n", "\n", "In the canonical ensemble the partition function is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Z=\\sum_{i=1}^{M}e^{-\\beta E_i},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the sum extends over all microstates $M$. \n", "\n", "\n", "## Helmoltz free energy, what does it mean?\n", "\n", "The potential of interest in this case is Helmholtz' free energy. It\n", "relates the expectation value of the energy at a given temperatur $T$\n", "to the entropy at the same temperature via" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "F=-k_{B}TlnZ=\\langle E \\rangle-TS.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Helmholtz' free energy expresses the\n", "struggle between two important principles in physics, namely the\n", "strive towards an energy minimum and the drive towards higher entropy\n", "as the temperature increases. A higher entropy may be interpreted as a\n", "larger degree of disorder. When equilibrium is reached at a given\n", "temperature, we have a balance between these two principles. The\n", "numerical expression is Helmholtz' free energy.\n", "\n", "## Thermodynamical quantities\n", "\n", "In the canonical ensemble the entropy is given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "S =k_{B}lnZ\n", "+k_{B}T\\left(\\frac{\\partial lnZ}{\\partial T}\\right)_{N, V},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the pressure by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p=k_{B}T\\left(\\frac{\\partial lnZ}{\\partial V}\\right)_{N, T}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly we can compute the chemical potential as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\mu =-k_{B}T\\left(\\frac{\\partial lnZ}{\\partial N}\\right)_{V, T}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thermodynamical quantities, the energy in the canonical ensemble\n", "\n", "For a system described by the canonical ensemble, the energy is an\n", "expectation value since we allow energy to be exchanged with the surroundings\n", "(a heat bath with temperature $T$). \n", "\n", "This expectation value, the mean energy,\n", "can be calculated using" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\langle E\\rangle =k_{B}T^{2}\\left(\\frac{\\partial lnZ}{\\partial T}\\right)_{V, N}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or \n", "using the probability distribution\n", "$P_i$ as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\langle E \\rangle = \\sum_{i=1}^M E_i P_i(\\beta)= \n", " \\frac{1}{Z}\\sum_{i=1}^M E_ie^{-\\beta E_i}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Energy and specific heat in the canonical ensemble\n", "\n", "\n", "The energy is proportional to the first derivative of the potential,\n", "Helmholtz' free energy. The corresponding variance is defined as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sigma_E^2=\\langle E^2 \\rangle-\\langle E \\rangle^2=\n", " \\frac{1}{Z}\\sum_{i=1}^M E_i^2e^{-\\beta E_i}-\n", " \\left(\\frac{1}{Z}\\sum_{i=1}^M E_ie^{-\\beta E_i}\\right)^2.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we divide the latter quantity with\n", "$kT^2$ we obtain the specific heat at constant volume" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "C_V= \\frac{1}{k_BT^2}\\left(\\langle E^2 \\rangle-\\langle E \\rangle^2\\right),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which again can be related to the second derivative of Helmholtz' free energy.\n", "\n", "\n", "## Magnetic moments and susceptibility in the canonical ensemble\n", "\n", "Using the same prescription, we can also evaluate the mean magnetization\n", "through" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\langle {\\cal M} \\rangle = \\sum_i^M {\\cal M}_i P_i(\\beta)= \n", " \\frac{1}{Z}\\sum_i^M {\\cal M}_ie^{-\\beta E_i},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the corresponding variance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sigma_{{\\cal M}}^2=\\langle {\\cal M}^2 \\rangle-\\langle {\\cal M} \\rangle^2=\n", " \\frac{1}{Z}\\sum_{i=1}^M {\\cal M}_i^2e^{-\\beta E_i}-\n", " \\left(\\frac{1}{Z}\\sum_{i=1}^M {\\cal M}_ie^{-\\beta E_i}\\right)^2.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This quantity defines also the susceptibility \n", "$\\chi$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\chi=\\frac{1}{k_BT}\\left(\\langle {\\cal M}^2 \\rangle-\\langle {\\cal M} \\rangle^2\\right).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Our model, the Ising model in one and two dimensions\n", "\n", "The model we will employ in our studies of phase transitions at finite temperature for \n", "magnetic systems is the so-called Ising model. In its simplest form\n", "the energy is expressed as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=-J\\sum_{}^{N}s_ks_l-{\\cal B}\\sum_k^Ns_k,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $s_k=\\pm 1$, $N$ is the total number of spins, \n", "$J$ is a coupling constant expressing the strength of the interaction\n", "between neighboring spins and \n", "${\\cal B}$ is an external magnetic field interacting with the magnetic\n", "moment set up by the spins.\n", "\n", "The symbol $$ indicates that we sum over nearest\n", "neighbors only. \n", "Notice that for $J>0$ it is energetically favorable for neighboring spins \n", "to be aligned. This feature leads to, at low enough temperatures,\n", "a cooperative phenomenon called spontaneous magnetization. That is, \n", "through interactions between nearest neighbors, a given magnetic\n", "moment can influence the alignment of spins that are separated \n", "from the given spin by a macroscopic distance. These long range correlations\n", "between spins are associated with a long-range order in which\n", "the lattice has a net magnetization in the absence of a magnetic field. \n", "\n", "## Boltzmann distribution\n", "\n", "In order to calculate expectation values such as the mean energy\n", "$\\langle E \\rangle $ or\n", "magnetization $\\langle {\\cal M} \\rangle $\n", "in statistical physics\n", "at a given temperature, we need a probability distribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "P_i(\\beta) = \\frac{e^{-\\beta E_i}}{Z}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\beta=1/kT$ being the inverse temperature, $k$ the \n", "Boltzmann constant, $E_i$ is the energy of a state $i$ while \n", "$Z$ is the partition function for the canonical ensemble\n", "defined as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Z=\\sum_{i=1}^{M}e^{-\\beta E_i},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the sum extends over all microstates\n", "$M$. \n", "$P_i$ expresses the probability of finding the system in a given \n", "configuration $i$.\n", "\n", "\n", "\n", "## Energy for a specific configuration\n", "\n", "The energy for a specific configuration $i$\n", "is given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_i =-J\\sum_{}^{N}s_ks_l.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configurations\n", "To better understand what is meant with a configuration, \n", "consider first the case of the one-dimensional Ising model\n", "with ${\\cal B}=0$. \n", "In general, a given configuration of $N$ spins in one\n", "dimension may look like" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{cccccccccc}\n", "\\uparrow&\\uparrow&\\uparrow&\\dots&\\uparrow&\\downarrow&\\uparrow&\\dots&\\uparrow&\\downarrow\\\\\n", "1&2&3&\\dots& i-1&i&i+1&\\dots&N-1&N\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to illustrate these features \n", "let us further specialize to\n", "just two spins.\n", "\n", "With two spins, since each spin takes two values only,\n", "we have $2^2=4$ possible arrangements of the two spins. \n", "These four possibilities are" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "1= \\uparrow\\uparrow\\hspace{1cm}\n", " 2= \\uparrow\\downarrow\\hspace{1cm}\n", " 3= \\downarrow\\uparrow\\hspace{1cm}\n", " 4=\\downarrow\\downarrow\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Boundary conditions, free ends\n", "\n", "What is the energy of each of these configurations? \n", "\n", "For small systems, the way we treat the ends matters. Two cases are\n", "often used.\n", "\n", "In the first case we employ what is called free ends. This means that there is no contribution from points to the right or left of the endpoints. For the one-dimensional case, the energy is then written as a sum over a single index" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_i =-J\\sum_{j=1}^{N-1}s_js_{j+1},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Free ends and the energy\n", "\n", "If we label the first spin as $s_1$ and the second as $s_2$ \n", "we obtain the following \n", "expression for the energy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=-Js_1s_2.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The calculation of the energy for the one-dimensional lattice\n", "with free ends for one specific spin-configuration \n", "can easily be implemented in the following lines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " for ( j=1; j < N; j++) {\n", " energy += spin[j]*spin[j+1]; \n", " }\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the vector $spin[]$ contains the spin value $s_k=\\pm 1$. \n", "\n", "## Free ends and energy\n", "\n", "For the specific state $E_1$, we have chosen all spins up. The energy of\n", "this configuration becomes then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_1=E_{\\uparrow\\uparrow}=-J.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The other configurations give" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_2=E_{\\uparrow\\downarrow}=+J,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_3=E_{\\downarrow\\uparrow}=+J,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_4=E_{\\downarrow\\downarrow}=-J.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Periodic boundary conditions\n", "\n", "We can also choose so-called periodic boundary conditions. This means\n", "that the neighbour to the right of $s_N$ is assumed to take the value\n", "of $s_1$. Similarly, the neighbour to the left of $s_1$ takes the\n", "value $s_N$. In this case the energy for the one-dimensional lattice\n", "reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_i =-J\\sum_{j=1}^{N}s_js_{j+1},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and we obtain the following expression for the\n", "two-spin case" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=-J(s_1s_2+s_2s_1).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Energy with PBC\n", "\n", "In this case the energy for $E_1$ is different, we obtain namely" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_1=E_{\\uparrow\\uparrow}=-2J.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The other cases do also differ and we have" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_2=E_{\\uparrow\\downarrow}=+2J,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_3=E_{\\downarrow\\uparrow}=+2J,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E_4=E_{\\downarrow\\downarrow}=-2J.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simple code for PBC\n", "\n", "\n", "If we choose to use periodic boundary conditions we can code the above\n", "expression as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " jm=N;\n", " for ( j=1; j <=N ; j++) {\n", " energy += spin[j]*spin[jm]; \n", " jm = j ;\n", " }\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The magnetization is however the same, defined as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "{\\cal M}_i=\\sum_{j=1}^N s_j,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where we sum over all spins for a given configuration $i$. \n", "\n", "\n", "## Summing up\n", "\n", "The table lists the energy and magnetization for both free ends\n", "and periodic boundary conditions. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
State Energy (FE) Energy (PBC) Magnetization
$1= \\uparrow\\uparrow$ $-J$ $-2J$ 2
$2=\\uparrow\\downarrow$ $J$ $2J$ 0
$ 3=\\downarrow\\uparrow$ $J$ $2J$ 0
$ 4=\\downarrow\\downarrow$ $-J$ $-2J$ -2
\n", "\n", "\n", "## Reorganizing\n", "\n", "We can reorganize according to the number of spins pointing up, as shown in the table here\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Number spins up Degeneracy Energy (FE) Energy (PBC) Magnetization
2 1 $-J$ $-2J$ 2
1 2 $J$ $2J$ 0
0 1 $-J$ $-2J$ -2
\n", "\n", "\n", "## Our model, the Ising model in one and two dimensions\n", "\n", "It is worth noting that for small dimensions of the lattice,\n", "the energy differs depending on whether we use\n", "periodic boundary conditions or free ends. This means also\n", "that the partition functions will be different, as discussed\n", "below. In the thermodynamic limit we have $N\\rightarrow \\infty$,\n", "and the final results do not depend on the kind of boundary conditions\n", "we choose. \n", "\n", "\n", "For a one-dimensional lattice with periodic boundary conditions, \n", "each spin sees two neighbors. For a\n", "two-dimensional lattice each spin sees four neighboring spins. \n", "How many neighbors does a spin see in three dimensions?\n", "\n", "## Ising model in one and two dimensions\n", "\n", "\n", "In a similar way, we could enumerate the number of states for\n", "a two-dimensional system consisting of two spins, i.e., \n", "a $2\\times 2$ Ising model on a square lattice with {\\em periodic\n", "boundary conditions}. In this case we have a total of \n", "$2^4=16$ states. \n", "Some \n", "examples of configurations with their respective energies are \n", "listed here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=-8J\\hspace{1cm}\\begin{array}{cc}\\uparrow & \\uparrow \\\\\n", " \\uparrow & \\uparrow\\end{array}\n", "\\hspace{0.5cm}\n", " E=0\\hspace{1cm}\\begin{array}{cc}\\uparrow & \\uparrow \\\\\n", " \\uparrow & \\downarrow\\end{array}\n", "\\hspace{0.5cm}\n", " E=0\\hspace{1cm}\\begin{array}{cc}\\downarrow & \\downarrow \\\\\n", " \\uparrow & \\downarrow\\end{array}\n", "\\hspace{0.5cm}\n", " E=-8J\\hspace{1cm}\\begin{array}{cc}\\downarrow & \\downarrow \\\\\n", " \\downarrow & \\downarrow\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## List of configurations with energies and magnetic moment\n", "\n", "In the table here we group these configurations\n", "according to their total energy and magnetization.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
Number spins up Degeneracy Energy Magnetization
4 1 $-8J$ 4
3 4 $0$ 2
2 4 $0$ 0
2 2 $8J$ 0
1 4 $0$ -2
0 1 $-8J$ -4
\n", "\n", "## Phase Transitions and Critical Phenomena\n", "\n", "A phase transition is marked by abrupt macroscopic changes as external\n", "parameters are changed, such as an increase of temperature. The point\n", "where a phase transition takes place is called a critical point.\n", "\n", "We distinguish normally between two types of phase transitions;\n", "first-order transitions and second-order transitions. An important\n", "quantity in studies of phase transitions is the so-called correlation\n", "length $\\xi$ and various correlations functions like spin-spin\n", "correlations. For the Ising model we shall show below that the\n", "correlation length is related to the spin-correlation function, which\n", "again defines the magnetic susceptibility. The spin-correlation\n", "function is nothing but the covariance and expresses the degree of\n", "correlation between spins.\n", "\n", "\n", "## Phase Transitions and Critical Phenomena, correlation length\n", "\n", "The correlation length defines the length scale at which the overall\n", "properties of a material start to differ from its bulk properties. It\n", "is the distance over which the fluctuations of the microscopic degrees\n", "of freedom (for example the position of atoms) are significantly\n", "correlated with each other. Usually it is of the order of few\n", "interatomic spacings for a solid. The correlation length $\\xi$\n", "depends however on external conditions such as pressure and\n", "temperature.\n", "\n", "\n", "\n", "\n", "## Classification of phase transitions\n", "\n", "First order/discontinuous phase transitions are characterized by two or more\n", "states on either side of the critical point that can coexist at the\n", "critical point. As we pass through the critical point we observe a\n", "discontinuous behavior of thermodynamical functions. The correlation\n", "length is normally finite at the critical point. Phenomena such as\n", "hysteris occur, viz. there is a continuation of state below the\n", "critical point into one above the critical point. This continuation is\n", "metastable so that the system may take a macroscopically long time to\n", "readjust. A classical example is the melting of ice. It takes a\n", "specific amount of time before all the ice has melted. The temperature\n", "remains constant and water and ice can coexist for a macroscopic\n", "time. The energy shows a discontinuity at the critical point,\n", "reflecting the fact that a certain amount of heat is needed in order\n", "to melt all the ice\n", "\n", "\n", "## Second-order phase Transitions\n", "\n", "Second order or continuous transitions are different and in general\n", "much difficult to understand and model. The correlation length\n", "diverges at the critical point, fluctuations are correlated over all\n", "distance scales, which forces the system to be in a unique critical\n", "phase. The two phases on either side of the critical point become\n", "identical. The disappearance of a spontaneous magnetization is a\n", "classical example of a second-order phase transitions. Structural\n", "transitions in solids are other types of second-order phase\n", "transitions. \n", "\n", "\n", "## Phase Transitions and Critical Phenomena\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
System Transition Order Parameter
Liquid-gas Condensation/evaporation Density difference $\\Delta\\rho=\\rho_{liquid}-\\rho_{gas}$
Binary liquid mixture/Unmixing Composition difference
Quantum liquid Normal fluid/superfluid $<\\phi>$, $\\psi$ = wavefunction
Liquid-solid Melting/crystallisation Reciprocal lattice vector
Magnetic solid Ferromagnetic Spontaneous magnetisation $M$
Antiferromagnetic Sublattice magnetisation $M$
Dielectric solid Ferroelectric Polarization $P$
Antiferroelectric Sublattice polarisation $P$
\n", "\n", "\n", "## Eherenfest definition of phase Transitions\n", "\n", "Using Ehrenfest's definition of the order of a phase transition we can\n", "relate the behavior around the critical point to various derivatives\n", "of the thermodynamical potential. In the canonical ensemble we are\n", "using, the thermodynamical potential is Helmholtz' free energy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "F= \\langle E\\rangle -TS = -kTln Z\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "meaning $ lnZ = -F/kT = -F\\beta$. The energy is given as the first derivative of $F$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\langle E \\rangle=-\\frac{\\partial lnZ}{\\partial \\beta} =\\frac{\\partial (\\beta F)}{\\partial \\beta}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the specific heat is defined via the second derivative of $F$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "C_V=-\\frac{1}{kT^2}\\frac{\\partial^2 (\\beta F)}{\\partial\\beta^2}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Phase Transitions and Critical Phenomena\n", "\n", "We can relate observables to various derivatives of the partition\n", "function and the free energy. When a given derivative of the free\n", "energy or the partition function is discontinuous or diverges\n", "(logarithmic divergence for the heat capacity from the Ising model) we\n", "talk of a phase transition of order of the derivative. A first-order\n", "phase transition is recognized in a discontinuity of the energy, or\n", "the first derivative of $F$. The Ising model exhibits a second-order\n", "phase transition since the heat capacity diverges. The susceptibility\n", "is given by the second derivative of $F$ with respect to external\n", "magnetic field. Both these quantities diverge.\n", "\n", "## The Ising Model and Phase Transitions\n", "\n", "The Ising model in two dimensions with ${\\cal B} = 0$ undergoes a\n", "phase transition of second order. What it actually means is that below\n", "a given critical temperature $T_C$, the Ising model exhibits a\n", "spontaneous magnetization with $\\langle {\\cal M} \\rangle\\ne 0$. Above\n", "$T_C$ the average magnetization is zero. The mean magnetization\n", "approaches zero at $T_C$ with an infinite slope. Such a behavior is\n", "an example of what are called critical phenomena. A critical\n", "phenomenon is normally marked by one or more thermodynamical variables\n", "which vanish above a critical point. In our case this is the mean\n", "magnetization $\\langle {\\cal M} \\rangle\\ne 0$. Such a parameter is\n", "normally called the order parameter.\n", "\n", "\n", "## The Ising Model and Phase Transitions, mean magnetization\n", "\n", "It is possible to show that the mean magnetization is given by\n", "(for temperature below $T_C$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\langle {\\cal M}(T) \\rangle \\sim \\left(T-T_C\\right)^{\\beta},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\beta=1/8$ is a so-called critical exponent. A similar relation\n", "applies to the heat capacity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "C_V(T) \\sim \\left|T_C-T\\right|^{-\\alpha},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the susceptibility" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\chi(T) \\sim \\left|T_C-T\\right|^{-\\gamma},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\alpha = 0$ and $\\gamma = -7/4$.\n", "\n", "\n", "## The Ising Model and Phase Transitions, correlation length\n", "\n", "Another important quantity is the correlation length, which is expected\n", "to be of the order of the lattice spacing for $T$ is close to $T_C$. Because the spins\n", "become more and more correlated as $T$ approaches $T_C$, the correlation\n", "length increases as we get closer to the critical temperature. The discontinuous \n", "behavior of the correlation $\\xi$ near $T_C$ is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\xi(T) \\sim \\left|T_C-T\\right|^{-\\nu}.\n", "\\label{eq:xi} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Ising Model and Phase Transitions, critical behavior\n", "\n", "A second-order phase transition is characterized by a correlation\n", "length which spans the whole system. The correlation length is\n", "typically of the order of some few interatomic distances. The fact\n", "that a system like the Ising model, whose energy is described by the\n", "interaction between neighboring spins only, can yield correlation\n", "lengths of macroscopic size at a critical point is still a feature\n", "which is not properly understood. \n", "\n", "## The Ising Model and Phase Transitions, critical temperature\n", "\n", "In our actual calculations of the two-dimensional Ising model, we are however \n", "always limited to a finite lattice and $\\xi$ will\n", "be proportional with the size of the lattice at the critical point. \n", "Through finite size scaling relations\n", "it is possible to relate the behavior at finite lattices with the \n", "results for an infinitely large lattice.\n", "The critical temperature scales then as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " T_C(L)-T_C(L=\\infty) \\propto aL^{-1/\\nu},\n", "\\label{eq:tc} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $a$ a constant and $\\nu$ defined in Eq. ([1](#eq:xi)).\n", "\n", "## The Ising Model and Phase Transitions, correlation length\n", "\n", "The correlation length for a finite lattice size can then be shown to be proportional to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\xi(T) \\propto L\\sim \\left|T_C-T\\right|^{-\\nu}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and if we set $T=T_C$ one can obtain the following relations for the \n", "magnetization, energy and susceptibility for $T \\le T_C$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\langle {\\cal M}(T) \\rangle \\sim \\left(T-T_C\\right)^{\\beta}\n", " \\propto L^{-\\beta/\\nu},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "C_V(T) \\sim \\left|T_C-T\\right|^{-\\gamma} \\propto L^{\\alpha/\\nu},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\chi(T) \\sim \\left|T_C-T\\right|^{-\\alpha} \\propto L^{\\gamma/\\nu}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Metropolis Algorithm and the Two-dimensional Ising Model\n", "\n", "In our case we have as the Monte Carlo sampling function the probability\n", "for finding the system in a state $s$ given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "P_s=\\frac{e^{-(\\beta E_s)}}{Z},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with energy $E_s$, $\\beta=1/kT$ and $Z$ is a normalization constant which\n", "defines the partition function in the canonical ensemble. As discussed\n", "above" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Z(\\beta)=\\sum_se^{-(\\beta E_s)}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "is difficult to compute since we need all states. \n", "\n", "## The Metropolis Algorithm and the Two-dimensional Ising Model\n", "\n", "In a calculation of the Ising model in two dimensions, the number of\n", "configurations is given by $2^N$ with $N=L\\times L$ the number of\n", "spins for a lattice of length $L$. Fortunately, the Metropolis\n", "algorithm considers only ratios between probabilities and we do not\n", "need to compute the partition function at all. The algorithm goes as\n", "follows\n", "\n", "1. Establish an initial state with energy $E_b$ by positioning yourself at a random configuration in the lattice\n", "\n", "2. Change the initial configuration by flipping e.g., one spin only. Compute the energy of this trial state $E_t$.\n", "\n", "3. Calculate $\\Delta E=E_t-E_b$. The number of values $\\Delta E$ is limited to five for the Ising model in two dimensions, see the discussion below.\n", "\n", "4. If $\\Delta E \\le 0$ we accept the new configuration, meaning that the energy is lowered and we are hopefully moving towards the energy minimum at a given temperature. Go to step 7.\n", "\n", "5. If $\\Delta E > 0$, calculate $w=e^{-(\\beta \\Delta E)}$.\n", "\n", "6. Compare $w$ with a random number $r$. If $r \\le w$, then accept the new configuration, else we keep the old configuration.\n", "\n", "7. The next step is to update various expectations values.\n", "\n", "8. The steps (2)-(7) are then repeated in order to obtain a sufficently good representation of states.\n", "\n", "9. Each time you sweep through the lattice, i.e., when you have summed over all spins, constitutes what is called a Monte Carlo cycle. You could think of one such cycle as a measurement. At the end, you should divide the various expectation values with the total number of cycles. You can choose whether you wish to divide by the number of spins or not. If you divide with the number of spins as well, your result for e.g., the energy is now the energy per spin.\n", "\n", "## The Metropolis Algorithm and the Two-dimensional Ising Model, practical issues\n", "\n", "The crucial step is the calculation of the energy difference and the\n", "change in magnetization. This part needs to be coded in an as\n", "efficient as possible way since the change in energy is computed many\n", "times. In the calculation of the energy difference from one spin\n", "configuration to the other, we will limit the change to the flipping\n", "of one spin only. For the Ising model in two dimensions it means that\n", "there will only be a limited set of values for $\\Delta E$. Actually,\n", "there are only five possible values. \n", "\n", "## Five possible energy differences\n", "\n", "To see this, select first a\n", "random spin position $x,y$ and assume that this spin and its nearest\n", "neighbors are all pointing up. The energy for this configuration is\n", "$E=-4J$. Now we flip this spin as shown below. The energy of the new\n", "configuration is $E=4J$, yielding $\\Delta E=8J$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=-4J\\hspace{1cm}\\begin{array}{ccc} & \\uparrow & \\\\\n", " \\uparrow & \\uparrow & \\uparrow\\\\\n", " & \\uparrow & \\end{array}\n", "\\hspace{1cm}\\Longrightarrow\\hspace{1cm}\n", " E=4J\\hspace{1cm}\\begin{array}{ccc} & \\uparrow & \\\\\n", " \\uparrow & \\downarrow & \\uparrow\\\\\n", " & \\uparrow & \\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The four other possibilities are as follows" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=-2J\\hspace{1cm}\\begin{array}{ccc} & \\uparrow & \\\\\n", " \\downarrow & \\uparrow & \\uparrow\\\\\n", " & \\uparrow & \\end{array}\n", "\\hspace{1cm}\\Longrightarrow\\hspace{1cm}\n", " E=2J\\hspace{1cm}\\begin{array}{ccc} & \\uparrow & \\\\\n", " \\downarrow & \\downarrow & \\uparrow\\\\\n", " & \\uparrow & \\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\Delta E=4J$," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=0\\hspace{1cm}\\begin{array}{ccc} & \\uparrow & \\\\\n", " \\downarrow & \\uparrow & \\uparrow\\\\\n", " & \\downarrow & \\end{array}\n", "\\hspace{1cm}\\Longrightarrow\\hspace{1cm}\n", " E=0\\hspace{1cm}\\begin{array}{ccc} & \\uparrow & \\\\\n", " \\downarrow & \\downarrow & \\uparrow\\\\\n", " & \\downarrow & \\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\Delta E=0$," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=2J\\hspace{1cm}\\begin{array}{ccc} & \\downarrow & \\\\\n", " \\downarrow & \\uparrow & \\uparrow\\\\\n", " & \\downarrow & \\end{array}\n", "\\hspace{1cm}\\Longrightarrow\\hspace{1cm}\n", " E=-2J\\hspace{1cm}\\begin{array}{ccc} & \\downarrow & \\\\\n", " \\downarrow & \\downarrow & \\uparrow\\\\\n", " & \\downarrow & \\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\Delta E=-4J$ and finally" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E=4J\\hspace{1cm}\\begin{array}{ccc} & \\downarrow & \\\\\n", " \\downarrow & \\uparrow & \\downarrow\\\\\n", " & \\downarrow & \\end{array}\n", "\\hspace{1cm}\\Longrightarrow\\hspace{1cm}\n", " E=-4J\\hspace{1cm}\\begin{array}{ccc} & \\downarrow & \\\\\n", " \\downarrow & \\downarrow & \\downarrow\\\\\n", " & \\downarrow & \\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with $\\Delta E=-8J$.\n", "\n", "\n", "## [The Metropolis Algorithm and the Two-dimensional Ising Model, elements of program](https://github.com/CompPhysics/ComputationalPhysics/blob/master/doc/Programs/ParallelizationMPI/MPIising.cpp)\n", "\n", "This means in turn that we could construct an array which contains all values\n", "of $e^{\\beta \\Delta E}$ before doing the Metropolis sampling. Else, we\n", "would have to evaluate the exponential at each Monte Carlo sampling. \n", "For the two-dimensional Ising model there are only five possible values. It is rather easy\n", "to convice oneself that for the one-dimensional Ising model we have only three possible values.\n", "The main part of the Ising model program is shown here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " /* \n", " Program to solve the two-dimensional Ising model \n", " The coupling constant J = 1\n", " Boltzmann's constant = 1, temperature has thus dimension energy\n", " Metropolis sampling is used. Periodic boundary conditions.\n", " */\n", " #include \n", " #include \n", " #include \n", " #include \"lib.h\"\n", " using namespace std;\n", " ofstream ofile;\n", " // inline function for periodic boundary conditions\n", " inline int periodic(int i, int limit, int add) { \n", " return (i+limit+add) % (limit);\n", " }\n", " // Function to read in data from screen \n", " void read_input(int&, int&, double&, double&, double&);\n", " // Function to initialise energy and magnetization\n", " void initialize(int, double, int **, double&, double&);\n", " // The metropolis algorithm \n", " void Metropolis(int, long&, int **, double&, double&, double *);\n", " // prints to file the results of the calculations \n", " void output(int, int, double, double *);\n", " \n", " // main program\n", " int main(int argc, char* argv[])\n", " {\n", " char *outfilename;\n", " long idum;\n", " int **spin_matrix, n_spins, mcs;\n", " double w[17], average[5], initial_temp, final_temp, E, M, temp_step;\n", " \n", " // Read in output file, abort if there are too few command-line arguments\n", " if( argc <= 1 ){\n", " cout << \"Bad Usage: \" << argv[0] << \n", " \" read also output file on same line\" << endl;\n", " exit(1);\n", " }\n", " else{\n", " outfilename=argv[1];\n", " }\n", " ofile.open(outfilename);\n", " // Read in initial values such as size of lattice, temp and cycles\n", " read_input(n_spins, mcs, initial_temp, final_temp, temp_step);\n", " spin_matrix = (int**) matrix(n_spins, n_spins, sizeof(int));\n", " idum = -1; // random starting point\n", " for ( double temp = initial_temp; temp <= final_temp; temp+=temp_step){\n", " // initialise energy and magnetization \n", " E = M = 0.;\n", " // setup array for possible energy changes\n", " for( int de =-8; de <= 8; de++) w[de+8] = 0;\n", " for( int de =-8; de <= 8; de+=4) w[de+8] = exp(-de/temp);\n", " // initialise array for expectation values\n", " for( int i = 0; i < 5; i++) average[i] = 0.;\n", " initialize(n_spins, double temp, spin_matrix, E, M);\n", " // start Monte Carlo computation\n", " for (int cycles = 1; cycles <= mcs; cycles++){\n", " Metropolis(n_spins, idum, spin_matrix, E, M, w);\n", " // update expectation values\n", " average[0] += E; average[1] += E*E;\n", " average[2] += M; average[3] += M*M; average[4] += fabs(M);\n", " }\n", " // print results\n", " output(n_spins, mcs, temp, average);\n", " }\n", " free_matrix((void **) spin_matrix); // free memory\n", " ofile.close(); // close output file\n", " return 0;\n", " }\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Coding energy differences\n", "\n", "The array $w[17]$ contains values of $\\Delta E$ spanning from $-8J$ to\n", "$8J$ and it is precalculated in the main part for every new\n", "temperature. The program takes as input the initial temperature, final\n", "temperature, a temperature step, the number of spins in one direction\n", "(we force the lattice to be a square lattice, meaning that we have the\n", "same number of spins in the $x$ and the $y$ directions) and the number\n", "of Monte Carlo cycles. \n", "\n", "\n", "## Efficient expression for energy differences\n", "\n", "For every Monte Carlo cycle we run through all\n", "spins in the lattice in the function metropolis and flip one spin at\n", "the time and perform the Metropolis test. However, every time we flip\n", "a spin we need to compute the actual energy difference $\\Delta E$ in\n", "order to access the right element of the array which stores $e^{\\beta\n", "\\Delta E}$. This is easily done in the Ising model since we can\n", "exploit the fact that only one spin is flipped, meaning in turn that\n", "all the remaining spins keep their values fixed. The energy\n", "difference between a state $E_1$ and a state $E_2$ with zero external\n", "magnetic field is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\Delta E = E_2-E_1 =J\\sum_{}^{N}s_k^1s_{l}^1-J\\sum_{}^{N}s_k^2s_{l}^2,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which we can rewrite as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\Delta E = -J \\sum_{}^{N}s_k^2(s_l^2-s_{l}^1),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the sum now runs only over the nearest neighbors $k$. \n", "\n", "\n", "## Final energy difference\n", "\n", "Since the spin to be flipped takes only two values, $s_l^1=\\pm 1$ and $s_l^2=\\pm 1$, it means that if\n", "$s_l^1= 1$, then $s_l^2=-1$ and if $s_l^1= -1$, then $s_l^2=1$. \n", "The other spins keep their values, meaning that\n", "$s_k^1=s_k^2$.\n", "If $s_l^1= 1$ we must have $s_l^1-s_{l}^2=2$, and \n", "if $s_l^1= -1$ we must have $s_l^1-s_{l}^2=-2$. From these results we see that the energy difference\n", "can be coded efficiently as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\Delta E = 2Js_l^1\\sum_{}^{N}s_k,\n", "\\label{eq:deltaeising} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where the sum runs only over the nearest neighbors $k$ of spin $l$.\n", "We can compute the change in magnetisation by flipping one spin as well.\n", "Since only spin $l$ is flipped, all the surrounding spins remain unchanged.\n", "\n", "## Change in magnetization\n", "\n", "The difference in magnetisation is therefore only given by the difference \n", "$s_l^1-s_{l}^2=\\pm 2$, or in a more compact way as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "M_2 = M_1+2s_l^2,\n", "\\label{eq:deltamising} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $M_1$ and $M_2$ are the magnetizations before and after the spin flip, respectively. \n", "Eqs. ([3](#eq:deltaeising)) and ([4](#eq:deltamising)) are implemented in the function **metropolis** shown here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " void Metropolis(int n_spins, long& idum, int **spin_matrix, double& E, double&M, double *w)\n", " {\n", " // loop over all spins\n", " for(int y =0; y < n_spins; y++) {\n", " for (int x= 0; x < n_spins; x++){\n", " // Find random position\n", " int ix = (int) (ran1(&idum)*(double)n_spins);\n", " int iy = (int) (ran1(&idum)*(double)n_spins);\n", " int deltaE = 2*spin_matrix[iy][ix]*\n", " \t(spin_matrix[iy][periodic(ix,n_spins,-1)]+\n", " \t spin_matrix[periodic(iy,n_spins,-1)][ix] +\n", " \t spin_matrix[iy][periodic(ix,n_spins,1)] +\n", " \t spin_matrix[periodic(iy,n_spins,1)][ix]);\n", " // Here we perform the Metropolis test\n", " if ( ran1(&idum) <= w[deltaE+8] ) {\n", " \tspin_matrix[iy][ix] *= -1; // flip one spin and accept new spin config\n", " // update energy and magnetization\n", " M += (double) 2*spin_matrix[iy][ix];\n", " E += (double) deltaE;\n", " }\n", " }\n", " }\n", " } // end of Metropolis sampling over spins\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A small note\n", "\n", "Note that we loop over all spins but that we choose the lattice positions $x$ and $y$ randomly.\n", "If the move is accepted after performing the Metropolis test, we update the energy and the magnetisation.\n", "The new values are used to update the averages computed in the main function.\n", "\n", "\n", "## Initialization\n", "\n", "We need also to initialize various variables.\n", "This is done in the function here." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " // function to initialise energy, spin matrix and magnetization\n", " void initialize(int n_spins, double temp, int **spin_matrix, \n", " \t\tdouble& E, double& M)\n", " {\n", " // setup spin matrix and intial magnetization\n", " for(int y =0; y < n_spins; y++) {\n", " for (int x= 0; x < n_spins; x++){\n", " if (temp < 1.5) spin_matrix[y][x] = 1; // spin orientation for the ground state\n", " M += (double) spin_matrix[y][x];\n", " }\n", " }\n", " // setup initial energy\n", " for(int y =0; y < n_spins; y++) {\n", " for (int x= 0; x < n_spins; x++){\n", " E -= (double) spin_matrix[y][x]*\n", " \t(spin_matrix[periodic(y,n_spins,-1)][x] +\n", " \t spin_matrix[y][periodic(x,n_spins,-1)]);\n", " }\n", " }\n", " }// end function initialise\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## [The Metropolis Algorithm and the Two-dimensional Ising Model, elements of program](https://github.com/CompPhysics/ComputationalPhysics/blob/master/doc/Programs/ParallelizationMPI/IsingModel.cpp)\n", "\n", "Here follows an alternative Ising model code using the [Mersenne twister engine](http://www.cplusplus.com/reference/random/mt19937_64/) as described in the c++ \"random class\":\" \"." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " /* \n", " Program to solve the two-dimensional Ising model \n", " with zero external field and no parallelization using the Mersenne twister engine for generating random \n", " numbers.\n", " The coupling constant J = 1\n", " Boltzmann's constant = 1, temperature has thus dimension energy\n", " Metropolis sampling is used. Periodic boundary conditions.\n", " The code needs an output file on the command line and the variables mcs, nspins,\n", " initial temp, final temp and temp step.\n", " Run as\n", " ./executable Outputfile numberof spins number of MC cycles initial temp final temp tempstep\n", " ./test.x Lattice 100 10000000 2.1 2.4 0.01\n", " Compile and link as \n", " c++ -O3 -std=c++11 -Rpass=loop-vectorize -o Ising.x IsingModel.cpp -larmadillo\n", " */\n", " \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " #include \n", " using namespace std;\n", " using namespace arma;\n", " // output file\n", " ofstream ofile;\n", " \n", " // inline function for periodic boundary conditions\n", " inline int periodic(int i, int limit, int add) { \n", " return (i+limit+add) % (limit);\n", " }\n", " // Function to initialise energy and magnetization\n", " void InitializeLattice(int, mat &, double&, double&);\n", " // The metropolis algorithm including the loop over Monte Carlo cycles\n", " void MetropolisSampling(int, int, double, vec &);\n", " // prints to file the results of the calculations \n", " void output(int, int, double, vec);\n", " \n", " // Main program begins here\n", " \n", " int main(int argc, char* argv[])\n", " {\n", " string filename;\n", " int NSpins, MCcycles;\n", " double InitialTemp, FinalTemp, TempStep;\n", " if (argc <= 5) {\n", " cout << \"Bad Usage: \" << argv[0] << \n", " \" read output file, Number of spins, MC cycles, initial and final temperature and tempurate step\" << endl;\n", " exit(1);\n", " }\n", " if (argc > 1) {\n", " filename=argv[1];\n", " NSpins = atoi(argv[2]);\n", " MCcycles = atoi(argv[3]); \n", " InitialTemp = atof(argv[4]);\n", " FinalTemp = atof(argv[5]);\n", " TempStep = atof(argv[6]);\n", " }\n", " // Declare new file name and add lattice size to file name\n", " string fileout = filename;\n", " string argument = to_string(NSpins);\n", " fileout.append(argument);\n", " ofile.open(fileout);\n", " // Start Monte Carlo sampling by looping over T first\n", " for (double Temperature = InitialTemp; Temperature <= FinalTemp; Temperature+=TempStep){\n", " vec ExpectationValues = zeros(5);\n", " // start Monte Carlo computation\n", " MetropolisSampling(NSpins, MCcycles, Temperature, ExpectationValues);\n", " output(NSpins, MCcycles, Temperature, ExpectationValues);\n", " }\n", " ofile.close(); // close output file\n", " return 0;\n", " }\n", " \n", " // function to initialise energy, spin matrix and magnetization\n", " void InitializeLattice(int NSpins, mat &SpinMatrix, double& Energy, double& MagneticMoment)\n", " {\n", " // setup spin matrix and initial magnetization\n", " for(int x =0; x < NSpins; x++) {\n", " for (int y= 0; y < NSpins; y++){\n", " SpinMatrix(x,y) = 1.0; // spin orientation for the ground state\n", " MagneticMoment += (double) SpinMatrix(x,y);\n", " }\n", " }\n", " // setup initial energy\n", " for(int x =0; x < NSpins; x++) {\n", " for (int y= 0; y < NSpins; y++){\n", " Energy -= (double) SpinMatrix(x,y)*\n", " \t(SpinMatrix(periodic(x,NSpins,-1),y) +\n", " \t SpinMatrix(x,periodic(y,NSpins,-1)));\n", " }\n", " }\n", " }// end function initialise\n", " \n", " \n", " // The Monte Carlo part with the Metropolis algo with sweeps over the lattice\n", " void MetropolisSampling(int NSpins, int MCcycles, double Temperature, vec &ExpectationValues)\n", " {\n", " // Initialize the seed and call the Mersenne algo\n", " std::random_device rd;\n", " std::mt19937_64 gen(rd());\n", " // Then set up the uniform distribution for x \\in [[0, 1]\n", " std::uniform_real_distribution distribution(0.0,1.0);\n", " // Allocate memory for spin matrix\n", " mat SpinMatrix = zeros(NSpins,NSpins);\n", " // initialise energy and magnetization \n", " double Energy = 0.; double MagneticMoment = 0.;\n", " // initialize array for expectation values\n", " InitializeLattice(NSpins, SpinMatrix, Energy, MagneticMoment);\n", " // setup array for possible energy changes\n", " vec EnergyDifference = zeros(17); \n", " for( int de =-8; de <= 8; de+=4) EnergyDifference(de+8) = exp(-de/Temperature);\n", " for (int cycles = 1; cycles <= MCcycles; cycles++){\n", " // The sweep over the lattice, looping over all spin sites\n", " for(int x =0; x < NSpins; x++) {\n", " for (int y= 0; y < NSpins; y++){\n", " \tint ix = (int) (distribution(gen)*(double)NSpins);\n", " \tint iy = (int) (distribution(gen)*(double)NSpins);\n", " \tint deltaE = 2*SpinMatrix(ix,iy)*\n", " \t (SpinMatrix(ix,periodic(iy,NSpins,-1))+\n", " \t SpinMatrix(periodic(ix,NSpins,-1),iy) +\n", " \t SpinMatrix(ix,periodic(iy,NSpins,1)) +\n", " \t SpinMatrix(periodic(ix,NSpins,1),iy));\n", " \tif ( distribution(gen) <= EnergyDifference(deltaE+8) ) {\n", " \t SpinMatrix(ix,iy) *= -1.0; // flip one spin and accept new spin config\n", " \t MagneticMoment += (double) 2*SpinMatrix(ix,iy);\n", " \t Energy += (double) deltaE;\n", " \t}\n", " }\n", " }\n", " // update expectation values for local node\n", " ExpectationValues(0) += Energy; ExpectationValues(1) += Energy*Energy;\n", " ExpectationValues(2) += MagneticMoment; \n", " ExpectationValues(3) += MagneticMoment*MagneticMoment; \n", " ExpectationValues(4) += fabs(MagneticMoment);\n", " }\n", " \n", " } // end of Metropolis sampling over spins\n", " \n", " void output(int NSpins, int MCcycles, double temperature, vec ExpectationValues)\n", " {\n", " double norm = 1.0/((double) (MCcycles)); // divided by number of cycles \n", " double E_ExpectationValues = ExpectationValues(0)*norm;\n", " double E2_ExpectationValues = ExpectationValues(1)*norm;\n", " double M_ExpectationValues = ExpectationValues(2)*norm;\n", " double M2_ExpectationValues = ExpectationValues(3)*norm;\n", " double Mabs_ExpectationValues = ExpectationValues(4)*norm;\n", " // all expectation values are per spin, divide by 1/NSpins/NSpins\n", " double Evariance = (E2_ExpectationValues- E_ExpectationValues*E_ExpectationValues)/NSpins/NSpins;\n", " double Mvariance = (M2_ExpectationValues - Mabs_ExpectationValues*Mabs_ExpectationValues)/NSpins/NSpins;\n", " ofile << setiosflags(ios::showpoint | ios::uppercase);\n", " ofile << setw(15) << setprecision(8) << temperature;\n", " ofile << setw(15) << setprecision(8) << E_ExpectationValues/NSpins/NSpins;\n", " ofile << setw(15) << setprecision(8) << Evariance/temperature/temperature;\n", " ofile << setw(15) << setprecision(8) << M_ExpectationValues/NSpins/NSpins;\n", " ofile << setw(15) << setprecision(8) << Mvariance/temperature;\n", " ofile << setw(15) << setprecision(8) << Mabs_ExpectationValues/NSpins/NSpins << endl;\n", " } // end output function\n", " \n", " \n", " \n", " \n", " \n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-dimensional Ising Model, energy per spin and specific heat\n", "\n", "The following Python program, based on the above C++ codes, plots the expectation value of the energy and its fluctuation, that is the specific heat. Both quantities are plotted per spin and genererated for a $20\\times 20$ lattice." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "# Code for the two-dimensional Ising model with periodic boundary conditions\n", "import numpy, sys, math\n", "from matplotlib import pyplot as plt\n", "import numpy as np\n", "\n", "\n", "def periodic (i, limit, add):\n", " \"\"\"\n", " Choose correct matrix index with periodic\n", " boundary conditions\n", "\n", " Input:\n", " - i: Base index\n", " - limit: Highest \\\"legal\\\" index\n", " - add: Number to add or subtract from i\n", " \"\"\"\n", " return (i+limit+add) % limit\n", "\n", "def monteCarlo(temp, size, trials):\n", " \"\"\"\n", " Calculate the energy and magnetization\n", " (\\\"straight\\\" and squared) for a given temperature\n", "\n", " Input:\n", " - temp: Temperature to calculate for\n", " - size: dimension of square matrix\n", " - trials: Monte-carlo trials (how many times do we\n", " flip the matrix?)\n", "\n", " Output:\n", " - E_av: Energy of matrix averaged over trials, normalized to spins**2\n", " - E_variance: Variance of energy, same normalization * temp**2\n", " \"\"\"\n", "\n", " #Setup spin matrix, initialize to ground state\n", " spin_matrix = numpy.zeros( (size,size), numpy.int8) + 1\n", "\n", " #Create and initialize variables\n", " E = 0\n", " E_av = E2_av = 0\n", " \n", " #Setup array for possible energy changes\n", " w = numpy.zeros(17,numpy.float64)\n", " for de in xrange(-8,9,4): #include +8\n", " w[de+8] = math.exp(-de/temp)\n", " \n", " #Calculate initial energy\n", " for j in xrange(size): \n", " for i in xrange(size):\n", " E -= spin_matrix.item(i,j)*\\\n", " (spin_matrix.item(periodic(i,size,-1),j) + spin_matrix.item(i,periodic(j,size,1)))\n", "\n", " #Start metropolis MonteCarlo computation \n", " for i in xrange(trials):\n", " #Metropolis\n", " #Loop over all spins, pick a random spin each time\n", " for s in xrange(size**2):\n", " x = int(numpy.random.random()*size)\n", " y = int(numpy.random.random()*size)\n", " deltaE = 2*spin_matrix.item(x,y)*\\\n", " (spin_matrix.item(periodic(x,size,-1), y) +\\\n", " spin_matrix.item(periodic(x,size,1), y) +\\\n", " spin_matrix.item(x, periodic(y,size,-1)) +\\\n", " spin_matrix.item(x, periodic(y,size,1)))\n", " if numpy.random.random() <= w[deltaE+8]:\n", " #Accept!\n", " spin_matrix[x,y] *= -1\n", " E += deltaE\n", " \n", " #Update expectation values\n", " E_av += E\n", " E2_av += E**2\n", "\n", " E_av /= float(trials);\n", " E2_av /= float(trials);\n", " #Calculate variance and normalize to per-point and temp\n", " E_variance = (E2_av-E_av*E_av)/float(size*size*temp*temp);\n", " #Normalize returned averages to per-point\n", " E_av /= float(size*size);\n", "\n", " return (E_av, E_variance)\n", " \n", " \n", "# Main program\n", "\n", "# values of the lattice, number of Monte Carlo cycles and temperature domain\n", "size = 20\n", "trials = 100000\n", "temp_init = 1.8\n", "temp_end = 2.6\n", "temp_step = 0.1\n", "\n", "\n", "temps = numpy.arange(temp_init,temp_end+temp_step/2,temp_step,float)\n", "Dim = np.size(temps)\n", "energy = np.zeros(Dim)\n", "heatcapacity = np.zeros(Dim) \n", "temperature = np.zeros(Dim)\n", "for temp in temps:\n", " (E_av, E_variance) = monteCarlo(temp,size,trials)\n", " temperature[temp] = temp\n", " energy[temp] = E_av\n", " heatcapacity[temp] = E_variance\n", "plt.figure(1)\n", "plt.subplot(211)\n", "plt.axis([1.8,2.6,-2.0, -1.0])\n", "plt.xlabel(r'Temperature $J/(k_B)$')\n", "plt.ylabel(r'Average energy per spin $E/N$')\n", "plt.plot(temperature, energy, 'b-')\n", "plt.subplot(212)\n", "plt.axis([1.8,2.6, 0.0, 2.0])\n", "plt.plot(temperature, heatcapacity, 'r-')\n", "plt.xlabel(r'Temperature $J/(k_B)$')\n", "plt.ylabel(r'Heat capacity per spin $C_V/N$')\n", "plt.savefig('energycv.pdf')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-dimensional Ising Model and analysis of spin values\n", "\n", "The following python code displays the values of the spins as function of temperature.\n", "The blue color corresponds to spin up states while red represents spin down states. Increasing the temperature as input parameter, see the parameters below, results in a a net magnetization which becomes zero. At low temperatures, the system is highly ordered with essentially only one specific spin value." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# coding=utf-8\n", "#2-dimensional ising model with visualization\n", "\n", "import numpy, sys, math\n", "import pygame\n", "\n", "#Needed for visualize when using SDL\n", "screen = None;\n", "font = None;\n", "BLOCKSIZE = 10\n", "\n", "def periodic (i, limit, add):\n", " \"\"\"\n", " Choose correct matrix index with periodic\n", " boundary conditions\n", "\n", " Input:\n", " - i: Base index\n", " - limit: Highest \\\"legal\\\" index\n", " - add: Number to add or subtract from i\n", " \"\"\"\n", " return (i+limit+add) % limit\n", "\n", "def visualize(spin_matrix, temp, E, M, method):\n", " \"\"\"\n", " Visualize the spin matrix\n", "\n", " Methods:\n", " method = -1:No visualization (testing)\n", " method = 0: Just print it to the terminal\n", " method = 1: Pretty-print to terminal\n", " method = 2: SDL/pygame single-pixel\n", " method = 3: SDL/pygame rectangle\n", " \"\"\"\n", "\n", " #Simple terminal dump\n", " if method == 0:\n", " print \"temp:\", temp, \"E:\", E, \"M:\", M\n", " print spin_matrix\n", " #Pretty-print to terminal\n", " elif method == 1:\n", " out = \"\"\n", " size = len(spin_matrix)\n", " for y in xrange(size):\n", " for x in xrange(size):\n", " if spin_matrix.item(x,y) == 1:\n", " out += \"X\"\n", " else:\n", " out += \" \"\n", " out += \"\\n\"\n", " print \"temp:\", temp, \"E:\", E, \"M:\", M\n", " print out + \"\\n\"\n", " #SDL single-pixel (useful for large arrays)\n", " elif method == 2:\n", " size = len(spin_matrix)\n", " screen.lock()\n", " for y in xrange(size):\n", " for x in xrange(size):\n", " if spin_matrix.item(x,y) == 1:\n", " screen.set_at((x,y),(0,0,255))\n", " else:\n", " screen.set_at((x,y),(255,0,0))\n", " screen.unlock()\n", " pygame.display.flip()\n", " #SDL block (usefull for smaller arrays)\n", " elif method == 3:\n", " size = len(spin_matrix)\n", " screen.lock()\n", " for y in xrange(size):\n", " for x in xrange(size):\n", " if spin_matrix.item(x,y) == 1:\n", " rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)\n", " pygame.draw.rect(screen,(0,0,255),rect)\n", " else:\n", " rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)\n", " pygame.draw.rect(screen,(255,0,0),rect)\n", " screen.unlock()\n", " pygame.display.flip()\n", " #SDL block w/ data-display\n", " elif method == 4:\n", " size = len(spin_matrix)\n", " screen.lock()\n", " for y in xrange(size):\n", " for x in xrange(size):\n", " if spin_matrix.item(x,y) == 1:\n", " rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)\n", " pygame.draw.rect(screen,(255,255,255),rect)\n", " else:\n", " rect = pygame.Rect(x*BLOCKSIZE,y*BLOCKSIZE,BLOCKSIZE,BLOCKSIZE)\n", " pygame.draw.rect(screen,(0,0,0),rect)\n", " s = font.render(\" = %5.3E; = %5.3E\" % E,M,False,(255,0,0))\n", " screen.blit(s,(0,0))\n", " \n", " screen.unlock()\n", " pygame.display.flip() \n", "\n", " \n", "\n", "def monteCarlo(temp, size, trials, visual_method):\n", " \"\"\"\n", " Calculate the energy and magnetization\n", " (\\\"straight\\\" and squared) for a given temperature\n", "\n", " Input:\n", " - temp: Temperature to calculate for\n", " - size: dimension of square matrix\n", " - trials: Monte-carlo trials (how many times do we\n", " flip the matrix?)\n", " - visual_method: What method should we use to visualize?\n", "\n", " Output:\n", " - E_av: Energy of matrix averaged over trials, normalized to spins**2\n", " - E_variance: Variance of energy, same normalization * temp**2\n", " - M_av: Magnetic field of matrix, averaged over trials, normalized to spins**2\n", " - M_variance: Variance of magnetic field, same normalization * temp\n", " - Mabs: Absolute value of magnetic field, averaged over trials\n", " \"\"\"\n", "\n", " #Setup spin matrix, initialize to ground state\n", " spin_matrix = numpy.zeros( (size,size), numpy.int8) + 1\n", "\n", " #Create and initialize variables\n", " E = M = 0\n", " E_av = E2_av = M_av = M2_av = Mabs_av = 0\n", " \n", " #Setup array for possible energy changes\n", " w = numpy.zeros(17,numpy.float64)\n", " for de in xrange(-8,9,4): #include +8\n", " w[de+8] = math.exp(-de/temp)\n", " \n", " #Calculate initial magnetization:\n", " M = spin_matrix.sum()\n", " #Calculate initial energy\n", " for j in xrange(size): \n", " for i in xrange(size):\n", " E -= spin_matrix.item(i,j)*\\\n", " (spin_matrix.item(periodic(i,size,-1),j) + spin_matrix.item(i,periodic(j,size,1)))\n", "\n", " #Start metropolis MonteCarlo computation \n", " for i in xrange(trials):\n", " #Metropolis\n", " #Loop over all spins, pick a random spin each time\n", " for s in xrange(size**2):\n", " x = int(numpy.random.random()*size)\n", " y = int(numpy.random.random()*size)\n", " deltaE = 2*spin_matrix.item(x,y)*\\\n", " (spin_matrix.item(periodic(x,size,-1), y) +\\\n", " spin_matrix.item(periodic(x,size,1), y) +\\\n", " spin_matrix.item(x, periodic(y,size,-1)) +\\\n", " spin_matrix.item(x, periodic(y,size,1)))\n", " if numpy.random.random() <= w[deltaE+8]:\n", " #Accept!\n", " spin_matrix[x,y] *= -1\n", " M += 2*spin_matrix[x,y]\n", " E += deltaE\n", " \n", " #Update expectation values\n", " E_av += E\n", " E2_av += E**2\n", " M_av += M\n", " M2_av += M**2\n", " Mabs_av += int(math.fabs(M))\n", "\n", " visualize(spin_matrix, temp,E/float(size**2),M/float(size**2), method);\n", "\n", " #Normalize average values\n", " E_av /= float(trials);\n", " E2_av /= float(trials);\n", " M_av /= float(trials);\n", " M2_av /= float(trials);\n", " Mabs_av /= float(trials);\n", " #Calculate variance and normalize to per-point and temp\n", " E_variance = (E2_av-E_av*E_av)/float(size*size*temp*temp);\n", " M_variance = (M2_av-M_av*M_av)/float(size*size*temp);\n", " #Normalize returned averages to per-point\n", " E_av /= float(size*size);\n", " M_av /= float(size*size);\n", " Mabs_av /= float(size*size);\n", " \n", " return (E_av, E_variance, M_av, M_variance, Mabs_av)\n", " \n", " \n", "# Main program\n", "\n", "size = 100\n", "trials = 100000\n", "temp = 2.1\n", "method = 3\n", "\n", "#Initialize pygame\n", "if method == 2 or method == 3 or method == 4:\n", " pygame.init()\n", " if method == 2:\n", " screen = pygame.display.set_mode((size,size))\n", " elif method == 3:\n", " screen = pygame.display.set_mode((size*10,size*10))\n", " elif method == 4:\n", " screen = pygame.display.set_mode((size*10,size*10))\n", " font = pygame.font.Font(None,12)\n", " \n", "\n", "(E_av, E_variance, M_av, M_variance, Mabs_av) = monteCarlo(temp,size,trials, method)\n", "print \"%15.8E %15.8E %15.8E %15.8E %15.8E %15.8E\\n\" % (temp, E_av, E_variance, M_av, M_variance, Mabs_av)\n", "\n", "pygame.quit();" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 4 }