====== 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 ^ | &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 [[https://github.com/electronic-structure/SIRIUS/blob/develop/apps/upf/upf_to_json.py|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 [[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 Energy -------------------------------------------------------------------------------- valence_eval_sum : -5.08815701 : -30.23739234 : -43.72437220 : 0.00000000 : 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. {{: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.