计算化学公社

标题: DFT+U 的 CeO2 表面催化体系做频率计算时，可动原子区域如何选择？ [打印本页]

作者
Author: zhj45318403 时间: 3 day ago
标题: DFT+U 的 CeO2 表面催化体系做频率计算时，可动原子区域如何选择？
各位老师、同学好，我想请教一个关于 VASP 表面催化体系频率计算中 partial Hessian 区域选择的问题。

我的体系是 Pd/CeO2 催化甲烷 C-H 活化相关模型，包含表面氧空位和次表层氧空位。体系总原子数约 140-150 个。前期已经完成了 IS 和 FS 的结构优化，现在想做频率计算，主要目的是确认 IS/FS 是否为局域极小态，同时获得与吸附物和表面反应区域相关的振动信息。

由于体系中含 CeO2，计算中考虑了 Ce 4f 电子的强关联效应，采用 DFT+U 处理 Ce 4f。主要设置如下：

LDAU    = .TRUE.
LDAUTYPE = 2
LDAUL = -1  3 -1 -1 -1
LDAUU = 0.0 5.0 0.0 0.0 0.0
LDAUJ = 0.0 0.0 0.0 0.0 0.0
LMAXMIX  = 6
LASPH = .TRUE.
ISPIN = 2

元素顺序对应  O Ce PD C H，即只对 Ce 的 4f 轨道加 U，Ueff = 5.0 eV。结构优化阶段使用了相同的 DFT+U 设置，并根据体系情况设置过 MAGMOM；频率计算中继承结构优化 INCAR 中的 MAGMOM 设置，不重新自动生成。由于涉及表面/次表层氧空位，部分 Ce 可能存在 Ce3+/Ce4+ 局域化和磁矩分布差异。

频率计算设置大致如下：

IBRION = 5
NSW = 1
NFREE  = 2
POTIM  = 0.015
EDIFF  = 1E-06
ISYM = 0

我使用 POSCAR 的 Selective dynamics 做 partial Hessian，只放开部分原子，其余 slab 原子固定。根据 VASP Wiki，IBRION=5 下 Selective dynamics 会只计算标记为 T 的 Hessian 分量；NFREE=2 时每个可动原子大约需要 6 个位移计算。因此可动原子数对计算量影响很大。

目前我的处理方式是按 z 方向分层，只放开最上方 5 层原子。用类似 VASPKIT 402 的分层方式判断，阈值 0.02 时整个 slab 大约可以识别出 8-12 层。放开 top 5 层后，可动原子数大致为：

- IS：约 34-67 个可动原子；
- FS：约 53-99 个可动原子。

这个区域通常包括 CH4 或 CH3+H、Pd、表面 O/Ce，以及最靠近吸附物的表面反应区域。我的反应构型中，不论氧空位位于表层还是次表层，C-H 活化时 C-H 基本都是对准表面 O，因此我初步认为主要振动自由度应集中在吸附物、Pd 和表面 O 附近。

但我担心的问题是：有些模型涉及次表层氧空位。如果只放开 top 5 层，实际可能主要包括吸附物、Pd、表面 O/Ce 等区域，而没有充分放开次表层氧空位附近的 O/Ce。如果再增加一层 O/Ce，可动原子数会显著增加，频率计算量也会大幅增加。

因此想请教大家：

1. 对这种 DFT+U 的 Pd/CeO2 表面甲烷活化体系，用 partial Hessian 只放开吸附物 + Pd + 表面反应区域 + 表面几层原子，是否可以作为判断 IS/FS 局域极小态的合理近似？

2. 如果存在次表层氧空位，是否有必要把氧空位附近的 Ce/O 也加入可动区域？还是说只要结构优化时这些原子已经充分弛豫，并且频率计算保持相同的 DFT+U 与 MAGMOM 设置，主要关注吸附物和表面反应 O/Pd 区域通常就可以接受？

3. 从经验上看，对于表面反应频率验证，是更推荐：
- 按 z 方向统一放开 top N 层；
- 还是手动选择“吸附物 + Pd + 成键表面 O/Ce + 氧空位近邻 Ce/O”作为可动区域？

我目前的想法是先用 top 5 层（只包含CH4和上面一层O和一层Ce批量计算，后续如果某些含次表层氧空位模型出现虚频，或者关键构型比较敏感，再扩大可动区域复算。想听听大家是否认为这个策略合理，或者是否有更规范的处理方式。

感谢各位指教！

作者
Author: nuonuo123456 时间: 3 day ago
1. 个人认为频率分析即便是找极小态放开五层也太多了，你看文献几乎做甲烷吸附很少有找极小态的吧
2. 感觉不太有必要
3. 手动选择毕竟好
另外，如果是甲烷吸附活化你肯定是以表面反应为主，我认为可以结构优化之后直接找过渡态，因为你过渡态CLNEB的过程中如果初末态不是极小态也可以再去重新优化。如果每一个吸附物种都这样做频率，付出和收获完全不成正比。

作者
Author: zhj45318403 时间: 3 day ago

nuonuo123456 发表于 2026-5-18 20:38
1. 个人认为频率分析即便是找极小态放开五层也太多了，你看文献几乎做甲烷吸附很少有找极小态的吧
2. 感觉 ...

谢谢老师回复！我理解您的意思是：其实结构优化后振动分析没有太大的必要，如果CI-NEB能跑通说明初、末态就是局域极小，如果跑不通就再优化就是了，没有单独去找极小态的必要。

仔细想来阅读文献确实也没有看到谁专门去找极小态的，谢谢您的指导！

作者
Author: nuonuo123456 时间: 3 day ago

zhj45318403 发表于 2026-5-18 20:57
谢谢老师回复！我理解您的意思是：其实结构优化后振动分析没有太大的必要，如果CI-NEB能跑通说明初、末态 ...

对的，我就是这个意思，而且可以clneb和dimer联用一般效率会更高。而且做neb的过程中要仔细观察结构，有时候它可能第一个点的能量比初态还低但是结构变化很小，或者是再次优化完能量或者结构差距很小这时候基本上可以默认是极小态，要灵活运用，多看结构。

作者
Author: zhj45318403 时间: 前天 14:31
本帖最后由 zhj45318403 于 2026-5-19 14:33 编辑

nuonuo123456 发表于 2026-5-18 21:15
对的，我就是这个意思，而且可以clneb和dimer联用一般效率会更高。而且做neb的过程中要仔细观察结构，有 ...

十分感谢老师。如果可以的话还想咨询一个与CI-NEB有关的问题。我目前计算200多步始终没有收敛到我设定的值，并且有越来越跑飞的趋势。我看到您说的和dimer联用，其它帖子也多次见到这个方法，我想额外的问题是：
1、是否可以罗列出所有轮次的力收敛结果，找到切线力接近于0 能量最高最大力比较小的几个点，拿出来做dimer，而不用等全部ci-neb跑完
2、以下是我进行CI-NEB其中一个模型的INCAR ，如您有时间，能否帮我见到看看我一下关于CI-NEB的参数的设置是否合理？十分感激！

# --- Global Parameters ---
SYSTEM = Pd
ISTART = 0
ISPIN  = 2
PREC = Normal
ENCUT  = 450
LWAVE  = .FALSE.
LCHARG = .TRUE.
LAECHG = .FALSE.
ICHARG = 2
LREAL  = Auto
IVDW = 11
ISYM = 0 # 取消对称性

# --- Electronic Relaxation ---
ALGO = Normal
ISMEAR = 0
SIGMA  = 0.05
NELM = 150
EDIFF  = 1E-06

# --- Ionic Relaxation ---
NSW = 200
IBRION = 3
POTIM = 0
ISIF = 2
EDIFFG = -0.05

# VTST Transitional state search
IOPT = 1
ICHAIN = 0
LCLIMB = .TRUE.
IMAGES = 5
SPRING = -5

# --- DFT+U Settings (The Soul of CeO2) ---
LDAU    = .TRUE.
LDAUTYPE  = 2
LDAUL    = -1  3 -1 -1 -1
LDAUU    = 0.0 5.0 0.0 0.0 0.0
LDAUJ    = 0.0 0.0 0.0 0.0 0.0
LMAXMIX = 6
LASPH  = .TRUE.
MAGMOM = 94*0.0 27*0.0 1.2 4*0.0 1.2 8*0.0 1.2 2*0.0 1.2 2*0.0 0.0 0.0 4*0.0

# --- Parallel Performance ---
KPAR = 4
NCORE  = 8

作者
Author: nuonuo123456 时间: 前天 17:57

zhj45318403 发表于 2026-5-19 14:31
十分感谢老师。如果可以的话还想咨询一个与CI-NEB有关的问题。我目前计算200多步始终没有收敛到我设定的 ...

1. 你可以自己在插点的文件夹里面看，也可以自己写脚本，vtst自带的有没有我没具体了解过。如果是把所有插点的能量一起罗列vtst有自带的脚本可以实现。能不能等dimer跑完要看你所有插点的收敛情况，比如如果是-0.1的收敛限，最高点到了-0.1，其他点0.2，0.3的样子也可以停掉，如果是-0.5的收敛限我建议等他跑完。如果跑飞了那就是没找到合适的路径，如果过渡态跑不出来大概率是初末态的机构有问题，或者是精度不够。
2. 参数我没看出来有什么明显的问题，EDIFF可以调到-7试试，不一定能跑出来。

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