Inertial and overdamped Brownian motion
Free Brownian particles with inertia, m x'' = -gamma x' + sqrt(2 kT gamma) W, or without it.
- Name:
brownian_motion - Field:
physics/soft_matter - Formalism: stochastic differential equations
- References: Argun et al. 2021; Li et al. 2010; Volpe and Volpe 2013
Factory
@register_model(
field='physics/soft_matter',
tags=('brownian motion', 'sde', 'diffusion', 'inertia'),
references=('Argun et al. 2021', 'Li et al. 2010', 'Volpe and Volpe 2013'),
defaults={'until': 10.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def brownian_motion(
n_particles: int = 500,
*,
inertial: bool = True,
m: float = 1.0,
gamma: float = 1.0,
kT: float = 1.0,
n_dim: int = 1,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return free Brownian particles with inertia, `m x'' = -gamma x' + sqrt(2 kT gamma) W`, or without it.
With inertia the motion is ballistic, `<x^2> = kT t^2 / m`, below the momentum relaxation time `m / gamma`, and
diffusive, `<x^2> = 2 kT t / gamma` per axis, above it. Without inertia it is diffusive at every time.
Args:
n_particles:
The number of particles.
inertial:
Whether the particles have inertia.
m:
The mass, used with inertia.
gamma:
The friction coefficient.
kT:
The thermal energy.
n_dim:
The number of dimensions.
random_state:
The seed or generator of the thermal initial velocities.
Returns:
system:
The model, with the positions `beads.x`.
Examples:
>>> from skcomplex.models.physics.soft_matter import brownian_motion
>>> brownian_motion(10, inertial=False).elements[0].element.__name__
'BrownianParticleElement'
"""
space = ContinuousSpace(n_dim=n_dim)
if not inertial:
beads = GroupElement(
'beads',
element=BrownianParticleElement,
size=n_particles,
values={'x': np.zeros((n_particles, n_dim)), 'gamma': gamma, 'temperature': kT},
)
return System(elements=[beads], interactions=[ThermalNoiseInteraction('bath', target='beads')], space=space)
velocities = np.random.default_rng(random_state).normal(0.0, math.sqrt(kT / m), (n_particles, n_dim))
beads = GroupElement(
'beads',
element=MassPointElement,
size=n_particles,
values={'x': np.zeros((n_particles, n_dim)), 'x_d1': velocities, 'm': m},
visual=POINTS,
)
bath = LangevinThermostatInteraction('bath', target='beads', temperature=kT, gamma=gamma)
return System(elements=[beads], interactions=[bath], space=space)
Components
class MassPointElement(PointElement):
"""Point mass, a point whose position `x` has order 2, with its velocity `x_d1` and a mass `m`.
The forces on the point divide by the mass through their divisor `'m'`.
Args:
label:
The label of the point.
x:
The initial position. A number applies to every component.
x_d1:
The initial velocity. A number applies to every component.
visual:
How the point is drawn, or `None` for `default_visual_`.
m:
The mass.
Attributes:
default_visual_:
Points at the positions.
Examples:
>>> from skcomplex.models.physics.common.particles import MassPointElement
>>> MassPointElement('ball', x=[0.0, 1.0], m=2.0).m
2.0
"""
def __init__(
self,
label: str,
x: npt.ArrayLike = 0.0,
x_d1: npt.ArrayLike = 0.0,
visual: VisualSpec | None = None,
*,
m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, x=x, x_d1=x_d1, visual=visual)
self.m = m
class LangevinThermostatInteraction(RelaxationInteraction):
"""Langevin thermostat, which adds a drag `-gamma * v` and a random force that hold particles at a temperature.
Every velocity component follows `dv = -gamma * v / m * dt + sqrt(2 * gamma * temperature) / m * dW`, with the
temperature in energy units, a relaxation of the velocity towards rest over the time `m / gamma` with additive
noise. The velocity is the first derivative of the position, a state of order 2, and the mass is the divisor. A
system with a thermostat is stochastic, and its velocities relax to the Maxwell-Boltzmann distribution at that
temperature.
Args:
label:
The label of the interaction.
target:
The path of the particle or particle group.
temperature:
The temperature in energy units.
gamma:
The friction coefficient.
state:
The name of the vector state of order 2 that holds the positions of the particles.
divisor:
The name of the state or parameter that holds the masses of the particles.
Examples:
>>> from skcomplex.models.physics.common.particles import LangevinThermostatInteraction
>>> LangevinThermostatInteraction('bath', target='gas', temperature=2.0).temperature
2.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
temperature: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
state: str = 'x',
divisor: str = 'm',
) -> None:
super().__init__(label, target=target, state=state)
self.temperature = temperature
self.gamma = gamma
self.divisor = divisor
def _get_time(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation time `1 / gamma` of the velocity, before the division by the mass.
Args:
namespace:
The namespace of the symbols the time can use.
Returns:
time:
The relaxation time.
"""
return 1 / namespace.gamma
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the amplitude `sqrt(2 * gamma * temperature)` of the noise, before the division by the mass.
Args:
namespace:
The namespace of the symbols the amplitude can use.
Returns:
amplitude:
The noise amplitude.
"""
return sqrt(2 * namespace.gamma * namespace.temperature)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation `-v / time` of the velocity `v`, the first derivative of the position.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The relaxation of the velocity towards rest, before the division by the mass.
"""
return -velocity_of(getattr(namespace.target, self.state)) / self._get_time(namespace)
Equations
beads.x[i, 0]: diff(x_beads_0(t), t, t) == -gamma_bath*diff(x_beads_0(t), t)/m_beads