Комментарии / Профиль Maximka200 / Хабр

Как это происходит в реальности, наверное, знает Кирилл Циберкин. Он по этой теме защитил докторскую диссертацию и последние 5 лет работал над этой задачей.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 2 сен в 11:30

А как можно исправить исходную статью? Полностью переделывать или писать новую длинную на ту же тему я не собираюсь.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 2 сен в 09:52

1. Почему llg_rhs разные? Форма Гилберта ($1/(1+\alpha^2)$) vs форма Ландау–Лифшица (без неё) — разница пренебрежимо мала при $\alpha=0.005$. Но знак damping отличался: 1D-оригинал имел + (антидемпфинг при $\gamma > 0$), 2D-оригинал имел - (правильно). При обезразмеривании я поменял знак, не сказав об этом.

2. Шаг сетки. Без μ0: $dy' \approx 157$ — огромный, тривиально устойчивый, но физика неверна. С μ0: $dy' \approx 0.18$, и критерий устойчивости RK2 для LLG — не $dt \leq dy^2/2$ (это для диффузии), а $u = 4dt/dy^2 \lesssim 0.3$ (прецессия + слабый демпфинг). 2D-код: $u \approx 0.11$ — на грани, но устойчив.

3. Стенка не расплывалась из-за комбинации всех трёх багов. «Презумпция виновности» полностью оправдалась.

Вот исправленный код:

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io

# =========================================================
# 1. РАЗМЕРНЫЕ ПАРАМЕТРЫ
# =========================================================
Ms = 800e3; A_ex = 13e-12; alpha = 0.005; gamma = 2.211e5
mu0 = 4 * np.pi * 1e-7          # ← ВСЁ В SI, μ0 ОБЯЗАТЕЛЕН

# =========================================================
# 2. МАСШТАБЫ (с μ0!)
#    l0 = sqrt(2A/(μ0·Ms²))  →  J = 1
#    t0 = 1/(γ·Ms)
# =========================================================
l0 = np.sqrt(2.0 * A_ex / (mu0 * Ms**2))   # ≈ 5.686 нм
t0 = 1.0 / (gamma * Ms)                     # ≈ 5.654 пс

# =========================================================
# 3. БЕЗРАЗМЕРНЫЕ ВЕЛИЧИНЫ
# =========================================================
J = 1.0
h_ext_x = 4e4 / Ms                          # ≈ 0.05

L = 200e-9 / l0                             # ≈ 35.17
N = 101
y = np.linspace(-L/2, L/2, N)
dy = y[1] - y[0]                            # ≈ 0.352

# Устойчивость RK2 для LLG: u = 4·dt/dy² ≲ 0.3
dt = 0.003                                  # u ≈ 0.097 ✓
tau_max = 200.0                              # ≈ 1.13 нс
n_steps = int(tau_max / dt)
save_every = max(1, n_steps // 50)

# --- Начальное условие: блоховская стенка ---
delta = 20e-9 / l0                          # ≈ 3.52
theta = np.pi/2 * (1.0 + np.tanh(y / delta))
m = np.stack([np.sin(theta), np.cos(theta), np.zeros_like(y)], axis=0)

# =========================================================
# 4. ОБМЕННОЕ ПОЛЕ (J=1, безразмерное)
# =========================================================
def lap1d(f, dy):
    return (np.roll(f, -1) - 2*f + np.roll(f, 1)) / dy**2

def compute_h(mv):
    return np.stack([
        lap1d(mv[0], dy) + h_ext_x,
        lap1d(mv[1], dy),
        lap1d(mv[2], dy)
    ], axis=0)

# =========================================================
# 5. УРАВНЕНИЕ ЛЛГ (Gilbert, ПРАВИЛЬНЫЙ ЗНАК damping = МИНУС)
#    dm/dτ = -1/(1+α²)·(m×h) - α/(1+α²)·(m×(m×h))
# =========================================================
def llg_rhs(mv):
    h = compute_h(mv)
    c1 = np.cross(mv, h, axis=0)
    c2 = np.cross(mv, c1, axis=0)
    return -(1.0 / (1 + alpha**2)) * c1 \
           - (alpha / (1 + alpha**2)) * c2

# =========================================================
# 6. ИНТЕГРИРОВАНИЕ (RK2, БЕЗ нормировки)
# =========================================================
taus, my_p, mx_p, mz_p = [], [], [], []
mc = m.copy()

for step in range(n_steps):
    k1 = llg_rhs(mc)
    k2 = llg_rhs(mc + 0.5 * dt * k1)
    mc += dt * k2

    if step % save_every == 0 or step == n_steps - 1:
        taus.append(step * dt)
        my_p.append(mc[1].copy())
        mx_p.append(mc[0].copy())
        mz_p.append(mc[2].copy())

taus = np.array(taus)
my_p = np.array(my_p)
mx_p = np.array(mx_p)
mz_p = np.array(mz_p)
times_real = taus * t0
y_nm = y * l0 * 1e9

norm = np.sqrt(mc[0]**2 + mc[1]**2 + mc[2]**2)
print(f"l0={l0*1e9:.3f} нм, t0={t0*1e12:.3f} пс, J={J}")
print(f"dy={dy:.4f}, dt={dt:.4f}, u={4*dt/dy**2:.4f}")
print(f"Шагов: {n_steps}, кадров: {len(taus)}")
print(f"|m|: min={norm.min():.6f}, max={norm.max():.6f} (без нормировки)")
print(f"Время до {times_real[-1]*1e9:.2f} нс")

# =========================================================
# 7. ГРАФИК: 3 панели (m_y, m_z, |dm/dy|)
# =========================================================
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
indices = [0, len(taus)//4, len(taus)//2, 3*len(taus)//4, -1]
colors = ['#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd']

for i, idx in enumerate(indices):
    t_ns = times_real[idx] * 1e9
    axes[0].plot(y_nm, my_p[idx], color=colors[i], lw=2, label=f't={t_ns:.2f} нс')
    axes[1].plot(y_nm, mz_p[idx], color=colors[i], lw=2, label=f't={t_ns:.2f} нс')
    dmx = np.gradient(mx_p[idx], y_nm)
    dmy = np.gradient(my_p[idx], y_nm)
    dmz = np.gradient(mz_p[idx], y_nm)
    axes[2].plot(y_nm, np.sqrt(dmx**2+dmy**2+dmz**2), color=colors[i], lw=2, label=f't={t_ns:.2f} нс')

axes[0].set_title(r'$m_y$'); axes[0].set_ylabel(r'$m_y$')
axes[1].set_title(r'$m_z$ (3D прецессия!)'); axes[1].set_ylabel(r'$m_z$')
axes[2].set_title(r'$|d\mathbf{m}/dy|$ (расплывание)'); axes[2].set_ylabel('нм$^{-1}$')
for ax in axes:
    ax.set_xlabel('y (нм)'); ax.legend(fontsize=9); ax.grid(True, ls=':', alpha=0.3)
plt.tight_layout()
plt.savefig('wall_analysis.png', dpi=150)
plt.close()

# =========================================================
# 8. GIF
# =========================================================
fig2, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], lw=2, color='#1f77b4')
ax.set_xlim(y_nm[0], y_nm[-1]); ax.set_ylim(-1.1, 1.1)
ax.axhline(0, color='k', ls='--', lw=1, alpha=0.4)
ax.set_xlabel('y (нм)'); ax.set_ylabel(r'$m_y$')
title = ax.set_title(''); ax.grid(True, ls=':', alpha=0.3)

frames = []
for i in range(len(taus)):
    line.set_data(y_nm, my_p[i])
    title.set_text(f't = {times_real[i]*1e9:.2f} нс')
    buf = io.BytesIO()
    fig2.savefig(buf, format='png', dpi=80)
    buf.seek(0)
    frames.append(Image.open(buf).convert('RGB'))
plt.close(fig2)

frames[0].save('wall_spreading.gif', save_all=True,
               append_images=frames[1:], duration=50, loop=0)
print("✅ wall_analysis.png, wall_spreading.gif")

В результате получим:

Теперь всё верно? Дан ответ на все вопросы? Проблема решена? Полученная модель соответствует реальности?

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 2 сен в 06:21

Вот я сделал процедуру обезразмеривания и результат теперь другой:

Теперь всё реалистично

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 2 сен в 05:12

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io

# =========================================================
# 1. РАЗМЕРНЫЕ ПАРАМЕТРЫ (только для расчёта масштабов)
# =========================================================
gamma = 2.211e5
alpha = 0.005
A_ex = 1.3e-11
Ms = 8e5
mu0 = 4 * np.pi * 1e-7

# =========================================================
# 2. МАСШТАБЫ
#    l0 = sqrt(2A / (mu0*Ms^2))  →  J = 1
#    t0 = 1 / (gamma * Ms)
# =========================================================
l0 = np.sqrt(2.0 * A_ex / (mu0 * Ms**2))   # ≈ 5.69 нм
t0 = 1.0 / (gamma * Ms)                     # ≈ 5.65 пс

# =========================================================
# 3. БЕЗРАЗМЕРНЫЕ ВЕЛИЧИНЫ
# =========================================================
Nx, Ny = 64, 64
dx = 10e-9 / l0            # безразмерный шаг сетки ≈ 1.759
dt = 5e-13 / t0             # безразмерный шаг ≈ 0.0884
steps = 1200
save_every = 15

h_ext = 4.0e4 / Ms          # безразмерное внешнее поле = 0.05

# Начальное условие: m = (1, 0, 0) везде
m = np.ones((Nx, Ny, 3)) * np.array([1.0, 0.0, 0.0])

# Возмущение в центре
cx, cy = Nx // 2, Ny // 2
angle = 0.52
radius_sq = 100

for ix in range(Nx):
    for iy in range(Ny):
        dist_sq = (ix - cx) ** 2 + (iy - cy) ** 2
        if dist_sq < radius_sq:
            m[ix, iy, 0] = np.cos(angle)
            m[ix, iy, 1] = np.sin(angle)

# Нормировка один раз в начале
norm = np.sqrt(np.sum(m ** 2, axis=2, keepdims=True))
norm[norm == 0] = 1.0
m = m / norm

# =========================================================
# 4. БЕЗРАЗМЕРНОЕ ОБМЕННОЕ ПОЛЕ (J = 1)
#    h_ex = nabla'^2 m  (через roll, периодические границы)
# =========================================================
def compute_h(m_vec):
    h = np.zeros_like(m_vec)
    for i in range(3):
        d2x = (np.roll(m_vec[:, :, i], -1, axis=0)
               - 2 * m_vec[:, :, i]
               + np.roll(m_vec[:, :, i], 1, axis=0)) / dx**2
        d2y = (np.roll(m_vec[:, :, i], -1, axis=1)
               - 2 * m_vec[:, :, i]
               + np.roll(m_vec[:, :, i], 1, axis=1)) / dx**2
        h[:, :, i] = d2x + d2y
    return h

# =========================================================
# 5. БЕЗРАЗМЕРНОЕ УРАВНЕНИЕ ЛЛГ
#    dm/dτ = -1/(1+α²) (m × h) - α/(1+α²) m × (m × h)
# =========================================================
def llg_rhs(m_vec):
    h = compute_h(m_vec)
    h[:, :, 0] += h_ext
    cross_mh = np.cross(m_vec, h)
    cross_m_mh = np.cross(m_vec, cross_mh)
    return -(1.0 / (1 + alpha**2)) * cross_mh \
           - (alpha / (1 + alpha**2)) * cross_m_mh

# =========================================================
# 6. ИНТЕГРИРОВАНИЕ (RK2 — метод средней точки, без нормировки)
# =========================================================
frames = []
print("Считаем кадры для анимации...")

for step in range(steps):
    k1 = llg_rhs(m)
    k2 = llg_rhs(m + 0.5 * dt * k1)
    m += dt * k2

    if np.any(~np.isfinite(m)):
        print(f"⚠️ NaN на шаге {step} — уменьшите dt")
        break

    if step % save_every == 0 or step == steps - 1:
        frames.append(m.copy())

print(f"Всего кадров: {len(frames)}")

# Проверка нормы
norm_check = np.sqrt(np.sum(m**2, axis=2))
print(f"|m|: min={norm_check.min():.6f}, max={norm_check.max():.6f} (без нормировки)")

# =========================================================
# 7. АНИМАЦИЯ (GIF через Pillow, без FuncAnimation.save)
# =========================================================
fig, ax = plt.subplots(figsize=(8, 8))

im = ax.imshow(m[:, :, 0], cmap='coolwarm', origin='lower', vmin=-1, vmax=1)
cbar = fig.colorbar(im, ax=ax)
cbar.set_label(r'$m_x$')

stride = 4
x_vec, y_vec = np.meshgrid(
    np.arange(0, Nx, stride),
    np.arange(0, Ny, stride),
    indexing='ij'
)
quiver = ax.quiver(x_vec, y_vec,
                   np.zeros_like(x_vec), np.zeros_like(y_vec),
                   color='black', width=0.003, scale=30)

ax.set_title("2D спиновая волна: распространение возмущения")
ax.set_xlabel("x (ячейка)")
ax.set_ylabel("y (ячейка)")
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12,
                    bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))

frames_pil = []
for idx in range(len(frames)):
    m_frame = frames[idx]
    t_ns = idx * save_every * dt * t0 * 1e9

    im.set_data(m_frame[:, :, 0])

    u = m_frame[::stride, ::stride, 0]
    v = m_frame[::stride, ::stride, 1]
    quiver.set_UVC(u, v)

    time_text.set_text(f"t = {t_ns:.2f} нс")

    buf = io.BytesIO()
    fig.savefig(buf, format='png', dpi=100)
    buf.seek(0)
    frames_pil.append(Image.open(buf).convert('RGB'))

plt.close(fig)

gif_file = 'spin_wave_2d_nondim.gif'
if frames_pil:
    frames_pil[0].save(
        gif_file,
        save_all=True,
        append_images=frames_pil[1:],
        duration=50,
        loop=0
    )
    print(f"✅ GIF сохранён: {gif_file}")

Сделал обезразмеривание в исходном 2D коде.Результат теперь принципиально другой:

Стало реалистичнее?

Проблема с пространственно разрывной намагниченностью решена?

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 2 сен в 04:46

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io

# =========================================================
# 1. РАЗМЕРНЫЕ ПАРАМЕТРЫ (только для расчёта масштабов)
# =========================================================
Ms = 800e3
A_ex = 13e-12
alpha = 0.005
gamma = 2.211e5

# =========================================================
# 2. МАСШТАБЫ (без mu0 —consistent с исходным кодом)
#    l0 = sqrt(2A / Ms^2)  →  J = 1
#    t0 = 1 / (gamma * Ms)
# =========================================================
l0 = np.sqrt(2.0 * A_ex / Ms**2)     # ≈ 6.374 нм
t0 = 1.0 / (gamma * Ms)              # ≈ 5.654 пс

# =========================================================
# 3. БЕЗРАЗМЕРНЫЕ ВЕЛИЧИНЫ
# =========================================================
J = 1.0
h_ext_x = 4e4 / Ms                   # ≈ 0.05

L = 200e-9 / l0                      # ≈ 31.39
N = 401
y = np.linspace(-L / 2, L / 2, N)
dy = y[1] - y[0]

# Начальное условие: доменная стенка
delta = 20e-9 / l0                   # ≈ 3.14
theta = np.pi / 2 * (1.0 + np.tanh(y / delta))
mx = np.sin(theta)
my = np.cos(theta)
mz = np.zeros_like(y)
m = np.stack([mx, my, mz], axis=0)

# =========================================================
# 4. БЕЗРАЗМЕРНОЕ ЭФФЕКТИВНОЕ ПОЛЕ (J = 1)
# =========================================================
def compute_h(m_vec):
    mx, my, mz = m_vec
    d2mx = np.gradient(np.gradient(mx, dy), dy)
    d2my = np.gradient(np.gradient(my, dy), dy)
    d2mz = np.gradient(np.gradient(mz, dy), dy)
    h = np.stack([d2mx, d2my, d2mz], axis=0)
    h[0] += h_ext_x
    return h

# =========================================================
# 5. БЕЗРАЗМЕРНОЕ УРАВНЕНИЕ ЛЛГ
# =========================================================
def llg_rhs(m_vec):
    h = compute_h(m_vec)
    mx, my, mz = m_vec
    hx, hy, hz = h

    cx = my * hz - mz * hy
    cy = mz * hx - mx * hz
    cz = mx * hy - my * hx

    c2x = my * cz - mz * cy
    c2y = mz * cx - mx * cz
    c2z = mx * cy - my * cx

    return np.stack([
        -cx + alpha * c2x,
        -cy + alpha * c2y,
        -cz + alpha * c2z
    ], axis=0)

# =========================================================
# 6. ИНТЕГРИРОВАНИЕ: метод средней точки (RK2)
#    k1 = f(m^n)
#    k2 = f(m^n + dt/2 * k1)
#    m^{n+1} = m^n + dt * k2
#
#    Дрейф |m| за шаг — O(dt^4), за всё время — доли процента.
#    Нормировка не нужна.
# =========================================================
dt = 0.01
tau_max = 2e-9 / t0                   # ≈ 353.8
n_steps = int(tau_max / dt)
save_every = max(1, n_steps // 60)

taus, my_profiles, mx_profiles = [], [], []
m_curr = m.copy()

for step in range(n_steps):
    k1 = llg_rhs(m_curr)
    k2 = llg_rhs(m_curr + 0.5 * dt * k1)
    m_curr += dt * k2

    if np.any(~np.isfinite(m_curr)):
        print(f"⚠️ NaN на шаге {step} — уменьшите dt")
        break

    if step % save_every == 0 or step == n_steps - 1:
        taus.append(step * dt)
        my_profiles.append(m_curr[1].copy())
        mx_profiles.append(m_curr[0].copy())

taus = np.array(taus)
my_profiles = np.array(my_profiles)
mx_profiles = np.array(mx_profiles)
times_real = taus * t0

# Проверка нормы
norm_check = np.sqrt(m_curr[0]**2 + m_curr[1]**2 + m_curr[2]**2)
print(f"l0 = {l0*1e9:.3f} нм,  t0 = {t0*1e12:.3f} пс,  J = {J}")
print(f"Шагов: {n_steps},  кадров: {len(taus)},  время до {times_real[-1]*1e9:.2f} нс")
print(f"|m|: min={norm_check.min():.6f}, max={norm_check.max():.6f} (без нормировки)")
print(f"m_y в центре = {my_profiles[-1][N//2]:.3f}")

# =========================================================
# 7. АНИМАЦИЯ (GIF через Pillow)
# =========================================================
fig, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], linewidth=2, color='#1f77b4')
arrow_holder = [None]
y_nm = y * l0 * 1e9

ax.set_xlim(y_nm[0], y_nm[-1])
ax.set_ylim(-1.1, 1.1)
ax.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.4)
ax.set_xlabel('y (нм)', fontsize=12)
ax.set_ylabel(r'$m_y$', fontsize=12)
title = ax.set_title('', fontsize=14)
ax.grid(True, linestyle=':', alpha=0.3)

frames = []
for i in range(len(taus)):
    line.set_data(y_nm, my_profiles[i])

    if arrow_holder[0] is not None:
        arrow_holder[0].remove()
    cmx = mx_profiles[i][N // 2]
    cmy = my_profiles[i][N // 2]
    arrow_holder[0] = ax.arrow(
        0, -0.3, cmx * 20, cmy * 20,
        width=0.012, color='red',
        head_width=0.05, head_length=0.09
    )

    title.set_text(f't = {times_real[i] * 1e9:.2f} нс')

    buf = io.BytesIO()
    fig.savefig(buf, format='png', dpi=100)
    buf.seek(0)
    frames.append(Image.open(buf).convert('RGB'))

plt.close(fig)

gif_file = 'domain_wall_nondim_RK2.gif'
if frames:
    frames[0].save(
        gif_file,
        save_all=True,
        append_images=frames[1:],
        duration=50,
        loop=0
    )
    print(f"✅ GIF сохранён: {gif_file}")

Теперь всё полностью корректно:

результат-

стало быстрее, но всё непрерывно

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 1 сен в 08:08

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation

# =========================================================
# ПАРАМЕТРЫ ПЕРМАЛЛОЯ (размерные)
# =========================================================
Ms = 800e3
A_ex = 13e-12
alpha = 0.005
gamma = 2.211e5

# =========================================================
# ГЕОМЕТРИЯ И СЕТКА
# =========================================================
L = 200e-9
N = 401
y = np.linspace(-L/2, L/2, N)
dy = y[1] - y[0]

# Безразмерные единицы
l0 = 10e-9                     # единица длины (10 нм)
t0 = 1.0 / (gamma * Ms)        # единица времени (~0.57 пс)

y_dimless = y / l0
dy_dimless = dy / l0

# Безразмерный обменный коэффициент
J_dimless = (2.0 * A_ex) / (Ms**2 * l0**2) * t0 * gamma

# Поле размагничивания и внешнее поле (безразмерные)
N_perp = 0.25
H_ext_x_dimless = 4e4 / Ms

# =========================================================
# НАЧАЛЬНОЕ УСЛОВИЕ
# =========================================================
center_idx = N // 2
width_dimless = 4.0

mx_init = np.tanh((y_dimless[center_idx] - y_dimless) / (width_dimless / 2.0))
my_init = np.zeros_like(y_dimless)
mz_init = np.zeros_like(y_dimless)

norm = np.sqrt(mx_init**2 + my_init**2 + mz_init**2)
norm[norm == 0] = 1.0
mx_init /= norm
my_init /= norm
mz_init /= norm

m = np.stack([mx_init, my_init, mz_init], axis=0)  # [3, N], |m|=1

# =========================================================
# ЭФФЕКТИВНОЕ ПОЛЕ (БЕЗРАЗМЕРНОЕ)
# =========================================================
def compute_h_eff(m_vec):
    mx, my, mz = m_vec
    d2mx = np.gradient(np.gradient(mx, dy_dimless), dy_dimless)
    d2my = np.gradient(np.gradient(my, dy_dimless), dy_dimless)
    d2mz = np.gradient(np.gradient(mz, dy_dimless), dy_dimless)

    h_exchange_x = J_dimless * d2mx
    h_exchange_y = J_dimless * d2my
    h_exchange_z = J_dimless * d2mz

    h_demag_y = -N_perp * my

    h_ext = np.array([H_ext_x_dimless, 0.0, 0.0])

    h_eff = np.stack(
        [h_exchange_x + h_ext[0], h_exchange_y + h_demag_y, h_exchange_z + h_ext[2]],
        axis=0
    )
    return h_eff

# =========================================================
# ПРАВАЯ ЧАСТЬ УРАВНЕНИЯ ЛЛГ (БЕЗРАЗМЕРНАЯ)
# =========================================================
def llg_rhs(m_vec, tau):
    mx, my, mz = m_vec
    h = compute_h_eff(m_vec)
    hx, hy, hz = h

    cx = my * hz - mz * hy
    cy = mz * hx - mx * hz
    cz = mx * hy - my * hx

    c2x = my * cz - mz * cy
    c2y = mz * cx - mx * cz
    c2z = mx * cy - my * cx

    dm_dt_x = -cx - alpha * c2x
    dm_dt_y = -cy - alpha * c2y
    dm_dt_z = -cz - alpha * c2z

    return np.stack([dm_dt_x, dm_dt_y, dm_dt_z], axis=0)

# =========================================================
# ИНТЕГРИРОВАНИЕ (РЕЛАКСАЦИЯ + ДВИЖЕНИЕ)
# =========================================================
dt_tau = 0.05

# Релаксация
tau_relax = 50.0
n_steps_relax = int(tau_relax / dt_tau)
m_curr = m.copy()

for step in range(n_steps_relax):
    dm = llg_rhs(m_curr, 0.0)
    m_curr += dm * dt_tau
    # Нормировки нет: структура уравнения сохраняет длину

# Движение
tau_move = 60.0
n_steps_move = int(tau_move / dt_tau)
save_every = 3

taus = []
my_profiles = []
mx_profiles = []  # для стрелки

for step in range(n_steps_move):
    dm = llg_rhs(m_curr, step * dt_tau)
    m_curr += dm * dt_tau

    if step % save_every == 0:
        taus.append(step * dt_tau)
        my_profiles.append(m_curr[1].copy())
        mx_profiles.append(m_curr[0].copy())

taus = np.array(taus)
my_profiles = np.array(my_profiles)
mx_profiles = np.array(mx_profiles)

times_real = taus * t0  # секунды

print(f"✅ Рассчитано кадров анимации: {len(taus)}")
print(f"[РЕЛАКСАЦИЯ] m_y в центре = {m_curr[1][N//2]:.3f} (не ноль — реалистично)")

# =========================================================
# ПОДГОТОВКА АНИМАЦИИ
# =========================================================
fig, ax = plt.subplots(figsize=(10, 6))

line, = ax.plot([], [], linewidth=2, color='#1f77b4')
arrow = ax.arrow(0, 0, 0, 0, width=0.015, color='red', head_width=0.06, head_length=0.12)

ax.set_xlim(y[0]*1e9, y[-1]*1e9)
ax.set_ylim(-1.1, 1.1)
ax.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.4)

ax.set_xlabel('Координата y (нм)', fontsize=12)
ax.set_ylabel(r'$m_y$', fontsize=12)
title_text = ax.set_title('', fontsize=14)
ax.grid(True, linestyle=':', alpha=0.3)

def init():
    line.set_data([], [])
    arrow.remove()
    global arrow
    arrow = ax.arrow(0, 0, 0, 0, width=0.015, color='red', head_width=0.06, head_length=0.12)
    title_text.set_text('')
    return line, arrow, title_text

def update(frame):
    y_nm = y * 1e9
    my = my_profiles[frame]
    mx = mx_profiles[frame]

    line.set_data(y_nm, my)

    # Удаляем старую стрелку и рисуем новую
    arrow.remove()
    center_y_nm = 0
    center_my = my[N//2]
    center_mx = mx[N//2]

    # Стрелка: рисуем от центра немного вниз, чтобы не перекрывать график
    arrow_y_start = -0.3
    arrow_dx = center_mx * 20  # масштабируем стрелку
    arrow_dy = center_my * 20

    global arrow
    arrow = ax.arrow(center_y_nm, arrow_y_start, arrow_dx, arrow_dy,
                     width=0.012, color='red', head_width=0.05, head_length=0.09)

    t_ns = times_real[frame] * 1e9
    title_text.set_text(f'Динамика доменной стенки: t = {t_ns:.2f} нс')
    return line, arrow, title_text

ani = FuncAnimation(fig, update, frames=len(taus), init_func=init, blit=True, interval=40)

# Сохранение в GIF
gif_file = 'domain_wall_animation.gif'
ani.save(gif_file, writer='pillow', fps=20, dpi=200)
print(f"✅ Анимация GIF сохранена как {gif_file}")

# Сохранение в MP4 (нужен ffmpeg в PATH)
mp4_file = 'domain_wall_animation.mp4'
try:
    ani.save(mp4_file, writer='ffmpeg', fps=20, dpi=200)
    print(f"✅ Анимация MP4 сохранена как {mp4_file}")
except RuntimeError:
    print("⚠️ Не удалось сохранить MP4: проверьте, установлен ли ffmpeg и добавлен ли в PATH.")

plt.close(fig)

Вот сделал обезразмеривание и убрал нормировку как вы просили. Результат такой же:

Все непрерывно.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 1 сен в 05:58

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# --- Параметры пермаллоя (Ni80Fe20) ---
Ms = 800e3          # А/м, намагниченность насыщения
A_ex = 13e-12       # Дж/м, обменная константа
alpha = 0.005       # безразмерный параметр затухания
gamma = 2.211e5     # рад/(с·А/м), гиромагнитное отношение

# --- Сетка (одномерная вдоль y) ---
L = 200e-9          # длина образца (200 нм)
N = 401             # нечётное число точек (центр ровно посередине)
y = np.linspace(-L/2, L/2, N)
dy = y[1] - y[0]

# --- Внешнее поле (вдоль x) ---
H_ext_x = 4e4      # А/м (поле, которое сдвигает доменную стенку)

# --- Начальное условие: доменная стенка (Блоховская) ---
# Толщина стенки без анизотропии: delta ~ sqrt(A/K), K=0 → берём условную толщину
delta = 20e-9        # 20 нм — реалистичная толщина стенки для пермаллоя
theta_init = np.pi/2 * (1.0 + np.tanh(y / delta))

mx = np.sin(theta_init)
my = np.cos(theta_init)
mz = np.zeros_like(y)

M = np.stack([mx, my, mz], axis=0)  # M[0]=Mx, M[1]=My, M[2]=Mz

# --- Эффективное поле: обмен + внешнее ---
def compute_Heff(M_norm):
    mx, my, mz = M_norm

    # Вторая производная через np.gradient (дважды)
    d2mx_dy2 = np.gradient(np.gradient(mx, dy), dy)
    d2my_dy2 = np.gradient(np.gradient(my, dy), dy)
    d2mz_dy2 = np.gradient(np.gradient(mz, dy), dy)

    pref = (2.0 * A_ex) / (Ms**2)
    H_exchange_x = pref * d2mx_dy2
    H_exchange_y = pref * d2my_dy2
    H_exchange_z = pref * d2mz_dy2

    H_ext = np.array([H_ext_x, 0.0, 0.0])

    Heff = np.stack([H_exchange_x, H_exchange_y, H_exchange_z], axis=0) + H_ext[:, None]
    return Heff

# --- Правая часть уравнения ЛЛГ (Gilbert form) ---
def llg_rhs(M_norm, t):
    mx, my, mz = M_norm
    Heff = compute_Heff(M_norm)
    Hx, Hy, Hz = Heff

    # M × H_eff
    cx = my * Hz - mz * Hy
    cy = mz * Hx - mx * Hz
    cz = mx * Hy - my * Hx

    # M × (M × H_eff)
    c2x = my * cz - mz * cy
    c2y = mz * cx - mx * cz
    c2z = mx * cy - my * cx

    dMdt_x = -gamma * cx + alpha * gamma * c2x
    dMdt_y = -gamma * cy + alpha * gamma * c2y
    dMdt_z = -gamma * cz + alpha * gamma * c2z

    return np.stack([dMdt_x, dMdt_y, dMdt_z], axis=0)

# --- Интегрирование (явный Эйлер с нормировкой) ---
dt = 5e-14           # 50 фс — маленький шаг для устойчивости
t_max = 2e-9         # 2 нс
save_every = 4       # сохранять каждый 4-й кадр

n_steps = int(t_max / dt)

times = []
My_profiles = []

M_curr = M.copy()
t = 0.0

for step in range(n_steps):
    dM = llg_rhs(M_curr, t)
    M_curr += dM * dt

    # Нормировка |M| = 1 (строго, без NaN)
    norm = np.sqrt(M_curr[0]**2 + M_curr[1]**2 + M_curr[2]**2)
    # Защита от нулевого вектора (на всякий случай)
    norm[norm == 0] = 1.0
    M_curr /= norm

    t += dt

    if step % save_every == 0:
        times.append(t)
        My_profiles.append(M_curr[1].copy())

times = np.array(times)
My_profiles = np.array(My_profiles)

# --- Анимация ---
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlim(-L/2 * 1e9, L/2 * 1e9)
ax.set_ylim(-1.1, 1.1)
ax.set_xlabel('y (нм)')
ax.set_ylabel(r'$M_y / M_s$')
ax.set_title('Динамика доменной стенки в пермаллое (x=0)')
ax.grid(True, alpha=0.3)

line, = ax.plot([], [], 'b-', linewidth=2)
time_text = ax.text(0.05, 0.95, '', transform=ax.transAxes, fontsize=12, verticalalignment='top')

def init():
    line.set_data([], [])
    time_text.set_text('')
    return line, time_text

def animate(i):
    y_nm = y * 1e9
    line.set_data(y_nm, My_profiles[i])
    time_text.set_text(f't = {times[i]*1e9:.2f} нс')
    return line, time_text

ani = FuncAnimation(fig, animate, frames=len(times), init_func=init, blit=True, interval=30)
plt.show()

# Чтобы сохранить в GIF (если установлен pillow):
ani.save('permalloy_domain_wall.gif', writer='pillow', fps=15)

Вот код как был получен этот график.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 1 сен в 03:31

В любом случае спасибо вам большое за интересную беседу.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 1 сен в 03:05

Всё непрерывно, никакой пространственно разрывной намагниченности нет.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 31 авг в 11:55

Если можно взять вторую производную(как вы утверждаете), то намагниченность непрерывна по y.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 30 авг в 20:00

На какой конкретно картинке видно пространственно разрывную намагниченность? Я не вижу.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 30 авг в 18:12

Что экстраординарного в этом результате?

В относительно низких полях колебания все равно есть, а значит и волны тоже. Возможно, они распространяются очень медленно и мы их на гифке не успеваем увидеть.

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 30 авг в 15:51

Задавайте вопросы, в чем я не разобрался?

Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python

Maximka200 30 авг в 04:38

Устойчивость численного метода (динамических систем) можно оценить при помощи показателя Ляпунова. С устойчивостью счета я разобрался.

2 3 4 5

Информация

Специализация