diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9741c52c..2267b6ea 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -26,7 +26,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.9", "3.10", "3.11", "3.12"] mpi: ["nompi", "mpich", "openmpi"] steps: - name: Checkout code @@ -68,7 +68,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.8", "3.9", "3.10", "3.11"] + python-version: ["3.9", "3.10", "3.11"] hoomd-version: ["2.9.7", "3.11.0", "4.7.0"] include: - python-version: "3.12" @@ -113,7 +113,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.8", "3.9", "3.10", "3.11"] + python-version: ["3.9", "3.10", "3.11"] lammps-version: ["2022.06.23", "2023.08.02"] mpi: ["mpich", "openmpi"] steps: diff --git a/.readthedocs.yml b/.readthedocs.yml index 41ac196e..b1f0c682 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -3,7 +3,7 @@ version: 2 build: os: ubuntu-22.04 tools: - python: "3.8" + python: "3.9" sphinx: configuration: doc/source/conf.py diff --git a/requirements.txt b/requirements.txt index aac4d42a..6578b7f3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,6 @@ freud-analysis>=2 -gsd<3.3.0; python_version < '3.9' -gsd; python_version >= '3.9' -lammpsio>=0.6.1 +gsd +lammpsio>=0.7.0 networkx>=2.5 numpy packaging diff --git a/setup.cfg b/setup.cfg index c5464bd7..5294968d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -31,7 +31,6 @@ classifiers = License :: OSI Approved :: BSD License Topic :: Scientific/Engineering :: Chemistry Topic :: Scientific/Engineering :: Physics - Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 Programming Language :: Python :: 3.10 Programming Language :: Python :: 3.11 @@ -42,12 +41,11 @@ packages = find: package_dir = = src include_package_data = True -python_requires = >=3.8 +python_requires = >=3.9 install_requires = freud-analysis>=2 - gsd<3.3.0; python_version < '3.9' - gsd; python_version >= '3.9' - lammpsio>=0.6.1 + gsd + lammpsio>=0.7.0 networkx>=2.5 numpy packaging diff --git a/src/relentless/model/potential/__init__.py b/src/relentless/model/potential/__init__.py index fbccc505..7231ccf9 100644 --- a/src/relentless/model/potential/__init__.py +++ b/src/relentless/model/potential/__init__.py @@ -16,6 +16,17 @@ PairSpline Yukawa +Bond potentials +=============== + +.. autosummary:: + :toctree: generated/ + + HarmonicBond + FENEWCA + BondSpline + + Developer classes ================= @@ -24,11 +35,15 @@ Potential Parameters + BondedPotential + BondPotential + BondParameters PairPotential PairParameters """ +from .bond import FENEWCA, BondParameters, BondPotential, BondSpline, HarmonicBond from .pair import ( Depletion, LennardJones, @@ -37,4 +52,4 @@ PairSpline, Yukawa, ) -from .potential import Parameters, Potential +from .potential import BondedPotential, Parameters, Potential diff --git a/src/relentless/model/potential/bond.py b/src/relentless/model/potential/bond.py new file mode 100644 index 00000000..763b4b9d --- /dev/null +++ b/src/relentless/model/potential/bond.py @@ -0,0 +1,560 @@ +import numpy + +from relentless import math +from relentless.model import variable + +from . import potential + + +class BondParameters(potential.Parameters): + """Parameters for a bond potential.""" + + pass + + +class BondPotential(potential.BondedPotential): + r"""Abstract base class for a bond potential.""" + + pass + + +class HarmonicBond(BondPotential): + r"""Harmonic bond potential. + + .. math:: + + u(r) = \frac{k}{2} (r - r_0)^2 + + where :math:`r` is the distance between two bonded particles. The parameters + for each type are: + + +-------------+--------------------------------------------------+ + | Parameter | Description | + +=============+==================================================+ + | ``k`` | Spring constant :math:`k`. | + +-------------+--------------------------------------------------+ + | ``r0`` | Minimum-energy length :math:`r_0`. | + +-------------+--------------------------------------------------+ + + Parameters + ---------- + types : tuple[str] + Types. + name : str + Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the + unique integer ID of the potential. + + Attributes + ---------- + coeff : :class:`BondParameters` + Parameters of the potential for each type. + + Examples + -------- + Harmonic Bond:: + + >>> u = relentless.potential.bond.Harmonic(("A",)) + >>> u.coeff["A"].update({'k': 1000, 'r0': 1}) + + """ + + def __init__(self, types, name=None): + super().__init__(keys=types, params=("k", "r0"), name=name) + + def energy(self, type_, r): + """Evaluate bond energy.""" + params = self.coeff.evaluate(type_) + k = params["k"] + r0 = params["r0"] + + r, u, s = self._zeros(r) + + u = 0.5 * k * (r - r0) ** 2 + + if s: + u = u.item() + return u + + def force(self, type_, r): + """Evaluate bond force.""" + params = self.coeff.evaluate(type_) + k = params["k"] + r0 = params["r0"] + + r, f, s = self._zeros(r) + + f = -k * (r - r0) + + if s: + f = f.item() + return f + + def _derivative(self, param, r, k, r0): + r"""Evaluate bond derivative with respect to a variable.""" + r, d, s = self._zeros(r) + + if param == "k": + d = (r - r0) ** 2 / 2 + elif param == "r0": + d = -k * (r - r0) + else: + raise ValueError("Unknown parameter") + + if s: + d = d.item() + return d + + +class FENEWCA(BondPotential): + r"""Finitely extensible nonlinear elastic (FENE) + Weeks-Chandler-Andersen + (WCA) bond interaction. + + .. math:: + + u(r) = - \frac{1}{2} k r_0^2 \ln \left[ 1 - \left( \frac{r}{r_0} + \right)^2 \right] + u_{\rm{WCA}}(r) + + where + + .. math:: + u_{\rm{WCA}}(r) = + \begin{cases} 4 \varepsilon \left[ \left( \frac{\sigma}{r} + \right)^{12} - \left( \frac{\sigma}{r} + \right)^6 \right] + \varepsilon + & r < 2^{1/6}\sigma\\ + 0 & r \ge 2^{1/6}\sigma + \end{cases} + + where :math:`r` is the distance between two bonded particles. The parameters + for each type are: + + +-------------+--------------------------------------------------+ + | Parameter | Description | + +=============+==================================================+ + | ``k`` | Spring constant :math:`k`. | + +-------------+--------------------------------------------------+ + | ``r0`` | Minimum-energy length :math:`r_0`. | + +-------------+--------------------------------------------------+ + | ``epsilon`` | Interaction energy :math:`\varepsilon`. | + +-------------+--------------------------------------------------+ + | ``sigma`` | Interaction length :math:`\sigma`. | + +-------------+--------------------------------------------------+ + + Parameters + ---------- + types : tuple[str] + Types. + name : str + Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the + unique integer ID of the potential. + + Attributes + ---------- + coeff : :class:`BondParameters` + Parameters of the potential for each type. + + Examples + -------- + Kremer-Grest bond (no repulsion included):: + + >>> u = relentless.potential.bond.FENEWCA(("A")) + >>> u.coeff["A"].update({'k': 30, 'r0': 1.5}) + + Kremer-Grest bond (repulsion included):: + + >>> u = relentless.potential.bond.FENEWCA(("A")) + >>> u.coeff["A"].update({'k': 30, 'r0': 1.5, "epsilon":1, "sigma": 1}) + """ + + def __init__(self, types, name=None): + super().__init__(keys=types, params=("k", "r0", "epsilon", "sigma"), name=name) + + def energy(self, type_, r): + """Evaluate bond energy.""" + params = self.coeff.evaluate(type_) + k = params["k"] + r0 = params["r0"] + epsilon = params["epsilon"] + sigma = params["sigma"] + + # initialize arrays + r, u_fene, s = self._zeros(r) + u_wca = u_fene.copy() + + # set flags for FENE potential + fene_flag = numpy.less(r, r0) + + # evaluate FENE potential + u_fene[fene_flag] = -0.5 * k * r0**2 * numpy.log(1 - (r[fene_flag] / r0) ** 2) + u_fene[~fene_flag] = numpy.inf + + # evaluate WCA potential + nonzero_flags = ~numpy.isclose(r, 0) + wca_flags = numpy.logical_and(nonzero_flags, r < 2 ** (1 / 6) * sigma) + r6_inv = numpy.power(sigma / r[wca_flags], 6) + u_wca[wca_flags] = 4.0 * epsilon * (r6_inv**2 - r6_inv) + epsilon + u_wca[~nonzero_flags] = numpy.inf + + if s: + u_fene = u_fene.item() + u_wca = u_wca.item() + + return u_fene + u_wca + + def force(self, type_, r): + """Evaluate bond force.""" + params = self.coeff.evaluate(type_) + k = params["k"] + r0 = params["r0"] + epsilon = params["epsilon"] + sigma = params["sigma"] + + # initialize arrays + r, f_fene, s = self._zeros(r) + f_wca = f_fene.copy() + + # set flags for FENE potential + fene_flag = numpy.less(r, r0) + + # evaluate FENE potential + f_fene[fene_flag] = -(k * r[fene_flag]) / (1 - (r[fene_flag] / r0) ** 2) + f_fene[~fene_flag] = numpy.inf + + # evaluate WCA potential + nonzero_flags = ~numpy.isclose(r, 0) + wca_flags = numpy.logical_and(nonzero_flags, r < 2 ** (1 / 6) * sigma) + rinv = 1.0 / r[wca_flags] + r6_inv = numpy.power(sigma * rinv, 6) + f_wca[wca_flags] = (48.0 * epsilon * rinv) * (r6_inv**2 - 0.5 * r6_inv) + f_wca[~nonzero_flags] = numpy.inf + + if s: + f_fene = f_fene.item() + f_wca = f_wca.item() + + return f_fene + f_wca + + def _derivative(self, param, r, k, r0, epsilon, sigma): + r"""Evaluate bond derivative with respect to a variable.""" + # initialize arrays + r, d, s = self._zeros(r) + + # set flags for FENE potential + fene_flag = numpy.less(r, r0) + + # set flags for WCA potential + nonzero_flags = ~numpy.isclose(r, 0) + wca_flags = numpy.logical_and(nonzero_flags, r < 2 ** (1 / 6) * sigma) + + # set r**6 for WCA potential + r6_inv = numpy.power(sigma / r[wca_flags], 6) + + if param == "k": + d[fene_flag] = 0.5 * r0**2 * numpy.log(1 - (r[fene_flag] / r0) ** 2) + d[~fene_flag] = numpy.inf + elif param == "r0": + d[fene_flag] = (k * r[fene_flag] ** 2) / ( + (1 - (r[fene_flag] / r0) ** 2) * r0 + ) + k * r0 * numpy.log(1 - (r[fene_flag] / r0) ** 2) + d[~fene_flag] = numpy.inf + elif param == "epsilon": + d[wca_flags] = 4 * (r6_inv**2 - r6_inv) + 1 + d[~nonzero_flags] = numpy.inf + elif param == "sigma": + d[wca_flags] = (48.0 * epsilon / sigma) * (r6_inv**2 - 0.5 * r6_inv) + d[~nonzero_flags] = numpy.inf + else: + raise ValueError("Unknown parameter") + if s: + d = d.item() + return d + + +class BondSpline(BondPotential): + """Spline bond potential. + + The bond potential is defined by interpolation through a set of knot points. + The interpolation scheme uses Akima splines. + + Parameters + ---------- + types : tuple[str] + Types. + num_knots : int + Number of knots. + mode : str + Mode for storing the values of the knots in + :class:`~relentless.variable.Variable` that can be optimized. If + ``mode='value'``, the knot amplitudes are manipulated directly. + If ``mode='diff'``, the amplitude of the *last* knot is fixed, and + differences between neighboring knots are manipulated for all other + knots. Defaults to ``'diff'``. + name : str + Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the + unique integer ID of the potential. + + Examples + -------- + The spline potential is setup from a tabulated potential instead of + specifying knot parameters directly:: + + spline = relentless.potential.bond.BondSpline(types=("bondA",), num_knots=3) + spline.from_array(("bondA"),[0,1,2],[4,2,0]) + + However, the knot variables can be iterated over and manipulated directly:: + + for r,k in spline.knots("bondA"): + k.value = 1.0 + + """ + + valid_modes = ("value", "diff") + + def __init__(self, types, num_knots, mode="diff", name=None): + if isinstance(num_knots, int) and num_knots >= 2: + self._num_knots = num_knots + else: + raise ValueError("Number of spline knots must be an integer >= 2.") + + if mode in self.valid_modes: + self._mode = mode + else: + raise ValueError("Invalid mode, choose from: " + ",".join(self.valid_modes)) + + params = [] + for i in range(self.num_knots): + ri, ki = self.knot_params(i) + params.append(ri) + params.append(ki) + super().__init__(keys=types, params=params, name=name) + + @classmethod + def from_json(cls, data, name=None): + u = super().from_json(data, name) + # reset the knot values as variables since they were set as floats + for type in u.coeff: + for i, (r, k) in enumerate(u.knots(type)): + u._set_knot(type, i, r, k) + + return u + + def from_array(self, types, r, u): + r"""Set up the potential from knot points. + + Parameters + ---------- + types : tuple[str] + The type for which to set up the potential. + r : list + Position of each knot. + u : list + Potential energy of each knot. + + Raises + ------ + ValueError + If the number of ``r`` values is not the same as the number of knots. + ValueError + If the number of ``u`` values is not the same as the number of knots. + + """ + # check that r and u have the right shape + if len(r) != self.num_knots: + raise ValueError("r must have the same length as the number of knots") + if len(u) != self.num_knots: + raise ValueError("u must have the same length as the number of knots") + + # convert to r,knot form given the mode + rs = numpy.asarray(r, dtype=float) + ks = numpy.asarray(u, dtype=float) + if self.mode == "diff": + # difference is next knot minus m y knot, + # with last knot fixed at its current value + ks[:-1] -= ks[1:] + + # convert knot positions to differences + drs = numpy.zeros_like(rs) + drs[0] = rs[0] + drs[1:] = rs[1:] - rs[:-1] + + for i in range(self.num_knots): + self._set_knot(types, i, drs[i], ks[i]) + + def to_json(self): + data = super().to_json() + data["num_knots"] = self.num_knots + data["mode"] = self.mode + return data + + def knot_params(self, i): + r"""Get the parameter names for a given knot. + + Parameters + ---------- + i : int + Key for the knot variable. + + Returns + ------- + str + The parameter name of the :math:`r` value. + str + The parameter name of the knot value. + + Raises + ------ + TypeError + If the knot key is not an integer. + + """ + if not isinstance(i, int): + raise TypeError("Knots are keyed by integers") + return f"dr-{i}", f"{self.mode}-{i}" + + def _set_knot(self, types, i, dr, k): + """Set the value of knot variables. + + The meaning of the value of the knot variable is defined by the ``mode``. + This method is mostly meant to coerce the knot variable types. + + Parameters + ---------- + types : tuple[str] + The type for which to set up the potential. + i : int + Index of the knot. + r : float + Relative position of each knot from previous one. + u : float + Value of the knot variable. + + """ + if i > 0 and dr <= 0: + raise ValueError("Knot spacings must be positive") + + dri, ki = self.knot_params(i) + if isinstance(self.coeff[types][dri], variable.IndependentVariable): + self.coeff[types][dri].value = dr + else: + self.coeff[types][dri] = variable.IndependentVariable(dr) + + if isinstance(self.coeff[types][ki], variable.IndependentVariable): + self.coeff[types][ki].value = k + else: + self.coeff[types][ki] = variable.IndependentVariable(k) + + def energy(self, type_, r): + params = self.coeff.evaluate(type_) + r, u, s = self._zeros(r) + u = self._interpolate(params)(r) + if s: + u = u.item() + return u + + def force(self, type_, r): + params = self.coeff.evaluate(type_) + r, f, s = self._zeros(r) + f = -self._interpolate(params).derivative(r, 1) + if s: + f = f.item() + return f + + def _derivative(self, param, r, **params): + r, d, s = self._zeros(r) + h = 0.001 + + if "dr-" in param: + f_low = self._interpolate(params)(r) + knot_p = params[param] + params[param] = knot_p + h + f_high = self._interpolate(params)(r) + params[param] = knot_p + d = (f_high - f_low) / h + elif self.mode + "-" in param: + # perturb knot param value + knot_p = params[param] + params[param] = knot_p + h + f_high = self._interpolate(params)(r) + params[param] = knot_p - h + f_low = self._interpolate(params)(r) + params[param] = knot_p + d = (f_high - f_low) / (2 * h) + else: + raise ValueError("Parameter cannot be differentiated") + + if s: + d = d.item() + return d + + def _interpolate(self, params): + """Interpolate the knot points into a spline potential. + + Parameters + ---------- + params : dict + The knot parameters + + Returns + ------- + :class:`~relentless.math.Interpolator` + The interpolated spline potential. + + """ + dr = numpy.zeros(self.num_knots) + u = numpy.zeros(self.num_knots) + for i in range(self.num_knots): + dri, ki = self.knot_params(i) + dr[i] = params[dri] + u[i] = params[ki] + + if numpy.any(dr[1:] <= 0): + raise ValueError("Knot differences must be positive") + + # reconstruct r from differences + r = numpy.cumsum(dr) + + # reconstruct the energies from differences, starting from the end + if self.mode == "diff": + u = numpy.flip(numpy.cumsum(numpy.flip(u))) + + return math.AkimaSpline(x=r, y=u) + + @property + def num_knots(self): + """int: Number of knots.""" + return self._num_knots + + @property + def mode(self): + """str: Spline construction mode.""" + return self._mode + + def knots(self, types): + r"""Generator for knot points. + + Parameters + ---------- + types : tuple[str] + The types for which to retrieve the knot points. + + Yields + ------ + :class:`~relentless.variable.Variable` + The next :math:`dr` variable in the parameters. + :class:`~relentless.variable.Variable` + The next knot variable in the parameters. + + """ + for i in range(self.num_knots): + dri, ki = self.knot_params(i) + yield self.coeff[types][dri], self.coeff[types][ki] + + @property + def design_variables(self): + """tuple: Designable variables of the spline.""" + dvars = [] + for types in self.coeff: + for i, (dr, k) in enumerate(self.knots(types)): + if i != self.num_knots - 1: + dvars.append(k) + return tuple(dvars) diff --git a/src/relentless/model/potential/potential.py b/src/relentless/model/potential/potential.py index dc4a5026..1cc7b9e9 100644 --- a/src/relentless/model/potential/potential.py +++ b/src/relentless/model/potential/potential.py @@ -375,3 +375,118 @@ def save(self, filename): data = self.to_json() with open(filename, "w") as f: json.dump(data, f, sort_keys=True, indent=4) + + +class BondedPotential(Potential): + r"""Abstract base class for bonded interaction potential. + + This class can be extended to evaluate the energy, force, and parameter + derivatives of a bonded potential with a given functional form. + :meth:`energy` specifies the potential energy :math:`u_0(r)`, :meth:`force` + specifies the force :math:`f_0(r) = -\partial u_0/\partial r`, and + :meth:`_derivative` specifies the derivative + :math:`u_{0,\lambda} = \partial u_0/\partial \lambda` with respect to + parameter :math:`\lambda`. + """ + + def derivative(self, type_, var, r): + r"""Evaluate bond derivative with respect to a variable. + + The derivative is evaluated using the :meth:`_derivative` function for all + :math:`u_{0,\lambda}(r)`. + + The derivative will be carried out with respect to ``var`` for all + :class:`~relentless.variable.Variable` parameters. The appropriate chain + rules are handled automatically. If the potential does not depend on + ``var``, the derivative will be zero by definition. + + Parameters + ---------- + _type : tuple[str] + The type for which to calculate the derivative. + var : :class:`~relentless.variable.Variable` + The variable with respect to which the derivative is calculated. + r : float or list + The bond distance(s) at which to evaluate the derivative. + + Returns + ------- + float or numpy.ndarray + The bond derivative evaluated at ``r``. The return type is consistent + with ``r``. + + Raises + ------ + ValueError + If any value in ``r`` is negative. + TypeError + If the parameter with respect to which to take the derivative + is not a :class:`~relentless.variable.Variable`. + ValueError + If the potential is shifted without setting ``rmax``. + + """ + params = self.coeff.evaluate(type_) + r, deriv, scalar_r = self._zeros(r) + if any(r < 0): + raise ValueError("r cannot be negative") + if not isinstance(var, variable.Variable): + raise TypeError( + "Parameter with respect to which to take the derivative" + " must be a Variable." + ) + + flags = numpy.ones(r.shape[0], dtype=bool) + + for p in self.coeff.params: + # try to take chain rule w.r.t. variable first + p_obj = self.coeff[type_][p] + if isinstance(p_obj, variable.DependentVariable): + dp_dvar = p_obj.derivative(var) + elif isinstance(p_obj, variable.IndependentVariable) and var is p_obj: + dp_dvar = 1.0 + else: + dp_dvar = 0.0 + + # skip when dp_dvar is exactly zero, since this does not contribute + if dp_dvar == 0.0: + continue + + # now take the parameter derivative + below = numpy.zeros(r.shape[0], dtype=bool) + above = numpy.zeros(r.shape[0], dtype=bool) + + flags = numpy.logical_and(~below, ~above) + deriv[flags] += self._derivative(p, r[flags], **params) * dp_dvar + + # coerce derivative back into shape of the input + if scalar_r: + deriv = deriv.item() + return deriv + + @abc.abstractmethod + def _derivative(self, param, r, **params): + """Implementation of the parameter derivative function. + + This abstract method defines the interface for computing the parameter + derivative of a bond interaction. ``**params`` will include all the + parameters from :class:`PairParameters`. The derivative should be + consistent with :meth:`_energy`. + + Parameters + ---------- + param : str + Name of the parameter. + r : float or list + The bond distance(s) at which to evaluate the derivative. + **params : kwargs + Named parameters of the potential. + + Returns + ------- + float or numpy.ndarray + The bond derivative evaluated at ``r``. The return type is consistent + with ``r``. + + """ + pass diff --git a/src/relentless/simulate/__init__.py b/src/relentless/simulate/__init__.py index dc40e3f7..b89b4f50 100644 --- a/src/relentless/simulate/__init__.py +++ b/src/relentless/simulate/__init__.py @@ -122,6 +122,7 @@ Potentials PotentialTabulator + BondPotentialTabulator PairPotentialTabulator Developer classes @@ -157,6 +158,7 @@ ) from .simulate import ( AnalysisOperation, + BondPotentialTabulator, PairPotentialTabulator, Potentials, PotentialTabulator, diff --git a/src/relentless/simulate/analyze.py b/src/relentless/simulate/analyze.py index f9538f22..b315f6ae 100644 --- a/src/relentless/simulate/analyze.py +++ b/src/relentless/simulate/analyze.py @@ -73,10 +73,15 @@ def _get_rdf_params(self, sim): else: rdf_every = self.every + if "exclude" in self.rdf: + rdf_exclude = bool(self.rdf["exclude"]) + else: + rdf_exclude = False rdf_params = { "bins": self.rdf["num"], "stop": self.rdf["stop"], "every": rdf_every, + "exclude": rdf_exclude, } else: rdf_params = None diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py index 33175336..018fbf12 100644 --- a/src/relentless/simulate/hoomd.py +++ b/src/relentless/simulate/hoomd.py @@ -66,7 +66,15 @@ def _call_v3(self, sim): sim.masses = self._get_masses_from_snapshot(sim, snap) self._assert_dimension_safe(sim, snap) # create the potentials, defer attaching until later - neighbor_list = hoomd.md.nlist.Tree(buffer=sim.potentials.pair.neighbor_buffer) + exclusion = sim.potentials.pair.exclusions + if exclusion is None: + exclusion = () + elif "1-2" in exclusion: + exclusion = ["bond" if ex == "1-2" else ex for ex in exclusion] + neighbor_list = hoomd.md.nlist.Tree( + buffer=sim.potentials.pair.neighbor_buffer, + exclusions=exclusion, + ) pair_potential = hoomd.md.pair.Table(nlist=neighbor_list) r, u, f = sim.potentials.pair.pairwise_energy_and_force( sim.types, tight=True, minimum_num=2 @@ -79,7 +87,22 @@ def _call_v3(self, sim): ) pair_potential.r_cut[(i, j)] = r[-1] sim[self]["_potentials"] = [pair_potential] - sim[self]["_potentials_rmax"] = r[-1] + + sim[self]["_bonds"] = self._get_bonds_from_snapshot(sim, snap) + if sim.potentials.bond is not None: + sim.bond_types = sim["engine"]["_hoomd"].state.bond_types + bond_potential = hoomd.md.bond.Table(width=sim.potentials.bond.num) + + for i in sim.bond_types: + r, u, f = ( + sim.potentials.bond.linear_space, + sim.potentials.bond.energy(key=i), + sim.potentials.bond.force(key=i), + ) + if numpy.any(numpy.isinf(u)): + raise ValueError("Bond potential/force is infinite at evaluated r") + bond_potential.params[i] = dict(r_min=r[0], r_max=r[-1], U=u[:], F=f[:]) + sim[self]["_potentials"].append(bond_potential) def _call_v2(self, sim): # initialize @@ -122,7 +145,6 @@ def _table_eval(r_i, rmin, rmax, **coeff): rmax=r[-1], coeff=dict(r=r, u=u[i, j], f=f[i, j]), ) - sim[self]["_potentials_rmax"] = r[-1] def _initialize_v3(self, sim): raise NotImplementedError( @@ -152,6 +174,12 @@ def _get_masses_from_snapshot(self, sim, snap): masses.update(masses_) return masses + def _get_bonds_from_snapshot(self, sim, snap): + if mpi.world.rank_is_root: + return snap.bonds.group + else: + return None + def _assert_dimension_safe(self, sim, snap): if sim.dimension == 3: dim_safe = True @@ -987,6 +1015,8 @@ def _pre_run_v3(self, sim, sim_op): system=sim["engine"]["_hoomd"].state, rdf_params=self._get_rdf_params(sim), constraints=self._get_constrained_quantities(sim, sim_op), + exclusions=sim.potentials.pair.exclusions, + bonds=sim[sim.initializer]["_bonds"], ) sim[self]["_hoomd_thermo_callback"] = hoomd.write.CustomWriter( trigger=self.every, @@ -1178,6 +1208,8 @@ def __init__( system, rdf_params=None, constraints=None, + exclusions=None, + bonds=None, ): if dimension not in (2, 3): raise ValueError("Only 2 or 3 dimensions supported") @@ -1188,6 +1220,8 @@ def __init__( self.system = system self.rdf_params = rdf_params self.constraints = constraints if constraints is not None else {} + self.exclusion = exclusions + self.bonds = bonds # this method handles all the initialization self.reset() @@ -1282,6 +1316,24 @@ def act(self, timestep): ), ).toNeighborList() + # filter bonds from the neighbor list if they are present + # bond exclusions apply regardless of order, so + # consider both (i,j) and (j,i) permutations + if ( + self.rdf_params["exclude"] + and snap.bonds.N != 0 + and len(neighbors[:]) > 0 + ): + bonds = numpy.vstack( + [self.bonds, numpy.flip(self.bonds, axis=1)], + ) + + bond_exclusion_filter = EnsembleAverage._cantor_pairing( + self, bonds, neighbors + ) + + neighbors.filter(bond_exclusion_filter) + for i in self.types: for j in self.types: filter_ij = numpy.logical_and( diff --git a/src/relentless/simulate/lammps.py b/src/relentless/simulate/lammps.py index 703f4cd3..8253758f 100644 --- a/src/relentless/simulate/lammps.py +++ b/src/relentless/simulate/lammps.py @@ -181,8 +181,19 @@ def _call_commands(self, sim): sim.masses = collections.FixedKeyDict(sim.types) masses = {} dim_safe = None + has_bonds = None + has_angles = None + has_dihedrals = None + has_impropers = None + bond_type_label = None if mpi.world.rank_is_root: snap = lammpsio.DataFile(sim[self]["_datafile"]).read() + has_bonds = snap.has_bonds() + if has_bonds: + bond_type_label = sim[self]["_bonds"].type_label + has_angles = snap.has_angles() + has_dihedrals = snap.has_dihedrals() + has_impropers = snap.has_impropers() for i in sim.types: mi = snap.mass[snap.typeid == sim["engine"]["types"][i]] if len(mi) == 0: @@ -205,6 +216,12 @@ def _call_commands(self, sim): if not dim_safe: raise ValueError("Simulation initialized inconsistent with dimension") + has_bonds = mpi.world.bcast(has_bonds) + has_angles = mpi.world.bcast(has_angles) + has_dihedrals = mpi.world.bcast(has_dihedrals) + has_impropers = mpi.world.bcast(has_impropers) + bond_type_label = mpi.world.bcast(bond_type_label) + # attach the potentials if sim.potentials.pair.start == 0: raise ValueError("LAMMPS requires start > 0 for pair potentials") @@ -212,7 +229,19 @@ def _call_commands(self, sim): sim.types, x=sim.potentials.pair.squared_space, tight=True, minimum_num=2 ) Nr = len(r) - sim[self]["_potentials_rmax"] = r[-1] + + if has_bonds: + uB = collections.FixedKeyDict(bond_type_label.types) + fB = collections.FixedKeyDict(bond_type_label.types) + for i in bond_type_label.types: + rB, uB[i], fB[i] = ( + sim.potentials.bond.linear_space, + sim.potentials.bond.energy(key=i), + sim.potentials.bond.force(key=i), + ) + if numpy.any(numpy.isinf(uB[i])): + raise ValueError("Bond potential/force is infinite at evaluated r") + NrB = len(rB) def pair_map(sim, pair): # Map lammps type indexes as a pair, lowest type first @@ -257,6 +286,27 @@ def pair_map(sim, pair): fw.write( "{idx} {r} {u} {f}\n".format(idx=idx, r=r_, u=u_, f=f_) ) + + # write bond potentials into the file + if has_bonds: + fw.write("# LAMMPS tabulated bond potentials\n") + for i in bond_type_label.types: + fw.write("# bond BOND_TABLE_{}\n\n".format(i)) + fw.write("BOND_TABLE_{}\n".format(i)) + fw.write( + "N {N} FP {f_low} {f_high} \n\n".format( + N=NrB, f_low=fB[i][0], f_high=fB[i][-1] + ) + ) + for idx, (rB_, uB_, fB_) in enumerate( + zip(rB, uB[i], fB[i]), start=1 + ): + fw.write( + "{id} {rB} {uB} {fB}\n".format( + id=idx, rB=rB_, uB=uB_, fB=fB_ + ) + ) + else: file_ = None file_ = mpi.world.bcast(file_) @@ -268,6 +318,17 @@ def pair_map(sim, pair): ] cmds += ["pair_style table linear {N}".format(N=Nr)] + # set exclusions if snap has topology + if has_bonds or has_angles or has_dihedrals: + excl_12, excl_13, excl_14 = 1.0, 1.0, 1.0 + if sim.potentials.pair.exclusions is not None: + if "1-2" in sim.potentials.pair.exclusions: + excl_12 = 0.0 + cmds += [f"special_bonds lj/coul {excl_12} {excl_13} {excl_14}"] + + if has_bonds: + cmds += ["bond_style table linear {Nb}".format(Nb=NrB)] + for i, j in sim.pairs: # get lammps type indexes, lowest type first id_i, id_j = pair_map(sim, (i, j)) @@ -277,6 +338,15 @@ def pair_map(sim, pair): ) ] + if has_bonds: + for id in bond_type_label.typeid: + cmds += [ + ("bond_coeff {typeid} {filename}" " BOND_TABLE_{id}").format( + typeid=id, + filename=file_, + id=bond_type_label.__getitem__(id), + ) + ] return cmds @abc.abstractmethod @@ -312,6 +382,19 @@ def _convert_to_data_file(self, sim): snap, type_map = lammpsio.Snapshot.from_hoomd_gsd(frame) type_map = {v: k for k, v in type_map.items()} + # check atom style is compatible with topology data if present + if ( + snap.has_bonds() + or snap.has_angles() + or snap.has_dihedrals() + or snap.has_impropers() + ) and sim["engine"]["atom_style"] == "atomic": + raise ValueError( + "Atomic atom style is not compatible with topology data." + ) + # store topology data + if snap.has_bonds(): + sim[self]["_bonds"] = snap.bonds # figure out dimensions dimension = self.dimension if dimension is None: @@ -330,6 +413,7 @@ def _convert_to_data_file(self, sim): data_filename, snap, sim["engine"]["atom_style"] ) else: + sim[self]["_bonds"] = None data_filename = None dimension = None type_map = None @@ -348,6 +432,20 @@ def _convert_to_data_file(self, sim): ).read() typeids = numpy.unique(snap.typeid) type_map = {str(typeid): typeid for typeid in typeids} + # store topology data + if snap.has_bonds(): + sim[self]["_bonds"] = snap.bonds + unique_bond_types = set() + for pot in sim.potentials.bond.potentials: + for i in pot.coeff.types: + unique_bond_types.add(i) + + bond_type_map = { + i + 1: t for i, t in enumerate(unique_bond_types) + } + sim[self]["_bonds"].type_label = lammpsio.LabelMap( + bond_type_map + ) else: type_map = None type_map = mpi.world.bcast(type_map) @@ -1176,6 +1274,31 @@ def process(self, sim, sim_op): exclude_ii=True, ), ).toNeighborList() + # filter bonds from the neighbor list if they are present + # bond exclusions apply regardless of order, so + # consider both (i,j) and (j,i) permutations + if ( + sim[self]["_rdf_params"]["exclude"] + and sim[sim.initializer]["_bonds"].N != 0 + and len(neighbors[:]) > 0 + ): + bonds = numpy.vstack( + [ + sim[sim.initializer]["_bonds"].members, + numpy.flip( + sim[sim.initializer]["_bonds"].members, axis=1 + ), + ], + ) + # Zero index bonds + bonds -= 1 + + bond_exclusion_filter = EnsembleAverage._cantor_pairing( + self, bonds, neighbors + ) + + neighbors.filter(bond_exclusion_filter) + for i in sim.types: _rdf_density[i] += N[i] / box.volume _rdf_num_origins[i] += N[i] @@ -1511,6 +1634,7 @@ def __init__( quiet=True, types=None, executable=None, + atom_style="atomic", ): # test executable if it is specified if executable is not None: @@ -1573,6 +1697,7 @@ def __init__( super().__init__(initializer, operations) self.quiet = quiet self.types = types + self.atom_style = atom_style def _post_run(self, sim): # force all the lammps commands to execute, since the operations did @@ -1627,7 +1752,7 @@ def _initialize_engine(self, sim): sim["engine"]["types"] = self.types sim["engine"]["units"] = "lj" - sim["engine"]["atom_style"] = "atomic" + sim["engine"]["atom_style"] = self.atom_style # initialize _InitializeFromFile = InitializeFromFile diff --git a/src/relentless/simulate/simulate.py b/src/relentless/simulate/simulate.py index bd33850f..a28d558c 100644 --- a/src/relentless/simulate/simulate.py +++ b/src/relentless/simulate/simulate.py @@ -78,6 +78,35 @@ def post_run(self, sim, sim_op): def process(self, sim, sim_op): pass + def _cantor_pairing(self, connection, neighbor): + """Compute the list intersect using the Cantor pairing function to + filter neighborlist. + + https://en.wikipedia.org/wiki/Pairing_function + + Parameters + ---------- + connection : numpy.ndarray + (N, 2) array of typeids to be filtered from neighborlist. + neighbor : numpy.ndarray + Neighborlist to be filtered. + + Returns + ------- + numpy.ndarray + Boolean array of length of neighborlist, True if neighbor is not in + connection. + """ + + pi_connection = (connection[:, 0] + connection[:, 1]) * ( + connection[:, 0] + connection[:, 1] + 1 + ) / 2 + connection[:, 1] + pi_neighbor = (neighbor[:, 0] + neighbor[:, 1]) * ( + neighbor[:, 0] + neighbor[:, 1] + 1 + ) / 2 + neighbor[:, 1] + + return ~numpy.isin(pi_neighbor, pi_connection) + class DelegatedInitializationOperation(InitializationOperation): def __call__(self, sim): @@ -350,6 +379,7 @@ class Potentials: def __init__(self, kB=1.0): self.kB = kB self.pair = None + self.bond = None @property def pair(self): @@ -362,6 +392,17 @@ def pair(self, val): raise TypeError("Pair potential must be tabulated") self._pair = val + @property + def bond(self): + """:class:`BondPotentialTabulator`: Bond potentials.""" + return self._bond + + @bond.setter + def bond(self, val): + if val is not None and not isinstance(val, BondPotentialTabulator): + raise TypeError("Bond potential must be tabulated") + self._bond = val + class PotentialTabulator: r"""Tabulator for an interaction potential. @@ -566,12 +607,18 @@ class PairPotentialTabulator(PotentialTabulator): potential. neighbor_buffer : float Buffer radius used in computing the neighbor list. + exclusions : list + The neighborlist nominally includes all pairs within ``rmax`` of each other. + This option allows for exclusions of pairs that should not be included in the + neighbor list. The string should be formatted as a tuple of strings. Allowed + values are '1-2', '1-3', and '1-4'. """ - def __init__(self, potentials, start, stop, num, neighbor_buffer): + def __init__(self, potentials, start, stop, num, neighbor_buffer, exclusions=None): super().__init__(potentials, start, stop, num) self.neighbor_buffer = neighbor_buffer + self.exclusions = exclusions @property def neighbor_buffer(self): @@ -583,6 +630,26 @@ def neighbor_buffer(self): def neighbor_buffer(self, val): self._neighbor_buffer = val + @property + def exclusions(self): + r"""tuple: The pairs to exclude from the neighbor list. + + Exclusions are formatted as a tuple of strings. Allowed values are: + + - ``'1-2'``: Exclude pairs separated by one bond. + - ``'1-3'``: Exclude pairs separated by two bonds. + - ``'1-4'``: Exclude pairs separated by three bonds. + """ + return self._exclusions + + @exclusions.setter + def exclusions(self, val): + if val is not None: + allowed = ["1-2", "1-3", "1-4"] + if not all([ex in allowed for ex in val]): + raise ValueError("Exclusions must be '1-2', '1-3', or '1-4'") + self._exclusions = val + def energy(self, pair, x=None): """Evaluates and accumulates energy for all potentials. @@ -746,3 +813,27 @@ def pairwise_energy_and_force(self, types, x=None, tight=False, minimum_num=2): f[pair] = f[pair].item() return x, u, f + + +class BondPotentialTabulator(PotentialTabulator): + r"""Tabulate one or more bond potentials. + + Enables evaluation of energy, force, and derivative at different + positional values (i.e. ``r``). + + Parameters + ---------- + potentials : :class:`~relentless.potential.bond.BondPotential` or array_like + The bond potential(s) to be tabulated. If array_like, all elements must + be :class:`~relentless.potential.bond.BondPotential`\s. + start : float + The minimum value of ``r`` at which to tabulate. + stop : float + The maximum value of ``r`` at which to tabulate. + num : int + The number of points (value of ``r``) at which to tabulate and evaluate the + potential. + + """ + + pass diff --git a/tests/model/potential/test_bond.py b/tests/model/potential/test_bond.py new file mode 100644 index 00000000..c6280787 --- /dev/null +++ b/tests/model/potential/test_bond.py @@ -0,0 +1,617 @@ +"""Unit tests for bond module.""" + +import tempfile +import unittest + +import numpy + +import relentless + + +class LinPotBond(relentless.model.potential.bond.BondPotential): + """Linear bond potential function""" + + def __init__(self, types, params): + super().__init__(types, params) + + def to_json(self): + data = super().to_json() + data["params"] = self.coeff.params + return data + + def energy(self, types, r): + m = self.coeff[types]["m"] + + r, u, s = self._zeros(r) + u[:] = m * r + if s: + u = u.item() + return u + + def force(self, types, r): + m = self.coeff[types]["m"] + + r, f, s = self._zeros(r) + f[:] = -m + if s: + f = f.item() + return f + + def _derivative(self, param, r, **params): + r, d, s = self._zeros(r) + if param == "m": + d[:] = r + if s: + d = d.item() + return d + + +class TwoVarPot(relentless.model.potential.bond.BondPotential): + """Mock potential function""" + + def __init__(self, types, params): + super().__init__(types, params) + + @classmethod + def from_json(cls, data): + raise NotImplementedError() + + def energy(self, r, x, y, **params): + pass + + def force(self, r, x, y, **params): + pass + + def _derivative(self, param, r, **params): + # not real derivative, just used to test functionality + r, d, s = self._zeros(r) + if param == "x": + d[:] = 2 * r + elif param == "y": + d[:] = 3 * r + if s: + d = d.item() + return d + + +class test_BondPotential(unittest.TestCase): + """Unit tests for relentless.model.bond.BondPotential""" + + def test_init(self): + """Test creation from data""" + # test creation with only m + p = LinPotBond(types=("1",), params=("m",)) + p.coeff["1"]["m"] = 3.5 + + coeff = relentless.model.potential.bond.BondParameters( + types=("1",), params=("m") + ) + coeff["1"]["m"] = 3.5 + + self.assertCountEqual(p.coeff.types, coeff.types) + self.assertCountEqual(p.coeff.params, coeff.params) + self.assertEqual(p.coeff.evaluate(("1")), coeff.evaluate(("1"))) + + def test_energy(self): + """Test energy method""" + p = LinPotBond(types=("1",), params=("m",)) + p.coeff["1"]["m"] = 2.0 + + # test with no cutoffs + u = p.energy(types=("1"), r=0.5) + self.assertAlmostEqual(u, 1.0) + u = p.energy(types=("1"), r=[0.25, 0.75]) + numpy.testing.assert_allclose(u, [0.5, 1.5]) + + def test_force(self): + """Test force method""" + p = LinPotBond(types=("1",), params=("m",)) + p.coeff["1"]["m"] = 2.0 + + # test with no cutoffs + f = p.force(types=("1"), r=0.5) + self.assertAlmostEqual(f, -2.0) + f = p.force(types=("1"), r=[0.25, 0.75]) + numpy.testing.assert_allclose(f, [-2.0, -2.0]) + + def test_derivative_values(self): + """Test derivative method with different param values""" + p = LinPotBond(types=("1",), params=("m",)) + x = relentless.model.IndependentVariable(value=2.0) + p.coeff["1"]["m"] = x + + # test with no cutoffs + d = p.derivative(type_=("1"), var=x, r=0.5) + self.assertAlmostEqual(d, 0.5) + d = p.derivative(type_=("1"), var=x, r=[0.25, 0.75]) + numpy.testing.assert_allclose(d, [0.25, 0.75]) + + def test_derivative_types(self): + """Test derivative method with different param types.""" + q = LinPotBond(types=("1",), params=("m",)) + x = relentless.model.IndependentVariable(value=4.0) + y = relentless.model.IndependentVariable(value=64.0) + z = relentless.model.GeometricMean(x, y) + q.coeff["1"]["m"] = z + + # test with respect to dependent variable parameter + d = q.derivative(type_=("1"), var=z, r=2.0) + self.assertAlmostEqual(d, 2.0) + + # test with respect to independent variable on which parameter is dependent + d = q.derivative(type_=("1"), var=x, r=1.5) + self.assertAlmostEqual(d, 3.0) + d = q.derivative(type_=("1"), var=y, r=4.0) + self.assertAlmostEqual(d, 0.5) + + # test invalid derivative w.r.t. scalar + a = 2.5 + q.coeff["1"]["m"] = a + with self.assertRaises(TypeError): + d = q.derivative(type_=("1"), var=a, r=2.0) + + # test with respect to independent variable which is + # related to a SameAs variable + r = TwoVarPot(types=("1",), params=("x", "y")) + + r.coeff["1"]["x"] = x + r.coeff["1"]["y"] = relentless.model.variable.SameAs(x) + d = r.derivative(type_=("1"), var=x, r=4.0) + self.assertAlmostEqual(d, 20.0) + + r.coeff["1"]["y"] = x + r.coeff["1"]["x"] = relentless.model.variable.SameAs(x) + d = r.derivative(type_=("1"), var=x, r=4.0) + self.assertAlmostEqual(d, 20.0) + + def test_iteration(self): + """Test iteration on typesPotential object""" + p = LinPotBond(types=("1",), params=("m",)) + for types in p.coeff: + p.coeff[types]["m"] = 2.0 + + self.assertEqual(dict(p.coeff["1"]), {"m": 2.0}) + + def test_json(self): + """Test saving to file""" + p = LinPotBond(types=("1",), params=("m")) + p.coeff["1"]["m"] = 2.0 + + data = p.to_json() + self.assertEqual(data["id"], p.id) + self.assertEqual(data["name"], p.name) + + p2 = LinPotBond.from_json(data) + self.assertEqual(p2.coeff["1"]["m"], 2.0) + + def test_save(self): + """Test saving to file""" + temp = tempfile.NamedTemporaryFile() + p = LinPotBond(types=("1",), params=("m")) + p.coeff["1"]["m"] = 2.0 + p.save(temp.name) + + p2 = LinPotBond.from_file(temp.name) + self.assertEqual(p2.coeff["1"]["m"], 2.0) + + temp.close() + + +class test_HarmonicBond(unittest.TestCase): + """Unit tests for relentless.model.potential.HarmonicBond""" + + def test_init(self): + """Test creation from data""" + harmonic_bond = relentless.model.potential.HarmonicBond(types=("1",)) + coeff = relentless.model.potential.BondParameters( + types=("1",), + params=( + "k", + "r0", + ), + ) + self.assertCountEqual(harmonic_bond.coeff.types, coeff.types) + self.assertCountEqual(harmonic_bond.coeff.params, coeff.params) + + def test_energy(self): + """Test _energy method""" + harmonic_bond = relentless.model.potential.HarmonicBond(types=("1",)) + harmonic_bond.coeff["1"].update(k=1000, r0=1.0) + # test scalar r + r_input = 0.5 + u_actual = 125.0 + u = harmonic_bond.energy(type_=("1"), r=r_input) + self.assertAlmostEqual(u, u_actual) + + # test array r + r_input = numpy.array([0.0, 0.5, 1.0]) + u_actual = numpy.array([500.0, 125.0, 0.0]) + u = harmonic_bond.energy(type_=("1"), r=r_input) + numpy.testing.assert_allclose(u, u_actual) + + def test_force(self): + """Test _force method""" + harmonic_bond = relentless.model.potential.HarmonicBond(types=("1",)) + harmonic_bond.coeff["1"].update(k=1000, r0=1.0) + + # test scalar r + r_input = 0.5 + f_actual = 500 + f = harmonic_bond.force(type_=("1"), r=r_input) + self.assertAlmostEqual(f, f_actual) + + # test array r + r_input = numpy.array([0.0, 0.5, 1.0]) + f_actual = numpy.array([1000, 500, 0]) + f = harmonic_bond.force(type_=("1"), r=r_input) + numpy.testing.assert_allclose(f, f_actual) + + def test_derivative(self): + """Test _derivative method""" + harmonic_bond = relentless.model.potential.HarmonicBond(types=("1",)) + + # w.r.t. k + # test scalar r + r_input = 0.5 + d_actual = 0.125 + d = harmonic_bond._derivative(param="k", r=r_input, k=1000, r0=1.0) + self.assertAlmostEqual(d, d_actual) + + # test array r + r_input = numpy.array([0, 1, 1.5]) + d_actual = numpy.array([0.50, 0.0, 0.125]) + d = harmonic_bond._derivative(param="k", r=r_input, k=1000, r0=1.0) + numpy.testing.assert_allclose(d, d_actual) + + # w.r.t. r0 + # test scalar r + r_input = 0.5 + d_actual = 500 + d = harmonic_bond._derivative(param="r0", r=r_input, k=1000.0, r0=1.0) + self.assertAlmostEqual(d, d_actual) + + # test array r + r_input = numpy.array([0, 1, 1.5]) + d_actual = numpy.array([1000, 0, -500]) + d = harmonic_bond._derivative(param="r0", r=r_input, k=1000.0, r0=1.0) + numpy.testing.assert_allclose(d, d_actual) + + # test invalid param + with self.assertRaises(ValueError): + harmonic_bond._derivative(param="ro", r=r_input, k=1000.0, r0=1.0) + + +class test_FENEWCA(unittest.TestCase): + """Unit tests for relentless.model.potential.FENEWCA""" + + def test_init(self): + """Test creation from data""" + FENEWCA = relentless.model.potential.FENEWCA(types=("1",)) + coeff = relentless.model.potential.BondParameters( + types=("1",), + params=( + "k", + "r0", + "epsilon", + "sigma", + ), + ) + self.assertCountEqual(FENEWCA.coeff.types, coeff.types) + self.assertCountEqual(FENEWCA.coeff.params, coeff.params) + + def test_energy(self): + """Test _energy method""" + FENEWCA = relentless.model.potential.FENEWCA(types=("1",)) + FENEWCA.coeff["1"].update(k=30, r0=1.5, epsilon=1.0, sigma=1.0) + + # test scalar r + r_input = 0.95 + u_actual = 20.2638974009 + u = FENEWCA.energy(type_=("1"), r=r_input) + self.assertAlmostEqual(u, u_actual) + + # test array r + r_input = numpy.array([0, 0.9, 1.2, 2]) + u_actual = numpy.array([numpy.inf, 22.698308667, 34.4807296042, numpy.inf]) + u = FENEWCA.energy(type_=("1"), r=r_input) + numpy.testing.assert_allclose(u, u_actual) + + def test_force(self): + """Test _force method""" + FENEWCA = relentless.model.potential.FENEWCA(types=("1",)) + FENEWCA.coeff["1"].update(k=30, r0=1.5, epsilon=1.0, sigma=1.0) + + # test scalar r + r_input = 0.95 + f_actual = 11.549426778 + f = FENEWCA.force(type_=("1"), r=r_input) + numpy.testing.assert_allclose(f, f_actual) + + # test array r + r_input = numpy.array([0, 0.9, 1.2, 2]) + f_actual = numpy.array([numpy.inf, 96.4721239943, -100, numpy.inf]) + f = FENEWCA.force(type_=("1"), r=r_input) + numpy.testing.assert_allclose(f, f_actual) + + def test_derivative(self): + """Test _derivative method""" + FENEWCA = relentless.model.potential.FENEWCA(types=("1",)) + + # w.r.t. k + # test scalar r + r_input = 0.95 + d_actual = -0.576764091467 + d = FENEWCA._derivative(param="k", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1) + self.assertAlmostEqual(d, d_actual) + + # test array r + r_input = numpy.array([0, 1, 1.6]) + d_actual = numpy.array([0, -0.661259998015, numpy.inf]) + d = FENEWCA._derivative(param="k", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1) + numpy.testing.assert_allclose(d, d_actual) + + # w.r.t. r0 + # test scalar r + r_input = 0.95 + d_actual = 7.06858290903 + d = FENEWCA._derivative(param="r0", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1) + self.assertAlmostEqual(d, d_actual) + + # test array r + r_input = numpy.array([0, 1, 1.6]) + d_actual = numpy.array([0, 9.5496000794, numpy.inf]) + d = FENEWCA._derivative(param="r0", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1) + numpy.testing.assert_allclose(d, d_actual) + + # w.r.t. epsilon + # test scalar r + r_input = 0.95 + d_actual = 2.96097465689 + d = FENEWCA._derivative( + param="epsilon", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1 + ) + self.assertAlmostEqual(d, d_actual) + + # test array r + r_input = numpy.array([0, 1.0, 1.1, 2.0]) + d_actual = numpy.array([numpy.inf, 1.0, 0.0166275506263, 0]) + d = FENEWCA._derivative( + param="epsilon", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1 + ) + + # test array r + r_input = numpy.array([0, 1, 1.6]) + d_actual = numpy.array([0, 9.5496000794, numpy.inf]) + d = FENEWCA._derivative(param="r0", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1) + numpy.testing.assert_allclose(d, d_actual) + + # w.r.t. sigma + # test scalar r + r_input = 0.95 + d_actual = 56.1806752906 + d = FENEWCA._derivative( + param="sigma", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1 + ) + self.assertAlmostEqual(d, d_actual) + + # test array r + r_input = numpy.array([0, 1.0, 1.1, 2.0]) + d_actual = numpy.array([numpy.inf, 24, 1.74690492881, 0]) + d = FENEWCA._derivative( + param="sigma", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1 + ) + numpy.testing.assert_allclose(d, d_actual) + + # test invalid param + with self.assertRaises(ValueError): + FENEWCA._derivative( + param="sgima", r=r_input, k=30, r0=1.5, epsilon=1, sigma=1 + ) + + def test_json(self): + FENEWCA = relentless.model.potential.FENEWCA(types=("1",)) + FENEWCA.coeff["1"].update(k=1000, r0=1.0, epsilon=1.0, sigma=1.0) + data = FENEWCA.to_json() + + FENEWCA2 = relentless.model.potential.FENEWCA.from_json(data) + self.assertEqual(FENEWCA2.coeff["1"]["k"], 1000) + self.assertEqual(FENEWCA2.coeff["1"]["r0"], 1.0) + + +class test_BondSpline(unittest.TestCase): + """Unit tests for relentless.model.potential.BondSpline""" + + def test_init(self): + """Test creation from data""" + # test diff mode + s = relentless.model.potential.BondSpline(types=("1",), num_knots=3) + self.assertEqual(s.num_knots, 3) + self.assertEqual(s.mode, "diff") + coeff = relentless.model.potential.BondParameters( + types=("1",), + params=( + "dr-0", + "dr-1", + "dr-2", + "diff-0", + "diff-1", + "diff-2", + ), + ) + self.assertCountEqual(s.coeff.types, coeff.types) + self.assertCountEqual(s.coeff.params, coeff.params) + + # test value mode + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=3, mode="value" + ) + self.assertEqual(s.num_knots, 3) + self.assertEqual(s.mode, "value") + coeff = relentless.model.potential.BondParameters( + types=("1",), + params=( + "dr-0", + "dr-1", + "dr-2", + "value-0", + "value-1", + "value-2", + ), + ) + self.assertCountEqual(s.coeff.types, coeff.types) + self.assertCountEqual(s.coeff.params, coeff.params) + + # test invalid number of knots + with self.assertRaises(ValueError): + s = relentless.model.potential.BondSpline(types=("1",), num_knots=1) + + # test invalid mode + with self.assertRaises(ValueError): + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=3, mode="val" + ) + + def test_from_array(self): + """Test from_array method and knots generator""" + r_arr = [1, 2, 3] + u_arr = [9, 4, 1] + u_arr_diff = [5, 3, 1] + + # test diff mode + s = relentless.model.potential.BondSpline(types=("1",), num_knots=3) + s.from_array(types=("1"), r=r_arr, u=u_arr) + + dvars = [] + for i, (r, k) in enumerate(s.knots(types=("1"))): + self.assertAlmostEqual(r.value, 1.0) + self.assertAlmostEqual(k.value, u_arr_diff[i]) + self.assertIsInstance(r, relentless.model.IndependentVariable) + self.assertIsInstance(k, relentless.model.IndependentVariable) + if i < s.num_knots - 1: + dvars.append(k) + self.assertCountEqual(s.design_variables, dvars) + + # test value mode + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=3, mode="value" + ) + s.from_array(types=("1"), r=r_arr, u=u_arr) + + dvars = [] + for i, (r, k) in enumerate(s.knots(types=("1"))): + self.assertAlmostEqual(r.value, 1.0) + self.assertAlmostEqual(k.value, u_arr[i]) + self.assertIsInstance(r, relentless.model.IndependentVariable) + self.assertIsInstance(k, relentless.model.IndependentVariable) + if i != s.num_knots - 1: + dvars.append(k) + self.assertCountEqual(s.design_variables, dvars) + + # test invalid r and u shapes + r_arr = [2, 3] + with self.assertRaises(ValueError): + s.from_array(types=("1"), r=r_arr, u=u_arr) + + r_arr = [1, 2, 3] + u_arr = [1, 2] + with self.assertRaises(ValueError): + s.from_array(types=("1"), r=r_arr, u=u_arr) + + def test_energy(self): + """Test energy method""" + r_arr = [1, 2, 3] + u_arr = [9, 4, 1] + + # test diff mode + s = relentless.model.potential.BondSpline(types=("1",), num_knots=3) + s.from_array(types=("1"), r=r_arr, u=u_arr) + u_actual = numpy.array([6.25, 2.25, 1]) + u = s.energy(type_=("1"), r=[1.5, 2.5, 3.5]) + numpy.testing.assert_allclose(u, u_actual) + + # test value mode + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=3, mode="value" + ) + s.from_array(types=("1"), r=r_arr, u=u_arr) + u_actual = numpy.array([6.25, 2.25, 1]) + u = s.energy(type_=("1"), r=[1.5, 2.5, 3.5]) + numpy.testing.assert_allclose(u, u_actual) + + # test BondSpline with 2 knots + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=2, mode="value" + ) + s.from_array(types=("1"), r=[1, 2], u=[4, 2]) + u = s.energy(type_=("1"), r=1.5) + self.assertAlmostEqual(u, 3) + + def test_force(self): + """Test force method""" + r_arr = [1, 2, 3] + u_arr = [9, 4, 1] + + # test diff mode + s = relentless.model.potential.BondSpline(types=("1",), num_knots=3) + s.from_array(types=("1"), r=r_arr, u=u_arr) + f_actual = numpy.array([5, 3, 0]) + f = s.force(type_=("1"), r=[1.5, 2.5, 3.5]) + numpy.testing.assert_allclose(f, f_actual) + + # test value mode + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=3, mode="value" + ) + s.from_array(types=("1"), r=r_arr, u=u_arr) + f_actual = numpy.array([5, 3, 0]) + f = s.force(type_=("1"), r=[1.5, 2.5, 3.5]) + numpy.testing.assert_allclose(f, f_actual) + + # test BondSpline with 2 knots + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=2, mode="value" + ) + s.from_array(types=("1"), r=[1, 2], u=[4, 2]) + f = s.force(type_=("1"), r=1.5) + self.assertAlmostEqual(f, 2) + + def test_derivative(self): + """Test derivative method""" + r_arr = [1, 2, 3] + u_arr = [9, 4, 1] + + # test diff mode + s = relentless.model.potential.BondSpline(types=("1",), num_knots=3) + s.from_array(types=("1"), r=r_arr, u=u_arr) + d_actual = numpy.array([1.125, 0.625, 0]) + param = list(s.knots(("1")))[1][1] + d = s.derivative(type_=("1"), var=param, r=[1.5, 2.5, 3.5]) + numpy.testing.assert_allclose(d, d_actual) + + # test value mode + s = relentless.model.potential.BondSpline( + types=("1",), num_knots=3, mode="value" + ) + s.from_array(types=("1"), r=r_arr, u=u_arr) + d_actual = numpy.array([0.75, 0.75, 0]) + param = list(s.knots(("1")))[1][1] + d = s.derivative(type_=("1"), var=param, r=[1.5, 2.5, 3.5]) + numpy.testing.assert_allclose(d, d_actual) + + def test_json(self): + r_arr = [1, 2, 3] + u_arr = [9, 4, 1] + u_arr_diff = [5, 3, 1] + + s = relentless.model.potential.BondSpline(types=("1",), num_knots=3) + s.from_array(types=("1"), r=r_arr, u=u_arr) + data = s.to_json() + + s2 = relentless.model.potential.BondSpline.from_json(data) + for i, (r, k) in enumerate(s2.knots(types=("1"))): + self.assertAlmostEqual(r.value, 1.0) + self.assertAlmostEqual(k.value, u_arr_diff[i]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/simulate/test_hoomd.py b/tests/simulate/test_hoomd.py index 604d0539..2dad12a1 100644 --- a/tests/simulate/test_hoomd.py +++ b/tests/simulate/test_hoomd.py @@ -14,6 +14,7 @@ from parameterized import parameterized_class import relentless +from tests.model.potential.test_bond import LinPotBond from tests.model.potential.test_pair import LinPot _has_modules = relentless.simulate.hoomd._hoomd_found @@ -77,6 +78,38 @@ def ens_pot(self): return (ens, pots) + def ens_pot_bonds(self): + if self.dim == 3: + ens = relentless.model.Ensemble( + T=2.0, V=relentless.model.Cube(L=20.0), N={"A": 2, "B": 2} + ) + elif self.dim == 2: + ens = relentless.model.Ensemble( + T=2.0, V=relentless.model.Square(L=20.0), N={"A": 2, "B": 2} + ) + else: + raise ValueError("HOOMD supports 2d and 3d simulations") + + # setup potentials + pot = LinPot(ens.types, params=("m",)) + for pair in pot.coeff: + pot.coeff[pair].update({"m": -2.0, "rmax": 1.0}) + pots = relentless.simulate.Potentials() + pots.pair = relentless.simulate.PairPotentialTabulator( + pot, start=0.0, stop=2.0, num=3, neighbor_buffer=0.4 + ) + bond = LinPotBond(("bondA", "bondB"), params=("m",)) + for type_ in bond.coeff: + bond.coeff[type_].update({"m": -2.0}) + pots.bond = relentless.simulate.BondPotentialTabulator( + bond, + start=0.0, + stop=2.0, + num=3, + ) + + return (ens, pots) + # mock gsd file for testing def create_gsd_file(self): filename = self.directory.file("test.gsd") @@ -103,6 +136,41 @@ def create_gsd_file(self): relentless.mpi.world.barrier() return filename + # mock gsd file with topology for testing + def create_gsd_file_bonds(self): + filename = self.directory.file("test.gsd") + if relentless.mpi.world.rank_is_root: + with gsd.hoomd.open(name=filename, mode=gsd_write_mode) as f: + s = HOOMDFrame() + s.particles.N = 4 + s.particles.types = ["A", "B"] + s.particles.typeid = [0, 0, 1, 1] + s.bonds.N = 2 + s.bonds.types = ["bondA", "bondB"] + s.bonds.group = [(0, 1), (2, 3)] + s.bonds.typeid = [0, 1] + if self.dim == 3: + s.particles.position = [ + [0, 0, 0], + [1, 0, 0], + [0, 2, 1], + [1, 2, 1], + ] + s.configuration.box = [20, 20, 20, 0, 0, 0] + elif self.dim == 2: + s.particles.position = [ + [0, 0, 0], + [1, 0, 0], + [0, 2, 0], + [1, 2, 0], + ] + s.configuration.box = [20, 20, 0, 0, 0, 0] + else: + raise ValueError("HOOMD supports 2d and 3d simulations") + f.append(s) + relentless.mpi.world.barrier() + return filename + def create_lammps_file(self): file_ = self.directory.file("test.data") @@ -120,6 +188,44 @@ def create_lammps_file(self): return file_ + def create_lammps_file_bonds(self): + file_ = self.directory.file("test.data") + + if relentless.mpi.world.rank_is_root: + low = [-5, -5, -5 if self.dim == 3 else -0.1] + high = [5, 5, 5 if self.dim == 3 else 0.1] + snap = lammpsio.Snapshot(N=4, box=lammpsio.Box(low, high)) + snap.position = [ + [0, 0, 0], + [1, 0, 0], + [0, 2, 0], + [1, 2, 0], + ] + if self.dim == 3: + snap.position = [ + [0, 0, 0], + [1, 0, 0], + [0, 2, 1], + [1, 2, 1], + ] + snap.typeid = [ + 1, + 1, + 2, + 2, + ] + snap.mass = [0.3, 0.3, 0.1, 0.1] + snap.bonds = lammpsio.topology.Bonds(N=2, num_types=2) + snap.bonds.id = [1, 2] + snap.bonds.typeid = [1, 2] + snap.bonds.members = [(1, 2), (3, 4)] + snap.bonds.type_label = lammpsio.topology.LabelMap({1: "bondA", 2: "bondB"}) + + lammpsio.DataFile.create(file_, snap) + relentless.mpi.world.barrier() + + return file_ + def test_initialize_from_gsd_file(self): ens, pot = self.ens_pot() f = self.create_gsd_file() @@ -137,6 +243,23 @@ def test_initialize_from_gsd_file(self): h.initializer = relentless.simulate.InitializeFromFile(pathlib.Path(f).name) h.run(pot, d) + def test_initialize_from_gsd_file_topology(self): + ens, pot = self.ens_pot_bonds() + f = self.create_gsd_file_bonds() + op = relentless.simulate.InitializeFromFile(filename=f) + h = relentless.simulate.HOOMD(op) + h.run(pot, self.directory) + + # Run in a different directory + with self.directory: + d = self.directory.directory( + "run", create=relentless.mpi.world.rank_is_root + ) + relentless.mpi.world.barrier() + op.filename = pathlib.Path(f).name + h.initializer = relentless.simulate.InitializeFromFile(pathlib.Path(f).name) + h.run(pot, d) + def test_initialize_from_lammps_file(self): """Test running initialization simulation operations.""" pot = LinPot(("1", "2"), params=("m",)) @@ -152,6 +275,21 @@ def test_initialize_from_lammps_file(self): h = relentless.simulate.HOOMD(op) h.run(pots, self.directory) + def test_initialize_from_lammps_file_topology(self): + """Test running initialization simulation operations.""" + pot = LinPot(("1", "2"), params=("m",)) + for pair in pot.coeff: + pot.coeff[pair].update({"m": -2.0, "rmax": 1.0}) + pots = relentless.simulate.Potentials() + pots.pair = relentless.simulate.PairPotentialTabulator( + pot, start=1e-6, stop=2.0, num=10, neighbor_buffer=0.1 + ) + + f = self.create_lammps_file_bonds() + op = relentless.simulate.InitializeFromFile(filename=f) + h = relentless.simulate.HOOMD(op) + h.run(pots, self.directory) + def test_initialize_randomly(self): ens, pot = self.ens_pot() op = relentless.simulate.InitializeRandomly(seed=1, N=ens.N, V=ens.V, T=ens.T) @@ -313,6 +451,19 @@ def test_molecular_dynamics(self): vrl.barostat = relentless.simulate.MTKBarostat(P=1, tau=0.5) h.run(pot, self.directory) + def test_topology_run(self): + ens, pot = self.ens_pot_bonds() + f = self.create_gsd_file_bonds() + init = relentless.simulate.InitializeFromFile(filename=f) + lgv = relentless.simulate.RunLangevinDynamics( + steps=1, timestep=1.0e-3, T=ens.T, friction=1.0, seed=2 + ) + h = relentless.simulate.HOOMD(init, lgv) + h.run(pot, self.directory) + + lgv.friction = {"A": 1.5, "B": 2.5} + h.run(pot, self.directory) + def test_temperature_ramp(self): if self.dim == 3: V = relentless.model.Cube(100.0) diff --git a/tests/simulate/test_lammps.py b/tests/simulate/test_lammps.py index 5e385d83..7a46d043 100644 --- a/tests/simulate/test_lammps.py +++ b/tests/simulate/test_lammps.py @@ -34,6 +34,7 @@ from packaging import version import relentless +from tests.model.potential.test_bond import LinPotBond from tests.model.potential.test_pair import LinPot # silence warnings about Snapshot being deprecated @@ -106,6 +107,38 @@ def ens_pot(self): return (ens, pots) + def ens_pot_bonds(self): + if self.dim == 3: + ens = relentless.model.Ensemble( + T=2.0, V=relentless.model.Cube(L=20.0), N={"1": 2, "2": 2} + ) + elif self.dim == 2: + ens = relentless.model.Ensemble( + T=2.0, V=relentless.model.Square(L=20.0), N={"1": 2, "2": 2} + ) + else: + raise ValueError("HOOMD supports 2d and 3d simulations") + + # setup potentials + pot = LinPot(ens.types, params=("m",)) + for pair in pot.coeff: + pot.coeff[pair].update({"m": -2.0, "rmax": 1.0}) + pots = relentless.simulate.Potentials() + pots.pair = relentless.simulate.PairPotentialTabulator( + pot, start=1e-6, stop=2.0, num=3, neighbor_buffer=0.4 + ) + bond = LinPotBond(("bond1", "bond2"), params=("m",)) + for type_ in bond.coeff: + bond.coeff[type_].update({"m": -1.0}) + pots.bond = relentless.simulate.BondPotentialTabulator( + bond, + start=0.0, + stop=4.0, + num=3, + ) + + return (ens, pots) + def create_gsd_file(self): filename = self.directory.file("test.gsd") if relentless.mpi.world.rank_is_root: @@ -131,6 +164,41 @@ def create_gsd_file(self): relentless.mpi.world.barrier() return filename + # mock gsd file with topology for testing + def create_gsd_file_bonds(self): + filename = self.directory.file("test.gsd") + if relentless.mpi.world.rank_is_root: + with gsd.hoomd.open(name=filename, mode=gsd_write_mode) as f: + s = HOOMDFrame() + s.particles.N = 4 + s.particles.types = ["A", "B"] + s.particles.typeid = [0, 0, 1, 1] + s.bonds.N = 2 + s.bonds.types = ["bondA", "bondB"] + s.bonds.group = [(0, 1), (2, 3)] + s.bonds.typeid = [0, 1] + if self.dim == 3: + s.particles.position = [ + [0, 0, 0], + [1, 0, 0], + [0, 2, 1], + [1, 2, 1], + ] + s.configuration.box = [20, 20, 20, 0, 0, 0] + elif self.dim == 2: + s.particles.position = [ + [0, 0, 0], + [1, 0, 0], + [0, 2, 0], + [1, 2, 0], + ] + s.configuration.box = [20, 20, 0, 0, 0, 0] + else: + raise ValueError("HOOMD supports 2d and 3d simulations") + f.append(s) + relentless.mpi.world.barrier() + return filename + def create_lammps_file(self): file_ = self.directory.file("test.data") @@ -148,6 +216,28 @@ def create_lammps_file(self): return file_ + def create_lammps_file_bonds(self): + file_ = self.directory.file("test.data") + + if relentless.mpi.world.rank_is_root: + low = [-5, -5, -5 if self.dim == 3 else -0.1] + high = [5, 5, 5 if self.dim == 3 else 0.1] + snap = lammpsio.Snapshot(N=5, box=lammpsio.Box(low, high)) + snap.position[:, :2] = [[-4, -4], [-2, -2], [0, 0], [2, 2], [4, 4]] + if self.dim == 3: + snap.position[:, 2] = [-4, -2, 0, 2, 4] + snap.typeid = [1, 1, 2, 2, 2] + snap.mass = [0.3, 0.3, 0.1, 0.1, 0.1] + + snap.bonds = lammpsio.topology.Bonds(N=2, num_types=2) + snap.bonds.id = [1, 2] + snap.bonds.typeid = [1, 2] + snap.bonds.members = [(1, 2), (3, 4)] + lammpsio.DataFile.create(file_, snap) + relentless.mpi.world.barrier() + + return file_ + def test_initialize_from_lammps_file(self): """Test running initialization simulation operations.""" # InitializeFromFile @@ -179,6 +269,33 @@ def test_initialize_from_lammps_file(self): lmp = relentless.simulate.LAMMPS(op, executable=self.executable) lmp.run(potentials=pot, directory=self.directory) + def test_initialize_from_lammps_file_topology(self): + """Test running initialization simulation operations.""" + # InitializeFromFile + ens, pot = self.ens_pot_bonds() + f = self.create_lammps_file_bonds() + op = relentless.simulate.InitializeFromFile(filename=f) + h = relentless.simulate.LAMMPS(op, atom_style="molecular") + h.run(pot, self.directory) + + lmp = relentless.simulate.LAMMPS( + op, + executable=self.executable, + atom_style="molecular", + ) + # Run in a different directory + with self.directory: + d = self.directory.directory( + "run", create=relentless.mpi.world.rank_is_root + ) + relentless.mpi.world.barrier() + op.filename = pathlib.Path(f).name + lmp.initializer = relentless.simulate.InitializeFromFile( + pathlib.Path(f).name + ) + + lmp.run(potentials=pot, directory=d) + def test_initialize_from_gsd_file(self): """Test running initialization simulation operations.""" _, pot = self.ens_pot() @@ -191,6 +308,31 @@ def test_initialize_from_gsd_file(self): ) lmp.run(potentials=pot, directory=self.directory) + def test_initialize_from_gsd_file_topology(self): + """Test running initialization simulation operations.""" + # setup potentials + pot = LinPot(("A", "B"), params=("m",)) + for pair in pot.coeff: + pot.coeff[pair].update({"m": -2.0, "rmax": 1.0}) + pots = relentless.simulate.Potentials() + pots.pair = relentless.simulate.PairPotentialTabulator( + pot, start=1e-6, stop=2.0, num=3, neighbor_buffer=0.4 + ) + bond = LinPotBond(("bondA", "bondB"), params=("m",)) + for type_ in bond.coeff: + bond.coeff[type_].update({"m": -2.0}) + pots.bond = relentless.simulate.BondPotentialTabulator( + bond, + start=0.0, + stop=2.0, + num=3, + ) + + f = self.create_gsd_file_bonds() + op = relentless.simulate.InitializeFromFile(filename=f) + h = relentless.simulate.LAMMPS(op, atom_style="molecular") + h.run(pots, self.directory) + def test_random_initialize_options(self): # no T ens, pot = self.ens_pot() @@ -418,6 +560,17 @@ def test_molecular_dynamics(self): vrl.thermostat = None lmp.run(pot, self.directory) + def test_topology_run(self): + ens, pot = self.ens_pot_bonds() + f = self.create_lammps_file_bonds() + init = relentless.simulate.InitializeFromFile(filename=f) + lmp = relentless.simulate.LAMMPS(init, atom_style="molecular") + + # VerletIntegrator + vrl = relentless.simulate.RunMolecularDynamics(steps=1, timestep=1e-3) + lmp.operations = vrl + lmp.run(pot, self.directory) + def test_temperature_ramp(self): if self.dim == 3: V = relentless.model.Cube(100.0)