CP2K教程系列之靜態計算(Pymatflow篇)

本系列CP2K教程是《CP2K菜根譚》的升級版，在舊版基礎上添加了如何結合Pymatflow工具簡化計算流程的內容。話不多說，本文將為您帶來CP2K系列教程中的靜態計算部分。

靜態計算設置

靜態計算是大多數人接觸第一性原理計算後第一次運行的計算類型。很多其它類型的計算的參數設置與靜態計算參數的設置緊密相關，比如結構優化，AIMD，過渡態搜索等。

只有靜態電子結構計算的參數設置合理才能夠保證後續的結構優化等步驟有一個堅實的基礎。比如有時候，很急衝衝地就想要對體系進行結構優化，但是忽略了電子步DFT參數的設置的合理性，就有可能無法得到好的優化結果。

可見掌握良好合理的單點能計算參數的設置對於進行後續的各種類型的模擬也是至關重要的。

單點能計算的核心也就是進行Kohn-Sham方程的一個自洽求解過程，通常又被稱為自洽場(SCF)過程。在CP2K中，這部分的設置主要通過FORCE_EVAL輸入模塊進行。

完整輸入文件

老傳統，先給出完整輸入文件，再一一進行解釋。

<code>&GLOBAL
    PROJECT NaCl-Static
    RUN_TYPE ENERGY_FORCE
    PRINT_LEVEL LOW
    FFTW_PLAN_TYPE ESTIMATE
&END GLOBAL
​
&FORCE_EVAL
    METHOD QS
    &SUBSYS
        &KIND Na
            BASIS_SET DZVP-MOLOPT-SR-GTH
            POTENTIAL GTH-PBE
        &END KIND
        &KIND Cl
            BASIS_SET DZVP-MOLOPT-SR-GTH
            POTENTIAL GTH-PBE
        &END KIND
        &CELL
            A 5.691694000 0.000000000 0.000000000
            B 0.000000000 5.691694000 0.000000000
            C 0.000000000 0.000000000 5.691694000
            PERIODIC XYZ
        &END CELL
        &TOPOLOGY
            COORD_FILE_FORMAT xyz
            COORD_FILE_NAME nacl.xyz
        &END TOPOLOGY
    &END SUBSYS
    &DFT
        BASIS_SET_FILE_NAME BASIS_MOLOPT
        POTENTIAL_FILE_NAME GTH_POTENTIALS
        &QS
            METHOD gpw
            EPS_DEFAULT 1e-14
        &END QS
        &POISSON
            PERIODIC XYZ
            POISSON_SOLVER PERIODIC
        &END POISSON
        &MGRID
            CUTOFF 100
            REL_CUTOFF 60
            NGRIDS 4
        &END MGRID
        &XC 

            &XC_FUNCTIONAL pbe
            &END XC_FUNCTIONAL
        &END XC
        &KPOINTS
            SCHEME monkhorst-pack 3 3 3
        &END KPOINTS
        &SCF
            ADDED_MOS 0
            SCF_GUESS ATOMIC
            EPS_SCF 1e-06
            MAX_SCF 50
            &DIAGONALIZATION TRUE
                ALGORITHM STANDARD
            &END DIAGONALIZATION
            &MIXING TRUE
                ALPHA 0.4
                METHOD BROYDEN_MIXING
            &END MIXING
        &END SCF
        &PRINT
        &END PRINT
    &END DFT
    &PROPERTIES
    &END PROPERTIES
&END FORCE_EVAL/<code>

其中nacl.xyz文件的內容為

<code>8
a comment line
Na  0.000000    0.000000    0.000000
Na  0.000000    2.845847    2.845847
Na  2.845847    0.000000    2.845847
Na  2.845847    2.845847    0.000000
Cl  0.000000    0.000000    2.845847
Cl  0.000000    2.845847    0.000000
Cl  2.845847    0.000000    0.000000
Cl  2.845847    2.845847    2.845847/<code>

(將輸入參數文件與座標文件nacl.xyz放在同一個路徑下就可以運行示例)

整個輸入文件分為兩個大的部分，即&GLOBAL與&FORCE_EVAL。其中GLOBAL比較簡單，參數PROJECT用於指定項目名，其會被用來輔助cp2k輸出文件的命名，這裡不用過多關心。我們需要關心的是設置運行類型，即RUN_TYPE為ENERGY_FORCE，或者也可以設置為ENERGY，兩者都是靜態計算，不過後者不會計算原子受到的力。另外PRINT_LEVEL決定了程序輸出內容的多少，可以為DEBUG、HIGH、LOW、MEDIUM、SILENT五個選項之一。

好了，該輪到重頭戲&FORCE_EVAL部分了。下面我們按照出現的順序依次進行解釋。

首先就是&FORCE_EVAL模塊的參數，常見的設置就是QS或者叫QUICKSTEP了。Quickstep是CP2K的DFT計算的核心方法實現。Quickstep基於高斯平面波(GPW)或高斯綴加平面波(GAPW)基組。關於這兩種基組的優點可以見文章CP2K 菜根譚之簡介。

然後就需要通過&SUBSYS模塊來設置計算體系。通過&KIND &END_KIND組合來為體系中所有的類型的元素定義DFT計算中使用的贋勢和基組類型，你可以進入你的CP2K安裝目錄下的data中查看贋勢和基組文件。需要注意的是這裡通過BASIS_SET和POTENTIAL定義的基組和贋勢需要分別出現在後面的通過&DFT中的BASIS_SET_FILE_NAME和POTENTIAL_FILE_NAME定義的基組文件和贋勢文件中。

&CELL是用來定義超胞的，而&TOPOLOG被用於定義體系的座標幾何。通過COORD_FILE_FORMAT定義座標文件的格式，如這裡的xyz格式，通過COORD_FILE_NAME定義座標文件名。得到支持的文件格式包括: cif、crd、g96、pdb、xyz等常見格式。另外除了通過這種引入外部文件的方式來定義體系座標外，也可以直接通過COORD來進行座標的定義。比如可以將上面的&TOPOLOGY部分替換為下面的內容:

<code>&COORD
Na  0.000000    0.000000    0.000000
Na  0.000000    2.845847    2.845847
Na  2.845847    0.000000    2.845847
Na  2.845847    2.845847    0.000000 

Cl  0.000000    0.000000    2.845847
Cl  0.000000    2.845847    0.000000
Cl  2.845847    0.000000    0.000000
Cl  2.845847    2.845847    2.845847
&END COORD/<code>

解析來就是定義DFT計算的一系列參數，其中BASIS_SET_FILE_NAME``和POTENTIAL_FILE_NAME在前面已經有提到過，用來定義基組和贋勢文件名。

<code>&QS
METHOD GPW
EPS_DEFAULT 1.0E-10
&END QS/<code>

這部分定義QS方法的一些參數，包括使用的方法，如這裡的GPW。以及為所有的能量收斂閾值設置一個默認值EPS_DEFAULT。

<code>&MGRID
  CUTOFF 100
  REL_CUTOFF 60
  NGRIDS 4
&END MGRID
​/<code>

MGRID定義了實空間mutltigrid的一些參數，包括實空間網格的截斷值CUTOFF(可以類比為截斷能)。使用到的multigrid的數量NGRIDS

<code>&XC
  &XC_FUNCTIONAL PBE
  &END XC_FUNCTIONAL
&END XC/<code>

這部分定義了使用的交換關聯泛函，如這裡使用的是PBE泛函。

<code>&SCF
  ADDED_MOS 0
  SCF_GUESS ATOMIC
  EPS_SCF 1e-05 

  MAX_SCF 50
  &DIAGONALIZATION TRUE
  ALGORITHM STANDARD
  &END DIAGONALIZATION
  &MIXING TRUE
      ALPHA 0.4
      METHOD BROYDEN_MIXING
  &END MIXING
  &PRINT
  &END PRINT
&END SCF/<code>

這裡定義了SCF自洽過程的相關參數。如ADDED_MOS定義了額外添加到計算中的軌道的數量。SCF_GUESS定義了初始波函數的構造方案。EPS_SCF定義了SCF的收斂閾值。MAX_SCF定義了最大的自洽循環步數。

&DIAGONALIZATION部分表明使用常規的正交化求解KS矩陣。至於CP2K中替代傳統正交化方法的OT方法，我們將在後面的文章中進行介紹。

&MIXING部分定義了SCF迭代過程中對密度矩陣進行混合的參數。ALPHA定義了密度矩陣混合的比例，METHOD定義了採取的混合方法，常見的有BROYDEN_MIXING和DIRECT_P_MIXING、PULAY_MIXING。

基本參數的解釋到這裡就結束了，至於偷偷潛入其中的&PRINT部分，後面在許多地方比如涉及到需要靜態計算獲取一些物理性質的時候需要用到它們。現在我們暫時忽略。

使用Pymatflow高效完成任務

為了完成上述計算任務，只需使用Pymatflow一條命令即可:

<code>~$ matflow cp2k -r 0 --cutoff 100 --kpoints-scheme "monkhorst-pack 3 3 3" -d nacl-static/<code>

上述命令會在當前路徑創建nacl-static目錄，然後在其下自動生成輸入文件並進行作業任務的提交。默認使用PBS作業調度器。

你也可以通過--server指定其它服務器類型，比如呂梁天河二號超算系統請使用--server llhpc參數。如果你只是想在本地單節點運行，請使用一下參數組合: --runopt genrun、--auto 0。

等待計算任務結束後，你可以通過postflow命令來進行部分計算結構的後處理:

<code>~$ postflow cp2k -r 0 -d nacl-static/<code>

上圖為在nacl-static/post-processing目錄下得到的靜態自洽過程能量的變化。

往期文章

Pymatflow是一個用於第一性原理模擬的工作流自動化管理軟件，目前支持的計算程序有VASP、Qutum ESPRESSO、Abinit、SIESTA、CP2K。

項目文檔地址:

http://pymatflow.readthedocs.org/

中文手冊:

https://pymatflow.readthedocs.io/en/latest/manual_cn.html

歡迎關注以獲取更多關於Pymatflow和第一性原理計算的精彩內容。

閱讀更多 開源科學計算 的文章

關鍵字: 系列 Scheme 靜態