计算化学公社
标题: 求助：高压条件下的同位素动力学效应 [打印本页]
作者Author: dnlx 时间: 2022-6-25 20:40
标题: 求助：高压条件下的同位素动力学效应
本帖最后由 dnlx 于 2024-5-26 07:41 编辑
各位，我一直在模拟乙烷裂解这个地质反应。气体压力大约是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）。但是我不知道它会不会带来百分之零点几的差别，以及能否进行理论模拟。
总之现在对模拟观察到的数据依然没有头绪，不知道是不是有哪里没考虑到或者没考虑对，请各位指教！
作者Author: wzkchem5 时间: 2022-6-25 21:36
我觉得更可能的原因是机理猜错了，比如可能是乙烷通过协同机理往活性炭（或者FeS2等）上面转移两个氢变成乙烯。百万年时间尺度下，常常会有实验室条件下观察不到的新机理出现，尤其对于乙烷这样的实验室条件下比较惰性的分子，可能完全走的是实验室条件下极难看到的非传统途径也说不定。
作者Author: dnlx 时间: 2022-6-26 20:13
本帖最后由 dnlx 于 2022-6-26 20:37 编辑
虽然有可能存在漏掉的机理，但是可能性不大，自由基抽氢应该是最优势的反应，自由能垒很低。其速度只是局限于自由基的低浓度。
乙烷参与的协同反应的能垒远高于抽氢反应，合适活性位的浓度更低，而且不能保持活性位的浓度（这个和自由基机理不一样）。
我估计了一下高压对KIE的影响，如果下面的公式（Basukevsky et al. 1985）可以用，
(, 下载次数 Times of downloads: 2)
那么高压（恒压）下12k/13k = exp(p ΔΔV/RT)，ΔΔV是同位素替换引起的活化体积（volume of activation）的变化。即使用非常夸大的键长变化（同位素替换之后，过渡态C-H键长降低0.005埃）和夸大的分子动力学直径（1 nm），12k/13k 是0.996，和1差得不远，不能弥补S带来的KIE（12k/13k = 1.03）。
不知道上面的计算是否合适。另外，由于长链烷烃也观察到类似的现象，所以这个反应不会是在非常局限的小孔中发生的，所以孔隙的位阻可以不考虑。
现在不知道非理想性的影响该怎么估算一下。
作者Author: wzkchem5 时间: 2022-6-26 20:59
你计算的时候有没有考虑13C取代的乙烷的symmetry number的影响？单个13C取代的乙烷的symmetry number和普通乙烷是不一样的，有的程序可能不能正确处理
作者Author: dnlx 时间: 2022-6-26 21:04
考虑了。如果对称数不对，结果会错得离谱。
作者Author: wzkchem5 时间: 2022-6-26 21:09
那我感觉还是挺奇怪的，如果从reaction complex算是inverse KIE，从单独反应物算是normal KIE，那就意味着从单独反应物生成reaction complex这一步有很大的equilibrium isotope effect。但是这一步就是一个非共价吸附而已，不应该有很大的同位素效应才对。搞不好可能和低频振动的处理有关？
作者Author: dnlx 时间: 2022-6-26 21:28
本帖最后由 dnlx 于 2022-6-26 22:22 编辑
是很奇怪，我可以保证数据处理没错，因为计算了不同的同位素取代（两个位置的13C，以及各个位置的氘），互相可以印证。
这个反应复合物的单点能低于反应物的加和，但是自由能高于反应物加和（几个kcal/mol）。确实如你所说，EIE不应该这么大。如果是低频振动处理带来的问题，请问哪个误差更小，是从反应复合物开始算，还是按照单独反应物加和？
=====================
补充：有可能又回到这个问题了：
http://bbs.keinsci.com/forum.php ... 2&fromuid=32447
作者Author: wzkchem5 时间: 2022-6-27 02:54
该不该从反应复合物开始算和有没有低频没关系，应该总是取反应物和反应复合物里面自由能较低的那个。低频问题可以用Shermo算一下试试
作者Author: dnlx 时间: 2022-6-27 04:11
本帖最后由 dnlx 于 2022-6-27 05:38 编辑
Shermo（ilowfreq设为2，Grimme's entropy interpolation）跟ORCA的结果是一样的。比如过渡态，熵只差了0.004 cal/mol/K，最后一位有效数字的区别。
我现在的疑问在这里：如你前面指出的，“单独反应物加和”与“反应复合物”之间是个很弱的吸附，不应该有明显的同位素分馏。所以我得到的二者的熵，应该至少有一个不对。在排除数据处理出错的情况下，是不是因为模拟精度问题或者 计算振动熵的模型的问题。
附件是单独反应物、反应复合物、过渡态的输出文件。如果需要，我可以上传同位素替换结果（这个简单而罗嗦）。
【补充：明白你的意思了，是用Shermo里面Truhlar的方法去计算频率100以下的部分】
(, 下载次数 Times of downloads: 2)
作者Author: wzkchem5 时间: 2022-6-27 15:11
哦，你本来用的就是orca是吧，那就不是Shermo能解决的问题了，Grimme的QRRHO一般比Truhlar的方法好。
可以试一下构象搜索。有一个可能性是你观察到的EIE是因为你只把乙烷里的一个C换成了13C，如果考虑到两个C都可能是13C、概率基本相等，可能EIE就抵消了。
作者Author: dnlx 时间: 2022-7-1 21:26
问题解决了。计算KIE需要知道熵变，但是熵既不能用ORCA默认值（Grimme方法），更不能用Truhlar方法。二者都存在近似，精度不够。应该自己用频率去计算振动熵，同位素异构体的频率都是一一对应的。这个本来是统计热力学的基础知识，但是之前一切交给软件，就把这些细节忽略了。
作者Author: wzkchem5 时间: 2022-7-1 21:43
不对吧，你是怎么用频率计算熵的，套的是什么公式？如果你手算的方法确实比Grimme、Truhlar的方法好，那问题就来了，为什么他们没有把你的手算方法写到程序里？
作者Author: dnlx 时间: 2022-7-1 21:59
本帖最后由 dnlx 于 2022-7-1 23:08 编辑
同样的基本公式（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。】
作者Author: wzkchem5 时间: 2022-7-1 23:06
没明白你是什么意思，所以你的方法和Truhlar的方法是等价的，只是后者是程序自己算的，前者是你手算的？那你上一个帖子为什么说“更不能用Truhlar方法”？
作者Author: dnlx 时间: 2022-7-1 23:12
根据得到的频率，用EXCEL算，公式推导见10.1002/9780470143476.ch2，实例见10.1016/j.gca.2008.12.022。
作者Author: wzkchem5 时间: 2022-7-2 03:28
本帖最后由 wzkchem5 于 2022-7-1 21:30 编辑
你在用这个公式的时候，到底有没有把小于100cm-1的频率改成100cm-1？如果有的话，你的方法等价于Truhlar的方法；如果没有的话，你的方法等价于RRHO。那么现在问题就来了，要么你的方法虽然本应该等价于Truhlar的方法，但是和程序用Truhlar的方法计算出来结果不一样；要么对于你这个体系RRHO反倒比quasi-RRHO好。不管是哪种情况，感觉都比较奇怪
作者Author: dnlx 时间: 2022-7-3 21:22
好的，我现在从头细说一下：
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也是正的。
我做到的就是这些。文献前面都给过了。
作者Author: wzkchem5 时间: 2022-7-4 02:17
我看了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的公式更差。
作者Author: dnlx 时间: 2022-7-9 04:36
恰巧2021年一篇JACS文章做了评价（J. Am. Chem. Soc. 2021, 143, 21079−21099），见截屏。虽然它引用的还是上面提到的1999年文章。
终于找到了合适的参考文献。
作者Author: wzkchem5 时间: 2022-7-9 05:21
不知道有没有文献比过B-M公式和QRRHO的自由能？这个截图提到B-M和自由能在理论上等价，但是事实上只有当用RRHO计算自由能时才等价，不知道QRRHO相比RRHO的改善是否比B-M避免使用小频率对精度的贡献更大。
作者Author: dnlx 时间: 2022-7-9 05:49
我没碰到这样的文献。KIE的准确程度不光在于自由能模型，文献提到误差的一个重要来源是数值计算过程。我算过不同取代位置的次级同位素效应，用自由能方法（ORCA默认的，也就是QRRHO）得到的不同位置的次级KIE差别很明显，所以不很可靠。但是B-M公式得到的不同位置的次级KIE差别不大。估计这就是网格精度造成的。
作者Author: wzkchem5 时间: 2022-7-9 15:26
那其实可以尝试一下提高格点，看看自由能作差的精度是否会有所改善。我还是觉得格点不会影响这么大，因为不同同位素用的是同一个Hessian，误差应该基本能抵消掉才对。
作者Author: wth1219 时间: 2022-7-10 12:18
老师您好，下图是一篇文章中关于Bigeleisen-Mayer公式的评述，使用B-M公式相比于直接求自由能之差似可避免数值误差。
作者Author: wth1219 时间: 2022-7-10 12:42
老师您好，Bigeleisen-Mayer公式中有一个容易被忽视的表征同位素替换前后对称性变化的项（因为我们通常研究的同位素交换体系并不涉及对称性的改变。），不知道您在计算中是否予以考虑了？
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |