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的困惑 发表于 2021-11-21 14:26
谢谢发发老师逐条解答！学到了很多！
关于“表面能”的讨论，我觉得问题在于两个方面：
1、表面能的定 ...

1、相同的体系表面积因为一样可以比较，要讨论的主要是一个块体劈开无论怎样都会形成两个表面，从这个角度说这两个表面都弛豫才是符合真实物理的，以足够的层高全部弛豫才是最合理的。反而去做固定是人为施加的，或者可以看出某些数学技巧，但这种方式并不见得能完整复现客观的约束。
2、9楼的意思应该并不是弛豫完整表面会与实际有较大偏离。我对书中的说法并不认同，我看法是这样：
（1）约束优化
把结构固定在体相的坐标不一定能保证这个部分的性质与体相一致，你固定的只是核坐标，从原理上说你其实只是能保证原子核所产生的核-核排斥与处于块体时的贡献是一致的，而电子的分布你没办法固定，因此还是会对这部分区域电子密度、静电势等产生扰动。如果要减少这部分的扰动，就得足够高的层高或者使用等效的电荷模型。非要结合固定底层原子的方法也是在上述基础上要保留额外足够高度被固定的原子，这样才是合理的，计算量上也未必有优势。
（2）对称表面
只有部分体系能保证切出的表面既存在对称性还能保证符合计量比。强行符合对称的情况而放弃计量比，处理不当不只是复现不了实际化学环境，甚至SCF迭代的稳定性都会出问题，所以这种情况我建议还是足够的真空层高并对静电势计算上面引入修正更为合适。
上述观点未必与主流看法一致，所以如何选择你可以考虑一下，或者我们可以继续讨论。

noodles的困惑

卡开发发 发表于 2021-11-20 17:43
1、仅未固定的原子。
2、VASP会出现构型优化的结果和自洽对不上的情况（即便不变晶格），迄今为止我没弄 ...

谢谢发发老师逐条解答！学到了很多！
关于“表面能”的讨论，我觉得问题在于两个方面：
1、表面能的定义。我之前帖子所说的“slab模型弛豫前后能量的变化是表面能”确实有不妥，还应加上“单位面积”这一条件，不过对于同一研究体系，不知是否可以直接比较能量的差异？
2、对称slab模型和非对称slab模型。对于这两种模型有不同的研究方式：对称slab模型是关于模型中心某个面对称，也就是slab上下面都是研究对象——“表面”，故在计算表面能时“固定模型中心处几层原子模拟体相，上下方的原子允许弛豫”，而非对称slab模型——也就是我用的结构——需要“固定底层原子，上层原子允许弛豫”。所以我的理解是，用slab模型模拟表面然后研究吸附时，选择性弛豫原子（固定底层和弛豫上层）是必要的，只不过根据所用的slab模型的对称性来选择相应的方式，九楼帖子中老师说的“弛豫整个slab”会导致弛豫后整个slab模型与实际的表面会有较大偏离，当然前面所说的slab模型也不是和实际情况完全相符，但应该是比弛豫整个slab来的合理一些？想听一听发发老师和九楼老师的看法~

谢谢!

下图是《Density functional theory A practical introduction》中文版中有关非对称slab模型的部分。

卡开发发

noodles的困惑 发表于 2021-11-20 15:44
感谢回复！
关于第二点，我可以确定，除了测试的NWRITE，其他的输入文件以及参数设置完全相同，这个问题 ...

1、仅未固定的原子。
2、VASP会出现构型优化的结果和自洽对不上的情况（即便不变晶格），迄今为止我没弄清楚全部的原因，虽然有时候单看能量差异并不大。建议精度充分的情况下多cp CONTCAR POSCAR优化几遍再自洽可能能拿到自洽和构型优化一致的结果。
继11楼：

当然我也想过您的这种方式，不过我觉得”固定底层原子，弛豫上层原子”的做法比较合理一些，而且弛豫前后的能量就是“表面能”吧。

这点我不能完全认可，表面能按计算上的定义应该是一个块体劈开为两个板后单位面积所发生的能量变化，此时产生的两个面（即上下面）都是弛豫的。相对来说9楼的说法稍显合理一些，但真正更合理的方式肯定是测试到收敛。且IDIPOLE确实容易导致SCF不收敛一些，建议先把结构优化做出了再加。
2、未必可靠，可能出现的情况是这样的结果自洽和你构型优化能量差很多，甚至可能自洽还是不收敛。所以一般最好建议第一个离子步中的电子步尽可能能做收敛。
3、真空层高度稳妥起见其实也可以进行收敛性测试，不过平面波计算量会与真空层尺寸正相关，更多因素还是用不起。isif=2不会收缩，此时根本就不允许晶格尺寸弛豫了。

noodles的困惑

Epikur 发表于 2021-11-20 15:23
未固定的原子，固定的原子的受力就不管了

根据实际计算的情况，确实是如此。
谢谢老师！

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%，这个差异在可以接受范围之内嘛？或者说这么一个比较的逻辑有意义吗？

感谢老师解答！

Epikur

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

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

noodles的困惑

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

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

Epikur

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

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

noodles的困惑

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

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

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

卡开发发

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，也比较简单，就不详细说了。

Epikur

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

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

问题有点多，谢谢老师！

卡开发发

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

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

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”来显示，还是因为我的计算是错误的呀？