计算化学公社

标题: 用CP2K做AlN(001)晶体表面吸附N2分子，表面发生形变怎么解决？ [打印本页]

作者
Author: THY 时间: 2025-7-17 22:55
标题: 用CP2K做AlN(001)晶体表面吸附N2分子，表面发生形变怎么解决？
本帖最后由 THY 于 2025-7-27 16:38 编辑

各位老师，我用AlN(001)面吸附N2分子，但是优化后发现晶体表面的某些原子的位置发生了形变，

如果说是与N2接触的表面原子发生形变可以理解，但是一些远离吸附位点的原子，如下图画圈的原子，同样发生了变化

请问这是为什么，而且奇怪的是，远离吸附位点的原子发生形变的概率是随机的，也就是说每次吸附都会有不同位子的原子发生形变（这里指远离吸附位点的表面原子）

                  优化前
(, 下载次数 Times of downloads: 1)
               优化后
(, 下载次数 Times of downloads: 1)

本人附上inp文件和out文件，inp文件中没有添加k点，但是本人实测添加k点也没法解决上诉问题

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

原胞文件：
(, 下载次数 Times of downloads: 3)

作者
Author: sobereva 时间: 2025-7-18 07:25
XY方向扩成更大晶胞再看情况。过小的晶胞往往没法正确展现表面化学吸附导致的附近的表面结构变化
当前用smearing完全多余

作者
Author: THY 时间: 2025-7-18 14:46

sobereva 发表于 2025-7-18 07:25
XY方向扩成更大晶胞再看情况。过小的晶胞往往没法正确展现表面化学吸附导致的附近的表面结构变化
当前用sm ...

sob老师，您说的这个办法是可行的，但是我之前是4x4的超胞，一共96个原子，现在改成5x5，一共150个原子，计算时间直接翻了三倍左右，之前吸附N2的计算时间大概1个小时（64核），现在要3个小时左右，时间成本上代价巨大，请问还有其他办法吗？

作者
Author: Stardust0831 时间: 2025-7-25 02:05

THY 发表于 2025-7-18 14:46
sob老师，您说的这个办法是可行的，但是我之前是4x4的超胞，一共96个原子，现在改成5x5，一共150个原子， ...

试试看对原胞结构优化，获得更合理的晶格常数。
我见过由于cif不可靠导致直接拿来切表面的时候表面会变形的情况。

作者
Author: THY 时间: 2025-7-25 09:04

Stardust0831 发表于 2025-7-25 02:05
试试看对原胞结构优化，获得更合理的晶格常数。
我见过由于cif不可靠导致直接拿来切表面的时候表面会变 ...

老师您好，试过了，好像就是不行，原胞切面，结构优化没问题，但是加上N2后就出现结构变形的问题了

作者
Author: Stardust0831 时间: 2025-7-25 12:34

THY 发表于 2025-7-25 09:04
老师您好，试过了，好像就是不行，原胞切面，结构优化没问题，但是加上N2后就出现结构变形的问题了

给出原胞文件，我来试试看

作者
Author: sobereva 时间: 2025-7-25 15:00

THY 发表于 2025-7-18 14:46
sob老师，您说的这个办法是可行的，但是我之前是4x4的超胞，一共96个原子，现在改成5x5，一共150个原子， ...

150个原子对于CP2K来说根本不大
而且3个小时优化出一个结构根本不算耗时过高

作者
Author: THY 时间: 2025-7-27 16:39

Stardust0831 发表于 2025-7-25 12:34
给出原胞文件，我来试试看

老师您好，原胞文件上传至1楼了

作者
Author: Stardust0831 时间: 2025-7-28 00:27
本帖最后由 Stardust0831 于 2025-7-28 00:31 编辑

@THY 我花了点心思摸索了一下，首先我解释一下为什么会出现严重的表面变形，原帖里的图其实是把下面几层全冻住了只留了最上面一层，但是如果直接尝试弛豫这个slab，会发现出现了比较大的变形，N和Al变得更平：
原始文件： (, 下载次数 Times of downloads: 2)
弛豫前： (, 下载次数 Times of downloads: 1) 弛豫后： (, 下载次数 Times of downloads: 1) （冻结最底层，放开了其余3层）
这与表面变形的趋势也解释了为什么后来N比Al高。
我归因于slab可能太薄，导致整个体系都出现了这种只有表面才可能有的“变平”，所以加厚了slab的层数以及冻结的层数（冻结2层放开4层）
同时注意到这个体系用PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH直接算会遇到SCF难收敛问题。尝试提高CUTOFF到1000Ry也还是一样的情况（其实根据经验，应该是400Ry的时候就已经可以CUTOFF收敛了），根据《CP2K中遇到SCF难收敛时的解决方法》做了一些尝试，最后发现换成PBE-D3(BJ)/pob-TZVP-rev2以后能正常收敛。
在完成以上处理后，slab板的弛豫变得正常了：
(, 下载次数 Times of downloads: 1)
(, 下载次数 Times of downloads: 2)
最后，乘胜追击，加上N2，吸附以后结构基本符合经验预期：
(, 下载次数 Times of downloads: 4)
(, 下载次数 Times of downloads: 1)

最后，社长@sobereva 提到150个原子对于CP2K来说根本不大，我深表认同。这是常见情况，我天天用cp2k算上千原子的情况，计算成本完全可以接受。但是在本例中，由于PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH时跑OT和对角化都遇到了scf难收敛，导致每一步的scf收敛步数远比正常情况高的多（甚至可能最后也无法迭代到正确的/稳定的波函数），带来了耗时的不合理变长，最终导致了耗时其实难以接受。（哪怕最后扩胞+加厚，因为scf收敛步数很少，所以每步的平均耗时也少于了最开始那个很薄的情况）（我怀疑最开始模型收敛出来的波函数也未必是“稳定”的波函数，不过目前只听说wsr老师贡献了代码使pyscf可以用类似“TDA计算没有负激发能”的原理来实现周期性体系的波函数稳定性判断，而cp2k应该没法判断波函数稳定性，这一怀疑只好搁置）

作者
Author: THY 时间: 2025-7-28 09:23

Stardust0831 发表于 2025-7-28 00:27
@THY 我花了点心思摸索了一下，首先我解释一下为什么会出现严重的表面变形，原帖里的图其实是把下面几层全 ...

超级感谢老师，我认真学习一下

作者
Author: THY 时间: 2025-7-28 09:34
本帖最后由 THY 于 2025-7-28 11:07 编辑

Stardust0831 发表于 2025-7-28 00:27
@THY 我花了点心思摸索了一下，首先我解释一下为什么会出现严重的表面变形，原帖里的图其实是把下面几层全 ...

老师，另外我还想咨询一下啊，您计算这些模型，slab或者slab+N2，用的是多少核数？

另外，老师，我仔细看了您的模型，发现在吸附N2后，AlN的表面同样发生了较大的形变（这里指远离N2吸附位点的表面），我个人觉得这些形变才是真正导致几何优化时间过长的主要原因

而且这些形变的位置是随机的，也就是说对于同样一个inp文件，同时进行多次计算，其计算结果不是唯一

作者
Author: Stardust0831 时间: 2025-7-28 12:25
本帖最后由 Stardust0831 于 2025-7-28 12:28 编辑

THY 发表于 2025-7-28 09:34
老师，另外我还想咨询一下啊，您计算这些模型，slab或者slab+N2，用的是多少核数？

另外，老师，我 ...

计算时使用的配置你可以查看压缩包里的输出文件（.log文件），cp2k在开头就输出了。

形变导致的是几何优化需要的步数变长。我上文提到的是scf所需步数的不合理变长，scf不停震荡而导致耗时大幅增加。换言之，scf顺利收敛的情况下，哪怕花了50多步才结构优化收敛，也比一开始scf难收敛的时候优化10步的耗时短，更何况后者只用10步没法结构优化收敛、波函数也未必稳定。
”AlN的表面同样发生了较大的形变“很正常，表层原子受到的束缚更弱，表面点群被打破后出现这种程度的弛豫很合理。要是跑一段AIMD，表面很可能也会出现这种程度的变形。这个弛豫不应该波及到内层的原子，这也是我加厚slab的原因，加厚slab以后只有表层出现的较大形变，这是合理的。

表层原子都发生弛豫和"这些形变的位置是随机的"不存在必然联系。

作者
Author: 卡开发发 时间: 2025-7-28 13:19

Stardust0831 发表于 2025-7-28 00:27
@THY 我花了点心思摸索了一下，首先我解释一下为什么会出现严重的表面变形，原帖里的图其实是把下面几层全 ...

这问题其实http://bbs.keinsci.com/forum.php?mod=viewthread&tid=46368这里有讨论过。可能不需要上升到波函数稳定性的讨论，有一个不完全相同的视角，其实可以来看：沿着c轴取一列原子，即...Al-N-Al-N...，实际结构或是块体的情形，Al和N原子间距离是长短交替的，但算出来是距离是等长的，大家是否能联想到这个很像是论坛热门的体系C18？那么反之是否换一个EXX比例更高的杂化泛函或是合适的自作用修正方式，可能会导致可能并不需要那么厚的层高其实就能够呈现原子距离的长短交替？造成这种缺陷可能是GGA等局域/半局域泛函自身的问题？

作者
Author: Stardust0831 时间: 2025-7-28 13:38
本帖最后由 Stardust0831 于 2025-7-28 13:40 编辑

卡开发发发表于 2025-7-28 13:19
这问题其实http://bbs.keinsci.com/forum.php?mod=viewthread&tid=46368这里有讨论过。可能不需要上升到 ...

有道理的。这个材料是共价晶体，所以我确实没有第一时间想到可能有磁性。
根据老师您的摸索，您觉得这个材料的磁性应该是怎样的呢？

另外，DFT+U是否有希望实现正确的“长短键交替”呢？

作者
Author: 卡开发发 时间: 2025-7-28 13:43
本帖最后由卡开发发于 2025-7-28 13:54 编辑

Stardust0831 发表于 2025-7-28 13:38
有道理的。这个材料是共价晶体，所以我确实没有第一时间想到可能有磁性。
根据老师您的摸索，您觉得这个 ...

之前测试下来铁磁能量低，但是否能够反铁磁我不太确定，因为反铁磁的排列可能性要很多种。直觉上我认为这可能和Peierls机制相关，或许能存在局域磁矩也与此有关，DFT+U有可能不行，有可能DFT+U+V可以，有时间我可以考虑QE测试下看看。

作者
Author: get-it 时间: 2025-7-28 22:17

卡开发发发表于 2025-7-28 13:19
这问题其实http://bbs.keinsci.com/forum.php?mod=viewthread&tid=46368这里有讨论过。可能不需要上升到 ...

补充一下，很早之前的文章就研究过这个问题：https://link.aps.org/doi/10.1103/PhysRevLett.96.066102

作者
Author: 卡开发发 时间: 2025-7-29 03:13

get-it 发表于 2025-7-28 22:17
补充一下，很早之前的文章就研究过这个问题：https://link.aps.org/doi/10.1103/PhysRevLett.96.066102

我想表达的问题可能更类似于10.1103/PhysRevLett.111.036601这样。

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