0%

【CP2K】02.计算力和能量

输入文件

输入文件需要准备:

  • BASIS_SET和GTH_POTENTIALS:该文件包含CP2K可用于此计算的基组和赝势的参数,在cp2k/data/中有相应的基组和赝势,它们包含大多数常用的元素,用户需要针对给定的计算选择自己需要的文件
  • Si_bulk8.inp中主要包含&GLOBAL&FORCE_EVAL两部分:·
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
&GLOBAL
PROJECT Si_bulk8 #修改任务名字
RUN_TYPE ENERGY_FORCE #修改任务类型
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&SUBSYS
&KIND Si
ELEMENT Si
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q4
&END KIND
&CELL
A 5.430697500 0.000000000 0.000000000
B 0.000000000 5.430697500 0.000000000
C 0.000000000 0.000000000 5.430697500
&END CELL
&COORD
Si 0.000000000 0.000000000 0.000000000
Si 0.000000000 2.715348700 2.715348700
Si 2.715348700 2.715348700 0.000000000
Si 2.715348700 0.000000000 2.715348700
Si 4.073023100 1.357674400 4.073023100
Si 1.357674400 1.357674400 1.357674400
Si 1.357674400 4.073023100 4.073023100
Si 4.073023100 4.073023100 1.357674400
&END COORD
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-10
&END QS
&MGRID
NGRIDS 4
CUTOFF 300
REL_CUTOFF 60
&END MGRID
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-7
MAX_SCF 300
&DIAGONALIZATION ON
ALGORITHM STANDARD
&END DIAGONALIZATION
&MIXING T
METHOD BROYDEN_MIXING
ALPHA 0.4
NBROYDEN 8
&END MIXING
&END SCF
&END DFT
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL

&FORCE_EVAL

  • METHOD:使用哪种方法来计算力
    • QS:默认,同QUICKSTEP,即使用高斯和平面波(GPW)方法的密度泛函理论。与&DFT一起使用。
    • FIST:分子力学

&SUBSYS

1
2
3
&SUBSYS
……
&END SUBSYS

定义模拟晶胞和计算中原子的初始坐标。

&KIND 修改基组

定义计算中的元素,每个元素必须有一个&KIND

1
2
3
4
5
&KIND Si
ELEMENT Si
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q4
&END KIND
  • BASIS_SET:定义基组,如DZVP-GTH-PADE(双 ζ 带有极化基础的ζ,针对Geodecker-Teter-Hutter PADE LDA赝势进行了优化),同BASIS_SET_FILE_NAME
  • POTENTIAL:定义赝势GTH-PADE-q4(具有4个价电子的Geodecker-Teter-Hutter PADE LDA赝势),同POTENTIAL_FILE_NAME

可以在 BASIS_SET和GTH_POTENTIALS文件中查找到Si DZVP-GTH-PADESi GTH-PADE-q4 GTH-LDA-q4

&CELL

定义了计算中使用的模拟单位。Angstrom是默认的晶格向量单位, ABC是第一,第二和第三晶格(单元)矢量。

1
2
3
4
5
&CELL
A 5.430697500 0.000000000 0.000000000
B 0.000000000 5.430697500 0.000000000
C 0.000000000 0.000000000 5.430697500
&END CELL
  • ABC:指定晶格向量A,B和C的长度,它们定义正交晶胞的h矩阵的对角元素。对于非斜方细胞,可以通过ALPHA_BETA_GAMMA关键字指定角度ALPHA,BETA,GAMMA,也可以使用关键字ABC。约定是A沿X轴放置,B在XY平面上。
  • PERIODIC:指定将应用周期性边界条件(PBC)的方向。默认是XYZ
  • MULTIPLE_UNIT_CELL:指定已定义单元格在空间(X,Y,Z)中的重复次数(假设其为单位单元格)。此关键字仅影响晶胞。默认为1 1 1

&COORD

1
2
3
4
5
6
7
8
9
10
&COORD
Si 0.000000000 0.000000000 0.000000000
Si 0.000000000 2.715348700 2.715348700
Si 2.715348700 2.715348700 0.000000000
Si 2.715348700 0.000000000 2.715348700
Si 4.073023100 1.357674400 4.073023100
Si 1.357674400 1.357674400 1.357674400
Si 1.357674400 4.073023100 4.073023100
Si 4.073023100 4.073023100 1.357674400
&END COORD
  • 定义初始原子单位,默认输入格式为:

    1
    <ATOM_KIND> X Y Z
    • XYZ是埃的笛卡尔坐标
    • <ATOM_KIND>应该是与&KIND小节中元素的定义相对应的标签。
  • SCALED:使坐标输入相对于晶格矢量为分数,默认是.FALSE.,可以改为.TRUE.

  • UNIT:更改笛卡尔坐标的单位,默认是 ANGSTROM

可以用@INCLUDE 'coord.inc'读取文件简化输出:

1
2
3
&COORD
@INCLUDE 'coord.inc'
&END COORD

&TOPOLOGY

  • COORD_FILE_NAME:指定包含坐标的文件名。
  • COORDINATE:等同于COORD_FILE_FORMAT,设置读取坐标的方式:
    • off:在输入文件的&COORD部分中读取的坐标
    • XYZ:通过XYZ文件格式提供的坐标
    • PDB:通过PDB文件格式提供的坐标
  • CONN_FILE_NAME:指定包含分子连接性的文件名。
  • CONNECTIVITY:等同CONN_FILE_FORMAT,确定和生成分子的方法:
    • GENERATE:默认。使用简单的距离标准。
    • OFF:不产生分子。用于QS或定义不明确的系统。
  • MULTIPLE_UNIT_CELL:指定已定义单元格在空间(X,Y,Z)中的重复次数(假设其为单位晶胞)。此关键字仅影响坐标规范。默认为1 1 1

&PRINT

1
2
3
4
5
6
7
8
&PRINT
&ATOMIC_COORDINATES
&EACH
JUST_ENERGY 1
GEO_OPT 1
&END EACH
&END ATOMIC_COORDINATES
&END PRINT

&DFT

1
2
3
4
5
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
……
&END DFT

该部分控制自洽的Kohn-Sham密度泛函理论计算的所有方面。仅当且仅当FORCE_EVAL中的METHOD关键字设置为QUICKSTEP时,&DFT才有意义。

  • BASIS_SET_FILE_NAME设定基组文件
  • POTENTIAL_FILE_NAME设定赝势文件
  • WFN_RESTART_FILE_NAME:wavefunction重新启动文件的名称和路径,如果未指定文件,则默认为打开由wfn restart打印键生成的文件。

  • CHARGE:设定系统的总电荷,默认为0

  • MULTIPLICITY:与MULTIP相同,为2S+1,偶数个电子默认是1,奇数个默认为2

&QS

包含QUICKSTEP使用的常规控制参数。

1
2
3
&QS
EPS_DEFAULT 1.0E-10
&END QS
  • EPS_DEFAULT:QUICKSTEP中使用的所有公差设置默认值。默认为。在进行比较高精度计算时,需要将EPS_DEFAULT 设为1.0E-14。
  • EPS_ *:设置各类公差,覆盖EPS_DEFAULT值。
    • EPS_CORE_CHARGE:核电荷映射精度,默认为EPS_DEFAULT/100.0
    • EPS_GVG_RSPACE:实空间KS 矩阵元积分精度。默认为SQRT(EPS_DEFAULT)
    • EPS_PGF_ORB:重叠矩阵元精度,默认为SQRT(EPS_DEFAULT)
    • EPS_KG_ORB:使用Kim-Gordon 方法时的精度,默认为SQRT(EPS_DEFAULT)
  • METHOD:指定应采用的电子结构方法
    • GBW:默认。高斯和平面波方法。
  • EXTRAPOLATION:默认是ASPC。波函数的外推策略,并非所有的选项都适用于所有的模拟方法。推荐使用PSASPC,参见EXTRAPOLATION_ORDER
  • EXTRAPOLATION_ORDER:默认是3。PS或ASPC推断的阶数(通常是2-4)。更高的阶数可能带来更多的精度,但对于大型系统来说,也要付出一定的代价。在某些情况下,高阶外推法并不稳定,需要降低阶数。

&MGRID

1
2
3
4
5
&MGRID
NGRIDS 4
CUTOFF 300
REL_CUTOFF 60
&END MGRID

定义了QUICKSTEP计算中使用的积分网格。设置4个级别的多重网格。详见【CP2K】03.&MGRID中的CUTOFF和RELEL_CUTOFF

&XC 修改泛函

使用的交换相关密度泛函的参数

1
2
3
&XC
……
&END XC
&XC_FUNCTIONAL

注意选择的泛函要与所选择的基组和赝势一致

1
2
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
  • PADE:使用的PADE,即LDA泛函
  • PBE:使用PBE泛函
  • BLYP
&XC_GRID
  • XC_SMOOTH_RHO:用于xc计算的密度平滑,默认为NONE。还有NN10NN4NN50NN6SPLINE2SPLINE3
  • XC_DERIV:用于计算导数的方法。默认为PW,还有COLLOCATENN10_SMOOTHNN4_SMOOTHNN50_SMOOTHNN6_SMOOTHSPLINE2SPLINE2_SMOOTHSPLINE3SPLINE3_SMOOTH
&VDW_POTENTIAL 色散校正
1
2
3
4
5
6
7
8
9
&VDW_POTENTIAL
POTENTIAL_TYPE PAIR_POTENTIAL
&PAIR_POTENTIAL
PARAMETER_FILE_NAME dftd3.dat
TYPE DFTD3
REFERENCE_FUNCTIONAL PBE
R_CUTOFF 15
&END PAIR_POTENTIAL
&END VDW_POTENTIAL
  • POTENTIAL_TYPE:同名DISPERSION_FUNCTIONAL
    • NONE:默认,不使用色散
    • NON_LOCAL:非局域范德华密度泛函
    • PAIR_POTENTIAL:对势范德华密度泛函
  • &PAIR_POTENTIAL:关于计算色散的对电势的信息:
    • PARAMETER_FILE_NAME:参数文件的名称,可能包括一个路径,默认是DISPERSION_PARAMETERS
    • TYPE:类型,有DFTD2Grimme D2 方法、DFTD3Grimme D3 (zero damping) 方法和DFTD3(BJ)Grimme D3 (Becke-Johnson damping)方法。
    • REFERENCE_FUNCTIONAL:使用这个特定密度泛函的参数。
    • R_CUTOFF:势的范围。cutoff是这个值的2倍,默认是10.5835442Å。

&SCF

定义自洽场方法的设置

1
2
3
4
5
6
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-7
MAX_SCF 300
……
&END SCF
  • SCF_GUESS设置初始试验电子密度泛函。自洽场迭代中,初猜对于快速收敛比较重要。

    • ATOMIC:将使用原子电荷密度的重叠来生成初始密度。
    • RESTART:将RESTART文件用作初始猜测(如果不存在,则使用ATOMIC
  • EPS_SCF:设置电荷密度剩余量的公差。这将覆盖&QS小节中设置的EPS_DEFAULT值。 默认为。如果要
    进行比较高精度的计算,可以将SCF 收敛精度设置为1E-6。对于频率计算等精度要求更高的计算,SCF 收敛精度可以设置为1E-7
  • MAX_SCF:设置允许QUICKSTEP执行的每次基态能量计算的最大自洽循环数。
  • ADDED_MOS:每个自旋附加的MOS数量,默认为0,在使用对角化方法&DIAGONALIZATION的时候必须使用。
&DIAGONALIZATION

告诉代码使用传统的对角化方法来查找基态Kohn-Sham能量和电子密度。如果该节中有参数,就应当设置为ON.TRUE.T。对角化的替代方法是使用轨道变换(OT)方法,在这种情况下,用户应删除&DIAGONALIZATION块或改为OFF.FALSE.,然后添加&OT代替。

1
2
3
&DIAGONALIZATION  ON
ALGORITHM STANDARD
&END DIAGONALIZATION
  • ALGORITHM用于对角化的算法
    • STANDARD:标准对角化:LAPACK/SCALAPACK方法或Jacobi。
&MIXING

包含自一致性计算中与电荷混合相关的所有参数。可以通过.TRUE.T)或者.FALSE.F)来打开或者关闭,默认是打开。仅适用于传统的对角线化方法, OT方法使用不同的方法进行电荷混合。

1
2
3
4
5
&MIXING  T
METHOD BROYDEN_MIXING
ALPHA 0.4
NBROYDEN 8
&END MIXING
  • METHOD:设置混合方法

    • BROYDEN_MIXING:使用Broyden混合
    • PULAY_MIXING:使用Pulay混合法
  • ALPHA:设置混合参数,默认是0.4。比如0.4是将0.4的输出密度与0.6的输入密度混合,以在下一个SCF迭代中形成新的输入密度。

  • NBUFFER:实际混合方案存储的先前步骤数,等同于NBROYDENNPULAYNMULTISECANT,设置在Broyden混合算法、Pulay混合法中使用的历史记录数。

  • &SMEAR:不使用Smearing对于金属或间隙很小的系统,这可能会导致计算不稳定,并且由于电子占据函数的不连续性,自洽迭代可能永远不会收敛。因此可以添加

    1
    2
    3
    4
    &SMEAR ON
    METHOD FERMI_DIRAC
    ELECTRONIC_TEMPERATURE [K] 300
    &END SMEA
    • METHOD:使用的Smearing方法,FERMI_DIRAC即Fermi-Dirac smearing function
    • ELECTRONIC_TEMPERATURE:电子温度
  • ADDED_MOS:如果不使用该关键词,仅使用&SMEAR,会导致导带中的分子轨道被占用,必须告诉CP2K在计算中包括多余的,空的分子轨道,否则为了以减少计算成本将被省略。

    • 比如在计算中不要忽略10个最低的空分子轨道。应该注意的是,给定选定的基组,会有一个分子轨道最大数量,即Hamiltonian的特征向量数量。理论上,最大值应为在计算中生成的Hamiltonian的范围内。

      1
      ADDED_MOS 10
&PRINT

在SCF期间打印信息:

  • &RESTART在SCF期间控制MO重新启动文件的转储。默认情况下,保留简短的三次重新启动的历史记录。
&OT
1
2
3
4
5
6
&OT T
MINIMIZER DIIS
#MINIMIZER CG
#LINESEARCH 2PNT
PRECONDITIONER FULL_SINGL
&END OT

如果体系有较大带隙的,如为半导体或者绝缘体等,推荐使用OT 算法,收敛速度比较快。但是对预处理敏感。好的预处理器可能很昂贵。不能smearing或进行高级SCF mixing ,HOMO-LUMO 带隙很小或者几乎没有的系统,如金属系统,收敛性较差(此时应该使用&DIAGONALIZATION&Smear)。设置轨道变换(OT)方法的各种选项:

  • PRECONDITIONER
    • FULL_KINETIC:默认。S和T的Cholesky反转,快速构建,强大且相对较好,可用于超大型系统。
    • FULL_ALL:基于对角化的最有效的状态选择性预处理器,要求ENERGY_GAP参数低估HOMO-LUMO间隙。对于几乎所有系统,都建议使用此预处理器,但在大型系统中,make_preconditioner将主导总计算成本。
  • MINIMIZER
    • CG:默认,共轭梯度:最可靠,适用于困难的系统。如果合适的话,总能量应该在每个OT CG步骤减少。
    • DIIS:迭代子空间中的直接反演,也称为Pulay混合,是一种外推技术。由Peter Pulay在计算量子化学领域开发,旨在加速和稳定Hartree-Fock自洽场方法的收敛。不如CG可靠,但有时快约50%。
    • BROYDEN:如果CGDIIS都有问题,可以使用这个方法。
  • ENERGY_GAP:默认是-1。应该是能隙[a.u.](HOMO-LUMO)的估计值,并应在预处理中使用,尤其是对FULL_ALL预处理器有效,在这种情况下,它应该比能隙小。(可以为小数,例如0.002) 。 FULL_SINGLE_INVERSE将其作为下限(值小于0.05可能会导致稳定性问题)。通常,在初始猜测很差的情况下,heigher值会驯服预处理器。负值将取决于预调节器的类型,使选择留给CP2K。
&OUTER_SCF
1
2
3
4
&OUTER_SCF ON
MAX_SCF 10
EPS_SCF 1.0E-05
&END OUTER_SCF

OUTER_SCF 为加速收敛的一种方法,如果SCF 经过100 次优化依然没有收敛,则进入OUTER_SCF 过程,对前一次计算的波函数进行调整,重新进行SCF 迭代。每次OUTER_SCF 中优化的次数依然是100 次,比如最多可以进行5次OUTER_SCF,那么最多可以进行500 次SCF 计算。

  • EPS_SCF: outer SCF变量的目标梯度。默认为。内部循环的EPS_SCF还确定了外部循环可以达到的值,通常,外部循环的EPS_SCF必须小于或等于内部循环的EPS_SCF。
  • MAX_SCF:最大循环数

&PRINT

  • &FORCES,在主要的输出中输出原子受力。GLOBAL 部分的RUN_TYPE 必须设置为ENERGY_FORCE 或者GEO_OPT。如果要将受力信息存储到文件中,则需要设定文件的名称;否则受力信息将会打印到out 文件中。

  • &TOTAL_NUMBERS:控制打印的原子总数,种类……

输出文件

.out

主要的输出文件

  • 自洽的Kohn-Sham基态计算的典型输出。它说明正在使用对角化方法和Broyden电荷混合,并且花了10个Broyden混合步骤(每个步骤都包含对角化过程以获得波函数),以达到所需的自洽公差。

    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
    Number of electrons:                                                         32
    Number of occupied orbitals: 16
    Number of molecular orbitals: 16

    Number of orbital functions: 104
    Number of independent orbital functions: 104

    Extrapolation method: initial_guess


    SCF WAVEFUNCTION OPTIMIZATION

    Step Update method Time Convergence Total energy Change
    ------------------------------------------------------------------------------
    1 NoMix/Diag. 0.40E+00 0.6 0.75558724 -32.2320848878 -3.22E+01
    2 Broy./Diag. 0.40E+00 1.1 0.05667976 -31.1418135481 1.09E+00
    3 Broy./Diag. 0.40E+00 1.1 0.09691469 -31.1974003416 -5.56E-02
    4 Broy./Diag. 0.40E+00 1.1 0.00245608 -31.3378474040 -1.40E-01
    5 Broy./Diag. 0.40E+00 1.1 0.00235460 -31.3009654398 3.69E-02
    6 Broy./Diag. 0.40E+00 1.1 0.00007565 -31.2972158934 3.75E-03
    7 Broy./Diag. 0.40E+00 1.1 0.00009004 -31.2977293749 -5.13E-04
    8 Broy./Diag. 0.40E+00 1.1 0.00000186 -31.2978454163 -1.16E-04
    9 Broy./Diag. 0.40E+00 1.1 0.00000252 -31.2978835492 -3.81E-05
    10 Broy./Diag. 0.40E+00 1.1 5.6405E-09 -31.2978852054 -1.66E-06

    *** SCF run converged in 10 steps ***
    • 如果使用了&SMEARADDED_MOS,分子轨道数目为26个,在计算中总共使用了104个基函数(由笛卡尔高斯分布的原子中心轨道),26在计算范围之内。
  • 最终能量和力的结果,应当检查最终电子密度计算出的电子总数Electronic density on regular grids是否正确。比如此例中为-31.9999999889。

    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
    Electronic density on regular grids:        -31.9999999889        0.0000000111
    Core density on regular grids: 31.9999999939 -0.0000000061
    Total charge density on r-space grids: 0.0000000051
    Total charge density g-space grids: 0.0000000051

    Overlap energy of the core charge distribution: 0.00000000005320
    Self energy of the core charge distribution: -82.06393942512820
    Core Hamiltonian energy: 18.06858429706010
    Hartree energy: 42.41172824581682
    Exchange-correlation energy: -9.71425832315952

    Total energy: -31.29788520535761

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


    ATOMIC FORCES in [a.u.]

    # Atom Kind Element X Y Z
    1 1 Si 0.00000000 0.00000000 0.00000000
    2 1 Si 0.00000000 0.00000001 0.00000001
    3 1 Si 0.00000001 0.00000001 0.00000000
    4 1 Si 0.00000001 0.00000000 0.00000001
    5 1 Si -0.00000001 -0.00000001 -0.00000001
    6 1 Si -0.00000001 -0.00000001 -0.00000001
    7 1 Si -0.00000001 -0.00000001 -0.00000001
    8 1 Si -0.00000001 -0.00000001 -0.00000001
    SUM OF ATOMIC FORCES -0.00000000 -0.00000000 -0.00000000 0.00000000
    • 如果使用了&SMEARADDED_MOS,在Exchange-correlation energy后面会多出

      • 熵entrop(TS)项:为了使计算成为零电子温度结果的可靠近似值,该值应该很小。

        1
        2
        Electronic entropic energy:                                  -0.00001687947145 
        Fermi energy: 0.20867150262130
      • 最终的自由能是DFT总能量与熵能量之和。

        1
        Total energy:                                               -31.29788686349247
      • 外推到TS→0,最终自由能为:(这是最终结果所引用的能量)

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

.wfn

  • 计算得出的最终Kohn-Sham波函数。

  • .wfn.bak-n文件记录了从第n个先前的SCF步骤获得的Kohn-Sham波函数,.wfn.bak-1包含从上一个SCF步骤获得的波函数。如果希望从先前计算的波函数开始新的计算,则需要:

    • FORCE_EVAL部分的SCF子部分中的SCF_GUESS关键字更改为RESTART
    • 新计算与生成波函数的计算的PROJECT_NAME相同,否则需要重命名重新启动wavefunction文件以与新计算的项目名称相对应。
有用可戳(●ˇ∀ˇ●)