API Reference

Geometry Generation

cdl_string_to_geometry

Generate geometry directly from a CDL string.

from crystal_geometry import cdl_string_to_geometry

geom = cdl_string_to_geometry("cubic[m3m]:{111}")

crystal_geometry.cdl_string_to_geometry(cdl, c_ratio=1.0)

Convenience function to convert CDL string directly to geometry.

Parameters:

Name	Type	Description	Default
cdl	str	CDL string like "cubic[m3m]:{111}@1.0 + {100}@1.3"	required
c_ratio	float	c/a ratio for non-cubic systems	1.0

Returns:

Type	Description
CrystalGeometry	CrystalGeometry

Source code in src/crystal_geometry/geometry.py

def cdl_string_to_geometry(cdl: str, c_ratio: float = 1.0) -> CrystalGeometry:
    """Convenience function to convert CDL string directly to geometry.

    Args:
        cdl: CDL string like "cubic[m3m]:{111}@1.0 + {100}@1.3"
        c_ratio: c/a ratio for non-cubic systems

    Returns:
        CrystalGeometry
    """
    desc = parse_cdl(cdl)
    return cdl_to_geometry(desc, c_ratio)

cdl_to_geometry

Generate geometry from a parsed CrystalDescription.

from cdl_parser import parse_cdl
from crystal_geometry import cdl_to_geometry

desc = parse_cdl("cubic[m3m]:{111}@1.0 + {100}@1.3")
geom = cdl_to_geometry(desc, c_ratio=1.0)

crystal_geometry.cdl_to_geometry(desc, c_ratio=1.0)

Convert CDL description to 3D geometry.

Supports crystalline descriptions, amorphous descriptions, nested growth, and aggregates.

Parameters:

Name	Type	Description	Default
desc	CrystalDescription | AmorphousDescription	Parsed CDL description (CrystalDescription or AmorphousDescription)	required
c_ratio	float	c/a ratio for non-cubic systems	1.0

Returns:

Type	Description
CrystalGeometry	CrystalGeometry with vertices and faces

Source code in src/crystal_geometry/geometry.py

def cdl_to_geometry(
    desc: CrystalDescription | AmorphousDescription, c_ratio: float = 1.0
) -> CrystalGeometry:
    """Convert CDL description to 3D geometry.

    Supports crystalline descriptions, amorphous descriptions,
    nested growth, and aggregates.

    Args:
        desc: Parsed CDL description (CrystalDescription or AmorphousDescription)
        c_ratio: c/a ratio for non-cubic systems

    Returns:
        CrystalGeometry with vertices and faces
    """
    # Handle amorphous descriptions
    if isinstance(desc, AmorphousDescription):
        from .amorphous import generate_amorphous_shape

        shape = desc.shapes[0] if desc.shapes else "massive"
        return generate_amorphous_shape(shape)

    lattice = get_lattice_for_system(desc.system, c_ratio)

    # Check for NestedGrowth or AggregateSpec nodes in the form tree
    # (before flattening, since flat_forms() strips these)
    nested_nodes: list[NestedGrowth] = []
    aggregate_nodes: list[AggregateSpec] = []
    regular_forms: list[FormNode] = []

    for form_node in desc.forms:
        if isinstance(form_node, NestedGrowth):
            nested_nodes.append(form_node)
        elif isinstance(form_node, AggregateSpec):
            aggregate_nodes.append(form_node)
        else:
            regular_forms.append(form_node)

    # Handle nested growth as primary geometry
    if nested_nodes:
        geometry = _generate_nested_growth(desc, nested_nodes[0], lattice, c_ratio)
        if desc.modifications:
            from .modifications import apply_modifications

            geometry = CrystalGeometry(
                vertices=apply_modifications(geometry.vertices, desc.modifications),
                faces=geometry.faces,
                face_normals=geometry.face_normals,
                face_forms=geometry.face_forms,
                face_millers=geometry.face_millers,
                forms=geometry.forms,
                component_ids=geometry.component_ids,
            )
        return geometry

    # Handle aggregates
    if aggregate_nodes:
        agg = aggregate_nodes[0]
        # Build geometry for the aggregate's base form
        from cdl_parser.models import CrystalForm, FormGroup

        if isinstance(agg.form, (CrystalForm, FormGroup)):
            agg_desc = CrystalDescription(
                system=desc.system,
                point_group=desc.point_group,
                forms=[agg.form],
            )
        else:
            agg_desc = CrystalDescription(
                system=desc.system,
                point_group=desc.point_group,
                forms=desc.forms[:1] if desc.forms else [],
            )

        base_geom = cdl_to_geometry(agg_desc, c_ratio)

        from .aggregates import generate_aggregate

        spacing_val: float | None = None
        if agg.spacing:
            # Parse spacing string (e.g. "0.5mm" -> 0.5)
            try:
                spacing_val = float("".join(c for c in agg.spacing if c in "0123456789."))
            except ValueError:
                spacing_val = None

        geometry = generate_aggregate(
            base_geometry=base_geom,
            arrangement=agg.arrangement,
            count=agg.count,
            spacing=spacing_val,
            orientation=agg.orientation,
        )
        return geometry

    # Standard crystalline geometry
    normals, distances, face_form_indices, face_millers = _build_halfspaces(desc, lattice)

    # Generate geometry (twinned or base)
    # For twins with deferred transforms (Japan Law), pass modifications
    # so they can be applied before the transform
    if desc.twin is not None:
        geometry = _generate_twinned_geometry(
            desc,
            normals,
            distances,
            face_form_indices,
            face_millers,
            modifications=desc.modifications,
        )
        # Modifications were applied inside for deferred-transform twins
        # For others, apply here
        if desc.modifications and not _uses_deferred_transform(desc):
            from .modifications import apply_modifications

            geometry = CrystalGeometry(
                vertices=apply_modifications(geometry.vertices, desc.modifications),
                faces=geometry.faces,
                face_normals=geometry.face_normals,
                face_forms=geometry.face_forms,
                face_millers=geometry.face_millers,
                forms=geometry.forms,
                component_ids=geometry.component_ids,
                twin_metadata=geometry.twin_metadata,
            )
    else:
        geometry = _generate_base_geometry(
            normals, distances, face_form_indices, face_millers, desc.flat_forms()
        )
        # Apply modifications if present
        if desc.modifications:
            from .modifications import apply_modifications

            geometry = CrystalGeometry(
                vertices=apply_modifications(geometry.vertices, desc.modifications),
                faces=geometry.faces,
                face_normals=geometry.face_normals,
                face_forms=geometry.face_forms,
                face_millers=geometry.face_millers,
                forms=geometry.forms,
                component_ids=geometry.component_ids,
                twin_metadata=geometry.twin_metadata,
            )

    return geometry

Convenience Constructors

create_octahedron

from crystal_geometry import create_octahedron

octahedron = create_octahedron(scale=1.0)

crystal_geometry.create_octahedron(scale=1.0)

Create a regular octahedron.

Parameters:

Name	Type	Description	Default
scale	float	Distance from origin to vertices	1.0

Returns:

Type	Description
CrystalGeometry	CrystalGeometry for octahedron

Source code in src/crystal_geometry/geometry.py

def create_octahedron(scale: float = 1.0) -> CrystalGeometry:
    """Create a regular octahedron.

    Args:
        scale: Distance from origin to vertices

    Returns:
        CrystalGeometry for octahedron
    """
    return cdl_string_to_geometry(f"cubic[m3m]:{{111}}@{scale}")

create_cube

from crystal_geometry import create_cube

cube = create_cube(scale=1.0)

crystal_geometry.create_cube(scale=1.0)

Create a cube.

Parameters:

Name	Type	Description	Default
scale	float	Distance from origin to face centers	1.0

Returns:

Type	Description
CrystalGeometry	CrystalGeometry for cube

Source code in src/crystal_geometry/geometry.py

def create_cube(scale: float = 1.0) -> CrystalGeometry:
    """Create a cube.

    Args:
        scale: Distance from origin to face centers

    Returns:
        CrystalGeometry for cube
    """
    return cdl_string_to_geometry(f"cubic[m3m]:{{100}}@{scale}")

create_dodecahedron

from crystal_geometry import create_dodecahedron

dodecahedron = create_dodecahedron(scale=1.0)

crystal_geometry.create_dodecahedron(scale=1.0)

Create a rhombic dodecahedron.

Parameters:

Name	Type	Description	Default
scale	float	Distance from origin to face centers	1.0

Returns:

Type	Description
CrystalGeometry	CrystalGeometry for dodecahedron

Source code in src/crystal_geometry/geometry.py

def create_dodecahedron(scale: float = 1.0) -> CrystalGeometry:
    """Create a rhombic dodecahedron.

    Args:
        scale: Distance from origin to face centers

    Returns:
        CrystalGeometry for dodecahedron
    """
    return cdl_string_to_geometry(f"cubic[m3m]:{{110}}@{scale}")

create_truncated_octahedron

from crystal_geometry import create_truncated_octahedron

truncated = create_truncated_octahedron(octahedron_scale=1.0, cube_scale=1.3)

crystal_geometry.create_truncated_octahedron(octahedron_scale=1.0, cube_scale=1.3)

Create a truncated octahedron (cuboctahedron-like).

Parameters:

Name	Type	Description	Default
octahedron_scale	float	Scale for octahedron faces	1.0
cube_scale	float	Scale for cube truncation	1.3

Returns:

Type	Description
CrystalGeometry	CrystalGeometry

Source code in src/crystal_geometry/geometry.py

def create_truncated_octahedron(
    octahedron_scale: float = 1.0, cube_scale: float = 1.3
) -> CrystalGeometry:
    """Create a truncated octahedron (cuboctahedron-like).

    Args:
        octahedron_scale: Scale for octahedron faces
        cube_scale: Scale for cube truncation

    Returns:
        CrystalGeometry
    """
    return cdl_string_to_geometry(f"cubic[m3m]:{{111}}@{octahedron_scale} + {{100}}@{cube_scale}")

CrystalGeometry Class

The main geometry container class.

from crystal_geometry import CrystalGeometry

geom = cdl_string_to_geometry("cubic[m3m]:{111}")

# Properties
geom.vertices      # Nx3 numpy array of vertex positions
geom.faces         # List of face vertex indices
geom.face_normals  # List of unit normal vectors
geom.face_forms    # Form index for each face
geom.face_millers  # Miller indices for each face

# Methods
geom.get_edges()           # Get unique edges as vertex pairs
geom.center()              # Compute centroid
geom.scale_to_unit()       # Scale to fit unit sphere
geom.translate(offset)     # Translate by vector
geom.euler_characteristic()  # V - E + F (should be 2)
geom.is_valid()            # Verify geometry integrity
geom.to_dict()             # Export to dictionary

crystal_geometry.CrystalGeometry dataclass

3D crystal geometry with vertices and faces.

Attributes:

Name	Type	Description
vertices	ndarray	Nx3 array of vertex positions
faces	list[list[int]]	List of faces, each face is list of vertex indices (counter-clockwise)
face_normals	list[ndarray]	Normal vector for each face
face_forms	list[int]	Which form each face belongs to (index into forms list)
face_millers	list[tuple[int, int, int]]	Miller index (h, k, l) for each face
forms	list[CrystalForm]	Original form definitions from CDL
component_ids	list[int] | None	Optional array of component IDs for each face (for twins)
twin_metadata	TwinMetadata | None	Optional metadata for twinned crystals
is_amorphous	bool	Whether this is an amorphous (non-crystalline) geometry
aggregate_metadata	AggregateMetadata | None	Optional metadata for aggregate geometries

Source code in src/crystal_geometry/models.py

@dataclass
class CrystalGeometry:
    """3D crystal geometry with vertices and faces.

    Attributes:
        vertices: Nx3 array of vertex positions
        faces: List of faces, each face is list of vertex indices (counter-clockwise)
        face_normals: Normal vector for each face
        face_forms: Which form each face belongs to (index into forms list)
        face_millers: Miller index (h, k, l) for each face
        forms: Original form definitions from CDL
        component_ids: Optional array of component IDs for each face (for twins)
        twin_metadata: Optional metadata for twinned crystals
        is_amorphous: Whether this is an amorphous (non-crystalline) geometry
        aggregate_metadata: Optional metadata for aggregate geometries
    """

    vertices: np.ndarray  # Nx3 array of vertex positions
    faces: list[list[int]]  # List of faces, each face is list of vertex indices
    face_normals: list[np.ndarray]  # Normal vector for each face
    face_forms: list[int]  # Which form each face belongs to (index into forms list)
    face_millers: list[tuple[int, int, int]]  # Miller index for each face
    forms: list[CrystalForm] = field(default_factory=list)  # Original form definitions
    component_ids: list[int] | None = None  # Component ID for each face (for twins)
    twin_metadata: TwinMetadata | None = None  # Metadata for twinned crystals
    is_amorphous: bool = False  # Whether this is amorphous geometry
    aggregate_metadata: AggregateMetadata | None = None  # Metadata for aggregates

    def get_edges(self) -> list[tuple[int, int]]:
        """Get all unique edges as vertex index pairs.

        Returns:
            List of (v1, v2) tuples where v1 < v2
        """
        edges: set[tuple[int, int]] = set()
        for face in self.faces:
            n = len(face)
            for i in range(n):
                v1, v2 = face[i], face[(i + 1) % n]
                edge = (min(v1, v2), max(v1, v2))
                edges.add(edge)
        return list(edges)

    def center(self) -> np.ndarray:
        """Get center of geometry (mean of all vertices)."""
        return np.mean(self.vertices, axis=0)

    def scale_to_unit(self) -> "CrystalGeometry":
        """Scale geometry to fit in unit sphere.

        Returns:
            New CrystalGeometry with vertices scaled to max distance of 1.0
        """
        max_dist = np.max(np.linalg.norm(self.vertices, axis=1))
        if max_dist > 0:
            new_verts = self.vertices / max_dist
        else:
            new_verts = self.vertices.copy()

        return CrystalGeometry(
            vertices=new_verts,
            faces=self.faces,
            face_normals=self.face_normals,
            face_forms=self.face_forms,
            face_millers=self.face_millers,
            forms=self.forms,
        )

    def translate(self, offset: np.ndarray) -> "CrystalGeometry":
        """Translate geometry by offset.

        Args:
            offset: 3D translation vector

        Returns:
            New translated CrystalGeometry
        """
        return CrystalGeometry(
            vertices=self.vertices + offset,
            faces=self.faces,
            face_normals=self.face_normals,
            face_forms=self.face_forms,
            face_millers=self.face_millers,
            forms=self.forms,
        )

    def rotate(self, matrix: np.ndarray) -> "CrystalGeometry":
        """Rotate geometry by rotation matrix.

        Args:
            matrix: 3x3 rotation matrix

        Returns:
            New rotated CrystalGeometry
        """
        new_verts = self.vertices @ matrix.T
        new_normals = [n @ matrix.T for n in self.face_normals]

        return CrystalGeometry(
            vertices=new_verts,
            faces=self.faces,
            face_normals=new_normals,
            face_forms=self.face_forms,
            face_millers=self.face_millers,
            forms=self.forms,
        )

    def euler_characteristic(self) -> int:
        """Compute Euler characteristic V - E + F.

        Returns:
            Euler characteristic (should be 2 for convex polyhedra)
        """
        V = len(self.vertices)
        E = len(self.get_edges())
        F = len(self.faces)
        return V - E + F

    def is_valid(self) -> bool:
        """Check if geometry is valid (Euler's formula holds).

        Returns:
            True if V - E + F = 2
        """
        return self.euler_characteristic() == 2

    def to_dict(self) -> dict[str, Any]:
        """Convert to dictionary representation."""
        result: dict[str, Any] = {
            "vertices": self.vertices.tolist(),
            "faces": self.faces,
            "face_normals": [n.tolist() for n in self.face_normals],
            "face_forms": self.face_forms,
            "face_millers": self.face_millers,
            "is_amorphous": self.is_amorphous,
        }
        if self.component_ids is not None:
            result["component_ids"] = self.component_ids
        if self.twin_metadata is not None:
            result["twin_metadata"] = self.twin_metadata.to_dict()
        if self.aggregate_metadata is not None:
            result["aggregate_metadata"] = self.aggregate_metadata.to_dict()
        return result

center()

Get center of geometry (mean of all vertices).

Source code in src/crystal_geometry/models.py

def center(self) -> np.ndarray:
    """Get center of geometry (mean of all vertices)."""
    return np.mean(self.vertices, axis=0)

euler_characteristic()

Compute Euler characteristic V - E + F.

Returns:

Type	Description
int	Euler characteristic (should be 2 for convex polyhedra)

Source code in src/crystal_geometry/models.py

def euler_characteristic(self) -> int:
    """Compute Euler characteristic V - E + F.

    Returns:
        Euler characteristic (should be 2 for convex polyhedra)
    """
    V = len(self.vertices)
    E = len(self.get_edges())
    F = len(self.faces)
    return V - E + F

get_edges()

Get all unique edges as vertex index pairs.

Returns:

Type	Description
list[tuple[int, int]]	List of (v1, v2) tuples where v1 < v2

Source code in src/crystal_geometry/models.py

def get_edges(self) -> list[tuple[int, int]]:
    """Get all unique edges as vertex index pairs.

    Returns:
        List of (v1, v2) tuples where v1 < v2
    """
    edges: set[tuple[int, int]] = set()
    for face in self.faces:
        n = len(face)
        for i in range(n):
            v1, v2 = face[i], face[(i + 1) % n]
            edge = (min(v1, v2), max(v1, v2))
            edges.add(edge)
    return list(edges)

is_valid()

Check if geometry is valid (Euler's formula holds).

Returns:

Type	Description
bool	True if V - E + F = 2

Source code in src/crystal_geometry/models.py

def is_valid(self) -> bool:
    """Check if geometry is valid (Euler's formula holds).

    Returns:
        True if V - E + F = 2
    """
    return self.euler_characteristic() == 2

rotate(matrix)

Rotate geometry by rotation matrix.

Parameters:

Name	Type	Description	Default
matrix	ndarray	3x3 rotation matrix	required

Returns:

Type	Description
CrystalGeometry	New rotated CrystalGeometry

Source code in src/crystal_geometry/models.py

def rotate(self, matrix: np.ndarray) -> "CrystalGeometry":
    """Rotate geometry by rotation matrix.

    Args:
        matrix: 3x3 rotation matrix

    Returns:
        New rotated CrystalGeometry
    """
    new_verts = self.vertices @ matrix.T
    new_normals = [n @ matrix.T for n in self.face_normals]

    return CrystalGeometry(
        vertices=new_verts,
        faces=self.faces,
        face_normals=new_normals,
        face_forms=self.face_forms,
        face_millers=self.face_millers,
        forms=self.forms,
    )

scale_to_unit()

Scale geometry to fit in unit sphere.

Returns:

Type	Description
CrystalGeometry	New CrystalGeometry with vertices scaled to max distance of 1.0

Source code in src/crystal_geometry/models.py

def scale_to_unit(self) -> "CrystalGeometry":
    """Scale geometry to fit in unit sphere.

    Returns:
        New CrystalGeometry with vertices scaled to max distance of 1.0
    """
    max_dist = np.max(np.linalg.norm(self.vertices, axis=1))
    if max_dist > 0:
        new_verts = self.vertices / max_dist
    else:
        new_verts = self.vertices.copy()

    return CrystalGeometry(
        vertices=new_verts,
        faces=self.faces,
        face_normals=self.face_normals,
        face_forms=self.face_forms,
        face_millers=self.face_millers,
        forms=self.forms,
    )

to_dict()

Convert to dictionary representation.

Source code in src/crystal_geometry/models.py

def to_dict(self) -> dict[str, Any]:
    """Convert to dictionary representation."""
    result: dict[str, Any] = {
        "vertices": self.vertices.tolist(),
        "faces": self.faces,
        "face_normals": [n.tolist() for n in self.face_normals],
        "face_forms": self.face_forms,
        "face_millers": self.face_millers,
        "is_amorphous": self.is_amorphous,
    }
    if self.component_ids is not None:
        result["component_ids"] = self.component_ids
    if self.twin_metadata is not None:
        result["twin_metadata"] = self.twin_metadata.to_dict()
    if self.aggregate_metadata is not None:
        result["aggregate_metadata"] = self.aggregate_metadata.to_dict()
    return result

translate(offset)

Translate geometry by offset.

Parameters:

Name	Type	Description	Default
offset	ndarray	3D translation vector	required

Returns:

Type	Description
CrystalGeometry	New translated CrystalGeometry

Source code in src/crystal_geometry/models.py

def translate(self, offset: np.ndarray) -> "CrystalGeometry":
    """Translate geometry by offset.

    Args:
        offset: 3D translation vector

    Returns:
        New translated CrystalGeometry
    """
    return CrystalGeometry(
        vertices=self.vertices + offset,
        faces=self.faces,
        face_normals=self.face_normals,
        face_forms=self.face_forms,
        face_millers=self.face_millers,
        forms=self.forms,
    )

Symmetry Operations

get_point_group_operations

Get symmetry operations for a crystallographic point group.

from crystal_geometry import get_point_group_operations

ops = get_point_group_operations('m3m')  # 48 operations
ops = get_point_group_operations('6/mmm')  # 24 operations

crystal_geometry.get_point_group_operations(point_group)

Get symmetry operations for a point group.

Parameters:

Name	Type	Description	Default
point_group	str	Hermann-Mauguin symbol (e.g., 'm3m', '6/mmm')	required

Returns:

Type	Description
list[ndarray]	List of 3x3 rotation/reflection matrices

Source code in src/crystal_geometry/symmetry.py

def get_point_group_operations(point_group: str) -> list[np.ndarray]:
    """Get symmetry operations for a point group.

    Args:
        point_group: Hermann-Mauguin symbol (e.g., 'm3m', '6/mmm')

    Returns:
        List of 3x3 rotation/reflection matrices
    """
    if point_group in _POINT_GROUP_CACHE:
        return _POINT_GROUP_CACHE[point_group]

    # Define generators for each point group
    generators_map = {
        # Cubic
        "m3m": [C4z, C3_111, I],
        "432": [C4z, C3_111],
        "-43m": [C4z @ I, C3_111],  # S4
        "m-3": [C2z, C3_111, I],
        "23": [C2z, C3_111],
        # Hexagonal
        "6/mmm": [C6z, C2x, Mxy],
        "622": [C6z, C2x],
        "6mm": [C6z, Mxz],
        "-6m2": [C3z, Mxy, Mxz],
        "6/m": [C6z, Mxy],
        "-6": [C3z, Mxy],
        "6": [C6z],
        # Trigonal
        "-3m": [C3z, C2x, I],
        "32": [C3z, C2x],
        "3m": [C3z, Mxz],
        "-3": [C3z, I],
        "3": [C3z],
        # Tetragonal
        "4/mmm": [C4z, C2x, Mxy],
        "422": [C4z, C2x],
        "4mm": [C4z, Mxz],
        "-42m": [C4z @ I, C2x],  # S4
        "4/m": [C4z, Mxy],
        "-4": [C4z @ I],  # S4
        "4": [C4z],
        # Orthorhombic
        "mmm": [C2z, C2x, I],
        "222": [C2z, C2x],
        "mm2": [C2z, Mxz],
        # Monoclinic
        "2/m": [C2y, Mxz],
        "2": [C2y],
        "m": [Mxz],
        # Triclinic
        "-1": [I],
        "1": [],
    }

    if point_group not in generators_map:
        raise ValueError(f"Unknown point group: {point_group}")

    generators = generators_map[point_group]
    if not generators:
        operations = [E.copy()]
    else:
        operations = _generate_group(generators)

    _POINT_GROUP_CACHE[point_group] = operations
    return operations

generate_equivalent_faces

Generate equivalent faces from one Miller index using point group symmetry.

from crystal_geometry import generate_equivalent_faces

faces = generate_equivalent_faces(1, 1, 1, 'm3m')  # 8 faces for {111}
faces = generate_equivalent_faces(1, 0, 0, 'm3m')  # 6 faces for {100}

crystal_geometry.generate_equivalent_faces(h, k, l, point_group, lattice=DEFAULT_LATTICE)

Generate all symmetry-equivalent Miller indices.

Parameters:

Name	Type	Description	Default
h, k, l		Miller indices	required
point_group	str	Hermann-Mauguin point group symbol	required
lattice	LatticeParams	Lattice parameters	DEFAULT_LATTICE

Returns:

Type	Description
list[tuple[int, int, int]]	List of unique (h, k, l) tuples for equivalent faces

Source code in src/crystal_geometry/symmetry.py

def generate_equivalent_faces(
    h: int, k: int, l: int, point_group: str, lattice: LatticeParams = DEFAULT_LATTICE
) -> list[tuple[int, int, int]]:
    """Generate all symmetry-equivalent Miller indices.

    Args:
        h, k, l: Miller indices
        point_group: Hermann-Mauguin point group symbol
        lattice: Lattice parameters

    Returns:
        List of unique (h, k, l) tuples for equivalent faces
    """
    # Use special Miller index transformations for hexagonal/trigonal systems
    if _is_hexagonal_trigonal(point_group):
        operations = _get_hex_point_group_operations(point_group)
        operations_arr = np.array(operations)
    else:
        # For other systems, use cached Cartesian operations
        if point_group not in _POINT_GROUP_ARRAY_CACHE:
            operations = get_point_group_operations(point_group)
            _POINT_GROUP_ARRAY_CACHE[point_group] = np.array(operations)
        operations_arr = _POINT_GROUP_ARRAY_CACHE[point_group]

    miller = np.array([h, k, l], dtype=np.float64)

    # Vectorized: apply all operations at once
    all_millers = operations_arr @ miller  # Shape: (N, 3)

    # Vectorized rounding to integers
    rounded = np.rint(all_millers).astype(np.int32)

    # Convert to set of tuples for uniqueness
    equivalent = {tuple(m) for m in rounded}

    return list(equivalent)

miller_to_normal

Convert Miller indices to a unit normal vector.

from crystal_geometry import miller_to_normal, LatticeParams

lattice = LatticeParams.cubic()
normal = miller_to_normal(1, 1, 1, lattice)

crystal_geometry.miller_to_normal(h, k, l, lattice=DEFAULT_LATTICE)

Convert Miller index to face normal vector.

For non-cubic systems, uses the reciprocal lattice (cached for performance).

Parameters:

Name	Type	Description	Default
h, k, l		Miller indices	required
lattice	LatticeParams	Lattice parameters	DEFAULT_LATTICE

Returns:

Type	Description
ndarray	Unit normal vector

Source code in src/crystal_geometry/symmetry.py

def miller_to_normal(
    h: int, k: int, l: int, lattice: LatticeParams = DEFAULT_LATTICE
) -> np.ndarray:
    """Convert Miller index to face normal vector.

    For non-cubic systems, uses the reciprocal lattice (cached for performance).

    Args:
        h, k, l: Miller indices
        lattice: Lattice parameters

    Returns:
        Unit normal vector
    """
    # For cubic systems, Miller indices directly give the normal
    # For other systems, we need the reciprocal lattice
    if (
        lattice.a == lattice.b == lattice.c
        and lattice.alpha == lattice.beta == lattice.gamma == np.pi / 2
    ):
        # Cubic
        normal = np.array([h, k, l], dtype=float)
    else:
        # Non-cubic: use cached reciprocal lattice vectors
        a_star, b_star, c_star = _get_reciprocal_lattice(
            lattice.a, lattice.b, lattice.c, lattice.alpha, lattice.beta, lattice.gamma
        )

        # Normal in Cartesian coordinates
        normal = h * a_star + k * b_star + l * c_star

    # Normalize
    norm = np.linalg.norm(normal)
    if norm > 0:
        return normal / norm
    return np.array([0, 0, 1])

Lattice Parameters

LatticeParams

Class for defining crystal lattice parameters.

from crystal_geometry import LatticeParams

# Factory methods for different crystal systems
cubic = LatticeParams.cubic()
tetragonal = LatticeParams.tetragonal(c_ratio=1.5)
hexagonal = LatticeParams.hexagonal(c_ratio=1.2)
orthorhombic = LatticeParams.orthorhombic(a=1.0, b=1.2, c=1.5)
monoclinic = LatticeParams.monoclinic(a=1.0, b=1.2, c=1.5, beta=100.0)
triclinic = LatticeParams.triclinic(a=1.0, b=1.2, c=1.5, alpha=80.0, beta=100.0, gamma=110.0)

crystal_geometry.LatticeParams dataclass

Crystal lattice parameters.

Attributes:

Name	Type	Description
a,	(b, c)	Lattice vector lengths
alpha,	(beta, gamma)	Angles between lattice vectors (in radians)

Source code in src/crystal_geometry/models.py

@dataclass
class LatticeParams:
    """Crystal lattice parameters.

    Attributes:
        a, b, c: Lattice vector lengths
        alpha, beta, gamma: Angles between lattice vectors (in radians)
    """

    a: float = 1.0
    b: float = 1.0
    c: float = 1.0
    alpha: float = np.pi / 2  # 90 degrees
    beta: float = np.pi / 2
    gamma: float = np.pi / 2

    @classmethod
    def cubic(cls) -> "LatticeParams":
        """Create cubic lattice (a = b = c, all angles 90°)."""
        return cls(1.0, 1.0, 1.0, np.pi / 2, np.pi / 2, np.pi / 2)

    @classmethod
    def hexagonal(cls, c_ratio: float = 1.633) -> "LatticeParams":
        """Create hexagonal lattice (a = b, gamma = 120°)."""
        return cls(1.0, 1.0, c_ratio, np.pi / 2, np.pi / 2, 2 * np.pi / 3)

    @classmethod
    def tetragonal(cls, c_ratio: float = 1.0) -> "LatticeParams":
        """Create tetragonal lattice (a = b, all angles 90°)."""
        return cls(1.0, 1.0, c_ratio, np.pi / 2, np.pi / 2, np.pi / 2)

cubic() classmethod

Create cubic lattice (a = b = c, all angles 90°).

Source code in src/crystal_geometry/models.py

@classmethod
def cubic(cls) -> "LatticeParams":
    """Create cubic lattice (a = b = c, all angles 90°)."""
    return cls(1.0, 1.0, 1.0, np.pi / 2, np.pi / 2, np.pi / 2)

hexagonal(c_ratio=1.633) classmethod

Create hexagonal lattice (a = b, gamma = 120°).

Source code in src/crystal_geometry/models.py

@classmethod
def hexagonal(cls, c_ratio: float = 1.633) -> "LatticeParams":
    """Create hexagonal lattice (a = b, gamma = 120°)."""
    return cls(1.0, 1.0, c_ratio, np.pi / 2, np.pi / 2, 2 * np.pi / 3)

tetragonal(c_ratio=1.0) classmethod

Create tetragonal lattice (a = b, all angles 90°).

Source code in src/crystal_geometry/models.py

@classmethod
def tetragonal(cls, c_ratio: float = 1.0) -> "LatticeParams":
    """Create tetragonal lattice (a = b, all angles 90°)."""
    return cls(1.0, 1.0, c_ratio, np.pi / 2, np.pi / 2, np.pi / 2)

get_lattice_for_system

Get default lattice parameters for a crystal system.

from crystal_geometry import get_lattice_for_system

lattice = get_lattice_for_system('cubic')
lattice = get_lattice_for_system('hexagonal', c_ratio=1.2)

crystal_geometry.get_lattice_for_system(system, c_ratio=1.0)

Get appropriate lattice parameters for a crystal system.

Parameters:

Name	Type	Description	Default
system	str	Crystal system name	required
c_ratio	float	c/a ratio for non-cubic systems	1.0

Returns:

Type	Description
LatticeParams	LatticeParams for the system

Source code in src/crystal_geometry/symmetry.py

def get_lattice_for_system(system: str, c_ratio: float = 1.0) -> LatticeParams:
    """Get appropriate lattice parameters for a crystal system.

    Args:
        system: Crystal system name
        c_ratio: c/a ratio for non-cubic systems

    Returns:
        LatticeParams for the system
    """
    if system == "cubic":
        return LatticeParams.cubic()
    elif system == "hexagonal" or system == "trigonal":
        return LatticeParams.hexagonal(c_ratio)
    elif system == "tetragonal":
        return LatticeParams.tetragonal(c_ratio)
    else:
        # Orthorhombic, monoclinic, triclinic - use defaults
        return LatticeParams()

Low-Level Functions

halfspace_intersection_3d

Compute the convex polyhedron from half-space intersection.

from crystal_geometry import halfspace_intersection_3d
import numpy as np

# Half-spaces as Ax <= b
A = np.array([...])  # Normals
b = np.array([...])  # Distances

vertices, faces = halfspace_intersection_3d(A, b)

crystal_geometry.halfspace_intersection_3d(normals, distances, interior_point=None)

Compute intersection of half-spaces in 3D.

Each half-space is defined by: normal . x <= distance

Parameters:

Name	Type	Description	Default
normals	list[ndarray] | ndarray	List or array of unit normal vectors pointing outward	required
distances	list[float] | ndarray	List or array of distances from origin to each plane	required
interior_point	ndarray | None	A point known to be inside the intersection	None

Returns:

Type	Description
ndarray | None	Array of vertices, or None if intersection is empty/unbounded

Source code in src/crystal_geometry/geometry.py

@prefer_native
def halfspace_intersection_3d(
    normals: list[np.ndarray] | np.ndarray,
    distances: list[float] | np.ndarray,
    interior_point: np.ndarray | None = None,
) -> np.ndarray | None:
    """Compute intersection of half-spaces in 3D.

    Each half-space is defined by: normal . x <= distance

    Args:
        normals: List or array of unit normal vectors pointing outward
        distances: List or array of distances from origin to each plane
        interior_point: A point known to be inside the intersection

    Returns:
        Array of vertices, or None if intersection is empty/unbounded
    """
    # Ensure contiguous arrays for optimal performance
    normals_arr = np.ascontiguousarray(normals, dtype=np.float64)
    distances_arr = np.ascontiguousarray(distances, dtype=np.float64)

    if interior_point is None:
        # Try Chebyshev center first (most robust)
        interior_point = _find_interior_point(normals_arr, distances_arr)

        if interior_point is None:
            # Fallback to iterative method
            interior_point = _iterative_interior_point(normals_arr, distances_arr)

        if interior_point is None:
            # Last resort: try origin
            interior_point = np.array([0.0, 0.0, 0.0])

    # Build halfspace matrix for scipy
    # Format: [A | -b] where Ax <= b becomes Ax - b <= 0
    halfspaces = np.hstack([normals_arr, -distances_arr.reshape(-1, 1)])

    try:
        hs = HalfspaceIntersection(halfspaces, interior_point)
        return hs.intersections
    except Exception:
        # Try with scaled interior point
        try:
            scaled_point = interior_point * 0.5
            hs = HalfspaceIntersection(halfspaces, scaled_point)
            return hs.intersections
        except Exception:
            return None

compute_face_vertices

Compute ordered vertices for each face of a polyhedron.

from crystal_geometry import compute_face_vertices

face_vertices = compute_face_vertices(vertices, faces)

crystal_geometry.compute_face_vertices(vertices, normal, distance, tolerance=1e-06)

Find vertices that lie on a face plane.

Uses vectorized operations for improved performance.

Parameters:

Name	Type	Description	Default
vertices	ndarray	All vertices (Nx3 array)	required
normal	ndarray	Face normal (3-element array)	required
distance	float	Distance from origin to face plane	required
tolerance	float	Numerical tolerance	1e-06

Returns:

Type	Description
list[int]	List of vertex indices on this face, ordered counter-clockwise

Source code in src/crystal_geometry/geometry.py

@prefer_native
def compute_face_vertices(
    vertices: np.ndarray, normal: np.ndarray, distance: float, tolerance: float = 1e-6
) -> list[int]:
    """Find vertices that lie on a face plane.

    Uses vectorized operations for improved performance.

    Args:
        vertices: All vertices (Nx3 array)
        normal: Face normal (3-element array)
        distance: Distance from origin to face plane
        tolerance: Numerical tolerance

    Returns:
        List of vertex indices on this face, ordered counter-clockwise
    """
    # Ensure contiguous arrays
    vertices = np.ascontiguousarray(vertices, dtype=np.float64)
    normal = np.ascontiguousarray(normal, dtype=np.float64)

    # Vectorized distance computation for all vertices
    projections = vertices @ normal
    on_face_mask = np.abs(projections - distance) < tolerance
    on_face_indices = np.where(on_face_mask)[0]

    if len(on_face_indices) < 3:
        return []

    # Get vertices on face
    on_face_verts = vertices[on_face_indices]

    # Order vertices counter-clockwise when viewed from outside
    center = np.mean(on_face_verts, axis=0)

    # Create local coordinate system on face
    u = on_face_verts[0] - center
    u = u - np.dot(u, normal) * normal
    u_norm = u @ u  # squared norm
    if u_norm < tolerance * tolerance:
        if len(on_face_indices) > 1:
            u = on_face_verts[1] - center
            u = u - np.dot(u, normal) * normal
            u_norm = u @ u
    u = u / (np.sqrt(u_norm) + 1e-10)
    v_axis = np.cross(normal, u)

    # Vectorized angle computation
    vecs = on_face_verts - center
    angles = np.arctan2(vecs @ v_axis, vecs @ u)

    # Sort by angle
    sorted_order = np.argsort(angles)
    return [int(on_face_indices[i]) for i in sorted_order]