3 Dimensional Ising Model for Dr. Engelhardt's meeting (sep 1) .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "e75a5740",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "from scipy import random\n",
    "import matplotlib.pyplot as plt\n",
    "import numba\n",
    "from numba import njit\n",
    "from scipy.ndimage import convolve, generate_binary_structure\n",
    "import numpy as np\n",
    "from mpl_toolkits import mplot3d\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-poster')\n",
    "\n",
    "def energy_ising_2d(spin_arr, J):\n",
    "    # applies the nearest neighbours summation\n",
    "    kern = generate_binary_structure(3, 1) \n",
    "    kern[1][1][1] = False\n",
    "    arr = -J *spin_arr * convolve(spin_arr, kern, mode='constant', cval=0)\n",
    "    return arr.sum()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "c2c27345",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_226723/840150636.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    118\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    119\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 120\u001b[0;31m \u001b[0maverage_spins\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcurrent_energys\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconfigurationq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmetropolis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlattice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_lattice_sites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_lattice_sites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0menergy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0menergy_ising_2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspin_arr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlattice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_lattice_sites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mJ\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdebug\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    121\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"-------------------------------------------------------------------------------------\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    122\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_226723/840150636.py\u001b[0m in \u001b[0;36mmetropolis\u001b[0;34m(spin_arr, n_steps, n_lattice_sites, beta, energy, debug, save_freq)\u001b[0m\n\u001b[1;32m     38\u001b[0m         \u001b[0;31m# 3 / 4. change state with designated probabilities\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     39\u001b[0m         \u001b[0mdE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mE_f\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mE_i\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 40\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdE\u001b[0m\u001b[0;34m>\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mbeta\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdE\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     41\u001b[0m             \u001b[0mspin_arr\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mspin_f\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     42\u001b[0m             \u001b[0menergy\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mdE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
     ]
    }
   ],
   "source": [
    "random_seed = 1\n",
    "np.random.seed(random_seed)\n",
    "def metropolis(spin_arr, n_steps, n_lattice_sites, beta, energy, debug=False,save_freq=10):\n",
    "    average_spins = []\n",
    "    current_energys = []\n",
    "    spin_arr = spin_arr.copy()\n",
    "    net_spins = np.zeros(n_steps-1)\n",
    "    net_energy = np.zeros(n_steps-1)\n",
    " \n",
    "    for t in range(0,n_steps-1):\n",
    "        # 2. pick random point on array and flip spin\n",
    "        x = np.random.randint(0,n_lattice_sites)\n",
    "        y = np.random.randint(0,n_lattice_sites)\n",
    "        z = np.random.randint(0,n_lattice_sites)\n",
    "        spin_i = spin_arr[x,y,z] #initial spin\n",
    "        spin_f = spin_i*-1 #proposed spin flip\n",
    "        # compute change in energy\n",
    "        E_i = 0\n",
    "        E_f = 0\n",
    "        if x>0:\n",
    "            E_i += -spin_i*spin_arr[x-1,y,z]\n",
    "            E_f += -spin_f*spin_arr[x-1,y,z]\n",
    "        if x<n_lattice_sites-1:\n",
    "            E_i += -spin_i*spin_arr[x+1,y,z]\n",
    "            E_f += -spin_f*spin_arr[x+1,y]\n",
    "        if y>0:\n",
    "            E_i += -spin_i*spin_arr[x,y-1,z]\n",
    "            E_f += -spin_f*spin_arr[x,y-1,z]\n",
    "        if y<n_lattice_sites-1:\n",
    "            E_i += -spin_i*spin_arr[x,y+1,z]\n",
    "            E_f += -spin_f*spin_arr[x,y+1,z]\n",
    "        if z>0:\n",
    "            E_i += -spin_i*spin_arr[x,y,z-1]\n",
    "            E_f += -spin_f*spin_arr[x,y,z-1]\n",
    "        if z<n_lattice_sites-1:\n",
    "            E_i += -spin_i*spin_arr[x,y,z+1]\n",
    "            E_f += -spin_f*spin_arr[x,y,z+1]\n",
    "        # 3 / 4. change state with designated probabilities\n",
    "        dE = E_f-E_i\n",
    "        if (dE>0)*(np.random.random() < np.exp(-beta*dE)):\n",
    "            spin_arr[x,y]=spin_f\n",
    "            energy += dE\n",
    "        elif dE<=0:\n",
    "            spin_arr[x,y]=spin_f\n",
    "            energy += dE\n",
    "        net_spins = spin_arr.mean()  \n",
    "        net_energy = energy\n",
    "        if t%save_freq == 0:\n",
    "            average_spins.append(net_spins) \n",
    "            current_energys.append(net_energy)\n",
    "            average_spin = spin_arr.mean()  \n",
    "    return average_spins, current_energys, configurationq\n",
    "\n",
    "\n",
    "\n",
    "####################\n",
    "\n",
    "###########\n",
    "random_seed = 1\n",
    "np.random.seed(random_seed)\n",
    "def metropoliss(spin_arr, n_steps, n_lattice_sites, beta, energy, debug=False):\n",
    "    average_spins = []\n",
    "    current_energys = []\n",
    "    spin_arr = spin_arr.copy()\n",
    "    net_spins = np.zeros(n_steps-1)\n",
    "    net_energy = np.zeros(n_steps-1)\n",
    "    \n",
    "\n",
    "    for t in range(0,n_steps-1):\n",
    "        # 2. pick random point on array and flip spin\n",
    "        x = np.random.randint(0,n_lattice_sites)\n",
    "        y = np.random.randint(0,n_lattice_sites)\n",
    "        z = np.random.randint(0,n_lattice_sites)\n",
    "        spin_i = spin_arr[x,y,z] #initial spin\n",
    "        spin_f = spin_i*-1 #proposed spin flip\n",
    "        # compute change in energy\n",
    "        E_i = 0\n",
    "        E_f = 0\n",
    "        if x>0:\n",
    "            E_i += -spin_i*spin_arr[x-1,y,z]\n",
    "            E_f += -spin_f*spin_arr[x-1,y,z]\n",
    "        if x<n_lattice_sites-1:\n",
    "            E_i += -spin_i*spin_arr[x+1,y,z]\n",
    "            E_f += -spin_f*spin_arr[x+1,y]\n",
    "        if y>0:\n",
    "            E_i += -spin_i*spin_arr[x,y-1,z]\n",
    "            E_f += -spin_f*spin_arr[x,y-1,z]\n",
    "        if y<n_lattice_sites-1:\n",
    "            E_i += -spin_i*spin_arr[x,y+1,z]\n",
    "            E_f += -spin_f*spin_arr[x,y+1,z]\n",
    "        if z>0:\n",
    "            E_i += -spin_i*spin_arr[x,y,z-1]\n",
    "            E_f += -spin_f*spin_arr[x,y,z-1]\n",
    "        if z<n_lattice_sites-1:\n",
    "            E_i += -spin_i*spin_arr[x,y,z+1]\n",
    "            E_f += -spin_f*spin_arr[x,y,z+1]\n",
    "        # 3 / 4. change state with designated probabilities\n",
    "        dE = E_f-E_i\n",
    "        if (dE>0)*(np.random.random() < np.exp(-beta*dE)):\n",
    "            spin_arr[x,y]=spin_f\n",
    "            energy += dE\n",
    "        elif dE<=0:\n",
    "            spin_arr[x,y]=spin_f\n",
    "            energy += dE\n",
    "        net_spins = spin_arr.mean()  \n",
    "        net_energy = energy\n",
    "        average_spins.append(net_spins) \n",
    "        current_energys.append(net_energy)\n",
    "        average_spin = spin_arr.mean()  \n",
    "        configurationq = spin_arr    \n",
    "    return average_spins, current_energys, configurationq\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "###########################################################################################################\n",
    "\n",
    "\n",
    "average_spins, current_energys, configurationq = metropolis(lattice(n_lattice_sites=10), n_steps=100, n_lattice_sites=10, beta=0.1, energy=energy_ising_2d(spin_arr=lattice(n_lattice_sites=10), J=1), debug=True)\n",
    "print(\"-------------------------------------------------------------------------------------\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "58506fd9",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_226723/4189027857.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0mtest_n_steps\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m100000\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0maverage_spins\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcurrent_energys\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconfigurationq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmetropoliss\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlattice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_lattice_sites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtest_n_lattice_sites\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtest_n_steps\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_lattice_sites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtest_n_lattice_sites\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtest_beta\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0menergy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0menergy_ising_2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspin_arr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlattice\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_lattice_sites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtest_n_lattice_sites\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mJ\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtest_J\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdebug\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfigurationq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"gray\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtitle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Lattice configuration in %d th MC step\"\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtest_n_steps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/tmp/ipykernel_226723/840150636.py\u001b[0m in \u001b[0;36mmetropoliss\u001b[0;34m(spin_arr, n_steps, n_lattice_sites, beta, energy, debug)\u001b[0m\n\u001b[1;32m     97\u001b[0m         \u001b[0;31m# 3 / 4. change state with designated probabilities\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     98\u001b[0m         \u001b[0mdE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mE_f\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mE_i\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdE\u001b[0m\u001b[0;34m>\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mbeta\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdE\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m             \u001b[0mspin_arr\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mspin_f\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m             \u001b[0menergy\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mdE\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "test_n_lattice_sites = 50\n",
    "test_beta = 10\n",
    "test_J = 1\n",
    "\n",
    "\n",
    "\n",
    "test_n_steps = 100000\n",
    "\n",
    "average_spins, current_energys, configurationq = metropoliss(lattice(n_lattice_sites=test_n_lattice_sites), n_steps=test_n_steps, n_lattice_sites=test_n_lattice_sites, beta=test_beta, energy=energy_ising_2d(spin_arr=lattice(n_lattice_sites=test_n_lattice_sites), J=test_J), debug=True)\n",
    "plt.imshow(configurationq,cmap=\"gray\")\n",
    "plt.title(\"Lattice configuration in %d th MC step\"%(test_n_steps))\n",
    "plt.show()\n",
    "print(\"-------------------------------------------------------------------------------------\")\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(average_spins)\n",
    "plt.ylabel(\"$<\\sigma> $ Spin Ave.\")\n",
    "plt.xlabel(\"MC Step\")\n",
    "plt.title(\"$\\\\beta=%.2f$,  $J=%.2f$\"%(test_beta,test_J))\n",
    "plt.show()\n",
    "\n",
    "print(\"-------------------------------------------------------------------------------------\")\n",
    "\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(current_energys)\n",
    "plt.ylabel(\"Energy\")\n",
    "plt.xlabel(\"MC Step\")\n",
    "plt.title(\"$\\\\beta=%.2f$,  $J=%.2f$\"%(test_beta,test_J))\n",
    "plt.show()\n",
    "print(\"-------------------------------------------------------------------------------------\")\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a84c765c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f13df6e3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a73fa82-f005-4af0-9409-ef796073257b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "NERSC Python",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}