# Copyright (c) 2026, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""This module defines the class for interfacing LAMMPS.
Important classes defined here
------------------------------
LAMMPSEngineSteps (:py:class:`.LAMMPSEngineSteps`)
The class responsible for interfacing with LAMMPS by relaunching the
executable once per order-parameter point (the slower, relaunch
implementation; the canonical streaming engine lives in
:py:mod:`pyretis.engines.lammps`).
"""
import logging
import os
import shlex
import numpy as np
from pyretis.core.common import counter
from pyretis.core.engine_time import validate_engine_timestep
from pyretis.core.velocity import is_aimless_velocity_setting
from pyretis.engines.external import ExternalMDEngine
from pyretis.inout.fileio import FileIO
from pyretis.inout.settings import look_for_input_files
logger = logging.getLogger(__name__) # pylint: disable=invalid-name
logger.addHandler(logging.NullHandler())
LAMMPS_ORDER_MODE = 'lammps'
PYRETIS_ORDER_MODE = 'pyretis'
ORDER_MODES = {LAMMPS_ORDER_MODE, PYRETIS_ORDER_MODE}
# Let us define the LAMMPS commands PyRETIS will make use of:
# 1) The commands related to the order parameter.
# This first command is to ensure that the order parameter is
# calculated also for cases where we do not run any MD steps. This
# can for instance be when we are generating velocities:
ORDER_HACK = """run 0 every 1 "print '{}' file {} screen no" # :-)"""
# The parameter is here the variables to print.
# The next command is outputting the order parameter:
ORDER_FIX = 'fix order_output all print {} "{}" append {} screen no'
# The parameters are: frequency, variables to print and the file name.
# 2) The commands related to thermodynamic output:
# Define a thermo style for LAMMPS:
THERMO_STYLE = 'thermo_style custom step temp press pe ke etotal'
# 3) The commands related to trajectory output and intput:
TRAJ_DUMP = 'dump {} all custom {} {} id {}'
# The three parameters are here: dump-id, frequency, output file name
# and the format for what to print. What to print depends on the number
# dimensions we consider and this is defined below:
TRAJ_OUT = {
1: 'x vx ix',
2: 'x y vx vy ix iy',
3: 'x y z vx vy vz ix iy iz',
}
# Modify the format for the output by increasing the number of digits
# used for positions and velocities:
TRAJ_FMT = 'dump_modify {} format {} %23.16g'
# The parameter is here the dump-id and the column to apply it to.
# Add a command to undump the trajectory
TRAJ_UNDUMP = 'undump traj_pyretis'
# We also need to read in trajectories and/or snapshots:
READ_DUMP = 'read_dump {} {} {} box yes'
# The parameters are: the filename and the index in that file.
# The next command will read in LAMMPS restart files:
READ_RESTART = 'read_restart {}'
# The parameter is here the file name.
# 4) The commands for generating velocities:
GENERATE_VEL = 'velocity all create ${SET_TEMP} ${VEL_SEED} dist gaussian'
# The parameters are here the temperature and the random seed, which
# are actually set elsewhere in the LAMMPS input script.
# The next command sets the seed to use:
GENERATE_VEL_SEED = 'variable VEL_SEED equal {}'
# The parameter is here the seed to use.
# 5) The commands for stopping a simulation based on the order parameter
# and the given interface positions:
INTERFACE_RIGHT_FIX = 'fix op_stop_right all halt {} v_op_1 > {}'
# The parameters are here the frequency and the right interface.
INTERFACE_LEFT_FIX = 'fix op_stop_left all halt {} v_op_1 < {}'
# The parameters are here the frequency and the left interface.
# 6) The commands for rescaling energies if this is requested:
RESCALE_ENERGY = [
'variable ke_rescale equal "v_SET_ENERGY - pe"',
'variable alpha equal "sqrt(v_ke_rescale / ke)"',
]
# Here we need different commands, depending on the number
# of dimensions we are simulating. These are hard-coded below:
RESCALE_DIM = {
3: [
'variable vx atom "v_alpha * vx"',
'variable vy atom "v_alpha * vy"',
'variable vz atom "v_alpha * vz"',
'velocity all set v_vx v_vy v_vz',
],
2: [
'variable vx atom "v_alpha * vx"',
'variable vy atom "v_alpha * vy"',
'velocity all set v_vx v_vy 0.0',
],
1: [
'variable vx atom "v_alpha * vx"',
'velocity all set v_vx 0.0 0.0',
],
}
# Map LAMMPS thermo column labels (matching THERMO_STYLE
# ``step temp press pe ke etotal``) to the PyRETIS-canonical thermo
# keys used by EnergyFormatter and ThermoTableFormatter. Without this
# rename the formatters look up ``vpot``/``ekin``/etc. and fall back
# to NaN, even though LAMMPS reported real values.
THERMO_KEY_MAP = {
'PotEng': 'vpot',
'KinEng': 'ekin',
'TotEng': 'etot',
'Temp': 'temp',
'Press': 'press',
'Step': 'step',
}
[docs]def rename_lammps_energy_keys(energy):
"""Rename LAMMPS thermo column labels to PyRETIS conventions."""
return {THERMO_KEY_MAP.get(key, key): val for key, val in energy.items()}
# 7) Commands for reversing velocities:
# First the command for reversing:
VEL_REV = 'variable v{0} atom -v{0}'
# The argument is here to dimension to apply this to, e.g. "x".
# Then, the command to actually set the velocities.
VEL_SET = 'velocity all set {} {} {}'
# The arguments are the velocities to set, either v_vx, v_vy, v_vz
# or "NULL".
[docs]def read_lammps_log(filename):
"""Read some info from a LAMMPS log file.
In particular, this method is used to read the thermodynamic
output from a simulation (e.g. potential and kinetic energies).
Parameters
----------
filename : string
The path to the LAMMPS log file.
Returns
-------
out : dict
A dict containing the data we found in the file.
"""
energy_keys = []
energy_data = {}
read_energy = False
with open(filename, 'r', encoding='utf-8') as logfile:
for lines in logfile:
if lines.lstrip().startswith('Step'):
# Assume that this is the start of the thermo output.
# Discard every previously read section, including keys
# absent from this header, so only the final thermo
# section survives. By construction of the LAMMPS input
# file that final section is the MD run; rebuilding the
# dict from scratch prevents stray columns from earlier
# setup-time ``run 0`` blocks (e.g. a one-line Density
# readout) from lingering as ragged, shorter-than-frame
# arrays.
energy_keys = [i.strip() for i in lines.strip().split()]
energy_data = {key: [] for key in energy_keys}
read_energy = True
continue
if lines.startswith('Loop time'):
# Assume this marks the end of the thermo output.
read_energy = False
energy_keys = []
continue
if read_energy and energy_keys:
# Assume that we are reading energies.
try:
data = [float(i) for i in lines.strip().split()]
except ValueError:
# Some text has snuck into the thermo output,
# ignore it:
continue
for key, value in zip(energy_keys, data):
if key == 'Step':
energy_data[key].append(int(value))
else:
energy_data[key].append(value)
for key, val in energy_data.items():
energy_data[key] = np.array(val)
out = {'energy': energy_data}
return out
[docs]def _parse_dump_box(bounds):
"""Return orthogonal box lengths from LAMMPS dump bounds."""
box = []
for bound in bounds:
if len(bound) < 2:
raise ValueError('LAMMPS dump box bounds must contain lo and hi')
box.append(bound[1] - bound[0])
return np.array(box, dtype=float)
[docs]def read_lammps_dump(filename, dimension):
"""Read custom LAMMPS dump frames written by the PyRETIS engine.
The parser expects the custom dump produced by :data:`TRAJ_DUMP`, i.e.
one atom table per frame with an ``id`` column, position columns, velocity
columns, and optional image flags. Atoms are returned sorted by LAMMPS id.
Parameters
----------
filename : string
Path to a LAMMPS dump trajectory.
dimension : integer
Number of spatial dimensions to extract.
Returns
-------
frames : list of dict
Each frame contains ``step``, ``box``, ``pos`` and ``vel`` entries.
"""
with open(filename, 'r', encoding='utf-8') as dump:
return _parse_lammps_dump_lines(iter(dump), dimension, source=filename)
[docs]def _parse_lammps_dump_lines(lines, dimension, source='<lines>'):
"""Parse LAMMPS custom dump frames from an iterator of text lines.
Split out from :func:`read_lammps_dump` so the streaming engine
(:class:`pyretis.engines.lammps.LAMMPSRunner`) can feed it only the
already-completed frames of a trajectory that LAMMPS is still
writing. The parser is strict: a frame truncated mid-way raises, so
callers reading a growing file must hand it complete frames only.
Parameters
----------
lines : iterator of strings
The dump lines to parse.
dimension : integer
Number of spatial dimensions to extract.
source : string, optional
Name used in error messages (the originating file).
Returns
-------
frames : list of dict
Each frame contains ``step``, ``box``, ``pos`` and ``vel``.
"""
dimension = int(dimension)
dim_labels = ('x', 'y', 'z')[:dimension]
frames = []
for line in lines:
if not line.startswith('ITEM: TIMESTEP'):
continue
try:
step = int(next(lines).strip())
marker = next(lines).strip()
if marker != 'ITEM: NUMBER OF ATOMS':
raise ValueError
natoms = int(next(lines).strip())
marker = next(lines).strip()
if not marker.startswith('ITEM: BOX BOUNDS'):
raise ValueError
bounds = [
[float(value) for value in next(lines).split()]
for _ in range(3)
]
marker = next(lines).strip()
if not marker.startswith('ITEM: ATOMS'):
raise ValueError
except (StopIteration, ValueError) as exc:
msg = f'Could not parse LAMMPS dump frame in "{source}"'
raise ValueError(msg) from exc
columns = marker.split()[2:]
index = {name: i for i, name in enumerate(columns)}
if 'id' not in index:
raise ValueError('LAMMPS dump atom table must contain id')
box = _parse_dump_box(bounds)
atoms = []
for _ in range(natoms):
values = next(lines).split()
atom_id = int(values[index['id']])
pos = []
vel = []
for dim, label in enumerate(dim_labels):
coord_key = f'{label}u'
if coord_key in index:
coord = float(values[index[coord_key]])
elif label in index:
coord = float(values[index[label]])
image_key = f'i{label}'
if image_key in index:
coord += float(values[index[image_key]]) * box[dim]
else:
raise ValueError(
f'Missing position column for {label}'
)
vel_key = f'v{label}'
velocity = (
float(values[index[vel_key]])
if vel_key in index else 0.0
)
pos.append(coord)
vel.append(velocity)
atoms.append((atom_id, pos, vel))
atoms.sort(key=lambda item: item[0])
frames.append({
'step': step,
'box': box[:dimension],
'pos': np.array([atom[1] for atom in atoms], dtype=float),
'vel': np.array([atom[2] for atom in atoms], dtype=float),
})
return frames
[docs]def _add_order_fix(settings):
"""Add the fix for printing the order parameter for LAMMPS."""
lines = []
if 'order-fix' in settings:
logger.debug('Adding fix for order parameter to LAMMPS input.')
lines.append('# Order parameter fix:')
for line in settings['order-fix']:
lines.append(line)
return lines
[docs]def _add_traj_dump(settings):
"""Add the dump commands for storing the LAMMPS trajectory."""
lines = [
'# PyRETIS trajectory settings:',
'# PyRETIS requested the following trajectory output:',
]
lines.append(
TRAJ_DUMP.format(
'traj_pyretis',
settings['subcycles'],
settings['traj'],
TRAJ_OUT[settings['dimension']],
)
)
# Format for positions:
idx = 0
for i in range(int(settings['dimension'])):
idx = i + 2
lines.append(TRAJ_FMT.format('traj_pyretis', idx))
# Format for velocities:
for i in range(int(settings['dimension'])):
lines.append(TRAJ_FMT.format('traj_pyretis', idx + 1 + i))
return lines
[docs]def _add_stopping_condition(settings):
"""Add the interface stopping condition to the LAMMPS input."""
lines = []
if (
'interfaces' in settings and
settings.get('order_mode', LAMMPS_ORDER_MODE) == LAMMPS_ORDER_MODE
):
lines.append('# Stopping condition fix:')
lines.append(
INTERFACE_RIGHT_FIX.format(
settings['subcycles'],
settings['interfaces'][-1],
)
)
if abs(settings['interfaces'][0]) != float('inf'):
lines.append(
INTERFACE_LEFT_FIX.format(
settings['subcycles'],
settings['interfaces'][0],
)
)
return lines
[docs]def _add_generate_vel(settings):
"""Add generation of velocities to the LAMMPS input."""
lines = []
if 'generate_vel' in settings:
seed = settings['generate_vel'].get('seed', 1)
lines.append('# PyRETIS requested random velocities')
lines.append(GENERATE_VEL_SEED.format(seed))
lines.append(GENERATE_VEL)
lines.append('run 0')
if settings['generate_vel'].get('rescale', None):
for i in RESCALE_ENERGY:
lines.append(i)
for i in RESCALE_DIM[settings['dimension']]:
lines.append(i)
return lines
[docs]def system_to_lammps(system, reverse_velocities, dimension):
"""Convert a LAMMPS system into LAMMPS commands for loading it.
This method will convert a system created by LAMMPS into
LAMMPS commands, so that LAMMPS can read the given snapshot
and use it.
Parameters
----------
system : object like :py:class:`.System`
The system we are converting.
reverse_velocities : boolean
True if the velocities in the system are to be reversed.
dimension : integer
The number of dimensions used in the LAMMPS simulation.
Returns
-------
out : list of strings
The commands needed by LAMMPS to read the configuration.
"""
config = system.config
lammps = []
if config[0].endswith('.restart') and config[1] == 0:
# Assume that this is a LAMMPS restart file.
lammps.append('# PyRETIS requested a LAMMPS restart file:')
lammps.append(READ_RESTART.format(config[0]))
elif config[0].endswith('.data') and config[1] in (0, None):
# Assume that this is a LAMMPS data file AND that there
# is already a line in the LAMMPS input file reading this
# file. The motivation behind this is that we define this
# file type to be the initial configuration and that data files
# will never be used for something else.
lammps.append(f'# PyRETIS requested the LAMMPS data file: {config[0]}')
else:
# Assume that this is a LAMMPS trajectory.
lammps.append('# PyRETIS requested the following snapshot:')
lammps.append(
READ_DUMP.format(config[0], config[1], TRAJ_OUT[dimension])
)
lammps.append('# Reset time step:')
lammps.append('reset_timestep 0')
if reverse_velocities:
# We add commands for reversing the velocities:
lammps.append('# PyRETIS requested reversing of velocities:')
vel = ['NULL', 'NULL', 'NULL']
for i, dim in enumerate(['x', 'y', 'z'][:dimension]):
lammps.append(VEL_REV.format(dim))
vel[i] = f'v_v{dim}'
lammps.append(VEL_SET.format(*vel))
lammps.append('run 0')
return lammps
[docs]class LAMMPSEngineSteps(ExternalMDEngine):
"""
A class for interfacing LAMMPS.
Attributes
----------
lmp : string
The command for executing LAMMPS
input_path : string
The directory where the input files are stored.
input_files : dict of strings
The names of the input files.
"""
needs_order = False
default_units = None
[docs] @classmethod
def get_default_units(cls, settings=None):
"""Return the default unit system for the LAMMPS engine.
This is read from the ``units`` command in the LAMMPS input
template. If the template can not be inspected, or if it does
not define ``units``, an exception is raised.
Parameters
----------
settings : dict, optional
Full simulation settings or engine settings.
Returns
-------
out : string
The default unit system for the LAMMPS engine.
Raises
------
ValueError
If the LAMMPS input template can not be read, or if no
``units`` command is found in that file.
"""
if settings is None:
return super().get_default_units(settings)
engine_settings = settings.get('engine', settings)
input_path = engine_settings.get('input_path', '.')
exe_path = engine_settings.get('exe_path')
if 'simulation' in settings:
exe_path = settings['simulation'].get('exe_path', exe_path)
if exe_path is None:
exe_path = os.path.abspath('.')
template = os.path.join(exe_path, input_path, 'lammps.in')
if not os.path.isfile(template):
msg = f'Could not determine LAMMPS units: Missing "{template}"'
raise ValueError(msg)
for key, value in read_lammps_input(template):
if key.lower() == 'units':
return value.strip()
msg = f'Could not determine LAMMPS units from "{template}"'
raise ValueError(msg)
[docs] def __init__(self, lmp, input_path, subcycles, extra_files=None,
order_mode=LAMMPS_ORDER_MODE, timestep=None,
velocity_generation='engine', temperature=None):
"""Set up the LAMMPS engine.
Parameters
----------
lmp : string
The LAMMPS executable.
input_path : string
The absolute path to where the input files are stored.
subcycles : integer
The frequency of output of data by LAMMPS.
extra_files : list of strings, optional
Additional files needed to run the LAMMPS simulation.
order_mode : string, optional
Select where order parameters are evaluated. ``"lammps"`` keeps
the traditional ``order.in``/``v_op_1`` mode. ``"pyretis"`` reads
dumped LAMMPS frames and evaluates the PyRETIS order parameter and
collective variables in Python.
timestep : float, optional
The LAMMPS time step is read from the LAMMPS input file; this
optional ``[engine]`` value is only used as a cross-check and
must equal the value in the input file (otherwise an error is
raised).
velocity_generation : string, optional
How shooting-move velocities are generated. ``"engine"``
(default) delegates to LAMMPS (``velocity create``);
``"maxwell"`` draws them in Python from a Maxwell-Boltzmann
distribution via the injected random generator (one reproducible
stream; ``'real'`` units; requires ``temperature``).
temperature : float, optional
The temperature (K) used to draw Maxwell-Boltzmann velocities
when ``velocity_generation == "maxwell"``.
"""
order_mode = str(order_mode).lower()
if order_mode not in ORDER_MODES:
msg = (
f'Unknown LAMMPS order mode "{order_mode}". '
f'Valid modes are: {", ".join(sorted(ORDER_MODES))}'
)
raise ValueError(msg)
# Note: We set time step to zero for now. This will be updated
# to the correct value when the timestep has been read from
# the input script:
super().__init__('LAMMPS Engine', 0.0, subcycles)
# Cumulative MD-step counter read by the infinite-swapping
# scheduler for subcycle accounting (unused by native runs).
self.steps = 0
self.ext = 'data'
self.lmp = lmp
self.input_path = os.path.abspath(input_path)
self.order_mode = order_mode
# Velocity-generation policy (4b-2b #4): "engine" delegates the
# shooting-velocity draw to LAMMPS (`velocity create`); "maxwell"
# draws it in Python from a Maxwell-Boltzmann distribution via the
# injected random generator (one reproducible stream, 'real' units).
# Default "engine" is unchanged.
if velocity_generation not in ('engine', 'maxwell'):
msg = ('Unknown velocity_generation "%s" '
'(expected "engine" or "maxwell")')
logger.error(msg, velocity_generation)
raise ValueError(msg % velocity_generation)
self.velocity_generation = velocity_generation
# k_B in LAMMPS 'real' units (kcal/(mol*K)); beta = 1 / (k_B T).
self.kb = 1.987204259e-3
self.temperature = temperature
self.beta = None
if temperature is not None:
self.beta = 1.0 / (self.kb * self.temperature)
if velocity_generation == 'maxwell' and temperature is None:
msg = ('velocity_generation = "maxwell" requires a '
'"temperature" ([engine] setting).')
logger.error(msg)
raise ValueError(msg)
# Define the files we require:
input_files = {
'conf': 'system.data',
'template': 'lammps.in', # The template file for LAMMPS.
}
if self.order_mode == LAMMPS_ORDER_MODE:
input_files['order'] = 'order.in'
# The user may choose to structure the LAMMPS input into
# several files. This is something we will not try to figure
# out, and we assume that the user knows what he/she is doing.
# We assume that these extra files will also be present and
# stored in the given input path. Now, we look for such files:
self.input_files = look_for_input_files(
self.input_path,
input_files,
)
self.run_files = {
'conf': self.input_files['conf'],
}
if 'order' in self.input_files:
self.run_files['order'] = self.input_files['order']
elif self.order_mode == PYRETIS_ORDER_MODE:
order_file = os.path.join(self.input_path, 'order.in')
if os.path.isfile(order_file):
self.input_files['order'] = order_file
self.run_files['order'] = order_file
if extra_files is not None:
extra = look_for_input_files(
self.input_path,
{f'file-{i}': val for i, val in enumerate(extra_files)}
)
for key, val in extra.items():
self.run_files[key] = val
# Read the LAMMPS input file:
self.settings = {
'dimension': 3,
'timestep': 0.0
}
for (key, val) in read_lammps_input(self.input_files['template']):
self.settings[key] = val
# Also set the timestep explicitly in case it is needed:
self.timestep = self.settings['timestep']
# The LAMMPS input owns the time step; if one was also given in the
# [engine] settings, it must agree with it.
validate_engine_timestep(timestep, self.timestep, 'LAMMPS')
[docs] def _make_order_fix(self, ordername):
"""Create a LAMMPS fix for the order parameter.
Note that we here make LAMMPS also output for step zero.
"""
order_settings = read_lammps_input(
os.path.join(self.input_path, self.input_files['order'])
)
ops = ['$(step)']
for key, val in order_settings:
if key == 'variable':
ops.append(f'${{{val.strip().split()[0]}}}')
txt = [
ORDER_HACK.format(' '.join(ops), ordername),
ORDER_FIX.format(self.subcycles, ' '.join(ops), ordername),
]
return txt
[docs] def read_order_parameters(self, ordername):
"""Read order parameters as calculated by LAMMPS.
We assume here that these can be found in the current execute
directory in a file named `ordername`, and further that they
contain the step number in the first column, followed by the
order parameters in following columns.
"""
orderfile = os.path.join(self.exe_dir, ordername)
order = np.loadtxt(orderfile, ndmin=2)
if order.shape[0] > 1 and order[0, 0] == order[1, 0]:
order = np.delete(order, 0, 0)
# Removing the first column, as this is not an order parameter,
# but the step number in the simulation. The usage of deleting
# along axis=1 is the motivation for forcing ndim=2 above.
order = np.delete(order, 0, 1)
return order
[docs] def read_energies(self, name):
"""Read energies obtained in a LAMMPS run.
Here, we assume that the energies can be found in a log file
in the current execute directory.
"""
logfile = os.path.join(self.exe_dir, f'{name}.log')
data = read_lammps_log(logfile)
return data['energy']
[docs] def _propagate_from(self, name, path, system, ens_set, msg_file,
reverse=False): # pragma: no cover
"""
Propagate with LAMMPS from the current system configuration.
Parameters
----------
name : string
A name to use for the trajectory we are generating.
path : object like :py:class:`.PathBase`
This is the path we use to fill in phase-space points.
system : object like :py:class:`.System`
The system object gives the initial state. The order
function is injected on the engine (``self.order_function``).
ens_set : dict
Ensemble settings. It contains:
* `interfaces`: list of floats
These interfaces define the stopping criterion.
msg_file : object like :py:class:`.FileIO`
An object we use for writing out messages that are useful
for inspecting the status of the current propagation.
reverse : boolean, optional
If True, the system will be propagated backward in time.
Returns
-------
success : boolean
This is True if we generated an acceptable path.
status : string
A text description of the current status of the propagation.
"""
return
[docs] def step(self, system, name):
"""Perform a single step with LAMMPS.
Parameters
----------
system : object like :py:class:`.System`
The system we are integrating.
name : string
To name the output files from the LAMMPS step.
Returns
-------
out : string
The name of the output configuration, obtained after
completing the step.
"""
data_file = os.path.join(self.exe_dir, f'{name}.data')
settings = {
'steps_subcycles': self.subcycles,
'reverse_velocities': system.vel_rev,
'write_data': data_file,
'collect_order': self.order_mode == LAMMPS_ORDER_MODE,
}
self.run_lammps(system, settings, name)
system.config = (data_file, None)
system.vel_rev = False
return os.path.basename(data_file)
@staticmethod
def _read_configuration(filename): # pragma: no cover
return
[docs] def _reverse_velocities(self, filename, outfile): # pragma: no cover
"""Reverse the velocities for a given configuration."""
return
# Atom-style -> (atom-type column, first position column) in the data
# file "Atoms" section (0-based). Image flags, when present, are the
# three columns after z.
_DATA_ATOM_COLS = {
'full': (2, 4),
'atomic': (1, 2),
'charge': (1, 3),
'molecular': (2, 3),
}
[docs] @classmethod
def _read_lammps_data_frame(cls, data_file):
"""Parse a LAMMPS data file for the Python velocity draw.
Parameters
----------
data_file : string
The LAMMPS data file to read (as written by ``dump_frame``).
Returns
-------
ids : list of int
The atom ids, sorted ascending.
positions : numpy.ndarray
The per-atom positions, shape ``(N, 3)``.
images : numpy.ndarray
The per-atom image flags, shape ``(N, 3)``.
masses : numpy.ndarray
The per-atom masses, shape ``(N, 1)``.
box : list of tuple
The ``(lo, hi)`` box bounds, one per spatial dimension.
"""
with open(data_file, encoding='utf-8') as readfile:
lines = readfile.read().splitlines()
box = []
type_mass = {}
atoms = {}
i = 0
while i < len(lines):
stripped = lines[i].strip()
tokens = stripped.split()
if tokens[-2:] in (['xlo', 'xhi'], ['ylo', 'yhi'],
['zlo', 'zhi']):
box.append((float(tokens[0]), float(tokens[1])))
elif stripped == 'Masses':
i += 2
while i < len(lines) and lines[i].strip():
parts = lines[i].split()
type_mass[int(parts[0])] = float(parts[1])
i += 1
continue
elif stripped.startswith('Atoms'):
style = (stripped.split('#')[1].strip()
if '#' in stripped else 'atomic')
if style not in cls._DATA_ATOM_COLS:
raise NotImplementedError(
'velocity_generation = "maxwell" does not support '
f'the LAMMPS atom_style "{style}" yet.'
)
tcol, xcol = cls._DATA_ATOM_COLS[style]
i += 2
while i < len(lines) and lines[i].strip():
parts = lines[i].split()
aid = int(parts[0])
pos = [float(parts[xcol]), float(parts[xcol + 1]),
float(parts[xcol + 2])]
if len(parts) >= xcol + 6:
img = [int(parts[xcol + 3]), int(parts[xcol + 4]),
int(parts[xcol + 5])]
else:
img = [0, 0, 0]
atoms[aid] = (int(parts[tcol]), pos, img)
i += 1
continue
i += 1
ids = sorted(atoms)
positions = np.array([atoms[a][1] for a in ids])
images = np.array([atoms[a][2] for a in ids])
masses = np.array([[type_mass[atoms[a][0]]] for a in ids])
return ids, positions, images, masses, box
[docs] @staticmethod
def _write_velocity_dump(dump_file, ids, pos, vel, images, box):
"""Write a LAMMPS dump frame (positions, velocities, image flags).
The columns ``id x y z vx vy vz ix iy iz`` match what the native
engine's ``read_dump`` overlay reads back.
"""
with open(dump_file, 'w', encoding='utf-8') as writefile:
writefile.write('ITEM: TIMESTEP\n0\n')
writefile.write(f'ITEM: NUMBER OF ATOMS\n{len(ids)}\n')
writefile.write('ITEM: BOX BOUNDS pp pp pp\n')
for low, high in box:
writefile.write(f'{low:.16e} {high:.16e}\n')
writefile.write('ITEM: ATOMS id x y z vx vy vz ix iy iz\n')
for k, aid in enumerate(ids):
writefile.write(
f'{aid} {pos[k, 0]:.16e} {pos[k, 1]:.16e} '
f'{pos[k, 2]:.16e} {vel[k, 0]:.16e} {vel[k, 1]:.16e} '
f'{vel[k, 2]:.16e} '
f'{int(images[k, 0])} {int(images[k, 1])} '
f'{int(images[k, 2])}\n'
)
[docs] def _draw_velocities_maxwell(self, system, vel_settings, rgen):
"""Draw shooting velocities in Python (Maxwell-Boltzmann).
Mirrors the engine path's contract -- sets ``system.config`` to a
dump frame the next propagation reads via ``read_dump`` and updates
``system.ekin`` -- but draws the velocities from the injected random
generator instead of LAMMPS's ``velocity create``.
Parameters
----------
system : object like :py:class:`.System`
The state whose velocities we redraw.
vel_settings : dict
``momentum`` (or the inf-style ``zero_momentum``) requests a
zero linear-momentum reset.
rgen : object like :py:class:`.RandomGenerator`
The velocity random generator, injected explicitly.
Returns
-------
float
The new kinetic energy (in ``(g/mol)(Angstrom/fs)^2``; transient
-- the propagation overwrites it with LAMMPS energies).
"""
# Lazy import: pyretis.engines._parts pulls in pyretis.setup at
# module load -> a circular import here (lammps_steps is imported
# while pyretis.engines is still initialising). Runtime path only.
# pylint: disable=import-outside-toplevel
from pyretis.engines._parts import kinetic_energy, reset_momentum
data_file = self.dump_frame(system)
ids, pos, images, mass, box = self._read_lammps_data_frame(data_file)
# The draw gives velocities in sqrt(kcal/g); convert to the LAMMPS
# 'real' velocity unit (Angstrom/fs) by dividing by this factor:
# uconvert((u"kcal/g")^0.5, 1u"Angstrom/fs") = 48.88821290839617
scale = 48.88821290839617
dim = int(self.settings['dimension'])
vel = np.zeros((len(ids), 3))
drawn, _ = self.draw_maxwellian_velocities(
vel[:, :dim], mass, self.beta, rgen
)
vel[:, :dim] = drawn / scale
if vel_settings.get('momentum',
vel_settings.get('zero_momentum', False)):
vel = reset_momentum(vel, mass)
dump_file = os.path.join(self.exe_dir, 'genvel.lammpstrj')
self._write_velocity_dump(dump_file, ids, pos, vel, images, box)
system.config = (dump_file, 0)
system.vel_rev = False
kin_new = kinetic_energy(vel, mass)[0]
system.ekin = kin_new
return kin_new
[docs] def modify_velocities(self, system, vel_settings, rgen):
"""Modify velocities.
Parameters
----------
system : object like :py:class:`.System`
The state whose velocities we modify.
vel_settings : dict
It contains all the info for the velocity modification.
For LAMMPS only aimless settings are supported; the
``rescale`` key is forwarded to LAMMPS.
rgen : object like :py:class:`.RandomGenerator`
The velocity random generator, injected explicitly per call.
In the default ``velocity_generation = "engine"`` mode it draws
the LAMMPS ``velocity create`` seed; in ``"maxwell"`` mode it
draws the velocities directly.
Returns
-------
dek : float
The change in the kinetic energy.
kin_new : float
The new kinetic energy.
"""
dek = None
kin_old = system.ekin
kin_new = None
if is_aimless_velocity_setting(vel_settings):
if self.velocity_generation == 'maxwell':
kin_new = self._draw_velocities_maxwell(
system, vel_settings, rgen
)
else:
logger.debug('Modifying velocities with LAMMPS.')
name = 'generate_vel'
if rgen:
seed = rgen.random_integers(1, 2147483647 - 2)
# Note: `random_integers` is inclusive for both
# numbers. We don't like to have > 2**31 - 1 as a
# possibility. For some reason, LAMMPS might also
# die/hang for a seed equal to 2**31 - 1, so this is
# why 2 is subtracted here, resulting in the range:
# [1, 2147483646]
else:
seed = 1
settings = {
'steps_subcycles': 0,
'reverse_velocities': False,
'generate_vel': {
'seed': seed,
'rescale': vel_settings.get('rescale'),
},
}
self.run_lammps(system, settings, name)
kin_new = system.ekin
if kin_old is None:
dek = float('inf')
else:
dek = kin_new - kin_old
return dek, kin_new
raise ValueError(
'LAMMPS only support the aimless velocity modification.'
)
[docs] def _calculate_pyretis_order(self, system, traj, reverse=False):
"""Calculate PyRETIS order parameters from a LAMMPS trajectory."""
order_function = self.order_function
if order_function is None:
raise ValueError(
'LAMMPS order_mode = pyretis requires an orderparameter '
'section in the PyRETIS input.'
)
frames = read_lammps_dump(traj, self.settings['dimension'])
work_system = system.copy()
orders = []
for frame in frames:
work_system.pos = frame['pos']
if reverse:
work_system.vel = -1.0 * frame['vel']
else:
work_system.vel = frame['vel']
work_system.update_box(frame['box'])
orders.append(order_function.calculate(work_system))
return np.array(orders, dtype=float)
[docs] @staticmethod
def _trim_to_common_length(order, energy):
"""Trim LAMMPS run data to the common sampled-frame length."""
if order is None:
return order, energy
lengths = [len(order)]
lengths.extend(len(values) for values in energy.values())
nframes = min(lengths)
if len(order) != nframes:
order = order[:nframes]
for key, values in energy.items():
if len(values) != nframes:
energy[key] = values[:nframes]
return order, energy
[docs] def _prepare_lammps_run(self, system, settings, name):
"""Write the LAMMPS input script and return the run command.
This is the input-preparation half of :meth:`run_lammps`,
factored out so the streaming engine
(:class:`pyretis.engines.lammps.LAMMPSEngine`) builds the
*identical* LAMMPS input -- the relaunch and streaming engines
must differ only in how the trajectory is consumed, not in what
LAMMPS is asked to do.
Parameters
----------
system : object like :py:class:`.System`
Defines the initial configuration to start from.
settings : dict
LAMMPS input settings; mutated in place with the trajectory,
order, subcycle, timestep, dimension and order-fix entries.
name : string
Base name for the LAMMPS input script and output files.
Returns
-------
cmd : list of strings
The command (with arguments) to launch LAMMPS.
"""
logger.debug('Adding input files for LAMMPS.')
self.add_input_files(self.exe_dir)
config = system.config
if config[0].endswith('.data') and config[1] in (0, None):
dest = os.path.join(
self.exe_dir,
os.path.basename(self.input_files['conf'])
)
source = config[0]
if not os.path.isabs(source):
source = os.path.join(self.exe_dir, source)
if os.path.abspath(source) != os.path.abspath(dest):
self._copyfile(source, dest)
# Name the trajectory output:
settings['traj'] = os.path.join(self.exe_dir, f'{name}.lammpstrj')
settings['order'] = f'order_{name}.txt'
# Add the subcycles & timestep settings:
settings['subcycles'] = self.subcycles
settings['timestep'] = self.timestep
settings['dimension'] = int(self.settings['dimension'])
settings['order_mode'] = self.order_mode
# Add the order-fix:
if self.order_mode == LAMMPS_ORDER_MODE:
settings['order-fix'] = self._make_order_fix(settings['order'])
# Name the input script for LAMMPS:
input_file = f'{name}.in'
script_file = os.path.join(self.exe_dir, input_file)
# Create the input file we will use or running LAMMPS:
create_lammps_md_input(
system,
self.input_files['template'],
script_file,
settings,
)
cmd = shlex.split(self.lmp)
cmd.extend([
'-in', input_file,
'-l', f'{name}.log',
'-screen', f'{name}.screen',
])
# Inject the per-run thermostat seed from the engine's random
# generator. Defining it on the command line makes ``${pyretis_seed}``
# available before the force-field template runs, so a stochastic
# thermostat (e.g. ``fix langevin ... ${pyretis_seed}``) draws a fresh,
# reproducible seed each run. Templates that do not reference it simply
# ignore the extra ``-var``.
if getattr(self, 'rgen', None) is not None:
seed = int(self.rgen.integers(1, 900000000))
cmd.extend(['-var', 'pyretis_seed', str(seed)])
return cmd
[docs] def run_lammps(self, system, settings, name):
"""Execute LAMMPS.
This method will handle input files, run LAMMPS and
return some data after the run.
Parameters
----------
system : object like :py:class:`.System`
The system defines the initial state we are using. When
``order_mode`` is ``"pyretis"`` it is also the template for
the per-frame order-parameter evaluation; the configured
order function is injected on the engine
(``self.order_function``).
settings : dict
This dict contains settings for creating the LAMMPS
input file.
name : string
This string is used as a base name for some of the
LAMMPS input scripts and output files.
Returns
-------
out[0] : numpy.array
The order parameters obtained during the run.
out[1] : dict of numpy.arrays
The energies, each dictionary key corresponds to a
specific energy term.
out[2] : string
The name of the trajectory created by LAMMPS in this run.
"""
cmd = self._prepare_lammps_run(system, settings, name)
logger.debug('Executing LAMMPS "%s".', name)
self.execute_command(cmd, cwd=self.exe_dir)
logger.debug('Reading LAMMPS order parameters & energies after run.')
energy = self.read_energies(name)
order = None
if settings.get('collect_order', True):
if self.order_mode == PYRETIS_ORDER_MODE:
order = self._calculate_pyretis_order(
system,
settings['traj'],
reverse=settings.get('reverse_order_velocities', False),
)
else:
order = self.read_order_parameters(settings['order'])
order, energy = self._trim_to_common_length(order, energy)
# Set the state of the system to the last point. ``order`` is
# None only when collect_order is False (e.g. dump_phasepoint);
# otherwise it was computed/read and an empty result is an error
# we want to surface, not silence -- so index unconditionally.
if order is not None:
system.order = order[-1]
system.ekin = energy['KinEng'][-1]
system.vpot = energy['PotEng'][-1]
system.etot = energy['TotEng'][-1]
system.temp = energy['Temp'][-1]
system.config = (settings['traj'], settings['steps_subcycles'])
return order, energy, settings['traj']
[docs] def integrate(self, system, order_function, steps, thermo='full'):
"""Propagate several integration steps with LAMMPS.
Parameters
----------
system : object like :py:class:`.System`
The system to integrate.
order_function : object like :py:class:`.OrderParameter` or None
An order function can be specified if we want to calculate
the order parameter along with the simulation.
steps : integer
The number of integration steps to perform. The LAMMPS run
length is ``max(steps - 1, 0)`` subcycles so that the first
frame is the initial configuration.
thermo : string, optional
Truthy values include a ``'thermo'`` key in each returned
step dict with the energy terms read from the LAMMPS log.
Returns
-------
out : list of dict
A list of step dicts each containing the integrated
``'system'`` (and ``'order'`` / ``'thermo'`` when requested).
"""
logger.debug('Integrating with LAMMPS.')
intervals = max(steps - 1, 0)
# Add settings for this run:
settings = {
'steps_subcycles': intervals * self.subcycles,
'reverse_velocities': system.vel_rev,
}
# Execute LAMMPS:
order, energy, traj = self.run_lammps(
system, settings, 'pyretis_md'
)
nresults = min(steps, len(order))
results = []
for i in range(nresults):
phase_point = system.copy()
phase_point.order = order[i]
phase_point.ekin = energy['KinEng'][i]
phase_point.vpot = energy['PotEng'][i]
phase_point.etot = energy['TotEng'][i]
phase_point.temp = energy['Temp'][i]
phase_point.config = (traj, i * self.subcycles)
step = {'system': phase_point}
if order_function is not None:
step['order'] = order[i]
if thermo:
step['thermo'] = rename_lammps_energy_keys(
{key: val[i] for key, val in energy.items()}
)
results.append(step)
return results
[docs] def propagate(self, path, ens_set, system, reverse=False):
"""
Propagate the equations of motion with the external code.
This is a LAMMPS-specific override: unlike the generic
:py:meth:`.ExternalMDEngine.propagate`, LAMMPS reverses the
velocities through its own input file rather than via
:py:meth:`._reverse_velocities`, so the common set-up is
handled here directly.
Parameters
----------
path : object like :py:class:`.PathBase`
This is the path we use to fill in phase-space points.
We are here not returning a new path - this since we want
to delegate the creation of the path to the method
that is running `propagate`.
ens_set : dict
Ensemble settings. It contains:
* `interfaces`: list of floats
These interfaces define the stopping criterion.
system : object like :py:class:`.System`
The system object gives the initial state for the
integration. The initial state is stored and the system is
reset to the initial state when the integration is done.
The order function is injected on the engine
(``self.order_function``).
reverse : boolean, optional
If True, the system will be propagated backward in time.
Returns
-------
success : boolean
This is True if we generated an acceptable path.
status : string
A text description of the current status of the propagation.
"""
initial_state = system
interfaces = ens_set['interfaces']
logger.debug('Running propagate with: "%s"', self.description)
prefix = str(counter())
if reverse:
logger.debug('Running backward in time.')
name = prefix + '_trajB'
else:
logger.debug('Running forward in time.')
name = prefix + '_trajF'
logger.debug('Trajectory name: "%s"', name)
# Also create a message file for inspecting progress:
msg_file_name = os.path.join(self.exe_dir, f'msg-{name}.txt')
logger.debug('Writing propagation progress to: %s', msg_file_name)
msg_file = FileIO(msg_file_name, 'w', None, backup=False)
msg_file.open()
msg_file.write(f'# Preparing propagation with {self.description}')
msg_file.write(f'# Trajectory label: {name}')
system = initial_state.copy()
logger.debug('Initial state: %s', system)
# For storing LAMMPS settings:
settings = {
'steps_subcycles': path.maxlen * self.subcycles,
'interfaces': interfaces,
'reverse_velocities': reverse != system.vel_rev,
'reverse_order_velocities': reverse,
}
if settings['reverse_velocities']:
logger.debug('Reversing velocities in initial config.')
system.vel_rev = reverse
# Propagate from the current point:
msg_file.write(f'# Interfaces: {interfaces}')
msg_file.write(f'# Running LAMMPS in folder: {self.exe_dir}')
order, energy, traj = self.run_lammps(
system, settings, name
)
msg_file.write('# Done running LAMMPS, adding points to path.')
for i, orderi in enumerate(order):
phase_point = system.copy()
phase_point.order = orderi
phase_point.ekin = energy['KinEng'][i]
phase_point.vpot = energy['PotEng'][i]
phase_point.etot = energy['TotEng'][i]
phase_point.temp = energy['Temp'][i]
phase_point.config = (traj, i * self.subcycles)
phase_point.vel_rev = reverse
status, success, stop, _ = self.add_to_path(
path,
phase_point,
interfaces[0],
interfaces[-1]
)
if stop:
break
msg_file.write('# Propagation with LAMMPS all done.')
msg_file.close()
# Accumulate the MD steps done this segment for the
# infinite-swapping scheduler's subcycle accounting.
self.steps += path.length
return success, status
[docs] def dump_phasepoint(self, phasepoint, deffnm='conf'):
"""Dump a phase point to a new file.
Note
----
We do not reverse the velocities here.
"""
settings = {
'steps_subcycles': 0,
'reverse_velocities': False,
'collect_order': False,
}
self.run_lammps(phasepoint, settings, deffnm)
[docs] def calculate_order(self, system, xyz=None, vel=None, box=None):
"""Return the last seen order parameter.
Parameters
----------
system : object like :py:class:`.System`
The state whose last seen order parameter we return.
xyz : numpy.array, optional
Unused; present to match the base signature.
vel : numpy.array, optional
Unused; present to match the base signature.
box : numpy.array, optional
Unused; present to match the base signature.
Returns
-------
out : list of floats
The last seen order parameter for the system.
"""
# xyz/vel/box unused: LAMMPS records the order during the run and
# this method only returns the already-computed last value.
_ = (xyz, vel, box)
return system.order