VASP基于线性响应近似的方法计算DFT+U的U值（有脚本，三分钟学会）
这里先介绍一些计算大神的想法，网址http://bbs.keinsci.com/thread-6299-1-1.html
具体内容：
一般有4种方式确定U：
1.凑实验带隙。
2.凑杂化泛函（如HSE）或者GW计算带隙。
3.线性响应方法。★★★
如果体系只有一个Hubbard site，而且盒子也比较大的情况，大体流程应该是这样的：
(1) 计算一次SCF，然后保存该任务的电荷密度。
(2) 你需要指定响应势α，对于vasp设置LDAUTYPE=3的情形，原来的LDAUU就变成了α。
(3) 做分别施加响应势α（比如-0.8 0.6...-0.2 0.2 ...0.8）的计算，每次计算需要读取无响应情形（α=0）下的电荷密度。
(4) 收集施加响应势后的on site占据数，对于每个α势响应的计算结果都有两个部分，一个是基于α势响应但电荷密度未自洽的on site占据数n0，另一个是基于α势响应但电荷密度自洽后的on site占据数n。
(5) 线性拟合α值和n0以及n的关系，你就可以得到U=X0^-1-X^-1=dα/dn0-dα/dn。
进一步做以下讨论：
(1)如果体系有N个Hubbard site的话，就得依次计算第J个位点施加响应以及第I位点的占据数，此时响应系数dα/dn和dα/dn0都是NxN的矩阵，最后取对角化的U.
(2)对于周期性计算，因为α势施加后也是周期性的，与计算声子谱类似，需要构建超胞来消除镜像上的影响，可以使用超胞外推的方式来使结果收敛，如PRB 71, 035105 (2005)所讨论。当然构建超胞后的Hubbard site的数目也会成倍增加，为了缩减计算量还需要找出等价的响应系数避免重复计算。
(3)DFT直接响应出来的结果也未必正确，实际上随着U增加，响应得到的Uout应与实际输入的Uin呈现一定自洽的关系，PRL 97, 103001 (2006)表明经验上U较大的时候两者构成线性的关系。和木虫帖子里说的不同，该方法并不是原始构建超胞线性响应的平行版本，使用这种方法还是得构建超胞才能用于周期性体系。PRL 106, 118501 (2011)的SI中提供了一种别的方案，值得参考一下。这两种自洽方法必须要求程序同时在+U的情况下进行α势响应才能做，因此vasp也实现不了。
4 constrainted RPA
响应系数和U之间的关系可以写成Dyson方程的形式，所以可以采用constrainted RPA的方法进行求解，不过我没有专门研究过这个，你可以PRB 74, 125106(2006)，目前我不知道有哪些程序支持。
上面的不太懂也没关系，下面我们介绍今天主要讲的方法
★★★基于线性响应近似的方法计算自洽Hubbard U值的方法
帖子原址：http://mp.weixin.qq.com/s?__biz=MzIwNTQzMTk5Mw==&mid=2247484324&idx=1&sn=bfe1474f2b332abbab95a9d40ff1e9d4&chksm=9731b672a0463f64207b4cc6c1cb65311882fb45edd82c917d346aeede579f4945be8e055f06&mpshare=1&scene=23&srcid=0328xEWuhZa5Hk7KnMj3lCd0&sharer_sharetime=1648453222322&sharer_shareid=44721a85a1edb3657384cdcc5f1c5fe4#rd
要使用到知乎和微信公众号苏理士多云分享的脚本，具体也是将VASP官方的计算方法进行了集成，省去了对POTCAR、INCAR的修改，一键直接进行U值的计算并进行线性拟合。
具体操作如下：
1. 准备必要文件
这里必要的文件是：POSCAR (你也可以自己生成KPOINTS、POTCAR、INCAR文件，也可以通过本脚本自动生成，注意自动生成INCAR后需要手动设置其中的MAGMOM值)
2. 在程序中输入本机运行VASP的代码
运行脚本后输入0
再输入0后，输入服务器运行VASP的代码。
随后即生成run.wsy文件，如果不想通过程序生成，也可事先自己写好run.wsy文件。（PS：像在超算上计算的时候，可以把作业提交脚本复制成run.wsy就可以用了，但貌似是在节点以外运行的。。。虽然是这样，但勉强可以算）
3. 计算DFT基态
运行程序后后输入1，这里会要求你输入需要计算U值原子在POSCAR中的位置以及U值加入的轨道。
例如像我算的是锰原子，就选择2，给d轨道加U，之后程序自动创建文件，调用vaspkit生成新的POTCAR文件，并调VASP进行计算 (此时会生成input.wsy文件，不要修改)
4. 进行+U的自洽与非自洽计算
运行脚本后输入2，程序自动调VASP进行自洽与非自洽计算。
5. 计算U值
上述计算全部完成之后，运行脚本后输入3，计算U值，并生成output.wsy文件
打开output.wsy，最后一行即为该原子所对应的U值。
也可如VASP网中通过上述的第一列和最后两列数据进行线性拟合得到。
(由于赝势选取的不同，对于NiO计算结果与官网有些许差别)
官网链接：https://www.vasp.at/wiki/index.php/Calculate_U_for_LSDA%2BU
官网给出的结果：
这是脚本初代版本，后续还会加入新功能，大家喜欢的话可以关注一下。
最新程序下载地址：https://github.com/Code-WSY/Code-WSY
up在使用过程中遇到了两个问题，给大家讲一下
1.在INCAR中不要设置ISTART和ICHARG参数，设置了之后，脚本无法在后续计算中对其修改，进而不能完成自洽和非自洽计算。如果不会写INCAR的话，就让它自动生成。它生成的脚本和官网例子是一样的。
2.用不同的K点会得到不同的U值结果，我用111K点和448K点算β-MnO2晶胞给出的U结果是不同的
K采用448算出来的算是比较接近的，K点的影响会有多大呢，我也不清楚，希望评论区有大神可以指点一下。
3.在VASPKIT群里咨询过一个问题，“大家看下这个帖子，脚本里的LDYUTYPE=3这句是不是写错了”，因为VASP说明书给出的选项只有1，2，4. 然后大神给出的回答是“没问题。3是内置线性响应算U值用的，官网教程也是3。算是一个小trick ”.
同样的，除了线性响应的方法计算U值，还可以在文献中查找与自己相关体系采用的U值，其中，up主计算老司机有一个视频分享了常见加U原子的U值视频，详细内容可以看他的视频DFT+U的U值哪里来
4.最后想说，无论是哪种方法得到的U值，都要对你关注的性能能够很好的描述才能去使用
整理不易，希望大家可以对我的工作关注，点赞，投币，您的关注是我不懈更新的动力！！！
