mc — Monte Carlo Placement

The Monte Carlo engine behind MC Placement & Chain Growth, also used by bead-spring generation.

Collision detection

AutoPoly.mc.collision

Collision Detection Module for Monte Carlo Placement

This module provides efficient 3D collision detection using monomer-level spheres and cell-linked list data structure for O(N) expected complexity.

Key Features: - Hard-sphere collision model at monomer level - Cell-linked list for spatial hashing - Efficient add/remove operations - Configurable tolerance for overlap detection

Created on 2026-01-29 @author: zwu

Classes:

Name	Description
MonomerSphere	Represents a monomer as a single sphere for collision detection.
CollisionDetector	Efficient collision detection using cell-linked list spatial hashing.

Functions:

Name	Description
calculate_box_size	Calculate cubic box size from monomer count and target density.

MonomerSphere dataclass

Represents a monomer as a single sphere for collision detection.

Attributes:

Name	Type	Description
center	ndarray	3D coordinates of the sphere center (monomer geometric center)
radius	float	Collision radius of the monomer
monomer_id	int	Unique identifier for this monomer

Source code in AutoPoly/mc/collision.py

@dataclass
class MonomerSphere:
    """
    Represents a monomer as a single sphere for collision detection.

    Attributes:
        center: 3D coordinates of the sphere center (monomer geometric center)
        radius: Collision radius of the monomer
        monomer_id: Unique identifier for this monomer
    """
    center: np.ndarray
    radius: float
    monomer_id: int

    def __post_init__(self):
        """Ensure center is a numpy array."""
        if not isinstance(self.center, np.ndarray):
            self.center = np.array(self.center, dtype=np.float64)

CollisionDetector

Efficient collision detection using cell-linked list spatial hashing.

This class implements a 3D grid-based spatial hashing scheme for fast collision queries. Each cell in the grid contains references to spheres whose centers fall within that cell.

The algorithm achieves O(N) expected time for both insertion and collision queries when the cell size is properly tuned to the sphere sizes.

Attributes:

Name	Type	Description
box_bounds		Tuple of ((xmin, xmax), (ymin, ymax), (zmin, zmax))
cell_size		Size of each cell in the spatial grid
cells	Dict[Tuple[int, int, int], Set[int]]	Dictionary mapping cell indices to sets of monomer IDs
spheres	Dict[int, MonomerSphere]	Dictionary mapping monomer IDs to MonomerSphere objects

Example

detector = CollisionDetector( ... box_bounds=((-50, 50), (-50, 50), (-50, 50)), ... cell_size=5.0 ... ) detector.add_monomer(0, np.array([0, 0, 0]), 2.0) detector.check_collision(np.array([1, 0, 0]), 2.0) True

Methods:

Name	Description
estimate_radius	Estimate collision radius for a monomer from connection atom positions.
estimate_radius_from_coords	Estimate collision radius from atom coordinates using bounding sphere.
add_monomer	Add a monomer sphere to the collision detector.
remove_monomer	Remove a monomer sphere from the collision detector.
check_collision	Check if a sphere at the given position would collide with existing spheres.
check_collision_detailed	Check collision with detailed information about the colliding sphere.
check_bounds	Check if a sphere is within the simulation box bounds.
get_sphere_count	Return the total number of spheres in the detector.
clear	Remove all spheres from the detector.

Source code in AutoPoly/mc/collision.py

class CollisionDetector:
    """
    Efficient collision detection using cell-linked list spatial hashing.

    This class implements a 3D grid-based spatial hashing scheme for fast
    collision queries. Each cell in the grid contains references to spheres
    whose centers fall within that cell.

    The algorithm achieves O(N) expected time for both insertion and collision
    queries when the cell size is properly tuned to the sphere sizes.

    Attributes:
        box_bounds: Tuple of ((xmin, xmax), (ymin, ymax), (zmin, zmax))
        cell_size: Size of each cell in the spatial grid
        cells: Dictionary mapping cell indices to sets of monomer IDs
        spheres: Dictionary mapping monomer IDs to MonomerSphere objects

    Example:
        >>> detector = CollisionDetector(
        ...     box_bounds=((-50, 50), (-50, 50), (-50, 50)),
        ...     cell_size=5.0
        ... )
        >>> detector.add_monomer(0, np.array([0, 0, 0]), 2.0)
        >>> detector.check_collision(np.array([1, 0, 0]), 2.0)
        True
    """

    def __init__(
        self,
        box_bounds: Tuple[Tuple[float, float], Tuple[float, float], Tuple[float, float]],
        cell_size: float = 5.0
    ):
        """
        Initialize the collision detector.

        Args:
            box_bounds: Simulation box boundaries as ((xmin, xmax), (ymin, ymax), (zmin, zmax))
            cell_size: Size of each cell in the spatial grid. Should be at least
                      as large as the largest expected sphere diameter for optimal
                      performance.
        """
        self.box_bounds = box_bounds
        self.cell_size = cell_size

        # Calculate grid dimensions
        self._xmin, self._xmax = box_bounds[0]
        self._ymin, self._ymax = box_bounds[1]
        self._zmin, self._zmax = box_bounds[2]

        self._nx = max(1, int(np.ceil((self._xmax - self._xmin) / cell_size)))
        self._ny = max(1, int(np.ceil((self._ymax - self._ymin) / cell_size)))
        self._nz = max(1, int(np.ceil((self._zmax - self._zmin) / cell_size)))

        # Cell storage: maps cell index tuple to set of monomer IDs
        self.cells: Dict[Tuple[int, int, int], Set[int]] = defaultdict(set)

        # Sphere storage: maps monomer ID to MonomerSphere
        self.spheres: Dict[int, MonomerSphere] = {}

    def _get_cell_index(self, position: np.ndarray) -> Tuple[int, int, int]:
        """
        Convert a 3D position to cell indices.

        Args:
            position: 3D coordinates

        Returns:
            Tuple of (ix, iy, iz) cell indices
        """
        ix = int((position[0] - self._xmin) / self.cell_size)
        iy = int((position[1] - self._ymin) / self.cell_size)
        iz = int((position[2] - self._zmin) / self.cell_size)

        # Clamp to valid range
        ix = max(0, min(ix, self._nx - 1))
        iy = max(0, min(iy, self._ny - 1))
        iz = max(0, min(iz, self._nz - 1))

        return (ix, iy, iz)

    def _get_neighbor_cells(
        self,
        cell_index: Tuple[int, int, int],
        search_radius: int = 1
    ) -> List[Tuple[int, int, int]]:
        """
        Get indices of neighboring cells within search radius.

        Args:
            cell_index: Center cell index
            search_radius: Number of cells to search in each direction

        Returns:
            List of valid neighboring cell indices
        """
        ix, iy, iz = cell_index
        neighbors = []

        for di in range(-search_radius, search_radius + 1):
            for dj in range(-search_radius, search_radius + 1):
                for dk in range(-search_radius, search_radius + 1):
                    ni = ix + di
                    nj = iy + dj
                    nk = iz + dk

                    if 0 <= ni < self._nx and 0 <= nj < self._ny and 0 <= nk < self._nz:
                        neighbors.append((ni, nj, nk))

        return neighbors

    @staticmethod
    def estimate_radius(
        left_conn: np.ndarray,
        right_conn: np.ndarray,
        buffer: float = 1.5
    ) -> float:
        """
        Estimate collision radius for a monomer from connection atom positions.

        The radius is estimated as half the distance between connection atoms
        plus a buffer to account for atoms extending beyond the backbone.

        Args:
            left_conn: Coordinates of left connection atom
            right_conn: Coordinates of right connection atom
            buffer: Additional padding in Angstroms (default: 1.5)

        Returns:
            Estimated collision radius in Angstroms
        """
        conn_distance = np.linalg.norm(right_conn - left_conn)
        return conn_distance / 2 + buffer

    @staticmethod
    def estimate_radius_from_coords(
        atom_coords: np.ndarray,
        buffer: float = 0.5
    ) -> float:
        """
        Estimate collision radius from atom coordinates using bounding sphere.

        Computes the center of mass and finds the maximum distance from center
        to any atom, then adds a buffer.

        Args:
            atom_coords: Array of shape (N, 3) containing atom coordinates
            buffer: Additional padding in Angstroms (default: 0.5)

        Returns:
            Estimated collision radius in Angstroms
        """
        if len(atom_coords) == 0:
            return buffer

        center = np.mean(atom_coords, axis=0)
        distances = np.linalg.norm(atom_coords - center, axis=1)
        max_distance = np.max(distances) if len(distances) > 0 else 0

        return max_distance + buffer

    def add_monomer(
        self,
        mon_id: int,
        center: np.ndarray,
        radius: float
    ) -> None:
        """
        Add a monomer sphere to the collision detector.

        Args:
            mon_id: Unique identifier for the monomer
            center: 3D coordinates of the sphere center
            radius: Collision radius of the sphere
        """
        sphere = MonomerSphere(center=center.copy(), radius=radius, monomer_id=mon_id)
        self.spheres[mon_id] = sphere

        # Add to cell
        cell_idx = self._get_cell_index(center)
        self.cells[cell_idx].add(mon_id)

    def remove_monomer(self, mon_id: int) -> None:
        """
        Remove a monomer sphere from the collision detector.

        Args:
            mon_id: Unique identifier of the monomer to remove
        """
        if mon_id not in self.spheres:
            return

        sphere = self.spheres[mon_id]
        cell_idx = self._get_cell_index(sphere.center)

        # Remove from cell
        if cell_idx in self.cells:
            self.cells[cell_idx].discard(mon_id)
            if not self.cells[cell_idx]:
                del self.cells[cell_idx]

        # Remove from spheres
        del self.spheres[mon_id]

    def check_collision(
        self,
        center: np.ndarray,
        radius: float,
        exclude_ids: Optional[Set[int]] = None,
        tolerance: float = 0.1
    ) -> bool:
        """
        Check if a sphere at the given position would collide with existing spheres.

        Args:
            center: 3D coordinates of the test sphere center
            radius: Radius of the test sphere
            exclude_ids: Set of monomer IDs to exclude from collision check
            tolerance: Overlap tolerance in Angstroms. Negative values allow slight
                      penetration before reporting collision. (default: 0.1)

        Returns:
            True if collision detected, False otherwise
        """
        if exclude_ids is None:
            exclude_ids = set()

        # Calculate search radius in cells
        max_search_distance = radius + max(
            s.radius for s in self.spheres.values()
        ) if self.spheres else radius
        search_cells = max(1, int(np.ceil(max_search_distance / self.cell_size)))

        # Get cell and neighbors
        cell_idx = self._get_cell_index(center)
        neighbor_cells = self._get_neighbor_cells(cell_idx, search_cells)

        # Check all spheres in neighboring cells
        for ncell in neighbor_cells:
            if ncell not in self.cells:
                continue

            for mon_id in self.cells[ncell]:
                if mon_id in exclude_ids:
                    continue

                sphere = self.spheres[mon_id]
                distance = np.linalg.norm(center - sphere.center)
                min_distance = radius + sphere.radius - tolerance

                if distance < min_distance:
                    return True

        return False

    def check_collision_detailed(
        self,
        center: np.ndarray,
        radius: float,
        exclude_ids: Optional[Set[int]] = None,
        tolerance: float = 0.1
    ) -> Tuple[bool, Optional[int], Optional[float]]:
        """
        Check collision with detailed information about the colliding sphere.

        Args:
            center: 3D coordinates of the test sphere center
            radius: Radius of the test sphere
            exclude_ids: Set of monomer IDs to exclude from collision check
            tolerance: Overlap tolerance in Angstroms

        Returns:
            Tuple of (collision_detected, colliding_monomer_id, overlap_distance)
        """
        if exclude_ids is None:
            exclude_ids = set()

        max_search_distance = radius + max(
            s.radius for s in self.spheres.values()
        ) if self.spheres else radius
        search_cells = max(1, int(np.ceil(max_search_distance / self.cell_size)))

        cell_idx = self._get_cell_index(center)
        neighbor_cells = self._get_neighbor_cells(cell_idx, search_cells)

        for ncell in neighbor_cells:
            if ncell not in self.cells:
                continue

            for mon_id in self.cells[ncell]:
                if mon_id in exclude_ids:
                    continue

                sphere = self.spheres[mon_id]
                distance = np.linalg.norm(center - sphere.center)
                min_distance = radius + sphere.radius - tolerance

                if distance < min_distance:
                    overlap = min_distance - distance
                    return True, mon_id, overlap

        return False, None, None

    def check_bounds(self, center: np.ndarray, radius: float) -> bool:
        """
        Check if a sphere is within the simulation box bounds.

        Args:
            center: 3D coordinates of the sphere center
            radius: Radius of the sphere

        Returns:
            True if the sphere is completely within bounds, False otherwise
        """
        return bool(
            center[0] - radius >= self._xmin and center[0] + radius <= self._xmax and
            center[1] - radius >= self._ymin and center[1] + radius <= self._ymax and
            center[2] - radius >= self._zmin and center[2] + radius <= self._zmax
        )

    def get_sphere_count(self) -> int:
        """Return the total number of spheres in the detector."""
        return len(self.spheres)

    def clear(self) -> None:
        """Remove all spheres from the detector."""
        self.cells.clear()
        self.spheres.clear()

estimate_radius(left_conn: np.ndarray, right_conn: np.ndarray, buffer: float = 1.5) -> float staticmethod

Estimate collision radius for a monomer from connection atom positions.

The radius is estimated as half the distance between connection atoms plus a buffer to account for atoms extending beyond the backbone.

Parameters:

Name	Type	Description	Default
left_conn	ndarray	Coordinates of left connection atom	required
right_conn	ndarray	Coordinates of right connection atom	required
buffer	float	Additional padding in Angstroms (default: 1.5)	1.5

Returns:

Type	Description
float	Estimated collision radius in Angstroms

Source code in AutoPoly/mc/collision.py

@staticmethod
def estimate_radius(
    left_conn: np.ndarray,
    right_conn: np.ndarray,
    buffer: float = 1.5
) -> float:
    """
    Estimate collision radius for a monomer from connection atom positions.

    The radius is estimated as half the distance between connection atoms
    plus a buffer to account for atoms extending beyond the backbone.

    Args:
        left_conn: Coordinates of left connection atom
        right_conn: Coordinates of right connection atom
        buffer: Additional padding in Angstroms (default: 1.5)

    Returns:
        Estimated collision radius in Angstroms
    """
    conn_distance = np.linalg.norm(right_conn - left_conn)
    return conn_distance / 2 + buffer

estimate_radius_from_coords(atom_coords: np.ndarray, buffer: float = 0.5) -> float staticmethod

Estimate collision radius from atom coordinates using bounding sphere.

Computes the center of mass and finds the maximum distance from center to any atom, then adds a buffer.

Parameters:

Name	Type	Description	Default
atom_coords	ndarray	Array of shape (N, 3) containing atom coordinates	required
buffer	float	Additional padding in Angstroms (default: 0.5)	0.5

Returns:

Type	Description
float	Estimated collision radius in Angstroms

Source code in AutoPoly/mc/collision.py

@staticmethod
def estimate_radius_from_coords(
    atom_coords: np.ndarray,
    buffer: float = 0.5
) -> float:
    """
    Estimate collision radius from atom coordinates using bounding sphere.

    Computes the center of mass and finds the maximum distance from center
    to any atom, then adds a buffer.

    Args:
        atom_coords: Array of shape (N, 3) containing atom coordinates
        buffer: Additional padding in Angstroms (default: 0.5)

    Returns:
        Estimated collision radius in Angstroms
    """
    if len(atom_coords) == 0:
        return buffer

    center = np.mean(atom_coords, axis=0)
    distances = np.linalg.norm(atom_coords - center, axis=1)
    max_distance = np.max(distances) if len(distances) > 0 else 0

    return max_distance + buffer

add_monomer(mon_id: int, center: np.ndarray, radius: float) -> None

Add a monomer sphere to the collision detector.

Parameters:

Name	Type	Description	Default
mon_id	int	Unique identifier for the monomer	required
center	ndarray	3D coordinates of the sphere center	required
radius	float	Collision radius of the sphere	required

Source code in AutoPoly/mc/collision.py

def add_monomer(
    self,
    mon_id: int,
    center: np.ndarray,
    radius: float
) -> None:
    """
    Add a monomer sphere to the collision detector.

    Args:
        mon_id: Unique identifier for the monomer
        center: 3D coordinates of the sphere center
        radius: Collision radius of the sphere
    """
    sphere = MonomerSphere(center=center.copy(), radius=radius, monomer_id=mon_id)
    self.spheres[mon_id] = sphere

    # Add to cell
    cell_idx = self._get_cell_index(center)
    self.cells[cell_idx].add(mon_id)

remove_monomer(mon_id: int) -> None

Remove a monomer sphere from the collision detector.

Parameters:

Name	Type	Description	Default
mon_id	int	Unique identifier of the monomer to remove	required

Source code in AutoPoly/mc/collision.py

def remove_monomer(self, mon_id: int) -> None:
    """
    Remove a monomer sphere from the collision detector.

    Args:
        mon_id: Unique identifier of the monomer to remove
    """
    if mon_id not in self.spheres:
        return

    sphere = self.spheres[mon_id]
    cell_idx = self._get_cell_index(sphere.center)

    # Remove from cell
    if cell_idx in self.cells:
        self.cells[cell_idx].discard(mon_id)
        if not self.cells[cell_idx]:
            del self.cells[cell_idx]

    # Remove from spheres
    del self.spheres[mon_id]

check_collision(center: np.ndarray, radius: float, exclude_ids: Optional[Set[int]] = None, tolerance: float = 0.1) -> bool

Check if a sphere at the given position would collide with existing spheres.

Parameters:

Name	Type	Description	Default
center	ndarray	3D coordinates of the test sphere center	required
radius	float	Radius of the test sphere	required
exclude_ids	Optional[Set[int]]	Set of monomer IDs to exclude from collision check	None
tolerance	float	Overlap tolerance in Angstroms. Negative values allow slight penetration before reporting collision. (default: 0.1)	0.1

Returns:

Type	Description
bool	True if collision detected, False otherwise

Source code in AutoPoly/mc/collision.py

def check_collision(
    self,
    center: np.ndarray,
    radius: float,
    exclude_ids: Optional[Set[int]] = None,
    tolerance: float = 0.1
) -> bool:
    """
    Check if a sphere at the given position would collide with existing spheres.

    Args:
        center: 3D coordinates of the test sphere center
        radius: Radius of the test sphere
        exclude_ids: Set of monomer IDs to exclude from collision check
        tolerance: Overlap tolerance in Angstroms. Negative values allow slight
                  penetration before reporting collision. (default: 0.1)

    Returns:
        True if collision detected, False otherwise
    """
    if exclude_ids is None:
        exclude_ids = set()

    # Calculate search radius in cells
    max_search_distance = radius + max(
        s.radius for s in self.spheres.values()
    ) if self.spheres else radius
    search_cells = max(1, int(np.ceil(max_search_distance / self.cell_size)))

    # Get cell and neighbors
    cell_idx = self._get_cell_index(center)
    neighbor_cells = self._get_neighbor_cells(cell_idx, search_cells)

    # Check all spheres in neighboring cells
    for ncell in neighbor_cells:
        if ncell not in self.cells:
            continue

        for mon_id in self.cells[ncell]:
            if mon_id in exclude_ids:
                continue

            sphere = self.spheres[mon_id]
            distance = np.linalg.norm(center - sphere.center)
            min_distance = radius + sphere.radius - tolerance

            if distance < min_distance:
                return True

    return False

check_collision_detailed(center: np.ndarray, radius: float, exclude_ids: Optional[Set[int]] = None, tolerance: float = 0.1) -> Tuple[bool, Optional[int], Optional[float]]

Check collision with detailed information about the colliding sphere.

Parameters:

Name	Type	Description	Default
center	ndarray	3D coordinates of the test sphere center	required
radius	float	Radius of the test sphere	required
exclude_ids	Optional[Set[int]]	Set of monomer IDs to exclude from collision check	None
tolerance	float	Overlap tolerance in Angstroms	0.1

Returns:

Type	Description
Tuple[bool, Optional[int], Optional[float]]	Tuple of (collision_detected, colliding_monomer_id, overlap_distance)

Source code in AutoPoly/mc/collision.py

def check_collision_detailed(
    self,
    center: np.ndarray,
    radius: float,
    exclude_ids: Optional[Set[int]] = None,
    tolerance: float = 0.1
) -> Tuple[bool, Optional[int], Optional[float]]:
    """
    Check collision with detailed information about the colliding sphere.

    Args:
        center: 3D coordinates of the test sphere center
        radius: Radius of the test sphere
        exclude_ids: Set of monomer IDs to exclude from collision check
        tolerance: Overlap tolerance in Angstroms

    Returns:
        Tuple of (collision_detected, colliding_monomer_id, overlap_distance)
    """
    if exclude_ids is None:
        exclude_ids = set()

    max_search_distance = radius + max(
        s.radius for s in self.spheres.values()
    ) if self.spheres else radius
    search_cells = max(1, int(np.ceil(max_search_distance / self.cell_size)))

    cell_idx = self._get_cell_index(center)
    neighbor_cells = self._get_neighbor_cells(cell_idx, search_cells)

    for ncell in neighbor_cells:
        if ncell not in self.cells:
            continue

        for mon_id in self.cells[ncell]:
            if mon_id in exclude_ids:
                continue

            sphere = self.spheres[mon_id]
            distance = np.linalg.norm(center - sphere.center)
            min_distance = radius + sphere.radius - tolerance

            if distance < min_distance:
                overlap = min_distance - distance
                return True, mon_id, overlap

    return False, None, None

check_bounds(center: np.ndarray, radius: float) -> bool

Check if a sphere is within the simulation box bounds.

Parameters:

Name	Type	Description	Default
center	ndarray	3D coordinates of the sphere center	required
radius	float	Radius of the sphere	required

Returns:

Type	Description
bool	True if the sphere is completely within bounds, False otherwise

Source code in AutoPoly/mc/collision.py

def check_bounds(self, center: np.ndarray, radius: float) -> bool:
    """
    Check if a sphere is within the simulation box bounds.

    Args:
        center: 3D coordinates of the sphere center
        radius: Radius of the sphere

    Returns:
        True if the sphere is completely within bounds, False otherwise
    """
    return bool(
        center[0] - radius >= self._xmin and center[0] + radius <= self._xmax and
        center[1] - radius >= self._ymin and center[1] + radius <= self._ymax and
        center[2] - radius >= self._zmin and center[2] + radius <= self._zmax
    )

get_sphere_count() -> int

Return the total number of spheres in the detector.

Source code in AutoPoly/mc/collision.py

def get_sphere_count(self) -> int:
    """Return the total number of spheres in the detector."""
    return len(self.spheres)

clear() -> None

Remove all spheres from the detector.

Source code in AutoPoly/mc/collision.py

def clear(self) -> None:
    """Remove all spheres from the detector."""
    self.cells.clear()
    self.spheres.clear()

calculate_box_size(total_monomers: int, monomer_density: float = 0.1) -> float

Calculate cubic box size from monomer count and target density.

A low default density (0.1 monomers/A^3) ensures sufficient space between monomers for successful MC placement.

Parameters:

Name	Type	Description	Default
total_monomers	int	Total number of monomers in the system	required
monomer_density	float	Target number density (monomers/A^3), default 0.1	0.1

Returns:

Type	Description
float	Box side length (Angstroms) for a cubic box

Example

calculate_box_size(100, monomer_density=0.1) 21.544... # 100 monomers at density 0.1 -> volume = 1000 A^3 -> box ~ 10 A

Source code in AutoPoly/mc/collision.py

def calculate_box_size(
    total_monomers: int,
    monomer_density: float = 0.1
) -> float:
    """
    Calculate cubic box size from monomer count and target density.

    A low default density (0.1 monomers/A^3) ensures sufficient space
    between monomers for successful MC placement.

    Args:
        total_monomers: Total number of monomers in the system
        monomer_density: Target number density (monomers/A^3), default 0.1

    Returns:
        Box side length (Angstroms) for a cubic box

    Example:
        >>> calculate_box_size(100, monomer_density=0.1)
        21.544...  # 100 monomers at density 0.1 -> volume = 1000 A^3 -> box ~ 10 A
    """
    if monomer_density <= 0:
        raise ValueError("monomer_density must be positive")
    if total_monomers <= 0:
        raise ValueError("total_monomers must be positive")

    volume = total_monomers / monomer_density
    box_size = volume ** (1/3)
    return box_size

Chain growth

AutoPoly.mc.chain_growth

Chain Growth Monte Carlo Module for Self-Avoiding Random Walk

This module implements self-avoiding random walk (SAW) for building polymer chains monomer-by-monomer with collision detection. The chain grows by aligning each new monomer's left connection point to the previous monomer's right connection point.

Key Features: - Parse monomer templates from .lt files - 3D transformation (rotation + translation) for monomer alignment - Self-avoiding random walk with dihedral angle sampling - Generation of moltemplate .rot().move() commands

Connection Point Convention (from .lt files): - First atom in "Data Atoms": Left connection point (C1) - Second atom in "Data Atoms": Right connection point (C2)

Created on 2026-01-29 @author: zwu

Classes:

Name	Description
AtomData	Data for a single atom parsed from .lt file.
MonomerTemplate	Template for a monomer parsed from .lt file.
MonomerPlacement	Represents a placed monomer with its transformation.
ChainGrowthMC	Self-Avoiding Random Walk chain growth using Monte Carlo.

Functions:

Name	Description
parse_lt_file	Parse a monomer .lt file to extract atom coordinates and connection points.
rotation_matrix_from_axis_angle	Create a 3x3 rotation matrix from axis-angle representation.
rotation_matrix_align_vectors	Create rotation matrix that aligns vector v1 to vector v2.
rotation_matrix_to_axis_angle	Convert a 3x3 rotation matrix to axis-angle representation.
random_rotation_matrix	Generate a uniformly distributed random rotation matrix.

AtomData dataclass

Data for a single atom parsed from .lt file.

Attributes:

Name	Type	Description
atom_id	str	Atom identifier (e.g., "C1", "H3")
element	str	Element symbol (e.g., "C", "H")
coords	ndarray	3D coordinates as numpy array
atom_type	str	Force field atom type
charge	float	Partial charge

Source code in AutoPoly/mc/chain_growth.py

@dataclass
class AtomData:
    """
    Data for a single atom parsed from .lt file.

    Attributes:
        atom_id: Atom identifier (e.g., "C1", "H3")
        element: Element symbol (e.g., "C", "H")
        coords: 3D coordinates as numpy array
        atom_type: Force field atom type
        charge: Partial charge
    """
    atom_id: str
    element: str
    coords: np.ndarray
    atom_type: str
    charge: float

MonomerTemplate dataclass

Template for a monomer parsed from .lt file.

Attributes:

Name	Type	Description
lt_file	str	Path to the .lt file
monomer_name	str	Name of the monomer (e.g., "monomer_0_1i")
monomer_type	str	Type of monomer ("first", "middle", "last", or "ring")
atoms	List[AtomData]	List of all atoms with local coordinates
left_conn_coords	Optional[ndarray]	Coordinates of left connection atom (first atom, C1)
right_conn_coords	Optional[ndarray]	Coordinates of right connection atom (second atom, C2)
left_conn_id	Optional[str]	Atom ID of left connection (e.g., "C1") or None for first monomer
right_conn_id	Optional[str]	Atom ID of right connection (e.g., "C2") or None for last monomer
bond_vector	Optional[ndarray]	Vector from left to right connection (for middle monomers)

Methods:

Name	Description
get_center	Calculate the geometric center of all atoms.
get_all_coords	Get array of all atom coordinates.

Source code in AutoPoly/mc/chain_growth.py

@dataclass
class MonomerTemplate:
    """
    Template for a monomer parsed from .lt file.

    Attributes:
        lt_file: Path to the .lt file
        monomer_name: Name of the monomer (e.g., "monomer_0_1i")
        monomer_type: Type of monomer ("first", "middle", "last", or "ring")
        atoms: List of all atoms with local coordinates
        left_conn_coords: Coordinates of left connection atom (first atom, C1)
        right_conn_coords: Coordinates of right connection atom (second atom, C2)
        left_conn_id: Atom ID of left connection (e.g., "C1") or None for first monomer
        right_conn_id: Atom ID of right connection (e.g., "C2") or None for last monomer
        bond_vector: Vector from left to right connection (for middle monomers)
    """
    lt_file: str
    monomer_name: str
    monomer_type: str  # "first", "middle", "last", "ring"
    atoms: List[AtomData]
    left_conn_coords: Optional[np.ndarray]
    right_conn_coords: Optional[np.ndarray]
    left_conn_id: Optional[str]
    right_conn_id: Optional[str]
    bond_vector: Optional[np.ndarray] = None

    def __post_init__(self):
        """Calculate bond vector if both connection points exist."""
        if self.left_conn_coords is not None and self.right_conn_coords is not None:
            self.bond_vector = self.right_conn_coords - self.left_conn_coords

    def get_center(self) -> np.ndarray:
        """Calculate the geometric center of all atoms."""
        if not self.atoms:
            return np.zeros(3)
        coords = np.array([a.coords for a in self.atoms])
        return np.mean(coords, axis=0)

    def get_all_coords(self) -> np.ndarray:
        """Get array of all atom coordinates."""
        return np.array([a.coords for a in self.atoms])

get_center() -> np.ndarray

Calculate the geometric center of all atoms.

Source code in AutoPoly/mc/chain_growth.py

def get_center(self) -> np.ndarray:
    """Calculate the geometric center of all atoms."""
    if not self.atoms:
        return np.zeros(3)
    coords = np.array([a.coords for a in self.atoms])
    return np.mean(coords, axis=0)

get_all_coords() -> np.ndarray

Get array of all atom coordinates.

Source code in AutoPoly/mc/chain_growth.py

def get_all_coords(self) -> np.ndarray:
    """Get array of all atom coordinates."""
    return np.array([a.coords for a in self.atoms])

MonomerPlacement dataclass

Represents a placed monomer with its transformation.

Attributes:

Name	Type	Description
template	MonomerTemplate	The monomer template used
position	ndarray	Translation vector applied
rotation_matrix	ndarray	3x3 rotation matrix applied
rotation_axis_angle	Tuple[float, float, float, float]	(angle_deg, ax, ay, az) for moltemplate .rot() command
monomer_index	int	Index in the chain
world_coords	ndarray	Atom coordinates after transformation
world_left_conn	Optional[ndarray]	Left connection point after transformation
world_right_conn	Optional[ndarray]	Right connection point after transformation

Source code in AutoPoly/mc/chain_growth.py

@dataclass
class MonomerPlacement:
    """
    Represents a placed monomer with its transformation.

    Attributes:
        template: The monomer template used
        position: Translation vector applied
        rotation_matrix: 3x3 rotation matrix applied
        rotation_axis_angle: (angle_deg, ax, ay, az) for moltemplate .rot() command
        monomer_index: Index in the chain
        world_coords: Atom coordinates after transformation
        world_left_conn: Left connection point after transformation
        world_right_conn: Right connection point after transformation
    """
    template: MonomerTemplate
    position: np.ndarray
    rotation_matrix: np.ndarray
    rotation_axis_angle: Tuple[float, float, float, float]
    monomer_index: int
    world_coords: np.ndarray = field(default_factory=lambda: np.array([]))
    world_left_conn: Optional[np.ndarray] = None
    world_right_conn: Optional[np.ndarray] = None

ChainGrowthMC

Self-Avoiding Random Walk chain growth using Monte Carlo.

This class builds polymer chains by placing monomers one at a time, using collision detection to ensure self-avoiding conformations.

The algorithm: 1. Place first monomer at origin with random orientation 2. For each subsequent monomer: a. Align left connection to previous monomer's right connection b. Sample random dihedral angles c. Check for collisions d. Accept if no collision, retry otherwise

Attributes:

Name	Type	Description
collision_detector		CollisionDetector for checking overlaps
max_attempts		Maximum placement attempts per monomer
templates	Dict[str, MonomerTemplate]	Cache of loaded monomer templates

Methods:

Name	Description
load_monomer_template	Load and cache a monomer template from .lt file.
align_monomer_to_connection	Align monomer so its left connection is at target position.
sample_bond_angle	Sample a bond angle randomly within the configured range.
grow_chain	Build a polymer chain using self-avoiding random walk.
grow_chain_from_templates	Build a polymer chain from pre-built monomer templates.
generate_lt_commands	Generate moltemplate instantiation commands from placements.

Source code in AutoPoly/mc/chain_growth.py

class ChainGrowthMC:
    """
    Self-Avoiding Random Walk chain growth using Monte Carlo.

    This class builds polymer chains by placing monomers one at a time,
    using collision detection to ensure self-avoiding conformations.

    The algorithm:
    1. Place first monomer at origin with random orientation
    2. For each subsequent monomer:
       a. Align left connection to previous monomer's right connection
       b. Sample random dihedral angles
       c. Check for collisions
       d. Accept if no collision, retry otherwise

    Attributes:
        collision_detector: CollisionDetector for checking overlaps
        max_attempts: Maximum placement attempts per monomer
        templates: Cache of loaded monomer templates
    """

    def __init__(
        self,
        collision_detector: CollisionDetector,
        max_attempts: int = 1000,
        bond_angle_min: float = 95.0,
        bond_angle_max: float = 150.0,
        intrachain_exclude_neighbors: int = 2,
        junction_bond_length: float = 1.54
    ):
        """
        Initialize the chain growth Monte Carlo sampler.

        Args:
            collision_detector: CollisionDetector instance for checking overlaps
            max_attempts: Maximum number of placement attempts per monomer
            bond_angle_min: Minimum bond angle in degrees (C-C-C angle between consecutive monomers)
            bond_angle_max: Maximum bond angle in degrees (C-C-C angle between consecutive monomers)
            intrachain_exclude_neighbors: Number of neighbors to exclude from intra-chain collision detection.
                                         0 = exclude only current monomer (maximum collision checking)
                                         1 = exclude i-1, i, i+1 (immediate neighbors)
                                         2 = exclude i-2, i-1, i, i+1, i+2 (default, recommended)
                                         3 = exclude i-3 through i+3
                                         Higher values allow tighter packing but must maintain zero self-intersections.
            junction_bond_length: Equilibrium bond length (Angstrom) of the
                                         inter-monomer bond formed at each junction.
                                         The incoming monomer's left connection is
                                         placed this far beyond the previous right
                                         connection (previously it was placed exactly
                                         ON it, giving zero-length junction bonds).
        """
        self.collision_detector = collision_detector
        self.max_attempts = max_attempts
        self.bond_angle_min = bond_angle_min
        self.bond_angle_max = bond_angle_max
        self.intrachain_exclude_neighbors = intrachain_exclude_neighbors
        self.junction_bond_length = junction_bond_length
        self.templates: Dict[str, MonomerTemplate] = {}

    def load_monomer_template(self, lt_file: str) -> MonomerTemplate:
        """
        Load and cache a monomer template from .lt file.

        Args:
            lt_file: Path to the .lt file

        Returns:
            MonomerTemplate with parsed data
        """
        if lt_file not in self.templates:
            self.templates[lt_file] = parse_lt_file(lt_file)
        return self.templates[lt_file]

    def align_monomer_to_connection(
        self,
        template: MonomerTemplate,
        target_position: np.ndarray,
        incoming_direction: np.ndarray,
        dihedral_angle: float = 0.0
    ) -> Tuple[np.ndarray, np.ndarray]:
        """
        Align monomer so its left connection is at target position.

        Note: In the new bond angle implementation, incoming_direction already
        includes both bond angle bend and dihedral rotation. The dihedral_angle
        parameter is kept for backward compatibility with legacy mode.

        Steps:
        1. Translate monomer so left_conn is at origin
        2. Rotate to align bond_vector with incoming_direction
        3. (Legacy) Apply dihedral rotation around backbone axis
        4. Translate so left_conn is at target_position

        Args:
            template: MonomerTemplate to align
            target_position: Where left_conn should be placed
            incoming_direction: Direction the backbone is coming from
                              (already includes bond angle and dihedral)
            dihedral_angle: (Deprecated) Dihedral angle in radians

        Returns:
            Tuple of (3x3 rotation matrix, translation vector)
        """
        if template.left_conn_coords is None:
            raise ValueError("Cannot align monomer without left connection point")

        # Get the bond vector (direction from left to right connection)
        if template.bond_vector is not None:
            bond_vec = template.bond_vector
        else:
            # Use a default direction if no bond vector
            bond_vec = np.array([1.0, 0.0, 0.0])

        # Normalize incoming direction
        incoming_dir = incoming_direction / np.linalg.norm(incoming_direction)

        # Align bond_vector with incoming_direction
        # (incoming_dir already includes bond angle and dihedral from _apply_bond_angle_bend)
        R_align = rotation_matrix_align_vectors(bond_vec, incoming_dir)

        # Apply additional dihedral rotation only if provided (legacy compatibility)
        if dihedral_angle != 0.0:
            R_dihedral = rotation_matrix_from_axis_angle(incoming_dir, dihedral_angle)
            R_total = R_dihedral @ R_align
        else:
            R_total = R_align

        # Calculate translation
        # After rotation, left_conn moves to: R_total @ left_conn_coords
        rotated_left_conn = R_total @ template.left_conn_coords
        translation = target_position - rotated_left_conn

        return R_total, translation

    def sample_bond_angle(self) -> float:
        """
        Sample a bond angle randomly within the configured range.

        Returns:
            Bond angle in radians (C-C-C angle between consecutive monomers)
        """
        # Uniform random sampling within range
        angle_deg = np.random.uniform(self.bond_angle_min, self.bond_angle_max)
        return np.radians(angle_deg)

    def _apply_bond_angle_bend(
        self,
        prev_bond_dir: np.ndarray,
        bond_angle: float,
        dihedral: float
    ) -> np.ndarray:
        """
        Apply bond angle bend to previous bond direction using spherical coordinates.

        Constructs direction on a cone around prev_bond_dir:
        - Cone half-angle = bond_angle (angle between consecutive bonds)
        - Azimuthal position = dihedral

        Args:
            prev_bond_dir: Normalized direction of previous bond
            bond_angle: C-C-C bond angle in radians (e.g., 109.5° for tetrahedral)
            dihedral: Dihedral rotation angle in radians

        Returns:
            Normalized incoming direction vector for new monomer
        """
        # The angle between consecutive bond vectors equals the bond angle
        # (not the supplement)

        # Find two orthogonal vectors perpendicular to prev_bond_dir
        if abs(prev_bond_dir[2]) < 0.9:
            u = np.cross(prev_bond_dir, np.array([0, 0, 1]))
        else:
            u = np.cross(prev_bond_dir, np.array([1, 0, 0]))
        u = u / np.linalg.norm(u)

        v = np.cross(prev_bond_dir, u)
        v = v / np.linalg.norm(v)

        # Construct direction on cone using spherical coordinates
        incoming_dir = (
            np.cos(bond_angle) * prev_bond_dir +
            np.sin(bond_angle) * (np.cos(dihedral) * u + np.sin(dihedral) * v)
        )

        return incoming_dir / np.linalg.norm(incoming_dir)

    def _transform_coords(
        self,
        coords: np.ndarray,
        rotation: np.ndarray,
        translation: np.ndarray
    ) -> np.ndarray:
        """Apply rotation and translation to coordinates."""
        return (rotation @ coords.T).T + translation

    def _estimate_monomer_radius(self, template: MonomerTemplate) -> float:
        """Estimate collision radius using backbone atoms only (C, O, N, S).

        Using backbone-only atoms produces a smaller radius (~1.0-1.5 Å) than
        the full bounding sphere (~3-5 Å). This keeps the radius below the
        monomer connection distance (~2.35 Å for PEO), so the auto-calibration
        condition ``collision_diameter > conn_dist`` stays false and intra-chain
        collision detection remains active.
        """
        backbone_elements = {'C', 'O', 'N', 'S'}
        backbone_coords = [a.coords for a in template.atoms if a.element in backbone_elements]

        if not backbone_coords:
            # Fallback to all heavy atoms
            backbone_coords = [a.coords for a in template.atoms if a.element != 'H']

        if not backbone_coords:
            return 1.5  # Minimum default

        coords = np.array(backbone_coords)
        center = np.mean(coords, axis=0)
        distances = np.linalg.norm(coords - center, axis=1)
        max_distance = np.max(distances) if len(distances) > 0 else 0

        # Small buffer (0.3 Å) — backbone-only radius should be smaller than
        # the connection distance to keep intra-chain collision active
        radius = max_distance + 0.3
        logger.debug(f"Monomer {template.monomer_name}: backbone-only radius = {radius:.2f} Å")
        return radius

    def grow_chain(
        self,
        monomer_lt_files: List[str],
        chain_id: int = 0,
        start_position: Optional[np.ndarray] = None
    ) -> List[MonomerPlacement]:
        """
        Build a polymer chain using self-avoiding random walk.

        Algorithm:
        1. Place first monomer at start_position with random 3D orientation
        2. For each subsequent monomer:
           - Get right_conn of previous monomer as target for left_conn
           - Sample random dihedral angles
           - Check collision with all placed atoms
           - Accept if no collision, retry otherwise

        Args:
            monomer_lt_files: List of .lt file paths for each monomer in sequence
            chain_id: Unique identifier for this chain
            start_position: Starting position (default: origin)

        Returns:
            List of MonomerPlacement objects for each placed monomer

        Raises:
            RuntimeError: If chain growth fails after max_attempts
        """
        if start_position is None:
            start_position = np.zeros(3)

        templates = [self.load_monomer_template(f) for f in monomer_lt_files]
        return self.grow_chain_from_templates(
            templates, chain_id=chain_id, start_position=start_position
        )

    def grow_chain_from_templates(
        self,
        templates: List[MonomerTemplate],
        chain_id: int = 0,
        start_position: Optional[np.ndarray] = None
    ) -> List[MonomerPlacement]:
        """
        Build a polymer chain from pre-built monomer templates.

        Same algorithm as grow_chain(), but takes MonomerTemplate objects
        directly instead of parsing them from .lt files. This lets the
        force-field-agnostic geometry stage grow chains before any typed
        .lt file exists.

        Args:
            templates: MonomerTemplate for each monomer in sequence
            chain_id: Unique identifier for this chain
            start_position: Starting position (default: origin)

        Returns:
            List of MonomerPlacement objects for each placed monomer

        Raises:
            RuntimeError: If chain growth fails after max_attempts
        """
        if start_position is None:
            start_position = np.zeros(3)

        placements = []

        for i, template in enumerate(templates):
            monomer_radius = self._estimate_monomer_radius(template)

            placement = None

            if i == 0:
                # First monomer: random orientation at start position
                placement = self._place_first_monomer(
                    template, start_position, chain_id, monomer_radius
                )
            else:
                # Subsequent monomers: align to previous
                prev_placement = placements[-1]
                placement = self._place_subsequent_monomer(
                    template, prev_placement, i, chain_id, monomer_radius
                )

            if placement is None:
                raise RuntimeError(
                    f"Failed to place monomer {i} ({template.monomer_name}) "
                    f"after {self.max_attempts} attempts"
                )

            placements.append(placement)

        return placements

    def _place_first_monomer(
        self,
        template: MonomerTemplate,
        position: np.ndarray,
        chain_id: int,
        radius: float
    ) -> Optional[MonomerPlacement]:
        """Place the first monomer with random orientation."""
        for _ in range(self.max_attempts):
            # Random orientation
            R = random_rotation_matrix()

            # Transform coordinates
            coords = template.get_all_coords()
            center = template.get_center()

            # Rotate around center, then translate to position
            rotated_coords = (R @ (coords - center).T).T + position

            # Calculate world connection points
            world_left_conn = None
            world_right_conn = None
            if template.left_conn_coords is not None:
                world_left_conn = R @ (template.left_conn_coords - center) + position
            if template.right_conn_coords is not None:
                world_right_conn = R @ (template.right_conn_coords - center) + position

            # Check collision
            monomer_center = np.mean(rotated_coords, axis=0)
            if not self.collision_detector.check_collision(monomer_center, radius):
                # Add to collision detector
                mon_id = chain_id * 10000  # Unique ID scheme
                self.collision_detector.add_monomer(mon_id, monomer_center, radius)

                axis_angle = rotation_matrix_to_axis_angle(R)

                return MonomerPlacement(
                    template=template,
                    position=position - R @ center,  # Offset for moltemplate
                    rotation_matrix=R,
                    rotation_axis_angle=axis_angle,
                    monomer_index=0,
                    world_coords=rotated_coords,
                    world_left_conn=world_left_conn,
                    world_right_conn=world_right_conn
                )

        return None

    def _place_subsequent_monomer(
        self,
        template: MonomerTemplate,
        prev_placement: MonomerPlacement,
        monomer_index: int,
        chain_id: int,
        radius: float
    ) -> Optional[MonomerPlacement]:
        """Place a subsequent monomer aligned to the previous one."""
        if prev_placement.world_right_conn is None:
            raise ValueError("Previous monomer has no right connection point")

        # Target position for this monomer's left connection
        target_pos = prev_placement.world_right_conn

        # Previous bond direction (to apply bend angle to)
        if prev_placement.world_left_conn is not None:
            prev_bond_dir = prev_placement.world_right_conn - prev_placement.world_left_conn
        else:
            prev_bond_dir = prev_placement.world_right_conn - np.mean(prev_placement.world_coords, axis=0)

        prev_bond_dir = prev_bond_dir / np.linalg.norm(prev_bond_dir)

        # Build exclusion set for selective intra-chain collision detection
        # Exclude only nearby neighbors (±intrachain_exclude_neighbors) - this prevents
        # the chain from folding back on itself while allowing realistic bond angles
        exclude_ids = set()
        for offset in range(-self.intrachain_exclude_neighbors,
                           self.intrachain_exclude_neighbors + 1):
            neighbor_idx = monomer_index + offset
            if 0 <= neighbor_idx < monomer_index:  # Only exclude already-placed monomers
                exclude_ids.add(chain_id * 10000 + neighbor_idx)

        # Always exclude current monomer
        exclude_ids.add(chain_id * 10000 + monomer_index)

        for attempt in range(self.max_attempts):
            # 1. Sample bond angle
            bond_angle = self.sample_bond_angle()

            # 2. Sample dihedral angle — uniform sampling for generic SAW.
            # Uniform dihedral + excluded volume naturally produces the correct
            # Flory exponent (ν ≈ 0.588) for self-avoiding walks.
            dihedral = np.random.uniform(0, 2 * np.pi)

            # 3. Construct incoming direction with bond angle bend
            incoming_dir = self._apply_bond_angle_bend(prev_bond_dir, bond_angle, dihedral)

            # 4. Align monomer (dihedral=0 since already applied in incoming_dir)
            #    The left connection is placed one equilibrium bond length
            #    BEYOND the previous right connection, along the growth
            #    direction — otherwise the junction atoms coincide exactly
            #    and the inter-monomer bond written into the .lt has zero length.
            junction_target = target_pos + self.junction_bond_length * incoming_dir
            R, translation = self.align_monomer_to_connection(
                template, junction_target, incoming_dir, 0.0
            )

            # Transform all coordinates
            coords = template.get_all_coords()
            world_coords = self._transform_coords(coords, R, translation)

            # Calculate world connection points
            world_left_conn = None
            world_right_conn = None
            if template.left_conn_coords is not None:
                world_left_conn = R @ template.left_conn_coords + translation
            if template.right_conn_coords is not None:
                world_right_conn = R @ template.right_conn_coords + translation

            # Check collision
            monomer_center = np.mean(world_coords, axis=0)

            # Check bounds
            if not self.collision_detector.check_bounds(monomer_center, radius):
                continue

            if not self.collision_detector.check_collision(
                monomer_center, radius, exclude_ids=exclude_ids
            ):
                # Add to collision detector
                mon_id = chain_id * 10000 + monomer_index
                self.collision_detector.add_monomer(mon_id, monomer_center, radius)

                axis_angle = rotation_matrix_to_axis_angle(R)

                # Enhanced logging with collision info
                logger.debug(
                    f"Chain {chain_id} monomer {monomer_index}: placed after {attempt + 1} attempts, "
                    f"bond_angle = {np.degrees(bond_angle):.1f}°, dihedral = {np.degrees(dihedral):.1f}°, "
                    f"excluded {len(exclude_ids)} neighbors"
                )

                return MonomerPlacement(
                    template=template,
                    position=translation,
                    rotation_matrix=R,
                    rotation_axis_angle=axis_angle,
                    monomer_index=monomer_index,
                    world_coords=world_coords,
                    world_left_conn=world_left_conn,
                    world_right_conn=world_right_conn
                )

        # Log failure for diagnostic purposes
        logger.debug(
            f"Chain {chain_id} monomer {monomer_index}: FAILED after {self.max_attempts} attempts"
        )
        return None

    def generate_lt_commands(
        self,
        placements: List[MonomerPlacement],
        use_rotation: bool = True
    ) -> List[str]:
        """
        Generate moltemplate instantiation commands from placements.

        Output format:
        monomer[0] = new MonomerA.rot(angle,ax,ay,az).move(x,y,z)

        Args:
            placements: List of MonomerPlacement from grow_chain()
            use_rotation: Whether to include .rot() commands

        Returns:
            List of moltemplate command strings
        """
        commands = []

        for placement in placements:
            monomer_name = placement.template.monomer_name
            idx = placement.monomer_index
            pos = placement.position

            if use_rotation and not np.allclose(placement.rotation_matrix, np.eye(3)):
                angle, ax, ay, az = placement.rotation_axis_angle
                cmd = (
                    f"    monomer[{idx}] = new {monomer_name}"
                    f".rot({angle:.4f},{ax:.4f},{ay:.4f},{az:.4f})"
                    f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
                )
            else:
                cmd = (
                    f"    monomer[{idx}] = new {monomer_name}"
                    f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
                )

            commands.append(cmd)

        return commands

load_monomer_template(lt_file: str) -> MonomerTemplate

Load and cache a monomer template from .lt file.

Parameters:

Name	Type	Description	Default
lt_file	str	Path to the .lt file	required

Returns:

Type	Description
MonomerTemplate	MonomerTemplate with parsed data

Source code in AutoPoly/mc/chain_growth.py

def load_monomer_template(self, lt_file: str) -> MonomerTemplate:
    """
    Load and cache a monomer template from .lt file.

    Args:
        lt_file: Path to the .lt file

    Returns:
        MonomerTemplate with parsed data
    """
    if lt_file not in self.templates:
        self.templates[lt_file] = parse_lt_file(lt_file)
    return self.templates[lt_file]

align_monomer_to_connection(template: MonomerTemplate, target_position: np.ndarray, incoming_direction: np.ndarray, dihedral_angle: float = 0.0) -> Tuple[np.ndarray, np.ndarray]

Align monomer so its left connection is at target position.

Note: In the new bond angle implementation, incoming_direction already includes both bond angle bend and dihedral rotation. The dihedral_angle parameter is kept for backward compatibility with legacy mode.

Steps: 1. Translate monomer so left_conn is at origin 2. Rotate to align bond_vector with incoming_direction 3. (Legacy) Apply dihedral rotation around backbone axis 4. Translate so left_conn is at target_position

Parameters:

Name	Type	Description	Default
template	MonomerTemplate	MonomerTemplate to align	required
target_position	ndarray	Where left_conn should be placed	required
incoming_direction	ndarray	Direction the backbone is coming from (already includes bond angle and dihedral)	required
dihedral_angle	float	(Deprecated) Dihedral angle in radians	0.0

Returns:

Type	Description
Tuple[ndarray, ndarray]	Tuple of (3x3 rotation matrix, translation vector)

Source code in AutoPoly/mc/chain_growth.py

def align_monomer_to_connection(
    self,
    template: MonomerTemplate,
    target_position: np.ndarray,
    incoming_direction: np.ndarray,
    dihedral_angle: float = 0.0
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Align monomer so its left connection is at target position.

    Note: In the new bond angle implementation, incoming_direction already
    includes both bond angle bend and dihedral rotation. The dihedral_angle
    parameter is kept for backward compatibility with legacy mode.

    Steps:
    1. Translate monomer so left_conn is at origin
    2. Rotate to align bond_vector with incoming_direction
    3. (Legacy) Apply dihedral rotation around backbone axis
    4. Translate so left_conn is at target_position

    Args:
        template: MonomerTemplate to align
        target_position: Where left_conn should be placed
        incoming_direction: Direction the backbone is coming from
                          (already includes bond angle and dihedral)
        dihedral_angle: (Deprecated) Dihedral angle in radians

    Returns:
        Tuple of (3x3 rotation matrix, translation vector)
    """
    if template.left_conn_coords is None:
        raise ValueError("Cannot align monomer without left connection point")

    # Get the bond vector (direction from left to right connection)
    if template.bond_vector is not None:
        bond_vec = template.bond_vector
    else:
        # Use a default direction if no bond vector
        bond_vec = np.array([1.0, 0.0, 0.0])

    # Normalize incoming direction
    incoming_dir = incoming_direction / np.linalg.norm(incoming_direction)

    # Align bond_vector with incoming_direction
    # (incoming_dir already includes bond angle and dihedral from _apply_bond_angle_bend)
    R_align = rotation_matrix_align_vectors(bond_vec, incoming_dir)

    # Apply additional dihedral rotation only if provided (legacy compatibility)
    if dihedral_angle != 0.0:
        R_dihedral = rotation_matrix_from_axis_angle(incoming_dir, dihedral_angle)
        R_total = R_dihedral @ R_align
    else:
        R_total = R_align

    # Calculate translation
    # After rotation, left_conn moves to: R_total @ left_conn_coords
    rotated_left_conn = R_total @ template.left_conn_coords
    translation = target_position - rotated_left_conn

    return R_total, translation

sample_bond_angle() -> float

Sample a bond angle randomly within the configured range.

Returns:

Type	Description
float	Bond angle in radians (C-C-C angle between consecutive monomers)

Source code in AutoPoly/mc/chain_growth.py

def sample_bond_angle(self) -> float:
    """
    Sample a bond angle randomly within the configured range.

    Returns:
        Bond angle in radians (C-C-C angle between consecutive monomers)
    """
    # Uniform random sampling within range
    angle_deg = np.random.uniform(self.bond_angle_min, self.bond_angle_max)
    return np.radians(angle_deg)

grow_chain(monomer_lt_files: List[str], chain_id: int = 0, start_position: Optional[np.ndarray] = None) -> List[MonomerPlacement]

Build a polymer chain using self-avoiding random walk.

Algorithm: 1. Place first monomer at start_position with random 3D orientation 2. For each subsequent monomer: - Get right_conn of previous monomer as target for left_conn - Sample random dihedral angles - Check collision with all placed atoms - Accept if no collision, retry otherwise

Parameters:

Name	Type	Description	Default
monomer_lt_files	List[str]	List of .lt file paths for each monomer in sequence	required
chain_id	int	Unique identifier for this chain	0
start_position	Optional[ndarray]	Starting position (default: origin)	None

Returns:

Type	Description
List[MonomerPlacement]	List of MonomerPlacement objects for each placed monomer

Raises:

Type	Description
RuntimeError	If chain growth fails after max_attempts

Source code in AutoPoly/mc/chain_growth.py

def grow_chain(
    self,
    monomer_lt_files: List[str],
    chain_id: int = 0,
    start_position: Optional[np.ndarray] = None
) -> List[MonomerPlacement]:
    """
    Build a polymer chain using self-avoiding random walk.

    Algorithm:
    1. Place first monomer at start_position with random 3D orientation
    2. For each subsequent monomer:
       - Get right_conn of previous monomer as target for left_conn
       - Sample random dihedral angles
       - Check collision with all placed atoms
       - Accept if no collision, retry otherwise

    Args:
        monomer_lt_files: List of .lt file paths for each monomer in sequence
        chain_id: Unique identifier for this chain
        start_position: Starting position (default: origin)

    Returns:
        List of MonomerPlacement objects for each placed monomer

    Raises:
        RuntimeError: If chain growth fails after max_attempts
    """
    if start_position is None:
        start_position = np.zeros(3)

    templates = [self.load_monomer_template(f) for f in monomer_lt_files]
    return self.grow_chain_from_templates(
        templates, chain_id=chain_id, start_position=start_position
    )

grow_chain_from_templates(templates: List[MonomerTemplate], chain_id: int = 0, start_position: Optional[np.ndarray] = None) -> List[MonomerPlacement]

Build a polymer chain from pre-built monomer templates.

Same algorithm as grow_chain(), but takes MonomerTemplate objects directly instead of parsing them from .lt files. This lets the force-field-agnostic geometry stage grow chains before any typed .lt file exists.

Parameters:

Name	Type	Description	Default
templates	List[MonomerTemplate]	MonomerTemplate for each monomer in sequence	required
chain_id	int	Unique identifier for this chain	0
start_position	Optional[ndarray]	Starting position (default: origin)	None

Returns:

Type	Description
List[MonomerPlacement]	List of MonomerPlacement objects for each placed monomer

Raises:

Type	Description
RuntimeError	If chain growth fails after max_attempts

Source code in AutoPoly/mc/chain_growth.py

def grow_chain_from_templates(
    self,
    templates: List[MonomerTemplate],
    chain_id: int = 0,
    start_position: Optional[np.ndarray] = None
) -> List[MonomerPlacement]:
    """
    Build a polymer chain from pre-built monomer templates.

    Same algorithm as grow_chain(), but takes MonomerTemplate objects
    directly instead of parsing them from .lt files. This lets the
    force-field-agnostic geometry stage grow chains before any typed
    .lt file exists.

    Args:
        templates: MonomerTemplate for each monomer in sequence
        chain_id: Unique identifier for this chain
        start_position: Starting position (default: origin)

    Returns:
        List of MonomerPlacement objects for each placed monomer

    Raises:
        RuntimeError: If chain growth fails after max_attempts
    """
    if start_position is None:
        start_position = np.zeros(3)

    placements = []

    for i, template in enumerate(templates):
        monomer_radius = self._estimate_monomer_radius(template)

        placement = None

        if i == 0:
            # First monomer: random orientation at start position
            placement = self._place_first_monomer(
                template, start_position, chain_id, monomer_radius
            )
        else:
            # Subsequent monomers: align to previous
            prev_placement = placements[-1]
            placement = self._place_subsequent_monomer(
                template, prev_placement, i, chain_id, monomer_radius
            )

        if placement is None:
            raise RuntimeError(
                f"Failed to place monomer {i} ({template.monomer_name}) "
                f"after {self.max_attempts} attempts"
            )

        placements.append(placement)

    return placements

generate_lt_commands(placements: List[MonomerPlacement], use_rotation: bool = True) -> List[str]

Generate moltemplate instantiation commands from placements.

Output format: monomer[0] = new MonomerA.rot(angle,ax,ay,az).move(x,y,z)

Parameters:

Name	Type	Description	Default
placements	List[MonomerPlacement]	List of MonomerPlacement from grow_chain()	required
use_rotation	bool	Whether to include .rot() commands	True

Returns:

Type	Description
List[str]	List of moltemplate command strings

Source code in AutoPoly/mc/chain_growth.py

def generate_lt_commands(
    self,
    placements: List[MonomerPlacement],
    use_rotation: bool = True
) -> List[str]:
    """
    Generate moltemplate instantiation commands from placements.

    Output format:
    monomer[0] = new MonomerA.rot(angle,ax,ay,az).move(x,y,z)

    Args:
        placements: List of MonomerPlacement from grow_chain()
        use_rotation: Whether to include .rot() commands

    Returns:
        List of moltemplate command strings
    """
    commands = []

    for placement in placements:
        monomer_name = placement.template.monomer_name
        idx = placement.monomer_index
        pos = placement.position

        if use_rotation and not np.allclose(placement.rotation_matrix, np.eye(3)):
            angle, ax, ay, az = placement.rotation_axis_angle
            cmd = (
                f"    monomer[{idx}] = new {monomer_name}"
                f".rot({angle:.4f},{ax:.4f},{ay:.4f},{az:.4f})"
                f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
            )
        else:
            cmd = (
                f"    monomer[{idx}] = new {monomer_name}"
                f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
            )

        commands.append(cmd)

    return commands

parse_lt_file(lt_file: str) -> MonomerTemplate

Parse a monomer .lt file to extract atom coordinates and connection points.

The .lt file format has atoms in the "Data Atoms" block with the format: $atom:ID $mol:... @atom:TYPE CHARGE X Y Z

Connection point convention: - First atom (C1): Left connection point - Second atom (C2): Right connection point

Parameters:

Name	Type	Description	Default
lt_file	str	Path to the .lt file	required

Returns:

Type	Description
MonomerTemplate	MonomerTemplate with parsed data

Raises:

Type	Description
FileNotFoundError	If the .lt file doesn't exist
ValueError	If the file format is invalid

Source code in AutoPoly/mc/chain_growth.py

def parse_lt_file(lt_file: str) -> MonomerTemplate:
    """
    Parse a monomer .lt file to extract atom coordinates and connection points.

    The .lt file format has atoms in the "Data Atoms" block with the format:
    $atom:ID $mol:... @atom:TYPE CHARGE X Y Z

    Connection point convention:
    - First atom (C1): Left connection point
    - Second atom (C2): Right connection point

    Args:
        lt_file: Path to the .lt file

    Returns:
        MonomerTemplate with parsed data

    Raises:
        FileNotFoundError: If the .lt file doesn't exist
        ValueError: If the file format is invalid
    """
    lt_path = Path(lt_file)
    if not lt_path.exists():
        raise FileNotFoundError(f"Monomer .lt file not found: {lt_file}")

    # Extract monomer name from file
    monomer_name = lt_path.stem

    atoms = []
    in_atoms_block = False

    with open(lt_path, 'r') as f:
        for line in f:
            line_stripped = line.strip()

            # Detect start of Data Atoms block
            if line_stripped == 'write("Data Atoms") {':
                in_atoms_block = True
                continue
            elif line_stripped == '}':
                if in_atoms_block:
                    in_atoms_block = False
                continue

            if in_atoms_block and line_stripped:
                atom = _parse_atom_line(line_stripped)
                if atom:
                    atoms.append(atom)

    if len(atoms) < 2:
        raise ValueError(f"Expected at least 2 atoms in {lt_file}, found {len(atoms)}")

    # Determine monomer type from filename
    monomer_type = _determine_monomer_type(monomer_name)

    # Connection points are first two atoms
    left_conn_coords = atoms[0].coords
    right_conn_coords = atoms[1].coords

    left_conn_id = atoms[0].atom_id
    right_conn_id = atoms[1].atom_id

    # For first/last monomers, one connection is "virtual" (used only for chain bonding)
    if monomer_type == "first":
        # First monomer has no left connection to previous
        left_conn_id = None
        left_conn_coords = None
    elif monomer_type == "last":
        # Last monomer has no right connection to next
        right_conn_id = None
        right_conn_coords = None

    return MonomerTemplate(
        lt_file=str(lt_path),
        monomer_name=monomer_name,
        monomer_type=monomer_type,
        atoms=atoms,
        left_conn_coords=left_conn_coords,
        right_conn_coords=right_conn_coords,
        left_conn_id=left_conn_id,
        right_conn_id=right_conn_id
    )

rotation_matrix_from_axis_angle(axis: np.ndarray, angle_rad: float) -> np.ndarray

Create a 3x3 rotation matrix from axis-angle representation.

Uses Rodrigues' rotation formula.

Parameters:

Name	Type	Description	Default
axis	ndarray	Unit vector defining rotation axis	required
angle_rad	float	Rotation angle in radians	required

Returns:

Type	Description
ndarray	3x3 rotation matrix

Source code in AutoPoly/mc/chain_growth.py

def rotation_matrix_from_axis_angle(axis: np.ndarray, angle_rad: float) -> np.ndarray:
    """
    Create a 3x3 rotation matrix from axis-angle representation.

    Uses Rodrigues' rotation formula.

    Args:
        axis: Unit vector defining rotation axis
        angle_rad: Rotation angle in radians

    Returns:
        3x3 rotation matrix
    """
    axis = axis / np.linalg.norm(axis)
    c = np.cos(angle_rad)
    s = np.sin(angle_rad)
    t = 1 - c

    x, y, z = axis

    return np.array([
        [t*x*x + c,    t*x*y - s*z,  t*x*z + s*y],
        [t*x*y + s*z,  t*y*y + c,    t*y*z - s*x],
        [t*x*z - s*y,  t*y*z + s*x,  t*z*z + c]
    ])

rotation_matrix_align_vectors(v1: np.ndarray, v2: np.ndarray) -> np.ndarray

Create rotation matrix that aligns vector v1 to vector v2.

Parameters:

Name	Type	Description	Default
v1	ndarray	Source vector	required
v2	ndarray	Target vector	required

Returns:

Type	Description
ndarray	3x3 rotation matrix R such that R @ v1 is parallel to v2

Source code in AutoPoly/mc/chain_growth.py

def rotation_matrix_align_vectors(v1: np.ndarray, v2: np.ndarray) -> np.ndarray:
    """
    Create rotation matrix that aligns vector v1 to vector v2.

    Args:
        v1: Source vector
        v2: Target vector

    Returns:
        3x3 rotation matrix R such that R @ v1 is parallel to v2
    """
    v1 = v1 / np.linalg.norm(v1)
    v2 = v2 / np.linalg.norm(v2)

    # Check if vectors are already aligned or opposite
    dot = np.dot(v1, v2)

    # For nearly parallel vectors, add small perturbation for numerical stability
    if dot > 0.9999:
        # Add tiny random perturbation to break symmetry and avoid numerical issues
        perturbation = np.random.randn(3) * 1e-6
        v2_perturbed = v2 + perturbation
        v2_perturbed = v2_perturbed / np.linalg.norm(v2_perturbed)
        dot = np.dot(v1, v2_perturbed)
        if dot > 0.9999:
            return np.eye(3)
        v2 = v2_perturbed
    elif dot < -0.9999:
        # Find a perpendicular axis
        perp = np.array([1, 0, 0]) if abs(v1[0]) < 0.9 else np.array([0, 1, 0])
        axis = np.cross(v1, perp)
        axis = axis / np.linalg.norm(axis)
        return rotation_matrix_from_axis_angle(axis, np.pi)

    # Rotation axis is cross product
    axis = np.cross(v1, v2)
    axis_norm = np.linalg.norm(axis)

    # Handle edge case where cross product is very small
    if axis_norm < 1e-10:
        return np.eye(3)

    axis = axis / axis_norm

    # Rotation angle from dot product
    angle = np.arccos(np.clip(dot, -1, 1))

    return rotation_matrix_from_axis_angle(axis, angle)

rotation_matrix_to_axis_angle(R: np.ndarray) -> Tuple[float, float, float, float]

Convert a 3x3 rotation matrix to axis-angle representation.

Returns format suitable for moltemplate: (angle_degrees, ax, ay, az)

Parameters:

Name	Type	Description	Default
R	ndarray	3x3 rotation matrix	required

Returns:

Type	Description
Tuple[float, float, float, float]	Tuple of (angle_degrees, axis_x, axis_y, axis_z)

Source code in AutoPoly/mc/chain_growth.py

def rotation_matrix_to_axis_angle(R: np.ndarray) -> Tuple[float, float, float, float]:
    """
    Convert a 3x3 rotation matrix to axis-angle representation.

    Returns format suitable for moltemplate: (angle_degrees, ax, ay, az)

    Args:
        R: 3x3 rotation matrix

    Returns:
        Tuple of (angle_degrees, axis_x, axis_y, axis_z)
    """
    # Handle identity matrix
    trace = np.trace(R)
    if np.isclose(trace, 3.0):
        return (0.0, 1.0, 0.0, 0.0)

    # Handle 180 degree rotation
    if np.isclose(trace, -1.0):
        # Find the column with largest diagonal element
        diag = np.diag(R)
        i = np.argmax(diag)
        axis = np.zeros(3)
        axis[i] = 1.0
        return (180.0, axis[0], axis[1], axis[2])

    # General case
    angle = np.arccos(np.clip((trace - 1) / 2, -1, 1))

    # Axis from skew-symmetric part
    axis = np.array([
        R[2, 1] - R[1, 2],
        R[0, 2] - R[2, 0],
        R[1, 0] - R[0, 1]
    ])

    axis_norm = np.linalg.norm(axis)
    if axis_norm < 1e-10:
        return (0.0, 1.0, 0.0, 0.0)

    axis = axis / axis_norm
    angle_deg = np.degrees(angle)

    return (angle_deg, axis[0], axis[1], axis[2])

random_rotation_matrix() -> np.ndarray

Generate a uniformly distributed random rotation matrix.

Uses the algorithm from Graphics Gems III for uniform random rotations based on quaternion sampling.

Returns:

Type	Description
ndarray	3x3 rotation matrix

Source code in AutoPoly/mc/chain_growth.py

def random_rotation_matrix() -> np.ndarray:
    """
    Generate a uniformly distributed random rotation matrix.

    Uses the algorithm from Graphics Gems III for uniform random rotations
    based on quaternion sampling.

    Returns:
        3x3 rotation matrix
    """
    # Generate uniform random quaternion
    u1, u2, u3 = np.random.random(3)

    q = np.array([
        np.sqrt(1 - u1) * np.sin(2 * np.pi * u2),
        np.sqrt(1 - u1) * np.cos(2 * np.pi * u2),
        np.sqrt(u1) * np.sin(2 * np.pi * u3),
        np.sqrt(u1) * np.cos(2 * np.pi * u3)
    ])

    # Convert quaternion to rotation matrix
    q0, q1, q2, q3 = q
    return np.array([
        [1 - 2*(q2**2 + q3**2), 2*(q1*q2 - q0*q3), 2*(q1*q3 + q0*q2)],
        [2*(q1*q2 + q0*q3), 1 - 2*(q1**2 + q3**2), 2*(q2*q3 - q0*q1)],
        [2*(q1*q3 - q0*q2), 2*(q2*q3 + q0*q1), 1 - 2*(q1**2 + q2**2)]
    ])

Molecular placement

AutoPoly.mc.placement

Molecular Placement Module for Monte Carlo Random Placement

This module provides Monte Carlo random placement of polymers and molecules in a simulation box with collision detection.

Key Features: - Random position generation within box bounds - Uniform random orientation using quaternion sampling - Collision-aware sequential placement - Generation of moltemplate .rot().move() commands

Created on 2026-01-29 @author: zwu

Classes:

Name	Description
PolymerPlacement	Represents a placed polymer chain with its transformation.
MoleculePlacement	Represents a placed molecule with its transformation.
MolecularPlacementMC	Monte Carlo random placement of polymers and molecules.

PolymerPlacement dataclass

Represents a placed polymer chain with its transformation.

Attributes:

Name	Type	Description
polymer_id	int	Unique identifier for this polymer
poly_name	str	Name of the polymer (e.g., "poly_1")
position	ndarray	Translation vector for the entire polymer
rotation_matrix	ndarray	3x3 rotation matrix for the entire polymer
rotation_axis_angle	Tuple[float, float, float, float]	(angle_deg, ax, ay, az) for moltemplate
monomer_placements	Optional[List[MonomerPlacement]]	List of MonomerPlacement for chain growth (optional)
center	Optional[ndarray]	Center of mass of the polymer
radius	float	Bounding sphere radius

Source code in AutoPoly/mc/placement.py

@dataclass
class PolymerPlacement:
    """
    Represents a placed polymer chain with its transformation.

    Attributes:
        polymer_id: Unique identifier for this polymer
        poly_name: Name of the polymer (e.g., "poly_1")
        position: Translation vector for the entire polymer
        rotation_matrix: 3x3 rotation matrix for the entire polymer
        rotation_axis_angle: (angle_deg, ax, ay, az) for moltemplate
        monomer_placements: List of MonomerPlacement for chain growth (optional)
        center: Center of mass of the polymer
        radius: Bounding sphere radius
    """
    polymer_id: int
    poly_name: str
    position: np.ndarray
    rotation_matrix: np.ndarray
    rotation_axis_angle: Tuple[float, float, float, float]
    monomer_placements: Optional[List[MonomerPlacement]] = None
    center: Optional[np.ndarray] = None
    radius: float = 0.0

MoleculePlacement dataclass

Represents a placed molecule with its transformation.

Attributes:

Name	Type	Description
molecule_id	int	Unique identifier for this molecule
molecule_name	str	Name of the molecule type (e.g., "water")
instance_name	str	Instance name in moltemplate (e.g., "molecule_1")
position	ndarray	Translation vector
rotation_matrix	ndarray	3x3 rotation matrix
rotation_axis_angle	Tuple[float, float, float, float]	(angle_deg, ax, ay, az) for moltemplate
center	Optional[ndarray]	Center of the molecule after placement
radius	float	Collision radius

Source code in AutoPoly/mc/placement.py

@dataclass
class MoleculePlacement:
    """
    Represents a placed molecule with its transformation.

    Attributes:
        molecule_id: Unique identifier for this molecule
        molecule_name: Name of the molecule type (e.g., "water")
        instance_name: Instance name in moltemplate (e.g., "molecule_1")
        position: Translation vector
        rotation_matrix: 3x3 rotation matrix
        rotation_axis_angle: (angle_deg, ax, ay, az) for moltemplate
        center: Center of the molecule after placement
        radius: Collision radius
    """
    molecule_id: int
    molecule_name: str
    instance_name: str
    position: np.ndarray
    rotation_matrix: np.ndarray
    rotation_axis_angle: Tuple[float, float, float, float]
    center: Optional[np.ndarray] = None
    radius: float = 0.0

MolecularPlacementMC

Monte Carlo random placement of polymers and molecules.

This class handles placing multiple polymers and molecules in a simulation box using random positions and orientations while avoiding collisions.

Attributes:

Name	Type	Description
box_bounds		Simulation box boundaries
collision_detector		CollisionDetector for overlap checking
max_attempts		Maximum placement attempts per molecule

Methods:

Name	Description
random_position	Generate a uniform random position within the box bounds.
random_orientation	Generate a uniform random orientation.
estimate_polymer_radius	Estimate bounding sphere for a polymer chain.
estimate_polymer_radius_from_templates	Estimate polymer radius from template files and spacing.
place_polymer	Place a polymer with random position and orientation.
place_molecule	Place a molecule with random position and orientation.
place_all_polymers	Place multiple polymers sequentially with collision avoidance.
place_all_molecules	Place multiple molecules sequentially with collision avoidance.
generate_polymer_lt_commands	Generate moltemplate instantiation commands for polymers.
generate_molecule_lt_commands	Generate moltemplate instantiation commands for molecules.
get_placement_stats	Get statistics about placed entities.
clear	Reset the placement state.

Source code in AutoPoly/mc/placement.py

class MolecularPlacementMC:
    """
    Monte Carlo random placement of polymers and molecules.

    This class handles placing multiple polymers and molecules in a simulation
    box using random positions and orientations while avoiding collisions.

    Attributes:
        box_bounds: Simulation box boundaries
        collision_detector: CollisionDetector for overlap checking
        max_attempts: Maximum placement attempts per molecule
    """

    def __init__(
        self,
        box_bounds: Tuple[Tuple[float, float], Tuple[float, float], Tuple[float, float]],
        collision_detector: Optional[CollisionDetector] = None,
        max_attempts: int = 10000
    ):
        """
        Initialize the molecular placement MC sampler.

        Args:
            box_bounds: Simulation box as ((xmin, xmax), (ymin, ymax), (zmin, zmax))
            collision_detector: Optional CollisionDetector. If None, creates one.
            max_attempts: Maximum placement attempts per molecule/polymer
        """
        self.box_bounds = box_bounds
        self.max_attempts = max_attempts

        if collision_detector is None:
            # Estimate cell size from box dimensions
            box_size = min(
                box_bounds[0][1] - box_bounds[0][0],
                box_bounds[1][1] - box_bounds[1][0],
                box_bounds[2][1] - box_bounds[2][0]
            )
            cell_size = max(5.0, box_size / 20)
            self.collision_detector = CollisionDetector(box_bounds, cell_size)
        else:
            self.collision_detector = collision_detector

        # Counters for unique IDs
        self._polymer_counter = 0
        self._molecule_counter = 0

    def random_position(self, margin: float = 0.0) -> np.ndarray:
        """
        Generate a uniform random position within the box bounds.

        Args:
            margin: Minimum distance from box walls

        Returns:
            3D position as numpy array
        """
        xmin, xmax = self.box_bounds[0]
        ymin, ymax = self.box_bounds[1]
        zmin, zmax = self.box_bounds[2]

        x = np.random.uniform(xmin + margin, xmax - margin)
        y = np.random.uniform(ymin + margin, ymax - margin)
        z = np.random.uniform(zmin + margin, zmax - margin)

        return np.array([x, y, z])

    def random_orientation(self) -> Tuple[np.ndarray, Tuple[float, float, float, float]]:
        """
        Generate a uniform random orientation.

        Uses quaternion sampling for uniform distribution on SO(3).

        Returns:
            Tuple of (3x3 rotation matrix, (angle_deg, ax, ay, az))
        """
        R = random_rotation_matrix()
        axis_angle = rotation_matrix_to_axis_angle(R)
        return R, axis_angle

    @staticmethod
    def estimate_polymer_radius(
        monomer_placements: List[MonomerPlacement]
    ) -> Tuple[np.ndarray, float]:
        """
        Estimate bounding sphere for a polymer chain.

        Computes the center of mass and maximum distance from center
        to any atom in the chain.

        Args:
            monomer_placements: List of MonomerPlacement from chain growth

        Returns:
            Tuple of (center position, bounding sphere radius)
        """
        if not monomer_placements:
            return np.zeros(3), 0.0

        # Collect all atom coordinates
        all_coords = []
        for placement in monomer_placements:
            all_coords.append(placement.world_coords)

        all_coords = np.vstack(all_coords)
        center = np.mean(all_coords, axis=0)

        # Find maximum distance from center
        distances = np.linalg.norm(all_coords - center, axis=1)
        radius = np.max(distances)

        return center, radius

    @staticmethod
    def estimate_polymer_radius_from_templates(
        lt_files: List[str],
        spacing: float = 3.5
    ) -> float:
        """
        Estimate polymer radius from template files and spacing.

        Uses a simple linear model: radius = N * spacing / 2

        Args:
            lt_files: List of .lt file paths
            spacing: Expected monomer spacing in Angstroms

        Returns:
            Estimated radius in Angstroms
        """
        n_monomers = len(lt_files)
        # Approximate as extended chain divided by 2
        return n_monomers * spacing / 2 + 2.0  # Extra buffer

    def place_polymer(
        self,
        poly_name: str,
        radius: float,
        center: Optional[np.ndarray] = None
    ) -> Optional[PolymerPlacement]:
        """
        Place a polymer with random position and orientation.

        Args:
            poly_name: Name of the polymer (e.g., "poly_1")
            radius: Bounding sphere radius for collision detection
            center: Optional center offset within polymer coordinates

        Returns:
            PolymerPlacement if successful, None if failed after max_attempts
        """
        if center is None:
            center = np.zeros(3)

        for attempt in range(self.max_attempts):
            # Random position with margin for polymer radius
            position = self.random_position(margin=radius)

            # Random orientation
            R, axis_angle = self.random_orientation()

            # World center after transformation
            world_center = R @ center + position

            # Check collision
            if not self.collision_detector.check_collision(world_center, radius):
                # Add to collision detector
                polymer_id = self._polymer_counter
                self._polymer_counter += 1

                self.collision_detector.add_monomer(
                    polymer_id + 100000,  # Offset to avoid collision with monomer IDs
                    world_center,
                    radius
                )

                return PolymerPlacement(
                    polymer_id=polymer_id,
                    poly_name=poly_name,
                    position=position,
                    rotation_matrix=R,
                    rotation_axis_angle=axis_angle,
                    center=world_center,
                    radius=radius
                )

        return None

    def place_molecule(
        self,
        molecule_name: str,
        radius: float,
        instance_name: Optional[str] = None
    ) -> Optional[MoleculePlacement]:
        """
        Place a molecule with random position and orientation.

        Args:
            molecule_name: Type name of the molecule (e.g., "water")
            radius: Collision radius for the molecule
            instance_name: Optional instance name. If None, auto-generated.

        Returns:
            MoleculePlacement if successful, None if failed after max_attempts
        """
        for attempt in range(self.max_attempts):
            # Random position with margin
            position = self.random_position(margin=radius)

            # Random orientation
            R, axis_angle = self.random_orientation()

            # Check collision
            if not self.collision_detector.check_collision(position, radius):
                # Add to collision detector
                molecule_id = self._molecule_counter
                self._molecule_counter += 1

                if instance_name is None:
                    instance_name = f"molecule_{molecule_id + 1}"

                self.collision_detector.add_monomer(
                    molecule_id + 200000,  # Offset to avoid collision with other IDs
                    position,
                    radius
                )

                return MoleculePlacement(
                    molecule_id=molecule_id,
                    molecule_name=molecule_name,
                    instance_name=instance_name,
                    position=position,
                    rotation_matrix=R,
                    rotation_axis_angle=axis_angle,
                    center=position,
                    radius=radius
                )

        return None

    def place_all_polymers(
        self,
        polymer_specs: List[Dict[str, Any]]
    ) -> List[PolymerPlacement]:
        """
        Place multiple polymers sequentially with collision avoidance.

        Args:
            polymer_specs: List of dicts with keys:
                - poly_name: Name of the polymer (e.g., "poly_1")
                - radius: Bounding sphere radius
                - center: Optional center offset (default: origin)

        Returns:
            List of PolymerPlacement objects

        Raises:
            RuntimeError: If any polymer fails to place
        """
        placements = []

        for spec in polymer_specs:
            placement = self.place_polymer(
                poly_name=spec['poly_name'],
                radius=spec['radius'],
                center=spec.get('center')
            )

            if placement is None:
                raise RuntimeError(
                    f"Failed to place polymer {spec['poly_name']} "
                    f"after {self.max_attempts} attempts"
                )

            placements.append(placement)

        return placements

    def place_all_molecules(
        self,
        molecule_specs: List[Dict[str, Any]]
    ) -> List[MoleculePlacement]:
        """
        Place multiple molecules sequentially with collision avoidance.

        Args:
            molecule_specs: List of dicts with keys:
                - molecule_name: Type name of the molecule
                - radius: Collision radius
                - instance_name: Optional instance name

        Returns:
            List of MoleculePlacement objects

        Raises:
            RuntimeError: If any molecule fails to place
        """
        placements = []

        for spec in molecule_specs:
            placement = self.place_molecule(
                molecule_name=spec['molecule_name'],
                radius=spec['radius'],
                instance_name=spec.get('instance_name')
            )

            if placement is None:
                raise RuntimeError(
                    f"Failed to place molecule {spec['molecule_name']} "
                    f"after {self.max_attempts} attempts"
                )

            placements.append(placement)

        return placements

    def generate_polymer_lt_commands(
        self,
        placements: List[PolymerPlacement],
        use_rotation: bool = True
    ) -> List[str]:
        """
        Generate moltemplate instantiation commands for polymers.

        Output format:
        polymer_1 = new poly_1.rot(angle,ax,ay,az).move(x,y,z)

        Args:
            placements: List of PolymerPlacement
            use_rotation: Whether to include .rot() commands

        Returns:
            List of moltemplate command strings
        """
        commands = []

        for placement in placements:
            poly_name = placement.poly_name
            pos = placement.position

            # Derive instance name from poly_name (e.g., "poly_1" -> "polymer_1")
            # This ensures correct naming even when some placements fail
            if poly_name.startswith("poly_"):
                poly_num = poly_name[5:]  # Extract number after "poly_"
                instance_name = f"polymer_{poly_num}"
            else:
                instance_name = f"polymer_{placement.polymer_id + 1}"

            if use_rotation and not np.allclose(placement.rotation_matrix, np.eye(3)):
                angle, ax, ay, az = placement.rotation_axis_angle
                cmd = (
                    f"{instance_name} = new {poly_name}"
                    f".rot({angle:.4f},{ax:.4f},{ay:.4f},{az:.4f})"
                    f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
                )
            else:
                cmd = (
                    f"{instance_name} = new {poly_name}"
                    f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
                )

            commands.append(cmd)

        return commands

    def generate_molecule_lt_commands(
        self,
        placements: List[MoleculePlacement],
        use_rotation: bool = True
    ) -> List[str]:
        """
        Generate moltemplate instantiation commands for molecules.

        Output format:
        molecule_1 = new water.rot(angle,ax,ay,az).move(x,y,z)

        Args:
            placements: List of MoleculePlacement
            use_rotation: Whether to include .rot() commands

        Returns:
            List of moltemplate command strings
        """
        commands = []

        for placement in placements:
            mol_name = placement.molecule_name
            instance = placement.instance_name
            pos = placement.position

            if use_rotation and not np.allclose(placement.rotation_matrix, np.eye(3)):
                angle, ax, ay, az = placement.rotation_axis_angle
                cmd = (
                    f"{instance} = new {mol_name}"
                    f".rot({angle:.4f},{ax:.4f},{ay:.4f},{az:.4f})"
                    f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
                )
            else:
                cmd = (
                    f"{instance} = new {mol_name}"
                    f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
                )

            commands.append(cmd)

        return commands

    def get_placement_stats(self) -> Dict[str, int]:
        """
        Get statistics about placed entities.

        Returns:
            Dict with counts of placed polymers and molecules
        """
        return {
            'polymers': self._polymer_counter,
            'molecules': self._molecule_counter,
            'total_spheres': self.collision_detector.get_sphere_count()
        }

    def clear(self) -> None:
        """Reset the placement state."""
        self.collision_detector.clear()
        self._polymer_counter = 0
        self._molecule_counter = 0

random_position(margin: float = 0.0) -> np.ndarray

Generate a uniform random position within the box bounds.

Parameters:

Name	Type	Description	Default
margin	float	Minimum distance from box walls	0.0

Returns:

Type	Description
ndarray	3D position as numpy array

Source code in AutoPoly/mc/placement.py

def random_position(self, margin: float = 0.0) -> np.ndarray:
    """
    Generate a uniform random position within the box bounds.

    Args:
        margin: Minimum distance from box walls

    Returns:
        3D position as numpy array
    """
    xmin, xmax = self.box_bounds[0]
    ymin, ymax = self.box_bounds[1]
    zmin, zmax = self.box_bounds[2]

    x = np.random.uniform(xmin + margin, xmax - margin)
    y = np.random.uniform(ymin + margin, ymax - margin)
    z = np.random.uniform(zmin + margin, zmax - margin)

    return np.array([x, y, z])

random_orientation() -> Tuple[np.ndarray, Tuple[float, float, float, float]]

Generate a uniform random orientation.

Uses quaternion sampling for uniform distribution on SO(3).

Returns:

Type	Description
Tuple[ndarray, Tuple[float, float, float, float]]	Tuple of (3x3 rotation matrix, (angle_deg, ax, ay, az))

Source code in AutoPoly/mc/placement.py

def random_orientation(self) -> Tuple[np.ndarray, Tuple[float, float, float, float]]:
    """
    Generate a uniform random orientation.

    Uses quaternion sampling for uniform distribution on SO(3).

    Returns:
        Tuple of (3x3 rotation matrix, (angle_deg, ax, ay, az))
    """
    R = random_rotation_matrix()
    axis_angle = rotation_matrix_to_axis_angle(R)
    return R, axis_angle

estimate_polymer_radius(monomer_placements: List[MonomerPlacement]) -> Tuple[np.ndarray, float] staticmethod

Estimate bounding sphere for a polymer chain.

Computes the center of mass and maximum distance from center to any atom in the chain.

Parameters:

Name	Type	Description	Default
monomer_placements	List[MonomerPlacement]	List of MonomerPlacement from chain growth	required

Returns:

Type	Description
Tuple[ndarray, float]	Tuple of (center position, bounding sphere radius)

Source code in AutoPoly/mc/placement.py

@staticmethod
def estimate_polymer_radius(
    monomer_placements: List[MonomerPlacement]
) -> Tuple[np.ndarray, float]:
    """
    Estimate bounding sphere for a polymer chain.

    Computes the center of mass and maximum distance from center
    to any atom in the chain.

    Args:
        monomer_placements: List of MonomerPlacement from chain growth

    Returns:
        Tuple of (center position, bounding sphere radius)
    """
    if not monomer_placements:
        return np.zeros(3), 0.0

    # Collect all atom coordinates
    all_coords = []
    for placement in monomer_placements:
        all_coords.append(placement.world_coords)

    all_coords = np.vstack(all_coords)
    center = np.mean(all_coords, axis=0)

    # Find maximum distance from center
    distances = np.linalg.norm(all_coords - center, axis=1)
    radius = np.max(distances)

    return center, radius

estimate_polymer_radius_from_templates(lt_files: List[str], spacing: float = 3.5) -> float staticmethod

Estimate polymer radius from template files and spacing.

Uses a simple linear model: radius = N * spacing / 2

Parameters:

Name	Type	Description	Default
lt_files	List[str]	List of .lt file paths	required
spacing	float	Expected monomer spacing in Angstroms	3.5

Returns:

Type	Description
float	Estimated radius in Angstroms

Source code in AutoPoly/mc/placement.py

@staticmethod
def estimate_polymer_radius_from_templates(
    lt_files: List[str],
    spacing: float = 3.5
) -> float:
    """
    Estimate polymer radius from template files and spacing.

    Uses a simple linear model: radius = N * spacing / 2

    Args:
        lt_files: List of .lt file paths
        spacing: Expected monomer spacing in Angstroms

    Returns:
        Estimated radius in Angstroms
    """
    n_monomers = len(lt_files)
    # Approximate as extended chain divided by 2
    return n_monomers * spacing / 2 + 2.0  # Extra buffer

place_polymer(poly_name: str, radius: float, center: Optional[np.ndarray] = None) -> Optional[PolymerPlacement]

Place a polymer with random position and orientation.

Parameters:

Name	Type	Description	Default
poly_name	str	Name of the polymer (e.g., "poly_1")	required
radius	float	Bounding sphere radius for collision detection	required
center	Optional[ndarray]	Optional center offset within polymer coordinates	None

Returns:

Type	Description
Optional[PolymerPlacement]	PolymerPlacement if successful, None if failed after max_attempts

Source code in AutoPoly/mc/placement.py

def place_polymer(
    self,
    poly_name: str,
    radius: float,
    center: Optional[np.ndarray] = None
) -> Optional[PolymerPlacement]:
    """
    Place a polymer with random position and orientation.

    Args:
        poly_name: Name of the polymer (e.g., "poly_1")
        radius: Bounding sphere radius for collision detection
        center: Optional center offset within polymer coordinates

    Returns:
        PolymerPlacement if successful, None if failed after max_attempts
    """
    if center is None:
        center = np.zeros(3)

    for attempt in range(self.max_attempts):
        # Random position with margin for polymer radius
        position = self.random_position(margin=radius)

        # Random orientation
        R, axis_angle = self.random_orientation()

        # World center after transformation
        world_center = R @ center + position

        # Check collision
        if not self.collision_detector.check_collision(world_center, radius):
            # Add to collision detector
            polymer_id = self._polymer_counter
            self._polymer_counter += 1

            self.collision_detector.add_monomer(
                polymer_id + 100000,  # Offset to avoid collision with monomer IDs
                world_center,
                radius
            )

            return PolymerPlacement(
                polymer_id=polymer_id,
                poly_name=poly_name,
                position=position,
                rotation_matrix=R,
                rotation_axis_angle=axis_angle,
                center=world_center,
                radius=radius
            )

    return None

place_molecule(molecule_name: str, radius: float, instance_name: Optional[str] = None) -> Optional[MoleculePlacement]

Place a molecule with random position and orientation.

Parameters:

Name	Type	Description	Default
molecule_name	str	Type name of the molecule (e.g., "water")	required
radius	float	Collision radius for the molecule	required
instance_name	Optional[str]	Optional instance name. If None, auto-generated.	None

Returns:

Type	Description
Optional[MoleculePlacement]	MoleculePlacement if successful, None if failed after max_attempts

Source code in AutoPoly/mc/placement.py

def place_molecule(
    self,
    molecule_name: str,
    radius: float,
    instance_name: Optional[str] = None
) -> Optional[MoleculePlacement]:
    """
    Place a molecule with random position and orientation.

    Args:
        molecule_name: Type name of the molecule (e.g., "water")
        radius: Collision radius for the molecule
        instance_name: Optional instance name. If None, auto-generated.

    Returns:
        MoleculePlacement if successful, None if failed after max_attempts
    """
    for attempt in range(self.max_attempts):
        # Random position with margin
        position = self.random_position(margin=radius)

        # Random orientation
        R, axis_angle = self.random_orientation()

        # Check collision
        if not self.collision_detector.check_collision(position, radius):
            # Add to collision detector
            molecule_id = self._molecule_counter
            self._molecule_counter += 1

            if instance_name is None:
                instance_name = f"molecule_{molecule_id + 1}"

            self.collision_detector.add_monomer(
                molecule_id + 200000,  # Offset to avoid collision with other IDs
                position,
                radius
            )

            return MoleculePlacement(
                molecule_id=molecule_id,
                molecule_name=molecule_name,
                instance_name=instance_name,
                position=position,
                rotation_matrix=R,
                rotation_axis_angle=axis_angle,
                center=position,
                radius=radius
            )

    return None

place_all_polymers(polymer_specs: List[Dict[str, Any]]) -> List[PolymerPlacement]

Place multiple polymers sequentially with collision avoidance.

Parameters:

Name	Type	Description	Default
polymer_specs	List[Dict[str, Any]]	List of dicts with keys: - poly_name: Name of the polymer (e.g., "poly_1") - radius: Bounding sphere radius - center: Optional center offset (default: origin)	required

Returns:

Type	Description
List[PolymerPlacement]	List of PolymerPlacement objects

Raises:

Type	Description
RuntimeError	If any polymer fails to place

Source code in AutoPoly/mc/placement.py

def place_all_polymers(
    self,
    polymer_specs: List[Dict[str, Any]]
) -> List[PolymerPlacement]:
    """
    Place multiple polymers sequentially with collision avoidance.

    Args:
        polymer_specs: List of dicts with keys:
            - poly_name: Name of the polymer (e.g., "poly_1")
            - radius: Bounding sphere radius
            - center: Optional center offset (default: origin)

    Returns:
        List of PolymerPlacement objects

    Raises:
        RuntimeError: If any polymer fails to place
    """
    placements = []

    for spec in polymer_specs:
        placement = self.place_polymer(
            poly_name=spec['poly_name'],
            radius=spec['radius'],
            center=spec.get('center')
        )

        if placement is None:
            raise RuntimeError(
                f"Failed to place polymer {spec['poly_name']} "
                f"after {self.max_attempts} attempts"
            )

        placements.append(placement)

    return placements

place_all_molecules(molecule_specs: List[Dict[str, Any]]) -> List[MoleculePlacement]

Place multiple molecules sequentially with collision avoidance.

Parameters:

Name	Type	Description	Default
molecule_specs	List[Dict[str, Any]]	List of dicts with keys: - molecule_name: Type name of the molecule - radius: Collision radius - instance_name: Optional instance name	required

Returns:

Type	Description
List[MoleculePlacement]	List of MoleculePlacement objects

Raises:

Type	Description
RuntimeError	If any molecule fails to place

Source code in AutoPoly/mc/placement.py

def place_all_molecules(
    self,
    molecule_specs: List[Dict[str, Any]]
) -> List[MoleculePlacement]:
    """
    Place multiple molecules sequentially with collision avoidance.

    Args:
        molecule_specs: List of dicts with keys:
            - molecule_name: Type name of the molecule
            - radius: Collision radius
            - instance_name: Optional instance name

    Returns:
        List of MoleculePlacement objects

    Raises:
        RuntimeError: If any molecule fails to place
    """
    placements = []

    for spec in molecule_specs:
        placement = self.place_molecule(
            molecule_name=spec['molecule_name'],
            radius=spec['radius'],
            instance_name=spec.get('instance_name')
        )

        if placement is None:
            raise RuntimeError(
                f"Failed to place molecule {spec['molecule_name']} "
                f"after {self.max_attempts} attempts"
            )

        placements.append(placement)

    return placements

generate_polymer_lt_commands(placements: List[PolymerPlacement], use_rotation: bool = True) -> List[str]

Generate moltemplate instantiation commands for polymers.

Output format: polymer_1 = new poly_1.rot(angle,ax,ay,az).move(x,y,z)

Parameters:

Name	Type	Description	Default
placements	List[PolymerPlacement]	List of PolymerPlacement	required
use_rotation	bool	Whether to include .rot() commands	True

Returns:

Type	Description
List[str]	List of moltemplate command strings

Source code in AutoPoly/mc/placement.py

def generate_polymer_lt_commands(
    self,
    placements: List[PolymerPlacement],
    use_rotation: bool = True
) -> List[str]:
    """
    Generate moltemplate instantiation commands for polymers.

    Output format:
    polymer_1 = new poly_1.rot(angle,ax,ay,az).move(x,y,z)

    Args:
        placements: List of PolymerPlacement
        use_rotation: Whether to include .rot() commands

    Returns:
        List of moltemplate command strings
    """
    commands = []

    for placement in placements:
        poly_name = placement.poly_name
        pos = placement.position

        # Derive instance name from poly_name (e.g., "poly_1" -> "polymer_1")
        # This ensures correct naming even when some placements fail
        if poly_name.startswith("poly_"):
            poly_num = poly_name[5:]  # Extract number after "poly_"
            instance_name = f"polymer_{poly_num}"
        else:
            instance_name = f"polymer_{placement.polymer_id + 1}"

        if use_rotation and not np.allclose(placement.rotation_matrix, np.eye(3)):
            angle, ax, ay, az = placement.rotation_axis_angle
            cmd = (
                f"{instance_name} = new {poly_name}"
                f".rot({angle:.4f},{ax:.4f},{ay:.4f},{az:.4f})"
                f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
            )
        else:
            cmd = (
                f"{instance_name} = new {poly_name}"
                f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
            )

        commands.append(cmd)

    return commands

generate_molecule_lt_commands(placements: List[MoleculePlacement], use_rotation: bool = True) -> List[str]

Generate moltemplate instantiation commands for molecules.

Output format: molecule_1 = new water.rot(angle,ax,ay,az).move(x,y,z)

Parameters:

Name	Type	Description	Default
placements	List[MoleculePlacement]	List of MoleculePlacement	required
use_rotation	bool	Whether to include .rot() commands	True

Returns:

Type	Description
List[str]	List of moltemplate command strings

Source code in AutoPoly/mc/placement.py

def generate_molecule_lt_commands(
    self,
    placements: List[MoleculePlacement],
    use_rotation: bool = True
) -> List[str]:
    """
    Generate moltemplate instantiation commands for molecules.

    Output format:
    molecule_1 = new water.rot(angle,ax,ay,az).move(x,y,z)

    Args:
        placements: List of MoleculePlacement
        use_rotation: Whether to include .rot() commands

    Returns:
        List of moltemplate command strings
    """
    commands = []

    for placement in placements:
        mol_name = placement.molecule_name
        instance = placement.instance_name
        pos = placement.position

        if use_rotation and not np.allclose(placement.rotation_matrix, np.eye(3)):
            angle, ax, ay, az = placement.rotation_axis_angle
            cmd = (
                f"{instance} = new {mol_name}"
                f".rot({angle:.4f},{ax:.4f},{ay:.4f},{az:.4f})"
                f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
            )
        else:
            cmd = (
                f"{instance} = new {mol_name}"
                f".move({pos[0]:.4f},{pos[1]:.4f},{pos[2]:.4f})"
            )

        commands.append(cmd)

    return commands

get_placement_stats() -> Dict[str, int]

Get statistics about placed entities.

Returns:

Type	Description
Dict[str, int]	Dict with counts of placed polymers and molecules

Source code in AutoPoly/mc/placement.py

def get_placement_stats(self) -> Dict[str, int]:
    """
    Get statistics about placed entities.

    Returns:
        Dict with counts of placed polymers and molecules
    """
    return {
        'polymers': self._polymer_counter,
        'molecules': self._molecule_counter,
        'total_spheres': self.collision_detector.get_sphere_count()
    }

clear() -> None

Reset the placement state.

Source code in AutoPoly/mc/placement.py

def clear(self) -> None:
    """Reset the placement state."""
    self.collision_detector.clear()
    self._polymer_counter = 0
    self._molecule_counter = 0