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编译的话，可以参考
下面，介绍使用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中找到
使用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
每次移动原子时的步长
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#!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。