{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "a01bc831", "metadata": { "tags": [ "remove-input" ] }, "outputs": [], "source": [ "%config InlineBackend.figure_formats = ['svg']" ] }, { "cell_type": "markdown", "id": "01981494", "metadata": {}, "source": [ "(annealing)=\n", "\n", "# Thermal annealing\n", "\n", "## Domain size in square ASI\n", "\n", "In section V. II. of the \"flatspin: A Large-Scale Artificial Spin Ice Simulator\" paper, we briefly describe the use of flatspin to perform an annealing procedure to match an experimental setup.\n", "The following is the setup for flatspin's parameters and a guide to how the experiment was performed.\n", "This example demonstrates a large part of flatspin's features, and most notably a novel use of its temperature implementation. \n", "\n", "The main goal of this experiment is to reproduce the annealing of square artificial spin ice, as in {cite}`Zhang2013`\n", "\n", "[“Crystallites of magnetic charges in artificial spin ice” Nature, vol. 500, no. 7464, 553–557 (2013)](https://doi.org/10.1038/nature12399).\n", "\n", "There are six main parts to this experiment:\n", "\n", "1. Define the systems temperature dependent parameters\n", "2. Define the probablities of flipping once (and more) for a given simulation step\n", "3. Use the probabilities and temperature dependent parameters to identify the nucleation volume\n", "4. Derive an annealing profile and generate the parameters.\n", "5. Produce the flatspin run command and input files\n", "6. Results" ] }, { "cell_type": "markdown", "id": "fa1d5ad8", "metadata": {}, "source": [ "## 1. Temperature dependent parameters\n", "\n", "At high temperatures, the physical magnetic materials constituting the nanomagnet dipoles change their physical properties. While these material properties are not directly modeled by flatspin, we can change relevant parameters of the simulation to fit the temperature induced change. \n", "\n", "**Saturation magnetization, $M_\\text{S}(T)$** \n", "One of the most important changes to a magnetic material as you heat it is the reduction of saturation magnetization. We use values obtained for a 25 nm permalloy film to find an expression for $M_\\text{S}(T)$.\n", "The {download}`experimental values ` are gathered from {cite}`Zhang2019`." ] }, { "cell_type": "code", "execution_count": null, "id": "c09830a0", "metadata": {}, "outputs": [], "source": [ "import json\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy import optimize\n", "from functools import partial\n", "\n", "# Some constant we will use later\n", "mu0 = 1.256637e-6 # kg m s^-2 A^-2\n", "kB = 1.38064852e-23 # m^2 kg s^-2 K^-1 \n", "\n", "# Open the data from Zhang2019\n", "with open('Zhang-Py-25nm-msat-vs-T.json') as f:\n", " msat_data = json.load(f) \n", "dataset_zfc = [dataset['data'] for dataset in msat_data['datasetColl'] if dataset['name']=='ZFC'][0]\n", "data = [p['value'] for p in dataset_zfc] \n", "x, y = list(zip(*data))\n", "y = [val*1e3 for val in y] # Convert to units of A/m\n", "# Plot it\n", "plt.figure(dpi=100)\n", "plt.scatter(x, y, c = 'r', s = 3, label='Extracted data')\n", "\n", "# Fit the following function to the msat curve\n", "def f1(x, a, b, c): \n", " return -a*np.log(-b/(x-c))\n", "a1, b1, c1 = 150000, 10, 820 # Some rough fitting parameters as initial guess\n", "(a0, b0, c0), _ = optimize.curve_fit(f1, x, y, p0=(a1, b1, c1))\n", "x_ax = np.linspace(min(x), max(x))\n", "plt.plot(x_ax, f1(x_ax, a0, b0, c0), label = '$M_S(T)$ Curve fit')\n", "\n", "plt.xlabel('Temperatue [K]')\n", "plt.ylabel('Saturation magnetization [A/m]')\n", "\n", "# Create the msat(T) function, units A/m (not kA/m, as in Zhang et al.)\n", "msat_vs_T = partial(f1, a=a0, b=b0, c=c0)\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "e48ab15a", "metadata": {}, "source": [ "**Coupling strength, $\\alpha$** \n", "While $M_\\text{S}(T)$ is not a direct part of flatspin, it influences other parameters relevant for flatspin simulations. One of these are the coupling parameter, $\\alpha$, which we can define given an $M_\\text{S}(T)$, a volume and a lattice spacing. \n", "\n", "**Switching threshold, $h_c$** \n", "As the saturation magnetization varies, so does the switching threshold, $h_c$. \n", "We create a function to estimate the switching threshold from the temperature, using our `msat_vs_T` function. \n", "*The coefficients of the scaling is taken from micromagnetic simulations.*" ] }, { "cell_type": "code", "execution_count": null, "id": "2f4d6947", "metadata": {}, "outputs": [], "source": [ "# Calculate alpha (in Tesla)\n", "def alpha_calc(msat, V, a):\n", " return mu0 * msat * V / (4 * np.pi * a**3)\n", "\n", "# Calculate switching threshold \n", "def hc_vs_T(T, msat=None):\n", " coeff = np.array([2.36501465e-07, -8.52767576e-03]) \n", " if msat is None:\n", " msat = msat_vs_T(T)\n", " return np.polyval(coeff, msat)" ] }, { "cell_type": "markdown", "id": "8097d8ac", "metadata": {}, "source": [ "Below is a plot of the calculated $\\alpha$ values and $h_c$ values for magnets of size $220\\times80\\times25$ nm and a separation of $320$ nm." ] }, { "cell_type": "code", "execution_count": null, "id": "16528f19", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "T_range = np.linspace(300,820,500)\n", "volume = 220e-9 * 80e-9 * 25e-9\n", "separation = 320e-9\n", "\n", "# Plot alpha\n", "plt.figure(dpi=100)\n", "plt.plot(T_range, alpha_calc(msat_vs_T(T_range), volume, separation), c='C0', ls='-', linewidth=2, label='alpha')\n", "plt.ylabel('Coupling coefficient, alpha, [T]')\n", "plt.xlabel('Temperature [K]')\n", "plt.legend(loc='upper right', bbox_to_anchor=(1, 1.0))\n", "\n", "# Plot hc\n", "ax2 = plt.gca().twinx()\n", "ax2.plot(T_range, hc_vs_T(T_range), c='C1', ls=':', lw=4, label='hc')\n", "plt.ylabel('Switching threshold, hc [T]')\n", "plt.legend(loc='upper right', bbox_to_anchor=(1, 0.9));" ] }, { "cell_type": "markdown", "id": "79d34bfc", "metadata": {}, "source": [ "## 2. Flipping probabilities\n", "\n", "We now use the temperature dependent parameters to define functions returning the probability of switching in a given scenario. \n", "The probability is evaluated for an unbiased magnet, i.e., we consider the shortest distance to the astroid edge from the origin.\n", "To calculate this, we use an example model instance with the same parameters as in the simulation we will be doing, and use its `_zero_to_astroid` attribute. \n", "We use the Poisson formula to calculate the probability.\n", "\n", "The `_zero_to_astroid` attribute returns the smallest field required to flip a magnet, **in units of $h_c$**, which we can use to multiply by `hc_vs_T(T)` to get the smallest field required to flip an unbiased magnet. Thus, only the astroid parameters, `sw_b`, `sw_c`, `sw_beta`, and `sw_gamma`, are relevant for the example model instance. \n", "\n", "We define the probability of flipping *at least once*, `P_flip`, and the probability of flipping *more than once*, `P_multiflip`." ] }, { "cell_type": "code", "execution_count": null, "id": "be86881f", "metadata": {}, "outputs": [], "source": [ "from flatspin.model import SquareSpinIceClosed\n", "\n", "# Example model instance\n", "model = SquareSpinIceClosed(sw_b=0.38, sw_c=1.0, sw_beta=1.3, sw_gamma=3.6)\n", "\n", "def P_flip(T, delta_t, volume_nuc):\n", " threshold_zero = model._zero_to_astroid * hc_vs_T(T)\n", " attempt_freq = model.attempt_freq\n", " m_therm = msat_vs_T(T) * volume_nuc\n", " \n", " # Poisson rate\n", " f = attempt_freq * np.exp(-m_therm * threshold_zero / (kB * T))\n", " P_flip_zero = np.exp(-f * delta_t)\n", " \n", " return 1 - P_flip_zero\n", "\n", "def P_multiflip(T, delta_t, volume_nuc):\n", " threshold_zero = model._zero_to_astroid * hc_vs_T(T)\n", " attempt_freq = model.attempt_freq\n", " m_therm = msat_vs_T(T) * volume_nuc\n", " \n", " # Poisson rate\n", " f = attempt_freq * np.exp(-m_therm * threshold_zero / (kB * T))\n", " P_flip_zero = np.exp(-f * delta_t)\n", " P_flip_once = f * delta_t * np.exp(-f * delta_t)\n", " P_multiflip = 1 - P_flip_zero - P_flip_once\n", " \n", " return P_multiflip" ] }, { "cell_type": "markdown", "id": "f2863b57", "metadata": {}, "source": [ "### The challenge of choosing delta_t\n", "\n", "Below is a visualization of the probability of flipping once, P_flip, and more than once, P_multiflip, as a function of temperature and simulation time interval.\n", "The example is for a magnet of size $220\\times80\\times25$ nm, assuming a 5% nucleation volume.\n", "\n", "At high temperatures (~800 K), the maximum simulation time interval that would keep P_multiflip below an acceptable value is about 10e-8 seconds.\n", "We see that continuous simulation from high temperatures to room temperature is difficult to do without changing the simulation time interval.\n", "A naive approach could be to stay at a low delta_t throughout all temperatures, but this would hardly simulate any physical time at all (each step accounting for 10e-8 seconds).\n", "Instead we can increase the time intervals as we decrease temperature, as long as we stay below an accepted value of P_multiflip.\n", "We call this process the \"annealing profile\" and in [](annealing-profile) we will have a closer look at how we go about choosing the right one for this experiment." ] }, { "cell_type": "code", "execution_count": null, "id": "93718260", "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# === First plot: heatmap of P_flip(T, delta_t) ===\n", "plt.figure(dpi=100)\n", "plt.subplot(121)\n", "# Range over which to plot P_flip and P_multiflip\n", "delta_t_range = np.geomspace(1e-10,1e10,500)\n", "T_range = np.linspace(500,820,500)\n", "volume = 220e-9 * 80e-9 * 25e-9\n", "\n", "# Calculate the probability heatmap of T and delta_t\n", "Tv, dtv = np.meshgrid(T_range, delta_t_range)\n", "T_v_dt_flip = P_flip(Tv, dtv, volume_nuc=volume*0.05)\n", "T_v_dt_multiflip = P_multiflip(Tv, dtv, volume_nuc=volume*0.05)\n", "\n", "# Show the probablity heatmap\n", "plt.pcolormesh(T_range, np.log10(delta_t_range), T_v_dt_flip, shading='auto', rasterized=True)\n", "plt.ylabel(r'$\\log(\\Delta t~[s])$')\n", "plt.xlabel(r'$T [K]$')\n", "plt.title('P_flip')\n", "\n", "# === Second plot: heatmap of P_flip(T, delta_t) ===\n", "plt.subplot(122)\n", "\n", "# Range over which to plot P_flip and P_multiflip\n", "delta_t_range = np.geomspace(1e-10,1e10,500)\n", "T_range = np.linspace(500,820,500)\n", "\n", "# Calculate the probability heatmap of T and delta_t\n", "Tv, dtv = np.meshgrid(T_range, delta_t_range)\n", "T_v_dt_multiflip = P_multiflip(Tv, dtv, volume_nuc=volume*0.05)\n", "\n", "# Show the probablity heatmap\n", "plt.pcolormesh(T_range, np.log10(delta_t_range), T_v_dt_multiflip, shading='auto', rasterized=True)\n", "plt.xlabel(r'$T [K]$')\n", "plt.title('P_multiflip')\n", "plt.colorbar();" ] }, { "cell_type": "markdown", "id": "8357ed55", "metadata": {}, "source": [ "## 3. Nucleation moment, $M_{th}$\n", "\n", "Now that we have the probabilities of flipping we can tune the final parameter, the nucleation moment, to fit with an expected (lack of) activity as described in {cite}`Zhang2019`.\n", "We expect the magnets to not flip at all below $T=673$ K, and we use this to define the nucleation volume (and thus the nucleation moment, $M_{th}$, to be used by flatspin). Thus we set `P_desired_flip = 0.00001` at timestep length `dt_freeze = 60` and temperature `T_freeze = 673`." ] }, { "cell_type": "code", "execution_count": null, "id": "e2b10024", "metadata": {}, "outputs": [], "source": [ "T_freeze = 673 # kelvin\n", "dt_freeze = 60 # If nothing happens in 60 seconds, we say the magnets are frozen\n", "P_desired_flip = 0.00001 # At freeze out, we want flips to be unlikely\n", "\n", "# ==== Find nucleation volume ====\n", "# First, create a f(nuc_volume) that we can evaluate against P_desired_flip\n", "volume = 220e-9*80e-9*25e-9 # The actual volume of the magnets\n", "\n", "def P_flip_from_V_nuc(V_nuc):\n", " return P_flip(T_freeze, dt_freeze, volume_nuc=V_nuc)\n", "\n", "# Plot P(nuc_volume) over the relevant nuc_volume range\n", "plt.figure(dpi=100)\n", "# We assume here that the interesting part will be for volumes less than half the total volume\n", "nuc_fractions = np.linspace(0.01,0.5,1000)\n", "plt.plot(nuc_fractions, P_flip_from_V_nuc(nuc_fractions*volume), label='P_flip = f(nuc_volume)')\n", "plt.xlabel('Nucleation volume fraction')\n", "plt.ylabel('P_flip')\n", "\n", "# Plot the target, P_desired_flip\n", "plt.plot(nuc_fractions, [P_desired_flip]*len(nuc_fractions), label='P_flip, target')\n", "\n", "# ==== Do bisection of P_flip_from_V_nuc to find the desired value ====\n", "# First create the function to search for roots\n", "def P_flip_min_desired(nuc_fraction):\n", " return (P_flip_from_V_nuc(nuc_fraction*volume)-P_desired_flip)\n", "plt.plot(nuc_fractions, P_flip_min_desired(nuc_fractions), label='Root searching function')\n", "\n", "# Search for roots\n", "nuc_vol_fraction = optimize.bisect(P_flip_min_desired, 0, 1, disp=True)\n", "volume_nuc = nuc_vol_fraction*volume\n", "\n", "# Print results\n", "print(f'Found nuc volume:\\n {volume_nuc} out of {volume},\\nnuc_fraction = {nuc_vol_fraction}')\n", "plt.scatter(nuc_vol_fraction, P_flip_from_V_nuc(volume_nuc), c='r', label='Identified nucleation volume')\n", "plt.legend()\n", "plt.ylim((-2*P_desired_flip,P_desired_flip*50))\n", "plt.title('T = 673 K, delta_t = 60 s');" ] }, { "cell_type": "markdown", "id": "0c15f698", "metadata": {}, "source": [ "(annealing-profile)=\n", "## 4. Annealing profile\n", "\n", "We are now ready to look at the effects of changing the timestep length and temperature, and use this to find a suitable annealing protocol.\n", "\n", "Ideally, we would do just like in the experiment {cite}`Zhang2013`, decreasing by 1 minute per kelvin, starting at `T_max = 800` kelvin, simulated as a single step per temperature value.T\n", "If `P_multiflip` rises above the accepted `multiflip_threshold`, we must reduce the timestep length and do more, shorter steps at these temperatures.\n", "\n", "We will see that this approach is a bit too naive." ] }, { "cell_type": "code", "execution_count": null, "id": "3adbe253", "metadata": {}, "outputs": [], "source": [ "# === The \"ideal\" annealing profile ===\n", "multiflip_threshold = 0.001\n", "dt_cap = 60\n", "\n", "# === Plot: heatmap of P_flip(T, delta_t) ===\n", "plt.figure(dpi=100)\n", "\n", "# Range over which to plot P_flip and P_multiflip\n", "delta_t_range = np.geomspace(1e-10,1e10,500)\n", "T_range = np.linspace(500,820,500)\n", "\n", "# Calculate the probability heatmap of T and delta_t\n", "Tv, dtv = np.meshgrid(T_range, delta_t_range)\n", "T_v_dt_multiflip = P_multiflip(Tv, dtv, volume_nuc=volume_nuc)\n", "\n", "# Show the probablity heatmap\n", "plt.pcolormesh(T_range, np.log10(delta_t_range), T_v_dt_multiflip, shading='auto', rasterized=True)\n", "plt.ylabel(r'$\\log(\\Delta t)$')\n", "plt.xlabel(r'$T [K]$') \n", "plt.title('P_multiflip')\n", "plt.colorbar()\n", "\n", "# Plot \"freeze out cursor\", i.e., where we expect the magnets to be stable\n", "plt.axvline(T_freeze, c='c', label=f'T = {T_freeze} K (blocking temp)')\n", "plt.axhline(np.log10(dt_cap), c='C1', label=f'delta_t = {dt_cap} s')\n", "\n", "# Find T and delta_t where P_multiflip == multiflip_threshold\n", "# For now, we create a delta_t(T) function from this contour object,\n", "# in the next section we will derive it properly.\n", "cs = plt.contour(T_range, np.log10(delta_t_range), T_v_dt_multiflip, [multiflip_threshold], colors='r')\n", "cs.collections[0].set_label(f'P_multiflip = {multiflip_threshold}')\n", "p = cs.collections[0].get_paths()[0]\n", "v = p.vertices\n", "\n", "# Plot the \"ideal\" temperature profile\n", "T_ideal = np.arange(800, T_freeze, -1)\n", "dt_ideal = np.minimum(np.repeat(np.log10(dt_freeze), len(T_ideal)), np.interp(T_ideal, v[:,0], v[:,1]))\n", "plt.plot(T_ideal, dt_ideal, c='C2', label = '\"Ideal\" temperature profile')\n", "plt.scatter(*[(min(T_ideal), max(T_ideal)),(max(dt_ideal), min(dt_ideal))], c='C2')\n", "\n", "plt.legend(loc='upper left', bbox_to_anchor=(1.2, 1.0));" ] }, { "cell_type": "markdown", "id": "dffba6c5", "metadata": {}, "source": [ "Implementing the above temperature profile, with the goal of simulating one minute at each degree kelvin would require about 10^10 steps, just for the $T = 800$ K step. \n", "To make our simulation feasible, we abandon the hope of simulating a full minute at all temperatures, and settle on a fixed number of simulation step per temperature. \n", "For the lower temperatures, we end up simulating at the desired 1 minute per kelvin, but at higher temperatures this compromise simulates less actual time.\n", "However, with enough simulation steps and high switching frequencies we assume that we do not miss important dynamics. \n", "\n", "To create this compromise temperature profile, it will be useful to calculate the maximum timestep length we can do without violating a certain `P_flip` or `P_multiflip` requirement. \n", "Below we create a function to yield the maximum timestep length and calculate a temperature profile." ] }, { "cell_type": "code", "execution_count": null, "id": "3c3c18c4", "metadata": {}, "outputs": [], "source": [ "from scipy.special import lambertw\n", "from numpy import real\n", "\n", "# === First we create a way to calculate the maximum delta_t ===\n", "def calc_delta_t(T, m_therm, delta_h, P_flip=None, P_multiflip=None, f0=1e9):\n", " \"\"\" Calculate maximum delta_t for target P_flip or P_multiflip.\n", " Returns the delta_t that would result in a probability that is\n", " smaller than (or equal to) the one supplied.\n", " \"\"\"\n", " f = f0 * np.exp(-m_therm * delta_h / (kB * T))\n", " if P_flip is not None:\n", " assert P_multiflip is None\n", " # Using the equation for of flipping at least once,\n", " # 1 - np.exp(-f * delta_t),\n", " # and solving for delta_t is straight forward\n", " return np.log(1 - P_flip)/-f\n", " \n", " elif P_multiflip is not None:\n", " assert P_flip is None\n", " # Using the expression for the likelihood of multiple switching events,\n", " # 1 - np.exp(-f * delta_t) - f * delta_t * np.exp(-f * delta_t),\n", " # to solve for delta_t is not trivial.\n", " # We do so here by using the Lambert W function from scipy.\n", " return real((1/f)*(-lambertw((P_multiflip - 1)/np.exp(1), -1) - 1))\n", " \n", " raise ValueError(\"Need either P_flip or P_multiflip\")\n", " \n", "def max_delta_t(T, P_flip=None, P_multiflip=None):\n", " delta_h = model._zero_to_astroid * hc_vs_T(T)\n", " f0 = model.attempt_freq\n", " msat = msat_vs_T(T)\n", " m_therm = msat * volume_nuc\n", " return calc_delta_t(T, m_therm, delta_h, P_flip, P_multiflip, f0)\n", "\n", "# === Return the temperatures, delta_ts, number of timesteps, \n", "# and temperature where 1 min/K is reached (T_breakoff)\n", "def make_fixed_temperature_profile(T_max, T_min, timesteps_per_K, multiflip_threshold, max_tpk=60):\n", " Ts = np.arange(T_max, T_min-.1, -1)\n", " max_dt = max_tpk / timesteps_per_K\n", " dts = max_delta_t(Ts, P_multiflip=multiflip_threshold)\n", " T_breakoff = Ts[np.argmax(dts > max_dt)]\n", " dts[dts>max_dt] = max_dt\n", " ns = np.repeat(timesteps_per_K, len(Ts))\n", "\n", " return Ts, dts, ns, T_breakoff\n", " \n", "T_max = 800\n", "timesteps_per_K = 50\n", "Ts, dts, ns, T_breakoff = make_fixed_temperature_profile(T_max, T_freeze, timesteps_per_K, multiflip_threshold)\n", "time_per_kelvin = ns * dts\n", "\n", "print(f\"Temperature range to be simulated at less than {np.max(time_per_kelvin)} seconds: {T_max} - {T_breakoff}\")\n", "print(f\"Time per kelvin range: {np.min(time_per_kelvin)} - {np.max(time_per_kelvin)}\")\n", "print(f\"Total timesteps: {sum(ns)}\")\n", "# print(\"Ts =\", Ts)\n", "# print(\"ns =\", ns)\n", "# print(\"dts =\", dts)\n", "# print(\"ns * dts=\", time_per_kelvin)" ] }, { "cell_type": "code", "execution_count": null, "id": "ea8fbfe4", "metadata": {}, "outputs": [], "source": [ "# === Plot everything ===\n", "# Plot the temperature profile in the P_flip heatmap\n", "plt.figure(dpi=100)\n", "plt.subplot(1,1,1)\n", "T_v_dt = P_flip(Tv, dtv, volume_nuc=volume_nuc)\n", "plt.pcolormesh(T_range, np.log10(delta_t_range), T_v_dt, shading='auto', rasterized=True)\n", "\n", "# T_freeze line\n", "plt.axvline(T_freeze, label='T_freeze', color='c')\n", "\n", "#for T_max in np.arange(740, 821, 10):\n", "Ts, dts, ns, T_breakoff = make_fixed_temperature_profile(T_max, T_freeze, timesteps_per_K, multiflip_threshold)\n", "tpk = ns * dts\n", "plt.plot(Ts, np.log10(dts), c='C2', label='Timestep length')\n", "plt.scatter([Ts[0], Ts[-1]], [np.log10(dts[0]), np.log10(dts[-1])], c='C2')\n", "\n", "plt.ylabel(r'$\\log10(\\Delta t)$')\n", "plt.xlabel(r'$T$ [K]') \n", "plt.title('P_flip')\n", "plt.colorbar()\n", "plt.plot(Ts, np.log10(ns*dts), c='C1', label='Total time per $T$')\n", "plt.axhline(np.log10(60), c='C1', ls=':', label=f'{dt_freeze} seconds')\n", "plt.legend(loc=(0.04,0.04));" ] }, { "cell_type": "markdown", "id": "ad93e6d7", "metadata": {}, "source": [ "In the above plot, we show the timestep length of the annealing profile (green). \n", "Starting at $T = 800$ K, and cooling down, the timestep length is gradually increased until the total simulated time at each temperature reaches 60 seconds. The timestep length is then kept constant (as the total simulated time per temperature stays at 60 seconds), until `T_freeze`." ] }, { "cell_type": "markdown", "id": "e15fe10f", "metadata": {}, "source": [ "## 5. flatspin run command and input files\n", "\n", "Now that we have our parameter arrays, we can generate a command to run the experiment.\n", "We will dump the parameters to CSV files and create the command in accordance with the other parameters outlined in the experiment {cite}`Zhang2013`." ] }, { "cell_type": "code", "execution_count": null, "id": "d5a95872", "metadata": {}, "outputs": [], "source": [ "import os\n", "from flatspin.data import save_table\n", "\n", "# === Save params.csv ===\n", "# Output files \n", "root_dir = \"params_folder\"\n", "params_path = root_dir + \"/params.csv\"\n", "if not os.path.exists(root_dir):\n", " os.makedirs(root_dir)\n", " print(f\"Creating directory {root_dir}\")\n", "\n", "# === Define the relevant model parameters (not temperature dependent) ===\n", "# Commented out parameters are normally supplied, but are here replaced by sweeping parameters.\n", "params = dict(\n", " size=(50, 50),\n", "# alpha=alpha_calc(msat = msat_vs_T(300), V = volume, a = 1e-9),\n", " neighbor_distance=10,\n", "# hc=hc0,\n", " disorder=0.05,\n", " sw_b=0.38, sw_c=1.0, sw_beta=1.3, sw_gamma=3.6,\n", "# m_therm=msat * volume_nuc,\n", " use_opencl=1,\n", " H=0, spp=1, timesteps=1,\n", ")\n", "\n", "print(f\"Saving static params to {params_path}\")\n", "save_table(params, params_path)" ] }, { "cell_type": "code", "execution_count": null, "id": "2a578d93", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "# === Save temperature dependent parameters and genereate command with all lattice spacings ===\n", "lattice_spacings = [320, 360, 400, 440, 480, 560, 680, 880]\n", "\n", "# Temperature dependent parameters\n", "# All are technically a function of time, hence the '_t' subscript\n", "Ts_out = np.repeat(Ts, ns)\n", "dts_out = np.repeat(dts, ns)\n", "msats = msat_vs_T(Ts_out)\n", "alphas = alpha_calc(msats, volume, 1e-9)\n", "m_therms = msats*volume_nuc\n", "hcs_out = hc_vs_T(Ts_out)\n", "\n", "csvs = [\n", " ('temperature_t', Ts_out),\n", " ('therm_timescale_t', dts_out),\n", " ('alpha_t', alphas),\n", " ('m_therm_t', m_therms),\n", " ('hc_t', hcs_out),\n", "]\n", "\n", "index = dict()\n", "for k, v in csvs:\n", " filename = f\"{k}.csv\"\n", " filepath = root_dir + \"/\" + filename\n", " print(f\"Creating {filepath}\")\n", " np.savetxt(filepath, v)\n", " index[k] = filepath\n", "\n", "# Parameters it is nice to keep track off, but which are not part of the model are also saved in the index\n", "index[\"T_max\"] = T_max\n", "index[\"T_breakoff\"] = T_breakoff\n", "index[\"T_min\"] = T_freeze\n", "index[\"multiflip_threshold\"] = multiflip_threshold\n", "index[\"periods\"] = len(Ts_out)\n", "\n", "index = pd.DataFrame([index])\n", "index = index.reindex(index.index.repeat(len(lattice_spacings))).reset_index(drop=True)\n", "index.insert(loc=0, column='lattice_spacing', value=lattice_spacings)\n", "\n", "print(f\"Generated {len(index)} runs\")\n", "display(index)\n", "\n", "index_file = \"index.csv\"\n", "index_path = root_dir + \"/\" + index_file\n", "print(f\"Saving index to {index_path}\")\n", "save_table(index, index_path)\n", "\n", "# Generate flatspin-run-sweep command\n", "print(\"\\nGENERATED COMMAND: \\n\")\n", "print(f\"flatspin-run-sweep -r dist -m SquareSpinIceClosed -e Sine \", end='')\n", "print(f\"--params-file={params_path} --sweep-file={index_path} \", end='')\n", "print(\"-o CHANGE_TO_OUTPUT_FOLDER\") # The folder to create the " ] }, { "cell_type": "markdown", "id": "a94205ab", "metadata": {}, "source": [ "**We are now ready to run the above command!**\n", "\n", "Remember to upload the `params-file` (`params.csv`) and the `sweep-file` (`index.csv`) to the file location you run flatspin from. \n", "\n", "Notes: \n", "\n", " - The `-r dist` flag tells flatspin to run distributed on a cluster. Skip this if you just run this locally. \n", " - The `-e Sine` flag signifies the use of a sine encoder for the external field, but there is no input in this experiment. This flag is there to satisfy the requirement of providing an encoder to flatspin runs, but will have noe effect since there's no input and the field is set to zero with the `H=0` parameter." ] }, { "cell_type": "markdown", "id": "8c657053", "metadata": {}, "source": [ "### 6. Results\n", "\n", "We have our results. Now, let's have a look! \n", "Below we define a function to read out the spin state of any timestep from our resulting files. We can use this to animate the timeseries of the annealing process." ] }, { "cell_type": "code", "execution_count": null, "id": "58c56c07", "metadata": {}, "outputs": [], "source": [ "from flatspin.data import Dataset, read_table, read_geometry, read_vectors, vector_grid\n", "from flatspin.grid import Grid\n", "from flatspin.vertices import find_vertices, vertex_pos, vertex_mag\n", "from flatspin.plotting import plot_vectors, plot_vector_image\n", "\n", "\n", "def read_data_vertices_grid(dataset, t=[0]):\n", " ''' Return the positions and magnetization of the dataset's vertices '''\n", " grid_size = None \n", " crop_width = ((1,1),(1,1))\n", " win_shape, win_step = ((3,3), (2, 2))\n", " \n", " df = read_table(dataset.tablefile('mag'), index_col='t')\n", " UV = np.array(df)\n", " UV = UV.reshape((UV.shape[0], -1, 2))\n", " UV = UV[t]\n", " \n", " pos, angle = read_geometry(dataset.tablefile('geometry'))\n", "\n", " XY, UV = vector_grid(pos, UV, grid_size, crop_width, win_shape, win_step, normalize=False)\n", " return XY, UV/2" ] }, { "cell_type": "markdown", "id": "c7156a70", "metadata": {}, "source": [ "#### Read the datasets\n", "\n", "This might take some time." ] }, { "cell_type": "code", "execution_count": null, "id": "97318c68", "metadata": {}, "outputs": [], "source": [ "# Read the dataset\n", "dataset = Dataset.read('/data/flatspin/annealing-example-run')\n", "dataset = dataset.filter(multiflip_threshold=multiflip_threshold, periods=len(Ts_out))\n", "\n", "# Read the states of the vertices at specified timesteps\n", "# (read at the end of each temperature step)\n", "periods = len(Ts_out)\n", "ts = np.arange(timesteps_per_K -1, periods, timesteps_per_K)\n", "vertex_states = dict()\n", "for ls, ds in dataset.groupby('lattice_spacing'):\n", " vertex_states[ls] = read_data_vertices_grid(ds, t=ts)" ] }, { "cell_type": "markdown", "id": "911a82f5", "metadata": {}, "source": [ "#### Plot the final frame" ] }, { "cell_type": "code", "execution_count": null, "id": "b39e9a31", "metadata": {}, "outputs": [], "source": [ "def plot_vertices(XY, UVi, ax):\n", " ax.axis('off')\n", " plot_vector_image(XY, UVi, ax=ax, replace=True)\n", "\n", "fig, axs = plt.subplots(2, 4, figsize=(12,6), dpi=200)\n", "for ax, (ls, vstate) in zip(axs.flat, vertex_states.items()):\n", " ax.set_title(f'lattice spacing = {ls} nm')\n", " XY, UVi = vstate\n", " mag = np.amax(np.linalg.norm(UVi[-1], axis=-1))\n", " plot_vertices(XY, UVi[-1]/mag, ax)\n", "fig.suptitle('Final timestep', fontsize='x-large');" ] }, { "cell_type": "markdown", "id": "eb1bd01f", "metadata": {}, "source": [ "#### Animate over all temperatures" ] }, { "cell_type": "code", "execution_count": null, "id": "5aabf4ea", "metadata": {}, "outputs": [], "source": [ "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML\n", "\n", "num_frames = periods // timesteps_per_K\n", "temperatures = list(read_table(dataset.tablefiles('params_t')[0][0])['temperature'][timesteps_per_K-1::timesteps_per_K])\n", "def animate_annealing(vstates):\n", " # Set up figure and axes\n", " fig, axs = plt.subplots(2, 4, figsize=(12,6), dpi=200)\n", " for ax, (ls, vstate) in zip(axs.flat, vstates.items()):\n", " ax.set_title(f'lattice spacing = {ls} nm')\n", " \n", " def do_plot(t):\n", " fig.suptitle(f'T = {temperatures[t]} K', y=0.0, va='bottom', fontsize='xx-large')\n", " for ax, (ls, vstate) in zip(axs.flat, vstates.items()):\n", "# ax.set_title(f'lattice spacing = {ls} nm')\n", " XY, UVi = vstate\n", " mag = np.amax(np.linalg.norm(UVi[-1], axis=-1))\n", " plot_vertices(XY, UVi[t]/mag, ax)\n", " \n", " anim = FuncAnimation(fig, do_plot, frames=num_frames, interval=150, blit=False)\n", " plt.close() # Only show the animation\n", " #anim.save(\"astroid.gif\")\n", " return HTML(anim.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "id": "71a9dc25", "metadata": {}, "outputs": [], "source": [ "animate_annealing(vertex_states)" ] }, { "cell_type": "markdown", "id": "31418a50", "metadata": {}, "source": [ "## References\n", "\n", "```{bibliography}\n", "```" ] } ], "metadata": { "jupytext": { "formats": "ipynb,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.8.10" }, "source_map": [ 15, 19, 44, 54, 90, 100, 111, 115, 134, 147, 176, 189, 227, 234, 273, 285, 327, 337, 398, 422, 427, 434, 464, 518, 529, 536, 558, 564, 576, 580, 592, 596, 622, 624 ] }, "nbformat": 4, "nbformat_minor": 5 }