分享一个用过渡态理论TST和Wigner tunneling 校正算反应速率的代码

Shannon

分享一个从Kisthelp源代码中提取出来的 用TST+wigner校正算反应速率的代码。 代码照抄自Kisthelp的 源代码，固结果是一样的。输入 只需要提供反应自由能变化和过渡态 的 那个虚频 即可。因此，除了高斯外，Orca/Mopac 等软件的输出也能用来算反应速率了。温度是默认的298K, 其他的温度需要改动下代码，因此没有 写。
阿伦尼乌斯公式 K=A*exp(Ea/RT) 和Erying 公式的对应关系， 反应的焓ΔH@298K 即是对应的活化能Ea，反应的熵变化对速率的影响包含在 指数前系数 A里面了。


Sob 老师写的关于Kisthelp软件用法的帖子：http://sobereva.com/246


另外，Kisthelp计算wigner校正时 频率单位转换似乎写错了，已经向作者反应，正在等回信中。

原来的kisthelp 不能计算 分子对称性对 反应自由能 的贡献，这个提取出来的代码因为不需要向Kisthelp那样从新计算delta G ，而是直接使用高斯、orca，mopac，等等的输出 因此反倒比kisthelp更准确点。


代码如下，用MATLAB写的。 只是些简单的加减乘除平方，各位如果想将它改成fortran, python,c++或者 Excel 都挺简单。

1. function [wigner_K,tunnelingFactor,K] = Calc_TST_bimolecular (deltaG,negmode_from_gauss)
   
2. % below is the code for bi-molecular reaction,
   
3. % result is consistent with Kisthelp
   
4. % reprogrammed from Kisthelp source code
   
5. deltaG2Joule=deltaG*1000;
   
6. h=6.6260696E-34 ;
   
7. kb=1.380649E-23;
   
8. T=298;
   
9. P0 = 1E05;
   
10. NA = 6.0221413E23;
    
11. c = 2.99792458E8;
    
12. R = NA * kb;
    
13. unitCoeff =  ((R*T) /P0)*(1.0e6/NA);
    
14. convertCm_1ToM_1=1.0E2;
    
15. convertCm_1ToKelvin= 100*c*h/kb;
    
16. 
    
17. K=(kb*T/h)*exp(-deltaG2Joule/(R*T));
    
18. K=K*unitCoeff;
    
19. freqImag=-negmode_from_gauss*(convertCm_1ToM_1/convertCm_1ToKelvin);%这个好像原作者错了
    
20. factor =  (  h * freqImag * ( c /( kb * T) ))^2;
    
21. tunnelingFactor = 1.0 + (1.0/24.0) * factor;
    
22. wigner_K= K*tunnelingFactor;
    
23. end
    

复制代码

1. function [wigner_K,tunnelingFactor,K] = Calc_TST_unimolecular (deltaG,negmode_from_gauss)
   
2. % below is the code for uni-molecular reaction,
   
3. % result is consistent with Kisthelp
   
4. % reprogrammed from Kisthelp source code
   
5. deltaG2Joule=deltaG*1000;
   
6. h=6.6260696E-34 ;
   
7. kb=1.380649E-23;
   
8. P0 = 1E05;
   
9. NA = 6.0221413E23;
   
10. c = 2.99792458E8;
    
11. R = NA * kb;
    
12. T=298;
    
13. convertCm_1ToM_1=1.0E2;
    
14. convertCm_1ToKelvin= 100*c*h/kb;
    
15. 
    
16. K=(kb*T/h)*exp(-deltaG2Joule/(R*T));
    
17. 
    
18. freqImag=-negmode_from_gauss*(convertCm_1ToM_1/convertCm_1ToKelvin);%这个好像原作者错了
    
19. factor =  (  h * freqImag * ( c /( kb * T) ))^2;
    
20. tunnelingFactor = 1.0 + (1.0/24.0) * factor;
    
21. 
    
22. wigner_K=K*tunnelingFactor;
    
23. end

复制代码

风飞

您好！您是说用这个代码在matlab中计算速率常数吗？

嘤嘤嘤

”Kisthelp计算wigner校正时 频率单位转换似乎写错了，已经向作者反应，正在等回信中“请问有回信嘛？目前想要用该软件计算Wigner correction。