? 這里先介紹一些計(jì)算大神的想法，網(wǎng)址http://bbs.keinsci.com/thread-6299-1-1.html

具體內(nèi)容：

一般有4種方式確定U：

1.湊實(shí)驗(yàn)帶隙。

2.湊雜化泛函（如HSE）或者GW計(jì)算帶隙。

3.線性響應(yīng)方法。★★★

如果體系只有一個Hubbard site，而且盒子也比較大的情況，大體流程應(yīng)該是這樣的：

? ? (1) 計(jì)算一次SCF，然后保存該任務(wù)的電荷密度。
? ? (2) 你需要指定響應(yīng)勢α，對于vasp設(shè)置LDAUTYPE=3的情形，原來的LDAUU就變成了α。
? ? (3) 做分別施加響應(yīng)勢α（比如-0.8 0.6...-0.2 0.2 ...0.8）的計(jì)算，每次計(jì)算需要讀取無響應(yīng)情形（α=0）下的電荷密度。
? ? (4) 收集施加響應(yīng)勢后的on site占據(jù)數(shù)，對于每個α勢響應(yīng)的計(jì)算結(jié)果都有兩個部分，一個是基于α勢響應(yīng)但電荷密度未自洽的on site占據(jù)數(shù)n0，另一個是基于α勢響應(yīng)但電荷密度自洽后的on site占據(jù)數(shù)n。
? ? (5) 線性擬合α值和n0以及n的關(guān)系，你就可以得到U=X0^-1-X^-1=dα/dn0-dα/dn。

進(jìn)一步做以下討論：

? ? (1)如果體系有N個Hubbard site的話，就得依次計(jì)算第J個位點(diǎn)施加響應(yīng)以及第I位點(diǎn)的占據(jù)數(shù)，此時響應(yīng)系數(shù)dα/dn和dα/dn0都是NxN的矩陣，最后取對角化的U.
? ? (2)對于周期性計(jì)算，因?yàn)棣羷菔┘雍笠彩侵芷谛缘模c計(jì)算聲子譜類似，需要構(gòu)建超胞來消除鏡像上的影響，可以使用超胞外推的方式來使結(jié)果收斂，如PRB 71, 035105 (2005)所討論。當(dāng)然構(gòu)建超胞后的Hubbard site的數(shù)目也會成倍增加，為了縮減計(jì)算量還需要找出等價的響應(yīng)系數(shù)避免重復(fù)計(jì)算。
? ? (3)DFT直接響應(yīng)出來的結(jié)果也未必正確，實(shí)際上隨著U增加，響應(yīng)得到的Uout應(yīng)與實(shí)際輸入的Uin呈現(xiàn)一定自洽的關(guān)系，PRL 97, 103001 (2006)表明經(jīng)驗(yàn)上U較大的時候兩者構(gòu)成線性的關(guān)系。和木蟲帖子里說的不同，該方法并不是原始構(gòu)建超胞線性響應(yīng)的平行版本，使用這種方法還是得構(gòu)建超胞才能用于周期性體系。PRL 106, 118501 (2011)的SI中提供了一種別的方案，值得參考一下。這兩種自洽方法必須要求程序同時在+U的情況下進(jìn)行α勢響應(yīng)才能做，因此vasp也實(shí)現(xiàn)不了。

4? constrainted RPA

? ? 響應(yīng)系數(shù)和U之間的關(guān)系可以寫成Dyson方程的形式，所以可以采用constrainted RPA的方法進(jìn)行求解，不過我沒有專門研究過這個，你可以PRB 74, 125106(2006)，目前我不知道有哪些程序支持。

上面的不太懂也沒關(guān)系，下面我們介紹今天主要講的方法

★★★基于線性響應(yīng)近似的方法計(jì)算自洽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官方的計(jì)算方法進(jìn)行了集成，省去了對POTCAR、INCAR的修改，一鍵直接進(jìn)行U值的計(jì)算并進(jìn)行線性擬合。

具體操作如下：

1. 準(zhǔn)備必要文件

這里必要的文件是：POSCAR (你也可以自己生成KPOINTS、POTCAR、INCAR文件，也可以通過本腳本自動生成，注意自動生成INCAR后需要手動設(shè)置其中的MAGMOM值)

2. 在程序中輸入本機(jī)運(yùn)行VASP的代碼

運(yùn)行腳本后輸入0

再輸入0后，輸入服務(wù)器運(yùn)行VASP的代碼。

隨后即生成run.wsy文件，如果不想通過程序生成，也可事先自己寫好run.wsy文件。（PS：像在超算上計(jì)算的時候，可以把作業(yè)提交腳本復(fù)制成run.wsy就可以用了，但貌似是在節(jié)點(diǎn)以外運(yùn)行的。。。雖然是這樣，但勉強(qiáng)可以算）

3. 計(jì)算DFT基態(tài)

運(yùn)行程序后后輸入1，這里會要求你輸入需要計(jì)算U值原子在POSCAR中的位置以及U值加入的軌道。

例如像我算的是錳原子，就選擇2，給d軌道加U，之后程序自動創(chuàng)建文件，調(diào)用vaspkit生成新的POTCAR文件，并調(diào)VASP進(jìn)行計(jì)算?(此時會生成input.wsy文件，不要修改)

4. 進(jìn)行+U的自洽與非自洽計(jì)算

運(yùn)行腳本后輸入2，程序自動調(diào)VASP進(jìn)行自洽與非自洽計(jì)算。

5.?計(jì)算U值

上述計(jì)算全部完成之后，運(yùn)行腳本后輸入3，計(jì)算U值，并生成output.wsy文件

打開output.wsy，最后一行即為該原子所對應(yīng)的U值。

也可如VASP網(wǎng)中通過上述的第一列和最后兩列數(shù)據(jù)進(jìn)行線性擬合得到。

(由于贗勢選取的不同，對于NiO計(jì)算結(jié)果與官網(wǎng)有些許差別)

官網(wǎng)鏈接：https://www.vasp.at/wiki/index.php/Calculate_U_for_LSDA%2BU

官網(wǎng)給出的結(jié)果：

這是腳本初代版本，后續(xù)還會加入新功能，大家喜歡的話可以關(guān)注一下。

最新程序下載地址：https://github.com/Code-WSY/Code-WSY

up在使用過程中遇到了兩個問題，給大家講一下

? 1.在INCAR中不要設(shè)置ISTART和ICHARG參數(shù)，設(shè)置了之后，腳本無法在后續(xù)計(jì)算中對其修改，進(jìn)而不能完成自洽和非自洽計(jì)算。如果不會寫INCAR的話，就讓它自動生成。它生成的腳本和官網(wǎng)例子是一樣的。

? 2.用不同的K點(diǎn)會得到不同的U值結(jié)果，我用111K點(diǎn)和448K點(diǎn)算β-MnO2晶胞給出的U結(jié)果是不同的

? K采用448算出來的算是比較接近的，K點(diǎn)的影響會有多大呢，我也不清楚，希望評論區(qū)有大神可以指點(diǎn)一下。

? ? 3.在VASPKIT群里咨詢過一個問題，“大家看下這個帖子，腳本里的LDYUTYPE=3這句是不是寫錯了”，因?yàn)閂ASP說明書給出的選項(xiàng)只有1，2，4. 然后大神給出的回答是“沒問題。3是內(nèi)置線性響應(yīng)算U值用的，官網(wǎng)教程也是3。算是一個小trick ”.

? 同樣的，除了線性響應(yīng)的方法計(jì)算U值，還可以在文獻(xiàn)中查找與自己相關(guān)體系采用的U值，其中，up主計(jì)算老司機(jī)有一個視頻分享了常見加U原子的U值視頻，詳細(xì)內(nèi)容可以看他的視頻DFT+U的U值哪里來

? ? 4.最后想說，無論是哪種方法得到的U值，都要對你關(guān)注的性能能夠很好的描述才能去使用

? ? ?整理不易，希望大家可以對我的工作關(guān)注，點(diǎn)贊，投幣，您的關(guān)注是我不懈更新的動力！??！