使用cp2k进行Ni(OH)2结构优化过程中SCF/几何优化无法收敛的问题

Eudaimonia

各位老师好，我使用cp2k优化Ni(OH)2表面，参考卡开发发老师的回答进行计算   [综合讨论] 新手刚开始做固体表面吸附，应该怎么样入手呢1.先对下载得到的Ni(OH)2单胞进行了结构优化，使用了5x5x1的k点(现在回顾此处有问题，应使用5x5x5或者类似的k 点更为恰当)，其他参数均为Multiwfn生成cp2k输入文件的默认参数，成功收敛；

2.取出优化好的结构，在切(001)面，模型有两层Ni(OH)2，构建10Å的真空层，使用5x5x1的k点再次进行优化，成功收敛；
3.由于吸附的分子较大，对优化好的(001)表面进行括胞，使用3x3x1的超晶胞进行计算，k点设置为2x2x1，发现SCF在0.9左右的水平上反复震荡，无法收敛(大致循环了80次左右就手动停止了)；之后尝试了设置Ni磁矩为3，整体自旋为7再次进行优化，仍然震荡；增加Cutoff至550eV也不行；
4.尝试将Diag方法更换为OT方法进行计算，SCF可以收敛但是下降速度极其缓慢，打算先粗略进行优化，将SCF收敛限降低至10E-4，几何优化收敛；
6.放置吸附分子至真空层发现真空取的太小，吸附分子距离顶端距离大致为3Å，将晶胞Z轴增加5Å以构建更大的真空层(此处直接修改cif文件中的cell的c坐标，也重新切面构建过，结果基本一致)，重新优化发现几何优化又无法收敛；
7.删去一层Ni(OH)2后再次优化，几何优化收敛；

此时学生非常迷茫，有以下几个问题：
1.很多帖子都表明，OT算法不适用于小晶胞，此处我的模型晶胞大小为12.6x12.6x18.5，感觉仍然是小晶胞，不确定使用OT算法得到的结果是否合理？
2.之前使用Diag方法进行计算时SCF反复震荡是否还有其他方法可以解决，感觉Diag+K点应该更为合适？
3.晶胞Z轴大小和几何优化是否有关系？
4.对于Ni(OH)2这样的结构，单层Ni(OH)2可以用于后续吸附计算吗？
5.卡卡开发老师之前提到，如果使用优化好的带切面的晶胞，理论来说扩胞后几何优化直接收敛，但是学生目前手头的一些结构，每一次扩胞后仍然需要多次(一般>20)几何优化才能够成功收敛，这样的情况是正常的吗？

感谢各位老师的回答！
由于实际操作时删去了很多中间得到的数据文件，暂时无法上传全部的输入输出文件，不好意思！
2021.12.26 UpDate
根据卡开发发老师的建议，重新对Ni猜测磁矩为3，增加了Smear进行了单胞及001切面的优化，还是可以正常完成。441晶面仍然在初始SCF计算时震荡，221、551、552的k点都进行了尝试，均震荡，相，请问老师还有什么可以进行的尝试？ undefined
Note：k点设置为552的文件不小心被覆盖，进行了约60步的SCF计算，震荡，实际电脑不太算的动这个等级的k点
目前打算更改一下超胞大小再进行尝试

2021.12.28 UpDate
根据卡开发发老师的进一步指点，该体系为反铁磁性(AFM)体系，模型中Ni原子的磁矩设定应为一层+3，一层-3，且可以打开对称性以节省kpoint；
同时对此进行更改以后SCF成功收敛，但是后续还是会输出一些警告，但是不影响计算结果的输出，有兴趣的老师可以移步至后续楼层进行讨论
再次感谢各位老师的帮助

卡开发发

0、（1）我估计有可能大概5x5x1或者5x5x2可能就可以收敛，你最好测试看看。
（3）另外Ni(OH)2应该是有磁性的，从你的结构看，你找的很可能是那个反铁磁的，按道理应该在优化晶格的时候考虑。如果你做的是多层的，这个多层也应当是反铁磁性的。
（4）不推荐做所谓的粗略优化，能量和力的稳定计算是结构优化最基本的保障，否则精度不足引起的力噪声导致优化出来更奇怪的结构只会更浪费时间。
（7）删掉层高能够收敛可能是因为你一共只有两层？我猜测还是磁矩猜错引起的问题。
最后再次重申猜测磁矩的重要性。
1、OT方法不适合小晶格，小晶格需要足够的k点甚至需要足够的展宽，不止如此，不建议OT在窄带体系尤其是金属性体系（比如纯金属块体或表面）使用，虽然你的体系还不涉及到这个问题。
2、按道理使用常规的对角化和多k点应该是最稳定的方式，只是标度不如OT，效率要差一些。
3、有关。你仔细思考一下周期性边界的概念，真空小的话，你的分子就被夹在两块很近的表面当中，并不符合实际的情况；即便假定Poisson方程引入了边界修正，那么也应该留有足够的真空区域保证密度和静电势求解有接近正确的边界行为。
4、可以用于后续吸附，但是你得思考你实际研究的问题到底是单层还是多层的。
5、你用标准对角化+多k点还有充分展宽的方式优化一下1x1的表面，然后再扩胞试试，这是原理上应当成立的结论。

丁越

卡开发发老师表示很无奈：“卡卡开发老师”
另外你上传一下你的cif文件呀

Eudaimonia

卡开发发 发表于 2021-12-25 18:39
0、（1）我估计有可能大概5x5x1或者5x5x2可能就可以收敛，你最好测试看看。
（3）另外Ni(OH)2应该是有磁性 ...

非常感谢卡开发发老师的帮助！
(1)目前打算先猜测磁矩重新计算，后续有空更换一下k点设置再测试一下
(3)学生不是特别确定应该在什么情况下优化晶格 在 该帖子中 sob老师认为只要计算不出现崩溃晶胞优化并不重要，请问卡开发发老师您怎么看？
(4)谢谢老师的教导，以后不会这么做了
(7)同(1),再次感谢老师的教导

1.学生有一点好奇为什么会这样，老师有空可以推荐一下相关学习资料吗？
3.这里学生的疑问主要在于平时一般认为原子之间距离>5Å，产生的相互作用就比较局限，我猜测真空层厚度从10Å提高到15Å对体系的影响会非常小，所以有这样的疑问。现在还是不确定究竟多少的真空区域能够保证密度和静电势求解有接近正确的边界行为。或者可以考虑计算时使用xy方向的周期性？
4.嗯嗯，根据表征结果应该还是一个多层Ni(OH)2的模型更为合理。此处的疑问在于我们认为“多厚”的体相可以近似认为模型是一个合理的表面。目前我觉得6层以上的原子可以认为大致是一个合理的表面，请问老师可以这样认为吗？
5.正在尝试，在VASP里展宽似乎是通过ISMEAR来实现，在cp2k里是通过&SMEAR字段来设置的吗？

再请问老师：
1.之前有帖子说设置原子磁矩后体系的自旋设置已经无意义，请问是这样的吗？因为小晶胞中金属原子数较少，设置自旋数较小没有疑问，可是大晶胞中金属原子数增加需要设置几十甚至上百的自旋，感觉非常的不靠谱。
再次感谢老师的帮助！

Eudaimonia

丁越 发表于 2021-12-26 10:14
卡开发发老师表示很无奈：“卡卡开发老师”
另外你上传一下你的cif文件呀

啊，我没确认居然犯了这么低级的错误（跪）
经过卡开发发老师的指点，目前打算重新进行一轮计算，旧的cif文件先不进行上传了。
PS. inp文件可以直接用Multiwfn打开，里面也储存了相关信息

卡开发发

Eudaimonia 发表于 2021-12-26 12:26
非常感谢卡开发发老师的帮助！
(1)目前打算先猜测磁矩重新计算，后续有空更换一下k点设置再测试一下
(3 ...

0、（3）有点搞事情啊。你要看清楚研究的问题，该贴中后续是要进行一个MD的过程，如果后续是NPT，晶格在指定的温度压力下会进一步被弛豫。但如果后续要构造表面，建议还是优化一下比较好。
1、smearing和多k点情形没办法用OT，算法的具体实现我也没来得及研究，不能完全说清cp2k没做完善是因为该方法不可被扩展还是他们就没去做。
3、指的其实就是考虑xy周期性的时候也得留有足够的真空，实际情况5A可能不够，最好测试下。比如一些修正方法（如Markov-Payne）可能建议是12A以上，不过这个确实看体系，有些体系库伦作用衰减快一些，也许要的比较短，你可以测试真空层高对能量和力的影响。
4、几层这个我不好说，得根据你们实验表征的情况，建议找实验的同事好好讨论一下。如果是比较厚的情况也是设置到足够就可以了，也可以通过测试确定下来。
5、对。

其他：
1、这个要看情况，最主要的原因是对于固体这种准连续能级的情况自旋多重度本身可能没什么意义。比如O2分子的计算，你给两个O分别做磁矩为1的猜测，得到的也应该是三重态。另外，PBC下的磁矩应该按照per cell或per site来讨论，总的磁矩应该是个广延量。

还行，没叫我开开发卡就不错了，ID不重要，但是科学问题不能讨论那么马虎。

丁越

Eudaimonia 发表于 2021-12-26 12:30
啊，我没确认居然犯了这么低级的错误（跪）
经过卡开发发老师的指点，目前打算重新进行一轮计算，旧的ci ...

我的意思是只有知道了Ni(OH)2是哪种结构，然后我才能去materials project看看这个材料的磁性之类的相关性质，然后才能考虑inp文件是否设置恰当。

Eudaimonia

卡开发发 发表于 2021-12-26 14:19
0、（3）有点搞事情啊。你要看清楚研究的问题，该贴中后续是要进行一个MD的过程，如果后续是NPT，晶格在 ...

（3）.哦哦，不好意思犯了这样的错误
我尝试了一下，对原始的NI(OH)2进行了晶胞优化，但是发现晶胞参数发生了非常夸张的变化，从3.16x3.16x9.35变化为7.56x7.7x4.95(还未优化完)，但大体结构还在，感觉也不合理，请问老师可能是哪里可能有问题？对Ni猜测磁矩为3，增加了D3校正和Smear，其他参数使用Multiwfn默认参数
1.3.4.5 感谢老师的指点，麻烦老师了

其他：
重新对Ni猜测磁矩为3，增加了Smear进行了单胞及001切面的优化，还是可以正常完成。441晶面仍然在初始SCF计算时震荡，221、551、552的k点都进行了尝试，均震荡，相关文件更新在1楼，请问老师还有什么可以进行的尝试？
Note：k点设置为552的文件不小心被覆盖，进行了约60步的SCF计算，震荡，实际电脑不太算的动这个等级的k点
目前打算更改一下超胞大小再进行尝试

卡开发发

Eudaimonia 发表于 2021-12-26 21:57
（3）.哦哦，不好意思犯了这样的错误
我尝试了一下，对原始的NI(OH)2进行了晶胞优化，但是发现晶胞 ...

（3）反铁磁（AFM）的情况Ni有磁矩为+3的，也有为-3的，z方向的层和层之间Ni的磁矩应该是相反的。前面专门提了反铁磁的事情，你可能没注意。
后续表面也是相同的问题。
表面体系应该5 5 1就行，5 5 2是针对块体的，你可以好好测试下单点，看看在能量受力等差异不大的情况可以使用稍小一些的k点，另外&KPOINTS这里，SYMMETRY F完全没必要，如果晶体本身有对称性可以用这个对称性来使用更少的（不可约）k点进行计算，有些k的贡献是一样的；如果能确定对称性，晶胞优化也尽可能使用对称性。这体系尽量猜测原子磁矩并且自动优化磁矩。

Eudaimonia

手滑恢复了错误的数据。。。。重新再写一次感谢卡开发发老师的进一步指点，反铁磁的情况我之前确实没有注意到，更改后SCF成功收敛，再一次感谢老师的指点和帮助！
考虑到吸附分子大小及计算时间，将441超晶胞改为331，对331超晶胞进行了k点测试，结果如下

k点设置为221时与551的能量差距约1.4meV，体系原子数160，k点选用221以节省计算时间
虽然体系成功收敛，但是还是有几个WARNING出现在输出文件中，如果老师有时间可以指点一下，
1.out文件输出警告 WARNING in cryssym.F:163 :: Symmetry library SPGLIB not available，应该指对称性虽然打开但是并没有被实际应用，请问应该是否应该在cif文件导出前寻找对称性或者其他什么操作找到对称性以成功应用？
Note：再次编辑到此处意识到一个问题，我所谓的打开对称性值得是将SYMMETRY F改为SYMMETRY T，程序可能将T识别为一种对称性出现此报错，应该将此处改为其他；
2.输出警告 WARNING in qs_environment.F:1163 :: Kpoints: Different number of MOs requested. The number of beta MOs will be set to the number alpha MOs ；此处应该指ADDED_MOS字段后给予alpha beta的设定不相等，但实际我设定的值均为100，不知道为什么没有被程序识别；
3.输出警告WARNING in qs_scf_post_gpw.F:2190 :: Spin contamination estimate not implemented for k-points.  似乎和k points + Broyden mixing + Fermi smearing共同设置引起的，具体讨论可以 看此处的讨论(需fq)，由于不影响计算结果且看不太懂，此处将原文黏贴如下，作为记录：
Your system is a bit more advanced than a simple gas phase O2 molecule. Moreover, you use k points with Broyden mixing and Fermi smearing and thus you won’t see the printout “Ideal and single determinant S**2” as indicated by the warning:
*** WARNING in qs_scf_post_gpw.F:1717 :: Spin contamination estimate not ***
*** implemented for k-points
Alternatively, you may check the final line of the Mulliken population analysis

4. 在后续441超晶胞的晶面优化中，未看到卡开发发老师12.27 22:07的回复，设置SYMMETRY F进行计算，此时MULTIPLICITY为49(Multiwfn程序自动输出，我在设置磁矩时只通过Multiwfn设置了一半Ni磁矩为3，另一半Ni磁矩手动更改)，SCF仍然在1附近反复震荡，但是改为SYMMETRY T，MULTIPLICITY随手设置为7(根据之前讨论此处设置可能无意义)后再计算就可收敛，但是根据此楼问题1.，对称性并未成功应用，为什么会出现这个差异呢？相关文件在本楼中；
Note：可能和MULTIPLICITY设置有关，暂时未进行测试
5.主楼问题5
"卡开发发老师之前提到，如果使用优化好的带切面的晶胞，理论来说扩胞后几何优化直接收敛，但是学生目前手头的一些结构，每一次扩胞后仍然需要多次(一般>20)几何优化才能够成功收敛，这样的情况是正常的吗？"
在优化时，441超晶胞设置SYMMETRY T，MULTIPLICITY设置为7，221 k点，Opt step4步即收敛，和括胞后直接收敛的描述比较接近；但是相同设置下331的超晶胞优化了100+步才收敛，还是存在这个问题，不明白是为什么；
Note：猜测可能是k点设置引起的，但是设置k点为111,331时也会在4步Opt step后突然大大远离收敛限


再次感谢各位老师尤其是卡开发发老师的指点！
由于主楼附件过多，此楼提到的数据文件上传在本楼

sobereva

以后发帖时别在贴子里出现fan qiang，过于敏感，给你改了

卡开发发

针对末尾的问题：1、对称性的问题实在不行你可以用一些程序约束磁矩差异或者临时将spin down的原子换成别的元素，重新找一下对称性试试，我不确定是不是程序没办法正确处理AFM的对称性。
2、不明，不知道程序内部怎么设计的。
3、那个Warning其实可以不用管，对窄带系统其实可以不要管自旋污染的问题，上面帖子提到过原因。
4、按道理两个Ni磁矩相反，总的多重度应该为1。
5、先看看纠正上述问题是否还有这个情况，不行再根据输出情况进行诊断。
另外，收敛性测试最好是在单点上做，然后看能量和力的收敛情况，这样能有效避免构型优化震荡。如果这种情况下还出现震荡，要么初始结构不合理，要么可以适当调整构型优化中算法的参数（比如步长控制或Hessian更新之类的算法）。

Eudaimonia

卡开发发 发表于 2021-12-28 21:08
针对末尾的问题：1、对称性的问题实在不行你可以用一些程序约束磁矩差异或者临时将spin down的原子换成别的 ...

谢谢卡开发发老师的指点，老师元旦快乐~
前几天去上海做实验了，没有及时回复和对问题进行测试，不好意思
1.对论坛中进行了检索，没有找到对磁性进行约束的讨论，请问一下老师能推荐一下相关程序吗？
4.嗯嗯
5.后续会对此处重新进行测试

谢谢老师的提示，论坛中有很多关于K点测试的讨论，10楼中的测试结果是单点的结果(不过取用的模型是已经优化过的331超胞(就是经过上百次Opt step重新收敛后的)，可能和单胞扩胞直接测试有一些差别，后续可能重新测试并更新再后面楼层，最近搬家+期末事情比较杂乱

卡开发发

Eudaimonia 发表于 2022-1-1 11:04
谢谢卡开发发老师的指点，老师元旦快乐~
前几天去上海做实验了，没有及时回复和对问题进行测试，不好意 ...

1、磁性约束指的是针对不等价磁矩设置的原子应当被识别为是不同类型的原子，目前能做到的例如MS或Pymatgen，在这个基础上找对称性看看获得的原胞如何。

欢迎讨论。

zmjsce

卡开发发 发表于 2021-12-28 21:08
针对末尾的问题：1、对称性的问题实在不行你可以用一些程序约束磁矩差异或者临时将spin down的原子换成别的 ...

卡开发发老师您好，请问使用CP2K优化Ni(OH)2时，在已经考虑磁性情况下，是否还需要考虑DFT+U呢？