计算化学公社

标题: 分享一个用过渡态理论TST和Wigner tunneling 校正算反应速率的代码 [打印本页]
作者
Author: Shannon    时间: 2015-2-2 06:04
标题: 分享一个用过渡态理论TST和Wigner tunneling 校正算反应速率的代码
分享一个从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;

  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;

  16. K=(kb*T/h)*exp(-deltaG2Joule/(R*T));

  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;

  22. wigner_K=K*tunnelingFactor;
  23. end
复制代码


作者
Author: 风飞    时间: 2018-8-5 15:47
您好！您是说用这个代码在matlab中计算速率常数吗？
作者
Author: 嘤嘤嘤    时间: 2022-7-30 23:53
”Kisthelp计算wigner校正时 频率单位转换似乎写错了，已经向作者反应，正在等回信中“请问有回信嘛？目前想要用该软件计算Wigner correction。




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