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

linqiaosong

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

一、简介
    最近无聊写了一个用于拟合双原子分子光谱常数的小工具，其目的是想在没有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工作区里面：

然后通过调用diatomic函数拟合16O2的光谱常数：

1. >> mol=diatomic(R,V,16,16);

复制代码

此时屏幕上会显示拟合的结果：

并打印拟合得到的双原子分子常数和振-转能级：

1. ==Diatomic Molecule Ro-vibration Calculation==
   
2. 
   
3. ----Ro-vibration Constants---
   
4. 
   
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
    
20. 
    
21. ----Ro-vibration Energy---
    
22. 
    
23. v=0        G(v) = 789.1823 cm-1           Bv = 1.4251 cm-1
    
24. 
    
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
    
38. 
    
39. v=1        G(v) = 2349.3391 cm-1           Bv = 1.4096 cm-1
    
40. 
    
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
    
54. 
    
55. v=2        G(v) = 3885.8640 cm-1           Bv = 1.3942 cm-1
    
56. 
    
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
    
70. 
    
71. v=3        G(v) = 5398.7571 cm-1           Bv = 1.3787 cm-1
    
72. 
    
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
    
86. 
    
87. v=4        G(v) = 6888.0182 cm-1           Bv = 1.3633 cm-1
    
88. 
    
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
     
102. 
     
103. v=5        G(v) = 8353.6475 cm-1           Bv = 1.3478 cm-1
     
104. 
     
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也欢迎提出。

qinjiu

冒昧的问一下，这个和LEVEL相比，两者的优缺点在什么地方呢？

linqiaosong

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

我没用过LEVEL，所以不知道，它只是我自己写的一个工具，谈不上有什么特别的优点

serenity

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

请问通过Gaussian得到的双原子分子势能曲线能否用LEVEL程序去计算它的振转能级？如果可以，这些(键长-能量)数据点如何处理写入到LEVEL程序的输入文件，从而计算出正确的振转能级？望回复。

qinjiu

serenity 发表于 2022-3-29 22:33
请问通过Gaussian得到的双原子分子势能曲线能否用LEVEL程序去计算它的振转能级？如果可以，这些(键长-能 ...

可以，LEVEL支持折线点和解析形式的两种输入模式，具体看手册和例子

量化菜鸡

，你好，我尝试了一下，然后发现这个值过于小了，换算后，好像应该是angstrom单位吧

linqiaosong

量化菜鸡 发表于 2023-6-7 10:17
，你好，我尝试了一下，然后发现这个值过于小了，换算后，好像应该是angstrom单位吧

应该是单位写错了