diff --git a/examples/2D-Ising.ipynb b/examples/2D-Ising.ipynb new file mode 100644 index 0000000..ea85963 --- /dev/null +++ b/examples/2D-Ising.ipynb @@ -0,0 +1,344 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Simulation a 2 dimensional Ising model on square lattice with or without frustrations\n", + "The simulation is a single spin flip Monte Carlo Metropolis algorithm.\n", + "To start the simulation call IsingSimulate() function\n", + "\n", + "MCM steps - Monte Carlo Metropolis steps\n", + "```\n", + "IsingSimulate() function\n", + "parameters:\n", + " L - linear lattice size\n", + " T - temperature value\n", + " sweeps - number of MCM steps\n", + " T_corr - number of MCM steps between configuration records in the output array\n", + " J - value of spins neigbours interaction\n", + " Jd - value of diagonal spins neigbours interaction\n", + " seed - seed for mc_lib rndm initialisation\n", + " rseed - seedSequence for mc_lib rndm initialisation\n", + "\n", + "outputs:\n", + " configs\n", + " e - mean energy for the whole simulation\n", + " m - mean magnetisation\n", + "```\n", + "Example [https://colab.research.google.com/drive/1Xwsx39wkBaopYTNswiy2OenMIx5lYBLZ?usp=sharing](https://colab.research.google.com/drive/1Xwsx39wkBaopYTNswiy2OenMIx5lYBLZ?usp=sharing)\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Check if mc_lib installed\n", + "try:\n", + " from mc_lib.rndm import RndmWrapper\n", + "except ImportError:\n", + " !pip install git+https://github.com/ev-br/mc_lib.git@master\n", + " pass" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%load_ext cython" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%%cython --cplus\n", + "import cython\n", + "import numpy as np\n", + "cimport numpy as np\n", + "from libc.math cimport exp\n", + "from mc_lib.rndm cimport RndmWrapper\n", + "from mc_lib.lattices import tabulate_neighbors\n", + "from mc_lib.observable cimport RealObservable\n", + "\n", + "@cython.boundscheck(False) # turn off bounds-checking for entire function\n", + "@cython.wraparound(False) # turn off negative index wrapping for entire function\n", + "cdef np.ndarray initState(long[::1] lattice,\n", + " RndmWrapper rndm):\n", + "\"\"\"\n", + "\n", + ":param lattice: unconfigured array of future lattice (empty array)\n", + ":param rndm: class with random methods from mc_lib\n", + ":return: change the 'Lattice' variable and return nothing\n", + "\"\"\"\n", + "for i in range(lattice.shape[0]):\n", + " lattice[i] = 1 if rndm.uniform() > 0.5 else -1\n", + "return\n", + "\n", + "\n", + "@cython.boundscheck(False)\n", + "@cython.wraparound(False)\n", + "cdef np.ndarray mcmove(long[::1] config,\n", + " double beta,\n", + " long[:, ::1] ngb,\n", + " RndmWrapper rndm,\n", + " const double[:,::1] Js):\n", + "\"\"\"\n", + "One flip attempt\n", + ":param config: Current configuration of lattice\n", + ":param beta: Inversed temperature of current configuration\n", + ":param L: Linear size 'L' of lattice\n", + ":param ngb: Array of neigbours\n", + "\"\"\"\n", + "cdef:\n", + "Py_ssize_t site = int(config.shape[0] * rndm.uniform())\n", + "Py_ssize_t site1 = 0\n", + "double dE = 0\n", + "long num_ngb = ngb[site, 0]\n", + "for n in range(1, num_ngb + 1):\n", + " site1 = ngb[site, n]\n", + " dE += config[site1] * config[site] * Js[site,site1]\n", + "cdef double ratio = exp(-2 * dE * beta)\n", + "if rndm.uniform() > ratio:\n", + " return\n", + "config[site] *= -1\n", + "return\n", + "\n", + "\n", + "@cython.boundscheck(False)\n", + "@cython.wraparound(False)\n", + "cdef double calcEnergy(long[::1] config,\n", + " long[:, ::1] ngb,\n", + " const double[:,::1] Js):\n", + "\"\"\"\n", + " Count current configuration energy\n", + " :param config: Current configuration of lattice\n", + " :param ngb: Array of neigbours\n", + " return:\n", + "\"\"\"\n", + "cdef:\n", + "Py_ssize_t site = 0\n", + "Py_ssize_t site1 = 0\n", + "double energy = 0\n", + "for site in range(config.shape[0]):\n", + " num_ngb = ngb[site, 0]\n", + " for n in range(1, num_ngb+1):\n", + " site1 = ngb[site, n]\n", + " energy += -1 * config[site] * config[site1] * Js[site, site1]\n", + "return energy / 4.\n", + "\n", + "\n", + "def get_J( double[:,::1] Js, double J, double Jd, int L1, int L2):\n", + "\"\"\"Tabulate the couplings per site.\"\"\"\n", + "cdef Py_ssize_t i\n", + "for i in range(L1*L2):\n", + " Js[i, ((i // L2 + 1) % L1 * L2 ) + (i + 1) % L2 ] = Jd\n", + " Js[i, ((i // L2 - 1) % L1 * L2 ) + (i - 1) % L2 ] = Jd\n", + " Js[i, (i // L2) * L2 + (i + 1) % L2] = J\n", + " Js[i, (i + L2) % (L1*L2)] = J\n", + " Js[i, (i // L2) * L2 + (i - 1) % L2] = J\n", + " Js[i, (i - L2) % (L1*L2)] = J\n", + "return\n", + "\n", + "\n", + "@cython.boundscheck(False)\n", + "@cython.wraparound(False)\n", + "def IsingSimulate(L, T, sweeps, int T_corr, double J=1, double Jd=0, int seed = np.random.randint(0, 1000), int rseed = 1234):\n", + "\"\"\"\n", + " L - linear_size\n", + " T - One temperature point\n", + " J - value of spins neigbours interaction\n", + " Jd - value of diagonal spins neigbours interaction\n", + " Sweeps - number of L^2 Metropolis Monte-Karlo steps\n", + "\"\"\"\n", + "cdef RndmWrapper rndm = RndmWrapper((rseed, seed))\n", + "cdef:\n", + "float beta = 1.0 / T\n", + "int sweep = 0\n", + "int steps_per_sweep = L * L\n", + "int num_therm = int(12 * L * L)\n", + "int i = 0\n", + "\n", + "cdef:\n", + "np.ndarray configs = np.empty([sweeps, L*L], dtype=int)\n", + "long[::1] config = np.empty(L*L, dtype=int)\n", + "int steps = 0\n", + "long[:, ::1] ngb\n", + "if Jd == 0:\n", + " ngb = tabulate_neighbors((L, L, 1), kind='sc')\n", + "else:\n", + " ngb = tabulate_neighbors(L, kind='triang')\n", + "\n", + "cdef double[:,::1] Js = np.zeros((L*L, L*L))\n", + "get_J(Js, J, Jd, L, L)\n", + "\n", + "cdef RealObservable m1 = RealObservable()\n", + "cdef RealObservable e1 = RealObservable()\n", + "\n", + "initState(config, rndm)\n", + "for sweep in range(num_therm):\n", + " mcmove(config, beta, ngb, rndm, Js)\n", + "\n", + "\n", + "for sweep in range(sweeps):\n", + " for i in range(T_corr):\n", + " steps += 1\n", + " for i in range(steps_per_sweep):\n", + " mcmove(config, beta, ngb, rndm, Js)\n", + " configs[sweep] = config.copy()\n", + " energ = calcEnergy(config, ngb, Js)\n", + " magn = np.sum(config)\n", + " e1.add_measurement(energ)\n", + " m1.add_measurement(magn)\n", + "\n", + "return configs, e1.mean, m1.mean" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "IsingSimulate(L=4, T=2.26, sweeps=100, T_corr=1, J=1, Jd=0, seed=10, rseed=10)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Energy and magnetisation graphs" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%%time\n", + "import numpy as np\n", + "ms = []\n", + "es = []\n", + "configs = []\n", + "L = 20\n", + "Temps = np.linspace(1.8, 3, 30)\n", + "for i, t in enumerate(Temps):\n", + " config, e, m = IsingSimulate(L=L, T=t, sweeps=10000, T_corr=50, J=1, Jd=0, seed=0)\n", + " es.append(e)\n", + " ms.append(m)\n", + " configs.append(config)\n", + " if i%10 == 0:\n", + " print(t)\n", + "es = np.array(es)\n", + "ms = np.array(ms)\n", + "configs = np.array(configs)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "try:\n", + " import matplotlib.pyplot as plt\n", + " fig = plt.figure(figsize=(7, 5))\n", + " fontsize = 16\n", + " T_c = 2 / (np.log(1 + 2 ** 0.5))\n", + " plt.scatter(Temps, es / L**2, marker = \"o\", linestyle=\"None\", label='E(T)')\n", + " plt.axvline(T_c, color=\"tab:orange\")\n", + " plt.legend(fontsize=fontsize)\n", + " plt.xlabel(\"T\", fontsize=fontsize)\n", + " plt.ylabel(r\"E(t) / $L^2$\", fontsize=fontsize)\n", + " plt.grid()\n", + " fig = plt.figure(figsize=(7, 5))\n", + " plt.scatter(Temps, np.abs(ms / L**2), marker = \"o\", linestyle=\"None\", label='M(T)')\n", + " plt.axvline(T_c, color=\"tab:orange\")\n", + " plt.legend(fontsize=fontsize)\n", + " plt.xlabel(\"T\", fontsize=fontsize)\n", + " plt.ylabel(r\"M(t) / $L^2$\", fontsize=fontsize)\n", + " plt.grid()\n", + "except ImportError:\n", + " print('You need matplotlib be installed')\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file