THEOS

The oxidation state of an atom is an empirical concept widely used in chemistry and in materials science - see e.g. Wikipedia for a definition: "The oxidation state, often called the oxidation number, is an indicator of the degree of oxidation of an atom in a chemical compound. The formal oxidation state is the hypothetical charge that an atom would have if all bonds to atoms of different elements were 100% ionic. Oxidation states are typically represented by integers, which can be positive, negative, or zero. In some cases, the average oxidation state of an element is a fraction, such as 8/3 for iron in magnetite (Fe3O4)."

Describing a system where the same ion appears more than once, and in different oxidation states, is a major challenge in approximate DFT. This is foremost a self-interaction problem - akin to the splitting of the hydrogen molecular ion H2+(when the two protons are very far apart, the Schroedinger equation tells us that the state with the single electron on the left proton should be degenerate with the state with the single electron of the right proton; and any linear combination of these two states is an admissible solution, and all these solution have the same energy (-1Ry). Move to LDA or GGA, and this degeneracy is lost, and the solution that is favoured is the one where half an electron sits on the left proton, and half an electron sits on the right proton. This per se is not a bad wave function (it's actually a perfectly admissible choice) but the LDA or GGA energy corresponding to it is much lower than -1Ry. In a different language, LDA or GGA favour fractional occupations too much. Hartree-Fock would get it right. A hybrid with xx% of exchange will get it xx% right, i.e. 25% right (this is 75% wrong...).

Going back to oxidations states, take e.g. a supercell with a two iron ions, very far apart, and only 13 valence electrons to share. It's now clear that LDA or GGA will fall in the same hydrogen trap, and provide two Fe ions with oxidation state 2.5+ each. Fluctuations in the environment would break degeneracy (both in the real world, and in approximate DFT), so that the electron would localise on one of the two sites (n the real world), and might help keeping it localised in the approximate DFT world (but not really).

In order to make progress, we introduced a few years ago a simple approach (find it here) based on the fact that approximate DFT gives very good results for atom at integer occupations. E.g., in the H2+ problem, where there is complete degeneracy between electron on the left and electron on the right, and approximate DFT splits it in two (and gets a solution with energy much lower than -1Ry), one could force the electronic ground state to be in the ~integer solution, i.e. force the electron to sit on the right, or to sit on the left, and then one would have a very good estimate of the energy (-1Ry) for states in this degenerate manifold.

Bottom line - in the presence of ions with different oxidation states, approximate DFT will split the electrons across all sites, with an energy that is much lower than the truth. Forcing each electron to be localised on a single site will instead select another admissible wave function, but a much better estimate of the energy. Now, how to localize electrons on an oxidation site? Well, you modify the functional to force that solution. Our choice was to add a penalty for any solution that didn't respect our target oxidation (measured with the matrix of projection on the d orbitals). Note that there is a conceptual basis for this - constrained DFT, as introduced by Peter Dederichs et al (see it here). Other groups (notably Van Voorhis and coworkers, see it here and Scheffler and Behler, here) introduced similar concepts, and applied it to a variety of molecular or solid-state systems, but in my biased opinion the cleanest case, and the one that is better defined, is that of 3d transition metals, or in general d or f ions.

It's fairly simple to implement these constraints or penalties in standard codes (especially those that work with DFT+U, since one has an occupation matrix ready). For Quantum-ESPRESSO, the Car-Parrinello code has still the penalty term implemented, albeit not tested recently (it's based on the DFT+U CP, implemented at the same time, and actually been recently made much more efficient by G. Mattioli), and so there is nothing that works automatically. Luckily, it's fairly simple to implement an exact constraint (less efficient than a penalty, but that's a different matter). Here below is a practical algorithm, in narrative form, kindly provided by David O'Regan (david.o.regan at tcd.ie), a former member now at Trinity College Dublin:

A home-made constrained DFT treatment can be used to enforce the required number of electrons exactly on a given atom in PWscf, for a given definition of that count. All that is needed for the cDFT potential, and hence ionic forces, is really already available using the alpha term in the DFT+U implementation, with the Hubbard U set to zero (or some very small number).

See input flags lda_plus_u and Hubbard_alpha.

This then gives the total energy, for one site and no spin-polarisation

E = E_DFT + alpha Tr[n]

where Tr[n] is the trace of the occupancy matrix defined by pseduoatomic orbitals (by default, other options are possible).

In cDFT you would want, however

E = E_DFT + alpha ( Tr[n] - N_target ) where alpha plays the role of a Lagrange multiplier.

Thus, you would need to take the total energy from PWscf and add to it

- alpha N_target

to get the cDFT energy for a Lagrange multiplier alpha and target occupancy N_target.

The “force” on the Lagrange multiplier, used to optimise it, is the occupancy deviation

\partial E / \partial alpha = Tr[n] - N_target

which allows you to optimise the alpha to make this occupancy deviation vanish (you *maximise* the energy with respect to alpha). To do this, iteratively call PWscf, updating the alpha, until the required Tr[n] - N_target is achieved within some tolerance and the cDFT term vanishes:

\partial E / \partial alpha = 0 --> alpha ( Tr[n] - N_target ) = 0

For a single constrained species, which I guess is the case here, the alpha is a scalar and simple steepest ascents on alpha is enough to optimise it efficiently.

If the Lagrange multiplier alpha is re-optimised for each structure, then the term proportional to d alpha / d R in the forces and stresses vanishes. So, no modification of the PWscf code is needed for cDFT geometry optimisation if both the total energy and Lagrange multiplier alpha are well converged. This implies nested loops of single restarted vc-relax steps, and alpha optimisation.

To summarise, what you wish to do requires no changes to PWscf whatsoever, just some careful bash/Python scripting or multiple manual restarts.