CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== 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. [[https://github.com/electronic-structure/SIRIUS|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 [[http://www.quantum-espresso.org|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 [[https://www.cp2k.org/_media/howto:si7ge.tar.gz|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 ^
|<code>
  &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
</code>|<code>
&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
</code>|
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 <code>CELL_PARAMETERS</code> we find
that the lattice vectors (in Cartesian coordinates) are given by (*Warning* the default unit is bohr)
^ QE input ^ CP2K input ^
|<code>
CELL_PARAMETERS
  10.26253566  0.00000000  0.00000000
  0.00000000  10.26253566  0.00000000
  0.00000000   0.00000000 10.26253566
</code>|<code>
&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
</code>|

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 ^ 
|<code>
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
</code>|<code>
&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
</code>|

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

<code>
&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
</code>

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 ^
|<code>
ATOMIC_SPECIES
  Ge 0.0 ge_lda_v1.4.uspp.F.UPF
  Si 0.0 si_lda_v1.uspp.F.UPF
</code>|<code>
&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
</code>|

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 [[https://github.com/electronic-structure/SIRIUS/blob/develop/apps/upf/upf_to_json.py|upf_to_json.py]] and use them as follow

<code>
upf_to_json.py ge_lda_v1.4.uspp.F.UPF
upf_to_json.py si_lda_v1.uspp.F.UPF
</code>

Regrouping everything together we obtain

<code>
&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
</code>

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 ^
|<code>
  K_POINTS (automatic)
    2 2 2  0 0 0
</code>|<code>
&FORCE_EVAL
  METHOD SIRIUS
  &PW_DFT
    &PARAMETERS
       NGRIDK 2 2 2
       SHIFTK 0 0 0
    &END PARAMETERS
  &END PW_DFT
&END FORCE_EVAL
</code>|
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 <code>GK_CUTOFF</code> 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 ^
|<code>
ecutwfc=25
ecutrho = 400
</code>|<code>
&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
</code>|

Regrouping again everything together, we obtain the following initial input file for cp2k.
<code>
&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

</code>

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 [[https://tddft.org/programs/libxc/|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

<code>
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
</code>

===== 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. 
 
<code>
&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
</code>
if you run this script using the command
<code>
PATH_TO_CP2K/cp2k.psmp Si7Ge-motion.inp
</code>

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. 
{{:howto:cp2k-sirius-md.png|}}
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.