Fast and Robust Algorithm for the Energy Minimization of Spin Systems Applied in an Analysis of High Temperature Spin Configurations in Terms of Skyrmion Density

TL;DR

This paper introduces a geometric energy minimization algorithm for magnetic systems, effectively analyzing thermal configurations and skyrmion densities, with improved convergence speed demonstrated in ferromagnetic and itinerant magnet models.

Contribution

A novel geometric algorithm for energy minimization that preserves magnetic moment length and accelerates convergence using LBFGS, applied to skyrmion analysis in magnetic systems.

Findings

Skyrmions appear within a narrow temperature range.

The algorithm accurately identifies local energy minima.

LBFGS significantly speeds up convergence.

Abstract

An algorithm for the minimization of the energy of magnetic systems is presented and applied to the analysis of thermal configurations of a ferromagnet to identify inherent structures, i.e. the nearest local energy minima, as a function of temperature. Over a rather narrow temperature interval, skyrmions appear and reach a high temperature limit for the skyrmion density. In addition, the performance of the algorithm is further demonstrated in a self-consistent field calculation of a skyrmion in an itinerant magnet. The algorithm is based on a geometric approach in which the curvature of the spherical domain is taken into account and as a result the length of the magnetic moments is preserved in every iteration. In the limit of infinitesimal rotations, the minimization path coincides with that obtained using damped spin dynamics while the use of limited-memory quasi-newton minimization algorithms, such as the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm, significantly accelerates the convergence.

Tables

Table 1. Table 1: Number of iterations required to reach various convergence criteria, defined as tolerance in the magnitude of the gradient of the energy on a spin. One iteration in dis-LL requires two evaluations of the gradient when the SIB algorithm is used.

Tolerance (meV)	dis-LL	VPO	CG	LBFGS
${10}^{-4}$	45300	34300	4560	715
$5\cdot {10}^{-5}$	49600	37600	4850	740
${10}^{-5}$	59750	45200	5750	750

Table 2. Table 2: Performance results averaged over 40 initial random configurations in terms of the number of energy and gradient evaluations (no. e/g), number of iterations (no. iter), and elapsed time in seconds required for various algorithms to converge on the energy minimum. The convergence is deemed to be achieved when the maximum torque acting on any one of the magnetic moments is less then 10 − 5 superscript 10 5 10^{-5} meV. All algorithms converge to the same energy minimum for these 40 initial configurations. The algorithms were implemented using Python and the Numpy package and the calculations performed on a single Intel Xeon Gold 6130 CPU (2.10GHz).

Method	no. e/g	no. iters	elapsed time
SIB-dis-LL	124309	62155	147.8
OSO-VPO	47041	47041	61.0
OSO-FR-CG	4827	4333	6.2
SN-PR-CG	2542	1826	4.0
OSO-LBFGS	724	687	1.0

Table 3. Table 3: Performance of various minimization algorithms in a calculation of a skyrmion in the NCAA model. The number of self-consistent field steps (no. SCF steps), number of energy and gradient evaluations (no. e/g), number of iterations (no. iters), and elapsed time in seconds are given. The algorithms were implemented using Python and the Numpy package and the calculations performed on 8 Intel Xeon Gold 6130 CPUs (2.10GHz).

Method	no. SCF steps	no. e/g	no. iters	elapsed time
SN-PR-CG	4566	423	285	816
OSO-FR-CG	4379	383	328	780
OSO-LBFGS	1589	126	120	281

Equations

$$E=E({\mathbf{e}}_1,{\mathbf{e}}_2..{\mathbf{e}}_N)$$

$${\mathbf{e}}_i=U^i{\mathbf{e}}_i^{\prime }.$$

$$U^i=e^{-A_i},$$

$$A_i=\left(\begin{array}{ccc}0& a_{12}^i& a_{13}^i\\ -a_{12}^i& 0& a_{23}^i\\ -a_{13}^i& -a_{23}^i& 0\end{array}\right).$$

$$E_{min}=\underset{\vec{a}\in {\mathbb{R}}^{3N}}{\min }F(\vec{a}),$$

$$\vec{a}=(a_{12}^1,a_{13}^1,a_{23}^1,...,a_{12}^i,a_{13}^i,a_{23}^i,...,a_{12}^N,a_{13}^N,a_{23}^N),$$

$$F(\vec{a})=E(e^{-A_1}{\mathbf{e}}_1^{\prime },e^{-A_2}{\mathbf{e}}_2^{\prime }\dots ,e^{-A_N}{\mathbf{e}}_N^{\prime }).$$

$$g_{\alpha \beta }^i:=\frac{\partial F}{\partial a_{\alpha \beta }^i}={\left({\int }_0^1{e}^{t{A}^i}{T}^i{e}^{-t{A}^i}dt\right)}_{\alpha \beta },$$

$$T^i=\left(\begin{array}{ccc}0& t_{iz}& -t_{iy}\\ -t_{iz}& 0& t_{ix}\\ t_{iy}& -t_{ix}& 0\end{array}\right),{\vec{t}}_i={\mathbf{e}}_i\times \frac{\partial E}{\partial {\mathbf{e}}_i}.$$

$${\mathbf{e}}_i^{(k+1)}=e^{-A{}_i^{(k)}}\cdots e^{-A{}_i^{(1)}}{e}^{-A{}_i^{(0)}}{\mathbf{e}}_i^{\prime },\forall i\in 1,2,..,N$$

$${\mathbf{e}}_i^{(k+1)}=e^{-(A{}_i^{(k)}+\cdots +A{}_i^{(1)}+A{}_i^{(0)})}{\mathbf{e}}_i^{\prime },\forall i\in 1,2,..,N.$$

$$g_{\alpha \beta }^i=T_{\alpha \beta }^i$$

$${\mathbf{e}}_i^{(k+1)}=e^{-A{}_i^{(k)}}{\mathbf{e}}_i^{(k)},\forall i\in 1,2,..,N.$$

$${\theta }_i=\sqrt{-\frac{1}{2}\mathrm{Tr}\left[{\left(A^i\right)}^2\right]}=\sqrt{{\left(a_{12}^i\right)}^2+{\left(a_{13}^i\right)}^2+{\left(a_{23}^i\right)}^2}$$

$${\mathbf{r}}_i=\frac{1}{{\theta }_i}(-a_{23}^i,a_{13}^i,-a_{12}^i{)}^T.$$

$$e^{A^i}=e^{\theta {A^{\prime }}^i}=I+\sin (\theta ){A^{\prime }}^i+\left(1-\cos (\theta )\right){\left({A^{\prime }}^i\right)}^2,$$

$${\theta }_i\leftarrow {\theta }_i\frac{{\theta }_{max}}{{\theta }_{rms}}\text{if }\frac{{\theta }_{max}}{{\theta }_{rms}}<1,$$

$$e^{-A^i}{\mathbf{e}}_i=\cos {\theta }_i{\mathbf{e}}_i+\sin {\theta }_i{\mathbf{r}}_i\times {\mathbf{e}}_i+(1-\cos \theta ){\mathbf{r}}_i\cdot {\mathbf{e}}_i{\mathbf{r}}_i.$$

$$\frac{d{\mathbf{e}}_i}{dt}=-\frac{\gamma }{M_i}{\mathbf{e}}_i\times {\vec{b}}_i^{\mathrm{eff}}-\alpha \frac{\gamma }{M_i}{\mathbf{e}}_i\times {\mathbf{e}}_i\times {\vec{b}}_i^{\mathrm{eff}},$$

$${\vec{b}}_i^{\mathrm{eff}}=-\frac{\partial E}{\partial {\mathbf{e}}_i}$$

$$\frac{d{\mathbf{e}}_i}{dt}=-\alpha \frac{\gamma }{M_i}{\mathbf{e}}_i\times {\mathbf{e}}_i\times {\vec{b}}_i^{\mathrm{eff}}.$$

$${\mathbf{e}}_i=e^{\lambda T^i}{\mathbf{e}}_i^{\prime }\approx (I+\lambda T^i){\mathbf{e}}_i^{\prime }={\mathbf{e}}_i^{\prime }+\lambda {\mathbf{e}}_i^{\prime }\times {\mathbf{e}}_i^{\prime }\times \frac{\partial E}{\partial {\mathbf{e}}_i^{\prime }}.$$

$$\frac{d{\mathbf{e}}_i^{\prime }}{dt}=\alpha \frac{\gamma }{M_i}{\mathbf{e}}_i^{\prime }\times {\mathbf{e}}_i^{\prime }\times \frac{\partial E}{\partial {\mathbf{e}}_i^{\prime }}$$

$$E=-\frac{1}{2}\sum _{<j,k>}(J_{jk}{\mathbf{e}}_j\cdot {\mathbf{e}}_k+{\vec{D}}_{jk}\cdot ({\mathbf{e}}_j\times {\mathbf{e}}_k))-K\sum _j{e}_{j,z}^2-M\vec{H}\sum _j{\mathbf{e}}_j,$$

$$\hat{H}=\sum _{i,\alpha }{E}_i^{\alpha }{\hat{n}}_{i\alpha }+\sum _{i,j,\alpha ,\beta }{V}_{ij}^{\alpha \beta }{\hat{d}}_{i\alpha }^{†}{\hat{d}}_{j\beta }-\frac{1}{4}\sum _i{U}_i\left(N_i^2-M_i^2\right),$$

$$E_i^{\alpha }=E_i^0+\frac{U_i}{2}\left(N_i-\alpha \cos {\theta }_i{M}_i\right),$$

$$V_{ij}^{\alpha \beta }=\frac{U_i}{2}\left({\delta }^{\alpha \beta }-1\right){\delta }_{ij}\exp \left(-\alpha \mathrm{i}{\varphi }_i\right)\sin {\theta }_i{M}_i+\left(1-{\delta }_{ij}\right){\delta }^{\alpha \beta }{V}_{ij}.$$

$$N_i=⟨{\hat{\tilde{d}}}_{i+}^{†}{\hat{\tilde{d}}}_{i+}⟩+⟨{\hat{\tilde{d}}}_{i-}^{†}{\hat{\tilde{d}}}_{i-}⟩=⟨{\hat{d}}_{i+}^{†}{\hat{d}}_{i+}⟩+⟨{\hat{d}}_{i-}^{†}{\hat{d}}_{i-}⟩,$$

$$\begin{array}{rl}{M}_i& =⟨{\hat{\tilde{d}}}_{i+}^{†}{\hat{\tilde{d}}}_{i+}⟩-⟨{\hat{\tilde{d}}}_{i-}^{†}{\hat{\tilde{d}}}_{i-}⟩=\cos {\theta }_i\left(⟨{\hat{d}}_{i+}^{†}{\hat{d}}_{i+}⟩-⟨{\hat{d}}_{i-}^{†}{\hat{d}}_{i-}⟩\right)\\ & +\sin {\theta }_i\left(e^{-\mathrm{i}{\varphi }_i}⟨{\hat{d}}_{i+}^{†}{\hat{d}}_{i-}⟩+e^{\mathrm{i}{\varphi }_i}⟨{\hat{d}}_{i-}^{†}{\hat{d}}_{i+}⟩\right).\end{array}$$

$${\mathbf{e}}_i={\left(\sin ({\theta }_i)\cos ({\varphi }_i),\sin ({\theta }_i)\sin ({\varphi }_i),\cos ({\theta }_i)\right)}^T$$

Full text

††journal: Computer Physics Communications

1 Introduction

Theoretical characterization of magnetic materials typically start with the identification of stable and/or metastable states of the system. These correspond to local minima on the energy surface, i.e. the energy as a function of the variables corresponding to the various degrees of freedom. At low temperature, a system in thermodynamic equilibrium is characterized by the global minimum on the energy surface. At a higher temperature, it explores regions of the energy surface favored more by entropy. For example, when a ferromagnetic material is heated up, thermal fluctuations in the spin configurations can bring the system into states with local non-collinear ordering such as skyrmions[1, 2, 3]. This has been proposed theoretically based on Heisenberg-type model [4] and demonstrated in experiments where short laser pulses briefly heat up the system followed by rapid cooling.[5]

High temperature configurations of the spins can be characterized by identifying the nearest local minima on the energy surface, i.e. minima that are connected to the high temperature configurations by energy minimization. The minimization provides a mapping between any configuration of the spins and a local minimum on the energy surface. The configuration at the local minimum is directly related to the given configuration in the sense that no energy barrier separates the two. Such local minima on the energy surface are referred to as ‘inherent structures’ of the system at a finite temperature. Analogous mapping is often used to characterize atomic configurations of liquids[6, 7].

In the present article, a fast method for minimizing the energy of magnetic systems is presented and applied in a study of the inherent structures of a ferromagnet over a wide range in temperature. A search for a local minimum on the energy surface requires an efficient algorithm for minimizing the energy. For magnetic systems, the minimization of the energy is challenging because it involves moving on a curved manifold. The energy of a system of $N$ spins is effectively a function of the orientation of the spin vectors

[TABLE]

since the length of the spin vectors is either fixed as in a Heisenberg-type Hamiltonian, or it is treated as a fast variable within an adiabatic approximation [9] in density functional theory (DFT) [10, 11] or semi-empirical model Hamiltonian calculations such as the non-collinear Alexander-Anderson (NCAA) model [12]. Well-established minimization algorithms that are widely used in atomic simulations to determine the location of atoms at an energy minimum cannot be applied to spin systems without significant modifications. It is possible to use Cartesian coordinates for the spin vectors, $\mathbf{e}_{i}=(e_{ix},e_{iy},e_{iz})$, but then either a normalization constraint needs to be added, $|\mathbf{e}_{i}|^{2}=1$ for $i=1...N$, by introducing Lagrange multipliers and thereby increasing the dimensionality of the problem, or a renormalization of the spin vectors needs to be introduced [13, 14]. Alternatively, one can use a spherical coordinate system or a stereographic coordinate system, but then a special treatment is required near the poles, $\theta=0$ and $\pi$ [16, 17, 18]. The choice of coordinate system, such as Cartesian, spherical, or stereographic, can influence both the numerical accuracy and the computational effort needed to reach convergence.

An alternative approach is to use an orthogonal optimization algorithm [19]. This approach has been successfully applied to a variety of problems, including signal processing [20] and electronic structure calculations. [21, 22, 23, 24] The advantage of the orthogonal spin optimisation is that it takes into account the curvature of the manifold on which the energy is defined and, therefore, can be combined with algorithms developed for unconstrained minimisation such as conjugate gradient algorithms and quasi-Newton methods. Here, we apply it to spin systems and show that orthogonal spin optimization (OSO) employed in combination with the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method [25] including inexact line search provides fast and stable convergence for large spin systems where the inherent structures reveal the presence of magnetic skyrmions. This OSO-LBFGS algorithm outperforms commonly used approaches based on the deterministic Landau-Lifshitz (LL) equation. [26]

The article is organised as follows. In section 2, the orthogonal optimization algorithm is presented. In section 3, results of numerical tests are presented where the orthogonal optimization is carried out in combination with the LBFGS method, the Fletcher–Reeves nonlinear conjugate gradient (CG) method [25] or the velocity projection optimization (VPO) method [27, 17] algorithms. Performance is also compared with two LL dynamics approaches, the damped LL and dissipative LL dynamics, as well as with the spin normalization (SN) algorithm of Cohen et al [13] combined with Polak-Ribière nonlinear conjugate gradient (CG) method. The SN method is commonly used in micromagnetics simulations [14, 15]. Finally, discussion and conclusions are presented in section 5. Details of the implementations of the orthogonal optimization algorithm are presented in the Appendices. The OSO-LBFGS algorithm has been implemented in the LAMMPS (https://lammps.sandia.gov) software [28, 29] and in the Spirit software [30]. The code as well as initial and final spin configurations of the calculations presented here are available in a GitLab repository. [31]

2 Methodology

First, the OSO algorithm is introduced and its implementation as a minimisation method. Then, two commonly used minimisation methods based on LL dynamics are discussed for comparison. The methods are illustrated and compared using calculations for a simple two-dimensional system using the energy surface shown in Fig. 1.

Consider a system consisting of N interacting spins and let $\{\mathbf{e}_{i}^{\prime}\}_{i=1}^{N}$ be a reference orientation of the spins. Then, any spin configuration, $\{\mathbf{e}_{i}\}_{i=1}^{N}$, can be obtained from this reference by applying an orthogonal transformation

[TABLE]

The objective is to find the set of orthogonal matrices that transform the reference spin configuration to the minimum energy configuration. Since orthogonal matrices must satisfy orthonormality constraints, ${U^{i}}^{T}U^{i}=I$, it is convenient to parametrise them using exponentials of skew-symmetric matrices, $A_{i}$:

[TABLE]

where

[TABLE]

Given that the reference spin vectors satisfy the normality constraints $|\mathbf{e}_{i}^{\prime}|^{2}=1$, the transformed spin vectors $\{\mathbf{e}_{i}\}_{i=1}^{N}$ also satisfy the normality constraints for any set of $\{A_{i}\}_{i=1}^{N}$.

Skew-symmetric matrices form a linear space and the minimum of the energy can be found as

[TABLE]

where $\vec{a}$ is a $3N-$dimensional vector

[TABLE]

$i$ refers to the spin site, and

[TABLE]

With this formulation, the energy can be minimized as a function of a $3N$-dimensional vector instead of a minimization with respect to the spin orientations subject to normality constraints. In practice, the reference spin orientation is an initial configuration, chosen in the present application from a finite temperature simulation. It can also represent a guess for the minimum energy spin configuration in order to obtain convergence with as few iterations as possible, or it can be chosen at random when local minima on the energy surface are being sampled.

The energy gradient with respect to $\vec{a}$ needs to be evaluated in order to minimize the energy in an efficient way. Within the OSO approach, the gradient is

[TABLE]

where the matrix $T^{i}$ is

[TABLE]

This equation can be obtained from the chain rule of differentiation and the definition of the first directional derivative of the matrix exponential [32].

An iterative minimization can be carried out in two ways, either

[TABLE]

or

[TABLE]

The advantage of Eq. (10) is that each rotation is small and the gradient vector at each iteration can be calculated as

[TABLE]

since the reference spins are updated at each iteration in Eq. (10)

[TABLE]

Both iteration processes have been implemented and tested. We have not found a significant difference in the number of energy and gradient evaluations while the evaluation of gradients in Eq. (10) requires less computational time. In order to reduce further the computational effort, the Cayley transformation could be used instead of the matrix exponential as it does not require calculation of trigonometric functions [42].

The skew-symmetric matrices used in the matrix exponential have a physical interpretation for rotations in spin systems. Namely, the quantity

[TABLE]

defines the rotation angle of the $i$th spin around the axis

[TABLE]

Using this interpretation, one can calculate the matrix exponential using Rodrigues’ formula [33] as

[TABLE]

where ${A^{\prime}}^{i}=A^{i}/\theta$ (see Appendix H for details). This opens the possibility of introducing a cutoff step length along the search direction in the minimization algorithm: if the root-mean-square angle $\theta_{rms}=\sqrt{\sum_{i=1}^{N}\theta_{i}^{2}/N}$ is larger than some threshold angle, $\theta_{max}$, then the rotation angles are rescaled but the rotation axes kept the same by carrying out the transformation:

[TABLE]

This makes the algorithm stable and can be particularly important in the beginning of a minimization if the line search procedure is not performed. In practical implementation, instead of calculating the matrix exponential (16) one can calculate the action of the matrix exponential on a vector as

[TABLE]

By defining the energy of the spin system as a function of variables in linear space and having the corresponding expression for the energy gradient, well-established algorithms can be used to carry out the energy minimization. We have tested three different minimization methods for finding the coefficients in the vector $\vec{a}$, namely the LBFGS, CG and VPO methods. The details of our implementation of the minimization algorithms and the calculation of the orthogonal matrix and gradient of the energy are given in A-H.

In simulations based on electronic structure calculations such as the DFT or NCAA methods as well as in calculations of large systems including many spins and long range interactions in a Heisnberg-type Hamiltonian, the evaluation of the energy and its gradient, i.e. the effective field, becomes the most time-consuming operation. Therefore, the computational effort in a minimization calculation is quantified here by the number of evaluations of the energy and its gradient.

A simple test problem involving two degrees of freedom is characterized by the energy surface shown in Fig. 1. The results of two types of OSO minimisation approaches are shown, the OSO-LBFGS and OSO-VPO algorithms.

The OSO-VPO calculation follows essentially a steepest descent path, but because of the acceleration that occurs when the gradient points in a similar direction in consecutive iterations, the steps size along the path varies. [27, 17] The method is quite stable but suffers from a slowing down as the magnitude of the gradient becomes small and the path curves near the minimum, leading to a small step size, see Fig. 1.

The combination of OSO with the LBFGS algorithm is more efficient. As can be seen in Fig. 1, the path taken from the initial point to the minimum is rather direct and involves much fewer iterations than for the OSO-VPO method.

In many cases the minimisation of the energy of a magnetic system is carried out using damped spin dynamics obtained by integrating the deterministic Landau-Lifshitz (LL) equation [26]

[TABLE]

where $\mathbf{e}_{i}$ is a unit vector defining the direction of the magnetic moment of the $i$th spin, $\gamma$ is the gyromagnetic ratio, $M_{i}$ is the length of the magnetic moment, $\vec{b}_{i}^{\rm eff}$ is the effective field

[TABLE]

and $\alpha$ is the damping parameter representing some coupling of the system to a heat bath. Along the trajectory obtained from this equation, the energy of the system leaks out through the damping term, eventually reaching a local minimum. We will refer to this as the ‘damp-LL’ minimisation method. As can be seen from Fig. 1, the calculation gives a long and winding path from the initial point to the minimum. This is clearly not an efficient minimisation method. Nevertheless, it is often used in computational studies of magnetic systems. Furthermore, since the energy leaks only gradually out of the system, it is possible that an energy barrier, lower than the energy of the initial point, is traversed during the minimisation, so the method may not provide a mapping to the inherent structure.

A more direct path to the minimum can be obtained by eliminating the oscillatory term in Eq. (19) and including only the dissipative term

[TABLE]

We will refer to this as the ‘dis-LL’ method and a calculated path for the two-dimensional test problem is shown in Fig. 1. The path obtained corresponds essentially to steepest descent on the energy surface, similar to OSO-VPO, but requires more iterations. Near the minimum especially, where the magnitude of the gradient is small, the method requires many evaluations for the final approach, thereby reducing its efficiency.

The orthogonal optimization method can be shown to reduce to the steepest descent and dis-LL in the limit of small rotations. For a small rotation, $\|A^{i}\|\ll 1$ for all $i=1\dots N$, the gradient becomes $g_{\alpha\beta}^{i}=T^{i}_{\alpha\beta}$ and a rotation in the steepest descent direction at $A^{i}=0$ is

[TABLE]

If one chooses $\lambda=\gamma\alpha\Delta t/M_{i}$ and $\Delta t\rightarrow 0$, then the equation above transforms to

[TABLE]

which is the same as Eq. (20). Thus, in the limit of small rotations, where the orthogonal optimisation becomes equivalent to a steepest descent algorithm, it corresponds to the dis-LL trajectory.

3 Performances tests

The performance tests are carried out using an extended Heisenberg-type Hamiltonian with short-range interaction as well as self-consistent mean field (SCF) calculations within NCAA model. The computational complexity of an extended Heisenberg-type Hamiltonian scales linearly with number of spins in the system, while the SCF calculations scale as the third power of the number of spins in conventional implementations.

3.1 Skyrmion lattice in an extended Heisenberg-type Hamiltonian.

We will use Heisenberg-like Hamiltonian with short-range interaction which describes magnetic moments arranged on a planar square lattice.

[TABLE]

where $\bf e_{j}$ is a three-dimensional vector of unit length along the magnetic moment on site j, $\it M$ is the magnitude of the magnetic moments, which is presupposed to be the same for all sites, $\it J_{jk}$ is the exchange parameter, which we assume to be non-zero ($\it J_{jk}=J$) only for the nearest-neighbour sites, $\vec{D}_{jk}$ is the Dzyaloshinsky-Moriya vector directed along the line connecting atomic sites i and j and K is the anisotropy parameter. Anisotropy vector points perpendicular to the lattice plane. The summation $<i,j>$ runs over all pairs of nearest neighbour sites.

The supercell consists of either 20$\times$20 or 40$\times$40 spins subject to periodic boundary conditions. The parameters in the Hamiltonian are in this case $J=10$ meV, $D=J/2$ (where D is the length of the DM vector), $K=0$, and the magnetic field is perpendicular to the monolayer with $M|\vec{H}|=J/5$. The ground state of the system corresponds to ferromagnetic ordering with energy of -22 meV per spin. The energy surface has multiple local minima corresponding to skyrmions of varying density. The skyrmionic states found in the 20x20 and 40x40 systems have energy -21.94 and -21.98 meV per spin, respectively.

The performance of four algorithms is compared. The OSO is combined with either VPO, CG or LBFGS algorithms using the gradient expression given in section 2. The CG and LBFGS algorithms include inexact line search. [25] The fourth algorithm is dis-LL, Eq. (20), where the SIB algorithm is used to generate the trajectory since it preserves the length of the spins even for large time steps. [34] The damping parameter is taken to be $\alpha$=0.1 and the time step set to 0.7 ps. For a larger time step of 0.8 ps the SIB algorithm diverges. The initial spin orientations are chosen from a random distribution on a unit sphere. Convergence is considered to be achieved when the maximum magnitude of the torque acting on the magnetic moments is less than $10^{-5}$ meV, that is $\max_{i}|\vec{t}_{i}|<10^{-5}$ meV.

The rate of convergence of the energy is shown in Fig.3. The OSO-LBFGS algorithm shows the best performance, requiring ca. 220 energy/gradient evaluations to converge on the skyrmion state. The OSO-CG requires about 50% more evaluations to reach the same level of convergence. The OSO-VPO is much slower requiring more than 3000 evaluations. The dis-LL algorithm has the worst performance requiring ca. 4000 iterations. As one iteration in the SIB algorithm requires two evaluations of the gradient for calculating the predictor and corrector, the total number of energy/gradient evaluations becomes ca. 8000 in the dis-LL calculation. This shows that large savings in computer time can be achieved by using the OSO-LBFGS in energy minimizations for magnetic systems.

While the LBFGS and CG algorithms in principle require only one evaluation of the gradient at each iteration, the energy must also be estimated in order to test the strong Wolfe conditions and if these conditions are violated a line search procedure needs to be applied (see Appendix A ). For LBFGS it is in general only necessary to perform the line search in the beginning phase of the minimization. After a few minimization steps have been carried out, the unity step length along the LBFGS search direction is guaranteed to satisfy the strong Wolfe conditions if earlier steps have satisfied the conditions. [25] Therefore, most iterations in the LBFGS algorithm require only one evaluation of the energy and the gradient.

When the number of degrees of freedom is increased by adding more spins to the simulated system, the number of energy/gradient evaluations needed to reach convergence also increases. This is illustrated by simulating a 40$\times$40 lattice where spins in a 20$\times$20 subsection of the lattice initially have random spin orientations, the same as in the calculation of the smaller system shown in Fig.2(a), while the rest of the spins are oriented as in the ferromagnetic phase, see Fig.4(a). The OSO-LBFGS calculation now requires 750 evaluations to converge to a tolerance of $10^{-5}$ meV, while OSO-CG and OSO-VPO require 5700 and 45000, respectively. The dis-LL minimization requires ca. 60000 iterations (120 000 energy/gradient evaluations). Interestingly, all four methods converge on the same local minimum energy spin configuration where the density of skyrmions is ca. half as large as in the smaller 20$\times$20 system. This calculation illustrates how the addition of more degrees of freedom, without introducing additional complexity (in that the number of randomly oriented spins is the same as in the calculation illustrated in Fig.2) increases the number of evaluations required to reach convergence. The damp-LL algorithm convergences to a different minimum and as a result is not included in the comparison. Furthermore, the number of iterations exceeds 70000 for convergence.

Table 1 shows how the number of energy/gradient evaluations depends on the convergence criterion. The LBFGS requires only slightly more iterations when the convergence criterion is decreased by an order of magnitude, while the other algorithms require ca. 20-30 $\%$ more iterations.

In order to reduce the dependence on the initial random configuration, 40 different random configurations were generated [49] for which all the tested algorithms converged to the same energy minimum. The average number of energy and gradient evaluations, number of iterations as well as elapsed time are presented in Table 2. The table also includes results obtained with the SN algorithm implemented with Polak-Rebière (PR) nonlinear conjugate gradient algorithm. The SN combined with PR algorithm is slower than OSO-LBFSG by a factor of 3-4 but performs better than the OSO-FR algorithm. Note that the ratio of energy/gradient evaluations and number of iteration in the LBFGS algorithm is close to one as the step size of 1 is the natural step-size for quasi-Newton algorithms [25] while nonlinear conjugate gradient algorithms perform a line search more often.

3.2 Skyrmions in NCAA model

In the NCAA model, the electronic structure of a 3d transition metal is approximated by two bands: One representing quasi-localized d-electrons and the other representing itinerant s(p)-electrons. In order to describe non-collinear magnetic states, we use a mean field approximation at each site $i$ where a local quantization axis, $z_{i}$, is chosen to be along the local magnetic moment associated with atom $i$ [12]. The Hamiltonian for the d-subsystem in mean field approximation can be written as

[TABLE]

where

[TABLE]

[TABLE]

Here, $E_{i}^{0}$ is a renormalized energy of unperturbed d-levels, $U_{i}$ is on site Coulomb repulsing and $V_{ij}$ are the hopping parameters that contain both a contribution due to direct exchange between d-states localized on sites $i$ and $j$ and a contribution from indirect d-d coupling through the conduction band. The d-electrons are included explicitly here, while the influence of the itinerant s(p)-electrons is indirectly taken into account via the renormalization of model parameters. The hybridization of s(p) and d bands leads to broadening of the d-band and the width, $\Gamma$, is assumed to be a parameter in the model. The polar angle $\theta_{i}$ and the azimuthal angle $\phi_{i}$ define the direction of the $i$th magnetic moment with respect to the laboratory quantization axis, $z$. In Eqns. (23)-(25), the number of d-electrons, $N_{i}$, and the magnitude of magnetic moment, $M_{i}$, associated with one of the five degenerate d-orbitals at atom $i$ are

[TABLE]

[TABLE]

For a given orientation of the magnetic moments, specified by the angles $\theta_{i}$ and $\phi_{i}$, the magnetic structure is described by a set of self-consistent values of $N_{i}$ and $M_{i}$. A detailed description of the self-consistent procedure can be found in [48].

In benchmark calculations the simulation cell consists of 21x21 spins placed on a square lattice subject to periodic boundary conditions. The initial guess for spin orientations is:

[TABLE]

with $A=5a$ where $a$ is the lattice constant. The magnetic moments were relaxed until the maximum torque dropped than $10^{-6}\Gamma$. At each minimization iteration $M_{i}$ and $N_{i}$ were converged with accuracy of $10^{-10}$. The final spin configuration is shown in Fig 5.

A comparison of the performance of the various minimization algorithms for this system is presented in Table 3. The most computationally demanding part is the SCF step which requires calculation of eigenvalues and eigenvectors of the NCAA Hamiltonian. The OSO-LBFGS algorithm converges faster than the other algorithms and requires around 1600 SCF steps which is ca. a third compared with the other algorithms. While SN-PR-CG carries out fewer outer steps then OSO-FR-CG, it is still slower then OSO-FR-CG because it performs more line searches. The total wall-time is proportional to the number of SCF steps for all algorithms which shows that it is indeed the most computational demanding operation and the performance of the algorithm can be quantified by counting the number of SCF steps.

4 Application: Inherent structures at high temperature

We now apply the OSO-LBFGS algorithm to study inherent structures of a ferromagnetic layer as a function of temperature. As mentioned in the introduction, both theoretical calculations based on a Heisenberg-type model [4] and experimental studies of thin films [5] have shown that high temperature can break the ferromagnetic order and introduce skyrmions which remain in the system if it is cooled quickly enough. The experiments were carried out for a system where skyrmions are large (micron size) and stable at room temperature. The simulation of such a large system on the atomic scale is a challenging task. Here, we have instead investigated the skyrmion creation due to local heating in a well known system that has been studied extensively, namely the Co/Pt(111) system. The values of the parameters of the Hamiltonian are chosen to describe Néel skyrmions: $J=29$ meV, $|\vec{D}_{ij}|=1.5$ meV, $K=0.4$ meV as has been done previously by Rohart and coworkers [35]. The external field is set to zero, $H=0$. The effect of the dipole-dipole interaction is accounted for effectively to a good approximation by reducing the value of the anisotropy constant to $K=0.293$ meV [35, 36]. The simulated system consists of 250000 spins in the simulation cell arranged on a triangular lattice subject to periodic boundary conditions. The initial configurations for the energy minimisations are prepared by heating spins in a region within a circle of radius 251$a$ (where $a$ is the lattice constant), including 167786 spins, using stochastic Landau-Lifshitz-Gilbert (LLG) dynamics over a period of 1.0 ns in order to reach thermal equilibrium. A total of 20 statistically uncorrelated configurations are prepared by recording the orientation of the spins at 0.2 ns intervals after the equilibration period. The other spins in the system remain frozen in the ferromagnetic arrangement during this preparation phase. This is repeated for a wide range in temperature. For each spin configuration generated in this way, an OSO-LBFGS minimization (see Appendix F) is carried out with all spin orientations being free degrees of freedom. The average number of skyrmions found at the minimum energy configuration is reported in Fig. 6 as a function of the temperature from which the initial configuration was obtained. The temperature is given in units of the exchange parameter, $J$. For temperature below 2.8 $J$ (corresponding to about 800 K for the Co/Pt(111) system) no skyrmion is found. Over a narrow temperature interval, between 2.8 $J$ and 3.2 $J$ (corresponding to 950 - 1050 K), there is a steep increase in the average number of skyrmions to a saturation value of 4 skyrmions within the simulation cell. At even higher temperature (5 J and 7 J) the average number of skyrmions in the inherent structures, i.e. after minimisation, is unchanged. Similar results are obtained by using the dis-LL algorithm but then the computational effort is more than an order of magnitude larger.

5 Discussion

By using the OSO algorithm, the task of minimising the energy on a curved manifold is transformed to a minimisation in linear space making it possible to use the highly efficient LBFGS algorithm and an inexact line search. This OSO-LBFGS algorithm is found to provide an efficient method for finding local energy minima on the energy surface of magnetic systems, outperforming in particular methods based on LL dynamics. The method is used to study inherent structures of a ferromagnetic layer mimicking the widely studied Co/Pt(111) system as a function of temperature. Over a remarkably narrow temperature range, the number of skyrmions in the inherent structures is found to rise abruptly to a high temperature saturation value. These skyrmion configurations represent local minima on the energy surface that are favored over the global ferromagnetic energy minimum by entropy. In this system, the skyrmion configurations are separated by a low energy barrier from the ferromagnetic state, as has been reported previously [36] so they are only expected to remain in the system after a fast quench to low temperature. This entropic preference for configurations involving skyrmions in the inherent structures is, however, likely a general property for spin systems and is consistent with the laser heating experiments where skyrmions were induced in a quite different magnetic system [5].

The number of skyrmions can, in principle, be identified by evaluating the total topological charge of the finite temperature spin configuration and dividing by the topological charge of a single skyrmion [8]. The topological charge can be estimated as [50]:

[TABLE]

where the sum is evaluated over triplets of neighboring spins. Such an approach can be useful at relatively low temperature where thermal fluctuations present only a small distortion of the minimum energy configuration. But, at high-temperature the magnetic configuration is highly disordered and can have an arbitrary value of the topological charge. Indeed, we have verified that the high temperature spin configurations generated here give essentially random values for the topological charge and are not related to the topological charge of the inherent structures at zero temperature.

Nearly 200 minimization calculations have been carried out in the tests and the Co/Pt(111) simulations using the OSO-LBFGS method and in all cases convergence to a local energy minimum has been obtained without problems. The number of iterations needed increases as one would expect the more degrees of freedom are involved and the higher the initial temperature is. The calculations presented here indicate that the method is robust, but additional calculations in the future with application to other types of systems will further test and hopefully demonstrate its robustness and wide applicability.

The method has furthermore been applied to a Hamiltonian that describes more accurately itinerant electrons. The interaction in such systems often includes non-quadratic spin terms and cannot, therefore, be described by a Heisenberg type model. At each step of the energy minimization one needs to solve self-consistent field equations. We have shown that the minimization method presented here drastically reduces the number of SCF steps and we expect the algorithm to be useful in other self-consistent field calculations, such as density functional theory calculations.

The performance of the OSO-LBFGS algorithm could likely be improved further by developing a preconditioner for the orthogonal optimization. The dis-LL and damp-LL methods could be accelerated especially towards the end of a minimization where the gradient (i.e. effective field) is small by using a larger time step in an adaptive time step algorithm. An adaptive time step propagation based on Barzilian-Borwein rule has been used in micromagnetic simulations and found to give an acceleration by a factor of 2  [37], but even with such acceleration, the dis-LL and damp-LL algorithms require significantly larger number of energy/gradient evaluations than the OSO-LBFGS algorithm.

Some implementations make use of higher order integration methods[45, 46], such as 5th order Runge-Kutta method with adaptive time step. The application of this algorithm to the test problems presented here shows that the adaptive time step reduces the number of iterations needed for convergence, but the 5th order algorithm requires 6 evaluations of the gradient per iteration and as a result the total number of energy/gradient evaluations needed for convergence is larger than for the SIB algorithm.

The OSO algorithm should also be useful in other spin configuration optimizations, such as minimum energy path calculations using the geodesic nudged elastic band method, which has previously been implemented with a VPO algorithm mimicking particle dynamics [17, 47]. Also, the algorithm could be used in searches for first order saddle points on energy surfaces for spin systems starting only from the initial state and without knowledge of any final state [38]. From the saddle points, the rate of magnetic transitions can be estimated using harmonic transition state theory [39] as has, for example, been done to estimate the lifetime of skyrmions as a function of temperature. [40, 41]

The focus in this article is on the use of the OSO approach in optimisation, but it could also be used in calculations of the dynamics of a spin system. If one defines the $\vec{t}_{i}$ vector as $\alpha\mathbf{e}_{i}\times\frac{\partial E}{\partial\mathbf{e}_{i}}-\frac{\partial E}{\partial\mathbf{e}_{i}}+\vec{f}_{i}$, where $\vec{f}_{i}$ is a random field, then the rotations can correspond to a trajectory obtained from the stochastic LL equation. In this way an algorithm for the integration of the stochastic LLG equation can be obtained. The use of exponential transformation for integration of the stochastic LLG equation has indeed been discussed by Lewis and Nigam. [42] An exact calculation of the rotation matrix has been used by Depondt and Mertens, [43] while an approximation of the exponential transformation using Cayley transform [42] is used in the method described by Mentik et al. [34]

6 Acknowledgments

We thank Julien Tranchida for helpful discussion and help with implementation of the algorithm in the LAMMPS software. Also, we thank Gideon Müller, Moritz Sallerman and Pavel Bessarab for helpful discussion and valuable comments on the manuscript. AVI thanks Filipp Rybakov for helpful discussions. This work was funded by the Icelandic Research Fund, the University of Iceland Research Fund, and the Russian Science Foundation (Grant 19-42-06302). AVI is supported by a doctoral fellowship from the University of Iceland.

Appendix A Orthogonal Spin Optimization A

Here the iteration process from Eq. (10) is described:

Set $k=0$ and choose initial spin orientations $\{\mathbf{e}\,^{(k)}_{i}\}_{i=1..N}$. Let $\vec{t}_{i}\,^{(k)}$ be the torque acting on $i$th spin at $k$th iteration:

[TABLE]

Calculate the gradient vector:

[TABLE] 2. 2.

Calculate initial search direction $\vec{p}\,^{(k)}$ according to the particular minimization algorithm chosen, see D-F for example. 3. 3.

Calculate the maximum magnitude of the torques:

[TABLE]

and set tolerance $\epsilon$ (for example, $\epsilon=10^{-5}$ meV). 4. 4.

While $\Delta\,^{(k)}>\epsilon$:

1. (a)

Compute $\lambda\,^{(k)}$ using a line search (for VPO algorithm $\lambda\,^{(k)}=1$) and calculate using (18):

[TABLE]

where skew-symmetric matrices

[TABLE]

constructed using the search direction vector:

[TABLE] 2. (b)

Calculate new gradient $\vec{g}\,^{(k+1)}$ and $\Delta\,^{(k+1)}$. Calculate new search direction, $\vec{p}\,^{(k+1)}$, according to the particular minimization algorithm chosen, see Appendix D-F . Set $k\leftarrow k+1$ 5. 5.

End

Appendix B Orthogonal Spin Optimization B

Here the iteration process from Eq. (11) is described. The task is to minimize the energy as a function of spin orientations, $E(\mathbf{e}_{1},\mathbf{e}_{2},..\mathbf{e}_{N})$. The $\mathbf{e}_{i}$ are parametrised using skew-symmetric matrices $\mathbf{e}_{i}=\mathrm{exp}(-A_{i})\mathbf{e}\,^{\prime}_{i}$

[TABLE]

It is enough to consider only the upper diagonal part of $A_{i}$, and the energy is therefore a function of a $3N$-dimensional vector

[TABLE]

where the upper index, i, refers to the spin site. Let $F$ be

[TABLE]

Algorithm:

Choose reference spin orientations $\{\mathbf{e}^{\prime}_{i}\}_{i=1..N}$, set initial skew-symmetric matrices to zero, $\vec{a}^{(0)}=0$, $\mathbf{e}_{i}^{(0)}=\mathbf{e}^{\prime}_{i}$ and calculate initial gradient $\vec{g}\,^{(0)}={\partial F}/{\partial\vec{a}\,^{(0)}}$. 2. 2.

Set $k=0$ and calculate initial search direction $\vec{p}\,^{(0)}$ according to the particular minimization algorithm chosen, see D-F. 3. 3.

Set integer $u$ that counts the number of steps before updating the reference spins. For example, $u=50$. 4. 4.

Calculate the maximum magnitude of the torques:

[TABLE]

and set tolerance $\epsilon$ (for example, $\epsilon=10^{-5}$ meV). 5. 5.

While $\Delta\,^{(k)}>\epsilon$:

1. (a)

If $k\,\mathrm{mod}\,u=0$ then update reference spins:

[TABLE] 2. (b)

Compute $\lambda\,^{(k)}$ using a line search (for VPO algorithm $\lambda\,^{(k)}=1$) and calculate:

[TABLE] 3. (c)

Calculate new gradient $\vec{g}\,^{(k+1)}$ and $\Delta\,^{(k+1)}$. Calculate new search direction, $\vec{p}\,^{(k+1)}$, according to the particular minimization algorithm chosen, see D-F . Set $k\leftarrow k+1$ 6. 6.

End

Appendix C Choice of the step length.

The step length parameter $\lambda$ is chosen in such a way that the strong Wolfe conditions [25] and/or approximate Wolfe conditions [44] are satisfied

[TABLE]

and/or

[TABLE]

with $0<c_{1}<c_{2}<1,\epsilon>0,\delta<\mathrm{min}{.5,\sigma},\sigma<1$. The parameters were chosen according to Ref. [25] and [44] as

[TABLE]

It can be shown that after several iterations, a step length of 1 guarantees satisfaction of the strong Wolfe conditions in the LBFGS algorithm, [25] and therefore, a trial step of $\lambda=1$ is always used first to test these conditions. If they are not satisfied, an inexact line search procedure based on the cubic interpolation is used. [25] Approximate Wolfe conditions [44] are always examined at the minimum of the cubic interpolation (see G).

The other approach for the choice of the step length is based on the maximum rotation of the spins as described in the main text.

Appendix D VPO algorithm

Values of two parameters need to be chosen, the time step $\Delta t$ and an effective (artificial) mass, $m$. Since the system accelerates when the gradient points in a similar direction in subsequent iterations, the time step can be chosen to be small, $\Delta t=0.005$. The effective mass was chosen to be $m=0.01$. At the $k$-th iteration the search direction is chosen according to the following:

if k = 0:

[TABLE]

else:

[TABLE]

Appendix E CG algorithm

The Fletcher–Reeves nonlinear conjugate gradient method [25] is given here for completeness of the paper. At the $k$-th iteration the search direction is chosen according to the following

if k = 0:

[TABLE]

else:

[TABLE]

Appendix F L-BFGS algorithm

Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [25] is given here for completeness.

Algorithm: Set $k=0$. Let $m$ be amount of previous step which are stored in memory.

1. (1)

Calculate $n=k\,\mathrm{mod}\,m$ 2. (2)

if k = 0:

[TABLE]

else:

[TABLE]

for $l$ in $(m-1),(m-2)\dots 0$:

[TABLE]

end for

[TABLE]

for l in $0,1,\dots(m-1)$:

[TABLE]

end for

end else 3. (3)

return -$\vec{p}\,^{(k)}$

Appendix G Cubic interpolation

Let $f$ be defined on $[0,r]$ and $f(0),f(r),f^{\prime}(0),f^{\prime}(r)$ be known. Then the cubic interpolation of $f$ is:

[TABLE]

and the minimum of this cubic function is at

[TABLE]

or at the boundaries of the interval $[0,r]$.

Appendix H Orthogonal matrix and gradient of the energy

Let $A$ be

[TABLE]

Let $\theta=\sqrt{a^{2}+b^{2}+c^{2}}$. Then the eigenvalues are:

[TABLE]

The eigenvectors are:

[TABLE]

if $a=0$ and $c=0$ then

[TABLE]

The matrix exponential can then be calculated as:

[TABLE]

where $L$ is a diagonal matrix, $\mathrm{diag}(L)=(1,e^{-i\theta},e^{i\theta})$, and $V$ is a unitary matrix, columns of which are eigenvectors of $A$. The gradient of the energy is then $g_{\alpha\beta}={\partial E}/{\partial a_{\alpha\beta}}$:

[TABLE]

where $\circ$ denotes Hadamard product,

[TABLE]

and

[TABLE]

*Rodrigues’ formula for rotations.

*Let $q=\cos\theta,w=1-\cos\theta$ and $x=a/\theta,y=b/\theta,z=c/\theta$, $s_{1}=-yzw$, $s_{2}=xzw$, $s_{3}=-xyw$, $p_{1}=x\sin\theta$, $p_{2}=y\sin\theta$, $p_{3}=z\sin\theta$ then the matrix exponential can be calculated as:

[TABLE]