1   NMAG User Manual (v0.2)¶

Importing nmag¶

First we need to import the nmag module, and any subpackages of nmag that we want to use. (In this basic example, this is just the SI module for dimensionful physical quantities).

import nmag
from nmag import SI

Creating the simulation object¶

Next, we need to create a simulation object. This will contain and provide information about our physical system.

sim = nmag.Simulation()

Defining (magnetic) materials¶

After importing the nmag module into Python’s workspace and creating the simulation object sim, we need to define a material using nmag.MagMaterial. We give it a name (as a Python string) which in this case we choose to be "Py" (a common abbreviation for PermAlloy) and we assign a saturation magnetisation and an exchange coupling strength.

Py = nmag.MagMaterial(name = 'Py',
                      Ms = SI(1e6, 'A/m'),
                      exchange_coupling = SI(13.0e-12, 'J/m'))

The name of the material is important, as we may want to simulate systems made up of multiple different materials, and the material name will be used as a postfix to the name of some Fields and subfields. The output files will also use that name to label output data. Names must be alphanumeric (i.e. formed exclusively out of the characters in the set 0-9\_a-zA-Z) here.

Rather than representing dimensionful physical quantities as numbers, nmag uses a special object class, the “SI object”. The underlying rationale is that this allows automated detection of mismatches of physical dimensions. If some physical parameter is given to nmag in a dimension different from the expected one, nmag will detect this and report an error. Also, any nmag output [e.g. a three-dimensional VTK visualisation file] will provide a sufficient amount of contextual information to clarify the physical meaning (i.e. dimensions) of numerical data.

We thus express the saturation magnetisation in Ampere per meter (Ms = SI(1e6,"A/m")) and the exchange coupling constant (often called A in micromagnetism) in Joules per meter (exchange_coupling = SI(13.0e-12, "J/m")). (Note that these are not the true physical parameters of PermAlloy, but have been chosen ad hoc for the sake of providing a simple example!)

Loading the mesh¶

The next step is to load the mesh.

sim.load_mesh('sphere1.nmesh.h5',
              [('sphere', Py)],
              unit_length = SI(1e-9, 'm'))

The first argument is the file name ("sphere1.nmesh.h5"). The second argument is a list of tuples which describe the domains (also called regions) within the mesh. In this example we have a one-element list containing the 2-tuple ("sphere", Py). The left element of this pair, "sphere", is a string (of the user’s choice) and this is the name given to mesh region 1 (i.e. the space occupied by all simplices that have the region id 1 in the mesh file).

[This information is currently only used for debugging purposes (such as when printing the simulation object).]

The second part of the tuple is the MagMaterial object that has been created in Defining (magnetic) materials and bound to the variable Py. This object determines the material properties of the material in this domain; in this example, we have specified the properties of PermAlloy.

The third argument to load_mesh is an SI object which defines what physical distance should be associated with the length 1.0 as given in the mesh file. In this example, the mesh has been created in nanometers, i.e. the distance 1.0 in the mesh file should correspond to 1 nanometer in the real world. We thus use a SI object representing 1 nm.

Setting the initial magnetisation¶

To set the initial magnetisation, we use the set_m method.

sim.set_m([1,0,0])

The field m describes the direction of magnetisation (as a field of normalised vectors) whereas the field M contains the magnetisation with its proper magnitude. So, |M| is the saturation magnetisation (in Amperes per meter), whereas m is dimensionless with |m|=1.0. There are different ways to set a particular magnetisation, in the simplest case of a homogeneously magnetised body, it is sufficient to provide the magnetisation vector. So, in this example, we provide a unit vector pointing in positive x-direction. (We could provide a vector with non-normalised magnitude, which would be normalised automatically. This is convenient for, say, setting an initial magnetisation in the x-y-plane with a 45 degree angle towards the x axis by specifying [1,1,0]).

Setting the external field¶

We can set the external field using the set_H_ext command

sim.set_H_ext([0,0,0], SI('A/m'))

In contrast to set_m, this method takes two arguments. The first defines numerical values for the direction and magnitude of the external field. The second determines the meaning of these numerical values using an SI object. Suppose we would like an external field of 1e6 A/m acting in the y-direction, then the command would read: sim.set_H_ext([0,1e6,0],SI(1,"A/m")). However, we could also use sim.set_H_ext([0,1,0],SI(1e6,"A/m")).

The default value for the external field is [0,0,0] A/m, so for this example, we could have omitted the set_H_ext command altogether.

Extracting and saving data¶

We have three different ways of extracting data from the simulation:

1. saving averaged values of fields (which can be analysed later)
2. saving spatially resolved fields (which can be analysed later)
3. extracting field values at arbitrary positions from within the running program

In this basic example, we demonstrate the use of all three methods:

sim.save_data()

H_demag = sim.probe_subfield_siv( 'H_demag', [x,0,0] )

x = -1e-08 : H_demag =  None
x = -9e-09 : H_demag =  [-329655.76203912671, 130.62999726469423, 194.84338557811344]
x = -8e-09 : H_demag =  [-329781.46587966662, 66.963624669268853, 137.47161381890737]
x = -7e-09 : H_demag =  [-329838.57852402801, 181.46249265908259, 160.61298054099865]
x = -6e-09 : H_demag =  [-329899.63327447395, 131.06488858715838, 71.383139326493094]
x = -5e-09 : H_demag =  [-329967.79622912291, 82.209856975234786, -16.893046828024836]
x = -4e-09 : H_demag =  [-329994.67306536058, 61.622521557150371, -34.433041910642359]
x = -3e-09 : H_demag =  [-329997.62759666931, 23.222244635691535, -65.991127111463769]
x = -2e-09 : H_demag =  [-330013.90370482224, 10.11035370824321, -61.358763616681067]
x = -1e-09 : H_demag =  [-330023.50844056415, -6.9714476825652287, -54.900260456937708]
x = 0.0   : H_demag =  [-330030.98847923806, -26.808832466764223, -48.465748009067141]
x = 1e-09 : H_demag =  [-330062.38479507214, -38.660812022013424, -42.83439139610747]
x = 2e-09 : H_demag =  [-330093.78111090627, -50.512791577262625, -37.2030347831478]
x = 3e-09 : H_demag =  [-330150.72580001026, -64.552170478617398, -23.120555702674721]
x = 4e-09 : H_demag =  [-330226.19050178828, -77.236085707456397, -5.5373829923226916]
x = 5e-09 : H_demag =  [-330304.59300913941, -90.584413821813229, 14.090609104026118]
x = 6e-09 : H_demag =  [-330380.1392610991, -115.83746059068679, 37.072085708324757]
x = 7e-09 : H_demag =  [-330418.85831447819, -122.47512022500726, 62.379121138009992]
x = 8e-09 : H_demag =  [-330476.40747455234, -110.84257225592108, 108.06217226524763]
x = 9e-09 : H_demag =  [-330500.20126762061, -68.175725285038382, 162.46166752217249]
x = 1e-08 : H_demag =  [-330517.86675206106, -24.351273685146875, 214.40344001233677]

Saving spatially resolved data¶

The command

sim.save_data(fields='all')

will save full spatially resolved data on all fields (see Fields and subfields) for the current configuration into a file with name sphere1_dat.h5. (It will also save the spatially averaged values as described in Saving averaged data.) Whenever the save_data function is called, it will write the averaged field values into the Data files (.ndt) file. This name is, by default, based on the name of the simulation script, but can be overridden with an optional argument to the Simulation constructor. The data in this file are kept in some compressed binary format (built on the `hdf5`_ standard) and can be extracted and converted later using the nmagpp tool.

For example, we can extract the magnetisation field from this file with the command:

$ nmagpp --dump sphere1

However, here we are interested in creating a vtk from the saved data file for visualisation. We use:

$ nmagpp --vtk sphere1.vtk sphere1

where sphere1.vtk is the base name of the vtk file that is to be generated.

In this manual, we use MayaVi as the visualisation tool for vtk files but there are others available (see vtk).

Starting MayaVi with the command mayavi -d sphere1-000000.vtk will load our simulation data. Using the pull-down menu Visualize -> Modules -> VelocityVector will then tell MayaVi to display the magnetisation vector field. (Likewise, we can use Visualize -> Modules -> Axes to add a 3d coordinate system to the visualization):

The magnetisation is pointing in positive x-direction because we initialised the magnetisation in this orientation by issuing the command sim.set_m([1,0,0]).

The Configure Data button in the DataVizManager section of MayaVi’s user interface allows to select:

• a vector field and
• a scalar field

which provide the data that is used for subsequent visualisation modules. Above, we have used the m_Py vector field.

The demagnetisation field should point in the opposite direction of the magnetisation. However, let’s first create a colour-coded plot of the scalar magnetic potential, phi, from which the demag field is computed (by taking its negative gradient):

We first need to select phi as the data source for ‘scalar’ visualisation modules: Through clicking on the Configure Data button in the DataVizManager section of MayaVi’s user interface, we can select phi<A> as the data source for scalar visualisations. (The <A> simply indicates that the units of the potential are Ampere).

To show the scalar potential, we use the Visualize->Module->SurfaceMap module.

We can see that the potential varies along the x-direction. The legend at the bottom of the figure shows the colour code used. We can also see from the legend title that the physical dimension of the potential phi is Ampere (this is the <A>).

Unless the user specifies a particular request for physical dimensions, the following rules apply for vtk files:

• position are given in the same coordinates as the mesh coordinates (that is why in this example, the x, y and z axis have values going from -10 to 10).
• all field data are given in SI units.

The next plot shows the demag field (the vectors) together with isosurfaces of the magnetic potential:

It can be seen that the isosurfaces are completely flat planes (i.e. the potential is changing only along x) and the demagnetisation field is perpendicular to the isosurfaces. The color bar on the left refers to the magnitude of the demagnetisation field which is expressed in Ampere per meter, as can be seen from the label <A/m>. (Note that all the H_demag arrows are colored red as they have identical length.)

Mesh generation¶

While it is down to the mesh generation software (see also Finite element mesh generation) to explain how to generate finite element meshes, we briefly summarize the steps necessary to create a mesh for this example in `Netgen`_, and how to convert it into an nmesh mesh.

1. The finite element method requires the domain of interest to be broken down into small regions. Such a subdivision of space is known as a mesh or grid. We use `Netgen`_ to create this mesh. Netgen reads a geometry file describing the three-dimensional structure. To create the mesh used here, we can start Netgen and load the geometry file by using the menu: File-> Load Geometry. We then tell Netgen that we like the edge length to be shorter than 3 by going to Mesh->Meshing Options->Mesh Size and enter 3.0 in the max mesh-size box. Then a click on the Generate Mesh button will generate the mesh. Finally, using File->Export will save the mesh as a “neutral” file (this is the default) under the name bar30_30_100.neutral. (We provide a gzipped version of this file for completeness.)

2. This neutral file needs to be converted into a nmesh file. We do this using the command:

   $ nmeshimport --netgen bar30_30_100.neutral bar30_30_100.nmesh.h5
   

   By providing the .h5 extension, we tell nmeshimport to write a compressed mesh file which is significantly smaller than an ascii file (see mesh file size).

The generated mesh looks like this:

We can examine the mesh using nmeshpp to obtain information about mesh quality, the statistical distribution of edge lengths, the overall number of points and elements etc.

If you like to script the mesh generation starting from a Netgen geometry file and ending with the nmesh file, you could use (for the example above), the following shell commands:

netgen -geofile=bar30_30_100.geo -meshfiletype="Neutral Format" -meshfile=bar30_30_100.neutral -batchmode
nmeshimport --netgen bar30_30_100.neutral bar30_30_100.nmesh.h5

The simulation¶

Having obtained the mesh file bar30_30_100.nmesh.h5, we can use the program bar30_30_100.py to run the simulation:

import nmag
from nmag import SI

mat_Py = nmag.MagMaterial(name="Py",
                          Ms=SI(0.86e6,"A/m"),
                          exchange_coupling=SI(13.0e-12, "J/m"),
                          llg_damping=0.5)

sim = nmag.Simulation("bar")

sim.load_mesh("bar30_30_100.nmesh.h5",
              [("Py", mat_Py)],
              unit_length=SI(1e-9,"m"))

sim.set_m([1,0,1])

dt = SI(5e-12, "s") 

for i in range(0, 61):
    sim.advance_time(dt*i)                  #compute time development

    if i % 10 == 0:                         #every 10 loop iterations, 
        sim.save_data(fields='all')         #save averages and all
                                            #fields spatially resolved
    else:
        sim.save_data()                     #otherwise just save averages

As in Example: demag field in uniformly magnetised sphere, we start by importing |nmag| and creating the material object.

import nmag
from nmag import SI

mat_Py = nmag.MagMaterial(name="Py",
                          Ms=SI(0.86e6,"A/m"),
                          exchange_coupling=SI(13.0e-12, "J/m"),
                          llg_damping=0.5)

We set the llg_damping parameter to 0.5. As this is a dimensionless parameter, we can pass a number. Alternatively, we may give it as SI(0.5). (Note that in this example, we give the appropriate physical value for the saturisation magnetisation of PermAlloy.)

The next line creates the simulation object:

sim = nmag.Simulation("bar")

Here, we provide a name for the simulation, which is bar. This will be used as the stem of the name of any data files that are being written. If this name is not specified (as in Example: demag field in uniformly magnetised sphere), it defaults to the name of the file that contains the script (but without the .py extension).

Next, we load the mesh file, and set the initial (normalised) magnetisation to point in the [1,0,1] direction, i.e. to have equal magnitude in the x- and z-direction and 0 in the y-direction.

sim.load_mesh("bar30_30_100.nmesh.h5",
              [("Py", mat_Py)],
              unit_length=SI(1e-9,"m"))

sim.set_m([1,0,1])

This vector will automatically be normalised within nmag, so that [1,0,1] is equivalent to the normalised vector [0.70710678,0,0.70710678].

In this example, we would like to study a dynamic process and will ask |nmag| to compute the time development over a certain amount of time dt. The line:

dt = SI(5e-12, "s")

simply creates a SI object which represents our timescale.

We then have a Python for-loop in which i will take integer values ranging from 0 to 60 for subsequent iterations. All indented lines are the body of the for-loop. (In the Python programming language, scoping is expressed through indentation rather than braces or other types of parentheses. Text editors such as Emacs come with built-in support for properly indenting Python code [by pressing the Tab key on a line to be indented].)

for i in range(0, 61):
    sim.advance_time(dt*i)

    if i % 10 == 0:
        sim.save_data(fields='all')
    else:
        sim.save_data()

In each iteration, we first call sim.advance_time(i*dt) which instructs |nmag| to carry on time integration up to the time i*dt.

The call to save_data will save the average data into the bar_dat.ndt file.

The last four lines contain an if statement which is used to save spatially resolved data every ten time steps only, and averaged data every time step. The percent operator % computes i modulo 10. This will be 0 when i takes values 0, 10, 20, 30, ... In this case, we call:

sim.save_data(fields='all')

which will save the (spatial) averages of all fields (going into the bar_dat.ndt file), and the spatially resolved data for all fields (that are saved to bar_dat.h5).

If i is not an integer multiple of 10, then the command:

sim.save_data()

is called, which only saves spatially averaged data.

Analysing the data¶

Time dependent averages¶

We first plot the average magnetisation vector against time. To see what data is available, we call ncol with just the name of the simulation (which here is bar):

$ ncol bar
 0:         #time           #<s>              0
 1:            id             <>              1
 2:          step             <>              0
 3:    stage_time            <s>              0
 4:    stage_step             <>              0
 5:         stage             <>              0
 6:    E_total_Py      <kg/ms^2>  -0.2603465789714
 7:           phi            <A>  0.0002507410390772
 8:      E_ext_Py      <kg/ms^2>              0
 9:     H_demag_0          <A/m>  -263661.6680783
10:     H_demag_1          <A/m>  -8.218106743355
11:     H_demag_2          <A/m>  -77027.641984
12:     dmdt_Py_0         <A/ms>  -8.250904652583e+15
13:     dmdt_Py_1         <A/ms>  2.333344983225e+16
14:     dmdt_Py_2         <A/ms>  8.250904652583e+15
15:   H_anis_Py_0          <A/m>              0
16:   H_anis_Py_1          <A/m>              0
17:   H_anis_Py_2          <A/m>              0
18:        m_Py_0             <>  0.7071067811865
19:        m_Py_1             <>              0
20:        m_Py_2             <>  0.7071067811865
21:        M_Py_0          <A/m>  608111.8318204
22:        M_Py_1          <A/m>              0
23:        M_Py_2          <A/m>  608111.8318204
24:     E_anis_Py      <kg/ms^2>              0
25:     E_exch_Py      <kg/ms^2>  5.046530179037e-17
26:           rho        <A/m^2>  0.03469702141876
27:       H_ext_0          <A/m>              0
28:       H_ext_1          <A/m>              0
29:       H_ext_2          <A/m>              0
30:  H_total_Py_0          <A/m>  -263661.6680783
31:  H_total_Py_1          <A/m>  -8.218106743352
32:  H_total_Py_2          <A/m>  -77027.641984
33:    E_demag_Py      <kg/ms^2>  -0.2603465789714
34:   H_exch_Py_0          <A/m>  3.301942533099e-11
35:   H_exch_Py_1          <A/m>              0
36:   H_exch_Py_2          <A/m>  3.301942533099e-11
37: maxangle_m_Py          <deg>              0
38:     localtime             <>  2007/08/15-11:16:19
39:      unixtime            <s>   1187172979.6

The meaning of the various entries is discussed in detail in section ncol. Here, we simply note that the column indices (given by the number at the beginning of every line) we are most interested in are:

• 0 for the time,
• 21 for M_Py_0 which is the x-component of the magnetisation of the Py material,
• 22 for M_Py_1 which is the y-component of the magnetisation of the Py material, and
• 23 for M_Py_2 which is the z-component of the magnetisation of the Py material,

We can use ncol to extract this data into a file data_M.dat which has the time for each time step in the first column and the x, y and z component of the magnetisation in columns 2, 3 and 4, respectively:

$ ncol bar 0 21 22 23 > data_M.txt

This creates a text file data_M.txt that can be read by other applications to create a plot.

Note, however, that the order of the entries in the ndt file is not guaranteed, i.e. the numbers corresponding to fields may change with different versions of the software, or different simulations (for example, the user may add extra fields). Therefore, the recommended approach is to directly specify the names of the columns that are to be extracted (i.e. time M_Py_0 M_Py_1 M_Py_2):

$ ncol bar time M_Py_0 M_Py_1 M_Py_2 > data_M.txt

We use the `xmgrace`_ command:

xmgrace -nxy data_M.txt

to create the following plot (manually adding the legend and axis labels):

Comparison with OOMMF and Magpar¶

We have carried out the same simulation with `Magpar`_ and `OOMMF`_. The following plot shows the corresponding OOMMF-curves (as spheres) together with |nmag|‘s results. (The Magpar curve, which is not shown here, follows the |nmag| data very closely.)

Spatially resolved fields¶

The command sim.save_data(fields='all') saves all fields into the file bar_dat.h5 (as explained, the filename is composed of the name of the simulation [here bar] and the extension _dat.h5). The code bar30_30_100.py above calls the save_data command every 10 iterations. As every dt corresponds to 0.5 picoseconds, the data hence is saved every 5 picoseconds.

We can confirm this by using the nmagpp command:

$ nmagpp --idlist bar

which produces the following output:

 id   stage   step      time  fields
 0->    1      0         0 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho
10->    1    312     5e-11 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho
20->    1    495     1e-10 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho
30->    1    603   1.5e-10 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho
40->    1    678     2e-10 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho
50->    1    726   2.5e-10 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho
60->    1    762     3e-10 E_anis E_demag E_exch E_ext E_total H_anis H_demag  ... phi pin rho

The first column is a unique identifier id for a configuration of the system. We can use the --range argument to select entries for further processing. The stage is only relevant for calculations of hysteresis curves (see Example: Simple hysteresis loop). The step is the time-stepper iteration counter for this calculation. The time is given in seconds (<s>). (Note the 5 pico-second interval between entries.) The stage, step and time data is provided for convenience. What follows is a list of fields that have been saved for each of these configurations.

We convert the first saved time step into a vtk file with base name bar_initial.vtk using

$ nmagpp --range 0 --vtk bar_initial.vtk bar

and we also convert the last saved time step at 300 picoseconds to a vtk file with base name bar_final.vtk using:

$ nmagpp --range 60 --vtk bar_final.vtk bar

The actual file names that are created by these two commands are bar_initial-000000.vtk and bar_final-000060.vtk. The appended number is the id of the saved configuration. This is useful if one wants to create vtk files for all saved configurations. For example:

$ nmagpp --vtk bar.vtk bar

will create the files:

bar-000000.vtk
bar-000010.vtk
bar-000020.vtk
bar-000030.vtk
bar-000040.vtk
bar-000050.vtk
bar-000060.vtk

Using MayaVi, we can display this data in a variety of ways. Remember that all field values are shown in SI units by default (see nmagpp), and positions are as provided in the mesh file. In this case, positions are expressed in nanometers (this comes from the unit_length=SI(1e-9,"m") expression in the sim.load_mesh() command.

This is the initial configuration with magnetisation pointing in the [1,0,1] direction:

The “final” configuration shows that the magnetisation aligns along the z-direction. The coloured surface shows the x-component of the magnetisation (and the colorbar provides the scale). It can be seen that the magnetisation at position z=100 nm goes into a flower state to minimise the overall energy. (Note that, strictly speaking, this system is not yet in a meta-stable state after 300 ps – but already quite close.):

Because we have saved all fields (not just the magnetisation), we can also study other properties. For example, the following image shows the demagnetisation field as vectors (and the legend refers to the magnitude of the vectors), as well as the magnetic scalar potential (as a stack of isosurfaces). Because the demagnetisation field is the (negative) gradient of the scalar potential, the vectors are perpendicular on the isosurfaces:

Higher level functions¶

We now have seen an overview over the fundamental commands used to set up a micromagnetic simulation and demonstrate how to advance the configuration of the system through time. In principle, this is all one would need to know to compute hysteresis loops and carry out most micromagnetic computations. However, there are more advanced functions that simplify and automatise the most frequent tasks, such as computing a hysteresis loop.

“Relaxing” the system¶

The relax command takes the current magnetisation configuration of a simulation and computes the time development until the torque on each mesh site is smaller than a certain threshold. This is useful for this particular example as we do not know for how long we need to integrate the system until it stops in a local energy minimum configuration. We can adjust the code of this example to make use of the relax command (modified source code):

import nmag
from nmag import SI, every, at

mat_Py = nmag.MagMaterial( name="Py",
                           Ms=SI(0.86e6, "A/m"),
                           exchange_coupling=SI(13.0e-12, "J/m"),
                           llg_damping=0.5)


sim = nmag.Simulation("bar_relax")

sim.load_mesh("bar30_30_100.nmesh.h5", [("Py", mat_Py)],
              unit_length=SI(1e-9, "m"))

sim.set_m([1, 0, 1])

ps = SI(1e-12,"s")
sim.relax(save = [('averages', every('time', 5*ps)),
                  ('fields', at('convergence'))])

(Note the additions to the import statement!)

The particular relax command employed here:

sim.relax(save = [('averages', every('time',5*ps)),
                  ('fields', at('convergence'))])

works as follows:

The argument save = [ ] tells relax to save data according to the instructions given in the form of a python list (i.e. enclosed by square brackets). The first relax instruction is this tuple:

('averages', every('time',5*ps)

and means that the averages should be saved every 5 picoseconds. The syntax used here breaks down into the following parts:

• 'averages' is just the keyword (a string) to say that the average data should be saved.
• every(...) is a special object which takes two parameters. They are here:
  • 'time' to indicate that something should be done every time a certain amount of simulated time has passed, and
  • 5*ps which is the amount of time after which the data should be saved again. This has to be a SI object, which we here obtain by multiplying a number (5) with the SI object ps which has been defined earlier in our example program to represent a pico-second.
  • We can provide further keywords to the every object (for example to save the data every 10 iteration steps we can use every('step', 10)).

Internally, the relax command uses the hysteresis command, so the documentation of hysteresis should be consulted for a more detailed explanation of parameters.

The second relax instruction is:

('fields', at('convergence'))

which means that the fields should be saved at convergence, i.e. when the relaxation process has finished and the magnetisation has converged to its (meta)stable configuration:

• 'fields' is a string that indicates that we would like to save all the defined fields.
• at('convergence') is a special object that indicates that this should happen exactly when the relaxation process has converged.

After running this program, we can use the ncol tool to look at the averages saved:

$ ncol bar_relax step time

gives output which starts like this:

  0              0
 82          5e-12
120          1e-11
146        1.5e-11
176          2e-11
201        2.5e-11
227          3e-11
248        3.5e-11

Here, we see the iterations on the left and the simulated time (in seconds) on the right. As requested, there is one data entry (i.e. line) every 5 picoseconds.

Note that it may happen that the system saves the data not exactly at the requested time, i.e.:

532        6.5e-11
580  7.047908066945e-11
620        7.5e-11

The middle line shows that the data has been saved when the simulated time was approximately 7.05e-11 seconds whereas we requested 7e-11 seconds. Such small deviations are tolerated by the system to improve performance [1].

From the data saved, we can obtain the following plot:

In summary, the relax function is useful to obtain a meta-stable configuration of the system. In particular, it will carry out the time integration until the remaining torque at any point in the system has dropped below a certain threshold.

“Relaxing” the system faster¶

If we are only interested in the final (meta-stable) configuration of a run, we can switch off the precession term in the Laundau Lifshitz and Gilbert equation. The MagMaterial definition in the following example shows how to do this:

import nmag
from nmag import SI, every, at

mat_Py = nmag.MagMaterial( name="Py",
                           Ms=SI(0.86e6, "A/m"),
                           exchange_coupling=SI(13.0e-12, "J/m"),
                           llg_damping=0.5,
                           do_precession=False )


sim = nmag.Simulation("bar_relax2")

sim.load_mesh("bar30_30_100.nmesh.h5", [("Py", mat_Py)],
              unit_length=SI(1e-9, "m"))

sim.set_m([1, 0, 1])

ps = SI(1e-12,"s")
sim.relax(save = [('averages', every('time', 5*ps)),
                  ('fields', at('convergence'))])

The new option is do_precession=False in the constructor of the PermAlloy material mat_Py. As a result, there will be no precession term in the equation of motion:

While the time-development of the system happens at the same time scale as for the system with the precession term (see “Relaxing” the system), the computation of the system without the precession is significantly faster (for this example, we needed about 3500 iterations with the precession term and 1500 without it, and the computation time scales similarly).

Note, that the ‘’dynamics’’ shown here are of course artificial and only used to obtain a meta-stable physical configuration more efficiently!

Decreasing execution time¶

Note that the execution time can generally be reduced significantly by decreasing the tolerances for the time integrator. In short, one has to use the set_params function (after set_m has been called). Decreasing the requested accuracy will of course make the simulation results less accurate but this is often acceptable. An example of how to use the set_m function and detailed discussion of the micromagnetic example shown in this section for a variety of tolerance values is given in the section Example timestepper tolerances.

Hysteresis simulation script¶

To compute the hysteresis loop for the ellipsoid, we use the script ellipsoid.py:

import nmag
from nmag import SI, at

#create simulation object
sim = nmag.Simulation()

# define magnetic material
Py = nmag.MagMaterial(name="Py",
                      Ms=SI(1e6,"A/m"),
                      exchange_coupling=SI(13.0e-12, "J/m"))

# load mesh: the mesh dimensions are scaled by 0.5 nm
sim.load_mesh("ellipsoid.nmesh.h5",
              [("ellipsoid", Py)],
              unit_length=SI(1e-9,"m"))

# set initial magnetisation
sim.set_m([1.,0.,0.])

Hs = nmag.vector_set(direction=[1.,0.01,0],
                     norm_list=[ 1.00,  0.95, [], -1.00,
                                -0.95, -0.90, [],  1.00],
                     units=1e6*SI('A/m'))

# loop over the applied fields Hs
sim.hysteresis(Hs, save=[('restart','fields', at('convergence'))])

As in the previous examples, we first need to import the modules necessary for the simulation. at('convergence') allows us to save the fields and the averages whenever convergence is reached. We then define the material of the magnetic object, load the mesh and set the initial configuration of the magnetisation as well as the external field.

Hysteresis loop computation¶

We apply the external magnetic fields in the x-direction with range of 1e6 A/m down to -1e6 A/m in steps of 0.05e6 A/m.

To convey this information efficiently to |nmag|, we use:

1. a direction for the applied field (here just [1,0.01,0]), (note that we have a small y-component of 1% in the applied field to break the symmetry)
2. a list of magnitudes of the field that will be multiplied with the direction vector,
3. another multiplier that defines the physical dimension of the applied fields (here 1000kA/m, given as 1e6*SI('A/m')).

Putting all this together, we obtain this command:

Hs = nmag.vector_set(direction=[1., 0.01, 0],
                     norm_list=[ 1.00,  0.95, [], -1.00,
                                -0.95, -0.90, [],  1.00],
                     units=1e6*SI('A/m'))

which computes a list of vectors Hs. Each entry in the list corresponds to one applied field.

The hysteresis command takes this list of applied fields Hs as one input parameter, and computes the hysteresis loop for these fields:

sim.hysteresis(Hs, save=[('restart', 'fields', at('convergence'))])

The save parameter is used to tell the hysteresis command what data to save, and how often. We have come across this notation when explaining the relax command in the section “Relaxing” the system of the previous example. In the example shown here, we request that the fields and the restart data should be saved at the point in time where we reach convergence. (The spatially averaged data is saved automatically to the Data files (.ndt) file when the fields are saved.) This is done in a compact notation shown above which is equivalent to this more explicit version:

sim.hysteresis(Hs,
               save=[('restart', at('convergence')),
                     ('fields',   at('convergence'))])

The compact notation can be used because here we want to save fields and restart data at the same time.

The hysteresis command computes the time development of the system for one applied field until a convergence criterion is met. It then proceeds to the next external field value provided in Hs.

We run the simulation as usual using:

$ nsim ellipsoid.py

If you have run the simulation before, we need to use the --clean switch to enforce overriding of existing data files:

$ nsim ellipsoid.py --clean

The simulation should take only a few minutes (for example 3 minutes on an Athlon64 3800+), and needs about 75MB of RAM.

If the simulation has been interrupted, it can be continued using

$ nsim ellipsoid.py –restart

Obtaining the hysteresis loop data¶

Once the calculation has finished, we can plot the graph of the magnetisation (projected along the direction of the applied field) as a function of the applied field.

We use the ncol command to extract the data into a text file named plot.dat:

$ ncol ellipsoid H_ext_0 m_Py_0 > plot.dat

Plotting the hysteresis loop with Gnuplot¶

In this example, rather than using `xmgrace`_, we show how to plot data using `Gnuplot`_:

$ gnuplot make_plot.gnu

The contents of the gnuplot script make_plot.gnu are:

set term postscript eps enhanced color
set out 'hysteresis.eps'
set xlabel 'Applied field H_x   (A/m)'
set ylabel 'M_x / M_s'
set xrange [-1050000:1050000]
set yrange [-1.2:1.2]
plot 'plot.dat' u 1:2 ti 'ellipsoid example' w lp 3

which generates the following hysteresis loop graph:

Plotting the hysteresis loop¶

To extract the data needed for plotting the hysteresis loop we proceed as explained in the previous example Example: Simple hysteresis loop. We use the ncol command and extract the data into a text file named plot.dat:

$ ncol ellipsoid H_ext_0 H_ext_1 H_ext_2 m_Py_0 m_Py_1 m_Py_2 > plot.dat

We then use `Gnuplot`_ to plot the loop:

$ gnuplot make_plot.gnu

The gnuplot script make_plot.gnu is:

set term postscript eps enhanced color
set out 'hysteresis.eps'
set xlabel 'Applied field (kA/m)'
set ylabel 'M / Ms'
versor_x = 1/sqrt(2)
versor_y = 1/sqrt(2)
versor_z = 0.0
scalar_prod(x1,x2,x3) = x1*versor_x + x2*versor_y + x3*versor_z

set mxtics 5            # minor tics and grid
set ytics 1
set mytics 5
set grid xtics ytics mxtics mytics lt -1 lw 0.5, lt 0
plot [-1050:1050] [-1.2:1.2] \
  'plot.dat' u (scalar_prod($1,$2,$3)/1000):(scalar_prod($4,$5,$6)) t 'Stoner-Wohlfarth' w lp 4

Note that within the gnuplot file, we project the magnetisation data in the [1,1,0] direction because the applied field was acting in this direction. We obtain this hysteresis loop:

The coercive field, which is located somewhere between 165 and 170 kA/m, can now be compared with the analytically known result for this particular system. To compute it, we need the demagnetizing factors Nx, Ny, Nz of the particle along the main axes. Since we deal with a prolate ellipsoid where two of the axes have the same dimension (y and z in this case), it is sufficient to compute the factor along the longest axis (x axis). The other two are easily derived from the relation Nx + Ny + Nz = 1. The expression to compute Nx is

where we call the length of the x semi-axis a, the length of the y (or z) semi-axis c, and take m to be the ratio m = a/c. Here, the value of Nx is therefore 0.1087, so we have Ny = Nz = 0.4456. With these values the shape anisotropy is easily computed according to the expression:

This gives Ha = 337 kA/m in the case of Ms = 1000 kA/m. The final step is to compute the coercive field hc using this analytical (Stoner-Wohlfarth) result:

Here, theta_0 is the angle between the easy-axis of the particle (x-axis in our case) and the direction of the applied field. Substituting theta_0 = 45 (degrees) in the formula, we obtain hc = 0.5, that is Hc = 0.5 * Ha = 168 kA/m. As we have seen before, the simulated hysteresis loop gives a value between 165 and 170 kA/m, which is in agreement with the analytical solution.

Note that this simulation is relatively slow due to a number of constraints: to get good Stoner-Wolfarth behaviour, we need to describe the shape of the ellipsoid well, and thus need a small edgelength when we generate the mesh. We further need uniform behaviour of the magnetisation, which limits the overall size of the ellipsoid. A general property of micromagnetic simulations is that the associated differential equations get stiffer if the edge lengths (or more generally: distances between neighbouring degrees of freedom) become smaller. Stiffer systems of differential equations are harder to intergrate, and thus take more time.

Thin disk hysteresis loop¶

Once the calculation has finished, we can plot the hysteresis loop, i.e. the graph of the magnetisation computed along the direction of the applied field as a function of the applied field strength.

We use the ncol command to extract the data into a text file plot.dat:

$ ncol nanodot1 H_ext_0 m_Py_0 > plot.dat

This file starts as follows:

      1000000  0.9995058139817 
      1000000  0.9995058139817 
       900000  0.9994226410102 
       900000  0.9994226410102 
       800000  0.9993139080655 

We use `Gnuplot`_ to plot the hysteresis loop:

$ gnuplot make_plot.gnu

using the gnuplot script make_plot.gnu:

set term postscript eps enhanced color
set out 'nanodot_hyst.eps'
set xlabel 'Applied field (A/m)'
set ylabel 'M / Ms'
set xrange [-1.2e6:1.2e6]
set yrange [-1.2:1.2]
plot 'plot.dat' u 1:2 ti 'nmag' with linespoints lw 3 pt 5

The resulting graph is:

and the comparison with the `Magpar`_ data, obtained with the script make_comparison_plot.gnu:

set term postscript eps enhanced color
set out 'nanodot_comparison_hyst.eps'
set xlabel 'Applied field (kA/m)'
set ylabel 'M / Ms'
set xrange [-0.2e3:0.2e3]
set yrange [-1.2:1.2]
plot 'plot.dat' u ($1/1000):2  ti 'nmag' w lp 3, 'magpar.dat' u 1:2  ti 'magpar' w p 4

is shown here (note that the `Magpar`_ computation only shows half of the hysteresis loop.):

Here we can see a slight difference between |nmag| and `Magpar`_ in the location of the switching point, probably due to different tolerances in both programs when determining time integrator convergence.

Modifying the magnetisation¶

In step 4, we use the get_subfield command. This will return a (NumPy) array that contains one 3d vector for every Site of the finite element mesh.

In step 5, we modify the first entry in this array (which has index 0), and set its value to [0,1,0]. Whereas the magnetisation is pointing everywhere in [1,0,0] (because we have used the set_m command in the very beginning of the program, it is now pointing in the [0,1,0] at site 0.

The information, which site corresponds to which entry in the data vector, that we have obtained using get_subfield, can be retrieved from get_subfield_sites. Correspondingly, the position of the sites can be obtained using get_subfield_positions.

We now need to set this modified magnetisation vector (Step 6) using the set_m command.

If we save the data again to the file (Step 7), we can subsequently convert this to a vtk file (using, for example, nmagpp --vtk data sphere_manipulate) and visualise with MayaVi:

We can see one blue cone in the centre of the sphere - this is the one site that he have modified to point in the y-direction (whereas all other cones point in the x-direction).

As before, we can probe the fields along a line through the center of the sphere (Step 8). For the demag field we obtain:

x = -1e-08 : H_demag =  None
x = -9e-09 : H_demag =  [-333816.99138074159, -1884.643376396662, 16.665519199152595]
x = -8e-09 : H_demag =  [-334670.87148225965, -2293.608410913705, -102.38526828192296]
x = -7e-09 : H_demag =  [-335258.77403632947, -3061.1708540342884, -532.73877752122235]
x = -6e-09 : H_demag =  [-339506.72150998382, -5316.1506383768137, -969.36630578549921]
x = -5e-09 : H_demag =  [-344177.83909963415, -8732.9787600552572, -1610.433091871927]
x = -4e-09 : H_demag =  [-344725.75257842313, -16708.164927667149, -5224.2484897904633]
x = -3e-09 : H_demag =  [-337963.49070659198, -24567.078937669514, -3321.016613832679]
x = -2e-09 : H_demag =  [-321612.85117992124, -30613.873989917105, -1385.6383061516099]
x = -1e-09 : H_demag =  [-298312.3363571504, -41265.117003123923, 636.60703829516081]
x = 0.0    : H_demag =  [-273449.78240732534, -52534.176864875568, 2793.5027588779139]
x = 1e-09  : H_demag =  [-293644.21931918303, -39844.049389551074, 4310.6449471266505]
x = 2e-09  : H_demag =  [-313838.65623104072, -27153.921914226579, 5827.7871353753881]
x = 3e-09  : H_demag =  [-330296.09687372146, -21814.293451835449, 5525.7290665358933]
x = 4e-09  : H_demag =  [-343611.94111195666, -18185.932406317523, 4931.5464761658959]
x = 5e-09  : H_demag =  [-348062.40814087034, -11029.603829202088, 3781.8263522408147]
x = 6e-09  : H_demag =  [-342272.36888512014, -6604.210117819096, 50.151907623841332]
x = 7e-09  : H_demag =  [-338716.66400897497, -3860.7761876767272, 485.90273674867018]
x = 8e-09  : H_demag =  [-335656.89887674141, -2610.0345208853882, 586.74812908870092]
x = 9e-09  : H_demag =  [-334985.59512328985, -2169.9546280837162, 542.76746044672041]
x = 1e-08  : H_demag =  [-334441.59096545313, -1634.8337299563193, 627.17874011463311]

The change of the magnetisation at position [0,0,0] from [1,0,0] to [0,1,0] has reduced the x-component of the demag field somewhat around x=0, and has introduced a significant demag field in the -y direction around x=0.

Looking at the exchange field (Step 9):

x = -1e-08 : H_exch_Py =  None
x = -9e-09 : H_exch_Py =  [-1.264324643856989e-09, 0.0, 0.0]
x = -8e-09 : H_exch_Py =  [-2.0419540595507732e-10, 0.0, 0.0]
x = -7e-09 : H_exch_Py =  [-1.4334754136843496e-09, 0.0, 0.0]
x = -6e-09 : H_exch_Py =  [-2.7214181426130964e-10, 0.0, 0.0]
x = -5e-09 : H_exch_Py =  [1.6323042074911775e-09, 0.0, 0.0]
x = -4e-09 : H_exch_Py =  [-153858.81305452777, 153858.81305452611, 0.0]
x = -3e-09 : H_exch_Py =  [-972420.67935341748, 972420.67935341166, 0.0]
x = -2e-09 : H_exch_Py =  [-2445371.8369108676, 2445371.8369108611, 0.0]
x = -1e-09 : H_exch_Py =  [5283169.701234119, -5283169.7012341227, 0.0]
x = 0.0    : H_exch_Py =  [15888993.991894867, -15888993.991894867, 0.0]
x = 1e-09 :  H_exch_Py =  [8434471.7912872285, -8434471.7912872266, 0.0]
x = 2e-09 :  H_exch_Py =  [979949.59067958547, -979949.59067958279, 0.0]
x = 3e-09 :  H_exch_Py =  [-1112837.3087986181, 1112837.3087986207, 0.0]
x = 4e-09 :  H_exch_Py =  [-193877.66176242317, 193877.6617624248, 0.0]
x = 5e-09 :  H_exch_Py =  [1.6602158583197139e-09, 0.0, 0.0]
x = 6e-09 :  H_exch_Py =  [1.8844573960991853e-09, 0.0, 0.0]
x = 7e-09 :  H_exch_Py =  [-6.2460015649740799e-09, 0.0, 0.0]
x = 8e-09 :  H_exch_Py =  [-1.1231714572170603e-08, 0.0, 0.0]
x = 9e-09 :  H_exch_Py =  [-7.3643182171284044e-09, 0.0, 0.0]
x = 1e-08 :  H_exch_Py =  [-3.4351784609779937e-09, 0.0, 0.0]

We can see that the exchange field is indeed very large around x=0.

Note that one of the fundamental problem of micromagnetic simulations is that the magnetisation must not vary significantly from one site to another. In this example, we have manually violated this requirement only to demonstrate how the magnetisation can be modified, and to see that this is reflected in the dependant fields (such as demag and exchange) immediately.

Pinning simulation script¶

import nmag
from nmag import SI, si

# Create simulation object
sim = nmag.Simulation()

# Define magnetic material: PermAlloy
Py = nmag.MagMaterial(name="Py",
                      Ms=SI(0.86e6, "A/m"),
                      exchange_coupling=SI(13.0e-12, "J/m"))

# Load mesh
sim.load_mesh("sphere1.nmesh.h5", [("sphere", Py)], unit_length=SI(1e-9, "m"))

# Set initial magnetisation to +x direction
sim.set_m([1, 0, 0])

# Pin magnetisation at center in radius of 4e-9m
def pin_at_center((x, y, z)):
  if (x*x + y*y + z*z) < (4e-9)*(4e-9):
    return 0.0 # Inside the 4nm sphere -> pin
  else:
    return 1.0 # Outside -> do not pin

sim.set_pinning(pin_at_center)

# Apply external field in +y direction
unit = 0.5*si.Tesla/si.mu0 # 500mT in A/m
sim.set_H_ext([0*unit, 1*unit, 0*unit])

# Relax the magnetisation
sim.relax()

Pinning magnetisation¶

In order to allow the user to fix the magnetisation, |nmag| provides a scalar field, the so-called pinning field: its value at each site is used as a scale factor for dm/dt, hence by setting it to 0 at certain locations of the mesh we can force magnetisation to remain constant at these locations for the entire simulation.

We set the pinning field using set_pinning (which is used like set_m and set_H_ext, except that it is a scalar field whereas the latter are vector fields) such that magnetisation is fixed at sites with distance less than 4 nm from the sphere’s center. First we define a Python function which we decide to call pin_at_center:

def pin_at_center((x, y, z)):
  if (x*x + y*y + z*z) < (4e-9)*(4e-9):
    return 0.0
  else:
    return 1.0

The function is called for each site of the mesh and receives the site position as an argument, a 3-tuple (x, y, z) containing the three components x, y and z (three floating point numbers), given in metres. The function returns either 0.0 (which means the magnetisation at this position is pinned) or 1.0 (in which case there is no pinning), for the given position vector.

The formula in the if statement simply evaluates the magnitude of the vector (x, y, z) by squaring each component. This number is then compared against (4nm)^2. As a result, the magnetisation is pinned at all the mesh nodes that are located within a sphere with center (0, 0, 0) and radius 4 nm. All the nodes that are located outside this sphere can change their magnetisation as usual.

Second, we need to tell |nmag| that it should use this function to decide where the magnetisation should be pinned:

sim.set_pinning(pin_at_center)

Note the slightly counterintuitive fact that value 1 means “no pinning”.

Finally we apply an external field of 0.5 T in +y direction, and use relax to compute the equilibrium configuration.

The relax command:

sim.relax()

will save the fields and averages at convergence (this is the default of the relax command).

Visualisation¶

After running the example via nsim sphere.py we convert the equilibrium data to VTK format:

$ nmagpp --vtk=sphere.vtk sphere

We would first like to verify that the pinning field has been set up properly. Hence we use MayaVi to visualise it by showing an isosurface of the pinning field (shown in blue), together with the magnetisation vector field.

The blue blob in the center of the sphere is the collection of those tetrahedra that have corners just inside the 4nm sphere. Because we have not generated the mesh to have nodes coinciding with the 4nm sphere, the shape of the blue region is not particularly spherical.

In the above diagram, we also see the magnetisation vectors of the final configuration. Their colour corresponds to the pinning field at their location. It can be seen that the blue magnetisation vectors emerging from the central region of the sphere are all pointing (strictly) in the x-direction. The magnetisation vectors outside the blue sphere are coloured red. The applied field drives these vectors to point into the y-direction. However, the magnetisation in the centre is pinned and the exchange interaction requires a gradual spatial change of magnetisation. This explains the spatial variation of the magnetisation.

The next figure shows the same data as the last figure but in addition a ScalarCutPlane (in MayaVi terminology) has been introduced which is coloured according to the x-component of the magnetisation. Red corresponds to 1.0 and blue corresponds to 0.73 (we have not shown the legend to provide a larger main plot). This demonstrates the gradual change from the pinned magnetisation in the centre to the outside.

Uniaxial anisotropy simulation script¶

import nmag
from nmag import SI, every, at
from numpy import array
import math

# Create simulation object (no demag field!)
sim = nmag.Simulation(do_demag=False)

# Define magnetic material (data from OOMMF materials file)
Co = nmag.MagMaterial(name="Co",
                      Ms=SI(1400e3, "A/m"),
                      exchange_coupling=SI(30e-12, "J/m"),
                      anisotropy=nmag.uniaxial_anisotropy(axis=[0, 0, 1], K1=SI(520e3, "J/m^3")))

# Load the mesh
sim.load_mesh("bar.nmesh.h5", [("bar", Co)], unit_length=SI(1e-9,"m") )

# Our bar is subdivided into 3 regions:
# - region A: for x < offset;
# - region B: for x between offset and offset+length
# - region C: for x > offset+length;
# The magnetisation is defined over all the three regions,
# but is pinned in region A and C.
offset = 2e-9   # m (meters)
length = 500e-9 # m

# Set initial magnetisation
def sample_m0((x, y, z)):
  # relative_position goes linearly from -1 to +1 in region B
  relative_position = -2*(x - offset)/length + 1
  mz = min(1.0, max(-1.0, relative_position))
  return [0, math.sqrt(1 - mz*mz), mz]

sim.set_m(sample_m0)

# Pin magnetisation outside region B
def sample_pinning((x, y, z)):
  return x >= offset and x <= offset + length

sim.set_pinning(sample_pinning)

# Save the magnetisation along the x-axis
def save_magnetisation_along_x(sim):
  f = open('bar_mag_x.dat', 'w')
  for i in range(0, 504):
    x = array([i+0.5, 0.5, 0.5]) * 1e-9
    M = sim.probe_subfield_siv('M_Co', x)
    print >>f, x[0], M[0], M[1], M[2]

# Relax the system
sim.relax(save=[(save_magnetisation_along_x, at('convergence'))])

We shall now discuss the bar.py script step-by-step:

In this particular example we are solely interested in energy terms resulting from exchange interaction and anisotropy. Hence we disable the demagnetisation field as follows:

sim = nmag.Simulation(do_demag=False)

We then create the material Co used for the bar, cobalt in this case, which exhibits uniaxial_anisotropy in z direction with phenomenological anisotropy constant K1 = SI(520e3, "J/m^3"):

Co = nmag.MagMaterial(name="Co",
                      Ms=SI(1400e3, "A/m"),
                      exchange_coupling=SI(30e-12, "J/m"),
                      anisotropy=nmag.uniaxial_anisotropy(axis=[0, 0, 1], K1=SI(520e3, "J/m^3")))

After loading the mesh, we set the initial magnetisation direction such that it rotates from +z to -z while staying in the plane normal to x direction (hence suggesting the development of a Bloch type domain wall):

def sample_m0((x, y, z)):
  # relative_position goes linearly from -1 to +1 in region B
  relative_position = -2*(x - offset)/length + 1
  mz = min(1.0, max(-1.0, relative_position))
  return [0, math.sqrt(1 - mz*mz), mz]

We further pin the magnetisation at the very left (x < offset = 2 nm) and right (x > offset + length = 502 nm) of the bar (note that the pinning function may also just return a python truth value rather than the number 0.0 or 1.0):

def sample_pinning((x, y, z)):
  return x >= offset and x <= offset + length

sim.set_pinning(sample_pinning)

Finally, we relax the system to find the equilibrium magnetisation configuration, which is saved to the file bar_mag_x.dat in a format understandable by `Gnuplot`_.

Visualization¶

We can then use the following `Gnuplot`_ script to visualize the equilibrium magnetisation:

set term png giant size 800, 600
set out 'bar_mag_x.png'
set xlabel 'x (nm)'
set ylabel 'M.z (millions of A/m)'

plot [0:504] [-1.5:1.5] \
  1.4 t "" w l 0, -1.4 t "" w l 0, \
  'bar_mag_x.dat' u ($1/1e-9):($4/1e6) t 'nmag' w l 2

The resulting plot clearly shows that a Bloch type domain wall has developed:

The figure shows also that the Bloch domain wall is well localized at the center of the bar, in the region where x goes from 200 to 300 nm.

Comparison¶

After simulating the same scenario with OOMMF (see oommf/bar.mif), we can compare results using another `Gnuplot`_ script:

#set term postscript enhanced eps color
set term png giant size 800, 600
set out 'bar_mag_x_compared.png'

set xlabel 'x (nm)'
set ylabel 'M.z (millions of A/m)'

Mz(x) = 1400e3 * cos(pi/2 + atan(sinh((x - 252e-9)/sqrt(30e-12/520e3))))

plot [220:280] [-1.5:1.5] \
  1.4 t "" w l 0, -1.4 t "" w l 0, \
  'bar_mag_x.dat' u ($1/1e-9):($4/1e6) t 'nmag' w lp 2, \
  'oommf/bar_mag_x.txt'u ($1/1e-9):($4/1e6) t 'oommf' w lp 1, \
  Mz(x*1e-9)/1e6 ti 'analytical' w l 3

which generates the following plot showing good agreement of both systems:

The plot shows also the known analytical solution:

Mz(x) = Ms * cos(pi/2 + atan(sinh((x - x_wall)/sqrt(A/K1))))

The plot shows only a restricted region located at the center of the bar, thus allowing an easier comparison between the three sets of data.

Cubic anisotropy simulation script¶

import nmag
from nmag import SI, si

# Create the simulation object
sim = nmag.Simulation()

# Define the magnetic material (data from OOMMF materials file)
Fe = nmag.MagMaterial(name="Fe",
                      Ms=SI(1700e3, "A/m"),
                      exchange_coupling=SI(21e-12, "J/m"),
                      anisotropy=nmag.cubic_anisotropy(axis1=[1, 0, 0],
                                                       axis2=[0, 1, 0],
                                                       K1=SI(48e3, "J/m^3")))

# Load the mesh
sim.load_mesh("cube.nmesh", [("cube", Fe)], unit_length=SI(1e-9, "m"))

# Set the initial magnetisation
sim.set_m([0, 0, 1])

# Launch the hysteresis loop
Hs = nmag.vector_set(direction=[1.0, 0, 0.0001],
                     norm_list=[0, 1, [], 19, 19.1, [], 21, 22, [], 50],
                     units=0.001*si.Tesla/si.mu0)
sim.hysteresis(Hs)

We will now discuss the cube.py script step-by-step:

After creating the simulation object we define a magnetic material Fe with cubic anisotropy representing iron:

Fe = nmag.MagMaterial(name="Fe",
                      Ms=SI(1700e3, "A/m"),
                      exchange_coupling=SI(21e-12, "J/m"),
                      anisotropy=nmag.cubic_anisotropy(axis1=[1,0,0],
                                                       axis2=[0,1,0],
                                                       K1=SI(48e3, "J/m^3")))

We load the mesh and set initial magnetisation pointing in +z direction (that is, in a local minimum of anisotropy energy density).

Finally, we use hysteresis to apply gradually stronger fields in +x direction (up to 50 mT):

Hs = nmag.vector_set(direction=[1.0, 0, 0.0001],
                     norm_list=[0, 1, [], 19, 19.1, [], 21, 22, [], 50],
                     units=0.001*si.Tesla/si.mu0)

Note that we sample more often the region between 19 and 21 mT where magnetisation direction changes rapidly due to having crossed the anisotropy energy “barrier” between +z and +x (as can be seen in the graph below).

Analyzing the result¶

First, we extract the magnitude of the applied field and the x component of magnetisation:

ncol cube H_ext_0 M_Fe_0 > cube_hext_vs_m.txt

Then we compare the result with OOMMF’s result (generated from the equivalent scene description oommf/cube.mif) using the following `Gnuplot`_ script:

set term png giant size 800,600
set out 'cube_hext_vs_m.png'
set xlabel 'H_ext.x (A/m)'
set ylabel 'M.x (A/m)'
plot 'cube_hext_vs_m.txt' t 'nmag' w l 2,\
 'oommf/cube_hext_vs_m.txt' u ($1*795.77471545947674):2 ti 'oommf' w p 1

which gives the following result:

|Nmag| provides advanced capabilities to conveniently handle arbitrary-order anisotropy energy functions. Details can be found in the documentation of the MagMaterial class.

Arbitrary anisotropy simulation script¶

import nmag
from nmag import SI, every, at
from nsim.si_units import si
import math

# Create simulation object (no demag field!)
sim = nmag.Simulation(do_demag=False)

# Function to compute the scalar product of the vectors a and b
def scalar_product(a, b): return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

# Here we define a function which returns the energy for a uniaxial
# anisotropy of order 4.
K1 = SI(43e3, "J/m^3")
K2 = SI(21e3, "J/m^3")
axis = [0, 0, 1]        # The (normalised) axis
def my_anisotropy(m):
    a = scalar_product(axis, m)
    return -K1*a**2 - K2*a**4

my_material = nmag.MagMaterial(name="MyMat",
                               Ms=SI(1e6, "A/m"),
                               exchange_coupling=SI(10e-12, "J/m"),
                               anisotropy=my_anisotropy,
                               anisotropy_order=4)

# Load the mesh
sim.load_mesh("coin.nmesh.h5", [("coin", my_material)], unit_length=SI(1e-9, "m"))

# Set the magnetization
sim.set_m([-1, 0, 0])

# Compute the hysteresis loop
Hs = nmag.vector_set(direction=[1.0, 0, 0.0001],
                     norm_list=[-0.4, -0.35, [], 0, 0.005, [], 0.15, 0.2, [], 0.4],
                     units=si.Tesla/si.mu0)

sim.hysteresis(Hs, save=[('fields', 'averages', at('convergence'))])

We simulate the hysteresis loop for a ferromagnetic thin disc, where the field is applied orthogonal to the axis of disc. This script includes one main element of novelty, which concerns the way the magnetic anisotropy is specified. In previous examples we found lines such as:

my_material = nmag.MagMaterial(name="MyMat",
                               Ms=SI(1e6, "A/m"),
                               exchange_coupling=SI(10e-12, "J/m"),
                               anisotropy=nmag.uniaxial_anisotropy(axis=[0, 0, 1],
                                                                   K1=SI(43e3, "J/m^3"),
                                                                   K2=SI(21e3, "J/m^3")))

where the material anisotropy was specified using the provided functions nmag.uniaxial_anisotropy (uniaxial_anisotropy) and nmag.cubic_anisotropy (cubic_anisotropy). In this example we are using a different approach to define the anisotropy. First we define the function my_anisotropy, which returns the energy density for the magnetic anisotropy:

# Here we define a function which returns the energy for a uniaxial
# anisotropy of order 4.
K1 = SI(43e3, "J/m^3")
K2 = SI(21e3, "J/m^3")
axis = [0, 0, 1]        # The (normalised) axis
def my_anisotropy(m):
    a = scalar_product(axis, m)
    return -K1*a**2 - K2*a**4

Note that the function returns a SI object with units “J/m^3” (energy density). The reader may have recognised the familiar expression for the uniaxial anisotropy: in fact the two code snippets we just presented are defining exactly the same anisotropy, they are just doing it in different ways. The function scalar_product, which we have used in the second code snippet just returns the scalar product of two three dimensional vectors a and b and is defined in the line above:

def scalar_product(a, b): return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

The function my_anisotropy has to be specified in the material definition: instead of passing anisotropy=nmag.uniaxial_anisotropy(...) we just pass anisotropy=my_anisotropy to the material constructor:

my_material = nmag.MagMaterial(name="MyMat",
                               Ms=SI(1e6, "A/m"),
                               exchange_coupling=SI(10e-12, "J/m"),
                               anisotropy=my_anisotropy,
                               anisotropy_order=4)

An important point to notice is that here we also provide an anisotropy order. To understand what this number is, we have to explain briefly what is going on behind the scenes. |nsim| calculates the values of the user provided function for an appropriately chosen set of normalised vectors, it then finds the polynomial in mx, my and mz (the components of the normalised magnetisation) of the specified order, which matches the sampled values.

The strength of this approach stands in the fact that the user has to provide just the energy density for the custom anisotropy. |nsim| is taking care of working out the other quantities which are needed for the simulation, such as the magnetic field resulting from the provided anisotropy energy, which would require a differentiation of the energy with respect to the normalised magnetisation.

However the user must be sure that the provided function can be expressed by a polynomial of the specified order in mx, my and mz. In the present case we are specifying anisotropy_order=4 because the energy for the uniaxial anisotropy can be expressed as a 4th-order polynomial in mx, my and mz.

In some cases the user may find useful to know that the functions nmag.uniaxial_anisotropy and nmag.cubic_anisotropy can be added: the resulting anisotropy will have as energy the sum of the energies of the original anisotropies.

The result¶

The steps involved to extract and plot the data for the simulation discussed in the previous section should be familiar to the user at this point of the manual. We then just show the graph obtained from the results of the script coin.py.

During the switching the system falls into an intermediate state, where the magnetisation is nearly aligned with the anisotropy easy axis.

Saving the state of the simulation¶

We re-consider the cubic anisotropy example (Cubic anisotropy simulation script) and replace the last line:

sim.hysteresis(Hs)

with the following lines:

from nmag import every, at
sim.hysteresis(Hs, save=[('averages', 'fields', at('stage_end')),
                         ('restart', at('stage_end') | every('step', 1000))])

The first two lines reproduce the default behaviour: the fields and their averages are saved at the end of each stage. The third line specifies that the restart file should be saved at the end of each stage and also every 1000 steps.

For convenience the modified script cube_restartable.py is shown below:

import nmag
from nmag import SI, si

# Create the simulation object
sim = nmag.Simulation()

# Define the magnetic material (data from OOMMF materials file)
Fe = nmag.MagMaterial(name="Fe",
                      Ms=SI(1700e3, "A/m"),
                      exchange_coupling=SI(21e-12, "J/m"),
                      anisotropy=nmag.cubic_anisotropy(axis1=[1, 0, 0],
                                                       axis2=[0, 1, 0],
                                                       K1=SI(48e3, "J/m^3")))

# Load the mesh
sim.load_mesh("cube.nmesh", [("cube", Fe)], unit_length=SI(1e-9, "m"))

# Set the initial magnetisation
sim.set_m([0, 0, 1])

# Launch the hysteresis loop
Hs = nmag.vector_set(direction=[1.0, 0, 0.0001],
                     norm_list=[0, 1, [], 19, 19.1, [], 21, 22, [], 50],
                     units=0.001*si.Tesla/si.mu0)
from nmag import every, at
sim.hysteresis(Hs, save=[('averages', 'fields', at('stage_end')),
                         ('restart', at('stage_end') | every('step', 1000))])

Starting and restarting the simulation¶

We will now demonstrate how the discussed nmag script can be restarted. To do that, we will have to interrupt it artificially. We start the simulation in the usual way:

$ nsim cube_restartable.py

We interrupt the execution after the hysteresis loop has started and several stages have been computed. Do this by pressing simultaneously the keys CTRL and C (in the same terminal window where nsim was started), thus simulating what could have been the result of a power cut. We then use the command:

$ ncol cube_restartable stage step time

to see at what point of the hysteresis loop the simulation was interrupted. We obtain (for this particular interruption):

 1            330  3.320127110062e-11
 2            480  5.042492488627e-10
 3            640  9.926580643272e-10
 4            805  1.464971830453e-09
 5            980  1.927649646634e-09
 6           1150  2.406521613682e-09
 7           1340  2.882400372552e-09
 8           1515  3.371522550051e-09
 9           1705  3.863380029345e-09
10           1920  4.365560120394e-09
11           2095  4.893234441813e-09
12           2295  5.436617525896e-09
13           2480  5.997866344586e-09
14           2680  6.570733097131e-09
15           2890  7.172534305054e-09
16           3100  7.803577637245e-09
17           3315  8.462827284047e-09

The simulation was interrupted at the seventeenth stage. We now try to run the simulation again with the command:

$ nsim cube_restartable.py

obtaining the following output:

<snip>
NmagUserError: Error: Found old file ./cube_restartable_dat.ndt -- cannot proceed.
To start a simulation script with old data files present you either need
to use '--clean' (and then the old files will be deleted), or use '--restart'
in which case the run will be continued.

nsim suggests the possible alternatives. We can start the simulation from scratch with the command (but this will override any data from the previous run):

$ nsim cube_restartable.py --clean

or we can continue from the configuration which was last saved:

$ nsim cube_restartable.py --restart

Here we choose the second possibility. After the simulation has finished we issue again the command ncol cube_restartable stage step time, obtaining the following output:

    1            330  3.320127110062e-11
    2            480  5.042492488627e-10
    3            640  9.926580643272e-10
    4            805  1.464971830453e-09
    5            980  1.927649646634e-09
    6           1150  2.406521613682e-09
    7           1340  2.882400372552e-09
    8           1515  3.371522550051e-09
    9           1705  3.863380029345e-09
   10           1920  4.365560120394e-09
   11           2095  4.893234441813e-09
   12           2295  5.436617525896e-09
   13           2480  5.997866344586e-09
   14           2680  6.570733097131e-09
   15           2890  7.172534305054e-09
   16           3100  7.803577637245e-09
   17           3315  8.462827284047e-09
stage           step          #time
   <>             <>           #<s>
   18           3715  8.519843629989e-09
   19           3975  9.300878866142e-09
   ...

The two lines between stage 17 and 18 stand as a reminder that the simulation was restarted at that point. (They need to be removed manually from the cube_restartable_dat.ndt file, before ncol can work in the usual way on the ndt file.)

Idea: pass simulation object to field-setting function¶

You can simulate an applied field which both changes in space and time: this may be useful to mimic the effect of a write head on the magnetic grains of an hard disk while the head is moving. The way we do this is by changing the applied field every delta_t picoseconds. This means that the applied field won’t change continuously in time: it will be piecewise constant in time (but, in general, it can be non uniform in space). You can do something like:

import math

def set_H(sim):
  width = 10.0        # nm
  v = 100.0           # nm/ns == m/s
  H_amplitude = 0.5e6 # A/m

  t = float(sim.time/SI(1e-9, 's')) # get the time in ns
  center = (v*t, 0, 0) # center of the applied field region
  def H(r):
    x, y, z = [ri/1e-9 - ci for ri, ci in zip(r, center)]
    factor = H_amplitude*math.exp(-(x*x + y*y + z*z)/(width*width))
    return [factor, factor, factor]

  sim.set_H_ext(H, unit=SI('A/m'))

sim.relax(do=[(set_H, every('time', SI(50e-12, 's'))),
              ('exit', at('time', SI(1000e-12, 's')))])

The function set_H is called every 50 ps and does the following: it sets a new field from the function H(r). This function sets a field which directed along the direction [1, 1, 1] and almost vanishes outside a sphere with radius ~ 30.0 nm. The center of this sphere moves along the direction [1, 0, 0] with velocity 100 nm/ns, thus simulating the motion of a write head in a hard disk. Obviously the piece of code is not complete, it shows only the technique in order to have a field changing in time and space. For a complete example see the next section.

Complete example: simple moving write-head example¶

Here is a simulation of five cubes made of cobalt and a write-head which moves on the top of the cubes and applies a time-varying field in order to change their magnetisation. At the beginning the magnetisation of all the cubes is pointing in the [0, 0, 1] direction. After the write-head has passed over the cubes, the magnetisation of cube 1, 3 and 5 are switched in the opposite direction, while cube 2 and 4 have unchanged magnetisation. This is possible because the write-head field, which is space-dependent (being intense only inside a sphere of radius 15-20 nm), changes also in time. It indeed translates in space, but also change in intensity, being directed in the [0, 0, -1] direction when the sphere is at the center of cube 1, 3 and 5 and in the [0, 0, 1] direction when the center of the sphere is in cube 2 and 4.

Here is the geo file used to generate the mesh (Netgen):

<pre>
algebraic3d

# cubes
solid cube1 = orthobrick (    0, 0, 0;  20.0, 20.0, 20.0) -maxh = 2;
solid cube2 = orthobrick ( 30.0, 0, 0;  50.0, 20.0, 20.0) -maxh = 2;
solid cube3 = orthobrick ( 60.0, 0, 0;  80.0, 20.0, 20.0) -maxh = 2;
solid cube4 = orthobrick ( 90.0, 0, 0; 110.0, 20.0, 20.0) -maxh = 2;
solid cube5 = orthobrick (120.0, 0, 0; 140.0, 20.0, 20.0) -maxh = 2;

tlo cube1;
tlo cube2;
tlo cube3;
tlo cube4;
tlo cube5;

And here is the full listing of the example:

from nmag.common import *
import math

# Define magnetic material (data from OOMMF materials file)
mat_Co = MagMaterial(name="Co",
                     Ms=SI(1400e3, "A/m"),
                     exchange_coupling=SI(30e-12, "J/m"),
                     anisotropy=uniaxial_anisotropy(axis=[0, 0, 1],
                                                    K1=SI(520e3, "J/m^3")))
sim = Simulation()
sim.load_mesh("cubes.nmesh.h5",
              [('cube1', mat_Co), ('cube2', mat_Co), ('cube3', mat_Co),
               ('cube4', mat_Co), ('cube5', mat_Co)],
              unit_length=SI(1e-9, 'm'))

sim.set_m([0, 0, 1])

sim.relax(save=[('fields', at('convergence'))])

t0 = [sim.time]

def set_H(sim):
    t = float((sim.time - t0[0])/SI(1e-9, 's'))  # get time in ns
    width = 10.0                                 # nm
    v = 25.0                                     # nm/ns = m/s
    H_amplitude = 4.0e6*math.sin(math.pi*t)      # A/m
    center = (v*t, 20, 10)
    print "CENTER IN", center
    def H(r):
        x, y, z = [ri/1e-9 - ci for ri, ci in zip(r, center)]
        factor = H_amplitude*math.exp(-(x*x + y*y + z*z)/(width*width))
        return [0, 0, -factor]

    sim.set_H_ext(H, unit=SI('A/m'))

set_H(sim)

sim.set_params(stopping_dm_dt=0*degrees_per_ns)

sim.relax(save=[('fields', every('time', SI(200e-12, 's'), first=t0[0]))],
          do=[(set_H, every('time', SI(50e-12, 's'), first=t0[0])),
              ('exit', at('time', SI(6000e-12, 's')))])

Here is the magnetisation at the beginning of the simulation, after the first relax command (whose purpose is just to find the zero field magnetisation configuration):

and here is the magnetisation after the write-head has passed over the cubes:

Introduction periodic boundary conditions (“macro geometry”)¶

Concerning the simulation of periodic magnetic structures, there are a few somewhat subtle issues to be taken into account, both with respect to the demagnetising and the exchange field.

The issue with the exchange field is that we may encounter situations where the magnetic material crosses the boundary of an elementary cell: a periodic array of non-touching spheres in a cubic lattice is fundamentally different from its complement, a cubic lattice made of spherical holes, insofar as that in the latter case, it is impossible to do a simulation using periodic boundary conditions without identifying degrees of freedom that live on boundaries of the simulation cell. |Nmag| can deal with this automatically, provided the mesh file contains periodicity information, i.e. data on how to identify nodes on exterior faces.

As for the demagnetising field, the most important problem is that one cannot ignore the effect of the faraway boundaries of the system: a 100 nm x 100 nm x 100 nm cell made of magnetic material in the center of a large (three-dimensional) periodic array will experience very different demagnetising fields depending on the shape of the outer boundaries of this array. Assuming spatially constant magnetisation, if these cells form a “macroscopic” (tree-dimensional) sphere, H_demag will be -1/3 M, while for a flat box, H_demag may be very close to -M. |Nmag| takes these “macro-geometry” effects into account by allowing the user to provide a geometrical layout for a finite number (say, 100-1000) of cells that approximates the shape of the faraway outer boundary of the system.

The macro geometry approach is described in [1] which may serve as a more detailed instruction to the concept.

1d example¶

In this example, we simulate a single cell in the middle of a long one-dimensional periodic array where for the purpose of computing the demagnetising field, we take three extra copies of this cell to the left and three copies to the right along the x axis. (For real applications, one would use more copies. The only effect of additional copies are to increase the setup time needed to compute an internal boundary/boundary interaction matrix.)

The next Example: 2D periodicity demonstrates the macro geometry concept for a thin film. This is followed by the Spin-waves example which includes exchange coupling between periodic copies (and is of more practical value).

The mesh of the central simulation cell used is described in cube.geo which reads:

algebraic3d

# prism
solid prism = orthobrick (-7.50, -7.50, -7.50; 7.50, 7.50, 7.50) -maxh = 1.8000;
tlo prism;

Note that the mesh is centered around the origin. This is recommended for periodic simulations. (We need to document this better.) The resulting mesh is this (the periodic copies are not shown):

The script periodic1.py reads:

import nmag
from nmag import SI

# define magnetic material
Py = nmag.MagMaterial(name="Py",
                      Ms=SI(1e6,"A/m"),
                      exchange_coupling=SI(13.0e-12, "J/m")
                      )

# size of simulation cell, plus extra spacing
# to avoid exchange interaction across interfaces
# between repeated copies of the simulation cell.
x_lattice = 15.01  # the spacing is 0.01
y_lattice = 0.0
z_lattice = 0.0


# list to store the lattice points where the periodic
# copies will be placed 
lattice_points = []

for xi in range(-3,4):
    lattice_points.append([xi*x_lattice,0.0*y_lattice,0.0*z_lattice])

#  create data structure pbc for this macro geometry
pbc = nmag.SetLatticePoints(vectorlist=lattice_points, scalefactor=SI(1e-9,'m'))

#create simulation object,  passing macro geometry data structure
sim = nmag.Simulation(periodic_bc=pbc.structure)

# load mesh
sim.load_mesh("cube1.nmesh.h5", [("repeated-cube-1D", Py)], unit_length=SI(1e-9,"m") )

# set initial magnetisation along the periodic axis
sim.set_m([1.0,0,0])

# compute the demagnetising field
sim.advance_time(SI(0,"s"))

# probe demag field at the centre of the cube, function
# returns an SI-Value ('siv')
H_demag = sim.probe_subfield_siv('H_demag', [0,0,0])

print "H_demag_x at centre of cube = ", H_demag[0]
print "H_demag_y at centre of cube = ", H_demag[1]
print "H_demag_z at centre of cube = ", H_demag[2]

Setup can be splitted into three steps. In the first step we set the x_lattice parameter to be slightly larger than the dimension of the unit cell (in order not to have any overlap between the cells) and set the y_lattice and z_lattice parameters to zero to indicate no periodidicity along these directions

x_lattice = 15.01 # the spacing is 0.01
y_lattice = 0.0
z_lattice = 0.0

In the second step we define the lattice points where we want the periodic copies to be:

for xi in range(-3,4):
    lattice_points.append([xi*x_lattice,0.0*y_lattice,0.0*z_lattice])

and in the third step we define the object whose structure attribute will be used as the parameter in the definition of the simulation object

pbc = nmag.SetLatticePoints(vectorlist=lattice_points, scalefactor=SI(1e-9,'m'))

#create simulation object
sim = nmag.Simulation(periodic_bc=pbc.structure)

The remaining part of the script computes the demagnetisation field at the center of the cube. This calculation can be carried out for a varying number of copies of the simulation cell. The next figures show components of demagnetising field in the center of the cube as a function of the number of periodic copies. As in the code above, we impose an uniform magnetisation along the periodic x-axis. The first figure shows the demagnetisation field along the x-axis, and the second figure along the y-axis. In both figures, we have added green crosses that have been obtained by computing the demagfield using OOMMF (where in OOMMF we have actually made the simulation cell larger and larger to represent the growing number of periodic copies).

Demagnetising field as a function of the number of periodic copies with the magnetisation aligned along the periodic axis.

Demagnetising field as a function of the number of periodic copies with the magnetisation aligned along an axis orthogonal to the periodic one.

Relaxation script¶

To see how the system relaxes, we use the following script (spinwaves.py):

import nmag
from nmag import SI
import math

# define magnetic material
Py = nmag.MagMaterial(name="Py",
                      Ms=SI(1e6,"A/m"),
                      exchange_coupling=SI(13.0e-12, "J/m"),
                      llg_damping = SI(0.02,"")
                      )

# lattice spacings along the main axes; 
# the value must be zero for no periodic copies,
# equal to the mesh dimension along the 
# given axis otherwise
x_lattice = 30.0
y_lattice = 0.0
z_lattice = 0.0


# list to store the lattice points where the periodic
# copies will be placed 
lattice_points = []

for xi in range(-1,2):
    lattice_points.append([xi*x_lattice,0.0*y_lattice,0.0*z_lattice])

# copies of the system along the x-axis
pbc = nmag.SetLatticePoints(vectorlist=lattice_points, scalefactor=SI(1e-9,'m'))

#create simulation object
sim = nmag.Simulation(periodic_bc=pbc.structure)

# load mesh
sim.load_mesh("periodic.nmesh", [("periodic-film", Py)], unit_length=SI(1e-9,"m") )

print ocaml.mesh_plotinfo_periodic_points_indices( sim.mesh.raw_mesh )

# function to set the magnetisation 
def perturbed_magnetisation(pos):
    x,y,z = pos
    newx = x*1e9
    newy = y*1e9
    if 8<newx<14 and -3<newy<3:
        # the magnetisation is twisted a bit 
        return [1.0, 5.*(math.cos(math.pi*((newx-11)/6.)))**3 *\
                        (math.cos(math.pi*(newy/6.)))**3, 0.0]
    else:
        return [1.0, 0.0, 0.0]


# set initial magnetisation
sim.set_m(perturbed_magnetisation)

# let the system relax generating spin waves
s = SI("s")
from nsim.when import every, at
sim.relax(save=[('averages','fields', every('time', 0.05e-12*s) | at('convergence'))],
          do=[('exit', at('time', 10e-12*s))])

To execute this script, we call the nsim executable, for example (on linux):

$ nsim spinwaves.py

As in the previous examples, we first need to import the modules necessary for the simulation, define the material of the magnetic object, load the mesh and set the initial configuration of the magnetisation. Here, we start from a spatially non-homogeneous configuration in order to excite spin waves. |Nmag| allows us to provide a function to be sampled on the mesh that defines a particular magnetisation configuration.

In our case, we use the function

def perturbed_magnetisation(pos):
    x,y,z = pos
    newx = x*1e9
    newy = y*1e9
    if 8<newx<14 and -3<newy<3:
        # the magnetisation is twisted a bit
        return [1.0, 5.*(math.cos(math.pi*((newx-11)/6.)))**3 *\
                        (math.cos(math.pi*(newy/6.)))**3, 0.0]
    else:
        return [1.0, 0.0, 0.0]

which is then passed on to set_m

# set initial magnetisation
sim.set_m(perturbed_magnetisation)

Visualising the magnetisation evolution¶

Once the calculation has finished, we can see how the system relaxed by means of snapshots of the magnetisation evolution.

The nmagpp command allows us to create vtk files from the data saved with the save option in the relax method:

nmagpp --vtk=fields spinwaves

The first few frames that show the evolution of the magnetic configuration are shown below.

Initial magnetisation configuration.

Magnetisation configuration after 0.15 ps. It is clearly visible that the spin waves travel from the center of the disturbance to the right and penetrate the system immediately from the left (due to the periodic boundary conditions in the x-direction).

Magnetisation configuration after 0.25 ps.