ENH: Enable option B and option F support for TDBs (#412)

This PR closes #202.

pycalphad/core/utils.py

@@ -428,3 +428,102 @@ def wrap_symbol(obj):
     else:
         return Symbol(obj)
 
+
+def recursive_tuplify(x):
+    """Recursively convert a nested list to a tuple"""
+    def _tuplify(y):
+        if isinstance(y, list) or isinstance(y, tuple):
+            return tuple(_tuplify(i) if isinstance(i, (list, tuple)) else i for i in y)
+        else:
+            return y
+    return tuple(map(_tuplify, x))
+
+
+def canonical_sort_key(x):
+    """
+    Wrap strings in tuples so they'll sort.
+
+    Parameters
+    ----------
+    x : list
+        List of strings to sort
+
+    Returns
+    -------
+    tuple
+        tuple of strings that can be sorted
+    """
+    return [tuple(i) if isinstance(i, (tuple, list)) else (i,) for i in x]
+
+
+def generate_symmetric_group(configuration, symmetry):
+    """
+    For a particular configuration and list of sublattices that are symmetric,
+    generate all the symmetrically equivalent configurations.
+
+    Parameters
+    ----------
+    configuration : Sequence[Any]
+        Typically a constituent array. The length should correspond to the number of
+        sublattices in the phase.
+    symmetry : Union[None, Sequence[Sequence[int]]]
+        A list of lists giving the indices of symmetrically equivalent sublattices.
+        For example: a symmetry of `[[0, 1, 2, 3]]` means that the first four
+        sublattices are symmetric to each other. If multiple sublattices are given, the
+        sublattices are internally equivalent and the sublattices themselves are assumed
+        interchangeble. That is, for a symmetry of `[[0, 1], [2, 3]]`, sublattices
+        0 and 1 are equivalent to each other (i.e. `[0, 1] == [1, 0]`) and similarly for
+        sublattices 2 and 3. It also implies that the sublattices are interchangeable,
+        (i.e. `[[0, 1], [2, 3]] == [[2, 3], [0, 1]]`), but note that constituents cannot
+        change sublattices (i.e. `[[0, 1], [2, 3]] != [[0, 3], [2, 1]]`).
+        If `symmetry=None` is given, no new configurations are generated.
+
+    Returns
+    -------
+    tuple
+        Tuple of configuration tuples that are all symmetrically equivalent.
+
+    Notes
+    -----
+    In the general case, equivalency between sublattices, for example
+    (`[[0, 1], [2, 3]] == [[2, 3], [0, 1]]`), is not necessarily required. It
+    could be that sublattices 0 and 1 represent equivalent substitutional
+    sublattices, while 2 and 3 represent equivalent interstitial sites.
+    Interchanging sublattices between substitutional sublattices is allowed, but
+    the substitutional sites would not be interchangeable with the interstitial
+    sites. To achieve this kind of effect with this function, you would need to
+    call it once with the equivalent substitutional sublattices, then for each
+    generated configuration, call this function again, giving the unique
+    configurations for symmetric interstitial sublattices.
+    """
+    # recursively casting sequences to tuples ensures that the generated configurations are hashable
+    configuration = recursive_tuplify(configuration)
+    sublattice_indices = list(range(len(configuration)))
+    if symmetry is None:
+        return [configuration]
+    seen_subl_indices = sorted([i for equiv_subl in symmetry for i in equiv_subl])
+    # fixed_subl_indices were not given, they are assumed to be inequivalent and constant
+    fixed_subl_indices = sorted(set(sublattice_indices) - set(seen_subl_indices))
+
+    # permute within each sublattice, i.e. [0, 1] -> [[0, 1], [1, 0]]
+    intra_sublattice_permutations = (itertools.permutations(equiv_subl) for equiv_subl in symmetry)
+    # product, combining all internal sublattice permutations, i.e.
+    # [[0, 1], [1, 0]] and [[2, 3], [3, 2]] become [ ([0, 1], [2, 3]), ... ]
+    sublattice_products = itertools.product(*intra_sublattice_permutations)
+    # finally, swap sets of equivalent sublattices, i.e.
+    # [ ([0, 1], [2, 3]), ... ] -> [[ ([0, 1], [2, 3]),  ([2, 3], [0, 1]) ], ... ]
+    inter_sublattice_permutations = (itertools.permutations(x) for x in sublattice_products)
+
+    symmetrically_distinct_configurations = set()
+    # chain.from_iterable calls flatten out nested permutation lists, i.e.
+    # ([0, 1], [2, 3]) -> [0, 1, 2, 3]
+    for proposed_distinct_indices in itertools.chain.from_iterable(inter_sublattice_permutations):
+        new_config = list(configuration[i] for i in itertools.chain.from_iterable(proposed_distinct_indices))
+        # The configuration only contains indices for symmetric sublattices. For the
+        # inequivalent sublattices, we need to insert them at their proper indices.
+        # Indices _must_ be in sorted order because we are changing the array size on insertion.
+        for fixed_idx in fixed_subl_indices:
+            new_config.insert(fixed_idx, configuration[fixed_idx])
+        symmetrically_distinct_configurations.add(tuple(new_config))
+    return sorted(symmetrically_distinct_configurations, key=canonical_sort_key)
+

pycalphad/io/database.py

@@ -10,18 +10,14 @@
 import os
 from pycalphad.variables import Species
 from pycalphad.core.cache import fhash
+from pycalphad.core.utils import recursive_tuplify
 
 
 class DatabaseExportError(Exception):
     """Raised when a database cannot be written."""
     pass
 
 
-def _to_tuple(lst):
-    "Convert nested list to nested tuple. Source: Martijn Pieters on StackOverflow"
-    return tuple(_to_tuple(i) if isinstance(i, list) else i for i in lst)
-
-
 class Phase(object): #pylint: disable=R0903
     """
     Phase in the database.
@@ -54,7 +50,7 @@ def __repr__(self):
         return 'Phase({0!r})'.format(self.__dict__)
     def __hash__(self):
         return hash((self.name, self.constituents, tuple(self.sublattices),
-                     tuple(sorted(_to_tuple(self.model_hints.items())))))
+                     tuple(sorted(recursive_tuplify(self.model_hints.items())))))
 
 DatabaseFormat = namedtuple('DatabaseFormat', ['read', 'write'])
 format_registry = {}

pycalphad/io/tdb.py

@@ -11,12 +11,14 @@
 from symengine.lib.symengine_wrapper import UniversalSet, Union, Complement
 from symengine import sympify, And, Or, Not, EmptySet, Interval, Piecewise, Add, Mul, Pow
 from symengine import Symbol, LessThan, StrictLessThan, S, E
+from tinydb import where
 from pycalphad import Database
 from pycalphad.io.database import DatabaseExportError
 from pycalphad.io.grammar import float_number, chemical_formula
 from pycalphad.variables import Species
 import pycalphad.variables as v
 from pycalphad.io.tdb_keywords import expand_keyword, TDB_PARAM_TYPES
+from pycalphad.core.utils import generate_symmetric_group
 from collections import defaultdict, namedtuple
 import ast
 import sys
@@ -572,6 +574,43 @@ def _apply_new_symbol_names(dbf, symbol_name_map):
         dbf._parameters.update({'parameter': S(p['parameter']).xreplace({Symbol(s): Symbol(v) for s, v in symbol_name_map.items()})}, doc_ids=[p.doc_id])
 
 
+KNOWN_SUBLATTICE_SYMMETRY_RELATIONS = {
+    # Keys should correspond to the model hints added via the `phase_options` dict
+    "symmetry_FCC_4SL": [[0, 1, 2, 3]],
+    "symmetry_BCC_4SL": [[0, 1], [2, 3]],
+}
+
+
+def add_phase_symmetry_ordering_parameters(dbf):
+    for phase_name, phase_obj in dbf.phases.items():
+        if phase_obj.model_hints.get("ordered_phase", "") == phase_name:
+            for symmetry_hint, symmetry in KNOWN_SUBLATTICE_SYMMETRY_RELATIONS.items():
+                if phase_obj.model_hints.get(symmetry_hint, False):
+                    for param in dbf.search(where("phase_name") == phase_name):
+                        const_array = param["constituent_array"]
+                        for symm_unique_const_array in set(generate_symmetric_group(const_array, symmetry)) - {const_array}:
+                            new_param = {key: val for key, val in param.items()}
+                            new_param["constituent_array"] = symm_unique_const_array
+                            new_param["_generated_by_symmetry_option"] = True  # flag to be able to remove it later and preserve the phase option
+                            dbf._parameters.insert(new_param)
+
+
+def _symmetry_added_parameter(dbf, param):
+    """
+    Return true if parameter belongs to a phase with an active symmetry
+    option and the parameter was added by a symmetry option.
+    """
+    phase_obj = dbf.phases.get(param["phase_name"])
+    if phase_obj is None:
+        # Phase isn't in the database at all, so it's impossible for this parameter
+        # to get added by symmetry
+        return False
+    for symm_hint in set(KNOWN_SUBLATTICE_SYMMETRY_RELATIONS.keys()).intersection(phase_obj.model_hints.keys()):
+            if phase_obj.model_hints[symm_hint] and param.get("_generated_by_symmetry_option", False):
+                return True
+    return False
+
+
 def write_tdb(dbf, fd, groupby='subsystem', if_incompatible='warn'):
     """
     Write a TDB file from a pycalphad Database object.
@@ -626,6 +665,7 @@ def write_tdb(dbf, fd, groupby='subsystem', if_incompatible='warn'):
             _apply_new_symbol_names(dbf, symbol_name_map)
         elif if_incompatible == 'warn':
             warnings.warn('Ignoring that the following function names are beyond the 8 character TDB limit: {}. Use the keyword argument \'if_incompatible\' to control this behavior.'.format(long_function_names))
+
     # Begin constructing the written database
     writetime = datetime.datetime.now()
     maxlen = 78
@@ -742,6 +782,8 @@ def write_tdb(dbf, fd, groupby='subsystem', if_incompatible='warn'):
     paramtuple = namedtuple('ParamTuple', ['phase_name', 'parameter_type', 'complexity', 'constituent_array',
                                            'parameter_order', 'diffusing_species', 'parameter', 'reference'])
     for param in dbf._parameters.all():
+        if _symmetry_added_parameter(dbf, param):
+            continue  # skip this parameter
         if groupby == 'subsystem':
             components = set()
             for subl in param['constituent_array']:
@@ -887,6 +929,9 @@ def read_tdb(dbf, fd):
 
     dbf.process_parameter_queue()
 
+    # Add phase option B/F parameters
+    # Must occur after adding model hints and parameters
+    add_phase_symmetry_ordering_parameters(dbf)
 
 
 Database.register_format("tdb", read=read_tdb, write=write_tdb)