计算化学公社

标题: 一个拟合双原子分子光谱常数的小工具 [打印本页]
作者
Author: linqiaosong    时间: 2022-1-30 18:15
标题: 一个拟合双原子分子光谱常数的小工具
本帖最后由 linqiaosong 于 2022-1-30 18:25 编辑

一个拟合双原子分子光谱常数的小工具


一、简介
    最近无聊写了一个用于拟合双原子分子光谱常数的小工具，其目的是想在没有Molcas或者Molpro程序的时候也能够拟合双原子分子光谱常数。需要的数据也比较简单，只需要两个原子的相对原子质量和键长扫描的势能曲线结果就可以了。相对原子质量可以从元素周期表上直接查到，势能曲线一般可以通过Gaussian或者ORCA等量子化学程序得到。

二、使用方法
    将diatomic.m文件放在MATLAB当前工作文件夹里，通过"help diatomic"命令可以查看函数的使用方法：
  1. >> help diatomic
复制代码
   通过工具栏上的“导入数据”可以把量子化学程序得到的势能曲线数据以“列向量”的形式导入，注意导入的时候键长的列向量和能量的列向量一定要对应，不要错位了。建议扫描势能面的时候尽可能的扫描超过20个点，如果想要拟合的解离能比较准，那至少得扫描相当长的范围，键长最好不要低于7埃，否则拟合的解离能基本上没有参考价值。函数的输入如下所示：
  1. >> mol=diatomic(R,V,m1,m2,n);
复制代码
    函数至少需要4个参数，其中R是以"Angstrom"为单位的键长，格式为向量，V是以"Hartree"为单位的能量，格式也为向量，m1和m2是两个原子的质量数，n是可选的参数，是拟合势能曲线的多项式阶数。函数的输出变量mol是一个结构体，包含了双原子分子光谱常数和振-转能级。此外，函数还会在屏幕上显示势能曲线的拟合结果，并打印双原子分子光谱常数和振-转能级的结果。请仔细检查拟合的结果是否合理，尤其是在平衡键长附近是否拟合的好，以及末端的势能曲线是否平整，如果出现震荡，那么应该更改多项式拟合的阶数。

三、测试结果
    为了和光谱学数据库可以比较，我使用Molpro在MRCI/cc-pVTZ的理论水平下扫描了O2三重态基态势能曲线，然后Molpro可以顺带输出双原子分子光谱常数。同时我也从NIST数据库和HITRAN数据库上下载了16O2三重态基态的双原子分子光谱常数实验结果用于对比。拟合的势能曲线数据如下：
  1.     0.8   -149.3612797
  2.     0.9   -149.7913697
  3.     1.0   -149.9951086
  4.     1.1   -150.0794592
  5.     1.2   -150.1021935
  6.     1.3   -150.0944188
  7.     1.4   -150.0725797
  8.     1.5   -150.0453371
  9.     1.6   -150.0174181
  10.     1.7   -149.9915713
  11.     1.8   -149.9694731
  12.     1.9   -149.9521309
  13.     2.0   -149.9399549
  14.     2.1   -149.9325257
  15.     2.2   -149.9285677
  16.     2.3   -149.9266310
  17.     2.4   -149.9256998
  18.     2.5   -149.9252285
  19.     2.6   -149.9249602
  20.     2.7   -149.9247824
  21.     2.8   -149.9246477
  22.     2.9   -149.9245367
  23.     3.0   -149.9244411
  24.     3.1   -149.9243570
  25.     3.2   -149.9242826
  26.     3.3   -149.9242169
  27.     3.4   -149.9241596
  28.     3.5   -149.9241105
  29.     3.6   -149.9240691
  30.     3.7   -149.9240350
  31.     3.8   -149.9240075
  32.     3.9   -149.9239856
  33.     4.0   -149.9239686
  34.     4.1   -149.9239555
  35.     4.2   -149.9239455
  36.     4.3   -149.9239380
  37.     4.4   -149.9239325
  38.     4.5   -149.9239283
  39.     4.6   -149.9239252
  40.     4.7   -149.9239228
  41.     4.8   -149.9239211
  42.     4.9   -149.9239197
  43.     5.0   -149.9239186
  44.     5.1   -149.9239178
  45.     5.2   -149.9239171
  46.     5.3   -149.9239165
  47.     5.4   -149.9239160
  48.     5.5   -149.9239156
  49.     5.6   -149.9239153
  50.     5.7   -149.9239150
  51.     5.8   -149.9239147
  52.     5.9   -149.9239145
  53.     6.0   -149.9239144
  54.     6.1   -149.9239142
  55.     6.2   -149.9239141
  56.     6.3   -149.9239140
  57.     6.4   -149.9239139
  58.     6.5   -149.9239138
  59.     6.6   -149.9239138
  60.     6.7   -149.9239137
  61.     6.8   -149.9239137
  62.     6.9   -149.9239137
  63.     7.0   -149.9239137
  64.     7.1   -149.9239136
  65.     7.2   -149.9239136
  66.     7.3   -149.9239136
  67.     7.4   -149.9239136
  68.     7.5   -149.9239136
  69.     7.6   -149.9239137
  70.     7.7   -149.9239137
  71.     7.8   -149.9239137
  72.     7.9   -149.9239137
  73.     8.0   -149.9239137
  74.     8.1   -149.9239137
  75.     8.2   -149.9239137
  76.     8.3   -149.9239138
  77.     8.4   -149.9239138
  78.     8.5   -149.9239138
  79.     8.6   -149.9239138
  80.     8.7   -149.9239138
  81.     8.8   -149.9239138
  82.     8.9   -149.9239139
  83.     9.0   -149.9239139
  84.     9.1   -149.9239139
  85.     9.2   -149.9239139
  86.     9.3   -149.9239139
  87.     9.4   -149.9239140
  88.     9.5   -149.9239140
  89.     9.6   -149.9239140
  90.     9.7   -149.9239140
  91.     9.8   -149.9239140
  92.     9.9   -149.9239141
  93.    10.0   -149.9239141
复制代码
将其以列向量的形式导入到MATLAB工作区里面：
(, 下载次数 Times of downloads: 33)
然后通过调用diatomic函数拟合16O2的光谱常数：
  1. >> mol=diatomic(R,V,16,16);
复制代码
此时屏幕上会显示拟合的结果：
(, 下载次数 Times of downloads: 31)
并打印拟合得到的双原子分子常数和振-转能级：
  1. ==Diatomic Molecule Ro-vibration Calculation==

  3. ----Ro-vibration Constants---

  5. Reduced mass:
  6.          u = 8.0000 amu
  7. Moment of inertia:
  8.          I = 11.7654 amu
  9. Balanced bond length:
  10.          Re = 1.2127 cm-1
  11. Rotation constant:
  12.          Be = 1.4328 cm-1                De = 4.6906E-06 cm-1
  13. Vibration constant:
  14.          ve = 1583.7887 cm-1                vexe = 11.8159 cm-1
  15.          ae = 0.0155 cm-1                Y00 = 0.2419 cm-1
  16. Dissociation energy:
  17.          D = 112.4785 kcal/mol
  18. Zero-point energy:
  19.          ZPE = 2.2564 kcal/mol

  21. ----Ro-vibration Energy---

  23. v=0        G(v) = 789.1823 cm-1           Bv = 1.4251 cm-1

  25.         J        F(J)/cm-1
  26.    ------------------------------
  27.         0        0.0000
  28.         1        2.8501
  29.         2        8.5503
  30.         3        17.1003
  31.         4        28.4998
  32.         5        42.7482
  33.         6        59.8452
  34.         7        79.7899
  35.         8        102.5816
  36.         9        128.2194
  37.         10        156.7023

  39. v=1        G(v) = 2349.3391 cm-1           Bv = 1.4096 cm-1

  41.         J        F(J)/cm-1
  42.    ------------------------------
  43.         0        0.0000
  44.         1        2.8192
  45.         2        8.4576
  46.         3        16.9149
  47.         4        28.1907
  48.         5        42.2846
  49.         6        59.1961
  50.         7        78.9245
  51.         8        101.4689
  52.         9        126.8286
  53.         10        155.0024

  55. v=2        G(v) = 3885.8640 cm-1           Bv = 1.3942 cm-1

  57.         J        F(J)/cm-1
  58.    ------------------------------
  59.         0        0.0000
  60.         1        2.7883
  61.         2        8.3649
  62.         3        16.7294
  63.         4        27.8816
  64.         5        41.8210
  65.         6        58.5471
  66.         7        78.0591
  67.         8        100.3563
  68.         9        125.4378
  69.         10        153.3025

  71. v=3        G(v) = 5398.7571 cm-1           Bv = 1.3787 cm-1

  73.         J        F(J)/cm-1
  74.    ------------------------------
  75.         0        0.0000
  76.         1        2.7574
  77.         2        8.2722
  78.         3        16.5440
  79.         4        27.5726
  80.         5        41.3574
  81.         6        57.8980
  82.         7        77.1937
  83.         8        99.2436
  84.         9        124.0469
  85.         10        151.6026

  87. v=4        G(v) = 6888.0182 cm-1           Bv = 1.3633 cm-1

  89.         J        F(J)/cm-1
  90.    ------------------------------
  91.         0        0.0000
  92.         1        2.7265
  93.         2        8.1794
  94.         3        16.3585
  95.         4        27.2635
  96.         5        40.8938
  97.         6        57.2490
  98.         7        76.3283
  99.         8        98.1310
  100.         9        122.6561
  101.         10        149.9027

  103. v=5        G(v) = 8353.6475 cm-1           Bv = 1.3478 cm-1

  105.         J        F(J)/cm-1
  106.    ------------------------------
  107.         0        0.0000
  108.         1        2.6956
  109.         2        8.0867
  110.         3        16.1731
  111.         4        26.9544
  112.         5        40.4302
  113.         6        56.5999
  114.         7        75.4629
  115.         8        97.0183
  116.         9        121.2653
  117.         10        148.2029
复制代码
这个地方由于veye我没有公式，所以就没有考虑了，离心畸变的Dv也没有考虑振动耦合的效应，一律用的De，所以对振-转能级可能会有一些微小的影响，但是它们都是高阶校正量，所以一般影响都比较小。    我们对比一下NIST库上的数据，拟合得到的结果和Molpro输出的结果：
  1. ========================================================================================================
  2. TITLE        ENRE     RE(AU)   RE(A)   BE       AE      WE     WEXE   WEYE   DE    D0      SIG
  3. ========================================================================================================
  4. MOLPRO     -150.1024270  2.2943  1.2141   1.430  0.0158  1558.11  11.25  0.05  4.86  4.76  0.4500955E-01
  5. FIT                             1.2127   1.433  0.0155  1583.79  11.82        4.69
  6. NIST                            1.2075   1.438  0.0159  1580.19  11.98  0.04  4.83   
  7. ========================================================================================================
复制代码
此外我还对比HITRAN数据库上16O2的1-0振动跃迁的Q(5)支和S(3)支的谱线：
  1. ===========================================
  2.                         1-0 Q(5)                1-0 S(3)
  3. FIT                        1559.6930                1585.3410
  4. HITRAN                1555.9121                1581.7485
  5. ===========================================
复制代码
目前来看结果都还是比较好的。

代码可以在附件下载，可以随意使用和修改，没有版权也不需要引用，如果有什么建议和bug也欢迎提出。




作者
Author: qinjiu    时间: 2022-3-15 11:50
冒昧的问一下，这个和LEVEL相比，两者的优缺点在什么地方呢？
作者
Author: linqiaosong    时间: 2022-3-15 20:37
qinjiu 发表于 2022-3-15 11:50
冒昧的问一下，这个和LEVEL相比，两者的优缺点在什么地方呢？

我没用过LEVEL，所以不知道，它只是我自己写的一个工具，谈不上有什么特别的优点
作者
Author: serenity    时间: 2022-3-29 22:33
qinjiu 发表于 2022-3-15 11:50
冒昧的问一下，这个和LEVEL相比，两者的优缺点在什么地方呢？

请问通过Gaussian得到的双原子分子势能曲线能否用LEVEL程序去计算它的振转能级？如果可以，这些(键长-能量)数据点如何处理写入到LEVEL程序的输入文件，从而计算出正确的振转能级？望回复。
作者
Author: qinjiu    时间: 2022-4-3 20:09
serenity 发表于 2022-3-29 22:33
请问通过Gaussian得到的双原子分子势能曲线能否用LEVEL程序去计算它的振转能级？如果可以，这些(键长-能 ...

可以，LEVEL支持折线点和解析形式的两种输入模式，具体看手册和例子
作者
Author: 量化菜鸡    时间: 2023-6-7 10:17
(, 下载次数 Times of downloads: 23) ，你好，我尝试了一下，然后发现这个值过于小了，换算后，好像应该是angstrom单位吧
作者
Author: linqiaosong    时间: 2023-12-12 10:16
量化菜鸡 发表于 2023-6-7 10:17
，你好，我尝试了一下，然后发现这个值过于小了，换算后，好像应该是angstrom单位吧

应该是单位写错了



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