diff --git a/gmso/core/topology.py b/gmso/core/topology.py
index 6cebb2772..2472361cb 100644
--- a/gmso/core/topology.py
+++ b/gmso/core/topology.py
@@ -30,7 +30,7 @@
     convert_params_units,
     convert_topology_expressions,
 )
-from gmso.utils.units import GMSO_UnitRegsitry as UnitReg
+from gmso.utils.units import GMSO_UnitRegistry as UnitReg
 
 scaling_interaction_idxes = {"12": 0, "13": 1, "14": 2}
 
@@ -170,6 +170,17 @@ def __init__(self, name="Topology", box=None):
         }
 
         self._unique_connections = {}
+        self._unit_system = None
+
+    @property
+    def unit_system(self):
+        """Return the unyt system of the topology."""
+        return self._unit_system
+
+    @unit_system.setter
+    def unit_system(self, unit_system):
+        """Set the unyt system of the topology."""
+        self._name = unit_system
 
     @property
     def name(self):
diff --git a/gmso/formats/lammpsdata.py b/gmso/formats/lammpsdata.py
index 97c865b8d..69ce71a64 100644
--- a/gmso/formats/lammpsdata.py
+++ b/gmso/formats/lammpsdata.py
@@ -11,7 +11,6 @@
 import numpy as np
 import unyt as u
 from sympy import Symbol
-from unyt import UnitRegistry
 from unyt.array import allclose_units
 
 import gmso
@@ -38,215 +37,7 @@
     convert_opls_to_ryckaert,
     convert_ryckaert_to_opls,
 )
-
-# TODO: move this to gmso.utils.units.py
-reg = UnitRegistry()
-dim = u.dimensions.current_mks * u.dimensions.time
-conversion = 1 * getattr(u.physical_constants, "elementary_charge").value
-reg.add(
-    "elementary_charge",
-    base_value=conversion,
-    dimensions=dim,
-    tex_repr=r"\rm{e}",
-)
-conversion = 1 * getattr(u.physical_constants, "boltzmann_constant_mks").value
-dim = u.dimensions.energy / u.dimensions.temperature
-reg.add(
-    "kb", base_value=conversion, dimensions=dim, tex_repr=r"\rm{kb}"
-)  # boltzmann temperature
-conversion = (
-    4
-    * np.pi
-    * getattr(u.physical_constants, "reduced_planck_constant").value ** 2
-    * getattr(u.physical_constants, "eps_0").value
-    / (
-        getattr(u.physical_constants, "electron_charge").value ** 2
-        * getattr(u.physical_constants, "electron_mass").value
-    )
-)
-dim = u.dimensions.length
-reg.add(
-    "a0", base_value=conversion, dimensions=dim, tex_repr=r"\rm{a0}"
-)  # bohr radius
-conversion = (
-    getattr(u.physical_constants, "reduced_planck_constant").value ** 2
-    / u.Unit("a0", registry=reg).base_value ** 2
-    / getattr(u.physical_constants, "electron_mass").value
-)
-dim = u.dimensions.energy
-reg.add(
-    "Ehartree", base_value=conversion, dimensions=dim, tex_repr=r"\rm{Ehartree}"
-)  # Hartree energy
-conversion = np.sqrt(
-    10**9 / (4 * np.pi * getattr(u.physical_constants, "eps_0").value)
-)
-dim = u.dimensions.charge
-reg.add(
-    "Statcoulomb_charge",
-    base_value=conversion,
-    dimensions=dim,
-    tex_repr=r"\rm{Statcoulomb_charge}",
-)  # Static charge
-
-
-def _unit_style_factory(style: str):
-    #  NOTE: the when an angle is measured in lammps is not straightforwards. It depends not on the unit_style, but on the
-    # angle_style, dihedral_style, or improper_style. For examples, harmonic angles, k is specificed in energy/radian, but the
-    # theta_eq is written in degrees. For fourier dihedrals, d_eq is specified in degrees. When adding new styles, make sure that
-    # this behavior is accounted for when converting the specific potential_type in the function
-    # _parameter_converted_to_float
-    if style == "real":
-        base_units = u.UnitSystem(
-            "lammps_real", "Å", "amu", "fs", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "kcal/mol"
-        base_units["charge"] = "elementary_charge"
-    elif style == "metal":
-        base_units = u.UnitSystem(
-            "lammps_metal", "Å", "amu", "picosecond", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "eV"
-        base_units["charge"] = "elementary_charge"
-    elif style == "si":
-        base_units = u.UnitSystem(
-            "lammps_si", "m", "kg", "s", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "joule"
-        base_units["charge"] = "coulomb"
-    elif style == "cgs":
-        base_units = u.UnitSystem(
-            "lammps_cgs", "cm", "g", "s", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "erg"
-        # Statcoulomb is strange. It is not a 1:1 correspondance to charge, with base units of
-        # mass**1/2*length**3/2*time**-1.
-        # However, assuming it is referring to a static charge and not a flux, it can be
-        # converted to coulomb units. See the registry for the unit conversion to Coulombs
-        base_units["charge"] = "Statcoulomb_charge"
-    elif style == "electron":
-        base_units = u.UnitSystem(
-            "lammps_electron", "a0", "amu", "s", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "Ehartree"
-        base_units["charge"] = "elementary_charge"
-    elif style == "micro":
-        base_units = u.UnitSystem(
-            "lammps_micro", "um", "picogram", "us", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "ug*um**2/us**2"
-        base_units["charge"] = "picocoulomb"
-    elif style == "nano":
-        base_units = u.UnitSystem(
-            "lammps_nano", "nm", "attogram", "ns", "K", "rad", registry=reg
-        )
-        base_units["energy"] = "attogram*nm**2/ns**2"
-        base_units["charge"] = "elementary_charge"
-    elif style == "lj":
-        base_units = ljUnitSystem()
-    else:
-        raise NotYetImplementedWarning
-
-    return base_units
-
-
-class ljUnitSystem:
-    """Use this so the empty unitsystem has getitem magic method."""
-
-    def __init__(self):
-        self.registry = reg
-        self.name = "lj"
-
-    def __getitem__(self, items):
-        """Return dimensionless units."""
-        return "dimensionless"
-
-
-def _parameter_converted_to_float(
-    parameter,
-    base_unyts,
-    conversion_factorDict=None,
-    n_decimals=3,
-    name="",
-):
-    """Take a given parameter, and return a float of the parameter in the given style.
-
-    This function will check the base_unyts, which is a unyt.UnitSystem object,
-    and convert the parameter to those units based on its dimensions. It can
-    also generate dimensionless units via normalization from conversion_factorsDict.
-    # TODO: move this to gmso.utils.units.py
-    """
-    # TODO: now I think phi_eq is what is really saved in the improper angle
-    if name in ["theta_eq", "chieq"]:  # eq angle are always in degrees
-        return round(float(parameter.to("degree").value), n_decimals)
-    new_dims = _dimensions_to_energy(parameter.units.dimensions)
-    new_dims = _dimensions_to_charge(new_dims)
-    if conversion_factorDict and isinstance(base_unyts, ljUnitSystem):
-        # multiply object -> split into length, mass, energy, charge -> grab conversion factor from dict
-        # first replace energy for (length)**2*(mass)/(time)**2 u.dimensions.energy. Then iterate through the free symbols
-        # and figure out a way how to add those to the overall conversion factor
-        dim_info = new_dims.as_terms()
-        conversion_factor = 1
-        for exponent, ind_dim in zip(dim_info[0][0][1][1], dim_info[1]):
-            factor = conversion_factorDict.get(
-                ind_dim.name[1:-1], 1
-            )  # replace () in name
-            conversion_factor *= float(factor) ** exponent
-        return float(
-            parameter / conversion_factor
-        )  # Assuming that conversion factor is in right units
-    new_dimStr = str(new_dims)
-    ind_units = re.sub("[^a-zA-Z]+", " ", new_dimStr).split()
-    for unit in ind_units:
-        new_dimStr = new_dimStr.replace(unit, str(base_unyts[unit]))
-
-    return round(
-        float(parameter.to(u.Unit(new_dimStr, registry=base_unyts.registry))),
-        n_decimals,
-    )
-
-
-def _dimensions_to_energy(dims):
-    """Take a set of dimensions and substitute in Symbol("energy") where possible."""
-    # TODO: move this to gmso.utils.units.py
-    symsStr = str(dims.free_symbols)
-    energy_inBool = np.all([dimStr in symsStr for dimStr in ["time", "mass"]])
-    if not energy_inBool:
-        return dims
-    energySym = Symbol("(energy)")  # create dummy symbol to replace in equation
-    dim_info = dims.as_terms()
-    time_idx = np.where(list(map(lambda x: x.name == "(time)", dim_info[1])))[
-        0
-    ][0]
-    energy_exp = (
-        dim_info[0][0][1][1][time_idx] // 2
-    )  # energy has 1/time**2 in it, so this is the hint of how many
-    return (
-        dims
-        * u.dimensions.energy**energy_exp
-        * energySym ** (-1 * energy_exp)
-    )
-
-
-def _dimensions_to_charge(dims):
-    """Take a set of dimensions and substitute in Symbol("charge") where possible."""
-    # TODO: move this to gmso.utils.units.py
-    symsStr = str(dims.free_symbols)
-    charge_inBool = np.all([dimStr in symsStr for dimStr in ["current_mks"]])
-    if not charge_inBool:
-        return dims
-    chargeSym = Symbol("(charge)")  # create dummy symbol to replace in equation
-    dim_info = dims.as_terms()
-    time_idx = np.where(
-        list(map(lambda x: x.name == "(current_mks)", dim_info[1]))
-    )[0][0]
-    charge_exp = dim_info[0][0][1][1][
-        time_idx
-    ]  # charge has (current_mks) in it, so this is the hint of how many
-    return (
-        dims
-        * u.dimensions.charge ** (-1 * charge_exp)
-        * chargeSym**charge_exp
-    )
+from gmso.utils.units import LAMMPS_UnitSystems, write_out_parameter_and_units
 
 
 @saves_as(".lammps", ".lammpsdata", ".data")
@@ -331,7 +122,7 @@ def write_lammpsdata(
         raise ValueError(
             "lj_cfactorsDict argument is only used if unit_style is lj."
         )
-    base_unyts = _unit_style_factory(unit_style)
+    base_unyts = LAMMPS_UnitSystems(unit_style)
     default_parameterMaps = {  # TODO: sites are not checked currently because gmso
         # doesn't store pair potential eqn the same way as the connections.
         "impropers": "HarmonicImproperPotential",
@@ -400,7 +191,9 @@ def write_lammpsdata(
 
 @loads_as(".lammps", ".lammpsdata", ".data")
 def read_lammpsdata(
-    filename, atom_style="full", unit_style="real", potential="lj"
+    filename,
+    atom_style="full",
+    unit_style="real",
 ):
     """Read in a lammps data file as a GMSO topology.
 
@@ -414,8 +207,6 @@ def read_lammpsdata(
     unit_style : str, optional, default='real
         LAMMPS unit style used for writing the datafile. Can be "real", "lj",
         "metal", "si", "cgs", "electron", "micro", "nano".
-    potential: str, optional, default='lj'
-        Potential type defined in data file. Only supporting lj as of now.
 
     Returns
     -------
@@ -434,7 +225,6 @@ def read_lammpsdata(
     "electron", "micro", "nano".
 
     Currently supporting the following potential styles: 'lj'
-
     Currently supporting the following bond styles: 'harmonic'
     Currently supporting the following angle styles: 'harmonic'
     Currently supporting the following dihedral styles: 'opls'
@@ -467,44 +257,44 @@ def read_lammpsdata(
                 unit_style
             )
         )
+    base_unyts = LAMMPS_UnitSystems(unit_style)
 
     # Parse box information
-    _get_box_coordinates(filename, unit_style, top)
+    _get_box_coordinates(filename, base_unyts, top)
     # Parse atom type information
-    top, type_list = _get_ff_information(filename, unit_style, top)
+    top, type_list = _get_ff_information(filename, base_unyts, top)
     # Parse atom information
-    _get_atoms(filename, top, unit_style, type_list)
+    _get_atoms(filename, top, base_unyts, type_list)
     # Parse connection (bonds, angles, dihedrals, impropers) information
     # TODO: Add more atom styles
     if atom_style in ["full"]:
-        _get_connection(filename, top, unit_style, connection_type="bond")
-        _get_connection(filename, top, unit_style, connection_type="angle")
-        _get_connection(filename, top, unit_style, connection_type="dihedral")
-        _get_connection(filename, top, unit_style, connection_type="improper")
+        _get_connection(filename, top, base_unyts, connection_type="bond")
+        _get_connection(filename, top, base_unyts, connection_type="angle")
+        _get_connection(filename, top, base_unyts, connection_type="dihedral")
+        _get_connection(filename, top, base_unyts, connection_type="improper")
 
     top.update_topology()
 
     return top
 
 
-def get_units(unit_style, dimension):
+def get_units(base_unyts, dimension):
     """Get u.Unit for specific LAMMPS unit style with given dimension."""
     # Need separate angle units for harmonic force constant and angle
-    if unit_style == "lj":
+    if base_unyts.usystem.name == "lj":
         if dimension == "angle":
             return u.radian
         return u.dimensionless
 
-    usystem = _unit_style_factory(unit_style)
     if dimension == "angle_eq":
         return (
             u.degree
         )  # LAMMPS specifies different units for some angles, such as equilibrium angles
 
-    return u.Unit(usystem[dimension], registry=reg)
+    return u.Unit(base_unyts.usystem[dimension], registry=base_unyts.reg)
 
 
-def _get_connection(filename, topology, unit_style, connection_type):
+def _get_connection(filename, topology, base_unyts, connection_type):
     """Parse connection types."""
     # TODO: check for other connection types besides the defaults
     with open(filename, "r") as lammps_file:
@@ -528,11 +318,11 @@ def _get_connection(filename, topology, unit_style, connection_type):
             # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term
             conn_params = {
                 "k": float(line.split()[1])
-                * get_units(unit_style, "energy")
-                / get_units(unit_style, "length") ** 2
+                * get_units(base_unyts, "energy")
+                / get_units(base_unyts, "length") ** 2
                 * 2,
                 "r_eq": float(line.split()[2])
-                * get_units(unit_style, "length"),
+                * get_units(base_unyts, "length"),
             }
             name = template_potential.name
             expression = template_potential.expression
@@ -548,11 +338,11 @@ def _get_connection(filename, topology, unit_style, connection_type):
             # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term
             conn_params = {
                 "k": float(line.split()[1])
-                * get_units(unit_style, "energy")
-                / get_units(unit_style, "angle") ** 2
+                * get_units(base_unyts, "energy")
+                / get_units(base_unyts, "angle") ** 2
                 * 2,
                 "theta_eq": float(line.split()[2])
-                * get_units(unit_style, "angle_eq"),
+                * get_units(base_unyts, "angle_eq"),
             }
             name = template_potential.name
             expression = template_potential.expression
@@ -566,10 +356,10 @@ def _get_connection(filename, topology, unit_style, connection_type):
         elif connection_type == "dihedral":
             template_potential = templates["OPLSTorsionPotential"]
             conn_params = {
-                "k1": float(line.split()[1]) * get_units(unit_style, "energy"),
-                "k2": float(line.split()[2]) * get_units(unit_style, "energy"),
-                "k3": float(line.split()[3]) * get_units(unit_style, "energy"),
-                "k4": float(line.split()[4]) * get_units(unit_style, "energy"),
+                "k1": float(line.split()[1]) * get_units(base_unyts, "energy"),
+                "k2": float(line.split()[2]) * get_units(base_unyts, "energy"),
+                "k3": float(line.split()[3]) * get_units(base_unyts, "energy"),
+                "k4": float(line.split()[4]) * get_units(base_unyts, "energy"),
             }
             name = template_potential.name
             expression = template_potential.expression
@@ -584,11 +374,11 @@ def _get_connection(filename, topology, unit_style, connection_type):
             template_potential = templates["HarmonicImproperPotential"]
             conn_params = {
                 "k": float(line.split()[2])
-                * get_units(unit_style, "energy")
-                / get_units(unit_style, "energy") ** 2
+                * get_units(base_unyts, "energy")
+                / get_units(base_unyts, "energy") ** 2
                 * 2,
                 "phi_eq": float(line.split()[3])
-                * get_units(unit_style, "angle_eq"),
+                * get_units(base_unyts, "angle_eq"),
             }
             name = template_potential.name
             expression = template_potential.expression
@@ -650,7 +440,7 @@ def _get_connection(filename, topology, unit_style, connection_type):
     return topology
 
 
-def _get_atoms(filename, topology, unit_style, type_list):
+def _get_atoms(filename, topology, base_unyts, type_list):
     """Parse the atom information in the LAMMPS data file."""
     with open(filename, "r") as lammps_file:
         for i, line in enumerate(lammps_file):
@@ -663,11 +453,11 @@ def _get_atoms(filename, topology, unit_style, type_list):
         atom_line = line.split()
         atom_type = atom_line[2]
         charge = u.unyt_quantity(
-            float(atom_line[3]), get_units(unit_style, "charge")
+            float(atom_line[3]), get_units(base_unyts, "charge")
         )
         coord = u.unyt_array(
             [float(atom_line[4]), float(atom_line[5]), float(atom_line[6])]
-        ) * get_units(unit_style, "length")
+        ) * get_units(base_unyts, "length")
         site = Atom(
             charge=charge,
             position=coord,
@@ -682,7 +472,7 @@ def _get_atoms(filename, topology, unit_style, type_list):
     return topology
 
 
-def _get_box_coordinates(filename, unit_style, topology):
+def _get_box_coordinates(filename, base_unyts, topology):
     """Parse box information."""
     with open(filename, "r") as lammps_file:
         for line in lammps_file:
@@ -723,20 +513,20 @@ def _get_box_coordinates(filename, unit_style, topology):
             gamma = np.arccos(xy / b)
 
             # Box Information
-            lengths = u.unyt_array([a, b, c], get_units(unit_style, "length"))
+            lengths = u.unyt_array([a, b, c], get_units(base_unyts, "length"))
             angles = u.unyt_array(
-                [alpha, beta, gamma], get_units(unit_style, "angle")
+                [alpha, beta, gamma], get_units(base_unyts, "angle")
             )
             topology.box = Box(lengths, angles)
         else:
             # Box Information
-            lengths = u.unyt_array([x, y, z], get_units(unit_style, "length"))
+            lengths = u.unyt_array([x, y, z], get_units(base_unyts, "length"))
             topology.box = Box(lengths)
 
         return topology
 
 
-def _get_ff_information(filename, unit_style, topology):
+def _get_ff_information(filename, base_unyts, topology):
     """Parse atom-type information."""
     with open(filename, "r") as lammps_file:
         types = False
@@ -753,7 +543,7 @@ def _get_ff_information(filename, unit_style, topology):
     for line in mass_lines:
         atom_type = AtomType(
             name=line.split()[0],
-            mass=float(line.split()[1]) * get_units(unit_style, "mass"),
+            mass=float(line.split()[1]) * get_units(base_unyts, "mass"),
         )
         type_list.append(atom_type)
 
@@ -769,10 +559,10 @@ def _get_ff_information(filename, unit_style, topology):
         if len(pair.split()) == 3:
             type_list[i].parameters["sigma"] = float(
                 pair.split()[2]
-            ) * get_units(unit_style, "length")
+            ) * get_units(base_unyts, "length")
             type_list[i].parameters["epsilon"] = float(
                 pair.split()[1]
-            ) * get_units(unit_style, "energy")
+            ) * get_units(base_unyts, "energy")
         elif len(pair.split()) == 4:
             warn_ljcutBool = True
 
@@ -825,12 +615,14 @@ def _validate_unit_compatibility(top, base_unyts):
         ]
         for parameter, name in parametersList:
             assert np.isclose(
-                _parameter_converted_to_float(
-                    parameter, base_unyts, n_decimals=6, name=name
+                float(
+                    base_unyts.convert_parameter(
+                        parameter, n_decimals=6, name=name
+                    )
                 ),
                 parameter.value,
                 atol=1e-3,
-            ), f"Units System {base_unyts} is not compatible with {atype} with value {parameter}"
+            ), f"Units System {base_unyts.usystem} is not compatible with {atype} with value {parameter}"
 
 
 def _write_header(out_file, top, atom_style):
@@ -886,8 +678,8 @@ def _write_box(out_file, top, base_unyts, cfactorsDict):
         atol=1e-8,
     ):
         box_lengths = [
-            _parameter_converted_to_float(
-                top.box.lengths[i], base_unyts, cfactorsDict
+            float(
+                base_unyts.convert_parameter(top.box.lengths[i], cfactorsDict)
             )
             for i in range(3)
         ]
@@ -898,8 +690,8 @@ def _write_box(out_file, top, base_unyts, cfactorsDict):
         out_file.write("0.000000 0.000000 0.000000 xy xz yz\n")
     else:
         box_lengths = [
-            _parameter_converted_to_float(
-                top.box.lengths[i], base_unyts, cfactorsDict
+            float(
+                base_unyts.convert_parameter(top.box.lengths[i], cfactorsDict)
             )
             for i in range(3)
         ]
@@ -941,15 +733,13 @@ def _write_box(out_file, top, base_unyts, cfactorsDict):
 def _write_atomtypes(out_file, top, base_unyts, cfactorsDict):
     """Write out atomtypes in GMSO topology to LAMMPS file."""
     out_file.write("\nMasses\n")
-    out_file.write(f"#\tmass ({base_unyts['mass']})\n")
+    out_file.write(f"#\tmass ({base_unyts.usystem['mass']})\n")
     atypesView = sorted(top.atom_types(filter_by=pfilter), key=lambda x: x.name)
     for atom_type in atypesView:
         out_file.write(
-            "{:d}\t{:.6f}\t# {}\n".format(
+            "{:d}\t{}\t# {}\n".format(
                 atypesView.index(atom_type) + 1,
-                _parameter_converted_to_float(
-                    atom_type.mass, base_unyts, cfactorsDict
-                ),
+                base_unyts.convert_parameter(atom_type.mass, cfactorsDict),
                 atom_type.name,
             )
         )
@@ -966,7 +756,7 @@ def _write_pairtypes(out_file, top, base_unyts, cfactorsDict):
         "sigma",
     )  # this will vary with new pair styles
     param_labels = [
-        _write_out_parameter_w_units(
+        write_out_parameter_and_units(
             key, test_atomtype.parameters[key], base_unyts
         )
         for key in nb_style_orderTuple
@@ -977,14 +767,11 @@ def _write_pairtypes(out_file, top, base_unyts, cfactorsDict):
     )
     for idx, param in enumerate(sorted_atomtypes):
         out_file.write(
-            "{}\t{:7.5f}\t\t{:7.5f}\t\t# {}\n".format(
+            "{}\t{:7}\t\t{:7}\t\t# {}\n".format(
                 idx + 1,
                 *[
-                    _parameter_converted_to_float(
-                        param.parameters[key],
-                        base_unyts,
-                        cfactorsDict,
-                        n_decimals=5,
+                    base_unyts.convert_parameter(
+                        param.parameters[key], cfactorsDict, n_decimals=5
                     )
                     for key in nb_style_orderTuple
                 ],
@@ -1000,7 +787,7 @@ def _write_bondtypes(out_file, top, base_unyts, cfactorsDict):
     out_file.write(f"\nBond Coeffs #{test_bondtype.name}\n")
     bond_style_orderTuple = ("k", "r_eq")
     param_labels = [
-        _write_out_parameter_w_units(
+        write_out_parameter_and_units(
             key, test_bondtype.parameters[key], base_unyts
         )
         for key in bond_style_orderTuple
@@ -1014,11 +801,11 @@ def _write_bondtypes(out_file, top, base_unyts, cfactorsDict):
             [bond_type.member_types[0], bond_type.member_types[1]]
         )
         out_file.write(
-            "{}\t{:7.5f}\t{:7.5f}\t\t# {}\t{}\n".format(
+            "{}\t{:7}\t{:7}\t\t# {}\t{}\n".format(
                 idx + 1,
                 *[
-                    _parameter_converted_to_float(
-                        bond_type.parameters[key], base_unyts, cfactorsDict
+                    base_unyts.convert_parameter(
+                        bond_type.parameters[key], cfactorsDict, n_decimals=6
                     )
                     for key in bond_style_orderTuple
                 ],
@@ -1037,7 +824,7 @@ def _write_angletypes(out_file, top, base_unyts, cfactorsDict):
         "theta_eq",
     )  # this will vary with new angle styles
     param_labels = [
-        _write_out_parameter_w_units(
+        write_out_parameter_and_units(
             key, test_angletype.parameters[key], base_unyts
         )
         for key in angle_style_orderTuple
@@ -1053,13 +840,13 @@ def _write_angletypes(out_file, top, base_unyts, cfactorsDict):
     )
     for idx, angle_type in enumerate(indexList):
         out_file.write(
-            "{}\t{:7.5f}\t{:7.5f}\t#{:11s}\t{:11s}\t{:11s}\n".format(
+            "{}\t{:7}\t{:7}\t#{:11s}\t{:11s}\t{:11s}\n".format(
                 idx + 1,
                 *[
-                    _parameter_converted_to_float(
+                    base_unyts.convert_parameter(
                         angle_type.parameters[key],
-                        base_unyts,
                         cfactorsDict,
+                        n_decimals=6,
                         name=key,
                     )
                     for key in angle_style_orderTuple
@@ -1080,7 +867,7 @@ def _write_dihedraltypes(out_file, top, base_unyts, cfactorsDict):
         "k4",
     )  # this will vary with new dihedral styles
     param_labels = [
-        _write_out_parameter_w_units(
+        write_out_parameter_and_units(
             key, test_dihedraltype.parameters[key], base_unyts
         )
         for key in dihedral_style_orderTuple
@@ -1094,13 +881,14 @@ def _write_dihedraltypes(out_file, top, base_unyts, cfactorsDict):
     index_membersList.sort(key=lambda x: ([x[1][i] for i in [1, 2, 0, 3]]))
     for idx, (dihedral_type, members) in enumerate(index_membersList):
         out_file.write(
-            "{}\t{:8.5f}\t{:8.5f}\t{:8.5f}\t{:8.5f}\t# {}\t{}\t{}\t{}\n".format(
+            "{}\t{:8}\t{:8}\t{:8}\t{:8}\t# {}\t{}\t{}\t{}\n".format(
                 idx + 1,
                 *[
-                    _parameter_converted_to_float(
+                    base_unyts.convert_parameter(
                         dihedral_type.parameters[parameterStr],
-                        base_unyts,
                         cfactorsDict,
+                        n_decimals=6,
+                        name=parameterStr,
                     )
                     for parameterStr in dihedral_style_orderTuple
                 ],
@@ -1119,7 +907,7 @@ def _write_impropertypes(out_file, top, base_unyts, cfactorsDict):
         "phi_eq",
     )  # this will vary with new improper styles
     param_labels = [
-        _write_out_parameter_w_units(
+        write_out_parameter_and_units(
             key, test_impropertype.parameters[key], base_unyts
         )
         for key in improper_style_orderTuple
@@ -1133,13 +921,13 @@ def _write_impropertypes(out_file, top, base_unyts, cfactorsDict):
     index_membersList.sort(key=lambda x: ([x[1][i] for i in [0, 1, 2, 3]]))
     for idx, (improper_type, members) in enumerate(index_membersList):
         out_file.write(
-            "{}\t{:7.5f}\t{:7.5f}\n".format(
+            "{}\t{:7}\t{:7}\n".format(
                 idx + 1,
                 *[
-                    _parameter_converted_to_float(
+                    base_unyts.convert_parameter(
                         improper_type.parameters[parameterStr],
-                        base_unyts,
                         cfactorsDict,
+                        n_decimals=6,
                         name=parameterStr,
                     )
                     for parameterStr in improper_style_orderTuple
@@ -1153,13 +941,15 @@ def _write_site_data(out_file, top, atom_style, base_unyts, cfactorsDict):
     """Write atomic positions and charges to LAMMPS file.."""
     out_file.write(f"\nAtoms #{atom_style}\n\n")
     if atom_style == "atomic":
-        atom_line = "{index:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
+        atom_line = "{index:d}\t{type_index:d}\t{x:.8}\t{y:.8}\t{z:.8}\n"
     elif atom_style == "charge":
-        atom_line = "{index:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
+        atom_line = (
+            "{index:d}\t{type_index:d}\t{charge:.8}\t{x:.8}\t{y:.8}\t{z:.8}\n"
+        )
     elif atom_style == "molecular":
-        atom_line = "{index:d}\t{moleculeid:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
+        atom_line = "{index:d}\t{moleculeid:d}\t{type_index:d}\t{x:.8}\t{y:.8}\t{z:.8}\n"
     elif atom_style == "full":
-        atom_line = "{index:d}\t{moleculeid:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n"
+        atom_line = "{index:d}\t{moleculeid:d}\t{type_index:d}\t{charge:.8}\t{x:.8}\t{y:.8}\t{z:.8}\n"
 
     unique_sorted_typesList = sorted(
         top.atom_types(filter_by=pfilter), key=lambda x: x.name
@@ -1170,17 +960,21 @@ def _write_site_data(out_file, top, atom_style, base_unyts, cfactorsDict):
                 index=i + 1,
                 moleculeid=site.molecule.number,
                 type_index=unique_sorted_typesList.index(site.atom_type) + 1,
-                charge=_parameter_converted_to_float(
-                    site.charge, base_unyts, cfactorsDict
+                charge=base_unyts.convert_parameter(
+                    site.charge,
+                    cfactorsDict,
+                    n_decimals=6,
                 ),
-                x=_parameter_converted_to_float(
-                    site.position[0], base_unyts, cfactorsDict, n_decimals=6
+                x=base_unyts.convert_parameter(
+                    site.position[0],
+                    cfactorsDict,
+                    n_decimals=6,
                 ),
-                y=_parameter_converted_to_float(
-                    site.position[1], base_unyts, cfactorsDict, n_decimals=6
+                y=base_unyts.convert_parameter(
+                    site.position[1], cfactorsDict, n_decimals=6
                 ),
-                z=_parameter_converted_to_float(
-                    site.position[2], base_unyts, cfactorsDict, n_decimals=6
+                z=base_unyts.convert_parameter(
+                    site.position[2], cfactorsDict, n_decimals=6
                 ),
             )
         )
@@ -1259,21 +1053,3 @@ def _default_lj_val(top, source):
         raise ValueError(
             f"Provided {source} for default LJ cannot be found in the topology."
         )
-
-
-def _write_out_parameter_w_units(parameter_name, parameter, base_unyts):
-    if parameter_name in ["theta_eq", "phi_eq"]:
-        return f"{parameter_name} ({'degrees'})"
-    if base_unyts.name == "lj":
-        return f"{parameter_name} ({'dimensionless'})"
-    new_dims = _dimensions_to_energy(parameter.units.dimensions)
-    new_dims = _dimensions_to_charge(new_dims)
-    new_dimStr = str(new_dims)
-    ind_units = re.sub("[^a-zA-Z]+", " ", new_dimStr).split()
-    for unit in ind_units:
-        new_dimStr = new_dimStr.replace(unit, str(base_unyts[unit]))
-
-    outputUnyt = str(
-        parameter.to(u.Unit(new_dimStr, registry=base_unyts.registry)).units
-    )
-    return f"{parameter_name} ({outputUnyt})"
diff --git a/gmso/tests/test_conversions.py b/gmso/tests/test_conversions.py
index cf1704949..539f708d0 100644
--- a/gmso/tests/test_conversions.py
+++ b/gmso/tests/test_conversions.py
@@ -189,13 +189,13 @@ def test_conversion_for_topology_sites(self, typed_ethane):
         ].units == u.Unit("kcal/mol")
 
     def test_lammps_dimensions_to_energy(self):
-        from gmso.formats.lammpsdata import _dimensions_to_energy
+        from gmso.utils.units import LAMMPS_UnitSystems
 
         units = u.Unit("kg")
-        outdims = _dimensions_to_energy(units.dimensions)
+        outdims = LAMMPS_UnitSystems._dimensions_to_energy(units.dimensions)
         assert outdims == units.dimensions == u.dimensions.mass
         units = u.Unit("J")
-        outdims = _dimensions_to_energy(units.dimensions)
+        outdims = LAMMPS_UnitSystems._dimensions_to_energy(units.dimensions)
         assert outdims == sympy.Symbol("(energy)")
         assert (
             units.dimensions
@@ -204,7 +204,7 @@ def test_lammps_dimensions_to_energy(self):
             / u.dimensions.time**2
         )
         units = u.Unit("kcal/nm")
-        outdims = _dimensions_to_energy(units.dimensions)
+        outdims = LAMMPS_UnitSystems._dimensions_to_energy(units.dimensions)
         assert outdims == sympy.Symbol("(energy)") / u.dimensions.length
         assert (
             units.dimensions
@@ -212,46 +212,40 @@ def test_lammps_dimensions_to_energy(self):
         )
 
     def test_lammps_conversion_parameters_base_units(self):
-        from gmso.formats.lammpsdata import (
-            _parameter_converted_to_float,
-            _unit_style_factory,
-        )
+        from gmso.utils.units import LAMMPS_UnitSystems
 
         parameter = 100 * u.Unit("kcal/mol*fs/Å")
-        base_unyts = _unit_style_factory(
+        base_unyts = LAMMPS_UnitSystems(
             "real"
         )  # "lammps_real", "Å", "amu", "fs", "K", "rad"
-        float_param = _parameter_converted_to_float(
-            parameter, base_unyts, conversion_factorDict=None
+        paramStr = base_unyts.convert_parameter(
+            parameter, conversion_factorDict=None
         )
-        assert float_param == 100
+        assert paramStr == "100.000"
         parameter = 100 * u.Unit("K*fs/amu/nm")
-        float_param = _parameter_converted_to_float(
-            parameter, base_unyts, conversion_factorDict=None
+        paramStr = base_unyts.convert_parameter(
+            parameter, conversion_factorDict=None
         )
-        assert float_param == 10
+        assert paramStr == "10.000"
         parameter = 100 * u.Unit("km*g*ms*kJ*degree")
-        base_unyts = _unit_style_factory(
+        base_unyts = LAMMPS_UnitSystems(
             "si"
         )  # "lammps_si", "m", "kg", "s", "K", "rad",
-        float_param = _parameter_converted_to_float(
-            parameter, base_unyts, conversion_factorDict=None, n_decimals=6
+        paramStr = base_unyts.convert_parameter(
+            parameter, conversion_factorDict=None, n_decimals=6
         )
-        assert float_param == round(100 * np.pi / 180, 6)
+        assert paramStr == str(round(100 * np.pi / 180, 6))
         parameter = 1 * u.Unit("g*kJ*Coulomb*m*degree")
-        base_unyts = _unit_style_factory(
+        base_unyts = LAMMPS_UnitSystems(
             "si"
         )  # "lammps_si", "m", "kg", "s", "K", "rad"
-        float_param = _parameter_converted_to_float(
-            parameter, base_unyts, conversion_factorDict=None, n_decimals=6
+        paramStr = base_unyts.convert_parameter(
+            parameter, conversion_factorDict=None, n_decimals=6
         )
-        assert np.isclose(float_param, np.pi / 180, 1e-3)
+        assert np.isclose(float(paramStr), np.pi / 180, 1e-3)
 
     def test_lammps_conversion_parameters_lj(self):
-        from gmso.formats.lammpsdata import (
-            _parameter_converted_to_float,
-            _unit_style_factory,
-        )
+        from gmso.utils.units import LAMMPS_UnitSystems
 
         parameter = 1 * u.Unit("g*kJ*Coulomb*m*degree")
         conversion_factorDict = {
@@ -260,11 +254,10 @@ def test_lammps_conversion_parameters_lj(self):
             "charge": 3 * u.Unit("Coulomb"),
             "length": 3 * u.Unit("m"),
         }
-        base_unyts = _unit_style_factory("lj")
-        float_param = _parameter_converted_to_float(
+        base_unyts = LAMMPS_UnitSystems("lj")
+        paramStr = base_unyts.convert_parameter(
             parameter,
-            base_unyts,
             conversion_factorDict=conversion_factorDict,
             n_decimals=6,
         )
-        assert np.isclose(float_param, 1 / 3**4, atol=1e-6)
+        assert np.isclose(float(paramStr), 1 / 3**4, atol=1e-6)
diff --git a/gmso/tests/test_lammps.py b/gmso/tests/test_lammps.py
index 625465e9b..882a1554e 100644
--- a/gmso/tests/test_lammps.py
+++ b/gmso/tests/test_lammps.py
@@ -14,7 +14,6 @@
 from gmso.exceptions import EngineIncompatibilityError
 from gmso.external import from_parmed, to_parmed
 from gmso.formats.formats_registry import UnsupportedFileFormatError
-from gmso.formats.lammpsdata import read_lammpsdata, write_lammpsdata
 from gmso.tests.base_test import BaseTest
 from gmso.tests.utils import get_path
 
@@ -346,7 +345,7 @@ def test_lammps_default_conversions(
         assert lines[38:41] == [
             "Dihedral Coeffs #OPLSTorsionPotential\n",
             "#\tk1 (kcal/mol)\tk2 (kcal/mol)\tk3 (kcal/mol)\tk4 (kcal/mol)\n",
-            "1\t 0.00000\t-0.00000\t 0.30000\t-0.00000\t# opls_140\topls_135\topls_135\topls_140\n",
+            "1\t0.000000\t-0.000000\t0.300000\t-0.000000\t# opls_140\topls_135\topls_135\topls_140\n",
         ]
 
         struc = harmonic_parmed_types_charmm
@@ -419,6 +418,9 @@ def test_lammps_units(self, typed_ethane, unit_style):
         """
         # check the initial set of units
         from gmso.formats.lammpsdata import get_units
+        from gmso.utils.units import LAMMPS_UnitSystems
+
+        base_unyts = LAMMPS_UnitSystems(unit_style)
 
         # real units should be in: [g/mol, angstroms, fs, kcal/mol, kelvin, electon charge, ...]
         mass_multiplierDict = {
@@ -434,10 +436,10 @@ def test_lammps_units(self, typed_ethane, unit_style):
                 atype.mass = 12 * mass_multiplierDict[unit_style] * u.kg
         typed_ethane.save("ethane.lammps", unit_style=unit_style)
         real_top = Topology().load("ethane.lammps", unit_style=unit_style)
-        energy_unit = get_units(unit_style, "energy")
-        angle_unit = get_units(unit_style, "angle_eq")
-        length_unit = get_units(unit_style, "length")
-        charge_unit = get_units(unit_style, "charge")
+        energy_unit = get_units(base_unyts, "energy")
+        angle_unit = get_units(base_unyts, "angle_eq")
+        length_unit = get_units(base_unyts, "length")
+        charge_unit = get_units(base_unyts, "charge")
         assert (
             real_top.dihedrals[0].dihedral_type.parameters["k1"].units
             == energy_unit
@@ -490,8 +492,6 @@ def test_lammps_units(self, typed_ethane, unit_style):
                 atol=1e-8,
             )
 
-    from gmso.exceptions import EngineIncompatibilityError
-
     def test_lammps_errors(self, typed_ethane):
         with pytest.raises(UnsupportedFileFormatError):
             typed_ethane.save("e.lammmps")
@@ -513,15 +513,15 @@ def test_lammps_errors(self, typed_ethane):
         with pytest.raises(ValueError):
             typed_ethane.save("error.lammps", unit_style="None")
 
-    def test_lammps_units(self, typed_methylnitroaniline):
+    def test_lammps_validate_units(self, typed_methylnitroaniline):
         from gmso.formats.lammpsdata import _validate_unit_compatibility
+        from gmso.utils.units import LAMMPS_UnitSystems
 
-        usys = u.unit_systems.mks_unit_system
+        base_unyts = LAMMPS_UnitSystems("si")
         with pytest.raises(AssertionError):
-            _validate_unit_compatibility(typed_methylnitroaniline, usys)
-        from gmso.formats.lammpsdata import _unit_style_factory
+            _validate_unit_compatibility(typed_methylnitroaniline, base_unyts)
 
-        usys = _unit_style_factory("real")
+        usys = LAMMPS_UnitSystems("real")
         _validate_unit_compatibility(typed_methylnitroaniline, usys)
 
     def test_units_in_headers(self, typed_ethane):
@@ -599,7 +599,7 @@ def test_lj_passed_units(self, typed_ethane):
         assert largest_eps_written == 0.5
 
     def test_unit_style_factor(self):
-        from gmso.formats.lammpsdata import _unit_style_factory
+        from gmso.utils.units import LAMMPS_UnitSystems
 
         for styleStr in [
             "real",
@@ -610,8 +610,11 @@ def test_unit_style_factor(self):
             "micro",
             "nano",
         ]:
-            assert _unit_style_factory(styleStr).name == "lammps_" + styleStr
+            assert (
+                LAMMPS_UnitSystems(styleStr).usystem.name
+                == "lammps_" + styleStr
+            )
         from gmso.exceptions import NotYetImplementedWarning
 
         with pytest.raises(NotYetImplementedWarning):
-            _unit_style_factory("None")
+            LAMMPS_UnitSystems("None")
diff --git a/gmso/tests/test_topology.py b/gmso/tests/test_topology.py
index 70ec8c14b..54fd29d7e 100644
--- a/gmso/tests/test_topology.py
+++ b/gmso/tests/test_topology.py
@@ -21,7 +21,7 @@
 from gmso.external.convert_parmed import from_parmed
 from gmso.tests.base_test import BaseTest
 from gmso.utils.io import get_fn, has_pandas, has_parmed, import_
-from gmso.utils.units import GMSO_UnitRegsitry as UnitReg
+from gmso.utils.units import GMSO_UnitRegistry as UnitReg
 
 if has_parmed:
     pmd = import_("parmed")
diff --git a/gmso/tests/test_units.py b/gmso/tests/test_units.py
new file mode 100644
index 000000000..ff2c71fe4
--- /dev/null
+++ b/gmso/tests/test_units.py
@@ -0,0 +1,243 @@
+import re
+
+import pytest
+import unyt as u
+
+from gmso.tests.base_test import BaseTest
+from gmso.utils.misc import unyt_to_hashable
+from gmso.utils.units import LAMMPS_UnitSystems
+
+
+class TestUnitHandling(BaseTest):
+    @pytest.fixture
+    def real_usys(self):
+        return LAMMPS_UnitSystems("real")
+
+    def test_unyt_to_hashable(self):
+        hash(unyt_to_hashable(None))
+        hash(unyt_to_hashable(1 * u.nm))
+        hash(unyt_to_hashable([4, 4] * u.nm))
+
+        assert hash(unyt_to_hashable(1 * u.nm)) == hash(
+            unyt_to_hashable(10 * u.angstrom)
+        )
+        assert hash(unyt_to_hashable(1 * u.kg)) == hash(
+            unyt_to_hashable(1000 * u.g)
+        )
+
+        assert hash(unyt_to_hashable(1 * u.nm)) != hash(
+            unyt_to_hashable(1.01 * u.nm)
+        )
+        assert hash(unyt_to_hashable(1 * u.nm)) != hash(
+            unyt_to_hashable(1.01 * u.second)
+        )
+        assert hash(unyt_to_hashable(1 * u.nm)) != hash(
+            unyt_to_hashable([1, 1] * u.nm)
+        )
+
+    """
+    Utilities to make
+    [a] register units needed for unit systems
+    [a] need to be able to generate unit systems
+    [a] take a unyt and check for energy
+    [a] take a unyt and check for electron volts
+    [?] take a unyt and look for thermal energy
+    [a] get all dimensions from a unit
+    [a] convert a unit using a base system
+    [a] return a rounded float for a unit
+    [] be able to write out a unit without the conversion
+    [] attach units to a float with a unit system and dimensions
+    [] should have associated somewhere to handle units outside of unit system
+    # units should be done as a top level conversion
+    # need to implement this module into the writers
+    # can probably subclass some functions out of LAMMPS_UnitSystems
+
+    Tests to make
+
+    """
+
+    def test_unit_conversion(self, real_usys):
+        parameter = 0.001 * u.Unit("mm")
+        n_decimals = 5
+        outStr = real_usys.convert_parameter(parameter, n_decimals=n_decimals)
+        assert float(outStr) == 10000.00000
+        assert outStr[::-1].find(".") == n_decimals
+
+        parameter = 1 * u.Unit("nm")
+        n_decimals = 5
+        outStr = real_usys.convert_parameter(parameter, n_decimals=n_decimals)
+        assert float(outStr) == 10.00000
+        assert outStr[::-1].find(".") == n_decimals
+
+    def test_unit_rounding(self, real_usys):
+        parameter = 0.001 * u.Unit("nm")
+        n_decimals = 5
+        outStr = real_usys.convert_parameter(parameter, n_decimals=n_decimals)
+        assert outStr[::-1].find(".") == n_decimals
+
+    def test_unitsystem_setup(self, real_usys):
+        assert real_usys.usystem.name == "lammps_real"
+
+        usys = LAMMPS_UnitSystems("lj", registry=u.UnitRegistry())
+        assert usys.usystem.name == "lj"
+
+    def test_dimensions_to_energy(self, real_usys):
+        real_usys = LAMMPS_UnitSystems("real")
+        parameter = 1 * u.kJ / u.nm * u.g
+        # Note: parameter cannot divide out mass or time from energy
+        out_parameter = real_usys._dimensions_to_energy(
+            parameter.units.dimensions
+        )
+        assert str(out_parameter) == "(energy)*(mass)/(length)"
+
+    def test_dimensions_to_charge(self, real_usys):
+        parameter = 1 * u.coulomb / u.nm
+        out_parameter = real_usys._dimensions_to_charge(
+            parameter.units.dimensions
+        )
+        assert str(out_parameter) == "(charge)/(length)"
+
+    def test_dimensions_thermal(self, real_usys):
+        parameter = 1 * u.K
+        out_parameter = real_usys._dimensions_from_thermal_to_energy(
+            parameter.units.dimensions
+        )
+        assert str(out_parameter) == "(energy)"
+
+    def test_get_dimensions(self):
+        usys = LAMMPS_UnitSystems("electron")
+        parametersList = list(
+            map(
+                lambda x: 1 * u.Unit(x, registry=usys.reg),
+                [
+                    "nm",
+                    "kJ",
+                    "kJ/mol",
+                    "K",
+                    "degree/angstrom",
+                    "elementary_charge/mm",
+                    "dimensionless",
+                    "kg*m**2/s**2",
+                    "coulomb",
+                    "kcal/nm**2",
+                    "K/nm",
+                ],
+            )
+        )
+
+        output_dimensionsList = [
+            "length",
+            "energy",
+            "energy",
+            "temperature",
+            "angle/length",
+            "charge/length",
+            "dimensionless",
+            "energy",
+            "charge",
+            "energy/length**2",
+            "energy/length",
+        ]
+        for parameter, dim in zip(parametersList, output_dimensionsList):
+            if str(parameter.units) == "K/nm":
+                thermalize = True
+            else:
+                thermalize = False
+            dimsStr = str(
+                usys._get_output_dimensions(
+                    parameter.units.dimensions, thermalize
+                )
+            )
+            remove_parStr = dimsStr.translate({ord(i): None for i in "()"})
+            assert remove_parStr == str(dim)
+
+    def test_convert_parameters(self, typed_ethane, real_usys):
+        parameter = typed_ethane.sites[0].atom_type.parameters["epsilon"]
+        assert real_usys.convert_parameter(parameter) == "0.066"
+        real_usys.usystem["energy"] = u.kJ / u.mol
+        assert (
+            real_usys.convert_parameter(parameter, n_decimals=6) == "0.276144"
+        )
+        usys = LAMMPS_UnitSystems("real")
+        parameter = typed_ethane.bonds[0].bond_type.parameters["k"]
+        assert usys.convert_parameter(parameter, n_decimals=0) == "680"
+        parameter = typed_ethane.angles[0].angle_type.parameters["theta_eq"]
+        assert usys.convert_parameter(parameter, name="theta_eq") == "110.700"
+        parameter = typed_ethane.dihedrals[0].dihedral_type.parameters["c0"]
+        assert usys.convert_parameter(parameter) == "0.150"
+
+    def test_get_parameter_dimension(self):
+        from gmso.utils.units import get_parameter_dimension
+
+        assert (
+            get_parameter_dimension(1 * u.kJ / u.mol / u.nm, "(energy)")
+            == u.kJ / u.mol
+        )
+        assert get_parameter_dimension(1 * u.kJ / u.nm, "(length)") == u.nm
+        assert (
+            get_parameter_dimension(1 * u.kJ / u.mol / u.nm, "(length)") == u.nm
+        )
+
+    def test_convert_to_unit_system(self):
+        # TODO: discuss if we would want a function to convert a whole
+        # topology at once.
+        # convert a whole topology to a specific unit system
+        pass
+
+    def test_generate_unit_styles(self):
+        # TODO: write all unit systems for these engines.
+        # look at libary of unit styles for lammps, gromacs, hoomd, gomc
+        pass
+
+    def test_lj_units(self, typed_ethane):
+        # write out unit styles from ljUnitSystem and a dictonary of non-dimesnional values
+        lj_usys = LAMMPS_UnitSystems("lj")
+        bond_parameter = typed_ethane.bonds[0].bond_type.parameters["k"]
+        errorStr = (
+            "Missing conversion_factorDict for a dimensionless unit system."
+        )
+        with pytest.raises(ValueError, match=errorStr):
+            lj_usys.convert_parameter(
+                bond_parameter, conversion_factorDict=None
+            )
+        cfactorDict = {"energy": 0.276144 * u.kJ / u.mol, "length": 0.35 * u.nm}
+        errorStr = f"Missing dimensionless constant in conversion_factorDict {cfactorDict}"
+        with pytest.raises(ValueError, match=re.escape(errorStr)):
+            lj_usys.convert_parameter(
+                bond_parameter, conversion_factorDict=cfactorDict
+            )
+        cfactorDict["charge"] = 1 * u.coulomb
+        cfactorDict["mass"] = 12.011 * u.amu
+        outStr = lj_usys.convert_parameter(
+            bond_parameter, conversion_factorDict=cfactorDict
+        )
+        assert outStr == str(round(284512.0 / 0.276144 * 0.35**2, 3))
+        outStr = lj_usys.convert_parameter(
+            typed_ethane.sites[0].atom_type.mass,
+            conversion_factorDict=cfactorDict,
+        )
+        assert outStr == f"{(round(1.000, 3)):.3f}"
+
+    def test_get_units(self, typed_ethane, real_usys):
+        # get the unit system used for a topology
+        typed_ethane.usystem = real_usys
+        assert typed_ethane.usystem == real_usys
+
+    def test_charmm_weighting_factors(self):
+        # write out dihedrals while taking into account weighting
+        pass
+
+    def test_parameter_and_units_writing(self, real_usys):
+        from gmso.utils.units import write_out_parameter_and_units
+
+        x = 1 * u.kJ / u.mol
+        outStr = write_out_parameter_and_units("x", x, real_usys)
+        assert outStr == "x (kcal/mol)"
+
+        x = 1 * u.rad
+        outStr = write_out_parameter_and_units("theta_eq", x, real_usys)
+        assert outStr == "theta_eq (degrees)"
+
+        lj_usys = LAMMPS_UnitSystems("lj")
+        outStr = write_out_parameter_and_units("x", x, lj_usys)
+        assert outStr == "x (dimensionless)"
diff --git a/gmso/utils/units.py b/gmso/utils/units.py
index 066d0b23d..6220038c5 100644
--- a/gmso/utils/units.py
+++ b/gmso/utils/units.py
@@ -1,10 +1,15 @@
 """Source of available units registered within GMSO."""
 
+import re
+
 import numpy as np
 import unyt as u
+from sympy import Symbol
+
+from gmso.exceptions import NotYetImplementedWarning
 
 
-class GMSO_UnitRegsitry(object):
+class GMSO_UnitRegistry(object):
     """A default unit registry class.
 
     The basic units that need to be added for various unit conversions done
@@ -18,15 +23,7 @@ class GMSO_UnitRegsitry(object):
 
     def __init__(self):
         self.reg_ = u.UnitRegistry()
-        conversion = (
-            1 * getattr(u.physical_constants, "elementary_charge").value
-        )
-        self.register_unit(
-            "elementary_charge",
-            conversion,
-            [u.dimensions.current_mks, u.dimensions.time],
-            r"\rm{e}",
-        )
+        register_general_units(self.reg)
 
     def register_unit(
         self,
@@ -79,12 +76,424 @@ def default_reg():
             A unyt registry with commonly used conversions defined.
         """
         reg = u.UnitRegistry()
-        conversion = (
-            1 * getattr(u.physical_constants, "elementary_charge").value
-        )
-        dimensionsList = [u.dimensions.current_mks, u.dimensions.time]
-        dim = np.prod(dimensionsList)
-        name = "elementary_charge"
-        symbol = r"\rm{e}"
-        reg.add(name, conversion, dim, symbol)
+        register_general_units(reg)
         return reg
+
+
+def register_general_units(reg: u.UnitRegistry):
+    """Register units that are generally useful to a basic unyt UnitSystem."""
+    conversion = 1 * getattr(u.physical_constants, "elementary_charge").value
+    dim = u.dimensions.current_mks * u.dimensions.time
+    reg.add(
+        "elementary_charge",
+        conversion,
+        dim,
+        r"\rm{e}",
+    )  # proton charge
+    conversion = (
+        1 * getattr(u.physical_constants, "boltzmann_constant_mks").value
+    )
+    dim = u.dimensions.energy / u.dimensions.temperature
+    reg.add(
+        "kb", base_value=conversion, dimensions=dim, tex_repr=r"\rm{kb}"
+    )  # boltzmann temperature
+    conversion = (
+        4
+        * np.pi
+        * getattr(u.physical_constants, "reduced_planck_constant").value ** 2
+        * getattr(u.physical_constants, "eps_0").value
+        / (
+            getattr(u.physical_constants, "electron_charge").value ** 2
+            * getattr(u.physical_constants, "electron_mass").value
+        )
+    )
+    dim = u.dimensions.length
+    reg.add(
+        "a0", base_value=conversion, dimensions=dim, tex_repr=r"\rm{a0}"
+    )  # bohr radius
+    conversion = (
+        getattr(u.physical_constants, "reduced_planck_constant").value ** 2
+        / u.Unit("a0", registry=reg).base_value ** 2
+        / getattr(u.physical_constants, "electron_mass").value
+    )
+    dim = u.dimensions.energy
+    reg.add(
+        "Ehartree",
+        base_value=conversion,
+        dimensions=dim,
+        tex_repr=r"\rm{Ehartree}",
+    )  # Hartree energy
+    conversion = np.sqrt(
+        10**9 / (4 * np.pi * getattr(u.physical_constants, "eps_0").value)
+    )
+    dim = u.dimensions.charge
+    reg.add(
+        "Statcoulomb_charge",
+        base_value=conversion,
+        dimensions=dim,
+        tex_repr=r"\rm{Statcoulomb_charge}",
+    )  # Static charge
+
+
+class LAMMPS_UnitSystems:
+    """Set of a unit systems distributed in LAMMPS (https://docs.lammps.org/units.html)."""
+
+    def __init__(self, style: str, registry=None):
+        if registry:
+            self.reg_ = registry
+        else:
+            self.reg_ = GMSO_UnitRegistry().reg
+        self.usystem_ = self.usystem_from_str(styleStr=style, reg=self.reg_)
+
+    @property
+    def usystem(self):
+        """Return the UnitSystem attribute for the class."""
+        return self.__dict__.get("usystem_")
+
+    @property
+    def reg(self):
+        """Return the UnytRegistry attribute for the class."""
+        return self.__dict__.get("reg_")
+
+    def usystem_from_str(self, styleStr: str, reg: u.UnitRegistry):
+        """Get systems for unit style."""
+        #  NOTE: the when an angle is measured in lammps is not straightforwards. It depends not on the unit_style, but on the
+        # angle_style, dihedral_style, or improper_style. For examples, harmonic angles, k is specificed in energy/radian, but the
+        # theta_eq is written in degrees. For fourier dihedrals, d_eq is specified in degrees. When adding new styles, make sure that
+        # this behavior is accounted for when converting the specific potential_type in the function
+        # _parameter_converted_to_float
+        if styleStr == "real":
+            base_units = u.UnitSystem(
+                "lammps_real",
+                length_unit="angstrom",
+                mass_unit="amu",
+                time_unit="fs",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "kcal/mol"
+            base_units["charge"] = "elementary_charge"
+        elif styleStr == "metal":
+            base_units = u.UnitSystem(
+                "lammps_metal",
+                length_unit="angstrom",
+                mass_unit="amu",
+                time_unit="picosecond",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "eV"
+            base_units["charge"] = "elementary_charge"
+        elif styleStr == "si":
+            base_units = u.UnitSystem(
+                "lammps_si",
+                length_unit="m",
+                mass_unit="kg",
+                time_unit="s",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "joule"
+            base_units["charge"] = "coulomb"
+        elif styleStr == "cgs":
+            base_units = u.UnitSystem(
+                "lammps_cgs",
+                length_unit="cm",
+                mass_unit="g",
+                time_unit="s",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "erg"
+            # Statcoulomb is strange. It is not a 1:1 correspondance to charge, with base units of
+            # mass**1/2*length**3/2*time**-1.
+            # However, assuming it is referring to a static charge and not a flux, it can be
+            # converted to coulomb units. See the registry for the unit conversion to Coulombs
+            base_units["charge"] = "Statcoulomb_charge"
+        elif styleStr == "electron":
+            base_units = u.UnitSystem(
+                "lammps_electron",
+                length_unit="a0",
+                mass_unit="amu",
+                time_unit="s",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "Ehartree"
+            base_units["charge"] = "elementary_charge"
+        elif styleStr == "micro":
+            base_units = u.UnitSystem(
+                "lammps_micro",
+                length_unit="um",
+                mass_unit="picogram",
+                time_unit="us",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "ug*um**2/us**2"
+            base_units["charge"] = "picocoulomb"
+        elif styleStr == "nano":
+            base_units = u.UnitSystem(
+                "lammps_nano",
+                length_unit="nm",
+                mass_unit="attogram",
+                time_unit="ns",
+                temperature_unit="K",
+                angle_unit="rad",
+                registry=reg,
+            )
+            base_units["energy"] = "attogram*nm**2/ns**2"
+            base_units["charge"] = "elementary_charge"
+        elif styleStr == "lj":
+            base_units = ljUnitSystem(reg)
+        else:
+            raise NotYetImplementedWarning
+
+        return base_units
+
+    def convert_parameter(
+        self,
+        parameter,
+        conversion_factorDict=None,
+        n_decimals=3,
+        name="",
+    ):
+        """Take a given parameter, and return a string of the parameter in the given style.
+
+        This function will check the base_unyts, which is a unyt.UnitSystem object,
+        and convert the parameter to those units based on its dimensions. It can
+        also generate dimensionless units via normalization from conversion_factorsDict.
+
+        Parameters
+        ----------
+        parameter : unyt.array.unyt_quantity
+            Any parameter to convert to a string in the dimensions of self.usystem
+        conversion_factorDict : dict, default=None
+            If the self.usystem is ljUnitSystem, handle conversion
+        n_decimals : int, default=3
+            The number of decimals used in string .f formatting
+        name : string, default=""
+            Additionally information about the parameter, required if handling a specific parameter
+            differently than the default self.usystem would.
+
+        Returns
+        -------
+        outStr : str
+            The parameter converted via self.usystem, and foramtted as a float string.
+        """
+        if name in [
+            "theta_eq",
+            "chieq",
+            "phi_eq",
+        ]:  # eq angle are always in degrees
+            return f"{round(float(parameter.to('degree').value), n_decimals):.{n_decimals}f}"
+        new_dims = self._get_output_dimensions(parameter.units.dimensions)
+        if isinstance(self.usystem, ljUnitSystem):
+            if not conversion_factorDict:
+                raise ValueError(
+                    "Missing conversion_factorDict for a dimensionless unit system."
+                )
+            elif not np.all(
+                [
+                    key in conversion_factorDict
+                    for key in ["energy", "length", "mass", "charge"]
+                ]
+            ):
+                raise ValueError(
+                    f"Missing dimensionless constant in conversion_factorDict {conversion_factorDict}"
+                )
+            # multiply object -> split into length, mass, energy, charge -> grab conversion factor from dict
+            # first replace energy for (length)**2*(mass)/(time)**2 u.dimensions.energy. Then iterate through the free symbols
+            # and figure out a way how to add those to the overall conversion factor
+            dim_info = new_dims.as_terms()
+            conversion_factor = 1
+            for exponent, ind_dim in zip(dim_info[0][0][1][1], dim_info[1]):
+                factor = conversion_factorDict.get(
+                    ind_dim.name[1:-1],
+                    1 * self.usystem[ind_dim.name[1:-1]],  # default value of 1
+                )  # replace () in name
+                current_unit = get_parameter_dimension(parameter, ind_dim.name)
+                factor = factor.to(
+                    current_unit
+                )  # convert factor to units of parameter
+                conversion_factor *= float(factor) ** (exponent)
+            return f"""{round(
+                float(parameter / conversion_factor),
+                n_decimals
+            ):.{n_decimals}f}"""  # Assuming that conversion factor is in right units
+        new_dimStr = str(new_dims)
+        ind_units = re.sub("[^a-zA-Z]+", " ", new_dimStr).split()
+        for unit in ind_units:
+            new_dimStr = new_dimStr.replace(unit, str(self.usystem[unit]))
+        outFloat = float(
+            parameter.to(u.Unit(new_dimStr, registry=self.usystem.registry))
+        )
+
+        return f"{outFloat:.{n_decimals}f}"
+
+    @staticmethod
+    def _dimensions_to_energy(dims):
+        """Take a set of dimensions and substitute in Symbol("energy") where possible."""
+        symsStr = str(dims.free_symbols)
+        energy_inBool = np.all(
+            [dimStr in symsStr for dimStr in ["time", "mass"]]
+        )  # TODO: this logic could be improved, it might fail on complex
+        # units where the units are energy/mass/time**2, or something where the
+        # dimensions are cancelled out
+        if not energy_inBool:
+            return dims
+        energySym = Symbol(
+            "(energy)"
+        )  # create dummy symbol to replace in equation
+        dim_info = dims.as_terms()
+        time_idx = np.where(
+            list(map(lambda x: x.name == "(time)", dim_info[1]))
+        )[0][0]
+        energy_exp = (
+            dim_info[0][0][1][1][time_idx] // 2
+        )  # energy has 1/time**2 in it, so this is the hint of how many
+        return (
+            dims
+            * u.dimensions.energy**energy_exp
+            * energySym ** (-1 * energy_exp)
+        )
+
+    @staticmethod
+    def _dimensions_to_charge(dims):
+        """Take a set of dimensions and substitute in Symbol("charge") where possible."""
+        symsStr = str(dims.free_symbols)
+        charge_inBool = np.all(
+            [dimStr in symsStr for dimStr in ["current_mks"]]
+        )
+        if not charge_inBool:
+            return dims
+        chargeSym = Symbol(
+            "(charge)"
+        )  # create dummy symbol to replace in equation
+        dim_info = dims.as_terms()
+        current_idx = np.where(
+            list(map(lambda x: x.name == "(current_mks)", dim_info[1]))
+        )[0][0]
+        charge_exp = dim_info[0][0][1][1][
+            current_idx
+        ]  # charge has (current_mks) in it, so this is the hint of how many
+        return (
+            dims
+            * u.dimensions.charge ** (-1 * charge_exp)
+            * chargeSym**charge_exp
+        )
+
+    @staticmethod
+    def _dimensions_from_thermal_to_energy(dims):
+        """Take a set of dimensions and substitute in Symbol("energy") to replace temperature."""
+        symsStr = str(dims.free_symbols)
+        temp_inBool = np.all([dimStr in symsStr for dimStr in ["temperature"]])
+        if not temp_inBool:
+            return dims
+        energySym = Symbol(
+            "(energy)"
+        )  # create dummy symbol to replace in equation
+        dim_info = dims.as_terms()
+        temp_idx = np.where(
+            list(map(lambda x: x.name == "(temperature)", dim_info[1]))
+        )[0][0]
+        temp_exp = dim_info[0][0][1][1][
+            temp_idx
+        ]  # energy has 1/time**2 in it, so this is the hint of how many
+        return (
+            dims
+            / u.dimensions.temperature**temp_exp
+            * energySym ** (temp_exp)
+        )
+
+    @classmethod
+    def _get_output_dimensions(cls, dims, thermal_equivalence=False):
+        if str(dims) == "1":  # use string as all dims can be converted
+            return u.dimensionless
+        dims = cls._dimensions_to_energy(dims)
+        dims = cls._dimensions_to_charge(dims)
+        if thermal_equivalence:
+            dims = cls._dimensions_from_thermal_to_energy(dims)
+        return dims
+
+
+class ljUnitSystem:
+    """Use this so the empty unitsystem has getitem magic method."""
+
+    def __init__(self, reg: u.UnitRegistry):
+        self.registry = reg
+        self.name = "lj"
+
+    def __getitem__(self, item):
+        """Return dimensionless units unless angle."""
+        if item == "angle":
+            return u.Unit("degree")
+        return u.Unit("dimensionless")
+
+
+def get_parameter_dimension(parameter, dimension):
+    """Return a unit from the parameter in a given dimension."""
+    param_terms = parameter.units.expr.as_terms()
+    uStr = ""
+    for symbol, exp in zip(param_terms[-1], param_terms[0][0][1][1]):
+        outputDim = LAMMPS_UnitSystems._get_output_dimensions(
+            u.Unit(symbol).dimensions
+        )
+        if str(outputDim) == dimension:
+            uStr += f"{symbol}*"
+        elif (
+            str(outputDim) == "dimensionless" and dimension == "(energy)"
+        ):  # add mol to units of energy
+            uStr += f"{symbol}**{exp}*"
+        elif (
+            str(outputDim) == "dimensionless" and dimension == "(mass)"
+        ):  # add mol to mass amu
+            uStr += f"{symbol}**{exp}*"
+    return u.Unit(uStr[:-1])
+
+
+def write_out_parameter_and_units(parameter_name, parameter, base_unyts=None):
+    """Take a parameter and return a heading for the parameter and the units it should be in.
+
+    Parameters
+    ----------
+    parameter_name : str
+        The name of the unyt parameter to  be written. The dict key of the
+        parameter associated with the GMSO object.
+    parameter : unyt.array.unyt_quantity
+        The unyt object with the units to be written out.
+    base_unyts : LAMMPS_UnitSystem or a more general GMSO_UnitSystem
+        The base units that house the relevant registry for
+        converting parameters into a specified system.
+
+
+    Returns
+    -------
+    output_parameter_units : str
+        parameter with name converted into output unyt system. Useful for
+        labeling parameters in output files, such as .data or .top files.
+    """
+    if not base_unyts:
+        return f"{parameter_name} ({parameter.units})"
+    if parameter_name in ["theta_eq", "phi_eq"]:
+        return f"{parameter_name} ({'degrees'})"  # default to always degrees
+    if base_unyts.usystem.name == "lj":
+        return f"{parameter_name} ({'dimensionless'})"
+    new_dims = LAMMPS_UnitSystems._get_output_dimensions(
+        parameter.units.dimensions
+    )
+    new_dimStr = str(new_dims)
+    ind_units = re.sub("[^a-zA-Z]+", " ", new_dimStr).split()
+    for unit in ind_units:
+        new_dimStr = new_dimStr.replace(unit, str(base_unyts.usystem[unit]))
+
+    outputUnyt = str(
+        parameter.to(u.Unit(new_dimStr, registry=base_unyts.reg)).units
+    )
+    return f"{parameter_name} ({outputUnyt})"