CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace: • running_qe_computation
howto:running_qe_computation

Table of Contents

How to run calculations like Quantum ESPRESSO

Introduction

Although CP2K uses gaussian as a main function basis, it also has the possibility to do computations with the plane wave basis. This functionality is provided by the SIRIUS library, which supports computation on both CPUs and GPUs. SIRIUS includes most of the core functionalities of QE (only ground state properties), such as non-colinear and colinear magnetism, spin-orbit coupling, on-site Hubbard corrections, forces and stresses.

An other popular program for plane wave DFT is Quantum Expresso (QE). A full description of QE functionalities is outside the scope of this tutorial, but one of the main functionalities is ground state minimization and geometry optimization of crystals. In this tutorial, we will show how to do QE calculations with CP2K + SIRIUS. We will consider the simple case of silicium doped with germanium and show how to convert a QE input file to cp2k input file. Files for this tutorial can be found here while the cp2k input files can be found in the cp2k tests directory.

Word of caution with the pseudo-potential files

Quantum Expresso uses the upf format to describe pseudo-potentials, while CP2K + SIRIUS can either use CP2K pseudo-potentials (they are tailored for the gaussian basis set but will work with the plane wave basis) or the same pseudo-potential files as QE (when the sirius backend is used). The latter case has drawbacks : the upf file should be converted to a json file with a small python script provided by the SIRIUS library, cp2k passing the file name directly to the SIRIUS backend. The SIRIUS backend supports norm-conserving and ultra-soft pseudo-potentials as well as full potential computation, something that is missing in QE.

Problems of units

QE uses atomic units by defaults for the distances but Rydberg unit for energy, while cp2k uses Angstroms by default for the unit of length. The outputs are however in atomic units (or specified otherwise). SIRIUS on the other hand works only in atomic unit. This issue only affect the section describing the crystal structure since Angstroms and Bohr can be mixed in the same input file. By *default* QE will use the Bohr units if nothing is specified while cp2k will use Angstroms.

Using cp2k for doing QE ground state minimization

We assume that the reader is familiar with the QE input format. The cp2k format is described in depth in the documentation but the purpose of this tutorial we will consider a simple example, a cluster of 8 atoms, namely silicium doped with germanium. The QE input file for this system looks like this

QE input CP2K input 
  &control
  calculation='scf',
  restart_mode='from_scratch',
  pseudo_dir = './',
  outdir='./',
  prefix = 'scf_',
  tstress = false,
  tprnfor = false,
  verbosity = 'high',
  wf_collect = false
  /
  &system
  ibrav=0, celldm(1)=1, ecutwfc=25, 
  ecutrho = 400, occupations = 'smearing', 
  smearing = 'gauss', degauss = 0.001,
  nat=8 ntyp=2
  /
  &electrons
  conv_thr =  1.0d-11,
  mixing_beta = 0.7,
  electron_maxstep = 100

  ATOMIC_SPECIES
  Ge 0.0 ge_lda_v1.uspp.F.UPF
  Si 0.0 si_lda_v1.uspp.F.UPF
  CELL_PARAMETERS
   10.26253566   0.00000000   0.00000000
    0.00000000  10.26253566   0.00000000
    0.00000000   0.00000000  10.26253566
  ATOMIC_POSITIONS (crystal)
  Ge  0.00000000  0.00000000  0.00000000
  Si  0.00000000  0.50000000  0.50000000
  Si  0.50000000  0.00000000  0.50000000
  Si  0.50000000  0.50000000  0.00000000
  Si  0.75000000  0.75000000  0.25000000
  Si  0.75000000  0.25000000  0.75000000
  Si  0.25000000  0.75000000  0.75000000
  Si  0.25000000  0.25000000  0.25000000
  K_POINTS (automatic)
  2 2 2  0 0 0
&FORCE_EVAL
  METHOD SIRIUS
  &PW_DFT
    &CONTROL
       PROCESSING_UNIT cpu
       STD_EVP_SOLVER_TYPE lapack
       GEN_EVP_SOLVER_TYPE lapack
    &END CONTROL
    &PARAMETERS
       ELECTRONIC_STRUCTURE_METHOD  pseudopotential
       SMEARING_WIDTH 0.01
       USE_SYMMETRY true
       NUM_MAG_DIMS 0
       GK_CUTOFF 5.0
       PW_CUTOFF 20.00
       ENERGY_TOL 1e-10
       POTENTIAL_TOL 1e-8
       NUM_DFT_ITER 100
       NGRIDK 2 2 2
    &END PARAMETERS
    &ITERATIVE_SOLVER
       ENERGY_TOLERANCE 1e-5
       NUM_STEPS 20
       SUBSPACE_SIZE 4
       TYPE davidson
       CONVERGE_BY_ENERGY 1
    &END ITERATIVE_SOLVER
    &MIXER
       TYPE broyden1
       MAX_HISTORY 8
    &END MIXER
  &END PW_DFT
    &DFT
      &XC
         &XC_FUNCTIONAL
            &LIBXC
               FUNCTIONAL XC_LDA_X
            &END LIBXC
            &LIBXC
               FUNCTIONAL XC_LDA_C_PZ
            &END LIBXC
         &END XC_FUNCTIONAL
      &END XC
    &END DFT
  &SUBSYS
    &CELL
      A [bohr] 10.26253566   0.00000000   0.00000000
      B [bohr]  0.00000000  10.26253566   0.00000000
      C [bohr]  0.00000000   0.00000000  10.26253566
    &END CELL
    &COORD
      SCALED
      Ge  0.00000000  0.00000000  0.00000000
      Si  0.00000000  0.50000000  0.50000000
      Si  0.50000000  0.00000000  0.50000000
      Si  0.50000000  0.50000000  0.00000000
      Si  0.75000000  0.75000000  0.25000000
      Si  0.75000000  0.25000000  0.75000000
      Si  0.25000000  0.75000000  0.75000000
      Si  0.25000000  0.25000000  0.25000000
    &END COORD
    &KIND Si
      POTENTIAL UPF si_lda_v1.uspp.F.UPF.json
    &END KIND
    &KIND Ge
      POTENTIAL UPF ge_lda_v1.4.uspp.F.UPF.json
    &END KIND
  &END SUBSYS
&END FORCE_EVAL
&GLOBAL
  PROJECT Si7Ge
  PRINT_LEVEL MEDIUM
&END GLOBAL

in the following we will explain how to arrive to the above script.

lattice and unit cell description

let us first start with the structure description. In QE, the coordinates of the atoms are often given in reduced coordinates (in lattice coordinates). So let us first start with the lattice vectors. Searching for the keyword 

CELL_PARAMETERS

we find that the lattice vectors (in Cartesian coordinates) are given by (*Warning* the default unit is bohr)

QE input CP2K input 
CELL_PARAMETERS
  10.26253566  0.00000000  0.00000000
  0.00000000  10.26253566  0.00000000
  0.00000000   0.00000000 10.26253566
&FORCE_EVAL
  &SUBSYS
    &CELL
      A [bohr] 10.26253566   0.00000000   0.00000000
      B [bohr]  0.00000000  10.26253566   0.00000000
      C [bohr]  0.00000000   0.00000000  10.26253566
    &END CELL
 &END SUSYS
&END FORCE_EVAL

Note the presence of the keyword [Bohr] which indicates that the lattice vectors are atomic unit. If the QE input file is in angstroms then change the keyword [Bohr] to [angstroms]. CP2K/SIRIUS will convert everything in atomic unit internally.

The atoms positions are given by the section ATOMIC_POSITIONS (crystal) in the QE input file. For this example, the coordinates are given relative to the lattice displacement.

QE input CP2K input 
ATOMIC_POSITIONS (crystal)
    Ge 0.00000000 0.00000000 0.00000000
    Si 0.00000000 0.50000000 0.50000000
    Si 0.50000000 0.00000000 0.50000000
    Si 0.50000000 0.50000000 0.00000000
    Si 0.75000000 0.75000000 0.25000000
    Si 0.75000000 0.25000000 0.75000000
    Si 0.25000000 0.75000000 0.75000000
    Si 0.25000000 0.25000000 0.25000000
&FORCE_EVAL
  &SUBSYS
    &COORD
      SCALED
    Ge 0.00000000 0.00000000 0.00000000
    Si 0.00000000 0.50000000 0.50000000
    Si 0.50000000 0.00000000 0.50000000
    Si 0.50000000 0.50000000 0.00000000
    Si 0.75000000 0.75000000 0.25000000
    Si 0.75000000 0.25000000 0.75000000
    Si 0.25000000 0.75000000 0.75000000
    Si 0.25000000 0.25000000 0.25000000
    &END COORD
 &END SUBSYS
&END FORCE_EVAL

for cp2k input file. Note the presence of the keyword SCALED in the section COORD which indicates that cp2k should treat the coordinates as given in the lattice basis. The coordinates do not have to be given in the lattice basis, any format supported by cp2k will work. Putting these two sections together, we have 

&FORCE_EVAL
  &SUBSYS
    &CELL
      A [bohr] 10.26253566        0.00000000        0.00000000
      B [bohr] 0.00000000       10.26253566        0.00000000
      C [bohr] 0.00000000        0.00000000       10.26253566
    &END CELL
    &COORD
      SCALED
      Ge          0.00000000         0.00000000         0.00000000
      Si          0.00000000         0.50000000         0.50000000
      Si          0.50000000         0.00000000         0.50000000
      Si          0.50000000         0.50000000         0.00000000
      Si          0.75000000         0.75000000         0.25000000
      Si          0.75000000         0.25000000         0.75000000
      Si          0.25000000         0.75000000         0.75000000
      Si          0.25000000         0.25000000         0.25000000
    &END COORD
  &END SUSYS
&END FORCE_EVAL

The information about the potential of each species can be found searching for the keyword ATOMIC_SPECIES in the QE input file.

QE input CP2K input 
ATOMIC_SPECIES
  Ge 0.0 ge_lda_v1.4.uspp.F.UPF
  Si 0.0 si_lda_v1.uspp.F.UPF
&FORCE_EVAL
  &SUBSYS
    KIND Si
       POTENTIAL upf si_lda_v1.uspp.F.UPF.json
    &KIND Si
    KIND Ge
       POTENTIAL upf ge_lda_v1.4.uspp.F.UPF.json
    &KIND
  &END SUSYS
&END FORCE_EVAL

Note the difference of naming for the potential files and the absence of 0.0 this cp2k input file (this parameter is the mass of the atom, it usually ignored in simple QE calculations). CP2K when using sirius as a back-end expects to find the file describing the atomic specie potential in json format. The SIRIUS library distribution contains a python script that can convert one format to the other.

To convert the upf file into json, we can use the script provided here upf_to_json.py and use them as follow 

upf_to_json.py ge_lda_v1.4.uspp.F.UPF
upf_to_json.py si_lda_v1.uspp.F.UPF

Regrouping everything together we obtain 

&FORCE_EVAL
  &SUBSYS
    &CELL
      A [bohr] 10.26253566        0.00000000        0.00000000
      B [bohr] 0.00000000       10.26253566        0.00000000
      C [bohr] 0.00000000        0.00000000       10.26253566
    &END CELL
    &COORD
      SCALED
      Ge          0.00000000         0.00000000         0.00000000
      Si          0.00000000         0.50000000         0.50000000
      Si          0.50000000         0.00000000         0.50000000
      Si          0.50000000         0.50000000         0.00000000
      Si          0.75000000         0.75000000         0.25000000
      Si          0.75000000         0.25000000         0.75000000
      Si          0.25000000         0.75000000         0.75000000
      Si          0.25000000         0.25000000         0.25000000
    &END COORD
        KIND Si
       POTENTIAL upf si_lda_v1.uspp.F.UPF.json
    &KIND Si
    KIND Ge
       POTENTIAL upf ge_lda_v1.4.uspp.F.UPF.json
    &KIND
  &END SUSYS
&END FORCE_EVAL

which is enough to describe the structure we want to study.

According to the QE documentation the default unit of length is the atomic unit (the Bohr) while cp2k uses angstroms as a default.

Energy cut-off and k-point grid

The next step is to convert the QE sections controlling the algorithms used to solve the KS equation. The three most important parameters are ecutwfc the plane wave cut-off for the wave functions, ecutrho the plane wave cut-off for the density, and K_POINTS which gives the dimension of the k-point grid. Let us start with the k point grid. In cp2k, the k-point grid is given in the section PW_DFT (section dedicated to SIRIUS back-end options) using the keyword ngridk.

QE input CP2K input 
  K_POINTS (automatic)
    2 2 2  0 0 0
&FORCE_EVAL
  METHOD SIRIUS
  &PW_DFT
    &PARAMETERS
       NGRIDK 2 2 2
       SHIFTK 0 0 0
    &END PARAMETERS
  &END PW_DFT
&END FORCE_EVAL

Note the presence of 0 0 0 after the k-grid size. Any value different from zero indicate a shift of the generated k-points compared to the gamma point. In some cases, shifting the k-point reduces the size of effective number of k-points to be calculated.

The section FORCE_EVAL now contains the keyword METHOD telling CP2K to use the sirius back-end. It also contains a subsection named *PW_DFT* containing all information necessary for the back-end to do computation. So most of the QE parameters contained in the section control, system and electrons will go in this section.

the parameter ecutwfc which fix the number of planes waves for computation of the wave functions is equivalent to the parameter 

GK_CUTOFF

in cp2k while ecutrho is the plane wave cutoff PW_CUTOFF for the density.

Both parameters in CP2K are the square root of the QE parameters.

QE input CP2K input 
ecutwfc=25
ecutrho = 400
&FORCE_EVAL
  METHOD SIRIUS
  &PW_DFT
    &PARAMETERS
       ELECTRONIC_STRUCTURE_METHOD  pseudopotential
       GK_CUTOFF 5.0
       PW_CUTOFF 20.00
    &END PARAMETERS
  &END PW_DFT
&END FORCE_EVAL

Regrouping again everything together, we obtain the following initial input file for cp2k. 

&FORCE_EVAL
  METHOD SIRIUS
  &PW_DFT
    &PARAMETERS
       ELECTRONIC_STRUCTURE_METHOD  pseudopotential
       GK_CUTOFF 5.0
       PW_CUTOFF 20.00
       NGRIDK 2 2 2
       SHIFTK 0 0 0
    &END PARAMETERS
  &END PW_DFT
  &SUBSYS
    &DFT
      &XC
         &XC_FUNCTIONAL
            &LIBXC
               FUNCTIONAL XC_LDA_X
            &END LIBXC
            &LIBXC
               FUNCTIONAL XC_LDA_C_PZ
            &END LIBXC
         &END XC_FUNCTIONAL
      &END XC
    &END DFT

    &CELL
      A [bohr] 10.26253566        0.00000000        0.00000000
      B [bohr] 0.00000000       10.26253566        0.00000000
      C [bohr] 0.00000000        0.00000000       10.26253566
    &END CELL
    &COORD
      SCALED
      Ge          0.00000000         0.00000000         0.00000000
      Si          0.00000000         0.50000000         0.50000000
      Si          0.50000000         0.00000000         0.50000000
      Si          0.50000000         0.50000000         0.00000000
      Si          0.75000000         0.75000000         0.25000000
      Si          0.75000000         0.25000000         0.75000000
      Si          0.25000000         0.75000000         0.75000000
      Si          0.25000000         0.25000000         0.25000000
    &END COORD
       KIND Si
       POTENTIAL upf si_lda_v1.uspp.F.UPF.json
    &KIND Si
    KIND Ge
       POTENTIAL upf ge_lda_v1.4.uspp.F.UPF.json
    &KIND
  &END SUSYS
&END FORCE_EVAL

&GLOBAL
  PROJECT Si7Ge
  PRINT_LEVEL MEDIUM
&END GLOBAL

The last section is specific to CP2K. The other default settings of the SIRIUS back-end should be enough to converge this simple structure.

By default, it will compute the ground state and return the total energy and the forces on each atom to CP2K. The input file contains a DFT section (which is absent by default in QE) that describes the type of exchange functional to use for the calculations. CP2K and SIRIUS support the libxc naming convention which might be different from the QE convention.

SIRIUS has more parameters to control its execution. Some of these parameters have no equivalent in QE other do. A full list of options that SIRIUS supports can be found in the cp2k documentation, but we will describe the most important one.

the following table summarizes the most common options of QE that have an equivalent in CP2K+SIRIUS.

QE keyword CP2K (SIRIUS) keyword section desciption
nspins = 1, 2 NUM_MAG_DIMS = 0, 1 PW_DFT/PARAMETERS computation in LDA and SLDA
nocolin = true NUM_MAG_DIMS = 3 PW_DFT/PARAMETERS non colinear magnetism
electron_maxstep NUM_DFT_ITER PW_DFT/PARAMETERS number of scf iterations
mixing_beta beta PW_DFT/MIXER mixing parameter for the broyden mixer
nosym = false USE_SYMMETRY = true PW_DFT/PARAMETERS use discret symmetries to build the k-point grid and the of plane waves
ecutrho PW_CUTOFF PW_DFT/PARAMETERS energy cutoff for the density. It is the square root of the QE parameter
ecutwfc GK_CUTOFF PW_DFT/PARAMETERS energy cutoff for the waavefunctions. It is the square root of the QE parameter

many of the other options can be dropped for the first cp2k input file.

If you run this input file with cp2k you should have the following output 

Energy
--------------------------------------------------------------------------------
valence_eval_sum          :        -5.08815701
<rho|V^{XC}>              :       -30.23739234
<rho|E^{XC}>              :       -43.72437220
<mag|B^{XC}>              :         0.00000000
<rho|V^{H}>               :       106.32392644
one-electron contribution :       -81.17469111 (Ha),      -162.34938222 (Ry)
hartree contribution      :        53.16196322
xc contribution           :       -43.72437220
ewald contribution        :       -68.41309922
PAW contribution          :         0.00000000
Total energy              :      -140.15019932 (Ha),      -280.30039863 (Ry)

band gap (eV) :         0.34564591
Efermi        :         0.24495087

iteration :  11, RMS 1.265254419560E-11, energy difference : 0.000000000000E+00

converged after 12 SCF iterations!

 ENERGY| Total FORCE_EVAL ( SIRIUS ) energy (a.u.):         -140.150199315127225

Molecular dynamics

Quantum expresso supports many types of calculations that cp2k supports natively. For instance, it is possible to do molecular dynamics in cp2k with the sirius backend, something that QE can do as well. To do this, we can start from the previous cp2k input file for Si7Ge and add a small section for setting the parameters relevant for the molecular dynamics. The symmetry parameter should be off during molecular dynamics configurations. 

&FORCE_EVAL
  METHOD SIRIUS
  &PW_DFT
    &CONTROL
       PROCESSING_UNIT cpu
       STD_EVP_SOLVER_TYPE lapack
       GEN_EVP_SOLVER_TYPE lapack
    &END CONTROL
    &PARAMETERS
       ELECTRONIC_STRUCTURE_METHOD  pseudopotential
       SMEARING_WIDTH 0.01
       USE_SYMMETRY false
       NUM_MAG_DIMS 0
       GK_CUTOFF 5.0
       PW_CUTOFF 20.00
       ENERGY_TOL 1e-10
       POTENTIAL_TOL 1e-8
       NUM_DFT_ITER 100
       NGRIDK 2 2 2
    &END PARAMETERS
    &ITERATIVE_SOLVER
       ENERGY_TOLERANCE 1e-5
       NUM_STEPS 20
       SUBSPACE_SIZE 4
       TYPE davidson
       CONVERGE_BY_ENERGY 1
    &END ITERATIVE_SOLVER
    &MIXER
       TYPE broyden1
       MAX_HISTORY 8
    &END MIXER
  &END PW_DFT
    &DFT
      &XC
         &XC_FUNCTIONAL
            &LIBXC
               FUNCTIONAL XC_LDA_X
            &END LIBXC
            &LIBXC
               FUNCTIONAL XC_LDA_C_PZ
            &END LIBXC
         &END XC_FUNCTIONAL
      &END XC
    &END DFT
  &SUBSYS
    &CELL
      A [bohr] 10.26253566        0.00000000        0.00000000
      B [bohr] 0.00000000       10.26253566        0.00000000
      C [bohr] 0.00000000        0.00000000       10.26253566
    &END CELL
    &COORD
      SCALED
      Ge          0.00000000         0.00000000         0.00000000
      Si          0.00000000         0.50000000         0.50000000
      Si          0.50000000         0.00000000         0.50000000
      Si          0.50000000         0.50000000         0.00000000
      Si          0.75000000         0.75000000         0.25000000
      Si          0.75000000         0.25000000         0.75000000
      Si          0.25000000         0.75000000         0.75000000
      Si          0.25000000         0.25000000         0.25000000
    &END COORD
    &KIND Si
      POTENTIAL UPF si_lda_v1.uspp.F.UPF.json
    &END KIND
    &KIND Ge
      POTENTIAL UPF ge_lda_v1.4.uspp.F.UPF.json
    &END KIND
  &END SUBSYS
&END FORCE_EVAL
&MOTION
   &GEO_OPT
      OPTIMIZER LBFGS
   &END GEO_OPT
   &MD
      ENSEMBLE NVE
      TIMESTEP 0.1 
      STEPS 125
      TEMPERATURE 300.0
   &END MD
&END MOTION
&GLOBAL
  PROJECT Si7Ge
  PRINT_LEVEL MEDIUM
  RUN_TYPE MD
&END GLOBAL

if you run this script using the command 

PATH_TO_CP2K/cp2k.psmp Si7Ge-motion.inp

you should obtain a file Si7Ge-pos-1.xyz that contains the positions of the atoms and Si7Ge-1.ener that contains the different energies and temperature as a function of time. The calculations are done in the micro-canonical ensemble which means that the total energy (the sum of the potential which is the total energy given by SIRIUS and the kinetic energy) is conserved during the integration process. It is possible to compare these results with QE, but instead of doing the molecular dynamics with QE, we take the positions of the atoms at each time step and compute the potential energy for this atom position with QE. The total energy given by QE and CP2K/SIRIUS are shifted by a small offset. More generally, converting the input file from QE to CP2K is not a warranty to obtain the same results for the total energy. The reasons for this multiple :

  • Indicating the same cutoff does not warranty the fft grid size.
  • QE and SIRIUS treat the radial integrals interpolation differently
  • The functional in QE do not use libxc by default, while CP2K/SIRIUS does.
  • The minimization method is sensitive to the initial states.

For all these reasons, we should not try to compare results at the binary level.

howto/running_qe_computation.txt · Last modified: 2021/12/08 10:18 by mtaillefumier