ORCA 点电荷嵌入实操

下载完整输入、真实结果、轨道导出、CSV与绘图脚本。解压进入资料包目录，先处理已有数据：

bash analyze.sh

图片在 reference/figures，数据表在 reference/analysis；原生Molden轨道文件在相应原子算例目录。想重新计算全部条件，在新解压目录执行：

bash submit.sh

脚本提交27个独立条件的数组任务，结束后导出6组原子轨道，再自动提取和绘图。新结果在 work，原参考保留。它们分属溶剂模型、点电荷环境和原子SOC三个小工作流，不是27种独立材料性质。

使用 orca/6.1.1 厂商预编译CPU程序，依赖专用OpenMPI4.1.8（GNU13.3环境）；单条件 nprocs=1。实际在P100节点使用CPU，不是ORCA GPU加速。附件Slurm模板每个任务申请1个CPU核心和1张P100，GPU不参与计算，资源按对应分区计费。不要在登录节点运行模拟。后处理使用ASE3.29.0环境中的NumPy/SciPy/Matplotlib。

1 先做一个能独立核对的小环境#

量子区是一份固定水分子，外部放两个相反的点电荷，位置分别为z=+4 Å和−4 Å。改变电荷大小和符号，观察量子电子的极化与能量变化。

这里是固定点电荷的静电嵌入示例，没有环境极化迭代，也不是完整的QM/MM力场、键切割或耦合动力学流程。

2 输入文件和外部电荷文件#

完整ORCA输入：

! PBE0 def2-SVP def2/J RIJCOSX DefGrid3 VeryTightSCF EnGrad
%pal nprocs 1 end
%maxcore 500
%pointcharges "charges.pc"
* xyz 0 1
O 0 0 0.119262
H 0 0.763239 -0.477047
H 0 -0.763239 -0.477047
*

charges.pc 实际内容：

2
0.1000000000 0 0 4.0000000000
-0.1000000000 0 0 -4.0000000000

第一行是点电荷数量。后面每行依次为：电荷q（单位e）、x、y、z（单位Å）。坐标不是Bohr；水分子量子区的输入电荷仍为0、多重度1，不要把外部电荷不加判断地计入量子区电子数。

单条件原生命令：

module load orca/6.1.1
"$ORCA_ROOT/orca" case.inp > case.out
grep 'FINAL SINGLE POINT ENERGY' case.out
cat case.pcgrad

附件 bash run.sh 22 对应该+0.1/−0.1 e、±4 Å条件；在已分配资源中运行。EnGrad 同时生成量子原子梯度与点电荷梯度。

3 外部电荷彼此的能量默认不在里面#

外部文件点电荷的相互作用通常交给外部MM程序处理，因此默认不会把点电荷彼此的经典库仑能加入这里的核排斥能。若独立模型需要这项，可明确加入 DoEQ：

! PBE0 def2-SVP def2/J RIJCOSX DefGrid3 VeryTightSCF EnGrad DoEQ
%pal nprocs 1 end
%maxcore 500
%pointcharges "charges.pc"
* xyz 0 1
O 0 0 0.119262
H 0 0.763239 -0.477047
H 0 -0.763239 -0.477047
*

同样的量子结构、电子设置和电荷排列，仅改变是否加入这项常数。不要在完整QM/MM程序已计算该能量时再重复加入。

本例两电荷相距8 Å，经典项为：

E_pc-pc = q1*q2 / r_bohr
        = -0.000661471513 Eh

实际 DoEQ 能量变化约−0.000661471510 Eh，与独立库仑公式相符。加入这项不改变本例水分子的电子极化和量子原子梯度，因为点电荷位置固定、该常数与量子坐标无关。

4 实际电荷扫描和位置对照#

下表能量相对零点电荷/气相条件，单位kJ/mol。q是+R处电荷，−R处为−q：

零电荷条件重现气相结果。正负电荷扫描显示方向相关的能量变化和电子极化；这是当前非均匀点电荷场中的响应，不能直接把偶极变化除以一个近似场强就称作精确孤立分子极化率。

5 梯度的方向和单位#

case.pcgrad 给出能量对外部点电荷位置的导数，单位Eh/Bohr；力是它的负号。不要把Å位移直接与Eh/Bohr数值相乘而漏掉单位换算。

本例把两电荷从±R一起移动到±(R±0.005 Å)，对最终能量作中心差分，并与原生梯度比较：

dE/dR = [E(R+h)-E(R-h)] / (2h)
      = [g_z(+R)-g_z(-R)] / Bohr_to_Angstrom

本次数值差分约0.00133195535 Eh/Å，原生梯度组合约0.00133194686 Eh/Å，相对差约6.4×10⁻⁶。它用于核对当前文件格式、单位和梯度方向，不是对任意外部环境给出的物理精度保证。

完整控制数据在 reference/analysis/point-charge-models.csv 和RESULTS.json。使用外部环境时，还需明确其电荷来源、固定或可极化、是否计入相互作用、量子区边界与总能量定义。

原生格式与DoEQ约定见 ORCA外部点电荷说明。