Quick Start of CP2K

本教程基于 HPC of zhenggroup 集群环境

一、CP2K介绍

关于CP2K

  1. CP2K是一种量子化学和固态物理第一性原理计算软件包，可以执行固态，液态，分子，和生物系统的计算模拟。
  2. CP2K是开源免费程序，使用混合高斯和平面波基组GPW和GAPW的DFT方法比平面波基组计算更快。
  3. CP2K用到好多非常新和前沿的算法
    OT(Orbital Transformation)，
    ASPC(Always stable predictor corrector)，
    ADMM(Auxiliary Density Matrix Methods)，
    PIMD(pathintegral Molecular Dynamics)，
    adaptive buffered QM/MM
  4. 相关文章 https://www.cp2k.org/science
  5. 催化：ACS Catal. 2019, 9 (9), 7876-7887

CP2K的应用和局限

应用: DFT,AIMD,(QM/MM),NEB,MD https://www.cp2k.org/features
局限: 1. 对导体计算较慢，OT算法算有带隙的材料可以很快收敛。导体一般用较慢的对角化方法。
2. 磁性体系提前指定自旋多重度。
3. 高斯基组带来的基组重叠误差（Basis Set Superposition Error, BSSE），不适合计算结合能。VASP使用平面波基组不存在此问题。
4. 能量计算使用VASP更准确

二、CP2K使用

输入文件+环境准备

*.inp,(coord.inc/xyz),CP2K.pbs

在提交脚本cp2k.pbs内添加以下内容

1
2
3
4
source /public/software/gcc-9.3.0/gcc-9.3.0_env.sh
source /public/software/intel2020u2//intel2020u2_env.sh
export OMP_NUM_THREADS=1
mpirun  -n 48 /public/software/cp2k-8.1/exe/local/cp2k.popt  -i cp2k.inp > cp2k.out

cp2k.inp 常用参数说明

CP2K手册https://manual.cp2k.org/trunk/index.html
关键词以section和subsection的形式套娃，每一个section都是以【&section_name】开头的，以【&END section_name 】结尾的。(顺序随意，嵌套不能乱)每一行只写一个关键词，关键词后面写参数，大小写和空格不敏感。

通用参数GLOBAL

1
2
3
4
5
&GLOBAL
  PROJECT cp2k                              # 任务名称
  RUN_TYPE GEO_OPT                          # 计算类型 ENERGY,ENERGY_FORCE,GEO_OPT,CELL_OPT,VIBRATIONAL_ANALYSIS,BAND,MD,PINT
  PRINT_LEVEL LOW                           # 输出等级
&END GLOBAL

体系的力和能量FORCE_EVAL

  1. SUBSYS(subsystem): coordinates, topology, molecules and cell
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    &FORCE_EVAL
      METHOD Quickstep                          # 使用DFT的GPW方法进行计算原子受力
      &SUBSYS
    # --------------- Structure --------------- #
        &TOPOLOGY
          &CENTER_COORDINATES                   # 结构放在盒子中间
          &END CENTER_COORDINATES
          COORD_FILE_FORMAT xyz
          COORD_FILE_NAME  ./coord.xyz
        &END TOPOLOGY
    # ------------ Lattice vectors ------------ #
        &CELL                                   # 也可输入 `ABC 10 11 13` 代表正交晶系，单位埃
          A    8.42    0.00    0.00
          B    4.21    7.29    0.00
          C    0.00    0.00   19.32
          PERIODIC XYZ                          # NONE,X,XY,XYZ,XZ,Y,YZ,Z  XYZ代表三个方向周期性计算
        &END CELL
    # ------------ Basis, Potential ----------- #
        &KIND Au                                # 元素种类
          ELEMENT Au                            
          BASIS_SET DZVP-MOLOPT-SR-GTH-q11      # 基组设置(对应文件里存在的条目)
          POTENTIAL GTH-PBE                     # 赝势设置
        &END KIND
    # ----------------------------------------- #
      &END SUBSYS
      ...
    &END FORCE_EVAL
    其他输入原子坐标方式
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    # Method 1：inc --------------------------- #
    &COORD  
        @INCLUDE 'coord.inc'  #xyz文件去掉前两行
    &END COORD

    # Method 2：Cartesian coordinates --------- #
    &COORD 
    C          -3.36178719   11.43419295    1.72360000
    H           3.11299203   11.00326708    1.72360000
    ...
    &END COORD
  2. DFT: Parameter needed by LCAO DFT programs (Quickstep) 关于自洽循环(电子步)的计算参数
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    &FORCE_EVAL
      &DFT
        BASIS_SET_FILE_NAME  BASIS_MOLOPT       # 基组文件路径
        POTENTIAL_FILE_NAME  POTENTIAL          # 赝势文件路径
    # --------------- Quickstep --------------- #
        &QS                                     # 
          EPS_DEFAULT 1.0E-12                   # 高精度计算时，EPS_DEFAULT 1.0E-14
        &END QS                                 
    # ----------------- Grid ------------------ #
        &MGRID                                  # 
          NGRIDS 4                              # 多重网格叠加做FFT, 值越大，算得越快，越不精确
          CUTOFF 500                            # 整体网格精度的最高值 https://manual.cp2k.org/trunk/methods/dft/cutoff.html，尽量>800
          REL_CUTOFF 50                         # 控制有多少网格点落到最精细的级别
        &END MGRID 
    # ----------------- SCF ------------------- #
        &SCF                                    
          SCF_GUESS ATOMIC                      # 初猜/可以改成RESTART
          EPS_SCF 1.0E-6                        # 能量收敛精度
          MAX_SCF 100                           # SCF循环圈数
          # ------ OT, DIAGONALIZATION -------- #
          &OT T                                 # 使用OT(Orbital Transformation)优化波函数
          MINIMIZER CG                          # 共轭梯度法(最可靠)(CG)/迭代子空间直接反演法(速度快)(DIIS)/布罗伊登法(BROYDEN)  
          LINESEARCH 3PNT                       # 搜索算法
          PRECONDITIONER FULL_SINGLE_INVERSE    # 适用于大型系统；难收敛改Full All
          &END OT 
          # &DIAGONALIZATION ON                 # 对角化方法，不如OT快
          #   ALGORITHM STANDARD 
          # &END DIAGONALIZATION 
          # ----------------------------------- #
        &END SCF                                
    # --- eXchange and Correlation potential -- #
        &XC                                     # 交换-关联密度泛函，要与基组和赝势的选择一致。
          &XC_FUNCTIONAL PBE 
          &END XC_FUNCTIONAL 
        &END XC
    # ----------------- Print ----------------- #
        &PRINT                                  # 输出参数 https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/PRINT.html 
          &FORCES ON                            # 输出力
          &END FORCES
        &END PRINT
    # ----------------------------------------- #
      &END DFT
      ...
    &END FORCE_EVAL
  3. 其他常用&DFT参数
&DFT 参数 设置 含义
BASIS_SET_FILE_NAME BASIS_MOLOPT, BASIS_MOLOPT_UCL 基组文件，前面可加路径${DATAPATH} 需要@SET DATAPATH ….
POTENTIAL_FILE_NAME GTH_POTENTIALS 赝势文件
UKS T 开壳层
ROK T 限制性开壳层
MULTIPLICITY 整数 闭壳层1，开壳层双重态2，三重态3
WFN_RESTART_FILE_NAME ./RESTART.wfn（波函数路径） 需要同时开启&DFT/&SCF中的SCF_GUESS RESTART
PLUS_U_METHOD 见手册 DFT+U （U值在 &KIND 里控制）
SURFACE_DIPOLE_CORRECTION T 表面，偶极校正 <=> IDIPOL = T
SURF_DIP_DIR Z Z方向 <=> IDIPOL = 3

  1. OT和Diagonalization比较
    OT：速度快，没有高斯展宽，对于金属系统收敛性差，不支持DFT+U，以及KPOINTS
    Dia：计算量大，要设置高斯展宽，支持DFT+U和KPOINTS

  2. PIRNT
    &FORCE_EVAL/ 中添加&PRINT/&FORCES ON/&END FORCES，输出力。
    &FORCE_EVAL/&DFT/&SCF 中添加&PRINT/&RESTART FILENAME =RESTART.wfn/&END RESTART，输出波函数。

  3. 基组
    计算精度，SZV < DZVP < TZVP < TZV2P < TZV2PX
    VP(Valence Polarized同时添加极化函数)
    VASP平面波基组(PW)，不同的体系设定相同的平面波截断能（ENCUT），那么用到的基组都是一样的。
    高斯基组计算AB体系和A体系使用的基组不一样，AB中B的部分基组也贡献于描述A原子。(基组重叠误差，基组越大，误差越小，SR相比与普通的基组，速度略有提升，BSSE误差增大)
    AB体系算出来的能量被低估。(计算吸附能、结合能等都会偏负，E(AB) – E(A) – E(B))

  4. 对角化K点设置

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    &FORCE_EVAL
      &DFT
    # ----------------KPOINTS------------------ #
        &KPOINTS
          SCHEME MONKHORST-PACK 3 3 3           
          SYMMETRY T                            # 识别晶胞对称性
          VERBOSE T                             # 输出
          FULL_GRID T
        &END KPOINTS
    # ----------------------------------------- #
      &END DFT
    &END FORCE_EVAL

原子运动参数MOTION

结构和晶胞优化

  1. 几何结构优化：GEO_OPT
    cp2k manual - geometry optimization
    ex:H2O
    &GLOBAL中RUN_TYPE设置GEO_OPT

    1
    2
    3
    4
    5
    6
    7
    8
    9
    &MOTION
    # ------------ Geometry Optimizer --------- #
      &GEO_OPT
        MAX_ITER 400
        OPTIMIZER BFGS                          # LBFGS大体系，BFGS针中小体系，CG更为Robust
        MAX_FORCE 6.0E-4
      &END GEO_OPT
    # ----------------------------------------- #
    &END MOTION

    输出文件：(以H2O为例) H2O.out, H2O-pos-1.xyz, H2O-1.restart, H2O-1.restart.bak-1,2,3
    H2O-1.restart.bak-*代表计算结束前*次迭代结果，..bak-1和H2O-1.restart相同
    终止后续算：使用H2O-1.restart作为input文件重新提交任务

  2. 晶胞优化：CELL_OPT
    cp2k manual - cell optimization
    &GLOBAL中RUN_TYPE设置CELL_OPT

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    &MOTION
    # -------------- Cell Optimizer ----------- #
      &CELL_OPT
        MAX_ITER 400
        OPTIMIZER CG                            # LBFGS大体系，BFGS针中小体系，CG更为Robust
        KEEP_ANGLES                             # 固定角度
        KEEP_SYMMETRY T                         # 固定对称性   
        # CONSTRAINT Z                          # 固定Z方向
        TYPE DIRECT_CELL_OPT                    # 同时优化晶胞和里面的位置(默认)
        MAX_FORCE 6.0E-4
      &END GEO_OPT
    # ----------------------------------------- #
    &END MOTION

  3. 坐标固定：CONSTRAINT

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    &MOTION
    # ------------- Constraints --------------- #
      &CONSTRAINT
        # ---------- Fixed atoms -------------- #
        &FIXED_ATOMS
          COMPONENTS_TO_FIX XYZ                 # 限制固定的方向，也可以设置X,XY等  
          LIST 1..54
          LIST 289..324
        &END FIXED_ATOMS
      &END CONSTRAINT
    # ----------------------------------------- #  
    &MOTION

频率计算

  1. &GLOBAL中RUN_TYPE设置VIBRATIONAL_ANALYSIS，输出信息精度调高
  2. SCF优化精度不变，同时添加&VIBRATIONAL_ANALYSIS
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    &VIBRATIONAL_ANALYSIS
      NPROC_REP 48                              # 总核数=节点数*核数
      DX 0.01                                   # 步长
      FULLY_PERIODIC F                          # 是否周期性
    # ----------------- Print ----------------- #
      &PRINT
        &MOLDEN_VIB                             # 生成振动模式的 MOLDEN 文件
        &END
        &CARTESIAN_EIGS                         # 输出输出频率计算中笛卡尔坐标系下的本征值和本征向量
        &END
        &PROGRAM_RUN_INFO                       # 记录程序运行的详细信息
          &EACH
            REPLICA_EVAL 1
          &END
        &END
      &END PRINT
    # ----------------------------------------- #
    &END VIBRATIONAL_ANALYSIS
  3. 自由能计算
    http://sobereva.com/552
    使用卢天的~/software/Shermo_2.6/Shermo计算ZPE，需要引用
    10.1016/j.comptc.2021.113249

过渡态计算

http://sobereva.com/44

  1. DIMER
    在势能放置两个相邻的点 Dimer 来寻找 TS，通过不断重复旋转移动dimer，直到找到能量最高点。
  2. CINEB
    在 IS 到 FS 之间插入一系列结构，共插入P-1个，反应为0，产物为P。优化不是对每个点孤立地优化，优化过程所有点一起运动。插点结构受到
    1.来自势能面的力（保留平行分量）；2.chain方向上的弹簧力（保留垂直分量）

CP2K-CINEB过程
&GLOBAL中RUN_TYPE设置BAND

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
&MOTION
# ----------------- Band ------------------ #
  &BAND
  NPROC_REP 28
  BAND_TYPE CI-NEB
  ALIGN_FRAMES F
  ROTATE_FRAMES F
  NUMBER_OF_REPLICA 6
  K_SPRING 0.08
    # ---------- Convergence -------------- #
    &CONVERGENCE_CONTROL
      MAX_DR 0.01
      MAX_FORCE 0.001
      RMS_DR 0.02
      RMS_FORCE 0.001
    &END
    &CI_NEB
      NSTEPS_IT 5
    &END
    &OPTIMIZE_BAND
      OPT_TYPE DIIS
      OPTIMIZE_END_POINTS F
    &END
    &REPLICA
      COORD_FILE_NAME ./1.xyz
    &END
    &REPLICA
      COORD_FILE_NAME ./2.xyz
    &END
    &REPLICA
      COORD_FILE_NAME ./3.xyz
    &END
    &REPLICA
      COORD_FILE_NAME ./4.xyz
    &END
    &REPLICA
      COORD_FILE_NAME ./5.xyz
    &END
    &PROGRAM_RUN_INFO
    &END
    &CONVERGENCE_INFO
    &END
  &END BAND
# ----------------------------------------- #
&END MOTION

CP2K-DIMER过程
http://bbs.keinsci.com/thread-23516-1-1.html

  1. 初猜选取dimer结构，得到coord.inc文件
  2. 设置优化方向 sed -i '1,2d' coord.inc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
METHOD DIMER
&DIMER
  DR 0.01
  ANGLE_TOLERANCE [deg] 4.0
  INTERPOLATE_GRADIENT
  &DIMER_VECTOR
     @INCLUDE 'DIMER_VECTOR'
  &END DIMER_VECTOR
  &ROT_OPT
    OPTIMIZER CG
    MAX_ITER 10

    MAX_DR 3.0E-3
    MAX_FORCE 6.0E-4
    &CG
      # MAX_STEEP_STEPS 15
      &LINE_SEARCH
        TYPE 2PNT
      &END LINE_SEARCH
    &END CG
  &END ROT_OPT
&END DIMER

杂化泛函（待更新）

单点计算（待更新）

AIMD

1
2
3
4
5
&QS
  EPS_DEFAULT 1.0E-14
  EXTRAPOLATION ASPC
  EXTRAPOLATION_ORDER 3
&END QS
  1. EXTRAPOLATION 控制波函数外推，在做结构优化/MD的时候，下一离子步的波函数用前几个离子步的波函数组合可以得到非常快的SCF收敛。建议用ASPC(Always stable predictor corrector)，注：不能和K点设置同时使用。
  2. EXTRAPOLATION_ORDER是ASPC extrapolation order。一般用默认值3。

这两个参数可以促进AIMD的时候的每一离子步的SCF收敛，默认打开。这也是cp2k算AIMD比较快的一个原因。

输出文件

(mermaid 模块还在测试)

graph TD
A[cp2k输入文件]-->B[GLOBAL]
A-->C[FORCE_EVAL]
A-->D[MOTION]
    C-->E[SUBSYS]
    C-->H[PRINT]
    C-->F[DFT]
        E-->G[CELL]
        E-->I[COORD]
        E-->J[TOPOLOGY]
        E-->K[KIND]
        F-->O[QS]
        F-->L[SCF]
        F-->N[MGRID]
        F-->P[XC]
        F-->T[PRINT]
            L-->M[OT]
            L-->Q[DIAGONALIZATION]

三、数据处理

脚本使用

*.out,(*.wfn.bak-*)wfn会不断覆盖，记录最后三步的波函数

四、学习资料

  1. CP2K的tutorials：
    https://www.cp2k.org/exercises
    https://www.cp2k.org/howto
  2. CP2K手册：http://manual.cp2k.org/trunk/
  3. CP2K的google group：https://groups.google.com/group/cp2k
  4. Zevan的博客：http://kenshin325.lofter.com/tag/cp2k
  5. 兰一的知乎专栏：https://zhuanlan.zhihu.com/cp2k-tutorial
  6. CP2K有用的工具： https://www.cp2k.org/tools
  7. python编写的cp2k输入文件生成器(最高支持5.1)：
    https://github.com/avishart/CP2K_Editor
  8. 难收敛http://sobereva.com/665
  9. http://wiki.cheng-group.net/wiki/software_usage/cp2k/cp2k/