Merge pull request #24 from ev-br/frustr_ising

Add the frustrated 2D Ising example

.github/workflows/python-package.yml

@@ -38,3 +38,7 @@ jobs:
         python setup.py build_ext --inplace
         python -c'from cy_ising import simulate; simulate(4, 1.25, 20000)'
 
+        # also check that the 2D triangular Ising example runs
+        pip install jupyter
+        ipython -c"%run frustrated_2d_ising.ipynb"
+

examples/frustrated_2d_ising.ipynb

@@ -0,0 +1,291 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "b92b380b",
+   "metadata": {},
+   "source": [
+    "# Simulate a frustrated 2D Ising model on the triangular lattice. \n",
+    "\n",
+    "Two exchange integrals: the nn coupling `J > 0` is ferromagnetic, and the coupling across the diagonal (across teh distance of $\\sqrt(2)$) $J_d < 0$ is antiferromagnetic.\n",
+    "\n",
+    "The code is written for the diploma work, https://www.hse.ru/edu/vkr/473082667\n",
+    "and the original code repository is https://github.com/vnsverchkova/VKR.\n",
+    "\n",
+    "\n",
+    "The simulation is a single spin flip Metropolis algorithm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "0bee3da3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext cython"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "da0196e3",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "In file included from /home/br/.virtualenvs/mc_lib/lib/python3.8/site-packages/numpy/core/include/numpy/ndarraytypes.h:1944,\n",
+      "                 from /home/br/.virtualenvs/mc_lib/lib/python3.8/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,\n",
+      "                 from /home/br/.virtualenvs/mc_lib/lib/python3.8/site-packages/numpy/core/include/numpy/arrayobject.h:4,\n",
+      "                 from /home/br/.cache/ipython/cython/_cython_magic_8c6439f2b10737e8c28c1e250f3a47b2.cpp:652:\n",
+      "/home/br/.virtualenvs/mc_lib/lib/python3.8/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning \"Using deprecated NumPy API, disable it with \" \"#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-Wcpp]\n",
+      "   17 | #warning \"Using deprecated NumPy API, disable it with \" \\\n",
+      "      |  ^~~~~~~\n",
+      "/home/br/.cache/ipython/cython/_cython_magic_8c6439f2b10737e8c28c1e250f3a47b2.cpp:2860:15: warning: ‘double __pyx_f_46_cython_magic_8c6439f2b10737e8c28c1e250f3a47b2_energy(__Pyx_memviewslice, __Pyx_memviewslice, __Pyx_memviewslice)’ defined but not used [-Wunused-function]\n",
+      " 2860 | static double __pyx_f_46_cython_magic_8c6439f2b10737e8c28c1e250f3a47b2_energy(__Pyx_memviewslice __pyx_v_spins, __Pyx_memviewslice __pyx_v_neighbors, __Pyx_memviewslice __pyx_v_Js) {\n",
+      "      |               ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%cython -+ \n",
+    "cimport cython\n",
+    "import numpy as np\n",
+    "cimport numpy as np\n",
+    "\n",
+    "from libc.math cimport exp, tanh\n",
+    "from mc_lib.rndm cimport RndmWrapper\n",
+    "from mc_lib.lattices import tabulate_neighbors\n",
+    "from mc_lib.observable cimport RealObservable\n",
+    "\n",
+    "\n",
+    "cdef void init_spins(long[::1] spins, RndmWrapper rndm): \n",
+    "    for j in range(spins.shape[0]):\n",
+    "        spins[j] = 1 if rndm.uniform() > 0.5 else -1\n",
+    "        \n",
+    "        \n",
+    "@cython.boundscheck(False)\n",
+    "@cython.wraparound(False)\n",
+    "cdef double energy(long[::1] spins, \n",
+    "                   long[:, ::1] neighbors,\n",
+    "                   const double[:,::1] Js):\n",
+    "\n",
+    "    cdef:\n",
+    "        double ene = 0.0\n",
+    "        Py_ssize_t site, site1, num_neighb\n",
+    "\n",
+    "    for site in range(spins.shape[0]):\n",
+    "        num_neighb = neighbors[site, 0]\n",
+    "        for j in range(1, num_neighb+1):\n",
+    "            site1 = neighbors[site, j]\n",
+    "            ene += -1 * Js[site, site1] * spins[site] * spins[site1] \n",
+    "    \n",
+    "    return ene / 2.0\n",
+    "\n",
+    "\n",
+    "@cython.boundscheck(False)\n",
+    "@cython.wraparound(False)\n",
+    "cdef void flip_spin(long[::1] spins, \n",
+    "                    const long[:, ::1] neighbors,\n",
+    "                    double beta,\n",
+    "                    const double[:,::1] Js,\n",
+    "                    RndmWrapper rndm): \n",
+    "    cdef:\n",
+    "        Py_ssize_t site = int(spins.shape[0] * rndm.uniform())\n",
+    "        Py_ssize_t site1\n",
+    "\n",
+    "    cdef long num_neighb = neighbors[site, 0]\n",
+    "    cdef double summ = 0.\n",
+    "    for j in range(1, num_neighb + 1):\n",
+    "        site1 = neighbors[site, j]\n",
+    "        summ += spins[site1] * spins[site] * Js[site,site1]\n",
+    "   \n",
+    "    cdef double ratio = exp(-2.0 * beta * summ )\n",
+    "    \n",
+    "    if rndm.uniform() > ratio:\n",
+    "        return\n",
+    "\n",
+    "    spins[site] = -spins[site]\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",
+    "def simulate(Py_ssize_t L,\n",
+    "             double T, double J, double Jd,\n",
+    "             Py_ssize_t num_sweeps, int seed, int rseed = 1234):\n",
+    "\n",
+    "    cdef:\n",
+    "        long[:, ::1] neighbors = tabulate_neighbors(L, kind='triang') \n",
+    "        double beta = 1./T\n",
+    "\n",
+    "    cdef:\n",
+    "        int num_therm = int(30 * L)\n",
+    "        int steps_per_sweep = L * L \n",
+    "        int sweep = 0\n",
+    "        int i\n",
+    "        double Z = 0., magn = 0., binder = 0., error = 0.\n",
+    "        \n",
+    "    \n",
+    "    cdef RndmWrapper rndm = RndmWrapper((rseed, seed)) \n",
+    "    cdef RealObservable m2 = RealObservable()\n",
+    "    cdef RealObservable m4 = RealObservable()\n",
+    "\n",
+    "    cdef long[::1] spins =  np.empty( L*L, dtype=int) \n",
+    "    init_spins(spins, rndm)\n",
+    "    \n",
+    "    cdef double[:,::1] Js = np.zeros((L*L, L*L)) \n",
+    "    get_J(Js, J, Jd, L, L)\n",
+    "\n",
+    "    for sweep in range(num_therm):\n",
+    "        for i in range(steps_per_sweep):\n",
+    "            flip_spin(spins, neighbors, beta, Js, rndm)\n",
+    "\n",
+    "    m = np.zeros(num_sweeps)\n",
+    "\n",
+    "    for sweep in range(num_sweeps):\n",
+    "        for i in range(steps_per_sweep):\n",
+    "            flip_spin(spins, neighbors, beta, Js, rndm)\n",
+    "            \n",
+    "        Z += 1\n",
+    "        magn = 0.\n",
+    "        for i in range(L*L):\n",
+    "            magn += spins[i]\n",
+    "            \n",
+    "        m2.add_measurement(magn**2)\n",
+    "        m4.add_measurement(magn**4)\n",
+    "        m[sweep] = magn\n",
+    "        \n",
+    "    binder = 1 - (m4.mean) / (3 * (m2.mean**2))\n",
+    "    \n",
+    "    error = np.sqrt( ( m4.errorbar/( 3*(m2.mean**2)) )**2 + ( 2*(m4.mean)*m2.errorbar/(3*(m2.mean**3) ) )**2 )\n",
+    "     \n",
+    "    return (binder)#, error, m)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "5628ca52",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.5110854799436115"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "simulate(L=4, T=5, J=1, Jd=1, num_sweeps=10, seed=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2a2d9104",
+   "metadata": {},
+   "source": [
+    "## Use multiprocessing for running several simulations in parallel\n",
+    "\n",
+    "(this can be streamlined, of course)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "id": "181e48c5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from multiprocessing import Pool"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "id": "94b126b6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "param = [\n",
+    "    # L, T, J, Jd, num_sweeps, seed\n",
+    "    (4, 1.0, 1, -0.1, 100, 0),\n",
+    "    (4, 1.5, 1, -0.1, 100, 0),\n",
+    "    (4, 2.0, 1, -0.1, 100, 0),\n",
+    "]\n",
+    "\n",
+    "with Pool(processes=2) as pool:\n",
+    "    res = pool.starmap(simulate, param)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "id": "0d5a40b7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res_sequential = [simulate(*par) for par in param ]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "id": "990cbb4f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from numpy.testing import assert_allclose\n",
+    "\n",
+    "assert_allclose(res,\n",
+    "                res_sequential,\n",
+    "                atol=1e-14)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fed870b8",
+   "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.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}