How to Converge the CUTOFF and REL_CUTOFF

How to Converge the CUTOFF and REL_CUTOFF

Introduction

QUICKSTEP, as with nearly all ab initio Density Functional Theory simulation packages, requires the use of a real-space (RS) integration grid to represent certain functions, such as the electron density and the product Gaussian functions.

QUICKSTEP uses a multi-grid system for mapping the product Gaussians onto the RS grid(s), so that wide and smooth Gaussian functions are mapped onto a coarser grid than narrow and sharp Gaussians.

The electron density is always mapped onto the finest grid.

Choosing a fine enough integration grid for a calculation is crucial in obtaining meaningful and accurate results.

In this tutorial, we will show the reader how to systematically find the correct settings for obtaining a sufficiently fine integration grid for his/her calculation.

This tutorial assumes the reader already has some knowledge of how to perform a simple energy calculation using QUICKSTEP (this can be found in tutorial: Calculating Energy and Forces using Quickstep).

A completed example from an earlier calculation can be obtained from the file converging_grid.tgz that comes with this tutorial.

The calculations were carried out using CP2K version 2.4.

‘’QUICKSTEP’’ Multi-Grid

Before we go through the input file, it is worthwhile to explain how the multi-grid is constructed in QUICKSTEP, and how the Gaussians are mapped onto the different grid levels.

Hopefully this will offer the reader a clear picture of how the key control parameters affect the grids, and thus the overall accuracy of a calculation.

All multi-grid related settings for a calculation is controlled via keywords in MULTIGRID subsection of DFT subsection in FORCE_EVAL.

The number of levels for the multi-grid is defined by NGRIDS, and by default this is set to 4.

The keyword CUTOFF defines the planewave cutoff (default unit is in Ry) for the finest level of the multi-grid.

The higher the planewave cutoff, the finer the grid.

The corresponding planewave cutoffs for the subsequent grid levels (from finer to coarser) are defined by the formula:

Eicut=E1cutα(i−1)

where α has a default value of 3.0, and since CP2K versions 2.0, can be configured by the keyword PROGRESSION_FACTOR.

Therefore, the higher the value of CUTOFF the finer grid for all multi-grid levels.

Having constructed the multi-grid, QUICKSTEP then needs to map the Gaussians onto the grids.

The keyword REL_CUTOFF controls which product Gaussians are mapped onto which level of the multi-grid.

CP2K tries to map each Gaussian onto a grid such that the number of grid points covered by the Gaussian—no matter how wide or narrow—are roughly the same.

REL_CUTOFF defines the planewave cutoff of a reference grid covered by a Gaussian with unit standard deviation (e|→r|2).

A Gaussian is mapped onto the coarsest level of the multi-grid, on which the function will cover number of grid points greater than or equal to the number of grid points e|→r|2 will cover on a reference grid defined by REL_CUTOFF.

Therefore, the two most important keywords effecting the integration grid and the accuracy of a calculation are CUTOFF and REL_CUTOFF.

If CUTOFF is too low, then all grids will be coarse and the calculation may become inaccurate; and if REL_CUTOFF is too low, then even if you have a high CUTOFF, all Gaussians will be mapped onto the coarsest level of the multi-grid, and thus the effective integration grid for the calculation may still be too coarse.

Example: Bulk Si with 8 atoms in a cubic cell

We demonstrate the process using an example based on Bulk Si with 8 atoms in a face centred cubic unit cell.

Template Input File

To systematically find the best CUTOFF and REL_CUTOFF values which are sufficient for a given accuracy (say, 10−6 Ry in total energy), we need to perform a series of single point energy calculations.

It is much easier to use a set of scripts that can automate this process.

To do this, we first write a template input file: template.inp, as shown below:

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
&GLOBAL
PROJECT Si_bulk8
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
NGRIDS 4
CUTOFF LT_cutoff
REL_CUTOFF LT_rel_cutoff
&END MGRID
&QS
EPS_DEFAULT 1.0E-10
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-6
MAX_SCF 1
ADDED_MOS 10
CHOLESKY INVERSE
&SMEAR ON
METHOD FERMI_DIRAC
ELECTRONIC_TEMPERATURE [K] 300
&END SMEAR
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION
&MIXING
METHOD BROYDEN_MIXING
ALPHA 0.4
BETA 0.5
NBROYDEN 8
&END MIXING
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&KIND Si
ELEMENT Si
BASIS_SET SZV-GTH-PADE
POTENTIAL GTH-PADE-q4
&END KIND
&CELL
SYMMETRY CUBIC
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
&PRINT
&TOTAL_NUMBERS ON
&END TOTAL_NUMBERS
&END PRINT
&END FORCE_EVAL