Skip to content
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

FIX: Fix assumptions for substitutional/interstitial sublattices in partitioned order/disorder phases #311

Merged
merged 26 commits into from Mar 25, 2021

Conversation

bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 9, 2021

Checklist

  • The documentation examples have been regenerated if the Jupyter notebooks in the examples/ have changed. To regenerate the documentation examples, run jupyter nbconvert --to rst --output-dir=docs/examples examples/*.ipynb from the top level directory)
  • If any dependencies have changed, the changes are reflected in the
    • setup.py
    • environment-dev.yml

This PR aims to close #310 by changing the heuristic for how pycalphad considers a sublattice in partitioned order/disorder and lays out in the code documentation how we treat the energetic contributions. Previously, the implementation was limited because it:

  1. only considered one interstitial sublattice
  2. that sublattice had to be the last sublattice in the ordered and disordered phase constituents
  3. the sublattice could only contain vacancies

To summarize the new changes that frees all these assumptions (copied from the new documentation):

This method assumes that the first sublattice of the disordered phase is the substitutional sublattice and all other sublattices are interstitial. In the ordered phase, all sublattices with constituents that match the disordered substitutional sublattice will be treated as disordered (with site fractions replaced by mole fractions in the ordered sublattices) and the interstitial sublattices will not have any site fractions substituted.

The only assumptions here is that there is one set of sublattices in the ordered phase that can disorder. Any number of sublattices may disorder and any number of interstitial sublattices are possible.

The main changes made in this PR:

  1. Refactor the Model.atomic_ordering_energy method to clearly split up the two responsibilities of the method: computing the ordering energy and replacing the ordered energy contributions with the disordered energy contributions
  2. Fix test case for Bug in order-disorder models when species other than VA are present in the "vacancy sublattice" #310 and a more complicated case where one of the substitutional elements is also an interstitial element
  3. Revamped documentation for Model.atomic_ordering_energy
  4. Make Model.mole_fraction a private method (Model._mole_fraction) to better reflect its purpose and to discourage users from using it, since it is not applicable outside the order-disorder model. It was not applicable outside the atomic_ordering_energy before, since it counts vacancies as regular species. This behavior is made more explicit by making Model._mole_fraction private. Since its usage is limited to Model.atomic_ordering_energy, I added a parameter to take in the substitutional sublattice indices. Strictly speaking, this could break backwards compatibility to anyone depending on mole_fraction, but mole_fraction does not normalize to mole fraction of atoms since it counts vacancies as atoms and anyone using it for the purpose of computing the mole fraction of atoms should be using Model.moles instead.

@bocklund bocklund changed the title Order disorder va subl FIX: Fix assumptions for substitutional/interstitial sublattices in partitioned order/disorder phases Mar 9, 2021
@codecov
Copy link

codecov bot commented Mar 9, 2021

Codecov Report

Merging #311 (c5d0852) into develop (5038950) will increase coverage by 0.05%.
The diff coverage is 91.66%.

Impacted file tree graph

@@             Coverage Diff             @@
##           develop     #311      +/-   ##
===========================================
+ Coverage    86.78%   86.83%   +0.05%     
===========================================
  Files           46       46              
  Lines         4464     4527      +63     
===========================================
+ Hits          3874     3931      +57     
- Misses         590      596       +6     
Impacted Files Coverage Δ
pycalphad/tests/test_energy.py 97.52% <82.75%> (-2.48%) ⬇️
pycalphad/model.py 91.50% <95.65%> (+0.22%) ⬆️
pycalphad/tests/datasets.py 100.00% <100.00%> (ø)
pycalphad/tests/test_model.py 100.00% <100.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 5038950...c5d0852. Read the comment docs.

Comment on lines 1000 to 1006
This method assumes that the first sublattice of the disordered phase is
the substitutional sublattice and all other sublattices are
interstitial. In the ordered phase, all sublattices with constituents
that match the disordered substitutional sublattice will be treated as
disordered (with site fractions replaced by mole fractions in the
ordered sublattices) and the interstitial sublattices will not have any
site fractions substituted.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know if I got it correctly, but for the ord-dis model, we should take the overall mole fraction and substitute in the ordered site fractions to get the ordering contribution, only if the phase can be totally disordered. But in the case of
A2: (A,B,Va)m+n(C,Va)p
B2: (A,B,Va)m(A,B,Va)n(C,Va)p
The substitution shouldn't be done with mole fraction but with the sub-lattice fraction of the disordered state that represents the disordered state of the sub-lattice sets that undergo disordering of the ordered state. Is this what you mean in this phrase, am I missing something? e. g.:

Gm = Gm_dis + ΔGm_ord
ΔGm_ord = Gm_ord (y1i, y2j, y3k) - Gm_ord (y1i = y2j, y3k)

One reference that used this model:
Acta Materialia 79 (2014) 1–15;
doi: http://dx.doi.org/10.1016/j.actamat.2014.07.006

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The substitution shouldn't be done with mole fraction but with the sub-lattice fraction of the disordered state that represents the disordered state of the sub-lattice sets that undergo disordering of the ordered state.

Thank you for the reference. Most of the original publications for the partitioned order/disorder model do not include the interstitial sublattices and leave that out of their derivations. I actually had trouble tracking down a reference for how the interstitial sublattices are intended to work, so thank you for the recent reference.

TL;DR

I believe this formulation is the same as the one we are using, just with slightly different notation. Please make any suggestions that would help clarify the documentation!

Discussion

That Equation (11) from the work of Phan, Paek, and Kang, Acta Materialia 79 (2014) 1–15 that you cited:

Screen Shot 2021-03-10 at 11 57 25 AM

treats the disordered sublattices the same as the one in Equation (3) of Servant and Ansara, Calphad 25, (2001) 509–525 that is cited in the documentation

Screen Shot 2021-03-10 at 12 02 27 PM

except that the work of Phan, Paek and Kang also specify that the interstitial sublattice should not be substituted (which is left out of Servant and Ansara's derivation, since they did not use an interstitial sublattice).

What the code here intends to do is as follows (you'll have to excuse the fast and loose notation here, since in pycalphad we don't make any assumptions about the number of sublattices so the notation from Phan, Paek and Kang which specifies the sublattices explicitly doesn't translate well):

image

where the "subs" superscript means "on substitutional sublattices" and the "inter" superscript means "on interstitial sublattices". The "mole fraction" from Model._mole_fraction that is used here only is really more of a "quasi mole fraction" rather than a real "mole fraction". It's named in the spirit of the original derivations. Model._mole_fraction replaces the site fractions on substitutional sublattices and does not normalize the mole fraction for vacancy site fractions as a true "mole fraction" would. This is the motivation for making the _mole_fraction method private: because it's not a real mole fraction.

Copy link
Collaborator Author

@bocklund bocklund Mar 11, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For another brief description in the published literature, see section 3.3 of Connetable et al. Calphad 32 (2008) 361-370 doi: 10.1016/j.calphad.2008.01.002, specifically equations 14-18.

I'm going to update the docstrings to be more consistent with how it's described in this paper.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you @bocklund. This was really helpful. I think I could understand it better!

@bocklund
Copy link
Collaborator Author

bocklund commented Mar 11, 2021

One of the TDB files I had uses the Al-C-Fe system from Connetable et al. as a subsystem. I wrote that out to a TDB file.
al-c-fe.tdb.txt

This is my test script, with the energies from Thermo-Calc (site fractions were plugged in via set-start-constitution)

# coding: utf-8
import pycalphad
print(pycalphad.__version__)
from pycalphad import Database, equilibrium, variables as v, calculate, Model
dbf = Database('al-c-fe.tdb')
phases = ['B2_BCC']

# Al-C
comps = ['AL', 'C', 'VA']
cr = calculate(dbf, comps, phases, N=1, P=101325, T=300, pdens=10)
mod = Model(dbf, comps, 'B2_BCC')
print(sorted(comps))
print([sorted([sp.name for sp in subl]) for subl in mod.constituents])
print(cr.isel(points=-1).Y.values.squeeze())
print(cr.isel(points=-1).GM.values.squeeze(), 'TC: 31911.163')

# Al-Fe
comps = ['AL', 'FE', 'VA']
cr = calculate(dbf, comps, phases, N=1, P=101325, T=300, pdens=10)
mod = Model(dbf, comps, 'B2_BCC')
print(sorted(comps))
print([sorted([sp.name for sp in subl]) for subl in mod.constituents])
print(cr.isel(points=-1).Y.values.squeeze())
print(cr.isel(points=-1).GM.values.squeeze(), 'TC: 86629.399')

# C-Fe
comps = ['FE', 'C', 'VA']
cr = calculate(dbf, comps, phases, N=1, P=101325, T=300, pdens=10)
mod = Model(dbf, comps, 'B2_BCC')
print(sorted(comps))
print([sorted([sp.name for sp in subl]) for subl in mod.constituents])
print(cr.isel(points=-1).Y.values.squeeze())
print(cr.isel(points=-1).GM.values.squeeze(), 'TC: 48974.223')

# Al-C-Fe
comps = ['AL', 'FE', 'C', 'VA']
cr = calculate(dbf, comps, phases, N=1, P=101325, T=300, pdens=10)
mod = Model(dbf, comps, 'B2_BCC')
print(sorted(comps))
print([sorted([sp.name for sp in subl]) for subl in mod.constituents])
print(cr.isel(points=-1).Y.values.squeeze())
print(cr.isel(points=-1).GM.values.squeeze(), 'TC: 34659.484')

with outputs:

0.8.4+26.gea5c234a
['AL', 'C', 'VA']
[['AL', 'VA'], ['AL', 'VA'], ['C', 'VA']]
[0.35136561 0.64863439 0.79724251 0.20275749 0.44382726 0.55617274]
31911.16739624921 TC: 31911.163
['AL', 'FE', 'VA']
[['AL', 'FE', 'VA'], ['AL', 'FE', 'VA'], ['VA']]
[0.46828033 0.00228176 0.52943791 0.43587845 0.01705016 0.54707139
 1.        ]
86629.40328004515 TC: 86629.399
['C', 'FE', 'VA']
[['FE', 'VA'], ['FE', 'VA'], ['C', 'VA']]
[0.35136561 0.64863439 0.79724251 0.20275749 0.44382726 0.55617274]
48974.22829354314 TC: 48974.223
['AL', 'C', 'FE', 'VA']
[['AL', 'FE', 'VA'], ['AL', 'FE', 'VA'], ['C', 'VA']]
[0.23632422 0.09387751 0.66979827 0.40269437 0.55906662 0.03823901
 0.12888967 0.87111033]
34707.15643943977 TC: 34659.484

all the binary cases agree, but the ternary does not. It's unclear why, but it's blocking this issue until it's figured out.

@bocklund
Copy link
Collaborator Author

bocklund commented Mar 11, 2021

This seems to be related to magnetism. I commented out the magnetic parameters in the ordered phase in the Al-C-Fe database and now the energies agree. Interestingly[1] the Curie temperature for the "non-magnetic ordered phase" is 301.8 K according to Thermo-Calc (this is also what pycalphad sees).

The atomic_ordering_energy method replaces the TC and BMAG in the ordered phase from the disordered phase. I need to look into it more, but my guess is that the TC and BMAG are parameters should be handled like the other energetic contributions, TC = TCdis + ΔTCord.

[1] This is actually really lucky that my test calculations were at 300 K, because if it were much higher than 300 K, the magnetic energy contribution would diminish and this bug would have gone un-noticed! It's actually unrelated to the issue this PR fixes and may have affected any database that includes magnetic parameters in the ordered phase)

@Ebert-Daniel
Copy link

I think you are right.
In H.L. Lukas, S.G. Fries, and B. Sundman. Computational Thermodynamics: the Calphad Method, volume 131. Cambridge University Press, Cambridge, United Kingdom, 2007. In section 5.8.6, it is said that the physical properties may also have their parameters partitioned between ordered and the disordered part.

@bocklund
Copy link
Collaborator Author

Thanks @Ebert-Daniel! I have read through 5.8.2.4 and 5.8.4 two or three times since investigating #310 and I missed 5.8.6! Definitely key information.

- Set beta = BMAG as a class property so it will be defined.
- Use _partitioned_expr to construct the ordering energy and to
  partition the physical properties (does not affect the Gibbs energy)
- Skip the failing test for the magnetic part for now with a TODO
  comment in the atomic_ordering_energy method
@bocklund
Copy link
Collaborator Author

Fixing the partitioning of the physical properties is going to be slightly more complicated. The changes here will partition the TC, BMAG and NT expressions and update the Model attributes appropriately (e.g. Model().TC will be correct). However, the AST expressions are different objects than the expressions in the attributes, so properly partitioning Model().TC in the ordered phase doesn't fix the Gibbs energy expression.

Using magnetic energy as an example: both the ordered energy and disordered energy need to have the partitioned TC, but there's currently no way for the ordered model to make the disordered magnetic energy use the partitioned TC because the values of the TC are inserted directly into the disordered magnetic energy expression and are not symbolically accessible later. This applies to all physical models, like the Einstein and two-state models.

The partitioning of the Gibbs energy contributions from physical models like the magnetic model is technically orthogonal to the issue this actually solves, so I propose to get this PR merged and raise a new issue to track the partitioning of physical properties. This PR still includes the failing test, but it is skipped for now.

That would make this ready for review.

magnetic property attributes

Add checking Curie temperature + BMAG to skipped test
@bocklund
Copy link
Collaborator Author

The failures are due to coverage upload failures. The tests passed.

pycalphad/model.py Outdated Show resolved Hide resolved
pycalphad/model.py Show resolved Hide resolved
@bocklund
Copy link
Collaborator Author

This should have addressed the outstanding review items and also added a warning from an out-of-band discussion to warn the user if they instantiate a Model that a physical Model.contribution that must be partitioned.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Bug in order-disorder models when species other than VA are present in the "vacancy sublattice"
3 participants