计算化学公社

标题: 求助Dmol3计算单原子催化剂NRR时，反应过程中如何计算自由能 [打印本页]

作者
Author: 462855 时间: 2022-9-12 10:19
标题: 求助Dmol3计算单原子催化剂NRR时，反应过程中如何计算自由能
如果整个结构计算频率，会把所有原子考虑进去，这样的自由能修正肯定不准确
但是选择partial Hessian计算时，out文件又不会输出最后的自由能分析。
求助计算吸附后分子自由能修正计算的方法。

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

(, 下载次数 Times of downloads: 22)

作者
Author: 卡开发发 时间: 2022-9-12 18:02
本帖最后由卡开发发于 2022-9-12 18:30 编辑

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

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

复制代码

输出结果：

Internal energy components at T = 298.15 K:
===============================
E_pot 0.000 eV
E_ZPE 0.568 eV
Cv_harm (0->T) 0.000 eV
-------------------------------
U 0.569 eV
===============================
Entropy components at T = 298.15 K:
=================================================
S T*S
S_harm 0.0000003 eV/K 0.000 eV
-------------------------------------------------
S 0.0000003 eV/K 0.000 eV
=================================================
Free energy components at T = 298.15 K:
=======================
U 0.569 eV
-T*S -0.000 eV
-----------------------
F 0.568 eV
=======================

复制代码

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

另外需要指明的问题是：

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

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

作者
Author: 462855 时间: 2022-9-12 20:05

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

输出结果：

好厉害，谢谢老师

作者
Author: 462855 时间: 2022-9-12 20:06

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

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

作者
Author: 卡开发发 时间: 2022-9-12 20:21
本帖最后由卡开发发于 2022-9-12 20:22 编辑

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单位代进去。

作者
Author: 含光君 时间: 2022-9-12 22:01

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

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

作者
Author: 卡开发发 时间: 2022-9-12 22:14

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

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

作者
Author: 含光君 时间: 2022-9-12 23:56

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

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

作者
Author: 462855 时间: 2022-9-13 08:21

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

索嘎，谢谢老师

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