In [1]:
import pycalphad.io.tdb
from pycalphad.io.tdb import expand_keyword

def _process_typedef(targetdb, typechar, line):
    """
    Process the TYPE_DEFINITION command.
    """
    # GES A_P_D BCC_A2 MAGNETIC  -1    0.4
    tokens = line.replace(',', '').split()
    if len(tokens) < 4:
        return
    keyword = expand_keyword(['DISORDERED_PART', 'MAGNETIC', 'NEVER_DISORDERS'], tokens[3].upper())[0]
    if len(keyword) == 0:
        raise ValueError('Unknown keyword: {}'.format(tokens[3]))
    if keyword == 'MAGNETIC':
        # magnetic model (IHJ model assumed by default)
        targetdb.tdbtypedefs[typechar] = {
            'ihj_magnetic':[float(tokens[4]), float(tokens[5])]
        }
    # GES A_P_D L12_FCC DIS_PART FCC_A1
    if keyword == 'DISORDERED_PART':
        # order-disorder model
        targetdb.tdbtypedefs[typechar] = {
            'disordered_phase': tokens[4].upper(),
            'ordered_phase': tokens[2].upper()
        }
        if tokens[2].upper() in targetdb.phases:
            # Since TDB files do not enforce any kind of ordering
            # on the specification of ordered and disordered phases,
            # we need to handle the case of when either phase is specified
            # first. In this case, we imagine the ordered phase is
            # specified first. If the disordered phase is specified
            # first, we will have to catch it in _process_phase().
            targetdb.phases[tokens[2].upper()].model_hints.update(
                targetdb.tdbtypedefs[typechar]
            )
    # GES AMEND_PHASE_DESCRIPTION SIGMA NEVER_DIS SIGMA_DIS
    if keyword == 'NEVER_DISORDERS':
        if len(tokens) < 4:
            return
        targetdb.tdbtypedefs[typechar] = {'never_disorders': True}
        if tokens[2].upper() in targetdb.phases:
            targetdb.phases[tokens[2].upper()].model_hints.update(
                targetdb.tdbtypedefs[typechar]
            )

def _process_phase(targetdb, name, typedefs, subls):
    """
    Process the PHASE command.
    """
    splitname = name.split(':')
    phase_name = splitname[0].upper()
    options = None
    if len(splitname) > 1:
        options = splitname[1]
    targetdb.add_structure_entry(phase_name, phase_name)
    model_hints = {}
    for typedef in list(typedefs):
        if typedef in targetdb.tdbtypedefs.keys():
            if 'ihj_magnetic' in targetdb.tdbtypedefs[typedef].keys():
                model_hints['ihj_magnetic_afm_factor'] = \
                    targetdb.tdbtypedefs[typedef]['ihj_magnetic'][0]
                model_hints['ihj_magnetic_structure_factor'] = \
                    targetdb.tdbtypedefs[typedef]['ihj_magnetic'][1]
            if 'ordered_phase' in targetdb.tdbtypedefs[typedef].keys():
                model_hints['ordered_phase'] = \
                    targetdb.tdbtypedefs[typedef]['ordered_phase']
                model_hints['disordered_phase'] = \
                    targetdb.tdbtypedefs[typedef]['disordered_phase']
                if model_hints['disordered_phase'] in targetdb.phases:
                    targetdb.phases[model_hints['disordered_phase']]\
                        .model_hints.update({'ordered_phase': model_hints['ordered_phase'],
                                             'disordered_phase': model_hints['disordered_phase']})
            if 'never_disorders' in targetdb.tdbtypedefs[typedef].keys():
                model_hints['never_disorders'] = \
                    targetdb.tdbtypedefs[typedef]['never_disorders']
    targetdb.add_phase(phase_name, model_hints, subls)

pycalphad.io.tdb._TDB_PROCESSOR['TYPE_DEFINITION'] = _process_typedef
pycalphad.io.tdb._TDB_PROCESSOR['PHASE'] = _process_phase
from pycalphad import Model
from sympy import S, Add
import pycalphad.variables as v

class NeverDisordersModel(Model):
    def atomic_ordering_energy(self, dbe):
        """
        Return the atomic ordering contribution in symbolic form.
        Description follows Servant and Ansara, Calphad, 2001.
        
        Also includes "never disordering" model following Lukas, Fries and Sundman, 2007, p. 145.
        """
        phase = dbe.phases[self.phase_name]
        ordered_phase_name = phase.model_hints.get('ordered_phase', None)
        disordered_phase_name = phase.model_hints.get('disordered_phase', None)
        never_disorders = phase.model_hints.get('never_disorders', False)
        if phase.name != ordered_phase_name:
            return S.Zero
        disordered_model = self.__class__(dbe, sorted(self.components),
                                          disordered_phase_name)
        constituents = [sorted(set(c).intersection(self.components)) \
                for c in dbe.phases[ordered_phase_name].constituents]

        # Fix variable names
        variable_rename_dict = {}
        for atom in disordered_model.energy.atoms(v.SiteFraction):
            # Replace disordered phase site fractions with mole fractions of
            # ordered phase site fractions.
            # Special case: Pure vacancy sublattices
            all_species_in_sublattice = \
                dbe.phases[disordered_phase_name].constituents[
                    atom.sublattice_index]
            if atom.species == 'VA' and len(all_species_in_sublattice) == 1:
                # Assume: Pure vacancy sublattices are always last
                vacancy_subl_index = \
                    len(dbe.phases[ordered_phase_name].constituents)-1
                variable_rename_dict[atom] = \
                    v.SiteFraction(
                        ordered_phase_name, vacancy_subl_index, atom.species)
            else:
                # All other cases: replace site fraction with mole fraction
                variable_rename_dict[atom] = \
                    self.mole_fraction(
                        atom.species,
                        ordered_phase_name,
                        constituents,
                        dbe.phases[ordered_phase_name].sublattices
                        )
        # Save all of the ordered energy contributions
        # This step is why this routine must be called _last_ in build_phase
        ordered_energy = Add(*list(self.models.values()))
        self.models.clear()
        # Copy the disordered energy contributions into the correct bins
        for name, value in disordered_model.models.items():
            self.models[name] = value.xreplace(variable_rename_dict)
        # All magnetic parameters will be defined in the disordered model
        self.TC = self.curie_temperature = disordered_model.TC
        self.TC = self.curie_temperature = self.TC.xreplace(variable_rename_dict)

        molefraction_dict = {}

        # Construct a dictionary that replaces every site fraction with its
        # corresponding mole fraction in the disordered state
        for sitefrac in ordered_energy.atoms(v.SiteFraction):
            all_species_in_sublattice = \
                dbe.phases[ordered_phase_name].constituents[
                    sitefrac.sublattice_index]
            if sitefrac.species == 'VA' and len(all_species_in_sublattice) == 1:
                # pure-vacancy sublattices should not be replaced
                # this handles cases like AL,NI,VA:AL,NI,VA:VA and
                # ensures the VA's don't get mixed up
                continue
            molefraction_dict[sitefrac] = \
                self.mole_fraction(sitefrac.species,
                                   ordered_phase_name, constituents,
                                   dbe.phases[ordered_phase_name].sublattices)

        if never_disorders:
            # Remove disordered model ideal mixing contribution
            self.models['idmix'] = 0.0
            return ordered_energy
        else:
            return ordered_energy - ordered_energy.subs(molefraction_dict,
                                                        simultaneous=True)

In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
from pycalphad import Database, ternplot


dbf = Database('MoNiRe_BEF.TDB')

In [ ]:
fig = plt.figure(figsize=(9,6))

conds = {v.T: 500, v.P:101325, v.X('NI'): (0,1,0.015), v.X('RE'): (0,1,0.015)}
phases = sorted(dbf.phases.keys())

ternplot(dbf, ['MO', 'NI', 'RE', 'VA'], phases, conds,
         x=v.X('NI'), y=v.X('RE'), eq_kwargs={'model': NeverDisordersModel})
plt.show()

In [ ]:
sorted(dbf.phases.keys())

In [ ]: