From ad53e1c9233966430240818410e3ac05e17945f5 Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Tue, 21 Nov 2023 11:26:56 -0500
Subject: [PATCH 1/6] Fix CASM::xtal::make_primitive, which was not copying
 unique_names. This also fixes the output of libcasm.xtal.make_primitive which
 was losing the occ_dof list as a result.

---
 src/casm/crystallography/BasicStructureTools.cc | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/casm/crystallography/BasicStructureTools.cc b/src/casm/crystallography/BasicStructureTools.cc
index afac129..7e6c9cc 100644
--- a/src/casm/crystallography/BasicStructureTools.cc
+++ b/src/casm/crystallography/BasicStructureTools.cc
@@ -300,14 +300,19 @@ BasicStructure make_primitive(const BasicStructure &non_primitive_struc,
 
   // Fill up the basis
   BasicStructure primitive_struc(primitive_lattice);
+  std::vector<std::vector<std::string>> _unique_names;
+  Index i_site = 0;
   for (Site site_for_prim : non_primitive_struc.basis()) {
     site_for_prim.set_lattice(primitive_struc.lattice(), CART);
     if (find_index(primitive_struc.basis(), site_for_prim, tol) ==
         primitive_struc.basis().size()) {
       site_for_prim.within();
       primitive_struc.set_basis().emplace_back(std::move(site_for_prim));
+      _unique_names.push_back(non_primitive_struc.unique_names()[i_site]);
+      i_site++;
     }
   }
+  primitive_struc.set_unique_names(_unique_names);
 
   // TODO: Do we want this?
   primitive_struc.set_title(non_primitive_struc.title());

From 94d095ce0b9254cf3df33715f8e473e4b71a2558 Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Tue, 21 Nov 2023 11:30:31 -0500
Subject: [PATCH 2/6] Add to libcasm.xtal: make_primitive_prim (equivalent to
 current make_primitive), make_primtive_structure, and
 make_canonical_structure. Changed make_primitive to act on either Prim or
 Structure. Add  options to the BCC and FCC structure factory functions in
 libcasm.xtal.structures.

---
 python/libcasm/xtal/__init__.py   |  5 +-
 python/libcasm/xtal/_methods.py   | 52 +++++++++++++++---
 python/libcasm/xtal/structures.py | 47 +++++++++++++++--
 python/src/xtal.cpp               | 87 +++++++++++++++++++++++++++++--
 python/tests/test_prim.py         |  2 +
 python/tests/test_structure.py    | 51 ++++++++++++++++++
 python/tests/test_structures.py   | 35 +++++++++++++
 7 files changed, 262 insertions(+), 17 deletions(-)

diff --git a/python/libcasm/xtal/__init__.py b/python/libcasm/xtal/__init__.py
index fcbfbde..b3624a6 100644
--- a/python/libcasm/xtal/__init__.py
+++ b/python/libcasm/xtal/__init__.py
@@ -3,6 +3,7 @@
     make_canonical,
     make_crystal_point_group,
     make_factor_group,
+    make_primitive,
     make_within,
 )
 from ._xtal import (
@@ -29,12 +30,14 @@
     make_atom,
     make_canonical_lattice,
     make_canonical_prim,
+    make_canonical_structure,
     make_equivalent_property_values,
     make_point_group,
     make_prim_crystal_point_group,
     make_prim_factor_group,
     make_prim_within,
-    make_primitive,
+    make_primitive_prim,
+    make_primitive_structure,
     make_structure_crystal_point_group,
     make_structure_factor_group,
     make_structure_within,
diff --git a/python/libcasm/xtal/_methods.py b/python/libcasm/xtal/_methods.py
index e1bcc66..d37a853 100644
--- a/python/libcasm/xtal/_methods.py
+++ b/python/libcasm/xtal/_methods.py
@@ -3,27 +3,63 @@
 import libcasm.xtal._xtal as _xtal
 
 
+def make_primitive(
+    obj: Union[_xtal.Prim, _xtal.Structure],
+) -> Any:
+    """Make the primitive cell of a Prim or atomic Structure
+
+    Notes
+    -----
+    Currently, for Structure this method only considers atom coordinates and types.
+    Molecular coordinates and types are not considered. Properties are not considered.
+    The default CASM tolerance is used for comparisons. To consider molecules
+    or properties, or to use a different tolerance, use a Prim.
+
+    Parameters
+    ----------
+    obj: Union[ _xtal.Prim, _xtal.Structure]
+        A Prim or an atomic Structure, which determines whether
+        :func:`~libcasm.xtal.make_primitive_prim`, or
+        :func:`~libcasm.xtal.make_primitive_structure` is called.
+
+    Returns
+    -------
+    canonical_obj : Union[_xtal.Prim, _xtal.Structure]
+        The primitive equivalent Prim or atomic Structure.
+    """
+    if isinstance(obj, _xtal.Prim):
+        return _xtal.make_primitive_prim(obj)
+    elif isinstance(obj, _xtal.Structure):
+        return _xtal.make_primitive_structure(obj)
+    else:
+        raise TypeError(f"TypeError in make_primitive: received {type(obj).__name__}")
+
+
 def make_canonical(
-    obj: Union[_xtal.Lattice, _xtal.Prim],
+    obj: Union[_xtal.Lattice, _xtal.Prim, _xtal.Structure],
 ) -> Any:
-    """Make the canonical form of a Lattice or Prim
+    """Make an equivalent Lattice, Prim, or Structure with the canonical form
+    of the lattice
 
     Parameters
     ----------
-    obj: Union[_xtal.Lattice, _xtal.Prim]
-        A Lattice or Prim, which determines whether
-        :func:`~libcasm.xtal.make_canonical_lattice` or
-        :func:`~libcasm.xtal.make_canonical_prim` is called.
+    obj: Union[_xtal.Lattice, _xtal.Prim, _xtal.Structure]
+        A Lattice, Prim, or Structure, which determines whether
+        :func:`~libcasm.xtal.make_canonical_lattice`, or
+        :func:`~libcasm.xtal.make_canonical_prim`,
+        :func:`~libcasm.xtal.make_canonical_structure` is called.
 
     Returns
     -------
-    canonical_obj : Union[_xtal.Lattice, _xtal.Prim]
-        The canonical equivalent Lattice or Prim.
+    canonical_obj : Union[_xtal.Lattice, _xtal.Prim, _xtal.Structure]
+        The equivalent Lattice, Prim, or Structure with canonical form of the lattice.
     """
     if isinstance(obj, _xtal.Prim):
         return _xtal.make_canonical_prim(obj)
     elif isinstance(obj, _xtal.Lattice):
         return _xtal.make_canonical_lattice(obj)
+    elif isinstance(obj, _xtal.Structure):
+        return _xtal.make_canonical_structure(obj)
     else:
         raise TypeError(f"TypeError in make_canonical: received {type(obj).__name__}")
 
diff --git a/python/libcasm/xtal/structures.py b/python/libcasm/xtal/structures.py
index c5d8365..35ed9ee 100644
--- a/python/libcasm/xtal/structures.py
+++ b/python/libcasm/xtal/structures.py
@@ -12,8 +12,9 @@ def BCC(
     atom_type: str = "A",
     atom_properties: dict[str, np.ndarray] = {},
     global_properties: dict[str, np.ndarray] = {},
+    conventional=False,
 ) -> xtal.Structure:
-    r"""Construct a primitive BCC structure
+    r"""Construct a BCC structure
 
     Parameters
     ----------
@@ -32,19 +33,36 @@ def BCC(
         Continuous properties associated with entire crystal, if present. Keys must be
         the name of a CASM-supported property type. Values are (m, 1) arrays with
         dimensions matching the standard dimension of the property type.
-
+    conventional : bool = False
+        If True, construct the 2-atom conventional BCC cell instead of the 1-atom
+        primitive cell. Default is False.
     Returns
     -------
     structure : xtal.Structure
         A primitive BCC structure
     """
-    return xtal.Structure(
+    structure = xtal.Structure(
         lattice=xtal_lattices.BCC(r=r, a=a),
         atom_coordinate_frac=np.array([0.0, 0.0, 0.0]),
         atom_type=[atom_type],
         atom_properties=atom_properties,
         global_properties=global_properties,
     )
+    if conventional is True:
+        T_bcc_conventional = np.array(
+            [
+                [0, 1, 1],
+                [1, 0, 1],
+                [1, 1, 0],
+            ],
+            dtype=int,
+        )
+        return xtal.make_superstructure(
+            T_bcc_conventional,
+            structure,
+        )
+    else:
+        return structure
 
 
 def FCC(
@@ -53,8 +71,9 @@ def FCC(
     atom_type: str = "A",
     atom_properties: dict[str, np.ndarray] = {},
     global_properties: dict[str, np.ndarray] = {},
+    conventional: bool = False,
 ) -> xtal.Structure:
-    r"""Construct a primitive FCC structure
+    r"""Construct a FCC structure
 
     Parameters
     ----------
@@ -73,19 +92,37 @@ def FCC(
         Continuous properties associated with entire crystal, if present. Keys must be
         the name of a CASM-supported property type. Values are (m, 1) arrays with
         dimensions matching the standard dimension of the property type.
+    conventional : bool = False
+        If True, construct the 4-atom conventional FCC cell instead of the 1-atom
+        primitive cell. Default is False.
 
     Returns
     -------
     structure : xtal.Structure
         A primitive FCC structure
     """
-    return xtal.Structure(
+    structure = xtal.Structure(
         lattice=xtal_lattices.FCC(r=r, a=a),
         atom_coordinate_frac=np.array([0.0, 0.0, 0.0]),
         atom_type=[atom_type],
         atom_properties=atom_properties,
         global_properties=global_properties,
     )
+    if conventional is True:
+        T_fcc_conventional = np.array(
+            [
+                [-1, 1, 1],
+                [1, -1, 1],
+                [1, 1, -1],
+            ],
+            dtype=int,
+        )
+        return xtal.make_superstructure(
+            T_fcc_conventional,
+            structure,
+        )
+    else:
+        return structure
 
 
 def HCP(
diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp
index aadeff5..d02f7b6 100644
--- a/python/src/xtal.cpp
+++ b/python/src/xtal.cpp
@@ -573,7 +573,7 @@ std::shared_ptr<xtal::BasicStructure const> make_within(
   return prim;
 }
 
-std::shared_ptr<xtal::BasicStructure const> make_primitive(
+std::shared_ptr<xtal::BasicStructure const> make_primitive_prim(
     std::shared_ptr<xtal::BasicStructure const> const &init_prim) {
   auto prim = std::make_shared<xtal::BasicStructure>(*init_prim);
   *prim = xtal::make_primitive(*prim, prim->lattice().tol());
@@ -818,6 +818,37 @@ xtal::SimpleStructure make_simplestructure_within(
   return structure;
 }
 
+xtal::SimpleStructure make_primitive_simplestructure(
+    xtal::SimpleStructure const &init_structure) {
+  std::vector<std::vector<std::string>> occ_dof;
+  for (std::string name : init_structure.atom_info.names) {
+    occ_dof.push_back({name});
+  }
+  std::shared_ptr<xtal::BasicStructure const> prim = make_prim(
+      get_simplestructure_lattice(init_structure, TOL),
+      get_simplestructure_atom_coordinate_frac(init_structure), occ_dof);
+  prim = make_primitive_prim(prim);
+  xtal::SimpleStructure structure;
+
+  std::vector<std::string> atom_type;
+  for (auto const &site_names : prim->unique_names()) {
+    atom_type.push_back(site_names[0]);
+  }
+
+  return make_simplestructure(prim->lattice(), get_prim_coordinate_cart(prim),
+                              atom_type);
+}
+
+xtal::SimpleStructure make_canonical_simplestructure(
+    xtal::SimpleStructure const &init_structure) {
+  xtal::SimpleStructure structure = init_structure;
+  xtal::Lattice lattice = get_simplestructure_lattice(init_structure, TOL);
+  lattice.make_right_handed();
+  lattice = xtal::canonical::equivalent(lattice);
+  structure.lat_column_mat = lattice.lat_column_mat();
+  return structure;
+}
+
 xtal::SimpleStructure make_superstructure(
     py::EigenDRef<Eigen::Matrix3l const> transformation_matrix_to_super,
     xtal::SimpleStructure const &simple) {
@@ -2022,7 +2053,8 @@ PYBIND11_MODULE(_xtal, m) {
   m.def("make_within", &make_within, py::arg("init_prim"),
         "Equivalent to :func:`make_prim_within`");
 
-  m.def("make_primitive", &make_primitive, py::arg("init_prim"), R"pbdoc(
+  m.def("make_primitive_prim", &make_primitive_prim, py::arg("init_prim"),
+        R"pbdoc(
             Returns a primitive equivalent Prim
 
             A :class:`Prim` object is not forced to be the primitive equivalent
@@ -2039,7 +2071,7 @@ PYBIND11_MODULE(_xtal, m) {
 
             Returns
             -------
-            prim : Lattice
+            prim : Prim
                 The primitive equivalent prim.
             )pbdoc");
 
@@ -2761,6 +2793,55 @@ PYBIND11_MODULE(_xtal, m) {
   m.def("make_within", &make_simplestructure_within, py::arg("init_structure"),
         "Equivalent to :func:`make_structure_within`");
 
+  m.def("make_primitive_structure", &make_primitive_simplestructure,
+        py::arg("init_structure"), R"pbdoc(
+        Returns a primitive equivalent atomic Structure
+
+        This function finds and returns the primitive equivalent cell by checking for
+        internal translations that map all atoms onto equivalent atoms.
+
+        Notes
+        -----
+        Currently this method only considers atom coordinates and types. Molecular
+        coordinates and types are not considered. Properties are not considered.
+        The default CASM tolerance is used for comparisons. To consider molecules
+        or properties, or to use a different tolerance, use a Prim.
+
+        Parameters
+        ----------
+        init_structure: _xtal.Structure
+            The initial Structure
+
+        Returns
+        -------
+        structure: _xtal.Structure
+            The primitive equivalent Structure
+        )pbdoc");
+
+  m.def("make_canonical_structure", &make_canonical_simplestructure,
+        py::arg("init_structure"), R"pbdoc(
+        Returns an equivalent Structure with canonical lattice
+
+        Finds the canonical right-handed Niggli cell of the lattice, applying
+        lattice point group operations to find the equivalent lattice in a
+        standardized orientation. The canonical orientation prefers lattice
+        vectors that form symmetric matrices with large positive values on the
+        diagonal and small values off the diagonal. See also `Lattice Canonical Form`_.
+
+        .. _`Lattice Canonical Form`: https://prisms-center.github.io/CASMcode_docs/formats/lattice_canonical_form/
+
+        Parameters
+        ----------
+        init_structure: _xtal.Structure
+            The initial Structure
+
+        Returns
+        -------
+        structure: _xtal.Structure
+            The structure with canonical lattice.
+
+        )pbdoc");
+
   m.def("make_superstructure", &make_superstructure,
         py::arg("transformation_matrix_to_super").noconvert(),
         py::arg("structure"),
diff --git a/python/tests/test_prim.py b/python/tests/test_prim.py
index dcbe94e..b32b66c 100644
--- a/python/tests/test_prim.py
+++ b/python/tests/test_prim.py
@@ -165,12 +165,14 @@ def test_prim_to_poscar_str(shared_datadir):
 def test_make_primitive_occ(nonprimitive_cubic_occ_prim):
     assert nonprimitive_cubic_occ_prim.coordinate_frac().shape[1] == 2
     prim = xtal.make_primitive(nonprimitive_cubic_occ_prim)
+    assert isinstance(prim, xtal.Prim)
     assert prim.coordinate_frac().shape[1] == 1
 
 
 def test_make_primitive_manydof(test_nonprimitive_manydof_prim):
     assert test_nonprimitive_manydof_prim.coordinate_frac().shape[1] == 2
     prim = xtal.make_primitive(test_nonprimitive_manydof_prim)
+    assert isinstance(prim, xtal.Prim)
     assert prim.coordinate_frac().shape[1] == 1
 
 
diff --git a/python/tests/test_structure.py b/python/tests/test_structure.py
index 0d34442..bf46042 100644
--- a/python/tests/test_structure.py
+++ b/python/tests/test_structure.py
@@ -442,6 +442,57 @@ def test_make_superstructure_2():
     )
 
 
+def test_make_primitive_structure_1():
+    struc = xtal_structures.BCC(r=1)
+    transformation_matrix = np.array([[0, 1, 1], [2, 0, 1], [1, 1, 0]], dtype=int).T
+    assert transformation_matrix.flags.f_contiguous
+    assert transformation_matrix.dtype is np.dtype(np.int64)
+    superstruc = xtal.make_superstructure(transformation_matrix, struc)
+    assert np.allclose(
+        superstruc.lattice().column_vector_matrix(),
+        struc.lattice().column_vector_matrix() @ transformation_matrix,
+    )
+
+    primitive_struc = xtal.make_primitive_structure(superstruc)
+    assert primitive_struc.is_equivalent_to(struc)
+
+
+def test_make_primitive_structure_2():
+    struc = xtal_structures.BCC(r=1)
+    conventional_struc = xtal_structures.BCC(r=1, conventional=True)
+    primitive_struc = xtal.make_primitive_structure(conventional_struc)
+    assert primitive_struc.is_equivalent_to(struc)
+
+
+def test_make_primitive_structure_3():
+    struc = xtal_structures.FCC(r=1)
+    conventional_struc = xtal_structures.FCC(r=1, conventional=True)
+    primitive_struc = xtal.make_primitive_structure(conventional_struc)
+    assert primitive_struc.is_equivalent_to(struc)
+
+
+def test_make_canonical_structure_1():
+    struc = xtal.Structure(
+        lattice=xtal.Lattice(
+            np.array(
+                [
+                    [0.0, 0.0, 1.0],  # z
+                    [1.0, 0.0, 0.0],  # x
+                    [0.0, 1.0, 0.0],  # y
+                ]
+            ).transpose()
+        ),
+        atom_coordinate_frac=np.array(
+            [
+                [0.0, 0.0, 0.0],
+            ]
+        ).transpose(),
+        atom_type=["A"],
+    )
+    canonical_struc = xtal.make_canonical_structure(struc)
+    assert np.allclose(canonical_struc.lattice().column_vector_matrix(), np.eye(3))
+
+
 def test_make_superstructure_and_rotate():
     struc = xtal_structures.BCC(r=1)
     assert len(struc.atom_type()) == 1
diff --git a/python/tests/test_structures.py b/python/tests/test_structures.py
index 63f5a2f..29812ae 100644
--- a/python/tests/test_structures.py
+++ b/python/tests/test_structures.py
@@ -1,14 +1,49 @@
+import numpy as np
+
 import libcasm.xtal.structures as structures
 from libcasm.xtal import Structure
 
 
 def test_construct_all_structures():
+    ## BCC ##
     assert isinstance(structures.BCC(r=1.0), Structure)
     assert isinstance(structures.BCC(a=1.0), Structure)
 
+    structure = structures.BCC(a=1.0, conventional=True)
+    assert isinstance(structure, Structure)
+    assert np.allclose(structure.lattice().column_vector_matrix(), np.eye(3))
+    assert structure.atom_type() == ["A", "A"]
+    assert np.allclose(
+        structure.atom_coordinate_frac(),
+        np.array(
+            [
+                [0.0, 0.0, 0.0],
+                [0.5, 0.5, 0.5],
+            ]
+        ).transpose(),
+    )
+
+    ## FCC ##
     assert isinstance(structures.FCC(r=1.0), Structure)
     assert isinstance(structures.FCC(a=1.0), Structure)
 
+    structure = structures.FCC(a=1.0, conventional=True)
+    assert isinstance(structure, Structure)
+    assert np.allclose(structure.lattice().column_vector_matrix(), np.eye(3))
+    assert structure.atom_type() == ["A", "A", "A", "A"]
+    assert np.allclose(
+        structure.atom_coordinate_frac(),
+        np.array(
+            [
+                [0.0, 0.0, 0.0],
+                [0.5, 0.0, 0.5],
+                [0.0, 0.5, 0.5],
+                [0.5, 0.5, 0.0],
+            ]
+        ).transpose(),
+    )
+
+    ## HCP ##
     assert isinstance(structures.HCP(r=1.0), Structure)
     assert isinstance(structures.HCP(r=1.0, c=2.0 * 1.8), Structure)
     assert isinstance(structures.HCP(a=1.0), Structure)

From 6fe565aad65a135c1b4e529e2d5ab52741727e5f Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Thu, 30 Nov 2023 01:09:41 -0500
Subject: [PATCH 3/6] add excluded_species and frac parameters to
 Structure.to_dict; add sort_structure methods

---
 CHANGELOG.md                    |  14 ++
 python/libcasm/xtal/__init__.py |   5 +
 python/libcasm/xtal/_methods.py | 249 +++++++++++++++++++++++++++++++-
 python/src/xtal.cpp             |  31 +++-
 python/tests/conftest.py        |   2 -
 python/tests/test_structure.py  | 164 +++++++++++++++++++++
 6 files changed, 457 insertions(+), 8 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index dc97148..951e5db 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,20 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [2.0a9] - X
+
+### Changed
+
+- Changed make_primitive to act on either Prim or Structure.
+- Added `excluded_species` and `frac` parameters to xtal.Structure.to_dict
+
+### Added
+
+- Add to libcasm.xtal: make_primitive_prim (equivalent to current make_primitive), make_primtive_structure, and make_canonical_structure. 
+- Add options to the BCC and FCC structure factory functions in libcasm.xtal.structures to make the conventional cubic unit cells.
+- Added StructureAtomInfo namedtuple, and methods sort_structure_by_atom_info, sort_structure_by_atom_type, sort_structure_by_atom_coordinate_frac, and sort_structure_by_atom_coordinate_cart
+
+
 ## [2.0a8] - 2023-11-15
 
 ### Changed
diff --git a/python/libcasm/xtal/__init__.py b/python/libcasm/xtal/__init__.py
index b3624a6..5b51167 100644
--- a/python/libcasm/xtal/__init__.py
+++ b/python/libcasm/xtal/__init__.py
@@ -1,10 +1,15 @@
 """CASM Crystallography"""
 from ._methods import (
+    StructureAtomInfo,
     make_canonical,
     make_crystal_point_group,
     make_factor_group,
     make_primitive,
     make_within,
+    sort_structure_by_atom_coordinate_cart,
+    sort_structure_by_atom_coordinate_frac,
+    sort_structure_by_atom_info,
+    sort_structure_by_atom_type,
 )
 from ._xtal import (
     AtomComponent,
diff --git a/python/libcasm/xtal/_methods.py b/python/libcasm/xtal/_methods.py
index d37a853..451eabd 100644
--- a/python/libcasm/xtal/_methods.py
+++ b/python/libcasm/xtal/_methods.py
@@ -1,5 +1,11 @@
-from typing import Any, Union
+import functools
+import math
+from collections import namedtuple
+from typing import Any, Callable, Union
 
+import numpy as np
+
+import libcasm.casmglobal
 import libcasm.xtal._xtal as _xtal
 
 
@@ -157,3 +163,244 @@ def make_within(
         return _xtal.make_structure_within(obj)
     else:
         raise TypeError(f"TypeError in make_within: received {type(obj).__name__}")
+
+
+@functools.total_ordering
+class ApproximateFloatArray:
+    def __init__(
+        self,
+        arr: np.ndarray,
+        abs_tol: float = libcasm.casmglobal.TOL,
+    ):
+        """Store an array that will be compared lexicographically up to a given
+        absolute tolerance using math.isclose
+
+        Parameters
+        ----------
+        arr: np.ndarray
+            The array to be compared
+
+        abs_tol: float = libcasm.casmglobal.TOL
+            The absolute tolerance
+        """
+        if not isinstance(arr, np.ndarray):
+            raise Exception("Error in ApproximateFloatArray: arr must be a np.ndarray")
+        self.arr = arr
+        self.abs_tol = abs_tol
+
+    def __eq__(self, other):
+        if len(self.arr) != len(other.arr):
+            return False
+        for i in range(len(self.arr)):
+            if not math.isclose(self.arr[i], other.arr[i], abs_tol=self.abs_tol):
+                return False
+        return True
+
+    def __lt__(self, other):
+        if len(self.arr) != len(other.arr):
+            return len(self.arr) < len(other.arr)
+        for i in range(len(self.arr)):
+            if not math.isclose(self.arr[i], other.arr[i], abs_tol=self.abs_tol):
+                return self.arr[i] < other.arr[i]
+        return False
+
+
+StructureAtomInfo = namedtuple(
+    "StructureAtomInfo",
+    ["atom_type", "atom_coordinate_frac", "atom_coordinate_cart", "atom_properties"],
+)
+
+
+def sort_structure_by_atom_info(
+    structure: _xtal.Structure,
+    key: Callable[[StructureAtomInfo], Any],
+    reverse: bool = False,
+) -> _xtal.Structure:
+    """Sort an atomic structure
+
+    Parameters
+    ----------
+    structure: _xtal.Structure
+        The structure to be sorted. Must be an atomic structure only. Raises if
+        ``len(structure.mol_type()) != 0``.
+    key: Callable[[StructureAtomInfo], Any]
+        The function used to return a value which is sorted. This is passed to the
+        `key` parameter of `list.sort()` to sort a `list[StructureAtomInfo]`.
+    reverse: bool = False
+        By default, sort in ascending order. If ``reverse==True``, then sort in
+        descending order.
+
+    Returns
+    -------
+    sorted_structure: _xtal.Structure
+        An equivalent structure with atoms sorted as specified.
+    """
+
+    if len(structure.mol_type()) != 0:
+        raise Exception(
+            "Error: only atomic structures may be sorted using sort_by_atom_info"
+        )
+    atom_type = structure.atom_type()
+    atom_coordinate_frac = structure.atom_coordinate_frac()
+    atom_coordinate_cart = structure.atom_coordinate_cart()
+    atom_properties = structure.atom_properties()
+
+    atoms = []
+    import copy
+
+    for i in range(len(atom_type)):
+        atoms.append(
+            StructureAtomInfo(
+                copy.copy(atom_type[i]),
+                atom_coordinate_frac[:, i].copy(),
+                atom_coordinate_cart[:, i].copy(),
+                {key: atom_properties[key][:, i].copy() for key in atom_properties},
+            )
+        )
+
+    atoms.sort(key=key, reverse=reverse)
+
+    for i, atom in enumerate(atoms):
+        atom_type[i] = atom[0]
+        atom_coordinate_frac[:, i] = atom[1]
+        for key in atom_properties:
+            atom_properties[key][:, i] = atom[2][key]
+
+    sorted_struc = _xtal.Structure(
+        lattice=structure.lattice(),
+        atom_type=atom_type,
+        atom_coordinate_frac=atom_coordinate_frac,
+        atom_properties=atom_properties,
+        global_properties=structure.global_properties(),
+    )
+
+    return sorted_struc
+
+
+def sort_structure_by_atom_type(
+    structure: _xtal.Structure,
+    reverse: bool = False,
+) -> _xtal.Structure:
+    """Sort an atomic structure by atom type
+
+    Parameters
+    ----------
+    structure: _xtal.Structure
+        The structure to be sorted. Must be an atomic structure only. Raises if
+        ``len(structure.mol_type()) != 0``.
+    reverse: bool = False
+        By default, sort in ascending order. If ``reverse==True``, then sort in
+        descending order.
+
+    Returns
+    -------
+    sorted_structure: _xtal.Structure
+        An equivalent structure with atoms sorted by atom type.
+    """
+    return sort_structure_by_atom_info(
+        structure,
+        key=lambda atom_info: atom_info.atom_type,
+        reverse=reverse,
+    )
+
+
+def sort_structure_by_atom_coordinate_frac(
+    structure: _xtal.Structure,
+    order: str = "cba",
+    abs_tol: float = libcasm.casmglobal.TOL,
+    reverse: bool = False,
+) -> _xtal.Structure:
+    """Sort an atomic structure by fractional coordinates
+
+    Parameters
+    ----------
+    structure: _xtal.Structure
+        The structure to be sorted. Must be an atomic structure only. Raises if
+        ``len(structure.mol_type()) != 0``.
+    order: str = "cba"
+        Sort order of fractional coordinate components. Default "cba" sorts by
+        fractional coordinate along the "c" (third) lattice vector first, "b" (second)
+        lattice vector second, and "a" (first) lattice vector third.
+    abs_tol: float = libcasm.casmglobal.TOL
+        Floating point tolerance for coordinate comparisons.
+    reverse: bool = False
+        By default, sort in ascending order. If ``reverse==True``, then sort in
+        descending order.
+
+    Returns
+    -------
+    sorted_structure: _xtal.Structure
+        An equivalent structure with atoms sorted by fractional coordinates.
+    """
+
+    def compare_f(atom_info):
+        values = []
+        for i in range(len(order)):
+            if order[i] == "a":
+                values.append(atom_info.atom_coordinate_frac[0])
+            elif order[i] == "b":
+                values.append(atom_info.atom_coordinate_frac[1])
+            elif order[i] == "c":
+                values.append(atom_info.atom_coordinate_frac[2])
+
+        return ApproximateFloatArray(
+            arr=np.array(values),
+            abs_tol=abs_tol,
+        )
+
+    return sort_structure_by_atom_info(
+        structure,
+        key=compare_f,
+        reverse=reverse,
+    )
+
+
+def sort_structure_by_atom_coordinate_cart(
+    structure: _xtal.Structure,
+    order: str = "zyx",
+    abs_tol: float = libcasm.casmglobal.TOL,
+    reverse: bool = False,
+) -> _xtal.Structure:
+    """Sort an atomic structure by Cartesian coordinates
+
+    Parameters
+    ----------
+    structure: _xtal.Structure
+        The structure to be sorted. Must be an atomic structure only. Raises if
+        ``len(structure.mol_type()) != 0``.
+    order: str = "zyx"
+        Sort order of Cartesian coordinate components. Default "zyx" sorts by
+        "z" Cartesian coordinate first, "y" Cartesian coordinate second, and "x"
+        Cartesian coordinate third.
+    abs_tol: float = libcasm.casmglobal.TOL
+        Floating point tolerance for coordinate comparisons.
+    reverse: bool = False
+        By default, sort in ascending order. If ``reverse==True``, then sort in
+        descending order.
+
+    Returns
+    -------
+    sorted_structure: _xtal.Structure
+        An equivalent structure with atoms sorted by Cartesian coordinates.
+    """
+
+    def compare_f(atom_info):
+        values = []
+        for i in range(len(order)):
+            if order[i] == "x":
+                values.append(atom_info.atom_coordinate_frac[0])
+            elif order[i] == "y":
+                values.append(atom_info.atom_coordinate_frac[1])
+            elif order[i] == "z":
+                values.append(atom_info.atom_coordinate_frac[2])
+
+        return ApproximateFloatArray(
+            arr=np.array(values),
+            abs_tol=abs_tol,
+        )
+
+    return sort_structure_by_atom_info(
+        structure,
+        key=compare_f,
+        reverse=reverse,
+    )
diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp
index d02f7b6..1dbe25d 100644
--- a/python/src/xtal.cpp
+++ b/python/src/xtal.cpp
@@ -2543,14 +2543,35 @@ PYBIND11_MODULE(_xtal, m) {
           py::arg("data"))
       .def(
           "to_dict",
-          [](xtal::SimpleStructure const &simple) {
+          [](xtal::SimpleStructure const &simple,
+             std::vector<std::string> const &excluded_species, bool frac) {
             jsonParser json;
-            to_json(simple, json);
+            COORD_TYPE mode = frac ? FRAC : CART;
+            std::set<std::string> _excluded_species(excluded_species.begin(),
+                                                    excluded_species.end());
+            to_json(simple, json, _excluded_species, mode);
             return static_cast<nlohmann::json>(json);
           },
-          "Represent the Structure as a Python dict. The `Structure reference "
-          "<https://prisms-center.github.io/CASMcode_docs/formats/casm/"
-          "crystallography/SimpleStructure/>`_ documents the format.")
+          py::arg("excluded_species") =
+              std::vector<std::string>({"Va", "VA", "va"}),
+          py::arg("frac") = false, R"pbdoc(
+          Represent the Structure as a Python dict.
+
+          Parameters
+          ----------
+          excluded_species : list[str] = ["Va", "VA", "va"]
+              The names of any molecular or atomic species that should not be included
+              in the output.
+          frac : boolean, default=False
+                If True, write coordinates using fractional coordinates
+                relative to the lattice vectors. If False, write coordinates using
+                Cartesian coordinates.
+
+          Returns
+          -------
+          data : json
+              The `Structure reference <https://prisms-center.github.io/CASMcode_docs/formats/casm/crystallography/SimpleStructure/>`_ documents the format.
+          )pbdoc")
       .def_static("from_json", &simplestructure_from_json, R"pbdoc(
           Construct a Structure from a JSON-formatted string.
 
diff --git a/python/tests/conftest.py b/python/tests/conftest.py
index db7d47e..e32ea13 100644
--- a/python/tests/conftest.py
+++ b/python/tests/conftest.py
@@ -407,7 +407,6 @@ def example_structure_1():
         ]
     ).transpose()
     atom_properties = {"disp": atom_disp}
-    print(atom_properties)
 
     # global properties
     np.array(
@@ -421,7 +420,6 @@ def example_structure_1():
     # Hstrain_vector = converter.from_F(F)
     Hstrain_vector = np.array([0.009950330853168087, 0.0, 0.0, 0.0, 0.0, 0.0])
     global_properties = {"Hstrain": Hstrain_vector}
-    print(global_properties)
 
     return xtal.Structure(
         lattice=lattice,
diff --git a/python/tests/test_structure.py b/python/tests/test_structure.py
index bf46042..eadeb17 100644
--- a/python/tests/test_structure.py
+++ b/python/tests/test_structure.py
@@ -69,15 +69,31 @@ def test_structure_to_dict(example_structure_1):
 
     assert "lattice_vectors" in data
     assert "coordinate_mode" in data
+    assert len(data["atom_type"]) == 4
     assert "atom_coords" in data
     assert "atom_properties" in data
     assert len(data["atom_properties"]) == 1
     assert "disp" in data["atom_properties"]
     assert len(data["global_properties"]) == 1
     assert "Hstrain" in data["global_properties"]
+    expected = np.array(
+        [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5], [0.0, 0.0, 1.0], [0.5, 0.5, 1.5]]
+    )
+    print(xtal.pretty_json(data["atom_coords"]))
+    assert np.allclose(np.array(data["atom_coords"]), expected)
 
     assert isinstance(xtal.pretty_json(data), str)
 
+    data = structure.to_dict(excluded_species=["B"])
+    assert data["atom_type"] == ["A", "A"]
+
+    data = structure.to_dict(frac=True)
+    expected = np.array(
+        [[0.0, 0.0, 0.0], [0.5, 0.5, 0.25], [0.0, 0.0, 0.5], [0.5, 0.5, 0.75]]
+    )
+    print(xtal.pretty_json(data["atom_coords"]))
+    assert np.allclose(np.array(data["atom_coords"]), expected)
+
 
 def test_structure_from_dict():
     data = {
@@ -531,3 +547,151 @@ def test_make_superstructure_and_rotate():
     S = rotated_superstruc.lattice().column_vector_matrix()
     assert np.allclose(S, rotation_matrix @ L @ transformation_matrix)
     assert len(rotated_superstruc.atom_type()) == 2
+
+
+def test_structure_sort_structure_by_atom_type():
+    struc = xtal_structures.BCC(a=1.0, conventional=True)
+    superstruc = xtal.make_superstructure(np.eye(3, dtype=int) * 3, struc)
+    unsorted_atom_type = ["A", "B"] * 27
+    unsorted_struc = xtal.Structure(
+        lattice=superstruc.lattice(),
+        atom_type=unsorted_atom_type,
+        atom_coordinate_frac=superstruc.atom_coordinate_frac(),
+    )
+    assert unsorted_struc.atom_type() == unsorted_atom_type
+    # print(xtal.pretty_json(unsorted_struc.to_dict()))
+
+    sorted_struc = xtal.sort_structure_by_atom_type(unsorted_struc)
+    # print(xtal.pretty_json(sorted_struc.to_dict()))
+    assert sorted_struc.atom_type() == ["A"] * 27 + ["B"] * 27
+
+
+def test_structure_sort_structure_by_atom_coordinate_frac():
+    struc = xtal_structures.BCC(r=1.0)
+    unsorted_struc = xtal.make_superstructure(np.eye(3, dtype=int) * 2, struc)
+
+    sorted_struc = xtal.sort_structure_by_atom_coordinate_frac(
+        unsorted_struc,
+        order="cba",
+    )
+    expected = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [0.5, 0.0, 0.0],
+            [0.0, 0.5, 0.0],
+            [0.5, 0.5, 0.0],
+            [0.0, 0.0, 0.5],
+            [0.5, 0.0, 0.5],
+            [0.0, 0.5, 0.5],
+            [0.5, 0.5, 0.5],
+        ]
+    ).transpose()
+    # print(xtal.pretty_json(unsorted_struc.to_dict(frac=True)))
+    assert np.allclose(sorted_struc.atom_coordinate_frac(), expected)
+
+    sorted_struc = xtal.sort_structure_by_atom_coordinate_frac(
+        unsorted_struc,
+        order="abc",
+    )
+    expected = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [0.0, 0.0, 0.5],
+            [0.0, 0.5, 0.0],
+            [0.0, 0.5, 0.5],
+            [0.5, 0.0, 0.0],
+            [0.5, 0.0, 0.5],
+            [0.5, 0.5, 0.0],
+            [0.5, 0.5, 0.5],
+        ]
+    ).transpose()
+    # print(xtal.pretty_json(sorted_struc.to_dict(frac=True)))
+    assert np.allclose(sorted_struc.atom_coordinate_frac(), expected)
+
+
+def test_structure_sort_structure_by_atom_coordinate_cart():
+    struc = xtal_structures.BCC(a=1.0, conventional=True)
+    unsorted_struc = xtal.make_superstructure(np.eye(3, dtype=int) * 2, struc)
+    # print(xtal.pretty_json(unsorted_struc.to_dict()))
+
+    sorted_struc = xtal.sort_structure_by_atom_coordinate_cart(
+        unsorted_struc,
+        order="zyx",
+    )
+    expected = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [1.0, 1.0, 0.0],
+            [0.5, 0.5, 0.5],
+            [1.5, 0.5, 0.5],
+            [0.5, 1.5, 0.5],
+            [1.5, 1.5, 0.5],
+            [0.0, 0.0, 1.0],
+            [1.0, 0.0, 1.0],
+            [0.0, 1.0, 1.0],
+            [1.0, 1.0, 1.0],
+            [0.5, 0.5, 1.5],
+            [1.5, 0.5, 1.5],
+            [0.5, 1.5, 1.5],
+            [1.5, 1.5, 1.5],
+        ]
+    ).transpose()
+    # print(xtal.pretty_json(sorted_struc.to_dict()))
+    assert np.allclose(sorted_struc.atom_coordinate_cart(), expected)
+
+    sorted_struc = xtal.sort_structure_by_atom_coordinate_cart(
+        unsorted_struc,
+        order="zyx",
+        reverse=True,
+    )
+    expected = np.array(
+        [
+            [1.5, 1.5, 1.5],
+            [0.5, 1.5, 1.5],
+            [1.5, 0.5, 1.5],
+            [0.5, 0.5, 1.5],
+            [1.0, 1.0, 1.0],
+            [0.0, 1.0, 1.0],
+            [1.0, 0.0, 1.0],
+            [0.0, 0.0, 1.0],
+            [1.5, 1.5, 0.5],
+            [0.5, 1.5, 0.5],
+            [1.5, 0.5, 0.5],
+            [0.5, 0.5, 0.5],
+            [1.0, 1.0, 0.0],
+            [0.0, 1.0, 0.0],
+            [1.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0],
+        ]
+    ).transpose()
+    # print(xtal.pretty_json(sorted_struc.to_dict()))
+    assert np.allclose(sorted_struc.atom_coordinate_cart(), expected)
+
+    sorted_struc = xtal.sort_structure_by_atom_coordinate_cart(
+        unsorted_struc,
+        order="xyz",
+    )
+    expected = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0],
+            [0.0, 1.0, 0.0],
+            [0.0, 1.0, 1.0],
+            [0.5, 0.5, 0.5],
+            [0.5, 0.5, 1.5],
+            [0.5, 1.5, 0.5],
+            [0.5, 1.5, 1.5],
+            [1.0, 0.0, 0.0],
+            [1.0, 0.0, 1.0],
+            [1.0, 1.0, 0.0],
+            [1.0, 1.0, 1.0],
+            [1.5, 0.5, 0.5],
+            [1.5, 0.5, 1.5],
+            [1.5, 1.5, 0.5],
+            [1.5, 1.5, 1.5],
+        ]
+    ).transpose()
+    # print(xtal.pretty_json(sorted_struc.to_dict()))
+    assert np.allclose(sorted_struc.atom_coordinate_cart(), expected)

From a2f4879c888db3c3f96c4de721e190d301a115a9 Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Thu, 30 Nov 2023 10:06:00 -0500
Subject: [PATCH 4/6] In xtal::make_primitive, do not copy
 xtal::BasicStructure::unique_names if they do not exist

---
 src/casm/crystallography/BasicStructureTools.cc | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/casm/crystallography/BasicStructureTools.cc b/src/casm/crystallography/BasicStructureTools.cc
index 7e6c9cc..3bc6d1a 100644
--- a/src/casm/crystallography/BasicStructureTools.cc
+++ b/src/casm/crystallography/BasicStructureTools.cc
@@ -308,9 +308,11 @@ BasicStructure make_primitive(const BasicStructure &non_primitive_struc,
         primitive_struc.basis().size()) {
       site_for_prim.within();
       primitive_struc.set_basis().emplace_back(std::move(site_for_prim));
-      _unique_names.push_back(non_primitive_struc.unique_names()[i_site]);
-      i_site++;
+      if (!non_primitive_struc.unique_names().empty()) {
+        _unique_names.push_back(non_primitive_struc.unique_names()[i_site]);
+      }
     }
+    i_site++;
   }
   primitive_struc.set_unique_names(_unique_names);
 

From e136659691cf4df1d286d6feaa8ee3652e21952b Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Thu, 30 Nov 2023 11:55:31 -0500
Subject: [PATCH 5/6] Add substitute_structure_species and update docs

---
 CHANGELOG.md                    |  5 +-
 python/libcasm/xtal/__init__.py |  1 +
 python/libcasm/xtal/_methods.py | 98 ++++++++++++++++++++++++++++-----
 python/src/xtal.cpp             | 17 +++---
 python/tests/test_structure.py  | 24 +++++++-
 5 files changed, 117 insertions(+), 28 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 951e5db..e4bffe3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,14 +11,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 
 - Changed make_primitive to act on either Prim or Structure.
+- Changed default of `xtal.Structure.to_dict` to output in fractional coordinates
 - Added `excluded_species` and `frac` parameters to xtal.Structure.to_dict
 
 ### Added
 
 - Add to libcasm.xtal: make_primitive_prim (equivalent to current make_primitive), make_primtive_structure, and make_canonical_structure. 
 - Add options to the BCC and FCC structure factory functions in libcasm.xtal.structures to make the conventional cubic unit cells.
-- Added StructureAtomInfo namedtuple, and methods sort_structure_by_atom_info, sort_structure_by_atom_type, sort_structure_by_atom_coordinate_frac, and sort_structure_by_atom_coordinate_cart
-
+- Add to libcasm.xtal: StructureAtomInfo namedtuple, and methods sort_structure_by_atom_info, sort_structure_by_atom_type, sort_structure_by_atom_coordinate_frac, and sort_structure_by_atom_coordinate_cart
+- Add to libcasm.xtal: substitute_structure_species 
 
 ## [2.0a8] - 2023-11-15
 
diff --git a/python/libcasm/xtal/__init__.py b/python/libcasm/xtal/__init__.py
index 5b51167..9c0c8c0 100644
--- a/python/libcasm/xtal/__init__.py
+++ b/python/libcasm/xtal/__init__.py
@@ -10,6 +10,7 @@
     sort_structure_by_atom_coordinate_frac,
     sort_structure_by_atom_info,
     sort_structure_by_atom_type,
+    substitute_structure_species,
 )
 from ._xtal import (
     AtomComponent,
diff --git a/python/libcasm/xtal/_methods.py b/python/libcasm/xtal/_methods.py
index 451eabd..03fbf60 100644
--- a/python/libcasm/xtal/_methods.py
+++ b/python/libcasm/xtal/_methods.py
@@ -177,14 +177,16 @@ def __init__(
 
         Parameters
         ----------
-        arr: np.ndarray
+        arr: numpy.ndarray
             The array to be compared
 
-        abs_tol: float = libcasm.casmglobal.TOL
+        abs_tol: float = :data:`~libcasm.casmglobal.TOL`
             The absolute tolerance
         """
         if not isinstance(arr, np.ndarray):
-            raise Exception("Error in ApproximateFloatArray: arr must be a np.ndarray")
+            raise TypeError(
+                "Error in ApproximateFloatArray: arr must be a numpy.ndarray"
+            )
         self.arr = arr
         self.abs_tol = abs_tol
 
@@ -209,6 +211,25 @@ def __lt__(self, other):
     "StructureAtomInfo",
     ["atom_type", "atom_coordinate_frac", "atom_coordinate_cart", "atom_properties"],
 )
+""" A namedtuple, used to hold atom info when sorting atoms in a 
+:class:`~_xtal.Structure`.
+
+.. rubric:: Constructor
+
+Parameters
+----------
+atom_type: str
+    The atom type, from :func:`~_xtal.Structure.atom_type`.
+atom_coordinate_frac: numpy.ndarray[numpy.float64[3]]
+    The fractional coordinate of the atom, from 
+    :func:`~_xtal.Structure.atom_type.atom_coordinate_frac`.
+atom_coordinate_cart: numpy.ndarray[numpy.float64[3]]
+    The Cartesian coordinate of the atom, from 
+    :func:`~_xtal.Structure.atom_type.atom_coordinate_cart`.
+atom_properties: dict[str, numpy.ndarray[numpy.float64[m]]]
+    The continuous properties associated with the atoms, if present, from 
+    :func:`~_xtal.Structure.atom_type.atom_coordinate_frac`.
+"""
 
 
 def sort_structure_by_atom_info(
@@ -221,8 +242,7 @@ def sort_structure_by_atom_info(
     Parameters
     ----------
     structure: _xtal.Structure
-        The structure to be sorted. Must be an atomic structure only. Raises if
-        ``len(structure.mol_type()) != 0``.
+        The structure to be sorted. Must be an atomic structure only.
     key: Callable[[StructureAtomInfo], Any]
         The function used to return a value which is sorted. This is passed to the
         `key` parameter of `list.sort()` to sort a `list[StructureAtomInfo]`.
@@ -234,10 +254,15 @@ def sort_structure_by_atom_info(
     -------
     sorted_structure: _xtal.Structure
         An equivalent structure with atoms sorted as specified.
+
+    Raises
+    ------
+    ValueError
+        For non-atomic structure, if ``len(structure.mol_type()) != 0``.
     """
 
     if len(structure.mol_type()) != 0:
-        raise Exception(
+        raise ValueError(
             "Error: only atomic structures may be sorted using sort_by_atom_info"
         )
     atom_type = structure.atom_type()
@@ -286,8 +311,7 @@ def sort_structure_by_atom_type(
     Parameters
     ----------
     structure: _xtal.Structure
-        The structure to be sorted. Must be an atomic structure only. Raises if
-        ``len(structure.mol_type()) != 0``.
+        The structure to be sorted. Must be an atomic structure only.
     reverse: bool = False
         By default, sort in ascending order. If ``reverse==True``, then sort in
         descending order.
@@ -296,6 +320,11 @@ def sort_structure_by_atom_type(
     -------
     sorted_structure: _xtal.Structure
         An equivalent structure with atoms sorted by atom type.
+
+    Raises
+    ------
+    ValueError
+        For non-atomic structure, if ``len(structure.mol_type()) != 0``.
     """
     return sort_structure_by_atom_info(
         structure,
@@ -315,13 +344,12 @@ def sort_structure_by_atom_coordinate_frac(
     Parameters
     ----------
     structure: _xtal.Structure
-        The structure to be sorted. Must be an atomic structure only. Raises if
-        ``len(structure.mol_type()) != 0``.
+        The structure to be sorted. Must be an atomic structure only.
     order: str = "cba"
         Sort order of fractional coordinate components. Default "cba" sorts by
         fractional coordinate along the "c" (third) lattice vector first, "b" (second)
         lattice vector second, and "a" (first) lattice vector third.
-    abs_tol: float = libcasm.casmglobal.TOL
+    abs_tol: float = :data:`~libcasm.casmglobal.TOL`
         Floating point tolerance for coordinate comparisons.
     reverse: bool = False
         By default, sort in ascending order. If ``reverse==True``, then sort in
@@ -331,6 +359,11 @@ def sort_structure_by_atom_coordinate_frac(
     -------
     sorted_structure: _xtal.Structure
         An equivalent structure with atoms sorted by fractional coordinates.
+
+    Raises
+    ------
+    ValueError
+        For non-atomic structure, if ``len(structure.mol_type()) != 0``.
     """
 
     def compare_f(atom_info):
@@ -366,13 +399,12 @@ def sort_structure_by_atom_coordinate_cart(
     Parameters
     ----------
     structure: _xtal.Structure
-        The structure to be sorted. Must be an atomic structure only. Raises if
-        ``len(structure.mol_type()) != 0``.
+        The structure to be sorted. Must be an atomic structure only.
     order: str = "zyx"
         Sort order of Cartesian coordinate components. Default "zyx" sorts by
         "z" Cartesian coordinate first, "y" Cartesian coordinate second, and "x"
         Cartesian coordinate third.
-    abs_tol: float = libcasm.casmglobal.TOL
+    abs_tol: float = :data:`~libcasm.casmglobal.TOL`
         Floating point tolerance for coordinate comparisons.
     reverse: bool = False
         By default, sort in ascending order. If ``reverse==True``, then sort in
@@ -382,6 +414,11 @@ def sort_structure_by_atom_coordinate_cart(
     -------
     sorted_structure: _xtal.Structure
         An equivalent structure with atoms sorted by Cartesian coordinates.
+
+    Raises
+    ------
+    ValueError
+        For non-atomic structure, if ``len(structure.mol_type()) != 0``.
     """
 
     def compare_f(atom_info):
@@ -404,3 +441,36 @@ def compare_f(atom_info):
         key=compare_f,
         reverse=reverse,
     )
+
+
+def substitute_structure_species(
+    structure: _xtal.Structure,
+    substitutions: dict[str, str],
+) -> _xtal.Structure:
+    """Create a copy of a structure with renamed atomic and molecular species
+
+    Parameters
+    ----------
+    structure: _xtal.Structure
+        The initial structure
+    substitutions: dict[str, str]
+        The substitutions to make, using the convention key->value. For example, using
+        ``substitutions = { "B": "C"}`` results in all `atom_type` and `mol_type`
+        equal to "B" in the input structure being changed to "C" in the output
+        structure.
+
+    Returns
+    -------
+    structure_with_substitutions: _xtal.Structure
+        A copy of `structure`, with substitutions of `atom_type` and `mol_type`.
+    """
+    return _xtal.Structure(
+        lattice=structure.lattice(),
+        atom_coordinate_frac=structure.atom_coordinate_frac(),
+        atom_type=[substitutions.get(x, x) for x in structure.atom_type()],
+        atom_properties=structure.atom_properties(),
+        mol_coordinate_frac=structure.mol_coordinate_frac(),
+        mol_type=[substitutions.get(x, x) for x in structure.mol_type()],
+        mol_properties=structure.mol_properties(),
+        global_properties=structure.global_properties(),
+    )
diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp
index 1dbe25d..b902099 100644
--- a/python/src/xtal.cpp
+++ b/python/src/xtal.cpp
@@ -1831,9 +1831,9 @@ PYBIND11_MODULE(_xtal, m) {
             Parameters
             ----------
             frac : boolean, default=True
-                If True, write basis site positions in fractional coordinates
-                relative to the lattice vectors. If False, write basis site positions
-                in Cartesian coordinates.
+                By default, basis site positions are written in fractional
+                coordinates relative to the lattice vectors. If False, write basis site
+                positions in Cartesian coordinates.
             include_va : boolean, default=False
                 If a basis site only allows vacancies, it is not printed by default.
                 If this is True, basis sites with only vacancies will be included.
@@ -1964,7 +1964,7 @@ PYBIND11_MODULE(_xtal, m) {
             Parameters
             ----------
             frac : boolean, default=True
-                If True, write basis site positions in fractional coordinates
+                By default, basis site positions are written in fractional coordinates
                 relative to the lattice vectors. If False, write basis site positions
                 in Cartesian coordinates.
             include_va : boolean, default=False
@@ -2554,7 +2554,7 @@ PYBIND11_MODULE(_xtal, m) {
           },
           py::arg("excluded_species") =
               std::vector<std::string>({"Va", "VA", "va"}),
-          py::arg("frac") = false, R"pbdoc(
+          py::arg("frac") = true, R"pbdoc(
           Represent the Structure as a Python dict.
 
           Parameters
@@ -2562,10 +2562,9 @@ PYBIND11_MODULE(_xtal, m) {
           excluded_species : list[str] = ["Va", "VA", "va"]
               The names of any molecular or atomic species that should not be included
               in the output.
-          frac : boolean, default=False
-                If True, write coordinates using fractional coordinates
-                relative to the lattice vectors. If False, write coordinates using
-                Cartesian coordinates.
+          frac : boolean, default=True
+              By default, coordinates are written in fractional coordinates relative to
+              the lattice vectors. If False, write coordinates in Cartesian coordinates.
 
           Returns
           -------
diff --git a/python/tests/test_structure.py b/python/tests/test_structure.py
index eadeb17..5e4bfa0 100644
--- a/python/tests/test_structure.py
+++ b/python/tests/test_structure.py
@@ -77,7 +77,7 @@ def test_structure_to_dict(example_structure_1):
     assert len(data["global_properties"]) == 1
     assert "Hstrain" in data["global_properties"]
     expected = np.array(
-        [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5], [0.0, 0.0, 1.0], [0.5, 0.5, 1.5]]
+        [[0.0, 0.0, 0.0], [0.5, 0.5, 0.25], [0.0, 0.0, 0.5], [0.5, 0.5, 0.75]]
     )
     print(xtal.pretty_json(data["atom_coords"]))
     assert np.allclose(np.array(data["atom_coords"]), expected)
@@ -87,9 +87,9 @@ def test_structure_to_dict(example_structure_1):
     data = structure.to_dict(excluded_species=["B"])
     assert data["atom_type"] == ["A", "A"]
 
-    data = structure.to_dict(frac=True)
+    data = structure.to_dict(frac=False)
     expected = np.array(
-        [[0.0, 0.0, 0.0], [0.5, 0.5, 0.25], [0.0, 0.0, 0.5], [0.5, 0.5, 0.75]]
+        [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5], [0.0, 0.0, 1.0], [0.5, 0.5, 1.5]]
     )
     print(xtal.pretty_json(data["atom_coords"]))
     assert np.allclose(np.array(data["atom_coords"]), expected)
@@ -695,3 +695,21 @@ def test_structure_sort_structure_by_atom_coordinate_cart():
     ).transpose()
     # print(xtal.pretty_json(sorted_struc.to_dict()))
     assert np.allclose(sorted_struc.atom_coordinate_cart(), expected)
+
+
+def test_substitute_structure_species_1(example_structure_1):
+    assert example_structure_1.atom_type() == ["A", "A", "B", "B"]
+    s2 = xtal.substitute_structure_species(
+        example_structure_1,
+        {"A": "C"},
+    )
+    assert s2.atom_type() == ["C", "C", "B", "B"]
+
+
+def test_substitute_structure_species_2(example_structure_1):
+    assert example_structure_1.atom_type() == ["A", "A", "B", "B"]
+    s2 = xtal.substitute_structure_species(
+        example_structure_1,
+        {"A": "C", "B": "D"},
+    )
+    assert s2.atom_type() == ["C", "C", "D", "D"]

From b983680a1aa62365abb2d3b3efc9e48fb600e42c Mon Sep 17 00:00:00 2001
From: bpuchala <bpuchala@users.noreply.github.com>
Date: Sat, 2 Dec 2023 10:32:58 -0500
Subject: [PATCH 6/6] Update CHANGELOG

---
 CHANGELOG.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e4bffe3..399cc4f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [2.0a9] - X
 
+### Fixed
+
+- Fix CASM::xtal::make_primitive, which was not copying unique_names. This also fixes the output of libcasm.xtal.make_primitive which was losing the occ_dof list as a result.
+
 ### Changed
 
 - Changed make_primitive to act on either Prim or Structure.