{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ae524ac",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_formats = ['svg']\n",
    "\n",
    "import numpy as np\n",
    "from numpy.linalg import norm\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['animation.frame_format'] = \"svg\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b6bec48",
   "metadata": {},
   "source": [
    "(dynamics)=\n",
    "\n",
    "# Dynamics\n",
    "\n",
    "As discussed in [](switching), flatspin employs a generalized Stoner-Wohlfarth (GSW) switching model to determine whether a spin will flip its magnetization.\n",
    "A spin may flip if the **total magnetic field** acting on the spin:\n",
    "\n",
    "1. is outside of the switching astroid (left hand side of the GSW equation is greater than 1) **and**\n",
    "2. is oriented in the opposite direction of the spin magnetization ($h_\\parallel < 0$)\n",
    "\n",
    "## Switching energy\n",
    "The **switching energy** captures both these conditions, and can be calculated with {func}`switching_energy() <flatspin.model.SpinIce.switching_energy>`.\n",
    "When the switching energy is positive, the total field is outside the astroid (1) and antiparallel to the spin magnetization (2).\n",
    "The magnitude of the switching energy corresponds to how far the total field is outside the astroid.\n",
    "A spin may flip if, and only if, the switching energy is positive.\n",
    "A list of flippable spins can be obtained by calling {func}`flippable() <flatspin.model.SpinIce.flippable>`.\n",
    "\n",
    "Below we color the spins in a square spin ice by their corresponding switching energies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9414bbeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from flatspin.model import SquareSpinIceClosed\n",
    "\n",
    "model = SquareSpinIceClosed()\n",
    "model.spin[[0, 4, 39]] = -1 # flip spins 0, 4 and 39\n",
    "\n",
    "print(\"Flippable spins:\", model.flippable())\n",
    "\n",
    "# Calculate the switching energy\n",
    "E = model.switching_energy()\n",
    "E_max = np.max(np.abs(E))\n",
    "quiv = model.plot(C=E, cmap='coolwarm', clim=(-E_max, E_max))\n",
    "plt.colorbar(quiv);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbca6902",
   "metadata": {},
   "source": [
    "In the example above, the switching energy was always negative, meaning no spins were flippable.\n",
    "If we increase the external field, we will see that some spins get a positive switching energy and become flippable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41983d61",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase the external field\n",
    "model.set_h_ext([0.2, 0])\n",
    "print(\"Flippable spins:\", model.flippable())\n",
    "\n",
    "# Calculate the switching energy\n",
    "E = model.switching_energy()\n",
    "E_max = np.max(np.abs(E))\n",
    "quiv = model.plot(C=E, cmap='coolwarm', clim=(-E_max, E_max))\n",
    "plt.colorbar(quiv);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ca6ae40",
   "metadata": {},
   "source": [
    "(dynamics-step)=\n",
    "## Flipping spins\n",
    "\n",
    "flatspin employs single spin flip dynamics, where only one spin is flipped at a time.\n",
    "Calling {func}`step() <flatspin.model.SpinIce.step>` advances the simulation one step by flipping a spin.\n",
    "The dipolar fields are always recalculated as part of a simulation step, **but the thermal fields are not**.\n",
    "\n",
    "If there are more than one flippable spin, **the spin with the highest switching energy will flip first**.\n",
    "This ensures that the spin flip order is completely *deterministic*.\n",
    "\n",
    "In the animation below we flip one spin at a time, until there are no more flippable spins.\n",
    "Notice how flipping a spin affects the switching energy of neighboring spins, due to the change in dipolar fields.\n",
    "Consequently, only 3 spins will actually flip, even though there were 6 flippable spins originally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a565f9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "\n",
    "def do_steps():\n",
    "    yield model\n",
    "    while model.step():\n",
    "        yield model\n",
    "\n",
    "def do_plot(model):\n",
    "    E = model.switching_energy()\n",
    "    E_max = np.max(np.abs(E))\n",
    "    model.plot(C=E, cmap='coolwarm', clim=(-E_max, E_max), ax=ax, replace=True)\n",
    "        \n",
    "fig, ax = plt.subplots()\n",
    "anim = FuncAnimation(fig, do_plot, frames=do_steps(), interval=500, blit=False)\n",
    "plt.close() # Only show the animation\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eea248c4",
   "metadata": {},
   "source": [
    "(dynamics-relax)=\n",
    "## Relaxation\n",
    "\n",
    "Flipping all spins until there are no more flippable spins is referred to as a *relaxation*, and is the core dynamical process of flatspin.\n",
    "A simulation is advanced by updating the external and/or thermal fields, and then flipping all flippable spins by calling {func}`relax() <flatspin.model.SpinIce.relax>`.\n",
    "{func}`relax() <flatspin.model.SpinIce.relax>` will re-sample the thermal fields *once*, before starting to flip spins.\n",
    "Because the external and thermal fields are held constant during relaxation, the process is guaranteed to terminate.\n",
    "{func}`relax() <flatspin.model.SpinIce.relax>` returns the number of spins that flipped.\n",
    "\n",
    "Below we use {func}`relax() <flatspin.model.SpinIce.relax>` to simulate the reversal of a square spin ice by an increasing external field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c410c68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Strength of external field\n",
    "H = np.linspace(0, 0.1, 100)\n",
    "# Direction of the external field\n",
    "phi = np.deg2rad(180+45)\n",
    "h_dir = np.array([np.cos(phi), np.sin(phi)])\n",
    "\n",
    "def do_reversal():\n",
    "    yield model\n",
    "    for h in H:\n",
    "        model.set_h_ext(h * h_dir)\n",
    "        if model.relax():\n",
    "            # Only show plot when there were spin flips\n",
    "            yield model\n",
    "\n",
    "def do_plot(model):\n",
    "    model.plot(ax=ax, replace=True)\n",
    "        \n",
    "fig, ax = plt.subplots()\n",
    "anim = FuncAnimation(fig, do_plot, frames=do_reversal(), interval=500, blit=False)\n",
    "plt.close() # Only show the animation\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75e70843",
   "metadata": {},
   "source": [
    "If the direction of the external field is fixed (*and there is no thermal field*), the relaxation process is invariant to the resolution of the external field.\n",
    "In other words, increasing the field strength gradually in 10 steps results in the same sequence of spin flips as increasing the strength in a single step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b431bc11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gradually increase strength of external field in 10 steps\n",
    "model.polarize()\n",
    "H = 0.058\n",
    "for h in np.linspace(0, H, 10):\n",
    "    model.set_h_ext([-h, -h])\n",
    "    model.relax()\n",
    "plt.figure()\n",
    "plt.title(\"After 10 field steps\")\n",
    "model.plot()\n",
    "\n",
    "# Set external field in a single step\n",
    "model.polarize()\n",
    "model.set_h_ext([-H, -H])\n",
    "model.relax()\n",
    "plt.figure()\n",
    "plt.title(\"After 1 field step\")\n",
    "model.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "214e657c",
   "metadata": {},
   "source": [
    "## Hysteresis loops\n",
    "\n",
    "Often we are interested in some average quantity of the ensemble, such as the total magnetization.\n",
    "\n",
    "A hysteresis loop is a plot of the total magnetization, projected along the direction of the external field. Below we gradually increase the strength `H` of the external field, applied at some angle `phi`, and record the total magnetization for each field value `H`. We also enable GPU acceleration to speed up the simulation with `use_opencl=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1058e4d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SquareSpinIceClosed(size=(10,10), use_opencl=True)\n",
    "\n",
    "# H decreases linearly from 0.1 to -0.1, then back to 0.1\n",
    "H = np.linspace(0.1, -0.1, 500)\n",
    "H = np.concatenate([H, -H])\n",
    "\n",
    "# Angle of the external field\n",
    "phi = np.deg2rad(45)\n",
    "h_dir = np.array([np.cos(phi), np.sin(phi)])\n",
    "\n",
    "M_H = []\n",
    "\n",
    "for h in H:\n",
    "    model.set_h_ext(h * h_dir)\n",
    "    model.relax()\n",
    "    # Magnetization projected along field direction\n",
    "    m = model.total_magnetization().dot(h_dir)\n",
    "    M_H.append(m)\n",
    "\n",
    "plt.plot(H, M_H)\n",
    "plt.xlabel(\"H (mT)\")\n",
    "plt.ylabel(r\"M_H (a.u.)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc0089ee",
   "metadata": {},
   "source": [
    "Here is an animation of the hysteresis loop, where changes to the state at each {func}`relax() <flatspin.model.SpinIce.relax>` are shown, omitting steps without change:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b12f5e0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset system back to the polarized state\n",
    "model.polarize()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "def do_hysteresis():\n",
    "    yield model\n",
    "    for h in H:\n",
    "        model.set_h_ext(h * h_dir)\n",
    "        if model.relax():\n",
    "            # Only show frames where any spins flipped\n",
    "            yield model\n",
    "\n",
    "def animate(model):\n",
    "    model.plot(ax=ax, replace=True)\n",
    "        \n",
    "anim = FuncAnimation(fig, animate, frames=do_hysteresis(), interval=200, blit=False)\n",
    "plt.close() # Only show the animation\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41353b54",
   "metadata": {},
   "source": [
    "## Disorder\n",
    "\n",
    "So far the coercive fields of all spins have been identical.\n",
    "If we introduce disorder in the coercive fields, the hysteresis loop changes.\n",
    "Below we can see how disorder causes the reversal process to become more gradual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27a80a9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "for disorder in [0, 0.05, 0.10]:\n",
    "    model = SquareSpinIceClosed(size=(10,10), use_opencl=True, disorder=disorder)\n",
    "\n",
    "    # H, phi and h_dir as before\n",
    "    M_H = []\n",
    "    for h in H:\n",
    "        model.set_h_ext(h * h_dir)\n",
    "        model.relax()\n",
    "        m = model.total_magnetization().dot(h_dir)\n",
    "        M_H.append(m)\n",
    "\n",
    "    plt.plot(H, M_H, label=f\"{disorder * 100}% disorder\")\n",
    "\n",
    "plt.xlabel(\"H (mT)\")\n",
    "plt.ylabel(r\"M_H (a.u.)\");\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91b9863e",
   "metadata": {},
   "source": [
    "## Temperature\n",
    "\n",
    "If we increase the temperature, spins flip more easily due to additional thermal energy. Below we can see how temperature affects the hysteresis loop from earlier. Notice how the hysteresis loop becomes narrower at higher temperatures, as the additional thermal energy causes spins to flip for weaker fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06d4680f",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = SquareSpinIceClosed(size=(10,10), m_therm=.5e-17, use_opencl=True)\n",
    "\n",
    "# H, phi and h_dir as before\n",
    "\n",
    "temperatures = [0, 200, 400, 600]\n",
    "\n",
    "for T in temperatures:\n",
    "    M_H = []\n",
    "\n",
    "    model.polarize()\n",
    "    model.set_temperature(T)\n",
    "    for h in H:\n",
    "        model.set_h_ext(h * h_dir)\n",
    "        model.relax()\n",
    "        m = model.total_magnetization().dot(h_dir)\n",
    "        M_H.append(m)\n",
    "\n",
    "    plt.plot(H, M_H, label=f\"{T}K\")\n",
    "\n",
    "plt.xlabel(\"H (mT)\")\n",
    "plt.ylabel(r\"M_H (a.u.)\");\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "md:myst",
   "text_representation": {
    "extension": ".md",
    "format_name": "myst",
    "format_version": 0.13,
    "jupytext_version": "1.11.5"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.6"
  },
  "source_map": [
   15,
   24,
   45,
   58,
   63,
   73,
   89,
   107,
   120,
   142,
   147,
   165,
   173,
   196,
   200,
   220,
   228,
   245,
   251
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 5
}