Skip to content

Soft matter

Soft matter: Brownian particles, optical traps, anomalous diffusion, and noise-driven transitions.

BrownianParticleElement(label, visual=POINTS, *, x=0.0, gamma=1.0, temperature=1.0)

Bases: Element

Overdamped Brownian particle at the position x, with the friction gamma and the thermal energy temperature.

The particle moves by dx = F / gamma dt + sqrt(2 temperature / gamma) dW, with the thermal noise of ThermalNoiseInteraction and the forces of other interactions, which divide by gamma.

Parameters:

Name Type Description Default
label str

The label of the particle or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 POINTS
x Annotated[ArrayLike, StateSpec(vector=True, role=position)]

The initial position.

 0.0
gamma Annotated[float, ParameterSpec(domain=open(0, oo))]

The friction coefficient.

 1.0
temperature Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The thermal energy.

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import BrownianParticleElement
>>> BrownianParticleElement('bead', temperature=2.0).temperature
2.0
Source code in src/skcomplex/models/physics/soft_matter/_brownian_particle.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.gamma = gamma
    self.temperature = temperature

ColoredNoiseParticleElement(label, visual=None, *, v=0.0, xi=0.0, gamma=1.0, D=1.0, tau=1.0)

Bases: Element

Particle with the velocity v, driven by the noise xi of correlation time tau.

Parameters:

Name Type Description Default
label str

The label of the particle.

 required
visual VisualSpec | None

How the particle is drawn, or None for none.

 None
v Annotated[float, StateSpec()]

The initial velocity.

 0.0
xi Annotated[float, StateSpec()]

The initial noise.

 0.0
gamma Annotated[float, ParameterSpec(domain=open(0, oo))]

The friction coefficient.

 1.0
D Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The noise intensity.

 1.0
tau Annotated[float, ParameterSpec(domain=open(0, oo))]

The correlation time of the noise.

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import ColoredNoiseParticleElement
>>> ColoredNoiseParticleElement('particle', tau=0.5).tau
0.5
Source code in src/skcomplex/models/physics/soft_matter/_colored_noise.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec()] = 0.0,
    xi: Annotated[float, StateSpec()] = 0.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.xi = xi
    self.gamma = gamma
    self.D = D
    self.tau = tau

ContinuousTimeWalkerElement(label, visual=POINTS, *, x=0.0, alpha=0.5, tau0=1.0, sigma=1.0)

Bases: Element

Walker that waits a power-law time tau0 r^(-1 / alpha) for r uniform, then jumps by a Gaussian step.

The jumps are a jump interaction of the model, which resets the position after every waiting time.

Parameters:

Name Type Description Default
label str

The label of the walker or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 POINTS
x Annotated[ArrayLike, StateSpec(vector=True, role=position)]

The initial position.

 0.0
alpha Annotated[float, ParameterSpec(domain=open(0, oo))]

The tail exponent of the waiting times, below 1 for waiting times without a mean.

 0.5
tau0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The shortest waiting time.

 1.0
sigma Annotated[float, ParameterSpec(domain=open(0, oo))]

The standard deviation of the jumps along every axis.

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import ContinuousTimeWalkerElement
>>> ContinuousTimeWalkerElement('walker', alpha=0.7).alpha
0.7
Source code in src/skcomplex/models/physics/soft_matter/_continuous_time_walker.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    tau0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.alpha = alpha
    self.tau0 = tau0
    self.sigma = sigma

DoubleWellElement(label, visual=None, *, x=-1.0, amplitude=0.1, frequency=0.01, D=0.12)

Bases: Element

Particle at the position x in the double well x^4 / 4 - x^2 / 2, forced periodically and by noise.

Parameters:

Name Type Description Default
label str

The label of the particle.

 required
visual VisualSpec | None

How the particle is drawn, or None for none.

 None
x Annotated[float, StateSpec(role=position)]

The initial position.

 -1.0
amplitude Annotated[float, ParameterSpec()]

The amplitude of the periodic forcing.

 0.1
frequency Annotated[float, ParameterSpec()]

The frequency of the forcing.

 0.01
D Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The noise intensity.

 0.12

Examples:

>>> from skcomplex.models.physics.soft_matter import DoubleWellElement
>>> DoubleWellElement('particle', D=0.2).D
0.2
Source code in src/skcomplex/models/physics/soft_matter/_double_well.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(role='position')] = -1.0,
    amplitude: Annotated[float, ParameterSpec()] = 0.1,
    frequency: Annotated[float, ParameterSpec()] = 0.01,
    D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.12,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.amplitude = amplitude
    self.frequency = frequency
    self.D = D

FlashingSawtoothInteraction(label, *, target, barrier=1.0, L=1.0, asymmetry=0.1, switches=())

Bases: RateInteraction

Sawtooth potential on overdamped particles on a line, turned on and off at given times.

Its force is the rate of the position, divided by the friction gamma of the target.

The potential of period L has a minimum at every multiple of L, rises linearly by barrier over the fraction asymmetry of the period, and falls back over the rest. It is on from time 0 to the first of the switches, off until the next one, and so on, and it keeps its state after the last one.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
target Annotated[str, RoleSpec()]

The path of the particle or the group.

 required
barrier Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The height of the sawtooth.

 1.0
L Annotated[float, ParameterSpec(domain=open(0, oo))]

The period of the sawtooth.

 1.0
asymmetry Annotated[float, ParameterSpec(domain=open(0, oo))]

The fraction of the period over which the potential rises, 0.5 for a symmetric potential.

 0.1
switches tuple[float, ...]

The increasing times at which the potential turns off and on in turn.

 ()

Examples:

>>> from skcomplex.models.physics.soft_matter import FlashingSawtoothInteraction
>>> FlashingSawtoothInteraction('ratchet', target='beads', switches=(1.0, 2.5)).switches
(1.0, 2.5)
Source code in src/skcomplex/models/physics/soft_matter/_flashing_sawtooth.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    barrier: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    asymmetry: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    switches: tuple[float, ...] = (),
) -> None:
    super().__init__(label, target=target, state='x', divisor='gamma')
    self.barrier = barrier
    self.L = L
    self.asymmetry = asymmetry
    self.switches = switches

FractionalWalkerElement(label, visual=POINTS, *, x=0.0, D=1.0, H=0.3)

Bases: Element

Walker driven by fractional Brownian motion, dx = sqrt(2 D) dB_H, whose increments are correlated in time.

Parameters:

Name Type Description Default
label str

The label of the walker or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 POINTS
x Annotated[ArrayLike, StateSpec(vector=True, role=position)]

The initial position.

 0.0
D Annotated[float, ParameterSpec(domain=open(0, oo))]

The diffusion coefficient.

 1.0
H Annotated[float, ParameterSpec(domain=open(0, 1))]

The Hurst index of the driving fractional Brownian motion.

 0.3

Examples:

>>> from skcomplex.models.physics.soft_matter import FractionalWalkerElement
>>> FractionalWalkerElement('walker', H=0.3).H
0.3
Source code in src/skcomplex/models/physics/soft_matter/_fractional_walker.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    H: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.D = D
    self.H = H

FrictionalBlockElement(label, visual=None, *, x=0.0, v=0.0, m=1.0, k=1.0, anchor=0.0, mu_s=1.0, mu_d=0.8, amplitude=2.0, shake_time=2.0, relax_time=8.0, v_stick=0.02)

Bases: Element

Block of a chain on a rough floor, at the position x with the velocity v, joined to its neighbours.

The springs of stiffness k and an anchor spring pull the block, Coulomb friction holds it while the force stays below mu_s and slows it by mu_d while it slides, and periodic shaking of amplitude amplitude kicks it during the first shake_time of every shake_time + relax_time.

Parameters:

Name Type Description Default
label str

The label of the block.

 required
visual VisualSpec | None

How the block is drawn, or None for none.

 None
x Annotated[float, StateSpec()]

The initial position.

 0.0
v Annotated[float, StateSpec()]

The initial velocity.

 0.0
m Annotated[float, ParameterSpec(domain=open(0, oo))]

The mass.

 1.0
k Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The stiffness of the springs between neighbours.

 1.0
anchor Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The stiffness of the spring to the rest position.

 0.0
mu_s Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The static friction force.

 1.0
mu_d Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The dynamic friction force.

 0.8
amplitude Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The amplitude of the shaking force.

 2.0
shake_time Annotated[float, ParameterSpec(domain=open(0, oo))]

The duration of every shaking.

 2.0
relax_time Annotated[float, ParameterSpec(domain=open(0, oo))]

The time between shakings.

 8.0
v_stick Annotated[float, ParameterSpec(domain=open(0, oo))]

The speed below which the block can stick.

 0.02

Examples:

>>> from skcomplex.models.physics.soft_matter import FrictionalBlockElement
>>> FrictionalBlockElement('block', mu_s=2.0).mu_s
2.0
Source code in src/skcomplex/models/physics/soft_matter/_friction.py 
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    v: Annotated[float, StateSpec()] = 0.0,
    m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    anchor: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    mu_s: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    mu_d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.8,
    amplitude: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    shake_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    relax_time: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 8.0,
    v_stick: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.02,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.m = m
    self.k = k
    self.anchor = anchor
    self.mu_s = mu_s
    self.mu_d = mu_d
    self.amplitude = amplitude
    self.shake_time = shake_time
    self.relax_time = relax_time
    self.v_stick = v_stick

LevyFlierElement(label, visual=POINTS, *, x=0.0, mu=1.2, l0=1.0)

Bases: Element

Walker that jumps at every step by a length of power-law tail l^(-1 - mu) in a random direction.

Parameters:

Name Type Description Default
label str

The label of the walker or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 POINTS
x Annotated[ArrayLike, StateSpec(vector=True, role=position)]

The initial position.

 0.0
mu Annotated[float, ParameterSpec(domain=open(0, oo))]

The tail exponent of the jump lengths, below 2 for jumps of infinite variance.

 1.2
l0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The shortest jump length.

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import LevyFlierElement
>>> LevyFlierElement('walker', mu=1.2).mu
1.2
Source code in src/skcomplex/models/physics/soft_matter/_levy_flier.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.2,
    l0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.mu = mu
    self.l0 = l0

LevyWalkerElement(label, visual=POINTS, *, x=0.0, v=0.0, sigma=1.5, tau0=1.0, speed=1.0)

Bases: Element

Walker that flies at constant speed in a random direction for a power-law flight time, then turns.

The flight is a rate interaction of the position, and the turns a jump interaction that resets the velocity after every flight time.

Parameters:

Name Type Description Default
label str

The label of the walker or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 POINTS
x Annotated[ArrayLike, StateSpec(vector=True, role=position)]

The initial position.

 0.0
v Annotated[ArrayLike, StateSpec(vector=True, role=velocity)]

The initial velocity.

 0.0
sigma Annotated[float, ParameterSpec(domain=open(0, oo))]

The tail exponent of the flight times, between 1 and 2 for superdiffusion.

 1.5
tau0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The shortest flight time.

 1.0
speed Annotated[float, ParameterSpec(domain=open(0, oo))]

The speed of every flight.

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import LevyWalkerElement
>>> LevyWalkerElement('walker', speed=2.0).speed
2.0
Source code in src/skcomplex/models/physics/soft_matter/_levy_walker.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = POINTS,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    v: Annotated[npt.ArrayLike, StateSpec(vector=True, role='velocity')] = 0.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5,
    tau0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    speed: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.sigma = sigma
    self.tau0 = tau0
    self.speed = speed

ModulatedParticleElement(label, visual=None, *, x=50.0, s0=1.0, ds=1.8, L=100.0)

Bases: Element

Particle at the position x in a box [0, L] whose walls reflect it.

The noise amplitude sigma = s0 + ds (x / L - 1/2) grows across the box. The walls reflect the particle at every frame, so frames must be short against the time to cross the box.

Parameters:

Name Type Description Default
label str

The label of the particle or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
x Annotated[float, StateSpec(role=position)]

The initial position.

 50.0
s0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The noise amplitude at the middle of the box.

 1.0
ds Annotated[float, ParameterSpec()]

The change of the noise amplitude across the box.

 1.8
L Annotated[float, ParameterSpec(domain=open(0, oo))]

The length of the box.

 100.0

Examples:

>>> from skcomplex.models.physics.soft_matter import ModulatedParticleElement
>>> ModulatedParticleElement('bead').ds
1.8
Source code in src/skcomplex/models/physics/soft_matter/_modulated_particle.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(role='position')] = 50.0,
    s0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    ds: Annotated[float, ParameterSpec()] = 1.8,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.s0 = s0
    self.ds = ds
    self.L = L

OpticalTrapInteraction(label, *, target, k=1.0, amplitude=0.0, frequency=0.0, blink_period=None, duty=0.5)

Bases: GradientInteraction

Harmonic optical trap on overdamped particles, the potential sum_d k_d (x_d - c_d(t))^2 / 2 of every axis.

The center oscillates as c(t) = amplitude sin(2 pi frequency t), and a blinking trap is on for the first duty fraction of every blink_period and off for the rest. The force divides by the friction gamma of the target.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
target Annotated[str, RoleSpec()]

The path of the particle or the group.

 required
k Annotated[ArrayLike, ParameterSpec(vector=True)]

The stiffness along every axis.

 1.0
amplitude Annotated[ArrayLike, ParameterSpec(vector=True)]

The amplitude of the oscillation of the center along every axis.

 0.0
frequency Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The frequency of the oscillation of the center.

 0.0
blink_period float | None

The period of a blinking trap, or None for a trap that is always on.

 None
duty float

The fraction of every blinking period during which the trap is on.

 0.5

Examples:

>>> from skcomplex.models.physics.soft_matter import OpticalTrapInteraction
>>> OpticalTrapInteraction('tweezers', target='bead', k=[1.0, 0.25]).k
[1.0, 0.25]
Source code in src/skcomplex/models/physics/soft_matter/_optical_trap.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    k: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 1.0,
    amplitude: Annotated[npt.ArrayLike, ParameterSpec(vector=True)] = 0.0,
    frequency: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    blink_period: float | None = None,
    duty: float = 0.5,
) -> None:
    super().__init__(label, target=target, states='x', divisor='gamma')
    self.k = k
    self.amplitude = amplitude
    self.frequency = frequency
    self.blink_period = blink_period
    self.duty = duty

ScaledWalkerElement(label, visual=None, *, x=0.0, K=1.0, alpha=0.5, t0=0.01)

Bases: Element

Walker at the position x whose diffusion coefficient K alpha (t + t0)^(alpha - 1) changes in time.

Parameters:

Name Type Description Default
label str

The label of the walker or the group.

 required
visual VisualSpec | None

How the walker is drawn, or None for none.

 None
x Annotated[float, StateSpec(role=position)]

The initial position.

 0.0
K Annotated[float, ParameterSpec(domain=open(0, oo))]

The generalized diffusion coefficient.

 1.0
alpha Annotated[float, ParameterSpec(domain=open(0, oo))]

The anomalous diffusion exponent.

 0.5
t0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The time offset of the diffusion coefficient.

 0.01

Examples:

>>> from skcomplex.models.physics.soft_matter import ScaledWalkerElement
>>> ScaledWalkerElement('walker', alpha=1.5).alpha
1.5
Source code in src/skcomplex/models/physics/soft_matter/_scaled_walker.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(role='position')] = 0.0,
    K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    t0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.K = K
    self.alpha = alpha
    self.t0 = t0

SpeckleFieldInteraction(label, *, target, strength=1.0, centers=(), intensities=(), lengths=())

Bases: GradientInteraction

Random optical speckle field of Gaussian hotspots, which drags particles up its gradient at u grad(I).

The potential is -u I(x), with the intensity I(x) = sum_i I_i exp(-|x - c_i|^2 / d_i^2) of hotspots whose centers, intensities relative to strength, and decay lengths are fixed when the field is created.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
target Annotated[str, RoleSpec()]

The path of the particle or the group.

 required
strength Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The product u I0 of the polarizability constant and the reference intensity.

 1.0
centers tuple[tuple[float, float], ...]

The centers of the hotspots.

 ()
intensities tuple[float, ...]

The intensities of the hotspots in units of the reference intensity.

 ()
lengths tuple[float, ...]

The decay lengths of the hotspots.

 ()

Examples:

>>> from skcomplex.models.physics.soft_matter import SpeckleFieldInteraction
>>> SpeckleFieldInteraction(
...     'speckle', target='bead', centers=((0.0, 0.0),), intensities=(1.0,), lengths=(1.0,)
... ).strength
1.0
Source code in src/skcomplex/models/physics/soft_matter/_speckle_field.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    strength: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    centers: tuple[tuple[float, float], ...] = (),
    intensities: tuple[float, ...] = (),
    lengths: tuple[float, ...] = (),
) -> None:
    super().__init__(label, target=target, states='x')
    self.strength = strength
    self.centers = centers
    self.intensities = intensities
    self.lengths = lengths

ThermalNoiseInteraction(label, *, target, state='x')

Bases: RateInteraction

Thermal noise of overdamped particles, the random force sqrt(2 gamma temperature) dW divided by gamma.

The friction gamma and the thermal energy temperature are parameters of the target, so every component of the position diffuses with the coefficient temperature / gamma.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
target Annotated[str, RoleSpec()]

The path of the particle or the group.

 required
state str

The name of the state that holds the positions.

 'x'

Examples:

>>> from skcomplex.models.physics.soft_matter import ThermalNoiseInteraction
>>> ThermalNoiseInteraction('bath', target='beads').divisor
'gamma'
Source code in src/skcomplex/models/physics/soft_matter/_thermal_noise.py 
37
38
def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], state: str = 'x') -> None:
    super().__init__(label, target=target, state=state, divisor='gamma')

TiltedWashboardInteraction(label, *, target, U0=1.0, L=1.0, f=0.0)

Bases: GradientInteraction

Tilted periodic potential U(x) - f x, with U(x) = U0 cos(2 pi x / L), on overdamped particles on a line.

The force divides by the friction gamma of the target.

Parameters:

Name Type Description Default
label str

The label of the interaction.

 required
target Annotated[str, RoleSpec()]

The path of the particle or the group.

 required
U0 Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The amplitude of the periodic potential, half its barrier height.

 1.0
L Annotated[float, ParameterSpec(domain=open(0, oo))]

The period of the potential.

 1.0
f Annotated[float, ParameterSpec()]

The constant non-conservative force that tilts the potential.

 0.0

Examples:

>>> from skcomplex.models.physics.soft_matter import TiltedWashboardInteraction
>>> TiltedWashboardInteraction('washboard', target='beads', f=0.5).f
0.5
Source code in src/skcomplex/models/physics/soft_matter/_tilted_washboard.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    U0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    L: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    f: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, target=target, states='x', divisor='gamma')
    self.U0 = U0
    self.L = L
    self.f = f

TransientWalkerElement(label, visual=VisualSpec('points', {'position': 'x', 'color_by': 'D', 'color_range': (0.0, 1.0)}), *, x=0.0, D=0.5, sigma=0.5, gamma=1.0)

Bases: Element

Brownian walker whose diffusion coefficient D is redrawn after waiting times of mean D^(-gamma).

The coefficient takes values of density sigma D^(sigma - 1) on (0, 1), the beta distribution B(sigma, 1), so slow walkers wait long before they change.

Parameters:

Name Type Description Default
label str

The label of the walker or the group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 VisualSpec('points', {'position': 'x', 'color_by': 'D', 'color_range': (0.0, 1.0)})
x Annotated[ArrayLike, StateSpec(vector=True, role=position)]

The initial position.

 0.0
D Annotated[float, StateSpec(domain=open(0, oo))]

The initial diffusion coefficient.

 0.5
sigma Annotated[float, ParameterSpec(domain=open(0, oo))]

The exponent of the density of diffusion coefficients near zero.

 0.5
gamma Annotated[float, ParameterSpec(domain=open(0, oo))]

The exponent of the mean waiting time D^(-gamma).

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import TransientWalkerElement
>>> TransientWalkerElement('walker', gamma=2.0).gamma
2.0
Source code in src/skcomplex/models/physics/soft_matter/_transient_walker.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('points', {'position': 'x', 'color_by': 'D', 'color_range': (0.0, 1.0)}),
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
    D: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 0.5,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.D = D
    self.sigma = sigma
    self.gamma = gamma

TrapSiteElement(label, visual=None, *, occupied=0.0, visited=0.0, clock=0.0, tau=1.0)

Bases: Element

Site of a line of traps, which holds a walker for its trapping time tau before the walker hops to the right.

The walker on a site counts the steps it has spent there with clock and hops to the next site after tau steps, and every site keeps a record of whether a walker has visited it. A walker hops only onto a site of positive trapping time, so on a lattice that reads -1 outside its edges a walker on the last site stays there.

Parameters:

Name Type Description Default
label str

The label of the site or the lattice group.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
occupied Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]

Whether a walker is on the site.

 0.0
visited Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]

Whether a walker has been on the site.

 0.0
clock Annotated[float, StateSpec(domain=NN)]

The number of steps the walker has spent on the site.

 0.0
tau Annotated[float, StateSpec(domain=open(0, oo))]

The trapping time of the site.

 1.0

Examples:

>>> from skcomplex.models.physics.soft_matter import TrapSiteElement
>>> TrapSiteElement('site', tau=3.0).tau
3.0
Source code in src/skcomplex/models/physics/soft_matter/_traps.py 
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    occupied: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    visited: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    clock: Annotated[float, StateSpec(domain=NN)] = 0.0,
    tau: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.occupied = occupied
    self.visited = visited
    self.clock = clock
    self.tau = tau

WallParticleElement(label, visual=None, *, z=1.0, D0=1.0, a=1.0, drift=0.2)

Bases: Element

Brownian particle at the distance z from a wall, which reflects it, and whose diffusion slows near the wall.

Parameters:

Name Type Description Default
label str

The label of the particle.

 required
visual VisualSpec | None

How the element is drawn, or None for default_visual_.

 None
z Annotated[float, StateSpec(domain=closed_open(0, oo))]

The distance from the wall.

 1.0
D0 Annotated[float, ParameterSpec(domain=open(0, oo))]

The diffusion coefficient far from the wall.

 1.0
a Annotated[float, ParameterSpec(domain=open(0, oo))]

The distance over which the wall slows diffusion.

 1.0
drift Annotated[float, ParameterSpec()]

The settling velocity far from the wall.

 0.2

Examples:

>>> from skcomplex.models.physics.soft_matter import WallParticleElement
>>> WallParticleElement('bead').a
1.0
Source code in src/skcomplex/models/physics/soft_matter/_wall_particle.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    z: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    D0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    drift: Annotated[float, ParameterSpec()] = 0.2,
) -> None:
    super().__init__(label, visual=visual)
    self.z = z
    self.D0 = D0
    self.a = a
    self.drift = drift

annealed_transient_time_motion(n_walkers=500, *, sigma=0.5, gamma=1.0, n_dim=2, random_state=0)

Return annealed transient time motion, Brownian walkers whose diffusion coefficient changes at random times.

Every walker diffuses with a coefficient D of density sigma D^(sigma - 1) on (0, 1) and keeps it for an exponential time of mean D^(-gamma), so slow walkers stay slow for long. For sigma < gamma the ensemble mean squared displacement grows as t^(sigma / gamma).

Parameters:

Name Type Description Default
n_walkers int

The number of independent walkers.

 500
sigma float

The exponent of the density of diffusion coefficients near zero.

 0.5
gamma float

The exponent of the mean waiting time D^(-gamma).

 1.0
n_dim int

The number of dimensions.

 2
random_state int | Generator | None

The seed or generator of the first diffusion coefficients.

 0

Returns:

Name Type Description
system System

The model, with the positions walkers.x, the coefficients walkers.D, and the change count changes.renewals.

Examples:

>>> from skcomplex.models.physics.soft_matter import annealed_transient_time_motion
>>> annealed_transient_time_motion(10, gamma=2.0).elements[0].values['gamma']
2.0
Source code in src/skcomplex/models/physics/soft_matter/_anomalous.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'renewal process', 'heterogeneity', 'sde'),
    references=('Argun et al. 2021', 'Massignan et al. 2014'),
    defaults={'until': 500.0, 'dt': 2.5, 'solver': EulerMaruyamaSolver(dt=0.5)},
)
def annealed_transient_time_motion(
    n_walkers: int = 500,
    *,
    sigma: float = 0.5,
    gamma: float = 1.0,
    n_dim: int = 2,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return annealed transient time motion, Brownian walkers whose diffusion coefficient changes at random times.

    Every walker diffuses with a coefficient `D` of density `sigma D^(sigma - 1)` on `(0, 1)` and keeps it for an
    exponential time of mean `D^(-gamma)`, so slow walkers stay slow for long. For `sigma < gamma` the ensemble mean
    squared displacement grows as `t^(sigma / gamma)`.

    Args:
        n_walkers:
            The number of independent walkers.

        sigma:
            The exponent of the density of diffusion coefficients near zero.

        gamma:
            The exponent of the mean waiting time `D^(-gamma)`.

        n_dim:
            The number of dimensions.

        random_state:
            The seed or generator of the first diffusion coefficients.

    Returns:
        system:
            The model, with the positions `walkers.x`, the coefficients `walkers.D`, and the change count
            `changes.renewals`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import annealed_transient_time_motion
        >>> annealed_transient_time_motion(10, gamma=2.0).elements[0].values['gamma']
        2.0
    """
    coefficients = np.random.default_rng(random_state).beta(sigma, 1.0, n_walkers)
    walkers = GroupElement(
        'walkers',
        element=TransientWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'D': coefficients, 'sigma': sigma, 'gamma': gamma},
    )
    noise = RateInteraction('noise', target='walkers', noise=lambda s: sqrt(2 * s.target.D))
    changes = JumpInteraction(
        'changes',
        target='walkers',
        waiting=lambda s: draw('wait', 'expon', scale=s.target.D ** (-s.target.gamma)),
        reset={'D': lambda s: draw('coefficient', 'beta', a=s.target.sigma, b=1)},
    )
    return System(elements=[walkers], interactions=[noise, changes], space=ContinuousSpace(n_dim=n_dim))

biased_trap_walk(n_walkers=100, *, n_sites=300, alpha=0.5, tau0=1.0, random_state=0)

Return walkers that hop always to the right through quenched traps with broadly distributed trapping times.

Every walker moves along its own line of sites and stays tau_i steps on site i before it hops to site i + 1. The trapping times are drawn once, as ceil(tau0 u^(-1 / alpha)) for uniform numbers u, so their distribution decays as tau0^alpha / tau^(1 + alpha). For alpha < 1 their mean diverges, the time to cross N sites grows as N^(1 / alpha), and the displacement grows subdiffusively as x ~ (t / tau0)^alpha, with the mean <x> = (t / tau0)^alpha sin(pi alpha) / (pi alpha) at long times. The walkers start on the first site.

Parameters:

Name Type Description Default
n_walkers int

The number of walkers, each in its own landscape.

 100
n_sites int

The number of sites of every landscape.

 300
alpha float

The exponent of the distribution of trapping times.

 0.5
tau0 float

The smallest trapping time.

 1.0
random_state int | Generator | None

The seed or generator of the trapping times.

 0

Returns:

Name Type Description
system System

The model, with the lattice traps of shape (n_walkers, n_sites) and the mean displacement position.

Examples:

>>> from skcomplex.models.physics.soft_matter import biased_trap_walk
>>> biased_trap_walk(4, n_sites=10).elements[0].check().prepare().size_
40
Source code in src/skcomplex/models/physics/soft_matter/_traps.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
@register_model(
    field='physics/soft_matter',
    tags=('lattice', 'anomalous diffusion', 'disorder', 'traps', 'subdiffusion'),
    references=('Bertin 2021', 'Bouchaud and Georges 1990'),
    defaults={'until': 2000, 'dt': 10},
    observables=_trap_observables,
)
def biased_trap_walk(
    n_walkers: int = 100,
    *,
    n_sites: int = 300,
    alpha: float = 0.5,
    tau0: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return walkers that hop always to the right through quenched traps with broadly distributed trapping times.

    Every walker moves along its own line of sites and stays `tau_i` steps on site `i` before it hops to site
    `i + 1`. The trapping times are drawn once, as `ceil(tau0 u^(-1 / alpha))` for uniform numbers `u`, so their
    distribution decays as `tau0^alpha / tau^(1 + alpha)`. For `alpha < 1` their mean diverges, the time to cross `N`
    sites grows as `N^(1 / alpha)`, and the displacement grows subdiffusively as `x ~ (t / tau0)^alpha`, with the
    mean `<x> = (t / tau0)^alpha sin(pi alpha) / (pi alpha)` at long times. The walkers start on the first site.

    Args:
        n_walkers:
            The number of walkers, each in its own landscape.

        n_sites:
            The number of sites of every landscape.

        alpha:
            The exponent of the distribution of trapping times.

        tau0:
            The smallest trapping time.

        random_state:
            The seed or generator of the trapping times.

    Returns:
        system:
            The model, with the lattice `traps` of shape `(n_walkers, n_sites)` and the mean displacement `position`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import biased_trap_walk
        >>> biased_trap_walk(4, n_sites=10).elements[0].check().prepare().size_
        40
    """
    shape = (n_walkers, n_sites)
    tau = np.ceil(tau0 * (1 - np.random.default_rng(random_state).random(shape)) ** (-1 / alpha))
    start = np.zeros(shape)
    start[:, 0] = 1.0
    traps = LatticeElement(
        'traps',
        element=TrapSiteElement,
        space=LatticeSpace(shape, boundary='fixed', fill=-1.0),
        values={'occupied': start, 'visited': start, 'tau': tau},
        visual=VISITED,
    )
    hops = LocalRuleInteraction(
        'hops',
        target='traps',
        state=('occupied', 'visited', 'clock'),
        rule=_hop,
        neighborhood=ALONG_ROW,
        reads=('occupied', 'visited', 'clock', 'tau'),
    )
    return System(elements=[traps], interactions=[hops])

brownian_motion(n_particles=500, *, inertial=True, m=1.0, gamma=1.0, kT=1.0, n_dim=1, random_state=0)

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.

Parameters:

Name Type Description Default
n_particles int

The number of particles.

 500
inertial bool

Whether the particles have inertia.

 True
m float

The mass, used with inertia.

 1.0
gamma float

The friction coefficient.

 1.0
kT float

The thermal energy.

 1.0
n_dim int

The number of dimensions.

 1
random_state int | Generator | None

The seed or generator of the thermal initial velocities.

 0

Returns:

Name Type Description
system 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'
Source code in src/skcomplex/models/physics/soft_matter/_diffusion.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
@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)

brownian_near_wall(z=1.0, *, D0=1.0, a=1.0, drift=0.2, convention='anti-ito')

Return a Brownian particle whose diffusion slows near a wall, D(z) = D0 z² / (z² + a²), as it settles on it.

The diffusion is a smooth form of the hydrodynamic hindrance that vanishes at the wall, and the drift is the settling velocity scaled by the same mobility. Diffusion that depends on the position makes the noise multiplicative, and the anti-Itô convention gives the equilibrium of physical particles. The laws are symmetric about the wall at z = 0, which reflects the particle at every frame.

Parameters:

Name Type Description Default
z float

The initial distance from the wall.

 1.0
D0 float

The diffusion coefficient far from the wall.

 1.0
a float

The distance over which the wall slows diffusion, such as the particle radius.

 1.0
drift float

The settling velocity far from the wall.

 0.2
convention str

The convention of the noise, 'ito', 'stratonovich', or 'anti-ito'.

 'anti-ito'

Returns:

Name Type Description
system System

The model, with the distance particle.z.

Examples:

>>> from skcomplex.models.physics.soft_matter import brownian_near_wall
>>> brownian_near_wall(convention='ito').interactions[0].convention
'ito'
Source code in src/skcomplex/models/physics/soft_matter/_brownian.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
@register_model(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'multiplicative noise', 'hydrodynamics'),
    references=('Brenner 1961', 'Volpe and Wehr 2016'),
    defaults={'until': 10.0, 'dt': 0.1},
)
def brownian_near_wall(
    z: float = 1.0,
    *,
    D0: float = 1.0,
    a: float = 1.0,
    drift: float = 0.2,
    convention: str = 'anti-ito',
) -> System:
    """Return a Brownian particle whose diffusion slows near a wall, `D(z) = D0 z² / (z² + a²)`, as it settles on it.

    The diffusion is a smooth form of the hydrodynamic hindrance that vanishes at the wall, and the drift is the
    settling velocity scaled by the same mobility. Diffusion that depends on the position makes the noise
    multiplicative, and the anti-Itô convention gives the equilibrium of physical particles. The laws are symmetric
    about the wall at `z = 0`, which reflects the particle at every frame.

    Args:
        z:
            The initial distance from the wall.

        D0:
            The diffusion coefficient far from the wall.

        a:
            The distance over which the wall slows diffusion, such as the particle radius.

        drift:
            The settling velocity far from the wall.

        convention:
            The convention of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.

    Returns:
        system:
            The model, with the distance `particle.z`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import brownian_near_wall
        >>> brownian_near_wall(convention='ito').interactions[0].convention
        'ito'
    """
    mobility = RateInteraction(
        'mobility',
        target='particle',
        state='z',
        rate=lambda s: -s.target.drift * s.target.z * abs(s.target.z) / (s.target.z**2 + s.target.a**2),
        noise=lambda s: sqrt(2 * s.target.D0) * s.target.z / sqrt(s.target.z**2 + s.target.a**2),
        convention=convention,
    )
    particle = WallParticleElement('particle', z=z, D0=D0, a=a, drift=drift)
    return System(elements=[particle], interactions=[mobility])

colored_noise_langevin(n_particles=500, *, gamma=1.0, D=1.0, tau=1.0, random_state=0)

Return particles whose velocity dv/dt = -gamma v + xi is driven by Ornstein-Uhlenbeck noise xi.

The noise follows tau dxi = -xi dt + sqrt(2 D) dW, so its correlation is <xi(t) xi(t')> = D / tau exp(-|t - t'| / tau), and it tends to the white noise of intensity 2 D as tau goes to 0. The stationary velocity variance is <v^2> = D / (gamma (1 + gamma tau)), which is D / gamma in the white noise limit. The velocities start at 0 and the noise starts in its stationary state. Solver steps must be short against tau.

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 500
gamma float

The friction coefficient.

 1.0
D float

The noise intensity.

 1.0
tau float

The correlation time of the noise.

 1.0
random_state int | Generator | None

The seed or generator of the initial noise values.

 0

Returns:

Name Type Description
system System

The model, with the velocities particles.v and the noise particles.xi.

Examples:

>>> from skcomplex.models.physics.soft_matter import colored_noise_langevin
>>> colored_noise_langevin(10, tau=0.5).elements[0].values['tau']
0.5
Source code in src/skcomplex/models/physics/soft_matter/_colored_noise.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
@register_model(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'colored noise', 'langevin equation'),
    references=('Bertin 2021', 'Hänggi and Jung 1995'),
    defaults={'until': 50.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def colored_noise_langevin(
    n_particles: int = 500,
    *,
    gamma: float = 1.0,
    D: float = 1.0,
    tau: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return particles whose velocity `dv/dt = -gamma v + xi` is driven by Ornstein-Uhlenbeck noise `xi`.

    The noise follows `tau dxi = -xi dt + sqrt(2 D) dW`, so its correlation is `<xi(t) xi(t')> = D / tau
    exp(-|t - t'| / tau)`, and it tends to the white noise of intensity `2 D` as `tau` goes to 0. The stationary
    velocity variance is `<v^2> = D / (gamma (1 + gamma tau))`, which is `D / gamma` in the white noise limit. The
    velocities start at 0 and the noise starts in its stationary state. Solver steps must be short against `tau`.

    Args:
        n_particles:
            The number of independent particles.

        gamma:
            The friction coefficient.

        D:
            The noise intensity.

        tau:
            The correlation time of the noise.

        random_state:
            The seed or generator of the initial noise values.

    Returns:
        system:
            The model, with the velocities `particles.v` and the noise `particles.xi`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import colored_noise_langevin
        >>> colored_noise_langevin(10, tau=0.5).elements[0].values['tau']
        0.5
    """
    noise = np.random.default_rng(random_state).normal(0.0, math.sqrt(D / tau), n_particles)
    particles = GroupElement(
        'particles',
        element=ColoredNoiseParticleElement,
        size=n_particles,
        values={'xi': noise, 'gamma': gamma, 'D': D, 'tau': tau},
    )
    drag = RateInteraction(
        'drag',
        target='particles',
        state='v',
        rate=lambda s: -s.target.gamma * s.target.v + s.target.xi,
    )
    noise = RelaxationInteraction(
        'noise',
        target='particles',
        state='xi',
        time=lambda s: s.target.tau,
        noise=lambda s: sqrt(2 * s.target.D) / s.target.tau,
    )
    return System(elements=[particles], interactions=[drag, noise])

continuous_time_random_walk(n_walkers=500, *, alpha=0.5, tau0=1.0, sigma=1.0, n_dim=2)

Return the continuous-time random walk, Gaussian jumps after waiting times with the tail tau^(-1 - alpha).

For alpha < 1 the waiting times have no mean, so the number of jumps grows as t^alpha and the ensemble mean squared displacement grows as t^alpha, slower than Brownian motion. The walkers renew on their own clocks.

Parameters:

Name Type Description Default
n_walkers int

The number of independent walkers.

 500
alpha float

The tail exponent of the waiting times.

 0.5
tau0 float

The shortest waiting time.

 1.0
sigma float

The standard deviation of the jumps along every axis.

 1.0
n_dim int

The number of dimensions.

 2

Returns:

Name Type Description
system System

The model, with the positions walkers.x and the jump count jumps.renewals.

Examples:

>>> from skcomplex.models.physics.soft_matter import continuous_time_random_walk
>>> continuous_time_random_walk(10, alpha=0.7).elements[0].values['alpha']
0.7
Source code in src/skcomplex/models/physics/soft_matter/_anomalous.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'renewal process', 'subdiffusion', 'random walk'),
    references=('Argun et al. 2021', 'Metzler and Klafter 2000', 'Montroll and Weiss 1965', 'Bertin 2021'),
    defaults={'until': 1000.0, 'dt': 5.0},
)
def continuous_time_random_walk(
    n_walkers: int = 500,
    *,
    alpha: float = 0.5,
    tau0: float = 1.0,
    sigma: float = 1.0,
    n_dim: int = 2,
) -> System:
    """Return the continuous-time random walk, Gaussian jumps after waiting times with the tail `tau^(-1 - alpha)`.

    For `alpha < 1` the waiting times have no mean, so the number of jumps grows as `t^alpha` and the ensemble mean
    squared displacement grows as `t^alpha`, slower than Brownian motion. The walkers renew on their own clocks.

    Args:
        n_walkers:
            The number of independent walkers.

        alpha:
            The tail exponent of the waiting times.

        tau0:
            The shortest waiting time.

        sigma:
            The standard deviation of the jumps along every axis.

        n_dim:
            The number of dimensions.

    Returns:
        system:
            The model, with the positions `walkers.x` and the jump count `jumps.renewals`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import continuous_time_random_walk
        >>> continuous_time_random_walk(10, alpha=0.7).elements[0].values['alpha']
        0.7
    """
    walkers = GroupElement(
        'walkers',
        element=ContinuousTimeWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'alpha': alpha, 'tau0': tau0, 'sigma': sigma},
    )
    jumps = JumpInteraction(
        'jumps',
        target='walkers',
        waiting=lambda s: draw('wait', 'pareto', b=s.target.alpha, scale=s.target.tau0),
        reset={
            'x': lambda s: [item + normal(f'jump{index}', 0, s.target.sigma) for index, item in enumerate(s.target.x)],
        },
    )
    return System(elements=[walkers], interactions=[jumps], space=ContinuousSpace(n_dim=n_dim))

diffusion_in_a_speckle_field(n_particles=100, *, n_hotspots=50, size=100.0, strength=30.0, length=10.0, D=1.0, random_state=0)

Return Brownian particles in a random speckle field, trapped at hotspots for longer as the field grows stronger.

The hotspots sit at random places in a square, with intensities r^(-1/2) and decay lengths length r^(-1/2) for independent uniform numbers r. The particles start at random places in the middle half of the square.

The default strength 30 draws the particles into the hotspots over the first ten time units and lets them escape now and then, so from time 10 to 60 their mean squared displacement grows as t^0.36, a tenth of free diffusion. A strength of 100 or more traps every particle within the first frames and leaves only the rattling inside one hotspot, with a flat mean squared displacement.

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 100
n_hotspots int

The number of hotspots.

 50
size float

The side of the square of the hotspots.

 100.0
strength float

The product u I0 of the polarizability constant and the reference intensity.

 30.0
length float

The reference decay length of the hotspots.

 10.0
D float

The diffusion coefficient.

 1.0
random_state int | Generator | None

The seed or generator of the hotspots and the initial positions.

 0

Returns:

Name Type Description
system System

The model, with the positions beads.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import diffusion_in_a_speckle_field
>>> len(diffusion_in_a_speckle_field(5, n_hotspots=3).interactions[0].centers)
3
Source code in src/skcomplex/models/physics/soft_matter/_diffusion.py 
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'sde', 'optical forces', 'disorder'),
    references=('Argun et al. 2021', 'Volpe, Volpe, and Gigan 2014'),
    defaults={'until': 60.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def diffusion_in_a_speckle_field(
    n_particles: int = 100,
    *,
    n_hotspots: int = 50,
    size: float = 100.0,
    strength: float = 30.0,
    length: float = 10.0,
    D: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Brownian particles in a random speckle field, trapped at hotspots for longer as the field grows stronger.

    The hotspots sit at random places in a square, with intensities `r^(-1/2)` and decay lengths `length r^(-1/2)` for
    independent uniform numbers `r`. The particles start at random places in the middle half of the square.

    The default strength 30 draws the particles into the hotspots over the first ten time units and lets them escape
    now and then, so from time 10 to 60 their mean squared displacement grows as `t^0.36`, a tenth of free diffusion.
    A strength of 100 or more traps every particle within the first frames and leaves only the rattling inside one
    hotspot, with a flat mean squared displacement.

    Args:
        n_particles:
            The number of independent particles.

        n_hotspots:
            The number of hotspots.

        size:
            The side of the square of the hotspots.

        strength:
            The product `u I0` of the polarizability constant and the reference intensity.

        length:
            The reference decay length of the hotspots.

        D:
            The diffusion coefficient.

        random_state:
            The seed or generator of the hotspots and the initial positions.

    Returns:
        system:
            The model, with the positions `beads.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import diffusion_in_a_speckle_field
        >>> len(diffusion_in_a_speckle_field(5, n_hotspots=3).interactions[0].centers)
        3
    """
    rng = np.random.default_rng(random_state)
    centers = rng.uniform(0.0, size, (n_hotspots, 2))
    intensities = (1 - rng.random(n_hotspots)) ** -0.5
    lengths = length * (1 - rng.random(n_hotspots)) ** -0.5
    speckle = SpeckleFieldInteraction(
        'speckle',
        target='beads',
        strength=strength,
        centers=tuple((float(x), float(y)) for x, y in centers),
        intensities=tuple(float(item) for item in intensities),
        lengths=tuple(float(item) for item in lengths),
    )
    positions = rng.uniform(size / 4, 3 * size / 4, (n_particles, 2))
    beads = GroupElement(
        'beads',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': positions, 'gamma': 1.0, 'temperature': D},
    )
    bath = ThermalNoiseInteraction('bath', target='beads')
    return System(elements=[beads], interactions=[speckle, bath], space=ContinuousSpace(n_dim=2))

flashing_ratchet(n_particles=500, *, barrier=1.0, asymmetry=0.1, L=1.0, rate=1.0, duration=100.0, gamma=1.0, kT=0.05, random_state=0)

Return the on-off flashing ratchet, Brownian particles in a sawtooth potential switched on and off at random.

The potential switches at the rate rate, so on and off phases last 1 / rate on average. While it is off the particles diffuse freely, and while it is on they slide to the minima of the sawtooth. A potential that rises over a short fraction asymmetry of its period sends more particles over its near maximum than over its far one, so the particles drift with a mean velocity of the sign of 0.5 - asymmetry, and a symmetric potential gives no current. The switching times are drawn once, as a Poisson process over duration shared by every particle. The particles start at the minimum at the origin.

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 500
barrier float

The height of the sawtooth.

 1.0
asymmetry float

The fraction of the period over which the potential rises, 0.5 for a symmetric potential.

 0.1
L float

The period of the sawtooth.

 1.0
rate float

The rate at which the potential switches.

 1.0
duration float

The time span of the drawn switching times.

 100.0
gamma float

The friction coefficient.

 1.0
kT float

The thermal energy.

 0.05
random_state int | Generator | None

The seed or generator of the switching times.

 0

Returns:

Name Type Description
system System

The model, with the positions particles.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import flashing_ratchet
>>> len(flashing_ratchet(10, rate=2.0, duration=5.0).interactions[0].switches) > 0
True
Source code in src/skcomplex/models/physics/soft_matter/_ratchets.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
@register_model(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'ratchet', 'nonequilibrium steady state', 'current'),
    references=('Bertin 2021', 'Ajdari and Prost 1992', 'Reimann 2002'),
    defaults={'until': 100.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.002)},
)
def flashing_ratchet(
    n_particles: int = 500,
    *,
    barrier: float = 1.0,
    asymmetry: float = 0.1,
    L: float = 1.0,
    rate: float = 1.0,
    duration: float = 100.0,
    gamma: float = 1.0,
    kT: float = 0.05,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the on-off flashing ratchet, Brownian particles in a sawtooth potential switched on and off at random.

    The potential switches at the rate `rate`, so on and off phases last `1 / rate` on average. While it is off the
    particles diffuse freely, and while it is on they slide to the minima of the sawtooth. A potential that rises over
    a short fraction `asymmetry` of its period sends more particles over its near maximum than over its far one, so
    the particles drift with a mean velocity of the sign of `0.5 - asymmetry`, and a symmetric potential gives no
    current. The switching times are drawn once, as a Poisson process over `duration` shared by every particle. The
    particles start at the minimum at the origin.

    Args:
        n_particles:
            The number of independent particles.

        barrier:
            The height of the sawtooth.

        asymmetry:
            The fraction of the period over which the potential rises, 0.5 for a symmetric potential.

        L:
            The period of the sawtooth.

        rate:
            The rate at which the potential switches.

        duration:
            The time span of the drawn switching times.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

        random_state:
            The seed or generator of the switching times.

    Returns:
        system:
            The model, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import flashing_ratchet
        >>> len(flashing_ratchet(10, rate=2.0, duration=5.0).interactions[0].switches) > 0
        True
    """
    rng = np.random.default_rng(random_state)
    times = np.cumsum(rng.exponential(1 / rate, int(2 * rate * duration) + 10))
    particles = GroupElement(
        'particles',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, 1)), 'gamma': gamma, 'temperature': kT},
    )
    sawtooth = FlashingSawtoothInteraction(
        'sawtooth',
        target='particles',
        barrier=barrier,
        L=L,
        asymmetry=asymmetry,
        switches=tuple(float(time) for time in times[times < duration]),
    )
    bath = ThermalNoiseInteraction('bath', target='particles')
    return System(elements=[particles], interactions=[sawtooth, bath], space=ContinuousSpace(n_dim=1))

fractional_brownian_motion(n_walkers=500, *, H=0.3, D=1.0, n_dim=2)

Return fractional Brownian motion, dx = sqrt(2 D) dB_H, with <x^2> = 2 D t^(2 H) along every axis.

The increments are correlated over the whole run: they anticorrelate for H < 1/2, which makes the walkers subdiffusive, and correlate for H > 1/2, which makes them superdiffusive.

Parameters:

Name Type Description Default
n_walkers int

The number of independent walkers.

 500
H float

The Hurst index.

 0.3
D float

The diffusion coefficient.

 1.0
n_dim int

The number of dimensions.

 2

Returns:

Name Type Description
system System

The model, with the positions walkers.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import fractional_brownian_motion
>>> fractional_brownian_motion(10, H=0.7).elements[0].values['H']
0.7
Source code in src/skcomplex/models/physics/soft_matter/_anomalous.py 
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'fractional noise', 'long memory', 'sde'),
    references=('Argun et al. 2021', 'Mandelbrot and Van Ness 1968'),
    defaults={'until': 10.0, 'dt': 0.05, 'solver': FractionalEulerMaruyamaSolver(dt=0.01)},
)
def fractional_brownian_motion(
    n_walkers: int = 500,
    *,
    H: float = 0.3,
    D: float = 1.0,
    n_dim: int = 2,
) -> System:
    """Return fractional Brownian motion, `dx = sqrt(2 D) dB_H`, with `<x^2> = 2 D t^(2 H)` along every axis.

    The increments are correlated over the whole run: they anticorrelate for `H < 1/2`, which makes the walkers
    subdiffusive, and correlate for `H > 1/2`, which makes them superdiffusive.

    Args:
        n_walkers:
            The number of independent walkers.

        H:
            The Hurst index.

        D:
            The diffusion coefficient.

        n_dim:
            The number of dimensions.

    Returns:
        system:
            The model, with the positions `walkers.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import fractional_brownian_motion
        >>> fractional_brownian_motion(10, H=0.7).elements[0].values['H']
        0.7
    """
    walkers = GroupElement(
        'walkers',
        element=FractionalWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'D': D, 'H': H},
    )
    noise = RateInteraction(
        'noise',
        target='walkers',
        noise=lambda s: sqrt(2 * s.target.D),
        hurst=lambda s: s.target.H,
    )
    return System(elements=[walkers], interactions=[noise], space=ContinuousSpace(n_dim=n_dim))

kramers_escape(n_particles=200, *, D=0.1, random_state=0)

Return particles that escape over the barrier of U = x^4 / 4 - x^2 / 2 by thermal noise of intensity D.

The particles start in the left well, and they hop between the wells at the Kramers rate sqrt(U''(-1) |U''(0)|) / (2 pi) exp(-1 / (4 D)), which is sqrt(2) / (2 pi) exp(-1 / (4 D)).

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 200
D float

The noise intensity.

 0.1
random_state int | Generator | None

The seed or generator of the small initial spread in the left well.

 0

Returns:

Name Type Description
system System

The model, with the positions particles.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import kramers_escape
>>> kramers_escape(10, D=0.2).elements[0].values['D']
0.2
Source code in src/skcomplex/models/physics/soft_matter/_diffusion.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
@register_model(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'escape', 'double well', 'rates'),
    references=('Kramers 1940', 'Hänggi, Talkner, and Borkovec 1990', 'Gros 2015'),
    defaults={'until': 100.0, 'dt': 1.0, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def kramers_escape(
    n_particles: int = 200,
    *,
    D: float = 0.1,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return particles that escape over the barrier of `U = x^4 / 4 - x^2 / 2` by thermal noise of intensity `D`.

    The particles start in the left well, and they hop between the wells at the Kramers rate
    `sqrt(U''(-1) |U''(0)|) / (2 pi) exp(-1 / (4 D))`, which is `sqrt(2) / (2 pi) exp(-1 / (4 D))`.

    Args:
        n_particles:
            The number of independent particles.

        D:
            The noise intensity.

        random_state:
            The seed or generator of the small initial spread in the left well.

    Returns:
        system:
            The model, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import kramers_escape
        >>> kramers_escape(10, D=0.2).elements[0].values['D']
        0.2
    """
    positions = -1.0 + np.random.default_rng(random_state).normal(0.0, 0.1, n_particles)
    group = GroupElement(
        'particles',
        element=DoubleWellElement,
        size=n_particles,
        values={'x': positions, 'amplitude': 0.0, 'D': D},
    )
    return System(elements=[group], interactions=_get_double_well_interactions('particles'))

levy_flight(n_walkers=500, *, mu=1.2, l0=1.0, n_dim=2)

Return the Lévy flight, a walk whose jumps have lengths with the tail l^(-1 - mu) in random directions.

For mu < 2 the jump lengths have no variance, so the walk is dominated by its longest jumps and the typical distance travelled after n steps grows as n^(1 / mu), faster than the n^(1/2) of a random walk.

Parameters:

Name Type Description Default
n_walkers int

The number of independent walkers.

 500
mu float

The tail exponent of the jump lengths.

 1.2
l0 float

The shortest jump length.

 1.0
n_dim int

The number of dimensions.

 2

Returns:

Name Type Description
system System

The model, with the positions walkers.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import levy_flight
>>> levy_flight(10, mu=1.5).elements[0].values['mu']
1.5
Source code in src/skcomplex/models/physics/soft_matter/_anomalous.py 
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'superdiffusion', 'random walk', 'stochastic map'),
    references=('Bertin 2021', 'Shlesinger, Zaslavsky, and Frisch 1995'),
    defaults={'until': 200.0, 'dt': 1.0},
)
def levy_flight(
    n_walkers: int = 500,
    *,
    mu: float = 1.2,
    l0: float = 1.0,
    n_dim: int = 2,
) -> System:
    """Return the Lévy flight, a walk whose jumps have lengths with the tail `l^(-1 - mu)` in random directions.

    For `mu < 2` the jump lengths have no variance, so the walk is dominated by its longest jumps and the typical
    distance travelled after `n` steps grows as `n^(1 / mu)`, faster than the `n^(1/2)` of a random walk.

    Args:
        n_walkers:
            The number of independent walkers.

        mu:
            The tail exponent of the jump lengths.

        l0:
            The shortest jump length.

        n_dim:
            The number of dimensions.

    Returns:
        system:
            The model, with the positions `walkers.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import levy_flight
        >>> levy_flight(10, mu=1.5).elements[0].values['mu']
        1.5
    """
    walkers = GroupElement(
        'walkers',
        element=LevyFlierElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'mu': mu, 'l0': l0},
    )
    flights = MapInteraction('flights', target='walkers', state='x', function=_get_flight)
    return System(elements=[walkers], interactions=[flights], space=ContinuousSpace(n_dim=n_dim))

levy_walk(n_walkers=500, *, sigma=1.5, tau0=1.0, speed=1.0, n_dim=2, random_state=0)

Return the Lévy walk, flights at constant speed in random directions for times with the tail tau^(-1 - sigma).

For 1 < sigma < 2 the flight times have a mean but no variance, so long flights dominate and the ensemble mean squared displacement grows as t^(3 - sigma), faster than Brownian motion and slower than ballistic motion.

Parameters:

Name Type Description Default
n_walkers int

The number of independent walkers.

 500
sigma float

The tail exponent of the flight times.

 1.5
tau0 float

The shortest flight time.

 1.0
speed float

The speed of every flight.

 1.0
n_dim int

The number of dimensions.

 2
random_state int | Generator | None

The seed or generator of the directions of the first flights.

 0

Returns:

Name Type Description
system System

The model, with the positions walkers.x, the velocities walkers.v, and the turn count turns.renewals.

Examples:

>>> from skcomplex.models.physics.soft_matter import levy_walk
>>> levy_walk(10, sigma=1.2).elements[0].values['sigma']
1.2
Source code in src/skcomplex/models/physics/soft_matter/_anomalous.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'renewal process', 'superdiffusion', 'random walk'),
    references=('Argun et al. 2021', 'Zaburdaev, Denisov, and Klafter 2015', 'Shlesinger, West, and Klafter 1987'),
    defaults={'until': 500.0, 'dt': 2.5, 'solver': RK4Solver(dt=0.5)},
)
def levy_walk(
    n_walkers: int = 500,
    *,
    sigma: float = 1.5,
    tau0: float = 1.0,
    speed: float = 1.0,
    n_dim: int = 2,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Lévy walk, flights at constant speed in random directions for times with the tail `tau^(-1 - sigma)`.

    For `1 < sigma < 2` the flight times have a mean but no variance, so long flights dominate and the ensemble mean
    squared displacement grows as `t^(3 - sigma)`, faster than Brownian motion and slower than ballistic motion.

    Args:
        n_walkers:
            The number of independent walkers.

        sigma:
            The tail exponent of the flight times.

        tau0:
            The shortest flight time.

        speed:
            The speed of every flight.

        n_dim:
            The number of dimensions.

        random_state:
            The seed or generator of the directions of the first flights.

    Returns:
        system:
            The model, with the positions `walkers.x`, the velocities `walkers.v`, and the turn count `turns.renewals`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import levy_walk
        >>> levy_walk(10, sigma=1.2).elements[0].values['sigma']
        1.2
    """
    directions = np.random.default_rng(random_state).standard_normal((n_walkers, n_dim))
    velocities = speed * directions / np.linalg.norm(directions, axis=1, keepdims=True)
    walkers = GroupElement(
        'walkers',
        element=LevyWalkerElement,
        size=n_walkers,
        values={'x': np.zeros((n_walkers, n_dim)), 'v': velocities, 'sigma': sigma, 'tau0': tau0, 'speed': speed},
    )
    flight = RateInteraction('flight', target='walkers', state='x', rate=lambda s: s.target.v)
    turns = JumpInteraction(
        'turns',
        target='walkers',
        waiting=lambda s: draw('flight', 'pareto', b=s.target.sigma, scale=s.target.tau0),
        reset={'v': lambda s: [s.target.speed * item for item in _direction('direction', len(list(s.target.v)))]},
    )
    return System(elements=[walkers], interactions=[flight, turns], space=ContinuousSpace(n_dim=n_dim))

multiplicative_noise(n_particles=500, *, s0=1.0, ds=1.8, L=100.0, convention='anti-ito')

Return free particles in a box with position-dependent noise, read in the Itô, Stratonovich, or anti-Itô sense.

Only the anti-Itô convention gives the flat stationary density of particles at thermal equilibrium without forces. The particles start at the middle of the box.

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 500
s0 float

The noise amplitude at the middle of the box.

 1.0
ds float

The change of the noise amplitude across the box.

 1.8
L float

The length of the box.

 100.0
convention str

The convention of the noise, 'ito', 'stratonovich', or 'anti-ito'.

 'anti-ito'

Returns:

Name Type Description
system System

The model, with the positions particles.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import multiplicative_noise
>>> multiplicative_noise(10, convention='ito').interactions[0].convention
'ito'
Source code in src/skcomplex/models/physics/soft_matter/_diffusion.py 
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
@register_model(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'multiplicative noise', 'spurious drift'),
    references=('Argun et al. 2021', 'Volpe and Wehr 2016'),
    defaults={'until': 100.0, 'dt': 0.05, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def multiplicative_noise(
    n_particles: int = 500,
    *,
    s0: float = 1.0,
    ds: float = 1.8,
    L: float = 100.0,
    convention: str = 'anti-ito',
) -> System:
    """Return free particles in a box with position-dependent noise, read in the Itô, Stratonovich, or anti-Itô sense.

    Only the anti-Itô convention gives the flat stationary density of particles at thermal equilibrium without
    forces. The particles start at the middle of the box.

    Args:
        n_particles:
            The number of independent particles.

        s0:
            The noise amplitude at the middle of the box.

        ds:
            The change of the noise amplitude across the box.

        L:
            The length of the box.

        convention:
            The convention of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.

    Returns:
        system:
            The model, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import multiplicative_noise
        >>> multiplicative_noise(10, convention='ito').interactions[0].convention
        'ito'
    """
    group = GroupElement(
        'particles',
        element=ModulatedParticleElement,
        size=n_particles,
        values={'x': L / 2, 's0': s0, 'ds': ds, 'L': L},
    )
    noise = RateInteraction(
        'noise',
        target='particles',
        noise=lambda s: s.target.s0 + s.target.ds * (s.target.x / s.target.L - 0.5),
        convention=convention,
    )
    return System(elements=[group], interactions=[noise])

optical_trap(n_particles=200, *, k=(1.0, 0.25), gamma=1.0, kT=1.0, amplitude=(0.0, 0.0), frequency=0.0, blink_period=None, duty=0.5)

Return overdamped particles in a harmonic optical trap, still, oscillating, or blinking.

In a still trap every coordinate is normal with variance kT / k and relaxes in the time gamma / k.

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 200
k Sequence[float]

The stiffness along each axis.

 (1.0, 0.25)
gamma float

The friction coefficient.

 1.0
kT float

The thermal energy.

 1.0
amplitude Sequence[float]

The amplitude of the oscillation of the center along each axis.

 (0.0, 0.0)
frequency float

The frequency of the oscillation of the center.

 0.0
blink_period float | None

The period of a blinking trap, or None for a trap that is always on.

 None
duty float

The fraction of every blinking period during which the trap is on.

 0.5

Returns:

Name Type Description
system System

The model, with the positions beads.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import optical_trap
>>> optical_trap(10, blink_period=2.0).interactions[0].blink_period
2.0
Source code in src/skcomplex/models/physics/soft_matter/_diffusion.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
@register_model(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'optical tweezers', 'forcing'),
    references=('Argun et al. 2021', 'Ashkin et al. 1986', 'Jones, Maragò, and Volpe 2015'),
    defaults={'until': 20.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def optical_trap(
    n_particles: int = 200,
    *,
    k: Sequence[float] = (1.0, 0.25),
    gamma: float = 1.0,
    kT: float = 1.0,
    amplitude: Sequence[float] = (0.0, 0.0),
    frequency: float = 0.0,
    blink_period: float | None = None,
    duty: float = 0.5,
) -> System:
    """Return overdamped particles in a harmonic optical trap, still, oscillating, or blinking.

    In a still trap every coordinate is normal with variance `kT / k` and relaxes in the time `gamma / k`.

    Args:
        n_particles:
            The number of independent particles.

        k:
            The stiffness along each axis.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

        amplitude:
            The amplitude of the oscillation of the center along each axis.

        frequency:
            The frequency of the oscillation of the center.

        blink_period:
            The period of a blinking trap, or `None` for a trap that is always on.

        duty:
            The fraction of every blinking period during which the trap is on.

    Returns:
        system:
            The model, with the positions `beads.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import optical_trap
        >>> optical_trap(10, blink_period=2.0).interactions[0].blink_period
        2.0
    """
    n_dim = len(k)
    beads = GroupElement(
        'beads',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, n_dim)), 'gamma': gamma, 'temperature': kT},
    )
    trap = OpticalTrapInteraction(
        'tweezers',
        target='beads',
        k=list(k),
        amplitude=list(amplitude),
        frequency=frequency,
        blink_period=blink_period,
        duty=duty,
    )
    bath = ThermalNoiseInteraction('bath', target='beads')
    return System(elements=[beads], interactions=[trap, bath], space=ContinuousSpace(n_dim=n_dim))

scaled_brownian_motion(n_particles=500, *, alpha=0.5, K=1.0, t0=0.01)

Return scaled Brownian motion, dx = sqrt(2 K alpha (t + t0)^(alpha - 1)) dW, with <x^2> ~ 2 K t^alpha.

The diffusion coefficient changes deterministically in time, which is the same as running Brownian motion on the rescaled clock t^alpha. The offset t0 keeps the coefficient finite at the start.

Parameters:

Name Type Description Default
n_particles int

The number of independent walkers.

 500
alpha float

The anomalous diffusion exponent, below 1 for subdiffusion and above 1 for superdiffusion.

 0.5
K float

The generalized diffusion coefficient.

 1.0
t0 float

The time offset of the diffusion coefficient.

 0.01

Returns:

Name Type Description
system System

The model, with the positions walkers.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import scaled_brownian_motion
>>> scaled_brownian_motion(10, alpha=1.5).elements[0].values['alpha']
1.5
Source code in src/skcomplex/models/physics/soft_matter/_diffusion.py 
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
@register_model(
    field='physics/soft_matter',
    tags=('anomalous diffusion', 'sde', 'non-ergodic'),
    references=('Argun et al. 2021', 'Lim and Muniandy 2002', 'Jeon et al. 2014'),
    defaults={'until': 10.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def scaled_brownian_motion(
    n_particles: int = 500,
    *,
    alpha: float = 0.5,
    K: float = 1.0,
    t0: float = 0.01,
) -> System:
    """Return scaled Brownian motion, `dx = sqrt(2 K alpha (t + t0)^(alpha - 1)) dW`, with `<x^2> ~ 2 K t^alpha`.

    The diffusion coefficient changes deterministically in time, which is the same as running Brownian motion on the
    rescaled clock `t^alpha`. The offset `t0` keeps the coefficient finite at the start.

    Args:
        n_particles:
            The number of independent walkers.

        alpha:
            The anomalous diffusion exponent, below 1 for subdiffusion and above 1 for superdiffusion.

        K:
            The generalized diffusion coefficient.

        t0:
            The time offset of the diffusion coefficient.

    Returns:
        system:
            The model, with the positions `walkers.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import scaled_brownian_motion
        >>> scaled_brownian_motion(10, alpha=1.5).elements[0].values['alpha']
        1.5
    """
    return System(
        elements=[
            GroupElement(
                'walkers',
                element=ScaledWalkerElement,
                size=n_particles,
                values={'K': K, 'alpha': alpha, 't0': t0},
            ),
        ],
        interactions=[
            RateInteraction(
                'noise',
                target='walkers',
                noise=lambda s: sqrt(2 * s.target.K * s.target.alpha * (s.t + s.target.t0) ** (s.target.alpha - 1)),
            ),
        ],
    )

shaken_spring_block(n_blocks=50, *, k=1.0, anchor=0.0, mu_s=1.0, mu_d=0.8, amplitude=2.0, shake_time=2.0, relax_time=8.0, v_stick=0.02)

Return a chain of blocks joined by springs on a frictional substrate, shaken and relaxed in turn.

Every block of unit weight feels the springs to its neighbours, k (x_{i+1} + x_{i-1} - 2 x_i), an optional spring of stiffness anchor to its rest position, and Coulomb friction. A sliding block feels the dynamic friction -mu_d sign(v), and a block at rest stays at rest while the spring force on it is below the static friction mu_s. In every cycle a random force of intensity amplitude^2 shakes the blocks for shake_time, and then the chain relaxes for relax_time to a blocked configuration, where the elongations xi_i = x_{i+1} - x_i satisfy k |xi_{i+1} - xi_i| < mu_s. Repeated cycles sample blocked configurations, whose energy per spring grows with the shaking, and whose elongations are correlated over a length that grows with that energy. Without springs between blocks and with an anchor, the blocks are independent, and their blocked positions satisfy |x_i| < mu_s / anchor. A block counts as at rest below the speed v_stick, which must exceed the change of speed (mu_s + mu_d) h of a solver step h. The ends of the chain are free, and the blocks start at rest.

Parameters:

Name Type Description Default
n_blocks int

The number of blocks.

 50
k float

The stiffness of the springs between neighbouring blocks.

 1.0
anchor float

The stiffness of the springs that hold every block to its rest position.

 0.0
mu_s float

The static friction coefficient.

 1.0
mu_d float

The dynamic friction coefficient.

 0.8
amplitude float

The amplitude of the random shaking force.

 2.0
shake_time float

The duration of the shaking in every cycle.

 2.0
relax_time float

The duration of the relaxation in every cycle.

 8.0
v_stick float

The speed below which a block counts as at rest.

 0.02

Returns:

Name Type Description
system System

The model, with the positions chain.x and velocities chain.v of the blocks.

Examples:

>>> from skcomplex.models.physics.soft_matter import shaken_spring_block
>>> shaken_spring_block(5).elements[0].check().prepare().size_
5
Source code in src/skcomplex/models/physics/soft_matter/_friction.py 
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
@register_model(
    field='physics/soft_matter',
    tags=('field', 'sde', 'friction', 'granular matter', 'blocked states'),
    references=('Bertin 2021', 'Edwards and Oakeshott 1989', 'Bertin, Dauchot, and Droz 2004'),
    defaults={'until': 200.0, 'dt': 10.0, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def shaken_spring_block(
    n_blocks: int = 50,
    *,
    k: float = 1.0,
    anchor: float = 0.0,
    mu_s: float = 1.0,
    mu_d: float = 0.8,
    amplitude: float = 2.0,
    shake_time: float = 2.0,
    relax_time: float = 8.0,
    v_stick: float = 0.02,
) -> System:
    """Return a chain of blocks joined by springs on a frictional substrate, shaken and relaxed in turn.

    Every block of unit weight feels the springs to its neighbours, `k (x_{i+1} + x_{i-1} - 2 x_i)`, an optional
    spring of stiffness `anchor` to its rest position, and Coulomb friction. A sliding block feels the dynamic friction
    `-mu_d sign(v)`, and a block at rest stays at rest while the spring force on it is below the static friction
    `mu_s`. In every cycle a random force of intensity `amplitude^2` shakes the blocks for `shake_time`, and then the
    chain relaxes for `relax_time` to a blocked configuration, where the elongations `xi_i = x_{i+1} - x_i` satisfy
    `k |xi_{i+1} - xi_i| < mu_s`. Repeated cycles sample blocked configurations, whose energy per spring grows with
    the shaking, and whose elongations are correlated over a length that grows with that energy. Without springs
    between blocks and with an anchor, the blocks are independent, and their blocked positions satisfy
    `|x_i| < mu_s / anchor`. A block counts as at rest below the speed `v_stick`, which must exceed the change of
    speed `(mu_s + mu_d) h` of a solver step `h`. The ends of the chain are free, and the blocks start at rest.

    Args:
        n_blocks:
            The number of blocks.

        k:
            The stiffness of the springs between neighbouring blocks.

        anchor:
            The stiffness of the springs that hold every block to its rest position.

        mu_s:
            The static friction coefficient.

        mu_d:
            The dynamic friction coefficient.

        amplitude:
            The amplitude of the random shaking force.

        shake_time:
            The duration of the shaking in every cycle.

        relax_time:
            The duration of the relaxation in every cycle.

        v_stick:
            The speed below which a block counts as at rest.

    Returns:
        system:
            The model, with the positions `chain.x` and velocities `chain.v` of the blocks.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import shaken_spring_block
        >>> shaken_spring_block(5).elements[0].check().prepare().size_
        5
    """
    chain = FieldElement(
        'chain',
        element=FrictionalBlockElement,
        space=LatticeSpace(n_blocks, boundary='reflecting'),
        boundary={'x': 'zero-flux', 'v': 'zero-flux'},
        values={
            'x': np.zeros(n_blocks),
            'k': k,
            'anchor': anchor,
            'mu_s': mu_s,
            'mu_d': mu_d,
            'amplitude': amplitude,
            'shake_time': shake_time,
            'relax_time': relax_time,
            'v_stick': v_stick,
        },
    )
    return System(elements=[chain])

stochastic_resonance(x=-1.0, *, amplitude=0.1, frequency=0.01, D=0.12)

Return a particle in a double well, dx = (x - x^3 + A cos(2 pi f t)) dt + sqrt(2 D) dW.

The forcing alone is too weak to move the particle between the wells, but noise of the right intensity makes the hops follow the forcing.

Parameters:

Name Type Description Default
x float

The initial position.

 -1.0
amplitude float

The amplitude of the periodic forcing.

 0.1
frequency float

The frequency of the forcing.

 0.01
D float

The noise intensity.

 0.12

Returns:

Name Type Description
system System

The model, with the position particle.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import stochastic_resonance
>>> stochastic_resonance(D=0.2).elements[0].D
0.2
Source code in src/skcomplex/models/physics/soft_matter/_brownian.py 
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
@register_model(
    field='physics/soft_matter',
    tags=('noise', 'sde', 'forcing', 'resonance'),
    references=('Benzi, Sutera, and Vulpiani 1981', 'Gammaitoni et al. 1998'),
    defaults={'until': 500.0, 'dt': 0.5},
)
def stochastic_resonance(
    x: float = -1.0,
    *,
    amplitude: float = 0.1,
    frequency: float = 0.01,
    D: float = 0.12,
) -> System:
    """Return a particle in a double well, `dx = (x - x^3 + A cos(2 pi f t)) dt + sqrt(2 D) dW`.

    The forcing alone is too weak to move the particle between the wells, but noise of the right intensity makes the
    hops follow the forcing.

    Args:
        x:
            The initial position.

        amplitude:
            The amplitude of the periodic forcing.

        frequency:
            The frequency of the forcing.

        D:
            The noise intensity.

    Returns:
        system:
            The model, with the position `particle.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import stochastic_resonance
        >>> stochastic_resonance(D=0.2).elements[0].D
        0.2
    """
    return System(
        elements=[DoubleWellElement('particle', x=x, amplitude=amplitude, frequency=frequency, D=D)],
        interactions=_get_double_well_interactions('particle'),
    )

tilted_periodic_potential(n_particles=500, *, f=1.0, U0=0.5, L=1.0, gamma=1.0, kT=0.5)

Return overdamped particles in the tilted periodic potential V(x) = U0 cos(2 pi x / L) - f x.

The force f breaks the left and right symmetry, so the particles reach a nonequilibrium steady state with a uniform probability current J. Their mean velocity is J L = L (1 - exp(-f L / kT)) / int_0^L dx int_0^L dy exp((V(x) - V(x - y)) / kT) kT / gamma, which is f / gamma without the potential and exponentially small in the barrier for a weak force. The particles start at the origin.

Parameters:

Name Type Description Default
n_particles int

The number of independent particles.

 500
f float

The constant force that tilts the potential.

 1.0
U0 float

The amplitude of the periodic potential.

 0.5
L float

The period of the potential.

 1.0
gamma float

The friction coefficient.

 1.0
kT float

The thermal energy.

 0.5

Returns:

Name Type Description
system System

The model, with the positions particles.x.

Examples:

>>> from skcomplex.models.physics.soft_matter import tilted_periodic_potential
>>> tilted_periodic_potential(10, f=2.0).interactions[0].f
2.0
Source code in src/skcomplex/models/physics/soft_matter/_ratchets.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
@register_model(
    field='physics/soft_matter',
    tags=('brownian motion', 'sde', 'nonequilibrium steady state', 'current'),
    references=('Bertin 2021', 'Risken 1996', 'Reimann 2002'),
    defaults={'until': 50.0, 'dt': 0.25, 'solver': EulerMaruyamaSolver(dt=0.005)},
)
def tilted_periodic_potential(
    n_particles: int = 500,
    *,
    f: float = 1.0,
    U0: float = 0.5,
    L: float = 1.0,
    gamma: float = 1.0,
    kT: float = 0.5,
) -> System:
    """Return overdamped particles in the tilted periodic potential `V(x) = U0 cos(2 pi x / L) - f x`.

    The force `f` breaks the left and right symmetry, so the particles reach a nonequilibrium steady state with a
    uniform probability current `J`. Their mean velocity is
    `J L = L (1 - exp(-f L / kT)) / int_0^L dx int_0^L dy exp((V(x) - V(x - y)) / kT) kT / gamma`, which is
    `f / gamma` without the potential and exponentially small in the barrier for a weak force. The particles start
    at the origin.

    Args:
        n_particles:
            The number of independent particles.

        f:
            The constant force that tilts the potential.

        U0:
            The amplitude of the periodic potential.

        L:
            The period of the potential.

        gamma:
            The friction coefficient.

        kT:
            The thermal energy.

    Returns:
        system:
            The model, with the positions `particles.x`.

    Examples:
        >>> from skcomplex.models.physics.soft_matter import tilted_periodic_potential
        >>> tilted_periodic_potential(10, f=2.0).interactions[0].f
        2.0
    """
    particles = GroupElement(
        'particles',
        element=BrownianParticleElement,
        size=n_particles,
        values={'x': np.zeros((n_particles, 1)), 'gamma': gamma, 'temperature': kT},
    )
    washboard = TiltedWashboardInteraction('washboard', target='particles', U0=U0, L=L, f=f)
    bath = ThermalNoiseInteraction('bath', target='particles')
    return System(elements=[particles], interactions=[washboard, bath], space=ContinuousSpace(n_dim=1))