Jing's Notes
0%

【CP2K】04.几何优化

发表于 更新于 分类于 ，
Valine： 本文字数： 11k 阅读时长 ≈ 10 分钟
以H2O（不更改单元尺寸）为例的几何优化，介绍&MOTION部分（&GEO_OPT、&CONSTRAINT、&CELL_OPT模块），&FORCE_EVAL部分（&MM分子动力学部分）

输入文件

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
&GLOBAL
  PROJECT H2O
  RUN_TYPE GEO_OPT	#修改任务类型
  PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
  METHOD QS
  &SUBSYS
    &CELL							#修改晶格向量单位
      ABC 12.4138 12.4138 12.4138
    &END CELL
    &COORD							#修改坐标
      O      12.235322       1.376642      10.869880
      H      12.415139       2.233125      11.257611
      H      11.922476       1.573799       9.986994
    &END COORD
    &KIND H							#设定基组
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q1
    &END KIND
    &KIND O							#设定基组
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q6
    &END KIND
  &END SUBSYS
  &DFT
    BASIS_SET_FILE_NAME ./BASIS_SET
    POTENTIAL_FILE_NAME ./POTENTIAL
    &QS
      EPS_DEFAULT 1.0E-7		#相比计算单点可以增大
    &END QS
    &MGRID						#修改MGRID
      CUTOFF 200			
      NGRIDS 4				
      REL_CUTOFF 30
    &END MGRID
    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-05		#修改EPS_SCF
      MAX_SCF 200			#修改MAX_SCF
      &DIAGONALIZATION T
        ALGORITHM STANDARD
      &END DIAGONALIZATION
      &MIXING T
        ALPHA 0.5
        METHOD PULAY_MIXING		#用Pulay混合法
        NPULAY 5
      &END MIXING
      &PRINT					#不重启SCF
        &RESTART OFF
        &END RESTART
      &END PRINT
    &END SCF
    &XC
      &XC_FUNCTIONAL PADE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
&END FORCE_EVAL
&MOTION
  &GEO_OPT
    TYPE MINIMIZATION
    MAX_DR    1.0E-03
    MAX_FORCE 1.0E-03
    RMS_DR    1.0E-03
    RMS_FORCE 1.0E-03
    MAX_ITER 200
    OPTIMIZER CG
    &CG
      MAX_STEEP_STEPS  0
      RESTART_LIMIT 9.0E-01
    &END CG
  &END GEO_OPT
  &CONSTRAINT
    &FIXED_ATOMS
      COMPONENTS_TO_FIX XYZ
      LIST 1
    &END FIXED_ATOMS
  &END CONSTRAINT
&END MOTION

&MOTION

定义了一组与原子核运动有关的工具

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
&MOTION
  &GEO_OPT
    TYPE MINIMIZATION
    MAX_DR    1.0E-03
    MAX_FORCE 1.0E-03
    RMS_DR    1.0E-03
    RMS_FORCE 1.0E-03
    MAX_ITER 200
    OPTIMIZER CG
    &CG
      MAX_STEEP_STEPS  0
      RESTART_LIMIT 9.0E-01
    &END CG
  &END GEO_OPT
  &CONSTRAINT
    &FIXED_ATOMS
      COMPONENTS_TO_FIX XYZ
      LIST 1
    &END FIXED_ATOMS
  &END CONSTRAINT
&END MOTION

&GEO_OPT

GEO_OPT小节仅适用于晶格尺寸不变的计算：

  • TYPE：指定要执行的几何优化类型
    • MINIMIZATION：执行几何最小化。默认。
    • TRANSITION_STATE：执行过渡状态优化。
  • 设置是否达到优化几何形状的标准：
    • MAX_DR：先前几何优化迭代中原子位移的最大公差。默认为Bohr
    • MAX_FORCE：对原子力的最大公差，默认为Bohr / Hartree
    • RMS_DR：先前几何优化迭代中原子位移的均方根的公差，默认为Bohr
    • RMS_FORCE：对原子力的均方根公差，默认为Bohr / Hartree
  • MAX_ITER：置最大几何优化迭代次数
  • OPTIMIZER：设置查找固定点的算法，同MINIMIZER
    • CG：共轭梯度，robust minimizer（取决于线搜索），也适用于大型系统
    • BFGS：默认。最有效的最小化器，但仅适用于“小型”系统，因为它依赖于完整的Hessian矩阵的对角化
    • LBFGS：BFGS的内存有限变体，适用于大型系统。微调效果不佳，但功能更强大。

&CG

设置共轭梯度算法的选项：

  • MAX_STEEP_STEPS：开始共轭梯度优化之前，最陡下降步骤的最大数量。默认是0，即不会执行最陡峭的下降步骤。
  • RESTART_LIMIT：两个连续搜索方向之间的角度的余弦值。如果CG优化过程中的角度小于对应于RESTART_LIMIT的角度，则会重置CG并执行最陡下降步骤。默认为0.9。

&CONSTRAINT

为原子运动添加约束

&FIXED_ATOMS

此部分用于约束整体原子位置（X，Y，Z）。如果指定了约束，则TARGET的值将被视为运行开始时的坐标值，或者是&FIX_ATOM_RESTART节中的相应值。

  • COMPONENTS_TO_FIX：指定将在本节中指定的原子的哪些成分（X，Y，Z或组合）受到约束/约束。默认是XYZ
  • LIST：指定要冻结的原子列表，对应于&FORCE_EVAL中的&SUBSYS的&COORD中的顺序

&CELL_OPT

设置用于优化模拟晶胞的环境。有两种可能的方案：（1）零温度优化； （2）有限温度优化。

  • KEEP_SYMMETRY：保持要求的初始细胞对称性（例如在细胞优化过程中）。初始对称性必须在＆CELL部分中指定。默认是.FALSE.

输出文件

会输出以下文件：

  • H2O.out：包含主要输出
  • H2O-pos-1.xyz：包含xyz文件格式中每个几何优化步骤的原子坐标轨迹。
  • H2O-1.restart：类似于输入文件，包含水分子的最新原子坐标。可以直接用于重启任务，也可以直接用其作为模板，使用其中的原子结构进一步计算
  • H2O-1.restart.bak-1、H2O-1.restart.bak-2、H2O-1.restart.bak-3：H2O-1.restart.bak-*是具有从先前的1、2和3几何优化迭代获得的原子坐标的备份重新启动文件。H2O-1.restart.bak-1应该与 H2O-1.restart相同。

.out文件

  • 每一步几何优化都会有如下信息：

    • 在几何优化步骤1结束时，系统的总能量为-17.1643447508（Ha），并且未达到优化几何的标准。因此，迭代将继续进行，直到所有条件变为“ YES”为止。
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    --------  Informations at step =     1 ------------
     Optimization Method        =                   SD
     Total Energy               =       -17.1643447508
     Real energy change         =        -0.0006776683
     Decrease in energy         =                  YES
     Used time                  =               90.837

     Convergence check :
     Max. step size             =         0.0336570168
     Conv. limit for step size  =         0.0010000000
     Convergence in step size   =                   NO
     RMS step size              =         0.0168136889
     Conv. limit for RMS step   =         0.0010000000
     Convergence in RMS step    =                   NO
     Max. gradient              =         0.0182785685
     Conv. limit for gradients  =         0.0010000000
     Conv. for gradients        =                   NO
     RMS gradient               =         0.0091312361
     Conv. limit for RMS grad.  =         0.0010000000
     Conv. for gradients        =                   NO
    ---------------------------------------------------
    • 结束时所有标准都已满足:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    --------  Informations at step =    11 ------------
     Optimization Method        =                   SD
     Total Energy               =       -17.1646204766
     Real energy change         =        -0.0000000529
     Decrease in energy         =                  YES
     Used time                  =               49.893

     Convergence check :
     Max. step size             =         0.0003393150
     Conv. limit for step size  =         0.0010000000
     Convergence in step size   =                  YES
     RMS step size              =         0.0001493298
     Conv. limit for RMS step   =         0.0010000000
     Convergence in RMS step    =                  YES
     Max. gradient              =         0.0001787448
     Conv. limit for gradients  =         0.0010000000
     Conv. in gradients         =                  YES
     RMS gradient               =         0.0000786642
     Conv. limit for RMS grad.  =         0.0010000000
     Conv. in RMS gradients     =                  YES

  • 最终的Kohn-Sham能量可以在输出的末尾获得：

    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
    46
    47
    48
    49
    50
    *******************************************************************************
    ***                    GEOMETRY OPTIMIZATION COMPLETED                      ***
    *******************************************************************************

                       Reevaluating energy at the minimum

    Number of electrons:                                                          8
    Number of occupied orbitals:                                                  4
    Number of molecular orbitals:                                                 4

    Number of orbital functions:                                                 23
    Number of independent orbital functions:                                     23

     Parameters for the always stable predictor-corrector (ASPC) method:

     ASPC order: 3

     B(1) =   3.000000
     B(2) =  -3.428571
     B(3) =   1.928571
     B(4) =  -0.571429
     B(5) =   0.071429

    Extrapolation method: ASPC


    SCF WAVEFUNCTION OPTIMIZATION

     Step     Update method      Time    Convergence         Total energy    Change
     ------------------------------------------------------------------------------
        1 Pulay/Diag. 0.50E+00    0.5     0.00005615       -17.1646204762 -1.72E+01
        2 Pulay/Diag. 0.50E+00    1.0     0.00000563       -17.1646347711 -1.43E-05

     *** SCF run converged in     2 steps ***


     Electronic density on regular grids:         -8.0000016293       -0.0000016293
     Core density on regular grids:                7.9999992554       -0.0000007446
     Total charge density on r-space grids:       -0.0000023739
     Total charge density g-space grids:          -0.0000023739

     Overlap energy of the core charge distribution:               0.00000004555422
     Self energy of the core charge distribution:                -43.83289054591484
     Core Hamiltonian energy:                                     12.82175605770555
     Hartree energy:                                              17.97395116120845
     Exchange-correlation energy:                                 -4.12745148966141

     Total energy:                                               -17.16463477110803

    ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):              -17.164634771108034

例子

NaCl团簇的几何优化

输入文件模板template.inp

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
@SET NREP MY_SUPERCELL
@SET OPTIMIZER LBFGS # BFGS

&FORCE_EVAL
  METHOD Fist
  &MM
    &FORCEFIELD
      &CHARGE
        ATOM Na
        CHARGE +1.000
      &END CHARGE
      &CHARGE
        ATOM Cl
        CHARGE -1.000
      &END CHARGE
      &NONBONDED
        &BMHFT
          map_atoms NA NA
          atoms NA NA
          RCUT 10.0
        &END BMHFT
        &BMHFT
          map_atoms NA CL
          atoms NA CL
          RCUT 10.0
        &END BMHFT
        &BMHFT
          map_atoms CL CL
          atoms CL CL
          RCUT 10.0
        &END BMHFT
      &END NONBONDED
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE spme
        ALPHA .35
        GMAX 12*${NREP}
        O_SPLINE 6
      &END EWALD
    &END POISSON
  &END MM
  &SUBSYS
    &CELL
      #ABC 5.620 5.620 5.620
      ABC 2*5.620 2*5.620 2*5.620
      MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME NaCl.pdb
      COORDINATE PDB
      CONN_FILE_FORMAT OFF
      MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
    &END TOPOLOGY
  &END SUBSYS
&END FORCE_EVAL

&GLOBAL
  PROJECT NaCl
  RUN_TYPE GEO_OPT
&END GLOBAL

&MOTION
  &GEO_OPT
    OPTIMIZER ${OPTIMIZER}
  &END
&END MOTION

  • 可以通过@SET key XX来设定变量

  • 通过脚本来取代MY_SUPERCELL设定一系列不同值的NREP的任务，并通过脚本提交任务

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    for ii in 2 4 6 8 10 12
    do
       mkdir ${ii}
       sed -e "s/MY_SUPERCELL/${ii}/g" template.inp > ${ii}/NaCl_supercell_${ii}.inp
       sed -e "s/filename/NaCl_supercell_${ii}/g" cp2k.sub > ${ii}/cp2k.sub
       cp NaCl.pdb ${ii}/NaCl.pdb
       cd ${ii}
       qsub cp2k.sub
       cd ..
    done

  • .pdb文件

1
2
3
4
5
6
7
8
9
10
11
12
REMARK   4 NaCl COMPLIES WITH FORMAT V. 2.0 
CRYST1    5.620    5.620    5.620  90.00  90.00  90.00 P 1         
HETATM    1 Na   NMO     0       0.000   0.000   0.000  1.00  0.00      NMO Na
HETATM    2 Cl   CMO     0       2.810   2.810   2.810  1.00  0.00      CMO Cl
HETATM    3 Na   NMO     0       0.000   2.810   2.810  1.00  0.00      NMO Na
HETATM    4 Cl   CMO     0       2.810   0.000   0.000  1.00  0.00      CMO Cl
HETATM    5 Na   NMO     0       2.810   0.000   2.810  1.00  0.00      NMO Na
HETATM    6 Cl   CMO     0       0.000   2.810   0.000  1.00  0.00      CMO Cl
HETATM    7 Na   NMO     0       2.810   2.810   0.000  1.00  0.00      NMO Na
HETATM    8 Cl   CMO     0       0.000   0.000   2.810  1.00  0.00      CMO Cl
TER       9      UNK     0 
END

&FORCE_EVAL

  • METHOD中使用FIST，使用分子动力学（MM）方法

&MM

  • 使用分子动力学的相关设定

&FORCEFIELD

指定有关如何为经典计算正确设置force_field的信息的部分。

&CHARGE

指定MM原子的电荷

  • ATOM：定义电荷的原子种类。
  • CHARGE：以电子电荷单位定义MM原子的电荷。
&NONBONDED

指定非键的输入参数

&BMHFT

指定BMHFT势能类型的输入参数。只适用于Na / Cl对。函数形式为

  • MAP_ATOMS：定义仅在Na和Cl时内部定义BMHFT非键势的种类。
  • ATOMS：定义BMHFT非键势涉及的原子种类
  • RCUT：定义BMHFT电位的截止参数，默认是4.12758223 angstrom。
&LENNARD-JONES

指定LENNARD-JONES是能类型的输入参数。函数形式:。

  • ATOMS：定义涉及非成键势的原子类型
  • EPSILON：定义LJ势的EPSILON参数，单位为是K_e
  • SIGMA：定义LJ电位的SIGMA参数，单位是Å
  • RCUT：定义LJ电位的截止参数，默认是10Å

&POISSON

设置泊松解析器。

&EWALD

控制CLASSICAL MM静电。

  • EWALD_TYPE：要执行的ewald类型
    • EWALD：默认，标准的基于非fft的ewald
    • NONE：无标准实空间库仑电势与非键合贡献一起计算
    • PME：使用fft插值的粒子网格
    • SPME：使用beta-Euler样条曲线的平滑粒子网格（推荐）

  • ALPHA：与Ewald相关的alpha参数（EWALD | PME | SPME）。默认为0.35 angstrom^-1。

    • 对于小型系统，建议为3.5/rcut。
    • 对于紧密键合，建议值为1.0。
    • 需要调整alpha，rcut和gmax以获得ewald的O（）缩放比例。

  • GMAX：网格点数（SPME和EWALD）。如果指定一个数字，则网格上所有三个方向都使用相同的点数。如果给出三个数字，则每个方向可以具有不同的点数。点数需要可FFTable（取决于使用的库），对于EWALD，使用奇数。最佳数量取决于Alpha和晶胞的大小。通常每埃1点。

  • O_SPLINE：beta-Euler样条的顺序（仅SPME）。默认为6。

NaCl的晶胞优化

输入文件模板template.inp

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
@SET SCALE_FACTOR MY_SCALING
@SET NREP 1
@SET OPTIMIZER BFGS # LBFGS
@SET CUTOFF MY_CUTOFF
@SET SAFTEY_CUTOFF 1.1

&FORCE_EVAL
  METHOD QS
  &DFT
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &MGRID
      CUTOFF ${CUTOFF}
      REL_CUTOFF 60
    &END MGRID
    &QS
      EPS_DEFAULT 1.0E-12
    &END QS
    &SCF
      SCF_GUESS RESTART
      &OT ON
        MINIMIZER DIIS
      &END OT
    &END SCF
    &XC
      &XC_FUNCTIONAL Pade
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC 5.620*${SCALE_FACTOR} 5.620*${SCALE_FACTOR} 5.620*${SCALE_FACTOR}
      MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
      #&CELL_REF
      #  ABC 5.620*${SAFETY_FACTOR} 5.620*${SAFETY_FACTOR} 5.620*${SAFETY_FACTOR}
      #  MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
      #&END
    &END CELL
    &COORD
		scaled
		Na 0.000   0.000   0.000
		Cl 0.500   0.500   0.500
		Na 0.000   0.500   0.500
		Cl 0.500   0.000   0.000
		Na 0.500   0.000   0.500
		Cl 0.000   0.500   0.000
		Na 0.500   0.500   0.000
		Cl 0.000   0.000   0.500
    &END
    &TOPOLOGY
      MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP}
    &END TOPOLOGY
    &KIND Na
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PADE-q9
    &END KIND
    &KIND Cl
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PADE-q7
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

&GLOBAL
  PROJECT NaCl
  RUN_TYPE ENERGY
&END GLOBAL

&MOTION
  &GEO_OPT
    OPTIMIZER ${OPTIMIZER}
  &END
  &CELL_OPT
    KEEP_SYMMETRY
  &END
&END MOTION
1
2
3
4
5
6
7
CUTOFF="280"

for ii in 0.800 0.825 0.850 0.875 0.900 0.925 0.950 0.975 1.000 1.025 1.050 1.075 1.100
do
   sed -e "s/MY_SCALING/${ii}/g" NaCl_QS.inp > temp.inp
   sed -e "s/MY_CUTOFF/${CUTOFF}/g" temp.inp > input_${ii}.inp
done
有用可戳(●ˇ∀ˇ●)