In this exercise, we are going to learn the basics of using CP2K by calculating the total energy, optimizing the geometry, and printing the projected density of states (PDOS) of a single H$_2$O molecule. We will also investigate how different parameters such as the DFT functional and basis sets affect the accuracy and duration of our calculations.

In the first part of the exercise, we will go through the input and steps needed to perform a single-point calculation of a water molecule. The first thing we need is to create an input file describing our system and choosing the methods we want to use. Below is a simple example input that we need.

water.inp

&GLOBAL
  PROJECT   h2o
  RUN_TYPE  ENERGY
  IOLEVEL   LOW
&END GLOBAL
&FORCE_EVAL
  METHOD QS
  &DFT
    BASIS_SET_FILE_NAME GTH_BASIS_SETS
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    LSD 0
    &QS      
      EPS_DEFAULT 1.0E-10
    &END QS
    &SCF
      MAX_SCF    300
      EPS_SCF    1.0E-06
      SCF_GUESS  ATOMIC
      &MIXING
        METHOD  DIRECT_P_MIXING
        ALPHA   0.6
      &END MIXING
      &DIAGONALIZATION
	    ALGORITHM STANDARD
      &END DIAGONALIZATION       
    &END SCF
    &MGRID
      NGRIDS        4
      CUTOFF        300
      REL_CUTOFF    60
    &END
    &XC
      &XC_FUNCTIONAL BLYP
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC    15.0  15.0  15.0
      ANGLES 90.00 90.00 90.00
    &END CELL
    &COORD
      O  0 0 0
      H  0.7 0.7 0
      H -0.7 0.7 0
    &END COORD
    &KIND H                 
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-BLYP-q1
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH
      POTENTIAL GTH-BLYP-q6
    &END KIND
    &PRINT
      &ATOMIC_COORDINATES
      &END
    &END PRINT 
  &END SUBSYS
&END FORCE_EVAL

The input file is structured in different sections, started by the “&SECTION_NAME” and ended by “&END SECTION NAME” where the including the section name after &END is optional but recommended. Each section may contain additional sections, and keywords used to control the parameters. There are three main input sections that we will use in this exercise:

1. “&GLOBAL” - contains the general options such as the kind of run (in this case ENERGY) and the name of the project.
2. “&FORCE_EVAL” - sets the parameters needed to evaluate the forces on the atoms, including the atomic positions, DFT functional, and basis set.
3. “&MOTION” - controls e.g. geometry optimizations and MD.

The &GLOBAL section is quite straightforward and controls the type of calculation, our project’s name, etc. but the &FORCE_EVAL section is more complicated and we will go through the most important parts more carefully. The keyword

METHOD Quickstep

chooses that we want to use the Quickstep code which means DFT with the Gaussian and Planewaves method (GPW). We then move in two the &DFT section where we have to specify in which files the program should search for the basis sets and potentials.

BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME GTH_POTENTIALS

In the &QS section, the EPS_DEFAULT keyword sets all tolerances within the Quickstep method. Individual tolerances can also be set, overwriting the EPS_DEFAULT value in the process. The &SCF section controls the self-consistent optimization method for the Kohn-Sham orbitals, with important keywords such as the maximum number of self-consistent loops, MAX_SCF, the initial guess of the wavefunction SCF_GUESS, and the convergence tolerance EPS_SCF. The &MIXING subsection

&MIXING
  METHOD  DIRECT_P_MIXING
  ALPHA   0.6
&END MIXING

sets how much of the new density that should be mixed with the old density in each step. In this case, 60 % of the new density will be combined with 40 % of the old density when moving on to the next step. A small mixing parameter ALPHA will slow down the calculation but can help with convergence in hard cases. Finally the &MGRID section

&MGRID
  NGRIDS        4
  CUTOFF        300
  REL_CUTOFF    60
&END

controls how to set up the integration grid used in the Quickstep method. NGRIDS chooses the number of multigrids to use, CUTOFF selects the number of plane-wave basis functions to describe the electron density, and the REL_CUTOFF determines the grid spacing of the Gaussian basis functions.

After the &DFT section, the other main part of &FORCE_EVAL is the &SUBSYS. Here we set up the system we would like to study including basis sets and potentials. The &CELL and the &COORD sections

&CELL
  ABC 15.0 15.0 15.0
  ANGLES 90.00 90.00 90.00
&END CELL
&COORD
 O  0.0 0.0 0.0
 H  0.7 0.7 0.0
 H -0.7 0.7 0.0
&END COORD 

sets up our simulation cell and atomic coordinates respectively. In the &KIND sections,

&KIND H                 
  BASIS_SET DZVP-GTH
  POTENTIAL GTH-BLYP-q1
&END KIND
&KIND O
  BASIS_SET DZVP-GTH
  POTENTIAL GTH-BLYP-q6
&END KIND

the basis sets and potentials associated with each atom is given. In this case, we use GTH pseudopotentials and corresponding basis sets optimized for them.

Remember that a lot of information about the different sections and the keywords can be found in the CP2K manual https://manual.cp2k.org/.

We are now ready to run the input in CP2K. To start the calculation, type

  ./cp2k.sopt -i water.inp -o water.out & 

or

  ./cp2k.sopt water.inp | tee water.out

if you want to follow the output as it runs. When the program is finished, you will have two new files in your working directory:

  water.out   
  h2o-RESTART.wfn

The first file ending with .out is the main output file of the calculation which contains information about the optimization