求助：高压条件下的同位素动力学效应

dnlx

各位，我一直在模拟乙烷裂解这个地质反应。气体压力大约是500大气压，气体中95%是甲烷，5%是乙烷。温度是低温、长时间（200度，百万年尺度）。存在类似活性炭这样的催化剂（多孔材料，纳米级孔隙，存在自由基）。反应第一步是乙烷和自由基反应变成乙基（可逆），部分乙基转变成为不可逆的产物。

实际观察到的数据是，带有13C的乙烷反应得更快(inverse KIE)，12k/13k大概是0.98。如果我从reaction complex（乙烷+一个芳香基，不论是苯基还是芘基，代表活性炭的活性位）出发，DFT模拟可以得到Inverse KIE。但是如果从单独的反应物开始计算，总是normal KIE。我试验了其它催化途径，比如乙烷质子化，比如乙烷和黄铁矿（FeS2）反应，都是normal KIE。

现在想到的原因是1）因为这是个多孔材料，可能是传质控制；2）高压下的steric isotope effect；3) 高压下热力学上的非理想性。

关于1），传质过程是有利于轻同位素分子参加反应，而不是重同位素分子，所以想不出如何解释后者反应更快。

关于2），有一些文章（尤其是酶催化反应的文章）提到这个概念，但是没发现这样的量子化学模拟。它的根据是，重同位素分子及过渡态的相关键长更短，所以在体积局限的环境下（高压，或者多孔材料），反应更快。也有人做出了这样的体系（https://www.osti.gov/servlets/purl/983282）。
但是读社长的http://bbs.keinsci.com/thread-10159-1-1.html，一来极难模拟同位素对键长的影响，二来社长举例说：“对于BH（氢化硼），DBOC效应对平衡键长的影响才区区0.0007埃……对于更重的原子影响会更微乎其微”，所以很难想象高压和泛泛的多孔材料会带来明显的同位素效应。

关于3），这篇文章讲到同位素混合物的非理想性（https://www.degruyter.com/document/doi/10.1351/pac200476010011/pdf），这种效应应该很微弱，否则不会长期经常被忽略（见文章introduction）。但是我不知道它会不会带来百分之零点几的差别，以及能否进行理论模拟。

总之现在对模拟观察到的数据依然没有头绪，不知道是不是有哪里没考虑到或者没考虑对，请各位指教！

wth1219

dnlx 发表于 2022-7-3 21:22
好的，我现在从头细说一下：

1.        利用频率来计算EIE和KIE的公式，Urey 1947 和Bigeleisen 1949 就已经 ...

老师您好，Bigeleisen-Mayer公式中有一个容易被忽视的表征同位素替换前后对称性变化的项（因为我们通常研究的同位素交换体系并不涉及对称性的改变。），不知道您在计算中是否予以考虑了？

wth1219

wzkchem5 发表于 2022-7-4 02:17
我看了Bigeleisen的公式的推导，一开始文章把总配分函数写成了平动项、rigid rotor转动项、harmonic osci ...

老师您好，下图是一篇文章中关于Bigeleisen-Mayer公式的评述，使用B-M公式相比于直接求自由能之差似可避免数值误差。

wzkchem5

dnlx 发表于 2022-7-8 22:49
我没碰到这样的文献。KIE的准确程度不光在于自由能模型，文献提到误差的一个重要来源是数值计算过程。我 ...

那其实可以尝试一下提高格点，看看自由能作差的精度是否会有所改善。我还是觉得格点不会影响这么大，因为不同同位素用的是同一个Hessian，误差应该基本能抵消掉才对。

dnlx

wzkchem5 发表于 2022-7-9 05:21
不知道有没有文献比过B-M公式和QRRHO的自由能？这个截图提到B-M和自由能在理论上等价，但是事实上只有当 ...

我没碰到这样的文献。KIE的准确程度不光在于自由能模型，文献提到误差的一个重要来源是数值计算过程。我算过不同取代位置的次级同位素效应，用自由能方法（ORCA默认的，也就是QRRHO）得到的不同位置的次级KIE差别很明显，所以不很可靠。但是B-M公式得到的不同位置的次级KIE差别不大。估计这就是网格精度造成的。

wzkchem5

dnlx 发表于 2022-7-8 21:36
恰巧2021年一篇JACS文章做了评价（J. Am. Chem. Soc. 2021, 143, 21079−21099），见截屏。虽然它引 ...

不知道有没有文献比过B-M公式和QRRHO的自由能？这个截图提到B-M和自由能在理论上等价，但是事实上只有当用RRHO计算自由能时才等价，不知道QRRHO相比RRHO的改善是否比B-M避免使用小频率对精度的贡献更大。

dnlx

wzkchem5 发表于 2022-7-4 02:17
我看了Bigeleisen的公式的推导，一开始文章把总配分函数写成了平动项、rigid rotor转动项、harmonic osci ...

恰巧2021年一篇JACS文章做了评价（J. Am. Chem. Soc. 2021, 143, 21079−21099），见截屏。虽然它引用的还是上面提到的1999年文章。

终于找到了合适的参考文献。

wzkchem5

dnlx 发表于 2022-7-3 14:22
好的，我现在从头细说一下：

1.        利用频率来计算EIE和KIE的公式，Urey 1947 和Bigeleisen 1949 就已经 ...

我看了Bigeleisen的公式的推导，一开始文章把总配分函数写成了平动项、rigid rotor转动项、harmonic oscillator振动项的乘积，到此为止公式和RRHO自由能之差严格等价。接下来文章用Teller-Redlich rule把平动、转动贡献消掉了（这一步是严格的），然后再近似认为同位素取代对频率的影响足够小，经过一系列数学化简，得到的Bigeleisen的公式。
所以可以看出，Bigeleisen的公式引入了如下近似：（1）RRHO；（2）同位素取代对频率的影响是小量，比RRHO自由能之差还要多一个近似；如果你用的是没有用“同位素取代对频率的影响足够小”之前的那个公式，那结果和RRHO自由能之差是等价的。而QRRHO几乎总是比RRHO更准（参见QRRHO原始文献10.1002/chem.201200497），我至少目前是没见过RRHO定性正确但QRRHO定性错误的情形。这也就是我为什么一直觉得你这个结果奇怪，因为由上可知，除非因为误差抵消，否则理论上QRRHO没有道理比Bigeleisen的公式更差。

dnlx

wzkchem5 发表于 2022-7-2 03:28
你在用这个公式的时候，到底有没有把小于100cm-1的频率改成100cm-1？如果有的话，你的方法等价于Truhlar ...

好的，我现在从头细说一下：

1.        利用频率来计算EIE和KIE的公式，Urey 1947 和Bigeleisen 1949 就已经提出了（二人都是曼哈顿计划的参加者，因此才做了大量同位素化学反应的研究），公式很简明，原因是借用了Redlich-Teller 规则。后来文章多次测试Bigeleisen 1949公式的准确程度，见Shaad et al. 1999 (Can. J. Chem. 77: 875–878)：“As Saunders et al. (4) point out, it has long been known that the Bigeleisen eq. [4] is less sensitive to low frequency vibrations and hence is to be used if these are inaccurately known.”

2.        我学了ORCA之后，一直用这个教程的方法，用自由能的差值来计算KIE： https://www.orcasoftware.de/tuto ... calculating-the-kie
然后才知道低频近似的问题（我没再去对比Trular 和 Grimme的方法，因为已经无关自己的计算）。

3.        自由能方法和Bigeleisen公式的计算结果存在差别。比如温度423.15 K乙烷被芘基抓氢的反应，乙烷上的碳的Primary 和 secondary 13C KIE，用Bigeleisen公式是k12/k13 = 1.0099和1.0005，而用自由能方法得到的是1.0279和1.0211。上面网站的内容不知道是不是ORCA 开发者编写的。如果开发者有兴趣，可以对比二者为什么有区别，哪个更合理。我觉得Bigeleisen 1949公式是严格推导的，没有近似，而且secondary 13C KIE合理（自由能算出的二级KIE太大），我的所有结果会改用Bigeleisen 1949公式。

4.        但是主贴的反向KIE 问题仍然没有解决。所有公式给出的都是normal KIE。后来发现文献说，断键时候KIE就应该是正的，而成键的KIE才可能是反的，所以现在在找乙烷合适的η2反应。

5.        你说的乙烷协同反应变乙烯我也试了一下，和环己烯反应，自由能垒是70kcal/mol，太大，而且KIE也是正的。

我做到的就是这些。文献前面都给过了。

wzkchem5

dnlx 发表于 2022-7-1 16:12
根据得到的频率，用EXCEL算，公式推导见10.1002/9780470143476.ch2，实例见10.1016/j.gca.2008.12.022。
...

你在用这个公式的时候，到底有没有把小于100cm-1的频率改成100cm-1？如果有的话，你的方法等价于Truhlar的方法；如果没有的话，你的方法等价于RRHO。那么现在问题就来了，要么你的方法虽然本应该等价于Truhlar的方法，但是和程序用Truhlar的方法计算出来结果不一样；要么对于你这个体系RRHO反倒比quasi-RRHO好。不管是哪种情况，感觉都比较奇怪

dnlx

wzkchem5 发表于 2022-7-1 23:06
没明白你是什么意思，所以你的方法和Truhlar的方法是等价的，只是后者是程序自己算的，前者是你手算的？ ...

根据得到的频率，用EXCEL算，公式推导见10.1002/9780470143476.ch2，实例见10.1016/j.gca.2008.12.022。

wzkchem5

dnlx 发表于 2022-7-1 14:59
同样的基本公式（Grimme文章里的公式3）。问题在于低频处理，Truhlar的方法是把低频设成100，这个当然是 ...

没明白你是什么意思，所以你的方法和Truhlar的方法是等价的，只是后者是程序自己算的，前者是你手算的？那你上一个帖子为什么说“更不能用Truhlar方法”？

dnlx

wzkchem5 发表于 2022-7-1 21:43
不对吧，你是怎么用频率计算熵的，套的是什么公式？如果你手算的方法确实比Grimme、Truhlar的方法好，那 ...

同样的基本公式（Grimme文章里的公式3）。问题在于低频处理，Truhlar的方法是把低频设成100，这个当然是经验性的。Grimme的方法包括了经验拟合和近似，见他文章公式5之后的讨论。这些近似在计算反应熵的时候或许足够准确，但是计算异构体的熵变之差（Delta Delta S）的时候，误差很可能就被放大了。

实际情况就是前面说的，从Grimme的方法出发，反应物到反应复合物之间有很大的同位素分馏，这个肯定不对。

所以还不如从基本公式开始计算。如果认为低频不应该用公式3，那就全部改成100，也就回到Truhlar的方法了（同位素异构体的频率一一对应，Truhlar方法的弊端，如你说的http://bbs.keinsci.com/forum.php ... 92259&fromuid=32447，在此没有影响）。
【补充：用量子化学的结果去计算KIE的时候，采用的方法 可能还得回到Bigeleisen 1950年前后的工作。一个例子是10.1016/j.gca.2008.12.022，公式7。】

wzkchem5

dnlx 发表于 2022-7-1 14:26
问题解决了。计算KIE需要知道熵变，但是熵既不能用ORCA默认值（Grimme方法），更不能用Truhlar方法。二者都 ...

不对吧，你是怎么用频率计算熵的，套的是什么公式？如果你手算的方法确实比Grimme、Truhlar的方法好，那问题就来了，为什么他们没有把你的手算方法写到程序里？

dnlx

问题解决了。计算KIE需要知道熵变，但是熵既不能用ORCA默认值（Grimme方法），更不能用Truhlar方法。二者都存在近似，精度不够。应该自己用频率去计算振动熵，同位素异构体的频率都是一一对应的。这个本来是统计热力学的基础知识，但是之前一切交给软件，就把这些细节忽略了。