自写python脚本实现热力学校正量的计算

ByTheWay

    大家好，这应该是楼主在本论坛的第一篇帖子。楼主平时是给实验组做有机机理计算的，偶尔也跑一些固体相关的计算（例如表面催化机理等等）。楼主大概在几年前就开始关注本论坛，很多工作也受到了论坛里的帖子和Sob老师博文的帮助，可惜楼主太懒，一直没有注册帐号。最近楼主终于想起来注册了一个帐号，希望借此机会和论坛里各位老师同学多多交流。

    咱们言归正传，楼主在前段时间的研究中使用了Sob老师编写的Shermo软件包，深切感受到了它的强大与易用。比较巧合的是楼主最近稍微了解了一些有关热力学计算的知识，因此就想到能否编写一些小程序，实现一些较为简单的分子体系的热力学性质计算，从而加深对于所学知识的理解。楼主本来想仿照Sob老师的Shermo软件来编写自己的代码，可惜楼主实在不懂得Fortran编程，也因此不很看得懂Shermo的源代码（尤其是对称性判断部分），所以楼主最终打算使用python来完成这个任务。可惜，通过检索楼主发现ASE软件包已经提供了不少热力学性质计算的功能，例如使用其中的数据类IdealGasThermo就能够在谐振近似下计算分子体系的各种热力学数据了。但ASE并不能很好地处理低频振动对于热力学量的贡献，并且不能自动判断分子的转动对称数。因此楼主的目标就变成了在ASE的代码中加入qRRHO方法，以及对于转动对称数进行识别的功能。最终楼主写出了这套称为HeatMeter的脚本，接下来楼主会简单介绍编写的思路，希望对有兴趣的同学们起到一定的启发作用。

    我们首先来看看Grimme的qRRHO方法是如何实现的。根据参考文献[1]中的公式(7)，不难发现在Grimme的qRRHO方法中，体系的振动熵由原本的仅仅使用谐振近似（公式(3)）计算，变为了由谐振近似下的振动熵与自由转子近似下的转动熵（公式(6)）共同贡献。由于当振动模式的频率较低时，它对于体系振动熵的贡献较为接近自由转子近似所求得的数值，反之，它对于体系振动熵的贡献较为接近谐振近似所求得的数值，因此对于每个振动模式，这两种不同的来源的熵对于总振动熵的贡献由切换函数（公式(8)）确定。随着振动模式频率的升高，由公式(8)求得的权重就越大，而公式(7)中的谐振近似熵所占的比例也就越大。公式(3)至公式(8)的实现非常容易，不过应当注意单位转换的问题：

1. import numpy as np
   
2. from ase import units
   
3. 
   
4. def vibrational_entropy_contribution_qRRHO_Grimme(vib_energies, temperature, low_freq):
   
5.     kT = units.kB*temperature
   
6.     B_av = 1e-44
   
7.     alpha = 4.0
   
8.     S_t = 0.0
   
9.     for e in vib_energies:
   
10.         # calculate S_v
    
11.         x = e/kT
    
12.         S_v = x/(np.exp(x) - 1.0) - np.log(1.0 - np.exp(-x))
    
13.         S_v *= units.kB
    
14.         # calculate S_r
    
15.         mu = units._hbar*units._hbar/(2.0 * e / units.J)  # J*s^2
    
16.         mu_av = mu * B_av / (mu + B_av)  # J*s^2
    
17.         tmp = 2.0 * np.pi * mu_av * kT / units.J / (units._hbar * units._hbar)
    
18.         S_r = units.kB * (0.5 + np.log(np.sqrt(tmp)))
    
19.         # Mix S_r and S_v with the Head-Gordon damping function.
    
20.         w = 1.0 / (1.0 + np.power((low_freq / (e / units.invcm)), alpha))
    
21.         S_t += w * S_v + (1.0 - w) * S_r

复制代码

   Truhlar的qRRHO方法则更加简单：将低于阈值的振动模式频率统一设置为阈值。不过应当注意的是，Truhlar的qRRHO方法不仅仅会影响到振动熵，同样也会影响体系的总零点能。这里楼主直接给出代码：

1. def vibrational_entropy_contribution_qRRHO_Truhlar(vib_energies, temperature, low_freq):
   
2.     kT = units.kB * temperature
   
3.     S_v = 0.0
   
4.     for e in vib_energies:
   
5.         if ((e/units.invcm) < low_freq):
   
6.             e = low_freq * units.invcm
   
7.         x = e/kT
   
8.         S_v += x/(np.exp(x) - 1.0) - np.log(1.0 - np.exp(-x))
   
9.     S_v *= units.kB
   
10.     return S_v
    
11. 
    
12. 
    
13. def vibrational_energy_contribution_qRRHO_Truhlar(vib_energies, temperature, low_freq):
    
14.     kT = units.kB * temperature
    
15.     dU = 0.0
    
16.     for e in vib_energies:
    
17.         if ((e/units.invcm) < low_freq):
    
18.             e = low_freq * units.invcm
    
19.         dU += e / (np.exp(e / kT) - 1.0)
    
20.     return dU

复制代码

   
    接下来我们讨论如何将实现了qRRHO方法的代码整合进ASE。通过阅读ASE软件包中的thermochemistry.py文件，我们发现用于计算理想气体近似下的热力学性质的类IdealGasThermo是类ThermoChem衍生出来的。而IdealGasThermo类中求解热力学性质的函数get_gibbs_energy则调用了类ThermoChem中计算熵的函数_vibrational_entropy_contribution与计算振动能的函数_vibrational_energy_contribution。因此我们只需要将这两个函数更新为我们前文中所定义的函数，即可使用qRRHO方法计算准谐振近似下的热力学量了。我们首先将ASE提供的热力学计算模块导入为origin。随后，我们把origin.IdealGasThermo变成一个我们自己的类，然后在这个类中更改_vibrational_entropy_contribution与_vibrational_energy_contribution两个函数的定义。这样，当我们调用自己的类中的get_gibbs_energy函数时，这一函数就能够调用我们定义的包含了qRRHO方法的_vibrational_entropy_contribution与_vibrational_energy_contribution两个函数了。

1. from ase import thermochemistry as origin
   
2. 
   
3. class IdealGasThermo(origin.IdealGasThermo):
   
4.     def __init__(self, vib_energies, geometry, potentialenergy=0.,
   
5.                          atoms=None, symmetrynumber=None, spin=None, natoms=None,
   
6.                          method=None, low_freq=100):
   
7.         self.method = method
   
8.         self.low_freq = low_freq
   
9.         super(IdealGasThermo, self).__init__(vib_energies=vib_energies,
   
10.                   geometry=geometry, potentialenergy=potentialenergy, atoms=atoms,
    
11.                   symmetrynumber=symmetrynumber, spin=spin)
    
12. 
    
13.     def _vibrational_entropy_contribution(self, temperature):
    
14.         if self.method == 'Grimme':
    
15.             return vibrational_entropy_contribution_qRRHO_Grimme(self.vib_energies,
    
16.                        temperature, self.low_freq)
    
17.         elif self.method == 'Truhlar':
    
18.             return vibrational_entropy_contribution_qRRHO_Truhlar(self.vib_energies,
    
19.                        temperature, self.low_freq)
    
20.         else:
    
21.             return origin.ThermoChem._vibrational_entropy_contribution(self, temperature)
    
22. 
    
23.     def _vibrational_energy_contribution(self, temperature):
    
24.         if self.method == 'Truhlar':
    
25.             return vibrational_energy_contribution_qRRHO_Truhlar(self.vib_energies,
    
26.                        temperature, self.low_freq)
    
27.         else:
    
28.             return origin.ThermoChem._vibrational_energy_contribution(self, temperature)

复制代码

    到目前为止我们的代码已经可以正常工作了。但当下我们的代码还没法自动判断分子的转动对称数，而是需要用户通过symmetrynumber关键词手动指定。而为了得到分子的转动对称数，我们首先需要判断分子的点群。在Sob老师的Shermo程序中，分子点群的判断似乎是通过SYVA这个库进行实现的（这里不得不赞扬一句Sob老师的Fortran水平太高，楼主根本没读懂是怎么调用库的），而Python语言没法像Fortran那样便捷地与这些常用的库相连（可能这也是Sob老师力荐计算化学工作者学习Fortran的原因之一吧）。但万幸的是，Python库PyMatGen直接提供了判断分子点群的工具，因此我们可以直接使用PyMatGen完成我们的任务。我们仅仅需要将ASE的ATOMS类转换为PyMatGen的Molecule类，随后就可以使用PyMatGen的PointGroupAnalyzer功能进行点群判断了。借助这些成熟的软件包，我们可以仅仅使用10行代码完成分子点群判断：

1. from pymatgen.core.structure import Molecule
   
2. from pymatgen.symmetry.analyzer import PointGroupAnalyzer
   
3. 
   
4. def get_sch_symbol(atoms, tolerance, eigen_tolerance, matrix_tol):
   
5.     symbols = atoms.get_chemical_symbols()
   
6.     positions = atoms.get_positions()
   
7.     mol = Molecule(symbols, positions)
   
8.     # Analyze the point group of a molecule
   
9.     pointgroup = PointGroupAnalyzer(mol, tolerance, eigen_tolerance, matrix_tol)
   
10.     # Get Schoenflies symbol
    
11.     sch_symbol = pointgroup.sch_symbol
    
12.     return sch_symbol

复制代码

    得到分子点群之后，我们还需要将点群转化为转动对称数。由于各种点群对应的转动对称数的不变的，因此我们可以直接将各个点群所对应的转动对称数[2]以字典的方式存储在程序当中，随后只要通过点群名称对转动对称数进行索引即可。但这样做会导致一个问题：相同的对称操作可能会对应无穷多个点群（例如C2，C3，...），我们可能需要无限长的字典才能够覆盖到所有可能的情况（尽管实际分子不会具有那么高的对称性，但这样做仍然会导致函数长度增加，不够简洁和优雅）。但我们又知道，点群名称与分子转动对称数间还存在着一系列直接的关系，例如CN点群分子的转动对称数就是N，因此我们可以利用正则表达式抓取点群名称中的数字，从而辅助检索转动对称数：

1. def get_symmetrynumber(sch_symbol):
   
2.     if sch_symbol in ['C1', 'Ci', 'Cs', 'C*v']:
   
3.         symmetrynumber = 1
   
4.     elif re.match('C\d+', sch_symbol):
   
5.         symmetrynumber = int(re.findall('C(\d+)', sch_symbol)[0])
   
6.     elif sch_symbol == 'D*h':
   
7.         symmetrynumber = 2
   
8.     elif re.match('D\d+', sch_symbol):
   
9.         symmetrynumber = int(re.findall('D(\d+)', sch_symbol)[0])*2
   
10.     elif sch_symbol in ['T', 'Td']:
    
11.         symmetrynumber = 12
    
12.     elif re.match('S\d+', sch_symbol):
    
13.         symmetrynumber = int(re.findall('S(\d+)', sch_symbol)[0])//2
    
14.     elif sch_symbol == 'Oh':
    
15.         symmetrynumber = 24
    
16.     elif sch_symbol == 'Ih':
    
17.         symmetrynumber = 60
    
18.     else:
    
19.         symmetrynumber = 1
    
20.     return symmetrynumber

复制代码

    而最后，我们只需要将所有的函数组合起来，就可以得到一套完整而易用的热力学性质计算代码了。成品代码文件已经上传到了附件当中。这里楼主将代码命名为了HeatMeter，取“量热器”之意。使用这段代码不需要进行任何引用，如果实在需要引文，引用ASE以及PyMatGen软件包即可。顺带一提，楼主的代码本体中并没有包含从各种量子化学软件的输出文件中获取频率信息的代码，因为各种软件格式并不统一，而且利用Python的正则表达式功能显然能够十分轻易地实现这样的读取。楼主仅仅给出一个从高斯log文件中读取频率的函数（包含在example.py当中）：

1. from ase import io
   
2. 
   
3. def read_gaussian_out(input_name):
   
4.     atoms = io.read(input_name, format='gaussian-out')
   
5.     f = open(input_name)
   
6.     lines = f.read()
   
7.     f.close()
   
8.     multiplicity = re.findall('Charge =\s*\d+ Multiplicity =\s*(\d+)', lines)
   
9.     multiplicity = int(multiplicity[-1])
   
10.     vib_energies = re.findall('Frequencies --(.*)', lines)
    
11.     vib_energies = np.array(' '.join(vib_energies).split(), dtype=float)
    
12.     low_freq = len(atoms)*3-len(vib_energies)
    
13.     vib_energies = np.append(np.zeros(low_freq, dtype=float), vib_energies)
    
14.     atoms.potential_energy = atoms.get_potential_energy()
    
15.     atoms.vib_energies = vib_energies
    
16.     atoms.multiplicity = multiplicity
    
17.     return atoms

复制代码

测试案例

    这里我们以Sob老师提供的乙基自由基为例，比较我们的代码的结果与Shemo的结果是否相同。当使用Grimme的qRRHO方法时，我们的代码与Shermo算得的体系的熵值分别为：254.933 J/mol/K与254.920 J/mol/K，相对误差为0.005%；而对体系自由能的总校正量分别为：93.577 kJ/mol 与93.581 kJ/mol，相对误差为-0.004%。当仅使用RRHO近似时，我们的代码与Shermo算得的体系的熵值分别为：255.512 J/mol/K与255.494 J/mol/K，相对误差为0.007%；而对体系自由能的总校正量分别为：93.404 kJ/mol 与93.410 kJ/mol，相对误差为0.006%。因此可见我们的代码应当是正确的（数值误差应当由物理常数选取的不同引起）。如果你需要运行这个例子，下载附件，解压，并在中终端输入python example.py即可。

结语

    在本文当中，楼主仿照Sob老师的Shermo代码的功能，简单实现了带有qRRHO方法的热力学计算代码，并实现了对于分子转动对称数的判断。但由于需要依赖于Python的软件包生态，楼主这段短短100行的代码的易用性确实不远远如Sob老师独立而稳定的Shermo程序，尽管它们所解决的问题是相同的（这里还是要感谢卢老师为众多计算化学工作者所耗费的大量心血，只有真正读过、写过一点代码才能够体会这背后所要付出的精力之大）。尽管如此，楼主还是认为本文的启发意义大于文中代码实用价值，也希望本文能够对想要实现一些简单脚本计算的同学带来一定的思考，可能有一些计算看起来困难，但是实际上手做起来并不真的复杂。

参考文献

[1] Grimme, Stefan. "Supramolecular binding thermodynamics by dispersion‐corrected density functional theory." Chemistry–A European Journal 18.32 (2012): 9955-9964.
[2] https://cccbdb.nist.gov/thermo.asp

风起~

pymatgen居然可以判断点群，好强大
谢谢楼主 学习了

WaterMirror

请问这个帖子是只是介绍如何写代码 还是 楼主的代码和Shermo相比有显著的优势呢？
因为我觉得Shermo已经很好了....

ByTheWay

WaterMirror 发表于 2022-5-31 10:45
**** 本内容被作者隐藏 ****

楼主当前积分不够，所以托组里师弟看了一眼你的留言。

请问这个帖子是只是介绍如何写代码 还是 楼主的代码和Shermo相比有显著的优势呢？
因为我觉得Shermo已经很好了....

楼主似乎从来没说过自己的代码比Shermo更好吧？这个帖子的目的只是和大家聊一聊楼主写码的一点点心得，希望能够给同样希望写一些代码，但是又不知道从何下手的坛友一点点灵感。如果你只是想把数据算出来，你当然应该去用shermo。

njfuzjs

WaterMirror 发表于 2022-5-31 10:45
**** 本内容被作者隐藏 ****

按照你的说法，Gaussian已经把SCF写了一遍，那么为什么还有必要写其他量子化学代码呢？

WaterMirror

ByTheWay 发表于 2022-5-31 15:14
楼主当前积分不够，所以托组里师弟看了一眼你的留言。

楼主似乎从来没说过自己的代码比Shermo更好吧？ ...

好的 谢谢您！
只是对代码效率的问题好奇，没有别的意思。

wzkchem5

WaterMirror 发表于 2022-5-31 11:20
好的 谢谢您！
只是对代码效率的问题好奇，没有别的意思。

这种代码没什么效率方面的问题，因为计算量极其小，不管写得多差都是瞬间算完。

WaterMirror

wzkchem5 发表于 2022-5-31 19:15
这种代码没什么效率方面的问题，因为计算量极其小，不管写得多差都是瞬间算完。

谢谢王老师！

卡开发发

wzkchem5 发表于 2022-5-31 19:15
这种代码没什么效率方面的问题，因为计算量极其小，不管写得多差都是瞬间算完。

计算量小只是一个方面，python这样的语言比较胶水，得看使用的库和底层实现的逻辑，比方说循环密集处及矢量计算都有相应的加速手段。反过来说，即便使用偏底层的语言也不应该不计数据结构和代码风格胡写一通，合适的场景选择合适的工具才是明智的。

卡开发发

例子中一处设计不当，实际上那个low_freq不应该放在read_gaussian_out下面，而应该单独放在外面，这样IdealGasThermo实例化就可以传入。不过这个想法确实很好，比方说我可以很简单实现一个vasp计算小分子热化学数据的接口：

1. from ase import io
   
2. from HeatMeter import IdealGasThermo, get_geometry
   
3. import numpy as np
   
4. import re
   
5. 
   
6. def read_vasp_out(input_name):
   
7.         atoms=io.read(input_name,format='vasp-out')
   
8.         atoms.potential_energy=atoms.get_potential_energy()
   
9.         magmom=atoms.get_calculator().results['magmom']
   
10.         atoms.multiplicity=int(np.round(magmom))+1
    
11.         f=open(input_name)
    
12.         lines=f.read()
    
13.         f.close()
    
14.         vib_energies=re.findall('PiTHz\s+(.*)\s+cm-1',lines)
    
15.         vib_energies=np.array(vib_energies,dtype=float)[::-1]
    
16.         atoms.vib_energies=vib_energies
    
17.         moltype=get_geometry(atoms)
    
18.         if moltype=='monatomic':
    
19.                 low_freq=0
    
20.         elif moltype=='linear':
    
21.                 low_freq=5
    
22.         else:
    
23.                 low_freq=6
    
24.         return atoms,low_freq
    
25. 
    
26. input_name='OUTCAR'
    
27. atoms,low_freq=read_vasp_out(input_name)
    
28. thermo=IdealGasThermo(atoms,method=None)
    
29. thermo.get_gibbs_energy(temperature=298.15,pressure=101325.)

复制代码

不过9～10行对多重度的处理以及21行对频率的读取我不能保证这是最合适的方案，如果大家有兴趣可以按照类似格式稍微调整下。其实其他例如CP2K或者DMol3的接口也可以简单实现，不过鉴于CP2K使用太复杂，这里没有提供相关方案。

hebrewsnabla

其实判断点群，得到转动对称数和热力学校正量这些东西在PySCF里面也有。只不过没有qRRHO。楼主有兴趣的话可以给他们提个PR贡献一下。