# Dynamics#

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\))

## Switching energy#

The **switching energy** captures both these conditions, and can be calculated with `switching_energy()`

.
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 `flippable()`

.

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]
```

## Flipping spins#

flatspin employs single spin flip dynamics, where only one spin is flipped at a time.
Calling `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, cache_frame_data=False)
plt.close() # Only show the animation
HTML(anim.to_jshtml())
```

## Relaxation#

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()`

.
`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, cache_frame_data=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();
```

## Hysteresis loops#

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 `use_opencl=True`

.

```
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, cache_frame_data=False)
plt.close() # Only show the animation
HTML(anim.to_jshtml())
```