求助高斯计算酞菁镍电催化CO2还原过程能量变化图

X教授

最近想研究用高斯计算酞菁镍电催化CO2还原过程（之前看到的文献都是用第一性原理计算的），看到一篇Nature Energy的文章（https://doi.org/10.1038/s41560-020-0667-9），是用高斯计算的，所以想先重复一下文献的工作。酞菁镍电催化CO2的还原过程主要有三步过程， 第一步是CO2在酞菁镍（*代替）表面转变成COOH（*+CO2+(H+)+(e-)=*COOH）, 第二步COOH转变成CO（*COOH+(H+)+(e-)=*CO+H2O）, 第三步CO脱附（*CO=*+CO），然后是计算方法，这是原文献使用的方法：Computational details. All DFT calculations on metal phthalocyanine molecules in this study were performed using the Gaussian 09 program. A PBE0 functional with D3 correction (Becke–Johnson damping) was adopted for its robustness and dispersion corrections, which make it widely accepted as the proper functional to study the reactions of transition metal complexes. The Stuttgart–Dresden pseudopotential and double-ξ valence basis set were used for transition metal atoms (nickel). For all other main group elements (H, C, N, O), the all-electron 6–31G* basis set was used. The geometric structures of all species were fully optimized. Harmonic vibrational frequencies were also computed, whose result showed that all reaction intermediates have no imaginary frequency. The gas-phase Gibbs free energies (G) of high- and low-spin forms of all intermediates, G, were calculated with the harmonic potential approximation at optimized structures at 298.15 K and 1 atm to determine the ground states in realistic conditions. These ground-state structures were then used in the mechanism study。
我的计算流程是：使用高斯16 A01优化酞菁镍（*），酞菁镍COOH（*COOH）以及酞菁镍CO（*CO）的结构，泛函使用pbe1pbe, C, H, O基组6–31G*，Ni使用SDD基组和赝势，同时使用empiricaldispersion=gd3bj，然后进行频率计算。结构优化过程中我发现如果不对酞菁镍骨架进行固定，酞菁镍COOH以及酞菁镍CO的结构会严重变形，因此我选择冻结酞菁镍骨架，最后得到优化后的结构以及Sum of electronic and thermal Free Energies：
E(*) = -1835.956069 Hartree
E(*COOH) = -2025.29839 Hartree
E(*CO) = -1949.5489 Hartree
E(CO2) = -188.39227 Hartree
E(H2O) = -76.32206 Hartree
E(H+) = -0.01 Hartree
E(CO) = -113.19593 Hartree
因此，起始物的能量E1 = E(*) + E(CO2) + 2E(H+) = -2025.47406 Hartree, 第一步能量E2 = E(*COOH) + E(H+) = -2025.28839 Hartree, 第三步能量E3 = E(*CO)  + E(H2O) = -2025.87096 Hartree, 最后一步能量E4 = E(*) + E(CO) + E(H2O) = -2025.47406 Hartree, 最后所有能量减去第一步能量再转换成eV得到能量变化图如下，结果非常奇怪，理论上E1应该是能量最低状态，而且第一步能量差应该很大（几十eV,太夸张）。
这是文献计算值，每一步能量差不到2eV, 而且E1与E4相同。所以想求教到底是什么原因导致差异如此之大的结果？

sobereva

贴图方式不对，其他人看不到。重新编辑帖子并仔细看置顶的新社员必读贴了解怎么正确贴图，此问题在这里还特意强调了：http://bbs.keinsci.com/thread-18961-1-1.html

X教授

sobereva 发表于 2020-10-21 02:49
贴图方式不对，其他人看不到。重新编辑帖子并仔细看置顶的新社员必读贴了解怎么正确贴图，此问题在这里还特 ...

多谢社长提醒，已更正

wzkchem5

卟啉、酞菁类的某些配合物配体严重变形是正常的，是二阶Jahn-Teller效应导致的。所以不要随意冻结骨架。另外注意各个物种的自旋多重度要写对。

X教授

好的，多谢，我试试不冻结骨架看看

radish-ding

我想问一下，那个电荷你设置的几

X教授

radish-ding 发表于 2020-11-14 18:58
我想问一下，那个电荷你设置的几

镍按+2价设置的话，酞菁镍就是0价，与COOH连接之后就是-1价，请问你也是做这个体系的么？

xiapin

你好！请问输入文件的自旋多重度你是怎么设置的？

jiangning198511

反应能量差十几个eV 这种结果明显就有问题，你需要调研一下其他类似的文献

大头攒毛

X教授 发表于 2020-11-17 07:05
镍按+2价设置的话，酞菁镍就是0价，与COOH连接之后就是-1价，请问你也是做这个体系的么？

*COOH是0价吧，因为是电催化，e-是电源给的，H+只能从溶液中来，一负一正，为0。

大头攒毛

请问一下，H+连电子都没有，它的能量可以用高斯计算吗？

sobereva

大头攒毛 发表于 2023-8-3 21:23
请问一下，H+连电子都没有，它的能量可以用高斯计算吗？

如果是计算真空下H+的电子能量，没有任何意义，直接就是0

熊琳

请问你这个问题解决了吗？我现在也是遇到这样的问题，直接算的每一步骤的自由能变，和相似文献中的能量曲线图相差极大

X教授

熊琳 发表于 2024-1-21 13:50
请问你这个问题解决了吗？我现在也是遇到这样的问题，直接算的每一步骤的自由能变，和相似文献中的能量曲线 ...

没解决，没继续搞这个体系了，如果严格按照热力学能来算能量呢？就是这个sob老师的这个帖子：http://sobereva.com/552，会有明显降低么？

熊琳

X教授 发表于 2024-1-21 23:49
没解决，没继续搞这个体系了，如果严格按照热力学能来算能量呢？就是这个sob老师的这个帖子：http://sobe ...

貌似需要用氢电极的概念对质子和电子的能量进行修正：e- + H+ → 1/2H2 ，即此时电子和质子的能量要用1/2H2的能量表示（见JPCB, 2004, 108,17886），修正之后反应自由能就正常了