计算化学公社
标题: 线性响应计算U值问题求助 [打印本页]
作者Author: physics_xw 时间: 2021-5-11 21:06
标题: 线性响应计算U值问题求助
问题背景
本人研究方向为二维本征铁磁性材料的第一性原理研究，目前的主要兴趣集中于单层或双层Cr基化合物(如CrI3、1T-CrSe2及1T-CrTe2等)。按照领域内的相关经验，含Cr的氧族或卤族化合物在处理时均有必要考虑电子的强关联效应，使用DFT+U进行计算。如果只是计算已知材料，完全可以参考一些已有文献进行U值的设置，但考虑到接下来会对一些新材料进行预测，那么自己测试U值将是必要的技能。所以本人近期学习了线性响应计算U值的方法，尝试重复相关文献的结果，但遇到了一些问题，不知道是否是自己操作上有不恰当的地方，还希望能有大佬提点一下。
本人学习线性响应测试U值参考的以下网站：
操作步骤
1.单点计算
不设置U值相关参数，设置LMAXMIX=11，将POSCAR中的Cr分成两部分，其中第一部分原子数为1，剩下的为第二部分，POTCAR中相应Cr的赝势也重复一次。之后仅进行常规的DFT单点能计算，写出WAVECAR和CHGCAR，记录OUTCAR文件末尾total charge中的单个Cr原子d电子占据。
2.非自洽计算(nscf)
添加U值相关参数，设置LDAUTYPE=3，并将所要施加扰动的单个Cr原子的LDAUU和LDAUJ设置为相应势α。读取1中电荷密度文件，并通过ICHARG=11进行非自洽计算。同样，通过OUTCAR末尾可以知道d电子占据数，通过施加不同的扰动可以得到相应变化的占据数。
3.自洽计算(scf)
U值相关参数设置与2相同，所不同的仅仅为取消ICHARG=11，即进行自洽计算。通过施加不同扰动也可以得到相应变化的占据数。
以1T-CrTe2为例，所得结果如下图所示。(一套典型的输入文件见附件)
第一列表示施加的扰动，第一行的111、221、331分别代表1×1×1、2×2×1、3×3×1的超胞规模，将d电子占据数与α线性拟合得到斜率如倒数第二行所示，nscf与scf相应斜率倒数之差即为有效U值。我一般取所计算的最大超胞得到的值为最后的有效U值。对于1T-CrSe2及CrI3也使用相同方法操作，得到类似的结果。有效U值均在5.5-6 eV之间。
遇到问题
1.文献[1][2]使用线性相应方法计算得到1T-CrSe2的有效U值分别为4eV和4.5eV，而我的计算结果为5.8 eV；文献[3]给出CrI3有效U值为2.8eV，而我的计算结果为6.0 eV。差异非常大，所以我想应该不能说是简单的误差，可能本人在这些步骤的操作上或是对这种方法的理解上有错误，只能尽可能将我所做的操作描述出来，求助于大佬们指出我其中的错误。
2.按照本人粗浅的理解，设置LDAUTYPE=2进行DFT+U计算时，起作用的是U-J，即有效U值。既然这样，那么设置J值还有意义吗?为什么相关文献中看到更多的是给了一个非零的J值，这些J值一般是通过什么方法获得的呢？这样设置非零J值和将J设置为0有什么本质的区别吗？
[1] Li, Bo, et al. “Van der Waals epitaxial growthof air-stable CrSe2 nanosheets with thickness-tunable magnetic order.” NatureMaterials (2021): (1-8).
[2] Wang, Cong, et al. “Bethe-Slater-curve-likebehavior and interlayer spin-exchange coupling mechanisms in two-dimensionalmagnetic bilayers.” Physical Reviewer B 102.2 (2020): 020402.
[3] Jiang, Peiheng, et al. “Stacking tunableinterlayer magnetism in bilayer CrI3.” Physical Review B 99.14 (2019): 144401.
作者Author: 卡开发发 时间: 2021-5-12 01:03
1、(1)我不确定你是如何响应的，看上去你可能只是做了一个位点的。实际应当考虑逐个扰动每个Hubbard位点，并且收集被扰动和其他涉及到的Hubbard位点的占据数，这样构造出来的slope其实是个矩阵（其实也就是原文的响应矩阵），然后应该是对bare（nscf）和screen(scf)情况得到的响应矩阵的逆进行做差，去对角元的值。如果你还是不太明确，可以参考Kulik组的博文（http://hjkgrp.mit.edu/content/ca ... -u-periodic-systems），里面的描述会更详细一些。当然你可以看到，实际上这样线性响应计算U的过程十分繁琐。
(2)如果晶胞足够大，k和动能截断足够大，那么可能没错/也许因为采用的赝势不同或者其他参数不同，会导致U出现不同。我私下测试了QE采用DFPT的方法来进行线性响应计算U，结果超软赝势、PAW和模守恒赝势得到的U值都不太一样。所以U我不认为是一个可以随意移植的参量。
2、LDAUTYPE=2的情况，Ueff=U-J才是有意义的，至于单独的U和J本身就没什么实际意义。但对于LDAUTYPE=1，J值就具备意义。例如QE中也提供了对J的线性响应方案（甚至包含双中心的DFT+U+V框架），但VASP可能缺乏此功能。
作者Author: jiangning198511 时间: 2021-5-12 10:40
之前看到一篇文章讲利用LDA+U 计算不同自旋态的能量差，他们用高自旋得到一个U值，然后在低自旋下得到一个U值，然后利用这两个值得到E(HS)-E(LS)值，对比不同泛函后发现这个结果误差最小。 然后发了JCTC. 我这个就有个疑问，你两个不同的U值相当于给出两个不同的LDA+U结果，这样的结果求差值和相同泛函下的的不同自旋差值能做对比吗？极端一点的话 我高自旋用一种泛函，低自旋用一种泛函，然后对比高精度方法看看哪个更小，岂不是也算一种计算方法的开放研究 
PS: 我用公式再表达一下
文章方法：E(HS-LS)=E(HS: LDA+U1, U1为高自旋下线性响应得到)-E(LS，LDA+U2, U2为低自旋下线性响应得到)
常规方法(以PBE0为例子)：E(HS-LS)=E(HS: PBE0)-E(LS，PBE0)
PPS: 文章：10.1021/acs.jctc.1c00034
作者Author: 卡开发发 时间: 2021-5-12 11:42
其实不同的化学环境下按道理不同U可以比较，还是以Kulik他们的帖子（http://hjkgrp.mit.edu/content/in ... -u-variations-dftur）为例子，U实际上是结构相关的。U所包含的局域排斥修正其实也应该隐含屏蔽效应。如果从cRPA响应U的角度说，甚至U还可以含频。
作者Author: physics_xw 时间: 2021-5-12 15:36
感谢卡老师的回复。
1.(1)确实如您所说，我只做了一个位点的计算。接下来我会再补齐其他位点的计算，按照您说的方法去计算。
随之而来有一个新问题，您发出来的链接博文中给出了博主计算的响应矩阵，我将所给的两个矩阵分别求逆，然后相减，取对角线元素，得到的值约为3.4 eV，而博文中给出的为2.9 eV。将相差的0.5 eV看作数值计算的误差合理吗？
(2)我在查阅相关文献时进行过筛选，只留下了PAW势PBE泛函的VASP计算结果。可能确实如您上面所说，近期我会再确认一下。
2.感谢您的指导。突然想尝试学习QE了
作者Author: 卡开发发 时间: 2021-5-12 16:58
1、（1）有可能，比如保留的小数等因素引起的，但我暂时不确定。
（2）可以复现文献试试看。
2、平面波方法基本原理相通，QE有QE的优势，VASP有VASP的优势，根据要做的事情选择擅长的程序即可（大前提是有版权）。
作者Author: nitro123 时间: 2021-8-23 16:13
你好，我最近也遇到了使用VASP线性响应计算U值偏大的问题。
如果愿意的话，方便加QQ或微信或者其他聊天软件聊一下这个问题吗？因为我也苦于这个问题😂
其实我看到一篇文章是标准线性响应计算U的改版，他是要求先加一个适当的U值，然后再去线性响应计算U
作者Author: 卡开发发 时间: 2021-8-23 23:34
是，因为直接响应其实还是用的plain DFT的密度，U肯定和真实密度下响应得到的有差别。
作者Author: ppy 时间: 2021-12-26 19:39
现在QE好像可以用hp.x去计算U值了：https://gitlab.com/QEF/q-e/-/tre ... 07974ba/HP/examples
作者Author: 荻花花花花 时间: 2022-11-4 22:19
你好，文我想请教您关于hp.x算U的问题可以吗？请问得出U首先是得出了两个矩阵。请问矩阵中的数据是从输出文件中的哪里得到的呢？矩阵的行数列数又由什么来决定呢？
作者Author: 卡开发发 时间: 2022-11-5 16:49
hp.x可以直接算出U值，不需要再去手动写脚本来得到响应矩阵求逆。
你可以看下hp.x下面的例子即可。
作者Author: 荻花花花花 时间: 2022-11-11 15:42
请问一下，响应矩阵的数据如何理解呢？我了解到使用vasp中线性响应的方法计算U值时，矩阵中的数据是每个Hubbard位点的d电子数目对于a势的响应。请问使用hp.x去计算U值的这种方法，矩阵中的数据如何理解呢？矩阵中的数据是怎么从输出文件中获取的呢？
作者Author: 卡开发发 时间: 2022-11-11 16:21
本帖最后由 卡开发发 于 2022-11-11 16:43 编辑
hp.x的输出不需要去读那个矩阵，程序会输出一个叫${prefix}.Hubbard_parameters.dat的文件，会按照你的pw.x设置的hubbard site来输出，例子里面大概是这样的格式：
- =-------------------------------------------------------------------=
- Hubbard U parameters:
- site n. type label spin new_type new_label Hubbard U (eV)
- 1 1 Co 1 1 Co 7.7514
- 2 2 O 1 2 O 8.0439
- 3 2 O 1 2 O 8.0439
- =-------------------------------------------------------------------=
QE用hp.x计算U要简单自动化的多，不需要自己手动构造响应矩阵（即局域占据数对响应势的导数），但是类似于原始实空间下做响应使用超胞，虽然现在的方法不需要指定超胞，但需要对响应的q点尺寸进行指定，这个原则上说也应当进行收敛性测试。你要是要用下面响应矩阵手动算，应该是用U=chi0的逆矩阵-chi的逆矩阵，然后再按照对角元取，实际上这两部分也被打印到那个.dat文件当中。新的VASP应该可以通过cRPA来计算U，不过我没做过相关测试。
这两种新的算法细节我现在还没有仔细研究过，所以底层实现我还不能和你说的很清楚。
作者Author: 鱼丸box 时间: 2023-3-18 01:26
大佬您好，看到这个输出文件，我想询问您一下，这意味着在LiCoO2体系中也需要对O原子+U吗？我记得在这个例子中也有计算V值的HUBBARD.dat文件，其中Co原子也加入了一个比较大的V值，这是什么原因呢？以及在后续计算中应该如何根据${prefix}.Hubbard_parameters.dat和HUBBARD.dat设置需要输入的U和V值呢？非常感谢
作者Author: 卡开发发 时间: 2023-3-18 08:46
简单说，就是把哪些原子作为Hubbard位点，然后考虑何种模型（+U还是+U+V）取决于进行hp.x前pw.x的计算参数，然后响应得到的参数是什么你再按照格式填回pw.x当中对应的参数进行后续计算即可。
作者Author: 鱼丸box 时间: 2023-3-18 14:42
谢谢卡老师回复，但是现在的我依然有两个疑问：
（1）关于老师回答中的"响应得到的参数是什么你再按照格式填回pw.x当中对应的参数进行后续计算即可"这句话，按照QE中的例10（LiCoO2）的输入文件在Co的3d轨道上+U，以及Co-3d对O-2p轨道作用上+V，为什么在LiCoO2.Hubbard_parameters.dat文件中会得出O原子的U值呢？也就是说，我需要在后续计算中也对O原子+U吗？
（2）如果只针对+U的话，我们进行线性响应计算U值，这个U值指的是Ueff吗？如果是，我们又该在哪里寻找J值呢？
作者Author: 卡开发发 时间: 2023-3-18 21:11
1、如果你要进行+U+V的计算，U应该参考LiCoO2.Hubbard_parameters.dat，V应该参考Hubbard.dat的数据。因为V当中设置了O-2p，那么U也会被算出来。
2、如果你按照DFT+U（也就是旋转不变形式）进行设置，此时U就是Ueff=U-J，但DFT+U+J目前看hp.x程序描述似乎没办法得到单独的J，有可能得通过pw.x在实空间下手动做，具体我不那么确定。
作者Author: 鱼丸box 时间: 2023-3-19 15:27
谢谢卡老师的回答，对于问题（1）我会先进行一些测试，对于（2），我的理解是pwscf的DFT+U相当于VASP中的LDAUTYPE=2，并设置了J=0。
作者Author: 卡开发发 时间: 2023-3-19 16:31
确实是这样。
作者Author: Tang-Tommy 时间: 2023-10-22 18:10
卡老师，请问一下，我使用Vasp对一个尖晶石体系进行计算，Zn、Fe、O元素。用DFT+U进行结构优化，然而结构优化提交任务后，每次得到的能量都不一样，这个是正常的吗？
作者Author: 卡开发发 时间: 2023-10-23 00:04
两次计算的OUTCAR和OSZICAR或STDOUT（也就是屏幕输出内容）发来看看，暂时不好说。特定的U值可能确实会导致计算不稳定。
作者Author: Tang-Tommy 时间: 2023-10-23 09:18
请问一下老师，怎么把这个文件发给您过目呢？
作者Author: 鱼丸box 时间: 2024-6-7 22:46
博主您好，请问最后您解决了这个问题吗？我使用了多位点计算，但是结果依然是5.6-5.9之间
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |