Shermo：计算气相分子配分函数和热力学数据的简单程序

sobereva · 发表于 Post on 2015-12-26 08:09:34

2020-May-19注：本文对应的是Shermo 1.0版，此程序已经完全没有意义了，因为后来笔者发布了Shermo 2.0版，介绍见《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）。2.0版比1.0版好用、强大、灵活得多得多得多，是日常量子化学研究者计算分子热力学数据离不开的工具。

Shermo：计算气相分子配分函数和热力学数据的简单程序

文/Sobereva @北京科音 2015-Dec-26

用过Gaussian的人都知道Freq任务会输出一堆热力学数据，但真正搞懂这些量到底怎么算出来的人不多。于是笔者开发了一个既有实用意义也有教学意义的Shermo程序（名字由来是Sob+Thermo）。

基于给定的谐振频率、惯性矩、温度、压力、原子质量、转动对称数等信息，Shermo程序可以输出分子配分函数和理想气体近似下的每mol的内能、焓、熵、自由能、热容，并且平动、转动、振动和电子贡献会独立输出，每种振动模式的贡献也能独立输出（这一点很有意义，能很方便直观地考察各个振动模式对热力学数据的影响）。

Shermo由Fortran编写，代码简洁易懂，也很适合学习热力学计算之用，以加深对概念的理解。程序附带了写得极为清楚的文档，里面有所有数据的计算公式，和源代码一对照就能透彻搞懂这些量是怎么算的了。

下载链接：http://sobereva.com/soft/Shermo.rar
压缩包里包含Windows版可执行程序、源代码、文档、示例输入文件以及与之对应的Gaussian freq任务的输出文件。

如果你的研究中使用了本程序，请这样引用：Tian Lu, Shermo program, http://sobereva.com/315 (accessed month day, year)

以下是输出例子：

Shermo: A utility to calculate various thermodynamic properties
Programmed by Sobereva (sobereva@sina.com)
First release: 2015-Dec-26

Temperature(K):    298.150
Pressure(Atm):    1.000
Rotational symmetry number: 6
Note: This is a non-linear molecule
Moments of inertia: 22.457510 90.665730 90.665730 amu*Bohr^2
Rotational constant: 80.362480 19.905440 19.905440 GHz
Rotational temperature: 3.856786 0.955309 0.955309 K
Spin multiplicity: 1
The number of atoms: 8

The number of frequencies: 18
Atom: 1 Mass:  12.000 amu
Atom: 2 Mass: 1.008 amu
Atom: 3 Mass: 1.008 amu
Atom: 4 Mass: 1.008 amu
Atom: 5 Mass:  12.000 amu
Atom: 6 Mass: 1.008 amu
Atom: 7 Mass: 1.008 amu
Atom: 8 Mass: 1.008 amu
Total mass: 30.046980 amu

Note: Only for translation motion, contribution to CV and U are different to CP
and H, respectively

                        ======= Translation =======

Translational q(T):       0.389858E+31
Translational q(T)/N:    0.647375E+07
Translational U(T): 0.888728 kcal/mol
Translational H(T): 1.481213 kcal/mol
Translational CV:    2.980807 cal/mol/K
Translational CP:    4.968012 cal/mol/K
Translational S(T): 36.133874 cal/mol/K

                        ========= Rotation ========

Rotational q(T):    0.810623E+03
Rotational U(T): 0.888728 kcal/mol
Rotational CV:    2.980807 cal/mol/K
Rotational S(T): 16.290714 cal/mol/K

                        ======== Vibration ========

  Mode  Wavenumber Freq    Vib. Temp. q(V=0)    q(BOT)
      cm^-1    GHz       K
1 312.37 0.9365E+04 449.43 1.284493 0.604507
2 827.24 0.2480E+05 1190.22 1.018810 0.138433
3 827.24 0.2480E+05 1190.22 1.018810 0.138432
4 1005.19 0.3013E+05 1446.24 1.007885 0.089144
5 1225.75 0.3675E+05 1763.58 1.002706 0.052087
6 1225.76 0.3675E+05 1763.59 1.002706 0.052087
7 1417.66 0.4250E+05 2039.69 1.001070 0.032728
8 1439.90 0.4317E+05 2071.69 1.000961 0.031015
9 1516.70 0.4547E+05 2182.19 1.000663 0.025761
10 1516.70 0.4547E+05 2182.20 1.000663 0.025761
11 1521.35 0.4561E+05 2188.88 1.000648 0.025473
12 1521.35 0.4561E+05 2188.89 1.000648 0.025473
13 3043.61 0.9125E+05 4379.08 1.000000 0.000647
14 3044.76 0.9128E+05 4380.72 1.000000 0.000645
15 3099.30 0.9291E+05 4459.20 1.000000 0.000565
16 3099.30 0.9291E+05 4459.20 1.000000 0.000565
17 3123.19 0.9363E+05 4493.57 1.000000 0.000534
18 3123.19 0.9363E+05 4493.57 1.000000 0.000534

  Mode  Wavenumber    ZPE    U(T)-U(0) U(T)    CV(T)    S(T)
      cm^-1    kcal/mol kcal/mol kcal/mol  cal/mol/K  cal/mol/K
1 312.373    0.447    0.254    0.701    1.650    1.350
2 827.245    1.183    0.044    1.227    0.607    0.186
3 827.245    1.183    0.044    1.227    0.607    0.186
4 1005.188    1.437    0.023    1.460    0.372    0.092
5 1225.753    1.752    0.009    1.762    0.189    0.037
6 1225.756    1.752    0.009    1.762    0.189    0.037
7 1417.658    2.027    0.004    2.031    0.100    0.017
8 1439.897    2.058    0.004    2.062    0.092    0.015
9 1516.703    2.168    0.003    2.171    0.071    0.011
10 1516.703    2.168    0.003    2.171    0.071    0.011
11 1521.351    2.175    0.003    2.178    0.070    0.011
12 1521.353    2.175    0.003    2.178    0.069    0.011
13 3043.615    4.351    0.000    4.351    0.000    0.000
14 3044.757    4.353    0.000    4.353    0.000    0.000
15 3099.303    4.431    0.000    4.431    0.000    0.000
16 3099.303    4.431    0.000    4.431    0.000    0.000
17 3123.188    4.465    0.000    4.465    0.000    0.000
18 3123.190    4.465    0.000    4.465    0.000    0.000

Vibrational q(V=0):       0.135737E+01
Vibrational q(BOT):       0.464776E-34
Vibrational ZPE: 0.074930 a.u.    47.019 kcal/mol    196.729 kJ/mol
Vibrational U(T)-U(0): 0.404393 kcal/mol
Vibrational U(T):       47.423838 kcal/mol
Vibrational CV(T):       4.085798 cal/mol/K
Vibrational S(T):       1.963524 cal/mol/K

                     ======== Electron spin ========

Note: Thermal excitation of electronic states is not taken into account, so ele
ctronic contribution to CV and U are zero
Electronic q:    1.000000
Electronic S:    0.000000 cal/mol/K

                        ===========================
                        ========== Total ==========
                        ===========================

Total q(V=0):          0.428966E+34
Total q(BOT):          0.146882E+00
Total q(V=0)/N:       0.712315E+10
Total q(BOT)/N:       0.243904E-24
Total CV(T):    10.047412 cal/mol/K
Total CP(T):    12.034617 cal/mol/K
Total S(T):    54.388112 cal/mol/K
Thermal correction to U(T):    49.201 kcal/mol 0.078407 a.u.
Thermal correction to H(T):    49.794 kcal/mol 0.079351 a.u.
Thermal correction to G(T):    33.578 kcal/mol 0.053510 a.u.

sobereva · 发表于 Post on 2018-3-29 16:08:15

清微发表于 2018-3-29 12:55
1.那么其实Gaussian用到的还是以势能面最低点作为零点求得的Ev(BOT),只是拆分成了ZPVE和Ev(V=0)两项，事 ...

1 姑且可以这么理解
2 仅在谐振近似下是这样

清微 · 发表于 Post on 2018-3-29 12:55:07

sobereva 发表于 2018-3-28 23:00
ZPE本身就算是振动产生的贡献的一部分
基于Q(V=0)计算升温导致的U的增加，将之与ZPE独立考虑，更便于 ...

1.那么其实Gaussian用到的还是以势能面最低点作为零点求得的Ev(BOT),只是拆分成了ZPVE和Ev(V=0)两项，事实上就只应该采用BOT作为零点，V=0为零点只是个辅助说明。
2.求解核的振动能级才能得到ZPVE，它的值是否严格等于Σi0.5hvi？

sobereva · 发表于 Post on 2018-3-28 23:00:36

清微发表于 2018-3-28 19:49
您在Shermo程序里，求振动对内能的贡献：第二个式子是以振动基态作为零点求得的Ev(V=0),而第一式加上第二 ...

ZPE本身就算是振动产生的贡献的一部分
基于Q(V=0)计算升温导致的U的增加，将之与ZPE独立考虑，更便于分别讨论两部分贡献。而且Q(V=0)的数量级显著大于Q(BOT)，后者太小

清微 · 发表于 Post on 2018-3-28 19:57:40

sobereva 发表于 2018-3-28 19:28
求解核的振动能级，即便是振动基态时能量也不为零，这个值就是ZPE

统计热力学公式里：直接用Ev(BOT)可以么？

清微 · 发表于 Post on 2018-3-28 19:49:31

sobereva 发表于 2018-3-28 19:28
求解核的振动能级，即便是振动基态时能量也不为零，这个值就是ZPE

您在Shermo程序里，求振动对内能的贡献：第二个式子是以振动基态作为零点求得的Ev(V=0),而第一式加上第二式相当于以势能面最低点作为零点求得的Ev(BOT),量化程序一般选用后者作为零点，得到Ev(BOT),然后从里面拿掉0.5hv这一项叫做零点能，而把其余部分叫做振动的贡献，这样做有什么必要？为什么不直接把整个Ev(BOT)再加上平动和转动叫做对内能的热力学校正量?

sobereva · 发表于 Post on 2018-3-28 19:28:28

清微发表于 2018-3-28 18:50
请问老师，零点振动能的物理意义是什么？
我们知道，在统计热力学中，计算振动配分函数时，有两种零点参考 ...

求解核的振动能级，即便是振动基态时能量也不为零，这个值就是ZPE

清微 · 发表于 Post on 2018-3-28 18:50:46

请问老师，零点振动能的物理意义是什么？
我们知道，在统计热力学中，计算振动配分函数时，有两种零点参考方法，一种以势能面最低点，另一种以振动基态，基于两者配分函数算出的内能就差了ZPE这一项。所以可以这么理解么：零点振动能仅仅是因为零点取得不同导致的，如果以振动基态作为零点，就不存在零点振动能了？

sobereva · 发表于 Post on 2018-1-27 00:50:36

清微发表于 2018-1-27 00:02
请问老师，我在windows 32位下测试，为什么出错了？错误信息为：

Shermo.ini不在当前目录

清微 · 发表于 Post on 2018-1-27 00:02:50

请问老师，我在windows 32位下测试，为什么出错了？错误信息为：

sobereva · 发表于 Post on 2017-12-5 23:45:35

让你变成回忆发表于 2017-12-5 20:21
Sob老师您好，跟着公式把这个程序写了一遍，仿照着您的程序在写。

程序的第53行，您写的是“rotcst3=h/( ...

是inert3，后来修改过这个bug，貌似当时忘了重新上传了。刚才已重新上传。

让你变成回忆 · 发表于 Post on 2017-12-5 20:21:01

Sob老师您好，跟着公式把这个程序写了一遍，仿照着您的程序在写。

程序的第53行，您写的是“rotcst3=h/(8D0*pi**2*inert2*amu2kg*(b2a*1D-10)**2)”

inert2是否是应该是inert3. 是您笔误还是就是这样？应该是笔误了吧？

我用其他的freq的任务，只有改了以后才能和Gaussian的结果对得上。

sobereva · 发表于 Post on 2017-8-14 08:16:50

beefly 发表于 2017-8-13 23:44
在程序里，电子对熵的贡献只取决于自旋多重度，这种近似是不是有问题？虽然Gaussian也是这么做的，但我估计 ...

是有这个问题，这只能靠有基础知识的用户自己把空间简并度等效归入到自旋多重度里体现了

beefly · 发表于 Post on 2017-8-13 23:44:52

本帖最后由 beefly 于 2017-8-13 23:45 编辑

在程序里，电子对熵的贡献只取决于自旋多重度，这种近似是不是有问题？虽然Gaussian也是这么做的，但我估计原因是Gaussian无法处理简并态。

按照Gaussian的说法，the first and higher excited states are assumed to be inaccessible at any temperature，一般情况下这是正确的。但是假如我们把激发态和基态无限拉近，也就是基态简并的情况，光考虑自旋多重度还不够，轨道角动量的简并度也要考虑。例如自由羟基OH，基态2Pi，总简并度是2*2=4，如果只取自旋多重度，得到的吉布斯自由能是不一样的。

flying1980 · 发表于 Post on 2016-11-30 14:19:06

非常感谢sob！

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[综合交流] Shermo：计算气相分子配分函数和热力学数据的简单程序

评分 Rate

浏览过的版块