Source code for pyretis.core._system_inf

"""The canonical PyRETIS ``System`` class.

Originally ported from ``infretis.classes.system`` as a lightweight
value type for phasepoint storage; A3.2c (2026-05-27) absorbed the
former pyretis-native "harness" System into this class so PyRETIS
has a single, unified type.

Both :mod:`pyretis.core.system` and :mod:`pyretis.core.snapshot`
re-export the class defined here. The underscore filename is
historical (the module pre-dated the unification); the class
itself is no longer "internal" or "inf-flavour."

Construction
~~~~~~~~~~~~

``System()`` -- bare phasepoint-storage shape (all the per-state
attributes None or empty until written).

``System(units=..., temperature=..., box=..., post_setup=...)`` --
live-System factory shape (matches what
:mod:`pyretis.setup.createsystem` expects).

Bridge properties
~~~~~~~~~~~~~~~~~

``pos`` / ``vel`` / ``vpot`` / ``ekin`` / ``mass`` / ``imass`` /
``force`` / ``virial`` / ``config`` / ``vel_rev`` route through
``self.particles`` when attached (harness-style behaviour for the
factory-built case) and through per-instance backing storage when
no particles are attached (pure-snapshot use as path phasepoints).

Behaviour
~~~~~~~~~

The methods on this class are thin wrappers around the canonical
free functions in :mod:`pyretis.core.forcefieldhelp`,
:mod:`pyretis.core.systemhelp`, and :mod:`pyretis.inout.restart`.
New code should call those directly; the methods exist for
backwards compatibility with the legacy ``system.X()`` shape.
"""

from __future__ import annotations

import logging
from typing import Any, Dict, List, Optional, Tuple

import numpy as np

# The forcefield/restart/system helpers below are imported lazily inside
# the methods that use them: ``System`` is the foundational object those
# helper modules (and the factory/restart/snapshot layers) build on, so
# importing them at module top would create import cycles. Keep the lazy
# pattern and silence import-outside-toplevel for the whole module.
# pylint: disable=import-outside-toplevel

logger = logging.getLogger(__name__)
logger.addHandler(logging.NullHandler())


[docs]class System: """Define the representation of System snapshots. The system is intended for sharing the current configuration (of a frame in a path) between different methods, but the same class also serves as the canonical live System (post A3.2c) -- it accepts the harness-style construction kwargs and carries the ``particles`` / ``forcefield`` / ``post_setup`` attributes the factory layer needs. """ # The live System aggregates the full phase-point state plus the # harness-compat fields, so it legitimately carries many attributes # and public accessors -- disable the size heuristics for the class. # pylint: disable=too-many-instance-attributes # pylint: disable=too-many-public-methods
[docs] def __init__( self, units: str = 'lj', temperature: Optional[float] = None, box: Any = None, post_setup: Optional[List[Tuple[str, tuple]]] = None, ) -> None: """Initiate class. Parameters ---------- units : str, optional Unit system identifier (e.g. ``"lj"``, ``"reduced"``, ``"openmm"``). Used by helpers such as :func:`pyretis.core.systemhelp.get_boltzmann`. temperature : float, optional Target temperature. Stored under the ``temperature['set']`` key for harness compatibility; ``dof`` and ``beta`` are initialised to ``None``. box : object, optional Either a :class:`pyretis.core.box.Box` instance (harness shape) or a 3x3 ``np.ndarray`` (snapshot shape). When ``None`` the legacy snapshot 3x3 zeros default is used. post_setup : list of (str, tuple), optional Deferred-setup callback list. Defaults to ``[]``. """ # --- phasepoint-storage core --- # pos / vel / vpot / ekin / mass / imass / force / virial / # config / vel_rev are exposed as bridge properties below # so they delegate to ``self.particles`` whenever it's set # (post-A3.2c: factory builds System then assigns particles; # engines integrate ``system.pos += ...`` and that must show # through to particles.pos for downstream consumers). # Backing storage holds the value when no particles attached # (legacy snapshot use as a pure phasepoint). self._config: Tuple[str, int] = ("", -1) # ``order`` default was ``[-float("nan")]`` in the pure # snapshot. Switched to ``None`` to match the harness's # construction-time default; the inf code path that read # ``frame.order[0]`` always overwrites order from disk # before reading. self.order: Optional[List[float]] = None # ``pos`` / ``vel`` backing storage default ``None`` (matches # the harness bridge surface where ``system.pos`` is ``None`` # when no particles are attached). The pure-snapshot inf # flow always sets these before reading. self._pos: Optional[np.ndarray] = None self._vel: Optional[np.ndarray] = None self._vel_rev: bool = False self._ekin: Optional[float] = None self._vpot: Optional[float] = None self.etot: Optional[float] = None self.temp: Optional[float] = None self._force: Optional[np.ndarray] = None self._virial: Optional[np.ndarray] = None self._mass: Optional[np.ndarray] = None self._imass: Optional[np.ndarray] = None # --- harness-compat construction surface --- self.units = units # Harness temperature is a 3-key dict; preserve that shape so # downstream callers reading ``temperature['beta']`` / # ``temperature['dof']`` don't break. self.temperature: Dict[str, Any] = { 'set': temperature, 'dof': None, 'beta': None } # Box: accept Box object, ndarray, or None (the harness # default; the factory may attach a Box later). self.box = box self.particles = None self.forcefield = None self.post_setup: List[Tuple[str, tuple]] = ( list(post_setup) if post_setup is not None else [] ) # Mirror the harness ``__init__``: when the box exposes a # ``periodic`` attribute (harness :class:`.Box`), subtract one # dof per periodic dimension. Silent no-op for the snapshot's # ndarray box (no ``periodic`` attribute). self._adjust_dof_according_to_box() # Match the harness: compute the thermodynamic beta at # construction time so downstream callers reading # ``system.temperature['beta']`` see a real value. self.temperature['beta'] = self.calculate_beta()
[docs] def extra_setup(self) -> None: """Run post-construction callbacks registered in ``post_setup``.""" for func_name, args in self.post_setup: func = getattr(self, func_name, None) if func is not None: func(*args)
# --- Convenience wrappers around the free-function helpers # (forcefieldhelp, systemhelp, inout.restart). These keep the # familiar ``system.X()`` shape for callers (tests, user scripts, # tutorials) while the actual behaviour lives in the helpers so # the snapshot stays free of forcefield/restart imports at # class-definition time.
[docs] def update_force(self): """See :func:`pyretis.core.forcefieldhelp.update_force`.""" from pyretis.core.forcefieldhelp import update_force return update_force(self, self.forcefield)
[docs] def potential(self): """See :func:`pyretis.core.forcefieldhelp.update_potential`.""" from pyretis.core.forcefieldhelp import update_potential return update_potential(self, self.forcefield)
[docs] def potential_and_force(self): """Update and return the system potential energy and forces. See :func:`pyretis.core.forcefieldhelp.update_potential_and_force`. """ from pyretis.core.forcefieldhelp import update_potential_and_force return update_potential_and_force(self, self.forcefield)
[docs] def evaluate_force(self): """See :func:`pyretis.core.forcefieldhelp.evaluate_force`.""" from pyretis.core.forcefieldhelp import evaluate_force return evaluate_force(self, self.forcefield)
[docs] def evaluate_potential(self): """See :func:`pyretis.core.forcefieldhelp.evaluate_potential`.""" from pyretis.core.forcefieldhelp import evaluate_potential return evaluate_potential(self, self.forcefield)
[docs] def evaluate_potential_and_force(self): """Evaluate and return the system potential energy and forces. See :func:`pyretis.core.forcefieldhelp.evaluate_potential_and_force`. """ from pyretis.core.forcefieldhelp import ( evaluate_potential_and_force, ) return evaluate_potential_and_force(self, self.forcefield)
[docs] def rescale_velocities(self, energy, external=False): """Rescale system velocities to the requested energy. See :func:`pyretis.core.forcefieldhelp.rescale_system_velocities`. """ from pyretis.core.forcefieldhelp import rescale_system_velocities rescale_system_velocities( self, energy, forcefield=self.forcefield, external=external, )
[docs] def update_box(self, length): """See :func:`pyretis.core.systemhelp.update_system_box`.""" from pyretis.core.systemhelp import update_system_box update_system_box(self, length)
[docs] def get_dim(self): """See :func:`pyretis.core.systemhelp.get_system_dim`.""" from pyretis.core.systemhelp import get_system_dim return get_system_dim(self)
[docs] def get_boltzmann(self): """See :func:`pyretis.core.systemhelp.get_boltzmann`.""" from pyretis.core.systemhelp import get_boltzmann return get_boltzmann(self)
[docs] def generate_velocities(self, rgen=None, seed=0, momentum=True, temperature=None, distribution='maxwell'): """See :func:`pyretis.core.systemhelp.generate_system_velocities`.""" # pylint: disable=too-many-arguments,too-many-positional-arguments from pyretis.core.systemhelp import generate_system_velocities generate_system_velocities( self, rgen=rgen, seed=seed, momentum=momentum, temperature=temperature, distribution=distribution, )
[docs] def calculate_temperature(self): """See :func:`pyretis.core.systemhelp.calculate_system_temperature`.""" from pyretis.core.systemhelp import calculate_system_temperature return calculate_system_temperature(self)
[docs] def calculate_beta(self, temperature=None): """Thermodynamic beta for this system.""" from pyretis.core.particlefunctions import calculate_beta if temperature is None: if self.temperature['set'] is None: return None temperature = self.temperature['set'] return calculate_beta(temperature, self.units)
[docs] def restart_info(self): """See :func:`pyretis.inout.restart.system_restart_info`.""" from pyretis.inout.restart import system_restart_info return system_restart_info(self)
[docs] def load_restart_info(self, info): """See :func:`pyretis.inout.restart.load_system_restart`.""" from pyretis.inout.restart import load_system_restart load_system_restart(self, info)
[docs] def reverse_velocities(self): """Flip the sign of the velocities. Mirrors the harness System's :meth:`.reverse_velocities`: when ``particles`` is attached (harness flavour) delegate through it (handles ``ParticlesExt.vel_rev`` flag toggle vs regular ``Particles.vel`` flip). Otherwise (pure snapshot, external-engine restart-on-disk style) toggle ``vel_rev``. """ if self.particles is not None: self.particles.reverse_velocities() else: self.vel_rev = not self.vel_rev
@property def particle_type(self): """Engine flavour ('internal'/'external') of the attached particles.""" if self.particles is None: return None return getattr(self.particles, 'particle_type', None)
[docs] def particles_restart_info(self): """Return restart info for the attached particles object.""" if self.particles is None: return {} return self.particles.restart_info()
[docs] def __eq__(self, other): """Compare two system objects. Mirrors the harness ``System.__eq__``: ignores the order parameter (depends on the order-function choice), compares units / post_setup / box / particles / forcefield, and the ``set`` / ``beta`` / ``dof`` keys of the temperature dict. """ from pyretis.core.common import compare_objects, numpy_allclose attrs = ('units', 'post_setup', 'box', 'particles', 'forcefield') check = compare_objects(self, other, attrs, numpy_attrs=None) # Temperature is a dict on both sides; both sides may have # had ``temperature`` deleted (unusual but exercised by # tests). If either side is missing it, the systems are not # equal. self_has_temp = hasattr(self, 'temperature') other_has_temp = hasattr(other, 'temperature') if self_has_temp != other_has_temp: return False if self_has_temp: check = check and len(self.temperature) == len( other.temperature ) for key in ('set', 'beta'): check = check and ( self.temperature[key] == other.temperature[key] ) check = check and numpy_allclose( self.temperature['dof'], other.temperature['dof'], ) return check
[docs] def __ne__(self, other): """Check if two systems are not equal.""" return not self == other
[docs] def __str__(self): """Print basic info about the system.""" msg = ['PyRETIS System', f'Order parameter: {self.order}', f'Box: {self.box}', f'Particles: {self.particles}'] return '\n'.join(msg)
__hash__ = None # mutable: not hashable
[docs] def add_particle(self, pos, vel=None, force=None, mass=1.0, name='?', ptype=0): """Add a particle. Requires ``self.particles`` to be set first.""" # pylint: disable=too-many-arguments,too-many-positional-arguments if vel is None: vel = np.zeros_like(pos) if force is None: force = np.zeros_like(pos) self.particles.add_particle(pos, vel, force, mass=mass, name=name, ptype=ptype)
@property def npart(self): """Number of particles (delegates to ``self.particles.npart``).""" if self.particles is None: return 0 return self.particles.npart # --- Bridge properties for the phasepoint state. # Read/write through particles when attached (harness-style # behaviour after the factory builds the System then assigns # particles); fall back to per-instance backing storage when no # particles are attached (inf-flavour pure-snapshot phasepoints # serialised from disk). @property def pos(self): """Particle positions (routed to ``particles`` when attached).""" if self.particles is not None: return self.particles.pos return self._pos @pos.setter def pos(self, value): if self.particles is not None: self.particles.pos = value else: self._pos = value @property def vel(self): """Particle velocities (routed to ``particles`` when attached).""" if self.particles is not None: return self.particles.vel return self._vel @vel.setter def vel(self, value): if self.particles is not None: self.particles.vel = value else: self._vel = value @property def vpot(self): """Potential energy (routed to ``particles`` when attached).""" if self.particles is not None: return self.particles.vpot return self._vpot @vpot.setter def vpot(self, value): if self.particles is not None: self.particles.vpot = value else: self._vpot = value @property def ekin(self): """Kinetic energy (routed to ``particles`` when attached).""" if self.particles is not None: return self.particles.ekin return self._ekin @ekin.setter def ekin(self, value): if self.particles is not None: self.particles.ekin = value else: self._ekin = value @property def force(self): """Per-particle forces (routed to ``particles`` when attached).""" if self.particles is not None: return getattr(self.particles, 'force', None) return self._force @force.setter def force(self, value): if self.particles is not None: self.particles.force = value else: self._force = value @property def virial(self): """Virial tensor (routed to ``particles`` when attached).""" if self.particles is not None: return getattr(self.particles, 'virial', None) return self._virial @virial.setter def virial(self, value): if self.particles is not None: self.particles.virial = value else: self._virial = value @property def mass(self): """Per-particle masses (routed to ``particles`` when attached).""" if self.particles is not None: return getattr(self.particles, 'mass', None) return self._mass @mass.setter def mass(self, value): if self.particles is not None: self.particles.mass = value else: self._mass = value @property def imass(self): """Inverse particle masses (routed to ``particles`` when attached).""" if self.particles is not None: return getattr(self.particles, 'imass', None) return self._imass @imass.setter def imass(self, value): if self.particles is not None: self.particles.imass = value else: self._imass = value @property def vel_rev(self): """Velocity-reversal flag (routed to ``particles`` when attached).""" if self.particles is not None: return getattr(self.particles, 'vel_rev', self._vel_rev) return self._vel_rev @vel_rev.setter def vel_rev(self, value): if self.particles is not None: # ParticlesExt has vel_rev; regular Particles doesn't. # setattr is safe (Python allows attribute creation). self.particles.vel_rev = value else: self._vel_rev = value @property def config(self): """External-engine trajectory pointer (file, frame_idx). When ``self.particles`` is a :class:`.ParticlesExt`, routes through its ``config`` attribute. For pure phasepoints (no particles) the backing ``self._config`` tuple is used. """ if self.particles is not None: return getattr(self.particles, 'config', None) return self._config @config.setter def config(self, value): if self.particles is not None: self.particles.config = value else: self._config = value
[docs] def adjust_dof(self, dof): """Subtract ``dof`` from the system's degree-of-freedom count. Mirrors the harness method: accumulates into ``self.temperature['dof']`` (initialised lazily). """ if self.temperature.get('dof') is None: self.temperature['dof'] = np.array(dof) else: self.temperature['dof'] += np.array(dof)
[docs] def _adjust_dof_according_to_box(self): """Walk the box's ``periodic`` attribute and subtract dofs. Subtracts one dof per periodic dimension (harness Box only). Snapshot ndarray boxes have no ``periodic`` -- silently no-ops via the AttributeError fallback. """ try: dof = [] all_false = True for peri in self.box.periodic: dof.append(1 if peri else 0) all_false = all_false and not peri if not all_false: self.adjust_dof(dof) return True except AttributeError: return False
[docs] def copy(self) -> System: # pylint: disable=too-many-branches """Return a copy of this system. Deep-copies the phasepoint-storage state (pos/vel/box/order/ energies/etc.) and the per-instance harness-compat fields (units/temperature/particles/post_setup). ``forcefield`` is shared by reference -- matches the harness ``System.copy()`` contract: if the force field changes, every system copy must see the change. """ # copy() legitimately seeds a freshly-built sibling System's # private phasepoint storage (_config/_pos/_vel/_force/...); # those attributes stay private by design, so silence the # protected-access check for the whole method. # pylint: disable=protected-access system_copy = self.__class__() # config is a bridge property -- setting it routes to # particles.config when particles is attached, else to # self._config. The particles aren't attached yet on the # fresh copy, so this writes to _config. self_config = self.config if self_config is not None: system_copy._config = tuple(self_config) else: system_copy._config = None if self.order is None: system_copy.order = None elif isinstance(self.order, (list, tuple)): system_copy.order = list(self.order) else: # Scalar or numpy array -- copy via the copy module to # preserve type (numpy.float64 etc.). from copy import copy as _copy system_copy.order = _copy(self.order) system_copy.etot = self.etot system_copy.temp = self.temp # Box may be None, ndarray, or a harness Box object with a # .copy() method. if self.box is None: system_copy.box = None elif isinstance(self.box, np.ndarray): system_copy.box = np.array(self.box, copy=True) else: box_copy = getattr(self.box, 'copy', None) system_copy.box = box_copy() if callable(box_copy) else self.box # ``temperature`` may have been deleted by an unusual caller # (some tests del-attribute it to exercise the "missing # attribute" code path). When that happens, leave the # default ``temperature`` dict that ``__init__`` already # placed on ``system_copy`` so the two systems compare # unequal (matches the harness copy contract). try: system_copy.temperature = dict(self.temperature) except AttributeError: logger.warning( 'Missing attribute "temperature" when copying system', ) # Harness-compat fields. system_copy.units = self.units system_copy.post_setup = list(self.post_setup) if self.particles is None: system_copy.particles = None else: particles_copy = getattr(self.particles, 'copy', None) if callable(particles_copy): # Guarded by callable() above; pylint can't infer the # getattr result is callable here. # pylint: disable-next=not-callable system_copy.particles = particles_copy() else: # Fallback (snapshot or unfamiliar type): share ref. system_copy.particles = self.particles system_copy.forcefield = self.forcefield # Phasepoint-state copy. # When particles is attached, ``particles.copy()`` above # already deep-copied pos/vel/force/virial/mass/imass. So we # ONLY copy the backing storage when no particles are # attached (pure-snapshot use). Avoids accidentally writing # bridge values back through the setter after the particles # copy already populated them. if self.particles is None: system_copy._pos = ( None if self._pos is None else np.array(self._pos, copy=True) ) system_copy._vel = ( None if self._vel is None else np.array(self._vel, copy=True) ) system_copy._vel_rev = self._vel_rev system_copy._ekin = self._ekin system_copy._vpot = self._vpot system_copy._force = ( None if self._force is None else np.array(self._force, copy=True) ) system_copy._virial = ( None if self._virial is None else np.array(self._virial, copy=True) ) system_copy._mass = ( None if self._mass is None else np.array(self._mass, copy=True) ) system_copy._imass = ( None if self._imass is None else np.array(self._imass, copy=True) ) return system_copy
[docs] def set_pos(self, pos: Tuple[str, int]) -> None: """Set positions for the particles.""" self.config = (pos[0], pos[1])