Standard problem 4#
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
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.
[1]:
import discretisedfield as df
import micromagneticmodel as mm
import oommfc as oc
Now, we can set all required geometry and material parameters.
[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.
[3]:
system = mm.System(name="stdprob4")
In order to completely define the micromagnetic system, we need to provide:
energy \(E\)
dynamics \(\text{d}\mathbf{m}/\text{d}t\)
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}\).
[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:
[5]:
mesh.mpl()
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.
[6]:
system.energy = mm.Exchange(A=A) + mm.Demag()
We can check what is the continuous model of system’s Hamiltonian.
[7]:
system.energy
[7]:
Dynamics: Similarly, the system’s dynamics is defined by providing precession and damping terms (LLG equation).
[8]:
system.dynamics = mm.Precession(gamma0=gamma0) + mm.Damping(alpha=alpha)
system.dynamics # check the dynamics equation
[8]:
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.
[9]:
system.m = df.Field(mesh, nvdim=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.
[10]:
md = oc.MinDriver() # create energy minimisation driver
md.drive(system) # minimise the system's energy
Running OOMMF (ExeOOMMFRunner)[2023/10/23 16:03]... (0.5 s)
The system is now relaxed. We can plot it in various ways. We start with simple plot using default settings:
[11]:
# x-component
system.m.x.sel("z").mpl()
[12]:
# y-component
system.m.y.sel("z").mpl()
[13]:
# vectors
system.m.sel("z").mpl()
With a little bit of additional code, we can fine tune the plot as appropriate for the relevant results:
[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.sel("z").resample((20, 5)).mpl.vector(ax=ax)
# add colouring for mx-component to this plot
system.m.x.sel("z").mpl.scalar(ax=ax, cmap="magma")
We can now obtain some numerical data characteristic to the magnetisation field as is useful for standard problem 4:
[15]:
print("The average magnetisation is {}.".format(system.m.mean()))
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 [ 7.73766211e+05 9.98568590e+04 -3.15994532e-03].
The magnetisation at the mesh centre [2.50e-07 6.25e-08 1.50e-09]
is [ 7.99979276e+05 -5.75829308e+03 -5.03057845e-03].
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.
[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.
[17]:
system.energy
[17]:
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.
[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)[2023/10/23 16:03]... (1.9 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:
[19]:
system.table.data.tail()
[19]:
E | E_calc_count | max_dm/dt | dE/dt | delta_E | E_exchange | max_spin_ang_exchange | stage_max_spin_ang_exchange | run_max_spin_ang_exchange | E_demag | E_zeeman | iteration | stage_iteration | stage | mx | my | mz | last_time_step | t | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
195 | -2.676342e-18 | 3736.0 | 1162.540569 | -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.312851 | -1.935225e-09 | -2.621610e-21 | 8.637446e-20 | 4.157988 | 4.157988 | 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.782329 | -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.611492 | -2.053171e-09 | -2.832075e-21 | 7.586011e-20 | 4.669400 | 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.086303 | -1.931645e-09 | -2.687311e-21 | 7.141300e-20 | 4.354488 | 4.669400 | 29.612838 | 9.291427e-19 | -3.716389e-18 | 800.0 | 3.0 | 199.0 | -0.983765 | 0.133793 | 0.042832 | 1.375197e-12 | 1.000000e-09 |
Finally, we want to plot the average magnetisation configuration my
as a function of time t
:
[20]:
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:
[21]:
# 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.sel("z").resample((20, 5)).mpl.vector(ax=ax)
# add colouring for mx-component to this plot
system.m.x.sel("z").mpl.scalar(ax=ax, cmap="magma")
References#
[1] µMAG Site Directory: http://www.ctcms.nist.gov/~rdm/mumag.org.html
[2] Pandas: http://pandas.pydata.org/