Calorine红外与拉曼谱实操

下载模型、真实数据、完整脚本与参考图片。相关教程：偶极与极化率 · 红外与拉曼谱 · NEP89微调 · 模型编辑 · GPUMD数据读写 · 完整接口索引

1 准备与运行#

先解压并重画已经完成的参考结果：

unzip calorine-response-workflows-wave266-*.zip
cd calorine-response-workflows-wave266
bash run.sh

六张参考图在figures，CSV在analysis。若要自己重新计算，在新解压目录以普通账号执行：

bash prepare.sh
bash calculate.sh
bash status.sh

calculate通过Slurm分配一张P100，自动完成本包所有计算；新结果在work，参考图与新算结果分开。已提交时继续查看status，不要删除记录重复提交。结束后：

cat work/COMPLETE.json
ls work/analysis

本例环境为Calorine 4.0、Python 3.11、GPUMD 5.9.1；GPUMD采用GCC 13.3/CUDA 12.9系列环境。可以用module show calorine/4.0和module show gpumd/5.9.1-p100查看当前模块。下面的Python片段在计算完成后的work目录执行；MD和训练由calculate提交。

2 这条工作流具体做什么#

先用能量势计算原子运动，再逐帧用偶极模型和极化率模型计算响应，最后把响应时间序列转成相关函数和频谱。分子响应模型不负责给出能量和原子力。

本例是3原子的孤立H2O、NVE、2 ps。初始内部动能0.01/0.03 eV，各用种子13579/24680，时间步0.25 fs；另用0.125 fs作同一初态的对照。每1 fs保存一次响应，各2001帧。这不是恒温300 K轨迹。较小时间步下能量漂移下降，完整数据均保留。

顺序是：准备模型与结构 → 分配资源 → 运行五条轨迹 → 逐帧预测μ和α → 计算相关函数与窗函数 → 得到IR/Raman CSV → 画图。

3 找到自己的时间序列#

ls work/MD
head -n 5 work/analysis/KE0.03-s13579-dt0.25-IR.csv
head -n 5 work/analysis/KE0.03-s13579-dt0.25-Raman.csv

每条轨迹的MD-RAW.json保留原始μ、α、能量、几何和时间；trajectory.extxyz保留真实结构。根目录reference/MD是已经完成的参考轨迹。

以下是提取和计算频谱的核心操作：

cd work
python - <<'CALORINE_SPECTRUM'
import json
import numpy as np
import pandas as pd
from calorine.tools import get_ir_spectrum, get_raman_spectrum
rows = json.load(open('MD/KE0.03-s13579-dt0.25/MD-RAW.json'))['rows']
t = [r['time_ps'] for r in rows]
mu = np.array([r['mu_eA'] for r in rows])
a = np.array([r['alpha_A3'] for r in rows])
signal = pd.DataFrame(dict(time=t, mu_x=mu[:,0], mu_y=mu[:,1], mu_z=mu[:,2]))
pol = pd.DataFrame(dict(time=t, xx=a[:,0,0], yy=a[:,1,1], zz=a[:,2,2],
                        xy=a[:,0,1], xz=a[:,0,2], yz=a[:,1,2]))
ir = get_ir_spectrum(signal, column='mu', t_sigma=0.2, window_size=1.0)
ra = get_raman_spectrum(pol, t_sigma=0.2, window_size=1.0)
ir.to_csv('my-IR.csv', index=False)
ra.to_csv('my-Raman.csv', index=False)
CALORINE_SPECTRUM

time、t_sigma和window_size均使用ps。这里截取1 ps相关窗口、Gaussian宽0.2 ps；谱的横坐标直接使用wavenumber_invcm，单位cm⁻¹。angular_frequency是角频率的10¹² s⁻¹表示，不能把它直接当普通频率THz。

以0.03 eV、种子13579、0.25 fs那条轨迹为例，红外两个区间峰在约1668.65与3988.08 cm⁻¹，拉曼强峰约3904.65 cm⁻¹。这些来自当前能量势和响应势，不是实验水谱；轨迹、训练范围和经典采样都会影响结果。有限窗口的小负尾没有被裁掉；展示时按区间内最大绝对值归一化，原始CSV保留原数值。

4 窗口、极化和量子修正#

更长轨迹、合适窗宽和多个初态有助于判断谱峰稳定性，不能只把线画得平滑。polarization可以指定IR方向；Raman的polarization_in/out分别指定入射和出射方向。apply_quantum_correction(ir, temperature=300, column='ir_intensity')添加ir_intensity_qm列：这里300 K是修正参数，不能因此把原NVE轨迹说成300 K恒温模拟。

5 介电函数和折射率接口怎么用#

分子线形首先保持不按体积归一化。资料包还给出了明确选定24³ ų盒和T=300的接口/单位演示，不能解释成体相水的介电常数。当前KK接口返回响应实部χ₁；转换时先加入真空基线1。弱吸收下κ的相减公式会损失小量，配套脚本用复数平方根保留符号：

from calorine.tools import (get_dielectric_function, apply_kramers_kronig,
    get_refractive_index, get_absorption_coefficient,
    get_dielectric_function_from_refractive_index)
eps = get_dielectric_function(signal, volume=24**3, temperature=300,
        column='mu', t_sigma=0.2, window_size=1.0, return_real_part=False)
kk = apply_kramers_kronig(eps)
iso = pd.DataFrame(dict(angular_frequency=kk.angular_frequency,
    epsilon_real=1 + sum(kk['epsilon_real_'+s] for s in ['xx','yy','zz'])/3,
    epsilon_imag=sum(kk['epsilon_imag_'+s] for s in ['xx','yy','zz'])/3))
nr = get_refractive_index(iso)
nr['refractive_index_imag'] = np.sqrt(
    iso.epsilon_real.to_numpy() + 1j*iso.epsilon_imag.to_numpy()).imag
optics = get_absorption_coefficient(nr)
back = get_dielectric_function_from_refractive_index(optics)

这里只对各向同性迹使用标量n+iκ代数，不把张量的每个非对角分量各自当成一个折射率。原始负谱值保留，不通过裁零制造正吸收；电子背景、取向平均和完整热力学采样也不在这条短轨迹中。

参考Calorine光谱与工具接口。