Source code for pycalphad.models.model_mqmqa

import itertools
from typing import List, Union
from collections import OrderedDict
from functools import partial
from symengine import log, S, Symbol
from tinydb import where
from pycalphad.model import _MAX_PARAM_NESTING
import pycalphad.variables as v
from pycalphad.core.utils import unpack_components, wrap_symbol
from pycalphad import Model
from pycalphad.core.errors import DofError


def _get_species(i, j, k, l) -> v.Species:
    """Return a Species for quadruplet given by constituents

    Canonicalizes the Species by sorting among the cation and anion sublattices.

    Parameters
    ----------
    constituents : List[Union[str, v.Species]]

    """
    constituents = [i, j, k, l]
    if all(isinstance(c, v.Species) for c in constituents):
        # if everything is a species, get the string names as constituents
        constituents = [c.name for c in constituents]
    if len(constituents) == 4:
        constituents = sorted(constituents[0:2]) + sorted(constituents[2:4])
    name = "".join(constituents)
    constituent_dict = {}
    # using get() will increment c instead of overriding if c is already defined
    for c in constituents:
        constituent_dict[c] = constituent_dict.get(c, 0.0) + 1.0
    return v.Species(name, constituents=constituent_dict)


[docs]class ModelMQMQA(Model): """ Symbolic implementation of the modified quasichemical model in the quadruplet approximation developed by Pelton _et al._ [1]_. The formulation here largely follows the derivation by Poschmann _et al._ [2]_. Following pycalphad convention, ``ModelMQMQA.G`` defines the energy for one "formula unit" of this phase. The MQMQA formulations by Pelton _et al._ and Poschmann _et al._ do not formally define the energy per formula unit in the same way as is done for the CEF. However, it is convienient to define one formula unit of an MQMQA phase to be the energy corresponding to one mole of quadruplets. In terms of the equations from Poschmann _et al._, which are largely followed in this class, that means :math:`\\sum_{ab/xy} n_{ab/xy} = 1` and therefore :math:`n_{ab/xy} = X_{ij/kl}` by definition. This class is only semantically a subclass of ``Model``. It implements the API expected for a Model in a self-contained way without any need to rely on the Model superclass, but it is created as a subclass to satisfy various ``isinstance`` checks in the codebase. The subclassing on Model should be removed once a suitable abstract base class or protocol is defined. References ---------- .. [1] Pelton, Chartrand, and Eriksson, The Modified Quasi-chemical Model: Part IV. Two-Sublattice Quadruplet Approximation, Metallurgical and Materials Transactions A, 32(6) (2001) 1409-1416 doi: `10.1007/s11661-001-0230-7 <https://doi.org/10.1007/s11661-001-0230-7>`_ .. [2] Poschmann, Bajpai, Fitzpatrick, and Piro, Recent developments for molten salt systems in Thermochimica, CALPHAD 75 (2021) 102341 doi: `10.1016/j.calphad.2021.102341 <https://doi.org/10.1016/j.calphad.2021.102341>`_ """ contributions = [ ("ref", "reference_energy"), ("idmix", "ideal_mixing_energy"), ("xsmix", "excess_mixing_energy"), ] def __init__(self, dbe, comps, phase_name, parameters=None): self._dbe = dbe self._reference_model = None self.components = set() self.constituents = [] self.phase_name = phase_name.upper() phase = dbe.phases[self.phase_name] self.site_ratios = tuple(list(phase.sublattices)) active_species = unpack_components(dbe, comps) constituents = [] for sublattice in dbe.phases[phase_name].constituents: sublattice_comps = set(sublattice).intersection(active_species) self.components |= sublattice_comps constituents.append(sublattice_comps) self.components = sorted(self.components) # create self.cations and self.anions properties to use instead of constituents if len(constituents) == 2: self.cations = sorted(constituents[0]) self.anions = sorted(constituents[1]) else: raise ValueError(f"Exactly two sublattices are required for the MQMQA. Got {len(constituents)} sublattices: {constituents}.") # In several places we use the assumption that the cation lattice and # anion lattice have no common species; we validate that assumption here shared_species = set(self.cations).intersection(set(self.anions)) assert len(shared_species) == 0, f"No species can be shared between the two MQMQA lattices, got {shared_species}" # Verify that this phase is still possible to build if len(self.cations) == 0: raise DofError(f"{self.phase_name}: Cation sublattice of {phase.constituents[0]} has no active species in {self.components}") if len(self.anions) == 0: raise DofError(f"{self.phase_name}: Anion sublattice of {phase.constituents[1]} has no active species in {self.components}") self._pairs = list(itertools.product(self.cations, self.anions)) quads = list(itertools.product( itertools.combinations_with_replacement(self.cations, 2), itertools.combinations_with_replacement(self.anions, 2), )) self._quadruplets = [(A, B, X, Y) for (A, B), (X, Y) in quads] quad_species = [_get_species(A, B, X, Y) for (A, B), (X, Y) in quads] self.constituents = [sorted(quad_species)] # Convert string symbol names to symengine Symbol objects # This makes xreplace work with the symbols dict symbols = {Symbol(s): val for s, val in dbe.symbols.items()} if parameters is not None: self._parameters_arg = parameters if isinstance(parameters, dict): symbols.update([(wrap_symbol(s), val) for s, val in parameters.items()]) else: # Lists of symbols that should remain symbolic for s in parameters: symbols.pop(wrap_symbol(s)) else: self._parameters_arg = None self._symbols = {wrap_symbol(key): value for key, value in symbols.items()} self.models = OrderedDict() self.build_phase(dbe) for name, value in self.models.items(): self.models[name] = self.symbol_replace(value, symbols) def __eq__(self, other): if self is other: return True elif type(self) != type(other): return False else: return self.__dict__ == other.__dict__ def __ne__(self, other): return not self.__eq__(other) def __hash__(self): return hash(repr(self)) # Default methods, don't need to be overridden unless you want to change the behavior @property def state_variables(self) -> List[v.StateVariable]: """Return a sorted list of state variables used in the ast which are not site fractions.""" return sorted((x for x in self.ast.free_symbols if not isinstance(x, v.SiteFraction) and isinstance(x, v.StateVariable)), key=str) @property def site_fractions(self) -> List[v.SiteFraction]: """Return a sorted list of site fractions used in the ast.""" return sorted((x for x in self.ast.free_symbols if isinstance(x, v.SiteFraction)), key=str) @property def ast(self): "Return the full abstract syntax tree of the model." return sum(self.models.values()) def _pair_test(self, constituent_array): """Return True if the constituent array represents a pair. Pairs have only one species in each sublattice. """ for subl in constituent_array: if len(subl) > 1: return False constituent = subl[0] if constituent not in self.components: return False return True def _array_validity(self, constituent_array): """Return True if the constituent array is satisfies all components.""" for subl in constituent_array: for constituent in subl: if constituent not in self.components: return False return True def _zeta_ik(self, i, k): pair_query = ( (where("phase_name") == self.phase_name) & \ (where("parameter_type") == "MQMG") & \ (where("constituent_array").test(lambda x: x == ((i,), (k,)))) ) params = list(self._dbe._parameters.search(pair_query)) assert len(params) == 1, f"Expected exactly one pair parameter for ({i, k}), got {len(params)}: {params}" param = params[0] return param["zeta"] def _X_ijkl(self, i, j, k, l) -> v.SiteFraction: """ Shorthand for creating a site fraction object v.Y for a quadruplet (ij/kl) """ return v.Y(self.phase_name, 0, _get_species(i, j, k, l)) def _n_ik(self, i, k): """ Return the amount of pair, n_i/k, for the pair i/k. Follows Poschmann Eq. 5 except that including zeta here (as by Poschmann) seems to break the tests when the coordination numbers vary across quadruplets. """ n_ik = S.Zero for a, b, x, y in self._quadruplets: n_ik += self._X_ijkl(a, b, x, y) * ((a == i) + (b == i)) * ((x == k) + (y == k)) return n_ik def _n_ik_star(self, i, k): """ Return the amount of pair, n^*_i/k, for the pair i/k following Lambotte JCT 2011 Eq. A.3 This term is only used in the pair part of the entropy. """ n_ik = S.Zero for a, b, x, y in self._quadruplets: n_ik += self._X_ijkl(a, b, x, y) * ((a == i) + (b == i)) * ((x == k) + (y == k)) n_ik /= self._zeta_ik(i, k) return n_ik def _X_ik(self, i: v.Species, k: v.Species): """ Return the endmember fraction, X_i/k, for a the pair i/k following Poschmann Eq. 6 """ X_ik = self._n_ik(i, k) / sum(self._n_ik(a, x) for a, x in self._pairs) return X_ik def _X_ik_star(self, i: v.Species, k: v.Species): """ Return the endmember fraction, X^*_i/k, for a the pair i/k. Follows Lambotte JCT 2011 Eq. A.5, but this is simply Poschmann Eq. 6 with n^*_ik instead of n_ik. This term is only used in the pair part of the entropy and as a sub-expression of F_i (also part of the pair part of the entropy). """ X_ik = self._n_ik_star(i, k) / sum(self._n_ik_star(a, x) for a, x in self._pairs) return X_ik def _n_i(self, dbe, species): """ Return the mass of the species following Poschmann Eq. 7 and 8. """ Z = partial(self.Z, dbe) # alias for notation n_i = S.Zero if species in self.cations: i = species for a, b, x, y in self._quadruplets: n_i += self._X_ijkl(a, b, x, y) * ((a == i) / Z(a, a, b, x, y) + (b == i) / Z(b, a, b, x, y)) else: assert species in self.anions k = species for a, b, x, y in self._quadruplets: n_i += self._X_ijkl(a, b, x, y) * ((x == k) / Z(x, a, b, x, y) + (y == k) / Z(y, a, b, x, y)) return n_i def _X_i(self, dbe, species: v.Species): """ Return the site fraction of species on it's sublattice. Poschmann Eq. 9 and 10. """ if species in self.cations: return self._n_i(dbe, species) / sum(self._n_i(dbe, a) for a in self.cations) else: assert species in self.anions return self._n_i(dbe, species) / sum(self._n_i(dbe, x) for x in self.anions) def _Y_i(self, species: v.Species): """ Return the site equivalent fraction of species following Poschmann Eq. 11 and 12. """ Y_i = S.Zero if species in self.cations: i = species for a, b, x, y in self._quadruplets: Y_i += self._X_ijkl(a, b, x, y) * ((a == i) + (b == i)) / 2 else: assert species in self.anions k = species for a, b, x, y in self._quadruplets: Y_i += self._X_ijkl(a, b, x, y) * ((x == k) + (y == k)) / 2 return Y_i def _F_i(self, species: v.Species): """ Return the coordination-equivalent fraction of species, Poschmann Eq. 13 and 14. """ soln_type = self._dbe.phases[self.phase_name].model_hints["mqmqa"]["type"] F_i = S.Zero if species in self.cations: i = species for a, x in self._pairs: if soln_type == "SUBG": F_i += int(a == i) * self._X_ik(a, x) elif soln_type == "SUBQ": F_i += int(a == i) * self._X_ik_star(a, x) else: assert species in self.anions k = species for a, x in self._pairs: if soln_type == "SUBG": F_i += int(x == k) * self._X_ik(a, x) elif soln_type == "SUBQ": F_i += int(x == k) * self._X_ik_star(a, x) return F_i def _chemical_group_filter(self, dbe, symmetric_species, asymmetric_species, sublattice): """ Return a function ``f(m)`` that returns ``True`` if m is symmetric with the symmetric_species and asymmetric with the asymmetric_species. """ # sublattice should be "cations" or "anions" chem_group_dict = dbe.phases[self.phase_name].model_hints["mqmqa"]["chemical_groups"][sublattice] def _f(species): if species == symmetric_species: return False elif species == asymmetric_species: return False elif chem_group_dict[species] == chem_group_dict[symmetric_species] and chem_group_dict[species] != chem_group_dict[asymmetric_species]: return True # This chemical group should be mixed else: return False return _f def _Chi_mix(self, dbe, i, j, k, l): """ (:math:`\\chi_{ij/k}`) following Poschmann Eq. 21 (SUBG-type model) or Eq. 22 (SUBQ-type models), respectively. """ cations = self.cations anions = self.anions X_ijkl = self._X_ijkl soln_type = dbe.phases[self.phase_name].model_hints["mqmqa"]["type"] mixing_term_numerator = S.Zero mixing_term_denominator = S.Zero if i == j and k == l: raise ValueError(f"Excess energies for pairs are not defined. Got quadruplet {(i, j, k, l)}") elif i != j and k == l: # Mixing on first sublattice # Chi_{ij/kk} type nu = list(filter(self._chemical_group_filter(dbe, i, j, "cations"), cations)) gamma = list(filter(self._chemical_group_filter(dbe, j, i, "cations"), cations)) for idx, a in enumerate([i] + nu): # enumerate to avoid double counting for b in ([i] + nu)[idx:]: if soln_type == "SUBG": mixing_term_numerator += X_ijkl(a, b, k, k) elif soln_type == "SUBQ": # Eq. 22 loop over anions for x in anions: for y in anions: mixing_term_numerator += X_ijkl(a, b, x, y) * ((k == x) + (k == y)) / 2 else: raise ValueError(f"Unknown solution type: {soln_type}") for idx, a in enumerate([i, j] + nu + gamma): # enumerate to avoid double counting for b in ([i, j] + nu + gamma)[idx:]: if soln_type == "SUBG": mixing_term_denominator += X_ijkl(a, b, k, k) elif soln_type == "SUBQ": # Eq. 22 loop over anions for x in anions: for y in anions: mixing_term_denominator += X_ijkl(a, b, x, y) * ((k == x) + (k == y)) / 2 else: raise ValueError(f"Unknown solution type: {soln_type}") return mixing_term_numerator / mixing_term_denominator elif i == j and k != l: # Mixing on second sublattice # Chi_{ii/kl} type nu = list(filter(self._chemical_group_filter(dbe, k, l, "anions"), anions)) gamma = list(filter(self._chemical_group_filter(dbe, l, k, "anions"), anions)) for idx, x in enumerate([k] + nu): # enumerate to avoid double counting for y in ([k] + nu)[idx:]: if soln_type == "SUBG": mixing_term_numerator += X_ijkl(i, i, x, y) elif soln_type == "SUBQ": # Eq. 22 loop over cations for a in cations: for b in cations: mixing_term_numerator += X_ijkl(a, b, x, y) * ((i == a) + (i == b)) / 2 else: raise ValueError(f"Unknown solution type: {soln_type}") for idx, x in enumerate([k, l] + nu + gamma): # enumerate to avoid double counting for y in ([k, l] + nu + gamma)[idx:]: if soln_type == "SUBG": mixing_term_denominator += X_ijkl(i, i, x, y) elif soln_type == "SUBQ": # Eq. 22 loop over cations for a in cations: for b in cations: mixing_term_denominator += X_ijkl(a, b, x, y) * ((i == a) + (i == b)) / 2 else: raise ValueError(f"Unknown solution type: {soln_type}") return mixing_term_numerator / mixing_term_denominator else: raise ValueError(f"Computing Chi_mix is not supported for reciprocal quadruplets. Got quadruplet {(i, j, k, l)}.") def _Y_ik(self, i, k): """ Poschmann Eq. 20 """ term = S.Zero for a, b, x, y in self._quadruplets: term += self._X_ijkl(a,b,x,y) * ((a == i) + (b == i)) * ((x == k) + (y == k)) / 4 return term def _Xi_mix(self, dbe, i, j, k, l): """ (:math:`\\xi_{ij/k}`) following Poschmann Eq. 19 """ # For mixing in cations (i != j, k == l), nu are cations (nu != i, nu != j) where i and nu have the same # chemical group and j has a different chemical group. mixing_term = S.Zero if i == j and k == l: raise ValueError(f"Computing Xi_mix is not supported for pair quadruplets (there must be mixing among cations or anions). Got quadruplet {(i, j, k, l)}.") elif i != j and k == l: # Mixing on first sublattice nu = list(filter(self._chemical_group_filter(dbe, i, j, "cations"), self.cations)) for a in [i] + nu: mixing_term += self._Y_ik(a, k) return mixing_term elif i == j and k != l: # Mixing on second sublattice nu = list(filter(self._chemical_group_filter(dbe, k, l, "anions"), self.anions)) for x in [k] + nu: mixing_term += self._Y_ik(i, x) return mixing_term else: raise ValueError(f"Computing Xi_mix is not supported for reciprocal quadruplets (there can only be mixing among cations _or_ anions). Got quadruplet {(i, j, k, l)}.") def _calc_Z(self, dbe, species, A, B, X, Y): # In derivations of the MQMQA, charges are written as if they have the same sign. The absolute values of the charges are used here. Z = partial(self.Z, dbe) if (species == A) or (species == B): species_is_cation = True elif (species == X) or (species == Y): species_is_cation = False else: raise ValueError(f"{species} is not A ({A}), B ({B}), X ({X}) or Y ({Y}).") if A == B and X == Y: raise ValueError(f"Z({species}, {A}{B}/{X}{Y}) is a pure pair and must be defined explictly") elif A != B and X != Y: # This is a reciprocal AB/XY quadruplet and needs to be calculated by eq 23 and 24 in Pelton et al. Met Trans B (2001) F = 1/8 * ( # eq. 24 abs(A.charge) / Z(A, A,A,X,Y) + abs(B.charge) / Z(B, B,B,X,Y) + abs(X.charge) / Z(X, A,B,X,X) + abs(Y.charge) / Z(Y, A,B,Y,Y) ) if species_is_cation: inv_Z = F * ( Z(X, A,B,X,X) / (abs(X.charge) * Z(species, A,B,X,X)) + Z(Y, A,B,Y,Y) / (abs(Y.charge) * Z(species, A,B,Y,Y)) ) else: inv_Z = F * ( Z(A, A,A,X,Y) / (abs(A.charge) * Z(species, A,A,X,Y)) + Z(B, B,B,X,Y) / (abs(B.charge) * Z(species, B,B,X,Y)) ) return 1 / inv_Z elif A != B: # X == Y # Need to calculate Z^i_AB/XX (Y = X). # We assume Z^A_ABXX = Z^A_AAXX = Z^A_AAYY # and Z^X_ABXX = (q_X + q_Y)/(q_A/Z^A_AAXX + q_B/Z^B_BBXX) # note: q_X = q_Y, etc. since Y = X # We don't know if these are correct, but that's what's implemented in Thermochimica if species_is_cation: return Z(species, species, species, X, X) else: return 2*abs(species.charge) / ( abs(A.charge) / Z(A, A,A,species,species) + abs(B.charge) / Z(B, B,B,species,species) ) elif X != Y: # A == B # These use the same equations as A != B case with the same assumptions if species_is_cation: # similarly, Z^A_AAXY = (q_A + q_B)/(q_X/Z^X_AAXX + q_Y/Z^Y_AAYY) return 2*abs(species.charge)/(abs(X.charge)/Z(X, species, species, X, X) + abs(Y.charge)/Z(Y, species, species, Y, Y)) else: return Z(species, A, A, species, species) raise ValueError("This should be unreachable")
[docs] def Z(self, dbe, species: v.Species, A: v.Species, B: v.Species, X: v.Species, Y: v.Species): # Canonicalize the order of cations and anions in alphabetical order # Note that the coordination number parameters (the constituents and # the values) _must_ obey this canonical sorted order. A, B = sorted((A, B)) X, Y = sorted((X, Y)) Zs = dbe._parameters.search( (where("phase_name") == self.phase_name) & \ (where("parameter_type") == "MQMZ") & \ (where("constituent_array").test(lambda x: x == ((A, B), (X, Y)))) ) if len(Zs) == 0: return self._calc_Z(dbe, species, A, B, X, Y) elif len(Zs) == 1: sp_idx = [A, B, X, Y].index(species) return Zs[0]["coordinations"][sp_idx] else: raise ValueError(f"Expected exactly one Z for {species} of {((A, B), (X, Y))}, got {len(Zs)}")
[docs] def get_internal_constraints(self): constraints = [] # Site fraction constraint site_frac_sum = S.Zero for a, b, x, y in self._quadruplets: site_frac_sum += self._X_ijkl(a, b, x, y) constraints.append(site_frac_sum - 1) # = 0 return constraints
@property def normalization(self): """ Total number of moles of this phase. Divide by this normalization factor to convert from J/mole-formula to J/mole-atoms """ return sum(self._n_i(self._dbe, i) for i in self.cations + self.anions if "VA" not in i.constituents)
[docs] def moles(self, species, per_formula_unit=False): "Number of moles of species or elements." species = v.Species(species) result = S.Zero for i in self.cations + self.anions: if list(species.constituents.keys())[0] in i.constituents: result += self._n_i(self._dbe, i) if per_formula_unit: return result else: return result / self.normalization
degree_of_ordering = DOO = S.Zero curie_temperature = TC = S.Zero beta = BMAG = S.Zero neel_temperature = NT = S.Zero # pylint: disable=C0103 # These are standard abbreviations from Thermo-Calc for these quantities GM = property(lambda self: self.ast / self.normalization) G = property(lambda self: self.ast) energy = GM entropy = SM = property(lambda self: -self.GM.diff(v.T)) enthalpy = HM = property(lambda self: self.GM - v.T * self.GM.diff(v.T)) heat_capacity = CPM = property(lambda self: -v.T * self.GM.diff(v.T, v.T)) # pylint: enable=C0103 mixing_energy = GM_MIX = property(lambda self: self.GM - self.reference_model.GM) mixing_enthalpy = HM_MIX = property(lambda self: self.GM_MIX - v.T * self.GM_MIX.diff(v.T)) mixing_entropy = SM_MIX = property(lambda self: -self.GM_MIX.diff(v.T)) mixing_heat_capacity = CPM_MIX = property(lambda self: -v.T * self.GM_MIX.diff(v.T, v.T))
[docs] def reference_energy(self, dbe): """ Returns energies contributed by pairs. """ pair_query = ( (where("phase_name") == self.phase_name) & \ (where("parameter_type") == "MQMG") & \ (where("constituent_array").test(self._pair_test)) ) params = dbe._parameters.search(pair_query) terms = S.Zero for param in params: a = i = param["constituent_array"][0][0] x = k = param["constituent_array"][1][0] X_ax = S.Zero for b in self.cations: for y in self.anions: X_ax += self._X_ijkl(a,b,x,y) * ((a == i) + (b == i)) * ((x == k) + (y == k)) / (2 * self.Z(dbe, a, a,b,x,y)) G_ax = param["parameter"] terms += X_ax * G_ax / param["stoichiometry"][0] return terms
[docs] def ideal_mixing_energy(self, dbe): # notational niceties n_i = partial(self._n_i, dbe) X_i = partial(self._X_i, dbe) X_ik = self._X_ik Y_i = self._Y_i F_i = self._F_i X_ijkl = self._X_ijkl soln_type = dbe.phases[self.phase_name].model_hints["mqmqa"]["type"] if soln_type == "SUBG": pow_X_ik = 1.0 pow_Y_i = 1.0 elif soln_type == "SUBQ": pow_X_ik = 0.75 pow_Y_i = 0.5 Sid = S.Zero # Individual constituents for A in self.cations: Sid += n_i(A) * log(X_i(A)) for X in self.anions: Sid += n_i(X) * log(X_i(X)) # Pairs if soln_type == "SUBG": for i, k in self._pairs: Sid += self._n_ik(i, k) * log(X_ik(i, k) / (F_i(i) * F_i(k))) elif soln_type == "SUBQ": for i, k in self._pairs: Sid += self._n_ik_star(i, k) * log(self._X_ik_star(i, k) / (F_i(i) * F_i(k))) # Quadruplets for i, j, k, l in self._quadruplets: C_ijkl = (2 - (i == j)) * (2 - (k == l)) Sid += X_ijkl(i, j, k, l) * log(X_ijkl(i, j, k, l) / (C_ijkl * (X_ik(i, k) * X_ik(i, l) * X_ik(j, k) * X_ik(j, l))**pow_X_ik / (Y_i(i) * Y_i(j) * Y_i(k) * Y_i(l))**pow_Y_i)) return Sid * v.T * v.R
[docs] def excess_mixing_energy(self, dbe): params = dbe._parameters.search( (where("phase_name") == self.phase_name) & (where("parameter_type") == "MQMX") & (where("constituent_array").test(self._array_validity)) ) cations = self.cations anions = self.anions X_ijkl = self._X_ijkl Z = partial(self.Z, dbe) energy = S.Zero for param in params: (A, B), (X, Y) = param["constituent_array"] exponents = param["exponents"] mixing_code = param["mixing_code"] m = param["additional_mixing_constituent"] p_alpha = exponents[0] q_alpha = exponents[1] r_alpha = param["additional_mixing_exponent"] # Poschmann Eq. 23-26 mixing_term = S.Zero if A != B and X == Y: if mixing_code == "G": # Poschmann Eq. 23 (cations mixing) mixing_term += self._Chi_mix(dbe, A, B, X, X)**p_alpha * self._Chi_mix(dbe, B, A, X, X)**q_alpha elif mixing_code == "Q": # Poschmann Eq. 19 and 20 (cations mixing) Xi_ijk = self._Xi_mix(dbe, A, B, X, X) Xi_jik = self._Xi_mix(dbe, B, A, X, X) # Poschmann Eq. 24 mixing_term += Xi_ijk**p_alpha * Xi_jik**q_alpha / (Xi_ijk + Xi_jik)**(p_alpha + q_alpha) else: raise ValueError(f"Unknown mixing code {mixing_code} for parameter {param}") if m != v.Species(None): # Poschmann Eq. 25 and 26 ternary term (same for both mixing codes) Xi_ijk = self._Xi_mix(dbe, A, B, X, X) Xi_jik = self._Xi_mix(dbe, B, A, X, X) Y_mk = self._Y_ik(m, X) nu = list(filter(self._chemical_group_filter(dbe, A, B, "cations"), cations)) gamma = list(filter(self._chemical_group_filter(dbe, B, A, "cations"), cations)) if m in gamma: mixing_term *= Y_mk / Xi_jik * (1 - self._Y_ik(B, X) / Xi_jik)**(r_alpha - 1) elif m in nu: mixing_term *= Y_mk / Xi_ijk * (1 - self._Y_ik(A, X) / Xi_ijk)**(r_alpha - 1) else: # not in nu or gamma mixing_term *= Y_mk * (1 - Xi_ijk - Xi_jik)**(r_alpha - 1) elif A == B and X != Y: if mixing_code == "G": # Poschmann Eq. 23 (anions mixing) mixing_term += self._Chi_mix(dbe, A, A, X, Y)**p_alpha * self._Chi_mix(dbe, A, A, Y, X)**q_alpha elif mixing_code == "Q": # Poschmann Eq. 19 and 20 (anions mixing) Xi_ikl = self._Xi_mix(dbe, A, A, X, Y) Xi_ilk = self._Xi_mix(dbe, A, A, Y, X) # Poschmann Eq. 24 mixing_term += Xi_ikl**p_alpha * Xi_ilk**q_alpha / (Xi_ikl + Xi_ilk)**(p_alpha + q_alpha) else: raise ValueError(f"Unknown mixing code {mixing_code} for parameter {param}") if m != v.Species(None): # Poschmann Eq. 25 and 26 ternary term (same for both mixing codes) Xi_ikl = self._Xi_mix(dbe, A, A, X, Y) Xi_ilk = self._Xi_mix(dbe, A, A, Y, X) Y_im = self._Y_ik(A, m) nu = list(filter(self._chemical_group_filter(dbe, X, Y, "anions"), anions)) gamma = list(filter(self._chemical_group_filter(dbe, Y, X, "anions"), anions)) if m in gamma: mixing_term *= Y_im / Xi_ilk * (1 - self._Y_ik(A, Y) / Xi_ilk)**(r_alpha - 1) elif m in nu: mixing_term *= Y_im / Xi_ikl * (1 - self._Y_ik(A, X) / Xi_ikl)**(r_alpha - 1) else: # not in nu or gamma mixing_term *= Y_im * (1 - Xi_ikl - Xi_ilk)**(r_alpha - 1) else: mixing_term = S.One # No mixing, this is a modification to the formation energy of this quadruplet g = param["parameter"] * mixing_term # Poschmann Eq. 17 cation_factor = S.Zero if A == B: for m in cations: if m != A: cation_factor += X_ijkl(A,m,X,Y) / Z(A, A,m,X,Y) cation_factor *= Z(A, A,B,X,Y) / 2 anion_factor = S.Zero if X == Y: for m in anions: if m != X: anion_factor += X_ijkl(A,B,X,m) / Z(X, A,B,X,m) anion_factor *= Z(X, A,B,X,Y) / 2 energy += 0.5 * g * (X_ijkl(A,B,X,Y) + cation_factor + anion_factor) return energy
[docs] def shift_reference_state(self, reference_states, dbe, contrib_mods=None, output=("GM", "HM", "SM", "CPM"), fmt_str="{}R"): raise NotImplementedError()
[docs] def build_phase(self, dbe): """ Generate the symbolic form of all the contributions to this phase. Parameters ---------- dbe : 'pycalphad.io.Database' """ self.models.clear() for key, value in self.__class__.contributions: self.models[key] = S(getattr(self, value)(dbe))
[docs] @staticmethod def symbol_replace(obj, symbols): """ Substitute values of symbols into 'obj'. Parameters ---------- obj : symengine object symbols : dict mapping symengine.Symbol to symengine object Returns ------- symengine object """ try: # Need to do more substitutions to catch symbols that are functions # of other symbols for iteration in range(_MAX_PARAM_NESTING): obj = obj.xreplace(symbols) undefs = [x for x in obj.free_symbols if not isinstance(x, v.StateVariable)] if len(undefs) == 0: break except AttributeError: # Can't use xreplace on a float pass return obj
def _build_reference_model(self, preserve_ideal=True): raise NotImplementedError() @property def reference_model(self): raise NotImplementedError("Endmember reference models do not have a physical meaning for MQMQA models")