Add support for hoomd.md.improper.Periodic in to_hoomd_forcefield (

#807)

* add dict entry for periodic impropers

* add sign factor (d) to Periodic Improper

* Change key name to match GAFF, hard code d = 1

* remove d param from template

* Remove d from expression

* change d from int to float

* update python version, remove ethanol from unit test

* if wild card in member classes, use member types instead

* use wildcard check for all forcefield writers

* fix typo in env files

* Fix variable name.

* use .get() for d parameter

* change box size and seed for packing in test

---------

Co-authored-by: Co Quach <43968221+daico007@users.noreply.github.com>

environment-dev.yml

@@ -12,7 +12,7 @@ dependencies:
   - pydantic>=2
   - networkx
   - pytest
-  - mbuild>=0.11.0
+  - mbuild>=0.17.0
   - openbabel>=3.0.0
   - foyer>=0.11.3
   - forcefield-utilities>=0.2.1

gmso/external/convert_hoomd.py

@@ -19,7 +19,11 @@
 from gmso.utils.conversions import convert_ryckaert_to_opls
 from gmso.utils.geometry import coord_shift
 from gmso.utils.io import has_gsd, has_hoomd
-from gmso.utils.sorting import sort_by_classes, sort_connection_members
+from gmso.utils.sorting import (
+    sort_by_classes,
+    sort_by_types,
+    sort_connection_members,
+)
 from gmso.utils.units import convert_params_units
 
 if has_gsd:
@@ -40,6 +44,8 @@
     "mass": u.g / u.mol,  # aka amu
 }
 
+hoomd_version = hoomd.version.version.split(".")
+
 
 def to_gsd_snapshot(
     top,
@@ -665,13 +671,15 @@ def _validate_compatibility(top):
     harmonic_bond_potential = templates["HarmonicBondPotential"]
     harmonic_angle_potential = templates["HarmonicAnglePotential"]
     periodic_torsion_potential = templates["PeriodicTorsionPotential"]
+    harmonic_torsion_potential = templates["HarmonicTorsionPotential"]
     opls_torsion_potential = templates["OPLSTorsionPotential"]
     rb_torsion_potential = templates["RyckaertBellemansTorsionPotential"]
     accepted_potentials = (
         lennard_jones_potential,
         harmonic_bond_potential,
         harmonic_angle_potential,
         periodic_torsion_potential,
+        harmonic_torsion_potential,
         opls_torsion_potential,
         rb_torsion_potential,
     )
@@ -968,8 +976,11 @@ def _parse_harmonic_bond(
 ):
     for btype in btypes:
         # TODO: Unit conversion
-        member_classes = sort_by_classes(btype)
-        container.params["-".join(member_classes)] = {
+        members = sort_by_classes(btype)
+        # If wild card in class, sort by types instead
+        if "*" in members:
+            members = sort_by_types(btype)
+        container.params["-".join(members)] = {
             "k": btype.parameters["k"],
             "r0": btype.parameters["r_eq"],
         }
@@ -1034,8 +1045,11 @@ def _parse_harmonic_angle(
     agtypes,
 ):
     for agtype in agtypes:
-        member_classes = sort_by_classes(agtype)
-        container.params["-".join(member_classes)] = {
+        members = sort_by_classes(agtype)
+        # If wild card in class, sort by types instead
+        if "*" in members:
+            members = sort_by_types(agtype)
+        container.params["-".join(members)] = {
             "k": agtype.parameters["k"],
             "t0": agtype.parameters["theta_eq"],
         }
@@ -1091,7 +1105,6 @@ def _parse_dihedral_forces(
         },
     }
 
-    hoomd_version = hoomd.version.version.split(".")
     if int(hoomd_version[0]) >= 4 or (
         int(hoomd_version[0]) == 3 and int(hoomd_version[1]) >= 8
     ):
@@ -1262,12 +1275,24 @@ def _parse_improper_forces(
             base_units,
         )
 
-    itype_group_map = {
-        "HarmonicImproperPotenial": {
-            "container": hoomd.md.improper.Harmonic,
-            "parser": _parse_harmonic_improper,
-        },
-    }
+    if int(hoomd_version[0]) >= 4 and int(hoomd_version[1]) >= 5:
+        itype_group_map = {
+            "HarmonicImproperPotential": {
+                "container": hoomd.md.improper.Harmonic,
+                "parser": _parse_harmonic_improper,
+            },
+            "PeriodicTorsionPotential": {
+                "container": hoomd.md.improper.Periodic,
+                "parser": _parse_periodic_improper,
+            },
+        }
+    else:
+        itype_group_map = {
+            "HarmonicImproperPotential": {
+                "container": hoomd.md.improper.Harmonic,
+                "parser": _parse_harmonic_improper,
+            },
+        }
 
     improper_forces = list()
     for group in groups:
@@ -1285,14 +1310,35 @@ def _parse_harmonic_improper(
     itypes,
 ):
     for itype in itypes:
-        member_types = sort_by_classes(itype)
-        container.params["-".join(member_types)] = {
+        members = sort_by_classes(itype)
+        # If wild card in class, sort by types instead
+        if "*" in members:
+            members = sort_by_types(itype)
+        container.params["-".join(members)] = {
             "k": itype.parameters["k"],
             "chi0": itype.parameters["phi_eq"],  # diff nomenclature?
         }
     return container
 
 
+def _parse_periodic_improper(
+    container,
+    itypes,
+):
+    for itype in itypes:
+        members = sort_by_classes(itype)
+        # If wild card in class, sort by types instead
+        if "*" in members:
+            members = sort_by_types(itype)
+        container.params["-".join(members)] = {
+            "k": itype.parameters["k"],
+            "chi0": itype.parameters["phi_eq"],
+            "n": itype.parameters["n"],
+            "d": itype.parameters.get("d", 1.0),
+        }
+    return container
+
+
 def _validate_base_units(base_units, top, auto_scale, potential_types=None):
     """Validate the provided base units, infer units (based on top's positions and masses) if none is provided."""
     if base_units and auto_scale:

gmso/tests/test_convert_mbuild.py

@@ -193,9 +193,9 @@ def test_custom_groups_from_compound(self):
         mb_cpd3 = mb.load("O", smiles=True)
         mb_cpd3.name = "O"
 
-        filled_box1 = mb.fill_box([mb_cpd1, mb_cpd2], n_compounds=[2, 2], box=[1, 1, 1])
+        filled_box1 = mb.fill_box([mb_cpd1, mb_cpd2], n_compounds=[2, 2], box=[5, 5, 5])
         filled_box1.name = "box1"
-        filled_box2 = mb.fill_box(mb_cpd3, n_compounds=2, box=[1, 1, 1])
+        filled_box2 = mb.fill_box(mb_cpd3, n_compounds=2, box=[5, 5, 5])
         filled_box2.name = "box2"
 
         top_box = mb.Compound()
@@ -222,7 +222,9 @@ def test_custom_groups_from_compound(self):
     def test_single_custom_group(self):
         mb_cpd1 = mb.Compound(name="_CH4")
         mb_cpd2 = mb.Compound(name="_CH3")
-        filled_box = mb.fill_box([mb_cpd1, mb_cpd2], n_compounds=[2, 2], box=[1, 1, 1])
+        filled_box = mb.fill_box(
+            [mb_cpd1, mb_cpd2], n_compounds=[2, 2], box=[5, 5, 5], seed=20
+        )
         filled_box.name = "box1"
 
         top = from_mbuild(filled_box, custom_groups=filled_box.name)
@@ -235,7 +237,9 @@ def test_single_custom_group(self):
     def test_bad_custom_groups_from_compound(self):
         mb_cpd1 = mb.Compound(name="_CH4")
         mb_cpd2 = mb.Compound(name="_CH3")
-        filled_box = mb.fill_box([mb_cpd1, mb_cpd2], n_compounds=[2, 2], box=[1, 1, 1])
+        filled_box = mb.fill_box(
+            [mb_cpd1, mb_cpd2], n_compounds=[2, 2], box=[5, 5, 5], seed=13
+        )
 
         with pytest.warns(Warning):
             from_mbuild(filled_box, custom_groups=["_CH4", "_CH3", "_CH5"])

gmso/tests/test_hoomd.py

@@ -396,9 +396,9 @@ def test_gaff_sim(self, gaff_forcefield):
             "length": u.nm,
             "energy": u.kJ / u.mol,
         }
-        ethanol = mb.load("CCO", smiles=True)
-        ethanol.box = mb.Box([5, 5, 5])
-        top = ethanol.to_gmso()
+        benzene = mb.load("c1ccccc1", smiles=True)
+        benzene.box = mb.Box([5, 5, 5])
+        top = benzene.to_gmso()
         parameterized_top = apply(top, gaff_forcefield, identify_connections=True)
         assert parameterized_top.is_fully_typed