DEPARTMENT OF SOLID STATE SCIENCES
DYNAMAT GROUP
MUMAX3-WORKSHOP
SESSION 1
Dr. Jonathan Leliaert,
Dr. Jeroen Mulkers
THE DYNAMAT TEAM
Research group of the department of solid-state sciences at Ghent University
Workshop instructors
Dr. Jonathan Leliaert Dr. Jeroen Mulkers Prof. Dr. Bartel Drs. Pieter Gypens
Van Waeyenberge
[Link] 2
WORKSHOP SCHEDULE
Monday 08/31, 6PM-8PM CET
Session 1: general introduction to micromagnetics in mumax3
Session 2: mumax3 ecosystem, workflow and a first simulation
Monday 09/07, 6PM-7:30PM CET
Session 3: basic examples
Homework
Monday 09/14, 6PM-7:30PM CET
Session 4: advanced features and more extensive examples
3
SESSION 1: GENERAL INTRODUCTION
Part 1: Micromagnetics recap Part 2: Design of mumax3
• Magnetization field • Mumax3 in a nutshell
• Magnetic energy • Scripting language
• Magnetization dynamics • Discretization
• Energy minimization • Shapes, geometry, regions
• Material parameters
• Initial magnetization
• Output
• Run/Relax/Minimize
• Effective field terms
• Spin transfer torques
4
MICROMAGNETISM
In ferromagnets, neighboring magnetic moments have the tendency to align
Continuum
approx.
Quasi-uniform state
Magnetization can be described by a continuous vector field
central quantity of interest Vortex
• picosecond time scale
• 1nm – 1µm length scale
Néel Skyrmion
5
MICROMAGNETISM
The physics is fully described by the total magnetic energy functional:
Crystal anisotropy
Dzyaloshinskii-Moriya
Exchange Zeeman Magneto-elasticity
Higher-order exchange
…
6
MICROMAGNETISM
Two main objectives
Time integration Energy minimization
dynamics statics
- Spin waves - Stable magnetic states
- Domain wall motion - Hysteresis curves
- Spin transfer torques - Phase diagrams
- Vortex excitation - Domain wall profiles
- …. - …
7
MICROMAGNETISM Time integration
dynamics
The magnetization dynamics is described by the Landau-Lifshitz-Gilbert (LLG) equation
𝛾 1 𝛿𝐸
𝒎̇ = − 𝒎×𝑯 +𝛼𝒎 × 𝒎×𝑯 with 𝑯 =−
1+𝛼 𝜇 𝑀 𝛿𝒎
precession damping
Computation of the magnetization dynamics:
LLG
𝐸 𝒎, 𝑡 → 𝐻 𝒓 𝒎(𝒓, 𝑡)
8
MICROMAGNETISM Time integration
dynamics
LLG adaptation 1
Spin transfer torques (Zhang-Li[1], Slonczewski[2,3] )
𝛾
𝒎̇ = − 𝒎×𝑯 +𝛼𝒎 × 𝒎×𝑯 +𝜏
1+𝛼
LLG adaptation 2
Additional random effective field term which scales with the temperature[4]
𝑯 → 𝑯 +𝑯 𝑯 (𝒓, 𝑡) = 0
2𝑘 𝑇𝛼
𝑯 (𝒓, 𝑡), 𝑯 (𝒓 , 𝑡′) = 𝛿 𝒓 − 𝒓 𝛿(𝑡 − 𝑡 )
𝑀𝛾
[1] Physical Review Letters 93(12), 127204 (2004). [3] Physical Review B 70(17), 172405 (2004).
9
[2] Journal of Magnetism and Magnetic Materials 159(1-2), L1–L7 (1996). [4] Journal of Applied Physics 112(2), (2012).
MICROMAGNETISM
Energy minimization
statics
Two approaches to minimize the free energy functional numerically:
1. Standard minimization scheme (e.g. steepest gradient)
2. LLG equation with strong damping (=removing the precession term)
𝒎̇ = −𝒎 × 𝒎 × 𝑯
10
WHAT IS MUMAX3?
• Free finite-difference based micromagnetic simulation package
• GPU-accelerated nvidia GPU required
• Developed at DyNaMat (Ugent) by Arne Vansteenkiste
• Latest official release mumax3.10 (Aug 13, 2020)
• Active community [Link]/forum/#!forum/mumax2
• Documented API [Link]
• Open source (GPLv3) [Link]/mumax/3
• Mainly written in Go
• CUDA C kernels for heavy lifting
• Scripting language + Web GUI
• Well tested (unit tests + NIST standard problems)
11
SCRIPTING LANGUAGE
Mumax3’s scripting language is a subset of golang
// saturation magnetization
Msat = 5e6
// declare new variable
Freq := 1e9
for i:=0; i<10; i++ {
print(i)
}
if 1+8 == 9 {
print(“Of course 1+8=9”)
}
12
DISCRETIZATION
• Rectangular simulation box (origin in the center)
• Single regular rectangular grid
• Uniform magnetization inside cell 𝒎 = 𝒎(𝑥 , 𝑦 , 𝑧 )
• Cell size < exchange length
setgridsize(256,64,1)
setcellsize(1e-9,1e-9,1e-9)
• PBC values are number of virtual repetitions of the
simulation box used to calculate dipolar interactions
setpbc(4,0,0)
13
DISCRETIZATION
The cuda fft library (used for the computation of the demag field) is
highly optimized for grid size dimensions with small prime factors.
Tip: try to use grid size dimensions which are ‘7-smooth’
Example of good and bad grid size dimensions:
190 = 2 ⋅ 5 ⋅ 19 191 = 191 192 = 2 ⋅ 3
14
SHAPES
• A shape can be considered as a function 𝑓: ℝ → 𝑡𝑟𝑢𝑒, 𝑓𝑎𝑙𝑠𝑒 where
𝑡𝑟𝑢𝑒 if (𝑥, 𝑦, 𝑧) in shape
𝑓 𝑥, 𝑦, 𝑧 =
𝑓𝑎𝑙𝑠𝑒 otherwise
• Most shapes do not depend on the grid
• Default location: center of universe (0,0,0)
• Large set of predefined basic shapes
• Define new shapes by combining and modifying the basic shapes
• Shapes are useful for different tasks:
The geometry
Defining regions
To set locally an initial magnetization
15
SHAPES
Shapes Shape methods // Rotated cheese example
Cell(j,k,l) d := 200e-9
Circle(diameter) Transl(dx,dy,dz) sq := square(d)
Cone(diameter,height) Scale(sx,sy,sz)
Cuboid(Lx,Ly,Lz) RotX(angle) h := 50e-9
Cylinder(diameter,height) RotY(angle) hole := cylinder(h, h)
Ellipse(a,b) RotZ(angle) hole1 := [Link](100e-9, 0, 0)
Ellipsoid(a,b,c) Repeat(dx,dy,dz) hole2 := [Link](0, -50e-9, 0)
ImageShape(filename)
Layer(i) Add(shape) cheese := [Link](hole1).sub(hole2)
Layers(i1,i2) Sub(shape) cheese = [Link](pi/6)
Rect(Lx,Ly) Inverse()
Square(L) Intersect(shape)
Xrange(xmin,xmax) Xor(shape)
Yrange(ymin,ymax)
Zrange(zmin,zmax)
rotated cheese
16
GEOMETRY
Optionally a magnet Shape other than the full simulation box can be specified
// Ring geometry example
setGridsize(100, 100, 10)
setCellsize(1e-9,1e-9,1e-9)
ring := circle(100e-9).sub(circle(50e-9))
setgeom(ring)
save(geom)
17
REGIONS
• 256 regions in total (index 0→ 255)
• Each cell is assigned to a single region (default region id is 0)
• Each region has its own set of material parameters
• Two ways to set the region id in cells:
1. Set region id of a single cell
2. Set region id of all cells in a shape
SetGridSize(32,32,1)
SetCellSize(1,1,1)
// Set region id of cells in a circle to 1
DefRegion(1, circle(30))
// Set region id of cell (10,10,0) to 2
DefRegionCell(2, 10, 10, 0)
Save(regions)
18
MATERIAL PARAMETERS
‒ Material parameters are assigned to the 256 regions.
‒ Material parameters can be functions of time
‒ There are vector and scalar material parameters
‒ Material parameters are predefined, they can not be created
// Assigning to a material parameter sets a value in all regions:
Msat = 800e3
AnisU = vector(1, 0, 0)
// When regions are defined, they can also be set region-wise:
[Link](0, 800e3)
[Link](1, 540e3)
// Material parameters can be functions of time as well:
f := 500e6
Ku1 = 500 * sin(2*pi*f*t)
19
MATERIAL PARAMETERS: EXCITATIONS
‒ An excitation is a regional material parameter
‒ Additionally, one can add an arbitrary number of time- and space-dependent
vector fields of the form :
𝑔 𝑥 ,𝑦 ,𝑧 ∗𝑓 𝑡
B_ext = vector(0,0,1)
B_ext.Add(LoadFile("[Link]"), sin(2*pi*f*t))
B_ext.removeExtraTerms()
20
INITIAL MAGNETIZATION
Different ways to set the magnetization:
m = config
[Link](filename)
[Link](regionId, config)
[Link](shape, config)
[Link](j,k,l, vector)
Here, a ‘config’ is an object which represents a magnetization configuration
21
INITIAL MAGNETIZATION
Config Config methods m = uniform(1,1,0) m=Vortex(1,-1).Add(0.1,randomMag())
Uniform(mx, my, mz)
RandomMag()
RandomMagSeed(seed)
TwoDomain( Transl(dx,dy,dz)
mx1, my1, mz1, Scale(sx,sy,sz) M = BlochSkyrmion(1, 1) m=Vortex(1,-1).transl(100e-9,50e-9,0)
my2, my2, mz2, Add(ratio, config)
mx3, my3, mz3) RotZ(angle)
Vortex(circ, pol)
AntiVortex(circ, pol)
VortexWall(mxLeft, mxRight, circ, pol)
NeelSkyrmion(charge, pol)
BlochSkyrmion(charge, pol) m = uniform(1, 1, 1)
Conical(kVec, coneDir, coneAngle) [Link](cylinder(400e-9,100e-9), vortex(1,-1))
Helical(kVec)
22
OUTPUT
3 output media:
• log file for input, logging and printing
• [Link] ( t, mx, my, mz, … )
tableadd(E_total)
tableaddvar(myVar,”myVar”,”unit”)
tablesave() // write single line
tableautosave(1e-12) // write periodically
• .ovf files for scalar and vector fields
save(Edens_total)
saveas(Edens_total,”[Link]”)
autosave(Edens_total, 1e-10) // write perdiodically
23
RUN/RELAX/MINIMIZE
• Solving the LLG equation (time integration) // WARNING: ADVANCED SETTINGS
// In most cases, these settings can be ignored
run(timeperiod)
steps(100) // Set the solver:
// 1:Euler, 2:Heun, 3:Bogaki-Shampine, 4: Runge-Kutta(RK45),
runWhile(condition)
// 5:Dormand-Prince(the default), 6:Fehlberg, -1:Backward Euler
SetSolver(5)
// set timestep
fixdt = 0 // if 0 (default): use adaptive timestep
• Minimizing the energy
// Advanced settings for adaptive timestep (default values are given)
relax() // LLG without precession Headroom = 0.8 // headroom dt correction
Minimize() // steepest descent[1] MaxDt = 0 // if 0, no maximal timestep
MinDt = 0 // if 0, no minimal timestep
MaxErr = 1e-5 // maximum allowed error/step
// Advanced settings for minimizer (default values are given)
MinimizerSamples = 10 //Number of max dM for convergence check
MinimizerStop = 1e-6 //Stopping max dM for Minimize
[1] JOURNAL OF APPLIED PHYSICS115, 17D118 (2014) 24
THE INTERACTIONS
Effective field terms Spin transfer torques
• Demagnetization • Zhang-Li STT
• Exchange • Slonczewski STT
• Anisotropy
• Dzyaloshinskii-Moriya
• External field
• Thermal field
• Custom field
25
DEMAGNETIZATION FIELD TERM
Regional Material Parameters
Demagnetization energy density
Msat Saturation magnetization (A/m)
𝜇 NoDemagSpins Disable magnetostatic interaction per region
𝜀= − 𝑀 𝒎⋅𝑯 (default=0, set to 1 to disable).
2
Output Quantities
B_demag Magnetostatic field (T)
Edens_demag Magnetostatic energy density (J/m3)
E_demag Magnetostatic energy (J)
Other functionalities
EnableDemag Enables/disables demag (default=true)
SetPBC Sets the number of repetitions in X,Y,Z to create periodic
boundary conditions. The number of repetitions
determines the cutoff range for the demagnetization.
DemagAccuracy Controls accuracy of demag kernel
26
EXCHANGE FIELD TERM
Regional Material Parameters
Exchange energy density
Aex Exchange stiffness (J/m)
Msat Saturation magnetization (A/m)
𝜀 = A ∇𝒎
Output Quantities
Harmonic mean for inter-region
B_exch Exchange field (T)
exchange coupling (default behavior)
Edens_exch Total exchange energy density (including DMI) (J/m3)
𝐴 𝐴 E_exch Total exchange energy (including DMI) (J)
𝐴 𝑀 𝑀 MaxAngle Maximum angle between exchanged coupled spins
=2 (rad)
𝑀 𝐴 𝐴
𝑀 + 𝑀
Other functionalities
1 1 1 2 2 2 Ext_InterExchange Sets exchange coupling between two regions
Ext_ScaleExchange Re-scales exchange coupling between two regions
27
ANISOTROPY FIELD TERM
Uniaxial anisotropy energy density Regional Material Parameters
anisU Uniaxial anisotropy direction
𝜀 =−𝐾 𝑢⋅𝒎 −𝐾 𝑢⋅𝒎 Ku1, Ku2 Uniaxial anisotropy constants (J/m3)
AnisC1, AnisC2 Cubic anisotropy directions
Kc1, Kc2, Kc3 Cubic anisotropy constants (J/m3)
Similar expression for cubic Msat Saturation magnetization (A/m)
anisotropy energy density [1]
Output Quantities
B_anis Anisotropy field (T)
Edens_anis Total anisotropy energy density (including DMI) (J/m3)
E_anis Total anisotropy energy (including DMI) (J)
[1] "The design and verification of mumax3", AIP Advances 4, 107133 (2014). 28
DZYALOSHINSKII-MORIYA FIELD TERM
Interfacially-induced DMI energy density
Regional Material Parameters
Dind Interfacially-induced DMI strength (J/m2)
𝜀 = 𝐷 [𝑚 ∇ ⋅ 𝒎 − 𝒎 ⋅ ∇ 𝑚 ] Dbulk Bulk DMI strength (J/m2)
Msat Saturation magnetization (A/m)
Bulk DMI energy density
Output Quantities
B_exch Exchange field (including DMI) (T)
𝜀 = 𝐷 𝒎 ⋅ (∇ × 𝒎)
Edens_exch Total exchange energy density (including DMI) (J/m3)
E_exch Total exchange energy (including DMI) (J)
Note:
Other functionalities
Only single DMI type allowed at once Ext_InterDind Sets Dind coupling between two regions
Ext_ScaleDind Re-scales Dind coupling between two regions
OpenBC Use open boundary conditions (default=false, use
Neumann BC). This setting is only relevant for DMI.
29
EXTERNAL FIELD TERM
Zeeman energy density Regional Material Parameters
Msat Saturation magnetization (A/m)
𝜀 = −𝑀 𝒎 ⋅ 𝑩 Excitation
B_ext Externally applied field(T)
Output Quantities
B_ext Externally applied field(T)
Edens_zeeman Zeeman energy density (J/m3)
E_zeeman Zeeman energy (J)
30
THERMAL FIELD TERM
Regional Material Parameters
Temp Temperature (K)
alpha Damping parameter
𝑯 (𝒓, 𝑡) = 0
Msat Saturation magnetization (A/m)
Output Quantities
𝑯 (𝒓, 𝑡), 𝑯 (𝒓 , 𝑡′)
B_therm Thermal field (T)
2𝑘 𝑇𝛼 Edens_therm Thermal energy density (J/m3)
= 𝛿 𝒓 − 𝒓 𝛿(𝑡 − 𝑡 )
𝑀𝛾 E_therm Thermal energy (J)
Other functionalities
ThermSeed Sets random seed for thermal noise
31
CUSTOM FIELD TERMS
Output Quantities
B_custom User-defined field (T)
Edens_custom Total energy density of user-defined field (J/m3)
E_custom Total energy of user-defined field (J)
Other functionalities
AddEdensTerm Add a custom energy density term
AddFieldTerm Add a custom effective field term
RemoveCustomFields Remove all custom fields
32
ZHANG-LI STT
Regional Material Parameters
Zhang-Li spin-transfer torque
Pol Electrical current polarization
1 + 𝜉𝛼 xi Non-adiabaticity of spin-transfer-torque
𝝉= 𝒎× 𝒎× 𝒖⋅∇ 𝐦
1+𝛼 alpha Damping parameter
𝜉−𝛼 Msat Saturation magnetization (A/m)
+ 𝒎× 𝒖⋅∇ 𝒎
1+𝛼 Excitation
J Electrical current density (A/m2)
with
𝜇 𝑃 Output Quantities
𝒖= 𝒋
2𝑒𝛾 𝑀 (1 + 𝜉 ) STTorque Spin-transfer torque/γ0 (T)
Other functionalities
DisableZhangLiTorque Disable Zhang-Li torque (default=false)
33
SLONCZEWSKI STT
Regional Material Parameters
Slonczewski spin-transfer torque Pol Electrical current polarization
Lambda Slonczewski Λ parameter
𝜖 + 𝛼𝜖′
𝝉=𝛽 𝒎× 𝒎 ×𝐦 EpsilonPrime Slonczewski secondairy STT term ε'
1+𝛼
alpha Damping parameter
𝜖′ − 𝛼𝜖
−𝛽 𝒎×𝒎 Msat Saturation magnetization (A/m)
1+𝛼
FreeLayerThickness Slonczewski free layer thickness (if set to zero (default),
with then the thickness will be deduced from the mesh size) (m)
Excitation
𝑗ℏ
𝛽= J Electrical current density (A/m2)
𝑀 𝑒𝑑
𝑃Λ FixedLayer Slonczewski fixed layer polarization
𝜖=
Λ + 1 + (Λ − 1)(𝒎 ⋅ 𝒎 )
Output Quantities
STTorque Spin-transfer torque/γ0 (T)
Other functionalities
DisableSlonczewskiTorque Disables Slonczewski torque (default=false)
FixedLayerPosition Position of the fixed layer: FIXEDLAYER_TOP,
FIXEDLAYER_BOTTOM (default=FIXEDLAYER_TOP) 34
[Link]
[Link]@[Link]
[Link]@[Link]