ORCA双原子势能曲线与局部振动

想从一组“键长—能量”数据估计双原子分子的平衡位置和局部振动，可以先做固定键长的单点计算，再用 orca_fitpes 处理。这篇采用鸢算科技集群实际完成的H₂ RHF/STO-3G扫描：从0.45到1.65 Å，间隔0.05 Å，共25点。

这里使用小基组帮助理解操作。数值属于这套计算模型，不作为高精度实验键长、解离能或谱线的预测。

1. 下载后先跑后处理#

下载双原子势能曲线实例包，在新目录执行：

python3 -m zipfile -e orca-diatomic-spline-wave528-*.zip orca-diatomic-spline
cd orca-diatomic-spline
bash run.sh

这个命令复用随包25个实际单点结果，重新运行H₂、D₂与bohr坐标三组原生后处理，然后生成表格和图片；它不重新做25次电子结构计算。

结果写到 work、tables 和 figures。只读取随包原始结果时执行 bash analyze.sh，重新分析本次后处理用 bash analyze.sh work。

ORCA6.1.1为官方CPU预编译发行包，配套OpenMPI4.1.8由GCC13.3构建。本例单进程、单线程，不使用GPU；分析使用已有ASE3.29.0环境中的Python、NumPy/SciPy与Matplotlib。

2. 原来的25个点怎么算#

每个点固定一个H–H距离，只计算该几何下的电子能量。例如0.75 Å的完整输入：

! RHF STO-3G VeryTightSCF NoAutoStart
%pal nprocs 1 end
%maxcore 300
* xyz 0 1
H 0 0 0
H 0 0 0.7500000000
*

RHF用于这个中性闭壳层示例，* xyz 0 1 表示电荷0、单重态。这里只变键长，不在每个点做结构优化；若加上 Opt 把每个点都放开优化，就不再是这里的固定距离扫描。

在已经分配的计算资源内，单独重做这个点：

module purge
module load orca/6.1.1
export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1
mkdir single-point
cp inputs/H2-single.inp single-point/case.inp
cd single-point
"$ORCA_ROOT/orca" case.inp > case.out
grep -E 'SCF CONVERGED AFTER|FINAL SINGLE POINT ENERGY|ORCA TERMINATED NORMALLY' case.out
cd ..

其余点保留在 reference/STO-3G-0.45/ 到 reference/STO-3G-1.65/，每个目录都带完整输入、实际输出和属性文本。自己的扫描应保持方法、基组、电荷、自旋和收敛设置一致，并检查每一个点正常完成。

3. fitpes输入怎么写#

本例H₂输入是 inputs/H2-A.inp，完整内容如下：

1 25 1 1.00782503223 1.00782503223
0.45000000000000 -0.98751314072000
0.50000000000000 -1.04299627675200
0.55000000000000 -1.07905073771900
0.60000000000000 -1.10112824324500
0.65000000000000 -1.11299654616500
0.70000000000000 -1.11734903508000
0.75000000000000 -1.11615144868100
0.80000000000000 -1.11085039692300
0.85000000000000 -1.10251055311700
0.90000000000000 -1.09191403999200
0.95000000000000 -1.07963692699000
1.00000000000000 -1.06610864792700
1.05000000000000 -1.05165675410400
1.10000000000000 -1.03653887335300
1.15000000000000 -1.02096414178600
1.20000000000000 -1.00510670467200
1.25000000000000 -0.98911381210700
1.30000000000000 -0.97311061372000
1.35000000000000 -0.95720317295800
1.40000000000000 -0.94148065254300
1.45000000000000 -0.92601717353800
1.50000000000000 -0.91087355237000
1.55000000000000 -0.89609895608600
1.60000000000000 -0.88173244770300
1.65000000000000 -0.86780438176100

首行的5个数依次表示：势能曲线条数、采样点数、坐标是否为Å、原子1质量、原子2质量。这里 1 25 1 表示一条曲线、25点、坐标用Å；两个质量单位为u。

后面每一行都是 同一个点的距离R、电子能量E。电子能量始终用Hartree。按距离递增排列，保留足够小数；不要先写一整列距离，再另写一整列能量，也不要把相对eV能量直接填进Hartree列。

手动处理这份文件：

cp inputs/H2-A.inp manual-fit.inp
"$ORCA_ROOT/orca_fitpes" manual-fit.inp > manual-fit.out
sed -n '1,20p' manual-fit.out

4. 先看SPLINE局部结果#

本篇读取输出前部的插值最低点、局部二阶导数，以及标有 based on SPLINE 的频率，不把后面其他模型的数值混在一起。

最低采样点与插值最低点不必重合：本例相邻采样距离是0.70和0.75 Å，插值在两者之间寻找局部最低点。最低点处的一阶导数接近0，二阶导数为正；曲率决定局部势阱的“硬度”。

这条小基组曲线只在采样范围内作参考。最低点附近需要更精确结果时，应选择更合适的电子结构方法和基组，并在最低点附近加密实际采样；插值不会自动补足电子结构方法本身的精度。

5. H₂换成D₂，能量面需要重算吗#

在这里的Born–Oppenheimer电子势能与谐近似下，替换核质量不重新改变这条电子能量曲线。随包D₂输入保留相同的25组距离和能量，只把两端质量改成 2.01410177812 u。

Å与bohr也有配套对照：H2-bohr.inp 的第三个数为0，距离已除以0.5291772083转换为bohr，电子能量和质量保持一致。

head -n 4 inputs/D2-A.inp
head -n 4 inputs/H2-bohr.inp
head work/fit-STO-3G-D2-A/native.log

局部谐振频率与约化质量μ的平方根成反比，μ=m₁m₂/(m₁+m₂)。同一条曲率下，D₂更重，振动波数降低；单位转换则不应改变物理结果。

图中的水平线使用本次局部谐模型的 (v+1/2)hν 绘制，展示质量与零点能的关系；它不是新的非谐计算，也不代表完整解离过程中的全部振动能级。

6. 最后怎样导出和画图#

head tables/scan-energies.csv
cat tables/local-minima.csv
bash analyze.sh work

scan-energies.csv 保留25个距离、输出能量和属性文件的高精度能量；local-minima.csv 保留三组实际原生插值结果及独立频率对照。figures 中是结构与扫描图、局部曲率图和同位素谐振能级图。

原生日志还会列出Chebyshev、数值导数和Morse模型等其他结果。不同近似不能只因出现在同一份输出里就互相替代。本例采用上面已经与原始采样点核对的SPLINE局部量；需要使用Morse解离深度或非谐参数时，应另行检查模型残差、采样区间和参数物理意义，不能把负深度、非数值输出或有限区间外推直接当作可信的解离能。

orca_fitpes 的用途见官方工具目录。长度与能量换算见ORCA单位说明。H/D质量取自NIST同位素质量表。