Repository: mosdef-hub/mbuild (language: python)

Issue:
Add warnings for deprecated file formats (protobuf, hoomdxml, Hoomd-Blue and LAMMPS simulation writers)
### PR Summary:

Per the last few dev meetings, we've been discussing what needs to be done before releasing 1.0. We decided to clean up and consolidate `conversion.py` which involves moving the backend to GMSO wherever possible, and removing `.hoomdxml` and `.pb2` readers and writers. Also, we decided to remove any `mbuild.format` writers that are now duplicated by GMSO such as `create_hoomd_simulation` and `create_hoomd_forcefield`, and `write_lammpsdata`. 


We decided that the first step was including the necessary deprecation warnings in a minor release, which this PR gets started. I think this covers everything that will be removed from mBuild, and any other changes to readers and writers will just involve changing the backend and won't require a deprecation warning. Let me know if I'm missing anything.

I can start making the changes needed in a separate PR.


Discussion / hints:
oops..I guess this includes a bunch of old commits. I'm not sure how that happened, as I started a new branch from `main`. Regardless, it looks like the only actual changs are ones I made for this PR specifically. 
Ok, I didn't look closely enough at the `deprecated` decorator in utils, I thought it just threw a warning message, but it actually deprecates the function its used on, which we don't want yet.  I'll make a new PR.
I found the issue..

The deprecated decorator in utils was missing a return statement, so it was just returning `None`.  I'll clean up my fork and branch and re-open this.

Fix patch:
diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml
index 864831ffd..d2ce820de 100644
--- a/.github/workflows/CI.yaml
+++ b/.github/workflows/CI.yaml
@@ -54,12 +54,6 @@ jobs:
       - name: Test (OS -> ${{ matrix.os }} / Python -> ${{ matrix.python-version }})
         run: python -m pytest -v --cov=mbuild --cov-report=xml --cov-append --cov-config=setup.cfg --color yes --pyargs mbuild
 
-      - name: Tests for HOOMD-blue (V3)
-        run: |
-          micromamba install -c conda-forge "hoomd>3"
-          python -m pytest -v --cov=mbuild --cov-report=xml --cov-append --cov-config=setup.cfg --color yes mbuild/tests/test_hoomd.py
-        if: runner.os == 'Linux'
-
       - name: Upload Coverage Report
         uses: codecov/codecov-action@v4
         with:
diff --git a/environment-dev-win.yml b/environment-dev-win.yml
index a95ac86a9..851d3e5b1 100644
--- a/environment-dev-win.yml
+++ b/environment-dev-win.yml
@@ -6,7 +6,7 @@ dependencies:
   - bump2version
   - codecov
   - ele
-  - gmso>=0.9.0
+  - gmso>=0.12.4
   - foyer>=0.11.0
   - garnett>=0.7.1
   - gsd>=2.9
@@ -20,7 +20,6 @@ dependencies:
   - parmed>=3.4.3
   - pip
   - pre-commit
-  - protobuf
   - py3Dmol
   - pycifrw
   - pytest
diff --git a/environment-dev.yml b/environment-dev.yml
index 0f151bc9e..42591309f 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -21,7 +21,6 @@ dependencies:
   - parmed>=3.4.3
   - pip
   - pre-commit
-  - protobuf
   - py3Dmol
   - pycifrw
   - pytest
diff --git a/mbuild/compound.py b/mbuild/compound.py
index 02eab970b..261f3b8ea 100644
--- a/mbuild/compound.py
+++ b/mbuild/compound.py
@@ -2007,7 +2007,6 @@ def _visualize_py3dmol(
             os.path.join(tmp_dir, "tmp.mol2"),
             include_ports=show_ports,
             overwrite=True,
-            parmed_kwargs={"infer_residues": False},
         )
 
         view = py3Dmol.view()
@@ -2979,15 +2978,9 @@ def save(
         self,
         filename,
         include_ports=False,
-        forcefield_name=None,
-        forcefield_files=None,
-        forcefield_debug=False,
         box=None,
         overwrite=False,
         residues=None,
-        combining_rule="lorentz",
-        foyer_kwargs=None,
-        parmed_kwargs=None,
         **kwargs,
     ):
         """Save the Compound to a file.
@@ -2997,22 +2990,11 @@ def save(
         filename : str
             Filesystem path in which to save the trajectory. The extension or
             prefix will be parsed and control the format. Supported extensions:
-            'hoomdxml', 'gsd', 'gro', 'top', 'lammps', 'lmp', 'mcf', 'pdb', 'xyz',
-            'json', 'mol2', 'sdf', 'psf'. See parmed/structure.py for more
-            information on savers.
+            'gsd', 'gro', 'top', 'mcf', 'pdb', 'xyz',
+            'json', 'mol2', 'sdf', 'psf'. See `mbuild.conversion.save()`
+            for more information about writer methods.
         include_ports : bool, optional, default=False
             Save ports contained within the compound.
-        forcefield_files : str, optional, default=None
-            Apply a forcefield to the output file using a forcefield provided
-            by the `foyer` package.
-        forcefield_name : str, optional, default=None
-            Apply a named forcefield to the output file using the `foyer`
-            package, e.g. 'oplsaa'. `Foyer forcefields
-            <https://github.com/mosdef-hub/foyer/tree/master/foyer/forcefields>`_
-        forcefield_debug : bool, optional, default=False
-            Choose verbosity level when applying a forcefield through `foyer`.
-            Specifically, when missing atom types in the forcefield xml file,
-            determine if the warning is condensed or verbose.
         box : mb.Box, optional, default=self.boundingbox (with buffer)
             Box information to be written to the output file. If 'None', a
             bounding box is used with 0.25nm buffers at each face to avoid
@@ -3022,50 +3004,25 @@ def save(
         residues : str of list of str
             Labels of residues in the Compound. Residues are assigned by
             checking against Compound.name.
-        combining_rule : str, optional, default='lorentz'
-            Specify the combining rule for nonbonded interactions. Only relevant
-            when the `foyer` package is used to apply a forcefield. Valid
-            options are 'lorentz' and 'geometric', specifying Lorentz-Berthelot
-            and geometric combining rules respectively.
-        foyer_kwargs : dict, optional, default=None
-            Keyword arguments to provide to `foyer.Forcefield.apply`.
-            Depending on the file extension these will be passed to either
-            `write_gsd`, `write_hoomdxml`, `write_lammpsdata`,
-            `write_mcf`, or `parmed.Structure.save`.
-            See `parmed structure documentation
-            <https://parmed.github.io/ParmEd/html/structobj/parmed.structure.Structure.html#parmed.structure.Structure.save>`_
-        parmed_kwargs : dict, optional, default=None
-            Keyword arguments to provide to :meth:`mbuild.Compound.to_parmed`
         **kwargs
+            See `mbuild.conversion.save()`.
             Depending on the file extension these will be passed to either
-            `write_gsd`, `write_hoomdxml`, `write_lammpsdata`, `write_mcf`, or
-            `parmed.Structure.save`.
+            Parmed or GMSO backend writers
             See https://parmed.github.io/ParmEd/html/structobj/parmed.structure.
-            Structure.html#parmed.structure.Structure.save
+            Structure.html#parmed.structure.Structure.save and
+            https://github.com/mosdef-hub/gmso/tree/main/gmso/formats
 
         Other Parameters
         ----------------
         ref_distance : float, optional, default=1.0
-            Normalization factor used when saving to .gsd and .hoomdxml formats
+            Normalization factor used when saving to the .gsd format
             for converting distance values to reduced units.
         ref_energy : float, optional, default=1.0
-            Normalization factor used when saving to .gsd and .hoomdxml formats
+            Normalization factor used when saving to the .gsd format
             for converting energy values to reduced units.
         ref_mass : float, optional, default=1.0
-            Normalization factor used when saving to .gsd and .hoomdxml formats
+            Normalization factor used when saving to the .gsd format
             for converting mass values to reduced units.
-        atom_style: str, default='full'
-            Defines the style of atoms to be saved in a LAMMPS data file. The
-            following atom styles are currently supported:
-            'full', 'atomic', 'charge', 'molecular'
-            See `LAMMPS atom style documentation
-            <https://lammps.sandia.gov/doc/atom_style.html>`_ for more
-            information.
-        unit_style: str, default='real'
-            Defines to unit style to be save in a LAMMPS data file.  Defaults
-            to 'real' units. Current styles are supported: 'real', 'lj'. See
-            `LAMMPS unit style documentation_
-            <https://lammps.sandia.gov/doc/units.html>`_ for more information.
 
         Notes
         -----
@@ -3074,29 +3031,24 @@ def save(
         * filename
         * include_ports
 
+        The savers used for each supported file type are:
+        GMSO: .gro, .gsd, .data, .xyz, .mcf, .top
+        Parmed: .mol2, .pdb, .prmtop, .cif, .crd
+        PyBel: .sdf
+
         See Also
         --------
-        conversion.save : Main saver logic
-        formats.gsdwrite.write_gsd : Write to GSD format
-        formats.hoomdxml.write_hoomdxml : Write to Hoomd XML format
-        formats.xyzwriter.write_xyz : Write to XYZ format
-        formats.lammpsdata.write_lammpsdata : Write to LAMMPS data format
-        formats.cassandramcf.write_mcf : Write to Cassandra MCF format
-        formats.json_formats.compound_to_json : Write to a json file
+        mbuild.conversion.save : Main saver logic
+        mbuild.formats.cassandramcf.write_mcf : Write to Cassandra MCF format
+        mbuild.formats.json_formats.compound_to_json : Write to a json file
         """
         conversion.save(
-            self,
-            filename,
-            include_ports,
-            forcefield_name,
-            forcefield_files,
-            forcefield_debug,
-            box,
-            overwrite,
-            residues,
-            combining_rule,
-            foyer_kwargs,
-            parmed_kwargs,
+            compound=self,
+            filename=filename,
+            include_ports=include_ports,
+            box=box,
+            overwrite=overwrite,
+            residues=residues,
             **kwargs,
         )
 
@@ -3246,6 +3198,21 @@ def to_gmso(self, **kwargs):
         """
         return conversion.to_gmso(self, **kwargs)
 
+    def to_hoomdsnapshot(self, **kwargs):
+        """Create a HOOMD-Blue snapshot from an mBuild Compound.
+
+        Parameters
+        ----------
+        compound : mb.Compound
+            The mb.Compound to be converted.
+
+        Returns
+        -------
+        snapshot : gsd.hoomd.Frame
+            HOOMD-Blue compatible topology.
+        """
+        return conversion.to_hoomdsnapshot(self, **kwargs)
+
     # Interface to Trajectory for reading/writing .pdb and .mol2 files.
     # -----------------------------------------------------------------
     def from_trajectory(
diff --git a/mbuild/conversion.py b/mbuild/conversion.py
index eb4361712..6d42091ab 100644
--- a/mbuild/conversion.py
+++ b/mbuild/conversion.py
@@ -7,6 +7,7 @@
 from pathlib import Path
 from warnings import warn
 
+import gmso
 import numpy as np
 import parmed as pmd
 from ele import (
@@ -19,12 +20,7 @@
 import mbuild as mb
 from mbuild.box import Box
 from mbuild.exceptions import MBuildError
-from mbuild.formats.gsdwriter import write_gsd
-from mbuild.formats.hoomdxml import write_hoomdxml
 from mbuild.formats.json_formats import compound_from_json, compound_to_json
-from mbuild.formats.lammpsdata import write_lammpsdata
-from mbuild.formats.par_writer import write_par
-from mbuild.formats.xyz import read_xyz, write_xyz
 from mbuild.utils.io import (
     has_gmso,
     has_mdtraj,
@@ -59,9 +55,9 @@ def load(
     Parameters
     ----------
     filename_or_object : str, mdtraj.Trajectory, parmed.Structure,
-        mbuild.Compound, pybel.Molecule,
+        mbuild.Compound, pybel.Molecule.
         Name of the file or topology from which to load atom and bond
-        information.
+        information, or a valid SMILES string.
     relative_to_module : str, optional, default=None
         Instead of looking in the current working directory, look for the file
         where this module is defined. This is typically used in Compound
@@ -80,7 +76,7 @@ def load(
     smiles: bool, optional, default=False
         Use RDKit or OpenBabel to parse filename as a SMILES string or file
         containing a SMILES string. If this is set to True, `rdkit` is the
-        default backend.
+        default backend and `filename_or_object` should be the SMILES string.
     infer_hierarchy : bool, optional, default=True
         If True, infer hierarchy from chains and residues
     ignore_box_warn : bool, optional, default=False
@@ -172,21 +168,20 @@ def load_object(
     mb.Compound
     """
     # Create type_dict type -> loading method
-    # Will need to add a gmso method soon
-    type_dict = {pmd.Structure: from_parmed}
-    if has_gmso:
-        gmso = import_("gmso")
-        type_dict.update({gmso.Topology: from_gmso})
+    type_dict = {
+        pmd.Structure: from_parmed,
+        gmso.Topology: from_gmso,
+    }
 
     if has_openbabel:
         pybel = import_("pybel")
         type_dict.update({pybel.Molecule: from_pybel})
-
     if has_mdtraj:
         md = import_("mdtraj")
         type_dict.update({md.Trajectory: from_trajectory})
 
     # Check if the given object is an mb.Compound
+    # TODO 1.0: What is the use case for this?
     if isinstance(obj, mb.Compound):
         if not compound:
             warn("Given object is already an mb.Compound, doing nothing.")
@@ -210,7 +205,7 @@ def load_object(
             return compound
 
     # If nothing is return raise an error
-    raise ValueError(f"Object of type {type(obj).__name__} is not supported")
+    raise ValueError(f"Object of type {type(obj).__name__} is not supported.")
 
 
 def load_smiles(
@@ -374,13 +369,11 @@ def load_file(
     # Need to come up with a different dict structure
     default_backends = {
         ".json": "internal",
-        ".xyz": "internal",
+        ".xyz": "gmso",
         ".sdf": "pybel",
-        ".hoomdxml": "mdtraj",
-        ".mol2": "mdtraj",
+        ".mol2": "gmso",
         ".pdb": "mdtraj",
     }
-
     # Handle mbuild *.py files containing a class that wraps a structure file
     # in its own folder. E.g., you build a system from ~/foo.py and it imports
     # from ~/bar/baz.py where baz.py loads ~/bar/baz.pdb.
@@ -391,11 +384,9 @@ def load_file(
     extension = Path(filename).suffix
 
     if not backend:
-        try:
-            # Try matching backend based on extension
+        try:  # Try matching backend based on extension
             backend = default_backends[extension]
-        except KeyError:
-            # Else use default backend
+        except KeyError:  # Else use default backend
             backend = "mdtraj" if has_mdtraj else "parmed"
 
     # First check internal readers
@@ -406,29 +397,8 @@ def load_file(
             compound = compound_from_json(filename)
             return compound
         # Handle xyz file
-        if extension == ".xyz" and not "top" in kwargs:
-            if coords_only:
-                tmp = read_xyz(filename)
-                if tmp.n_particles != compound.n_particles:
-                    raise ValueError(
-                        f"Number of atoms in {filename}"
-                        f"does not match {compound}"
-                    )
-                ref_and_compound = zip(
-                    tmp._particles(include_ports=False),
-                    compound.particles(include_ports=False),
-                )
-                for ref_particle, particle in ref_and_compound:
-                    particle.pos = ref_particle.pos
-            else:
-                compound = read_xyz(filename, compound=compound)
-        elif extension == ".xyz" and "top" in kwargs:
-            backend = "mdtraj"
-
     # Then gmso reader
     if backend == "gmso":
-        gmso = import_("gmso")
-
         top = gmso.Topology.load(filename=filename, **kwargs)
         compound = from_gmso(
             topology=top,
@@ -443,8 +413,7 @@ def load_file(
         if extension == ".sdf":
             pybel_mol = pybel.readfile("sdf", filename)
             # pybel returns a generator, so we grab the first molecule of a
-            # list of len 1.
-            # Raise ValueError if there are more molecules
+            # list of len 1. Raise ValueError if there are more molecules
             pybel_mol = [i for i in pybel_mol]
             if len(pybel_mol) == 1:
                 compound = from_pybel(
@@ -495,7 +464,6 @@ def load_file(
     # by the corresponding backend
     if rigid:
         compound.label_rigid_bodies()
-
     return compound
 
 
@@ -967,15 +935,9 @@ def save(
     compound,
     filename,
     include_ports=False,
-    forcefield_name=None,
-    forcefield_files=None,
-    forcefield_debug=False,
     box=None,
     overwrite=False,
     residues=None,
-    combining_rule=None,
-    foyer_kwargs=None,
-    parmed_kwargs=None,
     **kwargs,
 ):
     """Save the Compound to a file.
@@ -987,22 +949,10 @@ def save(
     filename : str
         Filesystem path in which to save the trajectory. The extension or prefix
         will be parsed and control the format. Supported extensions are:
-        'hoomdxml', 'gsd', 'gro', 'top', 'lammps', 'lmp', 'mcf', 'xyz', 'pdb',
-        'sdf', 'mol2', 'psf'. See parmed/structure.py for more information on
-        savers.
+        'gsd', 'gro', 'top', 'mcf', 'xyz', 'pdb', 'sdf', 'mol2', 'psf'.
+        See parmed/structure.py for more information on savers.
     include_ports : bool, optional, default=False
         Save ports contained within the compound.
-    forcefield_files : str, optional, default=None
-        Apply a forcefield to the output file using a forcefield provided by the
-        `foyer` package.
-    forcefield_name : str, optional, default=None
-        Apply a named forcefield to the output file using the `foyer` package,
-        e.g. 'oplsaa'. Forcefields listed here:
-        https://github.com/mosdef-hub/foyer/tree/master/foyer/forcefields
-    forcefield_debug : bool, optional, default=False
-        Choose level of verbosity when applying a forcefield through `foyer`.
-        Specifically, when missing atom types in the forcefield xml file,
-        determine if the warning is condensed or verbose.
     box : mb.Box, optional, default=compound.boundingbox (with buffer)
         Box information to be written to the output file. If 'None', a bounding
         box is used with 0.25nm buffers at each face to avoid overlapping atoms.
@@ -1011,44 +961,24 @@ def save(
     residues : str of list of str
         Labels of residues in the Compound. Residues are assigned by checking
         against Compound.name.
-    combining_rule : str or None, optional, default=None
-        Specify the combining rule for nonbonded interactions. Only relevant
-        when the `foyer` package is used to apply a forcefield. Valid options
-        are 'lorentz' and 'geometric' and None, specifying Lorentz-Berthelot and
-        geometric combining rules respectively. If None is provided, the combining_rule
-        specified from the forcefield will be used.
-    foyer_kwargs : dict, optional, default=None
-        Keyword arguments to provide to `foyer.Forcefield.apply`.
-    parmed_kwargs : dict, optional, default=None
-        Keyword arguments to provide to :meth:`mbuild.Compound.to_parmed`
-    **kwargs
+    **kwargs : dict, optional
         Depending on the file extension these will be passed to either
-        `write_gsd`, `write_hoomdxml`, `write_lammpsdata`, `write_mcf`, or
-        `parmed.Structure.save`.
+        GMSO's intenral writers or `parmed.Structure.save`.
+        See https://github.com/mosdef-hub/gmso/tree/main/gmso/formats
         See https://parmed.github.io/ParmEd/html/structobj/parmed.structure.
         Structure.html#parmed.structure.Structure.save
 
     Other Parameters
     ----------------
     ref_distance : float, optional, default=1.0
-        Normalization factor used when saving to .gsd and .hoomdxml formats for
+        Normalization factor used when saving to .gsd formats for
         converting distance values to reduced units.
     ref_energy : float, optional, default=1.0
-        Normalization factor used when saving to .gsd and .hoomdxml formats for
+        Normalization factor used when saving to .gsd formats for
         converting energy values to reduced units.
     ref_mass : float, optional, default=1.0
-        Normalization factor used when saving to .gsd and .hoomdxml formats for
+        Normalization factor used when saving to .gsd formats for
         converting mass values to reduced units.
-    atom_style: str, default='full'
-        Defines the style of atoms to be saved in a LAMMPS data file. The
-        following atom styles are currently supported: 'full', 'atomic',
-        'charge', 'molecular'. See http://lammps.sandia.gov/doc/atom_style.html
-        for more information on atom styles.
-    unit_style: str, default='real'
-        Defines to unit style to be save in a LAMMPS data file.  Defaults to
-        'real' units. Current styles are supported: 'real', 'lj'. See
-        https://lammps.sandia.gov/doc/99/units.html for more information on
-        unit styles
 
     Notes
     -----
@@ -1058,15 +988,19 @@ def save(
 
     See Also
     --------
-    formats.gsdwriter.write_gsd : Write to GSD format
-    formats.hoomdxml.write_hoomdxml : Write to Hoomd XML format
-    formats.xyzwriter.write_xyz : Write to XYZ format
-    formats.lammpsdata.write_lammpsdata : Write to LAMMPS data format
     formats.cassandramcf.write_mcf : Write to Cassandra MCF format
     formats.json_formats.compound_to_json : Write to a json file
     """
-    extension = os.path.splitext(filename)[-1]
+    if os.path.exists(filename) and not overwrite:
+        raise IOError(f"{filename} exists; not overwriting")
+    if compound.charge:
+        if round(compound.charge, 4) != 0.0:
+            warn(
+                f"System is not charge neutral. Total charge is {compound.charge}."
+            )
 
+    extension = os.path.splitext(filename)[-1]
+    # Keep json stuff with internal mbuild method
     if extension == ".json":
         compound_to_json(
             compound, file_path=filename, include_ports=include_ports
@@ -1075,98 +1009,60 @@ def save(
 
     # Savers supported by mbuild.formats
     savers = {
-        ".hoomdxml": write_hoomdxml,
-        ".gsd": write_gsd,
-        ".xyz": write_xyz,
-        ".lammps": write_lammpsdata,
-        ".lmp": write_lammpsdata,
-        ".par": write_par,
+        ".gro": save_in_gmso,
+        ".gsd": save_in_gmso,
+        ".data": save_in_gmso,
+        ".xyz": save_in_gmso,
+        # ".mol2": save_in_gmso,
+        ".mcf": save_in_gmso,
+        ".top": save_in_gmso,
     }
-    if has_networkx:
-        from mbuild.formats.cassandramcf import write_mcf
-
-        savers.update({".mcf": write_mcf})
 
     try:
         saver = savers[extension]
     except KeyError:
         saver = None
-
-    if os.path.exists(filename) and not overwrite:
-        raise IOError(f"{filename} exists; not overwriting")
-
-    if not parmed_kwargs:
-        parmed_kwargs = {}
-    structure = compound.to_parmed(
-        box=box,
-        residues=residues,
-        include_ports=include_ports,
-        **parmed_kwargs,
-    )
-    # Apply a force field with foyer if specified
-    if forcefield_name or forcefield_files:
-        foyer = import_("foyer")
-        ff = foyer.Forcefield(
-            forcefield_files=forcefield_files,
-            name=forcefield_name,
-            debug=forcefield_debug,
-        )
-
-        if not foyer_kwargs:
-            foyer_kwargs = {}
-        structure = ff.apply(structure, **foyer_kwargs)
-
-        if combining_rule is None:
-            combining_rule = structure.combining_rule
-
-        if structure.combining_rule != combining_rule:
-            warn(
-                f"Overwriting forcefield-specified combining rule ({combining_rule})"
-                f"to new combining rule ({combining_rule})."
-                "This can cause inconsistent between the 1-4 pair interactions,"
-                "calculated in foyer, and the new combining rule."
-                "Consider directly changing the metadata of the Forcefield."
-            )
-
-        structure.combining_rule = combining_rule
-        if structure.__dict__.get("defaults"):
-            structure.defaults.comb_rule = (
-                2 if combining_rule == "lorentz" else 3
-            )
-
-    total_charge = sum([atom.charge for atom in structure])
-    if round(total_charge, 4) != 0.0:
-        warn(f"System is not charge neutral. Total charge is {total_charge}.")
-
-    # Provide a warning if rigid_ids are not sequential from 0
-    if compound.contains_rigid:
-        unique_rigid_ids = sorted(
-            set([p.rigid_id for p in compound.rigid_particles()])
-        )
-        if max(unique_rigid_ids) != len(unique_rigid_ids) - 1:
-            warn("Unique rigid body IDs are not sequential starting from zero.")
-
     if saver:  # mBuild supported saver.
-        if extension in [".gsd", ".hoomdxml"]:
-            kwargs["rigid_bodies"] = [p.rigid_id for p in compound.particles()]
-        saver(filename=filename, structure=structure, **kwargs)
+        # Calling save_in_gmso
+        saver(
+            filename=filename,
+            compound=compound,
+            box=box,
+            overwrite=overwrite,
+            **kwargs,
+        )
 
     elif extension == ".sdf":
         pybel = import_("pybel")
-        new_compound = mb.Compound()
-        # Convert pmd.Structure to mb.Compound
-        new_compound.from_parmed(structure)
-        # Convert mb.Compound to pybel molecule
-        pybel_molecule = new_compound.to_pybel()
+        pybel_molecule = compound.to_pybel()
         # Write out pybel molecule to SDF file
         output_sdf = pybel.Outputfile("sdf", filename, overwrite=overwrite)
         output_sdf.write(pybel_molecule)
         output_sdf.close()
-
     else:  # ParmEd supported saver.
+        structure = compound.to_parmed()
         structure.save(filename, overwrite=overwrite, **kwargs)
 
 
+def save_in_gmso(compound, filename, box, overwrite, **kwargs):
+    """Convert to GMSO, call GMSO internal writers.
+
+    Parameters
+    ----------
+    compound : mbuild.compound.Compound, required.
+        The mBuild compound to convert to a GSMO topology.
+    box : mbuild.box.Box, required
+        The mBuild box to be converted to a gmso box, if different that compound.box
+    overwrite : bool, required
+        If `True` overwrite the file if it already exists.
+    **kwargs
+        Keyword arguments used in GMSO's savers.
+        See https://github.com/mosdef-hub/gmso/tree/main/gmso/formats
+    """
+    gmso_top = to_gmso(compound=compound, box=box)
+    gmso_top.save(filename=filename, overwrite=overwrite, **kwargs)
+
+
 def catalog_bondgraph_type(compound, bond_graph=None):
     """Identify type of subgraph found at this stage of the compound.
 
@@ -1298,6 +1194,69 @@ def pull_residues(
     return residuesList
 
 
+def to_hoomdsnapshot(
+    compound,
+    identify_connections=True,
+    ref_distance=1.0,
+    ref_mass=1.0,
+    rigid_bodies=None,
+    shift_coords=True,
+    write_special_pairs=True,
+    **kwargs,
+):
+    """Output a gsd.hoomd.Frame (HOOMD-Blue topology format).
+
+    Parameters
+    ----------
+    compound : mb.Compound
+        mBuild compound to save to the GSD format.
+    identify_connections : bool
+        If `True`, then infer angles and dihedrals in the topology.
+    ref_distance : float, optional, default=1.0
+        Reference distance for conversion to reduced units
+    ref_mass : float, optional, default=1.0
+        Reference mass for conversion to reduced units
+    rigid_bodies : list of int, optional, default=None
+        List of rigid body information. An integer value is required for each
+        atom corresponding to the index of the rigid body the particle is to be
+        associated with. A value of None indicates the atom is not part of a
+        rigid body.
+    shift_coords : bool, optional, default=True
+        Shift coordinates from (0, L) to (-L/2, L/2) if necessary.
+    write_special_pairs : bool, optional, default=True
+        Writes out special pair information necessary to correctly use the OPLS
+        fudged 1,4 interactions in HOOMD-Blue.
+
+    Notes
+    -----
+    See gmso.external.convert_hoomd on how to make sure the GSD file contains
+    forcefield information (e.g. atom types, angle types, etc.).
+    """
+    import unyt as u
+    from gmso.external import from_mbuild, to_gsd_snapshot
+
+    gmso_top = from_mbuild(compound=compound)
+
+    if identify_connections:
+        gmso_top.identify_connections()
+
+    base_units = {
+        "length": ref_distance * u.Unit("nm"),
+        "mass": ref_mass * u.Unit("amu"),
+        "energy": 1 * u.Unit("kJ/mol"),
+    }
+
+    snapshot, refs = to_gsd_snapshot(
+        top=gmso_top,
+        base_units=base_units,
+        rigid_bodies=rigid_bodies,
+        shift_coords=shift_coords,
+        parse_special_pairs=write_special_pairs,
+        auto_scale=False,
+    )
+    return snapshot
+
+
 def to_parmed(
     compound,
     box=None,
@@ -1811,7 +1770,6 @@ def to_rdkit(compound):
 
     for particle in compound.particles():
         temp_atom = Chem.Atom(particle.element.atomic_number)
-
         # this next line is necessary to prevent rdkit from adding hydrogens
         # this will also set the label to be the element with particle index
         temp_atom.SetProp(
@@ -1924,7 +1882,14 @@ def _iterate_children(compound, nodes, edges, names_only=False):
     return nodes, edges
 
 
-def to_gmso(compound, box=None, **kwargs):
+def to_gmso(
+    compound,
+    box=None,
+    parse_label=True,
+    custom_groups=None,
+    infer_elements=False,
+    **kwargs,
+):
     """Create a GMSO Topology from a mBuild Compound.
 
     Parameters
@@ -1941,7 +1906,13 @@ def to_gmso(compound, box=None, **kwargs):
     """
     from gmso.external.convert_mbuild import from_mbuild
 
-    return from_mbuild(compound, box=None, **kwargs)
+    return from_mbuild(
+        compound=compound,
+        box=box,
+        parse_label=parse_label,
+        custom_groups=custom_groups,
+        infer_elements=infer_elements,
+    )
 
 
 def to_intermol(compound, molecule_types=None):  # pragma: no cover
diff --git a/mbuild/formats/gsdwriter.py b/mbuild/formats/gsdwriter.py
deleted file mode 100644
index 2de3d5cd3..000000000
--- a/mbuild/formats/gsdwriter.py
+++ /dev/null
@@ -1,301 +0,0 @@
-"""GSD format.
-
-https://gsd.readthedocs.io/en/stable/
-"""
-
-import numpy as np
-
-from mbuild.utils.geometry import coord_shift
-from mbuild.utils.io import import_
-from mbuild.utils.sorting import natural_sort
-
-__all__ = ["write_gsd"]
-
-
-def write_gsd(
-    structure,
-    filename,
-    ref_distance=1.0,
-    ref_mass=1.0,
-    ref_energy=1.0,
-    rigid_bodies=None,
-    shift_coords=True,
-    write_special_pairs=True,
-    **kwargs
-):
-    """Output a GSD file (HOOMD v2 default data format).
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd Structure object
-    filename : str
-        Path of the output file.
-    ref_distance : float, optional, default=1.0
-        Reference distance for conversion to reduced units
-    ref_mass : float, optional, default=1.0
-        Reference mass for conversion to reduced units
-    ref_energy : float, optional, default=1.0
-        Reference energy for conversion to reduced units
-    rigid_bodies : list of int, optional, default=None
-        List of rigid body information. An integer value is required for each
-        atom corresponding to the index of the rigid body the particle is to be
-        associated with. A value of None indicates the atom is not part of a
-        rigid body.
-    shift_coords : bool, optional, default=True
-        Shift coordinates from (0, L) to (-L/2, L/2) if necessary.
-    write_special_pairs : bool, optional, default=True
-        Writes out special pair information necessary to correctly use the OPLS
-        fudged 1,4 interactions
-        in HOOMD.
-
-    Notes
-    -----
-    Force field parameters are not written to the GSD file and must be included
-    manually into a HOOMD input script.
-    """
-    import_("gsd")
-    import gsd.hoomd
-
-    xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms])
-    if shift_coords:
-        xyz = coord_shift(xyz, structure.box[:3])
-
-    gsd_snapshot = gsd.hoomd.Frame()
-
-    gsd_snapshot.configuration.step = 0
-    gsd_snapshot.configuration.dimensions = 3
-
-    # Write box information
-    if np.allclose(structure.box[3:6], np.array([90, 90, 90])):
-        gsd_snapshot.configuration.box = np.hstack(
-            (structure.box[:3] / ref_distance, np.zeros(3))
-        )
-    else:
-        a, b, c = structure.box[0:3] / ref_distance
-        alpha, beta, gamma = np.radians(structure.box[3:6])
-
-        lx = a
-        xy = b * np.cos(gamma)
-        xz = c * np.cos(beta)
-        ly = np.sqrt(b**2 - xy**2)
-        yz = (b * c * np.cos(alpha) - xy * xz) / ly
-        lz = np.sqrt(c**2 - xz**2 - yz**2)
-
-        gsd_snapshot.configuration.box = np.array([lx, ly, lz, xy, xz, yz])
-
-    _write_particle_information(
-        gsd_snapshot,
-        structure,
-        xyz,
-        ref_distance,
-        ref_mass,
-        ref_energy,
-        rigid_bodies,
-    )
-    if write_special_pairs:
-        _write_pair_information(gsd_snapshot, structure)
-    if structure.bonds:
-        _write_bond_information(gsd_snapshot, structure)
-    if structure.angles:
-        _write_angle_information(gsd_snapshot, structure)
-    if structure.rb_torsions:
-        _write_dihedral_information(gsd_snapshot, structure)
-
-    with gsd.hoomd.open(filename, mode="w") as gsd_file:
-        gsd_file.append(gsd_snapshot)
-
-
-def _write_particle_information(
-    gsd_snapshot,
-    structure,
-    xyz,
-    ref_distance,
-    ref_mass,
-    ref_energy,
-    rigid_bodies,
-):
-    """Write out the particle information."""
-    gsd_snapshot.particles.N = len(structure.atoms)
-    gsd_snapshot.particles.position = xyz / ref_distance
-
-    types = [a.name if a.type == "" else a.type for a in structure.atoms]
-
-    unique_types = list(set(types))
-    unique_types.sort(key=natural_sort)
-    gsd_snapshot.particles.types = unique_types
-
-    typeids = np.array([unique_types.index(t) for t in types])
-    gsd_snapshot.particles.typeid = typeids
-
-    masses = np.array([atom.mass for atom in structure.atoms])
-    masses[masses == 0] = 1.0
-    gsd_snapshot.particles.mass = masses / ref_mass
-
-    charges = np.array([atom.charge for atom in structure.atoms])
-    e0 = 2.396452e-04
-    """
-    Permittivity of free space = 2.396452e-04 e^2/((kcal/mol) Angstrom),
-    where e is the elementary charge
-    """
-    charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5
-    gsd_snapshot.particles.charge = charges / charge_factor
-
-    if rigid_bodies:
-        rigid_bodies = [-1 if body is None else body for body in rigid_bodies]
-    gsd_snapshot.particles.body = rigid_bodies
-
-
-def _write_pair_information(gsd_snapshot, structure):
-    """Write the special pairs in the system.
-
-    Parameters
-    ----------
-    gsd_snapshot : gsd.snapshot,
-        The file object of the GSD file being written
-    structure : parmed.Structure
-        Parmed structure object holding system information
-    """
-    pair_types = []
-    pair_typeid = []
-    pairs = []
-    for ai in structure.atoms:
-        for aj in ai.dihedral_partners:
-            # make sure we don't double add
-            if ai.idx > aj.idx:
-                ps = "-".join(sorted([ai.type, aj.type], key=natural_sort))
-                if ps not in pair_types:
-                    pair_types.append(ps)
-                pair_typeid.append(pair_types.index(ps))
-                pairs.append((ai.idx, aj.idx))
-    gsd_snapshot.pairs.types = pair_types
-    gsd_snapshot.pairs.typeid = pair_typeid
-    gsd_snapshot.pairs.group = pairs
-    gsd_snapshot.pairs.N = len(pairs)
-
-
-def _write_bond_information(gsd_snapshot, structure):
-    """Write the bonds in the system.
-
-    Parameters
-    ----------
-    gsd_snapshot :
-        The file object of the GSD file being written
-    structure : parmed.Structure
-        Parmed structure object holding system information
-    """
-    gsd_snapshot.bonds.N = len(structure.bonds)
-
-    unique_bond_types = set()
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        if t1 == "" or t2 == "":
-            t1, t2 = bond.atom1.name, bond.atom2.name
-        t1, t2 = sorted([t1, t2], key=natural_sort)
-        try:
-            bond_type = "-".join((t1, t2))
-        except AttributeError:  # no forcefield applied, bond.type is None
-            bond_type = ("-".join((t1, t2)), 0.0, 0.0)
-        unique_bond_types.add(bond_type)
-    unique_bond_types = sorted(list(unique_bond_types), key=natural_sort)
-    gsd_snapshot.bonds.types = unique_bond_types
-
-    bond_typeids = []
-    bond_groups = []
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        if t1 == "" or t2 == "":
-            t1, t2 = bond.atom1.name, bond.atom2.name
-        t1, t2 = sorted([t1, t2], key=natural_sort)
-        try:
-            bond_type = "-".join((t1, t2))
-        except AttributeError:  # no forcefield applied, bond.type is None
-            bond_type = ("-".join((t1, t2)), 0.0, 0.0)
-        bond_typeids.append(unique_bond_types.index(bond_type))
-        bond_groups.append((bond.atom1.idx, bond.atom2.idx))
-
-    gsd_snapshot.bonds.typeid = bond_typeids
-    gsd_snapshot.bonds.group = bond_groups
-
-
-def _write_angle_information(gsd_snapshot, structure):
-    """Write the angles in the system.
-
-    Parameters
-    ----------
-    gsd_snapshot :
-        The file object of the GSD file being written
-    structure : parmed.Structure
-        Parmed structure object holding system information
-    """
-    gsd_snapshot.angles.N = len(structure.angles)
-
-    unique_angle_types = set()
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3], key=natural_sort)
-        angle_type = "-".join((t1, t2, t3))
-        unique_angle_types.add(angle_type)
-    unique_angle_types = sorted(list(unique_angle_types), key=natural_sort)
-    gsd_snapshot.angles.types = unique_angle_types
-
-    angle_typeids = []
-    angle_groups = []
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3], key=natural_sort)
-        angle_type = "-".join((t1, t2, t3))
-        angle_typeids.append(unique_angle_types.index(angle_type))
-        angle_groups.append((angle.atom1.idx, angle.atom2.idx, angle.atom3.idx))
-
-    gsd_snapshot.angles.typeid = angle_typeids
-    gsd_snapshot.angles.group = angle_groups
-
-
-def _write_dihedral_information(gsd_snapshot, structure):
-    """Write the dihedrals in the system.
-
-    Parameters
-    ----------
-    gsd_snapshot :
-        The file object of the GSD file being written
-    structure : parmed.Structure
-        Parmed structure object holding system information
-    """
-    gsd_snapshot.dihedrals.N = len(structure.rb_torsions)
-
-    unique_dihedral_types = set()
-    for dihedral in structure.rb_torsions:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        unique_dihedral_types.add(dihedral_type)
-    unique_dihedral_types = sorted(
-        list(unique_dihedral_types), key=natural_sort
-    )
-    gsd_snapshot.dihedrals.types = unique_dihedral_types
-
-    dihedral_typeids = []
-    dihedral_groups = []
-    for dihedral in structure.rb_torsions:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        dihedral_typeids.append(unique_dihedral_types.index(dihedral_type))
-        dihedral_groups.append(
-            (
-                dihedral.atom1.idx,
-                dihedral.atom2.idx,
-                dihedral.atom3.idx,
-                dihedral.atom4.idx,
-            )
-        )
-
-    gsd_snapshot.dihedrals.typeid = dihedral_typeids
-    gsd_snapshot.dihedrals.group = dihedral_groups
diff --git a/mbuild/formats/hoomd_forcefield.py b/mbuild/formats/hoomd_forcefield.py
deleted file mode 100644
index 265947785..000000000
--- a/mbuild/formats/hoomd_forcefield.py
+++ /dev/null
@@ -1,466 +0,0 @@
-"""HOOMD v3 forcefield format."""
-
-import itertools
-import operator
-import warnings
-from collections import namedtuple
-
-import numpy as np
-import parmed as pmd
-
-import mbuild as mb
-from mbuild.utils.conversion import RB_to_OPLS
-from mbuild.utils.decorators import deprecated
-from mbuild.utils.io import import_
-from mbuild.utils.sorting import natural_sort
-
-from .hoomd_snapshot import _get_hoomd_version, to_hoomdsnapshot
-
-hoomd = import_("hoomd")
-
-
-dep_msg = """
-Support for Hoomd-Blue writers will be removed in mBuild 1.0.
-See GMSO (https://github.com/mosdef-hub/gmso/tree/main/gmso/external/convert_hoomd) for
-Hoomd-Blue 3.x and 4.x format support.
-To convert to GMSO, use the method `Compound.to_gmso()`.
-"""
-print(dep_msg)
-
-
-@deprecated(dep_msg)
-def create_hoomd_forcefield(
-    structure,
-    r_cut,
-    ref_distance=1.0,
-    ref_mass=1.0,
-    ref_energy=1.0,
-    auto_scale=False,
-    nlist_buffer=0.4,
-    snapshot_kwargs={},
-    pppm_kwargs={"Nx": 8, "Ny": 8, "Nz": 8, "order": 4},
-    init_snap=None,
-):
-    """Convert a parametrized pmd.Structure to a HOOMD snapshot and forces.
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd Structure object
-    r_cut : float
-        Cutoff radius in simulation units
-    ref_distance : float, optional, default=1.0
-        Reference distance for unit conversion ((Angstrom) / (desired units))
-    ref_mass : float, optional, default=1.0
-        Reference mass for unit conversion ((Dalton) / (desired units))
-    ref_energy : float, optional, default=1.0
-        Reference energy for unit conversion ((kcal/mol) / (desired units))
-    auto_scale : bool, optional, default=False
-        Scale to reduced units by automatically using the largest sigma value
-        as ref_distance, largest mass value as ref_mass, and largest epsilon
-        value as ref_energy
-    nlist_buffer : float, optional, default=True
-        buffer argument to pass to hoomd.md.nlist.Cell
-    snapshot_kwargs : dict
-        Keyword arguments to pass to to_hoomdsnapshot
-    pppm_kwargs : dict
-        Keyword arguments to pass to hoomd.md.long_range.pppm.make_pppm_coulomb_forces
-    init_snap : hoomd.Snapshot, optional, default=None
-        Initial snapshot to which to add the ParmEd structure object
-        (useful for rigid bodies)
-
-    Returns
-    -------
-    hoomd_snapshot : hoomd.Snapshot
-        HOOMD snapshot object to initialize the simulation
-    hoomd_forcefield : list[hoomd.md.force.Force]
-        List of hoomd force computes created during conversion
-    ReferenceValues : namedtuple
-        Values used in scaling
-
-    Notes
-    -----
-    If you pass a non-parametrized pmd.Structure, you will not have
-    angle, dihedral, or force field information. You may be better off
-    creating a hoomd.Snapshot.
-    See mbuild.formats.hoomd_snapshot.to_hoomdsnapshot()
-
-    About units: This method operates on a Parmed.Structure object
-    where the units used differ from those used in mBuild and Foyer
-    which may have been used when creating the typed Parmed.Structure.
-
-    The default units used when writing out the HOOMD Snapshot are:
-    Distance (Angstrom)
-    Mass (Dalton)
-    Energy (kcal/mol)
-
-    If you wish to convert this unit system to another, you can use the
-    reference parameters (ref_distance, ref_mass, ref_energy).
-    The values used here should be expected to convert from the Parmed
-    Structure units (above) to your desired units.
-    The Parmed.Structure values are divided by the reference values.
-
-    If you wish to used a reduced unit system, set auto_scale = True.
-    When auto_scale is True, the reference parameters won't be used.
-
-    Examples
-    --------
-    To convert the energy units from kcal/mol to kj/mol:
-        use ref_energy = 0.2390057 (kcal/kj)
-
-    To convert the distance units from Angstrom to nm:
-        use ref_distance = 10 (angstroms/nm)
-
-    To use a reduced unit system, where mass, sigma, and epsilon are
-    scaled by the largest value of each:
-        use auto_scale = True, ref_distance = ref_energy = ref_mass = 1
-
-    """
-    if isinstance(structure, mb.Compound):
-        raise ValueError(
-            "You passed mb.Compound to create_hoomd_simulation, there will be "
-            "no angles, dihedrals, or force field parameters. Please use "
-            "hoomd_snapshot.to_hoomdsnapshot to create a hoomd.Snapshot, then "
-            "create your own hoomd context and pass your hoomd.Snapshot to "
-            "hoomd.init.read_snapshot()"
-        )
-    elif not isinstance(structure, pmd.Structure):
-        raise ValueError(
-            "Please pass a parmed.Structure to create_hoomd_simulation"
-        )
-
-    hoomd_version = _get_hoomd_version()
-    if hoomd_version.major < 3:
-        raise RuntimeError(
-            "Unsupported HOOMD-blue version:", str(hoomd_version)
-        )
-
-    hoomd_forcefield = []
-
-    if auto_scale:
-        if not all([i == 1 for i in (ref_distance, ref_energy, ref_mass)]):
-            warnings.warn(
-                "Autoscale option selected--provided reference values will not "
-                "be used."
-            )
-
-        pair_coeffs = list(
-            set((a.type, a.epsilon, a.sigma) for a in structure.atoms)
-        )
-
-        ref_mass = max([atom.mass for atom in structure.atoms])
-        ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1]
-        ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2]
-
-    ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"])
-    ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy)
-
-    snapshot, _ = to_hoomdsnapshot(
-        structure,
-        ref_distance=ref_distance,
-        ref_mass=ref_mass,
-        ref_energy=ref_energy,
-        **snapshot_kwargs,
-        hoomd_snapshot=init_snap,
-    )
-
-    nl = hoomd.md.nlist.Cell(exclusions=["bond", "1-3"], buffer=nlist_buffer)
-
-    if structure.atoms[0].type != "":
-        print("Processing LJ and QQ")
-        lj = _init_hoomd_lj(
-            structure,
-            nl,
-            r_cut,
-            ref_distance=ref_distance,
-            ref_energy=ref_energy,
-        )
-        qq = _init_hoomd_qq(structure, nl, snapshot, r_cut, **pppm_kwargs)
-        hoomd_forcefield.append(lj)
-        if qq is not None:
-            hoomd_forcefield.extend(qq)
-    if structure.adjusts:
-        print("Processing 1-4 interactions, adjusting neighborlist exclusions")
-        lj_14, qq_14 = _init_hoomd_14_pairs(
-            structure,
-            nl,
-            snapshot,
-            r_cut,
-            ref_distance=ref_distance,
-            ref_energy=ref_energy,
-        )
-        hoomd_forcefield.append(lj_14)
-        hoomd_forcefield.append(qq_14)
-    if structure.bond_types:
-        print("Processing harmonic bonds")
-        harmonic_bond = _init_hoomd_bonds(
-            structure, ref_distance=ref_distance, ref_energy=ref_energy
-        )
-        hoomd_forcefield.append(harmonic_bond)
-    if structure.angle_types:
-        print("Processing harmonic angles")
-        harmonic_angle = _init_hoomd_angles(structure, ref_energy=ref_energy)
-        hoomd_forcefield.append(harmonic_angle)
-    if structure.dihedral_types:
-        print("Processing periodic torsions")
-        periodic_torsions = _init_hoomd_dihedrals(
-            structure, ref_energy=ref_energy
-        )
-        hoomd_forcefield.append(periodic_torsions)
-    if structure.rb_torsion_types:
-        print("Processing RB torsions")
-        rb_torsions = _init_hoomd_rb_torsions(structure, ref_energy=ref_energy)
-        hoomd_forcefield.append(rb_torsions)
-
-    return snapshot, hoomd_forcefield, ref_values
-
-
-def _init_hoomd_lj(structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0):
-    """LJ parameters."""
-    # Identify the unique atom types before setting
-    atom_type_params = {}
-    for atom in structure.atoms:
-        if atom.type not in atom_type_params:
-            atom_type_params[atom.type] = atom.atom_type
-
-    # Set the hoomd parameters for self-interactions
-    lj = hoomd.md.pair.LJ(nlist=nl)
-    for name, atom_type in atom_type_params.items():
-        lj.params[(name, name)] = dict(
-            sigma=atom_type.sigma / ref_distance,
-            epsilon=atom_type.epsilon / ref_energy,
-        )
-        if atom_type.epsilon / ref_energy == 0:
-            lj.r_cut[(name, name)] = 0
-        else:
-            lj.r_cut[(name, name)] = r_cut
-
-    # Cross interactions, mixing rules, NBfixes
-    all_atomtypes = sorted(atom_type_params.keys())
-    for a1, a2 in itertools.combinations_with_replacement(all_atomtypes, 2):
-        nb_fix_info = atom_type_params[a1].nbfix.get(a2, None)
-        # nb_fix_info = (rmin, eps, rmin14, eps14)
-        if nb_fix_info is None:
-            # No nbfix means use mixing rule to find cross-interaction
-            if structure.combining_rule == "lorentz":
-                sigma = (
-                    atom_type_params[a1].sigma + atom_type_params[a2].sigma
-                ) / (2 * ref_distance)
-                epsilon = (
-                    (
-                        atom_type_params[a1].epsilon
-                        * atom_type_params[a2].epsilon
-                    )
-                    / ref_energy**2
-                ) ** 0.5
-            elif structure.combining_rule == "geometric":
-                sigma = (
-                    (atom_type_params[a1].sigma * atom_type_params[a2].sigma)
-                    / ref_distance**2
-                ) ** 0.5
-                epsilon = (
-                    (
-                        atom_type_params[a1].epsilon
-                        * atom_type_params[a2].epsilon
-                    )
-                    / ref_energy**2
-                ) ** 0.5
-            else:
-                raise ValueError(
-                    f"Mixing rule {structure.combining_rule} not supported, "
-                    'use "lorentz" or "geometric"'
-                )
-        else:
-            # If we have nbfix info, use it
-            sigma = nb_fix_info[0] / (ref_distance * (2 ** (1 / 6)))
-            epsilon = nb_fix_info[1] / ref_energy
-        lj.params[(a1, a2)] = dict(sigma=sigma, epsilon=epsilon)
-        if epsilon == 0:
-            lj.r_cut[(a1, a2)] = 0
-        else:
-            lj.r_cut[(a1, a2)] = r_cut
-
-    return lj
-
-
-def _init_hoomd_qq(structure, nl, snapshot, r_cut, Nx=1, Ny=1, Nz=1, order=4):
-    """Charge interactions."""
-    num_charged = np.sum(snapshot.particles.charge[:] != 0)
-    if num_charged == 0:
-        print("No charged groups found, ignoring electrostatics")
-        return None
-    else:
-        qq = hoomd.md.long_range.pppm.make_pppm_coulomb_forces(
-            nlist=nl, resolution=(Nx, Ny, Nz), order=order, r_cut=r_cut
-        )
-        return qq
-
-
-def _init_hoomd_14_pairs(
-    structure, nl, snapshot, r_cut, ref_distance=1.0, ref_energy=1.0
-):
-    """Special_pairs to handle 14 scaling.
-
-    See discussion: https://groups.google.com/forum/
-    #!topic/hoomd-users/iZ9WCpHczg0
-    """
-    # Update neighborlist to exclude 1-4 interactions,
-    # but impose a special_pair force to handle these pairs
-    nl.exclusions = nl.exclusions + [
-        "1-4",
-    ]
-
-    if snapshot.pairs.N == 0:
-        print("No 1,4 pairs found in hoomd snapshot")
-        return None, None
-
-    lj_14 = hoomd.md.special_pair.LJ()
-    qq_14 = hoomd.md.special_pair.Coulomb()
-    params_14 = {}
-    # Identify unique 14 scalings
-    for adjust in structure.adjusts:
-        t1 = adjust.atom1.type
-        t2 = adjust.atom2.type
-        ps = "-".join(sorted([t1, t2]))
-        if ps not in params_14:
-            params_14[ps] = adjust.type
-    for name, adjust_type in params_14.items():
-        lj_14.params[name] = dict(
-            sigma=adjust_type.sigma / ref_distance,
-            # The adjust epsilon already carries the scaling
-            epsilon=adjust_type.epsilon / ref_energy,
-        )
-        if adjust_type.epsilon / ref_energy == 0:
-            lj_14.r_cut[name] = 0
-        else:
-            lj_14.r_cut[name] = r_cut
-        qq_14.params[name] = dict(alpha=adjust_type.chgscale)
-        qq_14.r_cut[name] = r_cut
-
-    return lj_14, qq_14
-
-
-def _init_hoomd_bonds(structure, ref_distance=1.0, ref_energy=1.0):
-    """Harmonic bonds."""
-    # Identify the unique bond types before setting
-    bond_type_params = {}
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        t1, t2 = sorted([t1, t2], key=natural_sort)
-        if t1 != "" and t2 != "":
-            bond_type = "-".join((t1, t2))
-            if bond_type not in bond_type_params:
-                bond_type_params[bond_type] = bond.type
-
-    # Set the hoomd parameters
-    harmonic_bond = hoomd.md.bond.Harmonic()
-    for name, bond_type in bond_type_params.items():
-        # A (paramerized) parmed structure with no bondtype
-        # is because of constraints
-        if bond_type is None:
-            print("Bond with no bondtype detected, setting coefficients to 0")
-            harmonic_bond.params[name] = dict(k=0, r0=0)
-        else:
-            harmonic_bond.params[name] = dict(
-                k=2 * bond_type.k * ref_distance**2 / ref_energy,
-                r0=bond_type.req / ref_distance,
-            )
-
-    return harmonic_bond
-
-
-def _init_hoomd_angles(structure, ref_energy=1.0):
-    """Harmonic angles."""
-    # Identify the unique angle types before setting
-    angle_type_params = {}
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3], key=natural_sort)
-        angle_type = "-".join((t1, t2, t3))
-        if angle_type not in angle_type_params:
-            angle_type_params[angle_type] = angle.type
-
-    # set the hoomd parameters
-    harmonic_angle = hoomd.md.angle.Harmonic()
-    for name, angle_type in angle_type_params.items():
-        harmonic_angle.params[name] = dict(
-            t0=np.deg2rad(angle_type.theteq),
-            k=2 * angle_type.k / ref_energy,
-        )
-
-    return harmonic_angle
-
-
-def _init_hoomd_dihedrals(structure, ref_energy=1.0):
-    """Periodic dihedrals (dubbed harmonic dihedrals in HOOMD)."""
-    dihedral_type_params = {}
-    for dihedral in structure.dihedrals:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        if dihedral_type not in dihedral_type_params:
-            if isinstance(dihedral.type, pmd.DihedralType):
-                dihedral_type_params[dihedral_type] = dihedral.type
-            elif isinstance(dihedral.type, pmd.DihedralTypeList):
-                if len(dihedral.type) > 1:
-                    warnings.warn(
-                        "Multiple dihedral types detected"
-                        + " for single dihedral, will ignore all except "
-                        + " first dihedral type."
-                        + "First dihedral type: {}".format(dihedral.type[0])
-                    )
-                dihedral_type_params[dihedral_type] = dihedral.type[0]
-
-    # Set the hoomd parameters
-    # These are periodic torsions
-    hoomd_version = _get_hoomd_version()
-    if (
-        hoomd_version.major == 3 and hoomd_version.minor >= 7
-    ) or hoomd_version.major == 4:
-        periodic_torsion = hoomd.md.dihedral.Periodic()
-    else:
-        periodic_torsion = hoomd.md.dihedral.Harmonic()
-    for name, dihedral_type in dihedral_type_params.items():
-        periodic_torsion.params[name] = dict(
-            k=2 * dihedral_type.phi_k / ref_energy,
-            d=1,
-            n=dihedral_type.per,
-            phi0=np.deg2rad(dihedral_type.phase),
-        )
-
-    return periodic_torsion
-
-
-def _init_hoomd_rb_torsions(structure, ref_energy=1.0):
-    """RB dihedrals (implemented as OPLS dihedrals in HOOMD)."""
-    # Identify the unique dihedral types before setting
-    dihedral_type_params = {}
-    for dihedral in structure.rb_torsions:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        if dihedral_type not in dihedral_type_params:
-            dihedral_type_params[dihedral_type] = dihedral.type
-
-    # Set the hoomd parameter
-    rb_torsion = hoomd.md.dihedral.OPLS()
-    for name, dihedral_type in dihedral_type_params.items():
-        F_coeffs = RB_to_OPLS(
-            dihedral_type.c0 / ref_energy,
-            dihedral_type.c1 / ref_energy,
-            dihedral_type.c2 / ref_energy,
-            dihedral_type.c3 / ref_energy,
-            dihedral_type.c4 / ref_energy,
-            dihedral_type.c5 / ref_energy,
-        )
-        rb_torsion.params[name] = dict(
-            k1=F_coeffs[1], k2=F_coeffs[2], k3=F_coeffs[3], k4=F_coeffs[4]
-        )
-
-    return rb_torsion
diff --git a/mbuild/formats/hoomd_simulation.py b/mbuild/formats/hoomd_simulation.py
deleted file mode 100644
index 12da75818..000000000
--- a/mbuild/formats/hoomd_simulation.py
+++ /dev/null
@@ -1,467 +0,0 @@
-"""HOOMD simulation format."""
-
-import itertools
-import operator
-import warnings
-from collections import namedtuple
-
-import numpy as np
-import parmed as pmd
-
-import mbuild as mb
-from mbuild.utils.conversion import RB_to_OPLS
-from mbuild.utils.decorators import deprecated
-from mbuild.utils.io import import_
-from mbuild.utils.sorting import natural_sort
-
-from .hoomd_snapshot import to_hoomdsnapshot
-
-gsd = import_("gsd")
-gsd.hoomd = import_("gsd.hoomd")
-hoomd = import_("hoomd")
-hoomd.md = import_("hoomd.md")
-
-
-dep_msg = """
-Support for Hoomd-Blue 2.0 will be removed in mBuild 1.0.
-See GMSO (https://github.com/mosdef-hub/gmso/tree/main/gmso/external/convert_hoomd) for
-Hoomd-Blue 3.x and 4.x format support.
-To convert to GMSO, use the method `Compound.to_gmso()`.
-"""
-print(dep_msg)
-
-
-@deprecated(dep_msg)
-def create_hoomd_simulation(
-    structure,
-    r_cut,
-    ref_distance=1.0,
-    ref_mass=1.0,
-    ref_energy=1.0,
-    auto_scale=False,
-    snapshot_kwargs={},
-    pppm_kwargs={"Nx": 8, "Ny": 8, "Nz": 8, "order": 4},
-    init_snap=None,
-    restart=None,
-    nlist=hoomd.md.nlist.cell,
-):
-    """Convert a parametrized pmd.Structure to hoomd.SimulationContext.
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd Structure object
-    r_cut : float
-        Cutoff radius in simulation units
-    ref_distance : float, optional, default=1.0
-        Reference distance for unit conversion (from Angstrom)
-    ref_mass : float, optional, default=1.0
-        Reference mass for unit conversion (from Dalton)
-    ref_energy : float, optional, default=1.0
-        Reference energy for unit conversion (from kcal/mol)
-    auto_scale : bool, optional, default=False
-        Scale to reduced units by automatically using the largest sigma value
-        as ref_distance, largest mass value as ref_mass, and largest epsilon
-        value as ref_energy
-    snapshot_kwargs : dict
-        Kwargs to pass to to_hoomdsnapshot
-    pppm_kwargs : dict
-        Kwargs to pass to hoomd's pppm function
-    init_snap : hoomd.data.SnapshotParticleData, optional, default=None
-        Initial snapshot to which to add the ParmEd structure object
-        (useful for rigid bodies)
-    restart : str, optional, default=None
-        Path to the gsd file from which to restart the simulation.
-        https://hoomd-blue.readthedocs.io/en/v2.9.4/restartable-jobs.html
-        Note: It is assumed that the ParmEd structure and the system in
-        restart.gsd contain the same types. The ParmEd structure is still used
-        to initialize the forces, but restart.gsd is used to initialize the
-        system state (e.g., particle positions, momenta, etc).
-    nlist : hoomd.md.nlist, default=hoomd.md.nlist.cell
-        Type of neighborlist to use, see
-        https://hoomd-blue.readthedocs.io/en/stable/module-md-nlist.html
-        for more information.
-
-    Returns
-    -------
-    hoomd_objects : list
-        List of hoomd objects created during conversion
-    ReferenceValues : namedtuple
-        Values used in scaling
-
-    Note
-    ----
-    While the hoomd objects are returned, the hoomd.SimulationContext is
-    accessible via `hoomd.context.current`. If you pass a non-parametrized
-    pmd.Structure, you will not have angle, dihedral, or force field
-    information. You may be better off creating a hoomd.Snapshot.
-    Reference units should be expected to convert parmed Structure units:
-
-    --- angstroms, kcal/mol, and daltons
-    """
-    if isinstance(structure, mb.Compound):
-        raise ValueError(
-            "You passed mb.Compound to create_hoomd_simulation, there will be "
-            "no angles, dihedrals, or force field parameters. Please use "
-            "hoomd_snapshot.to_hoomdsnapshot to create a hoomd.Snapshot, then "
-            "create your own hoomd context and pass your hoomd.Snapshot to "
-            "hoomd.init.read_snapshot()"
-        )
-    elif not isinstance(structure, pmd.Structure):
-        raise ValueError(
-            "Please pass a parmed.Structure to create_hoomd_simulation"
-        )
-
-    _check_hoomd_version()
-    version_numbers = _check_hoomd_version()
-    if float(version_numbers[0]) >= 3:
-        raise RuntimeError("This method does not support HOOMD-blue v3.x.")
-
-    hoomd_objects = []  # Potential adaptation for Hoomd v3 API
-
-    if auto_scale:
-        if not all([i == 1 for i in (ref_distance, ref_energy, ref_mass)]):
-            warnings.warn(
-                "Autoscale option selected--provided reference values will not "
-                "be used."
-            )
-        pair_coeffs = list(
-            set((a.type, a.epsilon, a.sigma) for a in structure.atoms)
-        )
-
-        ref_mass = max([atom.mass for atom in structure.atoms])
-        ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1]
-        ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2]
-
-    ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"])
-    ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy)
-
-    if not hoomd.context.current:
-        hoomd.context.initialize("")
-    if restart is None:
-        snapshot, _ = to_hoomdsnapshot(
-            structure,
-            ref_distance=ref_distance,
-            ref_mass=ref_mass,
-            ref_energy=ref_energy,
-            **snapshot_kwargs,
-            hoomd_snapshot=init_snap,
-        )
-        hoomd_objects.append(snapshot)
-        hoomd_system = hoomd.init.read_snapshot(snapshot)
-        hoomd_objects.append(hoomd_system)
-    else:
-        with gsd.hoomd.open(restart) as f:
-            snapshot = f[-1]
-        hoomd_objects.append(snapshot)
-        hoomd_system = hoomd.init.read_gsd(restart, restart=restart)
-        hoomd_objects.append(hoomd_system)
-        print("Simulation initialized from restart file")
-
-    nl = nlist()
-    nl.reset_exclusions(exclusions=["1-2", "1-3"])
-    hoomd_objects.append(nl)
-
-    if structure.atoms[0].type != "":
-        print("Processing LJ and QQ")
-        lj = _init_hoomd_lj(
-            structure,
-            nl,
-            r_cut,
-            ref_distance=ref_distance,
-            ref_energy=ref_energy,
-        )
-        qq = _init_hoomd_qq(structure, nl, r_cut, **pppm_kwargs)
-        hoomd_objects.append(lj)
-        hoomd_objects.append(qq)
-    if structure.adjusts:
-        print("Processing 1-4 interactions, adjusting neighborlist exclusions")
-        lj_14, qq_14 = _init_hoomd_14_pairs(
-            structure,
-            nl,
-            r_cut,
-            ref_distance=ref_distance,
-            ref_energy=ref_energy,
-        )
-        hoomd_objects.append(lj_14)
-        hoomd_objects.append(qq_14)
-    if structure.bond_types:
-        print("Processing harmonic bonds")
-        harmonic_bond = _init_hoomd_bonds(
-            structure, ref_distance=ref_distance, ref_energy=ref_energy
-        )
-        hoomd_objects.append(harmonic_bond)
-    if structure.angle_types:
-        print("Processing harmonic angles")
-        harmonic_angle = _init_hoomd_angles(structure, ref_energy=ref_energy)
-        hoomd_objects.append(harmonic_angle)
-    if structure.dihedral_types:
-        print("Processing periodic torsions")
-        periodic_torsions = _init_hoomd_dihedrals(
-            structure, ref_energy=ref_energy
-        )
-        hoomd_objects.append(periodic_torsions)
-    if structure.rb_torsion_types:
-        print("Processing RB torsions")
-        rb_torsions = _init_hoomd_rb_torsions(structure, ref_energy=ref_energy)
-        hoomd_objects.append(rb_torsions)
-    print("HOOMD SimulationContext updated from ParmEd Structure")
-
-    return hoomd_objects, ref_values
-
-
-def _init_hoomd_lj(structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0):
-    """LJ parameters."""
-    # Identify the unique atom types before setting
-    atom_type_params = {}
-    for atom in structure.atoms:
-        if atom.type not in atom_type_params:
-            atom_type_params[atom.type] = atom.atom_type
-
-    # Set the hoomd parameters for self-interactions
-    lj = hoomd.md.pair.lj(r_cut, nl)
-    for name, atom_type in atom_type_params.items():
-        lj.pair_coeff.set(
-            name,
-            name,
-            sigma=atom_type.sigma / ref_distance,
-            epsilon=atom_type.epsilon / ref_energy,
-        )
-        if atom_type.epsilon / ref_energy == 0:
-            lj.pair_coeff.set(name, name, r_cut=0)
-
-    # Cross interactions, mixing rules, NBfixes
-    all_atomtypes = sorted(atom_type_params.keys())
-    for a1, a2 in itertools.combinations_with_replacement(all_atomtypes, 2):
-        nb_fix_info = atom_type_params[a1].nbfix.get(a2, None)
-        # nb_fix_info = (rmin, eps, rmin14, eps14)
-        if nb_fix_info is None:
-            # No nbfix means use mixing rule to find cross-interaction
-            if structure.combining_rule == "lorentz":
-                sigma = (
-                    atom_type_params[a1].sigma + atom_type_params[a2].sigma
-                ) / (2 * ref_distance)
-                epsilon = (
-                    (
-                        atom_type_params[a1].epsilon
-                        * atom_type_params[a2].epsilon
-                    )
-                    / ref_energy**2
-                ) ** 0.5
-            elif structure.combining_rule == "geometric":
-                sigma = (
-                    (atom_type_params[a1].sigma * atom_type_params[a2].sigma)
-                    / ref_distance**2
-                ) ** 0.5
-                epsilon = (
-                    (
-                        atom_type_params[a1].epsilon
-                        * atom_type_params[a2].epsilon
-                    )
-                    / ref_energy**2
-                ) ** 0.5
-            else:
-                raise ValueError(
-                    f"Mixing rule {structure.combining_rule} not supported, "
-                    "use lorentz"
-                )
-        else:
-            # If we have nbfix info, use it
-            sigma = nb_fix_info[0] / (ref_distance * (2 ** (1 / 6)))
-            epsilon = nb_fix_info[1] / ref_energy
-        lj.pair_coeff.set(a1, a2, sigma=sigma, epsilon=epsilon)
-        if epsilon == 0:
-            lj.pair_coeff.set(a1, a2, r_cut=0)
-
-    return lj
-
-
-def _init_hoomd_qq(structure, nl, r_cut, Nx=1, Ny=1, Nz=1, order=4):
-    """Charge interactions."""
-    charged = hoomd.group.charged()
-    if len(charged) == 0:
-        print("No charged groups found, ignoring electrostatics")
-        return None
-    else:
-        qq = hoomd.md.charge.pppm(charged, nl)
-        qq.set_params(Nx, Ny, Nz, order, r_cut)
-        return qq
-
-
-def _init_hoomd_14_pairs(
-    structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0
-):
-    """Special_pairs to handle 14 scalings.
-
-    See discussion: https://groups.google.com/forum/#!topic/hoomd-users/
-    iZ9WCpHczg0
-    """
-    # Update neighborlist to exclude 1-4 interactions,
-    # but impose a special_pair force to handle these pairs
-    nl.exclusions.append("1-4")
-
-    if hoomd.context.current.system_definition.getPairData().getN() == 0:
-        print("No 1,4 pairs found in hoomd snapshot")
-        return None, None
-
-    lj_14 = hoomd.md.special_pair.lj()
-    qq_14 = hoomd.md.special_pair.coulomb()
-    params_14 = {}
-    # Identify unique 14 scalings
-    for adjust in structure.adjusts:
-        t1 = adjust.atom1.type
-        t2 = adjust.atom2.type
-        ps = "-".join(sorted([t1, t2]))
-        if ps not in params_14:
-            params_14[ps] = adjust.type
-    for name, adjust_type in params_14.items():
-        lj_14.pair_coeff.set(
-            name,
-            sigma=adjust_type.sigma / ref_distance,
-            # The adjust epsilon already carries the scaling
-            epsilon=adjust_type.epsilon / ref_energy,
-            # Do NOT use hoomd's alpha to modify any LJ terms
-            alpha=1,
-            r_cut=r_cut,
-        )
-        qq_14.pair_coeff.set(name, alpha=adjust_type.chgscale, r_cut=r_cut)
-
-    return lj_14, qq_14
-
-
-def _init_hoomd_bonds(structure, ref_distance=1.0, ref_energy=1.0):
-    """Harmonic bonds."""
-    # Identify the unique bond types before setting
-    bond_type_params = {}
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        t1, t2 = sorted([t1, t2], key=natural_sort)
-        if t1 != "" and t2 != "":
-            bond_type = "-".join((t1, t2))
-            if bond_type not in bond_type_params:
-                bond_type_params[bond_type] = bond.type
-
-    # Set the hoomd parameters
-    harmonic_bond = hoomd.md.bond.harmonic()
-    for name, bond_type in bond_type_params.items():
-        # A (paramerized) parmed structure with no bondtype
-        # is because of constraints
-        if bond_type is None:
-            print("Bond with no bondtype detected, setting coefficients to 0")
-            harmonic_bond.bond_coeff.set(name, k=0, r0=0)
-        else:
-            harmonic_bond.bond_coeff.set(
-                name,
-                k=2 * bond_type.k * ref_distance**2 / ref_energy,
-                r0=bond_type.req / ref_distance,
-            )
-
-    return harmonic_bond
-
-
-def _init_hoomd_angles(structure, ref_energy=1.0):
-    """Harmonic angles."""
-    # Identify the unique angle types before setting
-    angle_type_params = {}
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3], key=natural_sort)
-        angle_type = "-".join((t1, t2, t3))
-        if angle_type not in angle_type_params:
-            angle_type_params[angle_type] = angle.type
-
-    # set the hoomd parameters
-    harmonic_angle = hoomd.md.angle.harmonic()
-    for name, angle_type in angle_type_params.items():
-        harmonic_angle.angle_coeff.set(
-            name,
-            t0=np.deg2rad(angle_type.theteq),
-            k=2 * angle_type.k / ref_energy,
-        )
-
-    return harmonic_angle
-
-
-def _init_hoomd_dihedrals(structure, ref_energy=1.0):
-    """Periodic dihedrals (dubbed harmonic dihedrals in HOOMD)."""
-    # Identify the unique dihedral types before setting
-    # need Hoomd 2.8.0 to use proper dihedral implemtnation
-    # from this PR https://github.com/glotzerlab/hoomd-blue/pull/492
-    version_numbers = _check_hoomd_version()
-    if float(version_numbers[0]) < 2 or float(version_numbers[1]) < 8:
-        from mbuild.exceptions import MBuildError
-
-        raise MBuildError("Please upgrade Hoomd to at least 2.8.0")
-
-    dihedral_type_params = {}
-    for dihedral in structure.dihedrals:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        if dihedral_type not in dihedral_type_params:
-            if isinstance(dihedral.type, pmd.DihedralType):
-                dihedral_type_params[dihedral_type] = dihedral.type
-            elif isinstance(dihedral.type, pmd.DihedralTypeList):
-                if len(dihedral.type) > 1:
-                    warnings.warn(
-                        "Multiple dihedral types detected for single dihedral, "
-                        "will ignore all except  first dihedral type. First "
-                        "dihedral type: {}".format(dihedral.type[0])
-                    )
-                dihedral_type_params[dihedral_type] = dihedral.type[0]
-
-    # Set the hoomd parameters
-    # These are periodic torsions
-    periodic_torsion = hoomd.md.dihedral.harmonic()
-    for name, dihedral_type in dihedral_type_params.items():
-        periodic_torsion.dihedral_coeff.set(
-            name,
-            k=2 * dihedral_type.phi_k / ref_energy,
-            d=1,
-            n=dihedral_type.per,
-            phi_0=np.deg2rad(dihedral_type.phase),
-        )
-
-    return periodic_torsion
-
-
-def _init_hoomd_rb_torsions(structure, ref_energy=1.0):
-    """RB dihedrals (implemented as OPLS dihedrals in HOOMD)."""
-    # Identify the unique dihedral types before setting
-    dihedral_type_params = {}
-    for dihedral in structure.rb_torsions:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        if dihedral_type not in dihedral_type_params:
-            dihedral_type_params[dihedral_type] = dihedral.type
-
-    # Set the hoomd parameter
-    rb_torsion = hoomd.md.dihedral.opls()
-    for name, dihedral_type in dihedral_type_params.items():
-        F_coeffs = RB_to_OPLS(
-            dihedral_type.c0 / ref_energy,
-            dihedral_type.c1 / ref_energy,
-            dihedral_type.c2 / ref_energy,
-            dihedral_type.c3 / ref_energy,
-            dihedral_type.c4 / ref_energy,
-            dihedral_type.c5 / ref_energy,
-            error_if_outside_tolerance=False,
-        )
-        rb_torsion.dihedral_coeff.set(
-            name, k1=F_coeffs[1], k2=F_coeffs[2], k3=F_coeffs[3], k4=F_coeffs[4]
-        )
-
-    return rb_torsion
-
-
-def _check_hoomd_version():
-    version = hoomd.__version__
-    version_numbers = version.split(".")
-    return version_numbers
diff --git a/mbuild/formats/hoomd_snapshot.py b/mbuild/formats/hoomd_snapshot.py
deleted file mode 100644
index b52b9db72..000000000
--- a/mbuild/formats/hoomd_snapshot.py
+++ /dev/null
@@ -1,622 +0,0 @@
-"""HOOMD snapshot format."""
-
-import operator
-from collections import namedtuple
-
-import numpy as np
-import packaging.version
-import parmed as pmd
-
-from mbuild.box import Box
-from mbuild.compound import Compound, Particle
-from mbuild.utils.geometry import coord_shift
-from mbuild.utils.io import import_
-from mbuild.utils.sorting import natural_sort
-
-hoomd = import_("hoomd")
-
-__all__ = ["to_hoomdsnapshot", "from_snapshot"]
-
-
-def _get_hoomd_version():
-    if "version" in dir(hoomd):
-        return packaging.version.parse(hoomd.version.version)
-    else:
-        return packaging.version.parse(hoomd.__version__)
-
-
-def from_snapshot(snapshot, scale=1.0):
-    """Convert a Snapshot to a Compound.
-
-    Snapshot can be a hoomd.Snapshot or a gsd.hoomd.Frame.
-
-    Parameters
-    ----------
-    snapshot : hoomd.Snapshot or gsd.hoomd.Frame
-        Snapshot from which to build the mbuild Compound.
-    scale : float, optional, default 1.0
-        Value by which to scale the length values
-
-    Returns
-    -------
-    comp : Compound
-
-    Note
-    ----
-    GSD and HOOMD snapshots center their boxes on the origin (0,0,0), so the
-    compound is shifted by half the box lengths
-    """
-    comp = Compound()
-    bond_array = snapshot.bonds.group
-    n_atoms = snapshot.particles.N
-
-    if "SnapshotSystemData_float" in dir(hoomd._hoomd) and isinstance(
-        snapshot, hoomd._hoomd.SnapshotSystemData_float
-    ):
-        # hoomd v2
-        box = snapshot.box
-        comp.box = Box.from_lengths_tilt_factors(
-            lengths=np.array([box.Lx, box.Ly, box.Lz]) * scale,
-            tilt_factors=np.array([box.xy, box.xz, box.yz]),
-        )
-    else:
-        # gsd / hoomd v3
-        box = np.asarray(snapshot.configuration.box)
-        comp.box = Box.from_lengths_tilt_factors(
-            lengths=box[:3] * scale, tilt_factors=box[3:]
-        )
-
-    # GSD and HOOMD snapshots center their boxes on the origin (0,0,0)
-    shift = np.array(comp.box.lengths) / 2
-    # Add particles
-    for i in range(n_atoms):
-        name = snapshot.particles.types[snapshot.particles.typeid[i]]
-        xyz = snapshot.particles.position[i] * scale + shift
-        charge = snapshot.particles.charge[i]
-
-        atom = Particle(name=name, pos=xyz, charge=charge)
-        comp.add(atom, label=str(i))
-
-    # Add bonds
-    particle_dict = {idx: p for idx, p in enumerate(comp.particles())}
-    for i in range(bond_array.shape[0]):
-        atom1 = int(bond_array[i][0])
-        atom2 = int(bond_array[i][1])
-        comp.add_bond([particle_dict[atom1], particle_dict[atom2]])
-    return comp
-
-
-def to_hoomdsnapshot(
-    structure,
-    ref_distance=1.0,
-    ref_mass=1.0,
-    ref_energy=1.0,
-    rigid_bodies=None,
-    shift_coords=True,
-    write_special_pairs=True,
-    auto_scale=False,
-    parmed_kwargs={},
-    hoomd_snapshot=None,
-):
-    """Convert a Compound or parmed.Structure to hoomd.Snapshot.
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd Structure object
-        Reference distance for unit conversion ((Angstrom) / (desired units))
-    ref_mass : float, optional, default=1.0
-        Reference mass for unit conversion ((Dalton) / (desired units))
-    ref_energy : float, optional, default=1.0
-        Reference energy for unit conversion ((kcal/mol) / (desired units))
-    rigid_bodies : list of int, optional, default=None
-        List of rigid body information. An integer value is required for
-        each atom corresponding to the index of the rigid body the particle
-        is to be associated with. A value of None indicates the atom is not
-        part of a rigid body.
-    shift_coords : bool, optional, default=True
-        Shift coordinates from (0, L) to (-L/2, L/2) if necessary.
-    auto_scale : bool, optional, default=False
-        Automatically use largest sigma value as ref_distance,
-        largest mass value as ref_mass
-        and largest epsilon value as ref_energy
-    write_special_pairs : bool, optional, default=True
-        Writes out special pair information necessary to correctly use
-        the OPLS fudged 1,4 interactions in HOOMD.
-    hoomd_snapshot : hoomd.Snapshot, optional, default=None
-        Initial snapshot to which to add the ParmEd structure object.
-        The box information of the initial snapshot will be overwritten.
-        (useful for rigid bodies)
-
-    Returns
-    -------
-    hoomd_snapshot : hoomd.Snapshot
-    ReferenceValues : namedtuple
-        Values used in scaling
-
-    Notes
-    -----
-    This method does not create hoomd forcefield objects and
-    the snapshot returned does not store the forcefield parameters.
-    See mbuild.formats.hoomd_forcefield.create_hoomd_forcefield()
-
-    About units: This method operates on a Parmed.Structure object
-        where the units used differ from those used in mBuild and Foyer
-        which may have been used when creating the typed Parmed.Structure.
-
-    The default units used when writing out the HOOMD Snapshot are:
-    Distance (Angstrom)
-    Mass (Dalton)
-    Energy (kcal/mol)
-
-    If you wish to convert this unit system to another, you can use the
-    reference parameters (ref_distance, ref_mass, ref_energy).
-    The values used here should be expected to convert from the Parmed
-    Structure units (above) to your desired units.
-    The Parmed.Structure values are divided by the reference values.
-
-    If you wish to used a reduced unit system, set auto_scale = True.
-    When auto_scale is True, the reference parameters won't be used.
-
-    Examples
-    --------
-    To convert the energy units from kcal/mol to kj/mol:
-        use ref_energy = 0.2390057 (kcal/kj)
-
-    To convert the distance units from Angstrom to nm:
-        use ref_distance = 10 (angstroms/nm)
-
-    To use a reduced unit system, where mass, sigma, and epsilon are
-    scaled by the largest value of each:
-        use auto_scale = True, ref_distance = ref_energy = ref_mass = 1
-
-    """
-    if not isinstance(structure, (Compound, pmd.Structure)):
-        raise ValueError(
-            "You are trying to create a hoomd.Snapshot from a "
-            f"{type(structure)} please pass a Compound or pmd.Structure"
-        )
-    elif isinstance(structure, Compound):
-        structure = structure.to_parmed(**parmed_kwargs)
-
-    hoomd_version = _get_hoomd_version()
-    if hoomd_version.major == 2 and not hoomd.context.current:
-        hoomd.context.initialize("")
-
-    if auto_scale:
-        ref_mass = max([atom.mass for atom in structure.atoms])
-        pair_coeffs = list(
-            set(
-                (atom.type, atom.epsilon, atom.sigma)
-                for atom in structure.atoms
-            )
-        )
-        ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1]
-        ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2]
-
-    ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"])
-    ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy)
-
-    xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms])
-    if shift_coords:
-        xyz = coord_shift(xyz, structure.box[:3])
-
-    # Get box information
-    if np.allclose(structure.box[3:6], np.array([90, 90, 90])):
-        lx, ly, lz = structure.box[:3] / ref_distance
-        xy, xz, yz = 0, 0, 0
-    else:
-        a, b, c = structure.box[0:3] / ref_distance
-        alpha, beta, gamma = np.radians(structure.box[3:6])
-
-        lx = a
-        xy = b * np.cos(gamma)
-        xz = c * np.cos(beta)
-        ly = np.sqrt(b**2 - xy**2)
-        yz = (b * c * np.cos(alpha) - xy * xz) / ly
-        lz = np.sqrt(c**2 - xz**2 - yz**2)
-
-    (
-        n_particles,
-        scaled_positions,
-        unique_types,
-        typeids,
-        scaled_mass,
-        scaled_charges,
-        rigid_bodies,
-    ) = _parse_particle_information(
-        structure, xyz, ref_distance, ref_mass, ref_energy, rigid_bodies
-    )
-    (
-        n_bonds,
-        unique_bond_types,
-        bond_typeids,
-        bond_groups,
-    ) = _parse_bond_information(structure)
-    (
-        n_angles,
-        unique_angle_types,
-        angle_typeids,
-        angle_groups,
-    ) = _parse_angle_information(structure)
-    (
-        n_dihedrals,
-        unique_dihedral_types,
-        dihedral_typeids,
-        dihedral_groups,
-    ) = _parse_dihedral_information(structure)
-    (
-        n_impropers,
-        unique_improper_types,
-        improper_typeids,
-        improper_groups,
-    ) = _parse_improper_information(structure)
-    pair_types, pair_typeid, pairs, n_pairs = _parse_pair_information(structure)
-
-    if hoomd_snapshot is not None:
-        n_init = hoomd_snapshot.particles.N
-
-        if n_init > 0:
-            n_particles += n_init
-            # shift the typeids
-            typeids += len(set(hoomd_snapshot.particles.types))
-            hoomd_snapshot.particles.types += unique_types
-            # shift bond/angle/dihedral indices
-            if n_bonds > 0:
-                bond_groups = [
-                    tuple(row) for row in np.array(bond_groups) + n_init
-                ]
-                init_bondtypes = len(hoomd_snapshot.bonds.types)
-                hoomd_snapshot.bonds.types += unique_bond_types
-            if n_angles > 0:
-                angle_groups = [
-                    tuple(row) for row in np.array(angle_groups) + n_init
-                ]
-                init_angletypes = len(hoomd_snapshot.angles.types)
-                hoomd_snapshot.angles.types += unique_angle_types
-            if n_dihedrals > 0:
-                dihedral_groups = [
-                    tuple(row) for row in np.array(dihedral_groups) + n_init
-                ]
-                init_dihedraltypes = len(hoomd_snapshot.dihedrals.types)
-                hoomd_snapshot.dihedrals.types += unique_dihedral_types
-            if n_impropers > 0:
-                improper_groups = [
-                    tuple(row) for row in np.array(improper_groups) + n_init
-                ]
-                init_impropertypes = len(hoomd_snapshot.impropers.types)
-                hoomd_snapshot.impropers.types += unique_improper_types
-            if n_pairs > 0:
-                pairs = [tuple(row) for row in np.array(pairs) + n_init]
-                init_pairtypes = len(hoomd_snapshot.pairs.types)
-                hoomd_snapshot.pairs.types += pair_types
-        else:
-            raise RuntimeError(
-                "Initial snapshot provided, but it contains no particles"
-            )
-
-        if hoomd_version.major == 2:
-            hoomd_snapshot.box = hoomd.data.boxdim(
-                Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz
-            )
-
-            # save the box for later use when wrapping coordinates
-            box = hoomd_snapshot.box
-        elif hoomd_version.major == 3 or hoomd_version.major == 4:
-            hoomd_snapshot.configuration.box = [lx, ly, lz, xy, xz, yz]
-        else:
-            raise RuntimeError("Unsupported HOOMD version:", str(hoomd_version))
-
-        init_bonds = hoomd_snapshot.bonds.N
-        if init_bonds > 0:
-            n_bonds += init_bonds
-            bond_typeids = list(np.array(bond_typeids) + init_bondtypes)
-
-        init_angles = hoomd_snapshot.angles.N
-        if init_angles > 0:
-            n_angles += init_angles
-            angle_typeids = list(np.array(angle_typeids) + init_angletypes)
-
-        init_dihedrals = hoomd_snapshot.dihedrals.N
-        if init_dihedrals > 0:
-            n_dihedrals += init_dihedrals
-            dihedral_typeids = list(
-                np.array(dihedral_typeids) + init_dihedraltypes
-            )
-
-        init_impropers = hoomd_snapshot.impropers.N
-        if init_impropers > 0:
-            n_impropers += init_impropers
-            improper_typeids = list(
-                np.array(improper_typeids) + init_impropertypes
-            )
-
-        init_pairs = hoomd_snapshot.pairs.N
-        if init_pairs > 0:
-            n_pairs += init_pairs
-            pair_typeid = list(np.array(pair_typeid) + init_pairtypes)
-
-    else:
-        n_init = 0
-        init_bonds = 0
-        init_angles = 0
-        init_dihedrals = 0
-        init_impropers = 0
-        init_pairs = 0
-
-        if hoomd_version.major == 2:
-            hoomd_snapshot = hoomd.data.make_snapshot(
-                N=n_particles,
-                box=hoomd.data.boxdim(Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz),
-                particle_types=unique_types,
-                bond_types=unique_bond_types,
-                angle_types=unique_angle_types,
-                dihedral_types=unique_dihedral_types,
-                improper_types=unique_improper_types,
-                pair_types=pair_types,
-            )
-            box = hoomd.data.boxdim(Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz)
-        elif hoomd_version.major == 3 or hoomd_version.major == 4:
-            hoomd_snapshot = hoomd.Snapshot()
-            hoomd_snapshot.configuration.box = [lx, ly, lz, xy, xz, yz]
-            hoomd_snapshot.particles.types = unique_types
-            hoomd_snapshot.bonds.types = unique_bond_types
-            hoomd_snapshot.angles.types = unique_angle_types
-            hoomd_snapshot.dihedrals.types = unique_dihedral_types
-            hoomd_snapshot.impropers.types = unique_improper_types
-            hoomd_snapshot.pairs.types = pair_types
-            box = hoomd.Box(Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz)
-        else:
-            raise RuntimeError("Unsupported HOOMD version:", str(hoomd_version))
-
-    # wrap particles into the box manually for v2
-    if hoomd_version.major == 2:
-        scaled_positions = np.stack(
-            [box.wrap(xyz)[0] for xyz in scaled_positions]
-        )
-
-    def set_size(obj, n):
-        if hoomd_version.major == 2:
-            obj.resize(n)
-        elif hoomd_version.major == 3 or hoomd_version.major == 4:
-            obj.N = n
-        else:
-            raise RuntimeError("Unsupported HOOMD version:", str(hoomd_version))
-
-    set_size(hoomd_snapshot.particles, n_particles)
-    hoomd_snapshot.particles.position[n_init:] = scaled_positions
-    hoomd_snapshot.particles.typeid[n_init:] = typeids
-    hoomd_snapshot.particles.mass[n_init:] = scaled_mass
-    hoomd_snapshot.particles.charge[n_init:] = scaled_charges
-    hoomd_snapshot.particles.body[n_init:] = rigid_bodies
-
-    # wrap the particles into the box using the v3 API
-    if hoomd_version.major == 3 or hoomd_version.major == 4:
-        hoomd_snapshot.wrap()
-
-    if n_bonds > 0:
-        set_size(hoomd_snapshot.bonds, n_bonds)
-        hoomd_snapshot.bonds.typeid[init_bonds:] = bond_typeids
-        hoomd_snapshot.bonds.group[init_bonds:] = bond_groups
-
-    if n_angles > 0:
-        set_size(hoomd_snapshot.angles, n_angles)
-        hoomd_snapshot.angles.typeid[init_angles:] = angle_typeids
-        hoomd_snapshot.angles.group[init_angles:] = np.reshape(
-            angle_groups, (-1, 3)
-        )
-
-    if n_dihedrals > 0:
-        set_size(hoomd_snapshot.dihedrals, n_dihedrals)
-        hoomd_snapshot.dihedrals.typeid[init_dihedrals:] = dihedral_typeids
-        hoomd_snapshot.dihedrals.group[init_dihedrals:] = np.reshape(
-            dihedral_groups, (-1, 4)
-        )
-
-    if n_impropers > 0:
-        set_size(hoomd_snapshot.impropers, n_impropers)
-        hoomd_snapshot.impropers.typeid[init_impropers:] = improper_typeids
-        hoomd_snapshot.impropers.group[init_impropers:] = np.reshape(
-            improper_groups, (-1, 4)
-        )
-
-    if n_pairs > 0:
-        set_size(hoomd_snapshot.pairs, n_pairs)
-        hoomd_snapshot.pairs.typeid[init_pairs:] = pair_typeid
-        hoomd_snapshot.pairs.group[init_pairs:] = np.reshape(pairs, (-1, 2))
-
-    return hoomd_snapshot, ref_values
-
-
-def _parse_particle_information(
-    structure, xyz, ref_distance, ref_mass, ref_energy, rigid_bodies
-):
-    """Write out the particle information."""
-    n_particles = len(structure.atoms)
-    scaled_positions = xyz / ref_distance
-
-    types = [
-        atom.name if atom.type == "" else atom.type for atom in structure.atoms
-    ]
-
-    unique_types = list(set(types))
-    unique_types.sort(key=natural_sort)
-
-    typeids = np.array([unique_types.index(t) for t in types])
-
-    masses = np.array([atom.mass for atom in structure.atoms])
-    masses[masses == 0] = 1.0
-    scaled_mass = masses / ref_mass
-
-    charges = np.array([atom.charge for atom in structure.atoms])
-    e0 = 2.396452e-04
-    """
-    Permittivity of free space = 2.396452e-04 e^2/((kcal/mol) Angstrom),
-    where e is the elementary charge
-    """
-    charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5
-    scaled_charges = charges / charge_factor
-
-    if rigid_bodies:
-        rigid_bodies = [-1 if body is None else body for body in rigid_bodies]
-    else:
-        rigid_bodies = [-1 for _ in structure.atoms]
-
-    return (
-        n_particles,
-        scaled_positions,
-        unique_types,
-        typeids,
-        scaled_mass,
-        scaled_charges,
-        rigid_bodies,
-    )
-
-
-def _parse_pair_information(structure):
-    pair_types = []
-    pair_typeid = []
-    pairs = []
-    for ai in structure.atoms:
-        for aj in ai.dihedral_partners:
-            # make sure we don't double add
-            if ai.idx > aj.idx:
-                ps = "-".join(sorted([ai.type, aj.type]))
-                if ps not in pair_types:
-                    pair_types.append(ps)
-                pair_typeid.append(pair_types.index(ps))
-                pairs.append((ai.idx, aj.idx))
-    n_pairs = len(pairs)
-
-    return pair_types, pair_typeid, pairs, n_pairs
-
-
-def _parse_bond_information(structure):
-    n_bonds = len(structure.bonds)
-    unique_bond_types = set()
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        if t1 == "" or t2 == "":
-            t1, t2 = bond.atom1.name, bond.atom2.name
-        t1, t2 = sorted([t1, t2], key=natural_sort)
-        try:
-            bond_type = "-".join((t1, t2))
-        except AttributeError:  # no forcefield applied, bond.type is None
-            bond_type = ("-".join((t1, t2)), 0.0, 0.0)
-        unique_bond_types.add(bond_type)
-    unique_bond_types = sorted(list(unique_bond_types), key=natural_sort)
-
-    bond_typeids = []
-    bond_groups = []
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        if t1 == "" or t2 == "":
-            t1, t2 = bond.atom1.name, bond.atom2.name
-        t1, t2 = sorted([t1, t2], key=natural_sort)
-        try:
-            bond_type = "-".join((t1, t2))
-        except AttributeError:  # no forcefield applied, bond.type is None
-            bond_type = ("-".join((t1, t2)), 0.0, 0.0)
-        bond_typeids.append(unique_bond_types.index(bond_type))
-        bond_groups.append((bond.atom1.idx, bond.atom2.idx))
-
-    return n_bonds, unique_bond_types, bond_typeids, bond_groups
-
-
-def _parse_angle_information(structure):
-    n_angles = len(structure.angles)
-
-    unique_angle_types = set()
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3], key=natural_sort)
-        angle_type = "-".join((t1, t2, t3))
-        unique_angle_types.add(angle_type)
-    unique_angle_types = sorted(list(unique_angle_types), key=natural_sort)
-
-    angle_typeids = []
-    angle_groups = []
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3], key=natural_sort)
-        angle_type = "-".join((t1, t2, t3))
-        angle_typeids.append(unique_angle_types.index(angle_type))
-        angle_groups.append((angle.atom1.idx, angle.atom2.idx, angle.atom3.idx))
-
-    return n_angles, unique_angle_types, angle_typeids, angle_groups
-
-
-def _parse_dihedral_information(structure):
-    n_dihedrals = len(structure.rb_torsions + structure.dihedrals)
-
-    unique_dihedral_types = set()
-    for dihedral in structure.rb_torsions + structure.dihedrals:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        unique_dihedral_types.add(dihedral_type)
-    unique_dihedral_types = sorted(
-        list(unique_dihedral_types), key=natural_sort
-    )
-
-    dihedral_typeids = []
-    dihedral_groups = []
-    for dihedral in structure.rb_torsions + structure.dihedrals:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        dihedral_typeids.append(unique_dihedral_types.index(dihedral_type))
-        dihedral_groups.append(
-            (
-                dihedral.atom1.idx,
-                dihedral.atom2.idx,
-                dihedral.atom3.idx,
-                dihedral.atom4.idx,
-            )
-        )
-
-    return n_dihedrals, unique_dihedral_types, dihedral_typeids, dihedral_groups
-
-
-def _parse_improper_information(structure):
-    n_dihedrals = len(structure.impropers)
-
-    unique_dihedral_types = set()
-    for dihedral in structure.impropers:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        unique_dihedral_types.add(dihedral_type)
-    unique_dihedral_types = sorted(
-        list(unique_dihedral_types), key=natural_sort
-    )
-
-    dihedral_typeids = []
-    dihedral_groups = []
-    for dihedral in structure.impropers:
-        t1, t2 = dihedral.atom1.type, dihedral.atom2.type
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3], key=natural_sort):
-            dihedral_type = "-".join((t1, t2, t3, t4))
-        else:
-            dihedral_type = "-".join((t4, t3, t2, t1))
-        dihedral_typeids.append(unique_dihedral_types.index(dihedral_type))
-        dihedral_groups.append(
-            (
-                dihedral.atom1.idx,
-                dihedral.atom2.idx,
-                dihedral.atom3.idx,
-                dihedral.atom4.idx,
-            )
-        )
-
-    return n_dihedrals, unique_dihedral_types, dihedral_typeids, dihedral_groups
diff --git a/mbuild/formats/hoomdxml.py b/mbuild/formats/hoomdxml.py
deleted file mode 100644
index 1d026b653..000000000
--- a/mbuild/formats/hoomdxml.py
+++ /dev/null
@@ -1,433 +0,0 @@
-"""HOOMD xml format."""
-
-import operator
-from collections import namedtuple
-from math import radians
-
-import numpy as np
-
-from mbuild.utils.conversion import RB_to_OPLS
-from mbuild.utils.decorators import deprecated
-from mbuild.utils.geometry import coord_shift
-
-__all__ = ["write_hoomdxml"]
-
-
-@deprecated("The .hoomdxml file format will be removed in mBuild 1.0.")
-def write_hoomdxml(
-    structure,
-    filename,
-    ref_distance=1.0,
-    ref_mass=1.0,
-    ref_energy=1.0,
-    rigid_bodies=None,
-    shift_coords=True,
-    auto_scale=False,
-):
-    """Output a HOOMD XML file.
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd structure object
-    filename : str
-        Path of the output file.
-    ref_distance : float, optional, default=1.0, units=nanometers
-        Reference distance for conversion to reduced units
-    ref_mass : float, optional, default=1.0, units=amu
-        Reference mass for conversion to reduced units
-    ref_energy : float, optional, default=1.0, units=kJ/mol
-        Reference energy for conversion to reduced units
-    rigid_bodies : list
-        List of rigid body information. An integer value is required
-        for each particle corresponding to the number of the rigid body with
-        which the particle should be included. A value of None indicates the
-        particle is not part of any rigid body.
-    shift_coords : bool, optional, default=True
-        Shift coordinates from (0, L) to (-L/2, L/2) if necessary.
-    auto_scale : bool, optional, default=False
-        Automatically use largest sigma value as ref_distance, largest mass value
-        as ref_mass and largest epsilon value as ref_energy.
-
-    Returns
-    -------
-    ReferenceValues : namedtuple
-        Values used in scaling
-
-    Example
-    -------
-
-    >>> ref_values = ethane.save(
-    ...     filename='ethane-opls.hoomdxml',
-    ...     forcefield_name='oplsaa',
-    ...     auto_scale=True
-    ...     )
-    >>> print(ref_values.mass, ref_values.distance, ref_values.energy)
-
-    Notes
-    -----
-    The following elements are always written:
-
-    * **position** : particle positions
-    * **type** : particle types
-    * **mass** : particle masses (default 1.0)
-    * **charge** : particle charges
-
-    The following elements may be written if applicable:
-
-    * **pair_coeffs** : Pair coefficients for each particle type (assumes a
-        12-6 LJ pair style). The following information is written for each
-        particle type:
-
-        * type : particle type
-        * epsilon : LJ epsilon
-        * sigma : LJ sigma
-
-    * **bond_coeffs** : Coefficients for each bond type (assumes a harmonic bond
-        style). The following information is written for each bond type:
-
-        * type : bond type
-        * k : force constant (units of energy/distance^2)
-        * r0 : bond rest length (units of distance)
-
-    * **bond** : system bonds
-    * **angle_coeffs** : Coefficients for each angle type (assumes a harmonic
-        angle style). The following information is written for each angle type:
-
-        * type : angle type
-        * k : force constant (units of energy/radians^2)
-        * theta : rest angle (units of radians)
-
-    * **angle** : system angles
-    * **dihedral_coeffs** : Coefficients for each dihedral type (assumes an
-        OPLS dihedral style). The following information is written for each
-        dihedral type:
-
-        * type : dihedral type
-        * k1, k2, k3, k4 : force coefficients (units of energy)
-
-    * **dihedral** : system dihedrals
-
-    * **body** : ID of the rigid body to which each particle belongs
-    """
-    ref_distance *= 10  # Parmed unit hack
-    ref_energy /= 4.184  # Parmed unit hack
-    forcefield = True
-    if structure[0].type == "":
-        forcefield = False
-    if auto_scale and forcefield:
-        ref_mass = max([atom.mass for atom in structure.atoms])
-        pair_coeffs = list(
-            set(
-                (atom.type, atom.epsilon, atom.sigma)
-                for atom in structure.atoms
-            )
-        )
-        ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1]
-        ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2]
-
-    xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms])
-    if shift_coords:
-        xyz = coord_shift(xyz, structure.box[:3])
-
-    with open(filename, "w") as xml_file:
-        xml_file.write('<?xml version="1.2" encoding="UTF-8"?>\n')
-        xml_file.write('<hoomd_xml version="1.2">\n')
-        xml_file.write(
-            "<!-- ref_distance (nm) ref_mass (amu) ref_energy (kJ/mol) -->\n"
-        )
-        # We need to convert back into nm and kJ
-        xml_file.write(
-            "<!-- {} {} {} -->\n".format(
-                ref_distance / 10, ref_mass, ref_energy * 4.184
-            )
-        )
-        xml_file.write('<configuration time_step="0">\n')
-        _write_box_information(xml_file, structure, ref_distance)
-        _write_particle_information(
-            xml_file,
-            structure,
-            xyz,
-            forcefield,
-            ref_distance,
-            ref_mass,
-            ref_energy,
-        )
-        _write_bond_information(xml_file, structure, ref_distance, ref_energy)
-        _write_angle_information(xml_file, structure, ref_energy)
-        _write_dihedral_information(xml_file, structure, ref_energy)
-        _write_rigid_information(xml_file, rigid_bodies)
-        xml_file.write("</configuration>\n")
-        xml_file.write("</hoomd_xml>")
-
-    ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"])
-
-    return ReferenceValues(ref_distance, ref_mass, ref_energy)
-
-
-def _write_particle_information(
-    xml_file, structure, xyz, forcefield, ref_distance, ref_mass, ref_energy
-):
-    """Write out the particle information.
-
-    Parameters
-    ----------
-    xml_file : file object
-        The file object of the hoomdxml file being written
-    structure : parmed.Structure
-        Parmed structure object
-    xyz : np.ndarray, shape=(n,3), dtype=float
-        The particle positions to be written.
-    forcefield : bool
-        If True, write the particle "type". Write the particle "name" otherwise.
-    ref_distance : float, default=1.0
-        Reference distance for conversion to reduced units
-    ref_mass : float, default=1.0
-        Reference mass for conversion to reduced units
-    ref_energy : float, default=1.0
-        Reference energy for conversion to reduced units
-    """
-    xml_file.write('<position units="sigma" num="{}">\n'.format(xyz.shape[0]))
-    for pos in xyz:
-        xml_file.write("{}\t{}\t{}\n".format(*pos / ref_distance))
-    xml_file.write("</position>\n")
-    if forcefield:
-        types = [atom.type for atom in structure.atoms]
-    else:
-        types = [atom.name for atom in structure.atoms]
-
-    xml_file.write("<type>\n")
-    for atom_type in types:
-        xml_file.write("{}\n".format(atom_type))
-    xml_file.write("</type>\n")
-
-    masses = [atom.mass for atom in structure.atoms]
-    xml_file.write("<mass>\n")
-    for mass in masses:
-        if mass == 0:
-            mass = 1.0
-        xml_file.write("{}\n".format(mass / ref_mass))
-    xml_file.write("</mass>\n")
-
-    charges = [atom.charge for atom in structure.atoms]
-    xml_file.write("<charge>\n")
-    e0 = 2.396452e-04  # e^2 mol/(kcal A), permittivity of free space
-    charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5
-    for charge in charges:
-        xml_file.write("{}\n".format(charge / charge_factor))
-    xml_file.write("</charge>\n")
-    if forcefield:
-        pair_coeffs = list(
-            set(
-                (atom.type, atom.epsilon, atom.sigma)
-                for atom in structure.atoms
-            )
-        )
-        pair_coeffs.sort(key=lambda pair_type: pair_type[0])
-        xml_file.write("<pair_coeffs>\n")
-        for param_set in pair_coeffs:
-            xml_file.write(
-                "{}\t{:.4f}\t{:.4f}\n".format(
-                    param_set[0],
-                    param_set[1] / ref_energy,
-                    param_set[2] / ref_distance,
-                )
-            )
-        xml_file.write("</pair_coeffs>\n")
-
-
-def _write_bond_information(xml_file, structure, ref_distance, ref_energy):
-    """Write the bonds in the system.
-
-    Parameters
-    ----------
-    xml_file : file object
-        The file object of the hoomdxml file being written
-    structure : parmed.Structure
-        Parmed structure object
-    ref_distance : float, default=1.0
-        Reference distance for conversion to reduced units
-    ref_energy : float, default=1.0
-        Reference energy for conversion to reduced units
-    """
-    unique_bond_types = set()
-    xml_file.write("<bond>\n")
-    for bond in structure.bonds:
-        t1, t2 = bond.atom1.type, bond.atom2.type
-        if t1 == "" or t2 == "":
-            t1, t2 = bond.atom1.name, bond.atom2.name
-        t1, t2 = sorted([t1, t2])
-        try:
-            bond_type = ("-".join((t1, t2)), bond.type.k, bond.type.req)
-        except AttributeError:  # no forcefield applied, bond.type is None
-            bond_type = ("-".join((t1, t2)), 0.0, 0.0)
-        unique_bond_types.add(bond_type)
-        xml_file.write(
-            "{} {} {}\n".format(bond_type[0], bond.atom1.idx, bond.atom2.idx)
-        )
-    xml_file.write("</bond>\n")
-    xml_file.write("<bond_coeffs>\n")
-    xml_file.write("<!-- type k r_eq -->\n")
-    for bond_type, k, req in unique_bond_types:
-        xml_file.write(
-            "{} {} {}\n".format(
-                bond_type,
-                k * 2.0 / ref_energy * ref_distance**2.0,
-                req / ref_distance,
-            )
-        )
-    xml_file.write("</bond_coeffs>\n")
-
-
-def _write_angle_information(xml_file, structure, ref_energy):
-    """Write the angles in the system.
-
-    Parameters
-    ----------
-    xml_file : file object
-        The file object of the hoomdxml file being written
-    structure : parmed.Structure
-        Parmed structure object
-    ref_energy : float, default=1.0
-        Reference energy for conversion to reduced units
-    """
-    unique_angle_types = set()
-    xml_file.write("<angle>\n")
-    for angle in structure.angles:
-        t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type
-        t1, t3 = sorted([t1, t3])
-        angle_type = ("-".join((t1, t2, t3)), angle.type.k, angle.type.theteq)
-        unique_angle_types.add(angle_type)
-        xml_file.write(
-            "{} {} {} {}\n".format(
-                angle_type[0], angle.atom1.idx, angle.atom2.idx, angle.atom3.idx
-            )
-        )
-    xml_file.write("</angle>\n")
-    xml_file.write("<angle_coeffs>\n")
-    xml_file.write("<!-- type k theta_eq -->\n")
-    for angle_type, k, teq in unique_angle_types:
-        xml_file.write(
-            "{} {} {}\n".format(angle_type, k * 2.0 / ref_energy, radians(teq))
-        )
-    xml_file.write("</angle_coeffs>\n")
-
-
-def _write_dihedral_information(xml_file, structure, ref_energy):
-    """Write dihedrals in the system.
-
-    Parameters
-    ----------
-    xml_file : file object
-        The file object of the hoomdxml file being written
-    structure : parmed.Structure
-        Parmed structure object
-    ref_energy : float, default=1.0
-        Reference energy for conversion to reduced units
-    """
-    unique_dihedral_types = set()
-    xml_file.write("<dihedral>\n")
-    for dihedral in structure.rb_torsions:
-        t1, t2 = (dihedral.atom1.type, dihedral.atom2.type)
-        t3, t4 = dihedral.atom3.type, dihedral.atom4.type
-        if [t2, t3] == sorted([t2, t3]):
-            types_in_dihedral = "-".join((t1, t2, t3, t4))
-        else:
-            types_in_dihedral = "-".join((t4, t3, t2, t1))
-        dihedral_type = (
-            types_in_dihedral,
-            dihedral.type.c0,
-            dihedral.type.c1,
-            dihedral.type.c2,
-            dihedral.type.c3,
-            dihedral.type.c4,
-            dihedral.type.c5,
-            dihedral.type.scee,
-            dihedral.type.scnb,
-        )
-        unique_dihedral_types.add(dihedral_type)
-        xml_file.write(
-            "{} {} {} {} {}\n".format(
-                dihedral_type[0],
-                dihedral.atom1.idx,
-                dihedral.atom2.idx,
-                dihedral.atom3.idx,
-                dihedral.atom4.idx,
-            )
-        )
-    xml_file.write("</dihedral>\n")
-    xml_file.write("<dihedral_coeffs>\n")
-    xml_file.write("<!-- type k1 k2 k3 k4 -->\n")
-    for (
-        dihedral_type,
-        c0,
-        c1,
-        c2,
-        c3,
-        c4,
-        c5,
-        scee,
-        scnb,
-    ) in unique_dihedral_types:
-        opls_coeffs = RB_to_OPLS(
-            c0, c1, c2, c3, c4, c5, error_if_outside_tolerance=False
-        )
-        opls_coeffs /= ref_energy
-        xml_file.write(
-            "{} {:.5f} {:.5f} {:.5f} {:.5f}\n".format(
-                dihedral_type,
-                opls_coeffs[1],
-                opls_coeffs[2],
-                opls_coeffs[3],
-                opls_coeffs[4],
-            )
-        )
-    xml_file.write("</dihedral_coeffs>\n")
-
-
-def _write_rigid_information(xml_file, rigid_bodies):
-    """Write rigid body information.
-
-    Parameters
-    ----------
-    xml_file : file object
-        The file object of the hoomdxml file being written
-    rigid_bodies : list, len=n_particles
-        The rigid body that each particle belongs to (-1 for none)
-    """
-    if not all(body is None for body in rigid_bodies):
-        xml_file.write("<body>\n")
-        for body in rigid_bodies:
-            if body is None:
-                body = -1
-            xml_file.write("{}\n".format(int(body)))
-        xml_file.write("</body>\n")
-
-
-def _write_box_information(xml_file, structure, ref_distance):
-    """Write box information.
-
-    Parameters
-    ----------
-    xml_file : file object
-        The file object of the hoomdxml file being written
-    structure : parmed.Structure
-        Parmed structure object
-    ref_energy : float, default=1.0
-        Reference energy for conversion to reduced units
-    """
-    if np.allclose(structure.box[3:6], np.array([90, 90, 90])):
-        box_str = '<box units="sigma"  Lx="{}" Ly="{}" Lz="{}"/>\n'
-        xml_file.write(box_str.format(*structure.box[:3] / ref_distance))
-    else:
-        a, b, c = structure.box[0:3] / ref_distance
-        alpha, beta, gamma = np.radians(structure.box[3:6])
-
-        lx = a
-        xy = b * np.cos(gamma)
-        xz = c * np.cos(beta)
-        ly = np.sqrt(b**2 - xy**2)
-        yz = (b * c * np.cos(alpha) - xy * xz) / ly
-        lz = np.sqrt(c**2 - xz**2 - yz**2)
-        box_str = '<box units="sigma"  Lx="{}" Ly="{}" Lz="{}" xy="{}" xz="{}" yz="{}"/>\n'
-        xml_file.write(box_str.format(lx, ly, lz, xy, xz, yz))
diff --git a/mbuild/formats/lammpsdata.py b/mbuild/formats/lammpsdata.py
deleted file mode 100644
index 5cd5846fb..000000000
--- a/mbuild/formats/lammpsdata.py
+++ /dev/null
@@ -1,1522 +0,0 @@
-"""LAMMPS data format."""
-
-import itertools as it
-from collections import OrderedDict
-from warnings import warn
-
-import numpy as np
-from parmed import Structure
-from parmed.parameters import ParameterSet
-from scipy.constants import epsilon_0
-
-from mbuild import Box
-from mbuild.utils.conversion import RB_to_OPLS
-from mbuild.utils.decorators import deprecated
-from mbuild.utils.orderedset import OrderedSet
-from mbuild.utils.sorting import natural_sort
-
-__all__ = ["write_lammpsdata"]
-# Define constants for conversions
-KCAL_TO_KJ = 4.184
-NM2_TO_ANG2 = 100.0
-NM_TO_ANG = 10.0
-ELEM_TO_COUL = 1.602176634e-19
-
-
-dep_msg = """
-Support for writing out LAMMPS data files will be removed
-in mbuild 1.0.
-See GMSO (https://github.com/mosdef-hub/gmso/tree/main/gmso/formats/lammpsdata) for
-continued support for LAMMPS.
-"""
-print(dep_msg)
-
-
-@deprecated(dep_msg)
-def write_lammpsdata(
-    structure,
-    filename,
-    atom_style="full",
-    unit_style="real",
-    mins=None,
-    maxs=None,
-    pair_coeff_label=None,
-    detect_forcefield_style=True,
-    nbfix_in_data_file=True,
-    use_urey_bradleys=False,
-    use_rb_torsions=True,
-    use_dihedrals=False,
-    sigma_conversion_factor=None,
-    epsilon_conversion_factor=None,
-    mass_conversion_factor=None,
-    charge_conversion_factor=True,
-    zero_dihedral_weighting_factor=False,
-    moleculeID_offset=1,
-):
-    """Output a LAMMPS data file.
-
-    Outputs a LAMMPS data file in the 'full' atom style format. Default units
-    are 'real' units. See http://lammps.sandia.gov/doc/atom_style.html for
-    more information on atom styles.
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd structure object
-    filename : str
-        Path of the output file
-    atom_style: str, optional, default='full'
-        Defines the style of atoms to be saved in a LAMMPS data file. The
-        following atom styles are currently supported:
-        'full', 'atomic', 'charge', 'molecular'
-        see http://lammps.sandia.gov/doc/atom_style.html for more information
-        on atom styles.
-    unit_style: str, optional, default='real'
-        Defines to unit style to be save in a LAMMPS data file.  Defaults to
-        'real' units. Current styles are supported: 'real', 'lj', 'metal'
-        see https://lammps.sandia.gov/doc/99/units.html for more information
-        on unit styles
-    mins : list, optional, default=None
-        Minimum box dimension in x, y, z directions, nm
-    maxs : list, optional, default=None
-        Maximum box dimension in x, y, z directions, nm
-    pair_coeff_label : str, optional, default=None
-        Provide a custom label to the pair_coeffs section in the lammps data
-        file. A value of None means a suitable default will be chosen.
-    detect_forcefield_style : bool, optional, default=True
-        If True, format lammpsdata parameters based on the contents of
-        the parmed Structure
-    use_urey_bradleys : bool, optional, default=False
-        If True, will treat angles as CHARMM-style angles with urey bradley
-        terms while looking for `structure.urey_bradleys`
-    use_rb_torsions : bool, optional, default=True
-        If True, will treat dihedrals OPLS-style torsions while looking for
-        `structure.rb_torsions`
-    use_dihedrals : bool, optional, default=False
-        If True, will treat dihedrals as CHARMM-style dihedrals while looking
-        for `structure.dihedrals`
-    zero_dihedral_weighting_factor : bool, optional, default=False
-        If True, will set weighting parameter to zero in CHARMM-style dihedrals.
-        This should be True if the CHARMM dihedral style is used in non-CHARMM forcefields.
-    sigma_conversion_factor : float, optional, default=None
-        If unit_style is set to 'lj', then sigma conversion factor is used to non-dimensionalize.
-        Assume to be in units of nm. If None, will take the largest sigma value in
-        the structure.atoms.sigma values.
-    epsilon_conversion_factor : float, optional, default=None
-        If unit_style is set to 'lj', then epsilon conversion factor is used to non-dimensionalize.
-        Assume to be in units of kCal/mol. If None, will take the largest epsilon value in
-        the structure.atoms.epsilon values.
-    mass_conversion_factor : float, optional, default=None
-        If unit_style is set to 'lj', then mass conversion factor is used to non-dimensionalize.
-        Assume to be in units of amu. If None, will take the largest mass value in
-        the structure.atoms.mass values.
-    charge_conversion_factor : bool, optional, default=True
-        If unit_style is set to 'lj', then charge conversion factor may or may not be used to
-        non-dimensionalize. Assume to be in elementary charge units. If ``True``, the charges
-        are scaled by ``np.sqrt(4*np.pi()*eps_0*sigma_conversion_factor*epsilon_conversion_factor)``.
-        If ``False``, the charges are not scaled and the user must be wary to choose the dielectric
-        constant properly, which may be more convenient to implement an implicit solvent.
-    moleculeID_offset : int , optional, default=1
-        Since LAMMPS treats the MoleculeID as an additional set of information
-        to identify what molecule an atom belongs to, this currently
-        behaves as a residue id. This value needs to start at 1 to be
-        considered a real molecule.
-
-    Notes
-    -----
-    See http://lammps.sandia.gov/doc/2001/data_format.html for a full
-    description of the LAMMPS data format. Currently the following sections are
-    supported (in addition to the header): *Masses*, *Nonbond Coeffs*,
-    *Bond Coeffs*, *Angle Coeffs*, *Dihedral Coeffs*, *Atoms*, *Bonds*,
-    *Angles*, *Dihedrals*, *Impropers*
-    OPLS and CHARMM forcefield styles are supported, AMBER forcefield styles
-    are NOT
-
-    Some of this function has beed adopted from `mdtraj`'s support of the
-    LAMMPSTRJ trajectory format. See
-    https://github.com/mdtraj/mdtraj/blob/master/mdtraj/formats/lammpstrj.py
-    for details.
-
-    unique_types : a sorted list of unique atomtypes for all atoms in the structure.
-        Defined by:
-            atomtype : atom.type
-    unique_bond_types: an enumarated OrderedDict of unique bond types for all bonds in the structure.
-        Defined by bond parameters and component atomtypes, in order:
-        --- k : bond.type.k
-        --- req : bond.type.req
-        --- atomtypes : sorted((bond.atom1.type, bond.atom2.type))
-    unique_angle_types: an enumerated OrderedDict of unique angle types for all
-        angles in the structure.
-        Defined by angle parameters and component atomtypes, in order:
-        --- k : angle.type.k
-        --- theteq : angle.type.theteq
-        --- vertex atomtype: angle.atom2.type
-        --- atomtypes: sorted((bond.atom1.type, bond.atom3.type))
-
-    unique_dihedral_types: an enumerated OrderedDict of unique dihedrals type
-        for all dihedrals in the structure.
-        Defined by dihedral parameters and component atomtypes, in order:
-        --- c0 : dihedral.type.c0
-        --- c1 : dihedral.type.c1
-        --- c2 : dihedral.type.c2
-        --- c3 : dihedral.type.c3
-        --- c4 : dihedral.type.c4
-        --- c5 : dihedral.type.c5
-        --- scee : dihedral.type.scee
-        --- scnb : dihedral.type.scnb
-        --- atomtype 1 : dihedral.atom1.type
-        --- atomtype 2 : dihedral.atom2.type
-        --- atomtype 3 : dihedral.atom3.type
-        --- atomtype 4 : dihedral.atom4.type
-    """
-    if atom_style not in ["atomic", "charge", "molecular", "full"]:
-        raise ValueError(
-            'Atom style "{}" is invalid or is not currently supported'.format(
-                atom_style
-            )
-        )
-    if unit_style not in ["real", "lj", "metal"]:
-        raise ValueError(
-            'Unit style "{}" is invalid or is not currently supported'.format(
-                unit_style
-            )
-        )
-
-    forcefield = True
-    if structure[0].type == "":
-        forcefield = False
-    # copy structure so the input structure isn't modified in-place
-    structure = structure.copy(cls=Structure, split_dihedrals=True)
-
-    # units of angstroms
-    xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms])
-
-    if forcefield:
-        types = [atom.type for atom in structure.atoms]
-    else:
-        types = [atom.name for atom in structure.atoms]
-
-    unique_types = list(set(types))
-    unique_types.sort(key=natural_sort)
-
-    charges = np.array([atom.charge for atom in structure.atoms])
-
-    # lammps does not require the box to be centered at any a specific origin
-    # min and max dimensions are therefore needed to write the file in a
-    # consistent way the parmed structure only stores the box length.  It is
-    # not rigorous to assume bounds are 0 to L or -L/2 to L/2
-
-    # NOTE: 0 to L is current default, mins and maxs should be passed by user
-
-    if _check_minsmaxs(mins, maxs):
-        mins = np.array(mins) * NM_TO_ANG
-        maxs = np.array(maxs) * NM_TO_ANG
-        box = Box.from_mins_maxs_angles(
-            mins=mins, maxs=maxs, angles=structure.box[3:6]
-        )  # box lengths input in nm, so convert to angstrom
-    else:  # Parmed internally converts box to angstroms
-        box = Box(
-            lengths=np.array([val for val in structure.box[0:3]]),
-            angles=structure.box[3:6],
-        )
-
-        warn(
-            "Explicit box bounds (i.e., mins and maxs) were not provided. Box "
-            "bounds are assumed to be min = 0 and max = length in each "
-            "direction. This may not produce a system with the expected "
-            "spatial location and may cause non-periodic systems to fail. "
-            "Bounds can be defined explicitly by passing the them to the "
-            "write_lammpsdata function or by passing box info to the save "
-            "function."
-        )
-
-    # Lammps syntax depends on the functional form
-    # Infer functional form based on the properties of the structure
-    if detect_forcefield_style:
-        # Check angles
-        if len(structure.urey_bradleys) > 0:
-            print("Urey bradley terms detected, will use angle_style charmm")
-            use_urey_bradleys = True
-        else:
-            print(
-                "No urey bradley terms detected, will use angle_style harmonic"
-            )
-            use_urey_bradleys = False
-
-        # Check dihedrals
-        if len(structure.rb_torsions) > 0:
-            print("RB Torsions detected, will use dihedral_style opls")
-            use_rb_torsions = True
-        else:
-            use_rb_torsions = False
-        if len([d for d in structure.dihedrals if not d.improper]) > 0:
-            print("Charmm dihedrals detected, will use dihedral_style charmm")
-            use_dihedrals = True
-        else:
-            use_dihedrals = False
-    if use_rb_torsions and use_dihedrals:
-        raise ValueError(
-            "Multiple dihedral styles detected, check your "
-            "Forcefield XML and structure"
-        )
-
-    # save atom index information for all bonded params in structure
-    bonds = [[b.atom1.idx + 1, b.atom2.idx + 1] for b in structure.bonds]
-    angles = [
-        [angle.atom1.idx + 1, angle.atom2.idx + 1, angle.atom3.idx + 1]
-        for angle in structure.angles
-    ]
-    if use_rb_torsions:
-        dihedrals = [
-            [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1]
-            for d in structure.rb_torsions
-        ]
-    elif use_dihedrals:
-        dihedrals = [
-            [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1]
-            for d in structure.dihedrals
-            if not d.improper
-        ]
-    else:
-        dihedrals = []
-    impropers = [
-        [i.atom1.idx + 1, i.atom2.idx + 1, i.atom3.idx + 1, i.atom4.idx + 1]
-        for i in structure.impropers
-    ]
-    imp_dihedrals = [
-        [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1]
-        for d in structure.dihedrals
-        if d.improper
-    ]
-
-    if impropers and imp_dihedrals:
-        raise ValueError("Use of multiple improper styles is not supported")
-
-    # Get lj conversion factors if they exist and apply lj params to charges and box
-    # Params are in read in mbuild units of angstrom, kcal/mol, amu and converted to lammps units
-    if unit_style == "lj":
-        sigma_conversion_factor = _evaluate_lj_conversion_factors(
-            structure, "sigma", sigma_conversion_factor
-        )
-        epsilon_conversion_factor = _evaluate_lj_conversion_factors(
-            structure, "epsilon", epsilon_conversion_factor
-        )
-        mass_conversion_factor = _evaluate_lj_conversion_factors(
-            structure, "mass", mass_conversion_factor
-        )
-        # Convert coordinates and charges to LJ units
-        xyz = xyz / sigma_conversion_factor
-        if charge_conversion_factor:
-            charges = (charges * ELEM_TO_COUL) / np.sqrt(
-                4
-                * np.pi
-                * sigma_conversion_factor
-                * NM_TO_ANG**-1
-                * epsilon_conversion_factor
-                * KCAL_TO_KJ
-                * epsilon_0
-                * 10**-6
-            )
-        charges[np.isinf(charges)] = 0
-        eng_unit_str = " "
-
-    elif unit_style == "real":
-        sigma_conversion_factor = 1
-        epsilon_conversion_factor = 1
-        mass_conversion_factor = 1
-        eng_unit_str = "kcal/mol"
-
-    elif unit_style == "metal":
-        sigma_conversion_factor = 1  # unit in Angstrom
-        epsilon_conversion_factor = 1 / 0.043364115  #  kcal/mol to eV
-        mass_conversion_factor = 1
-        eng_unit_str = "eV"
-    else:
-        raise ValueError(
-            "Currently only support 'real', 'lj', and 'metal' unit_styles."
-        )
-
-    # Divide by conversion factor
-    Lx = box.Lx * (1 / sigma_conversion_factor)
-    Ly = box.Ly * (1 / sigma_conversion_factor)
-    Lz = box.Lz * (1 / sigma_conversion_factor)
-    box = Box(lengths=(Lx, Ly, Lz), angles=box.angles)
-
-    # Get bonded parameter information from structure
-    if bonds:
-        if len(structure.bond_types) == 0:
-            bond_types = np.ones(len(bonds), dtype=int)
-        else:
-            bond_types, unique_bond_types = _get_bond_types(
-                structure,
-                bonds,
-                unique_types,
-                sigma_conversion_factor,
-                epsilon_conversion_factor,
-            )
-
-    if angles:
-        angle_types, unique_angle_types = _get_angle_types(
-            structure,
-            unique_types,
-            use_urey_bradleys,
-            sigma_conversion_factor,
-            epsilon_conversion_factor,
-        )
-
-    if imp_dihedrals:
-        (
-            imp_dihedral_types,
-            unique_imp_dihedral_types,
-        ) = _get_improper_dihedral_types(
-            structure,
-            unique_types,
-            epsilon_conversion_factor,
-        )
-
-    if dihedrals:
-        dihedral_types, unique_dihedral_types = _get_dihedral_types(
-            structure,
-            unique_types,
-            use_rb_torsions,
-            use_dihedrals,
-            epsilon_conversion_factor,
-            zero_dihedral_weighting_factor,
-        )
-
-    if impropers:
-        improper_types, unique_improper_types = _get_impropers(
-            structure,
-            unique_types,
-            epsilon_conversion_factor,
-        )
-
-    # Write lammps data file https://docs.lammps.org/2001/data_format.html
-    with open(filename, "w") as data:
-        data.write(f"{filename} - created by mBuild; units = {unit_style}\n")
-        if unit_style == "lj":
-            data.write("#Normalization factors: ")
-            data.write(
-                "sigma - {:.3E} (angstrom), ".format(sigma_conversion_factor)
-            )
-            data.write(
-                "epsilon - {:.3E} (kcal/mol), ".format(
-                    epsilon_conversion_factor
-                )
-            )
-            data.write("mass - {:.3E} (amu)".format(mass_conversion_factor))
-        data.write("\n")
-
-        # Write counts of bonded interactions
-        data.write("{:d} atoms\n".format(len(structure.atoms)))
-        if atom_style in ["full", "molecular"]:
-            data.write("{:d} bonds\n".format(len(bonds)))
-            data.write("{:d} angles\n".format(len(angles)))
-            data.write("{:d} dihedrals\n".format(len(dihedrals)))
-            data.write(
-                "{:d} impropers\n\n".format(len(impropers) + len(imp_dihedrals))
-            )
-
-        # Write counts of unique bonded types
-        data.write("{:d} atom types\n".format(len(set(types))))
-        if atom_style in ["full", "molecular"]:
-            if bonds:
-                data.write("{:d} bond types\n".format(len(set(bond_types))))
-            if angles:
-                data.write("{:d} angle types\n".format(len(set(angle_types))))
-            if dihedrals:
-                data.write(
-                    "{:d} dihedral types\n".format(len(set(dihedral_types)))
-                )
-            if impropers:
-                data.write(
-                    "{:d} improper types\n".format(len(set(improper_types)))
-                )
-            elif imp_dihedrals:
-                data.write(
-                    "{:d} improper types\n".format(len(set(imp_dihedral_types)))
-                )
-
-        data.write("\n")
-        # Write box data
-        # NOTE: Needs better logic handling maxs and mins of a bounding box
-        # NOTE: CCC, "could be an option to grab mins and maxs of the structure boundingbox"
-        # Lammps supports any box origin, so it is okay if it matches the structure positions
-        _write_box_information(box, data, mins)
-
-        # Write mass data
-        _write_mass_information(
-            structure,
-            data,
-            mass_conversion_factor,
-            forcefield,
-            unique_types,
-            types,
-        )
-
-        if forcefield:
-            # Write pair coefficients data
-            _write_pair_information(
-                structure,
-                data,
-                forcefield,
-                sigma_conversion_factor,
-                epsilon_conversion_factor,
-                unique_types,
-                types,
-                pair_coeff_label,
-                unit_style,
-                nbfix_in_data_file,
-                eng_unit_str,
-            )
-
-            # Write bond coefficients
-            if bonds:
-                _write_bond_information(
-                    structure, data, unique_bond_types, unit_style, eng_unit_str
-                )
-
-            # Write angle coefficients
-            if angles:
-                _write_angle_information(
-                    structure,
-                    data,
-                    unique_angle_types,
-                    use_urey_bradleys,
-                    unit_style,
-                    eng_unit_str,
-                )
-
-            # Write dihedral coefficients
-            if dihedrals:
-                _write_dihedral_information(
-                    structure,
-                    data,
-                    unique_dihedral_types,
-                    unit_style,
-                    use_rb_torsions,
-                    use_dihedrals,
-                    eng_unit_str,
-                )
-
-            # Write improper coefficients
-            if impropers:
-                _write_improper_information(
-                    structure, data, unique_improper_types, unit_style
-                )
-                # Write improper dihedrals
-            elif imp_dihedrals:
-                _write_imp_dihedral_information(
-                    structure, data, unique_imp_dihedral_types, unit_style
-                )
-
-        # Write Atom data
-        data.write("\nAtoms # {}\n\n".format(atom_style))
-        if atom_style == "atomic":
-            atom_line = "{index:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
-        elif atom_style == "charge":
-            if unit_style in ["real", "metal"]:
-                atom_line = "{index:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
-            elif unit_style == "lj":
-                atom_line = "{index:d}\t{type_index:d}\t{charge:.4ef}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
-        elif atom_style == "molecular":
-            atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
-        elif atom_style == "full":
-            if unit_style in ["real", "metal"]:
-                atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
-            elif unit_style == "lj":
-                atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{charge:.4e}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
-
-        for i, coords in enumerate(xyz):
-            data.write(
-                atom_line.format(
-                    index=i + 1,
-                    type_index=unique_types.index(types[i]) + 1,
-                    zero=structure.atoms[i].residue.idx + moleculeID_offset,
-                    charge=charges[i],
-                    x=coords[0],
-                    y=coords[1],
-                    z=coords[2],
-                )
-            )
-
-        if atom_style in ["full", "molecular"]:
-            # Bond data
-            if bonds:
-                data.write("\nBonds\n\n")
-                for i, bond in enumerate(bonds):
-                    data.write(
-                        "{:d}\t{:d}\t{:d}\t{:d}\n".format(
-                            i + 1, bond_types[i], bond[0], bond[1]
-                        )
-                    )
-
-            # Angle data
-            if angles:
-                data.write("\nAngles\n\n")
-                for i, angle in enumerate(angles):
-                    data.write(
-                        "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format(
-                            i + 1, angle_types[i], angle[0], angle[1], angle[2]
-                        )
-                    )
-
-            # Dihedral data
-            if dihedrals:
-                data.write("\nDihedrals\n\n")
-                for i, dihedral in enumerate(dihedrals):
-                    data.write(
-                        "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format(
-                            i + 1,
-                            dihedral_types[i],
-                            dihedral[0],
-                            dihedral[1],
-                            dihedral[2],
-                            dihedral[3],
-                        )
-                    )
-            # Improper data for harmonic impropers
-            if impropers:
-                data.write("\nImpropers\n\n")
-                for i, improper in enumerate(impropers):
-                    # The atoms are written central-atom third in LAMMPS data file.
-                    # The first atom in structure.impropers is the central atom
-                    data.write(
-                        "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format(
-                            i + 1,
-                            improper_types[i],
-                            improper[2],
-                            improper[1],
-                            improper[0],
-                            improper[3],
-                        )
-                    )
-            elif imp_dihedrals:  # impropers in structure.dihedrals are cvff
-                data.write("\nImpropers\n\n")
-                for i, improper in enumerate(imp_dihedrals):
-                    # The atoms are written central-atom third in LAMMPS data file.
-                    # The third atom in structure.dihedral is the central atom
-                    # if dihedral.improper is True
-                    # This is correct for AMBER impropers even though
-                    # LAMMPS documentation implies central-atom-first.
-                    data.write(
-                        "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format(
-                            i + 1,
-                            imp_dihedral_types[i],
-                            improper[0],
-                            improper[1],
-                            improper[2],
-                            improper[3],
-                        )
-                    )
-
-
-def _evaluate_lj_conversion_factors(
-    structure, conversion_name, conversion_factor
-):
-    """Get Lennard Jones style conversion factors. `conversion_name`` can be sigma, epsilon, or mass."""
-    if conversion_factor is None:
-        # Check if structure is parametrized
-        if any([atom.sigma for atom in structure.atoms]) is None:
-            raise ValueError(
-                "LJ units specified but one or more atoms has undefined LJ "
-                "parameters."
-            )
-        else:
-            conversion_factor = np.max(
-                [getattr(atom, conversion_name) for atom in structure.atoms]
-            )
-            warn(
-                f"Assuming {conversion_name} conversion factor of "
-                + str(conversion_factor)
-            )
-        if conversion_factor == 0:
-            conversion_factor = 1
-            warn(
-                f"{conversion_name} conversion factor cannot be inferred from the maximum {conversion_name} value in the ParmEd Structure. "
-                "Setting the {conversion_name} conversion factor to 1"
-            )
-    elif conversion_factor <= 0:
-        raise ValueError(
-            f"The {conversion_name} conversion factor to convert to LJ units should be greater than 0."
-        )
-    else:
-        # assume conversion factor passed in mbuild units, convert to lammps real units
-        if conversion_name == "sigma":
-            conversion_factor *= NM_TO_ANG
-        elif conversion_name == "epsilon":
-            conversion_factor *= 1 / KCAL_TO_KJ
-        elif conversion_name == "mass":
-            pass
-    return conversion_factor
-
-
-def _check_minsmaxs(mins, maxs):
-    """Return True if both mins and maxs have been defined, and each have length 3 otherwise returns False."""
-    if mins and maxs:
-        if len(mins) == 3 and len(maxs) == 3:
-            return True
-        else:
-            warn(
-                "mins and maxs passed to write_lammpsdata, but list size is "
-                "incorrect. mins and maxs will be ignored."
-            )
-            return False
-    else:
-        return False
-
-
-def _get_bond_types(
-    structure,
-    bonds,
-    atom_types,
-    sigma_conversion_factor,
-    epsilon_conversion_factor,
-    bond_precision=3,
-):
-    """Will get the bond types from a parmed structure and convert them to lammps real units."""
-    unique_bond_types = OrderedDict(
-        enumerate(
-            OrderedSet(
-                *[
-                    (
-                        round(
-                            bond.type.k
-                            * (
-                                sigma_conversion_factor**2
-                                / epsilon_conversion_factor
-                            ),
-                            bond_precision,
-                        ),
-                        round(
-                            bond.type.req / sigma_conversion_factor,
-                            bond_precision,
-                        ),
-                        tuple(sorted((bond.atom1.type, bond.atom2.type))),
-                    )
-                    for bond in structure.bonds
-                ]
-            )
-        )
-    )
-
-    magnitude = np.ceil(np.log10(len(atom_types)))
-    bond_tuples = [x[2] for x in unique_bond_types.values()]
-    ordered_bond_tuples = bond_tuples.copy()
-    ordered_bond_tuples.sort(
-        key=lambda x: atom_types.index(x[0]) * 10**magnitude
-        + atom_types.index(x[1])
-    )
-
-    unique_bond_types = OrderedDict(
-        [
-            (unique_bond_types[bond_tuples.index(x)], i + 1)
-            for i, x in enumerate(ordered_bond_tuples)
-        ]
-    )
-
-    bond_types = [
-        unique_bond_types[
-            (
-                round(
-                    bond.type.k
-                    * (sigma_conversion_factor**2 / epsilon_conversion_factor),
-                    bond_precision,
-                ),
-                round(bond.type.req / sigma_conversion_factor, bond_precision),
-                tuple(sorted((bond.atom1.type, bond.atom2.type))),
-            )
-        ]
-        for bond in structure.bonds
-    ]
-    return bond_types, unique_bond_types
-
-
-def _get_angle_types(
-    structure,
-    atom_types,
-    use_urey_bradleys,
-    sigma_conversion_factor,
-    epsilon_conversion_factor,
-    angle_precision=3,
-):
-    """
-    Will get the angle types from a parmed structure and convert them to lammps real units.
-
-    Can get the parameters if urey bradleys or harmonic angles.
-    """
-    if use_urey_bradleys:
-        charmm_angle_types = []
-        for angle in structure.angles:
-            ub_k = 0
-            ub_req = 0
-            for ub in structure.urey_bradleys:
-                if (angle.atom1, angle.atom3) == (ub.atom1, ub.atom2):
-                    ub_k = ub.type.k
-                    ub_req = ub.type.req
-            charmm_angle_types.append(
-                (
-                    round(
-                        angle.type.k
-                        * (
-                            sigma_conversion_factor**2
-                            / epsilon_conversion_factor
-                        ),
-                        angle_precision,
-                    ),
-                    round(angle.type.theteq, angle_precision),
-                    round(ub_k / epsilon_conversion_factor, angle_precision),
-                    round(ub_req, angle_precision),
-                    tuple(sorted((angle.atom1.type, angle.atom3.type))),
-                )
-            )
-
-        unique_angle_types = OrderedDict(
-            enumerate(OrderedSet(*charmm_angle_types))
-        )
-
-        magnitude = np.ceil(np.log10(len(atom_types)))
-        angle_tuples = [x[-1] for x in unique_angle_types.values()]
-        ordered_angle_tuples = angle_tuples.copy()
-        ordered_angle_tuples.sort(
-            key=lambda x: atom_types.index(x[0]) * 10**magnitude
-            + atom_types.index(x[1])
-        )
-
-    else:
-        unique_angle_types = OrderedDict(
-            enumerate(
-                OrderedSet(
-                    *[
-                        (
-                            round(
-                                angle.type.k * (1 / epsilon_conversion_factor),
-                                angle_precision,
-                            ),
-                            round(angle.type.theteq, angle_precision),
-                            angle.atom2.type,
-                            tuple(sorted((angle.atom1.type, angle.atom3.type))),
-                        )
-                        for angle in structure.angles
-                    ]
-                )
-            )
-        )
-
-        magnitude = np.ceil(np.log10(len(atom_types)))
-        angle_tuples = [
-            (x[2], x[-1][0], x[-1][1]) for x in unique_angle_types.values()
-        ]
-        ordered_angle_tuples = angle_tuples.copy()
-        ordered_angle_tuples.sort(
-            key=lambda x: atom_types.index(x[0]) * 10 ** (2 * magnitude)
-            + atom_types.index(x[1]) * 10**magnitude
-            + atom_types.index(x[2])
-        )
-
-    unique_angle_types = OrderedDict(
-        [
-            (unique_angle_types[angle_tuples.index(x)], i + 1)
-            for i, x in enumerate(ordered_angle_tuples)
-        ]
-    )
-
-    if use_urey_bradleys:
-        angle_types = [
-            unique_angle_types[ub_info] for ub_info in charmm_angle_types
-        ]
-    else:
-        angle_types = [
-            unique_angle_types[
-                (
-                    round(
-                        angle.type.k * (1 / epsilon_conversion_factor),
-                        angle_precision,
-                    ),
-                    round(angle.type.theteq, angle_precision),
-                    angle.atom2.type,
-                    tuple(sorted((angle.atom1.type, angle.atom3.type))),
-                )
-            ]
-            for angle in structure.angles
-        ]
-
-    return angle_types, unique_angle_types
-
-
-typeLister = lambda x: (x.atom1.type, x.atom2.type, x.atom3.type, x.atom4.type)
-
-
-def _dihedral_sorter(dihedral):
-    if dihedral.atom2.type < dihedral.atom3.type or (
-        dihedral.atom2.type == dihedral.atom3.type
-        and dihedral.atom1.type < dihedral.atom4.type
-    ):
-        return typeLister(dihedral)
-    else:
-        return tuple(reversed(typeLister(dihedral)))
-
-
-def _get_dihedral_types(
-    structure,
-    atom_types,
-    use_rb_torsions,
-    use_dihedrals,
-    epsilon_conversion_factor,
-    zero_dihedral_weighting_factor,
-    dihedral_precision=5,
-):
-    """
-    Will get the dihedral types from a parmed structure and convert them to lammps real units.
-
-    Can be in the form of rb_torsions or charmm dihedrals.
-    """
-    lj_unit = 1.0 / epsilon_conversion_factor
-    if use_rb_torsions:
-        unique_dihedral_types = set(
-            [
-                (
-                    round(dihedral.type.c0 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c1 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c2 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c3 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c4 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c5 * lj_unit, dihedral_precision),
-                    round(dihedral.type.scee, 1),
-                    round(dihedral.type.scnb, 1),
-                    *_dihedral_sorter(dihedral),
-                )
-                for dihedral in structure.rb_torsions
-            ]
-        )
-        unique_dihedral_typesDict = {
-            v: k + 1
-            for k, v in enumerate(
-                sorted(
-                    unique_dihedral_types,
-                    key=lambda item: tuple(item[i] for i in [-3, -2, -4, -1]),
-                )
-            )
-        }
-
-    elif use_dihedrals:
-        structure.join_dihedrals()  # creates a list of dihedrals for dihedral.type
-        weight = 0.0
-        if not zero_dihedral_weighting_factor:
-            weight = 1.0
-        charmm_dihedralsList = []
-        for dihedral in structure.dihedrals:
-            if not dihedral.improper:
-                for dih_type in dihedral.type:
-                    charmm_dihedralsList.append(
-                        (
-                            dih_type.phi_k * lj_unit,
-                            int(round(dih_type.per, 0)),
-                            round(dih_type.phase, 3),
-                            round(weight / len(dihedral.type), 4),
-                            round(dih_type.scee, 1),
-                            round(dih_type.scnb, 1),
-                            *_dihedral_sorter(dihedral),
-                        )
-                    )
-        unique_dihedral_typesDict = {
-            v: k + 1
-            for k, v in enumerate(
-                sorted(
-                    charmm_dihedralsList,
-                    key=lambda item: tuple(item[i] for i in [-3, -2, -4, -1]),
-                )
-            )
-        }
-
-    if use_rb_torsions:
-        dihedral_types = [
-            unique_dihedral_typesDict[
-                (
-                    round(dihedral.type.c0 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c1 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c2 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c3 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c4 * lj_unit, dihedral_precision),
-                    round(dihedral.type.c5 * lj_unit, dihedral_precision),
-                    round(dihedral.type.scee, 1),
-                    round(dihedral.type.scnb, 1),
-                    *_dihedral_sorter(dihedral),
-                )
-            ]
-            for dihedral in structure.rb_torsions
-        ]
-    elif use_dihedrals:
-        dihedral_types = []
-        for dihedral in structure.dihedrals:
-            if not dihedral.improper:
-                for dih_type in dihedral.type:
-                    dihedral_types.append(
-                        unique_dihedral_typesDict[
-                            (
-                                dih_type.phi_k * lj_unit,
-                                int(round(dih_type.per, 0)),
-                                round(dih_type.phase, 3),
-                                round(weight / len(dihedral.type), 4),
-                                round(dih_type.scee, 1),
-                                round(dih_type.scnb, 1),
-                                *_dihedral_sorter(dihedral),
-                            )
-                        ]
-                    )
-
-    return dihedral_types, unique_dihedral_typesDict
-
-
-def _get_improper_dihedral_types(
-    structure, atom_types, epsilon_conversion_factor, imp_dih_precision=3
-):
-    """
-    Will get the improper types from a parmed structure and convert them to lammps real units.
-
-    Type cvff https://docs.lammps.org/improper_cvff.html
-    """
-    lj_unit = 1 / epsilon_conversion_factor
-    improper_dihedrals = []
-    for dihedral in structure.dihedrals:
-        if dihedral.improper:
-            dih_type = dihedral.type
-            phase = abs(int(round(dih_type.phase, 0)))
-            if not (phase == 0 or phase == 180):
-                raise ValueError("Improper dihedral phase must be 0 or 180")
-            if phase:
-                d = -1
-            else:
-                d = 1
-            improper_dihedrals.append(
-                (
-                    round(dih_type.phi_k * lj_unit, imp_dih_precision),
-                    d,
-                    int(round(dih_type.per, 0)),
-                    round(dih_type.scee, 1),
-                    round(dih_type.scnb, 1),
-                    dihedral.atom1.type,
-                    dihedral.atom2.type,
-                    dihedral.atom3.type,
-                    dihedral.atom4.type,
-                )
-            )
-    unique_imp_dihedral_types = dict(enumerate(OrderedSet(*improper_dihedrals)))
-    magnitude = np.ceil(np.log10(len(atom_types)))
-    dihedral_tuples = [
-        (x[1], x[5:]) for x in unique_imp_dihedral_types.values()
-    ]
-    ordered_dihedral_tuples = dihedral_tuples.copy()
-    ordered_dihedral_tuples.sort(
-        key=lambda x: atom_types.index(x[1][0]) * 10 ** (4 * magnitude)
-        + atom_types.index(x[1][1]) * 10 ** (3 * magnitude)
-        + atom_types.index(x[1][2]) * 10 ** (2 * magnitude)
-        + atom_types.index(x[1][3]) * 10**magnitude
-        + x[0]
-    )
-
-    unique_imp_dihedral_types = OrderedDict(
-        [
-            (unique_imp_dihedral_types[dihedral_tuples.index(x)], i + 1)
-            for i, x in enumerate(ordered_dihedral_tuples)
-        ]
-    )
-    imp_dihedral_types = [
-        unique_imp_dihedral_types[dihedral_info]
-        for dihedral_info in improper_dihedrals
-    ]
-
-    return imp_dihedral_types, unique_imp_dihedral_types
-
-
-def _get_impropers(
-    structure, atom_types, epsilon_conversion_factor, improper_precision=3
-):
-    """
-    Will get the improper types from a parmed structure and convert them to lammps real units.
-
-    Type harmonic https://docs.lammps.org/improper_harmonic.html.
-    """
-    lj_unit = 1 / epsilon_conversion_factor
-    unique_improper_types = dict(
-        enumerate(
-            OrderedSet(
-                *[
-                    (
-                        round(
-                            improper.type.psi_k * lj_unit, improper_precision
-                        ),
-                        round(improper.type.psi_eq, improper_precision),
-                        improper.atom3.type,
-                        improper.atom2.type,
-                        improper.atom1.type,
-                        improper.atom4.type,
-                    )
-                    for improper in structure.impropers
-                ]
-            )
-        )
-    )
-
-    magnitude = np.ceil(np.log10(len(atom_types)))
-    improper_tuples = [x[2:] for x in unique_improper_types.values()]
-    ordered_improper_tuples = improper_tuples.copy()
-    ordered_improper_tuples.sort(
-        key=lambda x: atom_types.index(x[0]) * 10 ** (4 * magnitude)
-        + atom_types.index(x[1]) * 10 ** (3 * magnitude)
-        + atom_types.index(x[2]) * 10 ** (2 * magnitude)
-        + atom_types.index(x[3]) * 10**magnitude
-    )
-
-    unique_improper_types = OrderedDict(
-        [
-            (unique_improper_types[improper_tuples.index(x)], i + 1)
-            for i, x in enumerate(ordered_improper_tuples)
-        ]
-    )
-
-    improper_types = [
-        unique_improper_types[
-            (
-                round(improper.type.psi_k * lj_unit, improper_precision),
-                round(improper.type.psi_eq, improper_precision),
-                improper.atom3.type,
-                improper.atom2.type,
-                improper.atom1.type,
-                improper.atom4.type,
-            )
-        ]
-        for improper in structure.impropers
-    ]
-
-    return improper_types, unique_improper_types
-
-
-def _write_box_information(box, data, mins):
-    """Write box information to lammps data file, can handle non-orthogonal boxes."""
-    if np.allclose(box.angles, 90.0, atol=1e-5) and (mins is None):
-        for i, dim in enumerate(["x", "y", "z"]):
-            data.write(
-                "{0:.6f} {1:.6f} {2}lo {2}hi\n".format(
-                    0.0,
-                    box.lengths[i],
-                    dim,
-                )
-            )
-    # NOTE:
-    # currently non-orthogonal bounding box translates
-    # Compound such that mins are new origin
-    else:
-        a = box.Lx
-        b = box.Ly
-        c = box.Lz
-        # alpha, beta, gamma = np.radians(box.angles)
-
-        xy = box.xy
-        xz = box.xz
-        yz = box.yz
-
-        # NOTE: using (0,0,0) as origin
-        xlo, ylo, zlo = (0.0, 0.0, 0.0)
-        xhi = xlo + a
-        yhi = ylo + b
-        zhi = zlo + c
-
-        xlo_bound = xlo + np.min([0.0, xy, xz, xy + xz])
-        xhi_bound = xhi + np.max([0.0, xy, xz, xy + xz])
-        ylo_bound = ylo + np.min([0.0, yz])
-        yhi_bound = yhi + np.max([0.0, yz])
-        zlo_bound = zlo
-        zhi_bound = zhi
-
-        data.write("{0:.6f} {1:.6f} xlo xhi\n".format(xlo_bound, xhi_bound))
-        data.write("{0:.6f} {1:.6f} ylo yhi\n".format(ylo_bound, yhi_bound))
-        data.write("{0:.6f} {1:.6f} zlo zhi\n".format(zlo_bound, zhi_bound))
-        data.write("{0:.6f} {1:.6f} {2:6f} xy xz yz\n".format(xy, xz, yz))
-
-
-def _write_mass_information(
-    structure, data, mass_conversion_factor, forcefield, unique_types, types
-):
-    """Write mass information from structure to lammps data file."""
-    if not forcefield:
-        masses = (
-            np.array([atom.mass for atom in structure.atoms])
-            / mass_conversion_factor
-        )
-    else:
-        tmp_masses = list()
-        for atom in structure.atoms:
-            # handle case where atomtype does not contain a mass
-            try:
-                tmp_masses.append(atom.atom_type.mass.in_units("g/mol").value)
-            except AttributeError:
-                warn(
-                    f"No mass or defined atomtype for atom: {atom}. Using atom mass of {atom.mass / mass_conversion_factor}"
-                )
-                tmp_masses.append(atom.mass)
-        masses = np.asarray(tmp_masses) / mass_conversion_factor
-
-    mass_dict = OrderedDict(
-        [
-            (unique_types.index(atom_type) + 1, mass)
-            for atom_type, mass in zip(types, masses)
-        ]
-    )
-    data.write("\nMasses\n\n")
-    for atom_type, mass in sorted(mass_dict.items()):
-        data.write(
-            "{:d}\t{:.6f}\t# {}\n".format(
-                atom_type, mass, unique_types[atom_type - 1]
-            )
-        )
-
-
-def _write_pair_information(
-    structure,
-    data,
-    forcefield,
-    sigma_conversion_factor,
-    epsilon_conversion_factor,
-    unique_types,
-    types,
-    pair_coeff_label,
-    unit_style,
-    nbfix_in_data_file,
-    eng_unit_str,
-):
-    """Write nonbonded pair information to lammps data file."""
-    epsilons = (
-        np.array([atom.epsilon for atom in structure.atoms])
-        / epsilon_conversion_factor
-    )
-    sigmas = (
-        np.array([atom.sigma for atom in structure.atoms])
-        / sigma_conversion_factor
-    )
-    forcefields = [atom.type for atom in structure.atoms]
-    epsilon_dict = dict(
-        [
-            (unique_types.index(atom_type) + 1, epsilon)
-            for atom_type, epsilon in zip(types, epsilons)
-        ]
-    )
-    sigma_dict = dict(
-        [
-            (unique_types.index(atom_type) + 1, sigma)
-            for atom_type, sigma in zip(types, sigmas)
-        ]
-    )
-    forcefield_dict = dict(
-        [
-            (unique_types.index(atom_type) + 1, forcefield)
-            for atom_type, forcefield in zip(types, forcefields)
-        ]
-    )
-
-    # Modified cross-interactions
-    if structure.has_NBFIX():
-        params = ParameterSet.from_structure(structure)
-        # Sort keys (maybe they should be sorted in ParmEd)
-        new_nbfix_types = OrderedDict()
-        for key in params.nbfix_types.keys():
-            sorted_key = tuple(sorted(key))
-            if sorted_key in new_nbfix_types:
-                warn("Sorted key matches an existing key")
-                if new_nbfix_types[sorted_key]:
-                    warn("nbfixes are not symmetric, overwriting old " "nbfix")
-            new_nbfix_types[sorted_key] = params.nbfix_types[key]
-        params.nbfix_types = new_nbfix_types
-        warn(
-            "Explicitly writing cross interactions using mixing rule: "
-            "{}".format(structure.combining_rule)
-        )
-        coeffs = OrderedDict()
-        for combo in it.combinations_with_replacement(unique_types, 2):
-            # Attempt to find pair coeffis in nbfixes
-            if combo in params.nbfix_types:
-                type1 = unique_types.index(combo[0]) + 1
-                type2 = unique_types.index(combo[1]) + 1
-                epsilon = params.nbfix_types[combo][0]  # kcal OR lj units
-                rmin = params.nbfix_types[combo][1]  # Angstrom OR lj units
-                sigma = rmin / 2 ** (1 / 6)
-                coeffs[(type1, type2)] = (
-                    round(sigma, 8),
-                    round(epsilon, 8),
-                )
-            else:
-                type1 = unique_types.index(combo[0]) + 1
-                type2 = unique_types.index(combo[1]) + 1
-                # Might not be necessary to be this explicit
-                if type1 == type2:
-                    sigma = sigma_dict[type1]
-                    epsilon = epsilon_dict[type1]
-                else:
-                    if structure.combining_rule == "lorentz":
-                        sigma = (sigma_dict[type1] + sigma_dict[type2]) * 0.5
-                    elif structure.combining_rule == "geometric":
-                        sigma = (sigma_dict[type1] * sigma_dict[type2]) ** 0.5
-                    else:
-                        raise ValueError(
-                            "Only lorentz and geometric combining "
-                            "rules are supported"
-                        )
-                    epsilon = (epsilon_dict[type1] * epsilon_dict[type2]) ** 0.5
-                coeffs[(type1, type2)] = (
-                    round(sigma, 8),
-                    round(epsilon, 8),
-                )
-        if nbfix_in_data_file:
-            if pair_coeff_label:
-                data.write("\nPairIJ Coeffs # {}\n".format(pair_coeff_label))
-            else:
-                data.write("\nPairIJ Coeffs # modified lj\n")
-
-            data.write(
-                "# type1 type2\tepsilon (%s)\tsigma (Angstrom)\n" % eng_unit_str
-            )
-
-            for (type1, type2), (sigma, epsilon) in coeffs.items():
-                data.write(
-                    "{0} \t{1} \t{2} \t\t{3}\t\t# {4}\t{5}\n".format(
-                        type1,
-                        type2,
-                        epsilon,
-                        sigma,
-                        forcefield_dict[type1],
-                        forcefield_dict[type2],
-                    )
-                )
-        else:
-            if pair_coeff_label:
-                data.write("\nPair Coeffs # {}\n".format(pair_coeff_label))
-            else:
-                data.write("\nPair Coeffs # lj\n")
-
-            for idx, epsilon in sorted(epsilon_dict.items()):
-                data.write(
-                    "{}\t{:.5f}\t{:.5f}\n".format(idx, epsilon, sigma_dict[idx])
-                )
-            print("Copy these commands into your input script:\n")
-            print(
-                "# type1 type2\tepsilon (%s)\tsigma (Angstrom)\n" % eng_unit_str
-            )
-            for (type1, type2), (sigma, epsilon) in coeffs.items():
-                print(
-                    "pair_coeff\t{0} \t{1} \t{2} \t\t{3} \t\t# {4} \t{5}".format(
-                        type1,
-                        type2,
-                        epsilon,
-                        sigma,
-                        forcefield_dict[type1],
-                        forcefield_dict[type2],
-                    )
-                )
-
-    # Pair coefficients
-    else:
-        if pair_coeff_label:
-            data.write("\nPair Coeffs # {}\n".format(pair_coeff_label))
-        else:
-            data.write("\nPair Coeffs # lj\n")
-
-        if unit_style in ["real", "metal"]:
-            data.write("#\tepsilon (%s)\t\tsigma (Angstrom)\n" % eng_unit_str)
-        elif unit_style == "lj":
-            data.write("#\treduced_epsilon \t\treduced_sigma \n")
-        for idx, epsilon in sorted(epsilon_dict.items()):
-            data.write(
-                "{}\t{:.5f}\t\t{:.5f}\t\t# {}\n".format(
-                    idx, epsilon, sigma_dict[idx], forcefield_dict[idx]
-                )
-            )
-
-
-def _write_bond_information(
-    structure, data, unique_bond_types, unit_style, eng_unit_str
-):
-    """Write Bond Coeffs section of lammps data file."""
-    data.write("\nBond Coeffs # harmonic\n")
-    if unit_style == "real" or unit_style == "metal":
-        data.write("#\tk(%s/angstrom^2)\t\treq(angstrom)\n" % eng_unit_str)
-    elif unit_style == "lj":
-        data.write("#\treduced_k\t\treduced_req\n")
-    sorted_bond_types = {
-        k: v
-        for k, v in sorted(unique_bond_types.items(), key=lambda item: item[1])
-    }
-    for params, idx in sorted_bond_types.items():
-        data.write(
-            "{}\t{}\t\t{}\t\t# {}\t{}\n".format(
-                idx,
-                params[0],
-                params[1],
-                params[2][0],
-                params[2][1],
-            )
-        )
-
-
-def _write_angle_information(
-    structure,
-    data,
-    unique_angle_types,
-    use_urey_bradleys,
-    unit_style,
-    eng_unit_str,
-):
-    """Write Angle Coeffs section of lammps data file."""
-    sorted_angle_types = {
-        k: v
-        for k, v in sorted(unique_angle_types.items(), key=lambda item: item[1])
-    }
-    if use_urey_bradleys:
-        data.write("\nAngle Coeffs # charmm\n")
-        data.write(
-            "#\tk(%s/rad^2)\t\ttheteq(deg)\tk(%s/angstrom^2)\treq(angstrom)\n"
-            % (eng_unit_str, eng_unit_str)
-        )
-        for params, idx in sorted_angle_types.items():
-            data.write("{}\t{}\t{:.5f}\t{:.5f}\t{:.5f}\n".format(idx, *params))
-
-    else:
-        data.write("\nAngle Coeffs # harmonic\n")
-
-        if unit_style == "lj":
-            data.write("#\treduced_k\t\ttheteq(deg)\n")
-        else:
-            data.write("#\tk(%s/rad^2)\t\ttheteq(deg)\n" % eng_unit_str)
-
-        for params, idx in sorted_angle_types.items():
-            data.write(
-                "{}\t{}\t\t{:.5f}\t# {}\t{}\t{}\n".format(
-                    idx,
-                    params[0],
-                    params[1],
-                    params[3][0],
-                    params[2],
-                    params[3][1],
-                )
-            )
-
-
-def _write_dihedral_information(
-    structure,
-    data,
-    unique_dihedral_types,
-    unit_style,
-    use_rb_torsions,
-    use_dihedrals,
-    eng_unit_str,
-):
-    """Write Dihedral Coeffs section of lammps data file."""
-    sorted_dihedral_types = {
-        k: v
-        for k, v in sorted(
-            unique_dihedral_types.items(),
-            key=lambda item: tuple(item[0][i] for i in [-3, -2, -4, -1]),
-        )
-    }
-    if use_rb_torsions:
-        data.write("\nDihedral Coeffs # opls\n")
-        if unit_style == "real" or unit_style == "metal":
-            data.write(
-                "#\tf1(%s)\tf2(%s)\tf3(%s)\tf4(%s)\n"
-                % (eng_unit_str, eng_unit_str, eng_unit_str, eng_unit_str)
-            )
-        elif unit_style == "lj":
-            data.write("#\tf1\tf2\tf3\tf4 (all lj reduced units)\n")
-        for idx, (params, _) in enumerate(sorted_dihedral_types.items()):
-            opls_coeffs = RB_to_OPLS(
-                params[0],
-                params[1],
-                params[2],
-                params[3],
-                params[4],
-                params[5],
-                error_if_outside_tolerance=False,
-            )
-            data.write(
-                "{}\t{:.5f}\t{:.5f}\t\t{:.5f}\t\t{:.5f}\t# {}\t{}\t{}\t{}\n".format(
-                    idx + 1,
-                    opls_coeffs[1],
-                    opls_coeffs[2],
-                    opls_coeffs[3],
-                    opls_coeffs[4],
-                    params[8],
-                    params[9],
-                    params[10],
-                    params[11],
-                )
-            )
-    elif use_dihedrals:
-        data.write("\nDihedral Coeffs # charmm\n")
-        data.write("#k, n, phi, weight\n")
-        for idx, (params, _) in enumerate(sorted_dihedral_types.items()):
-            data.write(
-                "{}\t{:.5f}\t{:d}\t{:d}\t{:.5f}\t# {}\t{}\t{}\t{}\n".format(
-                    idx + 1,
-                    params[0],
-                    params[1],
-                    int(params[2]),
-                    params[3],
-                    params[6],
-                    params[7],
-                    params[8],
-                    params[9],
-                )
-            )
-
-
-def _write_improper_information(
-    structure, data, unique_improper_types, unit_style
-):
-    """Write Impropers Coeffs section of lammps data file."""
-    sorted_improper_types = {
-        k: v
-        for k, v in sorted(
-            unique_improper_types.items(), key=lambda item: item[1]
-        )
-    }
-    data.write("\nImproper Coeffs # harmonic\n")
-    data.write("#k, phi\n")
-    for params, idx in sorted_improper_types.items():
-        data.write(
-            "{}\t{:.5f}\t{:.5f}\t# {}\t{}\t{}\t{}\n".format(
-                idx,
-                params[0],
-                params[1],
-                params[2],
-                params[3],
-                params[4],
-                params[5],
-            )
-        )
-
-
-def _write_imp_dihedral_information(
-    structure, data, unique_imp_dihedral_types, unit_style
-):
-    """Write Impropers Coeffs section of lammps data file."""
-    sorted_imp_dihedral_types = {
-        k: v
-        for k, v in sorted(
-            unique_imp_dihedral_types.items(),
-            key=lambda item: item[1],
-        )
-    }
-    data.write("\nImproper Coeffs # cvff\n")
-    data.write("#K, d, n\n")
-    for params, idx in sorted_imp_dihedral_types.items():
-        data.write(
-            "{}\t{:.5f}\t{:d}\t{:d}\t# {}\t{}\t{}\t{}\n".format(
-                idx,
-                params[0],
-                params[1],
-                params[2],
-                params[5],
-                params[6],
-                params[7],
-                params[8],
-            )
-        )
diff --git a/mbuild/formats/par_writer.py b/mbuild/formats/par_writer.py
deleted file mode 100644
index 26410efa9..000000000
--- a/mbuild/formats/par_writer.py
+++ /dev/null
@@ -1,201 +0,0 @@
-"""CHARMM Par format."""
-
-import warnings
-
-__all__ = ["write_par"]
-
-
-def write_par(structure, filename):
-    """Write CHARMM Par file given a parametrized structure.
-
-    Notes
-    -----
-    Follows format according to
-    https://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/
-    node25.html
-    Furthermore, ParmEd should support writing CHARMM par, rtf, str files
-    by converting the parmed.Structure into parmed.CharmmParameterSet
-
-    Parmed stores rmin/2 in "rmin"
-    """
-    # ATOMS
-    with open(filename, "w") as f:
-        f.write("ATOMS\n")
-        unique_atoms = set()
-        for atom in structure.atoms:
-            unique_atoms.add((atom.atom_type.name, atom.atom_type.mass))
-        for atom in unique_atoms:
-            f.write("MASS -1 {:8s} {:8.4f}\n".format(atom[0], atom[1]))
-
-        f.write("\nBONDS\n")
-        unique_bonds = set()
-        for bond in structure.bonds:
-            unique_bonds.add(
-                (
-                    bond.atom1.atom_type.name,
-                    bond.atom2.atom_type.name,
-                    bond.type,
-                )
-            )
-
-        for bond in unique_bonds:
-            f.write(
-                "{:8s} {:8s} {:.5f} {:.5f}\n".format(
-                    bond[0], bond[1], bond[2].k, bond[2].req
-                )
-            )
-
-        f.write("\nANGLES\n")
-        unique_angles = set()
-        unique_ubs = set()
-        for angle in structure.angles:
-            associated_ub = False
-            for ub in structure.urey_bradleys:
-                if ((angle.atom1, angle.atom3) == (ub.atom1, ub.atom2)) or (
-                    angle.atom3,
-                    angle.atom1,
-                ) == (ub.atom1, ub.atom2):
-                    unique_ubs.add(
-                        (
-                            angle.atom1.atom_type.name,
-                            angle.atom2.atom_type.name,
-                            angle.atom3.atom_type.name,
-                            angle.type,
-                            ub.type,
-                        )
-                    )
-                    associated_ub = True
-
-            if not associated_ub:
-                unique_angles.add(
-                    (
-                        angle.atom1.atom_type.name,
-                        angle.atom2.atom_type.name,
-                        angle.atom3.atom_type.name,
-                        angle.type,
-                    )
-                )
-
-        for ub in unique_ubs:
-            f.write(
-                "{:8s} {:8s} {:8s} {:.5f} {:.5f} {:.5f} {:.5f}\n".format(
-                    ub[0],
-                    ub[1],
-                    ub[2],
-                    ub[3].k,
-                    ub[3].theteq,
-                    ub[4].k,
-                    ub[4].req,
-                )
-            )
-        for angle in unique_angles:
-            f.write(
-                "{:8s} {:8s} {:8s} {:.5f} {:.5f}\n".format(
-                    angle[0], angle[1], angle[2], angle[3].k, angle[3].theteq
-                )
-            )
-
-        # These dihedrals need to be PeriodicTorsion Style (Charmm style)
-        if len(structure.rb_torsions) > 0:
-            warnings.warn("RB Torsions detected, but unsupported in par writer")
-        f.write("\nDIHEDRALS\n")
-        unique_dihedrals = set()
-        scnb = set()
-        for dihedral in structure.dihedrals:
-            if not dihedral.improper:
-                unique_dihedrals.add(
-                    (
-                        dihedral.atom1.atom_type.name,
-                        dihedral.atom2.atom_type.name,
-                        dihedral.atom3.atom_type.name,
-                        dihedral.atom4.atom_type.name,
-                        dihedral.type,
-                    )
-                )
-                scnb.add(dihedral.type.scnb)
-            else:
-                msg = (
-                    "AMBER-style improper detected between "
-                    + "{} {} {} {}".format(
-                        dihedral.atom1,
-                        dihedral.atom2,
-                        dihedral.atom3,
-                        dihedral.atom4,
-                    )
-                    + ", but unsupported in par writer"
-                )
-                warnings.warn(msg)
-
-        for dihedral in unique_dihedrals:
-            f.write(
-                "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format(
-                    dihedral[0],
-                    dihedral[1],
-                    dihedral[2],
-                    dihedral[3],
-                    dihedral[4].phi_k,
-                    dihedral[4].per,
-                    dihedral[4].phase,
-                )
-            )
-
-        f.write("\nIMPROPER\n")
-        unique_impropers = set()
-        for improper in structure.impropers:
-            unique_impropers.add(
-                (
-                    improper.atom1.atom_type.name,
-                    improper.atom2.atom_type.name,
-                    improper.atom3.atom_type.name,
-                    improper.atom4.atom_type.name,
-                    improper.type,
-                )
-            )
-        for improper in unique_impropers:
-            f.write(
-                "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format(
-                    improper[2],
-                    improper[0],
-                    improper[1],
-                    improper[3],
-                    improper[4].psi_k,
-                    0,
-                    improper[4].psi_eq,
-                )
-            )
-
-        sc_nb = [a for a in scnb]
-        if len(sc_nb) > 1:
-            warnings.warn(
-                "Multiple 1-4 LJ scalings were detected, "
-                "defaulting to first LJ scaling detected, {}".format(sc_nb[0])
-            )
-            sc_nb = sc_nb[0]
-        elif len(sc_nb) == 1:
-            sc_nb = sc_nb[0]
-        elif len(sc_nb) == 0:
-            warnings.warn("No 1-4 LJ scaling was detected, defaulting 1")
-            sc_nb = 1.0
-
-        f.write("\nNONBONDED\n")
-        unique_atypes = set()
-        for atom in structure.atoms:
-            unique_atypes.add(atom.atom_type)
-        for atype in unique_atypes:
-            # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4)
-            f.write(
-                "{:8s} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format(
-                    atype.name,
-                    0.0,
-                    -1 * atype.epsilon,
-                    atype.rmin,
-                    0.0,
-                    -1 * sc_nb * atype.epsilon,
-                    atype.rmin,
-                )
-            )
-
-        if structure.has_NBFIX():
-            warnings.warn("NBFixes detected but unsupported in par writer")
-
-        f.write("\nEND")
diff --git a/mbuild/formats/protobuf.py b/mbuild/formats/protobuf.py
deleted file mode 100644
index bc88198b7..000000000
--- a/mbuild/formats/protobuf.py
+++ /dev/null
@@ -1,195 +0,0 @@
-"""Protocol buffers (Protobuf) format.
-
-A language-agnostic data serialization format developed by Google
-https://developers.google.com/protocol-buffers
-"""
-
-import ele
-import numpy as np
-from google.protobuf.text_format import Merge, PrintMessage
-
-from mbuild import Box, Compound
-from mbuild.formats import compound_pb2
-from mbuild.utils.decorators import deprecated
-
-__all__ = ["write_pb2", "read_pb2"]
-
-
-dep_msg = """
-Support for the Protobuf file format will be removed in mBuild 1.0.
-"""
-print(dep_msg)
-
-
-@deprecated(dep_msg)
-def write_pb2(cmpd, filename, binary=True):
-    """Convert mb.Compound to Protobuf Message file.
-
-    Parameters
-    ----------
-    cmpd : mb.Compound
-    filename : str
-    binary: bool, default True
-        If True, will print a binary file
-        If False, will print to a text file
-    """
-    cmpd_to_proto = {}
-
-    root_proto = compound_pb2.Compound()
-    root_proto = _mb_to_proto(cmpd, root_proto)
-    cmpd_to_proto[cmpd] = root_proto
-
-    for sub_cmpd in cmpd.successors():
-        parent_cmpd = sub_cmpd.parent
-        sub_proto = cmpd_to_proto[parent_cmpd].children.add()
-        sub_proto = _mb_to_proto(sub_cmpd, sub_proto)
-        cmpd_to_proto[sub_cmpd] = sub_proto
-
-    _add_proto_bonds(cmpd, root_proto)
-
-    if binary:
-        with open(filename, "wb") as f:
-            f.write(root_proto.SerializeToString())
-    else:
-        with open(filename, "w") as f:
-            PrintMessage(root_proto, f)
-
-
-def read_pb2(filename, binary=True):
-    """Convert a Protobuf Message file into mb.Compound.
-
-    Parameters
-    ----------
-    filename : str
-    binary: bool, default True
-        If True, will print a binary file
-        If False, will print to a text file
-
-    Returns
-    -------
-    root_compound : mb.Compound
-    """
-    root_proto = compound_pb2.Compound()
-    if binary:
-        with open(filename, "rb") as f:
-            root_proto.ParseFromString(f.read())
-    else:
-        with open(filename, "r") as f:
-            Merge(f.read(), root_proto)
-
-    proto_to_cmpd = {}
-    root_compound = _proto_to_mb(root_proto)
-    proto_to_cmpd[root_proto.id] = root_compound
-
-    for sub_proto, parent_proto in _proto_successors(root_proto):
-        if parent_proto.id not in proto_to_cmpd:
-            parent_cmpd = _proto_to_mb(parent_proto)
-            proto_to_cmpd[parent_proto.id] = parent_cmpd
-        parent_cmpd = proto_to_cmpd[parent_proto.id]
-
-        if sub_proto.id not in proto_to_cmpd:
-            sub_cmpd = _proto_to_mb(sub_proto)
-            proto_to_cmpd[sub_proto.id] = sub_cmpd
-        sub_cmpd = proto_to_cmpd[sub_proto.id]
-
-        parent_cmpd.add(sub_cmpd)
-
-    _add_mb_bonds(root_proto, root_compound, proto_to_cmpd)
-    return root_compound
-
-
-def _mb_to_proto(cmpd, proto):
-    """Given mb.Compound, parse propertes into compound_pb2.Compound."""
-    proto.name = cmpd.name
-    proto.pos.x, proto.pos.y, proto.pos.z = cmpd.pos
-    proto.charge = cmpd.charge if cmpd.charge else 0.0
-    proto.id = id(cmpd)
-    (
-        proto.periodicity.x,
-        proto.periodicity.y,
-        proto.periodicity.z,
-    ) = cmpd.periodicity
-    if cmpd.element:
-        proto.element.name = cmpd.element.name
-        proto.element.symbol = cmpd.element.symbol
-        proto.element.atomic_number = cmpd.element.atomic_number
-        proto.element.mass = cmpd.element.mass
-
-    return proto
-
-
-def _add_proto_bonds(cmpd, proto):
-    """Parse the mb.Compound bonds, add to the proto bonds.
-
-    Parameters
-    ----------
-    cmpd : mb.Compound
-    proto : compound_pb2.Compound
-    """
-    for b in cmpd.bonds():
-        proto_bond = proto.bonds.add()
-        proto_bond.id1 = id(b[0])
-        proto_bond.id2 = id(b[1])
-
-
-def _proto_successors(proto):
-    """Recurisve method to look for a compound_pb2's children.
-
-    Parameters
-    ----------
-    proto : compound_pb2
-
-    Notes
-    -----
-    Base Case: there are no children to the proto, just return
-    Recursion: First look at proto's children and return these children
-        (sub_proto) Then make the recursive call to look at all the sub_proto's
-        successors. This is similar to mb.Compound().successors(). Unlike
-        mb.Compound(), we need to also keep track of parents in this recursion
-    """
-    if len(proto.children) == 0:
-        return
-    for sub_proto in proto.children:
-        yield (sub_proto, proto)
-        for sub_sub_proto, parent_proto in _proto_successors(sub_proto):
-            yield (sub_sub_proto, parent_proto)
-
-
-def _proto_to_mb(proto):
-    """Given compound_pb2.Compound, create mb.Compound.
-
-    Parameters
-    ----------
-    proto: compound_pb2.Compound()
-    """
-    if proto.element.symbol == "":
-        elem = None
-    else:
-        elem = ele.element_from_symbol(proto.element.symbol)
-    lengths = [proto.periodicity.x, proto.periodicity.y, proto.periodicity.z]
-    if np.all(np.array(lengths) != 0):
-        box = Box(lengths)
-    else:
-        box = None
-    return Compound(
-        name=proto.name,
-        pos=[proto.pos.x, proto.pos.y, proto.pos.z],
-        charge=proto.charge,
-        box=box,
-        element=elem,
-    )
-
-
-def _add_mb_bonds(proto, cmpd, proto_to_cmpd):
-    """Parse the compound_pb2.Compound bonds, add to mb.Compound.
-
-    Parameters
-    ----------
-    proto : compound_pb2.Compound
-    cmpd : mb.Compound
-    proto_to_cmpd : dict
-        keys : compound_pb2.Compound.id
-        value : mb.Compound
-    """
-    for bond in proto.bonds:
-        cmpd.add_bond([proto_to_cmpd[bond.id1], proto_to_cmpd[bond.id2]])
diff --git a/mbuild/formats/xyz.py b/mbuild/formats/xyz.py
deleted file mode 100644
index f4def4843..000000000
--- a/mbuild/formats/xyz.py
+++ /dev/null
@@ -1,134 +0,0 @@
-"""XYZ format."""
-
-from warnings import warn
-
-import numpy as np
-from ele import element_from_symbol
-from ele.exceptions import ElementError
-
-import mbuild as mb
-from mbuild.exceptions import MBuildError
-
-__all__ = ["read_xyz", "write_xyz"]
-
-
-def read_xyz(filename, compound=None):
-    """Read an XYZ file.
-
-    The expected format is as follows:
-    The first line contains the number of atoms in the file. The second line
-    contains a comment, which is not read. Remaining lines, one for each
-    atom in the file, include an elemental symbol followed by X, Y, and Z
-    coordinates in Angstroms. Columns are expected tbe separated by
-    whitespace. See https://openbabel.org/wiki/XYZ_(format).
-
-    Parameters
-    ----------
-    filename : str
-        Path of the input file
-
-    Returns
-    -------
-    compound : Compound
-
-    Notes
-    -----
-    The XYZ file format neglects many important details, including bonds,
-    residues, and box information.
-
-    There are some other flavors of the XYZ file format and not all are
-    guaranteed to be compatible with this reader. For example, the TINKER
-    XYZ format is not expected to be properly read.
-    """
-    if compound is None:
-        compound = mb.Compound()
-
-    guessed_elements = set()
-
-    compound_list = []
-    with open(filename, "r") as xyz_file:
-        n_atoms = int(xyz_file.readline())
-        xyz_file.readline()
-        coords = np.zeros(shape=(n_atoms, 3), dtype=np.float64)
-        for row, _ in enumerate(coords):
-            line = xyz_file.readline().split()
-            if not line:
-                msg = (
-                    "Incorrect number of lines in input file. Based on the "
-                    "number in the first line of the file, {} rows of atoms "
-                    "were expected, but at least one fewer was found."
-                )
-                raise MBuildError(msg.format(n_atoms))
-            coords[row] = line[1:4]
-            coords[row] *= 0.1
-            name = line[0]
-            try:
-                element = name.capitalize()
-                element = element_from_symbol(element)
-            except ElementError:
-                if name not in guessed_elements:
-                    warn(
-                        "No matching element found for {}; the particle will "
-                        "be added to the compound without an element "
-                        "attribute.".format(name)
-                    )
-                    guessed_elements.add(name)
-                element = None
-
-            particle = mb.Compound(pos=coords[row], name=name, element=element)
-            compound_list.append(particle)
-
-        # Verify we have read the last line by ensuring the next line in blank
-        line = xyz_file.readline().split()
-        if line:
-            msg = (
-                "Incorrect number of lines in input file. Based on the "
-                "number in the first line of the file, {} rows of atoms "
-                "were expected, but at least one more was found."
-            )
-            raise MBuildError(msg.format(n_atoms))
-
-    compound.add(compound_list)
-
-    return compound
-
-
-def write_xyz(structure, filename, write_atomnames=False, prec=6):
-    """Output an XYZ file.
-
-    Parameters
-    ----------
-    structure : parmed.Structure
-        ParmEd structure object
-    filename : str
-        Path of the output file
-    write_atomnames : bool
-        Write the `atom.name` attribute of the parmed structure
-        to the first column of the xyz file rather than the element
-    prec : int
-        The number of decimal places to write to the output file
-
-    Notes
-    -----
-    Coordinates are written in Angstroms. By default, the first column is the
-    element. These choices follow the conventions for the XYZ file format.
-
-    The XYZ file format neglects many important details, notably as bonds,
-    residues, and box information.
-    """
-    if isinstance(structure, mb.Compound):
-        raise ValueError("Expected a ParmEd structure, got an mbuild.Compound")
-
-    xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms])
-    if write_atomnames:
-        names = [atom.name for atom in structure.atoms]
-    else:
-        names = [atom.element_name for atom in structure.atoms]
-
-    with open(filename, "w") as xyz_file:
-        xyz_file.write(str(len(structure.atoms)))
-        xyz_file.write("\n" + filename + " - created by mBuild\n")
-        for name, (x, y, z) in zip(names, xyz):
-            xyz_file.write(
-                f"{name} {x:11.{prec}f} {y:11.{prec}f} {z:11.{prec}f}\n"
-            )
diff --git a/mbuild/utils/io.py b/mbuild/utils/io.py
index 1bc274c46..7c2f2e0f0 100644
--- a/mbuild/utils/io.py
+++ b/mbuild/utils/io.py
@@ -159,20 +159,6 @@ class DelayImportError(ImportError, SkipTest):
 # conda install -c conda-forge freud
 """
 
-MESSAGES[
-    "protobuf"
-] = """
-The code at {filename}:{line_number} requires the "protobuf" package
-
-protobuf can be installed using:
-
-# conda install -c conda-forge protobuf
-
-or
-
-# pip install protobuf
-"""
-
 
 def import_(module):
     """Import a module and issue a nice message to stderr if it isn't installed.
@@ -273,6 +259,14 @@ def import_(module):
         raise DelayImportError(m)
 
 
+try:
+    import hoomd
+
+    has_hoomd = True
+    del hoomd
+except ImportError:
+    has_hoomd = False
+
 try:
     import intermol
 
@@ -344,14 +338,6 @@ def import_(module):
 except ImportError:
     has_py3Dmol = False
 
-try:
-    from google import protobuf
-
-    has_protobuf = True
-    del protobuf
-except ImportError:
-    has_protobuf = False
-
 try:
     import garnett
 


Test patch:
diff --git a/mbuild/tests/base_test.py b/mbuild/tests/base_test.py
index 3437917c5..e58394dd1 100644
--- a/mbuild/tests/base_test.py
+++ b/mbuild/tests/base_test.py
@@ -26,6 +26,7 @@ def methane(self):
     @pytest.fixture
     def benzene(self):
         benzene = mb.load("C1=CC=CC=C1", smiles=True)
+        benzene.name = "Benzene"
 
         return benzene
 
@@ -121,12 +122,6 @@ def sixpoints(self):
         molecule.generate_bonds("C", "C", 0.9, 1.1)
         return molecule
 
-    @pytest.fixture
-    def benzene(self):
-        compound = mb.load(get_fn("benzene.mol2"))
-        compound.name = "Benzene"
-        return compound
-
     @pytest.fixture
     def rigid_benzene(self):
         compound = mb.load(get_fn("benzene.mol2"))
diff --git a/mbuild/tests/test_cassandramcf.py b/mbuild/tests/test_cassandramcf.py
deleted file mode 100644
index f54f4027a..000000000
--- a/mbuild/tests/test_cassandramcf.py
+++ /dev/null
@@ -1,467 +0,0 @@
-import pytest
-from numpy import isclose
-
-import mbuild as mb
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import has_foyer
-
-
-@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-class TestCassandraMCF(BaseTest):
-    def test_not_parameterized(self, ethane):
-        with pytest.raises(
-            ValueError, match=r"MCF writing not supported without"
-        ):
-            ethane.save(
-                filename="ethane.mcf",
-                angle_style="harmonic",
-                dihedral_style="opls",
-            )
-
-    def test_invalid_structure(self, ethane):
-        with pytest.raises(ValueError, match=r"requires parmed structure"):
-            mb.formats.cassandramcf.write_mcf(
-                ethane,
-                "ethane.mcf",
-                angle_style="harmonic",
-                dihedral_style="opls",
-            )
-
-    def test_invalid_angle_style(self, ethane):
-        with pytest.raises(
-            ValueError, match=r"Invalid selection for angle_style"
-        ):
-            ethane.save(
-                filename="ethane-opls.mcf",
-                forcefield_name="oplsaa",
-                angle_style="harm",
-                dihedral_style="opls",
-            )
-
-    def test_invalid_dihedral_style(self, ethane):
-        with pytest.raises(
-            ValueError, match=r"Invalid selection for dihedral_style"
-        ):
-            ethane.save(
-                filename="ethane-opls.mcf",
-                forcefield_name="oplsaa",
-                angle_style="harmonic",
-                dihedral_style="op",
-            )
-
-    def test_dihedral_style_none(self, ethane):
-        ethane.save(
-            filename="ethane-opls.mcf",
-            forcefield_name="oplsaa",
-            angle_style="harmonic",
-            dihedral_style="none",
-        )
-
-        mcf_data = []
-        with open("ethane-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Dihedral_Info":
-                    dihedral_section_start = idx
-
-        assert mcf_data[dihedral_section_start + 2][5] == "none"
-
-    def test_no_dihedrals(self, methane):
-        methane.save(
-            filename="methane-opls.mcf",
-            forcefield_name="oplsaa",
-            angle_style="harmonic",
-            dihedral_style="none",
-        )
-
-        mcf_data = []
-        with open("methane-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Dihedral_Info":
-                    dihedral_section_start = idx
-
-        assert mcf_data[dihedral_section_start + 1][0] == "0"
-
-    def test_no_14(self, methane):
-        methane.save(
-            filename="methane-opls.mcf",
-            forcefield_name="oplsaa",
-            angle_style="harmonic",
-            dihedral_style="none",
-        )
-
-        mcf_data = []
-        with open("methane-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Intra_Scaling":
-                    intrascaling_section_start = idx
-
-        assert isclose(float(mcf_data[intrascaling_section_start + 1][2]), 0.0)
-        assert isclose(float(mcf_data[intrascaling_section_start + 2][2]), 0.0)
-
-    def test_infer_14(self, ethane):
-        ethane.save(
-            filename="ethane-opls.mcf",
-            forcefield_name="oplsaa",
-            angle_style="harmonic",
-            dihedral_style="none",
-        )
-
-        mcf_data = []
-        with open("ethane-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Intra_Scaling":
-                    intrascaling_section_start = idx
-
-        assert isclose(float(mcf_data[intrascaling_section_start + 1][2]), 0.5)
-        assert isclose(float(mcf_data[intrascaling_section_start + 2][2]), 0.5)
-
-    def test_unmatched_dihedral_style(self, ethane):
-        with pytest.raises(ValueError, match=r"but RB torsions found"):
-            ethane.save(
-                filename="ethane-opls.mcf",
-                forcefield_name="oplsaa",
-                angle_style="harmonic",
-                dihedral_style="charmm",
-            )
-
-    def test_unreasonable_lj14(self, ethane):
-        import foyer
-
-        oplsaa = foyer.Forcefield(name="oplsaa")
-        ethane = oplsaa.apply(ethane)
-        with pytest.raises(ValueError, match=r"Unreasonable value"):
-            mb.formats.cassandramcf.write_mcf(
-                ethane,
-                "ethane.mcf",
-                angle_style="harmonic",
-                dihedral_style="opls",
-                lj14=2.0,
-            )
-
-    def test_unreasonable_coul14(self, ethane):
-        import foyer
-
-        oplsaa = foyer.Forcefield(name="oplsaa")
-        ethane = oplsaa.apply(ethane)
-        with pytest.raises(ValueError, match=r"Unreasonable value"):
-            mb.formats.cassandramcf.write_mcf(
-                ethane,
-                "ethane.mcf",
-                angle_style="harmonic",
-                dihedral_style="opls",
-                coul14=-1.0,
-            )
-
-    def test_multiple_molecules(self, ethane):
-        n_ethane = 2
-        ethane.name = "Ethane"
-        filled = mb.fill_box(
-            ethane, n_compounds=n_ethane, box=[0, 0, 0, 4, 4, 4]
-        )
-        with pytest.raises(
-            ValueError, match=r"Not all components of the molecule"
-        ):
-            filled.save(
-                filename="box-ethane-opls.mcf",
-                forcefield_name="oplsaa",
-                angle_style="harmonic",
-                dihedral_style="opls",
-            )
-
-    def test_save_forcefield(self, ethane):
-        ethane.save(
-            filename="ethane-opls.mcf",
-            forcefield_name="oplsaa",
-            angle_style="harmonic",
-            dihedral_style="opls",
-        )
-
-        mcf_data = []
-        with open("ethane-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Atom_Info":
-                    atom_section_start = idx
-                elif line[1] == "Bond_Info":
-                    bond_section_start = idx
-                elif line[1] == "Angle_Info":
-                    angle_section_start = idx
-                elif line[1] == "Dihedral_Info":
-                    dihedral_section_start = idx
-                elif line[1] == "Improper_Info":
-                    improper_section_start = idx
-                elif line[1] == "Fragment_Info":
-                    fragment_section_start = idx
-                elif line[1] == "Fragment_Connectivity":
-                    fragment_conn_start = idx
-
-        # Check a some atom info
-        assert mcf_data[atom_section_start + 1][0] == "8"
-        assert mcf_data[atom_section_start + 2][1] == "opls_135"
-        assert isclose(float(mcf_data[atom_section_start + 2][3]), 12.011)
-        assert isclose(float(mcf_data[atom_section_start + 2][4]), -0.18)
-        assert mcf_data[atom_section_start + 2][5] == "LJ"
-        assert isclose(float(mcf_data[atom_section_start + 2][6]), 33.21249)
-        assert isclose(float(mcf_data[atom_section_start + 2][7]), 3.500)
-
-        # Bond info
-        assert mcf_data[bond_section_start + 1][0] == "7"
-        passed_test = False
-        for line in mcf_data[bond_section_start + 2 : angle_section_start - 6]:
-            a1 = line[1]
-            a2 = line[2]
-            if (a1 == "1" and a2 == "2") or (a2 == "1" and a1 == "2"):
-                assert line[3] == "fixed"
-                assert isclose(float(line[4]), 1.090)
-                passed_test = True
-                break
-        assert passed_test
-
-        # Angle info
-        assert mcf_data[angle_section_start + 1][0] == "12"
-        passed_test = False
-        for line in mcf_data[
-            angle_section_start + 2 : dihedral_section_start - 8
-        ]:
-            if line[2] == "5":
-                a1 = line[1]
-                a3 = line[3]
-                if (a1 == "1" and a3 == "6") or (a3 == "1" and a1 == "6"):
-                    assert line[4] == "harmonic"
-                    assert isclose(float(line[5]), 18870.7)
-                    assert isclose(float(line[6]), 110.7)
-                    passed_test = True
-                    break
-        assert passed_test
-
-        # Dihedral info
-        assert mcf_data[dihedral_section_start + 1][0] == "9"
-        passed_test = False
-        for line in mcf_data[
-            dihedral_section_start + 2 : improper_section_start - 5
-        ]:
-            a1 = line[1]
-            a2 = line[2]
-            a3 = line[3]
-            a4 = line[4]
-            if (a1 == "2" and a2 == "1" and a3 == "5" and a4 == "6") or (
-                a4 == "2" and a3 == "1" and a2 == "5" and a1 == "6"
-            ):
-                assert line[5] == "OPLS"
-                assert isclose(float(line[6]), 0.000)
-                assert isclose(float(line[7]), 0.000)
-                assert isclose(float(line[8]), 0.000)
-                assert isclose(float(line[9]), 0.628)
-
-        assert mcf_data[improper_section_start + 1][0] == "0"
-        assert mcf_data[fragment_section_start + 1][0] == "2"
-
-        # Check fragment connectivity
-        assert mcf_data[fragment_conn_start + 1][0] == "1"
-        assert mcf_data[fragment_conn_start + 2][0] == "1"
-        assert mcf_data[fragment_conn_start + 2][1] == "1"
-        assert mcf_data[fragment_conn_start + 2][2] == "2"
-
-    def test_save_ring_forcefield(self, benzene):
-        benzene.save(
-            filename="benzene-opls.mcf",
-            forcefield_name="oplsaa",
-            angle_style="fixed",
-            dihedral_style="opls",
-        )
-
-        mcf_data = []
-        with open("benzene-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Atom_Info":
-                    atom_section_start = idx
-                elif line[1] == "Bond_Info":
-                    bond_section_start = idx
-                elif line[1] == "Angle_Info":
-                    angle_section_start = idx
-                elif line[1] == "Dihedral_Info":
-                    dihedral_section_start = idx
-                elif line[1] == "Improper_Info":
-                    improper_section_start = idx
-                elif line[1] == "Fragment_Info":
-                    fragment_section_start = idx
-                elif line[1] == "Fragment_Connectivity":
-                    fragment_conn_start = idx
-
-        # Check a some atom info
-        assert mcf_data[atom_section_start + 1][0] == "12"
-        assert mcf_data[atom_section_start + 2][1] == "opls_145"
-        assert isclose(float(mcf_data[atom_section_start + 2][3]), 12.011)
-        assert isclose(float(mcf_data[atom_section_start + 2][4]), -0.115)
-        assert mcf_data[atom_section_start + 2][5] == "LJ"
-        assert isclose(float(mcf_data[atom_section_start + 2][6]), 35.22537)
-        assert isclose(float(mcf_data[atom_section_start + 2][7]), 3.550)
-        assert mcf_data[atom_section_start + 2][8] == "ring"
-
-        # Bond info
-        assert mcf_data[bond_section_start + 1][0] == "12"
-        passed_test = False
-        for line in mcf_data[bond_section_start + 2 : angle_section_start - 6]:
-            a1 = line[1]
-            a2 = line[2]
-            if (a1 == "1" and a2 == "2") or (a2 == "1" and a1 == "2"):
-                assert line[3] == "fixed"
-                assert isclose(float(line[4]), 1.400)
-                passed_test = True
-                break
-        assert passed_test
-
-        # Angle info
-        assert mcf_data[angle_section_start + 1][0] == "18"
-        passed_test = False
-        for line in mcf_data[
-            angle_section_start + 2 : dihedral_section_start - 8
-        ]:
-            if line[2] == "2":
-                a1 = line[1]
-                a3 = line[3]
-                if (a1 == "1" and a3 == "3") or (a3 == "1" and a1 == "3"):
-                    assert line[4] == "fixed"
-                    assert isclose(float(line[5]), 120.00)
-                    passed_test = True
-                    break
-        assert passed_test
-
-        # Dihedral info
-        assert mcf_data[dihedral_section_start + 1][0] == "24"
-        passed_test = False
-        for line in mcf_data[
-            dihedral_section_start + 2 : improper_section_start - 5
-        ]:
-            a1 = line[1]
-            a2 = line[2]
-            a3 = line[3]
-            a4 = line[4]
-            if (a1 == "1" and a2 == "2" and a3 == "3" and a4 == "4") or (
-                a4 == "1" and a3 == "2" and a2 == "3" and a1 == "4"
-            ):
-                assert line[5] == "OPLS"
-                assert isclose(float(line[6]), 0.000)
-                assert isclose(float(line[7]), -0.000)
-                assert isclose(float(line[8]), 15.167)
-                assert isclose(float(line[9]), -0.000)
-
-        assert mcf_data[improper_section_start + 1][0] == "0"
-        assert mcf_data[fragment_section_start + 1][0] == "1"
-
-        # Check fragment connectivity
-        assert mcf_data[fragment_conn_start + 1][0] == "0"
-
-    def test_shorten_atomname(self, ethane):
-        import foyer
-
-        from mbuild.formats.cassandramcf import write_mcf
-
-        typed_ethane = foyer.forcefields.load_OPLSAA().apply(ethane)
-        typed_ethane[0].type = "C_very_very_very_extended"
-        write_mcf(
-            typed_ethane,
-            "ethane-opls.mcf",
-            angle_style="harmonic",
-            dihedral_style="opls",
-        )
-
-        mcf_data = []
-        with open("ethane-opls.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Atom_Info":
-                    atom_section_start = idx
-
-        assert mcf_data[atom_section_start + 2][1] == "y_very_very_extended"
-
-    def test_fused_rings(self):
-        import foyer
-
-        import mbuild
-        from mbuild.formats.cassandramcf import write_mcf
-
-        naph = mbuild.load("C1=CC=C2C=CC=CC2=C1", smiles=True)
-        # Note the atomtyping is wrong -- doesn't matter for test though
-        naph_ff = foyer.forcefields.load_OPLSAA().apply(naph)
-        write_mcf(
-            naph_ff, "naph.mcf", angle_style="harmonic", dihedral_style="opls"
-        )
-
-        mcf_data = []
-        with open("naph.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Fragment_Info":
-                    frag_section_start = idx
-
-        assert int(mcf_data[frag_section_start + 1][0]) == 1
-        assert int(mcf_data[frag_section_start + 2][1]) == 18
-
-        tmada = mbuild.load("C[C](C)(C)C12CC3CC(C1)CC(C3)C2", smiles=True)
-        tmada_ff = foyer.forcefields.load_OPLSAA().apply(tmada)
-        write_mcf(
-            tmada_ff, "tmada.mcf", angle_style="harmonic", dihedral_style="opls"
-        )
-
-        mcf_data = []
-        with open("tmada.mcf") as f:
-            for line in f:
-                mcf_data.append(line.strip().split())
-
-        for idx, line in enumerate(mcf_data):
-            if len(line) > 1:
-                if line[1] == "Fragment_Info":
-                    frag_section_start = idx
-
-        assert int(mcf_data[frag_section_start + 1][0]) == 5
-        assert int(mcf_data[frag_section_start + 2][1]) == 26
-        assert int(mcf_data[frag_section_start + 3][1]) == 5
-        assert int(mcf_data[frag_section_start + 4][1]) == 5
-        assert int(mcf_data[frag_section_start + 5][1]) == 5
-        assert int(mcf_data[frag_section_start + 6][1]) == 5
-
-    def test_infer_14_scaling_zero_eps(self):
-        import foyer
-
-        import mbuild
-        from mbuild.formats.cassandramcf import write_mcf
-
-        mol = mbuild.load("CO", smiles=True)
-        mol_ff = foyer.forcefields.load_OPLSAA().apply(mol)
-        with pytest.warns(UserWarning):
-            write_mcf(
-                mol_ff,
-                "test.mcf",
-                angle_style="harmonic",
-                dihedral_style="opls",
-            )
diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py
index 6d162c605..376f1eb6f 100644
--- a/mbuild/tests/test_compound.py
+++ b/mbuild/tests/test_compound.py
@@ -16,6 +16,7 @@
     get_fn,
     has_foyer,
     has_freud,
+    has_hoomd,
     has_intermol,
     has_mdtraj,
     has_networkx,
@@ -89,14 +90,11 @@ def test_load_conversion(self, ethane, h2o):
         assert np.allclose(test_converted1.xyz, test_converted2.xyz)
 
     def test_load_xyz(self):
-        class MyCompound(Compound):
-            def __init__(self):
-                super(MyCompound, self).__init__()
-
-                mb.load(get_fn("ethane.xyz"), compound=self)
-
-        myethane = MyCompound()
+        myethane = mb.load(get_fn("ethane.xyz"))
         assert myethane.n_particles == 8
+        assert myethane.n_bonds == 0
+        assert len(myethane.children) == 1
+        assert set([p.name for p in myethane.particles()]) == {"C", "H"}
 
     def test_update_from_file(self, ch3):
         ch3.update_coordinates(get_fn("methyl.pdb"))
@@ -373,12 +371,39 @@ def test_load_protein(self):
             # Chains info is lossed
             assert len(protein.children) == 393
 
-    def test_save_simple(self, ch3):
-        extensions = [".xyz", ".pdb", ".mol2", ".json", ".sdf"]
-        for ext in extensions:
-            outfile = "methyl_out" + ext
-            ch3.save(filename=outfile)
-            assert os.path.exists(outfile)
+    @pytest.mark.parametrize(
+        "extension",
+        [
+            ".xyz",
+            ".pdb",
+            ".mol2",
+            ".gro",
+            ".gsd",
+            ".json",
+            ".sdf",
+        ],
+    )
+    def test_save_simple(self, ch3, extension):
+        # Can't save gsd files with Windows
+        if extension == ".gsd" and not has_hoomd:
+            return True
+        outfile = "methyl_out" + extension
+        ch3.save(filename=outfile)
+        assert os.path.exists(outfile)
+
+    @pytest.mark.parametrize(
+        "extension",
+        [".xyz", ".pdb", ".mol2", ".json"],
+    )
+    def test_save_load_simple(self, ch3, extension):
+        outfile = "methyl_out" + extension
+        ch3.save(filename=outfile)
+        ch3_new = mb.load(outfile)
+        assert ch3_new.n_particles == ch3.n_particles
+        assert np.allclose(ch3_new.xyz, ch3.xyz)
+        ch3_names = [p.name for p in ch3.particles()]
+        ch3_new_names = [p.name for p in ch3_new.particles()]
+        assert set(ch3_names) == set(ch3_new_names)
 
     def test_save_json_loop(self, ethane):
         ethane.save("ethane.json", show_ports=True)
@@ -388,7 +413,7 @@ def test_save_json_loop(self, ethane):
         assert len(ethane.children) == len(ethane_copy.children)
 
     def test_save_box(self, ch3):
-        extensions = [".mol2", ".pdb", ".hoomdxml", ".gro", ".sdf"]
+        extensions = [".mol2", ".pdb", ".gro", ".sdf"]
         box_attributes = ["lengths"]
         custom_box = Box(lengths=[0.8, 0.8, 0.8], angles=[90, 90, 90])
         for ext in extensions:
@@ -404,7 +429,7 @@ def test_save_box(self, ch3):
                 assert np.array_equal(pad_attr, custom_attr)
 
     def test_save_overwrite(self, ch3):
-        extensions = [".gsd", ".hoomdxml", ".lammps", ".lmp", ".top", ".gro"]
+        extensions = [".mol2", ".xyz", ".gro"]
         for ext in extensions:
             outfile = "lyhtem" + ext
             ch3.save(filename=outfile)
@@ -412,86 +437,6 @@ def test_save_overwrite(self, ch3):
             with pytest.raises(IOError):
                 ch3.save(filename=outfile, overwrite=False)
 
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_save_forcefield(self, methane):
-        exts = [
-            ".gsd",
-            ".hoomdxml",
-            ".lammps",
-            ".lmp",
-            ".top",
-            ".gro",
-            ".mol2",
-            ".pdb",
-            ".xyz",
-            ".sdf",
-        ]
-        for ext in exts:
-            methane.save(
-                "lythem" + ext, forcefield_name="oplsaa", overwrite=True
-            )
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_save_forcefield_with_file(self, methane):
-        exts = [
-            ".gsd",
-            ".hoomdxml",
-            ".lammps",
-            ".lmp",
-            ".top",
-            ".gro",
-            ".mol2",
-            ".pdb",
-            ".xyz",
-            ".sdf",
-        ]
-        for ext in exts:
-            methane.save(
-                "lythem" + ext,
-                forcefield_files=get_fn("methane_oplssaa.xml"),
-                overwrite=True,
-            )
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    @pytest.mark.parametrize(
-        "ff_filename,foyer_kwargs",
-        [
-            ("ethane-angle-typo.xml", {"assert_angle_params": False}),
-            ("ethane-dihedral-typo.xml", {"assert_dihedral_params": False}),
-        ],
-    )
-    def test_save_missing_topo_params(self, ff_filename, foyer_kwargs):
-        """Test that the user is notified if not all topology parameters are found."""
-        from foyer.tests.utils import get_fn
-
-        ethane = mb.load(get_fn("ethane.mol2"))
-        with pytest.raises(Exception):
-            ethane.save("ethane.mol2", forcefield_files=get_fn(ff_filename))
-        with pytest.warns(UserWarning):
-            ethane.save(
-                "ethane.mol2",
-                forcefield_files=get_fn(ff_filename),
-                overwrite=True,
-                foyer_kwargs=foyer_kwargs,
-            )
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_save_forcefield_with_file_foyer_kwargs(self, methane):
-        foyer_kwargs = {"assert_improper_params": True}
-        with pytest.raises(Exception):
-            methane.save(
-                "lythem.hoomdxml",
-                forcefield_files=get_fn("methane_oplssaa.xml"),
-                overwrite=True,
-                foyer_kwargs=foyer_kwargs,
-            )
-        methane.save(
-            "lythem.hoomdxml",
-            forcefield_files=get_fn("methane_oplssaa.xml"),
-            overwrite=True,
-            foyer_kwargs={},
-        )
-
     def test_save_resnames(self, ch3, h2o):
         system = Compound([ch3, h2o])
         system.save("resnames.gro", residues=["CH3", "H2O"])
@@ -507,63 +452,6 @@ def test_save_resnames_single(self, c3, n4):
         assert struct.residues[0].number == 1
         assert struct.residues[1].number == 2
 
-    def test_save_residue_map(self, ethane):
-        filled = mb.fill_box(ethane, n_compounds=100, box=[0, 0, 0, 4, 4, 4])
-        t0 = time.time()
-        foyer_kwargs = {"use_residue_map": True}
-        filled.save(
-            "filled.mol2",
-            forcefield_name="oplsaa",
-            residues="Ethane",
-            foyer_kwargs=foyer_kwargs,
-        )
-        t1 = time.time()
-
-        foyer_kwargs = {"use_residue_map": False}
-        filled.save(
-            "filled.mol2",
-            forcefield_name="oplsaa",
-            overwrite=True,
-            residues="Ethane",
-            foyer_kwargs=foyer_kwargs,
-        )
-        t2 = time.time()
-        assert (t2 - t1) > (t1 - t0)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_save_references(self, methane):
-        foyer_kwargs = {"references_file": "methane.bib"}
-        methane.save(
-            "methyl.mol2", forcefield_name="oplsaa", foyer_kwargs=foyer_kwargs
-        )
-        assert os.path.isfile("methane.bib")
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_save_combining_rule(self, methane):
-        combining_rules = ["lorentz", "geometric"]
-        gmx_rules = {"lorentz": 2, "geometric": 3}
-        for combining_rule in combining_rules:
-            if combining_rule == "geometric":
-                methane.save(
-                    "methane.top",
-                    forcefield_name="oplsaa",
-                    combining_rule=combining_rule,
-                    overwrite=True,
-                )
-            else:
-                with pytest.warns(UserWarning):
-                    methane.save(
-                        "methane.top",
-                        forcefield_name="oplsaa",
-                        combining_rule=combining_rule,
-                        overwrite=True,
-                    )
-            with open("methane.top") as fp:
-                for i, line in enumerate(fp):
-                    if i == 18:
-                        gmx_rule = int(line.split()[1])
-                        assert gmx_rule == gmx_rules[combining_rule]
-
     def test_clone_with_box(self, ethane):
         ethane.box = ethane.get_boundingbox()
         ethane.periodicity = (True, True, False)
@@ -2672,9 +2560,9 @@ def test_elements_from_smiles(self, backend):
         for particle in mol.particles():
             assert particle.element is not None
 
-    def test_mins_maxs(self, benzene):
-        assert np.allclose(benzene.mins, [-0.2267, -0.15422, 0.0])
-        assert np.allclose(benzene.maxs, [0.20318, 0.34207, 0.0])
+    def test_mins_maxs(self, rigid_benzene):
+        assert np.allclose(rigid_benzene.mins, [-0.2267, -0.15422, 0.0])
+        assert np.allclose(rigid_benzene.maxs, [0.20318, 0.34207, 0.0])
 
     def test_periodicity_raises(self):
         with pytest.raises(ValueError):
diff --git a/mbuild/tests/test_gsd.py b/mbuild/tests/test_gsd.py
deleted file mode 100644
index 704a934a1..000000000
--- a/mbuild/tests/test_gsd.py
+++ /dev/null
@@ -1,320 +0,0 @@
-import numpy as np
-import pytest
-
-import mbuild as mb
-from mbuild import Box
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import has_foyer, has_gsd
-
-
-class TestGSD(BaseTest):
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_save(self, ethane):
-        ethane.save(filename="ethane.gsd")
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_save_forcefield(self, ethane):
-        ethane.save(filename="ethane-opls.gsd", forcefield_name="oplsaa")
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_save_box(self, ethane):
-        box = Box(lengths=[2.0, 2.0, 2.0], angles=[90.0, 90.0, 90.0])
-        ethane.save(
-            filename="ethane-box.gsd", forcefield_name="oplsaa", box=box
-        )
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_save_triclinic_box_(self, ethane):
-        box = Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[60, 70, 80])
-        ethane.save(
-            filename="triclinic-box.gsd", forcefield_name="oplsaa", box=box
-        )
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_particles(self, ethane):
-        from collections import OrderedDict
-
-        import gsd
-        import gsd.hoomd
-
-        from mbuild.utils.sorting import natural_sort
-
-        ethane.save(filename="ethane.gsd", forcefield_name="oplsaa")
-        with gsd.hoomd.open("ethane.gsd", mode="r") as f:
-            frame = f[0]
-
-        assert frame.configuration.step == 0
-        assert frame.configuration.dimensions == 3
-
-        mass_dict = {"C": 12.011, "H": 1.008}
-        masses = frame.particles.mass.astype(float)
-        for mass, particle in zip(masses, ethane.particles()):
-            assert round(mass, 3) == mass_dict[particle.name]
-
-        n_particles = frame.particles.N
-        assert n_particles == ethane.n_particles
-
-        positions = frame.particles.position.astype(float)
-        shift = positions[0] - (ethane[0].pos * 10)
-        shifted_xyz = (ethane.xyz * 10) + shift
-        assert np.array_equal(
-            np.round(positions, decimals=4), np.round(shifted_xyz, decimals=4)
-        )
-
-        opls_type_dict = OrderedDict([("C", "opls_135"), ("H", "opls_140")])
-        types_from_gsd = frame.particles.types
-        typeids_from_gsd = frame.particles.typeid.astype(int)
-        expected_types = [
-            opls_type_dict[particle.name] for particle in ethane.particles()
-        ]
-        unique_types = list(set(expected_types))
-        unique_types.sort(key=natural_sort)
-        expected_typeids = [
-            unique_types.index(atype) for atype in expected_types
-        ]
-        assert np.array_equal(types_from_gsd, unique_types)
-        assert np.array_equal(typeids_from_gsd, expected_typeids)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_box(self, ethane):
-        import gsd
-        import gsd.hoomd
-
-        lengths = [2.0, 3.0, 4.0]
-        ethane.box = Box(lengths=[2.0, 3.0, 4.0])
-
-        ethane.save(filename="ethane.gsd", forcefield_name="oplsaa")
-        with gsd.hoomd.open("ethane.gsd", mode="r") as f:
-            frame = f[0]
-
-        box_from_gsd = frame.configuration.box.astype(float)
-        (lx, ly, lz) = ethane.box.lengths
-        lx *= 10
-        ly *= 10
-        lz *= 10
-        assert np.array_equal(box_from_gsd[:3], [lx, ly, lz])
-        assert not np.any(box_from_gsd[3:])
-
-        ethane.periodicity = (True, True, True)
-        ethane.save(filename="ethane-periodicity.gsd", forcefield_name="oplsaa")
-        with gsd.hoomd.open("ethane-periodicity.gsd", mode="r") as f:
-            frame = f[0]
-        box_from_gsd_periodic = frame.configuration.box.astype(float)
-        assert np.array_equal(box_from_gsd, box_from_gsd_periodic)
-
-        box = Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[92, 104, 119])
-        # check that providing a box to save overwrites compound.box
-        ethane.save(
-            filename="triclinic-box.gsd", forcefield_name="oplsaa", box=box
-        )
-        with gsd.hoomd.open("triclinic-box.gsd", mode="r") as f:
-            frame = f[0]
-        lx, ly, lz, xy, xz, yz = frame.configuration.box
-
-        a = lx
-        b = np.sqrt(ly**2 + xy**2)
-        c = np.sqrt(lz**2 + xz**2 + yz**2)
-
-        assert np.isclose(np.cos(np.radians(92)), (xy * xz + ly * yz) / (b * c))
-        assert np.isclose(np.cos(np.radians(104)), xz / c)
-        assert np.isclose(np.cos(np.radians(119)), xy / b)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_rigid(self, benzene):
-        import gsd
-        import gsd.hoomd
-
-        benzene.label_rigid_bodies(rigid_particles="C")
-        benzene.save(filename="benzene.gsd", forcefield_name="oplsaa")
-        with gsd.hoomd.open("benzene.gsd", mode="r") as f:
-            frame = f[0]
-
-        rigid_bodies = frame.particles.body
-        expected_bodies = [
-            -1 if p.rigid_id is None else p.rigid_id
-            for p in benzene.particles()
-        ]
-        for gsd_body, expected_body in zip(rigid_bodies, expected_bodies):
-            assert gsd_body == expected_body
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_bonded(self, ethane):
-        import gsd
-        import gsd.hoomd
-        from foyer import Forcefield
-
-        ethane.save(filename="ethane.gsd", forcefield_name="oplsaa")
-        with gsd.hoomd.open("ethane.gsd", mode="r") as f:
-            frame = f[0]
-
-        structure = ethane.to_parmed()
-        forcefield = Forcefield(name="oplsaa")
-        structure = forcefield.apply(structure)
-
-        # Bonds
-        n_bonds = frame.bonds.N
-        assert n_bonds == len(structure.bonds)
-
-        expected_unique_bond_types = ["opls_135-opls_135", "opls_135-opls_140"]
-        bond_types = frame.bonds.types
-        assert np.array_equal(bond_types, expected_unique_bond_types)
-
-        bond_typeids = frame.bonds.typeid
-        bond_atoms = frame.bonds.group
-        expected_bond_atoms = [
-            [bond.atom1.idx, bond.atom2.idx] for bond in structure.bonds
-        ]
-        assert np.array_equal(bond_atoms, expected_bond_atoms)
-
-        bond_type_dict = {("C", "C"): 0, ("C", "H"): 1, ("H", "C"): 1}
-        expected_bond_typeids = []
-        for bond in structure.bonds:
-            expected_bond_typeids.append(
-                bond_type_dict[(bond.atom1.name, bond.atom2.name)]
-            )
-        assert np.array_equal(bond_typeids, expected_bond_typeids)
-
-        # Angles
-        n_angles = frame.angles.N
-        assert n_angles == len(structure.angles)
-
-        expected_unique_angle_types = [
-            "opls_135-opls_135-opls_140",
-            "opls_140-opls_135-opls_140",
-        ]
-        angle_types = frame.angles.types
-        assert np.array_equal(angle_types, expected_unique_angle_types)
-
-        angle_typeids = frame.angles.typeid
-        angle_atoms = frame.angles.group
-        expected_angle_atoms = [
-            [angle.atom1.idx, angle.atom2.idx, angle.atom3.idx]
-            for angle in structure.angles
-        ]
-        assert np.array_equal(angle_atoms, expected_angle_atoms)
-
-        angle_type_dict = {
-            ("C", "C", "H"): 0,
-            ("H", "C", "C"): 0,
-            ("H", "C", "H"): 1,
-        }
-        expected_angle_typeids = []
-        for angle in structure.angles:
-            expected_angle_typeids.append(
-                angle_type_dict[
-                    (angle.atom1.name, angle.atom2.name, angle.atom3.name)
-                ]
-            )
-        assert np.array_equal(angle_typeids, expected_angle_typeids)
-
-        # Dihedrals
-        n_dihedrals = frame.dihedrals.N
-        assert n_dihedrals == len(structure.rb_torsions)
-
-        expected_unique_dihedral_types = ["opls_140-opls_135-opls_135-opls_140"]
-        dihedral_types = frame.dihedrals.types
-        assert np.array_equal(dihedral_types, expected_unique_dihedral_types)
-
-        dihedral_typeids = frame.dihedrals.typeid
-        dihedral_atoms = frame.dihedrals.group
-        expected_dihedral_atoms = []
-        for dihedral in structure.rb_torsions:
-            expected_dihedral_atoms.append(
-                [
-                    dihedral.atom1.idx,
-                    dihedral.atom2.idx,
-                    dihedral.atom3.idx,
-                    dihedral.atom4.idx,
-                ]
-            )
-        assert np.array_equal(dihedral_atoms, expected_dihedral_atoms)
-        assert np.array_equal(dihedral_typeids, np.zeros(n_dihedrals))
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_pairs(self, benzene):
-        import gsd
-        import gsd.hoomd
-        from foyer import Forcefield
-
-        benzene.save(filename="benzene.gsd", forcefield_name="oplsaa")
-        with gsd.hoomd.open("benzene.gsd", mode="r") as f:
-            frame = f[0]
-
-        structure = benzene.to_parmed()
-        forcefield = Forcefield(name="oplsaa")
-        structure = forcefield.apply(structure)
-
-        # Pairs
-        assert len(frame.pairs.types) == 3
-        assert frame.pairs.N == 21
-
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_units(self, ethane):
-        import gsd
-        import gsd.hoomd
-
-        ref_distance = 3.5
-        ref_energy = 0.066
-        ref_mass = 12.011
-
-        box = Box(lengths=[2.0, 3.0, 4.0], angles=[90.0, 90.0, 90.0])
-        ethane.save(
-            filename="ethane.gsd",
-            forcefield_name="oplsaa",
-            ref_distance=ref_distance,
-            ref_energy=ref_energy,
-            ref_mass=ref_mass,
-            box=box,
-        )
-        with gsd.hoomd.open("ethane.gsd", mode="r") as f:
-            frame = f[0]
-
-        box_from_gsd = frame.configuration.box.astype(float)
-        assert np.array_equal(
-            np.round(box_from_gsd[:3], decimals=5),
-            np.round(np.asarray(box.lengths) * 10 / ref_distance, 5),
-        )
-
-        mass_dict = {"C": 12.011, "H": 1.008}
-        masses = frame.particles.mass.astype(float)
-        for mass, p in zip(masses, ethane.particles()):
-            assert round(mass, 3) == round(mass_dict[p.name] / ref_mass, 3)
-
-        charge_dict = {"C": -0.18, "H": 0.06}
-        charges = frame.particles.charge.astype(float)
-        e0 = 2.396452e-4
-        charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5
-        for charge, particle in zip(charges, ethane.particles()):
-            reduced_charge = charge_dict[particle.name] / charge_factor
-            assert round(charge, 3) == round(reduced_charge, 3)
-
-        positions = frame.particles.position.astype(float)
-        shift = positions[0] - (ethane[0].pos * 10 / ref_distance)
-        shifted_xyz = (ethane.xyz * 10 / ref_distance) + shift
-        assert np.array_equal(
-            np.round(positions, decimals=4), np.round(shifted_xyz, decimals=4)
-        )
-
-    @pytest.mark.skipif(not has_gsd, reason="GSD package not installed")
-    def test_box_dimensions(self, benzene):
-        import gsd
-        import gsd.hoomd
-
-        n_benzenes = 10
-        filled = mb.fill_box(
-            benzene, n_compounds=n_benzenes, box=[0, 0, 0, 4, 4, 4]
-        )
-        filled.save(filename="benzene.gsd")
-        with gsd.hoomd.open("benzene.gsd", mode="r") as f:
-            frame = f[0]
-        positions = frame.particles.position.astype(float)
-        for coords in positions:
-            assert coords.max() < 20
-            assert coords.min() > -20
diff --git a/mbuild/tests/test_hoomd.py b/mbuild/tests/test_hoomd.py
deleted file mode 100644
index 5dc243fc9..000000000
--- a/mbuild/tests/test_hoomd.py
+++ /dev/null
@@ -1,445 +0,0 @@
-import xml.etree.ElementTree
-
-import numpy as np
-import packaging.version
-import pytest
-
-import mbuild as mb
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import get_fn, has_foyer, has_gsd, has_hoomd, import_
-
-if has_hoomd:
-    import hoomd
-
-    if "version" in dir(hoomd):
-        hoomd_version = packaging.version.parse(hoomd.version.version)
-    else:
-        hoomd_version = packaging.version.parse(hoomd.__version__)
-
-
-@pytest.mark.skipif(not has_hoomd, reason="HOOMD is not installed")
-class TestHoomdAny(BaseTest):
-    def test_empty_initial_snapshot(self):
-        import hoomd
-
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        part = mb.Compound(name="Ar")
-        box = mb.Box(lengths=[5, 5, 5], angles=[90, 90, 90])
-        system = mb.fill_box(part, n_compounds=10, box=box)
-
-        if hoomd_version.major == 2:
-            hoomd.context.initialize("")
-            init_snap = hoomd.data.make_snapshot(
-                N=0, box=hoomd.data.boxdim(L=10)
-            )
-        else:
-            init_snap = hoomd.Snapshot()
-            init_snap.configuration.box = hoomd.Box.cube(L=10)
-
-        with pytest.raises(RuntimeError):
-            snap, _ = to_hoomdsnapshot(system, hoomd_snapshot=init_snap)
-
-    def test_compound_from_snapshot(self, ethane):
-        from mbuild.formats.hoomd_snapshot import (
-            from_snapshot,
-            to_hoomdsnapshot,
-        )
-
-        lengths = [5, 5, 5]
-        filled = mb.fill_box(ethane, n_compounds=5, box=mb.Box(lengths))
-        snap, _ = to_hoomdsnapshot(filled)
-        new_filled = from_snapshot(snap, scale=0.1)
-
-        assert filled.n_bonds == new_filled.n_bonds
-        assert filled.n_particles == new_filled.n_particles
-
-        assert np.array_equal(filled.box.angles, new_filled.box.angles)
-        assert np.array_equal(filled.box.lengths, new_filled.box.lengths)
-
-        for i in range(filled.n_particles):
-            assert np.allclose(filled[i].pos, new_filled[i].pos)
-
-    @pytest.mark.skipif(not has_gsd, reason="gsd is not installed")
-    def test_compound_from_gsdsnapshot(self, ethane):
-        import gsd.hoomd
-
-        from mbuild.formats.hoomd_snapshot import (
-            from_snapshot,
-            to_hoomdsnapshot,
-        )
-
-        lengths = [5, 5, 5]
-        filled = mb.fill_box(ethane, n_compounds=5, box=mb.Box(lengths))
-        snap, _ = to_hoomdsnapshot(filled)
-
-        # copy attributes from the snapshot to a gsd snapshot
-        gsd_snap = gsd.hoomd.Frame()
-        gsd_snap.particles.N = snap.particles.N
-        gsd_snap.particles.types = snap.particles.types
-        gsd_snap.particles.typeid = snap.particles.typeid
-        gsd_snap.particles.position = snap.particles.position
-        if hoomd_version.major == 2:
-            gsd_snap.configuration.box = np.array(
-                [
-                    snap.box.Lx,
-                    snap.box.Ly,
-                    snap.box.Lz,
-                    snap.box.xy,
-                    snap.box.xy,
-                    snap.box.yz,
-                ]
-            )
-        else:
-            gsd_snap.configuration.box = snap.configuration.box
-
-        gsd_snap.bonds.N = snap.bonds.N
-        gsd_snap.bonds.group = snap.bonds.group
-        gsd_snap.particles.charge = snap.particles.charge
-        gsd_snap.validate()
-
-        new_filled = from_snapshot(gsd_snap, scale=0.1)
-
-        assert filled.n_bonds == new_filled.n_bonds
-        assert filled.n_particles == new_filled.n_particles
-
-        assert np.array_equal(filled.box.angles, new_filled.box.angles)
-        assert np.array_equal(filled.box.lengths, new_filled.box.lengths)
-
-        for i in range(filled.n_particles):
-            assert np.allclose(filled[i].pos, new_filled[i].pos)
-
-    def test_compound_to_snapshot(self, ethane):
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        snap, _ = to_hoomdsnapshot(ethane)
-
-        assert snap.particles.N == 8
-        assert snap.bonds.N == 7
-        assert snap.angles.N == 0
-
-    def test_particles_to_snapshot(self):
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        part = mb.Compound(name="Ar")
-        box = mb.Box(lengths=[5, 5, 5], angles=[90, 90, 90])
-        system = mb.fill_box(part, n_compounds=10, box=box)
-        snap, _ = to_hoomdsnapshot(system)
-
-        assert snap.particles.N == 10
-        assert snap.bonds.N == 0
-        assert snap.angles.N == 0
-
-    def test_snapshot_from_initial(self):
-        import hoomd
-
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        part = mb.Compound(name="Ar")
-        box = mb.Box(lengths=[5, 5, 5], angles=[90, 90, 90])
-        system = mb.fill_box(part, n_compounds=10, box=box)
-        if hoomd_version.major == 2:
-            hoomd.context.initialize("")
-            init_snap = hoomd.data.make_snapshot(
-                N=10, box=hoomd.data.boxdim(L=10)
-            )
-        else:
-            init_snap = hoomd.Snapshot()
-            init_snap.particles.N = 10
-            init_snap.configuration.box = hoomd.Box.cube(L=10)
-
-        snap, _ = to_hoomdsnapshot(system, hoomd_snapshot=init_snap)
-
-        assert snap.particles.N == 20
-        assert snap.bonds.N == 0
-        assert snap.angles.N == 0
-        if hoomd_version.major == 2:
-            assert (snap.box.Lx, snap.box.Ly, snap.box.Lz) == (50, 50, 50)
-            assert (snap.box.xy, snap.box.xz, snap.box.yz) == (0, 0, 0)
-        else:
-            np.testing.assert_allclose(
-                snap.configuration.box, [50, 50, 50, 0, 0, 0]
-            )
-
-    def test_bad_input_to_snapshot(self):
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        with pytest.raises(ValueError):
-            to_hoomdsnapshot("fake_object")
-
-    def test_non_param_struc_to_snapshot(self, ethane):
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        structure = ethane.to_parmed()
-        snap, _ = to_hoomdsnapshot(structure)
-
-        assert snap.particles.N == 8
-        assert snap.bonds.N == 7
-        assert snap.angles.N == 0
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_param_structure_to_snapshot(self, ethane):
-        from foyer.forcefield import Forcefield
-
-        from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot
-
-        ff = Forcefield(name="oplsaa")
-        structure = ff.apply(ethane)
-        snap, _ = to_hoomdsnapshot(structure)
-
-        assert snap.particles.N == 8
-        assert snap.bonds.N == 7
-        assert snap.angles.N == 12
-        assert snap.dihedrals.N == 9
-        assert snap.pairs.N == 9
-
-
-@pytest.mark.skipif(
-    not has_hoomd or hoomd_version.major != 2,
-    reason="HOOMD is not installed or is the wrong version",
-)
-class TestHoomdSimulation(BaseTest):
-    def test_bad_input_to_hoomdsimulation(self):
-        from mbuild.formats.hoomd_simulation import create_hoomd_simulation
-
-        with pytest.raises(ValueError):
-            create_hoomd_simulation("fake_object", 2)
-
-    def test_compound_to_hoomdsimulation(self, ethane):
-        from mbuild.formats.hoomd_simulation import create_hoomd_simulation
-
-        with pytest.raises(ValueError):
-            create_hoomd_simulation(ethane, 2.5)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_structure_to_hoomdsimulation(self, ethane):
-        import hoomd
-        from foyer.forcefield import Forcefield
-
-        from mbuild.formats.hoomd_simulation import create_hoomd_simulation
-
-        ff = Forcefield(name="oplsaa")
-        structure = ff.apply(ethane)
-        sim = hoomd.context.SimulationContext()
-        with sim:
-            create_hoomd_simulation(structure, 2.5)
-
-            sim_forces = hoomd.context.current.forces
-            pair_force = import_("hoomd.md.pair")
-            charge_force = import_("hoomd.md.charge")
-            special_pair_force = import_("hoomd.md.special_pair")
-            bond_force = import_("hoomd.md.bond")
-            angle_force = import_("hoomd.md.angle")
-            dihedral_force = import_("hoomd.md.dihedral")
-
-            assert isinstance(sim_forces[0], pair_force.lj)
-            assert isinstance(sim_forces[1], charge_force.pppm)
-            assert isinstance(sim_forces[2], pair_force.ewald)
-            assert isinstance(sim_forces[3], special_pair_force.lj)
-            assert isinstance(sim_forces[4], special_pair_force.coulomb)
-            assert isinstance(sim_forces[5], bond_force.harmonic)
-            assert isinstance(sim_forces[6], angle_force.harmonic)
-            assert isinstance(sim_forces[7], dihedral_force.opls)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_lj_to_hoomdsimulation(self):
-        import hoomd
-        from foyer.forcefield import Forcefield
-
-        from mbuild.formats.hoomd_simulation import create_hoomd_simulation
-
-        box = mb.Compound()
-        box.add(mb.Compound(name="Ar", pos=[1, 1, 1]))
-        box.add(mb.Compound(name="Ar", pos=[1, 1, 1]))
-        ff = Forcefield(forcefield_files=get_fn("lj.xml"))
-        structure = ff.apply(box)
-        structure.box = [10, 10, 10, 90, 90, 90]
-        sim = hoomd.context.SimulationContext()
-        with sim:
-            create_hoomd_simulation(structure, 2.5)
-            sim_forces = hoomd.context.current.forces
-            pair_force = import_("hoomd.md.pair")
-
-            assert isinstance(sim_forces[0], pair_force.lj)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    @pytest.mark.skipif(not has_gsd, reason="GSD is not installed")
-    def test_hoomdsimulation_restart(self):
-        import gsd.hoomd
-        import hoomd
-        from foyer.forcefield import Forcefield
-
-        from mbuild.formats.hoomd_simulation import create_hoomd_simulation
-
-        box = mb.Compound()
-        box.add(mb.Compound(name="Ar", pos=[1, 1, 1]))
-        box.add(mb.Compound(name="Ar", pos=[1, 1, 1]))
-        ff = Forcefield(forcefield_files=get_fn("lj.xml"))
-        structure = ff.apply(box)
-        structure.box = [10, 10, 10, 90, 90, 90]
-        sim = hoomd.context.SimulationContext()
-        with sim:
-            hoomd_obj, ref_vals = create_hoomd_simulation(
-                structure, 2.5, restart=get_fn("restart.gsd")
-            )
-            sim_forces = hoomd.context.current.forces
-            pair_force = import_("hoomd.md.pair")
-
-            assert isinstance(sim_forces[0], pair_force.lj)
-
-        snap = hoomd_obj[0]
-        with gsd.hoomd.open(get_fn("restart.gsd")) as f:
-            rsnap = f[0]
-        assert np.array_equal(snap.particles.position, rsnap.particles.position)
-
-    def test_hoomdsimulation_nlist(self, ethane):
-        hoomd_simulation = import_("mbuild.formats.hoomd_simulation")
-        hoomd = import_("hoomd")
-        hoomd.md = import_("hoomd.md")
-
-        with pytest.raises(ValueError):
-            hoomd_simulation.create_hoomd_simulation(
-                ethane, 2.5, nlist=hoomd.md.nlist.tree
-            )
-
-
-@pytest.mark.skipif(
-    not has_hoomd or hoomd_version.major < 3,
-    reason="HOOMD is not installed or is the wrong version",
-)
-class TestHoomdForcefield(BaseTest):
-    def test_bad_input_to_hoomd_forcefield(self):
-        from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield
-
-        with pytest.raises(ValueError):
-            create_hoomd_forcefield("fake_object", 2.5)
-
-    def test_compound_to_hoomd_forcefield(self, ethane):
-        from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield
-
-        with pytest.raises(ValueError):
-            create_hoomd_forcefield(ethane, 2.5)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_structure_to_hoomd_forcefield(self, ethane):
-        import hoomd
-        from foyer.forcefield import Forcefield
-
-        from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield
-
-        ff = Forcefield(name="oplsaa")
-        structure = ff.apply(ethane)
-
-        snapshot, forces, ref_values = create_hoomd_forcefield(structure, 2.5)
-
-        assert isinstance(forces[0], hoomd.md.pair.LJ)
-        assert isinstance(forces[1], hoomd.md.pair.Ewald)
-        assert isinstance(forces[2], hoomd.md.long_range.pppm.Coulomb)
-        assert isinstance(forces[3], hoomd.md.special_pair.LJ)
-        assert isinstance(forces[4], hoomd.md.special_pair.Coulomb)
-        assert isinstance(forces[5], hoomd.md.bond.Harmonic)
-        assert isinstance(forces[6], hoomd.md.angle.Harmonic)
-        assert isinstance(forces[7], hoomd.md.dihedral.OPLS)
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed")
-    def test_lj_to_hoomd_forcefield(self):
-        import hoomd
-        from foyer.forcefield import Forcefield
-
-        from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield
-
-        box = mb.Compound()
-        box.add(mb.Compound(name="Ar", pos=[1, 1, 1]))
-        box.add(mb.Compound(name="Ar", pos=[1, 1, 1]))
-        ff = Forcefield(forcefield_files=get_fn("lj.xml"))
-        structure = ff.apply(box)
-        structure.box = [10, 10, 10, 90, 90, 90]
-
-        snapshot, forces, ref_values = create_hoomd_forcefield(structure, 2.5)
-
-        assert isinstance(forces[0], hoomd.md.pair.LJ)
-
-
-class TestHoomdXML(BaseTest):
-    def test_save(self, ethane):
-        ethane.save(filename="ethane.hoomdxml")
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_save_forcefield(self, ethane):
-        ethane.save(filename="ethane-opls.hoomdxml", forcefield_name="oplsaa")
-
-    def test_save_box(self, ethane):
-        box = mb.box.Box(lengths=[2.0, 2.0, 2.0], angles=[90, 90, 90])
-        ethane.save(filename="ethane-box.hoomdxml", box=box)
-
-    def test_save_triclinic_box_(self, ethane):
-        box = mb.Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[60, 70, 80])
-        ethane.save(filename="triclinic-box.hoomdxml", box=box)
-
-    def test_rigid(self, benzene):
-        n_benzenes = 10
-        benzene.name = "Benzene"
-        filled = mb.fill_box(
-            benzene, n_compounds=n_benzenes, box=[0, 0, 0, 4, 4, 4]
-        )
-        filled.label_rigid_bodies(
-            discrete_bodies="Benzene", rigid_particles="C"
-        )
-        filled.save(filename="benzene.hoomdxml")
-
-        xml_file = xml.etree.ElementTree.parse("benzene.hoomdxml").getroot()
-        body_text = xml_file[0].find("body").text
-        rigid_bodies = [int(body) for body in body_text.split("\n") if body]
-        for body_id in range(10):
-            assert rigid_bodies.count(body_id) == 6
-        assert rigid_bodies.count(-1) == n_benzenes * 6
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_number_in_each_section(self, box_of_benzenes):
-        box_of_benzenes.save(
-            filename="benzene.hoomdxml", forcefield_name="oplsaa"
-        )
-        xml_file = xml.etree.ElementTree.parse("benzene.hoomdxml").getroot()
-        for attribute in ["position", "type", "mass", "charge"]:
-            body_text = xml_file[0].find(attribute).text
-            list_of_things = [x for x in body_text.split("\n") if x]
-            assert len(list_of_things) == 12 * 10
-        for attribute, number in [
-            ("bond", 12),
-            ("angle", 18),
-            ("dihedral", 24),
-        ]:
-            body_text = xml_file[0].find(attribute).text
-            list_of_things = [x for x in body_text.split("\n") if x]
-            assert len(list_of_things) == number * 10
-
-    def test_box_dimensions(self, benzene):
-        n_benzenes = 10
-        filled = mb.fill_box(
-            benzene, n_compounds=n_benzenes, box=[0, 0, 0, 4, 4, 4]
-        )
-        filled.save(filename="benzene.hoomdxml")
-        for atom in mb.load("benzene.hoomdxml"):
-            assert atom.pos.max() < 20
-            assert atom.pos.min() > -20
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_auto_scale_forcefield(self, ethane):
-        ethane.save(
-            filename="ethane-opls.hoomdxml",
-            forcefield_name="oplsaa",
-            auto_scale=True,
-        )
-        xml_file = xml.etree.ElementTree.parse("ethane-opls.hoomdxml").getroot()
-        masses = xml_file[0].find("mass").text.splitlines()
-        # We use 1 and 5 since the first element of masses is empty
-        assert masses[1] == "1.0"
-        assert masses[5] == "1.0"
-        pair_coeffs = [
-            _.split("\t")
-            for _ in xml_file[0].find("pair_coeffs").text.splitlines()
-        ]
-        # The first element is empty, the next element should be
-        # ['opls_135', '1.0000', '1.0000']
-        assert pair_coeffs[1][1] == "1.0000"
-        assert pair_coeffs[1][2] == "1.0000"
diff --git a/mbuild/tests/test_lammpsdata.py b/mbuild/tests/test_lammpsdata.py
deleted file mode 100644
index 3b4711dcc..000000000
--- a/mbuild/tests/test_lammpsdata.py
+++ /dev/null
@@ -1,742 +0,0 @@
-from pathlib import Path
-
-import numpy as np
-import pytest
-from ele.element import element_from_symbol
-from pytest import FixtureRequest
-from scipy.constants import epsilon_0
-
-import mbuild as mb
-from mbuild.formats.lammpsdata import write_lammpsdata
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import get_fn, has_foyer
-
-KCAL_TO_KJ = 4.184
-ANG_TO_NM = 0.10
-ELEM_TO_COUL = 1.602176634e-19
-
-
-@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-class TestLammpsData(BaseTest):
-    def test_save(self, ethane):
-        ethane.save(filename="ethane.lammps")
-
-    @pytest.fixture(scope="session")
-    def lj_save(self, tmpdir_factory, sigma=None, epsilon=None, mass=None):
-        def _create_lammps(
-            ethane, tmp, sigma=sigma, epsilon=epsilon, mass=mass
-        ):
-            from foyer import Forcefield
-
-            OPLSAA = Forcefield(name="oplsaa")
-            structure = OPLSAA.apply(ethane)
-            fn = tmpdir_factory.mktemp("data").join("lj.lammps")
-            write_lammpsdata(
-                filename=str(fn),
-                structure=structure,
-                unit_style="lj",
-                sigma_conversion_factor=sigma,
-                epsilon_conversion_factor=epsilon,
-                mass_conversion_factor=mass,
-            )
-            return str(fn)
-
-        return _create_lammps
-
-    @pytest.fixture(scope="session")
-    def real_save(self, tmpdir_factory):
-        def _create_lammps(ethane, tmp):
-            from foyer import Forcefield
-
-            OPLSAA = Forcefield(name="oplsaa")
-            structure = OPLSAA.apply(ethane)
-            fn = tmpdir_factory.mktemp("data").join("lj.lammps")
-            write_lammpsdata(
-                filename=str(fn), structure=structure, unit_style="real"
-            )
-            return str(fn)
-
-        return _create_lammps
-
-    @pytest.mark.parametrize("unit_style", ["real", "lj"])
-    def test_save_forcefield(self, ethane, unit_style):
-        ethane.save(
-            filename="ethane-opls.lammps",
-            forcefield_name="oplsaa",
-            unit_style=unit_style,
-        )
-
-    @pytest.mark.parametrize("pair_coeff_label", [None, "lj", "lj/coul/cut"])
-    def test_save_pair_coeff_label(self, ethane, pair_coeff_label):
-        ethane.save(
-            filename="ethane-opls.lammps",
-            forcefield_name="oplsaa",
-            unit_style="real",
-            pair_coeff_label=pair_coeff_label,
-        )
-        with open("ethane-opls.lammps", "r") as fp:
-            for line in fp:
-                if "Pair Coeffs" in line:
-                    if pair_coeff_label is None:
-                        assert "# lj" in line
-                    else:
-                        assert "# {}".format(pair_coeff_label) in line
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_save_charmm(self):
-        from foyer import Forcefield
-
-        cmpd = mb.load(get_fn("charmm_dihedral.mol2"))
-        for i in cmpd.particles():
-            i.name = "_{}".format(i.name)
-        structure = cmpd.to_parmed(
-            box=cmpd.get_boundingbox(),
-            residues=set([p.parent.name for p in cmpd.particles()]),
-        )
-        ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")])
-        structure = ff.apply(structure, assert_dihedral_params=False)
-        write_lammpsdata(structure, "charmm_dihedral.lammps")
-        out_lammps = open("charmm_dihedral.lammps", "r").readlines()
-        found_angles = False
-        found_dihedrals = False
-        for i, line in enumerate(out_lammps):
-            if "Angle Coeffs" in line:
-                assert "# charmm" in line
-                assert (
-                    "#\tk(kcal/mol/rad^2)\t\ttheteq(deg)\tk(kcal/mol/angstrom^2)\treq(angstrom)\n"
-                    in out_lammps[i + 1]
-                )
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 5
-                found_angles = True
-            elif "Dihedral Coeffs" in line:
-                assert "# charmm" in line
-                assert "#k, n, phi, weight" in out_lammps[i + 1]
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 5
-                found_dihedrals = True
-            else:
-                pass
-        assert found_angles
-        assert found_dihedrals
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_singleterm_charmm(self):
-        from foyer import Forcefield
-
-        from mbuild.formats.lammpsdata import write_lammpsdata
-
-        cmpd = mb.load(get_fn("charmm_dihedral.mol2"))
-        for i in cmpd.particles():
-            i.name = "_{}".format(i.name)
-        structure = cmpd.to_parmed(
-            box=cmpd.get_boundingbox(),
-            residues=set([p.parent.name for p in cmpd.particles()]),
-        )
-        ff = Forcefield(
-            forcefield_files=[get_fn("charmm_truncated_singleterm.xml")]
-        )
-        structure = ff.apply(structure, assert_dihedral_params=False)
-        write_lammpsdata(structure, "charmm_dihedral_singleterm.lammps")
-        out_lammps = open("charmm_dihedral_singleterm.lammps", "r").readlines()
-        found_dihedrals = False
-        for i, line in enumerate(out_lammps):
-            if "Dihedral Coeffs" in line:
-                assert "# charmm" in line
-                assert "#k, n, phi, weight" in out_lammps[i + 1]
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 5
-                assert float(
-                    out_lammps[i + 2].split("#")[0].split()[4]
-                ) == float("1.0")
-                found_dihedrals = True
-            else:
-                pass
-        assert found_dihedrals
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_charmm_improper(self):
-        from foyer import Forcefield
-
-        import mbuild as mb
-        from mbuild.formats.lammpsdata import write_lammpsdata
-
-        system = mb.Compound()
-        first = mb.Particle(name="_CTL2", pos=[-1, 0, 0])
-        second = mb.Particle(name="_CL", pos=[0, 0, 0])
-        third = mb.Particle(name="_OBL", pos=[1, 0, 0])
-        fourth = mb.Particle(name="_OHL", pos=[0, 1, 0])
-
-        system.add([first, second, third, fourth])
-
-        system.add_bond((first, second))
-        system.add_bond((second, third))
-        system.add_bond((second, fourth))
-
-        ff = Forcefield(forcefield_files=[get_fn("charmm36_cooh.xml")])
-        struc = ff.apply(
-            system,
-            assert_angle_params=False,
-            assert_dihedral_params=False,
-            assert_improper_params=False,
-        )
-
-        write_lammpsdata(struc, "charmm_improper.lammps")
-        out_lammps = open("charmm_improper.lammps", "r").readlines()
-        found_impropers = False
-        for i, line in enumerate(out_lammps):
-            if "Improper Coeffs" in line:
-                assert "# harmonic" in line
-                assert "k, phi" in out_lammps[i + 1]
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 3
-                assert out_lammps[i + 2].split("#")[0].split()[0] == "1"
-                found_impropers = True
-        assert found_impropers
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_amber(self):
-        from foyer import Forcefield
-
-        from mbuild.formats.lammpsdata import write_lammpsdata
-
-        cmpd = mb.load("C1(=CC=CC=C1)F", smiles=True)
-
-        ff = Forcefield(forcefield_files=[get_fn("gaff_test.xml")])
-        structure = ff.apply(cmpd)
-
-        write_lammpsdata(
-            structure, "amber.lammps", zero_dihedral_weighting_factor=True
-        )
-        out_lammps = open("amber.lammps", "r").readlines()
-        found_angles = False
-        found_dihedrals = False
-        found_impropers = False
-        for i, line in enumerate(out_lammps):
-            if "Angle Coeffs" in line:
-                assert "# harmonic" in line
-                assert (
-                    "#\tk(kcal/mol/rad^2)\t\ttheteq(deg)" in out_lammps[i + 1]
-                )
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 3
-                assert out_lammps[i + 2].split("#")[0].split()[0] == "1"
-                found_angles = True
-            elif "Dihedral Coeffs" in line:
-                assert "# charmm" in line
-                assert "#k, n, phi, weight" in out_lammps[i + 1]
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 5
-                assert float(out_lammps[i + 2].split("#")[0].split()[4]) == 0.0
-                assert out_lammps[i + 2].split("#")[0].split()[0] == "1"
-                found_dihedrals = True
-            elif "Improper Coeffs" in line:
-                assert "# cvff" in line
-                assert "#K, d, n" in out_lammps[i + 1]
-                assert len(out_lammps[i + 2].split("#")[0].split()) == 4
-                assert out_lammps[i + 2].split("#")[0].split()[2] == "-1"
-                assert out_lammps[i + 2].split("#")[0].split()[0] == "1"
-                found_impropers = True
-            else:
-                pass
-        assert found_angles
-        assert found_dihedrals
-        assert found_impropers
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_save_forcefield_with_same_struct(self):
-        from foyer import Forcefield
-
-        from mbuild.formats.lammpsdata import write_lammpsdata
-
-        system = mb.load("C1(=CC=CC=C1)F", smiles=True)
-
-        ff = Forcefield(forcefield_files=[get_fn("gaff_test.xml")])
-        struc = ff.apply(
-            system,
-            assert_angle_params=False,
-            assert_dihedral_params=False,
-            assert_improper_params=False,
-        )
-        write_lammpsdata(
-            struc, "charmm_improper.lammps", zero_dihedral_weighting_factor=True
-        )
-        for i in range(3):
-            xyz = struc.coordinates
-            xyz = xyz + np.array([1, 1, 1])
-            struc.coordinates = xyz
-            write_lammpsdata(
-                struc,
-                f"charmm_improper{i}.lammps",
-                zero_dihedral_weighting_factor=True,
-            )
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    @pytest.mark.parametrize("unit_style", ["real", "lj"])
-    def test_save_box(self, ethane, unit_style):
-        box = mb.Box(
-            lengths=np.array([2.0, 2.0, 2.0]), angles=[90.0, 90.0, 90.0]
-        )
-        ethane.save(
-            filename="ethane-box.lammps",
-            forcefield_name="oplsaa",
-            box=box,
-            unit_style=unit_style,
-        )
-
-    def test_nbfix(self, ethane):
-        from foyer import Forcefield
-
-        OPLSAA = Forcefield(name="oplsaa")
-        structure = OPLSAA.apply(ethane)
-        structure.combining_rule = "lorentz"
-        # Add nbfixes
-        types = list(set([a.atom_type for a in structure.atoms]))
-        types[0].add_nbfix(types[1].name, 1.2, 2.1)
-        types[1].add_nbfix(types[0].name, 1.2, 2.1)
-        write_lammpsdata(filename="nbfix.lammps", structure=structure)
-
-        checked_section = False
-        with open("nbfix.lammps", "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "PairIJ Coeffs" in line:
-                    fi.readline()
-                    line = fi.readline().partition("#")[0]
-                    assert np.allclose(
-                        np.asarray(line.split(), dtype=float),
-                        [1, 1, 0.066, 3.5],
-                    )
-                    line = fi.readline().partition("#")[0]
-                    assert np.allclose(
-                        np.asarray(line.split(), dtype=float),
-                        [1, 2, 2.1, 1.06907846],
-                    )
-                    line = fi.readline().partition("#")[0]
-                    assert np.allclose(
-                        np.asarray(line.split(), dtype=float), [2, 2, 0.03, 2.5]
-                    )
-                    line = fi.readline()
-                    checked_section = True
-                # Break if PairIJ Coeffs is not found
-                if "Atoms" in line:
-                    break
-        structure.combining_rule = "geometric"
-        write_lammpsdata(filename="nbfix.lammps", structure=structure)
-        write_lammpsdata(
-            filename="nbfix.lammps",
-            structure=structure,
-            nbfix_in_data_file=False,
-        )
-        with pytest.raises(ValueError) as exc_info:
-            structure.combining_rule = "error"
-            write_lammpsdata(filename="nbfix.lammps", structure=structure)
-        assert str(exc_info.value).startswith("combining_rule must be")
-
-    def test_save_triclinic_box(self, ethane):
-        box = mb.Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[60, 70, 80])
-        ethane.save(
-            filename="triclinic-box.lammps", forcefield_name="oplsaa", box=box
-        )
-
-    def test_save_orthorhombic_box(self, ethane):
-        box = mb.Box(lengths=np.array([2.0, 2.0, 2.0]))
-        ethane.save(
-            filename="ortho-box.lammps", forcefield_name="oplsaa", box=box
-        )
-        with open("ortho-box.lammps", "r") as fi:
-            for line in fi:
-                assert "xy" not in line
-                assert "xz" not in line
-                assert "yz" not in line
-
-    @pytest.mark.parametrize(
-        "atom_style, n_columns",
-        [("full", 7), ("atomic", 5), ("molecular", 6), ("charge", 6)],
-    )
-    def test_writing_atom_styles(self, ethane, atom_style, n_columns):
-        ethane.save(filename="ethane.lammps", atom_style=atom_style)
-        with open("ethane.lammps", "r") as f:
-            for line in f:
-                if "Atoms" not in line:
-                    continue
-                    atoms_header = next(f)
-                    first_atom_line = next(f)
-                    columns = first_atom_line.split("\t")
-                    assert len(columns) == n_columns
-                else:
-                    assert "# " + atom_style in line
-
-    @pytest.mark.parametrize(
-        "offset, expected_value", [(0, ("0", "1")), (1, ("1", "2"))]
-    )
-    def test_resid(self, offset, expected_value, ethane, methane):
-        structure = ethane.to_parmed() + methane.to_parmed()
-        n_atoms = len(structure.atoms)
-        write_lammpsdata(structure, "compound.lammps", moleculeID_offset=offset)
-        res_list = list()
-        with open("compound.lammps", "r") as f:
-            for i, line in enumerate(f):
-                if "Atoms" in line:
-                    break
-        atom_lines = open("compound.lammps", "r").readlines()[
-            i + 2 : i + n_atoms + 2
-        ]
-        for line in atom_lines:
-            res_list.append(line.rstrip().split()[1])
-        assert set(res_list) == set(expected_value)
-
-    def test_box_bounds(self, ethane):
-        from foyer import Forcefield
-
-        OPLSAA = Forcefield(name="oplsaa")
-        structure = OPLSAA.apply(ethane)
-        box = mb.Box.from_mins_maxs_angles(
-            mins=np.array([-1.0, -2.0, -3.0]),
-            maxs=np.array([3.0, 2.0, 1.0]),
-            angles=[90.0, 90.0, 90.0],
-        )
-        # box = ethane.get_boundingbox()
-
-        write_lammpsdata(
-            filename="box.lammps",
-            structure=structure,
-            unit_style="real",
-            mins=[0.0, 0.0, 0.0],
-            maxs=[m for m in box.lengths],
-        )
-
-        checked_section = False
-        # NOTE, need to figure out how to handle lo and hi of a compound
-        with open("box.lammps", "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "xlo" in line:
-                    xlo = float(line.split()[0])
-                    xhi = float(line.split()[1])
-                    assert np.isclose(xlo, 0.0)
-                    assert np.isclose(xhi, 40.0)
-
-                    line = fi.readline()
-                    ylo = float(line.split()[0])
-                    yhi = float(line.split()[1])
-                    assert np.isclose(ylo, 0.0)
-                    assert np.isclose(yhi, 40.0)
-
-                    line = fi.readline()
-                    zlo = float(line.split()[0])
-                    zhi = float(line.split()[1])
-                    assert np.isclose(zlo, 0.0)
-                    assert np.isclose(zhi, 40.0)
-
-                    checked_section = True
-
-        write_lammpsdata(
-            filename="box.lammps", structure=structure, unit_style="real"
-        )
-
-        checked_section = False
-        with open("box.lammps", "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "xlo" in line:
-                    xlo = float(line.split()[0])
-                    xhi = float(line.split()[1])
-                    assert np.isclose(xlo, 0.0)
-                    assert np.isclose(xhi, 7.13999987)
-
-                    line = fi.readline()
-                    ylo = float(line.split()[0])
-                    yhi = float(line.split()[1])
-                    assert np.isclose(ylo, 0.0)
-                    assert np.isclose(yhi, 7.93800011)
-
-                    line = fi.readline()
-                    zlo = float(line.split()[0])
-                    zhi = float(line.split()[1])
-                    assert np.isclose(zlo, 0.0)
-                    assert np.isclose(zhi, 6.646)
-
-                    checked_section = True
-
-    def test_lj_box(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_sigma = 0.35  # nm
-        box_length = np.array(ethane.get_boundingbox().lengths) + [
-            0.5,
-            0.5,
-            0.5,
-        ]  # 0.5nm buffer around box
-        box_length /= ethane_sigma  # Unitless
-        with open(str(fn), "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "dihedral types" in line:  # line before box info
-                    fi.readline()
-                    for i in range(len(box_length)):
-                        line = float(fi.readline().split()[1])
-                        assert np.isclose(float(line), box_length[i])
-                    checked_section = True
-
-    def test_lj_masses(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Masses" in line:
-                    fi.readline()
-                    line = float(fi.readline().split()[1])
-                    assert np.isclose(line, 1.00)
-                    checked_section = True
-
-    def test_real_masses(self, ethane, real_save):
-        fn = real_save(ethane, Path.cwd())
-        checked_section = False
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Masses" in line:
-                    fi.readline()
-                    line = float(fi.readline().split()[1])
-                    assert np.isclose(line, 12.010780)
-                    checked_section = True
-
-    def test_lj_pairs(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Pair Coeffs" in line:
-                    fi.readline()
-                    line = fi.readline().split()
-                    epsilon = float(line[1])
-                    sigma = float(line[2])
-                    assert np.isclose(epsilon, 1.00)
-                    assert np.isclose(sigma, 1.00)
-                    checked_section = True
-
-    def test_lj_passed_pairs(self, ethane, lj_save):
-        fn = lj_save(
-            ethane,
-            Path.cwd(),
-            sigma=1,
-            epsilon=1,
-        )
-        checked_section = False
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Pair Coeffs" in line:
-                    fi.readline()
-                    line = fi.readline().split()
-                    epsilon = float(line[1])
-                    sigma = float(line[2])
-                    assert np.isclose(epsilon, 0.276144, atol=1e-5)
-                    assert np.isclose(sigma, 0.35)
-                    checked_section = True
-
-        with pytest.raises(ValueError) as exc_info:
-            lj_save(ethane, Path.cwd(), sigma=0)
-        exception_str = "The sigma conversion factor to convert to LJ"
-        assert str(exc_info.value).startswith(exception_str)
-
-    def test_real_pairs(self, ethane, real_save):
-        fn = real_save(ethane, Path.cwd())
-        checked_section = False
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Pair Coeffs" in line:
-                    fi.readline()
-                    line = fi.readline().split()
-                    epsilon = float(line[1])
-                    sigma = float(line[2])
-                    assert np.isclose(epsilon, 0.066)
-                    assert np.isclose(sigma, 3.5)
-                    checked_section = True
-
-    def test_lj_bonds(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_bondsk = (
-            np.array([224262.4, 284512.0]) * (0.35) ** 2 / (0.276144) / 2
-        )  # kj/mol/nm**2 to lammps
-        ethane_bondsreq = np.array([0.1529, 0.109]) * 0.35**-1  # unitless
-        ethane_lj_bonds = zip(ethane_bondsk, ethane_bondsreq)
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Bond Coeffs" in line:
-                    fi.readline()
-                    for i, bond_params in enumerate(ethane_lj_bonds):
-                        line = fi.readline()
-                        assert np.allclose(
-                            (float(line.split()[1]), float(line.split()[2])),
-                            bond_params,
-                            atol=1e-3,
-                        )
-                    assert int(line.split()[0]) == i + 1
-                    checked_section = True
-
-    def test_real_bonds(self, ethane, real_save):
-        fn = real_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_bondsk = (
-            np.array([224262.4, 284512.0]) / KCAL_TO_KJ * ANG_TO_NM**2 / 2
-        )  # kj/mol/nm**2 to lammps
-        ethane_bondsreq = np.array([0.1529, 0.109]) / ANG_TO_NM
-        ethane_real_bonds = zip(ethane_bondsk, ethane_bondsreq)
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Bond Coeffs" in line:
-                    fi.readline()
-                    for bond_params in ethane_real_bonds:
-                        line = fi.readline()
-                        assert np.allclose(
-                            (float(line.split()[1]), float(line.split()[2])),
-                            bond_params,
-                            atol=1e-3,
-                        )
-                    checked_section = True
-
-    def test_lj_angles(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_anglesk = (
-            np.array([313.8, 276.144]) / 0.276144 / 2
-        )  # kj/mol/nm**2 to lammps
-        ethane_anglesreq = (
-            np.array([1.93207948196, 1.88146493365]) * 180 / np.pi
-        )  # radians to lammps
-        ethane_lj_angles = zip(ethane_anglesk, ethane_anglesreq)
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Angle Coeffs" in line:
-                    fi.readline()
-                    for angle_params in ethane_lj_angles:
-                        line = fi.readline()
-                        assert np.allclose(
-                            (float(line.split()[1]), float(line.split()[2])),
-                            angle_params,
-                        )
-                    checked_section = True
-
-    def test_real_angles(self, ethane, real_save):
-        fn = real_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_anglesk = (
-            np.array([313.8, 276.144]) / 2 / KCAL_TO_KJ
-        )  # kj/mol to lammps kcal/mol
-        ethane_anglesreq = (
-            np.array([1.93207948196, 1.88146493365]) * 180 / np.pi
-        )  # radians to lammps
-        ethane_real_angles = zip(ethane_anglesk, ethane_anglesreq)
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Angle Coeffs" in line:
-                    fi.readline()
-                    for i, angle_params in enumerate(ethane_real_angles):
-                        line = fi.readline()
-                        assert np.allclose(
-                            (float(line.split()[1]), float(line.split()[2])),
-                            angle_params,
-                        )
-                        assert int(line.split()[0]) == i + 1
-                    checked_section = True
-
-    def test_lj_dihedrals(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_diheds = (
-            np.array([-0.00000, -0.00000, 0.30000, -0.00000]) / 0.066
-        )  # kcal/mol to lammps
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Dihedral Coeffs" in line:
-                    fi.readline()
-                    line = fi.readline()
-                    assert np.allclose(
-                        (
-                            float(line.split()[1]),
-                            float(line.split()[2]),
-                            float(line.split()[3]),
-                            float(line.split()[4]),
-                        ),
-                        ethane_diheds,
-                    )
-                    checked_section = True
-
-    def test_real_dihedrals(self, ethane, real_save):
-        fn = real_save(ethane, Path.cwd())
-        checked_section = False
-        ethane_diheds = np.array([-0.00000, -0.00000, 0.30000, -0.00000])
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Dihedral Coeffs" in line:
-                    fi.readline()
-                    line = fi.readline()
-                    assert np.allclose(
-                        (
-                            float(line.split()[1]),
-                            float(line.split()[2]),
-                            float(line.split()[3]),
-                            float(line.split()[4]),
-                        ),
-                        ethane_diheds,
-                    )
-                    assert int(line.split()[0]) == 1
-                    checked_section = True
-
-    def test_lj_charges(self, ethane, lj_save):
-        fn = lj_save(ethane, Path.cwd())
-        checked_section = False
-        charge_dict = {"C": -0.18, "H": 0.06}
-        ethane_charges = np.array(
-            [charge_dict[part.name] for part in ethane.particles()]
-        )
-        # in coulombs, convert to lj units
-        ethane_charges /= (
-            np.sqrt(4 * np.pi * 0.276144 * 0.35 * epsilon_0 * 10**-6)
-            / ELEM_TO_COUL
-        )
-        print(ethane_charges)
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Atoms " in line:
-                    fi.readline()
-                    assert np.allclose(
-                        float(fi.readline().split()[3]),
-                        ethane_charges[0],
-                        atol=1e-15,
-                    )
-                    assert np.allclose(
-                        float(fi.readline().split()[3]),
-                        ethane_charges[1],
-                        atol=1e-15,
-                    )
-                    checked_section = True
-
-    def test_real_charges(self, ethane, real_save):
-        fn = real_save(ethane, Path.cwd())
-        checked_section = False
-        charge_dict = {"C": -0.18, "H": 0.06}
-        ethane_charges = [charge_dict[part.name] for part in ethane.particles()]
-        print(ethane_charges)
-        # in coulombs
-        with open(fn, "r") as fi:
-            while not checked_section:
-                line = fi.readline()
-                if "Atoms " in line:
-                    fi.readline()
-                    assert np.allclose(
-                        float(fi.readline().split()[3]),
-                        ethane_charges[0],
-                    )
-                    assert np.allclose(
-                        float(fi.readline().split()[3]),
-                        ethane_charges[1],
-                    )
-                    checked_section = True
diff --git a/mbuild/tests/test_par.py b/mbuild/tests/test_par.py
deleted file mode 100644
index 732deb867..000000000
--- a/mbuild/tests/test_par.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import numpy as np
-import pytest
-
-import mbuild as mb
-from mbuild.formats.par_writer import write_par
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import get_fn, has_foyer
-
-
-@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-class TestPar(BaseTest):
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_save_charmm(self):
-        cmpd = mb.load(get_fn("charmm_dihedral.mol2"))
-        for i in cmpd.particles():
-            i.name = "_{}".format(i.name)
-        structure = cmpd.to_parmed(
-            box=cmpd.get_boundingbox(),
-            residues=set([p.parent.name for p in cmpd.particles()]),
-        )
-
-        from foyer import Forcefield
-
-        ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")])
-        structure = ff.apply(structure, assert_dihedral_params=False)
-
-        write_par(structure, "charmm_dihedral.par")
-
-    @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
-    def test_save_forcefield(self, ethane):
-        ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa")
-
-    def test_par_parameters(self, ethane):
-        ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa")
-        from parmed.charmm import CharmmParameterSet
-
-        pset = CharmmParameterSet.load_set(pfile="ethane-opls.par")
-        assert len(pset.bond_types) == 3
-        assert len(pset.angle_types) == 3
-        assert len(pset.atom_types) == 2
diff --git a/mbuild/tests/test_protobuf.py b/mbuild/tests/test_protobuf.py
deleted file mode 100644
index 001b8df55..000000000
--- a/mbuild/tests/test_protobuf.py
+++ /dev/null
@@ -1,18 +0,0 @@
-import pytest
-
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import has_protobuf
-
-
-class TestPB2(BaseTest):
-    @pytest.mark.skipif(
-        not has_protobuf, reason="Protobuf package not installed"
-    )
-    def test_loop(self, ethane):
-        from mbuild.formats.protobuf import read_pb2, write_pb2
-
-        write_pb2(ethane, "ethane.pb2")
-        proto = read_pb2("ethane.pb2")
-        assert ethane.n_particles == proto.n_particles
-        assert ethane.n_bonds == proto.n_bonds
-        assert len(ethane.children) == len(proto.children)
diff --git a/mbuild/tests/test_rigid.py b/mbuild/tests/test_rigid.py
index 41b61baec..d3739bd0a 100644
--- a/mbuild/tests/test_rigid.py
+++ b/mbuild/tests/test_rigid.py
@@ -3,7 +3,7 @@
 
 import mbuild as mb
 from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import get_fn
+from mbuild.utils.io import has_hoomd
 
 
 class TestRigid(BaseTest):
@@ -114,6 +114,7 @@ def test_label_rigid_bodies_single_partial(self, benzene):
         assert len(list(benzene.rigid_particles())) == 6
         assert len(list(benzene.rigid_particles(rigid_id=0))) == 6
 
+    @pytest.mark.skipif(not has_hoomd, reason="HOOMD not installed.")
     def test_save_non_sequential_rigid_ids(self, benzene):
         n_benzenes = 10
         filled = mb.fill_box(
@@ -122,7 +123,7 @@ def test_save_non_sequential_rigid_ids(self, benzene):
         filled.label_rigid_bodies(discrete_bodies="Benzene")
         filled.children[0]._increment_rigid_ids(increment=3)
         with pytest.warns(UserWarning):
-            filled.save("benzene-box.hoomdxml")
+            filled.save("benzene-box.gsd")
 
     def test_increment_rigid_id(self, rigid_benzene):
         compound = mb.Compound()
diff --git a/mbuild/tests/test_xyz.py b/mbuild/tests/test_xyz.py
deleted file mode 100644
index a567c29a7..000000000
--- a/mbuild/tests/test_xyz.py
+++ /dev/null
@@ -1,72 +0,0 @@
-import ele
-import numpy as np
-import pytest
-
-import mbuild as mb
-from mbuild.exceptions import MBuildError
-from mbuild.formats.xyz import write_xyz
-from mbuild.tests.base_test import BaseTest
-from mbuild.utils.io import get_fn
-
-
-class TestXYZ(BaseTest):
-    def test_load_no_top(self, ethane):
-        ethane.save(filename="ethane.xyz")
-        ethane_in = mb.load("ethane.xyz")
-        assert len(ethane_in.children) == 8
-        assert ethane_in.n_bonds == 0
-        assert set([child.name for child in ethane_in.children]) == {"C", "H"}
-        assert set([child.element for child in ethane_in.children]) == {
-            ele.Elements.C,
-            ele.Elements.H,
-        }
-
-    def test_wrong_n_atoms(self):
-        with pytest.raises(MBuildError):
-            mb.load(get_fn("too_few_atoms.xyz"))
-        with pytest.raises(MBuildError):
-            mb.load(get_fn("too_many_atoms.xyz"))
-
-    def test_bad_input(self, ethane):
-        with pytest.raises(ValueError):
-            assert isinstance(ethane, mb.Compound)
-            write_xyz(ethane, "compound.xyz")
-
-    def test_save(self, ethane):
-        ethane.save(filename="ethane.xyz")
-        ethane_in = mb.load("ethane.xyz")
-        assert len(ethane_in.children) == 8
-        assert ethane_in.n_bonds == 0
-        assert set([child.name for child in ethane_in.children]) == {"C", "H"}
-
-    def test_coordinates(self, ethane):
-        ethane.save(filename="ethane.xyz")
-        ethane_in = mb.load("ethane.xyz")
-        assert np.allclose(ethane.xyz, ethane_in.xyz)
-
-    def test_non_resolved_elements(self):
-        tip3p_water = mb.load(get_fn("tip3p_water.xyz"))
-        assert tip3p_water[0].element is None
-        assert tip3p_water[0].name == "opls_111"
-        assert tip3p_water[1].element is None
-        assert tip3p_water[1].name == "opls_112"
-        assert tip3p_water[2].element is None
-        assert tip3p_water[2].name == "opls_112"
-
-    def test_write_atomtypes(self):
-        tip3p_water = mb.load(get_fn("tip3p_water.xyz"))
-        tip3p_water.save(filename="test.xyz", write_atomnames=True)
-        lines = []
-        with open("test.xyz") as f:
-            for line in f:
-                lines.append(line.split())
-        assert lines[2][0] == "opls_111"
-        assert lines[3][0] == "opls_112"
-        assert lines[4][0] == "opls_112"
-        tip3p_water_in = mb.load("test.xyz")
-        assert tip3p_water_in[0].element is None
-        assert tip3p_water_in[0].name == "opls_111"
-        assert tip3p_water_in[1].element is None
-        assert tip3p_water_in[1].name == "opls_112"
-        assert tip3p_water_in[2].element is None
-        assert tip3p_water_in[2].name == "opls_112"