Lennard-Jones fluid

Lennard-Jones particles in a periodic square box, optionally held at a temperature by a thermostat.

• Name: lennard_jones_fluid
• Field: physics/statistical_mechanics
• Formalism: ordinary differential equations
• References: Lennard-Jones 1924; Verlet 1967

Factory

@register_model(
    field='physics/statistical_mechanics',
    tags=('molecular dynamics', 'pair potential', 'ode'),
    references=('Lennard-Jones 1924', 'Verlet 1967'),
    defaults={'until': 0.5, 'dt': 0.005, 'solver': VelocityVerletSolver()},
)
def lennard_jones_fluid(
    n_particles: int = 400,
    *,
    density: float = 0.5,
    temperature: float | None = None,
    cutoff: float = 2.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Lennard-Jones particles in a periodic square box, optionally held at a temperature by a thermostat.

    Args:
        n_particles:
            The number of particles.

        density:
            The number of particles per unit area.

        temperature:
            The temperature of the Langevin thermostat, or `None` for constant energy. A thermostat makes the model
            stochastic, so it runs with the Euler-Maruyama solver instead of the default velocity Verlet.

        cutoff:
            The distance beyond which pairs do not interact.

        random_state:
            The seed or generator of the initial jitter and velocities.

    Returns:
        system:
            The model, with the particle group `fluid` in two dimensions.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import lennard_jones_fluid
        >>> lennard_jones_fluid(100).elements[0].check().prepare().size_
        100
    """
    rng = np.random.default_rng(random_state)
    size = float(np.sqrt(n_particles / density))
    side = int(np.ceil(np.sqrt(n_particles)))
    spacing = size / side
    grid = np.stack(np.meshgrid(np.arange(side), np.arange(side)), axis=-1).reshape(-1, 2)[:n_particles]
    positions = np.mod((grid + 0.5) * spacing + rng.normal(0, 0.02 * spacing, (n_particles, 2)), size)
    velocities = rng.normal(0, np.sqrt(temperature or 1.0), (n_particles, 2))
    velocities -= velocities.mean(axis=0)
    interactions: list[BaseInteraction] = [
        LennardJonesInteraction('lj', first='fluid', second='fluid', cutoff=cutoff, divisor='m', period=size),
    ]
    if temperature is not None:
        interactions.append(
            LangevinThermostatInteraction('thermostat', target='fluid', temperature=temperature),
        )
    return System(
        elements=[
            GroupElement(
                'fluid',
                element=MassPointElement,
                size=n_particles,
                values={'x': positions, 'x_d1': velocities},
            ),
        ],
        interactions=interactions,
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
    )

Components

class MassPointElement(PointElement):
    """Point mass, a point whose position `x` has order 2, with its velocity `x_d1` and a mass `m`.

    The forces on the point divide by the mass through their divisor `'m'`.

    Args:
        label:
            The label of the point.

        x:
            The initial position. A number applies to every component.

        x_d1:
            The initial velocity. A number applies to every component.

        visual:
            How the point is drawn, or `None` for `default_visual_`.

        m:
            The mass.

    Attributes:
        default_visual_:
            Points at the positions.

    Examples:
        >>> from skcomplex.models.physics.common.particles import MassPointElement
        >>> MassPointElement('ball', x=[0.0, 1.0], m=2.0).m
        2.0
    """

    def __init__(
        self,
        label: str,
        x: npt.ArrayLike = 0.0,
        x_d1: npt.ArrayLike = 0.0,
        visual: VisualSpec | None = None,
        *,
        m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    ) -> None:
        super().__init__(label, x=x, x_d1=x_d1, visual=visual)
        self.m = m

class LennardJonesInteraction(PairCouplingInteraction):
    """Lennard-Jones pair potential, `4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)`, without a shift.

    Args:
        label:
            The label of the interaction.

        first:
            The path of the group or particle of the first members.

        second:
            The path of the group or particle of the second members.

        cutoff:
            The distance beyond which pairs do not interact.

        position:
            The name of the vector state that holds the positions of the members.

        write:
            The name of the state whose rate the force changes, the position of order 2 by default.

        divisor:
            The name of the parameter that holds the masses of the members, or `None` for unit masses.

        period:
            The side of a periodic box, for minimum image distances, or `None` for the period of the space.

        epsilon:
            The depth of the potential well.

        sigma:
            The distance at which the potential is zero.

    Examples:
        >>> from skcomplex.models.physics.common.particles import LennardJonesInteraction
        >>> LennardJonesInteraction('lj', first='gas', second='gas', cutoff=2.5).sigma
        1.0
    """

    def __init__(
        self,
        label: str,
        first: str,
        second: str,
        cutoff: float = 1.0,
        position: str = 'x',
        *,
        write: str = 'x',
        divisor: str | None = None,
        period: float | None = None,
        epsilon: 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,
            first,
            second,
            cutoff=cutoff,
            position=position,
            write=write,
            potential=self._get_lennard_jones_energy,
            divisor=divisor,
            period=period,
        )
        self.epsilon = epsilon
        self.sigma = sigma

    @staticmethod
    def _get_lennard_jones_energy(s: SimpleNamespace, r: Expression) -> Expression:
        """Return `4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)`."""
        return 4 * s.epsilon * ((s.sigma / r) ** 12 - (s.sigma / r) ** 6)