计算化学公社

标题: Shermo：计算气相分子配分函数和热力学数据的简单程序 [打印本页]

作者
Author: sobereva 时间: 2015-12-26 08:09
标题: Shermo：计算气相分子配分函数和热力学数据的简单程序
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.

作者
Author: psfan 时间: 2015-12-26 09:51
源代码可在linux编译运行么？

作者
Author: ORCA_in_TCC 时间: 2015-12-26 11:10
sobereva对“Science requires sacrifice”做了很好的诠释！

作者
Author: hlmkh 时间: 2015-12-26 12:19
Hi sobereva
tanks for you
please intro English ver.

作者
Author: 我本是个娃娃 时间: 2015-12-26 13:15
Gaussian拼错了

作者
Author: sobereva 时间: 2015-12-26 13:29

我本是个娃娃发表于 2015-12-26 13:15
Gaussian拼错了

已改

作者
Author: sobereva 时间: 2015-12-26 13:30

psfan 发表于 2015-12-26 09:51
源代码可在linux编译运行么？

可

作者
Author: 978142355 时间: 2015-12-26 16:54
我测试了一下CO运行和得出的结果是没问题的。
只是对sob老师intro.pdf中的“（5）原子质量”中的“比如 2、6 号原子是氘，那么紧跟着
在后面写上它们的质量”这句话没弄明白，它是如何填写在*.ini中的，我自己按照sob老师所说的，自己理解着填写的（已经贴出图了），哪一种是对的或者两种不对用什么形式是对的？（PS：不是我懒，我是不知道如何将氢改为氘………………顺便在这里请教老师如果真想将乙烷中的氢改为氘，如何书写输入文件？）

作者
Author: sobereva 时间: 2015-12-26 18:50

978142355 发表于 2015-12-26 16:54
我测试了一下CO运行和得出的结果是没问题的。
只是对sob老师intro.pdf中的“（5）原子质量”中的“比如 2 ...

先设定所有元素的默认值，然后再对单独的原子设定来覆盖前面元素的默认设定。
你写的两个都不对。
应该把第一个当中后两个H都去掉，只写原子序号就行了。

作者
Author: 978142355 时间: 2015-12-26 18:55

sobereva 发表于 2015-12-26 18:50
先设定所有元素的默认值，然后再对单独的原子设定来覆盖前面元素的默认设定。
你写的两个都不对。
应 ...

*.ini文件中原子序号和C与H对齐对吧？
氘的高斯输入文件当中是将原子坐标那一栏H改成H（2）是这个意思吗？

作者
Author: sobereva 时间: 2015-12-26 19:02

978142355 发表于 2015-12-26 18:55
*.ini文件中原子序号和C与H对齐对吧？
氘的高斯输入文件当中是将原子坐标那一栏H改成H（2）是这个意思吗 ...

C 12.0       //默认对所有C使用12.0
H 1.00783    //默认对所有H使用1.00783
2 2.01410    //把H2原子改成2.01410
6 2.01410    //把H6原子改成2.01410

从输出信息中你可以直接看到各个原子都被指认成了什么质量，弄没弄错立马知道。

作者
Author: 978142355 时间: 2015-12-26 19:12

sobereva 发表于 2015-12-26 19:02
C 12.0       //默认对所有C使用12.0
H 1.00783    //默认对所有H使用1.00783
2 2.01410    //把H2 ...

roger sob，我马上再做一个相关测试，谢谢sob老师

作者
Author: 978142355 时间: 2015-12-26 20:28
本帖最后由 978142355 于 2015-12-27 12:27 编辑

首先在这里感谢sob老师编写的程序以及在遇到问题时给予的热心解惑，我将使用sob老师程序的心得（我个人出错和不解的地方）写出来：（1）CO分子Shermo.ini输入部分（图1），对于一个线性分子，主要强调的东西都在下图中显示，其中rotational symmetry number给出的是高斯输出文件中对应数值的相反数（在高斯的输出文件中查询是1，则此处输入为-1）。
（2）CH2DT分子Shermo.ini输入部分（图2，其中2号原子为D，5号原子为T），即甲烷中两个H分别被氘和氚取代。强调的是后面填写原子的相对原子质量的格式，先写出相关原子对应的相对原子质量，再对特殊的如氘、氚等进行填写。
另外：
a.CH2DT输入文件也给出了（图3），对于大神们来说，这很简单了，但是我是今天刚刚知道，所以为了后面的人学习也一并给出。
b.对于初学者，尤其是未接触fortran或刚刚接触fortran的同学来说，可能不太清楚这些文件放在哪里，图4给出了它们应该放置的地方（VS+IVF编译器），此处大神们可忽略。
最后再次感谢在其中遇到困难的时候，sob给予的耐心解答。本人还是感觉能将东西整理出来才算是真正学会了。sob老师的pdf介绍的已经很好了，建议大家仔细阅读，自己动手操作一遍。在下只不过是将自己不太清楚的地方整理出来，还望sob大神别心生误会。

作者
Author: smutao 时间: 2015-12-27 01:14
学习了！谢谢sob！
高斯官网这个地方是shermo的理论背景
Thermochemistry in Gaussian
http://www.gaussian.com/g_whitepap/thermo.htm
要获得高斯freq输出的这些热力学信息，实际上需要以下的这些信息：
1. 笛卡尔坐标
2. 能量的二阶导矩阵Hessian（这也是为什么只有在freq计算中我们才能得到热力学数据）
3. 原子质量
4. 温度
5. 压强
有了这5个输入，其余的结果按照公式计算就可以得到

作者
Author: smutao 时间: 2015-12-27 02:34
附上2015年初做的一个ppt
内容就是热化学信息在量化中是如何被计算的

作者
Author: sobereva 时间: 2015-12-27 09:59

978142355 发表于 2015-12-26 20:28
首先在这里感谢sob老师编写的程序以及在遇到问题时给予的热心解惑，我将使用sob老师程序的心得（我个人出错 ...

有些地方你说得不准确。
1 线型分子转动对称数不是随便给个负值就行，要不然转动贡献肯定不对，如文档里公式所示转动配分函数是依赖于转动对称数的。线型分子转动对称数是多少就得输入多少，比如H2是2，HCN是1，你额外再在前面加个符号让程序知道这个体系是线型的，以使用对应的公式。

2 关于却怎么确定同位素的，看来这里还要更确切地说一下，比如
C 12.0
H 1.00783
2 2.01410
6 2.01410
程序首先看到C指定为12，程序就扫描所有原子，发现元素名是C，就给它设定为12
然后程序看到H指定为1.00783，程序就扫描所有原子，发现元素名是H，就给它设定为1.00783
之后程序看到原子2号指定为2.01410，就把2号原子质量设为2.01410覆盖前面设的默认值。然后6号原子也是这样。
所以建议先把所有元素进行指定，然后再对个别的需要更改的进行指定，就会覆盖前面元素设定。

3 对于线型分子，不要输入三个转动常数，而应当输入一个值（高斯输出文件里第一个是0，后两个一样，就把后两个之一添到ini里）。对于线型分子Shermo只会读取一个值（从文档里的转动配分函数公式上也可以看到只依赖于一个惯性矩值），所以你写成0.00000 31.66694 31.66694，相当于程序读取的惯性矩是0，这导致结果完全不对，转动配分函数成了0。

PS：所以通过使用Shermo，对于加深概念的理解很有益，能弄清楚细节。

我已经把程序包在1L进行了更新，对于你提到的问题在文档里描述得更详细了一点，并且附加了线型分子和自由基的示例文件。

作者
Author: 978142355 时间: 2015-12-27 11:04
本帖最后由 978142355 于 2015-12-27 12:24 编辑

sobereva 发表于 2015-12-27 09:59
有些地方你说得不准确。
1 线型分子转动对称数不是随便给个负值就行，要不然转动贡献肯定不对，如文档 ...

多谢sob老师提出批评与指正，相关错误我自己也更正了。
原来线性分子那个地方不是随意给的一个负数啊..........我改-3和-1得出的结果是一致的，认为就是随便给个负数就行，看来我是从根上就出错了。
相对原子质量那里我也理解错了，很sorry。
线性分子惯性矩那里是我是完全理解错了，也很sorry，很感谢sob老师及时改正，以免其他人看完出错。

作者
Author: sobereva 时间: 2015-12-27 13:18

978142355 发表于 2015-12-27 11:04
多谢sob老师提出批评与指正，相关错误我自己也更正了。
原来线性分子那个地方不是随意给的一个负数啊... ...

不是一致的啊，比如HCN的例子，用-1和-3得到的转动部分的数值明显不同

作者
Author: 虎王 时间: 2015-12-27 16:14
用了这个程序如何引用？

作者
Author: sobereva 时间: 2015-12-27 16:28

虎王发表于 2015-12-27 16:14
用了这个程序如何引用？

这样吧
Tian Lu, Shermo program, http://http://sobereva.com/315 (access 日期)

作者
Author: 978142355 时间: 2015-12-27 22:07

sobereva 发表于 2015-12-27 13:18
不是一致的啊，比如HCN的例子，用-1和-3得到的转动部分的数值明显不同

老师可能误会了，我之前之所以认为Rotational symmetry number写-1和-3（甚至是任何一个负数都行）是因为Moments of inertia地方写成了3个数值（第一个数值是0，第二和第三数值相等），我测试了一下，对于CO，如果Moments of inertia写成3个数值，那么Rotational symmetry number写-1和-3都相同。
最后，还是感谢老师热心而耐心的指导，使我对这个程序的功能有彻底的了解。

作者
Author: wqmartin 时间: 2016-1-5 08:34
要是用Linux的脚本程序编辑就更好了，不过还是谢谢Sob!

作者
Author: 北极天99 时间: 2016-3-4 12:25
我解压了之后是这样子shermo.f90不能正确显示运行exe文件就只是闪一下就没了我应该怎么改啊

作者
Author: sobereva 时间: 2016-3-4 12:50

北极天99 发表于 2016-3-4 12:25
我解压了之后是这样子shermo.f90不能正确显示运行exe文件就只是闪一下就没了我应该怎么改啊

.f90是源代码，显然你得用文本编辑器打开才能查看
shermo运行后闪了一下就关闭显然是你的输入文件有问题，仔细看使用文档，里面写得很清楚需要把输入文件命名为Shermo.ini。

作者
Author: 北极天99 时间: 2016-3-4 15:00

sobereva 发表于 2016-3-4 12:50
.f90是源代码，显然你得用文本编辑器打开才能查看
shermo运行后闪了一下就关闭显然是你的输入文件有问 ...

老师，我运行的是给出的例子也运行不了

作者
Author: sobereva 时间: 2016-3-4 21:41

北极天99 发表于 2016-3-4 15:00
老师，我运行的是给出的例子也运行不了

绝对能运行
给了三个例子，把你要运行的例子改名成Shermo.ini就完了，手册里很明确地说程序启动时要读取当前目录下的Shermo.ini。

作者
Author: flying1980 时间: 2016-11-30 14:19
非常感谢sob！

作者
Author: beefly 时间: 2017-8-13 23:44
本帖最后由 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，如果只取自旋多重度，得到的吉布斯自由能是不一样的。

作者
Author: sobereva 时间: 2017-8-14 08:16

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

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

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

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

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

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

作者
Author: sobereva 时间: 2017-12-5 23:45

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

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

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

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

作者
Author: sobereva 时间: 2018-1-27 00:50

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

Shermo.ini不在当前目录

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

作者
Author: sobereva 时间: 2018-3-28 19:28

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

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

作者
Author: 清微 时间: 2018-3-28 19:49

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

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

作者
Author: 清微 时间: 2018-3-28 19:57

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

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

作者
Author: sobereva 时间: 2018-3-28 23:00

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

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

作者
Author: 清微 时间: 2018-3-29 12:55

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？

作者
Author: sobereva 时间: 2018-3-29 16:08

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

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

欢迎光临计算化学公社 (http://bbs.keinsci.com/)