2.3.9 Voter method^2.6

The Voter method (also called ``Hyperdynamics of the infrequent events'') was proposed in Refs. [Voter 97b,Voter 97a] and used to calculate the diffusional processes of Ag atoms on a Ag(111) surface, achieving an acceleration of the calculation up to $8000$ times. For the first time we have applied the method to magnetic system and that is the motivation of the work presented in this section.

The method consists in modification of the external potential, basing on the Hessian energy matrix $\partial^2E/\partial x_i\partial x_j$, where $x_i$ are the system coordinates, so that the transition state (saddle point) remains unchanged. An additional external boost potential, $\Delta V_b$ (see Fig. 2.21) is slowly switched on at the minimum, rising its value, and is switched off near the transition surface, i.e. where the first eigenvalue of the Hessian matrix $\lambda_1$ becomes negative. This boost term is always positive making easier to escape from the potential minimum. The Langevin dynamics is then performed in this modified potential. The total time for the escape of the particle from the minimum can be evaluated as the sum of modified times at each timestep $\Delta t_b$, which could be computed from the Langevin dynamics timestep in the modified potential, $\Delta t_{LD}$, as the following:

$$\Delta t_b= \Delta t_{LD} \exp(\Delta V_b(t)/k_BT).$$ (2.64)

We have implemented the method for collection of non-interacting magnetic particles with external field applied at some angle to their anisotropy axis. According to the A.Voter's suggestion, we have tried two forms of the boost potentials to accelerate the stochastic dynamical calculations:

$$\Delta V_{b1}=a\theta(\lambda_1)\lambda_1^2,$$ (2.65)

$$\Delta V_{b2}=b\frac{a\theta(\lambda_1)\lambda_1^2}{1+2a\theta(\lambda_1)\lambda_1^2}.$$ (2.66)

Here $\theta$ is the standard Heaviside function and $a$ and $b$ are arbitrary parameters which could be tuned. On one hand, these parameters must be as large as possible to achieve the acceleration. On the other hand, the tuning of these parameters could be performed based on the fact that the equilibrium statistics should be achieved in the minimum, i.e. the number of LD time steps performed in the modified potential before the particle is escaped from the minimum must be large enough. The difference between the two potentials is that the second one, although more difficult to evaluate, produces much smoother modification and does not introduce a ``crest'' in the minimum, from which the particle could be scattered (see Fig. 2.21).

Figure: The initial potential for a magnetic moment, formula (2.43), and two possible boosting potentials (2.65) and (2.66), as function of magnetic moment angle $\theta$ and in the absence of external applied field.

The implementation of the method requires the constant evaluation of the lowest eigenvalue of the correspondent Hessian matrix. The direct normal mode analysis is a time consuming procedure which depends strongly on the system size and would limit the acceleration achieved by the method. The direct evaluation of the eigenvalue problem scales like $N^3$, being $N$ the number of moments in the system, whereas the number of operations in a Langevin dynamics steps scales like $N$. In order to obtain a good performance of the method, advanced methods that scale like $N$ are desirable. The iterative methods, such as the Gauss-Siegel, which could use the previous value as initial guess, could be very helpful. A.Voter [Voter 97a] also suggested to replace the direct evaluation of $\lambda_1$ by its approximate evaluation by means of the numerical minimization (with respect to the parameter s) of the following expression:

$$\lambda^{num}(\mathbf{s})=[E(\mathbf{x}+\eta \mathbf{s})+E(\mathbf{x}-\eta \mathbf{s})-2E(\mathbf{x})]/\eta^2$$ (2.67)

where $\eta$ is a small parameter.

Figure 2.22: Switching time obtained with direct LLG calculations and with the Voter's method for a magnetic particle at zero applied field.

The averaged time for the particle to escape from the minimum is calculated using the A.Voter method and compared to that obtained from the direct integration of the LLG equation with a random term representing temperature fluctuations. Fig. 2.22 presents results of the calculations for switching time for an ensemble of uniaxial particles averaged over many realizations. The computation is stopped when the standard deviation from the average value is below $1\%$. It is clear that Voter's method for reasonable computational time is much faster than the direct LLG integration. Remarkably, the method reproduces correctly all the features of the dynamics, including the precession. Fig. 2.23 presents the histogram for switching time of a particle for energy barrier value $KV/k_BT=4$. It shows that the accelerated dynamics correctly reproduces the form of a $t \exp (-t/\tau)$ distribution, where $t$ is the reversal time, although the general tendency of the Voter distribution is the displacement to larger values.

Figure: Histogram showing the distribution of reversal times for both LLG calculations and Voter's method with energy barrier value $KV/k_BT=4$.

We have checked the results for different values of the damping parameters and different angles between the applied field and the anisotropy directions. Fig. 2.24 presents the results for an angle between the anisotropy direction and applied field of $45$ degrees and for different values of the tuning parameter $a$. Therefore, unlike the TQMC method (see discussion in Section 2.3.5), the hyperdynamics method correctly reproduces the influence of the ellipticity of the precessional cone on the thermal switching statistics.

Figure: Average switching time for a magnetic moment with an applied field $(0.2 K/M_s, 0.2 K/M_s, 0)$ and the easy axis parallel to $z$ direction. Results for different values of parameter $a$ are represented, as well as LLG direct calculations.

The real acceleration of the method depends on the efficiency to calculate the lowest eigenvalue of a complex large system. Therefore, for the small barriers case, the direct integration of the LLG equation is faster. To compare, we present in Fig. 2.25 the ratio between the average CPU time used in the Voter method (using direct lowest eigenvalue evaluation in system described by Eq. (2.67)) and the average CPU time used in the LD dynamics. The acceleration in calculation appears for barrier values larger than $KV/k_BT = 7$. For this ``straightforward'' implementation, the acceleration up to $24$ times in CPU time has been reached. More sophisticated methods will improve this ratio. The method was also checked for a linear chain of 16 exchange coupled magnetic moments with open boundary condition without dipolar interactions. In this case the Hessian is a banded matrix, being a full matrix with dipolar interactions. The result is shown in Fig. 2.26. The disadvantages of a direct evaluation of the Hessian eigenvalue problem matrix become evident even in this relatively simple system.

Figure: Comparison of real time performance: average total CPU time to produce reversal in LLG calculations divided by the same quantity using Voter's method. The line $y=1$ represents equal performance and greater number reflects better efficiency of the tested method.

Figure: Switching time obtained with direct LLG calculations and with the Voter's method for a linear chain of 16 exchange coupled magnetic moments at zero applied field and parameters $J/KV= 0.4$ and $\alpha = 0.3$

In conclusion, we have implemented the method of ``Hyperdynamics of infrequent events'' to accelerate the molecular dynamics simulations in the case of magnetization dynamics achieving an acceleration up to $24$ times. Higher acceleration seems also possible if one uses more sophisticated modern methods to evaluate rapidly the lowest eigenvalue of the Hessian matrix. In comparison to the time-quantified Monte Carlo, the main advantage of the method is the correct description of the influence of the precession on the thermal switching process. In contrast to the Victora method this method does not suppose a priori the Arrhenius-Néel law with one unique temperature-independent attempt frequency. It may successfully be used when various reversal modes, with different attempt frequencies, coexist during the thermal magnetization process. The limitations of the method make it useful for intermediate timescale, up to hundreds of nanoseconds, for example, for the dynamic coercivity calculations. Higher time scale seems not reachable.