捐赠本论坛 (Donate this site)
开启辅助访问 切换到宽版

计算化学公社

 找回密码 Forget password
 注册 Register
“第13届北京科音中级量子化学培训班” 将于7月18至23日于北京举行，内容极为浩瀚，幻灯片超过2700页，干货满满，相较往届又增加了机器学习势等新内容。这是已经具备一定研究经验者利用暑期时间全面、系统提升量子化学研究水平和显著扩展知识的一年一度的非常难得的机会。报名正在进行中，请点击此链接查看，欢迎参加！
计算化学公社»论坛首页 › 理论与计算化学 (Theoretical and Computational Chemistry) › 第一性原理 (First-principles) › Quantum ESPRESSO吸附物振动自由能校正：DFPT线性响应声 ...
返回列表 Return 发新帖
Views: 402|回复 Reply: 3
打印 Print 上一主题 Last thread 下一主题 Next thread

[Quantum ESPRESSO] Quantum ESPRESSO吸附物振动自由能校正：DFPT线性响应声子法与有限位移法对比

[复制链接 Copy URL]
alkylaix

5

帖子

2

威望

513

eV
积分
558

Level 4 (黑子)
跳转到指定楼层 Go to specific reply
楼主
发表于 Post on 2026-6-20 18:10:59 | 只看该作者 Only view this author |只看大图 View big picture 回帖奖励 |倒序浏览 Reverse view |阅读模式 Reading model
本帖最后由 alkylaix 于 2026-6-24 02:09 编辑

在相同结构和相同低频处理规则下，Quantum ESPRESSO的DFPT声子法与有限位移法得到的CO2局域振动热校正非常接近：ZPE差异约0.0038 eV，G_vib差异约0.0058 eV。

按照帖子中的思路的辅助程序FDVIB已用C++实现。提供预编译的Linux x86_64版本，由于目前程序的输入输出模式，以及系统性的验证仍在讨论和完善中，待我们完善后，计划新发一个帖子系统的介绍使用方法，Github：文档  仓库

1. 目的

最近在用 Quantum ESPRESSO 处理气体分子在表面的物理吸附自由能校正。常规思路是用ph.x和dynmat.x做局域振动分析；这里我尝试了另一条路线：用有限位移法直接从QE受力构造吸附物局域 Hessian，然后与ph.x + dynmat.x的结果做一个对比。
测试体系是 CO2 在二维 Cu-MOF 表面的物理吸附。吸附构型已在相同计算参数下完成结构优化，优化收敛阈值为：
  1. etot_conv_thr =   1.000000000d-05
  2. forc_conv_thr =   1.000000000d-03
复制代码
在该优化结构基础上，本报告比较两种 Quantum ESPRESSO 振动分析流程得到的零点能（ZPE）和振动自由能校正（G_vib）：
  • 基于Quantum ESPRESSO ph.x + dynmat.x 的DFPT声子法；
  • 基于Quantum ESPRESSO Hellmann-Feynman力的有限位移法局域Hessian频率计算。
比较重点不是完整声子，而是吸附物局域振动对热力学校正的影响。两套后处理均采用298.15 K，并使用相同的低频、零频和虚频处理规则。

2. 频率处理规则

频率类型处理方式
虚频，ν < 0 cm^-1 跳过，不参与 ZPE 和 G_vib
零频跳过，不参与 ZPE 和 G_vib
正低频，0 < ν < 50 cm^-1设为 50 cm^-1
正频，ν >= 50 cm^-1保留原值

热校正公式为：

其中，vi为经过上述规则处理后参与热力学校正的正频率波数。严格来说，上式是谐振子近似下的振动Helmholtz自由能F_vib；对于固体表面/吸附体系，通常忽略pV项，并将其近似作为振动Gibbs自由能校正G_vib使用。本文后文沿用G_vib这一记号。

3. 计算设置说明
本次比较需要保持DFPT声子法和有限位移法两条路线的力常数来源尽量一致。
  • DFPT路线使用ph.in中的nat_todo = 3，只对CO2吸附物三个原子的9个扰动进行响应计算；dynmat.x仍按完整体系维度3N = 201输出频率，因此中间0.000 cm^-1 来自未计算的环境自由度，不应理解为真实物理零频。本报告只使用CO2局域振动块中按规则保留下来的正频率进行热力学校正。
  • 有限位移法通过对原子施加小位移，并由SCF计算得到的Hellmann-Feynman力进行有限差分，从而构造二阶力常数矩阵。该路线只选取吸附物的3个原子，因此共有9个局域自由度：

原子编号元素质量
65
C12.011
66O15.999
67O15.999


需要特别说明的是，这并不是因为pw.x / pwscf不支持DFT-D3的three-body term。相反，在Quantum ESPRESSO的DFT-D3总能量和力计算中，three-body term默认是开启的。这里显式设置 dftd3_threebody = .false.，是因为three-body term对应的二阶导数在ph.x / d3hess.x的声子或 Hessian 流程中没有实现。为了避免DFPT与有限位移法之间引入额外的色散Hessian差异，相关计算均关闭 three-body term。

赝势使用JTH v2.0 GGA-PBE PAW atomicdata。涉及的元素和赝势文件为：

元素赝势文件
CC.GGA_PBE-JTH.UPF
CuCu.GGA_PBE-JTH.UPF
HH.GGA_PBE-JTH.UPF
NN.GGA_PBE-JTH.UPF
OO.GGA_PBE-JTH.UPF

4. 输入文件
位移法实际计算中还包含由fd/scf.in生成的 18 个位移输入文件，它们仅改变outdir和被位移原子的一个坐标，不逐个展开。
4.1 ph_dynmat/scf.in
  1. &CONTROL
  2.   calculation = 'scf'
  3.   etot_conv_thr =   1.000000000d-05
  4.   forc_conv_thr =   1.000000000d-03
  5.   outdir = './out/'
  6.   prefix = 'ads'
  7.   tprnfor = .true.
  8.   pseudo_dir = '~/pseudo/JTH/JTH-v2.0-GGA-PBE-atomicdata-upf'
  9.   verbosity = 'low'
  10. /
  11. &SYSTEM
  12.   degauss =   3.0000000000d-03
  13.   ecutrho =   4.0000000000d+02
  14.   ecutwfc =   5.0000000000d+01
  15.   ibrav = 0
  16.   nat = 67
  17.   nosym = .true.
  18.   ntyp = 5
  19.   vdw_corr = 'dft-d3'
  20.   dftd3_version = 4
  21.   dftd3_threebody = .false.
  22.   occupations = 'smearing'
  23.   smearing = 'gauss'
  24.   nspin = 2
  25.   starting_magnetization(1) = 0
  26.   starting_magnetization(2) = 0
  27.   starting_magnetization(3) = 0
  28.   starting_magnetization(4) = 0
  29.   starting_magnetization(5) = 0
  30. /
  31. &ELECTRONS
  32.   conv_thr =   1.0000000d-09
  33.   electron_maxstep = 80
  34.   mixing_beta =   4.0000000000d-01
  35. /
  36. ATOMIC_SPECIES
  37. C   12.011  C.GGA_PBE-JTH.UPF
  38. Cu  63.546  Cu.GGA_PBE-JTH.UPF
  39. H    1.008  H.GGA_PBE-JTH.UPF
  40. N   14.007  N.GGA_PBE-JTH.UPF
  41. O   15.999  O.GGA_PBE-JTH.UPF
  42. ATOMIC_POSITIONS (angstrom)
  43. ! ... 64 atoms omitted ...
  44. C         7.9027208554      8.2451403636     18.1111719209
  45. O         8.7559409272      7.4441712661     18.0925952630
  46. O         7.0495099078      9.0464615449     18.1356223639
  47. K_POINTS automatic
  48. 1 1 1 0 0 0
  49. CELL_PARAMETERS (angstrom)
  50.      16.2978854140      0.0012415390      0.0000000000
  51.      -0.0018959120     16.4904555690      0.0000000000
  52.       0.0000000000      0.0000000000     30.0000000000
复制代码

4.2 ph_dynmat/d3hess.in
  1. &INPUT
  2.   prefix = 'ads'
  3.   outdir = './out'
复制代码
4.3 ph_dynmat/ph.in
  1. &INPUTPH
  2.   prefix = 'ads'
  3.   outdir = './out'
  4.   fildyn = 'ads.dynG'
  5.   tr2_ph = 1.0d-14
  6.   nat_todo = 3
  7.   nogg = .true.
  8.   dftd3_hess = 'ads.hess'
  9.   reduce_io = .true.
  10. /
  11. 0.0 0.0 0.0
  12. 65 66 67
复制代码

4.4 ph_dynmat/dynmat.in
  1. &INPUT
  2.   fildyn = 'ads.dynG'
  3.   filout = 'ads.freq.out'
  4.   asr = 'no'
  5.   remove_interaction_blocks = .true.
  6.   filxsf = 'ads_modes.axsf'
  7.   filmol = 'ads_modes.mold'
复制代码
4.5 fd/scf.in
有限位移法的基准输入与DFPT的SCF参数保持一致。18个位移输入由该基准结构沿原子65-67的 x/y/z 方向分别做 +/- 0.01 Angstrom 位移生成。
  1. &CONTROL
  2.   calculation = 'scf'
  3.   etot_conv_thr =   1.000000000d-05
  4.   forc_conv_thr =   1.000000000d-03
  5.   outdir = './out/'
  6.   prefix = 'ads'
  7.   tprnfor = .true.
  8.   pseudo_dir = '~/pseudo/JTH/JTH-v2.0-GGA-PBE-atomicdata-upf'
  9.   verbosity = 'low'
  10. /
  11. &SYSTEM
  12.   degauss =   3.0000000000d-03
  13.   ecutrho =   4.0000000000d+02
  14.   ecutwfc =   5.0000000000d+01
  15.   ibrav = 0
  16.   nat = 67
  17.   nosym = .true.
  18.   ntyp = 5
  19.   vdw_corr = 'dft-d3'
  20.   dftd3_version = 4
  21.   dftd3_threebody = .false.
  22.   occupations = 'smearing'
  23.   smearing = 'gauss'
  24.   nspin = 2
  25.   starting_magnetization(1) = 0
  26.   starting_magnetization(2) = 0
  27.   starting_magnetization(3) = 0
  28.   starting_magnetization(4) = 0
  29.   starting_magnetization(5) = 0
  30. /
  31. &ELECTRONS
  32.   conv_thr =   1.0000000d-09
  33.   electron_maxstep = 80
  34.   mixing_beta =   4.0000000000d-01
  35. /
  36. ATOMIC_SPECIES
  37. C   12.011  C.GGA_PBE-JTH.UPF
  38. Cu  63.546  Cu.GGA_PBE-JTH.UPF
  39. H    1.008  H.GGA_PBE-JTH.UPF
  40. N   14.007  N.GGA_PBE-JTH.UPF
  41. O   15.999  O.GGA_PBE-JTH.UPF
  42. ATOMIC_POSITIONS (angstrom)
  43. ! ... 64 atoms omitted ...
  44. C         7.9027208554      8.2451403636     18.1111719209
  45. O         8.7559409272      7.4441712661     18.0925952630
  46. O         7.0495099078      9.0464615449     18.1356223639
  47. K_POINTS automatic
  48. 1 1 1 0 0 0
  49. CELL_PARAMETERS (angstrom)
  50.      16.2978854140      0.0012415390      0.0000000000
  51.      -0.0018959120     16.4904555690      0.0000000000
  52.       0.0000000000      0.0000000000     30.0000000000
复制代码
有限位移输入只改变被位移原子的一个坐标。例如原子65沿x方向正位移+0.01 Angstrom时，CO2区域为：
  1. ATOMIC_POSITIONS (angstrom)
  2. ! ... 64 atoms omitted ...
  3. C         7.9127208554      8.2451403636     18.1111719209
  4. O         8.7559409272      7.4441712661     18.0925952630
  5. O         7.0495099078      9.0464615449     18.1356223639
复制代码

5. Quantum ESPRESSO ph.x + dynmat.x 结果
ph.x + dynmat.x 输出仍保留完整体系的3N = 201个频率位置。由于ph.in只计算原子65-67的局域扰动，表中的零频主要对应未计算的环境自由度或被移除的相互作用块。热力学校正实际使用的模式如下：
moderaw frequency / cm^-1used frequency / cm^-1处理方式
1-137.479skippedimag
2-72.324skippedimag
3-47.910skippedimag
4-1950.000skippednot computed / zero
19642.88050.000set to 50
19761.83661.836keep
198605.868605.868keep
199610.253610.253keep
2001326.8031326.803keep
2012369.0972369.097keep
得到：
项目数值
ZPE0.31143944 eV
ZPE7.1819 kcal/mol
ZPE30.0493 kJ/mol
G_vib0.23421412 eV
G_vib5.4011 kcal/mol
G_vib22.5982 kJ/mol

6. Quantum ESPRESSO 有限位移法局域Hessian结果

有限位移法频率来自对吸附物 3 个原子的 +/- 0.01 Angstrom 位移受力差分。理想情况下，由二阶导数构造的Hessian应满足对称性H_ij = H_ji。本次原始Hessian在对称化前的最大非对称项为 4.509226e+01 eV/Angstrom^2，说明有限位移得到的力常数矩阵存在明显数值不对称。我认为来源可能包括SCF力收敛、位移步长、结构残余力、局域Hessian截断近似等误差。本文最终使用对称化后的 Hessian 求频率。
不知道CO2的振动模式会不会对Hessian中的小误差是否敏感？

moderaw frequency / cm^-1used frequency / cm^-1处理方式
1-262.069skippedimag
2-60.048skippedimag
3-32.580skippedimag
450.28250.282keep
555.99155.991keep
6600.617600.617keep
7623.803623.803keep
81312.5941312.594keep
92319.3882319.388keep
得到：
项目数值
Raw modes9
Used modes6
ZPE0.30764659 eV
ZPE7.0945 kcal/mol
ZPE29.6834 kJ/mol
G_vib0.22839207 eV
G_vib5.2668 kcal/mol
G_vib22.0365 kJ/mol

7. 结果对比

项目Quantum ESPRESSO DFPT 声子法Quantum ESPRESSO 有限位移法有限位移法 - DFPT
Used modes660
ZPE / eV0.311439440.30764659-0.00379285
ZPE / kJ mol^-130.049329.6834-0.3659
G_vib / eV0.234214120.22839207-0.00582205
G_vib / kJ mol^-122.598222.0365-0.5617
两种方法的ZPE差异约为0.0038 eV，G_vib差异约为0.0058 eV。在常见吸附自由能分析尺度上，这属于较小差异；但有限位移法的G_vib仍比DFPT更低，主要来自低频模式对熵项的放大作用。

8. 讨论
8.1 高频区域
CO2 内部键振动对应较高的振动频率。两套方法在这些模式上比较接近：
模式DFPT 声子法 / cm^-1有限位移法 / cm^-1差值 / cm^-1
内部振动 1605.868600.617-5.251
内部振动 2610.253623.80313.550
内部振动 31326.8031312.594-14.209
内部振动 42369.0972319.388-49.709
8.2 低频区域
G_vib对低频模式更敏感。DFPT的两个最低正频率为42.880 cm^-1 和61.836 cm^-1，对应进入热校正的数值为50.000 cm^-1 和61.836 cm^-1；有限位移法的两个最低正频率为50.282 cm^-1和55.991 cm^-1。有限位移法的第二个低频更低，因此熵项贡献更负，使G_vib略低。
8.3 虚频模式
两套方法都出现3个虚频。对振动模式进行可视化检查后，这3个虚频可归属为整体相对表面的运动：
虚频模式模式特征
1吸附物整体水平平移
2吸附物整体水平平移
3吸附物绕吸附点摆动，运动方向几乎水平
需要注意的是，虚频被跳过只是本报告采用的一种热力学后处理约定。我也不确定这样处理是否最合理，欢迎大家进行讨论，因为这些吸附物整体平移/摆动软模式可能仍然对应低频受限平动/转动熵，对吸附自由能产生影响。

9. 总结
在统一的频率处理规则下，Quantum ESPRESSO的DFPT声子法与有限位移法局域Hessian给出的热力学校正非常接近：
方法ZPE / eVG_vib / eV
Quantum ESPRESSO DFPT 声子法0.311439440.23421412
Quantum ESPRESSO 有限位移法0.307646590.22839207

  • 两种方法的 ZPE 差异约0.0038eV。
  • 两种方法的G_vib 差异约0.0058eV。
  • 对该吸附物局域热校正而言，两种方法在当前精度下可以认为是一致的。

Python程序说明
我会把用到的Python程序一并上传。这里简单说明每个脚本的作用

脚本位置作用
make_fd_inputs.pyfd/从 fd/scf.in 读取未位移结构，对选定原子 65, 66, 67 沿 x/y/z 方向分别生成 +/- 0.01 Angstrom 位移输入文件。脚本会自动修改对应原子坐标和 outdir。
fd_thermo_qe.pyfd/读取有限位移计算得到的 Quantum ESPRESSO 输出力，构造 9 x 9 局域 Hessian；对 Hessian 对称化后求频率，并按本文规则处理虚频、零频和低频，最后输出 ZPE、G_vib、频率列表和 Hessian 文件。
make_fd_mode_xyz.pyfd/读取 fd_hessian_eV_A2.dat，重新求有限位移法的频率和本征矢，并把指定模式写成 .xyz 动画轨迹，方便用 OVITO、VESTA 或 Jmol 检查虚频/低频模式的振动方向。
ph_dynmat_thermo.pyph_dynmat/从 dynmat.x 输出的 ads.freq.out 中解析频率，按同样的虚频、零频和低频规则筛选模式，并计算 DFPT 路线对应的 ZPE 和 G_vib。

请大家指正！



本主题由 sobereva 于 2026-6-20 20:08 解除精华

benchmark_fd_ph.zip

817.52 KB, 下载次数 Times of downloads: 8

评分 Rate

参与人数
Participants 4		威望 +1 eV +7 收起 理由
Reason
1xiaohua1 + 1 赞！
丁越 + 5 赞！
asdf + 1 谢谢
sobereva + 1

查看全部评分 View all ratings

回复 Reply

举报 Report

sobereva

6万

帖子

99

威望

6万

eV
积分
128409

管理员

公社社长
2#
发表于 Post on 2026-6-20 20:08:26 | 只看该作者 Only view this author
其中对应公式的图片贴图方式有问题，请编辑帖子重新上传
北京科音自然科学研究中心（http://www.keinsci.com）致力于计算化学的发展和传播，长期开办极高质量的各种计算化学类培训：初级量子化学培训班、中级量子化学培训班、高级量子化学培训班、量子化学波函数分析与Multiwfn程序培训班、分子动力学与GROMACS培训班、CP2K第一性原理计算培训班，内容介绍以及往届资料购买请点击相应链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的高速路！培训各种常见问题见《北京科音办的培训班FAQ》。
欢迎加入北京科音微信公众号获取北京科音培训的最新消息，并避免错过网上有价值的计算化学文章！
欢迎加入人气极高、专业性特别强的理论与计算化学综合交流群思想家公社QQ群（群号见此链接），合计达一万多人。北京科音培训班的学员在群中可申请VIP头衔，提问将得到群主Sobereva的最优先解答。
思想家公社的门口Blog：http://sobereva.com（发布大量原创计算化学相关博文）
Multiwfn主页：http://sobereva.com/multiwfn（十分强大、极为流行的量子化学波函数分析程序）
Google Scholar：https://scholar.google.com/citations?user=tiKE0qkAAAAJ
ResearchGate：https://www.researchgate.net/profile/Tian_Lu
回复 Reply

举报 Report

alkylaix

5

帖子

2

威望

513

eV
积分
558

Level 4 (黑子)
3#
 楼主 Author| 发表于 Post on 2026-6-20 20:27:00 | 只看该作者 Only view this author
sobereva 发表于 2026-6-20 20:08
其中对应公式的图片贴图方式有问题，请编辑帖子重新上传

谢谢提醒，已重新上传
回复 Reply

举报 Report

1xiaohua1

39

帖子

0

威望

125

eV
积分
164

Level 3 能力者
4#
发表于 Post on 2026-6-23 10:50:33 | 只看该作者 Only view this author
赞！ 学习一下
回复 Reply

举报 Report

返回列表 Return 发新帖

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2026-7-17 07:03 , Processed in 0.225188 second(s), 28 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list