CP2K中遇到SCF難收斂時的解決方法
CP2K中遇到SCF難收斂時的解決方法
文/Sobereva@北京科音 2023-May-15
0 前言
經常有人問CP2K中遇到SCF難收斂怎么解決，這是很家常便飯的問題。這里就基于筆者講授的“北京科音CP2K第一性原理培訓班”（http://www.keinsci.com/workshop/KFP_content.html）里面相關部分的內容專門寫個博文，供CP2K用戶參考。這里假定讀者已經了解CP2K計算的基本的背景知識（如k點、CUTOFF、對角化、OT、&BS、smearing等等），如果不了解相關名詞、算法的話，十分推薦參加這個培訓全面、系統學習一遍。
筆者之前專門寫過Gaussian用戶遇到SCF不收斂的解決方法的文章《解決SCF不收斂問題的方法》（http://www.shanxitv.org/61），本文和那篇文章一樣，已經把所有可能有用的方法和要注意的問題說得盡可能全面、充分了，讀者遇到難收斂/不收斂時需要結合實際情況一個個嘗試。切勿問筆者“這文章我看了，但還是解決不了”之類的，我沒有任何可額外補充的（至多是給出當前計算的盡可能詳細的信息，比如給出輸入輸出文件，然后我能回答哪些問題需要優先考慮、哪些方法值得優先嘗試）。
下面第1節先說說遇到SCF不收斂時一般情況需要考慮的問題和解決策略，然后第2、3節再分別說使用對角化和OT時分別可以進一步考慮的。最后再說一下幾何優化、分子動力學時與SCF收斂有關的問題。在此之前先看培訓班ppt了解一下CP2K是怎么判斷SCF收斂的：
1 一般情況時考慮的
? 檢查結構合理性。SCF是否容易收斂和結構的合理性關系極大，化學意義越強的結構通常越容易收斂。因此遇到SCF很難收斂時，需要結合結構化學、固體物理常識檢查結構是否靠譜，如模型是否有意義、是否漏了原子、是否某些原子所處的化學環境和實際明顯不符（如過渡金屬本該與配體成鍵的地方卻光禿禿地懸著）等等。且要結合圖像肉眼檢查，確保不存在不該出現的結構嚴重扭曲、盒子邊緣原子和周期鏡像原子發生過近接觸等問題。若有這些問題，即便SCF收斂了，結果顯然也毫無意義。
? 確認關鍵詞、設置沒明顯不合理的。CP2K的關鍵詞很復雜，對使用者的知識水平要求很高，若用了明顯不當的關鍵詞和設置，SCF很難收斂或收斂到離譜結果是極其正常的事，然而很多要點是初學者難以注意到的。比如計算自旋極化體系的時候必須得在&DFT里寫UKS（或等價的LSD），否則即便自旋多重度設的是>1，程序也會愣當成是閉殼層算。再比如用wavelet的Poisson solver算孤立體系的時候，如果體系有明顯電子密度分布的區域跨越了盒子，結果就無意義。
CP2K屬于那種是完全手動擋還很個性的車。如果經驗不夠豐富的話，筆者建議總是用Multiwfn（http://www.shanxitv.org/multiwfn）產生其輸入文件，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）。Multiwfn的CP2K輸入文件創建功能將CP2K盡最大可能包裝成黑箱，不僅用起來極其方便省事，對初學者還能大大減小關鍵詞使用不當的可能性。比如將自旋多重度設成>1時自動就會添加UKS、用wavelet做非周期性計算時會自動設置足夠大又不太浪費的盒子尺寸并且加入自動讓體系在盒子里居中的設置。
? 確保體系的凈電荷、自旋多重度的設置合理。這倆如果設置不當的話，相當于計算了無意義或激發態的電子結構，SCF極難收斂也是很可能出現的事。特別是對于磁性體系，自旋多重度的設置對SCF收斂往往影響極大。
? 如果SCF有收斂趨勢，可增大SCF迭代次數上限再試，即&SCF里的MAX_SCF，其默認的50對于很多情況明顯太小。Multiwfn產生的CP2K輸入文件里默認用128（128是為了優雅，2的7次方。和Gaussian程序一致），絕大多數情況都夠了，但也可能碰上極個別體系仍需要加大。
筆者超級反感一碰到SCF不收斂，也不觀察收斂趨勢，上來就增加迭代次數上限的做法，因為這是最典型的菜鳥思維，絕大多數情況毫無用處而白浪費時間。因為如果SCF過程都發生顯著震蕩了、沒收斂的苗頭，再怎么迭代下去都是白搭，純粹是毫無意義地浪費時間和計算資源。
? 如果發現用其它泛函或基組下計算能收斂的話，可以讀取其收斂的波函數作為初猜，這至少比自動產生的初猜波函數要好、更容易收斂。也即&SCF里用SCF_GUESS RESTART，且&DFT里WFN_RESTART_FILE_NAME指定成之前任務的.wfn文件（若考慮了k點是.kp文件）。注意兩次的基組特征必須定性一致，比如之前計算用的大核贗勢，而現在用小核贗勢或全電子基組，這時候讀之前的波函數就沒意義了。
通常較小的基組下計算比大基組下計算更容易收斂，比如直接用TZV2P檔次的MOLOPT基組不收斂的話，可以嘗試用DZVP檔次的算，若收斂了，讀它收斂的波函數為初猜。
往往HF成份越大的泛函越容易收斂，很大程度是因為其HOMO-LUMO gap更大、SCF過程中占據和空軌道混合程度相對更緩和所致。因此遇上SCF難收斂時可以嘗試用HF成份比實際要用的泛函更大的泛函，若能收斂，讀取其收斂的波函數當初猜。HF成份可以參考《不同DFT泛函的HF成份一覽》（http://www.shanxitv.org/282）。
如果原本用的是純泛函，想嘗試雜化泛函但用不動，也可以用DFT+U試試，耗時和純泛函基本一致，但此時不支持k點，此外怎么獲得+U用的Ueff參數也是個問題。如果實際計算不打算+U，也可以嘗試3或4 eV的+U看看能否令SCF收斂，如果能，則不+U時讀取其收斂的波函數當初猜。
? 做雜化泛函計算時（特別是對于周期性體系）通常會在&DFT / &XC / &HF / &SCREENING里開啟SCREEN_ON_INITIAL_P從而基于初始的密度矩陣做雙電子積分的屏蔽來巨幅節約耗時，Multiwfn產生的雜化泛函計算的輸入文件里默認也開了這個。若此時沒有讀取之前收斂的波函數當初猜（通常是便宜得多的純泛函計算的），則CP2K會基于自動初猜的波函數對應的密度矩陣做積分屏蔽，此時結果通常無意義，SCF能收斂的概率也非常低！然而這個雜化泛函計算的關鍵要點卻經常被初學者所忽略！
? 確保恰當考慮了k點、k點數夠多。如果晶胞邊長不夠大，顯然必須考慮k點而不能只用默認的Gamma點。而且k點設得還得夠多，否則SCF難收斂極為正常。k點如果設置不當，就算SCF能收斂了也沒意義。眾所周知，算能量問題，原則上k點應當大到足以令能量隨k點數增加而基本收斂的程度。
? 嘗試增大&DFT / &MGRID里的CUTOFF、REL_CUTOFF。這倆截斷能設置控制格點精度，不夠大的話不僅結果不準，甚至可能造成SCF難收斂。如果當前明顯已經夠大了，就沒必要再盲目嘗試再增加了。注：還有種說法是嘗試&DFT / &XC / &XC_GRID里寫USE_FINER_GRID T將計算交換-相關泛函部分的格點精度提升為CUTOFF的四倍，但根據我的經驗這并沒有什么實際用處，而耗時暴增。
? CP2K默認的產生初猜波函數的方法是孤立原子密度的疊加，在整個任務一開始時會對每個KIND對應的原子在孤立狀態下做計算，用PRINT_LEVEL MEDIUM時可以看到具體細節。計算孤立原子所用的電子組態可以通過&KIND里的MAGNETIZATION（定義alpha電子數減beta電子數）或者&BS（可以具體控制各個主量子數的各個角量子數上的alpha和beta電子數）來設置。初猜波函數中的孤立原子的組態越接近整個體系收斂的波函數中原子的實際組態則SCF越容易收斂、SCF不收斂的可能性越低。
如果體系有磁性，遇到SCF不收斂時一定別忘了嘗試用MAGNETIZATION或&BS盡可能令初猜波函數中原子組態與實際情況接近，這對SCF是否能夠順利收斂有關鍵性影響。如果當前體系在文獻（可用google搜索來找）或高通量數據庫（如Material Project）中已經被計算過，可以參考那里面報道的原子自旋磁矩或原子組態來設。有時同一種元素在體系里處于不止一種化學環境，差異很大時（比如明顯處于不同價態）可分成多個&KIND分別設置。
即便體系沒有磁性，利用&BS使初猜波函數中的原子帶電狀態使之盡可能符合實際也往往是有益的。比如對于NaCl這樣高度離子性體系計算時SCF不收斂的話，可以嘗試令初猜中Na和Cl的原子組態分別對應Na+和Cl-的狀態（默認情況下分別對應中性的Na和Cl）。曾經我用PBEsol/pob-TZVP-rev2計算10層的NaCl板時通過這種做法成功解決了SCF不收斂問題。
? 對gap較小的體系可嘗試使用能級移動，即&SCF / LEVEL_SHIFT，設比如0.3、0.4左右，默認單位是Hartree。這使得SCF過程中人為加大空軌道與占據軌道之間的能級差，從而減輕彼此間混合的劇烈程度而有時令SCF更容易收斂。
? 嘗試將&QS / EPS_DEFAULT減小幾個數量級，默認為1E-10。EPS_DEFAULT這個參數控制一堆涉及數值計算精度的參數，如果設得太大，有可能造成SCF難收斂。若這個數值當前已經很小了，比如都1E-12了，嘗試繼續減小則沒必要。
? 改用其它泛函、基組。如果原先的計算級別就是極難收斂，而有其它級別性價比、精度相仿佛，且發現能收斂，那就別死磕當前的了。
? 在GPW與GAPW之間切換。如果之前用的GPW結合GTH贗勢基組不收斂，嘗試用GPAW結合一個全電子基組計算試試，反之亦然。
? 修改默認產生波函數的方式。參見手冊里&SCF / SCF_GUESS設置。但其實這沒什么用，除了默認的ATOMIC（即前述的基于孤立原子密度疊加產生的）以及RESTART（讀取之前的波函數當初猜），其它的一個能打的都沒有。
用ATOMIC產生初猜時，CP2K計算每個KIND對應的孤立原子也是需要SCF迭代的，上限是固定的50輪（沒法靠關鍵詞改），在PRINT_LEVEL MEDIUM時其迭代過程會在一開始輸出。有時候我發現到了50輪時收斂精度還很差，故有可能在某些情況下由于這個問題導致初猜質量爛，不利于實際體系的SCF收斂。我發現可以通過atom_kind_orbitals.F里atom%optimization%max_iter修改迭代次數上限，比如改成300然后重新編譯，一般孤立原子的SCF就都能收斂了。
? 如果加了外電場時SCF很難收斂，可能是外電場過大所致，可嘗試更小的數值。如果就是要用較大的外電場，可以試試讀取較小電場或不帶電場時收斂的波函數當初猜。有的體系在很強外電場下本來就不可能穩定存在，這時候SCF不收斂自然很正常，是物理上決定的，也別盲目去計算。
? 實在沒轍時干脆放寬SCF收斂限，即加大&SCF / EPS_SCF，默認是1E-5（這其實已經挺寬松了），這是走投無路時的做法。或者換程序，如免費且和CP2K同樣流行的Quantum ESPRESSO，或者嘗試CP2K里的純粹基于平面波的SIRIUS模塊。某些體系用CP2K的quickstep模塊就是很難收斂，如某些d族純金屬的slab，這和算法、基函數特征有關。
2 用對角化的情況可專門考慮的
? 對無gap（如金屬、石墨烯）或gap很小的體系（如非常窄帶半導體）使用smearing。即加入&SCF / &SMEAR / METHOD FERMI_DIRAC，且在&SCF里用ADDED_MOS要求求解一些空軌道從而能被熱激發的電子所占據。ADDED_MOS數目夠用就行，設太大了會增加耗時而且也沒用處。Multiwfn產生CP2K輸入文件時如果開了smearing會自動做ADDED_MOS設置，并根據體系原子數估計一個通常夠用的值。用較小電子溫度時若不收斂可嘗試再提升溫度，即&SMEAR / ELECTRONIC_TEMPERATURE設的。溫度設高時ADDED_MOS往往需要相應地增加，否則熱激發的電子沒法按Fermi-Dirac分布占該占的空軌道（此時計算時會有警告）。
? 修改&SCF / &MIXING / METHOD設置SCF過程中將舊密度矩陣與新密度矩陣混合得到實際每一輪密度矩陣的方法。默認的DIRECT_P_MIXING往往表現極差。通常用BROYDEN_MIXING，這也是Multiwfn產生的輸入文件里默認的。亦可嘗試PULAY_MIXING，但一般不比BROYDEN_MIXING好。
? 對于使用BROYDEN_MIXING的情況，嘗試修改&MIXING中的ALPHA（新密度矩陣混入舊的的比例，一般用默認的0.4，即混入40%。難收斂時可嘗試諸如0.1、0.2、0.3）和NBROYDEN（默認的4明顯太小，改大到8通常很有助于收斂，這是Multiwfn產生的輸入文件默認的。也可以嘗試進一步改大到比如12）
? 對于使用PULAY_MIXING的情況，嘗試修改&MIXING中的NPULAY （默認的4偏小，改大往往很有助于收斂）。注：NPULAY、NBROYDEN、NBUFFER三個關鍵詞實際上等價。
? 對于使用DIRECT_P_MIXING的情況，收斂到&SCF / EPS_DIIS值后會切換到加速收斂的DIIS算法。但默認的EPS_DIIS為0.1，有點偏小，可以嘗試改大以早點啟用DIIS（但借助DIIS通常仍不如BROYDEN_MIXING）。還可以把控制最大DIIS矢量的&SCF / MAX_DIIS從默認的4改大到諸如7。
? 在對角化框架下各種方法都不行的話，嘗試改用OT。如果OT能收斂的話，而且當前由于特殊目的又必須用對角化，也可以用OT收斂的波函數當對角化計算的初猜。由于OT不支持k點，因此如果之前對角化是結合k點算的較小晶胞情況，用OT時就只能靠超胞模型了。OT也不支持用ADDED_MOS設置在SCF計算過程中求解空軌道，因此必須用smearing的導體體系就別指望OT了。如果是用GFN1-xTB，總應當用OT，這也是Multiwfn默認的。
3 用OT的情況可專門考慮的
? 如果之前沒用outer SCF則開啟之，這對解決SCF不收斂非常有效，也即加入&SCF / &OUTER_SCF字段。由于其重要性，Multiwfn產生的輸入文件里開OT后默認就是啟用outer SCF的。此時每次到達&SCF / MAX_SCF所設的inner SCF的步數上限時，就會進入outer SCF過程更新OT的preconditioner，然后進入下一次的inner SCF迭代過程。每次更新preconditioner后經常會看到比之前以更快的速度趨近收斂限。原理上來說，OT每一步都更新preconditioner是最理想、SCF最容易收斂的，但由于這樣太昂貴，因此默認情況下（即不用outer SCF時）只會在最開始計算一次preconditioner。
開啟outer SCF時inner SCF的迭代次數上限應當設得比平時更小（否則經歷很多次SCF迭代才做一次outer SCF，無法有效達到用outer SCF的目的），15至35是比較合適的范圍。outer SCF迭代次數上限通過&SCF / &OUTER_SCF / MAX_SCF控制，總SCF迭代次數上限相當于inner SCF次數上限和(outer SCF次數上限+1)的乘積。
&SCF / EPS_SCF設的是inner SCF收斂限，&SCF / &OUTER_SCF / EPS_SCF設的是outer SCF收斂限，后者應當小于或等于前者（一般設為相同）。
? 嘗試其它的preconditioner。FULL_ALL通常最好，但最貴。大體系可以嘗試明顯更便宜的FULL_SINGLE_INVERSE和FULL_KINETIC。個別時候FULL_KINETIC反倒比FULL_ALL更容易收斂。GAPW計算總應當用FULL_ALL。
? 修改OT算法，把&OT里的ALGORITHM從默認的STRICT改為IRAC，SCF經常會收斂得更穩健，非常值得嘗試！特別是我發現用pob系列全電子基組的時候，用IRAC往往容易收斂得多，因此Multiwfn產生基于pob系列基組的輸入文件時若用戶選了OT，則默認用IRAC。而且用約束性DFT（CDFT）時，我發現用IRAC也會令SCF容易收斂得多得多。
? &OT / MINIMIZER控制OT用的極小化方法，若之前用的是較快的DIIS（也是Multiwfn產生的用OT的輸入文件里默認的），可改為更穩健但耗時更高的CG（共軛梯度法），也可以嘗試BROYDEN。
? MINIMIZER用CG時，嘗試用&OT / LINESEARCH選項把線搜索方式設為比默認2PNT更貴但也更穩的3PNT。
4 幾何優化和分子動力學過程中的SCF收斂問題
在幾何優化過程中，有可能中途偶爾出現幾步SCF不收斂。此時程序還會繼續算下去。只要最終幾何優化收斂時SCF也是收斂的狀態就行了。
有時候在幾何優化初期SCF一直不收斂，不要急于把任務停掉，可以先湊合跑跑看看情況。初始幾何結構相對于實際的極小點結構而言肯定是存在多多少少不合理性的，此時SCF也肯定比極小點結構處更難收斂。在SCF不收斂、受力也較糙的情況下湊合跑一定步數后，結構一般會變得比一開始更為合理，此時SCF也往往變得更容易收斂。
在幾何優化、分子動力學過程中，程序會自動利用之前一些結構的波函數外推出當前步的初猜波函數，這比起直接產生的初猜波函數質量通常好得多。外推方式可以通過&DFT / &QS / EXTRAPOLATION來具體設置，默認的ASPC通常是最好的選擇。由于這一點，當幾何優化、分子動力學跑一定步數之后，SCF達到收斂限所需的迭代次數明顯減小、出現SCF不收斂的可能性也大為降低。
注意當使用k點時CP2K不支持外推產生初猜波函數，因此會在幾何優化、分子動力學等涉及結構變化的任務中每一步重新產生初猜波函數，即EXTRAPOLATION USE_GUESS的情況。我一般建議改為EXTRAPOLATION USE_PREV_P，代表直接使用上一步SCF收斂的波函數當初猜。
如果想減小MD中途出現SCF不收斂的可能，可以嘗試減小步長（&MD / TIMESTEP）。如果想減小幾何優化中途出現SCF不收斂的可能，可減小優化的步長上限，即置信半徑（如通過&MOTION / &GEO_OPT / &BFGS / TRUST_RADIUS）。這樣做之后由于相鄰兩個結構間差異變小了，ASPC或USE_PREV_P給出的初始波函數對于當前結構而言就越理想，自然SCF收斂也會更容易，且有利于讓SCF能收斂的狀態保持下去。