Add `strict_anions` option to `MaterialsProject2020Compatibility`

[mkhorton]

Issue

This PR relates to the MaterialsProject2020Compatibility scheme, and cases in which structures might be heavily over-stabilized in error, leading to unrealistic structures potentially being placed on the convex hull.

In brief,

• Anion corrections are applied in cases where an element may not be an anion.
• This can happen either because (1) oxidation states are not known, and it is assumed to be an anion based on the element being the most electronegative element present or (2) oxi_state_guesses might give a very small fractional oxidation state in the range (-1, 0).
  • Case (1) was decided as a conservative choice due to unreliability of oxidation state guessing.
  • Case (2) I think was not considered, since this is not typically encountered in realistic materials.

During exploratory materials discovery, hypothetical structures might be generated that will encounter these cases. As an example, consider the composition CuI9, which results in the oxidation state guess of {'Cu': 3.0, 'I': -0.33333}. This becomes heavily over-stabilized by the correction scheme. This might lead to the hypothetical structure then being placed on the convex hull in error.

The experimental data used to fit the correction scheme did not consider any oxidation state in the domain (-1, 0), and so the corrections should be presumed invalid unless appropriately benchmarked. I want to acknowledge cases of anions which are > -1 such as suboxides, metal-rich pnictides, chacogenides, etc. do exist but as far as I can tell we did not have these in the experimental data used to fit the correction scheme.

Suggested Resolution

This PR adds a strict_anion keyword argument which is set to True by default. This is a breaking change. This will skip applying the anion correction if oxidation state is supplied and the oxidation state is > -1, i.e. the element is not present as an anion that is included in the experimental data used for fitting, and might not be acting as an anion at all.

However, this does not resolve the case where oxidation states could not be determined. Therefore, I also suggest we might consider a more significant breaking change to only apply the anion corrections where oxidation state is supplied.

Request for Comment

Notifying MP correction folks for option to comment: @awvio, @rkingsbury, @mattmcdermott, @computron, @shyamd

[mkhorton]

Notifying @munrojm too, since this would affect MP database builds.

[mkhorton]

To remind ourselves of the anions including in the experimental data used to fit the corrections, I re-generated this list using the same oxi_state_guess logic used within the correction scheme:

Element	Oxidation State	Count
O	-2	385
Cl	-1	100
F	-1	74
S	-2	73
Br	-1	52
I	-1	51
Si	-4	44
N	-3	35
C	-4	30
Se	-2	29
Pt	-2	28
B	-3	28
Te	-2	22
H	-1	22
Au	-1	16
P	-3	16
As	-3	8
Sb	-3	6
Sb	-2	5
C	-2	5
Os	-2	5
Os	-1	3
P	-1	2
S	-1	2
P	-2	1
Se	-1	1
C	-3	1

This list will not be complete.

The strictest application of the correction scheme would limit to just these oxidation states, excepting the oxides which have special handling for peroxides/superoxides.

[mkhorton]

Illustration of Consequences of Over-Stabilization

To better illustrate the issue, attached is a graph generated from a set of new structures that are calculated to be on the convex hull when using the correction scheme. The area of concern is this region between a reported -1 oxidation state and 0 oxidation state.

The algorithm used to generate the structures shown in this plot heavily samples this region because the correction scheme so heavily stabilizes the structures, but manual inspection suggests the structures are not physically reasonable.

The y-axis is not particularly meaningful since it is the total energy adjustment, but to give a sense of the consequence, here is the total number of "anions" as a fraction of total number of sites in the structure: structures with many purported anions are heavily over-represented. As a reminder, all of these structures are apparently on the convex hull as a result of the applied corrections.

Also tagging @esoteric-ephemera too, since you may have opinions -- trying to tag a lot of people in this issue because downstream effects of any change may be significant.

Further Thoughts on Resolution

I am thinking perhaps we should be even more strict, and only apply the empirical corrections when the oxidation state is positively identified as an oxidation state in the experimental data used for fitting the corrections. However, I am aware there are many edge cases in the MP database that will be affected by this change, so this will have a subtle knock-on effect in the convex hulls of many chemical systems. However, I think this is a clear bug that does need to be resolved in some way.

[rkingsbury]

Thanks for the thoughtful and detailed explanation @mkhorton . I am onboard with your suggestion to not apply corrections if the oxidation state is > -1 (though I would defer to @munrojm and other MP staff due to effects on database builds).

I am thinking perhaps we should be even more strict, and only apply the empirical corrections when the oxidation state is positively identified as an oxidation state in the experimental data used for fitting the corrections.

I see value in this. What if instead of a boolean, we make strict_anions a multi-level int where

• 0 = current behavior
• 1 (new default) = ignore anions with oxi state >-1
• 2 = ignore anions when oxi state is not in the experimental data

That would make it quite easy to test the effects of different choices on hulls during builds, and potentially to let users make different choices depending on their specific needs.

However, this does not resolve the case where oxidation states could not be determined. Therefore, I also suggest we might consider a more significant breaking change to only apply the anion corrections where oxidation state is supplied.

The current behavior when oxi states are missing is to apply anion corrections IFF the anion is the most electronegative element in the structure. Maybe we add an option 3 that disables this?

[mkhorton]

I concur @rkingsbury, we can have choose more explicit options rather than the int for readability but otherwise I agree. Maybe a Literal["no_check", "require_bound", "require_exact"] to be a bit clearer rather than the int, and maybe a separate flag require_oxidation_states (i.e., do not apply correction if no oxidation state provided).

The most difficult question is which should be the default I think.

[rkingsbury]

I concur @rkingsbury, we can have choose more explicit options rather than the int for readability but otherwise I agree. Maybe a Literal["no_check", "require_bound", "require_exact"] to be a bit clearer rather than the int, and maybe a separate flag require_oxidation_states (i.e., do not apply correction if no oxidation state provided).

Good ideas, no objections here

The most difficult question is which should be the default I think.

Here, I think we can be pragmatic and make the choice based on what will yield the most accurate hulls overall for the MP database. The formation energy benchmarking dataset we created (part of the r2SCAN workflow paper, also in matbench I think) could be a simple heuristic.

The other philosophy would be to use what I said in the previous paragraph to determine the default setting for MP builds only, but then set the default kwarg in pymatgen to either 1) preserve backward compatibility OR 2) for maximum correctness (i.e., no correction unless oxi states are provided).

[mkhorton]

Thanks Ryan.

For MP, the discriminator here is if it is mostly sourced from experimental crystal structures, or if the database also includes hypothetical structures too. I have seen discussion about MP accepting third party calculations for example, and if this happens it might pollute the hull considerably if this issue is not addressed. In other words, I think MP has "gotten away with" the correction scheme as it currently is, but the new default should work in the more general case.

All this to say, I am thinking maximum correctness should be the default, and maybe limiting to applying the correction only for identified, common oxidation states, and not applying the correction for, say, any assigned non-integer oxidation state.

However I would like MP folks to weigh in here. I am happy to merge this as a fix, but I do not want to cause problems for downstream usage.

[mattmcdermott]

@mkhorton I agree that with the increase in generative models we almost certainly need better handling of corrections and you presented a pretty compelling case for excluding anion oxidation states >-1.

I am in favor of @rkingsbury's suggestion to include the strict option as well. However is there any chance this could lead to weird changes in the hull if there is any instability in oxidation state determination? For example, are there any experimentally stable structures where we cannot reliably determine the oxidation states (or they are close to an integer but not exactly an integer)? Would this mean these phases were not compatible with the correction scheme and thus would be essentially "removed" from the hull? If there were even a few of these examples it would make me more hesitant to make the stricter options the default behavior

[JaGeo]

@mkhorton I agree that with the increase in generative models we almost certainly need better handling of corrections and you presented a pretty compelling case for excluding anion oxidation states >-1.

I am in favor of @rkingsbury's suggestion to include the strict option as well. However is there any chance this could lead to weird changes in the hull if there is any instability in oxidation state determination? For example, are there any experimentally stable structures where we cannot reliably determine the oxidation states (or they are close to an integer but not exactly an integer)? Would this mean these phases were not compatible with the correction scheme and thus would be essentially "removed" from the hull? If there were even a few of these examples it would make me more hesitant to make the stricter options the default behavior

Yes, there are. Al2O. Mostly because +1 is not an allowed Oxidation state. Surely, it is a suboxide but it does exist See https://next-gen.materialsproject.org/materials/mp-8022?formula=Al2O

There will be more cases like this.

[mkhorton]

Thanks @mattmcdermott, @JaGeo, welcome any suggestions here.

For example, are there any experimentally stable structures where we cannot reliably determine the oxidation states (or they are close to an integer but not exactly an integer)?

I would distinguish these cases. The former, certainly. The latter (close to integer but not exactly integer) I think more unlikely?

Yes, there are. Al2O. Mostly because +1 is not an allowed Oxidation state. Surely, it is a suboxide but it does exist

Agree, I mentioned suboxides in the original post, there are other examples too. The correction scheme already distinguishes between peroxides/superoxides/ozonides/regular oxides because there was enough experimental data for materials containing oxygen to fit these and because the correction is meaningfully different as the oxidation state changes. However, it does not make this distinction for other anions due to lack of experimental data.

The question really is if the correction would still be suitable for these more uncommon oxidation states. What is our null hypothesis? I think if a correction is fitted for an ion in one particular oxidation state, there is no reason to assume this correction would be transferrable to other oxidation states, and probably good physical reasons to suggest it should not be, as well as the example of e.g. oxygen that demonstrates this.

However is there any chance this could lead to weird changes in the hull if there is any instability in oxidation state determination? [...] Would this mean these phases were not compatible with the correction scheme and thus would be essentially "removed" from the hull?

Going back to @mattmcdermott's question, I'd be surprised if this didn't have some kind of unexpected side effect somewhere. It's definitely a concern. With that said, we are weighing a known error that is already over-stabilizing some compounds, against a possible error which we have not yet observed, so maybe we fix the known error first and then search for additional issues after that? Noting that the examples I provided in the plots above were not even from a generative model, but from a more traditional structure search.

[mkhorton]

Adding that there may be better ways of doing oxidation state assignment too, to minimize this risk.

Notably, three examples:

1. The improvement offered by doped and @kavanase
2. The BERTOS model
3. The electrochemical series method from Tim Mueller et al. See https://oxi.matr.io, they say a library will be available too but the link currently does not work (@montoyjh ?)

[JaGeo]

Adding that there may be better ways of doing oxidation state assignment too, to minimize this risk.

Notably, three examples:

1. The improvement offered by doped and @kavanase
2. The BERTOS model
3. The electrochemical series method from Tim Mueller et al. See https://oxi.matr.io, they say a library will be available too but the link currently does not work (@montoyjh ?)

Yes, I think it would be great to improve this.

[Andrew-S-Rosen]

Adding that there may be better ways of doing oxidation state assignment too, to minimize this risk.

An approach that I like, but have not seen applied to inorganic materials, is training an ML classifier based on the user-reported oxidation states in the Cambridge Structural Database (which now includes ICSD structures) as originally done by Jablonka et al.. As with all things metadata, it is by no means perfect, but it is pretty good (certainly better than BVS). Of course, this is more of a research project than a quick implementation. The model will also still be biased against properly handling very odd structures that may be generated by a diffusion model.

[JaGeo]

Adding that there may be better ways of doing oxidation state assignment too, to minimize this risk.

An approach that I like, but have not seen applied to inorganic materials, is training an ML classifier based on the user-reported oxidation states in the Cambridge Structural Database (which now includes ICSD structures) as originally done by Jablonka et al.. As with all things metadata, it is by no means perfect, but it is pretty good (certainly better than BVS). Of course, this is more of a research project than a quick implementation. The model will also still be biased against properly handling very odd structures that may be generated by a diffusion model.

Yup. It would be a nice project. (There is a todo list 😅.) But the most critical point would be the dataset.

[mkhorton]

@Andrew-S-Rosen are you sure this has not been done by some of what I linked above? It depends how broadly you interpret "classifier", but e.g. the BERTOS work is a transformer trained with annotated oxidation states from ICSD. I think I prefer Tim's method though, simply for its chemical interpretability.

If you mean structural-based assignment rather than composition-based I see your point. I'm not sure if the bond valence method in pymatgen is simply using an old dataset too, so perhaps performance can also be improved there.

[JaGeo]

I think we would certainly need a structure-based approach.

[Andrew-S-Rosen]

I'm biased by a world where composition is meaningless and the only thing that matters is structure. A MOF, after all, is basically all C, H, and O when your squint at it. 😉

That said, I wouldn't be surprised if it has been done. Thanks for the references too!

[mkhorton]

Point very well taken! :)

[montoyjh]

3. The electrochemical series method from Tim Mueller et al. See https://oxi.matr.io, they say a library will be available too but the link currently does not work (@montoyjh ?)

Will ask. Also @TimMueller-TRI-ctr may be able to answer.

[mkhorton]

Will merge this. I have added the new default, require_bound, as originally proposed. I have also added require_exact as a more conservative option that may be a more sensible default, but which I feel requires more due diligence before being made a default.

I have some concerns that regardless of the exact value of the correction, if it is oxidation state dependent, there will always be edge cases where the correction leads to distorted and/or inaccurate convex hulls.