CP2K入門教程-1：CP2K的安裝

CP2K的安裝?

1.1 直接使用二進制版本?

CP2K的安裝有很多種方法。最簡單的方法是直接使用預編譯版本的二進制可執行文件。

?用戶可以選擇從發行版所帶的軟件源安裝預編譯版本的CP2K：

Debian http://packages.debian.org/search?keywords=cp2k

Fedora https://apps.fedoraproject.org/packages/cp2k

Ubuntu http://packages.ubuntu.com/search?keywords=cp2k

?或者直接從CP2K官方網站下載預編譯可執行文件

?sourceforge.net/projects/cp2k/files/precompiled/

?預編譯版本的CP2K運行比較可靠，缺點是由于其采用了沒有優化過的blas庫和lapack庫，計算速度比較慢。

?或者從CP2K源碼（http://www.cp2k.org/download）進行編譯，獲得可執行文件。好處是可以使用MKL等速度更快的數學庫，計算速度比預編譯版本快很多。測試表明，使用Intel編譯器以及MKL數學庫編譯的CP2K執行速度是預編譯版本的3倍左右。

1.2 從源代碼編譯?

?從源碼進行編譯得到的CP2K可執行文件有更快的運行速度。CP2K的編譯比較復雜。希望進行CP2K編譯的話，可以參考

?www.cp2k.org/howto:compile

?下面，介紹使用Intel編譯器以及GNU編譯器編譯CP2K的步驟。下面的編譯均在Debian testing系統上完成。CPU型號為Intel(R) Core(TM) i5-4210M，內核版本為3.17-1。為了減少編譯的工作量，可以從Debian軟件源里安裝各種數學庫，如libint，elpa，fftw3，openblas, libxc等。以root權限運行：

1 apt-get install libxc1 libxc-dev libint-dev libint1 libelpa-dev libelpa0 libopenblas-base libopenblas-dev libfftw3-3 libfftw3-bin libfftw3-dev libfftw3-mpi3 openmpi-bin gcc gfortran g++

?安裝-dev包可以獲得數學庫的靜態鏈接庫，方便進行靜態編譯。

?系統安裝了4.9.1 版本的GNU編譯器，包括gcc以及gfortran，openmpi版本為1.6.5。

?使用Intel Fortran編譯器以及MKL編譯CP2K

?首先，安裝Intel Fortran編譯器以及MKL。我安裝的Intel Fortran編譯環境是composerxe-2011.3.174。其中包含了Fortran編譯器ifort以及MKL，沒有Intel MPI。ifort版本是12.0.3 20110309。Intel編譯器安裝到了/opt/intel目錄中。

?要使用Intel編譯器，打開終端，運行

1 source /opt/intel/bin/compilervars.sh intel64

?然后，編譯安裝MPI運行環境。如果你安裝的Intel編譯環境已經包含了impi，可以跳過該步驟。我安裝的MPI環境是openmpi 1.6.5，目前最新版本是1.8.3 。從http://www.open-mpi.org?上下載源碼，解壓縮后運行：

1 ./configure --prefix=/opt/openmpi-1.6.5 F77=ifort FC=ifort 2 3 make 4 5 make install

?注意，make install時需要root權限。

?從?http://sourceforge.net/projects/cp2k/files

?下載CP2K 2.5.1的源碼，解壓縮，修改arch目錄中的Linux-x86-64-intel.popt文件如下?：

1 CC = cc 2 3 CPP = 4 5 FC = /opt/openmpi-1.6.5/bin/mpif90 6 7 LD = /opt/openmpi-1.6.5/bin/mpif90 8 9 AR = ar -r 10 11 INTEL_MKL = /opt/intel/mkl 12 13 INTEL_INC = $(INTEL_MKL)/include 14 15 INTEL_LIB = $(INTEL_MKL)/lib/intel64 16 17 MKL_LIB = $(INTEL_MKL)/lib/intel64 18 19 FFTW3_INC = /usr/include/ 20 21 FFTW3_LIB = /usr/lib/x86_64-linux-gnu/ 22 23 DFLAGS = -D__INTEL -D__FFTSG -D__parallel -D__BLACS -D__SCALAPACK -D__FFTW3 -D__LIBINT 24 25 CPPFLAGS = -C -traditional $(DFLAGS) -I$(INTEL_INC) -I$(FFTW3_INC) 26 27 FCFLAGS = $(DFLAGS) -I$(INTEL_INC) -I$(FFTW3_INC) -O2 -xHost -heap-arrays 64 -funroll-loops -fpp -free 28 29 FCFLAGS2 = $(DFLAGS) -I$(INTEL_INC) -I$(FFTW3_INC) -O1 -xHost -heap-arrays 64 -fpp -free 30 31 LDFLAGS = $(FCFLAGS) -I$(INTEL_LIB) -L$(FFTW3_LIB) 32 33 LIBS = -L$(MKL_LIB) -lmkl_blas95_lp64 -lmkl_lapack95_lp64 -lmkl_scalapack_lp64 -lmkl_intel_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lmkl_blacs_openmpi_lp64 -lfftw3 -lpthread -lderiv -lint -lstdc++ 34 35 OBJECTS_ARCHITECTURE = machine_intel.o 36 37 graphcon.o: graphcon.F 38 39 $(FC) -c $(FCFLAGS2) $< 40 41 et_coupling.o: et_coupling.F 42 43 $(FC) -c $(FCFLAGS2) $< 44 45 qs_vxc_atom.o: qs_vxc_atom.F 46 47 $(FC) -c $(FCFLAGS2) $< 48 49 hfx_screening_methods.o: hfx_screening_methods.F 50 51 $(FC) -c $(FCFLAGS2) $<

?注意，如果編譯使用intelmpi，blacs庫就設置為-lmkl_blacs_intelmpi_lp64；如果使用openmpi，則設置為-lmkl_blacs_openmpi_lp64。在makefiles目錄中運行

1 make –j 4 ARCH=Linux-x86-64-intel VERSION=popt

?-j 4參數可以使編譯并行進行，以加快編譯速度。

?視CPU性能，一般15分鐘左右即可在exe/Linux-x86-64-intel目錄下獲得編譯好的cp2k.popt可執行文件。

?使用GNU Fortran編譯器,MKL以及ELPA庫編譯CP2K

?使用gfortran編譯器，并使用MKL數學庫和ELPA庫編譯CP2K，流程與上述類似。編輯arch目錄中的Linux-x86-64-gfortran_mkl_elpa.popt文件如下：

1 INTEL_MKL = /opt/intel/mkl 2 3 INTEL_INC = $(INTEL_MKL)/include 4 5 INTEL_MKL_LIB = $(INTEL_MKL)/lib/intel64 6 7 FFTW3_INC = /usr/include/ 8 9 FFTW3_LIB = /usr/lib/x86_64-linux-gnu/ 10 11 ELPA_INC = /usr/include/elpa/modules 12 13 CC = cc 14 15 CPP = 16 17 FC = /usr/bin/mpif90 18 19 LD = /usr/bin/mpif90 20 21 AR = ar -r 22 23 CPPFLAGS = 24 25 DFLAGS = -D__GFORTRAN -D__FFTSG -D__parallel -D__SCALAPACK -D__BLACS -D__LIBINT -D__LIBXC2 -D__FFTW3 -D__ELPA 26 27 FCFLAGS = -O3 -ffast-math -funroll-loops -ftree-vectorize -march=native -ffree-form $(DFLAGS) -g -I$(FFTW3_INC) -I$(ELPA_INC) -I$(INTEL_INC) 28 29 LDFLAGS = $(FCFLAGS) -L${FFTW3_LIB} -L${INTEL_MKL_LIB} 30 31 LIBS = \ 32 33 -lmkl_scalapack_lp64 -lmkl_blacs_openmpi_lp64 -lmkl_gf_lp64 -lmkl_sequential -lmkl_core \ 34 35 -lderiv -lint -lfftw3 -lxc -lelpa

在makefiles目錄中運行

1 make -j 4 ARCH=Linux-x86-64-gfortran_mkl_elpa VERSION=popt

在exe/Linux-x86-64-gfortran_mkl_elpa 目錄中即可獲得編譯好的cp2k.popt可執行文件。

經過測試，使用gfortran和intel編譯器編譯的CP2K可執行文件運行速度相當，后者稍微快一點點。這主要是因為它們均使用了MKL數學庫。

?

CP2K入門教程-2：如何獲取幫助？

CP2K一直在發展中，目前還沒有一個正式的手冊，所以學習CP2K并不是一個簡單的事，連自學的材料都不多。目前我所知道的途徑有以下幾個：

2.1 CP2K的Google Group

?這里是最好的獲得幫助的地方，CP2K的開發者經常在group里回復用戶的問題，非常專業。常見的問題在論壇上都能找到答案。很遺憾，拜國內惡劣的網絡環境所賜，訪問Google Group非常困難。可以使用Xskywalker瀏覽器等專用工具來訪問CP2K的Google Group。

2.2 CP2K 官方教程?

?目前已經舉辦了多次CP2K培訓教程，網站（http://www.cp2k.org/tutorials，

?http://www.cp2k.org/events）上有不少演示文稿可以下載，很多是程序開發人員做的報告。要了解CP2K，讀這些文稿是很好的入門材料。

2.3. CP2K的官方手冊?

CP2K的官方手冊（http://manual.cp2k.org/trunk/）實際上并不是“手冊”，因為這個網站只是解釋了各種關鍵詞的含義以及設置，并沒有教你如何使用CP2K。手冊本身是從CP2K的源碼直接生成的，只要下載了源碼就可以在本地生成CP2K的手冊。

2.4 CP2K源碼包中的測試文件?

?最直接的例子是源碼中的tests文件。CP2K源碼包中的tests目錄包含了各種方法的輸入文件。這些輸入文件并不是最適合計算的，其中測參數設置沒有經過優化。但這些輸入文件給了我們了解輸入文件結構的途徑。

?另外，學會使用grep命令。當你想了解CP2K的某個關鍵詞（keyword）時，不妨使用grep –iR keyword tests/?命令來查看使用了該關鍵詞的測試輸入文件。仔細閱讀這些輸入文件，就能知道這些關鍵詞的使用了。需要注意的是，tests目錄中的輸入文件主要是用來測試程序運行的正常與否，往往使用了不合理的參數，用戶需要參考手冊等其他資料自行進行調整。

2.5 CP2K相關的文獻?

?有關CP2K中使用的各種方法，CP2K的網站上放了一個參考文獻列表http://manual.cp2k.org/trunk/references.html

?要想真正了解CP2K的原理，可以閱讀這些文獻。

?

CP2K入門教程-3：CP2K的功能

? CP2K的功能很多，其輸入文件的結構也比Gaussian或者VASP等常見的計算化學軟件更加復雜。這里，我直接給出各種CP2K輸入文件作為例子，并說明其輸入文件中各種參數的設置。

3.1 CP2K的各種運行方式（RUN_TYPE）

• BAND使用band方法來進行最低能量路徑（MEP）以及反應過渡態的搜索，通常可以使用的方法有CI-NEB, IT-NEB等方法。

• CELL_OPT晶胞參數的優化。

• ENERGY計算當前結構的能量。對于DFT來說，就是一個SCF計算。

• ENERGY_FORCE計算當前結構的能量，同時計算出每個原子上的梯度（即受力）。

• GEO_OPT對當前結構進行幾何結構的優化，以獲得當前結構對應的局部極小點或者過渡態。GEO_OPT有兩種類型，一種是MINIMIZATION，即優化到極小點；一種是TRANSITION_STATE，即優化到過渡態。過渡態的優化采用Dimer方法來實現。

• MD分子動力學模擬。基于第一原理的分子動力學模擬是CP2K的一大優勢。此外，CP2K也支持使用力場以及DFTB等方法進行分子動力學模擬。

• VIBRATIONAL_ANALYSIS對當前結構進行振動分析，也就是頻率計算。頻率計算可以計算部分原子的頻率（INVOLVED_ATOMS），也可以計算某一頻率范圍內的頻率計算（MODE_SELECTIVE）。

?

CP2K入門教程-4：使用CP2K進行DFT計算

很多參數都可以影響CP2K的計算精度。下面我將介紹我所知道的參數。?

3.1.1 晶胞的大小?

CP2K只支持Gamma點的計算，沒有K點。因此，計算中必須使用足夠大的晶胞。如果晶胞太小，部分基組函數就會超過晶胞的邊界，導致重疊矩陣求逆過程出現問題，計算就會不可靠。不能直接使用VASP中的小晶胞來進行CP2K的計算。

3.1.2 CUTOFF以及REL_CUTOFF

CUTOFF和REL_CUTOFF兩個參數用來控制網格的精度。CP2K中的網格從粗糙到精細分為4個級別。CUTOFF參數控制整體網格精度的最高值，REL_CUTOFF參數控制有多少網格點落到最精細的級別。CUTOFF的設置取決于體系中元素的種類，默認值為280 Ry，但對有些原子需要設置到500 Ry甚至更高。

CP2K論壇中建議給不同的原子使用不同的CUTOFF。下圖給出了對不同原子建議使用的CUTOFF大小。可見，對包含Na，N，O，F，Ne，Ni，Ga等元素的計算，需要設置高達1000 Ry的CUTOFF來確保計算精度。在計算中包含這些元素時，需要額外小心。

REL_CUTOFF默認值為50 Ry，一般設置到60 Ry時精度就已經足夠了。

?除此以外，還有USE_FINER_GRID參數等，用于提高網格的精細程度。

3.1.3 基組和泛函?

CP2K中可以使用的泛函很多，但并非每個基組都為相應的泛函進行了優化。經常使用的泛函有LDA（PADE），BLYP以及PBE，在CP2K的tests目錄中有相應的優化基組。此外，CP2K中也可以使用B3LYP、HSE等雜化泛函，可以使用DFT-D3等色散校正。在CP2K中使用B3LYP泛函時，關鍵輸入文件如下：

1 2 3 &DFT 4 5 BASIS_SET_FILE_NAME ./BASIS_MOLOPT 6 7 POTENTIAL_FILE_NAME ./POTENTIAL 8 9 CHARGE 0 10 11 MULTIPLICITY 1 12 13 &SCF 14 15 SCF_GUESS ATOMIC 16 17 EPS_SCF 1.0E-6 18 19 MAX_SCF 50 20 21 &OUTER_SCF 22 23 MAX_SCF 10 24 25 &END OUTER_SCF 26 27 &OT 28 29 # My scheme 30 31 PRECONDITIONER FULL_SINGLE_INVERSE 32 33 MINIMIZER DIIS 34 35 N_DIIS 7 36 37 &END OT 38 39 &PRINT 40 41 &RESTART 42 43 &EACH 44 45 MD 20 46 47 &END EACH 48 49 &END RESTART 50 51 &RESTART_HISTORY OFF 52 53 &END RESTART_HISTORY 54 55 &END PRINT 56 57 &END SCF 58 59 60 61 &QS 62 63 METHOD GAPW 64 65 # My scheme 66 67 EPS_DEFAULT 1.0E-12 68 69 EPS_PGF_ORB 1.0E-32 70 71 EPS_FILTER_MATRIX 0.0E+0 72 73 &END QS 74 75 &MGRID 76 77 COMMENSURATE 78 79 CUTOFF 300 80 81 &END MGRID 82 83 &POISSON 84 85 POISSON_SOLVER MULTIPOLE 86 87 PERIODIC NONE 88 89 &MULTIPOLE 90 91 RCUT 40 92 93 &END MULTIPOLE 94 95 &END POISSON 96 97 &XC 98 99 #&XC_FUNCTIONAL BLYP 100 101 #&END XC_FUNCTIONAL 102 103 &XC_FUNCTIONAL 104 105 &LYP 106 107 SCALE_C 0.81 108 109 &END 110 111 &BECKE88 112 113 SCALE_X 0.72 114 115 &END 116 117 &VWN 118 119 FUNCTIONAL_TYPE VWN3 120 121 SCALE_C 0.19 122 123 &END 124 125 &XALPHA 126 127 SCALE_X 0.08 128 129 &END 130 131 &END XC_FUNCTIONAL 132 133 &HF 134 135 &SCREENING 136 137 EPS_SCHWARZ 1.0E-10 138 139 &END 140 141 &MEMORY 142 143 MAX_MEMORY 512 144 145 EPS_STORAGE_SCALING 1.0E-1 146 147 &END 148 149 FRACTION 0.20 150 151 &END 152 153 &XC_GRID 154 155 XC_SMOOTH_RHO NN10 156 157 XC_DERIV SPLINE2_SMOOTH 158 159 &END XC_GRID 160 161 &END XC 162 163 &END DFT 164 165 166 167

?完整的輸入文件可以在Google Group中找到

?https://groups.google.com/forum/?fromgroups#!searchin/cp2k/XC_SMOOTH_RHO|sort:date/cp2k/v32E3xXmgaY/-TygNY7t2msJ?

?使用HSE泛函的輸入文件例子：

1 &XC_FUNCTIONAL 2 3 &XWPBE 4 5 SCALE_X -0.25 6 7 SCALE_X0 1.0 8 9 OMEGA 0.11 10 11 &END 12 13 &PBE 14 15 SCALE_X 0.0 16 17 SCALE_C 1.0 18 19 &END PBE 20 21 &END XC_FUNCTIONAL 22 23 &HF 24 25 EPS_SCHWARZ 1.0E-10 26 27 MAX_MEMORY 10 28 29 FRACTION 0.25 30 31 SCREENING_TYPE SHORTRANGE 32 33 OMEGA 0.11 34 35 &END 36 37 38 39

?基組的大小也有很多種，如SZ，DZVP以及TZVP。一般計算中使用DZVP基組就足夠了。

3.1.4 SCF收斂精度?

?對于一般的測試性結算，SCF收斂到1E-5就可以得到相對準確的結果。如果要進行比較高精度的計算，可以將SCF收斂精度設置為1E-6。對于頻率計算等精度要求更高的計算，SCF收斂精度可以設置為1E-7。

3.1.5 收斂算法的選擇?

CP2K中主要有兩種SCF的收斂算法，一種是基于軌道變換（OT）的算法，一種是基于對角化(DIAG)的算法。如果體系有較大帶隙的，如為半導體或者絕緣體等，推薦使用OT算法，收斂速度比較快。如果體系中HOMO-LUMO 帶隙很小或者幾乎沒有，如金屬體系，則建議使用對角化的方法進行計算,并使用smear方法。

?下面是一個對Rh(1 1 1)表面進行計算使用的輸入文件的例子:

1 2 3 &SCF 4 5 SCF_GUESS RESTART 6 7 EPS_SCF 5.0E-7 8 9 MAX_SCF 500 10 11 ADDED_MOS 500 12 13 CHOLESKY INVERSE 14 15 &SMEAR ON 16 17 METHOD FERMI_DIRAC 18 19 ELECTRONIC_TEMPERATURE [K] 300 20 21 &END SMEAR 22 23 &DIAGONALIZATION 24 25 ALGORITHM STANDARD 26 27 &END DIAGONALIZATION 28 29 &MIXING 30 31 METHOD BROYDEN_MIXING 32 33 ALPHA 0.1 34 35 BETA 1.5 36 37 NBROYDEN 8 38 39 &END MIXING 40 41 &PRINT 42 43 &RESTART 44 45 &EACH 46 47 QS_SCF 50 48 49 &END 50 51 ADD_LAST NUMERIC 52 53 &END RESTART 54 55 &END PRINT 56 57 &END SCF 58 59

?注意使用對角化方法必須使用ADDED_MOS?關鍵詞。另外，設置正確的MIXING方案也是加速收斂的關鍵。

?實際上，即使是對于非金屬體系，有時候對角化方法也會比OT算法速度更快。?

?所以，在進行大規模的計算之前最好進行充分的測試。

?使用OT算法時，優化算法也有多重選擇。常用的有CG，DIIS以及BROYDEN。其中，CG算法是最為穩定的算法，一般的計算都可以使用CG算法。DIIS算法速度比較快，但不夠穩定。如果CG算法和DIIS算法收斂都有問題時，可以嘗試使用BROYDEN算法。下面是使用BROYDEN算法的輸入文件例子。

1 &OT T 2 3 MINIMIZER BROYDEN 4 5 N_HISTORY_VEC 4 6 7 BROYDEN_BETA 6.9999999999999996E-01 8 9 BROYDEN_SIGMA 1.4999999999999999E-01 10 11 LINESEARCH 2PNT 12 13 PRECONDITIONER FULL_SINGLE_INVERSE 14 15 &END OT 16 17 18 19

3.1.6 EPS_DEFAULT的設置?

QS部分使用的EPS_DEFAULT為1.0E-10。根據Google Group中的建議，在進行比較高精度計算時，需要將EPS_DEFAULT設為1.0E-14。

CP2K中有多個參數依賴EPS_DEFAULT，包括EPS_CORE_CHARGE，EPS_GVG_RSPACE，EPS_PGF_ORB，EPS_KG_ORB。各個參數的精度以及含義如下：

名稱

含義

默認精度

1 EPS_CORE_CHARGE

核電荷映射精度

1 EPS_DEFAULT/100.0 2 3 EPS_GVG_RSPACE

實空間KS矩陣元積分精度

1 SQRT(EPS_DEFAULT) 2 3 EPS_PGF_ORB

重疊矩陣元精度

1 SQRT(EPS_DEFAULT) 2 3 EPS_KG_ORB

使用Kim-Gordon方法時的精度

1 SQRT(EPS_DEFAULT)

?更多信息，可以參考：

?http://developer.berlios.de/forum/message.php?msg_id=35587?

?https://groups.google.com/forum/?fromgroups#!topic/cp2k/IIp-D6AA7ME?

3.2 使用CP2K進行能量計算 3.2.1 OUTER_SCF

?計算能量，需要設置RUN_TYPE為ENERGY；如果還要進行梯度（受力）計算，則設置為ENERGY_FORCE。

?使用DFT方法進行能量計算，使用OT算法，并開啟OUTER_SCF，示例如下：

1 &SCF 2 3 EPS_SCF 1.0E-6 4 5 SCF_GUESS RESTART 6 7 MAX_SCF 100 8 9 &OT T 10 11 PRECONDITIONER FULL_ALL 12 13 MINIMIZER DIIS 14 15 LINESEARCH 3PNT 16 17 &END OT 18 19 &OUTER_SCF ON 20 21 MAX_SCF 5 22 23 EPS_SCF 5.0E-6 24 25 &END OUTER_SCF 26 27 &END SC

?其中，OUTER_SCF為加速收斂的一種方法。以上面的輸入文件為例，計算過程中，如果SCF經過100次優化依然沒有收斂，則進入OUTER_SCF過程,對前一次計算的波函數進行調整，重新進行SCF迭代。每次OUTER_SCF中優化的次數依然是100次，最多可以進行5次OUTER_SCF。所以，最多可以進行500次SCF計算。

?更多的細節請參考：

?https://groups.google.com/forum/?fromgroups=#!topic/cp2k/6cikFLKeN34?

3.2.2 輸出每個原子上的受力

?要輸出每個原子上的受力，GLOBAL部分的RUN_TYPE必須設置為ENERGY_FORCE或者GEO_OPT。要輸出受力，需要在FORCE_EVAL部分開啟選項：

1 &PRINT 2 3 &FORCES ON 4 5 Filename ForceFileName 6 7 &END FORCES 8 9 &END PRINT

?如果要將受力信息存儲到文件中，則需要設定文件的名稱；否則受力信息將會打印到out文件中。輸出的受力格式如下：

1 ATOMIC FORCES in [a.u.] 2 3 4 5 # Atom Kind Element X Y Z 6 7 1 1 O 0.08722700 -0.04704030 0.08194080 8 9 2 2 H -0.07829459 0.00721899 -0.00996929 10 11 3 2 H -0.01049003 0.03981616 -0.06774948 12 13 SUM OF ATOMIC FORCES -0.00155761 -0.00000516 0.00422203 0.00450019

CP2K入門教程-5：使用CP2K進行幾何優化計算

要使用CP2K進行幾何優化計算，需要設置GLOBAL 部分的RUN_TYPE為GEO_OPT

1 &GLOBAL 2 3 PROJECT CP2K 4 5 RUN_TYPE GEO_OPT 6 7 &END GLOBAL

幾何優化的種類有兩種，一種是正常的能量極小化優化，一種是使用dimer算法進行過渡態優化。默認為能量極小化計算。下面將分別介紹這兩種計算的輸入文件。

3.3.1 能量極小化計算

使用CP2K進行常規的幾何優化，關鍵輸入文件如下：

1 &MOTION 2 3 &GEO_OPT 4 5 TYPE MINIMIZATION 6 7 MAX_ITER 400 8 9 OPTIMIZER LBFGS 10 11 MAX_FORCE 4.0E-4 12 13 &END GEO_OPT 14 15 &END MOTION

需要注意的是，幾何優化有三種算法，分別是CG、BFGS和LBFGS。其中，CG算法是最穩定的算法，但計算速度相對較慢；BFGS算法效率最高，計算中需要對Hessian矩陣進行對角化，如果初始結構不合理，BFGS算法容易出問題；LBFGS算法效率和BFGS類似，同時穩定性也很好。對于一般的幾何優化，推薦使用LBFGS算法。

如果在幾何優化過程中需要固定部分原子，可以在MOTION部分啟用CONSTRAINT選項。例子如下：

1 &CONSTRAINT 2 3 &FIXED_ATOMS 4 5 LIST 1 2 3 4 6 7 LIST 12 .. 43 8 9 LIST 76..91 10 11 &END FIXED_ATOMS 12 13 &END CONSTRAINT

3.3.2 使用Dimer算法進行過渡態優化

使用Dimer算法進行過渡態優化的關鍵輸入文件如下：

1 &MOTION 2 3 &GEO_OPT 4 5 TYPE TRANSITION_STATE 6 7 MAX_ITER 400 8 9 OPTIMIZER CG 10 11 MAX_FORCE 4.5E-4 12 13 &CG 14 15 &LINE_SEARCH 16 17 TYPE 2PNT 18 19 &END LINE_SEARCH 20 21 &END CG 22 23 &TRANSITION_STATE 24 25 METHOD DIMER 26 27 &DIMER 28 29 DR 0.01 30 31 ANGLE_TOLERANCE [deg] 4.0 32 33 INTERPOLATE_GRADIENT 34 35 &ROT_OPT 36 37 OPTIMIZER CG 38 39 MAX_ITER 10 40 41 #Loose these value. Just for faster converge. 42 43 MAX_DR 3.0E-3 44 45 MAX_FORCE 4.5E-4 46 47 &CG 48 49 # MAX_STEEP_STEPS 15 50 51 &LINE_SEARCH 52 53 TYPE 2PNT 54 55 &END LINE_SEARCH 56 57 &END CG 58 59 &END ROT_OPT 60 61 &END DIMER 62 63 &END TRANSITION_STATE 64 65 &END GEO_OPT 66 67 &END MOTION

注意，使用dimer方法進行過渡態搜索計算時，只能使用CG優化算法，不能用BFGS或者LBFGS。

使用dimer方法進行計算的時候，可以手工給定初始的Dimer Vector。Dimer Vector應該是指向過渡態方向的一個多維向量，使dimer方法在進行過渡態搜尋的時候沿著該方向進行搜索，可以提高搜索的效率。如果沒有給定初始Dimer Vector，CP2K程序中會隨機設定一個初始的搜索方向。

要理解dimer算法中的各種參數的含義，請參考

Henkelman, G; Jonsson, H. JOURNAL OF CHEMICAL PHYSICS, 111 (15), 7010-7022 (1999). http://dx.doi.org/10.1063/1.480097

A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives.

3.4 非周期性體系的計算

要使用CP2K對團簇等非周期性體系進行計算，需要額外進行一些設置：

首先，%FORCE_EVAL%DFT中POISSON部分要進行設置

1 &POISSON 2 3 PERIODIC none 4 5 POISSON_SOLVER wavelet 6 7 &END POISSON

如果使用MT作為POISSON_SOLVER，則設置為：

1 &POISSON 2 3 PERIODIC NONE 4 5 POISSON_SOLVER MT 6 7 &MT 8 9 ALPHA 7.0 10 11 REL_CUTOFF 1.2 12 13 &END MT 14 15 &END POISSON

對于周期性計算，POISSON_SOLVER設置為PERIODIC；對于非周期性計算，可以設置為MT或者 WAVELET，兩者略有不同。如果設置為MT，要保證計算使用的單胞體積足夠大，至少是電荷密度的兩倍。如果設置為WAVELET，不需要設置非常大的單胞，但分子必須處于單胞的中心，確保單胞的邊界處電子密度為0。可以使用TOPOLOGY參數來強制將分子置于單胞的中心。

1 &TOPOLOGY 2 3 &CENTER_COORDINATES 4 5 &END CENTER_COORDINATES 6 7 &END TOPOLOGY

此外，CELL部分周期性也要設置為NONE。

1 &CELL 2 3 ABC 20.0 20.0 20.0 4 5 PERIODIC NONE 6 7 &END CELL

3.5 使用CP2K進行晶胞參數的優化

CP2K可以對晶體的晶胞參數進行直接優化。

首先，GLOBAL部分需要設置RUN_TYPE為CELL_OPT

1 &GLOBAL 2 3 PROJECT cp2k 4 5 RUN_TYPE CELL_OPT 6 7 PRINT_LEVEL MEDIUM 8 9 &END GLOBAL

然后，在MOTION部分需要設置CELL_OPT的參數

1 &CELL_OPT 2 3 TYPE GEO_OPT 4 5 OPTIMIZER CG 6 7 MAX_ITER 200 8 9 EXTERNAL_PRESSURE 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 10 11 PRESSURE_TOLERANCE 0.1 12 13 KEEP_ANGLES 14 15 KEEP_SYMMETRY 16 17 &CG 18 19 &LINE_SEARCH 20 21 TYPE 2PNT 22 23 &2PNT 24 25 &END 2PNT 26 27 &END LINE_SEARCH 28 29 &END CG 30 31 &END CELL_OPT 32 33 &GEO_OPT 34 35 MAX_ITER 400 36 37 OPTIMIZER LBFGS 38 39 MAX_FORCE 4.0E-4 40 41 &END GEO_OPT

?

此外，在FORCE_EVAL部分需要設置STRESS_TENSOR的計算方法為ANALYTICAL。如果晶胞本身有特定的對稱性，可以在CELL部分增加晶胞的對稱性限制（SYMMETRY）。

更多細節可以參考CP2K源碼包中的tests/QS/regtest-gpw-4/cell-1.inp文件。

CP2K入門教程-6：多重度計算

當體系中有多個不成對電子的時候，計算的時候就需要考慮多重度問題。此時，計算必須是自旋非限制的，即使用LSD或者UKS參數。多重度的計算有兩種方案。

第一種方案是手工指定多重度，使用MUTIPLICITY關鍵詞設定多重度。這種方法最為簡單可靠，但如果體系自旋多重度的可能性很多，就需要進行多次計算。

第二種方案是自動進行多重度的猜測，使用RELAX_MULTIPLICITY關鍵詞。類似的方法在Dmol3程序中也有實現。下面是使用這種方法進行計算使用的輸入文件:

1 &DFT 2 3 LSD 4 5 BASIS_SET_FILE_NAME ./BASIS_MOLOPT 6 7 POTENTIAL_FILE_NAME ./rr_pot 8 9 WFN_RESTART_FILE_NAME ./cp2k-RESTART.wfn 10 11 CHARGE 0 12 13 RELAX_MULTIP 0.001 14 15 &MGRID 16 17 CUTOFF 300 18 19 &END MGRID 20 21 &QS 22 23 EPS_DEFAULT 1.0E-14 24 25 WF_INTERPOLATION ASPC 26 27 EXTRAPOLATION_ORDER 3 28 29 &END QS 30 31 &SCF 32 33 ADDED_MOS 50 50 34 35 EPS_SCF 5.0E-7 36 37 SCF_GUESS RESTART 38 39 MAX_SCF 200 40 41 CHOLESKY INVERSE 42 43 &DIAGONALIZATION 44 45 ALGORITHM STANDARD 46 47 &END DIAGONALIZATION 48 49 &MIXING 50 51 METHOD BROYDEN_MIXING 52 53 ALPHA 0.1 54 55 BETA 1.5 56 57 NBROYDEN 8 58 59 &END MIXING 60 61 &OUTER_SCF ON 62 63 MAX_SCF 5 64 65 EPS_SCF 1.0E-6 66 67 &END OUTER_SCF 68 69 &PRINT 70 71 &RESTART 72 73 &EACH 74 75 QS_SCF 100 76 77 &END EACH 78 79 ADD_LAST NUMERIC 80 81 &END RESTART 82 83 &END PRINT 84 85 &END SCF 86 87 &XC 88 89 &XC_FUNCTIONAL PBE 90 91 &END XC_FUNCTIONAL 92 93 &END XC 94 95 &END DFT

?

使用這種方法，需要注意以下幾點：

必須使用自旋非限制計算，即開啟UKS或者LSD。

必須使用對角化方法，不能使用OT算法。由于使用對角化方法，也必須使用ADDED_MOS參數。

不能使用SMEAR方法。

RELAX_MULTIP 設置為大于0的值，就開啟自旋優化模式。設置的值越大，自旋翻轉發生的概率越大。

綜上，盡管這種方法看似很誘人，但在使用中很受限制。

CP2K入門教程-7：使用NEB方法進行過渡態搜索

???? CP2K中可以使用NEB方法進行過渡態的搜索。

要使用NEB方法，首先在GLOBAL部分設置RUN_TYPE為BAND。然后，需要在MOTION部分設置NEB計算的參數。輸入文件例子如下：

1 &MOTION 2 3 &BAND 4 5 NPROC_REP 32 6 7 BAND_TYPE IT-NEB 8 9 NUMBER_OF_REPLICA 6 10 11 K_SPRING 0.02 12 13 &CONVERGENCE_CONTROL 14 15 MAX_DR 0.01 16 17 MAX_FORCE 0.001 18 19 RMS_DR 0.02 20 21 RMS_FORCE 0.0005 22 23 &END 24 25 ROTATE_FRAMES F 26 27 &CI_NEB 28 29 NSTEPS_IT 5 30 31 &END 32 33 &OPTIMIZE_BAND 34 35 OPT_TYPE MD 36 37 # OPTIMIZE_END_POINTS F 38 39 OPTIMIZE_END_POINTS T 40 41 &MD 42 43 TIMESTEP 0.5 44 45 TEMPERATURE 500.0 46 47 MAX_STEPS 300 48 49 &VEL_CONTROL 50 51 ANNEALING 0.99 52 53 PROJ_VELOCITY_VERLET T 54 55 &END 56 57 &END 58 59 &END 60 61 &REPLICA 62 63 COORD_FILE_NAME ./1.xyz 64 65 &END 66 67 &REPLICA 68 69 COORD_FILE_NAME ./2.xyz 70 71 &END 72 73 &REPLICA 74 75 COORD_FILE_NAME ./3.xyz 76 77 &END 78 79 &REPLICA 80 81 COORD_FILE_NAME ./4.xyz 82 83 &END 84 85 &REPLICA 86 87 COORD_FILE_NAME ./5.xyz 88 89 &END 90 91 &REPLICA 92 93 COORD_FILE_NAME ./6.xyz 94 95 &END 96 97 &PROGRAM_RUN_INFO 98 99 &END 100 101 &CONVERGENCE_INFO 102 103 &END 104 105 &END BAND 106 107 &CONSTRAINT 108 109 &FIXED_ATOMS 110 111 LIST 1.. 4 112 113 LIST 12 .. 43 114 115 LIST 76..91 116 117 &END FIXED_ATOMS 118 119 &END CONSTRAINT 120 121 &END MOTION

?

對于以上輸入文件中的參數，解釋如下：

關鍵詞

示例中的設置

解釋

NPROC_REP

32

進行BAND計算時，每個REPLICA使用的CPU數目

1 BAND_TYPE 2 3 IT-NEB

BAND計算方法類型。有IT-NEB，CI-NEB，B-NEB，D-NEB，EB，SM等多種。推薦使用IT-NEB以及CI-NEB

1 NUMBER_OF_REPLICA 2 3 6

BAND計算中使用的REPLICA的總數目。REPLICA的數目越多，計算越準確。如果使用較少的REPLICA無法得到正確的結果，可以考慮增加REPLIA的數目。CP2K使用的CPU總數目為NUMBER_OF_REPLICA*NPROC_REP。在本例中，就是32*6=192

1 K_SPRING 2 3 0.02

BAND計算中使用的彈簧勁度系數。K_SPRING越大，NEB計算收斂越快，但計算不準確。K_SPING越小，計算收斂越慢，但計算較為準確。在初步計算中，可以將K_SPRING設置為0.08左右，然后再放松至0.02以獲得精確結果。

?

CP2K入門教程-8：振動頻率分析

????? 用CP2K程序進行振動頻率分析，首先需要設置RUN_TYPE為VIBRATIONAL_ANALYSIS。輸入文件例子如下： 1 &GLOBAL 2 3 PROJECT cp2k 4 5 RUN_TYPE VIBRATIONAL_ANALYSIS 6 7 PRINT_LEVEL medium 8 9 &END GLOBAL

然后，設置頻率分析部分輸入文件

1 &VIBRATIONAL_ANALYSIS 2 3 DX 0.01 4 5 INTENSITIES F 6 7 NPROC_REP 128 8 9 FULLY_PERIODIC T 10 11 &END VIBRATIONAL_ANALYSIS

CP2K計算頻率使用的是數值算法，即對每個原子向+x, -x, +y, -y, +z, -z 6個方向分別進行移動，用數值的方法得到能量的二階導（即力常數），然后計算頻率。所以，如果有N個原子要進行移動，總共要進行6N+1次SCF收斂計算。

關鍵詞

設置示例

解釋

1 DX 2 3 0.01

每次移動原子時的步長

INTENSITIESF

是否計算紅外強度。如果設置為T，需要在DFT部分進行偶極矩的計算（關鍵詞MOMENTS）。

NPROC_REP128

并行計算頻率時，每個REPLICA使用的CPU數目

FULLY_PERIODICT

避免從Hessian矩陣中消除轉動模式。開啟該關鍵詞后，對于N個原子的體系會計算出3N-3個頻率，其中包含了3個轉動自由度。

要計算部分原子的振動頻率，有兩種辦法。一種是直接在MOTION中使用CONSTRAINT對不需要進行頻率分析的原子進行固定。一種是使用MODE_SELECTIVE模式。例子如下：

1 &VIBRATIONAL_ANALYSIS 2 3 NPROC_REP 16 4 5 DX 0.01 6 7 INTENSITIES T 8 9 &MODE_SELECTIVE 10 11 ATOMS 82 83 12 13 INITIAL_GUESS ATOMIC 14 15 EPS_NORM 1.0E-5 16 17 EPS_MAX_VAL 1.0E-6 18 19 &INVOLVED_ATOMS 20 21 INVOLVED_ATOMS 82 83 22 23 &END INVOLVED_ATOMS 24 25 &END &MODE_SELECTIVE 26 27 &END VIBRATIONAL_ANALYSIS

上面的例子中，對82 和83兩個原子進行了振動頻率分子。需要注意的是，使用這種方法計算頻率，使用的REPLICA數目不能太少。REPLICA的數目是這樣計算的：NREP=總CPU數目/NPROC_REP。上述輸入文件，如果使用的總CPU數目為64，則共有NREP=4，即共有4個REPLICA。如果只使用一個REPLICA，使用MODE_SELECTIVE算法計算頻率時，就會只跟蹤一個頻率，無法得到正確的結果。

另外，使用CP2K程序計算一個優化好的結構式的頻率時，也常會出現多個虛頻。這并非是幾何優化出現了問題，而是CP2K計算使用GTH贗勢時存在的一個問題。詳細內容請參考：

https://groups.google.com/forum/?fromgroups#!topic/cp2k/DVCV0epl7Wo

解決方案有四種：

使用NLCC贗勢。http://arxiv.org/abs/1212.6011 不過，NLCC贗勢很不完整，只有B-Cl的元素有，且只提供了PBE泛函的贗勢。

增大CUTOFF，使用600 Ry以上的CUTOFF。

在XC_GRID部分使用平滑參數SMOOTING。不推薦使用。

在XC_GRID部分使用USE_FINER_GRID。加上這個參數后，XC部分的格點的精度提高為4*CUTOFF。

轉載于:https://www.cnblogs.com/Shine-JK/p/10988556.html

總結

以上是生活随笔為你收集整理的CP2K入门教程分享的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。