As discussed in Switching, flatspin employs a generalized Stoner-Wohlfarth (GSW) switching model to determine whether a spin will flip its magnetization. A spin may flip if the total magnetic field acting on the spin:
is outside of the switching astroid (left hand side of the GSW equation is greater than 1) and
is oriented in the opposite direction of the spin magnetization (\(h_\parallel < 0\))
The switching energy captures both these conditions, and can be calculated with
When the switching energy is positive, the total field is outside the astroid (1) and antiparallel to the spin magnetization (2).
The magnitude of the switching energy corresponds to how far the total field is outside the astroid.
A spin may flip if, and only if, the switching energy is positive.
A list of flippable spins can be obtained by calling
Below we color the spins in a square spin ice by their corresponding switching energies:
from flatspin.model import SquareSpinIceClosed model = SquareSpinIceClosed() model.spin[[0, 4, 39]] = -1 # flip spins 0, 4 and 39 print("Flippable spins:", model.flippable()) # Calculate the switching energy E = model.switching_energy() E_max = np.max(np.abs(E)) quiv = model.plot(C=E, cmap='coolwarm', clim=(-E_max, E_max)) plt.colorbar(quiv);
Flippable spins: 
In the example above, the switching energy was always negative, meaning no spins were flippable. If we increase the external field, we will see that some spins get a positive switching energy and become flippable.
# Increase the external field model.set_h_ext([0.2, 0]) print("Flippable spins:", model.flippable()) # Calculate the switching energy E = model.switching_energy() E_max = np.max(np.abs(E)) quiv = model.plot(C=E, cmap='coolwarm', clim=(-E_max, E_max)) plt.colorbar(quiv);
Flippable spins: [ 0 4 5 13 34 39]
flatspin employs single spin flip dynamics, where only one spin is flipped at a time.
step() advances the simulation one step by flipping a spin.
The dipolar fields are always recalculated as part of a simulation step, but the thermal fields are not.
If there are more than one flippable spin, the spin with the highest switching energy will flip first. This ensures that the spin flip order is completely deterministic.
In the animation below we flip one spin at a time, until there are no more flippable spins. Notice how flipping a spin affects the switching energy of neighboring spins, due to the change in dipolar fields. Consequently, only 3 spins will actually flip, even though there were 6 flippable spins originally.
from matplotlib.animation import FuncAnimation from IPython.display import HTML def do_steps(): yield model while model.step(): yield model def do_plot(model): E = model.switching_energy() E_max = np.max(np.abs(E)) model.plot(C=E, cmap='coolwarm', clim=(-E_max, E_max), ax=ax, replace=True) fig, ax = plt.subplots() anim = FuncAnimation(fig, do_plot, frames=do_steps(), interval=500, blit=False) plt.close() # Only show the animation HTML(anim.to_jshtml())
Flipping all spins until there are no more flippable spins is referred to as a relaxation, and is the core dynamical process of flatspin.
A simulation is advanced by updating the external and/or thermal fields, and then flipping all flippable spins by calling
relax() will re-sample the thermal fields once, before starting to flip spins.
Because the external and thermal fields are held constant during relaxation, the process is guaranteed to terminate.
relax() returns the number of spins that flipped.
Below we use
relax() to simulate the reversal of a square spin ice by an increasing external field.
# Strength of external field H = np.linspace(0, 0.1, 100) # Direction of the external field phi = np.deg2rad(180+45) h_dir = np.array([np.cos(phi), np.sin(phi)]) def do_reversal(): yield model for h in H: model.set_h_ext(h * h_dir) if model.relax(): # Only show plot when there were spin flips yield model def do_plot(model): model.plot(ax=ax, replace=True) fig, ax = plt.subplots() anim = FuncAnimation(fig, do_plot, frames=do_reversal(), interval=500, blit=False) plt.close() # Only show the animation HTML(anim.to_jshtml())
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. 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:
# Gradually increase strength of external field in 10 steps model.polarize() H = 0.058 for h in np.linspace(0, H, 10): model.set_h_ext([-h, -h]) model.relax() plt.figure() plt.title("After 10 field steps") model.plot() # Set external field in a single step model.polarize() model.set_h_ext([-H, -H]) model.relax() plt.figure() plt.title("After 1 field step") model.plot();
Often we are interested in some average quantity of the ensemble, such as the total magnetization.
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
model = SquareSpinIceClosed(size=(10,10), use_opencl=True) # H decreases linearly from 0.1 to -0.1, then back to 0.1 H = np.linspace(0.1, -0.1, 500) H = np.concatenate([H, -H]) # Angle of the external field phi = np.deg2rad(45) h_dir = np.array([np.cos(phi), np.sin(phi)]) M_H =  for h in H: model.set_h_ext(h * h_dir) model.relax() # Magnetization projected along field direction m = model.total_magnetization().dot(h_dir) M_H.append(m) plt.plot(H, M_H) plt.xlabel("H (mT)") plt.ylabel(r"M_H (a.u.)");
Here is an animation of the hysteresis loop, where changes to the state at each
relax() are shown, omitting steps without change:
# Reset system back to the polarized state model.polarize() fig, ax = plt.subplots() def do_hysteresis(): yield model for h in H: model.set_h_ext(h * h_dir) if model.relax(): # Only show frames where any spins flipped yield model def animate(model): model.plot(ax=ax, replace=True) anim = FuncAnimation(fig, animate, frames=do_hysteresis(), interval=200, blit=False) plt.close() # Only show the animation HTML(anim.to_jshtml())