"""Facade for LAMMPS MD engine
This is a facade to ``PyLammps`` (added in 30th-Jul-2016 version), a convenience
wrapper for the LAMMPS Python interface i.e. where Python is extended with
LAMMPS.
Defining all interaction types requires that LAMMPS was built with the MOLECULE
package.
Notes
-----
When variables are either passed to or from ``PyLammps``, the ctypes conversion
can mean that they are unnecessarily cast, particularly from `float` to `int`.
This can cause issues as LAMMPS requires certain variables, e.g. number of
steps, to be `int`. Therefore it is always a good idea to be cast these
variables when they are read from ``PyLammps`` e.g.
``int(lmp.variables['steps'].value)``.
A minor bug in LAMMPS (Dec 2018 version) means that ``nangletypes`` returned
by ``PyLammps`` is incorrectly set to ``ndihedraltypes``. This was corrected in
Mar 2020 release.
"""
from __future__ import annotations
import logging
import os
import warnings
from collections import defaultdict, namedtuple
from contextlib import suppress
from copy import copy
from itertools import chain, combinations, count, product
from random import randint
from tempfile import NamedTemporaryFile
from typing import Any, Union
import numpy as np
try:
from lammps import PyLammps
except ModuleNotFoundError as err:
raise ModuleNotFoundError('The Python interface for LAMMPS (lammps.py) is'
' not in the PYTHONPATH. See LAMMPS documentation'
' on Python to rectify this.',
) from err
from MDMC.common import units
from MDMC.common.decorators import repr_decorator, unit_decorator, unit_decorator_getter
from MDMC.common.units import Unit
from MDMC.MD.engine_facades.facade import MDEngine, MDEngineError
from MDMC.MD.interactions import (
Bond,
BondAngle,
BondedInteraction,
Interaction,
NonBondedInteraction,
)
from MDMC.MD.simulation import ConstraintAlgorithm, KSpaceSolver, Universe
from MDMC.MD.structures import Atom
from MDMC.trajectory_analysis.compact_trajectory import CompactTrajectory
from MDMC.utilities.partitioning import partition, partition_interactions
LOGGER = logging.getLogger(__name__)
# pylint: disable=too-many-lines,possibly-used-before-assignment
[docs]
class PyLammpsAttribute:
"""
A class which has a ``PyLammps`` object as an attribute
It possesses attributes and methods relating to the ``PyLammps`` object
Parameters
----------
lmp : PyLammps, optional
Set the ``lmp`` attribute to a ``PyLammps`` object. Default is `None`,
which results in a new ``PyLammps`` object being initialised.
atom_style : str, optional
The LAMMPS ``atom_style``, which determines the properties that can be
associated which the atoms (e.g. ``charge``, ``bonds``). Default is
``full``.
Attributes
-----------
lmp : PyLammps
The ``PyLammps`` object owned by this class
"""
def __init__(self, lmp: PyLammps = None, atom_style: str = 'full'):
if lmp:
self.lmp = lmp
else:
self.lmp = PyLammps()
self.lmp.units('real')
self.lmp.atom_style(atom_style)
LOGGER.debug('%s: {lmp: %s, communicator: %s, atom_style: %s}. PyLammps'
' instance %s.',
self.__class__,
self.lmp,
None,
atom_style,
'added to class' if lmp else 'created by class')
@property
def system_state(self) -> namedtuple:
"""
Get the ``PyLammps`` wrapper system ``state`` `dict`
Returns
-------
``System``
Contains the properties of the simulation box.
"""
# Conversion from System class (which is a namedtuple) to
# ordered dict required as System cannot be pickled
system_state = self.lmp.system._asdict()
# Cast back to namedtuple to remain consist with LAMMPS system attribute
return namedtuple('System', system_state.keys())(*system_state.values())
@property
def fixes(self) -> "list[dict]":
"""
Get the ``PyLammps`` wrapper `list` of ``fixes``
Returns
-------
`list` of `dict`
Each `dict` states the ``group``, ``name`` and ``style`` of a LAMMPS
``fix` which is applied
"""
fixes = self.lmp.fixes
return fixes
@property
def fix_styles(self) -> "list[str]":
"""
Get the styles of the ``fixes`` applied in LAMMPS
Returns
-------
`list` of `str`
The styles of the ``fixes``
"""
return [fix['style'] for fix in self.fixes]
@property
def fix_names(self) -> "list[str]":
"""
Get the names of the ``fixes`` applied in LAMMPS
Returns
-------
`list` of `str`
The names of the ``fixes``
"""
return [fix['name'] for fix in self.fixes]
@property
def dumps(self) -> list:
"""
Get the PyLammps wrapper list of dumps
Dumps are LAMMPS commands which write atom quantities to file for
specified timesteps
"""
dumps = self.lmp.dumps
return dumps
[docs]
@repr_decorator('lmp', 'lmp_universe', 'lmp_simulation')
class LAMMPSEngine(PyLammpsAttribute, MDEngine):
"""
Facade for LAMMPS
One notable issue with creating the LAMMPS facade is that some
``pair_styles`` are combinations of both dispersion and coulombic
interactions. For example, although LAMMPS has ``lj/cut``, ``coul/cut``, and
``coul/long``, there is no ``lj/long`` pair_style; it only exists in
combination with ``coul/long`` (i.e. ``lj/long/coul/long``). This means that
the facade has to not just parse each nonbonded interaction, but also has to
determine if that ``pair_style`` should be combined with another
``pair_style`` (e.g. ``lj/long`` and ``coul/long`` have to be combined
before passing to ``pair_style`` and ``pair_coeff``). An additional
complication of this is that when setting the coefficients (``Parameters``)
of these interactions, they again have to be set together. So in theory a
coulombic interaction needs to know the dispersion interaction it has been
paired with an also pass those coefficients when calling ``pair_coeff``. In
practice a simplification can be made, at least in the case of currently
implemented pair styles: As coulombic ``pair_styles`` do no take any
coefficients, any coulombic ``pair_styles`` that comprise a combined
``pair_style`` can be ignored for the purposes of setting the coulombic
``pair_style``; this can be taken care of when ``pair_coeff`` is called from
the correspoding dispersion interaction, as this possesses the dispersion
coefficients that also need to be passed when the coefficients of this pair
style are set.
"""
def __init__(self):
super().__init__()
self.universe = None
self.lmp_universe = None
self._saved_config = None
self.trajectory_file = None
self.lmp_simulation = None
@property
def saved_config(self) -> np.ndarray:
"""
Get the saved configuration of the atomic positions
Returns
-------
numpy.ndarray
The configuration from the start of the run. Each row of the
``array`` corresponds to the LAMMPS atom ``ID - 1`` (offset is due
to array zero indexing) and the columns of the ``array`` are the
``x``, ``y``, ``z`` components of the ``position``, the ``mass`` and
the ``charge`` of each ``Atom``.
"""
return self._saved_config
@property
def temperature(self) -> float:
"""
Get or set the temperature of the simulation in ``K``
Returns
-------
`float`
Temperature in ``K``
"""
return self.lmp_simulation.temperature
@temperature.setter
@unit_decorator(unit=units.TEMPERATURE)
def temperature(self, value: float) -> None:
self.lmp_simulation.temperature = value
@property
def pressure(self) -> float:
"""
Get or set the pressure of the simulation in ``atm``
Returns
-------
`float`
Pressure in ``atm``
"""
return self.lmp_simulation.pressure
@pressure.setter
@unit_decorator(unit=units.PRESSURE)
def pressure(self, value: float) -> None:
self.lmp_simulation.pressure = value
@property
def ensemble(self) -> 'LAMMPSEnsemble':
"""
Get or set the ensemble object which applies a ``thermostat`` and/or
``barostat`` to LAMMPS
Returns
-------
LAMMPSEnsemble
The simulation thermodynamic ensemble
"""
return self.lmp_simulation.ensemble
@ensemble.setter
def ensemble(self, value: 'LAMMPSEnsemble') -> None:
self.lmp_simulation.ensemble = value
@property
def thermostat(self) -> str:
"""
Get or set the `str` which specifies the thermostat
Returns
-------
`str`
The ``thermostat`` name
"""
return self.ensemble.thermostat
@thermostat.setter
def thermostat(self, value: str) -> None:
self.ensemble.thermostat = value
@property
def barostat(self) -> str:
"""
Get or set the `str` which specifies the barostat
Returns
-------
`str`
The ``barostat`` name
"""
return self.ensemble.barostat
@barostat.setter
def barostat(self, value: str) -> None:
self.ensemble.barostat = value
[docs]
def setup_universe(self, universe: Universe, **settings: dict) -> None:
"""
Creates the simulation box, the atomic configuration, and the topology
in LAMMPS
Parameters
----------
universe : Universe
The MDMC ``Universe`` which will be setup in LAMMPS.
**settings
``atom_style`` (`str`)
A LAMMPS ``atom_style`` `str`. The default setting of ``real``
will generally be appropriate.
"""
super().__init__(atom_style=settings.get('atom_style', 'full'))
self.universe = universe
self.lmp_universe = LAMMPSUniverse(self.universe, self.lmp, **settings)
self._saved_config = None
[docs]
def setup_simulation(self, **settings: dict) -> None:
"""
Sets simulation parameters in LAMMPS, such as the thermodynamic
variables, thermostat/barostat parameters and trajectory settings
Parameters
----------
**settings
Passed to ``LAMMPSSimulation``
"""
self.lmp_simulation = LAMMPSSimulation(universe=self.universe,
time_step=self.time_step,
traj_step=self.traj_step,
lmp=self.lmp,
**settings)
[docs]
def minimize(self, n_steps: int, minimize_every: int = 10,
output_log: str = None,
work_dir: str = None, **settings: dict) -> None:
"""Moves the atoms towards a potential energy minimum,
by performing an MD simulation interrupted periodically
by a structure relaxation. In the end, the configuration
with the lowest potential energy reached during the run
is kept.
Parameters
----------
n_steps : int
Number of MD simulation steps
minimize_every : int, optional, default 10
The structure relaxation will be performed every
'minimize_every' steps of the MD simulation
output_log : str, optional
Not used in this facade
work_dir : str, optional
Not used in this facade
**settings
etol: float, energy tolerance criteria for energy minimisation
ftol: float, force tolerance criteria for force minimisation, active only if non-zero
maxiter: int, maximum number of steps in a single structure relaxation
maxeval: int, maximum number of force calculations in a single structure relaxation
"""
# Check fix styles for shake or rattle styles and remove them
if 'constrain' in self.fix_names:
LOGGER.debug('%s Remove %s constraint from fixes: %s',
self.__class__,
self.universe.constraint_algorithm,
self.fixes)
self.lmp.unfix('constrain')
etol = settings.get('etol', 1.e-4)
ftol = settings.get('ftol', 0.)
maxiter = settings.get('maxiter', 10000)
maxeval = settings.get('maxeval', 10000)
LOGGER.info('%s minimize: {n_steps: %s, minimize_every: %s, etol: %s, ftol: %s,'
' maxiter: %s, maxeval: %s}',
self.__class__,
n_steps,
minimize_every,
etol,
ftol,
maxiter,
maxeval)
self.lmp.minimize(etol,
ftol,
maxiter, # this is the number of relaxation steps
maxeval) # this is the number of force evaluations
best_energy = self.lmp.eval("pe")
self.save_config()
for _ in range(int(n_steps/minimize_every)):
if self.lmp_simulation.temperature is not None:
self.lmp.velocity(
'all', 'scale', convert_unit(self.lmp_simulation.temperature))
self.lmp.run(minimize_every)
self.lmp.minimize(etol,
ftol,
maxiter, # this is the number of relaxation steps
maxeval) # this is the number of force evaluations
energy = self.lmp.eval("pe")
if energy < best_energy:
self.save_config()
best_energy = energy
self.reset_config()
# Reapply the constraints
if self.universe.constraint_algorithm:
self.ensemble.remove_ensemble_fixes()
self.lmp_universe.apply_constraints()
self.ensemble.apply_ensemble_fixes()
[docs]
def run(self, n_steps: int, equilibration=False, output_log: str = None,
work_dir: str = None, **settings: dict):
if not equilibration:
# Remove previous dumps if they exist
if 'traj1' in [dump['name'] for dump in self.dumps]:
self.lmp.undump('traj1')
# Store the trajectory in a NamedTemporaryFile
# pylint: disable=consider-using-with
self.trajectory_file = NamedTemporaryFile()
# f_name = self.trajectory_file.name
# Custom trajectory output just saves the atom ID, type and
# positions
LOGGER.debug('%s set trajectory dump output to %s',
self.__class__,
self.trajectory_file)
self.lmp.dump('traj1', 'all', 'custom', self.traj_step,
self.trajectory_file.name, 'id', 'type', 'x', 'y',
'z')
else:
reset_to_nve = False
if self.thermostat is self.barostat is None:
# If NVE ensemble, add a berendsen thermostat for equilibration
reset_to_nve = True
self.thermostat = 'berendsen'
LOGGER.info('%s run: {n_steps: %s, equilibration: %s}',
self.__class__,
n_steps,
equilibration)
try:
self.lmp.run(n_steps)
except Exception as exc:
raise MDEngineError("There has been an error running the LAMMPS simulation.") from exc
if equilibration and reset_to_nve:
self.thermostat = None
[docs]
def convert_trajectory(self, start: int = 0, stop: int = None,
step: int = 1, **settings: dict) -> CompactTrajectory:
"""
Converts between a LAMMPS trajectory dump and an MDMC
``CompactTrajectory``
The LAMMPS dump must include at least ``id``, ``atom_type``, and ``xyz``
``positions``. The ``xyz`` ``positions`` must be consecutive and in that
order. The same is true of the ``xyz`` components of the ``velocity``, if
they are provided.
Parameters
----------
start : int
The index of the first trajectory, inclusive.
stop : int
The index of the last trajectory, exclusive.
step : int
The step size between trajectories.
**settings
``scaled_positions`` (`bool`)
If the ``trajectory_file`` has scaled ``positions``
``atom_IDs`` (`list`)
LAMMPS ``ID`` of the atoms which should be included. If not passed
then all atoms are included in the converted trajectory.
Returns
-------
``CompactTrajectory``
The MDMC ``CompactTrajectory`` corresponding to the LAMMPS
``trajectory_file``.
Raises
------
AssertionError
If ``universe`` is passed, and the number of atoms in the
``trajectory_file`` is not the same as in the ``universe``.
TypeError
If ``trajectory_file`` describes a triclinic universe.
"""
# Change expected position string if scaled positions are used
pos_string = 'xs' if settings.get('scaled_positions', False) else 'x'
traj_dimensions = np.zeros(3)
frame_n = start
# Use count to create range so that stop can be undefined
frame_indexes = count(start, step)
# next_frame_n next attribute is assigned dynamically
next_frame_n = next(frame_indexes) # pylint: disable=no-member
traj = CompactTrajectory() #the instance of our new trajectory object
def _make_gen(reader):
"""A support function for splitting a binary file into buffers
Args:
reader (file descriptor): an open file or a file-like object
Yields:
bytes: a byte string read from the file
"""
b = reader(1024 * 1024)
while b:
yield b
b = reader(1024*1024)
# here we check how long the trajectory really is
with open(self.trajectory_file.name, 'rb') as file_handler:
file_generator = _make_gen(file_handler.raw.read)
line_count = sum( buf.count(b'\n') for buf in file_generator )
if line_count == 0:
raise IOError("Trajectory file empty, run may have failed"
" or run length may be less than dump frequency.")
# And header_size will tell us how many lines per frame
# are added on top of the atom positions
header_size = 0
with open(self.trajectory_file.name, 'r', encoding='UTF-8') as file_handler:
line = file_handler.readline()
while line:
# LAMMPS TIMESTEP is the number of time steps that have elapsed. To
# avoid confusion with time_step (the amount of time that elapses in
# a single simulation step, i.e. dt), these are referred to as
# frames.
if 'ITEM: TIMESTEP' in line:
line = file_handler.readline()
frame = int(line.split()[0])
header_size += 2
if 'ITEM: NUMBER OF ATOMS' in line:
line = file_handler.readline()
n_atoms = int(line.split()[0])
header_size += 2
# Check that n_atoms is as expected, if a universe was passed
if self.universe:
assert n_atoms == len(self.universe.atoms)
if 'ITEM: BOX BOUNDS' in line:
header_size += 1
traj_dimensions *= 0.0
temp_dim = []
# CURRENTLY ASSUMES ORTHOGONAL SIMULATION BOX
if 'xy' in line:
raise TypeError('triclinic simulation boxes have not'
' been implemented')
# Test dimensions are as expected, if a universe was passed
# and we are not using an NPT or NPH ensemble
for i in range(3):
line = file_handler.readline()
header_size += 1
dmin, dmax = [float(splt) for splt in line.split()]
temp_dim.append((dmin, dmax))
traj_dimensions[i] = dmax-dmin
if self.universe and not ('npt' in self.fix_names or 'nph' in self.fix_names):
for i in range(3):
dmin, dmax = temp_dim[i]
assert dmin == 0.0
# unit is taken from universe dimensions (which is a
# UnitArray)
assert dmax == convert_unit(self.universe.dimensions[i],
self.universe.dimensions.unit)
if 'ITEM: ATOMS' in line:
header_size += 1
if frame_n == start:
# LAMMPS dump files contain order of LAMMPS atom properties,
# at each time step. As these should not change with time
# step only determine this order for first required time
# step. Assumes that position components (x y z) and
# velocity components (vx vy vz) are always adjacent and
# ordered as shown.
splt = line.split()
# Requires id, type and position to be defined, velocity is
# optional
i_id, i_type, i_pos = [splt.index(prop) - 2 for prop
in ['id', 'type', pos_string]]
i_vel = splt.index('vx') if 'vx' in splt else None
# now we try to get the correct number of frames in the trajectory
real_n_steps = 1 + line_count // (n_atoms + header_size)
traj.preAllocate(n_steps = real_n_steps,
n_atoms = n_atoms,
useVelocity = i_vel is not None)
traj.setDimensions(traj_dimensions)
if frame_n == next_frame_n:
# Reads all atom lines before creating any atoms. By
# creating a list of tuples of (LAMMPS_ID, atom), this
# allows the lines to be reordered based on LAMMPS_ID. This
# is required as by default LAMMPS does not sort by ID, so
# the same atom will not appear in the same place for each
# time step.
header_size = 0
lines = []
for _ in range(n_atoms):
split_line = file_handler.readline().split()
# convert id to int
split_line[i_id] = int(split_line[i_id])
lines.append(split_line)
# sort list of lists based on id
lines = sorted(lines, key=lambda x: x[i_id])
sorted_lines = np.array(lines, dtype = traj.dtype)
atom_types = sorted_lines[:,i_type].astype(np.int64)
if not traj.validateTypes(atom_types):
raise TypeError("CompactTrajectory received wrong atom type array")
if i_vel is not None:
traj.writeOneStep(step_num = frame_n,
time = frame * self.time_step,
positions = sorted_lines[:,i_pos:i_pos+3],
velocities = sorted_lines[:,i_vel:i_vel+3])
else:
traj.writeOneStep(step_num = frame_n,
time = frame * self.time_step,
positions = sorted_lines[:,i_pos:i_pos+3])
traj.setDimensions(traj_dimensions, step_num = frame_n)
# next_frame_n next attribute is assigned dynamically
# pylint: disable=no-member
next_frame_n = next(frame_indexes)
frame_n += 1
if stop is not None and frame_n >= stop:
break
line = file_handler.readline()
atom_symbols = {}
atom_masses = {}
for at_id in np.unique(traj.atom_types):
symbol, mass = self.lmp_universe.atom_type_properties[at_id-1]
atom_symbols[at_id] = symbol
atom_masses[at_id] = mass
traj.labelAtoms(atom_symbols, atom_masses)
# electric charges, for completeness
# (otherwise we cannot output an Atom from CompactTrajectory)
charge_list = []
for struc in self.parent_simulation.universe.structure_list:
if isinstance(struc, Atom):
charge_list.append([struc.ID, struc.charge])
charges = np.array(charge_list)
sequence = np.argsort(charges[:, 0])
charges = charges[sequence][:, 1]
traj.setCharge(charges)
# we conclude the creation of the trajectory
traj.postProcess()
return traj
[docs]
def update_parameters(self) -> None:
self.lmp_universe.update_parameters()
[docs]
def clear(self) -> None:
"""
Deletes all atoms of the MD engine, restores all settings to their default values,
and frees all memory in LAMMPS.
"""
self.lmp.clear()
[docs]
def save_config(self) -> None:
# It is not possible to deepcopy the LAMMPS wrapper atoms attribute,
# or the individual atoms, so instead this saves the x, y, z, mass
# and charge in a NumPy array with the indexes given by the atom ID
# (with a -1 offset due to zero index)
# The atoms attribute also is not iterable
n_atoms = self.system_state.natoms
LOGGER.info('%s save_config: {n_atoms: %s}. Config saved.',
self.__class__,
n_atoms)
atoms = np.zeros([n_atoms, 5])
tmp_mass = {}
for type_ID, atom_type_group in self.lmp_universe.atom_types.items():
tmp_mass[type_ID] = float(atom_type_group[0].mass)
for i in range(n_atoms):
atom = self.lmp.atoms[i]
atom_type = atom.type
# _, mass = self.lmp_universe.atom_type_properties[atom_type-1]
atoms[atom.id-1, :] = (list(atom.position) + [tmp_mass[atom_type], atom.charge])
saved_config = atoms
self._saved_config = saved_config
[docs]
def reset_config(self) -> None:
self.lmp_universe.set_config(self.saved_config)
[docs]
def eval(self, variable: str) -> Any:
return self.lmp.eval(variable)
[docs]
@repr_decorator('universe')
class LAMMPSUniverse(PyLammpsAttribute):
# Class has to maintain a lot of state (attributes) as PyLammps class does
# not
# pylint: disable=too-many-instance-attributes
"""
A class with what would be the equivalent in LAMMPS to the universe (i.e.
the configuration and topology)
Parameters
----------
universe : Universe
The MDMC ``Universe`` used to create the ``LAMMPSUniverse``
lmp : PyLammps, optional
Set the ``lmp`` attribute to a ``PyLammps`` object. Default is `None`,
which results in a new ``PyLammps`` object being initialised.
**settings
``atom_style`` (`str`)
A LAMMPS ``atom_style`` string. The default setting of ``real`` will
generally be appropriate.
``nonbonded_mix`` (`str`)
The name of the formula which determines non-bonded mixing
Attributes
----------
universe : Universe
The MDMC ``Universe`` which has been converted to this
``LAMMPSUniverse``.
atom_dict : dict
A `dict` with {``MDMC_atom``: ``LAMMPS_atom_id``}, where ``MDMC_atom``
is an MDMC ``Atom`` and ``LAMMPS_atom`` is the corresponding LAMMPS
``Atom``.
atom_types : dict
A `dict` with {``type_ID``: ``MDMC_atom_group``}, where the ``type_ID``
is a unique `int` and ``MDMC_atom_group`` is a `list` of ``Atom`` which
are identical in terms of element and interactions.
atom_type_properties : list of tuples
Each `tuple` is (``symbol``, ``mass``) for all ``atom_types`` (ordered)
by ``atom_type``, where ``symbol`` is a `str` specifying the element of
the ``atom_type`` and ``mass`` is a `float` specifying the ``mass`` of
the ``atom_type``.
bonds : list of Bonds
All ``Bond`` interactions in the MDMC ``Universe``.
angles : list of BondAngles
All ``BondAngle`` interactions in the MDMC ``Universe``.
couls : list of Coulombics
All ``Coulombic`` interactions in the MDMC ``Universe``.
disps : list of Dispersions
All ``Dispersion`` interactions in the MDMC ``Universe``.
bond_ID : dict
A `dict` of {``bond``: ``ID pairs``} relating each ``Bond`` to a LAMMPS
``ID``.
angle_ID : dict
A `dict` of {``angle``: ``ID pairs``} relating each ``BondAngle`` to a
LAMMPS ``ID``.
proper_ID : dict
A `dict` of {``proper``: ``ID pairs``} relating each ``DihedralAngle``
(with ``improper == False``) to a LAMMPS ``ID``.
improper_ID : dict
A `dict` of {``improper``: ``ID pairs``} relating each ``DihedralAngle``
(with ``improper == True``) to a LAMMPS ``ID``.
"""
def __init__(self, universe: Universe, lmp: PyLammps = None, **settings: dict):
super().__init__(lmp, settings.get('atom_style', 'full'))
self.universe = universe
self.atom_dict = {}
self.atom_types = {}
self.atom_type_properties = []
self.bonds = []
self.angles = []
self.dihedrals = []
self.disps = []
self.couls = []
self.propers = []
self.impropers = []
# ID is an acronym
# pylint: disable=invalid-name
self.bond_ID = {}
self.angle_ID = {}
self.proper_ID = {}
self.improper_ID = {}
self.nonbonded_mix = None
if "OMP_NUM_THREADS" in os.environ:
omp_num_threads_str = os.environ["OMP_NUM_THREADS"]
try:
omp_num_threads = int(omp_num_threads_str)
except ValueError:
pass
else:
if omp_num_threads > 1:
self.lmp.command("package omp 0") # use OMP_NUM_THREADS
self.lmp.command("suffix omp") # add /omp to relevant styles
self._define_simulation_box(self.universe)
self._build_config(self.universe)
self._add_topology(self.universe, **settings)
self.update_parameters()
@property
def nonbonded_mix(self) -> str:
"""
Get or set the formula used to calculate nonbonded interactions between
different ``atom_types``
Options are ``geometric``, ``arithmetic`` and ``sixthpower``, which are
defined in the LAMMPS documentation.
Returns
-------
`str`
The name of the formula which determines non-bonded mixing
Raises
------
ValueError
`str` specifies an unsupported mix name
"""
return self._nonbonded_mix
@nonbonded_mix.setter
def nonbonded_mix(self, value: str) -> None:
if not value:
self._nonbonded_mix = value
else:
supported = ['geometric', 'arithmetic', 'sixthpower']
if value.lower() not in supported:
raise ValueError('The supported mixes are: {0}, {1}, {2}'
''.format(*supported))
self._nonbonded_mix = value.lower()
[docs]
def update_parameters(self) -> None:
"""
Updates the LAMMPS force field parameters from the MDMC universe
"""
self._update_charges()
self._update_bonded_interactions('bond', self.bonds)
self._update_bonded_interactions('angle', self.angles)
self._update_bonded_interactions('dihedral', self.propers)
self._update_bonded_interactions('improper', self.impropers)
self._update_dispersions(self.universe)
def _define_simulation_box(self, universe: Universe) -> None:
"""
Defines a region and creates a simulation box that fills this region
Parameters
----------
universe : Universe
The MDMC ``Universe`` used to create the region and simulation box.
"""
# ID is an acronym
# pylint: disable=invalid-name
region_ID = 'universe'
self._create_lammps_region(universe, region_ID)
n_atom_types = len(universe.atom_types)
# Determine number of bond and angle types
bonded_interaction_types = [i.name for i in set(universe.interactions)
if issubclass(type(i), BondedInteraction)]
# Dihedrals are only separated into proper and improper at the attribute
# level
dihedrals = partition_interactions(set(universe.interactions),
['DihedralAngle'])[0]
dihedral_types = [dihedral.improper for dihedral in list(dihedrals)]
n_bond_types = bonded_interaction_types.count('Bond')
n_angle_types = bonded_interaction_types.count('BondAngle')
n_dihedral_types = dihedral_types.count(False)
n_improper_types = dihedral_types.count(True)
# Determine max number of bonds and angles per atom
atoms = universe.atoms
max_bonds_per_atom = self._max_n_interaction(atoms, 'Bond')
max_angles_per_atom = self._max_n_interaction(atoms, 'BondAngle')
max_dihedrals_per_atom = self._max_n_interaction(atoms, 'proper')
max_impropers_per_atom = self._max_n_interaction(atoms, 'improper')
LOGGER.debug('%s create LAMMPS box: {n_atom_types: %s'
' n_bond_types: %s, n_angle_types: %s,'
' n_dihedral_types: %s, n_improper_types: %s,'
' max_bonds_per_atom: %s, max_angles_per_atom: %s,'
' max_dihedrals_per_atom: %s, max_impropers_per_atom: %s}',
self.__class__,
n_atom_types,
n_bond_types,
n_angle_types,
n_dihedral_types,
n_improper_types,
max_bonds_per_atom,
max_angles_per_atom,
max_dihedrals_per_atom,
max_impropers_per_atom)
self.lmp.create_box(n_atom_types,
region_ID,
'bond/types', n_bond_types,
'angle/types', n_angle_types,
'dihedral/types', n_dihedral_types,
'improper/types', n_improper_types,
'extra/bond/per/atom', max_bonds_per_atom,
'extra/angle/per/atom', max_angles_per_atom,
'extra/dihedral/per/atom', max_dihedrals_per_atom,
'extra/improper/per/atom', max_impropers_per_atom,
)
# ID is an acronym
# pylint: disable=invalid-name
def _create_lammps_region(self, universe: Universe, region_ID: str) -> None:
"""
Create a geometry of the simulation box in LAMMPS
Parameters
----------
universe : Universe
The MDMC ``Universe`` used to create the region.
region_ID : str
The LAMMPS region ``ID``
"""
xlo = ylo = zlo = 0.
xhi, yhi, zhi = universe.dimensions
region_ID = 'universe'
# 'block' gives a cuboidal universe
self.lmp.region(region_ID, 'block', xlo, xhi, ylo, yhi, zlo, zhi,
units='box')
def _build_config(self, universe: Universe) -> None:
"""
Adds atoms to LAMMPS
Parameters
----------
universe : Universe
The MDMC ``Universe`` used to fill the LAMMPS box with atoms.
"""
self.atom_types = universe.atom_types
# Assume all atoms of the same type have the same element and mass
# Sort atoms based on atom_type (i.e. numerically starting at 1) so
# that it is possible to index into atom_type_properties with
# atom_type - 1 (where the -1 is to account for 0 indexing on lists)
self.atom_type_properties = [(atom[0].element, atom[0].mass)
for _, atom
in sorted(self.atom_types.items())]
LOGGER.info('%s Add atoms to LAMMPS',
self.__class__)
for type_ID, atom_type_group in self.atom_types.items():
# The mass below has associated units, which causes a segfault when
# it is passed to LAMMPS - cast to float to remove units
self.lmp.mass(type_ID, float(atom_type_group[0].mass))
for atom in atom_type_group:
self.lmp.create_atoms(type_ID, 'single', *atom.position)
# As PyLammps has a bug preventing getting atom id from
# self.lmp.atoms[index].id, use number of atoms as proxy for
# new atom id (as it is sequential)
lmp_atom_id = self.lmp.atoms.natoms
self.lmp.set('atom', lmp_atom_id,
'vx', atom.velocity[0],
'vy', atom.velocity[1],
'vz', atom.velocity[2])
self.atom_dict[atom] = lmp_atom_id
[docs]
def set_config(self, config: np.ndarray) -> None:
"""
Changes the positions of all of the atoms in the LAMMPS wrapper
Parameters
----------
config : numpy.ndarray
The ``positions``, ``mass`` and ``charge`` of the ``Atom`` objects,
used to set the LAMMPS configuration. Each row of the array must
correspond to the LAMMPS ``atom ID - 1`` (offset is due to ``array``
zero indexing) and the columns of the ``array`` must be the ``x``,
``y``, ``z`` components of the ``position``, the ``mass`` and the
``charge`` of each ``Atom``.
Raises
------
IndexError
If ``config`` does not contain the same number of atoms as LAMMPS
possesses.
"""
# Raise an IndexError if the config is not the correct size
n_atoms = self.system_state.natoms
if len(config) != n_atoms:
raise IndexError('the new configuration does not specify the'
' correct number of atoms')
# The LAMMPS wrapper does not allow the configuration to be updated
# simply by setting all atoms. Instead the position of the atoms must be
# reset.
index_components = list(enumerate(['x', 'y', 'z']))
# LAMMPS IDs start at 1, so are offset from config indexes
for id_offset in range(n_atoms):
for index, component in index_components:
self.lmp.set('atom', id_offset+1, component,
config[id_offset][index])
@staticmethod
def _max_n_interaction(atoms: list[Atom], name: str) -> int:
"""
Parameters
----------
atoms : list of Atoms
name : str
An ``Interaction`` type, for example ``Bond``.
Returns
-------
`int`
The maximum number of interactions with a given name that any
``Atom`` in ``atoms`` possesses
"""
max_inters = 0
for atom in atoms:
# Filter interactions by name
if name in ['proper', 'improper']:
inters = filter(lambda i: i.name == 'DihedralAngle' and
i.improper == bool(name == 'improper'), atom.interactions)
else:
inters = filter(lambda i: i.name == name, atom.interactions)
n_inters = len(list(inters))
max_inters = n_inters if n_inters > max_inters else max_inters
return max_inters
def _add_topology(self, universe: Universe, **settings: dict) -> None:
"""
Add the bonded and nonbonded interactions to LAMMPS
Parameters
----------
universe : Universe
The MDMC ``Universe`` used to define the topology.
**settings
``nonbonded_mix`` (`str`)
The name of the formula which determines non-bonded mixing
Raises
------
NotImplementedError
If ``universe`` contains an interaction type that has not been
implemented in the LAMMPS facade
"""
LOGGER.info('%s Add topology to LAMMPS',
self.__class__)
# *_ is for pylint as it does not know about the output of partition_interactions
bonds, angles, dihedrals, disps, couls, others, *_ = partition_interactions(
set(universe.interactions),
['Bond', 'BondAngle', 'DihedralAngle', 'Dispersion', 'Coulombic'],
unpartitioned=True,
lst=True)
if others:
raise NotImplementedError('This interaction type has not been'
' implemented in the LAMMPS facade')
self.bonds = bonds
self.angles = angles
self.dihedrals = dihedrals
self.disps = disps
self.couls = couls
pair_styles, pair_mods, pair_coeff_cmds = self._pair_commands(universe)
if pair_styles:
LOGGER.debug('%s Add nonbonded pair styles to LAMMPS: {pair_styles:'
' %s, pair_mods: %s, pair_coeff_cmds: %s}',
self.__class__,
pair_styles,
pair_mods,
pair_coeff_cmds)
# hybrid/overlay allows multiple pair_styles for same atom_type pair
self.lmp.pair_style('hybrid/overlay', *pair_styles)
self._update_charges()
self.nonbonded_mix = settings.get('nonbonded_mix')
self._update_dispersions(universe, pair_coeff_cmds)
# Apply LAMMPS modifications to nonbonded interactions
self._modify_nonbonded_styles(pair_mods)
if bonds:
LOGGER.debug('%s Add bonds to LAMMPS',
self.__class__)
# Set used to remove duplicate bond styles, which are not required
# to be (and in fact cannot) be passed to LAMMPS hybrid bond_style
self.lmp.bond_style('hybrid',
*set(tuple(parse_bonded_styles(b)
for b in list(bonds))))
self._create_bonded_interactions('bond', bonds)
if angles:
LOGGER.debug('%s Add angles to LAMMPS',
self.__class__)
# Set used to remove duplicate angle styles, which are not required
# to be (and in fact cannot) be passed to LAMMPS hybrid angle_style
# pylint: disable=not-an-iterable
# false positive raised here
self.lmp.angle_style('hybrid',
*set(tuple(parse_bonded_styles(a)
for a in angles)))
self._create_bonded_interactions('angle', angles)
if dihedrals:
LOGGER.debug('%s Add dihedrals to LAMMPS',
self.__class__)
# Split dihedrals into proper and impropers
impropers, propers = partition(dihedrals, lambda d: d.improper)
self.impropers, self.propers = list(impropers), list(propers)
# Set used to remove duplicate dihedral styles, which are not
# required to be (and in fact cannot) be passed to LAMMPS hybrid
# dihedral_style or improper_style
proper_styles = set(tuple(parse_bonded_styles(p) for p
in self.propers))
improper_styles = set(tuple(parse_bonded_styles(i) for i
in self.impropers))
if proper_styles:
self.lmp.dihedral_style('hybrid', *proper_styles)
self._create_bonded_interactions('dihedral', self.propers)
if improper_styles:
self.lmp.improper_style('hybrid', *improper_styles)
self._create_bonded_interactions('improper', self.impropers)
if self.universe.constraint_algorithm:
self.apply_constraints()
@staticmethod
def _pair_commands(universe: Universe) -> tuple:
"""
Parses all the ``NonBondedInteractions`` for every appropriate
combination of ``atom_type`` pairs in an MDMC ``Universe``, returning
the correctly formatted input for ``pair_style`` and ``pair_coeff``
LAMMPS commands.
THE PARSING OF PAIR MODIFIERS IS LIMITED - if there is more than one
interaction with the same ``pair_style``, and if these pair styles have
different modifiers, all of the modifiers will be applied to this
``pair_style``. The correct implementation would apply each modifier to
a copy of the ``pair_style``.
Parameters
----------
universe : Universe
The MDMC ``Universe`` containing the ``NonBondedInteraction``
objects to be parsed.
Returns
-------
`tuple`
(``pair_styles``, ``pair_modifiers``, ``pair_coeff_cmds``), where
``pair_styles`` is a flattened `list` of `str` for the
``pair_style`` commands to be set in the LAMMPS interface,
``pair_modifiers`` is a `list` of `list` of `str` for the
``pair_modify`` commands to be set in the LAMMPS interface, and
``pair_coeff_cmds`` is a `list` of `str` for the ``pair_coeff``
commands for each ``atom_type`` pair to be set individually in the
LAMMPS interface.
"""
pair_styles = set()
pair_modifiers = defaultdict(set)
pair_coeff_cmds = []
for pair, inters in universe.nbis_by_atom_type_pairs.items():
# Generates list of tuples containing styles and cutoffs
styles_mods = parse_all_nonbonded_styles(inters)
# One list of pair modifier for
for style, mod in styles_mods.items():
# Only append if there is a modifier on the pair style
if mod:
# Only the first style index is required as the rest contain
# cutoffs etc
pair_modifiers[style[0]].update(set(mod))
# Generates list of tuples that contain each pair_coeff command
parsed_coeffs = parse_dispersion_coefficients(inters,
styles_mods.keys())
coeff_cmds = [' '.join([str(a_type) for a_type in pair]) + ' '
+ coeff for coeff in parsed_coeffs]
pair_styles.update(styles_mods.keys())
pair_coeff_cmds += coeff_cmds
# Remove duplicates in pair_styles, create flattened list
pair_styles = list(chain.from_iterable(pair_styles))
pair_modifiers = [[style, *mod] for style, mod
in pair_modifiers.items()]
return pair_styles, pair_modifiers, pair_coeff_cmds
def _update_charges(self) -> None:
"""
Updates the ``charges`` in LAMMPS
Raises
------
AttributeError
If one or more ``Atom`` do not have a ``charge`` (or
``charge is None``)
"""
for atom, lmp_atom_id in self.atom_dict.items():
try:
self.lmp.set('atom',
lmp_atom_id,
'charge',
convert_unit(atom.charge))
except ValueError as error:
raise AttributeError('LAMMPS requires all atoms in the universe'
' to have a charge.') from error
def _update_dispersions(self, universe: Universe, pair_coeff_cmds: 'list[str]' = None) -> None:
"""
Updates ``Dispersion`` interactions in LAMMPS
Parameters
----------
universe : Universe
The MDMC ``Universe`` containing ``NonBondedInteractions`` used to
generate the ``pair_coeff`` commands if ``pair_coeff_cmds`` is not
passed.
pair_coeff_cmds : list of str
A `list` of ``pair_coeff`` commands for each ``atom_type`` pair to
be set individually in the LAMMPS interface.
NOTE: the ``Coulombic`` for the appropriate ``atom_type`` pairs are
also set in this method.
"""
if not pair_coeff_cmds:
_, _, pair_coeff_cmds = self._pair_commands(universe)
for cmd in pair_coeff_cmds:
self.lmp.pair_coeff(cmd)
def _modify_nonbonded_styles(self, pair_mods: list) -> None:
"""
Applies modifications to nonbonded ``pair_styles``, such as the VdW tail
correction or setting a mixing style for interactions acting on unlike
``atom_type`` pairs
Parameters
----------
pair_mods : list
A `list` of `list` which are valid ``pair_modify`` commands
"""
# MDMC only allows nonbonding mixing on all NonbondedInteractions or
# none
if self.nonbonded_mix:
self.lmp.pair_modify('mix', self.nonbonded_mix)
for mod in pair_mods:
self.lmp.pair_modify('pair', *mod)
def _create_bonded_interactions(self, lmp_name: str,
bonded_interactions: list[BondedInteraction]) -> None:
"""
Creates coefficients and new bonded interactions in LAMMPS, and fills
the relevant ``BondedInteraction`` ID (e.g. ``self.bond_ID`` for bonds,
``self.angle_ID`` for angles) `dict` with
``bonded_interaction``: ``ID pairs``
This generic method can be used for bonds, angles, dihedrals (proper),
and impropers. It can only be used for a single type of
``bonded_interactions`` (e.g. only bonds)
Parameters
----------
lmp_name : str
The name of the bonded interaction type used for setting coeffs in
LAMMPS. This must be one of: ``bond``, ``angle``, ``dihedral``,
``improper``. In the case of ``dihedral``, LAMMPS is just referring
to proper dihedral interactions.
bonded_interactions : list of bonded_interactions (or single type)
``BondedInteractions`` which will be created in LAMMPS. This list
must only contain a single type of ``bonded_interactions`` (e.g.
only bonds), which must correspond to ``lmp_name``.
"""
special = 'no'
ID_attr = getattr(self, '{0}_ID'.format(lmp_name)
if lmp_name != 'dihedral' else 'proper_ID')
coeff_function = getattr(self.lmp, f'{lmp_name}_coeff')
# If bonds already exist, new bond IDs are generated from lowest unused
# integer
start = max(ID_attr.values()) + 1 if ID_attr else 1
for ID, b_i in enumerate(bonded_interactions, start=start):
# Create the bonded interaction coefficients
coeff_function(ID, *parse_bonded_coefficients(b_i))
# Relate each bonded interaction with its ID
ID_attr[b_i] = ID
# Create the bonded_interactions
# Special triggers the internal interaction list in LAMMPS
# This must at least occur at the end, and is an expensive
# operation
if b_i is bonded_interactions[-1]:
special = 'yes'
# LAMMPS create_bonds is used for creating all types of bonded
# interactions, by appending the lmp_name to 'single/'
c_b_type = f'single/{lmp_name}'
for atom_tpl in b_i.atoms:
atom_IDs = [self.atom_dict[atom] for atom in atom_tpl]
self.lmp.create_bonds(c_b_type,
ID,
*atom_IDs,
'special',
special)
def _update_bonded_interactions(self, lmp_name: str,
bonded_interactions: list[BondedInteraction]) -> None:
"""
Updates the bonded interaction coefficients, which are then applied to
any bonded interactions which have previously been set
Parameters
----------
lmp_name : str
The name of the bonded interaction type used for setting coeffs in
LAMMPS. This must be one of: ``bond``, ``angle``, ``dihedral``,
``improper``. In the case of ``dihedral``, LAMMPS is just referring
to proper dihedral interactions.
bonded_interactions : list of BondedInteractions
``BondedInteractions`` which will be updated in LAMMPS. This list
must only contain a single type of ``bonded_interactions`` (e.g.
only bonds), which must correspond to ``lmp_name``.
"""
# Get LAMMPS function for setting bonded interaction attributes (e.g.
# bond_coeff)
coeff_function = getattr(self.lmp, f'{lmp_name}_coeff')
# Get ID dict attribute from self (e.g. bond_ID)
b_i_IDs = getattr(self, '{0}_ID'.format(lmp_name)
if lmp_name != 'dihedral' else 'proper_ID')
for b_i in bonded_interactions:
coeff_function(b_i_IDs[b_i], *parse_bonded_coefficients(b_i))
[docs]
def apply_constraints(self) -> None:
"""
Adds a constraint ``fix`` to LAMMPS for all bonds and bond angles which
are constrained
"""
# Sort bonded interactions in the Universe which are constrained into
# bonds and angles
b_inters = set(self.universe.bonded_interactions)
# *_ is for pylint as it does not know about the output of partition_interactions
bonds, angles, *_ = partition_interactions([inter for inter
in b_inters
if getattr(inter, 'constrained', False)],
['Bond', 'BondAngle'], lst=True)
algorithm = parse_constraint(self.universe.constraint_algorithm,
bonds=bonds, bond_ID_dict=self.bond_ID,
angles=angles, angle_ID_dict=self.angle_ID)
# Create a group from all of the atom types in the constrained bonds and
# angles - the fix will be applied to this group. Applying constraint
# fix only to this group improves performance.
# chain is used to flatten inter.atoms, which is a list of tuples
atom_types = {atom.atom_type for inter in bonds+angles
for atom in chain.from_iterable(inter.atoms)}
constrain_group = 'constrain_group'
self.lmp.group(constrain_group, 'type', *atom_types)
self.lmp.fix('constrain', constrain_group, *algorithm)
[docs]
@repr_decorator('universe', 'time_step', 'traj_step', 'skin', 'neighbor_steps',
'lin_momentum_steps', 'ang_momentum_steps', 'ensemble')
class LAMMPSSimulation(PyLammpsAttribute):
# Class has to maintain a lot of state (attributes) as PyLammps class does
# not
# pylint: disable=too-many-instance-attributes
"""
The attributes and methods related running a simulation in LAMMPS using a
``LAMMPSUniverse`` object
Parameters
----------
universe : Universe
The MDMC ``Universe`` used to create the ``LAMMPSUniverse``.
traj_step : int
How many steps the simulation should take between dumping each
``CompactTrajectory`` frame
time_step : float, optional
Simulation timestep in ``fs``, default is ``1.``
lmp : PyLammps, optional
Set the ``lmp`` attribute to a ``PyLammps`` object. Default is `None`,
which results in a new ``PyLammps`` object being initialised.
**settings
``temperature`` (`float`)
``skin`` (`float`)
``neighbor_steps`` (`int`)
``remove_linear_momentum`` (`int`)
``remove_angular_momentum`` (`int`)
``velocity_seed`` (`int`): The seed to be used by LAMMPS to create velocities.
Attributes
----------
universe : Universe
An MDMC ``Universe`` object.
traj_step : int
Number of simulation steps that elapse between the ``CompactTrajectory``
being stored.
ensemble : LAMMPSEnsemble
Simulation ensemble, which applies a ``thermostat`` and ``barostat``.
"""
def __init__(self, universe, traj_step: int, time_step: float = 1., lmp=None, **settings: dict):
super().__init__(lmp, settings.get('atom_style', 'full'))
self.universe = universe
self.ensemble = LAMMPSEnsemble(self.lmp, **settings)
self.velocity_seed = settings.get('velocity_seed', randint(1, 9999))
self.temperature = settings.get('temperature')
self.traj_step = traj_step
self.time_step = time_step
self.skin = settings.get('skin', 2.0)
self.neighbor_steps = settings.get('neighbor_steps', 1)
# Setting _lin_momentum_steps and _ang_momentum_steps allows
# _set_momentum_removers to be called when setting
# self.lin_momentum_steps and self.ang_momentum_steps
self._lin_momentum_steps = None
self._ang_momentum_steps = None
self.lin_momentum_steps = settings.get('remove_linear_momentum', 1)
self.ang_momentum_steps = settings.get('remove_angular_momentum')
self._set_kspace_solver()
@property
def time_step(self) -> float:
"""
Get or set the simulation time step in ``fs``
Returns
-------
`float`
Simulation time step in ``fs``
"""
return self._time_step
@time_step.setter
@unit_decorator(unit=units.TIME)
def time_step(self, value: float) -> None:
self._time_step = value
with suppress(ValueError):
# Set the timestep in LAMMPS wrapper
self.lmp.timestep(convert_unit(self._time_step))
@property
def temperature(self) -> float:
"""
Get or set the temperature of the simulation in ``K``
Returns
-------
`float`
Temperature in ``K``
"""
return self._temperature
@temperature.setter
@unit_decorator(unit=units.TEMPERATURE)
def temperature(self, value: float) -> None:
self._temperature = value
self.ensemble.temperature = value
try:
# Set the initial temperature in the LAMMPS wrapper
if self.system_state.natoms > 0:
zero_velocity = [np.array_equal(atom.velocity, (0, 0, 0))
for atom in self.universe.atoms]
if all(zero_velocity):
# If we have not set any velocities (they are all the default value of zero)
# then "create" a velocity for each atom according to user-specified or
# random seed
self.lmp.velocity('all', 'create',
convert_unit(self._temperature), self.velocity_seed)
else:
if any(zero_velocity):
msg = ('Some but not all atom velocities set. Atoms with non-zero velocity'
' will be re-scaled to match target temperature, atoms with zero '
'velocity will remain stationary.')
LOGGER.info('%s %s', self.__class__, msg)
print(msg)
# If we have set velocities then "scale" the velocities we have to the correct
# temperature
self.lmp.velocity(
'all', 'scale', convert_unit(self._temperature))
except ValueError:
pass
@property
def pressure(self) -> float:
"""
Get or set the pressure of the simulation in ``atm``
Returns
-------
`float`
Pressure in ``atm``
"""
return self.ensemble.pressure
@pressure.setter
@unit_decorator(unit=units.PRESSURE)
def pressure(self, value: float) -> None:
self.ensemble.pressure = value
@property
def thermostat(self) -> str:
"""
Get or set the string which specifies the thermostat
Returns
-------
`str`
The thermostat name
"""
return self.ensemble.thermostat
@thermostat.setter
def thermostat(self, value: str) -> None:
self.ensemble.thermostat = value
@property
def barostat(self) -> str:
"""
Get or set the string which specifies the barostat
Returns
-------
`str`
The barostat name
"""
return self.ensemble.barostat
@barostat.setter
def barostat(self, value: str) -> None:
self.ensemble.barostat = value
@property
def skin(self) -> float:
"""
Get or set the skin distance in ``Ang``
Returns
-------
`float`
The skin distance in ``Ang``. This distance plus the force
``cutoff`` distance determines which atom pairs are stored in the
neighbor list.
"""
return self._skin
@skin.setter
@unit_decorator(unit=units.LENGTH)
def skin(self, value: float) -> None:
self._skin = value
# Set the neighor list parameters in the LAMMPS wrapper
self.lmp.neighbor(convert_unit(self._skin), 'bin')
@property
def neighbor_steps(self) -> int:
"""
Get or set the number of steps between neighbor list updates
Returns
-------
`int`
Number of steps between neighbor list updates
"""
return self._neighbor_steps
@neighbor_steps.setter
def neighbor_steps(self, value: int) -> None:
self._neighbor_steps = value
with suppress(TypeError):
# Set the modifiers to the neighbor list in the LAMMPS wrapper
self.lmp.neigh_modify('every', int(self.neighbor_steps), 'delay', 0,
'check', 'yes')
@property
def lin_momentum_steps(self) -> int:
"""
Get or set the number of steps between resetting the linear momentum
Returns
-------
`int`
Number of steps between the linear momentum being removed
"""
return self._lin_momentum_steps
@lin_momentum_steps.setter
def lin_momentum_steps(self, value: int) -> None:
self._lin_momentum_steps = value
# Set the momentum removers in the LAMMPS wrapper
self._set_momentum_removers()
@property
def ang_momentum_steps(self) -> int:
"""
Get or set the number of steps between resetting the angular momentum
Returns
-------
`int`
Number of steps between the angular momentum being removed
"""
return self._ang_momentum_steps
@ang_momentum_steps.setter
def ang_momentum_steps(self, value: int) -> None:
self._ang_momentum_steps = value
# Set the momentum removers in the LAMMPS wrapper
self._set_momentum_removers()
def _set_momentum_removers(self) -> None:
"""
Creates the ``fixes`` in LAMMPS which remove the linear and angular
momentum of the simulation
Removes all pre-existing momentum remover ``fixes``
"""
for name in self.fix_names:
if name in ['RemoveMomentum', 'RemoveLinearMomentum',
'RemoveAngularMomentum']:
self.lmp.unfix(name)
if self.system_state.natoms > 0:
if self.lin_momentum_steps == self.ang_momentum_steps:
self.lmp.fix('RemoveMomentum', 'all', 'momentum',
self.lin_momentum_steps, 'linear', 1, 1, 1, 'angular')
else:
if self.lin_momentum_steps:
self.lmp.fix('RemoveLinearMomentum', 'all', 'momentum',
self.lin_momentum_steps, 'linear', 1, 1, 1)
if self.ang_momentum_steps:
self.lmp.fix('RemoveAngularMomentum', 'all', 'momentum',
self.ang_momentum_steps, 'angular')
def _set_kspace_solver(self) -> None:
"""
Creates a k-space solver in LAMMPS using ``kspace_style``, if one is
required
Uses either the ``kspace_solver``, the ``electrostatic_solver`` or both
the ``electrostatic_solver`` and ``dispersive_solver`` attribute of the
MDMC ``Universe`` to set the ``kspace_style``. Note that LAMMPS does not
support different electrostatic and dispersive solvers. Setting with
equivalent electrostatic and dispserive solvers is equivalent to setting
with ``kspace_solver``. LAMMPS also does not support just applying a
dispersive solver.
Raises
------
TypeError
If ``self.universe`` has both a ``kspace_solver`` and either an
``electrostatic_solver`` or a ``dispersive_solver``.
TypeError
If ``self.universe`` has a different ``electrostatic_solver`` and
``dispersive_solver``.
TypeError
If ``self.universe`` only has a ``dispersive_solver``.
"""
kspace = self.universe.kspace_solver
electrostatic = self.universe.electrostatic_solver
dispersive = self.universe.dispersive_solver
# LAMMPS supports a single kspace solver, which can be taken from kspace
# electrostatic or dispersive solvers of the universe. If any other
# solver is defined as well as kspace (which effecitvely defines both
# other solves), raise an error.
if kspace and (electrostatic or dispersive):
raise TypeError('kspace solver cannot be applied in conjunction'
' with either electrostatic or dispersive solvers')
if kspace:
self.lmp.kspace_style(*parse_kspace_solver(kspace))
# If either an electrostatic or a dispersive solver has been defined,
# use this. If both have been defined, this is valid as long as they are
# equivalent. If not, raise an error.
if electrostatic:
if dispersive and electrostatic != dispersive:
raise TypeError('LAMMPS only supports a single solver, so if'
' both dispersive and electrostatic solvers'
' have been defined then they must be'
' equivalent')
self.lmp.kspace_style(*parse_kspace_solver(electrostatic))
elif dispersive:
raise TypeError('LAMMPS does not support only applying a dispersive'
' solver - it must be applied in conjunction with'
' an electrostatic solver')
[docs]
@repr_decorator('temperature', 'pressure', 'thermostat', 'barostat')
class LAMMPSEnsemble(PyLammpsAttribute):
# Class has to maintain a lot of state (attributes) as PyLammps class does
# not
# pylint: disable=too-many-instance-attributes
"""
A thermodynamic ensemble determined by applying a thermostat and/or barostat
Parameters
----------
lmp : PyLammps
Set the ``lmp`` attribute to a ``PyLammps`` object.
temperature : float, optional
Thermostat temperature. Default is `None`, which is only valid if a
``thermostat`` is also `None`.
pressure : float, optional
Barostat pressure. Default is `None`, which is only valid if a
``barostat`` is also `None`.
thermostat : str
Name of a thermostat to be applied.
barostat : str
Name of a barostat to be applied.
**settings
``time_step`` (`float`)
``t_damp`` (`int`)
``p_damp`` (`int`)
``t_window`` (`float`)
``t_fraction`` (`float`)
``rescale_step`` (`int`)
Attributes
----------
rescale_step : int
Number of steps between applying temperature rescaling. This only
applies to rescale thermostats.
"""
def __init__(self, lmp, temperature: float = None, pressure: float = None,
thermostat: str = None, barostat: str = None, **settings: dict):
# Requires a lmp object as thermostats cannot be applied before
# configuration is defined
super().__init__(lmp)
# Setting _thermostat and _barostat allows apply_ensemble_fixes to be
# called when setting self.temperature and self.pressure
self._thermostat = None
self._barostat = None
self.temperature = temperature
self.pressure = pressure
self.time_step = settings.get('time_step', 1.0)
self.t_damp = settings.get('t_damp', 100)
self.p_damp = settings.get('p_damp', 1000)
self.t_window = settings.get('t_window')
self.t_fraction = settings.get('t_fraction')
self.rescale_step = settings.get('rescale_step')
self.thermostat = thermostat
self.barostat = barostat
@property
def time_step(self) -> float:
"""
Get or set the simulation time step in ``fs``
Returns
-------
`float`
Simulation time step in ``fs``
"""
return self._time_step
@time_step.setter
@unit_decorator(unit=units.TIME)
def time_step(self, value):
self._time_step = value
@property
def temperature(self) -> float:
"""
Get or set the temperature of the simulation in ``K``
"""
return self._temperature
@temperature.setter
@unit_decorator(unit=units.TEMPERATURE)
def temperature(self, value: float) -> None:
self._temperature = value
self.apply_ensemble_fixes()
@property
def pressure(self) -> float:
"""
Get or set the pressure of the simulation in ``atm``
"""
return self._pressure
@pressure.setter
@unit_decorator(unit=units.PRESSURE)
def pressure(self, value: float) -> None:
self._pressure = value
self.apply_ensemble_fixes()
# Unit has to be applied to getter due to operation in setter
@property
@unit_decorator_getter(unit=units.TIME)
def t_damp(self) -> int:
"""
Get or set the number of time steps over which the ``temperature`` is
relaxed
Required for ``Nose-Hoover``, ``Berendsen`` and ``Langevin``
thermostats. The ``time_step`` must have been set before ``t_damp``.
Returns
-------
`int`
Number of time steps
Raises
------
AttributeError
If ``self.time_step`` has not been set.
"""
# t_damp is stored in units of time - convert back to number of steps
# here
try:
return self._t_damp / self.time_step
except TypeError:
return self._t_damp
@t_damp.setter
def t_damp(self, value: int) -> None:
try:
# LAMMPS requires t_damp to be given in units of time, but it is
# more natural to give it in units of steps - convert between them
# here
self._t_damp = value * self.time_step
except TypeError as error:
if value is None:
self._t_damp = value
else:
raise AttributeError('the time_step attribute must be set'
' before t_damp') from error
# Unit has to be applied to getter due to operation in setter
@property
@unit_decorator_getter(unit=units.TIME)
def p_damp(self) -> int:
"""
Get or set the number of time steps over which the ``pressure`` is
relaxed
Required for ``Nose-Hoover`` or ``Berendsen`` barostats. The
``time_step`` must have been set before ``p_damp``.
Returns
-------
`int`
Number of time steps
Raises
------
AttributeError
If ``self.time_step`` has not been set.
"""
# p_damp is stored in units of time - convert back to number of steps
# here
try:
return self._p_damp / self.time_step
except TypeError:
return self._p_damp
@p_damp.setter
def p_damp(self, value: int) -> None:
try:
# LAMMPS requires p_damp to be given in units of time, but it is
# more natural to give it in units of steps - convert between them
# here
self._p_damp = value * self.time_step
except TypeError as error:
if value is None:
self._p_damp = value
else:
raise AttributeError('the time_step attribute must be set'
' before p_damp') from error
@property
def t_fraction(self) -> float:
"""
Get or set the fraction by which the ``temperature`` is rescaled to the
target temperature
This is required for the rescale ``thermostat``.
Returns
-------
`float`
Fraction (i.e. between 0.0 and 1.0 inclusive) by which the
``temperature`` is rescaled
Raises
------
ValueError
If set to a value outside of 0.0 and 1.0 inclusive.
"""
return self._t_fraction
@t_fraction.setter
def t_fraction(self, value: float) -> None:
# Must be a fraction
if value is None or 0. <= value <= 1.0:
self._t_fraction = value
else:
raise ValueError('the t_fraction must be between 0.0 and 1.0')
@property
def t_window(self) -> float:
"""
Get or set the ``temperature`` range in ``K`` in which the
``temperature`` is not rescaled
This only applies to rescale ``thermostats``.
Returns
-------
`float`
temperature range in ``K``
"""
return self._t_window
@t_window.setter
@unit_decorator(unit=units.TEMPERATURE)
def t_window(self, value: float) -> None:
self._t_window = value
@property
def thermostat(self) -> str:
"""
Get or set the `str` which specifies the thermostat
Raises
------
AttributeError
If ``self.temperature`` has not been set.
"""
return self._thermostat
@thermostat.setter
def thermostat(self, value: str) -> None:
if value and not self.temperature:
raise AttributeError('all ensembles with a thermostat must have a'
' temperature')
self._thermostat = value
# Set the thermostat and barostat in LAMMPS wrapper
self.apply_ensemble_fixes()
@property
def barostat(self) -> str:
"""
Get or set the `str` which specifies the barostat
Raises
------
AttributeError
If ``self.pressure`` has not been set.
"""
return self._barostat
@barostat.setter
def barostat(self, value: str) -> None:
if value and not self.pressure:
raise AttributeError('all ensembles with a barostat must have a'
' pressure')
self._barostat = value
# Set the thermostat and barostat in LAMMPS wrapper
self.apply_ensemble_fixes()
[docs]
def remove_ensemble_fixes(self) -> None:
"""
Removes all LAMMPS ``fixes`` relating to the ensemble i.e. removes all
thermostats and barostats
This must be done before thermostat and barostat fixes are added, so
that there is no conflict with existing thermostat and barostats
``fixes``. It is also required for Shake and Rattle ``fixes`` which
cannot be added after barostat fixes have been applied.
"""
for name in self.fix_names:
if name in ['nve', 'nvt', 'npt', 'nph', 't_berendsen',
'p_berendsen', 'langevin', 'rescale']:
LOGGER.debug('%s Remove %s fix.',
self.__class__,
name)
self.lmp.unfix(name)
[docs]
def apply_ensemble_fixes(self) -> None:
"""
Passes the required LAMMPS ``fixes`` to apply a specific thermodynamic
ensemble to the simulation
Removes all pre-existing thermostat and barostat fixes
"""
self.remove_ensemble_fixes()
LOGGER.debug('%s apply_ensemble_fixes before fixes applied:'
' {fixes: %s}',
self.__class__,
self.fixes)
if not self.thermostat and not self.barostat:
self.lmp.fix('nve', 'all', 'nve')
else:
thermo_parameters = []
press_parameters = []
if self.thermostat:
temp = convert_unit(self.temperature)
if self.thermostat != 'rescale':
t_damp = convert_unit(self.t_damp)
thermo_parameters = [temp, temp, t_damp]
if self.barostat:
press = convert_unit(self.pressure)
p_damp = convert_unit(self.p_damp)
press_parameters = ['iso', press, press, p_damp]
# Apply thermostat
# we create local funcs for each thermostat, and use a dict to get the correct one.
# think of this like a proto-factory pattern, or like a switch-case block in C++.
def nose():
if self.barostat == 'nose':
self.lmp.fix('npt', 'all', 'npt', 'temp',
*thermo_parameters + press_parameters)
else:
self.lmp.fix('nvt', 'all', 'nvt',
'temp', *thermo_parameters)
def berendsen():
# berendsen does not do time integration so also requires nve
self.lmp.fix('nve', 'all', 'nve')
self.lmp.fix('t_berendsen', 'all', 'temp/berendsen',
*thermo_parameters)
def langevin():
# langevin does not do time integration so also requires nve
self.lmp.fix('nve', 'all', 'nve')
self.lmp.fix('langevin', 'all', 'langevin',
*thermo_parameters + [randint(0, 9999)])
def rescale():
# temp/rescale does not do time integration so also requires nve
t_window = convert_unit(self.t_window)
self.lmp.fix('nve', 'all', 'nve')
self.lmp.fix('rescale', 'all', 'temp/rescale',
self.rescale_step, temp, temp, t_window,
self.t_fraction)
def csvr():
# csvr does not do time integration so also requires nve
self.lmp.fix('nve', 'all', 'nve')
self.lmp.fix('csvr', 'all', 'temp/csvr',
*thermo_parameters + [randint(0, 9999)])
thermostat = {'nose': nose,
'berendsen': berendsen,
'langevin': langevin,
'rescale': rescale,
'csvr': csvr}
try:
thermostat[self.thermostat]()
except KeyError:
warnings.warn("Thermostat type not recognised. No thermostat will be applied."
" Your thermostat type may be misspelt or not currently implemented.")
# Apply barostat
if self.barostat == 'berendsen':
self.lmp.fix('p_berendsen', 'all', 'press/berendsen',
*press_parameters)
elif self.barostat == 'nose' and self.thermostat != 'nose':
if 'nve' in self.fix_names:
self.lmp.unfix('nve')
self.lmp.fix('nph', 'all', 'nph', *press_parameters)
LOGGER.debug('%s apply_ensemble_fixes after fixes applied:'
' {fixes: %s}',
self.__class__,
self.fixes)
# Define the unit system used in LAMMPS
# NB: LAMMPS uses deg for angle but radian for derived quantities of angle:
# e.g. harmonic angle potential strength is in kcal / mol radian ^ 2
SYSTEM = {
'LENGTH': units.Unit('Ang'),
'TIME': units.Unit('fs'),
'MASS': units.Unit('g') / units.Unit('mol'),
'CHARGE': units.Unit('e'),
'ANGLE': units.Unit('deg'),
'TEMPERATURE': units.Unit('K'),
'ENERGY': units.Unit('kcal') / units.Unit('mol'),
'FORCE': units.Unit('kcal') / (units.Unit('Ang') * units.Unit('mol')),
'PRESSURE': units.Unit('atm'),
}
[docs]
def convert_unit(value: Union[np.ndarray, float],
unit: Unit = None, to_lammps: bool = True):
"""
Converts between MDMC units and LAMMPS real units
Parameters
----------
value : array_like or float_like
The value of the physical property to be converted, in MDMC units.
Must derive from either ``ndarray`` or `float`.
unit : Unit, optional
The unit of the ``value``. If `None`, the ``value`` must possess a
``unit`` attribute i.e. derive from ``UnitFloat`` or ``UnitArray``.
Default is `None`.
to_lammps : bool, optional
If `True` the conversion is from MDMC units to LAMMPS units. Default is
`True`.
Returns
-------
`float` or `numpy.ndarray`
Value in LAMMPS units if ``to_lammps`` is `True`, otherwise value in
MDMC units. Return type is same as ``value`` type.
"""
def expand_components(unit: Unit, system: dict) -> tuple:
"""
Expands out the ``components`` of a ``Unit``, so that the ``Unit`` is
expressed purely in terms of ``base`` ``Unit`` objects. The only
exception to this is ``Unit`` objects which occur in ``System``: these
are kept in the `list` of ``components``.
Parameters
----------
unit : Unit
The ``Unit`` to be expanded out
system : dict
A `dict` of {``Property``: ``Unit``} pairs, which is used to
substitute ``System`` ``Unit`` objects into the expanded
``components``.
Returns
-------
`tuple`
A `tuple` of (``num``, ``denom``), where ``num`` is a `list` of all
``base`` ``Unit`` objects in the numerator, and ``denom`` is a
`list` of all ``base`` ``Unit`` objects in the denominator
"""
def is_sublist_of_list(sub: list, lst: list) -> bool:
"""
Determines if all of the elements in a sublist are in a `list`,
including ensuring that any duplicates in the sublist have at least
the same number of duplicates in the `list`
Parameters
----------
sub : list
The sublist to be tested
lst : list
The `list` the ``sub`` is tested against
Returns
-------
`bool`
`True` if all of the elements of ``sub`` are in ``lst``, and
`False` if not
"""
return all(sub.count(x) <= lst.count(x) for x in set(sub))
def remove_components(remove_comps: list, comps: list) -> list:
"""
Removes all elements of a `list` of ``components`` from another
`list` of ``components``
Parameters
----------
remove_comps : list
The ``components`` to be removed
comps : list
The ``components`` from which ``remove_comps`` is removed
Returns
-------
`list`
A `list` of ``components`` from which ``remove_comps`` have been
removed
"""
for remove_comp in remove_comps:
comps.remove(remove_comp)
return comps
num, denom = [], []
# If unit is in system of units, don't break down to components, as
# a unit in the system of units always has a conversion, even if it is
# a compound unit.
if unit.base or unit in system.values():
num.append(unit)
else:
# tuple unpacking style below appends to num and denom lists
for comp in unit.components['numerator']:
num[len(num):], denom[len(denom):] = expand_components(comp,
system)
for comp in unit.components['denominator']:
denom[len(denom):], num[len(num):] = expand_components(comp,
system)
# Substitute in units rather than separated out components if a unit
# exists in the system of units. This is because the conversion can
# always be determined for units in the system of units. Base units
# are filtered as they can only replace themselves (as they are
# their only component).
for sys_unit in (unit for unit in system.values() if not unit.base):
# while is required for powers of sys_unit, as otherwise only
# a single power will be removed
while True:
sys_unit_num = sys_unit.components['numerator']
sys_unit_denom = sys_unit.components['denominator']
# Determine if all of the components of unit are in the
# numerator and denominator lists. If so, remove them and
# replace with the unit in the numerator.
if (is_sublist_of_list(sys_unit_num, num) and
is_sublist_of_list(sys_unit_denom, denom)):
num = remove_components(sys_unit_num, num)
denom = remove_components(sys_unit_denom, denom)
num.append(sys_unit)
# Do the same for the inverse (i.e. unit's numberator
# components in the denominator list and vice versa). If so,
# remove them and replace with the unit in the denominator.
elif (is_sublist_of_list(sys_unit_num, denom) and
is_sublist_of_list(sys_unit_denom, num)):
denom = remove_components(sys_unit_num, denom)
num = remove_components(sys_unit_denom, num)
denom.append(sys_unit)
# Breaks if the components of sys_unit are not found in num
# and denom
else:
break
return num, denom
# If no unit argument is passed, the value must possess a unit
if not unit:
# If value is unitless, no conversion is required
try:
unit = value.unit
except AttributeError as error:
if value is None:
raise ValueError('Cannot convert NoneType value') from error
return value
# Expand the unit in terms of its base units (for numerator and denominator)
if to_lammps:
# SYSTEM represents the LAMMPS system units, units.SYSTEM the MDMC
# system units
l_sys = copy(SYSTEM)
mdmc_sys = copy(units.SYSTEM)
# For angular potential strength LAMMPS uses units derived from 'rad',
# namely 'kcal / mol rad^2', but uses 'deg' when measuring angles
# themselves. As a result 'deg' is recorded as their system unit. In
# order to prevent us from erroneously expressing potential strength in
# 'kcal / mol deg^2', we must override the LAMMPS system unit for angle
# ONLY in the case where we are dealing with angular potential, namely
# when the MDMC unit is measuring 'energy / angle^2'.
if (unit in (mdmc_sys['ENERGY'] / mdmc_sys['ANGLE'] ** 2,
mdmc_sys['ENERGY'] / units.Unit('rad') ** 2)):
l_sys['ANGLE'] = units.Unit('rad')
expanded_unit = expand_components(unit, mdmc_sys)
# For each MDMC unit, determine the LAMMPS system unit that corresponds
# to that physical property (e.g. 'kJ / mol' -> 'ENERGY' -> 'kcal / mol')
l_unit_nums, l_unit_denoms = map(lambda unit_comps: [l_sys[comp.physical_property]
for comp in unit_comps],
expanded_unit)
# In general we may have an MDMC measurement that is not wholly
# expressed in MDMC system units, for example 'kJ / mol rad^2' uses
# 'rad' rather than 'deg'. To account for this, multiply by the
# conversion_factor of the original MDMC unit to get the value in terms
# of MDMC system units only.
value *= unit.conversion_factor
# Then, for each component of the LAMMPS unit, divide by the
# conversion_factor to go from MDMC system units to LAMMPS system units.
for component in l_unit_nums:
value /= component.conversion_factor
for component in l_unit_denoms:
value *= component.conversion_factor
else:
value *= unit.conversion_factor
return value
[docs]
def parse_bonded_styles(interaction: BondedInteraction) -> str:
"""
Converts MDMC ``InteractionFunction`` names for ``BondedInteractions`` to
LAMMP bond styles
Parameters
----------
interaction : BondedInteraction
``BondedInteraction`` to be parsed into LAMMPS bond style.
Returns
-------
`str`
LAMMPS bond style corresponding to ``interaction``
Raises
------
NotImplementedError
If ``interaction`` has a ``function`` that has not been implemented in
the LAMMPS facade.
"""
if interaction.function_name == 'HarmonicPotential':
if interaction.name == 'DihedralAngle' and not interaction.improper:
raise TypeError(
'LAMMPS does not support harmonic proper dihedrals')
return 'harmonic'
if interaction.function_name == 'Periodic':
if interaction.name != 'DihedralAngle':
raise TypeError('LAMMPS only supports periodic interaction function'
' for Dihedrals')
if interaction.improper:
# LAMMPS has an improper_style called cvff, which is a Simplified
# version of Periodic (it is only first order and d1=0). Return cvff
# here and deal with incompatible Periodic interactions (e.g. d1!=0)
# in parse_bonded_coefficients
return 'cvff'
return 'fourier'
raise NotImplementedError('This InteractionFunction has not been'
' implemented in the LAMMPS facade')
[docs]
def parse_nonbonded_styles(interaction: NonBondedInteraction) -> tuple:
"""
Converts MDMC ``InteractionFunction`` names for ``NonBondedInteractions`` to
LAMMPS pair styles
Parameters
----------
interaction : NonBondedInteraction
``NonBondedInteraction`` to be parsed into LAMMPS pair style.
Returns
-------
`tuple`
(``styles``, ``mods``) where ``styles`` is a `list` of `str` of LAMMPS
pair style corresponding to the ``interaction``, and ``mods`` is a
`list` of `str` of LAMMPS pair modifications corresponding to the
interaction
Raises
------
NotImplementedError
If ``interaction`` has a ``function`` that has not been implemented in
the LAMMPS facade.
"""
lmp_str = []
if interaction.function_name == 'Buckingham':
lmp_str.append('buck')
elif interaction.function_name == 'LennardJones':
lmp_str.append('lj')
elif interaction.function_name == 'Coulomb':
lmp_str.append('coul')
else:
raise NotImplementedError('This InteractionFunction has not been'
' implemented in the LAMMPS facade')
if interaction.cutoff:
cutoff = convert_unit(interaction.cutoff)
kspace = interaction.universe.kspace_solver
electrostatic = interaction.universe.electrostatic_solver
dispersive = interaction.universe.dispersive_solver
if (kspace or (electrostatic and interaction.name == 'Coulombic')
or (dispersive and interaction.name == 'Dispersion')):
lmp_str[-1] += '/long'
else:
if interaction.function_name != 'Buckingham':
lmp_str[-1] += '/cut'
lmp_str.append(cutoff)
else:
raise AttributeError(f'You have not specified a `cutoff` for the'
f'{interaction} InteractionFunction.')
return lmp_str, parse_nonbonded_modifications(interaction)
[docs]
def parse_nonbonded_modifications(interaction: Interaction) -> 'list[str]':
"""
Parses MDMC ``Interaction`` attributes into `list` that can be used with
LAMMPS ``pair_modify`` command
Attributes
----------
interaction : Interaction
The ``Interaction`` for which the attributes are parsed
Returns
-------
`list`
A `list` of `str` which are a valid ``pair_modify`` command, or an empty
`list` if no ``Interaction`` attributes require setting ``pair_modify``
"""
# LAMMPS pair_modify is of the following form:
# pair_modify('pair', 'lj/cut', 'mix', 'geometric', 'tail', 'yes')
# pair_modify('coul/long', 'mix', 'arithmetic')
# where each pair style (lj/cut, coul/long etc) has any mix or tail keywords
# defined after the pair style. If the same pair_style occurs multiple times
# but with different modifiers
# (e.g. 'lj/cut', 'mix', 'geometric' and 'lj/cut', 'mix', 'arithmetic')
# whichever pair_modify occurs last will be applied to all identical
# pair_styles.
mods = []
if interaction.name == 'Dispersion' and interaction.vdw_tail_correction:
mods.extend(['tail yes'])
return mods
[docs]
def parse_all_nonbonded_styles(interactions: NonBondedInteraction) -> 'dict[tuple, list[str]]':
"""
Converts all ``NonBondedInteractions`` to LAMMPS pair styles
This is required because LAMMPS frequently treats ``Coulombic`` and
``Disperion`` interactions together; these cases need to be dealt with to
generate the correct input to LAMMPS ``pair_styles``. For example, while
the ``pair_styles`` ``buck``, ``lj/cut``, ``coul/cut`` and ``coul/long`` can
all be passed separately, the dispersive and coulombic styles are combined
(dispersive always preceding coulombic) as this is dealt with more
efficiently in the LAMMPS engine. ``buck/long`` and ``lj/long`` only exist
as part of other pair styles, such as ``buck/long/coul/long`` and
``lj/long/coul/long``, so must be combined.
Parameters
----------
interactions : list of NonBondedInteractions
``NonBondedInteractions`` to be parsed into LAMMPS ``pair_styles``.
Returns
-------
`dict`
A `dict` with {``pair_styles``: ``pair_modifiers``} where
``pair_styles`` is a `tuple` of all of the combined LAMMPS pair styles
corresponding to ``interactions``. Each `tuple` contains the combined
``pair_style`` as one element, and the ``cutoffs`` as the second element
e.g. ``('lj/cut/coul/cut', 5.0 10.0)``. If additional arguments are
required when combining styles, they will be added and returned in the
`tuple` as the third and middle element
e.g. ``('buck/long/coul/long', 'long long', '10.0'``).
``pair_modifiers`` is a `list` of `str` for setting ``pair_modify`` for
the corresponding ``pair_style``.
Raises
------
ValueError
If ``Dispersion`` and ``Coulombic`` interactions have different long
range ``cutoff``, which is not implemented in LAMMPS.
ValueError
If long range ``Dispersion`` not defined in conjunction with a long
range ``Coulombic``.
"""
def check_validity(pair_style: str, cutoffs: 'list[float]' = None) -> 'list[str]':
"""
Tests the validity of a LAMMPS ``pair_style``.
Parameters
----------
pair_style : str
The LAMMPS ``pair_style`` to be checked.
cutoffs : list of float
The cutoff distances for each style in the ``pair_style``.
Returns
-------
`list` of `str`
Contains the inputted ``pair_style`` (if valid) as well as any extra
terms required in specific cases (i.e. ``'long long'`` if the pair
style ``'buck/long/coul/long'`` is passed).
Raises
------
ValueError
If the ``pair_style`` passed is not valid.
"""
if pair_style in ['buck/long/coul/cut', 'lj/long/coul/cut']:
raise ValueError('Invalid pair_style: ' + pair_style
+ '. LAMMPS requires long range Coulombics to be'
' defined in conjunction with long range Dispersion'
' interactions.')
if pair_style in ['buck/long/coul/long', 'lj/long/coul/long']:
if cutoffs[0] != cutoffs[1]:
raise ValueError('LAMMPS requires both cutoffs to be the same'
' for long range Dispersion and Coulombic pair'
' styles.')
return [pair_style, 'long long']
return [pair_style]
# Some pair styles cannot have vdw tail corrections modifiers - exclude
# these and warn about this exclusion
# This is dealth with here rather than in parse_nonbonded_modifications as
# it is dependent on combined pair_styles
excluded = ['lj/long/coul/long']
# Change parse_nonbonded_styles to return the parsed nonbonded style and the
# parsed_ modifiers from that interaction
# Create a dictionary of parsed_interactions (keys) and modifiers (values)
single_parsed_inters = {tuple(parsed[0]): parsed[1] for parsed in
[parse_nonbonded_styles(nb) for nb in interactions]}
combined_parsed_inters = copy(single_parsed_inters)
# Define all Dispersion and Coulombic styles currently supported
# Dispersion styles always precede coulombic styles in LAMMPS pair styles
disp_styles = ['buck', 'buck/long', 'lj/cut', 'lj/long']
coul_styles = ['coul/cut', 'coul/long']
for d_style, c_style in product(disp_styles, coul_styles):
dc_styles = [d_style, c_style]
for int1, int2 in combinations(single_parsed_inters.keys(), 2):
if (int1[0] in dc_styles and int2[0] in dc_styles):
indiv_cmd = check_validity('/'.join(dc_styles),
cutoffs=[int1[1], int2[1]])
# Format cutoffs so that disp cutoff precedes coul cutoff if
# they are different
if int1[1] == int2[1]:
cutoffs = str(int1[1])
else:
d_cut, c_cut = ((int1[1], int2[1]) if int1[0] == d_style
else (int2[1], int1[1]))
cutoffs = f'{d_cut} {c_cut}'
# Add indiv_cmd to parsed_interactions dict instead. Use
# modifier from parsed_interactions as value. set is used to
# remove duplicated modifiers which occur for both interactions.
# Delete int1 and int2 from parsed_interactions
mod = set(single_parsed_inters[int1]
+ single_parsed_inters[int2])
# Warn if incompatible pair_style and pair_modification have
# been used
if indiv_cmd in excluded and 'tail yes' in mod:
LOGGER.warning('parse_all_nonbonded_styles: %s pair style'
' cannot have a vdw tail correction applied',
indiv_cmd)
mod = set(md for md in mod if md != 'tail yes')
combined_parsed_inters[tuple(
indiv_cmd + [cutoffs])] = list(mod)
for key in [int1, int2]:
with suppress(KeyError):
del combined_parsed_inters[key]
return combined_parsed_inters
[docs]
def parse_bonded_coefficients(interaction: BondedInteraction) -> 'list[str]':
"""
Orders MDMC ``Parameter`` objects for input to LAMMPS ``bond_coeff`` and
``angle_coeff``
Parameters
----------
interaction : BondedInteraction
``BondedInteraction`` where its style and parameters will be parsed.
Returns
-------
`list` of `str`
Style and parameters converted to the input format for LAMMPS
``bond_coeff`` and ``angle_coeff``
Raises
------
NotImplementedError
If ``interaction`` has a ``function`` that has not been implemented in
the LAMMPS facade.
"""
parameters = {p.type: convert_unit(p.value)
for p in interaction.parameters.as_array}
style = parse_bonded_styles(interaction)
if style == 'harmonic':
ordered_parameters = [parameters['potential_strength'],
parameters['equilibrium_state']]
elif style == 'fourier':
# There are three parameters (K, n and d) for each order of the equation
fourier_order = int(len(interaction.parameters) / 3)
# LAMMPS requires parameters to be ordered K1, n1, d1, K2, n2, d2 etc
# So the letter sequence is:
ord_let = ('K', 'n', 'd')
# Sort by fourier_order first and then by letter sequence
ordered_p_names = sorted(parameters,
key=lambda p_type: (p_type[1],
ord_let.index(p_type[0])))
# Get parameters values and prepend the order of the equation
ordered_parameters = [fourier_order] + [parameters[p_name] for p_name
in ordered_p_names]
elif style == 'cvff':
if len(parameters) > 3:
raise TypeError('LAMMPS improper dihedrals can only have a Periodic'
' interaction function with first order'
' coefficients (e.g. K1, n1, d1)')
if parameters['d1'] != 0. and parameters['d1'] != 180.:
raise TypeError('LAMMPS improper dihedrals can only have a Periodic'
' interaction with d1 = 0 deg or d = 180 deg')
if not 0 <= parameters['n1'] <= 6:
raise TypeError('LAMMPS improper dihedrals can only have a Periodic'
' interaction 0 <= n1 <= 6')
# fourier and cvff have different d parameters. The definition of
# Periodic is the same as the fourier style, so need to convert to cvff
# d parameter (which can only take values of 1 or -1)
cvff_d = 1 if parameters['d1'] == 0. else -1
ordered_parameters = [parameters['K1'], cvff_d, parameters['n1']]
else:
raise NotImplementedError('This InteractionFunction has not been'
' implemented in the LAMMPS facade')
return [style] + ordered_parameters
[docs]
def parse_dispersion_coefficients(interactions: list[NonBondedInteraction],
nonbonded_styles: 'list[tuple]' = None) -> 'list[str]':
"""
Orders MDMC ``Parameter`` objects for input to LAMMPS ``pair_coeff``
Parameters
----------
interactions : list of NonBondedInteraction
``NonBondedInteractions`` where their style and parameters will be
parsed.
nonbonded_styles : list of tuple
The parsed ``NonBondedInteractions``, where the combined ``pair_styles``
have been created by the ``parse_all_nonbonded_styles`` function.
If `None` (default), the ``interactions`` are parsed manually.
Returns
-------
`list` of `str`
Each element is a partial ``pair_coeff`` command, containing a
``pair_style``, the ordered ``Dispersion`` ``Parameter`` objects
converted to the correct input format for LAMMPS ``pair_coeff``, and the
``cutoff``.
For example, if interactions contains a ``Buckingham`` with
``Parameter`` ``A, B, C = 4.184, 1.0, 4.184``, and ``cutoff = 20.0``,
and a ``Coulombic`` with a ``cutoff`` of ``10.0``, the `str` element of
the `list` will be::
'buck/coul/cut 1.0 1.0 1.0 20.0 10.0'
Raises
------
NotImplementedError
If ``interaction`` has a ``function`` that has not been implemented in
the LAMMPS facade.
ValueError
If the LAMMPS Buckingham parameter rho is less than or equal to zero.
"""
if not nonbonded_styles:
nonbonded_styles = parse_all_nonbonded_styles(interactions)
coul_styles = ['coul/cut', 'coul/long']
coeff_cmds = []
for style in nonbonded_styles:
pair_style = style[0]
cutoffs = style[-1]
if 'buck' in pair_style:
for inter in interactions:
if inter.function.name == 'Buckingham':
parameters = {p.type: convert_unit(p.value)
for p in inter.parameters.as_array}
ordered_parameters = [parameters['A'],
parameters['B'] ** -1,
parameters['C']]
try:
assert ordered_parameters[1] > 0
except AssertionError as error:
raise ValueError('LAMMPS Buckingham parameter rho (= 1 / B)'
' must be greater than 0') from error
coeff_cmd = (pair_style + ' '
+ ' '.join(str(p) for p in ordered_parameters) + ' '
+ cutoffs)
elif 'lj' in pair_style:
for inter in interactions:
if inter.function.name == 'LennardJones':
parameters = {p.type: convert_unit(p.value)
for p in inter.parameters.as_array}
ordered_parameters = [parameters['epsilon'],
parameters['sigma']]
coeff_cmd = (pair_style + ' '
+ ' '.join(str(p) for p in ordered_parameters) + ' '
+ cutoffs)
elif pair_style in coul_styles:
coeff_cmd = pair_style
else:
raise NotImplementedError('This InteractionFunction has not been'
' implemented in the LAMMPS facade')
coeff_cmds.append(coeff_cmd)
return coeff_cmds
[docs]
def parse_kspace_solver(solver: KSpaceSolver) -> list:
"""
Converts an MDMC ``KSpaceSolver`` for input to LAMMPS ``kspace_style``
Parameters
----------
solver : KSpaceSolver
A ``KSpaceSolver`` to be parsed.
Returns
-------
`list`
Style and parameters for input to LAMMPS ``kspace_style``
Raises
------
NotImplementedError
If ``solver`` type has has not been implemented in the LAMMPS facade.
"""
solver_name = solver.name.lower()
if solver_name not in ('ewald', 'pppm'):
raise NotImplementedError(f'This k-space solver ({solver_name}) has not been implemented'
' in the LAMMPS facade')
return [solver_name, solver.accuracy]
[docs]
def parse_constraint(constraint_algorithm: ConstraintAlgorithm,
bonds: list[Bond] = None, bond_ID_dict: dict = None,
angles: list[BondAngle] = None, angle_ID_dict: dict = None) -> list:
# ID is an acronym
# pylint: disable=invalid-name
"""
Converts an MDMC ``ConstraintAlgorithm`` for input to LAMMPS fix
At least one of bonds and angles must be passed.
Parameters
----------
constraint_algorithm : ConstraintAlgorithm
An object that derives from ``ConstraintAlgorithm`` to be parsed.
bonds : list of Bonds, optional
Constrained ``Bond`` interactions.
bond_ID_dict : dict, optional
`dict` with {``bond``: ``ID pairs``} relating each ``Bond`` object to a
LAMMPS ``ID``.
angles : list of BondAngles, optional
Constrained ``BondAngle`` interactions.
angle_ID_dict : dict, optional
`dict` with {``angle``: ``ID pairs``} relating each ``BondAngle`` object
to a LAMMPS ``ID``.
Returns
-------
`list`
Input parameters for LAMMPS ``fix``, not including the first two terms
(``fix ID``, ``group-ID``). The output list has a maximum length of 7,
where the last four entries are optional but a minimum of two is
required::
[algorithm name, accuracy, max iterations, 'b', bond IDs,
'a', angle IDs]
Raises
------
TypeError
If there is not at least one ``constrained`` ``Interaction`` passed.
NotImplementedError
If ``constraint_algorithm`` not been implemented in the LAMMPS facade.
"""
# Raise error if there is not at least one constrained interaction passed
if not (bonds or angles):
raise TypeError('A LAMMPS constraint fix must have constraints on at'
' least one bond or one bond angle')
lmp_str = []
# Add algorithm name
if constraint_algorithm.name.upper() == 'SHAKE':
lmp_str.append('shake')
elif constraint_algorithm.name.upper() == 'RATTLE':
lmp_str.append('rattle')
else:
raise NotImplementedError('This constraint is not implemented in the'
' LAMMPS facade')
# Add accuracy and max iterations
lmp_str.append(float(constraint_algorithm.accuracy))
lmp_str.append(int(constraint_algorithm.max_iterations))
# Never display the constraint statistics
lmp_str.append(0)
# Add bonds and their LAMMPS IDs and angles and their LAMMPS IDs
if bonds:
if not bond_ID_dict:
bond_ID_dict = {}
lmp_str.append('b')
lmp_str += [bond_ID_dict[bond] for bond in bonds if bond.constrained]
if angles:
if not angle_ID_dict:
angle_ID_dict = {}
lmp_str.append('a')
lmp_str += [angle_ID_dict[angle] for angle in angles
if angle.constrained]
return lmp_str