Source code for espressomd.system

#
# Copyright (C) 2013-2026 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

import numpy as np

# pylint: disable=unused-import
from . import accumulators
from . import analyze
from . import bond_breakage
from . import cell_system
from . import cuda_init
from . import collision_detection
from . import comfixed
from . import constraints
from . import electrokinetics
from . import electrostatics
from . import magnetostatics
from . import galilei
from . import interactions
from . import integrate
from . import lb
from . import lees_edwards
from . import particle_data
from . import thermostat
# pylint: enable=unused-import

from .code_features import has_features, assert_features
from .script_interface import script_interface_register, ScriptInterfaceHelper


[docs] @script_interface_register class System(ScriptInterfaceHelper): """ The ESPResSo system class. Attributes ---------- analysis: :class:`espressomd.analyze.Analysis` auto_update_accumulators: :class:`espressomd.accumulators.AutoUpdateAccumulators` bond_breakage: :class:`espressomd.bond_breakage.BreakageSpecs` bonded_inter: :class:`espressomd.interactions.BondedInteractions` cell_system: :class:`espressomd.cell_system.CellSystem` collision_detection: :class:`espressomd.collision_detection.CollisionDetection` comfixed: :class:`espressomd.comfixed.ComFixed` constraints: :class:`espressomd.constraints.Constraints` cuda_init_handle: :class:`espressomd.cuda_init.CudaInitHandle` galilei: :class:`espressomd.galilei.GalileiTransform` integrator: :class:`espressomd.integrate.IntegratorHandle` lees_edwards: :class:`espressomd.lees_edwards.LeesEdwards` non_bonded_inter: :class:`espressomd.interactions.NonBondedInteractions` part: :class:`espressomd.particle_data.ParticleList` thermostat: :class:`espressomd.thermostat.Thermostat` ekcontainer: :class:`espressomd.electrokinetics.EKContainer` EK system (diffusion-advection-reaction models). box_l: (3,) array_like of :obj:`float` Dimensions of the simulation box. periodicity: (3,) array_like of :obj:`bool` System periodicity in ``[x, y, z]``, ``False`` for no periodicity in this direction, ``True`` for periodicity min_global_cut : :obj:`float` Minimal interaction cutoff. Methods ------- number_of_particles() Count the number of particles of a given type. Parameters ---------- type : :obj:`int` (:attr:`~espressomd.particle_data.ParticleHandle.type`) Particle type to count the number for. Returns ------- :obj:`int` The number of particles which have the given type. rotate_system() Rotate the particles in the system about the center of mass. If ``ROTATION`` is activated, the internal rotation degrees of freedom are rotated accordingly. Parameters ---------- phi : :obj:`float` Angle between the z-axis and the rotation axis. theta : :obj:`float` Rotation of the axis around the y-axis. alpha : :obj:`float` How much to rotate """ _so_name = "System::System" _so_creation_policy = "GLOBAL" _so_bind_methods = ( "number_of_particles", "rotate_system") def __init__(self, **kwargs): if "sip" in kwargs: super().__init__(**kwargs) self._setup_atexit() else: super().__init__(_regular_constructor=True, **kwargs) if has_features("CUDA"): self.cuda_init_handle = cuda_init.CudaInitHandle() # lock class self.call_method("lock_system_creation") self._setup_atexit() def _setup_atexit(self): import atexit def session_shutdown(): self.call_method("session_shutdown") atexit.register(session_shutdown) def __reduce__(self): so_callback, so_callback_args = super().__reduce__() return (System._restore_object, (so_callback, so_callback_args, self.__getstate__())) @classmethod def _restore_object(cls, so_callback, so_callback_args, state): so = so_callback(*so_callback_args) so.__setstate__(state) return so def __getstate__(self): return {} def __setstate__(self, params): # initialize Python-only members for property_name in params.keys(): System.__setattr__(self, property_name, params[property_name]) self.call_method("lock_system_creation") @property def force_cap(self): """ If > 0, the magnitude of the force on the particles are capped to this value. Type: :obj:`float` """ return self.integrator.force_cap @force_cap.setter def force_cap(self, value): self.integrator.force_cap = value @property def time(self): """ Total simulation time. Type: :obj:`float` """ return self.integrator.time @time.setter def time(self, value): self.integrator.time = value @property def time_step(self): """ MD time step. Type: :obj:`float` """ return self.integrator.time_step @time_step.setter def time_step(self, value): self.integrator.time_step = value @property def max_cut_nonbonded(self): """ Maximal cutoff for non-bonded interactions. Type: :obj:`float` """ return self.cell_system.max_cut_nonbonded @property def max_cut_bonded(self): """ Maximal cutoff for bonded interactions. Type: :obj:`float` """ return self.cell_system.max_cut_bonded @property def lb(self): """ Lattice-Boltzmann method. Type: :class:`espressomd.lb.LBFluid` """ return self.lbcontainer.solver @lb.setter def lb(self, lb): self.lbcontainer.solver = lb
[docs] def change_volume_and_rescale_particles(self, d_new, dir="xyz"): """Change box size and rescale particle coordinates. Parameters ---------- d_new : :obj:`float` New box length dir : :obj:`str`, optional Coordinate to work on, ``"x"``, ``"y"``, ``"z"`` or ``"xyz"`` for isotropic. Isotropic makes the box cubic. """ coord = {"x": 0, "y": 1, "z": 2, 0: 0, 1: 1, 2: 2, "xyz": 3}.get(dir) if coord is None: raise ValueError( 'Usage: change_volume_and_rescale_particles(<L_new>, [{ "x" | "y" | "z" | "xyz" }])') self.call_method("rescale_boxl", length=d_new, coord=coord)
[docs] def volume(self): """Return volume of the cuboid box. """ return float(np.prod(self.box_l))
[docs] def distance(self, p1, p2): """Return the scalar distance between particles, between a particle and a point or between two points, respecting periodic boundaries. Parameters ---------- p1 : :class:`~espressomd.particle_data.ParticleHandle` or (3,) array_like of :obj:`float` First particle or position. p2 : :class:`~espressomd.particle_data.ParticleHandle` or (3,) array_like of :obj:`float` Second particle or position. """ return np.linalg.norm(self.distance_vec(p1, p2))
[docs] def distance_vec(self, p1, p2): """Return the distance vector between particles, between a particle and a point or between two points, respecting periodic boundaries. Parameters ---------- p1 : :class:`~espressomd.particle_data.ParticleHandle` or (3,) array_like of :obj:`float` First particle or position. p2 : :class:`~espressomd.particle_data.ParticleHandle` or (3,) array_like of :obj:`float` Second particle or position. """ if isinstance(p1, particle_data.ParticleHandle): pos1 = p1.pos_folded else: pos1 = p1 if isinstance(p2, particle_data.ParticleHandle): pos2 = p2.pos_folded else: pos2 = p2 return self.call_method("distance_vec", pos1=pos1, pos2=pos2)
[docs] def velocity_difference(self, p1, p2): """ Return the velocity difference between two particles, considering Lees-Edwards boundary conditions, if active. Parameters ---------- p1 : :class:`~espressomd.particle_data.ParticleHandle` p2 : :class:`~espressomd.particle_data.ParticleHandle` """ return self.call_method("velocity_difference", pos1=p1.pos_folded, pos2=p2.pos_folded, v1=p1.v, v2=p2.v)
[docs] def auto_exclusions(self, distance): """ Add exclusions between particles that are bonded. This only considers pair bonds. Requires feature ``EXCLUSIONS``. Parameters ---------- distance : :obj:`int` Bond distance up to which the exclusions should be added. """ assert_features("EXCLUSIONS") self.part.auto_exclusions(distance=distance)
[docs] def setup_type_map(self, *args, **kwargs): # pylint: disable=unused-argument raise NotImplementedError( "system.setup_type_map() was removed in release 5.1")