#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
1次元単原子チェーン（最近接相互作用）における非調和フォノン振動数の温度シフトを計算するモジュール。
関連リンク: :doc:`anharmonic3D_usage`
詳細説明:
モデル（相対変位 Δ_n = u_{n+1} - u_n）:
V = Σ_n [ (K2/2) Δ_n^2 + (K3/3!) Δ_n^3 + (K4/4!) Δ_n^4 ]
以下の摂動論的計算を行います:
ω_q(T) ≈ ω_q + Δω_q^(4)(T) + Δω_q^(3)(T)
ここで
Δω_q ≈ Re Σ(q, ω_q) / (2 ω_q)
- 4次 (K4) "Tadpole" シフト: ⟨|u_k|^2⟩ の単純な熱平均。
- 3次 (K3) "Bubble" シフト: kに対する主値積分（小さな eta で正則化）。
単位系:
- デフォルトでSI単位系を使用 (a: m, M: kg, K2: N/m, K3: N/m^2, K4: N/m^3)。
- 内部の振動数は rad/s で計算され、出力時に THz へ変換可能です。
k2_func/k3_func/k4_func を編集することで、K2/K3/K4 を (delta, T) の関数としてカスタマイズできます。
"""
from __future__ import annotations
import argparse
import math
from dataclasses import dataclass
from typing import Callable, Tuple, List
import numpy as np
# Physical constants (SI)
HBAR = 1.054_571_817e-34 # J*s
KB = 1.380_649e-23 # J/K
TWOPI = 2.0 * math.pi
# -----------------------------
# User-customizable force constants (as functions)
# -----------------------------
[ドキュメント]
def k2_func(delta: float = 0.0, T: float | None = None) -> float:
"""
調和力の定数K2を計算する。
詳細説明:
調和力の定数K2 [N/m] (無次元規則を使用する場合は一貫した単位)。
任意の関数に置き換え可能。
:param delta: float, 歪み (デフォルト: 0.0)
:param T: float | None, 温度 (デフォルト: None)
:returns: float, 調和力の定数K2
"""
return 10.0 # example
[ドキュメント]
def k3_func(delta: float = 0.0, T: float | None = None) -> float:
"""
3次の非調和定数K3を計算する。
詳細説明:
3次の非調和定数K3 [N/m^2]。
任意の関数に置き換え可能。ゼロ以外の値を設定するとBubbleシフトが有効になります。
:param delta: float, 歪み (デフォルト: 0.0)
:param T: float | None, 温度 (デフォルト: None)
:returns: float, 3次の非調和定数K3
"""
return 0.0 # example (set nonzero to enable bubble shift)
[ドキュメント]
def k4_func(delta: float = 0.0, T: float | None = None) -> float:
"""
4次の非調和定数K4を計算する。
詳細説明:
4次の非調和定数K4 [N/m^3]。
任意の関数に置き換え可能。ゼロ以外の値を設定するとTadpoleシフトが有効になります。
:param delta: float, 歪み (デフォルト: 0.0)
:param T: float | None, 温度 (デフォルト: None)
:returns: float, 4次の非調和定数K4
"""
return 1.0e-5 # example (set nonzero to enable tadpole shift)
# -----------------------------
# Core model helpers
# -----------------------------
[ドキュメント]
def C_of_q(q: np.ndarray, a: float) -> np.ndarray:
"""
関数 C(q) = e^{iqa} - 1 を計算する。
詳細説明:
主に |C(q)|^2 = 4 sin^2(qa/2) として利用される。
:param q: np.ndarray, 波数ベクトル
:param a: float, 格子定数
:returns: np.ndarray, 計算されたC(q)の複素数値
"""
return np.exp(1j * q * a) - 1.0
[ドキュメント]
def absC2(q: np.ndarray, a: float) -> np.ndarray:
"""
|C(q)|^2 = 4 sin^2(qa/2) を計算する。
:param q: np.ndarray, 波数ベクトル
:param a: float, 格子定数
:returns: np.ndarray, |C(q)|^2の実数値
"""
return 4.0 * np.sin(0.5 * q * a) ** 2
[ドキュメント]
def omega_harmonic(q: np.ndarray, a: float, M: float, K2: float) -> np.ndarray:
"""
調和分散（最近接相互作用）による振動数を計算する。
詳細説明:
ω_q^2 = (4 K2 / M) sin^2(qa/2) を用いて振動数を求める。
:param q: np.ndarray, 波数ベクトル
:param a: float, 格子定数
:param M: float, 質量
:param K2: float, 調和力の定数
:returns: np.ndarray, 計算された調和振動数
"""
return np.sqrt((4.0 * K2 / M) * np.sin(0.5 * q * a) ** 2)
[ドキュメント]
def bose_n(omega: np.ndarray, T: float, use_quantum: bool = True) -> np.ndarray:
"""
ボース分布関数 n(ω) = 1/(exp(β ħ ω) - 1) を計算する。
詳細説明:
use_quantum=False の場合は古典極限 n ~ kT/(ħ ω) を用いる。
:param omega: np.ndarray, 振動数
:param T: float, 温度
:param use_quantum: bool, 量子論的計算を使用するかどうか (デフォルト: True)
:returns: np.ndarray, ボース分布関数の値
"""
if T <= 0.0:
return np.zeros_like(omega)
if not use_quantum:
# classical limit
return (KB * T) / (HBAR * np.maximum(omega, 1e-300))
x = (HBAR * omega) / (KB * T)
# avoid overflow for large x
x = np.clip(x, 0.0, 700.0)
return 1.0 / (np.expm1(x) + 1e-300)
[ドキュメント]
def coth_halfbeta_hw(omega: np.ndarray, T: float, use_quantum: bool = True) -> np.ndarray:
"""
双曲線余接関数 coth(βħω/2) を計算する。
詳細説明:
T -> 0 で値は 1 に漸近する。古典極限では coth ~ 2kT/(ħω) を用いる。
:param omega: np.ndarray, 振動数
:param T: float, 温度
:param use_quantum: bool, 量子論的計算を使用するかどうか (デフォルト: True)
:returns: np.ndarray, coth(βħω/2) の値
"""
if T <= 0.0:
return np.ones_like(omega)
if not use_quantum:
return (2.0 * KB * T) / (HBAR * np.maximum(omega, 1e-300))
x = (HBAR * omega) / (2.0 * KB * T)
x = np.clip(x, 0.0, 350.0)
# coth(x) = cosh(x)/sinh(x) = (e^x + e^-x)/(e^x - e^-x)
ex = np.exp(x)
emx = 1.0 / ex
return (ex + emx) / (ex - emx + 1e-300)
# -----------------------------
# Quartic (K4) tadpole-like shift
# -----------------------------
[ドキュメント]
def delta_omega_quartic_tadpole(
q: float,
kgrid: np.ndarray,
a: float,
M: float,
K2: float,
K4: float,
T: float,
use_quantum: bool,
k_cut: float,
) -> float:
"""
4次（K4）のTadpoleシフト Δω_q^(4) を計算する。
詳細説明:
以下の構造を持つ1次元の実装:
Σ^(4)(q) ∝ K4 |C(q)|^2 * (1/N) Σ_k |C(k)|^2 <|u_k|^2>
<|u_k|^2> = ħ/(2 M ω_k) coth(βħω_k/2)
Δω_q^(4) = ReΣ/(2ω_q) を返す。
:param q: float, 波数q
:param kgrid: np.ndarray, 波数kのグリッド
:param a: float, 格子定数
:param M: float, 質量
:param K2: float, 調和力の定数
:param K4: float, 4次非調和定数
:param T: float, 温度
:param use_quantum: bool, 量子論的計算を使用するかどうか
:param k_cut: float, 赤外発散を防ぐためのkのカットオフ
:returns: float, 4次周波数シフト量 Δω_q^(4)
"""
omega_k = omega_harmonic(kgrid, a, M, K2)
omega_q = float(omega_harmonic(np.array([q]), a, M, K2)[0])
if omega_q <= 0.0:
return 0.0
mask = np.abs(kgrid) >= k_cut
kk = kgrid[mask]
omk = omega_k[mask]
# <|u_k|^2>
u2 = (HBAR / (2.0 * M * np.maximum(omk, 1e-300))) * coth_halfbeta_hw(omk, T, use_quantum=use_quantum)
S = np.mean(absC2(kk, a) * u2) # (1/N) Σ_k ...
# prefactor (model choice)
# Σ^(4) ~ (K4 / 2) |C(q)|^2 * S
Sigma4 = 0.5 * K4 * absC2(np.array([q]), a)[0] * S
d_omega = Sigma4 / (2.0 * omega_q)
return float(d_omega)
# -----------------------------
# Cubic (K3) bubble-like shift
# -----------------------------
[ドキュメント]
def v3_squared(
q: float,
k: np.ndarray,
a: float,
M: float,
K3: float,
omega_q: float,
omega_k: np.ndarray,
omega_qmk: np.ndarray,
) -> np.ndarray:
"""
3次相互作用の行列要素の2乗 |V3(q,k,q-k)|^2 を計算する。
詳細説明:
V3 ∝ K3 * C(q)C(k)C(q-k) / sqrt( (2M)^3 ω_q ω_k ω_{q-k} )
複素位相を避けるため |C|^2 = 4 sin^2(...) を使用する。
:param q: float, 波数q
:param k: np.ndarray, 波数kの配列
:param a: float, 格子定数
:param M: float, 質量
:param K3: float, 3次非調和定数
:param omega_q: float, 波数qにおける調和振動数
:param omega_k: np.ndarray, 波数kにおける調和振動数
:param omega_qmk: np.ndarray, 波数q-kにおける調和振動数
:returns: np.ndarray, |V3(q,k,q-k)|^2 の値
"""
# use |C| instead of C itself (phases cancel in |V|^2)
num = absC2(np.array([q]), a)[0] * absC2(k, a) * absC2((q - k), a)
den = (2.0 * M) ** 3 * np.maximum(omega_q * omega_k * omega_qmk, 1e-300)
return (K3 ** 2) * num / den
[ドキュメント]
def re_sigma_cubic_bubble(
q: float,
kgrid: np.ndarray,
a: float,
M: float,
K2: float,
K3: float,
T: float,
use_quantum: bool,
eta: float,
k_cut: float,
) -> float:
"""
Bubble構造による自己エネルギーの実部 Re Σ^(3)(q, ω_q) を計算する。
詳細説明:
etaを用いた正則化により主値積分を近似する:
Re[1/(x + iη)] = x/(x^2 + η^2)
:param q: float, 波数q
:param kgrid: np.ndarray, 波数kのグリッド
:param a: float, 格子定数
:param M: float, 質量
:param K2: float, 調和力の定数
:param K3: float, 3次非調和定数
:param T: float, 温度
:param use_quantum: bool, 量子論的計算を使用するかどうか
:param eta: float, 正則化パラメータ
:param k_cut: float, 赤外発散を防ぐためのkのカットオフ
:returns: float, 自己エネルギーの実部 Re Σ^(3)(q, ω_q)
"""
omega_q = float(omega_harmonic(np.array([q]), a, M, K2)[0])
if omega_q <= 0.0:
return 0.0
mask = np.abs(kgrid) >= k_cut
k = kgrid[mask]
omega_k = omega_harmonic(k, a, M, K2)
omega_qmk = omega_harmonic(q - k, a, M, K2)
nk = bose_n(omega_k, T, use_quantum=use_quantum)
nqm = bose_n(omega_qmk, T, use_quantum=use_quantum)
V2 = v3_squared(q, k, a, M, K3, omega_q, omega_k, omega_qmk)
def re_inv(x: np.ndarray) -> np.ndarray:
"""
実部の逆数近似を計算する内部関数。
詳細説明:
x / (x^2 + eta^2) を計算し、主値積分を正則化する。
:param x: np.ndarray, 分母の変数
:returns: np.ndarray, 正則化された値
"""
return x / (x * x + eta * eta)
w = omega_q
x1 = w - omega_k - omega_qmk
x2 = w - omega_k + omega_qmk
x3 = w + omega_k - omega_qmk
x4 = w + omega_k + omega_qmk
term1 = (nk + nqm + 1.0) * re_inv(x1)
term2 = (nk - nqm) * re_inv(x2)
term3 = (nqm - nk) * re_inv(x3)
term4 = (nk + nqm + 1.0) * re_inv(x4)
Sigma3 = np.sum(V2 * (term1 + term2 + term3 + term4)) / k.size
return float(Sigma3)
[ドキュメント]
def delta_omega_cubic_bubble(
q: float,
kgrid: np.ndarray,
a: float,
M: float,
K2: float,
K3: float,
T: float,
use_quantum: bool,
eta: float,
k_cut: float,
) -> float:
"""
3次（K3）のBubbleシフト Δω_q^(3) を計算する。
詳細説明:
Δω_q^(3) = Re Σ^(3)(q, ω_q) / (2 ω_q)
:param q: float, 波数q
:param kgrid: np.ndarray, 波数kのグリッド
:param a: float, 格子定数
:param M: float, 質量
:param K2: float, 調和力の定数
:param K3: float, 3次非調和定数
:param T: float, 温度
:param use_quantum: bool, 量子論的計算を使用するかどうか
:param eta: float, 正則化パラメータ
:param k_cut: float, 赤外発散を防ぐためのkのカットオフ
:returns: float, 3次周波数シフト量 Δω_q^(3)
"""
omega_q = float(omega_harmonic(np.array([q]), a, M, K2)[0])
if omega_q <= 0.0:
return 0.0
reS = re_sigma_cubic_bubble(q, kgrid, a, M, K2, K3, T, use_quantum, eta, k_cut)
return reS / (2.0 * omega_q)
# -----------------------------
# Main calculation
# -----------------------------
[ドキュメント]
@dataclass
class Params:
"""
計算に使用する物理パラメータや条件を保持するデータクラス。
詳細説明:
格子定数、質量、k点数、カットオフ、量子論フラグ、正則化係数などをまとめている。
:param a: float, 格子定数
:param M: float, 質量
:param Nk: int, k点数
:param k_cut: float, 波数kのカットオフ
:param use_quantum: bool, 量子論的計算を使用するかどうか
:param eta_rel: float, 正則化パラメータの相対値
"""
a: float
M: float
Nk: int
k_cut: float
use_quantum: bool
eta_rel: float
[ドキュメント]
def make_kgrid(Nk: int, a: float) -> np.ndarray:
"""
第1ブリュアンゾーン内の均一なk点グリッドを作成する。
詳細説明:
区間 (-pi/a, pi/a] における均一なk点グリッドを生成する。
:param Nk: int, k点グリッドの数
:param a: float, 格子定数
:returns: np.ndarray, k点グリッドの配列
"""
k = np.linspace(-math.pi / a, math.pi / a, Nk, endpoint=False)
return k
[ドキュメント]
def omega_T_for_q(
q: float,
T: float,
kgrid: np.ndarray,
p: Params,
k2f: Callable[[float, float | None], float],
k3f: Callable[[float, float | None], float],
k4f: Callable[[float, float | None], float],
) -> Tuple[float, float, float, float]:
"""
指定された波数qと温度Tに対する振動数および各シフト量を計算する。
:param q: float, 波数q
:param T: float, 温度
:param kgrid: np.ndarray, 波数kのグリッド
:param p: Params, 計算パラメータを保持するデータクラス
:param k2f: Callable, K2を計算する関数
:param k3f: Callable, K3を計算する関数
:param k4f: Callable, K4を計算する関数
:returns: Tuple[float, float, float, float], (初期振動数, 3次シフト量, 4次シフト量, シフト後振動数)
"""
K2 = float(k2f(0.0, T))
K3 = float(k3f(0.0, T))
K4 = float(k4f(0.0, T))
omega0 = float(omega_harmonic(np.array([q]), p.a, p.M, K2)[0])
# regularization scale for principal value
eta = max(p.eta_rel * omega0, 1e-12)
d3 = 0.0
d4 = 0.0
if abs(K3) > 0.0 and omega0 > 0.0:
d3 = delta_omega_cubic_bubble(q, kgrid, p.a, p.M, K2, K3, T, p.use_quantum, eta, p.k_cut)
if abs(K4) > 0.0 and omega0 > 0.0:
d4 = delta_omega_quartic_tadpole(q, kgrid, p.a, p.M, K2, K4, T, p.use_quantum, p.k_cut)
omegaT = omega0 + d3 + d4
return omega0, d3, d4, omegaT
[ドキュメント]
def parse_args() -> argparse.Namespace:
"""
コマンドライン引数をパースする。
:returns: argparse.Namespace, パースされた引数
"""
ap = argparse.ArgumentParser(description="Anharmonic phonon frequency shift vs T (1D chain, K2/K3/K4 as functions).")
ap.add_argument("--a", type=float, default=1.0, help="Lattice constant a [m] (or your unit).")
ap.add_argument("--M", type=float, default=1.0, help="Mass M [kg] (or your unit).")
ap.add_argument("--Nk", type=int, default=4096, help="Number of k-points in BZ.")
ap.add_argument("--kcut", type=float, default=0.0, help="Exclude |k| < kcut to avoid IR issues (1D). Unit: 1/a.")
ap.add_argument("--q", type=float, default=0.5, help="Target q as fraction of pi/a (q = frac * pi/a).")
ap.add_argument("--Tmin", type=float, default=0.0, help="Minimum temperature [K].")
ap.add_argument("--Tmax", type=float, default=300.0, help="Maximum temperature [K].")
ap.add_argument("--nT", type=int, default=31, help="Number of temperature points.")
ap.add_argument("--classical", action="store_true", help="Use classical (high-T) occupations instead of quantum.")
ap.add_argument("--eta_rel", type=float, default=1e-3, help="Regularization eta = eta_rel * omega0 for bubble principal value.")
ap.add_argument("--out", type=str, default="omega_vs_T.csv", help="Output CSV filename.")
ap.add_argument("--thz", action="store_true", help="Output frequencies in THz (cycles/s) instead of rad/s.")
return ap.parse_args()
[ドキュメント]
def main() -> None:
"""
メインの実行関数。
詳細説明:
コマンドライン引数を読み込み、計算パラメータを設定して各温度における
フォノン振動数のシフトを計算し、CSVファイルとして出力する。
:returns: None
"""
args = parse_args()
a = args.a
M = args.M
Nk = args.Nk
q = args.q * (math.pi / a) # q in 1/m (or 1/a)
k_cut = args.kcut * (1.0 / a) if args.kcut > 0 else 0.0
p = Params(
a=a,
M=M,
Nk=Nk,
k_cut=k_cut,
use_quantum=(not args.classical),
eta_rel=args.eta_rel,
)
kgrid = make_kgrid(Nk, a)
Ts = np.linspace(args.Tmin, args.Tmax, args.nT)
rows: List[List[float]] = []
for T in Ts:
omega0, d3, d4, omegaT = omega_T_for_q(q, float(T), kgrid, p, k2_func, k3_func, k4_func)
rows.append([float(T), omega0, d3, d4, omegaT])
arr = np.array(rows, dtype=float)
# unit conversion for output
# rad/s -> THz (cycles/s) = omega/(2π)/1e12
if args.thz:
conv = 1.0 / (TWOPI * 1e12)
arr[:, 1:] *= conv
header = "T,omega0,domega3_bubble,domega4_tadpole,omegaT"
np.savetxt(args.out, arr, delimiter=",", header=header, comments="")
print(f"Wrote: {args.out}")
print("Columns:", header)
if args.thz:
print("Frequencies are in THz.")
else:
print("Frequencies are in rad/s.")
if __name__ == "__main__":
main()