4. 速度自相关函数与振动谱

速度自相关函数与振动谱：从轨迹到频谱

速度自相关函数（VACF）描述粒子速度随时间的“记忆衰减”，其傅里叶变换即振动谱（态密度）。对LJ流体，VACF在短时间呈现振荡衰减，反映粒子在配位壳层内的准局域振动；长时间趋于零，反映扩散导致的关联丢失。振动谱中低频峰对应集体声子模式，高频肩峰对应局域振动。

关键参数与计算设置

VACF计算与谱分析实现

import numpy as np

# 假设已从轨迹读取速度数组 vel，形状 (n_frames, N, 3)
# 采样间隔 dt_sample = 0.02 (约化时间单位)

def compute_vacf(vel, max_lag=None):
    """时间原点平均计算归一化VACF"""
    n_frames, N, _ = vel.shape
    if max_lag is None:
        max_lag = n_frames // 2
    vacf = np.zeros(max_lag)
    for lag in range(max_lag):
        # 对所有时间原点 t 和所有粒子 i 求 v(t)·v(t+lag)
        dot = np.sum(vel[:n_frames-lag] * vel[lag:], axis=(1, 2))
        vacf[lag] = dot.mean()
    vacf /= vacf[0]          # 归一化 C(0)=1
    return vacf

def compute_spectrum(vacf, dt_sample):
    """加Hann窗后FFT，返回频率与功率谱"""
    n = len(vacf)
    window = np.hanning(n)
    spec = np.fft.rfft(vacf * window)
    freq = np.fft.rfftfreq(n, d=dt_sample)   # 单位: 1/τ
    power = np.abs(spec)**2
    return freq, power

# 主流程
vacf = compute_vacf(vel)
freq, power = compute_spectrum(vacf, dt_sample=0.02)

# 峰位识别
from scipy.signal import find_peaks
peaks, props = find_peaks(power, height=0.05*power.max(), distance=5)
print("峰频率 (1/τ):", freq[peaks])
print("峰强度:", props["peak_heights"])

结果判读

• VACF振荡周期：第一个负谷位置对应粒子在配位壳层内的反弹时间，约 0.3–0.5 τ。
• 低频峰（< 2 /τ）：声子模式，频率随密度增大而升高。
• 高频肩峰（约 8–12 /τ）：局域振动，对应粒子在势阱内的瞬时振荡。
• 谱线平滑度：关联长度不足会导致谱线毛刺，建议 max_lag ≥ 1000 帧。

验证与常见问题