Skip to content

Ginelli nematic particles

Polar particles with nematic alignment, which order along a common axis but move both ways along it.

  • Name: ginelli_nematics
  • Field: physics/active_matter
  • Formalism: discrete map
  • References: Ginelli et al. 2010; Chaté 2020

Select this system

from skcomplex.systems import create_simulation, create_system, describe

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

Factory

@register_system(
    field='physics/active_matter',
    tags=('active matter', 'flocking', 'nematic', 'discrete time', 'neighborhood'),
    references=('Ginelli et al. 2010', 'Chaté 2020'),
    defaults={'until': 300.0, 'dt': 1.0},
    observables=lambda system: _orders('flock'),
)
def ginelli_nematics(
    n_particles: int = 400,
    *,
    size: float = 10.0,
    v0: float = 0.3,
    eta: float = 0.08,
    seed: int = 0,
) -> System:
    """Return polar particles with nematic alignment, which order along a common axis but move both ways along it.

    At every step particle `i` takes the direction of `sum_j sign(cos(theta_j - theta_i)) exp(i theta_j)` over the
    particles within the unit distance, itself included, so it aligns with the axis of every neighbour on the side
    nearer its own heading, plus a uniform angle in `[-eta pi, eta pi]`, and moves `v0` along it. At low noise the
    particles share an axis, with as many moving one way as the other: the nematic order parameter
    `|<exp(2 i theta)>|` grows towards 1 while the polar order parameter stays small.

    Args:
        n_particles:
            The number of particles.

        size:
            The side of the periodic box, in units of the interaction radius.

        v0:
            The distance a particle moves in one step.

        eta:
            The amplitude of the noise, as a fraction of `pi`.

        seed:
            The seed of the uniform initial positions and headings.

    Returns:
        system:
            The system, with the positions `flock.x`, the headings `flock.theta`, and the observables `polarization`
            and `nematic_order`.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import ginelli_nematics
        >>> ginelli_nematics(10).elements[0].values['eta']
        0.08
    """
    alignment = HeadingAlignmentInteraction('alignment', target='flock', nematic=True)
    return _flock(NematicParticleElement, n_particles, size, seed, v0=v0, alignment=alignment, eta=eta)

Components

Elements:

  • NematicParticleElement, from Element

Interactions:

  • PropulsionStepInteraction, from MapInteraction
  • HeadingAlignmentInteraction, from CouplingInteraction
class NematicParticleElement(Element):
    """Polar particle with nematic alignment, which aligns with the axis of every neighbour on its nearer side.

    `HeadingAlignmentInteraction` with `nematic=True` sets the new heading to the circular mean of the headings within
    the unit distance, itself included, each flipped to the side nearer the particle's own heading, plus a uniform
    angle in `[-eta pi, eta pi]`. A propulsion step moves the particle along its heading.

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

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

        x:
            The initial position.

        theta:
            The initial heading.

        eta:
            The amplitude of the noise, as a fraction of `pi`.

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

    Examples:
        >>> from skcomplex.systems.physics.active_matter import NematicParticleElement
        >>> NematicParticleElement('rod', eta=0.2).eta
        0.2
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
        theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 0.0,
        eta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    ) -> None:
        super().__init__(label, visual=visual)
        self.x = x
        self.theta = theta
        self.eta = eta

class PropulsionStepInteraction(MapInteraction):
    """Step `x(t + 1) = x + speed (cos theta, sin theta)` of a position in the plane along a heading, in discrete time.

    The step reads the heading at the start of the step, as every map of the step does.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the particle or the group whose positions move.

        speed:
            The distance moved in one step.

        state:
            The name of the vector state that holds the position.

        heading:
            The name of the state that holds the heading.

    Examples:
        >>> from types import SimpleNamespace
        >>> from sage.all__sagemath_symbolics import var, vector
        >>> from skcomplex.systems.physics.common.particles import PropulsionStepInteraction
        >>> step = PropulsionStepInteraction('propulsion', target='flock').prepare()
        >>> target = SimpleNamespace(x=vector([var('x1'), var('x2')]), theta=var('theta'))
        >>> [law.value for law in step.run(SimpleNamespace(target=target, speed=var('v'))).laws_]
        [(v*cos(theta) + x1, v*sin(theta) + x2)]
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        speed: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
        state: str = 'x',
        heading: str = 'theta',
    ) -> None:
        super().__init__(label, target=target, state=state, function=self._step, mode='increment')
        self.speed = speed
        self.heading = heading

    def _step(self, namespace: SimpleNamespace) -> FreeModuleElement:
        """Return the step `speed (cos theta, sin theta)`."""
        heading = getattr(namespace.target, self.heading)
        return namespace.speed * vector(SR, [cos(heading), sin(heading)])

class HeadingAlignmentInteraction(CouplingInteraction):
    """Vicsek alignment, the turn of every particle to the circular mean heading of its neighbours plus noise.

    Every step the heading `theta` of every particle becomes the circular mean of the headings within the unit
    distance, itself included, plus a uniform angle in `[-eta pi, eta pi]` with the parameter `eta` of the particle.
    With `nematic`, every heading is flipped to the side nearer the particle's own heading first, as for particles
    with an axis rather than a direction. It is the map form of coupling with the circular mean over a radius. Another
    neighbourhood, such as the `k` nearest particles or a vision cone, replaces the unit disk.

    With a positive `latency` `h`, every particle reads the headings of its neighbours `h` steps ago. With a negative
    latency it reads their headings predicted `|h|` steps ahead, by a least-squares line through their last `memory`
    headings, unwrapped by taking every change between two steps as the smallest angle. A latent alignment needs
    particles that keep enough past headings, such as `DelayedVicsekParticleElement`.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the group of particles.

        nematic:
            Whether every heading is flipped to the side nearer the particle's own heading.

        neighborhood:
            The neighbourhood of every particle, itself included, or `None` for the unit disk.

        latency:
            The number of steps by which the headings of the neighbours are read late, or predicted ahead when
            negative.

        memory:
            The number of past headings a prediction fits, at least 2.

    Examples:
        >>> from skcomplex.systems.physics.active_matter import HeadingAlignmentInteraction
        >>> HeadingAlignmentInteraction('alignment', target='flock', nematic=True).reduce
        'circular-mean'
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        nematic: bool = False,
        neighborhood: BaseNeighborhood | None = None,
        latency: int = 0,
        memory: int = 5,
    ) -> None:
        super().__init__(
            label,
            target=target,
            read='theta',
            write='theta',
            function=self._get_heading,
            reduce='circular-mean',
            weights=self._get_side,
            form='map',
            combine=self._get_noisy,
            neighborhood=RadiusNeighborhood(1.0, center=True) if neighborhood is None else neighborhood,
        )
        self.nematic = nematic
        self.latency = latency
        self.memory = memory

    def _get_heading(self, namespace: SimpleNamespace) -> Expression:
        """Return the heading of the neighbour, delayed or predicted."""
        heading = namespace.other.theta
        if self.latency >= 0:
            return lag(heading, self.latency)
        unwrapped = [heading]
        for step in range(1, self.memory):
            change = lag(heading, step - 1) - lag(heading, step)
            unwrapped.append(unwrapped[-1] - atan2(sin(change), cos(change)))
        times = [-step for step in range(self.memory)]
        mean_time = sum(times) / self.memory
        mean_heading = sum(unwrapped) / self.memory
        spread = sum((time - mean_time) ** 2 for time in times)
        slope = sum((time - mean_time) * (value - mean_heading) for time, value in zip(times, unwrapped, strict=True))
        return mean_heading + slope / spread * (-self.latency - mean_time)

    def _get_side(self, namespace: SimpleNamespace) -> Expression:
        """Return -1 for a neighbour heading away from the particle when nematic, and 1 otherwise."""
        if not self.nematic:
            return SR(1)
        return where(cos(namespace.other.theta - namespace.member.theta) >= 0, 1, -1)

    @staticmethod
    def _get_noisy(namespace: SimpleNamespace, reduction: Target) -> Target:
        """Return the mean heading plus a uniform angle in `[-eta pi, eta pi]`."""
        eta = namespace.member.eta
        return reduction + uniform('noise', -eta * math.pi, eta * math.pi)

Equations

flock.x[i, 0]: x_flock_0(t + 1) == speed_propulsion*cos(theta_flock(t)) + x_flock_0(t)
flock.x[i, 1]: x_flock_1(t + 1) == speed_propulsion*sin(theta_flock(t)) + x_flock_1(t)
flock.theta[i]: theta_flock(t + 1) == 6.283185307179586*eta_flock*uniform_noise_alignment - 3.141592653589793*eta_flock + arctan2(_neighbours_1, _neighbours_0)