计算化学公社

 找回密码 Forget password
 注册 Register
Views: 18842|回复 Reply: 22
打印 Print 上一主题 Last thread 下一主题 Next thread

[VASP] 线性响应计算U值问题求助

[复制链接 Copy URL]

9

帖子

0

威望

1141

eV
积分
1150

Level 4 (黑子)

问题背景
    本人研究方向为二维本征铁磁性材料的第一性原理研究,目前的主要兴趣集中于单层或双层Cr基化合物(CrI31T-CrSe21T-CrTe2)。按照领域内的相关经验,含Cr的氧族或卤族化合物在处理时均有必要考虑电子的强关联效应,使用DFT+U进行计算。如果只是计算已知材料,完全可以参考一些已有文献进行U值的设置,但考虑到接下来会对一些新材料进行预测,那么自己测试U值将是必要的技能。所以本人近期学习了线性响应计算U值的方法,尝试重复相关文献的结果,但遇到了一些问题,不知道是否是自己操作上有不恰当的地方,还希望能有大佬提点一下。
    本人学习线性响应测试U值参考的以下网站:
操作步骤
1.单点计算
    不设置U值相关参数,设置LMAXMIX=11,将POSCAR中的Cr分成两部分,其中第一部分原子数为1,剩下的为第二部分,POTCAR中相应Cr的赝势也重复一次。之后仅进行常规的DFT单点能计算,写出WAVECARCHGCAR,记录OUTCAR文件末尾total charge中的单个Cr原子d电子占据。
2.非自洽计算(nscf)
    添加U值相关参数,设置LDAUTYPE=3,并将所要施加扰动的单个Cr原子的LDAUULDAUJ设置为相应势α。读取1中电荷密度文件,并通过ICHARG=11进行非自洽计算。同样,通过OUTCAR末尾可以知道d电子占据数,通过施加不同的扰动可以得到相应变化的占据数。
3.自洽计算(scf)
        U值相关参数设置与2相同,所不同的仅仅为取消ICHARG=11,即进行自洽计算。通过施加不同扰动也可以得到相应变化的占据数。
1T-CrTe2为例,所得结果如下图所示。(一套典型的输入文件见附件)

    第一列表示施加的扰动,第一行的111221331分别代表1×1×12×2×13×3×1的超胞规模,将d电子占据数与α线性拟合得到斜率如倒数第二行所示,nscfscf相应斜率倒数之差即为有效U值。我一般取所计算的最大超胞得到的值为最后的有效U值。对于1T-CrSe2CrI3也使用相同方法操作,得到类似的结果。有效U值均在5.5-6 eV之间。
遇到问题
1.文献[1][2]使用线性相应方法计算得到1T-CrSe2的有效U值分别为4eV4.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.

U-1T-CrTe2.jpg (100.65 KB, 下载次数 Times of downloads: 73)

1T-CrTe2计算结果

1T-CrTe2计算结果

InputFiles.rar

157.06 KB, 下载次数 Times of downloads: 67

InputFiles

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

2#
发表于 Post on 2021-5-12 01:03:10 | 只看该作者 Only view this author
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可能缺乏此功能。

评分 Rate

参与人数
Participants 2
eV +8 收起 理由
Reason
physics_xw + 5 谢谢
sobereva + 3

查看全部评分 View all ratings

日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

689

帖子

2

威望

4104

eV
积分
4833

Level 6 (一方通行)

3#
发表于 Post on 2021-5-12 10:40:52 | 只看该作者 Only view this author
卡开发发 发表于 2021-5-12 01:03
1、(1)我不确定你是如何响应的,看上去你可能只是做了一个位点的。实际应当考虑逐个扰动每个Hubbard位点, ...

之前看到一篇文章讲利用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

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

4#
发表于 Post on 2021-5-12 11:42:09 | 只看该作者 Only view this author
jiangning198511 发表于 2021-5-12 10:40
之前看到一篇文章讲利用LDA+U 计算不同自旋态的能量差,他们用高自旋得到一个U值,然后在低自旋下得到一 ...

其实不同的化学环境下按道理不同U可以比较,还是以Kulik他们的帖子(http://hjkgrp.mit.edu/content/in ... -u-variations-dftur)为例子,U实际上是结构相关的。U所包含的局域排斥修正其实也应该隐含屏蔽效应。如果从cRPA响应U的角度说,甚至U还可以含频。
日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

9

帖子

0

威望

1141

eV
积分
1150

Level 4 (黑子)

5#
 楼主 Author| 发表于 Post on 2021-5-12 15:36:26 | 只看该作者 Only view this author
卡开发发 发表于 2021-5-12 01:03
1、(1)我不确定你是如何响应的,看上去你可能只是做了一个位点的。实际应当考虑逐个扰动每个Hubbard位点, ...

感谢卡老师的回复。
1.(1)确实如您所说,我只做了一个位点的计算。接下来我会再补齐其他位点的计算,按照您说的方法去计算。
    随之而来有一个新问题,您发出来的链接博文中给出了博主计算的响应矩阵,我将所给的两个矩阵分别求逆,然后相减,取对角线元素,得到的值约为3.4 eV,而博文中给出的为2.9 eV。将相差的0.5 eV看作数值计算的误差合理吗?
(2)我在查阅相关文献时进行过筛选,只留下了PAW势PBE泛函的VASP计算结果。可能确实如您上面所说,近期我会再确认一下。
2.感谢您的指导。突然想尝试学习QE了

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

6#
发表于 Post on 2021-5-12 16:58:48 | 只看该作者 Only view this author
physics_xw 发表于 2021-5-12 15:36
感谢卡老师的回复。
1.(1)确实如您所说,我只做了一个位点的计算。接下来我会再补齐其他位点的计算,按 ...

1、(1)有可能,比如保留的小数等因素引起的,但我暂时不确定。
(2)可以复现文献试试看。
2、平面波方法基本原理相通,QE有QE的优势,VASP有VASP的优势,根据要做的事情选择擅长的程序即可(大前提是有版权)。

评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
physics_xw + 5 谢谢

查看全部评分 View all ratings

日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

2

帖子

0

威望

14

eV
积分
16

Level 1 能力者

7#
发表于 Post on 2021-8-23 16:13:54 | 只看该作者 Only view this author
你好,我最近也遇到了使用VASP线性响应计算U值偏大的问题。
如果愿意的话,方便加QQ或微信或者其他聊天软件聊一下这个问题吗?因为我也苦于这个问题😂
其实我看到一篇文章是标准线性响应计算U的改版,他是要求先加一个适当的U值,然后再去线性响应计算U

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

8#
发表于 Post on 2021-8-23 23:34:38 | 只看该作者 Only view this author
nitro123 发表于 2021-8-23 16:13
你好,我最近也遇到了使用VASP线性响应计算U值偏大的问题。
如果愿意的话,方便加QQ或微信或者其他聊天软 ...

是,因为直接响应其实还是用的plain DFT的密度,U肯定和真实密度下响应得到的有差别。
日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

1

帖子

0

威望

161

eV
积分
162

Level 3 能力者

9#
发表于 Post on 2021-12-26 19:39:42 | 只看该作者 Only view this author
nitro123 发表于 2021-8-23 16:13
你好,我最近也遇到了使用VASP线性响应计算U值偏大的问题。
如果愿意的话,方便加QQ或微信或者其他聊天软 ...

现在QE好像可以用hp.x去计算U值了:https://gitlab.com/QEF/q-e/-/tre ... 07974ba/HP/examples

5

帖子

0

威望

39

eV
积分
44

Level 2 能力者

10#
发表于 Post on 2022-11-4 22:19:58 | 只看该作者 Only view this author
ppy 发表于 2021-12-26 19:39
**** 作者被禁止或删除 内容自动屏蔽 ****

你好,文我想请教您关于hp.x算U的问题可以吗?请问得出U首先是得出了两个矩阵。请问矩阵中的数据是从输出文件中的哪里得到的呢?矩阵的行数列数又由什么来决定呢?

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

11#
发表于 Post on 2022-11-5 16:49:54 | 只看该作者 Only view this author
荻花花花花 发表于 2022-11-4 22:19
你好,文我想请教您关于hp.x算U的问题可以吗?请问得出U首先是得出了两个矩阵。请问矩阵中的数据是从输出 ...

hp.x可以直接算出U值,不需要再去手动写脚本来得到响应矩阵求逆。
你可以看下hp.x下面的例子即可。
日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

5

帖子

0

威望

39

eV
积分
44

Level 2 能力者

12#
发表于 Post on 2022-11-11 15:42:36 | 只看该作者 Only view this author
卡开发发 发表于 2022-11-5 16:49
hp.x可以直接算出U值,不需要再去手动写脚本来得到响应矩阵求逆。
你可以看下hp.x下面的例子即可。

请问一下,响应矩阵的数据如何理解呢?我了解到使用vasp中线性响应的方法计算U值时,矩阵中的数据是每个Hubbard位点的d电子数目对于a势的响应。请问使用hp.x去计算U值的这种方法,矩阵中的数据如何理解呢?矩阵中的数据是怎么从输出文件中获取的呢?

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

13#
发表于 Post on 2022-11-11 16:21:16 | 只看该作者 Only view this author
本帖最后由 卡开发发 于 2022-11-11 16:43 编辑
荻花花花花 发表于 2022-11-11 15:42
请问一下,响应矩阵的数据如何理解呢?我了解到使用vasp中线性响应的方法计算U值时,矩阵中的数据是每个H ...

hp.x的输出不需要去读那个矩阵,程序会输出一个叫${prefix}.Hubbard_parameters.dat的文件,会按照你的pw.x设置的hubbard site来输出,例子里面大概是这样的格式:
  1.   =-------------------------------------------------------------------=

  2.                            Hubbard U parameters:

  3.        site n.  type  label  spin  new_type  new_label  Hubbard U (eV)
  4.          1        1    Co      1      1         Co         7.7514
  5.          2        2    O       1      2         O          8.0439
  6.          3        2    O       1      2         O          8.0439

  7.   =-------------------------------------------------------------------=
复制代码

QE用hp.x计算U要简单自动化的多,不需要自己手动构造响应矩阵(即局域占据数对响应势的导数),但是类似于原始实空间下做响应使用超胞,虽然现在的方法不需要指定超胞,但需要对响应的q点尺寸进行指定,这个原则上说也应当进行收敛性测试。你要是要用下面响应矩阵手动算,应该是用U=chi0的逆矩阵-chi的逆矩阵,然后再按照对角元取,实际上这两部分也被打印到那个.dat文件当中。新的VASP应该可以通过cRPA来计算U,不过我没做过相关测试。
这两种新的算法细节我现在还没有仔细研究过,所以底层实现我还不能和你说的很清楚。

评分 Rate

参与人数
Participants 1
eV +3 收起 理由
Reason
荻花花花花 + 3 感谢您的回复

查看全部评分 View all ratings

日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

11

帖子

0

威望

215

eV
积分
226

Level 3 能力者

14#
发表于 Post on 2023-3-18 01:26:46 | 只看该作者 Only view this author
卡开发发 发表于 2022-11-11 16:21
hp.x的输出不需要去读那个矩阵,程序会输出一个叫${prefix}.Hubbard_parameters.dat的文件,会按照你的pw ...

大佬您好,看到这个输出文件,我想询问您一下,这意味着在LiCoO2体系中也需要对O原子+U吗?我记得在这个例子中也有计算V值的HUBBARD.dat文件,其中Co原子也加入了一个比较大的V值,这是什么原因呢?以及在后续计算中应该如何根据${prefix}.Hubbard_parameters.dat和HUBBARD.dat设置需要输入的U和V值呢?非常感谢

3621

帖子

3

威望

1万

eV
积分
18430

Level 6 (一方通行)

第一原理惨品小作坊

15#
发表于 Post on 2023-3-18 08:46:47 | 只看该作者 Only view this author
鱼丸box 发表于 2023-3-18 01:26
大佬您好,看到这个输出文件,我想询问您一下,这意味着在LiCoO2体系中也需要对O原子+U吗?我记得在这个 ...

简单说,就是把哪些原子作为Hubbard位点,然后考虑何种模型(+U还是+U+V)取决于进行hp.x前pw.x的计算参数,然后响应得到的参数是什么你再按照格式填回pw.x当中对应的参数进行后续计算即可。
日常打哑谜&&探寻更多可能。
原理问题不公开讨论,非商业性质讨论欢迎私聊。
本周忙

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2024-11-24 13:26 , Processed in 0.212720 second(s), 31 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list