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.
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.
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.
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.
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.
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
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 :
For all these reasons, we should not try to compare results at the binary level.