New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Compute metastable/unstable single phase driving forces in ZPF error #151
Compute metastable/unstable single phase driving forces in ZPF error #151
Conversation
This requires the following before it can be merged:
|
I constructed a super simple Cu-Ni comparison in the miscibility gap using two ZPF data points and no other data (thermochemical, ZPF, etc.). Below are the phase diagram comparisons of the optimal parameter set after 500 iterations. The "global min" version uses the |
a31df24
to
82930b0
Compare
Thanks Brandon for bringing up this solution! I am running an MCMC optimization with only this single-phase-ZPF data and no other data. The color of the circles (color map on the right of the above figure) represents the driving force for the individual data points at the start of the optimization. (With Still the MCMC using only the single phase ZPF data changes the PD in a manner that the experimental ZPF data lie in a single phase region of the phase diagram (see figure below). The two points which had previously a large driving force still have a very large driving force compared to the other points. The logarithmic probability went from -295 to -269. All of the above calculations were done with Should one avoid using |
Hi @tobiasspt, thanks for your question. I think what’s happening is that the approximate equilibrium target hyperplane is less accurate than the exact for these points. I’d have to have the database and data in hand to check for sure. Since you have a development version installed, you could try to increase the You give the log probability decrease, what’s the decrease in the driving force? Also, what’s the total probability for all the data? If the driving forces/probabilities of these high driving force points stay relatively constant compared to the other data, it’s probably safe to continue using approximate equilibrium. If I am correct that it’s just poor sampling and a poor target hyperplane in approximate compared to exact, then the driving force and probability will stay relatively constant and not influence whether parameters are accepted or rejected by the Metropolis criteria. |
Also, I’d love to see the code you are using to color code the driving forces like this. That would be an extremely useful debugging tool for me and many others. Have you been extracting the driving force data by hand? I’ve been toying with an idea to add an extra function that returns the driving forces for the data as a list instead of the log probability sum, which would make it easier to generate these type of figure programmatically. |
Thx for the help @bocklund. Increasing the The change in driving force is depicted in the figure below. It is clearly seen that the highest change is in the middle, where at the start of the run the miscibility gap was. This would explain, why the miscibility gap dissapears after the run. Also the driving force of the high driving force points stay relatively constant, suggesting that using The total probability of all the ZPF data at the beginning of the run is -1236 out of which -295 stem from the single-phase-ZPF data. |
I am extracting the driving forces from the ESPEI logfile (verbosity 3) with a simple python function (more or less by hand), which manipultes the text to read out the desired value. This function is not very general and works for the the start of a run. (In order to extract the driving forces for the end of the run, I simply start an MCMC with the optimised parameters with 1 iteration). See the small code below:
A function which writes out the driving forces into a list would indeed be very helpfull |
Huh - that is puzzling indeed. Can you share a TDB (or two: before & after) with me, along with some of those tricky points around 10% W? |
single_phase_zpf.zip |
Thanks for sending the data! This helped uncover the issue. The target hyperplane is well sampled, that's why increasing the point density of the target hyperplane doesn't do much. The problem lies in the sampling of the single phase site fraction space. TL;DRThe problem is that having vacancies in the first sublattice of BCC produces is causing the I would update my reccomendation to "proceed with caution" when using approximate equilibrium when there are independent vacancy site fractions. Moving forward, this is how I see the options:
ExplanationWithout vacancies in the first sublattice, there are no degrees of freedom to sample because the site fractions of Ti and W must be equal to the mole fractions. With vacancies in the first sublattice, it's more complicated. For one data point:
That is, the independent site fraction is the site fraction of W in the first sublattice. Notably, the expression for the site fraction of vacancies in the first sublattice is Right now, It's not clear to me how to achieve sampling around that solution. One idea is to try to change the variables to make all the vacancy-like site fractions independent. That would allow us to sample:
This sampling would catch edge cases like this and efficiently sample the site fractions of vacancy-like species. Unfortunately, SymPy doesn't seem like it lets you specify a "preference" for which variables should be independent in the solution. I spent some time last weekend when I was implementing this trying to change the variables on my own, but I wasn't successful in coming up with a way to do it. We might have to require exact equilibrium for only the driving force calculation in these cases. It will increase the cost a little bit, but the single phase driving force calculations with restricted points will still be much faster than the total cost for exact equilibrium everywhere, since the multi-phase calculations with global minimization are more expensive. (Note that this sampling issue does not affect target hyperplane determination because all the site fractions are sampled independently - the root cause is that the site fractions here are not all independent.) |
With a second look from @richardotis, we were able to come up with a way to prefer vacancy and vacancy-like species as independent symbols in the symbolic solution. After these last two commits, I now get much better agreement: Before the changes:
After the changes:
|
We use the ability to isolate the independent vacancy-like site fractions to be able to sample in logspace near the site fraction limits @tobiasspt can you try pulling the latest commits and trying again? I'd be curious to see how the driving forces of the starting points compare between the exact and approximate. |
I refactored the internals a little to make it easier to calculate the driving forces programmatically. This makes it a lot easier plot driving forces directly in Python. You can also provide your own import numpy as np
import matplotlib.pyplot as plt
from pycalphad import Database, binplot, variables as v
from pycalphad.core.utils import extract_parameters
from espei.datasets import load_datasets, recursive_glob
from espei.error_functions.zpf_error import get_zpf_data, calculate_zpf_driving_forces
# User input variables
TDB_PATH = '/Users/brandon/Downloads/single_phase_zpf/ti-w_end.tdb'
DATASETS_DIR = '/Users/brandon/Downloads/single_phase_zpf/input-data'
COMPS = ['TI', 'W', 'VA']
INDEP_COMP = v.X('W') # binary assumed
CONDS = {v.N: 1, v.P: 101325, v.T: (500, 3800, 20), INDEP_COMP: (0, 1, 0.01)}
CMAP = 'hot'
parameters = {} # e.g. {'VV0001': 10000.0}
approximate_equilibrium = True
# Script below:
dbf = Database(TDB_PATH)
phases = list(dbf.phases.keys())
# Get the datasets, construct ZPF data and compute driving forces
# Driving forces and weights are ragged 2D arrays of shape (len(zpf_data), len(vertices in each zpf_data))
ds = load_datasets(recursive_glob(DATASETS_DIR, '*.json'))
zpf_data = get_zpf_data(dbf, COMPS, phases, ds, parameters=parameters)
param_vec = extract_parameters(parameters)[1]
driving_forces, weights = calculate_zpf_driving_forces(zpf_data, param_vec, approximate_equilibrium=approximate_equilibrium)
# Construct the plotting compositions, temperatures and driving forces
# Each should have len() == (number of vertices)
# Driving forces already have the vertices unrolled so we can concatenate directly
Xs = []
Ts = []
dfs = np.concatenate(driving_forces)
for data in zpf_data:
for phase_region in data['phase_regions']:
for vertex_comp_cond in phase_region.comp_conds:
Ts.append(phase_region.potential_conds[v.T])
# Binary assumptions here
assert len(vertex_comp_cond) == 1
if INDEP_COMP in vertex_comp_cond:
Xs.append(vertex_comp_cond[INDEP_COMP])
else:
# Switch the dependent and independent component
Xs.append(1.0 - vertex_comp_cond[INDEP_COMP])
# Plot the phase diagram with driving forces
fig = plt.figure(dpi=100)
ax = fig.gca()
binplot(dbf, COMPS, phases, CONDS, plot_kwargs={'ax': ax}, eq_kwargs={'parameters': parameters})
sm = plt.cm.ScalarMappable(cmap=CMAP)
sm.set_array(dfs)
ax.scatter(Xs, Ts, c=dfs, cmap=CMAP, edgecolors='k')
fig.colorbar(sm, ax=ax, pad=0.25)
fig.show() Here's the ending point TDB that you sent with approximate and exact equilibrium, respecitvely: |
It's part of the stdlib in Python 3.3+
5d360b4
to
7766946
Compare
@tobiasspt I rebased this on |
Thanks a lot for the detailed explanation and the help with my question! Comparison before the changes:
Comparison after the changes:
|
Small question. It seems that with the new commits in my logfiles the driving force for the ZPF data is printed, but the probability is omitted. This is not much of a problem for me right now and I could have made some mistake while getting the new commits installed. |
Part of the refactoring was to move all the probabilities to be calculated in bulk from the driving forces, so they are now logged at the end of calculating the ZPF error in each step in something like:
I tried to move it out so it was a little cleaner in the code as far as where the driving forces were calculated and where the probabilities were calculted, but maybe it would make more sense to do them together like we used to. Instead of |
Acutally I think that the new behaviour of calculating the probabilities from the driving forces in bulk after the ZPF error is a good solution. Thx agaun for the answer! |
Thanks for the feedback 😄 - it helps to have more fresh eyes look at this. |
* Make `PhaseRegion` a dataclass * Introduce a `RegionVertex` dataclass instead of storing the vertex data in lists of `PhaseRegion` * Fix a regression in ZPF data introduced by #151 where prescribed phase compositions of stoichiometric phases that used to work no longer work because the phase composition of a stoichiometric phase may be unsatisfiable. Now there's a check for stoichiometric phases that will not try to solve for the points. * Fixes a bug where stoichiometric phases used `equilibrium` in driving force calculations which could give bad energies for exact equilibrium where mass balance could not be satisfied (approximate was not affected). Now exact and approximate equilibrium use driving force estimation via `calculate` which is always exact for stoichiometric phases because they have exactly one point.
Thanks to Tobias Spitaler for suggesting this and to @richardotis for brainstorming this solution concept.
This PR introduces two new functions in ZPF error,
_solve_sitefracs_composition
and_sample_solution_constitution
. Their purpose is to facilitate computing metastable or unstable single phase driving forces when a phase has a miscibility gap. This should improve the convergence for any phase that has a stable or metastable miscibility gap.Rationale
ESPEI currently computes the "single-phase hyperplane" at a vertex by performing an equilibrium calculate at a black point and then subtracting that from the target hyperplane energy at that composition. As illustrated in the figure Tobias constructed (below), this is problematic for phases with a miscibility gap because a "single-phase" equilibrium calculation in pycalphad will always compute the global minimum energy and give two composition sets.
What ESPEI should do is what Tobias illustrates by the orange
x
and the green driving force line. This solution ensures that minimizing the driving force will force the Gibbs energy curve to match the energy of the black points on the multi-phase target hyperplane.Historically, we didn't implement this because one would like to use
equilibrium
to minimize the internal degeres of freedom, but pycalphad always computes the global minimum energy, so it was not possible to do viaequilibrium
. More recently, ESPEI had introduced the idea ofapproximate_equilibrium
, which usesstarting_point
to more quickly determine a minimum energy solution from a discrete point smapling grid. Theapproximate_equilibrium
method we use still has the same problem as pycalphad'sequilibrium
becausestarting_point
will still give the global minimum solution for the discrete sampling.Solution
In an ideal world, pycalphad should be able to turn off global minimization (automatically introducing new composition sets) and enable a condition to be set for the composition of a phase, i.e.
X(BCC,B)
. In practice, being able to turn off global minimum and provide a valid starting point for only one composition set that has a global composition condition would simulate a phase composition condition. Unfortunately, neither turning off global minimization nor phase composition conditions are currently implemented. So we need to do a workaround.The two functions introduced here consider each single phase composition at a tie-vertex and construct a point grid that only contains points which satisfy the prescribed overall composition (and the internal phase constraints). This can be used in either approximate or exact equilibrium modes to find lowest energy starting point and then to pass that
equilibrium
with the constrained point grid so the global minimization step has no new composition sets to introduce (i.e. it cannot detect a miscibility gap).For perfomance, we pre-compute the grid of points for every phase composition in the ZPF datasets and re-use them to compute the grid, starting point and equilibrium at every parameter iteration (note that this would be invalid if a parameter changes the number of moles, like varying coordination number in the MQMQA).
To summarize the impact: