Optical trap including oscillating and blinking traps
Overdamped particles in a harmonic optical trap, still, oscillating, or blinking.
- Name:
optical_trap - Field:
physics/soft_matter - Formalism: stochastic differential equations
- References: Argun et al. 2021; Ashkin et al. 1986; Jones, Maragò, and Volpe 2015
Factory
@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))
Components
class BrownianParticleElement(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`.
Args:
label:
The label of the particle or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The initial position.
gamma:
The friction coefficient.
temperature:
The thermal energy.
Examples:
>>> from skcomplex.models.physics.soft_matter import BrownianParticleElement
>>> BrownianParticleElement('bead', temperature=2.0).temperature
2.0
"""
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
class OpticalTrapInteraction(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.
Args:
label:
The label of the interaction.
target:
The path of the particle or the group.
k:
The stiffness along every axis.
amplitude:
The amplitude of the oscillation of the center along every 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.
Examples:
>>> from skcomplex.models.physics.soft_matter import OpticalTrapInteraction
>>> OpticalTrapInteraction('tweezers', target='bead', k=[1.0, 0.25]).k
[1.0, 0.25]
"""
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
def _get_potential(self, namespace: SimpleNamespace) -> Expression:
"""Return the potential of the trap, switched by the blinking.
Args:
namespace:
The namespace of the symbols the potential can use.
Returns:
potential:
The potential.
"""
s = namespace
center = s.amplitude * sin(2 * math.pi * s.frequency * s.t)
on = 1 if self.blink_period is None else pulse_train(s.t, self.blink_period, self.duty * self.blink_period)
energy = sum(
(
stiffness * (item - middle) ** 2 / 2
for stiffness, item, middle in zip(s.k, s.target.x, center, strict=True)
),
SR(0),
)
return on * energy
class ThermalNoiseInteraction(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`.
Args:
label:
The label of the interaction.
target:
The path of the particle or the group.
state:
The name of the state that holds the positions.
Examples:
>>> from skcomplex.models.physics.soft_matter import ThermalNoiseInteraction
>>> ThermalNoiseInteraction('bath', target='beads').divisor
'gamma'
"""
def __init__(self, label: str, *, target: Annotated[str, RoleSpec()], state: str = 'x') -> None:
super().__init__(label, target=target, state=state, divisor='gamma')
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the amplitude `sqrt(2 gamma temperature)` of the random force.
Args:
namespace:
The namespace of the symbols the amplitude can use.
Returns:
amplitude:
The noise amplitude, before the division by the friction.
"""
target = namespace.target
return sqrt(2 * target.gamma * target.temperature)
Equations
beads.x[i, 0]: diff(x_beads_0(t), t) == (amplitude_tweezers_0*sin(6.283185307179586*frequency_tweezers*t) - x_beads_0(t))*k_tweezers_0/gamma_beads
beads.x[i, 1]: diff(x_beads_1(t), t) == (amplitude_tweezers_1*sin(6.283185307179586*frequency_tweezers*t) - x_beads_1(t))*k_tweezers_1/gamma_beads