Hi @johan28 , sorry for my slow response and thank you for the deep thought you’ve put into this. In general I agree with you; I think I misunderstood you a few posts back and caused some confusion about formation vs. free energy. After thinking about this for a while, I also think you have found an error in how this correction is implemented.
Let me re-state my understanding of how
E denote DFT (electronic energy)
dE denote an energy correction applied to DFT energy
g denote free energy (as in chemical potential, not formation free energy)
dG note the formation free energy change
First, we add room-temperature entropies to the corrected DFT energies of water and O2,
(E_H2O + dE_H2O) and
(E_O2+dE_O2) , to get g_H2O and g_O2
Next, we use these free energies to calculate the corrected free energy of H2, such that the formation free energy of water matches the experimental value
g_H2O - 1/2 g_O2 - 2 g_H2= MU_H2O = dG_H2O
Now for the hydrate adjustment. The purpose of this adjustment is that when you calculate the formation free energy of FeO.nH2O, the
H2O part “contributes” the correct amount of free energy, equal to
n dG_H2O or
n MU_H2O. Another way of thinking about this is that we want to treat the free energy (g_FeO.nH2O) as a combination of g_FeO (hyd) and n g_H2O:
g_FeO.nH2O = g_FeO + n g_H2O
But the solid compound FeO.nH2O may have had other energy corrections applied that are not consistent with this objective. For example, in the MP2020 compatibility scheme there are corrections for hydrides that really aren’t appropriate if the compound is a hydrate. The compound above may have had its energy adjusted by
n * hydride correction.
So we need to get rid of any previous energy corrections on the
O that are associated with water. After that, we need to ADD any energy corrections that were applied to a single water molecule, e.g.
g_FeO.nH2O = g_FeO.nH2O - n dE_O - 2n dE_H + n dE_H2O
Note that if there are no corrections applied to hydrides, oxides, or water, then this energy adjustment is zero. But that won’t be the case when using MP data because H2O will not receive a
hydride correction (see our corrections file).
Now that we have adjusted the free energy of the embedded H2O, we can verify that dG_FeO.nH2O = dG_FeO + n dG_H2O (as intended) as follows:
dG_FeO.nH2O = g_FeO.nH2O - g_Fe - (n+1)/2 g_O2 - 2n g_H2
dG_FeO.nH2O = g_FeO.nH2O - g_Fe - 1/2 g_O2 - n ( 1/2 g_O2 + 2 g_H2 )
( 1/2 g_O2 + 2 g_H2 ) = dG_H2O - g_H2O, from above, so
dG_FeO.nH2O = g_FeO.nH2O - g_Fe - 1/2 g_O2 - n ( dG_H2O - g_H2O )
Now, substitute the definition that we started with: g_FeO.nH2O = g_FeO (hyd) + n g_H2O
dG_FeO.nH2O = g_FeO (hyd) + n g_H2O - g_Fe - 1/2 g_O2 - n ( dG_H2O - g_H2O )
dG_FeO.nH2O = [ g_FeO (hyd) - g_Fe - 1/2 g_O2 ] + n [ g_H2O - n ( dG_H2O - g_H2O ) ]
dG_FeO.nH2O = dG_FeO (hyd) + n dG_H2O
Now this is not what currently happens in
pymatgen, and I think that’s an error. Specifically, I think this line
hydrate_adjustment = -1 * (self.h2o_adjustments * 3 + MU_H2O)
Should be something like
hydrate_adjustment = - dE_O - 2 dE_H + self.h2o_adjustments * 3
dE_H are the corrections applied to oxides or hydrides in the material, and
dE_H2O, and all values are in eV/atom.
For MP data,
dE_H2O = dE_O / 3 (i.e., a water molecule receives one oxide correction), so the line would simplify to
hydrate_adjustment = - 2 ( dE_H )
So basically, the hydrate correction boils down to just removing any hydride energy adjustments from the solid phase.
Does the above make sense to you? If so, would you be able to do some testing to see how correcting this adjustment affects some representative Pourbaix diagrams? Flagging @Joseph_Montoya in case he has any thoughts as well.
If this looks like it makes sense, I can open an issue in
pymatgen and we can start working on a fix.