ORCA静电势拟合电荷实操

需要用原子点电荷近似分子外部的静电势时，可以考虑CHELPG拟合。它与Mulliken、Löwdin电荷的定义不同。这篇用水分子和CF3Cl的实际结果，说明怎样运行、读取文件、改变拟合条件，以及为什么原子电荷不能代替完整的空间静电势。

1. 下载并运行#

下载ORCA电荷与片段轨道实例包，在已分配的CPU计算资源中执行：

python3 -m zipfile -e orca-charges-fragments-wave520-*.zip orca-charges-fragments
cd orca-charges-fragments
bash run.sh

完整流程计算四个起始分子、两组默认CHELPG电荷、四组拟合参数对照，并完成片段合并与整体SCF。再次运行完整流程时，请解压到新目录。

环境为 orca/6.1.1，采用官方预编译ORCA6.1.1和配套OpenMPI4.1.8（GCC13.3构建）。本例每次使用1个CPU进程，无需GPU；ORCA本体未重新编译。通常的环境与提交方法见ORCA使用指南。

inputs/water-HF.inp 使用HF/def2-SVP中性单重态水分子，inputs/CF3Cl.inp 使用PBE/def2-SVP中性单重态CF3Cl；均为固定几何单点。完整输入均在包内，水分子几何和输入也见局域轨道实例，CF3Cl完整输入见指定位置静电势。

2. 对已有结果运行CHELPG#

已有波函数时，不必重做单点。在相应目录中运行：

module purge
module load orca/6.1.1
cd work/charges-water-HF
"$ORCA_ROOT/orca_chelpg" case.gbw > chelpg-manual.out
cat case.pc_chelpg
cd ../..

同一计算的 case.gbw、case.densities、case.densitiesinfo 保留在一起。本次默认设置为网格间距 0.3 Å、截断距离 2.8 Å、COSMO半径，不约束偶极矩。

水分子的原生输出文件开头为：

3
# orca_chelpg point-charges stored as xyzfile
Q    -0.772054    -1.551007    -0.114520     0.000000
Q     0.386764    -1.934259     0.762503     0.000000
Q     0.385290    -0.599677     0.040712     0.000000

第一行是点数，第二行是说明，之后每行依次为 Q、电荷、x、y、z。电荷单位为e，坐标单位为Å，行序与原结构一致。这是点电荷列表，不是普通的“元素、x、y、z”分子XYZ格式；传给其他程序时，应按目标程序要求转换。

查看另一组实际输出：

cat work/charges-CF3Cl/case.pc_chelpg

本次CF3Cl的C、Cl、F1、F2、F3电荷分别为 0.276388、−0.053020、−0.074644、−0.074242、−0.074482 e，总电荷为0。有限网格下等价位置的拟合电荷可能有小差异，不应只凭最后几位数推断原子有不同的物理电荷。

3. 四个拟合参数怎么改#

对已经收敛、同几何同方法同基组的波函数，可以用 MOREAD NoIter CHELPG 读取轨道并重新拟合，无需再迭代优化轨道。下面实际演示加上偶极矩约束：

module purge
module load orca/6.1.1
mkdir manual-water
cp work/seed/water-HF/case.gbw manual-water/seed.gbw
cd manual-water
cat > charge.inp <<'EOF'
! HF def2-SVP TightSCF NoSym NoRI MOREAD NoIter CHELPG NoAutoStart
%pal nprocs 1 end
%maxcore 1024
%moinp "seed.gbw"
%chelpg
  GRID 0.3
  RMAX 2.8
  VDWRADII COSMO
  DIPOLE true
end
* xyz 0 1
O -1.551007000 -0.114520000 0.000000000
H -1.934259000  0.762503000 0.000000000
H -0.599677000  0.040712000 0.000000000
*
EOF
"$ORCA_ROOT/orca" charge.inp > charge.out
cat charge.pc_chelpg
cd ..

%moinp读取的是 seed.gbw，输出主名是 charge，两者明确分开。这里的 NoIter 仅用于保持这次已收敛波函数；改变体系、方法、基组或电子状态时，应按需要重新做SCF。

每组只改变一项，完整输入保存在对应目录。把上面参数块改成同样的设置即可操作，不需要额外安装插件。

保持同一电子密度时，默认、较密网格、较大截断、BW半径、偶极矩约束对应的氧电荷依次为 −0.772054、−0.772243、−0.770493、−0.769119、−0.755775 e。这说明拟合条件本身就会改变点电荷，不等于分子电子结构发生了变化。

4. 拟合电荷能代替空间静电势吗#

左图用与静电势点值教程相同的九个位置，比较量子静电势与本次CHELPG原子点电荷的库仑势。这些是独立比较点，不是拟合网格。两者存在明显差别，尤其是Cl轴向的正静电势区域。

所以，Cl的拟合原子电荷为负，并不排除它的外侧存在正静电势区域；一个原子中心的单点电荷无法完整表示各方向的电子分布。分析σ-hole等空间性质时，继续使用实际波函数和Critic2表面静电势方法；将电荷用于其他模型时，还应采用那个模型相应的参数和约定。

5. 核对并画图#

bash analyze.sh work
/public/toolkit/miniconda3/envs/ase-3.29.0/bin/python plot_results.py

第一条核对本次完整结果，包括总电荷、轨道密度、实际参数和偶极矩约束。第二条从包内参考数据生成 figures/ 下的三张PNG；换成自己的数据时，相应修改画图脚本的数据目录。

接口和参数参考ORCA CHELPG说明。