CP2K Polaron Geometry Optimisation with PBE+U¶
Author: Chris Ahart
In this tutorial I will show how to perform calculations of polaron structures in CP2K. Hopefully by the end of this tutorial you will be comfortable generating images such as the one above. The system chosen is the hole polaron in anatase TiO2, where DFT+U calculations are capable of reproducing polaron structures in qualitative agreement with hybrid calculations. To make this tutorial suitable for those without tier 1 HPC access, I will not include any hybrid calculations and will only consider small supercells. If you are interested to use hybrid functionals in CP2K, then please speak with me.
This tutorial is not intended as a beginners guide to DFT or CP2K. For those new to CP2K, I recommend watching the recorded lectures and working through the exercises from the CP2K Workshop 2019 at Ghent University. I have included some additional links at the at end of tutorial, as well as to some polaron papers that you can read if you are interested to learn more about polarons.
All files corresponding to the tutorial are available on our group Github page. I have organised the files such that the examples you can run are in the folder /examples, and my original files are also included for reference in a separate folder.
To begin, we perform PBE geometry optimisation of anatase TiO2. Refer to your HPC provider of choice for a job submission script, and run CP2K using the provided input file. CP2K will generate a number of output files, including a log file which contains all relevant information for the calculation.
For a CP2K log file 'cp2k_log.log' run command
For a geometry optimisation including 6 steps, there are 8 total energies. The last two energies should be equal within numerical error, as CP2K will re-evaluate the energy after completing the geometry optimisation. Because this is a optimisation for neutral geometry starting from a good starting geometry, the energy smoothly decreases to the minimum.
Run command
This is the HOMO LUMO gap, 2.07 eV compared to the experimental band gap 3.2 eV. This underestimation of the band gap is expected for PBE, and generally hybrid functionals are required to reproduce experimental band gaps.
Next we perform a geometry optimisation starting using the automatically generated restart file 'tio2-1.restart'. You should make a new folder and re-name 'tio2-1.restart' to an input filename of your choice such as 'input.inp'. Change lines 'CHARGE 0' to 'CHARGE 1' and 'MULTIPLICITY 1' to 'MULTIPLICITY 2' in 'input.inp' to create a charge of 1, equivalent to removing an electron from the system. The multiplicity keyword can just be removed and CP2K will automatically calculate it, but it is worthwhile to get in the habit of controlling the multiplicity yourself as it must be manually set if studying magnetic systems. Note that CP2K will automatically add a homogeneous background charge in order to keep the system neutral, as otherwise the electrostatic energy will diverge. This is standard in all DFT codes, but is worth understanding that performing calculations on charged systems introduces a number of challenges and approximations.
Run this new calculation, and run same commands
We find that that the geometry optimisation does not result in polaron formation. This can be seen by the negligible change in energy upon geometry optimisation, and by the HOMO-LUMO gap. CP2K will always remove electrons from the beta spin channel and add electrons to the alpha spin channel, as such the alpha spin channel HOMO-LUMO gap is unchanged as the number of alpha electrons is the same while we have removed an electron from the valence band of the beta spin channel. The HOMO-LUMO gap of 0.0 eV for the beta spin channel means that the electron hole is fully delocalised across the valence band. When a polaron forms, there will be a non-zero HOMO-LUMO gap that is referred to as a 'mid-gap state'.
Now, let us include a Hubbard U term. This is performed by adding a new section to the input file to the section that describes the electronic structure of oxygen atoms
CP2K automatically will default to units of Hartree for U_MINUS_J, so it is important to add units [eV] to tell CP2K to use eV. Running the file should result in a longer geometry optimisation, which for my calculation is 28 steps including the 6 steps from neutral geometry optimisation as we have restarted using the .restart file without resetting the number of geometry optimisation steps (keyword 'STEP_START_VAL').
Running the same commands:
We can calculate the difference in the first and last energy -5772.858935726821983 - -5772.898198909129860 = 0.039263182307877 and convert to eV = 1.07 eV. This is the trapping energy, also known as the polaron self-trapping energy, polaron formation energy or the polaron reorganisation energy. In addition, we see that there is a HOMO-LUMO gap of 2.29 eV in the beta spin channel. This means that a polaron mid-gap state has formed close to the conduction band. You can also look at the Hirshfeld charge analysis, which for my calculation is found at hirshfeld/hirshfeld-1_28.hirshfeld. '28' refers to the index of the geometry optimisation, so find your equivalent file. In this file I find line
To finish, you can plot the CP2K output .cube files using VESTA or VMD. The files 'SPIN_DENSITY' contain the spin density (difference in electron density for the alpha and beta spin channels) printed at each geometry optimisation step. You should see that the spin density is initially delocalised over all atoms, and as the geometry optimisation progresses the spin density localises onto a single oxygen atom. The final spin density should be consistent with the image presented at the beginning of this tutorial.
This tutorial continues in part 2, 'CP2K Polaron DFT-MD with PBE+U'.
Extensions:¶
- To extend this work, you could look at the change in trapping energy and band gap as a function of the Hubbard U. You will find that the trapping energy increases linearly with U, and that no value of oxygen Hubbard U can reproduce the band gap. To reproduce the band gap, you can add a Hubbard U to the Ti atoms. I find that a value of U = 6 eV for Ti and a value of U = 2 for O reproduces the experimental band gap.
- You could also look at the change in trapping energy with supercell size. Ideally, you should increase the supercell size until the trapping energy is converged. In practice, this may not always be possible and some compromise must be found between accuracy and computational cost.
- You can also run hybrid DFT calculations, which will allow you to get a much more accurate polaron trapping energy and related electron transfer parameters. This will require a few hundred CPUs, so I do not include such calculations here.
Regarding accuracy:¶
From my understanding of plane wave DFT calculations, the total energy can be converged with respect to the plane wave cutoff in a well behaved manner. This is not the case for Gaussian type orbital codes such as CP2K. You are never going to perform a CP2K calculation that is converged with respect to total energy. Instead, focus on energy differences relevant to properties that you are calculating. Generally, most structural and electronic properties are well described at double-zeta basis set level and it is rare to find condensed phase calculations that are at triple-zeta level or above. There is however some evidence that the widespread use of double-zeta basis sets is problematic. If you are interested, refer to the thesis of Tiziano Müller where he presents the convergence of total energy with respect to basis set size in CP2K versus the all-electron plane wave code Wien2k. He ultimately recommends use of triple-zeta basis sets.
Other remarks:¶
- CP2K is a challenging code to use, with many keywords and generally poor documentation. It can take a while to become experienced with CP2K, and get a 'feel' for appropriate keyword choices. The values in this work are a reasonable starting point, and should work for most systems that are not magnetic or strongly correlated.
- Polarons are generally only indirectly measured by experiments. For example in conductivity measurements polarons are characterised by an increase in conductivity with temperature, as there is an activation barrier for polaron hopping. This is different to band transport, characterised by a decrease in conductivity with temperature as a result of phonon scattering. As such, polaron calculations can be seen as problematic becasue the structure and dynamics of polarons are only indirectly measurable by experiment. There is often considerable disagreement regarding the structure and dynamics of polarons in different materials.
- Optimally tuned range-separated hybrid functionals such as HSE or PBE0-TC-LRC remain the only appropriate functionals for polaron structure and dynamics. As these are expensive, I am currently exploring how machine learning can be used to increase the accessible length-scales and time-scales while preserving hybrid level accuracy. I hope that in the future advances in density functional theory will provide more accurate and cheaper functionals, however this seems unlikely to happen in the near-future.
- There are a number of physics papers that present analytical modelling of polarons, typically based on model Hamiltonians (e.g. Luo et al, Lafuente-Bartolome et al and Sio et al ). I do not trust these approaches, as they neglect the most important factor in small polaron transport: the electronic coupling of the initial and final diabatic states.
CP2K website resources:¶
Other CP2K resources:¶
Other tutorials I have written:¶
- Potential control and current induced forces using CP2K+SMEAGOL, Imperial College London
- Converging magnetic systems in CP2K, Imperial College London
Other tutorials I have contributed to:¶
Polaron reference material:¶
- Small Polarons in Transition Metal Oxides, Handbook of Materials Modeling.
- Formation and dynamics of small polarons on the rutile TiO2(110) surface, Article.
- First-Principles Modeling of Polaron Formation in TiO2 Polymorphs, Article.
- Ferrotoroidic polarons in antiferrodistortive SrTiO3
- Machine Learning Small Polaron Dynamics
- Polarons on transition-metal oxide surfaces, Thesis.
- Charge transport in bulk hematite and at the hematite/water interface, Thesis.