駁網上流傳的對CP2K缺點的不實描述
駁網上流傳的對CP2K缺點的不實描述
Refuting the false description of shortcomings about CP2K circulating on the Internet
文/Sobereva@北京科音 2024-Nov-28
0 前言
CP2K是極其流行的第一性原理程序，在《2024年計算化學公社舉辦的計算化學程序和DFT泛函的流行程度投票結果》（http://www.shanxitv.org/706）的流行度得票統計中甚至已超過了VASP。筆者在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/KFP）里也積極推廣CP2K、全面介紹其相關理論背景知識和正確的使用方式。筆者之前在網上答疑時，就曾偶爾聽到一些關于CP2K的不實謠言，諸如什么算導體慢、算相互作用能不準之類。今天在群里看到有人發了一頁講CP2K缺點的幻燈片，里面列了四條，根本就沒有一條是對的，令我深感這種幻燈片實在太以訛傳訛了！估計已有很多對CP2K一無所知的人看了這種幻燈片對CP2K產生了嚴重錯誤印象，導致其放棄了其實極其適合他們的CP2K而轉向使用某些又貴又慢的程序。因此筆者忍不住寫本文，對幻燈片里的不實說法進行駁斥和斧正。實際上，如果讀者參加過北京科音CP2K第一性原理計算培訓班，正確、深入、系統學習過CP2K計算相關的理論知識，并結合培訓里的例子用CP2K做過各類計算，就自然而然就知道那幻燈片里對CP2K缺點的說法有多么離譜。我也希望讀者看到有人再給CP2K亂扣屎盆子時轉發本文鏈接。
1 對導體計算慢
那幻燈片里說CP2K“對導體計算較慢，CP2K之所以算的快是因為OT算法，對有帶隙的體系可以迅速SCF迭代收斂。但是此法原理上不適用于導體。所以導體只能用對角化方法”。
以上說法大錯特錯。首先說CP2K為什么快。對于純泛函的DFT計算，CP2K對中、大體系遠比VASP、Quantum ESPRESSO、Abinit等純粹基于平面波的程序快，這在于CP2K是原子中心基函數（具體來說是Gaussian型基函數）展開軌道波函數，同時用傳統平面波程序的方式計算經典庫侖作用。用原子中心基函數使得CP2K的得以利用矩陣的稀疏性巨幅節約算較大體系的時間，而且CP2K的開發者把程序效率視為重中之重，把CP2K的效率優化到了極致。由于相關公式較多，這里直接貼北京科音CP2K第一性原理計算培訓班講CP2K的GPW算法的兩頁ppt進行說明
此外，對于雜化泛函的計算，使用原子中心基函數的CP2K由于可以充分利用積分屏蔽技術大幅節約要計算的雙電子積分的量，使得CP2K在主流的雙路服務器上做雜化泛函計算算到上千原子都不是難事。
上述原因是CP2K做DFT特別快的最關鍵原因，而不是在于CP2K有OT。OT只不過是在這CP2K本來就特別高效的基礎上進一步令大體系每一輪的SCF耗時更低，因為OT避免了算大體系耗時很高的對角化巨大KS矩陣的步驟（體系越大時越重要，對小體系OT和對角化的差異不大）。另外，用OT時SCF往往比對角化更容易收斂。一些相關信息看《CP2K中遇到SCF難收斂時的解決方法》（http://www.shanxitv.org/665）。即便CP2K只用傳統對角化算法做SCF而不用OT，算中、大體系的速度依然吊打VASP等純平面波程序。
算金屬體系時需要開smearing，此時需要算空軌道，由于用OT做SCF過程中無法求解空軌道，因此OT一般來說不能用于金屬體系。金屬體系的KS矩陣不是很稀疏，因此前述的CP2K基于原子中心基函數帶來的優勢沒有非金屬的情況那么明顯。但即便如此，CP2K算不小的以金屬為主的體系也完全不比平面波程序慢！（這里只從gamma點的能量和受力計算來說，至于對稱性的利用與k點的考慮、SCF收斂性等其它方面另談，那是其它層面的事）
上述的北京科音CP2K第一性原理計算培訓班里的豐富的例子中就有不少是金屬為主的體系，比如“Pd(100)晶面對苯分子的吸附”、“Cu(001)表面的銀原子遷移”、“Au(111)表面H2分子解離”、“Ir(001)表面活化甲烷”等，學員親自用CP2K跑過一遍這些例子就會感受到CP2K對這些金屬為主的體系的計算效率也非常高，何來“對導體計算較慢”？
另外，前述網傳的ppt中說OT“對有帶隙的體系可以迅速SCF迭代收斂”也是錯誤說法。不少情況下OT達到收斂所需的SCF圈數反倒顯著多于對角化。OT也并不保證肯定能收斂，而且也經常遇到OT需要一百輪以上SCF才收斂的收斂緩慢的情況。OT遠沒有那么神乎其神。
2 磁性體系處理比較麻煩
那幻燈片里說CP2K“對于磁性體系處理比較麻煩，因為CP2K計算需要提前指定自旋多重度，不像VASP一樣可以自動尋找合適的磁矩”。
以上說法嚴重不實。為了保證收斂到正確的磁性狀態，CP2K需要用MAGNETIZATION指定初猜波函數中的原子自旋磁矩或者用&BS指定初始的原子電子組態，VASP也需要用MAGMOM來設定初猜波函數中的原子自旋磁矩，顯然在這點上CP2K根本沒有什么額外的更麻煩之處。
CP2K的計算根本就不是必須“提前指定自旋多重度”。北京科音CP2K第一性原理計算培訓班的下面這兩頁幻燈片里我把情況說得非常清楚：
3 k點功能不完善、小體系計算不準確
那幻燈片里說CP2K“K點功能不完善，只有用Gamma點計算比較好用，在添加的K點參數以后一些很好的加速算法就不能用了比如ASPC波函數外推，老版本的CP2K里只有Gamma點計算的功能。這意味著對大體系可以，但是小體系計算不準確”
以上說法極具誤導性。說CP2K的k點（不叫K點）的功能不完善我沒意見，確實尚有不少特性不支持k點，截止到2024版包括OT、NMR、SCCS溶劑模型、DFT+U等等，但若說成“只有用Gamma點計算比較好用”就完全是在開玩笑。CP2K早就能在純泛函計算的時候很好地考慮k點，而且效率也很好，在能量計算、（變胞）優化、CI-NEB產生反應路徑和搜索過渡態、能帶計算等等方面都表現得很好，北京科音CP2K第一性原理計算培訓班里大量例子都是開著k點做的，根本什么問題都沒有。而且如http://bbs.keinsci.com/thread-43683-1-1.html所說，CP2K從2024開始還支持了雜化泛函結合k點計算的RI-HFk算法。CP2K做DFT計算對k點考慮的一個不足（至少截止到撰文時最新的2024版）是不支持利用晶格的對稱性節約k點來降低耗時（而Quantum ESPRESSO、VASP支持），但這也不是什么大問題，因為本來CP2K就已經很快了。
CP2K支持波函數的外推，也就是在結構優化、動力學過程中用之前幾步的波函數外推出這一步的較好的初猜波函數，ASPC是其中一種。用k點時任何外推方法都沒法用，因此上面那段話提到的ASPC不支持k點不假。然而這在實際中僅僅是一個小問題，上面那段話把這問題的重要性完全夸大了。考慮k點時做幾何優化完全可以用USE_PREV_P直接使用上一步SCF收斂的密度矩陣當初猜，從而讓SCF收斂更快、降低幾何優化每一步的耗時。而分子動力學用的盒子尺寸通常不小，否則會有虛假的周期性效應導致動力學行為/現象不合理，因此此時一般都用不著考慮k點。
那幻燈片里說CP2K“小體系計算不準確”純屬無稽之談。CP2K考慮k點計算小晶胞體系不僅效率上足夠好，也不存在任何bug，何來不準確？北京科音CP2K第一性原理計算培訓班講“能量的計算及相關問題”的部分還專門給了個計算例子，令學員清楚地看到CP2K做n*n*n的超胞只考慮gamma點算的能量正好是用n*n*n k點算原胞的n*n*n倍，充分體現出CP2K考慮k點計算小體系的正確性沒有絲毫問題。
4 高斯基組帶來巨大的BSSE
那幻燈片里說CP2K用的“高斯基組帶來巨大的BSSE，基組不完備誤差，不適合計算結合能。VASP使用平面波基組就不存在這個問題”。
以上說法真是**！
BSSE問題我在《談談BSSE校正與Gaussian對它的處理》（http://www.shanxitv.org/46）有簡要介紹，北京科音中級量子化學培訓班（http://www.keinsci.com/KBQC）里“弱相互作用的計算與相關問題”一節有深入介紹。BSSE是計算弱相互作用能時候需要注意的問題，而不說具體基組談BSSE的大小根本毫無意義！！！CP2K能用的基組多了去了，最常用的是MOLOPT系列，按照此順序尺寸依次增大、BSSE依次減小：SZV-MOLOPT-GTH、DZVP-MOLOPT-SR-GTH、DZVP-MOLOPT-GTH、TZVP-MOLOPT-GTH、TZV2P-MOLOPT-GTH、TZV2PX-MOLOPT-GTH。在計算弱相互作用能方面，常用的DZVP-MOLOPT-SR-GTH有嚴重的BSSE問題，而到了TZV2P-MOLOPT-GTH檔次BSSE問題就已經相當小了。下面是MOLOPT基組原文JCP, 127, 114105 (2007)中的不同基組的對比測試，m-開頭的對應MOLOPT系列，可見TZV2P-MOLOPT-GTH算小分子相互作用能BSSE帶來的誤差僅僅有區區不到0.2 kcal/mol，遠遠小于相互作用能自身的數量級（水分子二聚體間相互作用能約5 kcal/mol），難道這能叫“巨大的BSSE”？更何況CP2K還直接支持做counterpoise校正，可以令本來就很小的BSSE問題更是完全忽略不計。雖然算大體系之間弱相互作用會比小體系間的BSSE更大，但只要用到TZV2P-MOLOPT-GTH結合counterpoise檔次依然誤差微乎其微。
有些人對CP2K使用的基組可謂近乎一無所知，就隨便網上看幾個零碎的（往往還是不適合自己情況或者有嚴重誤導性的）例子，拿比如TZVP-GTH、DZVP-MOLOPT-SR-GTH這種BSSE挺大的基組算弱相互作用能，也不知道counterpoise是何物，發現算出來結果差（一般會嚴重高估）就直接賴CP2K算不準，這些人太令人無語了。我不得不說，參加北京科音CP2K第一性原理計算培訓班把這些必備知識好好系統學學再做計算實在太重要了。以不恰當的方式做計算，不管什么程序都算不準，諸如VASP計算時用很小的平面波截斷能、Gaussian計算時用6-31G，結果能不爛么？
注意前面說的是弱相互作用能，對于涉及到成鍵的相互作用，如化學吸附的能量變化，相同的基組體現出的BSSE問題會更小，并且也不建議用counterpoise校正，見《計算化學鍵鍵能時以counterpoise方式考慮BSSE不僅是多余的甚至是有害的》（http://www.shanxitv.org/381），直接用足夠大的比如TZV2P-MOLOPT-GTH檔次做單點計算就夠了（優化和振動分析一般用DZVP-MOLOPT-SR-GTH即可）。
前述網上的幻燈片里的說法不僅嚴重誤導了CP2K的業余和潛在用戶，還令所有使用Gaussian基函數的程序中槍！諸如Gaussian、ORCA、NWChem、GAMESS-US等等等等，難道這些被廣為使用的量子化學程序算相互作用能都有巨大誤差不成！？
那幻燈片最后還來了句“VASP使用平面波基組就不存在這個問題”，顯得好像VASP、Quantum ESPRESSO平面波程序相對于CP2K就明顯很適合算相互作用似的。不僅前面說了，這種說法有嚴重誤導性，而且還忽略了純平面波程序在這方面的一個普遍劣勢，即必須用三維周期性。關于這點，挺值得一提的是在Mol. Phys., 117, 1298 (2019)文中，作者專門比較了Quantum ESPRESSO以三維周期性的方式和Q-Chem以孤立體系的方式計算的22種小分子相互作用能，結論是基于平面波計算的時候必須盒子取得很大，即便開偶極校正，也建議至少留15埃真空層，以確保令誤差小到不至于嚴重影響結果。CP2K相較于純平面波程序的一個優點是還可以以非周期性（像量子化學程序一樣）、一維周期性、二維周期性的方式做計算，因此在非周期性方向可以不考慮周期性，相應方向就可以用小得多的真空區（具體看Poission solver的選擇，在培訓班里會很詳細講），避免需要設很大的真空區浪費耗時。
值得一提的是，CP2K里有個專門做純平面波計算的SIRIUS模塊，可以得到和Quantum ESPRESSO極為相近的結果。如果你就是有特殊原因非要做平面波計算也可以用這個模塊。
5 CP2K真正有哪些缺點？
在我來看CP2K的缺點是以下這些
上面這張幻燈片里關于CP2K的k點的局限性前文已經說過一部分了，對一般應用并不構成問題。如果做的某些計算真是必須考慮k點但CP2K此時不支持的話，實際上擴胞到只需要考慮gamma點也可以解決，由于CP2K算大體系的速度非常快此時一般也都算得動。
上面這張幻燈片里說的第三方輔助工具相對有限，注意這只是目前跟VASP、Quantum ESPRESSO等歷史更長的程序相比，實際上CP2K已有大量輔助工具可以用。比如創建輸入文件可以非常便利地用Multiwfn實現，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）；Multiwfn基于CP2K的輸出文件可以做非常豐富的波函數分析以及繪制DOS和能帶圖；計算各種熱力學量CP2K可以結合Shermo，見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）；觀看振動模式可以用MfakeG+GaussView，見http://www.shanxitv.org/soft/MfakeG；計算聲子譜可以CP2K結合Phonopy；做晶體結構預測可以CPK結合USPEX；做非絕熱動力學可以CP2K結合Newton-X；ASE有CP2K的接口，等等。由于CP2K的巨大、無可取代的優點，無疑以后輔助工具還會越來越多！
上面這張幻燈片里說的算某些磁性的d族金屬表面體系SCF難收斂，典型例子就是諸如Ni(111)表面、Fe(001)表面這種，而且改各種SCF收斂相關設置也不太好解決。碰到這種情況可以用CP2K的SIRIUS模塊基于平面波計算（CP2K用戶做DFT默認用的是Quickstep模塊），或者改用免費的Quantum ESPRESSO程序。注意絕對絕對絕對、千萬千萬千萬不要以為CP2K算磁性體系就有SCF必然難收斂的問題！北京科音CP2K第一性原理計算培訓班里講了很多例子，諸如“Fe2O3的鐵磁性和反鐵磁性狀態的單點計算”、“Fe (bcc)原胞的單點計算”、“反鐵磁性的UO2的DFT+U的計算”、“順磁性物質CuCl2晶胞的優化”等等，都是順利收斂的。
上面這張幻燈片里說的沒有解析Hessian，因此只能基于解析受力做有限差分得到Hessian，是相對于Quantum ESPRESSO、VASP的一個缺點。但在實際中也不是大問題。有對稱性的小晶胞的振動分析結合Phonopy程序可以利用對稱性巨幅節約要做的有限差分的次數，再加上CP2K算其中要涉及的超胞的受力計算本來就很快，所以振動分析也花不了什么時間。而對于沒對稱性的原子較多的大晶胞，多數情況也可以通過恰當的凍結巨幅節約要做的有限差分的次數（如計算表面吸附，基底的不與被吸附分子接觸的區域都可以凍住）。
除了上述以外CP2K還有些次要缺點，例如SCCS溶劑模型下容易遇到SCF難收斂、考慮k點的變胞優化容易中途卡住（但完全可以解決，培訓里講了，計算化學公社論壇首頁google框搜也能找到我的回答）等等。
總的來說，CP2K的缺點相對于其優點是相當次要的，因此我十分推薦把極為高效、功能十分強大、特別流行還開源免費的CP2K作為第一性原理計算的默認選擇、當做主力程序，結合GaussView建模和Multiwfn產生輸入文件和做各種后處理分析（北京科音CP2K第一性原理計算培訓班里有非常全面系統的講解，零散的相關博文見http://www.shanxitv.org/category/CP2K/里的文章）用起來相當絲滑，有額外需求時再利用其它程序作為補充。
6 正確分辨網上關于CP2K的說法
最后我想再強調一下，大家切勿隨便聽信網上的可信度不明的關于CP2K的評價。比如今天看到群里有人說什么“CP2K對磁性系統太難收斂了”、“我的體系沒有磁性，也很慢”。這種說法非常不負責任。但凡是嚴謹的內行，首先就不會在沒有具體前提的情況下說這種話，而如果不是內行，又有什么能力客觀評價CP2K？像是上來就說CP2K磁性體系難收斂的，我要反問：具體是什么磁性體系（是否容易收斂都是case by case的事）？你確認初始原子磁矩合理設置了么？其它程序算這個體系就能容易收斂么、有參照物么？沒有這些最基本前提，就說“CP2K對磁性系統太難收斂了”毫無意義。至于那個說“我的體系沒有磁性，也很慢”的，我要問：慢具體是多慢、什么體系花了多少耗時？你的計算資源如何、什么CPU多少核并行？確保以合理的方式做并行計算了么（別告訴我說你用的是往往很慢的ssmp）？算的什么任務？用的什么計算級別和CUTOFF等直接影響耗時的設定？你監控SCF或幾何優化步數了么、是否發生震蕩？我這里說的這些因素都是嚴重關乎收斂性和耗時的，但凡有人沒有把這些細節準確交代出來就批CP2K，大概率其水平堪憂，因此切勿隨便聽信他們的說法，否則極大概率會誤導你。我在網上答疑時見到過太多外行人自己根本不懂怎么正確使用程序，遇到使用程序不順心的情況就直接賴程序有問題，他們應該自我檢討。