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 

 使用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 

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= 

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

每次移动原子时的步长

INTENSITIES

F

是否计算红外强度。如果设置为T，需要在DFT部分进行偶极矩的计算（关键词MOMENTS）。

NPROC_REP

128

并行计算频率时，每个REPLICA使用的CPU数目

FULLY_PERIODIC

T

避免从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

解决方案有四种：

使用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。