Occupation site basis function updates

CHANGELOG.md

@@ -17,13 +17,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `to_dict` methods to `ClusterFunctionsBuilder`, `OrbitMatrixRepBuilder` and `ClusterMatrixRepBuilder`
 - Added `to_dict` and `from_dict` methods to `ExponentSumConstraint`
 - Added `casm.bset.json_io` module
+- Added option to specify occ_site_basis_functions_specs of type occupation with a choice of reference occupant on each sublattice.
+- Added a check that prints an error message and raises an exception if the constant occupation site basis function mixes with other site basis functions.
 
 ### Changed
 
 - Updated `ClexulatorWriter.write` to write a "variables.json.gz" file for each Clexulator (including local Clexulator) which contains the variables used by the jinja2 templates as well as information like basis function formulas generated during the write process
 - Updated `ClexulatorWriter.write` to write a "cluster_functions.json.gz" file with the generated clusters, matrix reps, and functions.
 - Changed the name used for occupation variables. The name "\\phi" is the base, but the number and meaning of the indices can vary depending on the site basis functions. An "occ_site_functions_info" dict in the "variables.json" file contains the string value used as a template ("occ_var_name") and a description of the variable including its indices, if any ("occ_var_desc" and "occ_var_indices").
 
+### Fixed
+
+- Fixed `v1.basic` templates so that the Clexulator compiles for a prim with occupation DoF, but only the constant basis  function is included.
+- Fixed tests with anisotropic occupants and occupation site basis functions.
+
 
 ## [2.0a1] - 2024-08-15
 

casm/bset/_methods.py

@@ -1,7 +1,7 @@
 """Simple top-level casm.bset methods for the most common use cases"""
 
 import pathlib
-from typing import Optional, Union
+from typing import Any, Optional, Union
 
 import numpy as np
 
@@ -44,7 +44,7 @@ def make_clex_basis_specs(
     custom_generators: Optional[list[ClusterOrbitGenerator]] = [],
     phenomenal: Union[Cluster, OccEvent, None] = None,
     cutoff_radius: Optional[list[float]] = [],
-    occ_site_basis_functions_specs: Union[str, list[dict], None] = None,
+    occ_site_basis_functions_specs: Any = None,
     global_max_poly_order: Optional[int] = None,
     orbit_branch_max_poly_order: Optional[dict] = None,
 ) -> ClexBasisSpecs:
@@ -81,10 +81,9 @@ def make_clex_basis_specs(
         site to include in the local environment, by number of sites in the cluster.
         The null cluster value (element 0) is arbitrary.
 
-    occ_site_basis_functions_specs: Union[str, list[dict], None] = None
+    occ_site_basis_functions_specs: Any = None
         Provides instructions for constructing occupation site basis functions.
-        The accepted options are "chebychev", "occupation", or a `list[dict]`
-        a specifying sublattice-specific choice of site basis functions. This
+        The most common options are "chebychev" or "occupation". This
         parameter corresponds to the value of
 
         .. code-block:: Python

casm/bset/clexwriter/_write_v1_basic.py

@@ -111,7 +111,7 @@ def __init__(
             - `"occ_var_name"`: str, A variable name template for the site functions,
               which may be formated using `b` for sublattice index and `m` for site
               function index (i..e ``occ_var_name.format(b=0, m=1)``).
-            - `"occ_var_desc": str, A description of the occupation
+            - `"occ_var_desc"`: str, A description of the occupation
               variable, including a description of the subscript indices.
             - `"occ_var_indices"`: list[list[str, str]], A list of lists, where each
               sublist contains the variable name and description for each subscript
@@ -212,7 +212,7 @@ def __init__(
         - `"occ_var_name"`: str, A variable name template for the site functions,
             which may be formated using `b` for sublattice index and `m` for site 
             function index (i..e ``occ_var_name.format(b=0, m=1)``).
-        - `"occ_var_desc": str, A description of the occupation
+        - `"occ_var_desc"`: str, A description of the occupation
             variable, including a description of the subscript indices.
             
         """

casm/bset/cluster_functions/_cluster_functions.py

@@ -1,4 +1,4 @@
-from typing import Callable, Iterable, Optional, Union
+from typing import Any, Callable, Iterable, Optional
 
 import numpy as np
 
@@ -574,7 +574,7 @@ def __init__(
         dofs: Optional[Iterable[str]] = None,
         global_max_poly_order: Optional[int] = None,
         orbit_branch_max_poly_order: Optional[dict[int, int]] = None,
-        occ_site_basis_functions_specs: Union[str, list[dict], None] = None,
+        occ_site_basis_functions_specs: Any = None,
         prim_neighbor_list: Optional[PrimNeighborList] = None,
         make_equivalents: bool = True,
         make_all_local_basis_sets: bool = True,
@@ -620,10 +620,9 @@ def __init__(
             created. Higher order polynomials are requested either according to cluster
             size using `orbit_branch_max_poly_order` or globally using
             `global_max_poly_order`. The most specific level specified is used.
-        occ_site_basis_functions_specs: Union[str, list[dict], None] = None
+        occ_site_basis_functions_specs: Any = None
             Provides instructions for constructing occupation site basis functions.
-            The accepted options are "chebychev", "occupation", or a `list[dict]`
-            a specifying sublattice-specific choice of site basis functions. This
+            The most common options are "chebychev" or "occupation". This
             parameter corresponds to the value of
 
             .. code-block:: Python
@@ -818,7 +817,7 @@ class used is :class:`~casm.bset.cluster_functions.MakeVariableName`.
         - `"occ_var_name"`: str, A variable name template for the site functions,
           which may be formated using `b` for sublattice index and `m` for site function
           index (i..e ``occ_var_name.format(b=0, m=1)``).
-        - `"occ_var_desc": str, A description of the occupation
+        - `"occ_var_desc"`: str, A description of the occupation
           variable, including a description of the subscript indices.
         """
 

casm/bset/cluster_functions/_discrete_functions.py

@@ -1,4 +1,4 @@
-from typing import Any, Optional, Union
+from typing import Any, Optional
 
 import numpy as np
 
@@ -281,6 +281,11 @@ def make_chebychev_site_functions(
 
 
 def _is_occupation_site_functions(site_basis_functions_specs: Any):
+    if isinstance(site_basis_functions_specs, dict):
+        if "type" in site_basis_functions_specs:
+            return site_basis_functions_specs["type"].lower() == "occupation"
+        else:
+            return False
     if not isinstance(site_basis_functions_specs, str):
         return False
     return site_basis_functions_specs.lower() == "occupation"
@@ -289,6 +294,7 @@ def _is_occupation_site_functions(site_basis_functions_specs: Any):
 def make_occupation_site_functions(
     prim: casmconfig.Prim,
     abs_tol: float = 1e-10,
+    reference_occ: list[int] = None,
 ):
     """Make discrete occupation site functions using the "occupation" basis, which \
     expands about the default configuration with each sublattice occupied by the first \
@@ -302,6 +308,11 @@ def make_occupation_site_functions(
     abs_tol: float = 1e-10
         A absolute tolerance for comparing values.
 
+    reference_occ: list[str] = None
+        For each sublattice, the names of the occupant that should only give a non-zero
+        value for the constant site function. If not provided, the first occupant
+        listed in the prim is used.
+
     Returns
     -------
     occ_site_functions: list[dict]
@@ -318,6 +329,22 @@ def make_occupation_site_functions(
     site_rep = prim.integral_site_coordinate_symgroup_rep
     indicator_matrix_rep = prim.local_dof_matrix_rep("occ")
 
+    if reference_occ is None:
+        reference_occ = [site_occ_dofs[0] for site_occ_dofs in occ_dofs]
+
+    if len(reference_occ) != len(occ_dofs):
+        raise Exception(
+            "Error in make_occupation_site_functions: "
+            "reference_occ must have the same length as the number of sublattices."
+        )
+    for b in range(len(occ_dofs)):
+        if reference_occ[b] not in occ_dofs[b]:
+            raise Exception(
+                "Error in make_occupation_site_functions: "
+                f"for sublattice {b}, reference_occ[{b}]={reference_occ[b]} "
+                "is not an allowed occupant."
+            )
+
     phi = {}
     for unit_indices in all_indices:
         b_init = unit_indices[0]
@@ -327,7 +354,8 @@ def make_occupation_site_functions(
             continue
 
         occ_probs = np.zeros((n_allowed_occs,))
-        occ_probs[0] = 1.0
+        reference_index = occ_dofs[b_init].index(reference_occ[b_init])
+        occ_probs[reference_index] = 1.0
         phi[b_init] = make_orthonormal_discrete_functions(occ_probs, abs_tol)
 
         site_init = xtal.IntegralSiteCoordinate(sublattice=b_init, unitcell=[0, 0, 0])
@@ -601,7 +629,7 @@ def make_direct_site_functions(
 
 def make_occ_site_functions(
     prim: casmconfig.Prim,
-    occ_site_basis_functions_specs: Union[str, list[dict]],
+    occ_site_basis_functions_specs: Any,
     abs_tol: float = 1e-10,
 ) -> list[dict]:
     """Make discrete occupation site functions from `dof_specs` input
@@ -611,7 +639,7 @@ def make_occ_site_functions(
     prim: libcasm.configation.Prim
         The prim, with symmetry information.
 
-    occ_site_basis_functions_specs: Union[str, list[dict]]
+    occ_site_basis_functions_specs: Any
         Provides instructions for constructing occupation site basis functions. As
         described in detail in the section
         :ref:`DoF Specifications <sec-dof-specifications>`, the options are:
@@ -624,6 +652,30 @@ def make_occ_site_functions(
           correlation values all equal to `0`) about the default configuration where
           each site is occupied by the first allowed occupant in the prim
           :func:`~libcasm.xtal.Prim.occ_dof` list.
+        - `{"type": "occupation", ...}`: An alternative way to specify "occupation"
+          site basis functions that allows specifying a "reference_occ" list of occupant
+          names to use as the reference occupant for each sublattice. If
+          "reference_occ" is not provided, the first occupant listed in the prim is
+          used. The reference occupant must the same on all symmetrically equivalent
+          sites. For example, an FCC A-B-C binary could use any of the following:
+
+          .. code-block:: Python
+
+              A_ref_specs = {
+                  "type": "occupation",
+                  "reference_occ": ["A"],
+              }
+
+              B_ref_specs = {
+                  "type": "occupation",
+                  "reference_occ": ["B"],
+              }
+
+              C_ref_specs = {
+                  "type": "occupation",
+                  "reference_occ": ["C"],
+              }
+
         - `list[dict]`: A list of dictionaries in one of the following formats:
 
           - Composition-centered basis functions, which give an expansion (with
@@ -700,7 +752,14 @@ def make_occ_site_functions(
     if _is_chebychev_site_functions(x):
         return make_chebychev_site_functions(prim=prim, abs_tol=abs_tol)
     elif _is_occupation_site_functions(x):
-        return make_occupation_site_functions(prim=prim, abs_tol=abs_tol)
+        reference_occ = None
+        if isinstance(x, dict):
+            reference_occ = x.get("reference_occ", None)
+        return make_occupation_site_functions(
+            prim=prim,
+            abs_tol=abs_tol,
+            reference_occ=reference_occ,
+        )
     elif _is_composition_site_functions(x):
         return make_composition_site_functions(
             site_basis_functions_specs=x, prim=prim, abs_tol=abs_tol

casm/bset/cluster_functions/_matrix_rep.py

@@ -233,7 +233,7 @@ class MakeVariableName:
 
     - For continuous DoF, use the :class:`DoFSetBasis` axis names
     - For discrete DoF, use a template string. For example, the default is
-      ``"\\phi_{{{b},{m}}}"` where "{b}" is replaced by the sublattice index and
+      ``"\\phi_{{{b},{m}}}"`` where "{b}" is replaced by the sublattice index and
       "{m}" is replaced by the site function index, a number in the range 0 to M-1,
       where M is the number of occupants allowed on the site.
 
@@ -963,17 +963,42 @@ def make_occ_site_functions_matrix_rep(
     """
     site_rep = integral_site_coordinate_symgroup_rep
 
+    # P = B * Phi, P: indicator functions (as rows), Phi: site functions (as rows)
+
+    # Get Phi and Phi^-1 for each sublattice
     phi = []
     phi_inv = []
+    constant_function_index = []
     for sublattice_index, indicator_rep in enumerate(indicator_matrix_rep[0]):
         site_functions = get_occ_site_functions(
             occ_site_functions=occ_site_functions,
             sublattice_index=sublattice_index,
         )
+
+        # Check for constant site function
+        if site_functions.shape == (0, 0):
+            constant_function_index.append(None)
+        else:
+            found_constant_function_index = False
+            for function_index, site_function in enumerate(site_functions):
+                if (site_function == 1.0).all():
+                    found_constant_function_index = True
+                    constant_function_index.append(function_index)
+            if found_constant_function_index is False:
+                raise ValueError(
+                    "Error in make_occ_site_functions_matrix_rep: "
+                    "No constant function found for the occupation site basis "
+                    f"functions on sublattice {sublattice_index}."
+                )
+
         phi.append(site_functions)
         phi_inv.append(np.linalg.inv(site_functions))
 
+    # Generate the matrix rep for transforming phi, M_phi,
+    # from the matrix rep for transforming indicator functions, M_indicator:
+    #     M_phi[b_init] = phi[b_final] * M_indicator[b_init] * phi_inv[b_init]
     occ_site_functions_matrix_rep = []
+    mixing = False
     for prim_factor_group_index, indicator_op_rep in enumerate(indicator_matrix_rep):
         site_functions_op_rep = []
         for sublattice_index, indicator_rep in enumerate(indicator_op_rep):
@@ -990,9 +1015,56 @@ def make_occ_site_functions_matrix_rep(
                 site_functions_op_rep.append(np.zeros((0, 0)))
                 continue
             site_function_rep = phi[b_final] @ indicator_rep @ phi_inv[b]
+
+            x = site_function_rep[:, constant_function_index[b]]
+            if np.isclose(x, 0.0).sum() != len(x) - 1:
+                mixing = True
+
             site_functions_op_rep.append(site_function_rep)
         occ_site_functions_matrix_rep.append(site_functions_op_rep)
 
+    if mixing:
+        print()
+        print(
+            "\n"
+            "**** Error in make_occ_site_functions_matrix_rep: **********************\n"
+            "* Invalid choice of occupation DoF site basis functions.               *\n"
+            "* The constant function is mixing with non-constant functions.         *\n"
+            "* Try using chebychev site basis functions or an isotropic default     *\n"
+            "* occupant.                                                            *\n"
+            "************************************************************************\n"
+            "\n"
+        )
+
+        print("Occupation site basis functions:")
+        for sublattice_index, indicator_rep in enumerate(indicator_matrix_rep[0]):
+            print(f"~ sublattice_index: {sublattice_index} ~")
+            site_functions = get_occ_site_functions(
+                occ_site_functions=occ_site_functions,
+                sublattice_index=sublattice_index,
+            )
+            print(site_functions)
+            print()
+        print()
+
+        print("Site function matrix reps:")
+        for i_fg, site_functions_op_rep in enumerate(occ_site_functions_matrix_rep):
+            print(f"~~~ prim_factor_group_index: {i_fg} ~~~")
+            for b, site_function_rep in enumerate(site_functions_op_rep):
+                print(f"~ sublattice_index: {b} ~")
+                x = site_function_rep[:, constant_function_index[b]]
+                if np.isclose(x, 0.0).sum() != len(x) - 1:
+                    print("********** Mixing **********")
+                print("M:\n", site_function_rep)
+                print()
+
+        raise ValueError(
+            "Error in make_occ_site_functions_matrix_rep: "
+            "Invalid choice of occupation DoF site basis functions. "
+            "The constant function is mixing with non-constant functions. "
+            "Try using chebychev site basis functions."
+        )
+
     return occ_site_functions_matrix_rep
 
 
@@ -1125,6 +1197,15 @@ class used is :class:`~casm.bset.cluster_functions.MakeVariableName`.
                     "occ_site_functions is None for discrete DoF."
                 )
 
+            # # DEBUG:
+            # indicator_matrix_rep = prim.local_dof_matrix_rep(key)
+            # for i_fg, sublat_rep in enumerate(indicator_matrix_rep):
+            #     print(f"~~~ fg: {i_fg} ~~~")
+            #     for i_sublat, M in enumerate(sublat_rep):
+            #         print(f"~ b: {i_sublat} ~")
+            #         print(M)
+            #         print()
+
             # For discrete DoF, we need to construct the matrix rep for the
             # site basis functions actually being used from the matrix rep
             # for the indicator variables
@@ -1184,7 +1265,11 @@ class used is :class:`~casm.bset.cluster_functions.MakeVariableName`.
         ``x_final = M @ x_init``, where `x_init` are the site variables on an initial 
         site on the `i_sublat_init`-th sublattice and `x_final` are the site variables
         on the final site under the application of the `i_factor_group`-th prim factor
-        group operation."""
+        group operation.
+        
+        For example, for occupation DoF, these are the matrix representations 
+        that transform the site basis functions, not the indicator functions.
+        """
 
         self.cluster_group = cluster_group
         """libcasm.sym_info.SymGroup: The subgroup of the generating group that leaves \