Unable to use activity data in binary Fe-C with Graphite as reference state

[tobiasspt]

Hi,

We are currently trying to use activity data for Fe-C. Lobo1976 measured the activity of C in alpha-iron relative to Graphite as the standard state, but get erroneous results. (Lobo, Joseph A., and Gordon H. Geiger. "Thermodynamics and solubility of carbon in ferrite and ferritic Fe-Mo alloys." Metallurgical Transactions A 7.8 (1976): 1347-1357.)

I have added the input file below. With this input file, we get chemical potential difference: [nan] (verbosity 3 output).
Is the input file correct or are we missing something?
I have had a look at the value of ref_result within the activity_error.py and this does give only nan results for the specified reference state. Graphite only has C as a component. An equilibrium calculation of Graphite specifying x.V('C') gives an error as Number of dependent components different from one. Can this cause an error here as well?
Used versions: espei: 0.8.6 and pycalphad 0.9.2.
I have added a zip-file with the TDB file and espei input files which reproduces this behaviour.

Thank you for your help,
Tobias

{
        "components": ["FE", "C", "VA"],
        "phases": ["BCC_A2", "GRAPHITE"],
        "weight": 1000,
        "reference_state": {
                "phases": ["GRAPHITE"],
                "conditions": {
                        "P": 101325,
                        "T": 1056.15,
                        "X_C": 1

                }
        },
        "conditions": {
                "P": 101325,
                "T": 1056.15,
                "X_C": [0.00013017]
        },
        "output": "ACR_C",
        "values": [[[0.087]]
                ],
        "reference": "Lobo1976_1056K",
        "meta_data": {
                "DOI": "10.1007/BF02658820",
                "literature reference": "Thermodynamics and Solubility of Carbon in Ferrite and Ferritic Fe-Mo Alloys",
                "table/figure": "table 1",
                "measured data": "C-activity in Alpha-Iron",
                "experimental details": "not available",
                "weight": "default"
        }
}

minimal_example.zip

[bocklund]

Thanks for the great MWE, as always. I set up this calculation in pycalphad and can confirm that the solver doesn't converge.

from pycalphad import Database, equilibrium, variables as v
dbf = Database("Fe-C_matcalc.tdb")
comps = ["FE", "C", "VA"]
phases = ["GRAPHITE"]
eq = equilibrium(dbf, comps, phases, {v.P: 101235, v.T: 1056.15, v.X("C"): 1, v.N: 1}, verbose=True)
print(eq)

I suspect it's from doing this calculation at the edge of composition space, where there's no mass of Fe and the Fe chemical potentials become undefined. There are some rudimentary cases where pycalphad can nudge you off the edge of composition space, but that won't help here because there's no Fe present in any of the phases in the chosen reference state (just graphite).

There's undocumented support in other places in the codebase for defining the reference state per element as something like the following, but it's not supported for activity data yet.

    "reference_states": {
        "C": {"phase": "GRAPHITE", "T": 1056.15, "P": 101325},
        "FE": {"phase": "BCC_A2", "T": 1056.15, "P": 101315}
    },

I will try to add support for using these kind of reference states for activity data, such that the equivalent equilibrium calculation for the reference state would be:

from pycalphad import Database, equilibrium, variables as v
dbf = Database("Fe-C_matcalc.tdb")
phases = ["GRAPHITE"]
eq = equilibrium(dbf, ["C", "VA"], ["GRAPHITE"], {v.P: 101235, v.T: 1056.15, v.N: 1})
print(eq)