Active polymers
An active polymer, a Rouse chain of active Brownian beads that each swim along their own orientation.
- Name:
active_polymer - Field:
physics/active_matter - Formalism: stochastic differential equations
- References: Winkler and Gompper 2020; Osmanović and Rabin 2017; Eisenstecken et al. 2016
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('active_polymer')
entry = describe('active_polymer')
trajectory = create_simulation('active_polymer').prepare().run().trajectory_
Factory
@register_system(
field='physics/active_matter',
tags=('active matter', 'sde', 'polymer', 'self-propulsion', 'rouse chain'),
references=('Winkler and Gompper 2020', 'Osmanović and Rabin 2017', 'Eisenstecken et al. 2016'),
defaults={'until': 50.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def active_polymer(
n_beads: int = 20,
*,
v0: float = 1.0,
k: float = 10.0,
D_t: float = 0.1,
D_r: float = 1.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return an active polymer, a Rouse chain of active Brownian beads that each swim along their own orientation.
Bead `i` follows `r_i' = v0 n_i + k (r_{i+1} + r_{i-1} - 2 r_i) + sqrt(2 D_t) xi_i`, with its orientation
`n_i = (cos theta_i, sin theta_i)` diffusing with `D_r`, and the end beads bonded to one neighbour. The bonds cancel
in the motion of the centre of mass, whose mean squared displacement is
`(4 D_t t + 2 v0² (D_r t - 1 + exp(-D_r t)) / D_r²) / N`, that of an active Brownian particle divided by `N`. The
activity swells the chain: a passive chain has the mean squared end-to-end distance `2 D_t (N - 1) / k` in the
plane, and self-propulsion of persistence time `1 / D_r` adds to the bond fluctuations and stretches it further. The
beads start on a straight line with unit spacing and random orientations.
Args:
n_beads:
The number of beads.
v0:
The self-propulsion speed of every bead.
k:
The stiffness of the bonds.
D_t:
The translational diffusion coefficient.
D_r:
The rotational diffusion coefficient.
random_state:
The seed or generator of the initial orientations.
Returns:
system:
The system, with the coordinates `chain.x` and `chain.y` and the orientations `chain.theta` of the beads.
Examples:
>>> from skcomplex.systems.physics.active_matter import active_polymer
>>> active_polymer(5).elements[0].prepare().size_
5
"""
angles = np.random.default_rng(random_state).uniform(0.0, 2 * np.pi, n_beads)
chain = FieldElement(
'chain',
element=ActiveBeadElement,
space=LatticeSpace(n_beads, boundary='reflecting'),
boundary={'x': 'zero-flux', 'y': 'zero-flux'},
values={
'x': np.arange(n_beads, dtype=float),
'y': np.zeros(n_beads),
'theta': angles,
'v0': v0,
'k': k,
'D_t': D_t,
'D_r': D_r,
},
)
return System(elements=[chain])
Components
Elements:
ActiveBeadElement, fromElement
class ActiveBeadElement(Element):
"""Bead of an active polymer, an active Brownian particle bonded to its neighbours along a chain by springs.
The bead at `(x, y)` swims at the speed `v0` along its orientation `theta`, which diffuses with `D_r`, its position
diffuses with `D_t`, and harmonic springs of stiffness `k` and zero rest length pull it towards its neighbours
along the chain, `k (r_{i+1} + r_{i-1} - 2 r_i)`, the Rouse chain. The beads at the ends have one neighbour.
Args:
label:
The label of the bead or the chain.
visual:
How the bead is drawn, or `None` for `default_visual_`.
x:
The initial first coordinate.
y:
The initial second coordinate.
theta:
The initial orientation.
v0:
The self-propulsion speed.
k:
The stiffness of the bonds.
D_t:
The translational diffusion coefficient.
D_r:
The rotational diffusion coefficient.
Attributes:
default_visual_:
The visual of the bead when it sets none.
Examples:
>>> from skcomplex.systems.physics.active_matter import ActiveBeadElement
>>> ActiveBeadElement('bead').k
10.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec()] = 0.0,
y: Annotated[float, StateSpec()] = 0.0,
theta: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
v0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 10.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.x = x
self.y = y
self.theta = theta
self.v0 = v0
self.k = k
self.D_t = D_t
self.D_r = D_r
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the propulsion and the bonds of the bead, and the noises of its position and orientation.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [
RateLaw(s.x, s.v0 * cos(s.theta) + s.k * laplacian(s.x)),
RateLaw(s.y, s.v0 * sin(s.theta) + s.k * laplacian(s.y)),
NoiseLaw(s.x, sqrt(2 * s.D_t)),
NoiseLaw(s.y, sqrt(2 * s.D_t)),
NoiseLaw(s.theta, sqrt(2 * s.D_r)),
]
Equations
chain.x[i]: diff(x_chain(t), t) == v0_chain*cos(theta_chain(t)) + k_chain*laplacian(x_chain(t))
chain.y[i]: diff(y_chain(t), t) == k_chain*laplacian(y_chain(t)) + v0_chain*sin(theta_chain(t))
chain.theta[i]: diff(theta_chain(t), t) == 0