卡开发发 发表于 2022-9-12 20:21
上面输出G是correction的部分，实际总的Gtot可以这样算
Gtot=Gcorr+Eele，Eele是电子自洽场计算的能量。 ...

索嘎，谢谢老师

卡开发发 发表于 2022-9-12 22:14
没有，我说的其实比较零散且了解比较有限，如果你感兴趣可以看看Phonopy和ASE给的公式和资料，也许会有些 ...

老师自谦啦，谢谢您推荐的资料~

含光君 发表于 2022-9-12 22:01
发发老师讲得很明白，比学校统计热力学课清晰太多太多，非常受益！

没有，我说的其实比较零散且了解比较有限，如果你感兴趣可以看看Phonopy和ASE给的公式和资料，也许会有些帮助。

卡开发发 发表于 2022-9-12 20:21
上面输出G是correction的部分，实际总的Gtot可以这样算
Gtot=Gcorr+Eele，Eele是电子自洽场计算的能量。 ...

发发老师讲得很明白，比学校统计热力学课清晰太多太多，非常受益！

462855 发表于 2022-9-12 20:06
老师，一般G=E+EZPE-TS,这个H我不知道怎么加进去，能说明白一点吗？

上面输出G是correction的部分，实际总的Gtot可以这样算
Gtot=Gcorr+Eele，Eele是电子自洽场计算的能量。
G=Eele+ZPE-TS只是个近似公式，其中Eele+ZPE是近似的H(0)，其实少了H(0->T)的部分。
如果偷懒一点，上述代码中HarmonicThermo(vib_energies)当中其实可以改为HarmonicThermo(vib_energies,potential)，这个potential可以把电子的自洽场能量换算为eV单位代进去。

含光君 发表于 2022-9-12 13:12
把频率分析文件中的实频记录下来按下图公式计算，注意能量单位的转化。

老师，一般G=E+EZPE-TS,这个H我不知道怎么加进去，能说明白一点吗？

卡开发发 发表于 2022-9-12 18:02
用python的ase包就能很简单解决，给个简单的例子：

输出结果：

好厉害，谢谢老师

用python的ase包就能很简单解决，给个简单的例子：

1. import numpy as np
   
2. from ase.thermochemistry import HarmonicThermo
   
3. 
   
4. """
   
5. Ex:
   
6.      vibrational frequencies,                          intensities
   
7.   mode     au_amu        cm-1         meV         THz      km/mol
   
8.     1   -0.047478      -244.1      -30.26      -7.317      105.12
   
9.     2   -0.017344       -89.2      -11.05      -2.673       83.42
   
10.     3   -0.000089        -0.5       -0.06      -0.014        0.00
    
11.     4    0.000044         0.2        0.03       0.007        0.00
    
12.     5    0.000075         0.4        0.05       0.012        0.00
    
13.     6    0.017593        90.4       11.21       2.711       66.41
    
14.     7    0.314760      1618.0      200.61      48.507       76.57
    
15.     8    0.723107      3717.1      460.87     111.437       66.40
    
16.     9    0.745950      3834.6      475.43     114.957      457.09
    
17. """
    
18. # 读取meV那一列，并且去除前6个平动+转动自由度的贡献
    
19. vib_energies=[200.61,460.87,475.43]
    
20. # meV->eV
    
21. vib_energies=np.array(vib_energies)/1000
    
22. thermo=HarmonicThermo(vib_energies)
    
23. # 计算298.15K下的振动自由能（含ZPE）
    
24. FreeE=thermo.get_helmholtz_energy(298.15)

复制代码

输出结果：

1. Internal energy components at T = 298.15 K:
   
2. ===============================
   
3. E_pot                  0.000 eV
   
4. E_ZPE                  0.568 eV
   
5. Cv_harm (0->T)         0.000 eV
   
6. -------------------------------
   
7. U                      0.569 eV
   
8. ===============================
   
9. 
   
10. 
    
11. Entropy components at T = 298.15 K:
    
12. =================================================
    
13.                            S               T*S
    
14. S_harm             0.0000003 eV/K        0.000 eV
    
15. -------------------------------------------------
    
16. S                  0.0000003 eV/K        0.000 eV
    
17. =================================================
    
18. 
    
19. Free energy components at T = 298.15 K:
    
20. =======================
    
21.     U          0.569 eV
    
22. -T*S         -0.000 eV
    
23. -----------------------
    
24.     F          0.568 eV
    
25. =======================
    

复制代码

上面还有修改的余地，例如增设argparse来直接做成个外部读入文件和参数的程序，以及可以通过正则表达式re模块从outmol抓取频率信息。
另外涉及到的公式是2楼的公式，具体底层实现你可以看HarmonicThermo类其中的源码部分。

另外需要指明的问题是：

如果整个结构计算频率，会把所有原子考虑进去，这样的自由能修正肯定不准确

1、这个说法并不合理，如果假定整个体系是符合谐振近似的，那么应当计算热化学过程考虑整个体系的Hessian矩阵，即便不符合近似的条件，仅通过不考虑非活性区域部分对振动模式的贡献也不能得到正确的结果，甚至会影响到你吸附分子的振动模式（因为吸附分子会和表面的原子发生作用，振动模式很难解除耦合）。这在数学上比较容易理解，其实做partial Hessian只是相当于full Hessian的一个子空间，频率计算要对角化Hessian矩阵，若两者要得到相同的本征值则必须保证你关注的活性区域和其他非活性区域在Hessian矩阵中是可以分块的，比方说一个很弱的物理吸附也许能够满足这样的条件，否则可能需要考虑分子吸附时周围相互作用更多的原子，才有可能具备上述条件。很多文章这么算的原因不是准确与否这样的原因，而是振动频率计算代价确实很大。
2、上述公式成立的条件是振动态相对比较局域化，我们按照类似于分子那样去处理固体的振动，并且考虑其实固体表面是充分大。否则，整个计算也应当使用类似于超胞方法来进行完整的声子计算，在处理热化学数据的时候应该也要考虑布里渊区积分。

把频率分析文件中的实频记录下来按下图公式计算，注意能量单位的转化。