Most parts of this exercise are adapted from the cubic Si example on the CP2K “How to” website: How to Calculate Energy and Forces

As a first step, we perform a static self-consistent Kohn–Sham density functional theory (DFT) calculation to obtain the electronic energy. Our model system is the L-alanine crystal, consisting of 13 atoms per molecule and 4 molecules in the unit cell.

We first look at the files required for this calculation, that can be downloaded here:

• Lalanine.inp: the main input file for the calculation, which defines the system and the job parameters
• Lalanine.cif: the coordinate file when defined externally

Parameter files:

• BASIS_SET: file containing parameters for the basis sets
• GTH_POTENTIALS: file containing parameters for the pseudopotentials

A list of basis set and pseudopotential files may be found in cp2k/data/, which comes with a CP2K source release.

Let us look at the main input: Lalanine.inp in detail. The filename and extension can be chosen freely. The input file is structured into ordered blocks and keywords, the order of which are unimportant. Each input block is referred to as a 'section' in this tutorial, and some sections are 'subsections' of other sections. Full definitions of the input file format and keywords are available via the CP2K input reference manual.

The main sections in the input file are:

• GLOBAL: contains general options for the CP2K run, such as the name of the job, the type of run etc.
• FORCE_EVAL: contains all parameters associated with the evaluation of forces on atoms, including the initial atomic coordinates.

We now look at each section in detail, starting with the GLOBAL section from Lalanine.inp:

&GLOBAL
  PROJECT LALANINE
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL

Here we perform a static energy calculation; in this case, we must set RUN_TYPE to ENERGY. The keyword PROJECT sets the root name of the calculation, in this case Lalanine. Any output files automatically generated by CP2K will have the name prefixed by Lalanine. PRINT_LEVEL controls the default verbosity of the main output of CP2K, in this example it is set to 'medium'.

Next, we examine the FORCE_EVAL section line-by-line.

METHOD Quickstep

The keyword METHOD specifies QUICKSTEP as the approach for evaluating the forces on the atoms, i.e. density functional theory using the Gaussian and Plane Waves (GPW) method.

The subsection SUBSYS defines the simulation unit cell and the initial coordinates of the atoms in the calculation.

&SUBSYS
  &CELL
    ABC 5.94 12.274 5.806
    ALPHA_BETA_GAMMA 90.0 90.0 90.0
  &END CELL
  &COORD
    H 2.524500 3.041497 1.962428
    ...
  &END COORD
  &KIND H
    ELEMENT H
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q1
  &END KIND
  &KIND C
    ELEMENT C
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q4
  &END KIND
  &KIND N
    ELEMENT N
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q5
  &END KIND
  &KIND O
    ELEMENT O
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q6
  &END KIND
&END SUBSYS

The subsection CELL defines the simulation cell used in a calculation. In this example, we have a tetragonal cell with $a$=5.94 Å, $b$=12.274 Å, $c$=5.806 Å.

The initial atomic coordinates are specified in the COORD subsection. The default input format for atomic coordinates in CP2K is: <ATOM_KIND> X Y Z where $X$, $Y$, $Z$ are Cartesian coordinates in Å. This can be changed by configuring keyword SCALED to .TRUE. in the COORDsubsection, which makes the coordinate input $X$ $Y$ $Z$ fractional with respect to the lattice vectors. One can also change the unit for the Cartesian coordinates by setting the keyword UNIT within the subsection. <ATOM_KIND> should be a label that corresponds to the definition of the elements in the KIND subsections.

The subsection KIND gives definitions of elements in the calculation. There must be one KIND subsection per element. In this example, we have defined the basis set to be DZVP-GTH-PBE (double-ζ with polarisation basis optimised for Goedecker–Teter–Hutter PBE pseudopotentials); and the pseudopotentials such as GTH-PBE-q4 (Goedecker–Teter–Hutter PBE pseudopotential with 4 valence electrons).

The basis set and pseudopotential names must correspond to an existing entry in the corresponding basis set and pseudopotential files defined by the BASIS_SET_FILE_NAME and POTENTIAL_FILE_NAME keywords in the DFT subsection of the FORCE_EVAL section.

After the SUBSYS section in Lalanine.inp comes the DFT subsection, which controls all aspects of the self-consistent Kohn–Sham density functional theory calculation. This subsection is relevant only if the METHOD keyword in FORCE_EVAL is set to QUICKSTEP.

BASIS_SET_FILE_NAME  BASIS_SET
POTENTIAL_FILE_NAME  GTH_POTENTIALS

As already mentioned above, the keywords BASIS_SET_FILE_NAME and POTENTIAL_FILE_NAME set the files that contain basis set and pseudopotential parameters.

&QS
  EPS_DEFAULT 1.0E-10
&END QS

General control parameters for QUICKSTEP are specified in the QS subsection. The parameter EPS_DEFAULT defines the global default tolerance applied throughout the calculation. More specific tolerance parameters (EPS_*), if present, override this global value.

&MGRID
  NGRIDS 4
  CUTOFF 300
  REL_CUTOFF 60
&END MGRID

The MGRID subsection defines how the integration grid used in QUICKSTEP calculations should be set up. QUICKSTEP employs a multigrid scheme to represent Gaussian functions numerically: Narrow and sharp Gaussians are mapped onto finer grids than wider and smoother Gaussians. In this case, we are telling the code to set up 4 levels of multigrids, with the planewave cutoff of the finest grid set to 300 Ry, and with the grid spacing underneath any Gaussian functions to be finer than the equivalent planewave cutoff of 60 Ry. The users should read the tutorial How to Converge the CUTOFF and REL_CUTOFF for details on how these parameters affect the grid constructed, and how to define a sufficient grid for their calculation.

&XC
  &XC_FUNCTIONAL PBE
  &END XC_FUNCTIONAL
&END XC

This subsection defines which exchange–correlation density functional we want to use. In this case we choose the PBE functional, which is consistent with the basis set and pseudopotential we have chosen.

&SCF
  SCF_GUESS ATOMIC
  EPS_SCF 1.0E-7
  MAX_SCF 300
  &DIAGONALIZATION
    ALGORITHM STANDARD
  &END DIAGONALIZATION
  &MIXING
    METHOD BROYDEN_MIXING
    ALPHA 0.4
    NBROYDEN 8
  &END MIXING
&END SCF

The SCF subsection defines all the settings related to methods used to find a self-consistent solution of the Kohn–Sham DFT formalism.

SCF_GUESS specifies how the initial trial electron density function $\vec{r}$ is generated. In this example (ATOMIC), the initial density is obtained from overlapping atomic charge densities. A good starting point for the electron density in the self-consistency loop is important for achieving convergence efficiently. EPS_SCF sets the tolerance of the charge density residual. This overrides the EPS_DEFAULT value set in the QS subsection. MAX_SCF sets the maximum number of self-consistency loops QUICKSTEP is allowed to perform for each ground-state energy calculation.

&DIAGONALIZATION ON
  ALGORITHM STANDARD
&END DIAGONALIZATION

The DIAGONALIZATION subsection instructs CP2K to use the traditional diagonalisation method for obtaining the ground-state Kohn–Sham energy and electron density. The subsection heading also takes an argument; in this case it is set to ON, which is equivalent to .TRUE. or T, and enables diagonalisation. If the argument is omitted, the default is also .TRUE.. The alternative to diagonalisation is the Orbital Transform (OT) method. In that case, the user should either delete the DIAGONALIZATION block or set the heading to OFF (or .FALSE.), and instead add the OT subsection. The keyword ALGORITHM specifies the diagonalisation algorithm; here STANDARD means the standard LAPACK/SCALAPACK routines are used.

&MIXING T
  METHOD BROYDEN_MIXING
  ALPHA 0.4
  NBROYDEN 8
&END MIXING

The MIXING subsection defines the parameters for charge mixing in a self-consistency calculation. It can take a value of either .TRUE. (T) or .FALSE. (F), which switches charge mixing on or off. The default is .TRUE.. Note that this subsection only applies when using the diagonalisation method; the OT method uses a different approach to charge mixing, which is explained in other tutorials. The keyword ALPHA sets the mixing parameter: in this example, 40% of the output density is mixed with 60% of the input density to form the new input density in the next SCF iteration. The keyword METHOD specifies the mixing method; here Broyden mixing is used. The keyword NBROYDEN (an alias for NBUFFER) sets the number of previous iterations to be stored in the Broyden mixing algorithm.

To run the calculation, the user needs to have a working CP2K package compiled in the system PATH. The calculation can be started with the command:

cp2k.psmp -o Lalanine.out Lalanine.inp &

This launches CP2K in the background. The option -o redirects the CP2K output to the file Lalanine.out. Note that the -o option appends output from successive runs to the same file. If you want to start from scratch, you should delete Lalanine.out before launching a new calculation.

After the job has finished, the following files should be present:

1. Lalanine.out
2. Lalanine-RESTART.wfn

The file Lalanine.out contains the main output of the job. The file Lalanine-RESTART.wfn contains the final Kohn–Sham wavefunctions from the calculation. Files of the form Lalanine-RESTART.wfn.bak-<n> store the wavefunctions obtained from the $n$-th previous SCF step; in this case, Lalanine-RESTART.wfn.bak-1 contains the wavefunctions from the last SCF step.

If you want to start a new calculation from previously calculated wavefunctions, change the keyword SCF_GUESS in the SCF subsection of the FORCE_EVAL section to:

SCF_GUESS RESTART

This requires that the new calculation uses the same PROJECT_NAME as the one that generated the restart files. Otherwise, the restart wavefunction files need to be renamed accordingly.

In the output file, CP2K reports quantities such as the number of electrons, the number of occupied orbitals, the SCF convergence behaviour, the total energy, and the atomic forces. During the class we will examine these sections of the output in detail.