GaN slab模型结构优化后表面原子弛豫问题

noodles的困惑

最近在做GaN0001面吸附计算，建立slab模型，吸附前先进行clean surface结构优化。发现结构优化后2x2 GaN（0001）面最上层四个Ga原子：两个Ga原子上移两个Ga原子下移，具体的结果图在下面。
slab模型具体信息：四层GaN双原子层，表面是0001面（Ga面）；结构优化时最底层N原子用0.75H饱和，下两层GaN和赝氢固定，上面两层GaN弛豫；真空层15埃。

各位老师，我的问题是：
1、首先，在计算GaN0001面吸附单原子前，clean surface的slab模型的结构优化是必要的吗？

• 我的理解：吸附前是肯定要进行结构优化的——固定晶胞，固定下面两层原子，弛豫上面两层，slab模型结构优化后能量最低，应力释放，原子到达平衡位置，slab模型能更好的与GaN实际表面吻合，然后再进行后续的单原子吸附计算。
  

2、为什么优化后会出现两头原子上移，中间原子下移的结构呢？

• 我的理解：这和VASP结构优化算法有关——离子弛豫相关参数：IBRION、POTIM、EDIFFG，在结构优化时，第一步尝试的太大，离子偏离的太大，然后就沿着这个节奏一直优化下去，最终导致了最外层Ga原子z方向的高低，不过这种构型也能达到设置的收敛条件，所以计算收敛。我的理解比较笼统并且还有一些错误的地方，望各位老师指正。
  

3、这个计算是收敛了的，但是如何判断合理性呢？从力/能量/物理/实际表面的角度判断？对优化后（计算上都收敛）两种情况（1）最上层Ga原子整体上下移动，即表面平整的（2）最上层Ga原子z方向上会有高低，即表面不平整，我们该选择哪种结构进行后续的单原子吸附呢？

• 我的理解：计算上收敛但是是否符合物理上/实际表面的情况还要具体分析（但是如何分析我比较困惑）。关于两种结构选择，我认为是都可以，但是要结合自己研究的问题来判断选择哪种，我倾向于选择优化后平整的最外层Ga原子的模型。
  

下面是计算时INCAR、KPOINTS、POSCAR以及优化前后结构图

noodles的困惑

@卡开发发 请教一下卡老师，打扰了~

卡开发发

1、一般建议优化而已，也有因为吸附而偏离clean情况的构型，自然是手动调整的越接近实际情况所需要优化的步骤越少越好。通常势能面不大复杂的情况，先优化后优化可能得到的结果基本上是差不多的（数值误差范围内）。固定这问题见仁见智，尤其是固定的层数可多可少，如果因为结构稍微跑偏就使用固定那么我会怀疑计算的意义，还是说使用更充分的层高就能避免这样的问题？或是说需要更近一步的模型修正才能保证这种固定使用能复现足够层高的结果，我不好说这个问题。

2、导致这个问题的原因很多，算法只是一个环节，你说的构型优化算法参数选择确实是个问题，所以通常不建议一上来就IBRION=1，而是用IBRION=2，这样有线搜索的过程来保证步长的合理性，此时POTIM影响的是步长的初猜。类似的体系我可能恰好做过，我猜测是你的体系实际有磁性，可能有磁矩的原子在表面，但具体是上下哪个面你可以自己测试下，磁矩初猜通过magmom来设置，要开ispin=2的选项，你可以试试看这个方案，我不确定是不是这个原因，仅作为一个排查方向。

3、确实看研究问题而定（1）有些情况是针对稳定表面，可以看表面能等作为参考，这是从统计的角度看（2）从研究反应的角度看，稳定面不一定是活性面，真正能起到催化作用的根本不是稳定面（3）实验上通过某些技术观察、佐证了某些晶面。

noodles的困惑

卡开发发 发表于 2021-11-14 14:03
1、一般建议优化而已，也有因为吸附而偏离clean情况的构型，自然是手动调整的越接近实际情况所需要优化的步 ...

晚上好！感谢老师回复！由于我课题的主要部分就是吸附，无论是clean surface还是吸附后表面，结构优化非常重要，如果糊里糊涂过去心里难过这个坎，总有个疙瘩。。我先消化一下这几个问题还有老师您的建议与见解，后续会给您我的反馈！

卡开发发

noodles的困惑 发表于 2021-11-14 20:32
晚上好！感谢老师回复！由于我课题的主要部分就是吸附，无论是clean surface还是吸附后表面，结构优化非 ...

可以，虽然我不能保证我的观点完全正确，但是欢迎讨论，细节弄明白有助于大家对模型有更深刻的理解。

noodles的困惑

卡开发发 发表于 2021-11-14 21:22
可以，虽然我不能保证我的观点完全正确，但是欢迎讨论，细节弄明白有助于大家对模型有更深刻的理解。

老师，我想分析slab模型结构优化后（reached required accuracy）原子受力情况，为什么在OUTCAR尾部没有找到TOTAL-FORCE部分呀？只有“FORCE on cell =-STRESS in cart. coord.  units (eV):”这一栏，是因为slab模型OUTCAR中不用TOTAL-FORCE来表明力的收敛情况，而用“FORCE on cell”来显示，还是因为我的计算是错误的呀？

卡开发发

noodles的困惑 发表于 2021-11-16 16:03
老师，我想分析slab模型结构优化后（reached required accuracy）原子受力情况，为什么在OUTCAR尾部没有 ...

应该有，你可以再找找，实在不行把OUTCAR发上来看看。

noodles的困惑

卡开发发 发表于 2021-11-16 19:16
应该有，你可以再找找，实在不行把OUTCAR发上来看看。

卡卡老师，晚上好！感谢回复！关于“OUTCAR末尾找不到TOTAL-FORCE”已经解决啦，原因是我在INCAR中NWRITE设置成0了，现在用改为2或者直接不写（默认值也为2）就好啦——这真的是低级错误，汗颜。。。

为了确定NWRITE每个参数对OUTCAR中写入内容的作用，我做了个测试，其他参数和输入文件完全相同情况下，取不同NWRITE的值，出现几个问题：
1、vasp手册上说，设置NWRITE=0时，关于force这一项，会在OUTCAR中第一个和最后一个离子步写入，但是我设置NWRITE=0后，得到的OUTCAR中只有第一个离子步后有TOTAL-FORCE，最后一个离子步后面没有。这是为什么呢？参数为1、2、3的测试OUTCAR中显示内容和手册中相同。
2、除了NWRITE参数不同外，其他参数和输入文件完全相同，但是output末尾显示的内容不一样，这是为什么呢？具体不同如下图，NWRITE=0和3的output相同，1和2的output相同。3、对于图中output中的 ZBRENT: fatal error in bracketing please rerun with smaller EDIFF, or copy CONTCAR to POSCAR and continue“是怎么产生的，是如何解决的呢？我再查找这个问题时，这个帖子有讨论，我想听一听您对这个问题的看法或者对这个帖子里的说法的理解。
http://mu****g.com/t-11528282-1

问题有点多，谢谢老师！

Epikur

我觉得可以先完全弛豫整个slab，再固定做吸附能的计算
有几点可以探讨：
第一是对称性是否关掉了
第二是IDIPOLE对能量的影响很小，一些情况会导致电子步不收敛
第三个Slab的厚度是不是需要考虑，有的时候slab如果薄了，会发生收缩，审稿人有的也会问你的slab厚度对吸附能是否有影响

卡开发发

noodles的困惑 发表于 2021-11-19 17:31
卡卡老师，晚上好！感谢回复！关于“OUTCAR末尾找不到TOTAL-FORCE”已经解决啦，原因是我在INCAR中NWRITE ...

1、这问题我没啥印象，如果非要弄清楚可以看下源码的逻辑才知道，我印象中最后一步应该有。如果还不行，我建议你查看下vasprun.xml里面信息记录的情况，这样应该会保险一些，如果这里没有可能有必要一块仔细看下源码。
2、你确定参数完全一样？两个收敛的步数都不一样啊。
3、ZBRENT的提示来源于vasp的CG算法中的线搜索过程，使用的是修改的Brent算法，然后Brent算法找极小出问题给出的提示。这个原因是，Brent算法先要划定一个导数相反的区间（保证区间内有极小值，也就是所谓的bracketing），然后再进行区间内逼近找极小，如果Brent过程中一直找不到这个区间，就会一直延伸，直至程序达到内置的阈值而报错。当势能面噪声较大或者初始结构极差的情况，就有可能发生这种错误，所以要避免的也是结构太差精度不足这类情况，当然造成这个的原因也远不止EDIFF这一个参数的问题。

关于线搜索即Brent具体的一种实现可以看Numerical Receipts，也比较简单，就不详细说了。

noodles的困惑

Epikur 发表于 2021-11-19 17:55
我觉得可以先完全弛豫整个slab，再固定做吸附能的计算
有几点可以探讨：
第一是对称性是否关掉了

感谢回复！
关于slab模型弛豫，我看包括大师兄和VASP hand on session：surface的教程中是：固定晶胞，固定底层原子模拟体相，弛豫上面几层原子模拟表面来进行结构优化的，然后在优化后的表面上的位点上放一个原子进行吸附结构的优化。当然我也想过您的这种方式，不过我觉得”固定底层原子，弛豫上层原子”的做法比较合理一些，而且弛豫前后的能量就是“表面能”吧。
关于您的探讨：
第一：我用的是默认的对称性设置
第二：我记得电子步的默认值是60，检查了几个计算，并没有出现电子步没收敛的情况。不过关于电子步收敛我有个疑问，如果前几个离子步的电子步没有收敛但是最后达到了收敛条件，那么计算结果是可靠的吗？还是说在结构优化时，第一个离子步里的电子步就必须收敛？
第三：slab厚度问题，严格来说是要做测试，不过我目前计算是直接取了15埃，另外在结构优化时ISIF=2，应该不会出现真空层收缩的情况吧？

这是我的思考与一些疑惑，基础薄弱，还望老师多多指教多多交流~谢谢！

Epikur

noodles的困惑 发表于 2021-11-19 20:22
感谢回复！
关于slab模型弛豫，我看包括大师兄和VASP hand on session：surface的教程中是：固定晶胞， ...

互相学习哈
一般来说SCF的精度决定了力的精度，早期的电子步不收敛对最终结果的影响还真没考虑过，期待高手解惑
如果你的slab已经15 A了算是相当厚了，应该不用做啥测试了

noodles的困惑

Epikur 发表于 2021-11-19 21:10
互相学习哈
一般来说SCF的精度决定了力的精度，早期的电子步不收敛对最终结果的影响还真没考虑过，期待 ...

嗯嗯好的。
请问slab模型结构优化后判断力的收敛时，vasp程序是判断TOTAL-FORCE中所有原子的受力小于INCAR中设定的EDIFFG值停止计算还是未固定的原子受力小于EDIFFG时停止计算呀？

Epikur

noodles的困惑 发表于 2021-11-20 15:17
嗯嗯好的。
请问slab模型结构优化后判断力的收敛时，vasp程序是判断TOTAL-FORCE中所有原子的受力小于INC ...

未固定的原子，固定的原子的受力就不管了

noodles的困惑

卡开发发 发表于 2021-11-19 20:18
1、这问题我没啥印象，如果非要弄清楚可以看下源码的逻辑才知道，我印象中最后一步应该有。如果还不行， ...

感谢回复！
关于第二点，我可以确定，除了测试的NWRITE，其他的输入文件以及参数设置完全相同，这个问题我还会继续排查。

新的问题：
1、在slab模型结构优化时，固定了底层原子模拟bulk，上面几层原子可以弛豫模拟表面，力的收敛条件设置为-0.03，那么请问发发老师，VASP是当所有原子上的力都小于|EDIFFG|时还是未固定原子的力小于|EDIFFG|时认为力达到收敛了呢？
我的理解：我实际的计算结果是，固定原子未达到力收敛条件，未固定原子全部达到收敛条件，total-drift z方向未达到力收敛条件，所以我判断未固定原子受力都小于|EDIFFG|时VASP认为达到力的收敛条件。

2、在我判断slab模型结构优化计算的能量（EDIFF）、力（EDIFFG）均达到收敛条件后，有进行了一次静态计算（NSW=0），这两次计算结果OUTCAR中的未固定原子受力有不同，这是为什么呀？z方向上有个原子前后两次计算受力：（0.01416-0.009238）/0.009238=12%，这个差异在可以接受范围之内嘛？或者说这么一个比较的逻辑有意义吗？

感谢老师解答！