Handlers for template potentials that are different forms of the same…

… equation. (#855)

* Add support for modifying template expression via set expression, and testing for handling accepted_potentials as scaled versions of template potentials

* Remove abstract_potential set_expression method

* Fix typing and remove unused variables

* Improve type checking in conversions.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* update docstrings

* Fix check for importing types

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* Add test for failing conversion handling

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Chris Jones <[email protected]>

gmso/abc/abstract_potential.py

@@ -1,6 +1,5 @@
 """Abstract representation of a Potential object."""
 
-from abc import abstractmethod
 from typing import Any, Dict, Iterator, List
 
 import unyt as u
@@ -160,11 +159,6 @@ def validate_potential_expression(cls, v):
             v = PotentialExpression(**v)
         return v
 
-    @abstractmethod
-    def set_expression(self):
-        """Set the functional form of the expression."""
-        raise NotImplementedError
-
     def __setattr__(self, key: Any, value: Any) -> None:
         """Set attributes of the potential."""
         if key == "expression":

gmso/lib/potential_templates.py

@@ -1,6 +1,7 @@
 """Module supporting template potential objects."""
 
 import json
+from copy import deepcopy
 from pathlib import Path
 from typing import Dict
 
@@ -113,9 +114,11 @@ def expected_parameters_dimensions(self):
         """Return the expected dimensions of the parameters for this template."""
         return self.__dict__.get("expected_parameters_dimensions_")
 
-    def set_expression(self, *args, **kwargs):
+    def set_expression(self, new_expression):
         """Set the expression of the PotentialTemplate."""
-        raise NotImplementedError
+        copied_template = deepcopy(self)
+        copied_template.expression = new_expression
+        return copied_template
 
     def assert_can_parameterize_with(
         self, parameters: Dict[str, u.unyt_quantity]

gmso/tests/test_conversions.py

@@ -2,7 +2,9 @@
 import pytest
 import sympy
 import unyt as u
+from sympy import sympify
 
+from gmso.exceptions import EngineIncompatibilityError
 from gmso.tests.base_test import BaseTest
 from gmso.utils.conversions import (
     convert_kelvin_to_energy_units,
@@ -17,10 +19,21 @@ def _convert_potential_types(top, connStr, expected_units_dim, base_units):
     return potentials
 
 
-class TestKelvinToEnergy(BaseTest):
-    def test_convert_potential_styles(self, typed_ethane):
-        from sympy import sympify
+class TestConversions(BaseTest):
+    def test_rescale_potentials(self, typed_ethane):
+        from gmso.lib.potential_templates import PotentialTemplateLibrary
+
+        library = PotentialTemplateLibrary()
+        template = library["LennardJonesPotential"]
+        template = template.set_expression(
+            template.expression / 4
+        )  # use setter to not set in place
+        atype = list(typed_ethane.atom_types)[0]
+        assert atype.expression == sympify("4*epsilon*((sigma/r)**12 - (sigma/r)**6)")
+        typed_ethane.convert_potential_styles({"sites": template})
+        assert atype.expression == sympify("epsilon*((sigma/r)**12 - (sigma/r)**6)")
 
+    def test_convert_potential_styles(self, typed_ethane):
         rb_expr = sympify(
             "c0 * cos(phi)**0 + c1 * cos(phi)**1 + c2 * cos(phi)**2 + c3 * cos(phi)**3 + c4 * cos(phi)**4 + c5 * cos(phi)**5"
         )
@@ -227,3 +240,9 @@ def test_lammps_conversion_parameters_lj(self):
             n_decimals=6,
         )
         assert np.isclose(float(paramStr), 1 / 3**4, atol=1e-6)
+
+    def test_failed_conversion_method(self, typed_ethane):
+        with pytest.raises(
+            EngineIncompatibilityError,
+        ):
+            typed_ethane.convert_potential_styles({"bond": 2})

gmso/tests/test_hoomd.py

@@ -388,7 +388,7 @@ def test_zero_charges(self):
     @pytest.mark.skipif(not has_hoomd, reason="hoomd is not installed")
     @pytest.mark.skipif(not has_mbuild, reason="mbuild not installed")
     @pytest.mark.skipif(
-        int(hoomd_version[0]) <= 3.8, reason="Deprecated features in HOOMD 4"
+        int(hoomd_version[0]) < 4.5, reason="No periodic impropers in hoomd < 4.5"
     )
     def test_gaff_sim(self, gaff_forcefield):
         base_units = {

gmso/tests/test_template.py

@@ -34,8 +34,10 @@ def test_template_set_expression(self):
             independent_variables={"x"},
             expected_parameters_dimensions={"a": "length", "b": "length"},
         )
-        with pytest.raises(NotImplementedError):
-            template.set_expression(expression="a*y+b")
+        with pytest.raises(ValueError):
+            template.set_expression(new_expression="a*y+b")
+        template2 = template.set_expression(new_expression="3*a*x+b")
+        assert template2.expression != template.expression
 
     def test_parameterization_non_dict_expected_dimensions(self):
         template = PotentialTemplate(

gmso/utils/conversions.py

@@ -1,5 +1,7 @@
 """Module for standard conversions needed in molecular simulations."""
 
+from __future__ import annotations
+
 import re
 from functools import lru_cache
 
@@ -10,20 +12,20 @@
 from unyt.dimensions import length, mass, time
 
 import gmso
-from gmso.exceptions import GMSOError
-from gmso.lib.potential_templates import PotentialTemplateLibrary
+from gmso.exceptions import EngineIncompatibilityError, GMSOError
+from gmso.lib.potential_templates import (
+    PotentialTemplate,
+    PotentialTemplateLibrary,
+)
 
 templates = PotentialTemplateLibrary()
 
 
 @lru_cache(maxsize=128)
 def _constant_multiplier(pot1, pot2):
-    # TODO: Doc string
-    # TODO: Test outputs
-    # TODO: Check speed
+    """Attempt to convert from pot1 to pot2 using a scalar."""
     try:
         constant = symengine.expand(pot1.expression / pot2.expression)
-        # constant = sympy.simplify(pot1.expression / pot2.expression)
         if constant.is_number:
             for eq_term in pot1.expression.args:
                 if eq_term.is_symbol:
@@ -40,9 +42,7 @@ def _constant_multiplier(pot1, pot2):
 
 @lru_cache(maxsize=128)
 def _try_sympy_conversions(pot1, pot2):
-    # TODO: Doc string
-    # TODO: Test outputs
-    # TODO: Check speed
+    """Attempt to convert from pot1 to pot2 using any generalized conversion function."""
     convertersList = []
     for conversion in sympy_conversionsList:
         convertersList.append(conversion(pot1, pot2))
@@ -52,23 +52,12 @@ def _try_sympy_conversions(pot1, pot2):
     return None
 
 
-def convert_topology_expressions(top, expressionMap={}):
-    """Convert from one parameter form to another.
-
-    Parameters
-    ----------
-    expressionMap : dict, default={}
-        map with keys of the potential type and the potential to change to
-
-    Examples
-    --------
-    Convert from RB torsions to OPLS torsions
-    top.convert_expressions({"dihedrals": "OPLSTorsionPotential"})
-    """
-    # TODO: Raise errors
-
-    # Apply from predefined conversions or easy sympy conversions
-    conversions_map = {
+def _conversion_from_template_name(
+    top, connStr: str, conn_typeStr: str, convStr: str
+) -> "gmso.Topology":
+    """Use the name of convStr to identify function to convert sympy expressions."""
+    conversions_map = {  # these are predefined between template types
+        # More functions, and `(to, from)` key pairs added to this dictionary
         (
             "OPLSTorsionPotential",
             "RyckaertBellemansTorsionPotential",
@@ -83,40 +72,112 @@ def convert_topology_expressions(top, expressionMap={}):
         ): convert_ryckaert_to_opls,
     }  # map of all accessible conversions currently supported
 
-    for conv in expressionMap:
-        # check all connections with these types for compatibility
-        for conn in getattr(top, conv):
-            current_expression = getattr(conn, conv[:-1] + "_type")
-            if (
-                current_expression.name == expressionMap[conv]
-            ):  # check to see if we can skip this one
-                # TODO: Do something instead of just comparing the names
-                continue
-
-            # convert it using pre-defined conversion functions
-            conversion_from_conversion_toTuple = (
-                current_expression.name,
-                expressionMap[conv],
-            )
-            if (
-                conversion_from_conversion_toTuple in conversions_map
-            ):  # Try mapped conversions
-                new_conn_type = conversions_map.get(conversion_from_conversion_toTuple)(
-                    current_expression
-                )
-                setattr(conn, conv[:-1] + "_type", new_conn_type)
-                continue
-
-            # convert it using sympy expression conversion
-            new_potential = templates[expressionMap[conv]]
-            modified_connection_parametersDict = _try_sympy_conversions(
-                current_expression, new_potential
+    # check all connections with these types for compatibility
+    for conn in getattr(top, connStr):
+        current_expression = getattr(conn, conn_typeStr[:-1])  # strip off extra s
+        # convert it using pre-defined names with conversion functions
+        conversion_from_conversion_toTuple = (current_expression.name, convStr)
+        if (
+            conversion_from_conversion_toTuple in conversions_map
+        ):  # Try mapped conversions
+            new_conn_type = conversions_map.get(conversion_from_conversion_toTuple)(
+                current_expression
             )
-            if modified_connection_parametersDict:  # try sympy conversions
-                current_expression.name = new_potential.name
-                current_expression.expression = new_potential.expression
-                current_expression.parameters.update(modified_connection_parametersDict)
+            setattr(conn, conn_typeStr[:-1], new_conn_type)
+            continue
+
+        # convert it using sympy expression conversion (handles constant multipliers)
+        new_potential = templates[convStr]
+        modified_connection_parametersDict = _try_sympy_conversions(
+            current_expression, new_potential
+        )
+        if modified_connection_parametersDict:  # try sympy conversions
+            current_expression.name = new_potential.name
+            current_expression.expression = new_potential.expression
+            current_expression.parameters.update(modified_connection_parametersDict)
+
 
+def _conversion_from_template_obj(
+    top: "gmso.Topology",
+    connStr: str,
+    conn_typeStr: str,
+    potential_template: gmso.core.ParametricPotential,
+):
+    """Use a passed template object to identify conversion between expressions."""
+    for conn in getattr(top, connStr):
+        current_expression = getattr(conn, conn_typeStr[:-1])  # strip off extra s
+
+        # convert it using sympy expression conversion (handles constant multipliers)
+        modified_connection_parametersDict = _try_sympy_conversions(
+            current_expression, potential_template
+        )
+        if modified_connection_parametersDict:  # try sympy conversions
+            current_expression.name = potential_template.name
+            current_expression.expression = potential_template.expression
+            current_expression.parameters.update(modified_connection_parametersDict)
+
+
+def convert_topology_expressions(top, expressionMap={}):
+    """Convert from one parameter form to another.
+
+    Parameters
+    ----------
+    expressionMap : dict, default={}
+        map with keys of the potential type and the potential to change to.
+        Note that all of the possible template can be found in `gmso/lib/jsons`.
+        Use the string for the json property `name`, or load the template and pass
+        that as the value for conversion.
+
+    Examples
+    --------
+    Convert from RB torsions to OPLS torsions
+    top.convert_expressions({"dihedrals": "OPLSTorsionPotential"})
+
+    Convert from LJ sites to 1*epsilon LJ sites
+    ```
+    template = gmso.lib.potential_templates.PotentialTemplatesLibrary()["LennardJonesPotential"]
+    template = template.set_expression(template.expression / 4) # modify equation
+    top.convert_expressions({"sites":template})
+    # or
+    top = gmso.utils.conversion.convert_topology_expression(top, {"sites":template})
+    ```
+    """
+    # Apply from predefined conversions or easy sympy conversions
+    # handler for various keys passed to expressionMap for conversion
+    for connStr, conv in expressionMap.items():
+        possible_connections = ["bond", "angle", "dihedral", "improper"]
+        if connStr.lower() in [
+            "sites",
+            "site",
+            "atom",
+            "atoms",
+            "atom_types",
+            "atom_type",
+            "atomtype",
+            "atomtypes",
+        ]:
+            # handle renaming type names in relationship to the site or connection
+            conn_typeStr = "atom_types"
+            connStr = "sites"
+        for possible_connection in possible_connections:
+            if possible_connection in connStr.lower():
+                connStr = possible_connection + "s"
+                conn_typeStr = possible_connection + "_types"
+                break
+
+        if isinstance(conv, str):
+            _conversion_from_template_name(top, connStr, conn_typeStr, conv)
+        elif isinstance(conv, PotentialTemplate):
+            _conversion_from_template_obj(top, connStr, conn_typeStr, conv)
+        else:
+            connType = list(getattr(top, conn_typeStr))[0]
+            errormsg = f"""
+            Failed to convert {top} for {connStr} components, with conversion
+            of {connType.name}: Attempted to convert {connType} with style {conv}, which is a {type(conv)}.
+            Please use a template conversion from gmso.lib.potential_templates by passing the name: i.e.
+            HarmonicBondPotential, or by loading a template from gmso.lib.potential_template PotentialTemplateLibrary().
+            """
+            raise EngineIncompatibilityError(errormsg)
     return top
 
 
@@ -132,7 +193,6 @@ def convert_opls_to_ryckaert(opls_connection_type):
     for OPLS and RB torsions. OPLS torsions are defined with
     phi_cis = 0 while RB torsions are defined as phi_trans = 0.
     """
-    # TODO: this function really converts the fourier torsion to rb, not opls
     opls_torsion_potential = templates["FourierTorsionPotential"]
     valid_connection_type = False
     if (