MERRILL: Micromagnetic Earth Related Robust Interpreted Language Laboratory.

Complex magnetic domain structures and the energy barriers between them are responsible for pseudo-single-domain phenomena in rock magnetism and contribute significantly to the magnetic remanence of paleomagnetic samples. This article introduces MERRILL, an open source software package for three-dimensional micromagnetics optimized and designed for the calculation of such complex structures. MERRILL has a simple scripting user interface that requires little computational knowledge to use but provides research strength algorithms to model complex, inhomogeneous domain structures in magnetic materials. It uses a finite element/boundary element numerical method, optimally suited for calculating magnetization structures of local energy minima (LEM) in irregular grain geometries that are of interest to the rock and paleomagnetic community. MERRILL is able to simulate the magnetic characteristics of LEM states in both single grains, and small assemblies of interacting grains, including saddle-point paths between nearby LEMs. Here the numerical model is briefly described, and an overview of the scripting language and available commands is provided. The open source nature of the code encourages future development of the model by the scientific community.

Introduction

Paleomagnetic observations have contributed a wealth of information about the evolution of the Earth and other planetary bodies (Dunlop & Özdemir, 2001; Merrill, 1998), through the interpretation of remanent magnetization in naturally occurring magnetic minerals. The vast majority of the natural magnetic archives are recorded by nanosized particles whose magnetic properties initially could only be inferred from experimental observations, either on natural samples or on man‐made analogues of bulk particle arrays which have relatively broad particle size distributions. Rock magnetic interpretation is typically based on the assumption that paleomagnetic samples are dominated by single‐domain particles which are theoretically accessible by Néel's theoretical description (Néel, 1955).

With the advent of numerical micromagnetics and high‐resolution imaging techniques, it has been demonstrated that many of the remanence carriers in natural samples are inhomogeneously magnetized pseudo‐single‐domain (PSD) grains, which in their simplest form take on a single vortex (SV) magnetic structure. The first three‐dimensional micromagnetic models were introduced in the late 1980s (Schabes & Bertram, 1988; Williams & Dunlop, 1989), and over the last 30 years have significantly advanced our understanding of magnetic recording in interacting and noninteracting nonuniformly magnetized particles. In paleomagnetism, for example, such studies have demonstrated not only that small PSD grains, ubiquitous in rocks, primarily occupy SV states but that these states can provide reliable and stable recording of the ancient magnetic field that remain stable of over billions of years (Nagy et al., 2017).

Modern desktop computers and workstations are now powerful enough to perform micromagnetic simulations of a wide range of grain sizes, materials, and even clusters of grains that are of interest to rock‐magnetists, paleomagnetists, and environmental magnetists. The currently available open source programs (e.g. OOMMF (Donahue & Porter, 1999), Magpar (Scholz et al., 2003), and NMag (Fischbacher et al., 2007)) are developed for applications in material science and physics and require substantial specialist knowledge to install, maintain, and apply in Earth related contexts. The objective of this article is to introduce an open source micromagnetic model that we believe can significantly increase the ability of the rock and paleomagnetic community to define the magnetic properties of their samples through forward modeling of the behavior of their possible domain states.

We describe here a finite element micromagnetic modeling package called MERRILL (Micromagnetic Earth Related Robust Interpreted Language Laboratory) in honor of the early work of Ron T. Merrill on micromagnetics in rock magnetism (Merrill, 1977; Moon & Merrill, 1984, 1985). MERRILL is a script based modeling program designed for the Earth science community that needs no specialist computing knowledge or proprietary additional software to run the models and visualize the magnetic domain structures. Yet it is a fully tested research‐strength modeling platform that uniquely provides specific features relevant for natural samples and may even outperform some of the above mentioned software in certain applications. Both, precompiled binaries and the FORTRAN source code are freely available and run on LINUX, macOS, and Windows. They can be downloaded from the MERRILL homepage at http://www.rockmag.org.

Workflow

A typical workflow using MERRILL is outlined in Figure 1. A tetrahedral mesh of the geometry of interest must first be generated in an external program (although MERRILL contains routines for generating some commonly used geometries), representing the magnetic material. An “MScript” command file is then passed to MERRILL to drive the model, e.g., loading the mesh, setting material parameters, varying external fields, minimizing the micromagnetic energy, and outputting the magnetization to disk. When the magnetization has been output to disk, an external visualization program can be used to inspect the results. The visualizations of the geometries and magnetizations presented in this paper were generated with ParaView (Ahrens et al., 2005). This magnetization may also be used as a starting point in future models.

Figure 1

A typical workflow using MERRILL.

A productive workflow might involve running a micromagnetic model, inspecting the results visually, and then continuing the model using the previous results as the new starting point, until some desired result is achieved.

Micromagnetism

In micromagnetism, a real physical magnetic system is described in terms of a continuous vector valued function: where represents a mathematical magnetization vector at the mathematical point (Brown, 1963; Hubert & Schäfer, 1998). All physical energies and processes are then mathematically studied through this continuous function. To link the results to the physical situation, one can assume that the magnetization at a given mathematical point is constructed by averaging the discrete physical sources of magnetism, e.g., electron spins and orbitals, over a small volume centered at that point. This volume must be large enough that the behavior of individual atoms are averaged out, but small enough to resolve inhomogeneous magnetic structures such as domain walls, vortices, or flower states. It turns out that this continuum approximation represents real magnetization structures astonishingly well, even if the mathematical magnetization at a given point represents only the average over a few atoms. Quantum mechanical effects, which are essential in magnetism, are represented purely phenomenologically through material constants like the exchange constant or the magnetocrystalline anisotropy constant.

In micromagnetic models, the magnetic structure in a particle gives rise to various contributions to the total free magnetic energy functional , or to their associated effective fields. The effective field is given as a function of as

The dynamics of a micromagnetic system is described by the Landau‐Lifshitz‐Gilbert (LLG) equation (Gilbert, 2004): where γ is the electron gyromagnetic ratio and λ is a material‐dependent damping factor. This behaves like a dynamical system where a force acts on a system with a nonzero angular momentum vector pointing parallel to . A sufficient condition for an equilibrium magnetization is given by

This solution is local energy minimum (LEM), which we might associate with a remanent magnetization state, and so we will denote it . From equation (3) we can see two types of solutions. Either is parallel to or . In MERRILL, we focus on the second solution, that is, we solve for the LEM solutions by optimizing the expression for the total free magnetic energy of the system. This is generally much more efficient at finding stable domain structures that the full solution to the LLG equation, and we are not normally interested in the details of the domain transition dynamics that the LLG describes.

Effective Fields

The total effective field for a typical cubic ferromagnetic crystal involves four primary components: Zeeman, anisotropy, exchange and demagnetizing fields (Brown, 1963; Kittel, 1949).

Zeeman Field

The Zeeman field represents the interaction of external sources of magnetic field with the magnetic system under investigation. As such, it is often referred to as the “external field.” It is assumed that the system under investigation has no effect upon the external source.

Anisotropy Field

The anisotropy field couples the magnetization to the crystal lattice. It is the primary mechanism by which the symmetries of the lattice affect the magnetization. For a cubic ferromagnetic crystal, if the crystal axes are along the x, y, and z coordinate axes, it can be written with K 1 the anisotropy constant. In magnetite and iron higher‐order terms are much smaller and can be safely ignored for micromagnetic calculations. However, MERRILL does allow a value for K 2 to be set.

The vector here represents the directional cosines of the magnetization with respect to the crystal axes. For a crystal with cubic axes , and , the vector is defined: such that .

For other noncubic anisotropies, appropriate equations are used, e.g., for a crystal with a uniaxial symmetry, e.g., a tetragonal mineral like tetrataenite, a uniaxial anisotropy is commonly used, with a the anisotropy axis.

Exchange Field

The exchange field serves to align nearest neighbor magnetizations. Although this incorporates quantum mechanical spin coupling, in the continuum approximation a micromagnetic expression can be given as with A the exchange coupling constant.

Demagnetizing Field

The demagnetizing field represents the magnetic field generated by the magnetic material itself, derived from the “magnetic self‐energy.” It is called the demagnetizing field, because a higher magnetic self‐energy represents a higher energy configuration, so it typically acts on the magnetization in such a way that it minimizes external flux and thereby itself.

Outside the magnetic material, this is the “stray field.”

Minimum Energy Solutions Using Finite Elements

There are a variety of numerical micromagnetic approaches than can be used to solve for locally stable magnetic domain structures. Here we use a method that is both numerically efficient and robust while requiring the simplest input for representing the geometry of the grain. In most micromagnetic models, the primary consideration is the efficiency with which in the internal demagnetizing field can be computed. This field calculation scales as , where N is the number nodes, but if the geometry is meshed using a regular grid this scaling can be reduced to using FFT (Fabian et al., 1996; Wright et al., 1997). However, such regular grids do not easily account for arbitrary grain geometries, although some success has been achieved by relative scaling of surface elements (Witt et al., 2005), or irregular FFT techniques (Kritsikis et al., 2008).

In MERRILL, we employ the Finite Element Method (FEM), a standard technique for describing functions over a geometry and solving differential equations in terms of those functions (Davies, 2011). MERRILL uses arbitrarily shaped linear tetrahedral finite elements to describe the geometry of a particle and to solve for the LEM stable domain states. Some care is needed in the calculation of the demagnetizing field described by equations (8), (9), (10), which involves solving the Poisson equation for the magnetic scalar potential over an infinite space.

MERRILL makes use of a Boundary Element Method (BEM) technique (Fredkin & Koehler, 1990; Lindholm, 1984), which is a specialization of the FEM for homogeneous Poisson equations, particularly suited to problems defined over an infinite space. This method has significant advantage in that we need not create a mesh in the free‐space region outside the geometry of the magnetic particle (even for multiparticle solutions) but the method does, however, increase the memory requirements of the programme. In our experience, this is acceptable for single‐grain geometries requiring up to about a million elements. A more detailed account of different mciromagnetic methodologies can be found in Fidler and Schrefl (2000).

LEM states are found by solving for the minimum free magnetic energy of the system where the total effective field is given by

To include further phenomena, like stress fields or surface anisotropy, in theory one need only derive the corresponding effective field, typically by taking the derivative of the energy with respect to the magnetization, and add it to equation (11).

In order to determine minimum energy solutions, by default MERRILL makes use of an accelerated adaptive step‐size steepest descent algorithm across the energy landscape, optimized for micromagnetics, here called “Hubert Minimizer” (Berkov, 1998a; Ramstöck, 1997). A standard conjugate gradient optimizer is also available as an option, but in most cases is slightly less efficient in finding energy minima.

Mesh Generation

MERRILL requires that the geometry of the magnetic particle is described using a linear tetrahedral mesh (by default defined in micron units). The magnetization is thus specified only at the four vertices of each element and linearly interpolated at all other locations. MERRILL is able to generate suitable meshes for simple grain geometries such as cubes and spheres, however, more complex particle geometries require additional meshing software of which there are many free and commercial programmes. The only requirement is that the mesh is formatted according to ASCII text based PATRAN (.neu) standard. Most finite element meshing applications will support this and more detailed information about the PATRAN format can be obtained from the website of the MSC Software Corporation (https://simcompanion.mscsoftware.com).

Like all finite element models, the quality of the mesh will affect the convergence efficiency of the model. In most cases meshing software will take care to produce a qood quality mesh, but for highly irregular geometries problems may still arise. Discussion of mesh quality metrics can be found in many publications (e.g., Dai et al., 2014; Knupp, 2006)

For micromagnetic applications, it is important to ensure that that mesh is fine enough to resolve the expected spatial variation of the magnetization within the model geometry. The maximum element size is usually described in terms of the “exchange length,” (Rave et al., 1998), which is dependent on the magnetic material parameters:

For iron and magnetite, for example, the exchange length is around 3 and 9 nm, respectively, at 30°C and slightly larger near the Curie temperature.

Clearly, the bigger particle and finer the mesh the longer the model will take to converge to a LEM state. For large grains, the mesh size is usually set to the exchange length. However, small grain geometries will require a mesh that is finer than the exchange length in order adequately represent the grain shape.

Model Validation

Validation of micromagnetic models should ultimately be done against experimental observations. However, such direct validation has until recently been extremely difficult to achieve since the maximum grain size that could be modeled numerically was much smaller than that which could be directly observed experimentally. As a result, the earliest micromagnetic models could only be validated against bulk observations such as assemblies of sized particle fractions or magnetosome observations (Fabian et al., 1996; Williams & Dunlop, 1995; Witt et al., 2005) in the case of natural materials, or on thin film and particulate man‐made recording media (Labrune & Miltat, 1990; Silva & Bertram, 1990). However, results from MERRILL have been directly compared with nanoscale experimental data via electron holography with good agreement (Almeida et al., 2016).

With an increasing number of micromagnetic models being published, a number of standard tests were developed in order to ensure that the models were at least self‐consistent. Validation against these standard tests is now a prerequisite for any newly published micromagnetic code. One such test is μMAG Standard Problem 3 (see Appendix A), which tests for the critical edge length of a cube for transition between a flower state and a vortex state. Our solutions found the flower and vortex states had equal energies at an edge length of 8.47 when extrapolated to an infinitely fine mesh, which is in good agreement with other submissions to the μMAG problem.

MERRILL also tests the effective field components against some analytic solutions. For example, the demagnetizing field for a uniform sphere can be written and should be independent of the sphere size.

MScript: The MERRILL Scripting Language

MERRILL is run at the command line with an input text file containing a series of MScript commands for MERRILL to execute. These commands allow the user to interact with MERRILL in a number of ways, for example, setting material constants, loading meshes or finding a local energy minimum. In this section, we will outline a few typical computer experiments that can be used to probe the behavior of the magnetization of magnetic materials using MERRILL. This will also serve as a brief tutorial on the MScript language and introduce some of features in MERRILL.

The default units for MERRILL are microns and degrees Celsius, although different units may be used as specified in the user documentation. A full list of MERRILL commands can be found in Appendix B and on MERRILL download page (http://www.rockmag.org) where it will be updated as MERRILL is developed further.

Some typical commands include (but are not limited to)

Set material constants to magnetite at temperature degrees Celsius.

Set material constants to iron at temperature degrees Celsius.

Generate a cubic geometry of width width and with the average length between mesh nodes of length edgelen .

Set the magnetization of the material to

Run the minimizer to find the local energy minimum.

Write the magnetization to disk. Two files are written, filename .dat and filename _mult.tec. The file filename .dat contains a list of vertex points and the magnetization at that point. This is suitable for use with the ReadMagnetization command (see below). The filename _mult.tec file is the mesh and the magnetization in a TecPlot file format. This can be read by a number of visualization tools, including the free, open source and cross‐platform ParaView software (Ahrens et al., 2005).

Read a previously written file filename from the WriteMagnetization command, and set the current magnetization based on that.

Set the Zeeman field direction to

Set the Zeeman field strength to

MScript also includes supports for variables and loops. For example, the script

Print #myvalue

EndLoop

creates a loop with the local variable myvalue. The numbers in the loop command determine initial value, end value and step size, such that the above loop prints out the values 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, and 100. To distinguish between floating point, integer and string values, the variable myvalue can be referred to as %myvalue, #myvalue, and $myvalue$, depending on what is needed. For example, in

EndLoop

the filenames written to by the two WriteMagnetization commands are quite different. The argument SolnFileA_%myvalue produces filenames of the form SolnFileA_+0.10000E2 when myvalue is 10, while the argument SolnFileB_#myvalue produces filename of the form SolnFileB_10.

Minimization

A basic minimization from a random start, or a uniformly magnetized start can be a quick way to tell if a grain is in the SD or PSD size range. Moreover, it is the basic building block from which most other MERRILL modeling experiments come.

A simple minimization script to determine a LEM state from a random initial state for a meshed geometry can be given as follows:

The output of this script can be seen in Figure 2 using a mesh generated from experimental focused ion beam nanotomography data for an example iron inclusion in dusty olivine in a chondritic meteor sample (Einsle et al., 2016).

Figure 2

Minimization of an iron inclusion from a dusty olivine sample at 20 C with random initial magnetization. The shape is approximately a prolate ellipsoid with semimajor axes of 0.134 and 0.126 μm, and semiminor axis of 0.085 μm, with axes aligned along the cardinal directions. The LEM is a single vortex state with core aligned along the semiminor axis.

It is important to note that there may be several possible magnetization states the minimization might settle upon, representing multiple possible, and completely valid local energy minima states. However, from a given starting point, barring numerical noise, the minimization should always reach the same end point.

The LEM found can be highly dependent on the choice of initial state. A simple minimization script to determine a LEM state for a cubic grain of magnetite, starting with uniform magnetization, can be given as follows:

! easy axis

The output of this script can be seen in Figure 3, representing a flower state.

Figure 3

Minimization of a 0.08 μm cube of Magnetite to a flower state from an initial uniform [111] magnetization.

By starting from a rough vortex state rather than a uniform state, we can instead nucleate a vortex LEM. This can be accomplished by

.

The output of this new script (seen in Figure 4) highlights the “local” aspect of a local energy minimum.

Figure 4

Minimization of a 0.08 μm cube of Magnetite to a vortex state from an initial approximate [111] aligned vortex state, colored by helicity ( ).

Hysteresis Loops

Hysteresis loops are a useful tool when a magnetic material has several remanent magnetization states for a given set of parameters. From a hysteresis loop, it is possible to deterministically move from one state to another, and find the tipping point where the variation of a given parameter will cause one state to spontaneously switch to the other. This also provides information on the range of values for the given parameter where both states can coexist. In other words, it provides some information about the stability of a given state when multiple valid states exist.

Since the minimization should always reach the same solution when given the same starting point, a hysteresis loop run with the same parameters and the same changes in parameters should always return the same results in MERRILL.

Magnetic Field Hysteresis

A magnetic hysteresis loop has several uses. A single hysteresis loop can be useful to get a feel for the behavior of the magnetization of a system (i.e., its coercivity, if it is SD or PSD). The average of hysteresis loops in many directions can be used to compare simulations with experimental observations of magnetic characteristics such as coercivity.

A hysteresis loop is a quasi‐static thermodynamic process. That is to say, at each point of the hysteresis loop, the system is assumed to be in equilibrium and the energy at a local energy minimum. A change from one point to the next in a hysteresis loop represents the system moving from what was a local energy minimum in the previous step to the nearest, new local energy minimum in the current step.

The simulation of a hysteresis loop is accomplished by the following scheme:

Set Zeeman field to saturating value in a fixed direction: Z = Zmax.

Find LEM.

Update Zeeman field: Z = Z − a*Zmax.

If Z != −Zmax, go to 2.

For a saturating field , a hysteresis loop will run from to − in small increments. Here “small” means sufficiently small enough that any phenomenon sensitive to the external field, e.g., magnetization switching or nucleation of vortices, are accounted for. In practice, “small” means just small enough to resolve changes in the magnetization you are looking for, but as large as possible to reduce the total number of steps needed for the loop. Otherwise a large amount of work may be done unnecessarily for little progress.

A MERRILL script which accomplishes this is

Since hysteresis is generally symmetrical to forward and reverse fields, it is only necessary to run either the upper or lower branch of the hysteresis loop, not both.

A plot of the upper branch is shown in Figure 5. This curve is different to the usual SD curve, where the value remains near the saturation value, until passing the coercive field. Here the grain enters a vortex state, and the value can minimized while keeping the anisotropy energy relatively low by varying the shape of the vortex. An example of the vortex state at mT is shown in Figure 6.

Figure 5

The lower branch of a hysteresis loop with respect to an external magnetic field, varying from −100 to 100 mT for a 0.1 μm sphere of magnetite. The arrow denotes the direction of change of during the loop.

Figure 6

Magnetization of a 0.1 μm sphere of magnetite at mT during a field hysteresis from mT to H = 100 mT, colored by helicity ( ).

Size Hysteresis

A hysteresis loop where the size of the grain is varied rather than an external field can be a useful tool for determining SD and PSD ranges, and the evolution of magnetic domain states. Since several remanent states can exist for a grain, particularly in the early PSD size range, it can be difficult to pinpoint exactly where that regime begins.

In a size hysteresis loop, the grain is started in an SD magnetization state. The size of the grain is increased until it spontaneously switches to a vortex PSD state. The size of the grain is then scaled down until it is in a SD state again. The branch of increasing size can tell us what the largest grain size is that supports an SD state, and the branch of decreasing size can tell us what the smallest grain size is that supports a PSD state.

Note that the exact initial SD was not defined. Ideally, a size hysteresis should be run for each possible SD remanent state. For symmetric grains, symmetry should make many of these redundant. However, for asymmetric grains, the direction of the initial magnetization may greatly affect the stability of the SD solution.

When scaling sizes, it is important to consider the size of the elements of the mesh during scaling. Typically, the edges of a mesh should be around the exchange length. When the size of a mesh is increased, the average edge length should not exceed the exchange length. By starting with a larger mesh and scaling down, rather than the other way around, we can avoid this problem. However, finer meshes take longer to run. A compromise can be found using the “Remesh” command. A user might use a mesh that is suitable for scaling up to, e.g., 0.1 μm, and at that size, switch to a mesh suitable for up to, e.g., 0.2 μm.

So if we want a maximum node spacing of, e.g., 0.005 μm, we make two meshes of width 0.1 μm: one with node spacing 0.005 μm, and one with node spacing of 0.0025 μm. The first mesh will cover grains from 0 to 0.1 μm, while the second can cover grains from 0.1 to 0.2 μm. We make the meshes the same initial size so the Remesh command can interpolate the magnetization directly from one to the other.

In MERRILL, the “Resize” command can be used to scale the mesh. An example size hysteresis script incorporating all of this is

! Run the minimization

EndLoop

A graph of the saturation magnetization versus the grain size is shown in Figure 7, with three distinct magnetic phases marked: a flower state (FS), a hard‐aligned single vortex (HSV), and an easy‐aligned easy vortex state (ESV). A comparison of the upper and lower branches at 0.13 μm and at 0.15 μm, are shown in Figure 8, demonstrating the ESV and HSV states for the lower branch, and the corresponding FS states on the upper branch.

Figure 7

A full hysteresis loop of saturation magnetization versus size for an octahedral grain of magnetite. Three distinct magnetic phases are marked: a flower state (FS), a hard‐aligned single vortex state (HSV), and an easy‐aligned single vortex state (ESV).

Figure 8

Comparison of upper and lower branch results for a size hysteresis of an octahedral grain of magnetite. The upper branch goes from small to large, and the lower branch goes from large to small. (a) The upper branch at 0.13 μm, (b) the lower branch at 0.13 μm, (c) the upper branch at 0.15 μm, and (d) the lower branch at 0.15 μm. The coloring is helicity ( ).

Temperature Hysteresis

A temperature hysteresis can be used to observe the effects of heating, say iron, from room temperature to the Curie temperature and back. The approach MERRILL takes to this is to simply change the temperature dependant material parameters for each step. In this way, the model is of a “cold” material, since no other thermal effects are taken into account which might spontaneously and nondeterministically effect the magnetization. However, the changes in material properties can have significant impact on the shape of the energy landscape with respect to the depths, positions and even number of the local energy minima.

An example temperature hysteresis of iron is

.

Energy Barriers

MERRILL includes a new method for finding the minimum energy transition between two given LEM states (K. Fabian & V. P. Shcherbakov, Energy barriers in three‐dimensional micromagnetic models and the physics of thermo‐viscous magnetization in multidomain particles, arXiv 1702.00070v1, 2017). It uses a combination of the Nudged Elastic Band technique (Dittrich et al., 2002; Henkelman & Jónsson, 2000) and an action minimization method (Berkov, 1998a, 1998b; Fabian & Shcherbakov, arXiv 1702.00070v1, 2017). It finally constructs a prescribed number of intermediary states between the given LEM states L 1, L 2, such that their interpolation represents the average physical switching path for a fully damped (λ = 0) switching process from L 1 to L 2. The path with the lowest maximum energy barrier determines the thermal activation energy required to switch between these states. From that, the probability of switching and the relaxation time across the energy barrier can be obtained.

These quantities are of paramount importance in paleomagnetism, where small changes in size, shape, or material parameters can change the relaxation time from the order of seconds, to millions of years, and the corresponding paleomagnetic information is either completely unblocked, or can be regarded as stable or blocked. In contrast to, for example, the simulation of writing heads for magnetic hard drives, the nanosecond dynamics of the switching processes are of no interest in paleomagnetism, such that time‐consuming full‐fledged LLG models do not improve the result. On the contrary, the exact determination of the energy barrier through the saddle‐point path provides a better estimate of the thermal relaxation time. Note in that respect that the LLG equation does not include temperature, and that thermal activation has to be added in a phenomenological way, e.g., as a stochastic zero‐mean Gaussian fluctuation field (Torres et al., 2001), if it is considered at all.

If morphology and composition of a natural magnetic mineral system is known MERRILL can provide quantitative insight to whether the signal recorded has remained stable since it was recorded (Nagy et al., 2017).

An example script for a NEB calculation is

! predefined material (iron)

! point 1

.

Practical Considerations

There are a number of practical considerations when using MERRILL. Presently, MERRILL is not parallelized. However, many instances of MERRILL may be run in task‐farm parallel manner for reasonably small grains. This approach can be used to model, for example, large assemblies of noninteracting grains, and hysteresis loops about many axes in parallel. The maximum grain size that you will be able to model will be dependent on a number of factors, but primarily the number of elements in the model. MERRILL will comfortably cope with models of up to approximately 106 elements for single grain of simple geometry on a machine with access to 64 GB of RAM (see Figure 9). If the model consists of a cluster particles where the total grain surface to volume ratio is significantly greater than that of a single grain then the maximum number of elements in the model would need to be less than 106. It is not possible to formulate any hard rules that cover the combination of different number of particles, grain sizes and mesh sizes, but the user should be aware that the memory requirement increases as the square of the number of surface nodes (which is typically proportional to the number of surface elements).

Figure 9

A log‐log plot of the time taken for MERRILL to find a local energy minimum of a cube, starting from a random state, versus the number of elements in the mesh.

The time taken for a model to converge on a local energy minimum will again be largely dependent on the model size (see Figure 9), but also on how close your initial guess is to the final LEM solution. Thus solution times for a random initial guess and a uniform saturated state initial guess can have significantly different solution times. A similar situation occurs during a change in the magnetic phase, i.e., when a solution changes from SD to SV, or hard‐aligned SV to easy‐aligned SV.

In rare cases, the model might attempt a large number of iterations because of slow convergence to the final solution. In such cases it is better to place a limit in on the maximum number of allowed iterations (which has default value of 100,000). Often such slowly converging solutions can be avoided by trying a slightly different initial guess, slightly changing the model mesh, or adding a small random kick to the solution.

As mentioned in section 4, the mesh shape and quality also impact on convergence times. Regular tetrahedra are the most numerically stable elements, so the more regular the tetrahedra in a mesh, the quicker the model will converge. The single most important concern, however, is to avoid slivers. Slivers are elements that are so flat, their volumes approach zero (Alliez et al., 2017). These lead to extreme numerical instability and very slow convergence times, if the model converges at all. Slivers might occur during Delaunay triangulation, as the algorithm is well known to be blind to a class of slivers where the vertices all lie on the same circumcircle. An additional mesh optimization step like Lloyd optimization, or ODT smoothing is often performed during mesh generation to eliminate slivers. Most finite element meshing packages will have algorithms that attempt to detect and avoid such elements, but in all cases when the model fails to converge it is important to check the mesh quality (e.g., using the “Mesh Quality” filter in ParaView).

Also mentioned in section 4 is that mesh edge lengths should be smaller than the exchange length. In most meshing programs, it is unclear what the defined “mesh size” is. In some cases, it is the upper bound of the radius of the circumcircle of an element allowed before the element is decomposed into smaller elements. In some cases, it is the diameter of this circumcircle. For some software, it is the target average edge length of the elements in the mesh. In any case, the ReportEnergy function of MERRILL reports the average edge length of the mesh in nanometres, which is then directly comparable to the exchange length reported by this function in nanometres. Users are encouraged to use these values as a useful and direct report of whether the mesh is fine enough for the given material.

The solvers used in MERRILL can be susceptible to finding unstable equilibria. Gimbal locking, for instance, can occur for field hysteresis loops directly along the hard‐axis of the material. Two simple techniques can be used to avoid this. First, prefer not to perform hysteresis loops directly along the major axes. Add a slight perturbation, i.e., use the [1 1 1.0001] direction instead of [1 1 1]. Another technique is to add a small random “kick” to the magnetization. Typically, moving each magnetization randomly by around 5 or 10 degrees is not so large that it will move the state out of a deep, stable LEM, but it can move it off a local maximum, and also out of a shallow LEM, where the magnetization is around a SD/SV phase transition.

Discussion

The formulation of a micromagnetic model in an easy to use form can significantly enhance its application in paleomagnetic and rock magnetic investigations. MERRILL presents such a tool that is particularly focused on finding remanent states and studying their stability. The parallels and usefulness to paleomagnetic and rock magnetic studies should be clear. MERRILL has been used in a number of publications and talks using functionality not presented in this paper, e.g., behaviors of assemblages of interacting grains, large models, strongly anisotropic materials and detailed simulations of magnetostrictive effects (e.g., Almeida et al., 2016; Chang et al., 2012; Conbhuí et al., 2016; Einsle et al., 2016; Li et al., 2013; Nagy et al., 2017; Williams et al., 2010). Not all of these, however, used the scripting interface presented here.

The simple scripting language makes it particularly friendly to nontechnical users. The fast and efficient minimization scheme means simple computer experiments can be run quite quickly. The scripts presented here represent the sort of scripts run by the authors day to day. Some effort has been made to make these copy/pasteable, but a curious reader is recommended to look in the “demo” directory of the MERRILL package for more information and examples.

If the scripting language is not up to a particular task, MERRILL can also be used as a library and called from a Fortran program. In addition, a plugin interface has been included so that users can compile libraries that can be loaded from the MScript interface and hook into the MScript parser to add commands and variables, and also add effective field calculators to not originally shipped with MERRILL. This should allow MERRILL to be adopted to a wide range of problems.

For viewing solutions, the authors recommend ParaView. It is a free, open source 3D viewer, easily downloaded and installed from the ParaView website, which can open and visualize solutions generated by MERRILL from the WriteMagnetization command. A plugin for ParaView for opening MERRILL solutions with some preprocessing already done can be found in the demo directory.

Table A1

Extrapolated Partial Energies at the Critical Edge Length in Units of

	Flower	Vortex
Demag	0.2791 ± 0.0001	0.0775 ± 0.0006
Anisotropy	0.00563 ± 0.00002	0.0520 ± 0.0001
Exchange	0.01773 ± 0.00005	0.1729 ± 0.0007

Table A2

Extrapolated Magnetizations at the Critical Edge Length in Units of M

	Flower	Vortex
〈m_x〉	0.97086 ± 0.00008	−0.0000 ± 0.0008
〈m_y〉	−0.0000 ± 0.0002	0.0 ± 0.3
〈m_z〉	−0.0003 ± 0.0002	0.0000 ± 0.0007