计算化学公社

 找回密码 Forget password
 注册 Register
Views: 4576|回复 Reply: 4

[辅助/分析程序] 一个拟合双原子分子光谱常数的小工具

[复制链接 Copy URL]

66

帖子

4

威望

1469

eV
积分
1615

Level 5 (御坂)

大鸭鸭

发表于 Post on 2022-1-30 18:15:20 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 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工作区里面:
202201301759228734..png
然后通过调用diatomic函数拟合16O2的光谱常数:
  1. >> mol=diatomic(R,V,16,16);
复制代码
此时屏幕上会显示拟合的结果:
202201301803012147..png
并打印拟合得到的双原子分子常数和振-转能级:
  1. ==Diatomic Molecule Ro-vibration Calculation==

  2. ----Ro-vibration Constants---

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

  18. ----Ro-vibration Energy---

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

  20.         J        F(J)/cm-1
  21.    ------------------------------
  22.         0        0.0000
  23.         1        2.8501
  24.         2        8.5503
  25.         3        17.1003
  26.         4        28.4998
  27.         5        42.7482
  28.         6        59.8452
  29.         7        79.7899
  30.         8        102.5816
  31.         9        128.2194
  32.         10        156.7023

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

  34.         J        F(J)/cm-1
  35.    ------------------------------
  36.         0        0.0000
  37.         1        2.8192
  38.         2        8.4576
  39.         3        16.9149
  40.         4        28.1907
  41.         5        42.2846
  42.         6        59.1961
  43.         7        78.9245
  44.         8        101.4689
  45.         9        126.8286
  46.         10        155.0024

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

  48.         J        F(J)/cm-1
  49.    ------------------------------
  50.         0        0.0000
  51.         1        2.7883
  52.         2        8.3649
  53.         3        16.7294
  54.         4        27.8816
  55.         5        41.8210
  56.         6        58.5471
  57.         7        78.0591
  58.         8        100.3563
  59.         9        125.4378
  60.         10        153.3025

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

  62.         J        F(J)/cm-1
  63.    ------------------------------
  64.         0        0.0000
  65.         1        2.7574
  66.         2        8.2722
  67.         3        16.5440
  68.         4        27.5726
  69.         5        41.3574
  70.         6        57.8980
  71.         7        77.1937
  72.         8        99.2436
  73.         9        124.0469
  74.         10        151.6026

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

  76.         J        F(J)/cm-1
  77.    ------------------------------
  78.         0        0.0000
  79.         1        2.7265
  80.         2        8.1794
  81.         3        16.3585
  82.         4        27.2635
  83.         5        40.8938
  84.         6        57.2490
  85.         7        76.3283
  86.         8        98.1310
  87.         9        122.6561
  88.         10        149.9027

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

  90.         J        F(J)/cm-1
  91.    ------------------------------
  92.         0        0.0000
  93.         1        2.6956
  94.         2        8.0867
  95.         3        16.1731
  96.         4        26.9544
  97.         5        40.4302
  98.         6        56.5999
  99.         7        75.4629
  100.         8        97.0183
  101.         9        121.2653
  102.         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也欢迎提出。




diatomic.m

3.37 KB, 下载次数 Times of downloads: 21

MATLAB Function

评分 Rate

参与人数
Participants 9
威望 +1 eV +33 收起 理由
Reason
北大-陶豫 + 5 好物!
LEVELF + 1 好物!
ggdh + 5 不明觉厉
njfuzjs + 5 鸭鸭鸭鸭
Freeman + 2 好久不见
ene + 5 好!
卡开发发 + 5 хорошо!
hebrewsnabla + 5 好物!
sobereva + 1

查看全部评分 View all ratings

The leaping arc at your finger
is my faith and will eliminate never
Only the railgun will live forever

42

帖子

0

威望

1023

eV
积分
1065

Level 4 (黑子)

发表于 Post on 2022-3-15 11:50:48 | 显示全部楼层 Show all
冒昧的问一下,这个和LEVEL相比,两者的优缺点在什么地方呢?

66

帖子

4

威望

1469

eV
积分
1615

Level 5 (御坂)

大鸭鸭

 楼主 Author| 发表于 Post on 2022-3-15 20:37:55 | 显示全部楼层 Show all
qinjiu 发表于 2022-3-15 11:50
冒昧的问一下,这个和LEVEL相比,两者的优缺点在什么地方呢?

我没用过LEVEL,所以不知道,它只是我自己写的一个工具,谈不上有什么特别的优点
The leaping arc at your finger
is my faith and will eliminate never
Only the railgun will live forever

2

帖子

0

威望

53

eV
积分
55

Level 2 能力者

发表于 Post on 2022-3-29 22:33:15 | 显示全部楼层 Show all
qinjiu 发表于 2022-3-15 11:50
冒昧的问一下,这个和LEVEL相比,两者的优缺点在什么地方呢?

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

42

帖子

0

威望

1023

eV
积分
1065

Level 4 (黑子)

发表于 Post on 2022-4-3 20:09:48 | 显示全部楼层 Show all
serenity 发表于 2022-3-29 22:33
请问通过Gaussian得到的双原子分子势能曲线能否用LEVEL程序去计算它的振转能级?如果可以,这些(键长-能 ...

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

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 04:16 , Processed in 0.468595 second(s), 31 queries .

快速回复 返回顶部 返回列表 Return to list