CP2K入门

type

status

date

slug

summary

CP2K 的输入文件一般命名为 cp2k.inp，其中关键词以 section 和 subsection 的形式，一环套一环，如下图所示。

每一个 section 都是以 &section 开头，以 &END section 结尾，顺序随意，但是嵌套不能乱。

在 section 中，每行只写一个关键词，关键词后面接参数，CP2K 对大小写和空格不敏感。

GLOBAL

GLOBAL 控制 CP2K 的全局参数，GLOBAL 的设置如下所示：

PROJECT 指定计算任务名，个人习惯是统一设置成 cp2k，方便脚本后处理。

RUN_TYPE 指定 CP2K 计算任务的类型，包括

MD，分子动力学模拟；

ENERGY，单点能计算；

ENERGY_FORCE，计算能量和原子受力；

GEO_OPT，几何优化；

CELL_OPT，晶胞优化，常与 GEO_OPT 联用；

BAND，NEB计算；

PRINT_LEVEL 控制输出详略，一般设置成 LOW 或 MEDIUM 即可。

CP2K 中计算能量需要设置 GLOBAL/RUN_TYPE 为 ENERGY；如果还要计算受力，需要设置 GLOBAL/RUN_TYPE 为 ENERGY_FORCE，并设置 FORCE_EVAL/PRINT/FORCES 为 ON。如果要将受力信息存储到文件中，则需要设定 FROCE_EVAL/PRINT/FORCES/FILENAME 文件名；否则受力信息将会打印到 *.out 文件中。

注意，在 CP2K 中，默认的能量和距离单位分别是 Hartree (Hartree 原子单位) 和 Bohr (波尔半径)。如果需要将它们转换成常用的 eV（电子伏特）和 Å（埃）单位，可以使用以下换算关系：

反之：

CP2K 如果每步SCF警告 *** WARNING in fm/cp_fm_elpa.F:522 :: Setting real_kernel for ELPA failed *** 的问题，说明 CP2K 没有编译 ELPA 数学库，在 &GLOBAL 中设置 PREFERRED_DIAG_LIBRARY SL 即可。

FORCE_EVAL

FORCE_EVAL 控制能量和原子受力（类似于 VASP 中的电子步），是 CP2K 中最重要的 section。以 SiSb 的计算为例，一个典型的 FORCE_EVAL（使用 DIAG 算法）如下所示：

CP2K 通常默认使用单 Gamma 点计算，需要使用足够大的晶胞以减少周期性误差。如果晶胞太小，部分基组函数可能超出晶胞边界，导致重叠矩阵求逆过程出现数值问题，从而使得结果不可靠。可以考虑增大晶胞尺寸或在输入文件中调整边界条件（如 POISSON 模块中的 PERIODIC NONE），并适当提高平面波基组的截断能量（CUTOFF）和相对截断能量（REL_CUTOFF）。

切记，不能直接使用 VASP 中的小晶胞来进行 CP2K 的计算。需要使用多 k 点时，可以在 KPOINTS 中指定网格的 k 点采样，如下所示：

CP2K 计算中需要指定赝势和基组文件（类似于 VASP 中的 POTCAR）。赝势文件通常选择位于 GTH_POTENTIALS 目录下的 PBE 或 BLYP 泛函文件，并根据元素的价电子数选择不同的 -q 后缀，如 -q4 表示该赝势包含 4 个价电子。基组文件常用 BASIS_MOLOPT，它经过专门的优化，适用于大部分元素体系；而 BASIS_MOLOPT_UCL 是伦敦大学学院提供的扩展基组，适用于更多类型的元素。

在基组文件中，SZV（Single Zeta Valence）、DZVP（Double Zeta Valence with Polarization）、TZVP（Triple Zeta Valence with Polarization）分别代表基组的劈裂（多少层电子壳层）和极化函数的添加情况。其中，VP（Valence Polarized）表示基组中包含极化函数。通常，劈裂层数和极化函数越多，计算结果越精确，但计算量也会显著增加。计算量的顺序一般为：SZV < DZVP < TZVP < TZV2P < TZV2PX。

选择基组时应注意与赝势中的价电子数保持一致。基组的大小对基组重叠误差（Basis Set Superposition Error, BSSE）有显著影响。较大的基组（如 TZVP 或 TZV2P）通常可以有效降低 BSSE 误差，因此对于大体系或精度要求较高的计算，推荐使用 DZVP 或 TZVP。在计算过程中，可以使用 COUNTERPOISE 方法来校正 BSSE 误差，确保结果的精度。

FORCE_EVAL 中的 SCF 收敛算法

CP2K 提供了两种常用的 SCF（Self-Consistent Field）收敛算法：

对角化算法（DIAG）：基于直接对角化的方式求解 KS（Kohn-Sham）方程。对于体系的带隙很小（如半导体或金属）或几乎没有（如金属体系），推荐使用 DIAG 算法，并开启电子展宽（SMEAR）来帮助收敛。

注意：在使用 DIAG 算法时，必须设置 ADDED_MOS 参数来指定额外的虚拟轨道数量，以防止出现电子填充不完全的问题，并提升收敛效果。 ADDED_MOS 的值通常需要设置为体系中价带电子数的 5-10% 左右，但对于金属体系或需要处理激发态的计算，可能需要更多的虚拟轨道以确保所有电子可以正确填充。

轨道变换算法（OT）：一种基于轨道变换的 SCF 优化方法，通常比对角化算法更高效，尤其适用于具有较大带隙的绝缘体或介电材料体系。在体系中若带隙较大或者体系较为复杂时，OT 算法往往能显著加快 SCF 的收敛。

在使用 DIAG 或 OT 算法时，设置合适的混合参数（MIXING）是加速收敛的关键。通常，对于金属体系（小带隙或无带隙），较低的混合参数（如 MIXING 0.1）可以避免 SCF 振荡；而对于带隙较大的绝缘体，较高的混合参数（如 MIXING 0.4）则有助于加快收敛。

DIAG 算法的输入文件参考上一节。 OT 算法的 .inp 输入文件如下所示（只展示 FORCE_EVAL/DFT 下的 SCF），替换上一节中 SCF 即可：

OT/MINIMIZER 常用的设置有 CG、DIIS 以及 BROYDEN。其中，CG 是最为稳定的算法，一般的计算任务都可以使用 CG 算法，并设置 OT/LINESEARCH 为比默认 2PNT 更贵但也更稳健的 3PNT。 DIIS 算法速度比较快，但没有 CG 稳定。如果 CG 算法和 DIIS 算法收敛都有问题时，可以尝试使用 BROYDEN 算法。采用 OT 时，推荐开启 SCF/OUTER_SCF，这是加速收敛的一种方法。开启 OUTER_SCF 时 SCF/MAX_SCF 应当非常小，15 至 35 是比较合适的范围。 OUTER_SCF 迭代圈数通过 SCF/OUTER_SCF/MAX_SCF 控制，一般设为 5 即可。实际 SCF 迭代次数上限等于 SCF/MAX_SCF *（1+SCF/OUTER_SCF/MAX_SCF）。另外，SCF/EPS_SCF 设置的是 SCF 总的收敛标准（一般是 1E-6 ），而 SCF/OUTER_SCF/EPS_SCF 设置的是 OUTER_SCF 的收敛标准（一般是 1E-5），后者绝对值应当大于或等于前者。 OT/PRECONDITIONER 中，FULL_ALL 通常最稳定，但耗时最长（GAPW 计算只能使用 FULL_ALL 另说）；大体系可以尝试 FULL_SINGLE_INVERSE 和 FULL_KINETIC。 OT/ALGORITHM 可以从默认的 STRICT 改为 IRAC，SCF 收敛会更稳健。但实际上，即使是对于非金属体系，有时候 DIAG 算法也可能比 OT 算法速度更快。所以，在进行大规模的计算之前最好进行充分的测试。

SCF 收敛相关

当 CP2K 计算出现 SCF（Self-Consistent Field）难以收敛的情况时，首先应检查体系结构的合理性，以确保没有出现原子位置不佳或结构畸变的问题。以下是一些常用的收敛优化方法：

增大 SCF 迭代次数：如果 SCF 过程表现出一定的收敛趋势，但无法在当前迭代次数内收敛，可以通过增加 SCF/MAX_SCF 来允许更多的迭代次数继续计算，直至满足收敛标准。

初猜波函数优化：通常较小的基组（例如 DZVP）比大的基组（例如 TZVP）更容易收敛。可以先用较小基组进行一次计算，并将其收敛的波函数保存下来，作为大基组计算的初始猜测波函数（initial guess），这将有助于提高 SCF 的收敛性。

检查截断能设置： DFT 计算中，DFT/MGRID 中的 CUTOFF 和 REL_CUTOFF 参数用于控制格点的精度。若截断能量设置不够大，不仅会导致计算结果不准确，还会增加 SCF 的收敛难度。通常建议将 CUTOFF 设置为 400-600 Ry，并将 REL_CUTOFF 设置为 40-60 Ry，根据实际体系进行适当调整。

调整 SCF 混合策略（MIXING）：设置 SCF/MIXING/METHOD 为 BROYDEN_MIXING 来控制 SCF 过程中旧密度矩阵与新密度矩阵的混合比例，有利于加快 SCF 的收敛速度。

ALPHA：表示新密度矩阵混入旧密度矩阵的比例，默认为 0.4。当 SCF 难以收敛时，可尝试减小 ALPHA 值（如 0.1、0.2、0.3）以增加收敛稳定性。

NBROYDEN：控制 Broyden 方法中使用的历史矩阵数量。默认值为 4，通常过小，可将其增大到 8 或 12 以提升收敛效果。

电子温度控制：对于金属体系或半导体体系，可以在 SCF/SMEAR 中设置 METHOD 为 FERMI_DIRAC，并适当增大 ELECTRONIC_TEMPERATURE（如 300-1000 K），可以帮助电子填充更平滑，从而加快收敛。

总结：SCF 收敛的难度通常与体系结构、基组选择、计算精度以及 SCF 迭代策略密切相关。可综合利用上述方法进行逐步优化，以达到更稳定的 SCF 收敛。

MOTION

MOTION 控制原子/离子运动（类似于 VASP 中的离子步）。无论是几何优化、晶胞优化、AIMD 还是过渡态搜索，都需要仔细地设置 MOTION 板块。

几何优化

使用 CP2K 进行几何优化时，需要设置 GLOBAL/RUN_TYPE 为 GEO_OPT。 CP2K 中几何优化分为两种，一种是正常的能量极小化优化，一种是使用 dimer 算法进行过渡态优化，默认前者。能量极小化优化有三种算法，分别是 CG、BFGS 和 LBFGS。 CG 算法是最稳定的算法，但计算速度相对较慢；BFGS 算法效率最高也最常用，计算中需要对 Hessian 矩阵进行对角化，如果初始结构不合理，BFGS 算法容易出问题； LBFGS 算法效率和 BFGS 类似，同时稳定性也很好，适合超大体系。对于一般的几何优化，推荐使用 LBFGS 算法。

注意，使用 dimer 进行过渡态优化时，只能使用 CG 算法。

几何优化（能量极小化）的 MOTION 如下所示：

如果在几何优化过程中需要固定部分原子，可以在 MOTION/CONSTRAINT 设置 FIXED_ATOMS 选项。在几何优化过程中，有可能中途偶尔出现几步 SCF 不收敛。此时 CP2K 程序还会继续算下去，只要最终几何优化收敛时 SCF 是收敛的状态就行了。几何优化每一步输出的结构都写在 cp2k-pos-1.xyz 中，使用命令 grep = cp2k-pos-1.xyz 查看每一步的能量，使用命令 grep Max\\.\\ g cp2k.out 查看每一步的受力。

晶胞优化

晶胞优化是在几何结构优化的基础上，同时优化晶胞参数。首先，GLOBAL/RUN_TYPE 设置为 CELL_OPT，然后在 MOTION 部分需要同时设置 CELL_OPT 和 GEO_OPT 的参数，如下所示：

此外，在 FORCE_EVAL/STRESS_TENSOR 需要设置为 ANALYTICAL。使用命令 grep CELL cp2k.out 查看晶胞优化结果。

AIMD

计算 AIMD 是 CP2K 的优势项目，特别是 CP2K 自带的 CSVR 热浴能更好地控温。

CP2K 计算 AIMD 时，可以在 FORCE_EVAL/DFT/SCF/PRINT 中关掉波函数输出（即设置 RESTART OFF），避免I/O浪费时间。

在 VASP 和 LAMMPS 中可以设置初温与末温，而 CP2K 对于变温过程不是那么友好。 CP2K 变温的方式有两种，一种是准静态变温，另一种是使用 ANNEALING 参数控制变温。准静态变温是将某一温度下平衡的结构输出，用于下一温度的模拟。 ANNEALING 参数是通过设置温度标度因子，每步使用该因子重新进行温度标度。当 ANNEALING > 1 时为升温过程，反之为降温，其默认值为1。使用该命令时不能使用热浴，这意味着仅可选择 NVE 或 NPE 系综。一般操作是先在 NPT 或 NVT 系综下平衡结构，得到 .restart 文件后改为 NVE 系综通过 ANNEALING 参数变温，运行的步数与你设置的标度因子有关，最后再使用 .restart 文件在 NPT 或 NVT 系综下平衡。

NEB 过渡态搜索

准备初、末态的结构文件，分别使用 CP2K 进行固定晶格的结构优化（即 GEO_OPT）。优化好的初态和末态分别命名为 ini.xyz 和 fin.xyz。 GLOBAL/RUN_TYPE 设置为 BAND，FORCE_EVAL 和前述基本一致，其中 SUBSYS/TOPOLOGY 的结构文件用 ini.xyz 即可。然后在 MOTION 部分设置 BAND 参数，如下所示：

输出文件中，ener 文件输出 NEB 每步的能量和原子间距变化， -pos-Replica*.xyz 文件输出各点的原子坐标优化过程。

CP2K 输出结果分析

使用 Diag 算法时，grep "SCF run converged" cp2k.out 查看是否收敛；grep "/Diag." cp2k.out 查看自洽过程收敛趋势。

使用 OT 算法时，grep "outer SCF loop converged" cp2k.out 查看是否收敛；grep "OT DIIS" cp2k.out 查看自洽过程收敛趋势。

另外，无论是 Diag 还是 OT，grep "Total energy" cp2k.out | tail -1 查看总能； grep -A12 "ATOMIC FORCES" cp2k.out # -A 后面等于晶胞原子数加 4 查看原子受力。

使用 time.sh 提取计算任务用时：