Nudged Elastic Band Method (NEBM)

The NEBM is an algorithm to find minimum energy transitions between equilibrium states. This method was developed by Henkelman and Jónsson [1] to solve molecular problems in Chemistry. Later, Dittrich et al. [2] applied the NEBM to magnetic systems based on the micromagnetic theory. Recently, Bessarab et al. [3] have proposed an optimised version of the NEBM with improved control over the behaviour of the algorithm. Their review of the method allows to apply the technique to systems described by either a micromagnetic or atomistic model.

The algorithm is based on, firstly, generating \(N\) copies of a magnetic system that we describe by \(P\) spins or magnetic moments arranged in a lattice or mesh. Every copy of the system, is called an image \(\mathbf{Y}_{i}\), \(i\in \{0,\ldots,P-1\}\) and we specify an image using the spins directions in a given coordinate system. For example, in Cartesian coordinates we have

\[\mathbf{Y}_{i} = \left( s_{x,0}^{(i)}, s_{y,0}^{(i)}, s_{z,0}^{(i)}, s_{x,1}^{(i)}, s_{y,1}^{(i)},\ldots, s_{y,P-1}^{(i)}, s_{z,P-1}^{(i)} \right)\]

for a system described by a discrete spin model. Within the micromagnetic model we use \(\mathbf{m}\) rather than \(\mathbf{s}\), but both are unit directions.

This sequence of images defines a band, where every image is in a (ideally) different magnetic configuration. The energy of the system depends on the magnetic configuration (e.g. a skyrmion, vortex, uniform state,etc.), thus the energy is parametrised by the number of degrees of freedom, which is \(n\times P\), where \(n\) is the number of coordinates to describe a spin (e.g. \(n=3\) for Cartesian coordinates and \(n=2\) in Spherical coordinates). Accordingly, every image will have a specific energy, i.e. \(E=E(\mathbf{Y})\), which determines the position of an image in an energy landscape. The first and last images in a band are chosen as equilibrium states of the magnetic system and they are kept fixed during the NEBM evolution.

Secondly, to initiate the NEBM evolution, it is necessary to specify an initial guess for the images, which means generating different magnetic configurations between the extrema of the band. It is customary to use an interpolation of the spins directions to achieve this.

And thirdly, the band is relaxed to find a path in energy space that costs less energy. This path is characterised by passing through a saddle point, which is a maximum in energy along certain directions in phase space, and this point determines an energy barrier between the two equilibrium states. The energy barrier is the energy necessary to drive one equilibrium state towards the other. A first order saddle point is the one that is a maximum along a single direction in phase space and will usually in the energy path that costs less energy. It is possible that there are more than one energy paths.

It is also necessary to consider that:

• To distinguish different images in the energy landscape we need to define a distance

• To keep the images equally spaced to avoid clustering around saddle points or equilibrium states, we use a spring force between images.

• Spins or magnetic moments can be described in any coordinate system. The most commonly used are spherical and Cartesian coordinates. For the later we have to specify the constraint to fix the spin length.

For a thorough explanation of the method see references [3] and [4].

NEBM relaxation

Cartesian

When using Cartesian coordinates, every image of the band \(\mathbf{Y}_{i}\) is iterated with the following dynamical equation with a fictitious time \(\tau\)

\[\frac{\partial \mathbf{Y}_i}{\partial \tau} = -\gamma \mathbf{Y}_{i} \times \mathbf{Y}_{i} \times \mathbf{G}_{i} + c \sqrt{ \left( \frac{\partial \mathbf{Y}_{i}}{\partial \tau} \right)^{2} } \left( 1 - \mathbf{Y}_{i}^{2} \right) \mathbf{Y}_{i}\]

In this equation, \(\gamma\) is in units of \(\text{Hz T}^{-1}\) and it determines the time scale, which is irrelevant here so we set \(\gamma=1\). The second term to the right is necessary to keep the spins/magnetisation length equal to one, using an appropriate factor \(c\) that we set to 6. The \(\mathbf{G}\) is the total force on an image, which in the atomistic theory is defined as

\[\mathbf{G}_{i} = - \boldsymbol{\nabla}_{\boldsymbol{\mu}} E(\mathbf{Y}_{i})|_{\perp} + \mathbf{F}(\mathbf{Y}_{i})|_{\parallel}\]

and in the micromagnetic theory as

\[\mathbf{G}_{i} = - \boldsymbol{\nabla}_{\mathbf{M}} E(\mathbf{Y}_{i})|_{\perp} + \mathbf{F}(\mathbf{Y}_{i})|_{\parallel}\]

We use the magnetic effective field definition to evaluate the gradients, i.e. \(\boldsymbol{\nabla}_{\boldsymbol{\mu}}E=\partial E / (\mu_{s}\partial\mathbf{s})=-\mathbf{H}_{\text{eff}}\) or \(\boldsymbol{\nabla}_{\mathbf{M}}E=\partial E / (M_{s}\partial\mathbf{m})=-\mathbf{H}_{\text{eff}}\). The perpendicular component is with respect to the tangents \(\mathbf{t}\) to the energy band, thus for a vector \(\mathbf{A}\)

\[\mathbf{A}|_{\perp} = \mathbf{A} - (\mathbf{A}\cdot\mathbf{t})\mathbf{t}\]

The tangents are defined according to the energies of the neighbouring images [3]. The second term to the right hand side of the equation for \(\mathbf{G}\) is the spring force that tries to keep images at equal distance and is defined using the distance between neighbouring images

\[\mathbf{F}(\mathbf{Y}_{i})|_{\parallel}=k\left(|\mathbf{Y}_{i+1}-\mathbf{Y}_{i}|- |\mathbf{Y}_{i}-\mathbf{Y}_{i-1}|\right)\mathbf{t}_{i}\]

which is parallel to the band, i.e. in the direction of the tangent.

According to Bessarab et al. [3], the tangents and the total force \(\mathbf{G}\) must be projected into the spin/magnetisation tangent space.

Climbing Image NEBM

The climbing image technique is a modification of the NEBM where the forces of the image with largest energy in the band, are redefined so this image can climb up in energy along the band to get a better estimate of the saddle point energy [1] and [3]. The image with largest energy is chosen after relaxing the band with the usual NEBM algorithm. The total force on this climbing image is

\[\mathbf{G}^{\text{CI}}_{i} = - \boldsymbol{\nabla}_{\boldsymbol{\mu}} E(\mathbf{Y}_{i})|_{\perp} + \boldsymbol{\nabla}_{\boldsymbol{\mu}} E(\mathbf{Y}_{i})|_{\parallel}\]

where the spring force was removed. The climbing image is still allowed to climb down in energy in a direction perpendicular to the band, thus it is possible that the energy barrier magnitude decreases after applying this technique.

Vectors

Following the definition of an image in Cartesian coordinates, we mentioned that the number of degrees of freedom is \(n\times P\), where \(n\) is the number of coordinates to describe a spin. Accordingly, many of the vectors in the NEBM algorithm such as the tangents, total forces, etc. have the same number of components, which agree with the spin components of an image.

For instance, the total force (or tangents, spring forces, etc.) has three components in Cartesian coordinates, corresponding to every spin direction:

\[\mathbf{G}_{i} = \left( G_{x,0}^{(i)}, G_{y,0}^{(i)}, G_{z,0}^{(i)}, G_{x,1}^{(i)}, G_{y,1}^{(i)},\ldots, G_{y,P-1}^{(i)}, G_{z,P-1}^{(i)} \right)\]

Projections

The projection of a vector into the spin/magnetisation tangent space simply means projecting its components with the corresponding spin/magnetisation field components. For example, for a vector \(\mathbf{A}\) associated to the \(i\) image of the band (we will omit the \((i)\) superscripts in the spin directions \(\mathbf{s}\) and the \(\mathbf{A}\) vector components)

\[\mathbf{A} = \left( \mathbf{A}_{0}, \ldots \mathbf{A}_{P-1}\right) = \left( A_{x,0}, A_{y,0}, A_{z,0}, A_{x,1}, A_{y,1},\ldots, A_{y,P-1}, A_{z,P-1} \right)\]

the projection \(\mathcal{P}\) is defined as

\[\mathcal{P}\mathbf{A} = \left( \mathcal{P}_{\mathbf{s}_{0}}\mathbf{A}_{0}, \mathcal{P}_{\mathbf{s}_{1}}\mathbf{A}_{1}, \ldots \mathcal{P}_{\mathbf{s}_{1}}\mathbf{A}_{P-1}, \right)\]

where

\[\mathcal{P}_{\mathbf{s}_{j}}\mathbf{A}_{j} = \mathbf{A}_{j} - \left( \mathbf{A}_{j} \cdot \mathbf{s}_{j} \right) \mathbf{s}_{j}\]

with \(j\in\{0,\ldots,P-1 \}\), hence

\[\mathbf{A} = \left( \mathcal{P}A_{x,0}, \mathcal{P}A_{y,0}, \ldots, \mathcal{P}A_{y,P-1}, \mathcal{P}A_{z,P-1} \right)\]

Distances

There are different ways of defining the distance in phase space between two images, \(d_{j,k}=|\mathbf{Y}_{j} - \mathbf{Y}_{k}|\).

Geodesic

The optimised version of the NEBM [3] proposes a Geodesic distance based on Vicenty’s formulae:

\[d_{j,k} = \sqrt{ \left( \delta_{0}^{(j,k)} \right)^{2} + \left( \delta_{1}^{(j,k)} \right)^{2} + \ldots \left( \delta_{P-1}^{(j,k)} \right)^{2} }\]

where

\[\delta_{i}^{(j,k)} = \arctan2 \left( \left| \mathbf{m}_{i}^{(j)}\times \mathbf{m}_{i}^{(k)} \right|, \mathbf{m}_{i}^{(j)}\cdot \mathbf{m}_{i}^{(k)} \right)\]

This definition seems to work better with the NEBM since the spin directions are defined in a unit sphere.

Euclidean

The first versions of the method simply used an Euclidean distance based on the difference between corresponding spins. In Cartesian coordinates it reads

\[d_{j,k} = \frac{1}{3 P} \left\{ \sum_{j=0}^{P-1} \sum_{\alpha\in\{x,y,z\}} \left[ \left( s_{\alpha}^{(j)} - s_{\alpha}^{(k)} \right)^{2} \right] \right\}^{1/2}\]

where we have scaled the distance by the number of degrees of freedom of the system (or an image). In spherical coordinates we use the difference of the polar and azimuthal angles, but without the scaling, since the tangents are normalised in the same metric (see Spherical coordinates code).

Algorithm

The algorithm can be summarised as:

1. Define a magnetic system and find two equilibrium states for which we want to find a minimum energy transition.

2. Set up a band of images and an initial sequence between the extrema. We can use linear interpolations on the spherical angles that define the spin directions [4] or Rodrigues formulae [3].

3. Evolve the system using the NEBM dynamical equation, which depends on the chosen coordinate system. This equation involves:

   1. Compute the effective field for every image (they are in different magnetic configurations) and the total energy of every image

   2. Compute the tangents according to the energies of the images and project them into the spin/magnetisation tangent space

   3. Compute the total force for every image in the band using the tangents and distances between neighbouring images. This allows to calculate the gradient (which uses the effective field) and the spring forces on the images

   4. Project the total force into the spin/magnetisation tangent space

   5. Use the dynamical equation according to the coordinate system

Early versions of the NEBM did not project the vectors into the tangent space in steps I and II. This leads to an uncontrolled/poor behaviour of the band evolution since the vectors that are supposed to be perpendicular to the band still have a component along the band and interfere with the images movement in phase space.

Fidimag Code

The NEBM is not the only method that relaxes a chain of images, so the code is organised around what these methods have in common rather than around the NEBM alone. A ChainMethodBase class holds everything that does not depend on the coordinate system or on the method: the band and the arrays that go with it, the energies, the interpolations that generate the initial band, the integrators, the data saving and the relax loop. The methods themselves derive from it and supply the parts that do differ, which are the distance between two images, the tangents, the spring force and the step taken at every iteration.

Three classes derive from it:

1. NEBM_Geodesic: Using Cartesian coordinates for the spin directions and Geodesic distances, with vectors projected in tangent space. This is the optimised version of the NEBM, for which [3] is the main reference, and it appears to work well with every system we have tried so far. Cartesian coordinates have the advantage that they are well defined close to the poles of the spin directions. Fabian and Shcherbakov [6] have also applied geodesic distances, in a minimum action formulation. Prefer this class for most calculations.

2. NEBM_Spherical: Using spherical coordinates for the spin directions and Euclidean distances with no projections into spin space. The azimuthal and polar angles need to be redefined when performing differences or computing Euclidean distances, specially because the polar angle gets undefined when it is close to the north or south. This class finds the same energy barriers as NEBM_Geodesic, but it is slower when the spins pass close to the poles, see Spherical coordinates code.

3. StringMethod: Not a NEBM class. Instead of a spring force, the images are redistributed along the path by interpolating it after every step, which removes the spring constant from the problem altogether. See [5] for a description of both methods side by side.

An earlier NEBM_Cartesian class, using Cartesian coordinates with Euclidean distances and no projections, has been removed. It worked for simple systems, but as soon as the system was complex enough to stabilise a vortex or a skyrmion the spring force interfered with the convergence of the band, and avoiding that required an optimal spring constant that depends on the system size and on the interactions involved. NEBM_Geodesic supersedes it. Some scripts under examples/ still import it and have not been updated.

The following diagram shows how the code is structured:

There is a fidimag.common.chain_method_tools module (formerly nebm_tools) with common functions for these classes:

fidimag.common.chain_method_tools
|
--> cartesian2spherical
    spherical2cartesian
    compute_norm
    interpolation_Rodrigues_rotation
    linear_interpolation_spherical
    m_to_zero_nomaterial

and a fidimag.common.chain_method_integrators module with the integrators that do not come from Sundials:

fidimag.common.chain_method_integrators
|
--> StepIntegrator      # Euler and Runge-Kutta 4 steps
    VerletIntegrator    # Velocity projection, see below

Arrays

In fidimag we mainly use Numpy to define the NEBM vectors. When calling one of the NEBM classes, we have to pass a sim object with the specification of the magnetic system which has associated a mesh with n_spins. According to the coordinate system, we set the dof variable. For instance dof = 3 for Cartesian coordinates. Consequently, we define the number of degrees of freedom per image n_dofs_image = dof * n_spins, thus if the NEBM class was specified with n_images, the total number of degrees of freedom for the band is n_band = n_dofs_image * n_images.

As explained in our discussion about the NEBM, we set up band, negH, tangents and spring_force arrays whose length is n_band. The order is the same than how we defined the images, thus the Numpy array, when using Cartesian coordinates to describe the spins, looks like

band = [ s(0)_{x,0}  s(0)_{y,0}  ... s(0)_{z,n_dofs_image-1}
         s(1)_{x,0}  s(1)_{y,0}  ... s(1)_{z,n_dofs_image-1}
         ...
         s(n_images-1)_{x,0}  ... s(n_images-1)_{z,n_dofs_image-1}
       ]

and similarly for the other vectors since they follow the same order of the spins. This band array is passed to the Cython codes to compute the NEBM forces. Notice that we can easily redefine this array into a (n_images, n_dofs_image) shaped Numpy array using

band = band.reshape(-1, n_dofs_image)

so every row is a different image. We can even take the inner images (no extrema) and use the same piece of code since n_dofs_image does not change.

Cython Codes

Most of the calculations are made using C code through Cython. The files for these libraries are located in fidimag/common/neb_method/. Every library has a .c file and a .h header file, and they are built through CMake, which is driven by pyproject.toml (the old setup.py build is gone, see BUILD.md).

There used to be one Cython module per coordinate system. They are now a single nebm_clib.pyx, exposed as fidimag.extensions.nebm_clib, which wraps all of the C libraries below. The base one, with the functions that do not depend on the coordinate system, is nebm_lib.c

nebm_lib.c
nebm_lib.h
|
--> compute_tangents_C
    compute_spring_force_C
    compute_effective_force_C
    compute_dYdt_C
    compute_dYdt_nc_C
    compute_distance_cartesian
    compute_image_distances
    project_images_C
    project_vector_C
    normalise
    normalise_images_C
    normalise_spins_C
    compute_norm
    cross_product
    dot_product

and the two coordinate-system-dependent ones only have to define what a distance is:

nebm_geodesic_lib.c
nebm_geodesic_lib.h
|
--> compute_geodesic_GreatCircle
    compute_geodesic_Vincenty


nebm_spherical_lib.c
nebm_spherical_lib.h
|
--> compute_distance_spherical
    normalise_spherical
    normalise_images_spherical_C

There is also nebm_integrators.c, with step_Verlet_C, which is the velocity projection integrator described below.

Note that the Cartesian library was removed together with the NEBM_Cartesian class. The spring force is still computed by the single compute_spring_force_C of nebm_lib.c, which takes the distances as an argument, so it is the distance definition, and not the spring force, that changes with the coordinate system.

The functions that fidimag.extensions.nebm_clib exposes to the Python classes are

fidimag.extensions.nebm_clib
|
--> compute_tangents
    compute_spring_force
    compute_effective_force
    compute_dYdt
    compute_dYdt_nc
    project_images
    project_vector
    normalise_clib
    normalise_images
    normalise_spins
    geodesic_distance_GreatCircle
    geodesic_distance_Vincenty
    image_distances_GreatCircle
    image_distances_Spherical
    normalise_images_spherical

Geodesic distances code

Based on the aforementioned NEBM algorithm, the class initialise the NEBM calling the following methods in this order:

1. generate_initial_band   # Using linear interpolations or Rodrigues rotation formulae
|
--> chain_method_tools.cartesian2spherical
    --> chain_method_tools.linear_interpolation_spherical
    # or
    chain_method_tools.interpolation_Rodrigues_rotation

2. initialise_energies     # Fill the energies array
3. initialise_integrator   # Start the integrator
4. create_tablewriter      # To pass data into .ndt files per every iteration of the integrator

The linear interpolation function requires that the input array is in spherical coordinates.

The integrator of step 3. is chosen with the integrator argument of initialise_integrator. The default is 'sundials', i.e. the cvode.CvodeSolver (or cvode.CvodeSolver_OpenMP) integrator, and the alternatives are 'euler' and 'rk4', which use the StepIntegrator class, and 'verlet', which uses the VerletIntegrator velocity projection scheme. All of them need a function giving the right hand side of the \(\partial \mathbf{Y} / \partial \tau\) dynamical equation, which is Sundials_RHS for CVODE and step_RHS for the others; both of them call nebm_step, where the NEBM forces are calculated as

Sundials_RHS  /  step_RHS
|
--> nebm_step
    |
    --> compute_effective_field_and_energy  # Gradient = - Eff field
                                            # Which we compute for every image
                                            # using the sim class
        nebm_clib.project_images(negH)      # Project the gradient
        compute_tangents
        |
        --> nebm_clib.compute_tangents      #
            nebm_clib.project_images        # Project tangents
            nebm_clib.normalise_images      #
        compute_spring_force                # Using Geodesic distances
        nebm_clib.compute_effective_force

    G[extrema] = 0                          # The force on the extreme images
                                            # should already be zero

Notice that it is the energy gradient that is projected into tangent space, before the tangents and the spring force are computed, rather than the total force at the end. Projecting either one is equivalent, and the second form is left commented out in nebm_step.

What the two functions do with the resulting force differs, and this is where the spin length is kept equal to one. Sundials_RHS writes the total force \(\mathbf{G}\) into ydot and then calls

nebm_clib.compute_dYdt   # dY/dt = Y x Y x G - correction factor

so the band is evolved with the \(\mathbf{Y}\times\mathbf{Y}\times\mathbf{G}\) form of the equation together with the correction factor of the dynamical equation given above. step_RHS does not: the StepIntegrator normalises the spins after every Euler or Runge-Kutta step instead, and the VerletIntegrator descends along the total force \(\mathbf{G}\) directly.

If a climbing image was specified as an argument for the class, we compute its modified force in the compute_effective_force method.

The function that iterates the integrator is the relax method. On every iteration, we compute the difference with the previous step using a scaled Euclidean distance. The definitions of this process are specified in the compute_maximum_dYdt method from the ChainMethodBase class. According to the magnitude stopping criteria specified in the stopping_dYdt argument of relax, the iterations of the integrator will stop if the difference with the previous step is smaller than stopping_dYdt.

Convergence criteria and force units

relax can stop on either of two criteria. stopping_dYdt compares how much the coordinates of the band moved in the last step, which is convenient but depends on the step the integrator happened to take: it can become small simply because CVODE shortened its step, even far from equilibrium. stopping_max_force instead thresholds the largest residual force on the band, which does not have that weakness, and is off by default.

The forces are reported and thresholded in the raw units of the effective field, and the first line that relax writes to the log states which those are for the simulation at hand:

Relaxation parameters: stopping_dYdt=0.0001, time_step=1e-06 s,
max_iterations=2, stopping_max_force=0.01,
forces in [A/m] over material sites

• micromagnetic simulations: amperes per metre, as for any effective field. This is the same kind of quantity as OOMMF’s \(|m \times H \times m|\), so a threshold of around 1e-2 means here what it means there.

• atomistic simulations: tesla, since the atomistic effective field is an energy per magnetic moment (\(\mathbf{H} = -\partial E / \partial (\mu_s \mathbf{m})\)). Thresholds have to be chosen on that scale instead.

The three quantities in the per-step log line, max|G|, max|H| and max|F_k|, are all in those units, so they can be compared with one another: G is assembled from the energy gradient and the spring force, so the log shows directly whether the spring force is dominating the band. max|G| is exactly the value that stopping_max_force is compared against, so a threshold can be read straight off a previous run. The same values are collected in self.G_log and written to <name>_G_log.txt.

Two points are worth keeping in mind when reading these numbers.

Only the sites that carry material are included. In a patterned sample the surrounding cells have \(M_s = 0\), but they still see the stray field of the magnetised region: the exchange and anisotropy fields vanish there because they are scaled by the magnetisation, the demagnetising field does not. Nothing relaxes those cells, so their \(|G|\) stays fixed for the whole relaxation and would hold max|G| at a floor that stopping_max_force could never get below.

The gradient is projected onto the tangent space of each spin before it is used, in nebm_step. This matters most in the atomistic case, where the exchange field is the full molecular field, of the order of \(z J / \mu_s\) and so hundreds of tesla, but points essentially along \(\mathbf{m}\) and exerts no torque. What is left after the projection comes from the curvature of the texture, of order \((J/\mu_s)(a / \lambda)^2\) for a structure of width \(\lambda\), which is the part that actually competes with the anisotropy and Zeeman terms.

self.log_energy_scale divides the three logged values, for display only. It is 1 by default, and the log records the division in the units field when it is not, since the threshold then no longer matches the printed number.

CVODE linear solver and memory usage

The initialise_integrator method exposes the choice of linear solver that CVODE uses inside its Newton iteration, through the linear_solver, maxl and maxrs arguments:

• linear_solver='spgmr' (default) uses restarted GMRES, i.e. SUNDIALS’ SUNLinSol_SPGMR, with a Krylov subspace of dimension maxl and at most maxrs restarts.

• linear_solver='diag' uses CVODE’s diagonal approximate Jacobian solver, CVDiag, which allocates no Krylov basis at all but converges more slowly.

The NEBM classes do not supply an analytical Jacobian-times-vector function, so GMRES runs unpreconditioned and CVODE approximates the Jacobian-vector products by difference quotients.

Memory. GMRES keeps an orthonormal basis of the Krylov search space, which costs maxl + 1 copies of the whole band:

memory = (maxl + 1) * n_images * n_dofs_image * 8 bytes

SUNDIALS reserves this basis up front, on the first call to the integrator. The pages are only faulted in as GMRES actually uses the vectors, so an oversized maxl does not appear as a single large allocation: the resident memory of the process instead creeps upwards during a relaxation, in jumps of exactly one band, every time GMRES needs one more iteration than it has ever needed before. This is easily mistaken for a memory leak. The virtual size (VmSize) is the quantity that reveals the reservation, and it stays constant after the first integrator step.

For reference, a band of 26 images over a \(200 \times 200\) mesh has n_dofs_image = 120000, so each copy of the band is about 24 MB. The value maxl = 300 that was hardcoded in the Cython wrapper until 2026 therefore reserved about 7 GB for a single relaxation. The defaults are now maxl = 30 and maxrs = 10, which bound the basis to about 740 MB for that band while keeping the same total iteration reach of maxl * (1 + maxrs) = 330 iterations.

Choosing maxl. Restarting is not merely cheaper in memory. The orthogonalisation cost of a GMRES cycle of \(m\) iterations grows as \(O(m^2)\), so a long cycle is also expensive in CPU time and memory bandwidth; restarted GMRES(30) reaches 300 iterations at roughly a tenth of the orthogonalisation work of a single full GMRES(300) cycle. The cost of restarting is that GMRES loses the residual-minimisation property over the discarded subspace, so it may need more iterations overall, and each iteration is one evaluation of the effective field over every image of the band. When GMRES converges within maxl iterations no restart is performed and the results are identical to those of full GMRES.

If a system needs repeated restarts, and CVODE responds by cutting the step size, it is better to raise maxl than maxrs, since a larger single cycle converges better than more restarts. Note, however, that the underlying issue in that regime is the absence of a preconditioner rather than the size of the subspace.

Spherical coordinates code

This class follows the same process as the Geodesic distances code. The main difference is that in the nebm_step process, the projections are not performed and the distances are computed using the Euclidean distance of the angles.

For spherical coordinates, the vectors are smaller, with n_dofs_image = 2 * n_spins, where

band = [ theta(0)_{0}  phi(0)_{0} theta(0)_{1} ... phi(0)_{n_dofs_image-1}
         theta(1)_{0}  phi(1)_{0} theta(1)_{1} ... phi(1)_{n_dofs_image-1}
         ...
         theta(n_images-1)_{0} ... phi(n_images-1)_{n_dofs_image-1}
       ]

The Sundials_RHS function does not include the correction factor since spherical coordinates have implicit the constraint of fixed length for the magnetisation, thus every image is simply evolved with

\[\frac{\partial \mathbf{Y}_i}{\partial \tau} = \mathbf{G}_{i}\]

Gradient

The degrees of freedom are now the angles, so the gradient of the energy is taken with respect to them. Writing a spin as \(\mathbf{m}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)\) and using the chain rule with the effective field, we have, for every spin

\[\frac{\partial E}{\partial \theta} = - w\, \mathbf{H}_{\text{eff}} \cdot \frac{\partial \mathbf{m}}{\partial \theta} , \qquad \frac{\partial E}{\partial \phi} = - w\, \mathbf{H}_{\text{eff}} \cdot \frac{\partial \mathbf{m}}{\partial \phi} = - w \sin\theta\, \mathbf{H}_{\text{eff}} \cdot \hat{\boldsymbol{\phi}}\]

where \(w=\mu_{0}M_{s}\Delta V\) in the micromagnetic case and \(w=\mu_{s}\) in the atomistic case. The negH array stores these components divided by \(w\), i.e. in the units of the field, as in the other classes. It is worth noticing the \(\sin\theta\) factor in the \(\phi\) component, since \(\hat{\boldsymbol{\phi}}\) is a unit vector but \(\partial\mathbf{m}/\partial\phi\) is not. Sites without material have no energy, so their gradient is set to zero and their angles do not change.

Tangents and distances

The tangents follow the same rules of Henkelman and Jónsson [1] as in the other classes, but they are computed in Python, in compute_tangents, because the differences of the angles have to be redefined before they are weighted by the energy differences. When a spin crosses \(\phi=\pm\pi\) the raw difference of its azimuthal angle is close to \(2\pi\), although the spin barely moved, thus every difference is brought into the \([-\pi, \pi]\) range, keeping its sign. The tangents are then normalised.

The distance between two images is the Euclidean distance of the angles of the spins with material, now without any scaling,

\[d_{j,k} = \left\{ \sum_{a} \left[ \left( \theta_{a}^{(j)} - \theta_{a}^{(k)} \right)^{2} + \left( \phi_{a}^{(j)} - \phi_{a}^{(k)} \right)^{2} \right] \right\}^{1/2}\]

This is the same metric in which the tangents are normalised, so that the projection of the force, the spring force and the interpolation of the energy band, which uses \(\partial E/\partial x\) along the path, all agree.

When an angle leaves its range during the relaxation it is redefined in correct_angles: a polar angle outside \([0, \pi]\) is reflected back, and the azimuthal angle of that spin is turned by \(\pi\), otherwise the reflection would describe a different spin. The azimuthal angles are kept in \([-\pi, \pi]\).

Comparison with the Geodesic code

Both classes find the same energy barriers. We compared the number of effective field evaluations (calls to nebm_step) using the same CVODE integrator and tolerances, and stopping_dYdt=1e-4. The first number is the number of evaluations to get the energy barrier within a relative error of \(10^{-5}\), and the second one the total number until the relaxation stops:

System	Geodesic	Spherical
Two particles (this test system)	675 / 802	546 / 902
Wire of 40 spins, rotating through y	581 / 1028	780 / 1217
Wire of 40 spins, rotating through z	585 / 1016	231 / 2103

The wire is 40 cells of 2 nm with exchange and a uniaxial anisotropy along \(x\), reversed through a middle image pointing close to \(y\) or to \(z\). Notice that stopping_dYdt is measured in the angles in one class and in the spin components in the other, so the totals are not exactly the same criterion. In the last wire the spins pass close to \(\theta=0\), where the azimuthal angle is not defined: the barrier is found early, but the band takes twice as long to relax. Accordingly, spherical coordinates are only an option when the spins stay away from the poles, and even then there is no clear advantage over the Geodesic code, which is also faster per step, since its tangents and forces are computed in C.

We also tried using the metric of the sphere, \(\mathrm{d}\theta^{2}+\sin^{2}\theta\,\mathrm{d}\phi^{2}\), in this class, i.e. Geodesic distances and the force obtained by dividing the \(\phi\) component of the gradient by \(\sin^{2}\theta\), which is the Geodesic NEBM written in the angles. It gave the same barriers and practically the same number of evaluations away from the poles, and it was even slower close to them (about 3300 evaluations in total for the last wire), because of the \(1/\sin^{2}\theta\) factor. For this reason it was not kept.

[1] (1,2,3)

Henkelman, G. & Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. The J. Chem. Phys. 113, 9978–9985 (2000)

[2]

Dittrich, R. et al. A path method for finding energy barriers and minimum energy paths in complex micromagnetic systems. J. Magn. Magn. Mater. 250, 12–19 (2002).

[3] (1,2,3,4,5,6,7,8)

Bessarab, P. F., Uzdin, V. M. & Jónsson, H. Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation. Comput. Phys. Commun. 196, 1–37 (2015)

[4] (1,2)

Cortés-Ortuño, D. et al. Thermal stability and topological protection of skyrmions in nanotracks. Preprint at arXiv:1611.07079.

[5]

Abert, C. Micromagnetics and spintronics: models and numerical methods. Eur. Phys. J. B 92, 120 (2019)

[6]

Fabian, K. & Shcherbakov, V. P. Energy barriers in three-dimensional micromagnetic models and the physics of thermoviscous magnetization. Geophys. J. Int. 215, 314–324 (2018)