As we discussed, there are several challenges to modeling polarons. Simply using DFT+U or hybrid functionals does not guarantee polaron formation. A deliberate approach must be taken, or polaron formation may occur in a random or happenstance manner. A general strategy to modeling polarons can be summarized as: (1) Choose a method to overcome self-interaction errors (e.g. DFT+U or hybrid functionals). (2) Choose initial settings/input that help with convergence to a localized polaron. We discuss several approaches below that help with polaron convergence.

**Bond Distortion Method**

This is one of the simplest, easiest ways to enable charge localization. An analysis of this method and further discussion can be found in our recent paper^{1}. This method has several advantages, as it is easy to implement and does not rely on any particular modeling code. It can be used with many different codes, e.g. VASP, CP2K, Siesta, etc. The basic premise of this method is to create initial geometries that look like polaron geometries, and this helps ensure convergence to localized solutions. Rather than using a symmetrical crystal as the initial structure, in this method the initial geometries have lengthened bonds around the polaron site. Since the structure ‘looks’ like a polaron, convergence to a localized solution more readily occurs. We found that this method is very efficient (more than the electron attractor method) and usually converges to polarons just fine. Scripts and sample input files to create initial geometries with distorted bonds can be found here.

**Electron Attractor Method**

This approach aims to attract electrons to a particular atomic site to enable polaron localization. More positive charge is placed somewhere in the lattice, and since electrons are attracted to positive charge the electrons may localize at that site. The extra positive charge can for instance be obtained by carefully choosing the pseudopotential, which represents the nucleus and core electrons. For instance, Ti has an atomic number of 22, meaning that the nucleus has 22 protons. V has an atomic number of 23 and has a slightly more positive nucleus than Ti. When modeling TiO_{2}, the pseudopotential for one Ti atom can be replaced with a V pseudopotential. When the simulation is run, an excess electron should in principle be attracted to the V site because of its larger positive charge. Once a wavefunction and geometry with localized charge is obtained, then the original pseudopotential can be used (i.e. Ti instead of V). We found that this method can be tricky to use and is slower than the bond distortion method. Sometimes it works, while other times it does not. It also does not always converge to lowest energy states, so use it with care.

**Orbital Occupancy/Control**

If one can choose which orbitals are occupied by electrons, then one can choose to put an extra electron at a particular site and orbital. Again, if the initial wavefunction ‘looks’ like the localized wavefunction, then convergence to the desired localized solution could occur. A method to control orbital occupancy has been implemented in VASP^{2}^{,}^{3}. Other codes may also allow specification or control of the initial orbitals. Constrained DFT^{4}^{,}^{5} adds a constraining potential, and can also be used to control charge/spin.

- Pham, T. D., & Deskins, N. A. (2020). Efficient Method for Modeling Polarons Using Electronic Structure Methods. Journal of Chemical Theory and Computation, 16(8), 5264–5278.
- Allen, J. P., & Watson, G. W. (2014). Occupation matrix control of d- and f-electron localisations using DFT + U. Phys. Chem. Chem. Phys., 16(39), 21016–21031
- https://www.chemistry.tcd.ie/staff/people/gww/gw_new/research/methodology/
- Kaduk, B., Kowalczyk, T., & Van Voorhis, T. (2012). Constrained Density Functional Theory. Chemical Reviews, 112(1), 321–370.
- Ma, H., Wang, W., Kim, S., Cheng, M., Govoni, M., & Galli, G. (2020). PyCDFT: A Python package for constrained density functional theory. Journal of Computational Chemistry, 41(20), 1859–1867.