Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Handlers for template potentials that are different forms of the same equation. #855

Merged
merged 12 commits into from
Oct 15, 2024
Merged
6 changes: 0 additions & 6 deletions gmso/abc/abstract_potential.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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":
Expand Down
7 changes: 5 additions & 2 deletions gmso/lib/potential_templates.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
"""Module supporting template potential objects."""

import json
from copy import deepcopy
from pathlib import Path
from typing import Dict

Expand Down Expand Up @@ -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):
Fixed Show fixed Hide fixed
"""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]
Expand Down
25 changes: 22 additions & 3 deletions gmso/tests/test_conversions.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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"
)
Expand Down Expand Up @@ -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})
2 changes: 1 addition & 1 deletion gmso/tests/test_hoomd.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {
Expand Down
6 changes: 4 additions & 2 deletions gmso/tests/test_template.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
178 changes: 119 additions & 59 deletions gmso/utils/conversions.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
"""Module for standard conversions needed in molecular simulations."""

from __future__ import annotations

import re
from functools import lru_cache

Expand All @@ -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:
Expand All @@ -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))
Expand All @@ -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",
Expand All @@ -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


Expand All @@ -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 (
Expand Down
Loading