Pycalphad most efficient way to gather equilibirum phase data for set of alloy compositions?

Question

What's the most efficient way to get the equilibrium phases and phase fractions for a set of alloy compositions? For context, I'm coming from Thermo-Calc/TC-Python where I would just loop through the alloy compositions and do a single equilbrium calculation. However, when I attempt to do this with pycalphad it seems quite slow. For example, here is what my code for the "one at a time" looping:

from pycalphad import Database, calculate, equilibrium, variables as v
dbf = Database("Cr-Fe-Ni.tdb")
components = ['FE','CR', 'NI', 'VA'] 
phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'SIGMA']

#Example set of alloy compositions, normally I would want to sample more finely but I'll keep it simple for demonstration purposes
compositions = array([[0. , 0. ],
                      [0. , 0.2],
                      [0. , 0.4],
                      [0. , 0.6],
                      [0. , 0.8],
                      [0.2, 0. ],
                      [0.2, 0.2],
                      [0.2, 0.4],
                      [0.2, 0.6],
                      [0.2, 0.8],
                      [0.4, 0. ],
                      [0.4, 0.2],
                      [0.4, 0.4],
                      [0.4, 0.6],
                      [0.6, 0. ],
                      [0.6, 0.2],
                      [0.6, 0.4],
                      [0.8, 0. ],
                      [0.8, 0.2]])

phase_data = []
for i in tqdm.tqdm(range(len(compositions)), total=len(compositions), position=0, leave=True, smoothing=0.3):
    eq_result = equilibrium(dbf, components, phases, {v.X('CR'):compositions[i,0], v.X('NI'):compositions[i,1], v.T:1000, v.P:101325})
    phase_data.append([*eq_result.NP.values.squeeze().tolist(), *eq_result.Phase.values.squeeze().tolist()])

This takes about 15 seconds, or about 1.2 alloys evaluated per second. This seems quite slow for a ternary, and I suspect it's because I'm calling the equilibrium function each time instead of just changing the composition.

I tried feeding the whole column at once (see code below), but this seemed to take the Cartesian product of the compositions, which is not what I want.

#database, phase, and composition setup is same as in above code

eq_result = equilibrium(dbf, components, phases, {v.X('CR'):compositions[:,0], v.X('NI'):compositions[:,1], v.T:1000, v.P:101325})
 

So, please let me know what the better approach is. Thanks!

Answer 1

In your comment, you mentioned that you may be testing many alloys (100,000s). I will assume that interested in taking more of a screening approach, rather than specifying precisely 100,000 alloys compositions.

In that case, it is easy and faster to let pycalphad broadcast the conditions for you. It will skip any conditions that have independent mole fractions that sum to greater than unity, e.g. X(CR)=0.8 and X(NI)=0.8. The result that you get from equilibrium is an xarray Dataset object and uses padded NumPy arrays. Floating point values (like NP) will be padded with nan and string values (like Phase) will be padded with '' when no phases are present.

xarray and most of the tools in the ecosystem (like matplotlib) play nicely with these padded values. Below, I give an example where equilibrium calculations are performed on a dense grid and I plot the phase fractions of FCC over the whole grid and under a constraint where the composition of Cr is under 40 at%.

from pycalphad import Database, calculate, equilibrium, variables as v
import matplotlib.pyplot as plt

dbf = Database("Cr-Fe-Ni.tdb")
components = ['FE','CR', 'NI', 'VA'] 
phases = ['LIQUID', 'BCC_A2', 'FCC_A1', 'SIGMA']

eq_result = equilibrium(dbf, components, phases, {v.X('CR'): (0, 1, 0.1), v.X('NI'): (0, 1, 0.1), v.T:1000, v.P:101325})


# verify what the Phase and NP look like by printing
# note these are multi-dimensional arrays and are padded with '' and nan, respectively
print(eq_result.Phase.values)
print(eq_result.NP.values)


# Plot FCC phase fractions over the whole grid
plt.figure()
plt.subplot(projection="triangular")  # projection provided by pycalphad

mask = eq_result.Phase == "FCC_A1" # choose only FCC, others will be NaN
plt.scatter(
    eq_result.X_CR.where(mask).broadcast_like(mask),
    eq_result.X_NI.where(mask).broadcast_like(mask), 
    c=eq_result.NP.where(mask)
)
plt.title("Cr-Fe-Ni FCC Phase Fractions")
plt.xlabel("X(CR)")
plt.ylabel("X(NI)")
plt.colorbar(label="NP(FCC_A1)")
plt.show()


# Plot FCC phase fractions only at compositions under 40 at% Cr
plt.figure()
plt.subplot(projection='triangular')  # projection provided by pycalphad

# choose only FCC AND where the composition of Cr is less than 40%, others will be NaN
mask = (eq_result.Phase == 'FCC_A1') & (eq_result.X_CR < 0.4)
plt.scatter(
    eq_result.X_CR.where(mask).broadcast_like(mask),
    eq_result.X_NI.where(mask).broadcast_like(mask), 
    c=eq_result.NP.where(mask)
)
plt.title("Cr-Fe-Ni FCC Phase Fractions")
plt.xlabel("X(CR)")
plt.ylabel("X(NI)")
plt.colorbar(label="NP(FCC_A1)")
plt.show()

This produces the following two plots: