Brownian motion of Janus, ellipsoidal, and chiral particles in three dimensions
Overdamped axisymmetric Brownian particles in three dimensions that translate and rotate.
- Name:
anisotropic_brownian_particles - Field:
physics/soft_matter - Formalism: stochastic differential equations
- References: Perrin 1934; Han et al. 2006; Argun et al. 2021
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('anisotropic_brownian_particles')
entry = describe('anisotropic_brownian_particles')
trajectory = create_simulation('anisotropic_brownian_particles').prepare().run().trajectory_
Factory
@register_system(
field='physics/soft_matter',
tags=('brownian motion', 'sde', 'rotational diffusion', 'three dimensions', 'anisotropy', 'chirality'),
references=('Perrin 1934', 'Han et al. 2006', 'Argun et al. 2021'),
defaults={'until': 10.0, 'dt': 0.05, 'solver': EulerMaruyamaSolver(dt=0.001)},
)
def anisotropic_brownian_particles(
n_particles: int = 200,
*,
D_par: float = 1.0,
D_perp: float = 1.0,
D_r: float = 0.75,
D_spin: float = 0.75,
D_c: float = 0.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return overdamped axisymmetric Brownian particles in three dimensions that translate and rotate.
Every particle diffuses with `D_par` along its axis `u` and with `D_perp` across it, its axis diffuses on the unit
sphere with `D_r`, so that `<u(t) . u(0)> = exp(-2 D_r t)`, and it spins about its axis with `D_spin`. The defaults
are a Janus sphere of unit radius, `D_par = D_perp = D_t` and `D_r = 3 D_t / 4`, whose orientation is visible
although its diffusion is isotropic. A prolate ellipsoid has `D_par > D_perp` from Perrin's factors, so its
displacement over times shorter than `1 / D_r` is elongated along its axis. A chiral body, such as a spiral, has a
coupling `D_c` of the translation along its axis and the rotation about it, `<dr_par dpsi> = 2 D_c dt`. These are
Challenges 5.1 to 5.3 of Argun et al. The axis is held as a unit vector and rescaled to unit length after every
step, and the particles start at the origin with axes drawn uniformly on the sphere.
Args:
n_particles:
The number of particles.
D_par:
The translational diffusion coefficient along the axis.
D_perp:
The translational diffusion coefficient across the axis.
D_r:
The rotational diffusion coefficient of the axis.
D_spin:
The rotational diffusion coefficient about the axis.
D_c:
The coupling of the translation along the axis and the rotation about it.
random_state:
The seed or generator of the initial axes.
Returns:
system:
The system, with the positions `particles.x`, `particles.y`, `particles.z`, the axes `particles.ux`,
`particles.uy`, `particles.uz`, and the spins `particles.psi`.
Raises:
ParameterError:
The coupling exceeds `sqrt(D_par D_spin)`.
Examples:
>>> from skcomplex.systems.physics.soft_matter import anisotropic_brownian_particles
>>> anisotropic_brownian_particles(5, D_par=1.5, D_perp=0.8).elements[0].values['D_perp']
0.8
"""
if D_c**2 > D_par * D_spin:
error_msg = f'The coupling {D_c!r} exceeds sqrt(D_par D_spin), so the diffusion tensor is not positive.'
raise ParameterError(error_msg)
axes = np.random.default_rng(random_state).normal(size=(n_particles, 3))
axes /= np.linalg.norm(axes, axis=1, keepdims=True)
values = {
'ux': axes[:, 0],
'uy': axes[:, 1],
'uz': axes[:, 2],
'D_par': D_par,
'D_perp': D_perp,
'D_r': D_r,
'D_spin': D_spin,
'D_c': D_c,
}
particles = GroupElement('particles', element=AnisotropicParticleElement, size=n_particles, values=values)
return System(elements=[particles])
Components
Elements:
AnisotropicParticleElement, fromElement
class AnisotropicParticleElement(Element):
"""Axisymmetric Brownian particle in three dimensions, whose diffusion tensor turns with its axis.
The position `r = (x, y, z)` diffuses with `D_par` along the unit axis `u = (ux, uy, uz)` and with `D_perp` across
it, `dr = (sqrt(2 D_par) u u + sqrt(2 D_perp) (I - u u)) dW`. The axis diffuses on the unit sphere,
`du = -2 D_r u dt + sqrt(2 D_r) (I - u u) dW_r` in the Itô form, and a boundary rescales it to unit length after
every step against the drift of the discrete steps. The spin `psi` about the axis diffuses with `D_spin`, and its
noise is correlated with the displacement along the axis through the coupling `D_c` of a chiral body, so that
`<dr_par dpsi> = 2 D_c dt`, which needs `D_c² <= D_par D_spin`.
Args:
label:
The label of the particle or the group.
visual:
How the particle is drawn, or `None` for `default_visual_`.
x:
The initial first coordinate.
y:
The initial second coordinate.
z:
The initial third coordinate.
ux:
The initial first component of the axis.
uy:
The initial second component of the axis.
uz:
The initial third component of the axis.
psi:
The initial spin angle about the axis.
D_par:
The translational diffusion coefficient along the axis.
D_perp:
The translational diffusion coefficient across the axis.
D_r:
The rotational diffusion coefficient of the axis.
D_spin:
The rotational diffusion coefficient about the axis.
D_c:
The coupling of the translation along the axis and the rotation about it.
Attributes:
default_visual_:
The visual of the particle when it sets none.
Examples:
>>> from skcomplex.systems.physics.soft_matter import AnisotropicParticleElement
>>> AnisotropicParticleElement('janus').uz
1.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec()] = 0.0,
y: Annotated[float, StateSpec()] = 0.0,
z: Annotated[float, StateSpec()] = 0.0,
ux: Annotated[float, StateSpec()] = 0.0,
uy: Annotated[float, StateSpec()] = 0.0,
uz: Annotated[float, StateSpec()] = 1.0,
psi: Annotated[float, StateSpec()] = 0.0,
D_par: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
D_perp: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
D_r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
D_spin: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.75,
D_c: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.y = y
self.z = z
self.ux = ux
self.uy = uy
self.uz = uz
self.psi = psi
self.D_par = D_par
self.D_perp = D_perp
self.D_r = D_r
self.D_spin = D_spin
self.D_c = D_c
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the translation and the rotation with their correlated noises, and the boundary of the axis.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
axis = [getattr(s, name) for name in DIRECTIONS]
position = [getattr(s, name) for name in AXES]
length = sqrt(sum(component**2 for component in axis))
laws: list[BaseLaw] = []
for i, (coordinate, component) in enumerate(zip(position, axis, strict=True)):
laws.append(RateLaw(component, -2 * s.D_r * component))
for j, other in enumerate(axis):
projector = (1 if i == j else 0) - component * other
laws.append(NoiseLaw(coordinate, sqrt(2 * s.D_perp) * projector, process=f'translation_{j}'))
laws.append(NoiseLaw(coordinate, sqrt(2 * s.D_par) * component * other, process=f'translation_{j}'))
laws.append(NoiseLaw(component, sqrt(2 * s.D_r) * projector, process=f'rotation_{j}'))
laws.append(BoundaryLaw(component, component / length))
coupling = s.D_c / sqrt(s.D_par)
for j, other in enumerate(axis):
laws.append(NoiseLaw(s.psi, sqrt(2) * coupling * other, process=f'translation_{j}'))
laws.append(NoiseLaw(s.psi, sqrt(2 * (s.D_spin - coupling**2)), process='spin'))
return laws
Equations
particles.x[i]: diff(x_particles(t), t) == 0
particles.y[i]: diff(y_particles(t), t) == 0
particles.z[i]: diff(z_particles(t), t) == 0
particles.ux[i]: diff(ux_particles(t), t) == -2*D_r_particles*ux_particles(t)
particles.uy[i]: diff(uy_particles(t), t) == -2*D_r_particles*uy_particles(t)
particles.uz[i]: diff(uz_particles(t), t) == -2*D_r_particles*uz_particles(t)
particles.psi[i]: diff(psi_particles(t), t) == 0