Xf_ThermSpinXferEvolve

— 17 Comments

Download

The Xf_ThermSpinXferEvolve OOMMF extension is a combination of the Oxs_SpinXferEvolve evolver that comes default with OOMMF, and the thermal fluctuation model that is in the UHH_ThetaEvolve evolver. The Xf_ThermSpinXferEvolve extension was created to allow the simultaneous simulation of spin-transfer torque and thermal fluctuation effects on a magnetic ensemble.

The spin torque equation derived by Sloczewski [1] and modified by Xiao [2] is written as (as implemented in OOMMF):

$\frac{\partial\textbf{m}}{\partial t} = -|\gamma|\textbf{m}\times\textbf{H}_\text{EFF}+\alpha\textbf{m}\times\frac{\partial\textbf{m}}{\partial t}+\textbf{STT}$
$\textbf{STT} = \left|\gamma\right|\beta \left( \epsilon \left( \textbf{m}\times \textbf{m}_{p} \times \textbf{m} \right) - \epsilon'\left(\textbf{m}\times\textbf{m}_{p}\right) \right)$

According to the OOMMF user guide:

$\textbf{m}=\text{reduced magnetization,} \frac{\textbf{M}}{M_{sat}}$
$\gamma=\text{Gilbert gyromagnetic ratio}$
$\beta=\left|\frac{\hbar}{\mu_{0}e}\right|\frac{J}{tM_{sat}}$
$\textbf{m}_{p}=\text{(unit) electron polarization direction}$

$\gamma$ has a default value of $-2.211\times{}10^{5} \text{m/A}\cdot\text{s}$ . In the definition of $\beta$ , $\hbar$ is the reduced Planck’s constant, $e$ is the electronic charge in C, $J$ is the current density exerting spin-torque in $\text{A/m}^{2}$ , $t$ is the thickness of the free layer in meters in the direction which the current density is flowing, and $M_{sat}$ is the saturation magnetization in $\text{A/m}$ . Note that $\beta$ may be rewritten as,

$\beta=\left|\frac{\hbar}{\mu_{0}e}\right|\frac{I}{M_{sat}Vol}$

where $I$ is the current flowing homogeneously into the magnetic cell in A, and $Vol$ is the volume of the magnetic cell in $\text{m}^{3}$ .

$\epsilon$ and $\epsilon'$ gives the in-plane and out-of-plane spin-torque terms, respectively. As implemented in OOMMF,

$\epsilon=\frac{q_{+}}{A_{+}+A_{-}\left(\textbf{m}\cdot\textbf{m}_{p}\right)}+\frac{q_{-}}{A_{+}-A_{-}\left(\textbf{m}\cdot\textbf{m}_{p}\right)}$
$q_{\pm}=P_{fixed}\Lambda_{fixed}^{2}\sqrt{\frac{\Lambda_{free}^{2}+1}{\Lambda_{fixed}^{2}+1}}\pm P_{free}\Lambda_{free}^{2}\sqrt{\frac{\Lambda_{fixed}^{2}-1}{\Lambda_{free}^{2}-1}}$
$A_{\pm}=\sqrt{\left(\Lambda_{fixed}^{2}\pm 1\right)\left(\Lambda_{free}^{2}\pm 1\right)}$

In the case where $P_{fixed}=P_{free}$ and $\Lambda_{fixed}=\Lambda_{free}$ , $\epsilon$ reduces to

$\epsilon=\frac{P\Lambda^{2}}{\left(\Lambda^{2}+1\right)+\left(\Lambda^{2}-1\right)\left(\textbf{m}\cdot\textbf{m}_{p}\right)}$

The Specify block for the Xf_ThermSpinXferEvolve class has the form

$\text{Specify Xf\_ThermSpinXferEvolve:}\textit{name}~\{$
     $\text{alpha}~~~~~~~~~~~~\alpha$
     $\text{gamma\_LL}~~\bar{\gamma}$
     $\text{gamma\_G}~~~~\gamma$
     $\text{do\_precess}~~~~\textit{precess}$
     $\text{allow\_signed\_gamma}~~~~\textit{signed\_gamma}$
     $\text{min\_timestep}~~~~\textit{minimum\_stepsize}$
     $\text{min\_timestep}~~~~\textit{maximum\_stepsize}$
     $\text{fixed\_spins}~~~~\{$
             $\textit{atlas\_spec}$
             $\textit{region1 region2 }\cdots$
     $\}$

     $\text{start\_dm}~~~~\Delta\textbf{m}$
     $\text{start\_dt}~~~~\text{start\_timestep}$
     $\text{stage\_start}~~~~\textit{scontinuity}$
     $\text{error\_rate}~~~~\textit{rate}$
     $\text{absolute\_step\_error}~~~~\textit{abs\_error}$
     $\text{relative\_step\_error}~~~~\textit{rel\_error}$
     $\text{error\_precision}~~~~\textit{eprecision}$
     $\text{min\_step\_headroom}~~~~\textit{min\_headroom}$
     $\text{max\_step\_headroom}~~~~\textit{max\_headroom}$
     $\text{reject\_goal}~~~~\textit{reject\_proportion}$
     $\text{method}~~~~\textit{subtype}$
     $\text{P}~~~~\textit{polarization}$
     $\text{P\_fixed}~~~~\textit{p\_fixed\_layer}$
     $\text{P\_free}~~~~\textit{p\_free\_layer}$
     $\text{Lambda}~~~~\Lambda$
     $\text{Lambda\_fixed}~~~~\Lambda\textit{\_fixed\_layer}$
     $\text{Lambda\_free}~~~~\Lambda\textit{\_free\_layer}$
     $\text{eps\_prime}~~~~\textit{ep}$
     $\text{J}~~~~\textit{current\_density}$
     $\text{J\_direction}~~~~\textit{current\_flow\_direction}$
     $\text{J\_profile}~~~~\textit{Jprofile\_script}$
     $\text{J\_profile\_args}~~~~\textit{Jprofile\_script\_args}$
     $\text{mp}~~~~\textit{p\_direction}$
     $\text{propagate\_mp}~~~~\textit{prop\_mp}$
     $\text{temperature}~~~~\textit{cell\_temperature}$
     $\text{tempscript}~~~~\textit{temperature\_profile\_script}$
     $\text{tempscript\_args}~~~~\textit{temperature\_profile\_script\_arguments}$
     $\text{uniform\_seed}~~~~\textit{uniform\_seed}$
$\}$

Most of these options appear also in the Oxs_EulerEvolve class. The repeats have the same meaning as in that class, and the same default values except for relative_step_error and error_rate, which for Xf_ThermSpinXferEvolve have the default values of 0.01 and -1.0, respectively. Additionally, the alpha, gamma_LL and gamma_G options may be initialized using scalar field objects, to allow these material parameters to vary spatially.

The default values for P and Lambda are 0.4 and 2, respectively. If preferred, values for the fixed and free layers may be instead specified separately, through P_fixed, P_free, Lambda_fixed, and Lambda_free. Otherwise P_fixed = P_free = P and
Lambda_fixed = Lambda_free = Lambda. Lambda must be larger than or equal to 1; set Lambda=1 to remove the dependence of $\epsilon$ on $\textbf{m}\cdot\textbf{m}_{p}$ . If you want non-zero $\epsilon'$ , it is set directly as eps_prime.

Current density J and unit polarization direction mp are required. The units on J are $\text{A/m}^{2}$ . Positive J produces torque that tends to align $\textbf{m}$ towards $\textbf{m}_{p}$ .

Simulation of domain-wall dynamics under current-induced spin-torque is enabled by setting propogate_mp to 1. The setting propogate_mp is 0 (disabled) by default. When propogate_mp is enabled, mp is actually $\Delta_{x}\times\frac{\partial\textbf{m}}{\partial x}$ , where $x$ is the flow direction and $\Delta_{x}$ is the cell dimension in that direction. The flow direction may be set by setting J_direction as one of six options:
-z, +z, -y, +y, -x, +x. The default is -z. The direction changes the mp used to calculate the spin torque at each cell site.

Parameters J, mp, P, Lambda, and eps_prime may all be varied pointwise (by specifying them as Oxs_ScalarField objects), but are fixed with respect to time. However, J can be multiplied by a time varying “profile,” to model current rise times, pulses, etc. Use the J_profile and J_profile_args options to enable this feature. The Jprofile_script should be a Tcl script that returns a single scalar. Jprofile_script_args should be a subset of {stage stage_time total_time }, to specify arguments appended to Jprofile_script on each time step. Default is the entire set, in the order as listed.

Simulations with thermal fluctuations are enabled using the parameters temperature (in units of Kelvin). The parameter temperature sets a constant temperature of thermal fluctuations. Heating effects may be simulated using the tempscript, which behaves similarly to Jprofile, to provide a time dependent scaling of temperature. Only spatially uniform temperature is allowed in the simulation, and temperature should be given as a scalar constant. The Tcl script given to tempscript should also return a scalar. The arguments to tempscript may be given as a list defined through tempscript_args. The arguments in tempscript_args should come from the following list: stage, stage_time, total_time.

The parameter uniform_seed sets up the random number generator used to generate the thermal fluctuation field for thermal simulations. (Developmental note: thermal simulations seem to have random errors if this is not set during simulations with non-zero temperature.)

The allow_signed_gamma parameter is for simulation testing purposes, and is intended for advanced use only. There is some lack of consistency in the literature with respect to the sign of $\gamma$ . For this reason the Landau-Lifshitz-Gilbert equations are presented above using the absolute value of $\gamma$ . This is the interpretation used if allow_signed_gamma is 0 (the default). If instead allow_signed_gamma is set to 1, then the Landau-Lifshitz-Gilbert equations are interpreted without the absolute values and with a sign change on the $\gamma$ terms, i.e., the default value for $\gamma$ in this case is $-2.211\times 10^{5}$ (units are $\text{m/A}\cdot\text{s}$ ). In this setting, if $\gamma$ is set positive then the spins will precess backwards about the effective field, and the damping term will force the spins away from the effective field and increase the total energy. If you are experimenting with $\gamma > 0$ , you should either set $\alpha \leq 0$ to force spins back towards the effective field, or disable the energy precision control (discussed below).

The two controls min_step_headroom (default value 0.33) and max_step_headroom (default value 0.95) replace the single step_headroom option in Oxs_EulerEvolve. The effective step_headroom is automatically adjusted by the evolver between the min_headroom and max_headroom limits to make the observed reject proportion approach the reject_goal (default value 0.05).

The method entry selects a particular Runge-Kutta implementation. It should be set to one of rk2, rk2heun, rk4, rkf54, rkf54m, or rkf54s; the default value is rkf54. The rk2 and rk4 methods implement canonical second and fourth global order Runge-Kutta methods [1], respectively. For rk2, stepsize control is managed by comparing $\dot{{\textbf{m}}}$ at the middle and final points of the interval, similar to what is done for stepsize control for the Oxs_EulerEvolve class. One step of the rk2 method involves 2 evaluations of $\dot{{\textbf{m}}}$ . The rk2heun method implements Heun’s method and is essentially a modified version of the forward Euler method. Heun’s method calculates $\dot{{\textbf{m}}}$ at the initial and the predict step, and uses the average of the two as the actual $\dot{{\textbf{m}}}$ for calculating the next step.

In the rk4 method, two successive steps are taken at half the nominal step size, and the difference between that end point and that obtained with one full size step are compared. The error is estimated at 1/15th the maximum difference between these two states. One step of the rk4 method involves 11 evaluations of $\dot{{\textbf{m}}}$ , but the end result is that of the 2 half-sized steps.

The remaining methods, rkf54, rkf54m, and rkf54s, are closely related Runge-Kutta-Fehlberg methods derived by Dormand and Prince [2, 3]. In the nomenclature of these papers, rkf54 implements RK5(4)7FC, rkf54m implements RK5(4)7FM, and rkf54s implements RK5(4)7FS. All are 5th global order with an embedded 4th order method for stepsize control. Each step of these methods requires 6 evaluations of $\dot{{\textbf{m}}}$ if the step is accepted, 7 if rejected. The difference between the methods involves tradeoffs between stability and error minimization. The RK5(4)7FS method has the best stability, RK5(4)7FM the smallest error, and RK5(4)7FC represents a compromise between the two. The default method used by Oxs_RungeKuttaEvolve is RK5(4)7FC.

The remaining undiscussed entry in the Xf_ThermSpinXferEvolve Specify block is energy_precision. This should be set to an estimate of the expected relative accuracy of the energy calculation. After accounting for any change in the total energy arising from time-varying applied fields, the energy remainder should decrease from one step of the LLG ODE to the next. Xf_ThermSpinXferEvolve will reject a step if the energy remainder is found to increase by more than that allowed by eprecision. The default value for eprecision is 1e-10. This control may be disabled by setting eprecision to -1.

The Xf_ThermSpinXferEvolve module provides the same five scalar outputs and three vector outputs as Oxs_RungeKuttaEvolve, plus the scalar outputs “average J” and “Temperature,” and the vector field outputs “Spin torque” (which is $|\gamma|\beta\epsilon\left(\textbf{m}\times\textbf{m}_{p}\times\textbf{m}\right)$ ) and “J*mp.”

[1] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis (Springer, New York, 1993), 2nd edn.
[2] J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comp. Appl. Math., 6, 19-26 (1980)
[3] J. R. Dormand and P. J. Prince, “A reconsideration of some embedded Runge-Kutta formulae,” J. Comp. Appl. Math., 15, 203-211 (1986).

17 comments on “Xf_ThermSpinXferEvolve”

Jinjun Ding says:

June 7, 2014 at 12:27 am

Thank you very much for your sharing your source code. I’ve got one question about your code. And it lies in line 47 in xf_thermspinxferevolve_example1.mif. You set current_area as $width*$width. Why is it not $width*$lenth which means that the current is vertically injected, or $width*$thick which means that the current is injected in plane? I don’t quite understand the reason you set current_area as $width*$width.
Waiting for you reply.
Thank you again for your sharing your code.

Jinjun Ding
State Key Laboratory for Magnetism
Institute of Physics
Chinese Academy of Sciences
Beijing 100190, P.R. China
Phone: 0086-10-82648195
E-mail: jinjun.ding@gmail.com

LikeLike

Reply
xfong says:

June 8, 2014 at 2:38 pm

Thanks for trying out the source code Jinjun. The examples simulate the injection of current vertically into only a portion of the magnet. In the examples, current is only injected into a cross-sectional area that is $width*$width.

Regards,
Xuanyao (Kelvin) Fong

LikeLike

Reply
- Patrick says:
  
  October 1, 2015 at 8:41 pm
  
  Thank you very much for sharing this source code and I would like to follow up the question from Jinjun.
  
  As you have said, $width*$width is used to inject a current only into that particular cross-sectional area. If I wish to change the direction of current from vertically into in-plane and so on, how would you go by doing that?
  
  Thank you again for the code and Regards,
  Patrick
  
  LikeLike
  
  Reply
Brahim Messaoudi says:

September 24, 2015 at 9:36 pm

Could you please tell me whether you have had so far any experience in using temperature scripts together with your source code in calculations involving temperature gradients?. If so, would you please share your experience in terms of how did you structure your temperature script and so on…
I look forward to hearing from you.
Brahim@MERL
Cambridge MA

LikeLike

Reply
user says:

January 17, 2019 at 7:07 am

Hi there,

How would you distinguish STT and SOT in this module?
Thanks.

LikeLike

Reply
- Xuanyao (Kelvin) Fong says:
  
  August 8, 2019 at 2:56 am
  
  STT and SOT are not distinguished in this module. All the Xf_STT allows you to do is to include different spin torques (either STT or SOT) into the simulation (for example by injecting current through both the heavy metal and the MTJ in an SOT MTJ structure).
  
  This is achieved by mathematically rewriting expression of spin torque into a form that looks like a Zeeman field. Since Xf_STT is an energy class, it was written to not contribute to the energy of the system. OOMMF will include as many energy classes as declared and thus, the STT and SOT terms can be included as separate Xf_STT modules.
  
  LikeLike
  
  Reply
Yen-Lin says:

January 7, 2021 at 6:39 am

Hi Kevin,

Thank you very much for sharing this module, it has been very useful for me.
I have one question regarding the example:

proc dualPinnedLayers { xrel yrel zrel } {
if {$xrel 0.5, mp pointing up). Do you have any particular reason to set the mp has this spin polarization?

LikeLike

Reply
- Xuanyao (Kelvin) Fong says:
  
  January 9, 2021 at 9:46 am
  
  This example file is associated with the paper doi:10.1109/IMW.2013.6582105
  
  LikeLike
  
  Reply
  - Yen-Lin says:
    
    February 2, 2021 at 6:46 am
    
    Great thanks!
    
    LikeLike
Stanislav says:

November 24, 2021 at 11:10 am

Good day dear Kevin.
Thank you very much for this module. Can I ask you a question regarding your examples?
I tested your examples, and I think your code seriously overestimates thermal fluctuations.
Even in the absence of external torque, it produces large deviations of the magnetization vector. I checked your source code, and it seems I know the problem.
In your code, you calculated the standard deviation of the distribution like sqrt((2*alpha*kb*T)/(mu0^2*gamma*Ms*V*dt)). In other words, it is proportional to 1/sqrt(dt). As far as I know, the conventional model of the Brownian motion (or Weiner process) tells that the standard deviation is proportional to sqrt(dt). Do you have any particular reason to add the dt into the denominator?

LikeLike

Reply
- Xuanyao (Kelvin) Fong says:
  
  November 24, 2021 at 3:01 pm
  
  Hi Stanislav,
  
  I think you might be confused between dm and dm/dt. The formulation says dm is proportional to sqrt(dt) and hence, dm/dt is proportional to 1/sqrt(dt). When the time integrator calculates m(t) from the initial value problem, it will multiply dm/dt with dt to solve for dm at every step, starting at m(t=0). The multiplication will result in dm being proportional to sqrt(dt) which is what you expect.
  
  LikeLike
  
  Reply
  - Stanislav says:
    
    November 25, 2021 at 12:22 am
    
    Thank you very much for your kind reply.
    Excuse me, but I meant thermal field Hth, not dm/dt. As far as I know, the thermal fluctuations are represented by field Hth. And Hth should be proportional to sqrt(dt), but not dm/dt.The term sqrt((2*alpha*kb*T*dt)/(mu0^2*gamma*Ms*V))=dHth; therefore, Hth=dHth+Hth0. After calculations of all Fields, they will be integrated.
    In the absence of the external torques, the magnetization vector will move towards the new random position of the Heff=Hanis+Hdemag+Hexch+Hth.
    
    LikeLike
  - Xuanyao (Kelvin) Fong says:
    
    November 25, 2021 at 1:09 am
    
    Hth cannot be proportional to sqrt(dt)
    
    If it is, integrating over time will yield the “noise power” of Hth having a time dependence. This cannot be true. The only possibility is for Hth to have Hth to be proportional to 1/sqrt(dt)
    
    LikeLike
Stanislav says:

November 25, 2021 at 2:19 am

Excuse me for my mistake.
dHth is proportional to sqrt(dt), not Hth. Hth modeled using the stochastic Wiener process. Such that dHth=Hth-Hth0=N(0,1)*(sigma^2). Where N(0,1) is a random number with a normal distribution. Sigma^2 is proportional to sqrt(dt). In the Wiener process, increments are independent of the previous values. However, the motion itself is continuous.
If you don’t mind, can we discuss it in more detail via e-mail?
Thank you very much.

LikeLike

Reply
- Xuanyao (Kelvin) Fong says:
  
  November 25, 2021 at 3:27 am
  
  Sure.
  kelvin.xy.fong@nus.edu.sg
  
  LikeLike
  
  Reply
HARIKRISHNAN O K says:

May 12, 2022 at 5:44 pm

Hi Kevin,
I am trying to do SOT simulation using this module.But actually, I am suffering from an error which given as follows.

Error thrown from inside “Oxs_Director::GetAllScalarOutputs” — Error evaluating output “Max dm/dt” of object “Oxs_SpinXferEvolve:” — Oxs_Ext ERROR in object Oxs_SpinXferEvolve:: Oxs_SpinXferEvolve::UpdateDerivedOutputs: Can’t derive Delta E from single state.
(Oxs_Output)

————–

LikeLike

Reply
- Xuanyao (Kelvin) Fong says:
  
  May 22, 2022 at 12:38 pm
  
  Hi Harikrishnan,
  
  I recommend that you try using the 1.2a4 version of OOMMF with the module first. If the problems are still there, there is probably a compatibility issue with the module.
  
  Best,
  Kelvin FONG
  
  LikeLike
  
  Reply