计算化学公社

标题: 关于用Gaussian做水分子的OH键的柔性扫描遇到的能量突然跳变的问题 [打印本页]
作者
Author: gingaorb    时间: 2023-4-14 12:16
标题: 关于用Gaussian做水分子的OH键的柔性扫描遇到的能量突然跳变的问题
各位小哥哥，各位大佬们，请教一个问题呀~
是这样的，我最近想探索水分子中的一个氢氧键在断裂的过程中有没有可能性均裂生成羟基自由基和氢原子，于是做了一个水分子的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”中提到的原因一致呢？万分感谢大家了！

作者
Author: sobereva    时间: 2023-4-14 15:04
这么扫描没意义
涉及共价键断裂情况的扫描怎么做看此文，对于DFT必须考虑对称破缺
详谈使用Gaussian做势能面扫描
http://sobereva.com/474（http://bbs.keinsci.com/thread-12660-1-1.html）
作者
Author: gingaorb    时间: 2023-4-14 22:14
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~
作者
Author: wjc404    时间: 2023-4-18 03:22
本帖最后由 wjc404 于 2023-4-18 03:46 编辑

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

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

至于为什么是得到氢和氢氧自由基而不是两个离子，可以去算一下真空中H原子的电离能和氢氧自由基的电子亲和能，就知道了。水中的氢离子和氢氧根离子都是高度水合的，被充分稳定了，和真空不同。
作者
Author: gingaorb    时间: 2023-4-18 21:00
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)，但第一步就几何优化不收敛。
作者
Author: Kmetsch    时间: 2023-4-18 21:14
本帖最后由 Kmetsch 于 2023-4-18 21:16 编辑
gingaorb 发表于 2023-4-18 21:00
能请问下“真空中H原子的电离能和氢氧自由基的电子亲和能”这个用Gaussian要怎么计算呢？另外我也尝试了 ...

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

  5. %casscf
  6.   nel 8
  7.   norb 6
  8.   ptmethod fic_nevpt2
  9.   ptsettings
  10.     d3tpre 0
  11.     d4tpre 0
  12.   end
  13. end

  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

  21. %numgrad dx 0.001 end

  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. *

复制代码

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
INF0.974 N.A.0.0
作者
Author: gingaorb    时间: 2023-4-28 22:57
wjc404 发表于 2023-4-28 21:27
用ORCA在NEVPT2/def2-TZVPP水平下对水的一个O-H键做柔性扫描

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


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




欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3