CC BY
24
Truncated minimum energy path method for finding first order saddle points
I. S. Lobanov1, M.N. Potkina2'3, H. Jonsson3'4, V.M. Uzdin1'2
xITMO University, Kronverkskiy, 49, St. Petersburg, 197101, Russia Kronverkskiy, 49, St. Petersburg, 197101, Russia 2 St. Petersburg State University, St. Petersburg, 198504, Russia 3 Science Institute and Faculty of Physical Sciences, University of Iceland 107 Reykjavik, Iceland 4Center for Nonlinear Studies, Los Alamos, NM 87545, USA [email protected], [email protected], [email protected]
PACS 82.20.Pm, 34.10.+x, 75.40.Mg DOI 10.17586/2220-8054-2017-8-5-586-595
A method for finding a selected region of the minimum energy path between two local minima on an energy surface is presented. It can be used to find the highest saddle point and thereby estimate the activation energy for the corresponding transition when the shape of the path is known reasonably well and a good guess can be made of the approximate location of the saddle point. The computational effort is then reduced significantly as compared with a calculation of the full minimum energy path by focusing the images on the selected part of the path and making one of the images, the climbing image, converge rigorously on the saddle point. Unlike the commonly used implementation where a restraint is used to distribute the images along the path, the present implementation makes use of a constraint where the distance between images is controlled based on a predefined overall length of the path. A relatively even density of images on each side of the climbing image is maintained by allowing images to move from one side to the other. Applications to magnetic skyrmion annihilation and escape through boundary are used to illustrate the savings in computational effort as compared with full minimum energy path calculations. Keywords: saddle point, minimum energy path, constraint, rate theory, skyrmion. Received: 9 October 2017 Revised: 19 October 2017
1. Introduction
The mechanism of a transition between (meta)stable states of an atomic system, involving rearrangement of atomic coordinates and/or rotation of magnetic moments, can be characterized by the minimum energy path (MEP) connecting the corresponding (local) minima on the energy surface. The MEP represents the path of highest statistical weight and is a natural choice for a reaction coordinate. Typically, the primary task is to estimate the transition rate. Within the harmonic approximation to transition state theory (HTST) this requires the evaluation of the highest energy along the MEP, a first order saddle point (SP) on the energy surface. While the MEP may have more than one energy maximum, only the highest one is needed for the HTST rate estimate as it determines the activation energy. For this purpose, the MEP only needs to be known well enough to identify the point of highest energy.
The nudged elastic band (NEB) method is frequently used to find MEPs to determine transition mechanisms and at the same time estimate the activation energy for transitions [1,2]. Given initial and final state minima on the energy surface and some initial path between them, the NEB method can be used to iteratively converge on an MEP, typically the one nearest to the initial path. The MEP does not, however, need to be close to the initial path. An NEB calculation can reveal an unexpected transition mechanism and intermediate local (or even global) minima that were not known beforehand (see, for example, ref. [3]). If, on the other hand, the shape of the MEP is quite well known and the region where the (highest) SP lies can be identified beforehand, the calculation of the full MEP is not needed and the computational effort in estimating the activation energy can be reduced by launching a SP search starting from a best estimate of the SP. The NEB calculation can, for example, be carried out only to the extent that the rough shape of the MEP becomes evident, and then a minimum mode following (MMF) method [4-6] launched from the point of highest energy [7].
In the NEB method, the path is discretized using a set of replicas of the system, i.e. sets of values of the variables defining a system configuration. The replicas are referred to as images of the system. Starting with some initial path, an iterative procedure is used to bring the images to an MEP. In the climbing image NEB (CI-NEB) [8], the highest energy image is treated separately and made to converge on the highest SP. This modification of the NEB method is simple and easy to implement and is a popular alternative to a saddle point search method such as
MMF. In either case, the computational effort of converging accurately on the whole MEP can be avoided if the primary goal is to find the SP and estimate the activation energy for the transition.
Another option for reducing the computational effort compared with a full MEP calculation is to focus the convergence on a limited region of the MEP, one that is known to contain the (highest) SP. There have been several different implementations of this approach. In its simplest form, intermediate images of partly converged NEB calculation are chosen to be the fixed endpoints instead of the local minima (see, for example, ref. [9]), but a better approach is to make the new endpoints follow selected equipotential contours towards the MEP [10,11].
A key aspect of the NEB method is a clear separation of the adjustments in image positions that represent a change in the shape of the path and adjustments that affect the distribution of images along the path. This is referred to as 'nudging' and it requires a numerically stable estimate of the local tangent to the path [12]. Displacements representing changes in the shape of the path are only driven by the component of the energy gradient perpendicular to the local tangent. The distribution of the images along the path is then a separate issue and can be handled in various ways. The most common implementation makes use of a restraint approach, i.e. a harmonic penalty term in the energy which gives a spring force between adjacent images acting only in the direction of the tangent [1,2]. Alternatively, a constraint approach can be used where the total length of the path is estimated and the images placed along the path so as to obtain a desired (usually even) distribution [13]. This is one of many examples in computational methodology where a desired feature of the simulation can be implemented either with a restraint or a constraint (see, for example, ref. [14]). Tests indicate that the restraint approach can converge with fewer iterations when a small but typical number of images is included in the MEP calculation, while the constraint approach ensures better distribution of the images along the path when a large number of images is included (because the springs only act between adjacent images) [15]. While the restraint approach has become widely used, the constraint approach is used less frequently. The reason for this may partly be that useful features such as the CI and ways to focus the convergence onto a limited region of the MEP have not, as far as we know, been formulated for the constraint approach. In the present article, we describe how this can be done.
The applications we present involve annihilation and escape of a magnetic skyrmion. The NEB method with restraints has been adapted to magnetic systems, the so-called geodesic NEB (GNEB) [16] and applied to various types of magnetic systems, in particular skyrmions [16-18]. Magnetic skyrmions are topologically nontrivial states of magnetic systems. The states are localized in the sense that only in a finite domain (called the support) the vectors are not directed along the external magnetic field direction. The support can be large (containing millions of spins), making the simulation domain large, thus numerical modelling of the skyrmions requires significant computational effort. During annihilation in the interior of the domain, the skyrmion first shrinks before reaching the SP and then the spins rotate to form the ferromagnetic (FM) state. We are primarily interested here in the activation energy, hence the part of the MEP close to the SP is of primary importance. The shape of the MEP is simple, has a single maximum, and is well known from previous calculations of similar systems. In calculations of skyrmion annihilation, the support can be reduced by an order of magnitude by focusing on a limited region of the MEP, thereby reducing the simulation domain and the computational effort. The skyrmion escape through a free boundary of the simulated system has similar properties: The skyrmion is repelled by the boundary, so the local energy minimum corresponds to a skyrmion located far from the boundary. Large part of the MEP represents movement of the skyrmion towards the boundary while the interesting part of the MEP is in a small region near the boundary. By focusing on a limited region of the MEP, the calculation can be carried out for a smaller simulation domain without loss of precision in the estimate of the SP energy. Calculations of these two skyrmion transitions are used here as benchmarks to calibrate the efficiency.
In the following section, we review briefly the NEB method and then discuss the constraint approach for distributing the images along the path. The formulation of the CI variant of the method within the context of the constraint approach is then described and finally the method for focusing on only a limited part of the MEP, a method we will refer to as truncated MEP (TMEP). Applications to skyrmion annihilation and escape through the boundary of the simulated system are then presented and the computational effort compared with full MEP calculations.
2. Method
2.1. Nudging and image distribution
Let E denote the energy of the system, and pn denote a set of values of all the variables defining a configuration of the sysetm. This will be referred to as an image of the system. A path in configuration space is defined as a sequence of images p = (pn), n = 0... N. The path p is an MEP between two states if
(1) p0 and pN are local minima of E,
(2) pn for n =1... N — 1 is a local minimum of energy in the plane orthogonal to the path p.
The first condition implies that the gradient of energy at the first and the last image is zero: g0 = gN = 0, gn = VE[pn]. Given a unit tangent tn to the path at the image pn, the second condition implies that the projection of the gradient of energy on the space orthogonal to tn is zero:
Ptn(gn) = gn — tn(Tn • gn) = 0.
Care must be taken in the estimation of the the tangent. A simple and numerically stable way has been presented in ref. [12]:
V+, En+1 > En > En-1, t" En+1 < En < En-1,
n
where t" = pn+1 — pn is the right tangent, t- = pn — pn 1 is the left tangent. If En is a local maximum En-1 < En A En+1 < En or a local minimum En < En-1 A En < En+1 on the path, then the tangent should be a weighted sum of t" and t- that transforms continuously between the definition above as energy becomes monotone. The tangent vector fn should be normalized:
Given some initial path p", the images are moved towards the MEP using an iterative procedure. The projected gradient PTn (gn) can be used as a measure of convergence to the MEP, the smaller the norm, the closer the convergence. The images are moved in the direction of the antigradient,
w" = p" - ag"
where a > 0 controlls the length of an optimization step. This displacement of the images should only change the shape of the path, but not affect directly the distribution of the images along the path. Therefore, the gradient is projected on the orthogonal plane of the path (nudging):
w" = p" - ah", hn = Ptn(gn).
If the configuration space is constrained (for example the fixed norm of the spins in the applications presented here), then additional steps need to be carried out to satisfy the constrains. Let qn = n(gn) be projection of gn to the tangent space of the manifold defined by the constraints, t" = n(r„)/||n(Tn)|| be projection of the approximation of the tangent vector to the tangent space, and R be operator projecting arbitrary point to the closest state satisfying the constrains. Then, the gradient descent step is defined as:
w" = R(p" - ah"), h" = Pt?(q").
The procedure described above changes the location of the images so as to adjust the shape of the path. However, this does not specify how the images are distributed along the path. If this distribution is not controlled, the images will eventually slide down to the local energy minima. The nudging algorithm separates the distribution of the images along the path from the displacements of images affecting the shape of the path. While the distribution is most commonly controlled using a restraint method where a harmonic spring acts between adjacent images [1, 2], we will here use a constraint method based on an estimate of the total length of the path [13]. A reparametrization step is introduced to distribute the images. Given pre-images w", piecewise linear interpolation of the path is constructed
(t) m (t — dn)wn+1 — (t — rfn+1)wn\ [dn ,n+1i
where d" is the distance from the pre-image w" to the first image w0 in the metric of the phase space:
d"+1 = d" +dist(w"+1,w"), d0 =0.
Almost equidistant distribution of the images on the next step is obtained by selecting images with constant step in terms of the distance d along the path:
p" ^ w(dNn/N).
T = T
2.2. Climbing image
The procedure described above will bring the images to the nearest MEP and make the magnitude of the projected gradient arbitrarily small for a sufficiently small a. But, the most important point on the MEP, the maximum which corresponds to the SP, may turn out to be in between images and the estimate of the SP energy thus inaccurate. The algorithm can be modified by making one of the images, the climbing image (CI), converge rigorously on the SP [8]. The CI is moved along the path in the uphill direction, to higher energy:
wn = R(pn - ahn), hn = qn - • qn).
It is desirable to choose the CI to be the highest energy image of the path. Hence every image on the path is marked as climbing Mn = 2 (closest to a maximum), descending Mn = 0 (closest to a minimum) or holding Mn = 1 (all other images). The general equation for the update of the images is then as follows:
wn = R(pn - ahn), hn = qn - MV^« • qn).
Climbing and descending images should not be moved during the reparametrization, otherwise they will not reach stationary points. To save positions of climbing and descending images, every new image w(dNn/N) introduced during reparametrization that is closest to a climbing/descending image should be removed and the climbing/descending image should be used instead. This guaranties that climbing/descending images will not come too close to the holding images as that could cause problems in the estimation of the tangent. This approach allows for an image to move from one side of the climbing image to the other. This helps maintain similar resolution on both sides of the CI. (This feature has, as far as we know, not been implemented previously).
The norm of the gradient at the climbing and the descending images can be used to monitor convergence. If a climbing image reaches a SP, its gradient becomes zero. Hence another important measure of the convergence to an MEP is the maximum of norm of gradients at all extremes of the energy along the path.
2.3. Truncated MEP
Above, the calculation of a full MEP between local minima has been described. When enough information about the shape of the MEP has been obtained to identify the region where the (highest) SP is located, or it is known a priori that the shape is simple with only one maximum, computational effort can be reduced by focusing the calculation on the region near the SP. Only a part of the MEP is then calculated. We will refer to this as truncated MEP (TMEP). To specify the selected part, an image called an anchor is specified. It should always belong to the TMEP. The anchor can be a point of maximum energy, or a local minimum. To compute TMEP, the same steps are carried out as described above, namely compute wn and then w(t). But when new images are introduced, they are taken from the interval t G [tmin,tmax], the increment of t is still constant:
p w(tmin + (tmax tmin)n/N).
Let n = a be an anchor, then the part of the MEP used for reparametrization is chosen in such a way that image a lies in the middle of the TMEP and the length of the TMEP is a predefined constant p, that is
tmin d p/2, tmax d + p/2.
If da < p/2, then the interval should be moved to the right, to preserve the length of the TMEP:
tmin ° tmax min(p, d ).
If da + p/2 > dN, then the interval should be moved to the left:
tmin = max(°, d p): tmax = d .
If the whole path is shorter than p, then nothing is dropped and the interval is the same as for the computation of the full MEP:
t ■ = ° t = dN
°min '-'max ^
The length p of the TMEP can be chosen to equal the length of the initial path or as some multiple of that length.
During computation of a full MEP the end images are fixed at local minima, and can be marked as descending images. The same can be done for TMEP, but the gradient at the end images will in general no longer become zero, only the projected gradient. The tangent at the ends can be found as left and right tangents:
-0 n -N N T = T+, T = T_ .
During reparametrization the indices of images can change, that should be taken into account when defining an anchor image for the next step. To preserve the anchor image, it is excluded from the reparametrization stage
in the same way as a climbing or descending image. Namely, after computation of a new image pA, the image closest to the old anchor wa, is replaced and the new anchor value set to A.
The selection of the TMEP can be done in various ways. MEP often can be divided to two parts: from the initial state to the SP, or from SP to the final state. One of the possible uses of TMEP is the calculation of only one of the two parts. In calculations of only the first part of the path, the TMEP should contain the initial state and at least one image lying on the other side of the energy ridge, possibly close to the SP. Another choice of TMEP is to find only the highest SP, i.e. the highest energy point on the MEP between two local minima. In that case, two images lying on opposite sides of the energy ridge need to be identified and the initial path set to a geodesic between the states. The anchor image in that case is the image with maximal energy. It is important to place the images at the ends of TMEP on opposite sides of the energy ridge, otherwise all images collapse to a metastable state. If the length of TMEP is insufficient and the anchor is set to an image different from the image with maximal energy, all images may end up on the same side of the energy ridge, even if initially they were separated by the energy ridge. This all requires considerable knowledge of the MEP, but can be useful when repeated calculations are carried out on similar systems or only slight changes have been made in the energy landscape, so previous knowledge can be used to speed up the calculation.
3. Applications
The TMEP method is applied here in studies of a magnetic skyrmion. Calculations of two example transitions are described below, annihilation of the skyrmion and escape of the skyrmion through a free boundary of the simulated system.
3.1. Model system
The energetics of the system are described by a Heisenberg model. Every state of the system is given by the vectors of magnetic moments pSn of each atom n in the system, where p is value of magnetic moment and Sn is a unit direction vector. The energy of the system is given by
E[S] = - Y (JjkSk • Sk + Djik • (Sj x Sk)) - Y K(K0 • Sj)2 - Y PBj • Sj, (1)
(j,k) j j
where the sum over (j, k) is taken over all pairs of neighboring atoms, B is strength of external magnetic field, J is Heisenberg exchange constant, D is Dzyaloshinskii-Moriya (DM) vector, K and K0 are strength and direction of anisotropy.
In this case, the computation of an MEP involves a constraint, because the length of the magnetic moments is fixed, ||Sn|| = 1. The operator R projects to the manifold corresponding to this restriction and the projection gn
of the gradient on the tangent space of the manifold are
S
r(S)= m' n(gn) = gn - Sn(Sn •gn).
The system is represented by a square lattice in the x-y plane, 100 x 100. The parameters in the Heisenberg model are chosen as follows: Easy axis anisotropy constant of K = 0.004 J with the axis perpendicular to the plane of the system, and DM vector directed along pairs of interacting spins of length D = 0.05 J. An external magnetic field is not applied. For these parameters, the skyrmion is a metastable state with a radius of « 21 lattice constants.
3.2. Skyrmion annihilation
In the calculations of skyrmion annihilation, periodic boundary conditions are applied in the x-y plane. The initial path is generated as a geodesic between a single skyrmion and a ferromagnetic state where all spins are directed along z. The skyrmion annihilation process consists of two stages: First, skyrmion shrinkage and then collapse of the skyrmion core. The saddle point is located in the second stage where the radius has shrunk to about 3 lattice constants. A large part of the full MEP contains the first stage and, therefore, requires a large simulation domain. By using TMEP to focus on the part of the path between transition state and ferromagnetic state (FM) the calculation of the SP can be done quicker. To obtain the full MEP we first relaxed both ends of the path to the local minima and then applied the constraint method with steps a = 0.22. TMEP was obtained as described in the previous section with the same a. FM state was set as the anchor and as one of the ends. The other end was set to a skyrmion like state of radius 10 lattice constants.
The calculated results are shown in Fig. 1. The TMEP calculation gives the same activation energy for FM to skyrmion transition as a full MEP calculation up to computational precision. Convergence of TMEP is several
Fig. 1. MEP for skyrmion annihilation (100 x 100 lattice): (left) full MEP, (right) truncated path (TMEP) with anchor at ferromagnetic state (last image). (top) first images of MEP after 2000 iterations, (middle) MEP after 2000 iterations, (bottom) convergence history. Solid black line depicts energy, dash-dotted line is norm of gradient, and dashed line is norm of the gradient component perpendicular to the path. Down-pointing triangle, circle and up-pointing triangle markers show descending, climbing and holding images, respectively. images with small markers are subjects of reparametrization step. Support of all images in TMEP is 10 times smaller and the resolution of the path in TMEP is higher than when the full MEP is calculated. SP energy was found with a precision of 0.1% five times faster using TMEP
times faster, probably because of better tangent estimate due to higher resolution of the path around the saddle point. Both the full MEP and TMEP were computed with the same number of images, a total of 8 images. Because of the lower resolution, the full MEP calculation does not contain intermediate points between the saddle point and FM state. Since the images are situated closer to the SP in the TMEP calculation, details of the transition are better resolved. First images of the paths have largest support, and the support for the full MEP calculation is two orders of magnitude larger than for TMEP. As a result, the computational effort of TMEP can be two orders of magnitude smaller under appropriate correction of simulation domain size. TMEP also gives almost correct value of the activation energy from the beginning, while the estimate of the activation energy given by the full MEP calculation slowly increases as climbing image reaches the saddle point.
Fig. 2. MEP for skyrmion escape through a free bondary of the (100 x 100 lattice): (left) full MEP, (right) truncated path (TMEP) with anchor at ferromagnetic state (last image). (top) first images of MEP after 2000 iterations, (middle) MEP after 2000 iterations, (bottom) convergence history. Solid black line depicts energy, dash-dotted line is l™ norm of gradient, and dashed line is l™ norm of the gradient component perpendicular to the path. Down-pointing triangle, circle and up-pointing triangle markers show descending, climbing and holding images, respectively. Images with small markers are subjects of reparametrization step. All images of TMEP are closer to the boundary than those of the full MEP. Resolution in the part of interest is higher, convergence is faster and precision is higher for TMEP as compared with the full MEP calculation
3.3. Skyrmion escape
As a second application, the SP for a skyrmion escaping through free boundaries of the simulated system (periodic boundaries turned off). The initial path is generated from a skyrmion configuration having a radius of 30 lattice constants with center gradually moving from one radius toward center of the domain from the interface to one radius outside the domain. In the calculation of the full MEP, the endpoints of the path were relaxed, but for TMEP calculation the ends were not relaxed. The beginning of TMEP is relaxed to the single skyrmion state; since skyrmions are repelled by the boundary, the isolated skyrmion is located at the center of the simulated system. The final configuration is FM (all spins are directed along z axis except of spins near the boundary). The optimization methods and parameters were the same as in the previous application.
Fig. 3. MEP for skyrmion annihilation (1000 x 1000 lattice): (left) full MEP, (right) truncated path (TMEP) with anchor at ferromagnetic state (last image). (top) first images of MEP after 2000 iterations, (middle) MEP after 2000 iterations, (bottom) convergence history
The calculated results are shown in Fig. 2. Since the lowest energy configuration of the skyrmion is at the center of the simulated system, the first part of the MEP describes translation of the skyrmion with only a slight variation in the energy. This first part is of little interest, but can be long if the domain is large. Hence, most images of the full MEP belong to this first part. The second part of the MEP describes the escape of the skyrmion through one of the edges and is of main interest. In the TMEP calculation the focus is placed in the second part while the first part is not included. In this benchmark it is hard to estimate the length of the TMEP from the beginning, so the anchor is set to the image with maximal energy. TMEP again gives good resolution of the important part of the MEP. The precision of the calculation of the SP is several times higher using TMEP when the same number of images is used to represent the full MEP. The first image on the path is most distant from the interface, for the TMEP the image is located in the region about the interface of the width about three skyrmion radii, when for MEP the image occupies the region about half of the simulation domain. As before, the TMEP can be used to significantly decrease computational effort.
For comparison, the calculations described above are repeated for a denser lattice consisting of 1000 x 1000 atoms. All parameters are scaled accordingly to make the skyrmion have the same relative dimensions to the domain size as in the previous example. For a small number of images, the CI in this case does not converge on the SP (the gradient does not become zero), although the projected gradient at the CI becomes zero. This is a result of inaccuracy in the tangent estimate. To improve the tangent estimate, the number of images was increased
Fig. 4. MEP for skyrmion escape through a free boundary of the (1000 x 1000 lattice): (left) full MEP, (right) truncated path (TMEP) with anchor at ferromagnetic state (last image). (top) first images of MEP after 2000 iterations, (middle) MEP after 2000 iterations, (bottom) convergence history
to 32 and that enabled convergence of the CI to the SP. All optimization parameters are the same as in the previous example. The calculated results are shown in Fig. 3 and Fig. 4.
4. Acknowledgements
This work was supported by the Government of the Russian Federation (grant 074U01) and by grant 16-1110330 of Russian Science Foundation.
References
[1] Mills G., Jonsson H., Schenter G.K. Reversible work based transition state theory: Application to H2 dissociative adsorption. Surf. Sci., 1995, 324, P. 305-337.
[2] Jonsson H., Mills G., Jacobsen K.W., edited by Berne B.J., Ciccotti G., Coker D.F. Nudged elastic band method for finding minimum energy paths of transitions, in Classical and Quantum Dynamics in Condensed Phase Simulations. World Scientific, Singapore, 1998, P. 385-404.
[3] Henkelman G., Arnaldsson A., Jonsson H. Theoretical calculations of CH4 and H2 associative desorption from Ni(111): Could subsurface hydrogen play an important role? J. Chem. Phys., 2006, 124, P. 044706(9).
[4] Henkelman G., Jonsson H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys., 1999, 111(15), P. 7010-7022.
[5] Olsen R.A., Kroes G.J., Henkelman G., Arnaldsson A., Jonsson H. Comparison of methods for finding saddle points without knowledge of the final states. J. Chem. Phys., 2004, 121(20), P. 9776-9792.
[6] Gutierrez M.P., Argaez C., Jonsson H. Improved minimum mode following method for finding first order saddle points. J. Chem. Theo. Comput., 2017, 13(1), P. 125-134.
[7] Henkelman G., Johannesson G., Jonsson H. Methods for finding saddle points and minimum energy paths. in Theoretical Methods in Condensed Phase Chemistry, 2002, P. 269-302.
[8] Henkelman G., Uberuaga B.P., Jonsson H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys., 2000, 113(22), P. 9901-9904.
[9] Maragakis P., Andreev S.A., Brumer Y., Reichman D.R., Kaxiras E. Adaptive nudged elastic band approach for transition state calculation. J. Chem. Phys. 2002, 117, P. 4651-4658.
[10] Zhu T., Li J., Samanta A., Kim H.G., Suresh S. Interfacial plasticity governs strain rate sensitivity and ductility in nanostructured metals. PNAS, 2007, 104 (9), P. 3031-3036.
[11] Einarsdottir D.M., Arnaldsson A., Oskarsson F. and Jonsson H. Path optimization with application to tunneling. Lecture Notes in Computer Science, 2012, 7134, P. 45-55.
[12] Henkelman G., Jonsson H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys., 2000, 113(22), P. 9978-9985.
[13] Weinan E., Weiqing R., Vanden-Eijnden E. String method for the study of rare events. Phys. Rev. B, 2002, 66, P. 052301(4).
[14] Ravishenker, et al. Reviews in Computational Chemistry, Volume 11, edited by Kenny B. Lipkowitz, Donald B. Boyd, John Wiley & Sons, Sep 22, 2009, chapter 6, P. 346-347.
[15] Sheppard D., Terrell R., Henkelman G. Optimization methods for finding minimum energy paths. J. Chem. Phys., 2008,128, P. 134106(10).
[16] Bessarab P.F., Uzdin V.M. Jonsson H. Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation. Comp. Phys. Commun., 2015, 196, P. 335-347.
[17] Bessarab P.F. Comment on "Path to collapse for an isolated Neel skyrmion". Phys. Rev. B , 2017, 95(13), P. 136401(2).
[18] Malottki S.V., Dupé B., Bessarab P.F., Delin A., Heinze S. Enhanced skyrmion stability due to exchange frustration. Sci. Rep., 2017, 7, P. 12299(10).
