关于用Gaussian做水分子的OH键的柔性扫描遇到的能量突然跳变的问题

gingaorb

各位小哥哥，各位大佬们，请教一个问题呀~
是这样的，我最近想探索水分子中的一个氢氧键在断裂的过程中有没有可能性均裂生成羟基自由基和氢原子，于是做了一个水分子的OH键的柔性扫描，其中，将一个优化好的水分子其中一个OH键对齐Z轴作为初始结构，H不断远离O做柔性扫描，以下是我的具体输入文件：
%chk=H2O.chk
%nprocshared=48
%mem=60GB
# opt=(modred) nosymm ub3lyp/def2tzvp  empiricaldispersion=gd3bj pop=(always,nbo)

Title Card Required

0 1
O        -0.092869        0.000000       -0.070943
H        -0.092869        0.000000        0.891633
H         0.835824        0.000000       -0.324087

B 1 2 S 60 0.1

然后扫描出来的图在第34帧时能量突然跳变，构型也突然发生了一个巨大变化（如下动图所示），显示2H原子上的NBO电荷和自旋布居也在这时发生突变，从氢离子突然跳变成氢原子，这太奇怪了。我不能理解这种跳变现象发生的原因，有没有哪位小哥哥，大佬能为我解答一下这种变化发生的原因呢？是否和站长大大在帖子“详谈使用Gaussian做势能面扫描 - 思想家公社的门口：量子化学·分子模拟·二次元  http://sobereva.com/474”中提到的原因一致呢？万分感谢大家了！

sobereva

这么扫描没意义
涉及共价键断裂情况的扫描怎么做看此文，对于DFT必须考虑对称破缺
详谈使用Gaussian做势能面扫描
http://sobereva.com/474（http://bbs.keinsci.com/thread-12660-1-1.html）

gingaorb

sobereva 发表于 2023-4-14 15:04
这么扫描没意义
涉及共价键断裂情况的扫描怎么做看此文，对于DFT必须考虑对称破缺
详谈使用Gaussian做势 ...

我也尝试了，将输入文件改为如下（每一步都用guess=mix，以产生对称破缺）：
%chk=H2O.chk
%nprocshared=48
%mem=60GB
# opt=(modred) nosymm ub3lyp/def2tzvp  empiricaldispersion=gd3bj pop=(always,nbo) guess(mix,always)

Title Card Required

0 1
O        -0.092869        0.000000       -0.070943
H        -0.092869        0.000000        0.891633
H         0.835824        0.000000       -0.324087

B 1 2 S 60 0.1

但是结果却在第39个点发生了SCF不收敛的报错（Convergence failure -- run terminated. Error termination via Lnk1e in /public1/home/scf0192/Group/senxu/soft/g16/l502.exe），并且无论我如何按照教程（量子化学计算中帮助几何优化收敛的常用方法 - 思想家公社的门口：量子化学·分子模拟·二次元  http://sobereva.com/164、解决SCF不收敛问题的方法 - 思想家公社的门口：量子化学·分子模拟·二次元  http://sobereva.com/61）尝试修改输入文件关键词都无济于事，如下是我最后修改的版本（所用高斯版本为G16，我是真没办法了才用的SCF=(XQC)，结果还不收敛，我技穷了呜呜~）：
%chk=H2O.chk
%nprocshared=48
%mem=60GB
# opt=(modred, MaxStep=5, NoTrustUpdate, GDIIS, calcall) nosymm ub3lyp/def2tzvp  empiricaldispersion=gd3bj pop=(always,nbo) guess=(mix,always) SCF=(XQC)

Title Card Required

0 1
O        -0.092869        0.000000       -0.070943
H        -0.092869        0.000000        0.891633
H         0.835824        0.000000       -0.324087

B 1 2 S 60 0.1

报错：Convergence failure.
Error termination via Lnk1e in /public1/home/scf0192/Group/senxu/soft/g16/l508.exe

而且我还有个疑问，就是按理说，水的OH键异裂生成氢离子和氢氧根离子会比较容易，均裂产生羟基自由基和氢原子应该很困难，并且就算能产生，体系的能量相比于生成氢离子和氢氧根应该会更高，但1#楼我发的能量曲线里，能量却突然降低，显得很反常理，这是为什么呢？感觉好伤脑筋呢o(╥﹏╥)o~

wjc404

这种扫描应该最好从最长长度开始，而不是最短键长。因为最长长度时情况更复杂，更容易出现对称破缺波函数。先确保最长键长下得到了正确的稳定波函数，再（读取初猜）往短走更稳妥。

你这个体系很小，就一个水分子，基组可以用def2-TZVPP，方法可以用比UKS更靠谱的NEVPT2（ORCA程序支持，这么小体系用数值梯度优化的开销可以接受）。

至于为什么是得到氢和氢氧自由基而不是两个离子，可以去算一下真空中H原子的电离能和氢氧自由基的电子亲和能，就知道了。水中的氢离子和氢氧根离子都是高度水合的，被充分稳定了，和真空不同。

gingaorb

wjc404 发表于 2023-4-18 03:22
这种扫描应该最好从最长长度开始，而不是最短键长。因为最长长度时情况更复杂，更容易出现对称破缺波函数。 ...

能请问下“真空中H原子的电离能和氢氧自由基的电子亲和能”这个用Gaussian要怎么计算呢？另外我也尝试了从最长长度开始，用# opt=(modred, MaxStep=5, NoTrustUpdate, GDIIS, calcall) nosymm ub3lyp/def2tzvp  empiricaldispersion=gd3bj pop=(always,nbo) guess=(mix,always) SCF=(XQC)，但第一步就几何优化不收敛。

Kmetsch

gingaorb 发表于 2023-4-18 21:00
能请问下“真空中H原子的电离能和氢氧自由基的电子亲和能”这个用Gaussian要怎么计算呢？另外我也尝试了 ...

从双自由基状态开始向内扫描最好使用片段组合波函数http://sobereva.com/474，对两个片段分别赋1个alpha和1个beta自旋，计算单点能后，再让扫描的任务读取自旋对称破缺单重态的波函数作为初猜。
真空电离能可以单纯的在同计算水平下对自由基和离子的能量作差，电子亲和能同理。
最后，建议去掉XQC，99.99%没用

wjc404

用ORCA在NEVPT2/def2-TZVPP水平下对水的一个O-H键做柔性扫描

步骤一、先在平衡键长附近生成初始轨道：
(1) 计算MP2自然轨道
(2) 以MP2自然轨道为初猜，做CASSCF(8,6)计算，查看活性空间轨道
以上步骤的输入文件这里暂不示出。假设第二步生成的轨道文件为H2O_INITCASV.gbw。
这里采用冻芯近似，剩下8个价电子全部放入活性空间。价层轨道O原子贡献4个（2S、2P）而每个H原子贡献一个（1S），故活性空间轨道数取6。

步骤二、读取初始轨道做柔性扫描：
输入文件如下

1. ! Opt NumGrad def2-TZVPP def2/JK RIJK ExtremeSCF MORead
   
2. %moinp "H2O_INITCASV.gbw"
   
3. %maxcore 2000
   
4. 
   
5. %casscf
   
6.   nel 8
   
7.   norb 6
   
8.   ptmethod fic_nevpt2
   
9.   ptsettings
   
10.     d3tpre 0
    
11.     d4tpre 0
    
12.   end
    
13. end
    
14. 
    
15. %geom
    
16.   scan
    
17.     B 0 1 [ 0.96 0.98 1.00 1.03 1.06 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.1 2.15 2.2 2.25 2.3 2.35 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 ]
    
18.   end
    
19. end
    
20. 
    
21. %numgrad dx 0.001 end
    
22. 
    
23. * int 0 1
    
24. O 0 0 0 0.0 0.0 0.0
    
25. H 1 0 0 0.96 0.0 0.0
    
26. H 1 2 0 0.97 102.0 0.0
    
27. *
    
28. 
    

复制代码

ORCA暂时未实现NEVPT2的解析梯度，所以需要用数值梯度来优化。计算数值梯度时，为了减小噪声的干扰，这里提高了收敛限和积分精度。以上计算使用RI-JK来加速。

计算结果经过处理后见下表

B(O-H1) scan (Angs.)	B(O-H2) opt (Angs.)	A(H1-O-H2) (Deg.)	E (kJ/mol)
0.96	0.966	102.3	-507.8
0.98	0.966	102.0	-507.5
1.00	0.966	101.7	-505.4
1.03	0.966	101.2	-499.4
1.06	0.966	100.7	-490.6
1.10	0.967	100.1	-475.5
1.15	0.967	99.2	-452.6
1.20	0.967	98.3	-426.8
1.25	0.968	97.4	-399.3
1.30	0.968	96.6	-370.9
1.35	0.969	95.7	-342.4
1.40	0.969	94.9	-314.4
1.45	0.970	94.0	-287.0
1.50	0.970	93.2	-260.8
1.55	0.971	92.4	-235.8
1.60	0.971	91.7	-212.3
1.65	0.972	90.9	-190.2
1.70	0.972	90.3	-169.8
1.75	0.972	89.6	-150.9
1.80	0.973	89.0	-133.6
1.85	0.973	88.5	-117.8
1.90	0.973	88.0	-103.6
1.95	0.973	87.5	-90.7
2.00	0.974	87.1	-79.2
2.05	0.974	86.7	-69.0
2.10	0.974	86.3	-59.9
2.15	0.974	85.9	-51.9
2.20	0.974	85.6	-44.9
2.25	0.974	85.4	-38.7
2.30	0.974	85.0	-33.4
2.35	0.974	84.8	-28.7
2.40	0.974	84.6	-24.6
2.50	0.974	84.2	-18.1
2.60	0.974	83.9	-13.3
2.70	0.974	83.6	-9.7
2.80	0.974	83.2	-7.1
2.90	0.974	79.9	-5.2
3.00	0.974	76.7	-3.8
3.10	0.974	73.5	-2.7
3.20	0.974	70.3	-2.0
3.30	0.974	67.2	-1.4
3.40	0.974	64.1	-1.0
3.50	0.974	61.1	-0.8
3.60	0.974	58.1	-0.6
3.70	0.974	55.1	-0.4
3.80	0.974	52.2	-0.3
3.90	0.974	49.3	-0.2
4.00	0.974	46.5	-0.2
INF	0.974	N.A.	0.0

gingaorb

wjc404 发表于 2023-4-28 21:27
用ORCA在NEVPT2/def2-TZVPP水平下对水的一个O-H键做柔性扫描

步骤一、先在平衡键长附近生成初始轨道：

哇！实在是太感谢您啦！这份教程真是详细到不能再详细了！❥(^_^)❥