Tutorial 08: Standard problem 4

Interactive online tutorial: Binder

Problem specification

The sample is a thin film cuboid with dimensions:

  • length $l_{x} = 500 \,\text{nm}$,
  • width $l_{y} = 125 \,\text{nm}$, and
  • thickness $l_{z} = 3 \,\text{nm}$.

The material parameters (similar to permalloy) are:

  • exchange energy constant $A = 1.3 \times 10^{-11} \,\text{J/m}$,
  • magnetisation saturation $M_\text{s} = 8 \times 10^{5} \,\text{A/m}$.

Magnetisation dynamics are governed by the Landau-Lifshitz-Gilbert equation

$$\frac{d\mathbf{m}}{dt} = \underbrace{-\gamma_{0}(\mathbf{m} \times \mathbf{H}_\text{eff})}_\text{precession} + \underbrace{\alpha\left(\mathbf{m} \times \frac{d\mathbf{m}}{dt}\right)}_\text{damping}$$

where $\gamma_{0} = 2.211 \times 10^{5} \,\text{m}\,\text{A}^{-1}\,\text{s}^{-1}$ and Gilbert damping $\alpha=0.02$.

In the standard problem 4, the system is first relaxed at zero external magnetic field and then, starting from the obtained equlibrium configuration, the magnetisation dynamics are simulated for two external magnetic fields $\mathbf{B}_{1} = (-24.6, 4.3, 0.0) \,\text{mT}$ and $\mathbf{B}_{2} = (-35.5, -6.3, 0.0) \,\text{mT}$.

More detailed specification of Standard problem 4 can be found in Ref. 1.

Simulation

In the first step, we import the required discretisedfield and oommfc modules.

In [1]:
import discretisedfield as df
import micromagneticmodel as mm
import oommfc as oc

Now, we can set all required geometry and material parameters.

In [2]:
# Geometry
lx = 500e-9  # x dimension of the sample(m)
ly = 125e-9  # y dimension of the sample (m)
lz = 3e-9  # sample thickness (m)

# Material (permalloy) parameters
Ms = 8e5  # saturation magnetisation (A/m)
A = 1.3e-11  # exchange energy constant (J/m)

# Dynamics (LLG equation) parameters
gamma0 = 2.211e5  # gyromagnetic ratio (m/As)
alpha = 0.02  # Gilbert damping

First stage

In the first stage, we need to relax the system at zero external magnetic field.

We choose stdprob4 to be the name of the system. This name will be used to name all output files created by OOMMF.

In [3]:
system = mm.System(name='stdprob4')

In order to completely define the micromagnetic system, we need to provide:

  1. energy $E$
  2. dynamics $\text{d}\mathbf{m}/\text{d}t$
  3. magnetisation $\mathbf{m}$

The mesh is created by providing two points p1 and p2 between which the mesh domain spans and the size of a discretisation cell. We choose the discretisation to be $(5, 5, 3) \,\text{nm}$.

In [4]:
cell = (5e-9, 5e-9, 3e-9)  # mesh discretisation (m)
mesh = df.Mesh(p1=(0, 0, 0), p2=(lx, ly, lz), cell=cell)  # Create a mesh object.

We can visualise the mesh domain and a discretisation cell:

In [5]:
%matplotlib inline
mesh.k3d()

Hamiltonian: In the second step, we define the system's Hamiltonian. In this standard problem, the Hamiltonian contains only exchange and demagnetisation energy terms. Please note that in the first simulation stage, there is no applied external magnetic field. Therefore, we do not add Zeeman energy term to the Hamiltonian.

In [6]:
system.energy = mm.Exchange(A=A) + mm.Demag()

We can check what is the continuous model of system's Hamiltonian.

In [7]:
system.energy
Out[7]:
$A \left[ (\nabla \mathbf{m}_\text{x})^{2} + (\nabla \mathbf{m}_\text{y})^{2} + (\nabla \mathbf{m}_\text{z})^{2} \right]-\frac{1}{2}\mu_{0}M_\text{s}\mathbf{m} \cdot \mathbf{H}_\text{d}$

Dynamics: Similarly, the system's dynamics is defined by providing precession and damping terms (LLG equation).

In [8]:
system.dynamics = mm.Precession(gamma0=gamma0) + mm.Damping(alpha=alpha)

system.dynamics  # check the dynamics equation
Out[8]:
$-\frac{\gamma_{0}}{1 + \alpha^{2}} \mathbf{m} \times \mathbf{H}_\text{eff}-\frac{\gamma_{0} \alpha}{1 + \alpha^{2}} \mathbf{m} \times (\mathbf{m} \times \mathbf{H}_\text{eff})$

Magnetisation: Finally, we have to provide the magnetisation configuration that is going to be relaxed subsequently. We choose the uniform configuration in $(1, 0.25, 0.1)$ direction, and as norm (magnitude) we set the magnetisation saturation $M_\text{s}$. In order to create the magnetisation configuration, we create a Field object from the discretisedfield module.

In [9]:
system.m = df.Field(mesh, dim=3, value=(1, 0.25, 0.1), norm=Ms)

Now, the system is fully defined.

Energy minimisation: The system (its magnetisation) is evolved using a particular driver. In the first stage, we need to relax the system - minimise its energy. Therefore, we create MinDriver object and drive the system using its drive method.

In [10]:
md = oc.MinDriver()  # create energy minimisation driver
md.drive(system)  # minimise the system's energy
Running OOMMF (ExeOOMMFRunner) [2020/06/14 11:27]... (3.2 s)

The system is now relaxed. We can plot it in various ways. We start with simple plot using default settings:

In [11]:
# x-component
system.m.x.plane('z').mpl()
In [12]:
# y-component
system.m.y.plane('z').mpl()
In [13]:
# vectors
system.m.plane('z').mpl()

With a little bit of additional code, we can fine tune the plot as appropriate for the relevant results:

In [14]:
# vectors
import matplotlib.pyplot as plt

# make figure larger
fig, ax = plt.subplots(figsize=(20, 5))

# plot vectors on grid of 20 x 5 over the numerical resulotion
system.m.plane('z', n=(20, 5)).mpl_vector(ax=ax)

# add colouring for mx-component to this plot
system.m.x.plane('z').mpl_scalar(ax=ax)

We can now obtain some numerical data characteristic to the magnetisation field as is useful for standard problem 4:

In [15]:
print('The average magnetisation is {}.'.format(system.m.average))

print('The magnetisation at the mesh centre {}\nis {}.'.format(
        system.m.mesh.region.centre, system.m(system.m.mesh.region.centre)))
The average magnetisation is (773766.2112987704, 99856.85901465818, -0.003159960962047351).
The magnetisation at the mesh centre (2.5e-07, 6.25e-08, 1.5e-09)
is (799979.276019471, -5758.29309460144, -0.00502709656498145).

Second stage: field $\mathbf{B}_{1}$

In the second stage, we need to apply an external magnetic field $\mathbf{B}_{1} = (-24.6, 4.3, 0.0) \,\text{mT}$ to the system. In other words, we have to add Zeeman energy term to the Hamiltonian.

In [16]:
# Add Zeeman energy term to the Hamiltonian
H1 = (-24.6e-3/mm.consts.mu0, 4.3e-3/mm.consts.mu0, 0.0)
system.energy += mm.Zeeman(H=H1)

If we now inspect the Hamiltonian, we see that an additional Zeeman term is added.

In [17]:
system.energy
Out[17]:
$A \left[ (\nabla \mathbf{m}_\text{x})^{2} + (\nabla \mathbf{m}_\text{y})^{2} + (\nabla \mathbf{m}_\text{z})^{2} \right]-\frac{1}{2}\mu_{0}M_\text{s}\mathbf{m} \cdot \mathbf{H}_\text{d}-\mu_{0}M_\text{s} \mathbf{m} \cdot \mathbf{H}$

Finally, we can run the simulation using TimeDriver this time. We run the magnetisation evolution for $t=1 \,\text{ns}$, during which we save the system's state $n=200$ times.

In [18]:
t = 1e-9  # simulation time (s)
n = 200  # number of data saving steps

td = oc.TimeDriver()  # create time driver
td.drive(system, t=t, n=n)  # drive the system
Running OOMMF (ExeOOMMFRunner) [2020/06/14 11:27]... (7.1 s)

Postprocessing

When we drove the system using the TimeDriver, we specified that we want to save the magnetisation configuration $n=200$ times. A detailed table of all computed parameters from the last simulation can be shown from the datatable (system.dt), which is a pandas dataframe [2].

For instance, if we want to show the last 10 rows in the table, we run:

In [19]:
system.table.data.tail()
Out[19]:
E E_calc_count max_dm/dt dE/dt delta_E E_exchange max_spin_ang stage_max_spin_ang run_max_spin_ang E_demag E_zeeman iteration stage_iteration stage mx my mz last_time_step t
195 -2.676342e-18 3736.0 1162.539933 -1.697359e-09 -2.278389e-21 9.548545e-20 3.913558 4.281302 29.612838 8.530166e-19 -3.624844e-18 784.0 3.0 195.0 -0.984260 -0.010971 0.033987 1.374339e-12 9.800000e-10
196 -2.685463e-18 3755.0 1260.312762 -1.935225e-09 -2.621610e-21 8.637446e-20 4.157989 4.157989 29.612838 8.844369e-19 -3.656275e-18 788.0 3.0 196.0 -0.987086 0.021592 0.039286 1.373845e-12 9.850000e-10
197 -2.695498e-18 3774.0 1384.782389 -2.057789e-09 -2.813912e-21 8.039753e-20 4.727728 4.727728 29.612838 9.074607e-19 -3.683357e-18 792.0 3.0 197.0 -0.988092 0.057824 0.042604 1.373937e-12 9.900000e-10
198 -2.705827e-18 3793.0 1351.611377 -2.053171e-09 -2.832075e-21 7.586011e-20 4.669401 4.772024 29.612838 9.220444e-19 -3.703732e-18 796.0 3.0 198.0 -0.986964 0.095870 0.043794 1.374355e-12 9.950000e-10
199 -2.715833e-18 3812.0 1191.086615 -1.931645e-09 -2.687310e-21 7.141300e-20 4.354488 4.669401 29.612838 9.291427e-19 -3.716389e-18 800.0 3.0 199.0 -0.983765 0.133793 0.042832 1.375196e-12 1.000000e-09

Finally, we want to plot the average magnetisation configuration my as a function of time t:

In [21]:
myplot = system.table.data.plot('t', 'my')

We can plot another snap shot of the magnetisation at the end of the run, which shows the fluctuations in the simulation at this time:

In [22]:
# vectors
import matplotlib.pyplot as plt

# make figure larger
fig, ax = plt.subplots(figsize=(20, 5))

# plot vectors on grid of 20 x 5 over the numerical resulotion
system.m.plane('z', n=(20, 5)).mpl_vector(ax=ax)

# add colouring for mx-component to this plot
system.m.x.plane('z').mpl_scalar(ax=ax)

References

[1] ┬ÁMAG Site Directory: http://www.ctcms.nist.gov/~rdm/mumag.org.html

[2] Pandas: http://pandas.pydata.org/