使用FDVIB控制Quantum ESPRESSO高效计算局域振动，并用Shermo完成热化学后处理

alkylaix

在先前的《Quantum ESPRESSO吸附物振动自由能校正：DFPT线性响应声子法与有限位移法对比》一文中，我们用有限位移差分构造吸附物局域Hessian与ph.x结果进行了对比。本文在此基础上介绍FDVIB程序，展示如何把这一流程更简单方便的用于周期体系中的局域振动分析和热化学后处理。

FDVIB是一个用于局域振动分析的有限位移差分工具。它通过调用Quantum ESPRESSO[1,2]的pw.x计算正、负位移结构的原子力，然后由力的差分构造局域Hessian，并进一步生成QE兼容的Γ点动力学矩阵，可直接用dynmat.x得到频率（这样能够兼容ph.x+dynmat.x的结果进行后处理）。还能进一步将dynmat.x的结果转为Molden格式并结合MfakeG查看振动模式，生成 .shm文件使用Shermo[3]进行化学后处理

FDVIB最适合这类问题：周期体系，但真正关心的只是某个吸附物、反应中间体或局部结构的振动。例如表面催化、分子筛催化等。此时如果只需要局域振动频率进行热化学校正，就不一定要对整个周期体系做完整声子计算。此时用FDVIB能较快的得到频率，加上直接输出 .shm更是方便热化学后处理！

FDVIB支持两类计算：

• local：用于周期体系或嵌入模型。只位移selected_atoms指定的原子，并构造这些原子的局域Hessian。
• gas：用于孤立分子。位移所有原子，构造完整分子Hessian，并在后处理中去除平动和转动。
  

本文以CO吸附在H-LTA分子筛中的模型为例，介绍FDVIB的完整计算流程，并与QE的ph.x结果进行对照。

---

1.算例说明

本文以CO@H-LTA为例，计算吸附CO的局域振动。我们只选择第74、75号O、C原子，在完整周期环境中构造6×6局域Hessian，得到6个局域振动模式。其余分子筛骨架原子不参与位移。



计算时，FDVIB先完成未位移结构的参考SCF，再将收敛的电荷密度作为所有正、负位移任务的共同初猜。这样可以减少SCF迭代次数。
最后，我们将FDVIB结果与ph.x的nat_todo=2结果对照，并进一步导出振动模式文件和Shermo输入文件，用于模式可视化和局域热化学后处理。
本文使用Quantum ESPRESSO-7.5和Shermo-2.6.2。

---

2.软件准备

FDVIB可从GitHub Releases下载Linux x86_64预编译版本。下载并解压后，将fdvib放入PATH即可使用。同时需确保Quantum ESPRESSO的pw.x和dynmat.x也能正常调用。
仓库：https://github.com/wxia529/fdvib
Releases：https://github.com/wxia529/fdvib/releases
文档：https://wxia529.github.io/fdvib/

FDVIB常用命令如下：

1. fdvib init local
   
2. fdvib init gas
   
3. fdvib -inp fdvib.in
   
4. fdvib modes fdvib/results
   
5. fdvib shm fdvib/results
   
6. fdvib thermo fdvib/results -inp thermo.in

复制代码

一个典型流程是：

1. 1. 准备并测试QE的scf.in
   
2. 2. 生成或编写fdvib.in
   
3. 3. 运行fdvib -inp fdvib.in
   
4. 4. 根据需要导出Molden、Shermo或热化学结果

复制代码

FDVIB没有需要手动设置的restart参数。如果计算中断，重新运行同一个命令即可。程序会自动跳过已经完成的任务，从第一个未完成或失败的任务继续。

---

3.输入文件

运行FDVIB前，算例目录中至少需要两个文件：

1. case_directory/
   
2. |-- scf.in
   
3. `-- fdvib.in

复制代码

其中，scf.in是原始未位移结构的QE输入文件，fdvib.in是FDVIB控制文件。thermo.in不是主计算必需文件，只在后续运行fdvib thermo时使用。
FDVIB不会修改用户顶层目录中的scf.in。它会根据原始scf.in生成参考SCF和所有位移任务的QE输入文件。例如本例选择第74、75号原子后，运行目录类似：

1. fdvib/calculations/
   
2. |-- init_scf_001/scf.in
   
3. |-- disp_0074_x_p_001/pw.in
   
4. |-- disp_0074_x_m_001/pw.in
   
5. |-- disp_0074_y_p_001/pw.in
   
6. |-- disp_0074_y_m_001/pw.in
   
7. |-- disp_0074_z_p_001/pw.in
   
8. |-- disp_0074_z_m_001/pw.in
   
9. |-- disp_0075_x_p_001/pw.in
   
10. |-- disp_0075_x_m_001/pw.in
    
11. |-- disp_0075_y_p_001/pw.in
    
12. |-- disp_0075_y_m_001/pw.in
    
13. |-- disp_0075_z_p_001/pw.in
    
14. `-- disp_0075_z_m_001/pw.in

复制代码

这些pw.in才是实际带有位移坐标的QE输入。用户原始的case_directory/scf.in保持不变。

---

4.QE输入文件scf.in

scf.in应先单独用pw.x测试，确保能够正常收敛并输出完整原子力。FDVIB用有限差分力构造Hessian，因此力的精度会直接影响频率。
本例scf.in的关键部分如下，实际计算时必须提供完整的75个原子坐标：

1. &CONTROL
   
2.   calculation = 'scf'
   
3.   etot_conv_thr = 1.0d-5
   
4.   forc_conv_thr = 1.0d-3
   
5.   outdir = './out/'
   
6.   prefix = 'lta'
   
7.   tprnfor = .true.
   
8.   pseudo_dir = '~/pseudo/dojo/nc-sr-05_pbe_standard_upf'
   
9.   verbosity = 'low'
   
10. /
    
11. &SYSTEM
    
12.   ecutwfc = 100.0d0
    
13.   ibrav = 0
    
14.   nat = 75
    
15.   nosym = .true.
    
16.   ntyp = 5
    
17.   vdw_corr = 'dft-d3'
    
18.   dftd3_version = 4
    
19.   dftd3_threebody = .false.
    
20.   occupations = 'fixed'
    
21.   nspin = 1
    
22. /
    
23. &ELECTRONS
    
24.   conv_thr = 1.0d-9
    
25.   electron_maxstep = 50
    
26.   mixing_beta = 0.4
    
27. /
    
28. 
    
29. ATOMIC_SPECIES
    
30. Al  26.982  Al.upf
    
31. C   12.011  C.upf
    
32. H    1.008  H.upf
    
33. O   15.999  O.upf
    
34. Si  28.086  Si.upf
    
35. 
    
36. ATOMIC_POSITIONS angstrom
    
37. !...前72个分子筛骨架原子...
    
38. H   4.3422304235  0.2660025317  4.4486388564
    
39. O   6.5168625726  1.2666840637  6.1417271027
    
40. C   5.6654265957  0.8046175659  5.5658904674
    
41. 
    
42. K_POINTS gamma
    
43. 
    
44. CELL_PARAMETERS angstrom
    
45. 11.9190000000   0.0000000000   0.0000000000
    
46. 0.0000000000  11.9190000000   0.0000000000
    
47. 0.0000000000   0.0000000000  11.9190000000
    

复制代码

这里重要的是:
第一，calculation='scf'。FDVIB做的是固定离子有限差分，位移结构不应继续做几何优化。
第二，必须设置tprnfor=.true.。否则QE不会输出完整原子力，FDVIB无法提取力常数。
第三，原始scf.in中不要设置startingpot='file'。这是因为FDVIB内部已经专门处理了电荷密度读取流程：程序会先运行一次未位移结构的参考SCF，然后把参考SCF得到的收敛电荷密度复制到每个位移任务中，再自动在位移任务里加入startingpot='file'。（因为位移结构和原始结构只差0.01Å，参考电荷密度通常已经非常接近新结构的自洽密度，读取它可以明显减少SCF迭代次数）
此外，FDVIB会在位移计算中自动插入或替换：

1. disk_io = 'nowf'
   

复制代码

因为有限差分只需要原子力，不需要保存波函数文件。参考SCF则保留用户原始输入中的disk_io设置。

---

5.FDVIB输入文件fdvib.in

第一次准备输入时，可以运行：

1. fdvib init local
   

复制代码

该命令会生成一个fdvib.in模板。本文的实际输入如下：

1. scf_input = scf.in
   
2. outdir = fdvib
   
3. system_type = local
   
4. selected_atoms = 74,75
   
5. displacement_angstrom = 0.01
   
6. pw_command = mpirun -np 64 pw.x
   
7. prefix = system
   
8. run_dynmat = true
   
9. dynmat_command = dynmat.x
   

复制代码

fdvib.in是普通的key=value文本文件，不是QE的namelist。
本例最关键的是这三行：

1. system_type = local
   
2. selected_atoms = 74,75
   
3. displacement_angstrom = 0.01
   

复制代码

它们表示FDVIB只对第74、75号原子进行正、负位移，位移量为0.01Å。两个原子共有6个笛卡尔坐标，因此会产生12个位移SCF任务。
其余73个原子不会被位移，但它们仍保留在完整周期结构中，并参与每一次SCF计算。因此，FDVIB不是把CO单独拿出来算，而是在完整CO@H-LTA周期环境中计算CO的局域振动。
还要注意，prefix=system是FDVIB结果文件前缀，不是QE的prefix。QE的prefix='lta'用于命名lta.save目录；FDVIB的prefix=system用于命名system.dynG、system.freq.out和system.shm等结果文件。

---

6.运行FDVIB
准备好scf.in和fdvib.in后，在算例目录中运行：

1. fdvib -inp fdvib.in > fdvib.out

复制代码

本例fdvib.out的主要内容如下：

1. Running unperturbed reference SCF
   
2. Running disp_0074_x_p
   
3. Running disp_0074_x_m
   
4. Running disp_0074_y_p
   
5. Running disp_0074_y_m
   
6. Running disp_0074_z_p
   
7. Running disp_0074_z_m
   
8. Running disp_0075_x_p
   
9. Running disp_0075_x_m
   
10. Running disp_0075_y_p
    
11. Running disp_0075_y_m
    
12. Running disp_0075_z_p
    
13. Running disp_0075_z_m
    
14. Completed 12, preserved 0 displacement jobs
    
15. Hessian max asymmetry: 5.474338e-04 Ry/Bohr^2
    
16. Wrote fdvib/results/system.dynG and fdvib/results/dynmat.in
    
17. Running dynmat.x
    
18. FDVIB calculation completed
    

复制代码

日志表明参考SCF和12个位移任务均正常完成。FDVIB会检查SCF收敛、JOB DONE和完整力输出，而不是简单地把文件存在视为任务成功。
如果任务中断，重新运行同一个命令即可。FDVIB会跳过已经完成的任务，并从第一个未完成或失败的任务继续。

---

7.输出文件
主计算完成后，主要结果位于：

1. fdvib/results/
   
2. ├── metadata.dat
   
3. ├── system.dynG
   
4. ├── dynmat.in
   
5. ├── dynmat.out
   
6. └── system.freq.out
   

复制代码

其中：

• system.dynG是FDVIB生成的dynmat.x能读取Γ点动力学矩阵；
• dynmat.in是自动生成的dynmat.x输入；
• dynmat.out是dynmat.x输出；
• system.freq.out包含频率和振动向量；
• metadata.dat记录电子能量、活动原子和模式筛选信息。
  

本例metadata.dat为：

1. program = qe
   
2. electronic_energy_hartree = -9.2030325431500000e+02
   
3. multiplicity = 1
   
4. mode_selection = local
   
5. selected_atoms = 74,75

复制代码

---

8.导出Molden模式文件

运行：

1. fdvib modes fdvib/results

复制代码

程序输出：

1. Wrote 6 Molden modes to fdvib/results/system.mol

复制代码

生成的system.mol保留完整CO@H-LTA结构，但只有第74、75号原子具有非零振动位移。未选中的环境原子仍显示在结构中
本例system.mol中的频率为：

1. [Molden Format]
   
2. [FREQ]
   
3. 47.77492500
   
4. 55.86191100
   
5. 143.35994200
   
6. 189.98604800
   
7. 226.43112400
   
8. 2180.89792400

复制代码

FDVIB不计算红外强度，因此生成的.mol文件不包含[INT]部分。
导出cp2k风格的molden文件的好处不言而喻：可以方便的使用sob老师的MfakeG生成一个高斯View能打开的out文件查看振动

---

9.与ph.x结果对照

为了对比FDVIB有限位移差分方法和QE的DFPT方法，本文使用相同结构和电子结构参数，通过ph.x的nat_todo=2只计算第74、75号原子。
ph.in输入如下：

1. &INPUTPH
   
2.   prefix = 'lta'
   
3.   outdir = './out'
   
4.   fildyn = 'lta.dynG'
   
5.   tr2_ph = 1.0d-16
   
6.   nat_todo = 2
   
7.   nogg = .true.
   
8.   dftd3_hess = 'lta.hess'
   
9.   reduce_io = .true.
   
10. /
    
11. 0.0 0.0 0.0
    
12. 74 75
    

复制代码

两条路线随后都使用以下dynmat.x设置：

1. &INPUT
   
2.   asr = 'no'
   
3.   remove_interaction_blocks = .true.
   
4. /
   

复制代码

下表列出局域Hessian的6个非零模式。负值表示虚频。

方法	6个局域频率/cm-1
ph.x,tr2_ph=10d-12	−350.71,−283.56,−253.76,111.11,169.82,2166.94
ph.x,tr2_ph=10d-14	−46.72,−34.66,112.27,184.52,221.60,2180.24
ph.x,tr2_ph=10d-16	48.82,57.25,145.73,190.22,226.92,2180.71
ph.x,tr2_ph=10d-18	48.59,56.73,143.78,190.71,227.34,2180.81
ph.x,tr2_ph=10d-20	48.88,56.79,143.73,190.67,227.35,2180.80
FDVIB,Γ点	47.77,55.86,143.36,189.99,226.43,2180.90
FDVIB,1×1×1	47.63,55.64,143.58,190.06,226.39,2180.90

可以看到，tr2_ph=10d-12和10d-14时出现的大虚频会随着tr2_ph收紧而消失，因此主要来自数值未充分收敛。到tr2_ph=10d-16时，6个局域频率均为正；继续收紧到10d-18和10d-20后，结果进一步稳定。
以ph.x,tr2_ph=10d-20为基准，频率差异如下：

方法	MAD/cm-1	RMSD/cm-1	最大绝对差/cm-1
ph.x,tr2_ph=10d-16	0.580	0.876	2.001
ph.x,tr2_ph=10d-18	0.079	0.125	0.291
FDVIB,Γ点	0.683	0.768	1.110
FDVIB,1×1×1	0.703	0.839	1.258

这里的MAD和RMSD只描述不同计算结果之间的数值差距（实际上DFPT并没有进一步测试）。FDVIB和DFPT采用不同方法，因此这些数值不能简单理解为“哪一种方法绝对更准确”。

---

10.计算时间对比
所有测试都使用同一节点上的64个MPI进程。时间对比采用6个频率均为正的tr2_ph=10d-16。

方法	时间
FDVIB,K_POINTS gamma,总时间	3分10秒
FDVIB,1×1×1,总时间	6分13秒
DFPT,tr2_ph=10d-16,总时间	14分24秒

DFPT各阶段时间为：

DFPT阶段	时间
SCF	40秒
D3Hessian	1分29秒
ph.x,tr2_ph=10d-16	12分13秒
dynmat.x	2秒
合计	14分24秒

按总时间比较，本例中FDVIB Γ点约快4.5倍，FDVIB1×1×1约快2.3倍。
这些时间只代表本文硬件、并行方式和模型，不能视为所有体系上的固定加速比。实际成本还取决于能带数、赝势、SCF收敛、活动原子数量和并行效率。

---

11.导出Shermo输入

本文建议使用Shermo进行热化学后处理。FDVIB可以直接导出Shermo的.shm文件：

1. fdvib shm fdvib/results
   

复制代码

本例输出为：

1. Wrote fdvib/results/system.shm
   
2. SHM mode selection: local, retained 6, removed 219, largest removed |frequency|=0, smallest retained |frequency|=47.7749 cm^-1
   
3. Run Shermo with: Shermo fdvib/results/system.shm -imode 1 -PGlabel C1

复制代码

这里mode_selection=local表示只保留所选原子的局域模式。本例最终导出的正是6个局域频率：

1. *wavenum  //Wavenumbers (cm-1). Negative value means imaginary frequency
   
2. 47.7749250000
   
3. 55.8619110000
   
4. 143.3599420000
   
5. 189.9860480000
   
6. 226.4311240000
   
7. 2180.8979240000
   

复制代码

---

12.使用Shermo进行局域热化学后处理

对于吸附物局域振动，我们通常只需要振动贡献，不应把整个周期模型当成自由气相分子加入平动和转动。因此运行：

1. Shermo fdvib/results/system.shm -imode 1 -PGlabel C1
   

复制代码

其中，-imode 1表示忽略平动和转动，只计算振动及电子贡献；-PGlabel C1用于指定点群。
本文使用Shermo-2.6.2，温度为298.15K，振动频率缩放因子为1.0，并采用100cm-1阈值的Grimme低频熵插值。Shermo识别到6个频率：

1. There are 6 frequencies (cm^-1):
   
2.    47.8    55.9   143.4   190.0   226.4  2180.9
   
3. 
   
4. Translation contribution is ignored since imode=1
   
5. Rotation contribution is ignored since imode=1
   

复制代码

节选结果输出如下：

1.                            ===========================
   
2.                            ========== Total ==========
   
3.                            ===========================
   
4. Total q(V=0):       1.031806E+002
   
5. Total q(bot):       1.079181E-001
   
6. Total q(V=0)/NA:    1.713353E-022
   
7. Total q(bot)/NA:    1.792021E-025
   
8. Total CV:      39.848 J/mol/K       9.524 cal/mol/K
   
9. Total CP:      39.848 J/mol/K       9.524 cal/mol/K
   
10. Total S:       58.737 J/mol/K      14.038 cal/mol/K    -TS:    -4.186 kcal/mol
    
11. Zero point energy (ZPE):     17.013 kJ/mol      4.066 kcal/mol   0.006480 a.u.
    
12. Thermal correction to U:     25.977 kJ/mol      6.209 kcal/mol   0.009894 a.u.
    
13. Thermal correction to H:     25.977 kJ/mol      6.209 kcal/mol   0.009894 a.u.
    
14. Thermal correction to G:      8.465 kJ/mol      2.023 kcal/mol   0.003224 a.u.
    
15. Electronic energy:       -920.3032543 a.u.
    
16. Sum of electronic energy and ZPE, namely U/H/G at 0 K:       -920.2967745 a.u.
    
17. Sum of electronic energy and thermal correction to U:        -920.2933601 a.u.
    
18. Sum of electronic energy and thermal correction to H:        -920.2933601 a.u.
    
19. Sum of electronic energy and thermal correction to G:        -920.3000302 a.u.
    

复制代码

记录对比不同频率来源得到的热化学校正如下，单位均为kcal/mol。

方法	ZPE	U	H	G
ph.x,tr2_ph=10d-12	3.499	4.330	4.330	2.697
ph.x,tr2_ph=10d-14	3.858	5.003	5.003	2.835
ph.x,tr2_ph=10d-16	4.074	6.210	6.210	2.044
ph.x,tr2_ph=10d-18	4.071	6.210	6.210	2.037
ph.x,tr2_ph=10d-20	4.072	6.210	6.210	2.038
FDVIB,Γ点	4.066	6.209	6.209	2.023
FDVIB,1×1×1	4.066	6.209	6.209	2.022

tr2_ph=10d-12含3个虚频，tr2_ph=10d-14含2个虚频。在对应热化学计算中，这些虚频被忽略，因此这两行只用于展示频率未收敛对热化学结果的影响。其余5组结果均使用全部6个正频率。
由于-imode 1不计算平动和转动，表中的U与H相同。
计算吸附自由能时，应对吸附态、气相分子和参照态分别使用相适应的热化学处理。

---

13.气相CO分子的计算

FDVIB也支持孤立分子振动分析。对于孤立分子，应使用system_type=gas。此时FDVIB会位移所有原子，构造完整分子Hessian，并在后处理中去除平动和转动模式。
气相CO的fdvib.in如下：

1. scf_input = scf.in
   
2. outdir = fdvib
   
3. system_type = gas
   
4. selected_atoms = all
   
5. displacement_angstrom = 0.01
   
6. multiplicity = 1
   
7. pw_command = mpirun -np 64 pw.x
   
8. prefix = molecule
   
9. run_dynmat = true
   
10. dynmat_command = dynmat.x
    

复制代码

计算得到的6个原始频率为：

模式	频率/cm-1
1	−6.7468
2	−5.8031
3	−0.1468
4	21.7268
5	22.0869
6	2128.0609

前5个是数值残留的平动和转动模式，最后1个是CO伸缩振动。CO是线性分子，因此FDVIB会去除5个平动和转动模式，只保留1个内部振动。(气相计算中，FDVIB先得到完整的3N个分子模式，再根据分子类型扣除刚体平动和转动自由度：单原子无振动，线性分子保留 3N−5个振动自由度，非线性分子保留3N−6个振动自由度。)

运行：

1. fdvib shm fdvib/results
   

复制代码

得到：

1. Wrote fdvib/results/molecule.shm
   
2. SHM mode selection: linear, retained 1, removed 5, largest removed |frequency|=22.0869, smallest retained |frequency|=2128.06 cm^-1
   
3. Run Shermo with: Shermo fdvib/results/molecule.shm

复制代码

气相分子需要包含平动和转动贡献，因此这里直接运行Shermo：

1. Shermo fdvib/results/molecule.shm
   

复制代码

Shermo将CO识别为线性分子，并读取1个内部振动。298.15K、1atm下的结果为：

ZPE/(kcal/mol)	U/(kcal/mol)	H/(kcal/mol)	G/(kcal/mol)
3.042	4.524	5.116	−8.963

---

14.实际使用建议

第一，先保证QE力的精度。有限差分力常数对力噪声很敏感，正式计算前应测试ecutwfc、K点、SCF阈值、赝势和位移量。
第二，关注Hessian max asymmetry。如果该数值明显偏大，应优先检查SCF收敛、位移量和活动区域，而不是直接使用频率结果。
第三，低频处理必须保持一致。比较不同结构时，应统一温度、频率缩放因子和低频处理方法。
第四，局域振动不是完整固体热力学。FDVIB适合比较相同或相近骨架背景下的局域吸附构型，但不适合代替完整晶格振动自由能。
第五，新体系最好与ph.x或其他可靠方法做交叉对照，确认位移量和局域模式选择没有引入明显偏差。
第六，如果您在实际工作中使用了Shermo，请按照http://sobereva.com/soft/shermo中的要求引用参考文献

---

15.结语

在CO@H-LTA测试中，FDVIB与ph.x对相同两个原子得到的局域振动频率非常接近。相对于ph.x,tr2_ph=10d-20结果，FDVIBΓ点结果的最大频率差约为1.11cm-1。相对于采用tr2_ph=10d-16的DFPT流程，本文算例中FDVIBΓ点总时间由14分24秒降至3分10秒。
FDVIB的优势在于流程简单、局域目标明确，并且可以直接生成Molden模式文件和Shermo输入文件。对于吸附物局域振动和局域热化学校正，它是一条实用的计算路线。
但FDVIB结果应始终理解为局域谐振近似，而不是完整周期声子谱。对于完整晶体振动性质和完整固体热力学，仍需要使用完整声子计算。


更详细的解释请参考文档：https://wxia529.github.io/fdvib


[1] Paolo Giannozzi et al., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter, 21, 395502 (2009) DOI: 10.1088/0953-8984/21/39/395502
[2] Paolo Giannozzi et al., Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter, 29, 465901 (2017) DOI: 10.1088/1361-648X/aa8f79
[3] Tian Lu, Qinxue Chen, Shermo: A general code for calculating molecular thermodynamic properties, Comput. Theor. Chem., 1200, 113249 (2021) DOI: 10.1016/j.comptc.2021.113249