Living crystals
Active colloids of radius 1 with phoretic attraction, which gather into clusters that slowly coarsen.
- Name:
living_crystals - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Argun et al. 2021; Palacci et al. 2013
Factory
@register_model(
field='physics/active_matter',
tags=('active matter', 'sde', 'clusters', 'phoresis', 'pair interaction'),
references=('Argun et al. 2021', 'Palacci et al. 2013'),
defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.002)},
)
def living_crystals(
n_particles: int = 100,
*,
size: float = 50.0,
speed: float = 3.0,
D_r: float = 1.0,
v0: float = 20.0,
cutoff: float = 10.0,
k: float = 200.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return active colloids of radius 1 with phoretic attraction, which gather into clusters that slowly coarsen.
The colloids swim as active Brownian particles, drift towards each other at the speed `v0 / r^2` within the
cutoff, and a stiff soft repulsion stands in for the hard-sphere correction at contact.
Args:
n_particles:
The number of colloids.
size:
The side of the periodic box.
speed:
The self-propulsion speed.
D_r:
The rotational diffusion coefficient.
v0:
The phoretic drift speed at the distance of one radius.
cutoff:
The distance beyond which the attraction is neglected.
k:
The stiffness of the repulsion times the mobility.
random_state:
The seed or generator of the initial positions and orientations.
Returns:
system:
The model, with the positions `swimmers.x` and orientations `swimmers.theta`.
Examples:
>>> from skcomplex.models.physics.active_matter import living_crystals
>>> living_crystals(16).interactions[2].v0
20.0
"""
system = _active_disks(n_particles, size, speed, D_r, random_state)
repulsion = SoftSphereInteraction('contact', first='swimmers', second='swimmers', cutoff=2.0, period=size, k=k)
attraction = PhoreticAttractionInteraction(
'phoresis',
first='swimmers',
second='swimmers',
cutoff=cutoff,
period=size,
v0=v0,
R=1.0,
)
interactions = [*system.interactions, repulsion, attraction]
return System(elements=system.elements, interactions=interactions, space=system.space)
Components
class ActiveParticleElement(Element):
"""Active particle in the plane, whose position diffuses, `dx = sqrt(2 D_t) dW`, and whose orientation turns.
The orientation turns at the angular velocity `omega` and diffuses, `dtheta = omega dt + sqrt(2 D_r) dW`. A
`SelfPropulsionInteraction` adds the drift `v (cos theta, sin theta)` of the position. With `omega = 0` it is an
active Brownian particle, and otherwise a chiral active particle.
Args:
label:
The label of the particle or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The initial position.
theta:
The initial orientation.
D_t:
The translational diffusion coefficient.
D_r:
The rotational diffusion coefficient.
omega:
The angular velocity of a chiral particle.
Examples:
>>> from skcomplex.models.physics.active_matter import ActiveParticleElement
>>> ActiveParticleElement('swimmer', D_r=2.0).D_r
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = POINTS,
*,
x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
theta: Annotated[float, StateSpec(period=2 * pi, role='direction')] = 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,
omega: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.theta = theta
self.D_t = D_t
self.D_r = D_r
self.omega = omega
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the translational noise, and the turning and the rotational noise.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [
NoiseLaw(s.x, sqrt(2 * s.D_t)),
RateLaw(s.theta, s.omega),
NoiseLaw(s.theta, sqrt(2 * s.D_r)),
]
class SelfPropulsionInteraction(RateInteraction):
"""Self-propulsion `x' += speed (cos theta, sin theta)` of a position in the plane along a heading `theta`.
The heading is a state of the same target, which other laws turn and diffuse.
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 self-propulsion speed.
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.models.physics.active_matter import SelfPropulsionInteraction
>>> propulsion = SelfPropulsionInteraction('propulsion', target='swimmers').check().prepare()
>>> target = SimpleNamespace(x=vector([var('x1'), var('x2')]), theta=var('theta'))
>>> [law.value for law in propulsion.run(SimpleNamespace(target=target, speed=var('v'))).laws_]
[(v*cos(theta), v*sin(theta))]
"""
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)
self.speed = speed
self.heading = heading
def _get_rate(self, namespace: SimpleNamespace) -> FreeModuleElement:
"""Return the velocity `speed (cos theta, sin theta)`.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the position.
"""
heading = getattr(namespace.target, self.heading)
return namespace.speed * vector(SR, [cos(heading), sin(heading)])
class SoftSphereInteraction(PairCouplingInteraction):
"""Soft repulsion between overlapping spheres, `k * (cutoff - r) ** 2 / 2` closer than the cutoff.
Args:
label:
The label of the interaction.
first:
The path of the group or particle of the first members.
second:
The path of the group or particle of the second members.
cutoff:
The contact distance, the sum of the two radii.
position:
The name of the vector state that holds the positions of the members.
write:
The name of the state whose rate the force changes, the position of order 2 by default.
divisor:
The name of the parameter that holds the masses of the members, or `None` for unit masses.
period:
The side of a periodic box, for minimum image distances, or `None` for the period of the space.
k:
The stiffness.
Examples:
>>> from skcomplex.models.physics.common.particles import SoftSphereInteraction
>>> SoftSphereInteraction('contact', first='gas', second='gas', cutoff=0.1, k=50.0).k
50.0
"""
def __init__(
self,
label: str,
first: str,
second: str,
cutoff: float = 1.0,
position: str = 'x',
*,
write: str = 'x',
divisor: str | None = None,
period: float | None = None,
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
) -> None:
super().__init__(
label,
first,
second,
cutoff=cutoff,
position=position,
write=write,
potential=self._get_soft_sphere_energy,
divisor=divisor,
period=period,
)
self.k = k
@staticmethod
def _get_soft_sphere_energy(s: SimpleNamespace, r: Expression) -> Expression:
"""Return `k * (cutoff - r) ** 2 / 2`."""
return s.k * (s.cutoff - r) ** 2 / 2
class PhoreticAttractionInteraction(PairCouplingInteraction):
"""Phoretic attraction between active colloids, the pair potential `-v0 R^2 / r` of an overdamped position.
Every member of a pair drifts towards the other at the speed `v0 R^2 / r^2`.
Args:
label:
The label of the interaction, unique among the parts of its system.
first:
The path of the group or particle of the first members.
second:
The path of the group or particle of the second members.
cutoff:
The distance beyond which the attraction is neglected.
position:
The name of the vector state that holds the positions of the members.
v0:
The drift speed at the distance `R`.
R:
The particle radius.
period:
The side of a periodic box, for minimum image distances, or `None` for the period of the space.
Examples:
>>> from skcomplex.models.physics.active_matter import PhoreticAttractionInteraction
>>> PhoreticAttractionInteraction('phoresis', first='colloids', second='colloids', v0=20.0).v0
20.0
"""
def __init__(
self,
label: str,
first: str,
second: str,
cutoff: float = 1.0,
position: str = 'x',
*,
v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 20.0,
R: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
period: float | None = None,
) -> None:
super().__init__(
label,
first,
second,
cutoff=cutoff,
position=position,
write=position,
potential=self._attraction,
period=period,
)
self.v0 = v0
self.R = R
@staticmethod
def _attraction(namespace: SimpleNamespace, r: Expression) -> Expression:
"""Return the potential `-v0 R^2 / r` of a pair at the distance `r`."""
return -namespace.v0 * namespace.R**2 / r
Equations
swimmers.theta[i]: diff(theta_swimmers(t), t) == omega_swimmers