#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
:doc:`anharmonic_usage`
anharmonic_chain_Tshift.py
1次元単原子チェーン（最隣接）の非調和フォノン振動数シフトの温度依存性を計算するモジュール。
詳細説明:
モデル（相対変位 Δ_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は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を返す。
詳細説明:
デフォルトでは定数を返しますが、ユーザーが任意の関数に置き換えることができます。
単位は [N/m]、または無次元の慣例に従う場合は一貫した単位となります。
:param delta: float: 相対変位
:param T: float | None: 温度 [K]
:returns: float: 調和力定数 K2
"""
return 10.0 # example
[ドキュメント]
def k3_func(delta: float = 0.0, T: float | None = None) -> float:
"""
3次の非調和力定数K3を返す。
詳細説明:
デフォルトでは0を返しますが、バブル型シフトを有効にするためにゼロ以外の値に設定することができます。
:param delta: float: 相対変位
:param T: float | None: 温度 [K]
:returns: float: 3次の非調和定数 K3 [N/m^2]
"""
return 0.0 # example (set nonzero to enable bubble shift)
[ドキュメント]
def k4_func(delta: float = 0.0, T: float | None = None) -> float:
"""
4次の非調和力定数K4を返す。
詳細説明:
デフォルトの値を返しますが、オタマジャクシ型シフトの挙動を変更するために任意の関数に置き換え可能です。
:param delta: float: 相対変位
:param T: float | None: 温度 [K]
:returns: float: 4次の非調和定数 K4 [N/m^3]
"""
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:
"""
位相因子の絶対値の2乗 |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(ω) を計算する。
詳細説明:
量子的な場合は n(ω) = 1/(exp(β ħ ω) - 1) を計算します。
use_quantum=False の場合は、古典極限 n ~ kT/(ħ ω) を計算します。
:param omega: np.ndarray: 振動数の配列
:param T: float: 温度 [K]
: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: 温度 [K]
:param use_quantum: bool: 量子分布を使用するかどうか（デフォルトはTrue）
:returns: np.ndarray: 計算結果の配列
"""
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) のオタマジャクシ型シフト Δω_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: 対象の波数
:param kgrid: np.ndarray: 積分用のkグリッド配列
:param a: float: 格子定数
:param M: float: 質量
:param K2: float: 調和力定数
:param K4: float: 4次の非調和力定数
:param T: float: 温度 [K]
:param use_quantum: bool: 量子分布を使用するかどうか
:param k_cut: float: 赤外発散を避けるためのkのカットオフ値
:returns: float: オタマジャクシ型シフトによる振動数の補正量
"""
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次相互作用の要素 |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: 対象の波数
:param k: np.ndarray: 積分用の波数配列
: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|^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:
"""
3次 (K3) のバブル型シフトによる自己エネルギーの実部 Re Σ^(3)(q, ω_q) を計算する。
詳細説明:
バブル型の標準的な数式構造を用い、etaで正則化された主値積分 Re[1/(x + iη)] = x/(x^2 + η^2) を用いて評価します。
Σ^(3)(q, ω) = Σ_k |V3|^2 [ (n_k + n_{q-k} + 1)/(ω - ω_k - ω_{q-k} + i0)
+ (n_k - n_{q-k})/(ω - ω_k + ω_{q-k} + i0)
+ (n_{q-k} - n_k)/(ω + ω_k - ω_{q-k} + i0)
+ (n_k + n_{q-k} + 1)/(ω + ω_k + ω_{q-k} + i0) ]
:param q: float: 対象の波数
:param kgrid: np.ndarray: 積分用のkグリッド配列
:param a: float: 格子定数
:param M: float: 質量
:param K2: float: 調和力定数
:param K3: float: 3次の非調和力定数
:param T: float: 温度 [K]
:param use_quantum: bool: 量子分布を使用するかどうか
:param eta: float: 正則化のための微小パラメータ
:param k_cut: float: 赤外発散を避けるためのkのカットオフ値
:returns: float: 自己エネルギーの実部
"""
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:
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) のバブル型シフト Δω_q^(3) を計算する。
詳細説明:
自己エネルギーの実部を用いて、Δω_q^(3) = Re Σ^(3)(q, ω_q) / (2 ω_q) を算出します。
:param q: float: 対象の波数
:param kgrid: np.ndarray: 積分用のkグリッド配列
:param a: float: 格子定数
:param M: float: 質量
:param K2: float: 調和力定数
:param K3: float: 3次の非調和力定数
:param T: float: 温度 [K]
:param use_quantum: bool: 量子分布を使用するかどうか
:param eta: float: 正則化のための微小パラメータ
:param k_cut: float: 赤外発散を避けるためのkのカットオフ値
:returns: float: バブル型シフトによる振動数の補正量
"""
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:
"""
計算に使用するパラメータを保持するデータクラス。
: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] となります。
: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: 対象の波数
:param T: float: 温度 [K]
:param kgrid: np.ndarray: 積分用のkグリッド配列
:param p: Params: 計算パラメータのデータクラス
:param k2f: Callable[[float, float | None], float]: 調和力定数を返す関数
:param k3f: Callable[[float, float | None], float]: 3次の非調和力定数を返す関数
:param k4f: Callable[[float, float | None], float]: 4次の非調和力定数を返す関数
: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:
"""
スクリプトのエントリポイント。
詳細説明:
パラメータの設定、kグリッドの作成、各温度に対する振動数シフトの計算を実行し、
結果を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()