Skip to content

Rigid Brownian and active dimers

Colloidal molecules: pairs of Brownian or active spheres joined by rigid rods, free in the plane.

  • Name: rigid_dimers
  • Field: physics/soft_matter
  • Formalism: stochastic differential equations
  • References: Argun et al. 2021

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('rigid_dimers')
entry = describe('rigid_dimers')
trajectory = create_simulation('rigid_dimers').prepare().run().trajectory_

Factory

@register_system(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'colloidal molecule', 'rotational diffusion', 'active matter'),
    references=('Argun et al. 2021',),
    defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def rigid_dimers(
    n_dimers: int = 100,
    *,
    L: float = 1.0,
    v0: float = 0.0,
    D_t: float = 0.1,
    D_r: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return colloidal molecules: pairs of Brownian or active spheres joined by rigid rods, free in the plane.

    Each sphere diffuses with `D_t`, and an active sphere also swims at the speed `v0` along its own orientation,
    which diffuses with `D_r`. The rod keeps the spheres at the distance `L`. Diffusing both spheres freely and then
    restoring their distance, the scheme of Problem 5.1 of Argun et al., gives in the limit of small steps a centre
    that diffuses with `D_t / 2` and a rod angle that diffuses with `2 D_t / L²`, so the rotational diffusion of the
    molecule falls as `L^-2` in this free-draining model. Active spheres, Problem 9.3, swim together along the mean of
    their orientations and turn the rod with the part of their difference across it, which raises the long-time
    diffusion of the centre to `D_t / 2 + v0² / (4 D_r)`. The dimers start at the origin with random rod angles and
    orientations.

    Args:
        n_dimers:
            The number of dimers.

        L:
            The length of the rod.

        v0:
            The swimming speed of each sphere, 0 for Brownian spheres.

        D_t:
            The translational diffusion coefficient of each sphere.

        D_r:
            The rotational diffusion coefficient of the orientation of each sphere.

        random_state:
            The seed or generator of the initial angles and orientations.

    Returns:
        system:
            The system, with the centres `dimers.c`, the rod angles `dimers.phi`, and the orientations `dimers.theta1`
            and `dimers.theta2`. The spheres sit at `c ± (L / 2) (cos phi, sin phi)`.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import rigid_dimers
        >>> rigid_dimers(5, L=2.0).elements[0].values['L']
        2.0
    """
    rng = np.random.default_rng(random_state)
    angles = rng.uniform(0.0, 2 * np.pi, (3, n_dimers))
    values = {
        'c': np.zeros((n_dimers, 2)),
        'phi': angles[0],
        'theta1': angles[1],
        'theta2': angles[2],
        'L': L,
        'v0': v0,
        'D_t': D_t,
        'D_r': D_r,
    }
    dimers = GroupElement('dimers', element=DimerElement, size=n_dimers, values=values)
    return System(elements=[dimers], space=ContinuousSpace(n_dim=2))

Components

Elements:

  • DimerElement, from Element
class DimerElement(Element):
    """Two Brownian or active spheres joined by a rigid rod, described by the centre and the angle of the rod.

    Each sphere diffuses with `D_t` and swims at the speed `v0` along its own orientation, `theta1` or `theta2`, which
    diffuse with `D_r`, and the rod keeps them at the distance `L`. Moving both spheres freely and then restoring
    their distance moves the centre `c` by the mean of their displacements and turns the rod by the part of their
    relative displacement across it, so in the limit of small steps `dc = v0 (n1 + n2) / 2 dt + sqrt(D_t) dW` and
    `dphi = v0 (n1 - n2) . e / L dt + 2 sqrt(D_t) / L dW`, with `e = (-sin phi, cos phi)` across the rod.

    Args:
        label:
            The label of the dimer or the group.

        visual:
            How the dimer is drawn, or `None` for `default_visual_`.

        c:
            The initial centre.

        phi:
            The initial angle of the rod.

        theta1:
            The initial orientation of the first sphere.

        theta2:
            The initial orientation of the second sphere.

        L:
            The length of the rod.

        v0:
            The swimming speed of each sphere, 0 for Brownian spheres.

        D_t:
            The translational diffusion coefficient of each sphere.

        D_r:
            The rotational diffusion coefficient of the orientation of each sphere.

    Attributes:
        default_visual_:
            The visual of the dimer when it sets none.

    Examples:
        >>> from skcomplex.systems.physics.soft_matter import DimerElement
        >>> DimerElement('molecule').L
        1.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        c: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
        phi: Annotated[float, StateSpec(role='direction')] = 0.0,
        theta1: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
        theta2: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
        L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        D_t: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.c = c
        self.phi = phi
        self.theta1 = theta1
        self.theta2 = theta2
        self.L = L
        self.v0 = v0
        self.D_t = D_t
        self.D_r = D_r

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the motion of the centre and the turning of the rod, and the diffusion of the orientations.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        first = vector(SR, [cos(s.theta1), sin(s.theta1)])
        second = vector(SR, [cos(s.theta2), sin(s.theta2)])
        across = vector(SR, [-sin(s.phi), cos(s.phi)])
        return [
            RateLaw(s.c, s.v0 * (first + second) / 2),
            NoiseLaw(s.c, sqrt(s.D_t)),
            RateLaw(s.phi, s.v0 * (first - second).dot_product(across) / s.L),
            NoiseLaw(s.phi, 2 * sqrt(s.D_t) / s.L),
            NoiseLaw(s.theta1, sqrt(2 * s.D_r)),
            NoiseLaw(s.theta2, sqrt(2 * s.D_r)),
        ]

Equations

dimers.c[i, 0]: diff(c_dimers_0(t), t) == 1/2*v0_dimers*(cos(theta1_dimers(t)) + cos(theta2_dimers(t)))
dimers.c[i, 1]: diff(c_dimers_1(t), t) == 1/2*v0_dimers*(sin(theta1_dimers(t)) + sin(theta2_dimers(t)))
dimers.phi[i]: diff(phi_dimers(t), t) == ((sin(theta1_dimers(t)) - sin(theta2_dimers(t)))*cos(phi_dimers(t)) - (cos(theta1_dimers(t)) - cos(theta2_dimers(t)))*sin(phi_dimers(t)))*v0_dimers/L_dimers
dimers.theta1[i]: diff(theta1_dimers(t), t) == 0
dimers.theta2[i]: diff(theta2_dimers(t), t) == 0