BraggWilliams.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "from __future__ import division\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "from numba import njit\n",
    "# Simulating the Bragg-Williams model\n",
    "import numpy as np\n",
    "from numpy.random import rand\n",
    "import matplotlib.pyplot as plt\n",
    "# Interaction parameters\n",
    "waa = -1.\n",
    "wbb = -1.\n",
    "wab = 1.\n",
    "\n",
    "# System size\n",
    "N = 100  # sites in x and y-directions\n",
    "M = N*N # total number of sites\n",
    "\n",
    "#    ''' Simulating the Bragg-Williams model a la Martin Trulsson'''  \n",
    "# Code re-adated from https://github.com/rajeshrinet/compPhy and the Ising model\n",
    "#    ## monte carlo moves\n",
    "@njit\n",
    "def mcmove(config, N, beta):\n",
    "    ''' This is to execute the monte carlo moves using \n",
    "    Metropolis algorithm such that detailed\n",
    "    balance condition is satisified'''\n",
    "\n",
    "    for i in range(N):\n",
    "        for j in range(N):            \n",
    "                x1  = np.random.randint(0, N)\n",
    "                y1  = np.random.randint(0, N)\n",
    "                x2 = np.random.randint(0, N)\n",
    "                y2 = np.random.randint(0, N)\n",
    "                s1 = config[x1, y1]\n",
    "                s2 = config[x2, y2]\n",
    "                if s1!=s2: # only if 1 and 2 is different otherwise does not change state\n",
    "                    # number of a neighbours place 1\n",
    "                    nb1 = config[(x1+1)%N,y1] + config[x1,(y1+1)%N] + config[(x1-1)%N,y1] + config[x1,(y1-1)%N]\n",
    "                    # number of a neighbours place 2 -> What does %N do? \n",
    "                    nb2 = config[(x2+1)%N,y2] + config[x2,(y2+1)%N] + config[(x2-1)%N,y2] + config[x2,(y2-1)%N]\n",
    "                    na1 = 4.-nb1\n",
    "                    na2 = 4.-nb2\n",
    "                    if s1==1: # inital 1:b 2:a center # final 1:a 2:b\n",
    "                        de = -(waa+wbb-2.*wab)*(nb1-nb2)\n",
    "                    else: # a center\n",
    "                        de = (waa+wbb-2.*wab)*(nb1-nb2)\n",
    "                    if de < 0: # energy decreases -> always accept\t\n",
    "                        if s1==1: \n",
    "                            config[x1,y1]=0\n",
    "                            config[x2,y2]=1\n",
    "                        else:\n",
    "                            config[x1,y1]=1\n",
    "                            config[x2,y2]=0\n",
    "                    elif np.random.random() < np.exp(-de*beta): \n",
    "                        if s1==1:\n",
    "                            config[x1,y1]=0\n",
    "                            config[x2,y2]=1\n",
    "                        else:\n",
    "                            config[x1,y1]=1\n",
    "                            config[x2,y2]=0\n",
    "    return config\n",
    "    \n",
    "def simulate():   \n",
    "    ''' This module simulates the Bragg-Williams model'''\n",
    "    # Initialse the lattice (random configuration)\n",
    "    config = np.zeros((N,N))   # Matrix N*N\n",
    "    xmix = 0.5\n",
    "    # In the while loop, a fraction of the lattice sites are assigned to particle B (value set to 1)\n",
    "    i = 0\n",
    "    while i < M*xmix:\n",
    "        x  = np.random.randint(0, N)\n",
    "        y  = np.random.randint(0, N)\n",
    "        if config[x,y]==0:\n",
    "          config[x,y]=1\n",
    "          i=i+1\n",
    "\n",
    "    # Plot data\n",
    "    f = plt.figure(figsize=(15, 15), dpi=80);    \n",
    "    plt.suptitle(\"Temperature=\"+str(temp))\n",
    "    # Plot the first random configuartion\n",
    "    configPlot(f, config, 0, N, 1);\n",
    "    msrmnt = 501\n",
    "    intU = np.empty(0)\n",
    "    for i in range(msrmnt):\n",
    "        mcmove(config, N, 1.0/temp)\n",
    "        if i == 1:       configPlot(f, config, i, N, 2);\n",
    "        if i == 50:       configPlot(f, config, i, N, 3);\n",
    "        if i == 100:      configPlot(f, config, i, N, 4);\n",
    "        if i == 200:     configPlot(f, config, i, N, 5);\n",
    "        if i == 500:     configPlot(f, config, i, N, 6);\n",
    "        # Calculate the total energy\n",
    "        if i in np.arange(200,500,1):\n",
    "            u = 0\n",
    "            for i in range(N):\n",
    "                for j in range(N): \n",
    "                    # Find an expression for the total energy\n",
    "                    u=u+0.\n",
    "            intU = np.append(intU,u)         \n",
    "    return intU.mean()/(M),intU.std()/float(M)\n",
    "                 \n",
    "                    \n",
    "def configPlot(f, config, i, N, n_):\n",
    "    ''' This modules plts the configuration once passed to it along with time etc '''\n",
    "    X, Y = np.meshgrid(range(N), range(N))\n",
    "    sp =  f.add_subplot(3, 3, n_ )  \n",
    "    plt.setp(sp.get_yticklabels(), visible=False)\n",
    "    plt.setp(sp.get_xticklabels(), visible=False)      \n",
    "    plt.imshow(config, interpolation='none',cmap=plt.cm.RdBu)\n",
    "    plt.title('Time=%d'%i)     \n",
    "\n",
    "intUvsT = np.empty(0)\n",
    "intUerrvsT = np.empty(0)\n",
    "intU2vsT = np.empty(0)\n",
    "tempvec = np.arange(4.,6.5,0.5)\n",
    "for temp in tempvec:\n",
    "    # rm = BraggWilliams()\n",
    "    intU_mean, intU_std = simulate()\n",
    "    intUvsT = np.append(intUvsT, intU_mean)\n",
    "    intUerrvsT = np.append(intUerrvsT, intU_std)\n",
    "    \n",
    "plt.show()\n",
    "plt.errorbar(tempvec, intUvsT/tempvec, intUerrvsT/tempvec, lw=0, marker='o',color='black')\n",
    "plt.title('Internal Energy')\n",
    "plt.ylabel('$U$  / $kT$')\n",
    "plt.xlabel('$kT$')\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}