Smit-Beljers model for ferromagnetic resonance¶

In this notebook you can see how to run a FMR calculation with a theoretical Smit-Beljers model.

Compared to standard CMTJ here we work in Ms of the A/m units instead of T, and also spherical coordinates. Use VectorObj reference to create objects in convenient frame of reference for you.
Here's the reference

In [1]:

import numpy as np

from collections import defaultdict
from cmtj.models import LayerSB, VectorObj, Solver
from cmtj.utils import mu0
from tqdm import tqdm

Ms1 = (
    1.0 / mu0
)  # here we pass the saturation magnetisation in A/m, but in the dynamic model we use T!
Ms2 = 1.2 / mu0
layerA = LayerSB(
    _id=0,
    thickness=1e-9,
    Kv=VectorObj(
        np.deg2rad(0.0), np.deg2rad(0), 1e1
    ),  # for the Kv only phi angle counts !
    Ks=3e4,
    Ms=Ms1,
)
layerB = LayerSB(
    _id=1,
    thickness=1.3e-9,
    Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), 1e4),
    Ks=1e1,
    Ms=Ms2,
)

# we indicate the "guess" of the initial position
# it's generally good to align it with the field, but it's not necessary
current_position = [np.deg2rad(89), np.deg2rad(0.1), np.deg2rad(180), np.deg2rad(0.1)]
Hspace = np.linspace(-400e3, 400e3, 100)
result_dictionary = defaultdict(list)
# we perform a sweep over the field magnitude
for Hmag in tqdm(Hspace):
    solver = Solver(
        layers=[layerA, layerB],
        J1=[1e-4],
        J2=[0.0],
        H=VectorObj(np.deg2rad(89), np.deg2rad(0.1), Hmag),
    )
    # check for additional parameters in the solver
    # such as gradient convergence tolerance, max iterations, etc.
    # also in the solver there are root finding parameters
    (t1, p1, t2, p2), frequencies = solver.solve(init_position=current_position)
    # frequencies are already in GHz
    for frequency in frequencies:
        result_dictionary["frequency"].append(frequency)
        result_dictionary["Hmag"].append(Hmag)

    # we note the final position of the magnetisation in spherical coordinates
    result_dictionary["theta_1"].append(t1)
    result_dictionary["phi_1"].append(p1)
    result_dictionary["theta_2"].append(t2)
    result_dictionary["phi_2"].append(p2)
    # we reuse the previous solution as the initial guess for the next iteration
    current_position = [t1, p1, t2, p2]

import numpy as np from collections import defaultdict from cmtj.models import LayerSB, VectorObj, Solver from cmtj.utils import mu0 from tqdm import tqdm Ms1 = ( 1.0 / mu0 ) # here we pass the saturation magnetisation in A/m, but in the dynamic model we use T! Ms2 = 1.2 / mu0 layerA = LayerSB( _id=0, thickness=1e-9, Kv=VectorObj( np.deg2rad(0.0), np.deg2rad(0), 1e1 ), # for the Kv only phi angle counts ! Ks=3e4, Ms=Ms1, ) layerB = LayerSB( _id=1, thickness=1.3e-9, Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), 1e4), Ks=1e1, Ms=Ms2, ) # we indicate the "guess" of the initial position # it's generally good to align it with the field, but it's not necessary current_position = [np.deg2rad(89), np.deg2rad(0.1), np.deg2rad(180), np.deg2rad(0.1)] Hspace = np.linspace(-400e3, 400e3, 100) result_dictionary = defaultdict(list) # we perform a sweep over the field magnitude for Hmag in tqdm(Hspace): solver = Solver( layers=[layerA, layerB], J1=[1e-4], J2=[0.0], H=VectorObj(np.deg2rad(89), np.deg2rad(0.1), Hmag), ) # check for additional parameters in the solver # such as gradient convergence tolerance, max iterations, etc. # also in the solver there are root finding parameters (t1, p1, t2, p2), frequencies = solver.solve(init_position=current_position) # frequencies are already in GHz for frequency in frequencies: result_dictionary["frequency"].append(frequency) result_dictionary["Hmag"].append(Hmag) # we note the final position of the magnetisation in spherical coordinates result_dictionary["theta_1"].append(t1) result_dictionary["phi_1"].append(p1) result_dictionary["theta_2"].append(t2) result_dictionary["phi_2"].append(p2) # we reuse the previous solution as the initial guess for the next iteration current_position = [t1, p1, t2, p2]

100%|██████████| 100/100 [00:53<00:00,  1.85it/s]

In [2]:

import matplotlib.pyplot as plt

with plt.style.context(["science", "nature"]):
    w, h = plt.figaspect(0.5)
    fig, (ax1, ax2) = plt.subplots(1, 2, dpi=200, figsize=(w, h))
    Hvals = np.asarray(result_dictionary["Hmag"]) / 1e3
    ax1.plot(
        Hvals,
        result_dictionary["frequency"],
        "ro",
        label="frequency",
    )
    ax1.set_xlabel("H (kA/m)")
    ax1.set_ylabel("f (GHz)")

    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary["phi_1"]),
        color="forestgreen",
        label=r"$\phi_1$",
    )
    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary["theta_1"]),
        color="navy",
        label=r"$\theta_1$",
    )

    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary["phi_2"]),
        color="forestgreen",
        label=r"$\phi_2$",
        marker="o",
        markersize=1,
    )
    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary["theta_2"]),
        color="navy",
        label=r"$\theta_2$",
        marker="o",
        markersize=1,
    )
    ax2.legend()
    ax2.set_xlabel("H (kA/m)")
    ax2.set_ylabel("angle (deg)")
    fig.tight_layout()

import matplotlib.pyplot as plt with plt.style.context(["science", "nature"]): w, h = plt.figaspect(0.5) fig, (ax1, ax2) = plt.subplots(1, 2, dpi=200, figsize=(w, h)) Hvals = np.asarray(result_dictionary["Hmag"]) / 1e3 ax1.plot( Hvals, result_dictionary["frequency"], "ro", label="frequency", ) ax1.set_xlabel("H (kA/m)") ax1.set_ylabel("f (GHz)") ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary["phi_1"]), color="forestgreen", label=r"$\phi_1$", ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary["theta_1"]), color="navy", label=r"$\theta_1$", ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary["phi_2"]), color="forestgreen", label=r"$\phi_2$", marker="o", markersize=1, ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary["theta_2"]), color="navy", label=r"$\theta_2$", marker="o", markersize=1, ) ax2.legend() ax2.set_xlabel("H (kA/m)") ax2.set_ylabel("angle (deg)") fig.tight_layout()

Dynamic LLG from Smit-Beljers model, spherical model¶

Here, we can also return eigenvectors in addition to eigenvalues (frequencies). Under the hood the solver is a bit different, hence we use a slightly different input object.

In [3]:

import numpy as np

from collections import defaultdict
from cmtj.models import LayerDynamic, VectorObj, Solver
from cmtj.utils import mu0
from tqdm import tqdm

Ms1 = (
    1.0 / mu0
)  # here we pass the saturation magnetisation in A/m, but in the dynamic model we use T!
Ms2 = 1.2 / mu0
alpha = 1e-3
layerA = LayerDynamic(
    _id=0,
    thickness=1e-9,
    Kv=VectorObj(
        np.deg2rad(0.0), np.deg2rad(0), 1e1
    ),  # for the Kv only phi angle counts !
    Ks=3e4,
    Ms=Ms1,
    alpha=alpha,
)
layerB = LayerDynamic(
    _id=1,
    thickness=1.3e-9,
    Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), 1e4),
    Ks=1e1,
    Ms=Ms2,
    alpha=alpha,
)

# we indicate the "guess" of the initial position
# it's generally good to align it with the field, but it's not necessary
current_position = [np.deg2rad(89), np.deg2rad(0.1), np.deg2rad(180), np.deg2rad(0.1)]
Hspace = np.linspace(-400e3, 400e3, 100)
result_dictionary_dynamic = defaultdict(list)
# we perform a sweep over the field magnitude
for Hmag in tqdm(Hspace):
    solver = Solver(
        layers=[layerA, layerB],
        J1=[1e-4],
        J2=[0.0],
        H=VectorObj(np.deg2rad(89), np.deg2rad(0.1), Hmag),
    )
    # check for additional parameters in the solver
    # such as gradient convergence tolerance, max iterations, etc.
    # also in the solver there are root finding parameters
    # the last param here is an eigenvector
    (t1, p1, t2, p2), frequencies, _ = solver.solve(init_position=current_position)
    # frequencies are already in GHz
    for frequency in frequencies:
        result_dictionary_dynamic["frequency"].append(frequency)
        result_dictionary_dynamic["Hmag"].append(Hmag)

    # we note the final position of the magnetisation in spherical coordinates
    result_dictionary_dynamic["theta_1"].append(t1)
    result_dictionary_dynamic["phi_1"].append(p1)
    result_dictionary_dynamic["theta_2"].append(t2)
    result_dictionary_dynamic["phi_2"].append(p2)
    # we reuse the previous solution as the initial guess for the next iteration
    current_position = [t1, p1, t2, p2]

import numpy as np from collections import defaultdict from cmtj.models import LayerDynamic, VectorObj, Solver from cmtj.utils import mu0 from tqdm import tqdm Ms1 = ( 1.0 / mu0 ) # here we pass the saturation magnetisation in A/m, but in the dynamic model we use T! Ms2 = 1.2 / mu0 alpha = 1e-3 layerA = LayerDynamic( _id=0, thickness=1e-9, Kv=VectorObj( np.deg2rad(0.0), np.deg2rad(0), 1e1 ), # for the Kv only phi angle counts ! Ks=3e4, Ms=Ms1, alpha=alpha, ) layerB = LayerDynamic( _id=1, thickness=1.3e-9, Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), 1e4), Ks=1e1, Ms=Ms2, alpha=alpha, ) # we indicate the "guess" of the initial position # it's generally good to align it with the field, but it's not necessary current_position = [np.deg2rad(89), np.deg2rad(0.1), np.deg2rad(180), np.deg2rad(0.1)] Hspace = np.linspace(-400e3, 400e3, 100) result_dictionary_dynamic = defaultdict(list) # we perform a sweep over the field magnitude for Hmag in tqdm(Hspace): solver = Solver( layers=[layerA, layerB], J1=[1e-4], J2=[0.0], H=VectorObj(np.deg2rad(89), np.deg2rad(0.1), Hmag), ) # check for additional parameters in the solver # such as gradient convergence tolerance, max iterations, etc. # also in the solver there are root finding parameters # the last param here is an eigenvector (t1, p1, t2, p2), frequencies, _ = solver.solve(init_position=current_position) # frequencies are already in GHz for frequency in frequencies: result_dictionary_dynamic["frequency"].append(frequency) result_dictionary_dynamic["Hmag"].append(Hmag) # we note the final position of the magnetisation in spherical coordinates result_dictionary_dynamic["theta_1"].append(t1) result_dictionary_dynamic["phi_1"].append(p1) result_dictionary_dynamic["theta_2"].append(t2) result_dictionary_dynamic["phi_2"].append(p2) # we reuse the previous solution as the initial guess for the next iteration current_position = [t1, p1, t2, p2]

100%|██████████| 100/100 [00:25<00:00,  3.92it/s]

You see that we get the same result in the end -- even though the dynamic method seems to be a bit faster! There can be still some differences in the results due to the differences in the solver, let us know if and when that happens!

In [4]:

import matplotlib.pyplot as plt

with plt.style.context(["science", "nature"]):
    w, h = plt.figaspect(0.5)
    fig, (ax1, ax2) = plt.subplots(1, 2, dpi=200, figsize=(w, h))
    Hvals = np.asarray(result_dictionary_dynamic["Hmag"]) / 1e3
    ax1.plot(
        Hvals,
        result_dictionary_dynamic["frequency"],
        "bo",
        label="dynamic",
    )
    ax1.plot(
        Hvals,
        result_dictionary["frequency"],
        "ro",
        markersize=0.9,
        label="classical",
    )
    ax1.set_xlabel("H (kA/m)")
    ax1.set_ylabel("f (GHz)")

    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary_dynamic["phi_1"]),
        color="forestgreen",
        label=r"$\phi_1$",
    )
    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary_dynamic["theta_1"]),
        color="navy",
        label=r"$\theta_1$",
    )

    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary_dynamic["phi_2"]),
        color="forestgreen",
        label=r"$\phi_2$",
        marker="o",
        markersize=1,
    )
    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dictionary_dynamic["theta_2"]),
        color="navy",
        label=r"$\theta_2$",
        marker="o",
        markersize=1,
    )
    ax1.legend()
    ax2.legend()
    ax2.set_xlabel("H (kA/m)")
    ax2.set_ylabel("angle (deg)")
    fig.tight_layout()

import matplotlib.pyplot as plt with plt.style.context(["science", "nature"]): w, h = plt.figaspect(0.5) fig, (ax1, ax2) = plt.subplots(1, 2, dpi=200, figsize=(w, h)) Hvals = np.asarray(result_dictionary_dynamic["Hmag"]) / 1e3 ax1.plot( Hvals, result_dictionary_dynamic["frequency"], "bo", label="dynamic", ) ax1.plot( Hvals, result_dictionary["frequency"], "ro", markersize=0.9, label="classical", ) ax1.set_xlabel("H (kA/m)") ax1.set_ylabel("f (GHz)") ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary_dynamic["phi_1"]), color="forestgreen", label=r"$\phi_1$", ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary_dynamic["theta_1"]), color="navy", label=r"$\theta_1$", ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary_dynamic["phi_2"]), color="forestgreen", label=r"$\phi_2$", marker="o", markersize=1, ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dictionary_dynamic["theta_2"]), color="navy", label=r"$\theta_2$", marker="o", markersize=1, ) ax1.legend() ax2.legend() ax2.set_xlabel("H (kA/m)") ax2.set_ylabel("angle (deg)") fig.tight_layout()

Adding dipole interaction¶

You can also add dipole interaction to the model. This is done by setting the Ndipole parameter to a list of dipole tensors which are 3x3 matrices each. Each dipole tensor should describe the interaction between the corresponding layer and the next one. This is WIP feature and any feedback is appreciated.

In [5]:

import numpy as np

from collections import defaultdict
from cmtj.models import VectorObj, Solver
from cmtj.utils import mu0
from tqdm import tqdm


def compose_dipole_tensor(V, pos, other):
    r = np.asarray(pos) - np.asarray(other)  # vector from other to pos
    dr = np.linalg.norm(r)  # euclidean distance
    norm = V / (4 * np.pi * dr**5)  # normalization factor
    return [
        norm
        * VectorObj.from_cartesian(
            3 * r[0] ** 2 - dr**2, 3 * r[0] * r[1], 3 * r[0] * r[2]
        ),
        norm
        * VectorObj.from_cartesian(
            3 * r[0] * r[1], 3 * r[1] ** 2 - dr**2, 3 * r[1] * r[2]
        ),
        norm
        * VectorObj.from_cartesian(
            3 * r[0] * r[2], 3 * r[1] * r[2], 3 * r[2] ** 2 - dr**2
        ),
    ]


Ms1 = (
    1.0 / mu0
)  # here we pass the saturation magnetisation in A/m, but in the dynamic model we use T!
Ms2 = 1.2 / mu0
alpha = 1e-3
layerA = LayerSB(
    _id=0,
    thickness=1e-9,
    Kv=VectorObj(
        np.deg2rad(0.0), np.deg2rad(0), 1e1
    ),  # for the Kv only phi angle counts !
    Ks=3e4,
    Ms=Ms1,
)
layerB = LayerSB(
    _id=1,
    thickness=1.3e-9,
    Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), 1e4),
    Ks=1e1,
    Ms=Ms2,
)

# we indicate the "guess" of the initial position
# it's generally good to align it with the field, but it's not necessary
current_position = [np.deg2rad(89), np.deg2rad(0.1), np.deg2rad(180), np.deg2rad(0.1)]
Hspace = np.linspace(-400e3, 400e3, 30)
result_dict_dipole = defaultdict(list)
# we perform a sweep over the field magnitude
r = 300e-9
surf = r**2 * np.pi

dipoleA = compose_dipole_tensor(surf * 1e-9, [0, 0, 0], [100e-9, 0, 0])
for Hmag in tqdm(Hspace):
    solver = Solver(
        layers=[layerA, layerB],
        J1=[1e-4],
        J2=[0.0],
        H=VectorObj(np.deg2rad(89), np.deg2rad(0.1), Hmag),
        Ndipole=[dipoleA],
    )
    # check for additional parameters in the solver
    # such as gradient convergence tolerance, max iterations, etc.
    # also in the solver there are root finding parameters
    # the last param here is an eigenvector
    (t1, p1, t2, p2), frequencies = solver.solve(init_position=current_position)
    # frequencies are already in GHz
    for frequency in frequencies:
        result_dict_dipole["frequency"].append(frequency)
        result_dict_dipole["Hmag"].append(Hmag)

    # we note the final position of the magnetisation in spherical coordinates
    result_dict_dipole["theta_1"].append(t1)
    result_dict_dipole["phi_1"].append(p1)
    result_dict_dipole["theta_2"].append(t2)
    result_dict_dipole["phi_2"].append(p2)
    # we reuse the previous solution as the initial guess for the next iteration
    current_position = [t1, p1, t2, p2]

import numpy as np from collections import defaultdict from cmtj.models import VectorObj, Solver from cmtj.utils import mu0 from tqdm import tqdm def compose_dipole_tensor(V, pos, other): r = np.asarray(pos) - np.asarray(other) # vector from other to pos dr = np.linalg.norm(r) # euclidean distance norm = V / (4 * np.pi * dr**5) # normalization factor return [ norm * VectorObj.from_cartesian( 3 * r[0] ** 2 - dr**2, 3 * r[0] * r[1], 3 * r[0] * r[2] ), norm * VectorObj.from_cartesian( 3 * r[0] * r[1], 3 * r[1] ** 2 - dr**2, 3 * r[1] * r[2] ), norm * VectorObj.from_cartesian( 3 * r[0] * r[2], 3 * r[1] * r[2], 3 * r[2] ** 2 - dr**2 ), ] Ms1 = ( 1.0 / mu0 ) # here we pass the saturation magnetisation in A/m, but in the dynamic model we use T! Ms2 = 1.2 / mu0 alpha = 1e-3 layerA = LayerSB( _id=0, thickness=1e-9, Kv=VectorObj( np.deg2rad(0.0), np.deg2rad(0), 1e1 ), # for the Kv only phi angle counts ! Ks=3e4, Ms=Ms1, ) layerB = LayerSB( _id=1, thickness=1.3e-9, Kv=VectorObj(np.deg2rad(0.0), np.deg2rad(0), 1e4), Ks=1e1, Ms=Ms2, ) # we indicate the "guess" of the initial position # it's generally good to align it with the field, but it's not necessary current_position = [np.deg2rad(89), np.deg2rad(0.1), np.deg2rad(180), np.deg2rad(0.1)] Hspace = np.linspace(-400e3, 400e3, 30) result_dict_dipole = defaultdict(list) # we perform a sweep over the field magnitude r = 300e-9 surf = r**2 * np.pi dipoleA = compose_dipole_tensor(surf * 1e-9, [0, 0, 0], [100e-9, 0, 0]) for Hmag in tqdm(Hspace): solver = Solver( layers=[layerA, layerB], J1=[1e-4], J2=[0.0], H=VectorObj(np.deg2rad(89), np.deg2rad(0.1), Hmag), Ndipole=[dipoleA], ) # check for additional parameters in the solver # such as gradient convergence tolerance, max iterations, etc. # also in the solver there are root finding parameters # the last param here is an eigenvector (t1, p1, t2, p2), frequencies = solver.solve(init_position=current_position) # frequencies are already in GHz for frequency in frequencies: result_dict_dipole["frequency"].append(frequency) result_dict_dipole["Hmag"].append(Hmag) # we note the final position of the magnetisation in spherical coordinates result_dict_dipole["theta_1"].append(t1) result_dict_dipole["phi_1"].append(p1) result_dict_dipole["theta_2"].append(t2) result_dict_dipole["phi_2"].append(p2) # we reuse the previous solution as the initial guess for the next iteration current_position = [t1, p1, t2, p2]

100%|██████████| 30/30 [00:21<00:00,  1.42it/s]

In [6]:

import matplotlib.pyplot as plt

with plt.style.context(["science", "nature"]):
    w, h = plt.figaspect(0.5)
    fig, (ax1, ax2) = plt.subplots(1, 2, dpi=200, figsize=(w, h))
    ax1.plot(
        np.asarray(result_dict_dipole["Hmag"]) / 1e3,
        result_dict_dipole["frequency"],
        "bo",
        label="classical+dipole",
    )
    ax1.plot(
        np.asarray(result_dictionary["Hmag"]) / 1e3,
        result_dictionary["frequency"],
        "ro",
        markersize=0.9,
        label="classical",
    )
    ax1.set_xlabel("H (kA/m)")
    ax1.set_ylabel("f (GHz)")

    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dict_dipole["phi_1"]),
        color="forestgreen",
        label=r"$\phi_1$",
    )
    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dict_dipole["theta_1"]),
        color="navy",
        label=r"$\theta_1$",
    )

    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dict_dipole["phi_2"]),
        color="orange",
        label=r"$\phi_2$",
        marker="o",
        markersize=1,
    )
    ax2.plot(
        Hspace / 1e3,
        np.rad2deg(result_dict_dipole["theta_2"]),
        color="yellow",
        label=r"$\theta_2$",
        marker="o",
        markersize=1,
    )
    ax1.legend()
    ax2.legend()
    ax2.set_xlabel("H (kA/m)")
    ax2.set_ylabel("angle (deg)")
    fig.tight_layout()

import matplotlib.pyplot as plt with plt.style.context(["science", "nature"]): w, h = plt.figaspect(0.5) fig, (ax1, ax2) = plt.subplots(1, 2, dpi=200, figsize=(w, h)) ax1.plot( np.asarray(result_dict_dipole["Hmag"]) / 1e3, result_dict_dipole["frequency"], "bo", label="classical+dipole", ) ax1.plot( np.asarray(result_dictionary["Hmag"]) / 1e3, result_dictionary["frequency"], "ro", markersize=0.9, label="classical", ) ax1.set_xlabel("H (kA/m)") ax1.set_ylabel("f (GHz)") ax2.plot( Hspace / 1e3, np.rad2deg(result_dict_dipole["phi_1"]), color="forestgreen", label=r"$\phi_1$", ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dict_dipole["theta_1"]), color="navy", label=r"$\theta_1$", ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dict_dipole["phi_2"]), color="orange", label=r"$\phi_2$", marker="o", markersize=1, ) ax2.plot( Hspace / 1e3, np.rad2deg(result_dict_dipole["theta_2"]), color="yellow", label=r"$\theta_2$", marker="o", markersize=1, ) ax1.legend() ax2.legend() ax2.set_xlabel("H (kA/m)") ax2.set_ylabel("angle (deg)") fig.tight_layout()