From 172ba2ba90981b0b060c5087123349b05b186263 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Mon, 11 Sep 2023 15:47:06 -0600
Subject: [PATCH 01/41] create classmethods for two ways of initialization

---
 hoomd_organics/base/simulation.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 3e5c7ad2..4aae2fb1 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -52,6 +52,7 @@ def __init__(
         self,
         initial_state,
         forcefield=None,
+        reference_values=None,
         r_cut=2.5,
         dt=0.0001,
         device=hoomd.device.auto_select(),
@@ -83,12 +84,31 @@ def __init__(
         self.integrator = None
         self._dt = dt
         self._reference_values = dict()
+        if reference_values is not None:
+            self._reference_values = reference_values
         self._integrate_group = hoomd.filter.All()
         self._wall_forces = dict()
         self._create_state(self.initial_state)
         # Add a gsd and thermo props logger to sim operations
         self._add_hoomd_writers()
 
+    @classmethod
+    def from_system(cls, system, **kwargs):
+        """Initialize a simulation from a system object."""
+
+        return cls(
+            initial_state=system.state,
+            forcefield=system.hoomd_forcefield,
+            reference_values=system.reference_values,
+            **kwargs,
+        )
+
+    @classmethod
+    def from_init_state_forces(cls, init_state, forcefield, **kwargs):
+        """Initialize a simulation from an initial state object and a
+        list of HOOMD forces."""
+        return cls(initial_state=init_state, forcefield=forcefield, **kwargs)
+
     @property
     def forces(self):
         if self.integrator:

From d95de7969607147b2387e5c677a843f06a034b0f Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Mon, 11 Sep 2023 16:04:38 -0600
Subject: [PATCH 02/41] add density to run_update_volume

---
 hoomd_organics/base/simulation.py | 31 +++++++++++++++++++++----------
 1 file changed, 21 insertions(+), 10 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 4aae2fb1..c4b51964 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -52,7 +52,7 @@ def __init__(
         self,
         initial_state,
         forcefield=None,
-        reference_values=None,
+        reference_values=dict(),
         r_cut=2.5,
         dt=0.0001,
         device=hoomd.device.auto_select(),
@@ -84,8 +84,7 @@ def __init__(
         self.integrator = None
         self._dt = dt
         self._reference_values = dict()
-        if reference_values is not None:
-            self._reference_values = reference_values
+        self._reference_values = reference_values
         self._integrate_group = hoomd.filter.All()
         self._wall_forces = dict()
         self._create_state(self.initial_state)
@@ -399,7 +398,8 @@ def run_update_volume(
         period,
         kT,
         tau_kt,
-        final_box_lengths,
+        final_box_lengths=None,
+        final_density=None,
         thermalize_particles=True,
     ):
         """Runs an NVT simulation while shrinking or expanding
@@ -415,20 +415,31 @@ def run_update_volume(
             The temperature to use during shrinking.
         tau_kt : float; required
             Thermostat coupling period (in simulation time units)
-        final_box_lengths : np.ndarray, shape=(3,), dtype=float; required
+        final_box_lengths : np.ndarray, shape=(3,), dtype=float; optional
             The final box edge lengths in (x, y, z) order
+        final_density : float; optional
+            The final density of the simulation
 
         """
+        if final_box_lengths is None and final_density is None:
+            raise ValueError(
+                "Must provide either final_box_lengths or final_density"
+            )
+        if final_box_lengths:
+            final_box = hoomd.Box(
+                Lx=final_box_lengths[0],
+                Ly=final_box_lengths[1],
+                Lz=final_box_lengths[2],
+            )
+        else:
+            pass
+
         resize_trigger = hoomd.trigger.Periodic(period)
         box_ramp = hoomd.variant.Ramp(
             A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps)
         )
         initial_box = self.state.box
-        final_box = hoomd.Box(
-            Lx=final_box_lengths[0],
-            Ly=final_box_lengths[1],
-            Lz=final_box_lengths[2],
-        )
+
         box_resizer = hoomd.update.BoxResize(
             box1=initial_box,
             box2=final_box,

From eafe4a1018691d89e46f9712503247c281a28a42 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Mon, 11 Sep 2023 16:05:43 -0600
Subject: [PATCH 03/41] add target_box property and scale it if necessary

---
 hoomd_organics/base/system.py | 25 ++++++++++++++++++++++---
 1 file changed, 22 insertions(+), 3 deletions(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 24492227..41ce9f90 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -62,7 +62,7 @@ def __init__(
         self.remove_hydrogens = remove_hydrogens
         self.remove_charges = remove_charges
         self.scale_charges = scale_charges
-        self.target_box = None
+        self._target_box = None
         self.all_molecules = []
         self._hoomd_snapshot = None
         self._hoomd_forcefield = []
@@ -208,6 +208,10 @@ def hoomd_forcefield(self):
         self._hoomd_forcefield = self._create_hoomd_forcefield()
         return self._hoomd_forcefield
 
+    @property
+    def target_box(self):
+        return self._target_box / self.reference_length
+
     def _remove_hydrogens(self):
         """Call this method to remove hydrogen atoms from the system.
         The masses and charges of the hydrogens are absorbed into
@@ -344,7 +348,7 @@ def set_target_box(
             constraints[np.where(constraints is None)] = L
             Lx, Ly, Lz = constraints
 
-        self.target_box = np.array([Lx, Ly, Lz])
+        self._target_box = np.array([Lx, Ly, Lz])
 
     def visualize(self):
         if self.system:
@@ -389,6 +393,21 @@ class Pack(System):
     The box used for packing is expanded to allow PACKMOL
     to more easily place all the molecules.
 
+    Warnings:
+    ---------
+    Note that the default `packing_expand_factor` for pack is 5, which means
+    that the box density will not be the same as the specified density. This is
+    because in some cases PACKMOL will not be able to fit all the molecules
+    into the box if the target box is too small, therefore, we need to expand
+    the box by a factor (default:5) to allow PACKMOL to fit all the molecules.
+
+    In order to get the
+    specified density there are two options:
+    1) set the `packing_expand_factor` to 1, which will not expand the box,
+     however, this may result in PACKMOL errors if the box is too small.
+    2) Update the box volume after creating the simulation object to the target
+    box length. This property is called `target_box`.
+
     Parameters
     ----------
     packing_expand_factor : int; optional, default 5
@@ -428,7 +447,7 @@ def _build_system(self):
         system = mb.packing.fill_box(
             compound=self.all_molecules,
             n_compounds=[1 for i in self.all_molecules],
-            box=list(self.target_box * self.packing_expand_factor),
+            box=list(self._target_box * self.packing_expand_factor),
             overlap=0.2,
             edge=self.edge,
         )

From 607d4fbdd7f32da434d315dc8e33ad3e61bcd1d4 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Mon, 11 Sep 2023 16:06:14 -0600
Subject: [PATCH 04/41] calculate box length given mass and density

---
 hoomd_organics/utils/utils.py | 39 +++++++++++++++++++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/hoomd_organics/utils/utils.py b/hoomd_organics/utils/utils.py
index ead926e0..fa65be37 100644
--- a/hoomd_organics/utils/utils.py
+++ b/hoomd_organics/utils/utils.py
@@ -1,3 +1,4 @@
+import numpy as np
 import unyt as u
 
 from hoomd_organics.utils.exceptions import ReferenceUnitError
@@ -89,3 +90,41 @@ def _is_float(num):
             f"a reference value with unit of "
             f"{dimension} dimension."
         )
+
+
+def calculate_box_length(mass, density, fixed_L=None):
+    """Calculates the required box length(s) given the
+    mass of a sytem and the target density.
+
+    Box edge length constraints can be set by set_target_box().
+    If constraints are set, this will solve for the required
+    lengths of the remaining non-constrained edges to match
+    the target density.
+
+    Parameters
+    ----------
+    mass : float, required
+        Mass of the system
+    density : float, required
+        Target density of the system
+    fixed_L : np.array, optional, defualt=None
+        Array of fixed box lengths to be accounted for
+        when solving for L
+
+    Returns
+    -------
+    L : float
+        Box edge length
+    """
+    # Convert from amu to grams
+    M = mass * 1.66054e-24
+    vol = M / density  # cm^3
+    if fixed_L is None:
+        L = vol ** (1 / 3)
+    else:
+        L = vol / np.prod(fixed_L)
+        if len(fixed_L) == 1:  # L is cm^2
+            L = L ** (1 / 2)
+    # Convert from cm back to nm
+    L *= 1e7
+    return L

From e9c4d55ba2132e514bdebe4c26a46d1e7ca9a4b8 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Mon, 11 Sep 2023 16:14:16 -0600
Subject: [PATCH 05/41] call calculate_box_length

---
 hoomd_organics/base/simulation.py |  6 ++++--
 hoomd_organics/base/system.py     | 16 +++-------------
 hoomd_organics/utils/__init__.py  |  6 +++++-
 3 files changed, 12 insertions(+), 16 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index c4b51964..d90622b3 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -7,7 +7,7 @@
 import numpy as np
 import unyt as u
 
-from hoomd_organics.utils import StdOutLogger, UpdateWalls
+from hoomd_organics.utils import StdOutLogger, UpdateWalls, calculate_box_length
 from hoomd_organics.utils.exceptions import ReferenceUnitError
 
 
@@ -432,7 +432,9 @@ def run_update_volume(
                 Lz=final_box_lengths[2],
             )
         else:
-            pass
+            # TODO: should we convert the final_density to g/cm^3?
+            L = calculate_box_length(self.mass, final_density)
+            final_box = hoomd.Box(Lx=L, Ly=L, Lz=L)
 
         resize_trigger = hoomd.trigger.Periodic(period)
         box_ramp = hoomd.variant.Ramp(
diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 41ce9f90..9aa20bb6 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -18,6 +18,7 @@
 from hoomd_organics.base.molecule import Molecule
 from hoomd_organics.utils import (
     FF_Types,
+    calculate_box_length,
     check_return_iterable,
     validate_ref_value,
     xml_to_gmso_ff,
@@ -360,7 +361,7 @@ def visualize(self):
 
     def _calculate_L(self, fixed_L=None):
         """Calculates the required box length(s) given the
-        mass of a sytem and the target density.
+        mass of a system and the target density.
 
         Box edge length constraints can be set by set_target_box().
         If constraints are set, this will solve for the required
@@ -374,18 +375,7 @@ def _calculate_L(self, fixed_L=None):
             when solving for L
 
         """
-        # Convert from amu to grams
-        M = self.mass * 1.66054e-24
-        vol = M / self.density  # cm^3
-        if fixed_L is None:
-            L = vol ** (1 / 3)
-        else:
-            L = vol / np.prod(fixed_L)
-            if len(fixed_L) == 1:  # L is cm^2
-                L = L ** (1 / 2)
-        # Convert from cm back to nm
-        L *= 1e7
-        return L
+        return calculate_box_length(self.mass, self.density, fixed_L=fixed_L)
 
 
 class Pack(System):
diff --git a/hoomd_organics/utils/__init__.py b/hoomd_organics/utils/__init__.py
index 116f22ac..9558082f 100644
--- a/hoomd_organics/utils/__init__.py
+++ b/hoomd_organics/utils/__init__.py
@@ -1,4 +1,8 @@
 from .actions import *
 from .base_types import FF_Types
 from .ff_utils import xml_to_gmso_ff
-from .utils import check_return_iterable, validate_ref_value
+from .utils import (
+    calculate_box_length,
+    check_return_iterable,
+    validate_ref_value,
+)

From 4567dfc76026397e2eeea24e02524fe04ebce411 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Mon, 11 Sep 2023 16:18:18 -0600
Subject: [PATCH 06/41] get reference_length value

---
 hoomd_organics/base/system.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 9aa20bb6..0c1a33ba 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -211,7 +211,7 @@ def hoomd_forcefield(self):
 
     @property
     def target_box(self):
-        return self._target_box / self.reference_length
+        return self._target_box / self.reference_length.value
 
     def _remove_hydrogens(self):
         """Call this method to remove hydrogen atoms from the system.

From 4d73b12c506b725f0054b7042e6ce9953858f727 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 09:40:58 -0600
Subject: [PATCH 07/41] handle cases where system has no forcefield

---
 hoomd_organics/base/simulation.py | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index d90622b3..6cb261f5 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -95,15 +95,24 @@ def __init__(
     def from_system(cls, system, **kwargs):
         """Initialize a simulation from a system object."""
 
+        if system.hoomd_forcefield:
+            forcefield = system.hoomd_forcefield
+        elif kwargs["forcefield"]:
+            forcefield = kwargs["forcefield"]
+        else:
+            raise ValueError(
+                "No forcefield provided. Please provide a forcefield "
+                "or a system with a forcefield."
+            )
         return cls(
-            initial_state=system.state,
-            forcefield=system.hoomd_forcefield,
+            initial_state=system.hoomd_snapshot,
+            forcefield=forcefield,
             reference_values=system.reference_values,
             **kwargs,
         )
 
     @classmethod
-    def from_init_state_forces(cls, init_state, forcefield, **kwargs):
+    def from_init_state(cls, init_state, forcefield, **kwargs):
         """Initialize a simulation from an initial state object and a
         list of HOOMD forces."""
         return cls(initial_state=init_state, forcefield=forcefield, **kwargs)

From 2e466692ce0a010126faaff09fe8a74c36c019b9 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 10:27:21 -0600
Subject: [PATCH 08/41] unit tests for new sim interface

---
 hoomd_organics/base/simulation.py            |   6 +-
 hoomd_organics/tests/base/test_simulation.py | 101 +++++--------------
 2 files changed, 31 insertions(+), 76 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 6cb261f5..e5788fbf 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -112,10 +112,10 @@ def from_system(cls, system, **kwargs):
         )
 
     @classmethod
-    def from_init_state(cls, init_state, forcefield, **kwargs):
+    def from_state(cls, initial_state, forcefield, **kwargs):
         """Initialize a simulation from an initial state object and a
         list of HOOMD forces."""
-        return cls(initial_state=init_state, forcefield=forcefield, **kwargs)
+        return cls(initial_state=initial_state, forcefield=forcefield, **kwargs)
 
     @property
     def forces(self):
@@ -434,7 +434,7 @@ def run_update_volume(
             raise ValueError(
                 "Must provide either final_box_lengths or final_density"
             )
-        if final_box_lengths:
+        if final_box_lengths is not None:
             final_box = hoomd.Box(
                 Lx=final_box_lengths[0],
                 Ly=final_box_lengths[1],
diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 8b37be39..9f624161 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -9,14 +9,20 @@
 
 
 class TestSimulate(BaseTest):
-    def test_initialize_from_snap(self, benzene_system):
-        Simulation(
+    def test_initialize_from_system(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        assert len(sim.forces) == len(benzene_system.hoomd_forcefield)
+        assert sim.reference_values == benzene_system.reference_values
+
+    def test_initialize_from_state(self, benzene_system):
+        Simulation.from_state(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
+            reference_values=benzene_system.reference_values,
         )
 
     def test_no_reference_values(self, benzene_system):
-        sim = Simulation(
+        sim = Simulation.from_state(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
         )
@@ -29,24 +35,18 @@ def test_reference_values(self, benzene_system):
         sim = Simulation(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
+            reference_values=benzene_system.reference_values,
         )
-        sim.reference_values = benzene_system.reference_values
         assert np.isclose(float(sim.mass.value), benzene_system.mass, atol=1e-4)
         assert np.allclose(benzene_system.box.lengths, sim.box_lengths.value)
 
     def test_NVT(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500)
         assert isinstance(sim.method, hoomd.md.methods.NVT)
 
     def test_NPT(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_NPT(
             kT=1.0,
             n_steps=500,
@@ -57,34 +57,22 @@ def test_NPT(self, benzene_system):
         assert isinstance(sim.method, hoomd.md.methods.NPT)
 
     def test_langevin(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_langevin(n_steps=500, kT=1.0, alpha=0.5)
         assert isinstance(sim.method, hoomd.md.methods.Langevin)
 
     def test_NVE(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_NVE(n_steps=500)
         assert isinstance(sim.method, hoomd.md.methods.NVE)
 
     def test_displacement_cap(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_displacement_cap(n_steps=500, maximum_displacement=1e-4)
         assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped)
 
     def test_update_volume(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_update_volume(
             kT=1.0,
             tau_kt=0.01,
@@ -94,10 +82,7 @@ def test_update_volume(self, benzene_system):
         )
 
     def test_update_volume_walls(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.add_walls(wall_axis=(1, 0, 0), sigma=1.0, epsilon=1.0, r_cut=1.12)
         sim.run_update_volume(
             kT=1.0,
@@ -108,10 +93,7 @@ def test_update_volume_walls(self, benzene_system):
         )
 
     def test_change_methods(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)
         assert isinstance(sim.method, hoomd.md.methods.NVT)
         sim.run_NPT(
@@ -120,20 +102,14 @@ def test_change_methods(self, benzene_system):
         assert isinstance(sim.method, hoomd.md.methods.NPT)
 
     def test_change_dt(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)
         sim.dt = 0.003
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)
         assert sim.dt == 0.003
 
     def test_scale_epsilon(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         epsilons = []
         for param in sim._lj_force().params:
             epsilons.append(sim._lj_force().params[param]["epsilon"])
@@ -145,10 +121,7 @@ def test_scale_epsilon(self, benzene_system):
             assert np.allclose(i * 0.5, j, atol=1e-3)
 
     def test_shift_epsilon(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         epsilons = []
         for param in sim._lj_force().params:
             epsilons.append(sim._lj_force().params[param]["epsilon"])
@@ -160,10 +133,7 @@ def test_shift_epsilon(self, benzene_system):
             assert np.allclose(i + 1, j, atol=1e-3)
 
     def test_scale_sigma(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sigmas = []
         for param in sim._lj_force().params:
             sigmas.append(sim._lj_force().params[param]["sigma"])
@@ -175,10 +145,7 @@ def test_scale_sigma(self, benzene_system):
             assert np.allclose(i * 0.5, j, atol=1e-3)
 
     def test_shift_sigma(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sigmas = []
         for param in sim._lj_force().params:
             sigmas.append(sim._lj_force().params[param]["sigma"])
@@ -190,19 +157,13 @@ def test_shift_sigma(self, benzene_system):
             assert np.allclose(i + 1, j, atol=1e-3)
 
     def test_remove_force(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.remove_force(sim._lj_force())
         for i in sim.forces:
             assert not isinstance(i, hoomd.md.pair.LJ)
 
     def test_set_integrate_group(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         assert isinstance(sim.integrate_group, hoomd.filter.All)
         tag_filter = hoomd.filter.Tags([0, 1, 2, 3])
         sim.integrate_group = tag_filter
@@ -210,10 +171,7 @@ def test_set_integrate_group(self, benzene_system):
         sim.run_NVT(n_steps=200, kT=1.0, tau_kt=0.01)
 
     def test_pickle_ff(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.pickle_forcefield("forcefield.pickle")
         assert os.path.isfile("forcefield.pickle")
         f = open("forcefield.pickle", "rb")
@@ -224,15 +182,12 @@ def test_pickle_ff(self, benzene_system):
         os.remove("forcefield.pickle")
 
     def test_save_restart_gsd(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
+        sim = Simulation.from_system(benzene_system)
         sim.save_restart_gsd("restart.gsd")
         assert os.path.isfile("restart.gsd")
         sim.pickle_forcefield("forcefield.pickle")
         f = open("forcefield.pickle", "rb")
         hoomd_ff = pickle.load(f)
-        Simulation(initial_state="restart.gsd", forcefield=hoomd_ff)
+        Simulation.from_state(initial_state="restart.gsd", forcefield=hoomd_ff)
         os.remove("forcefield.pickle")
         os.remove("restart.gsd")

From 5873afb42204653f8d6d97fc42ef71cbdb78ee4b Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 10:27:36 -0600
Subject: [PATCH 09/41] add a cg fixture

---
 hoomd_organics/tests/base_test.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py
index a7f1c241..7f21ea01 100644
--- a/hoomd_organics/tests/base_test.py
+++ b/hoomd_organics/tests/base_test.py
@@ -121,6 +121,11 @@ def _benzene_molecule(n_mols):
 
         return _benzene_molecule
 
+    @pytest.fixture()
+    def benzene_cg(self, benzene_molecule, benzene_smiles):
+        benzene_molecule.coarse_grain(beads={"A": benzene_smiles})
+        return benzene_molecule
+
     @pytest.fixture()
     def ethane_molecule(self, ethane_smiles):
         def _ethane_molecule(n_mols):

From 0a01f1a05c8ba8279a853d13c8fc5aac45e61a3c Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:10:01 -0600
Subject: [PATCH 10/41] fix bug in kwargs[forcefield]

---
 hoomd_organics/base/simulation.py | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index e5788fbf..d8c83a1b 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -96,20 +96,23 @@ def from_system(cls, system, **kwargs):
         """Initialize a simulation from a system object."""
 
         if system.hoomd_forcefield:
-            forcefield = system.hoomd_forcefield
+            return cls(
+                initial_state=system.hoomd_snapshot,
+                forcefield=system.hoomd_forcefield,
+                reference_values=system.reference_values,
+                **kwargs,
+            )
         elif kwargs["forcefield"]:
-            forcefield = kwargs["forcefield"]
+            return cls(
+                initial_state=system.hoomd_snapshot,
+                reference_values=system.reference_values,
+                **kwargs,
+            )
         else:
             raise ValueError(
                 "No forcefield provided. Please provide a forcefield "
                 "or a system with a forcefield."
             )
-        return cls(
-            initial_state=system.hoomd_snapshot,
-            forcefield=forcefield,
-            reference_values=system.reference_values,
-            **kwargs,
-        )
 
     @classmethod
     def from_state(cls, initial_state, forcefield, **kwargs):

From 294e3ef627b67d6a50bac36d655af67fe0d5504f Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:10:22 -0600
Subject: [PATCH 11/41] change mol name when ff is provided

---
 hoomd_organics/base/system.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 0c1a33ba..8e8bc88a 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -75,7 +75,8 @@ def __init__(
         self.n_mol_types = 0
         for mol_item in self._molecules:
             if isinstance(mol_item, Molecule):
-                mol_item.assign_mol_name(str(self.n_mol_types))
+                if self._force_field:
+                    mol_item.assign_mol_name(str(self.n_mol_types))
                 self.all_molecules.extend(mol_item.molecules)
                 # if ff is provided in Molecule class
                 if mol_item.force_field:
@@ -206,8 +207,11 @@ def hoomd_snapshot(self):
 
     @property
     def hoomd_forcefield(self):
-        self._hoomd_forcefield = self._create_hoomd_forcefield()
-        return self._hoomd_forcefield
+        if self._force_field:
+            self._hoomd_forcefield = self._create_hoomd_forcefield()
+            return self._hoomd_forcefield
+        else:
+            return self._hoomd_forcefield
 
     @property
     def target_box(self):

From 330be4f2dc4bd9f6388939720811d08ab0944fd2 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:10:45 -0600
Subject: [PATCH 12/41] cg system fixture

---
 hoomd_organics/tests/base_test.py | 29 +++++++++++++++++++++++------
 1 file changed, 23 insertions(+), 6 deletions(-)

diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py
index 7f21ea01..9642bd08 100644
--- a/hoomd_organics/tests/base_test.py
+++ b/hoomd_organics/tests/base_test.py
@@ -6,7 +6,7 @@
 from gmso.external.convert_mbuild import from_mbuild
 
 from hoomd_organics import Molecule, Pack, Polymer, Simulation
-from hoomd_organics.library import OPLS_AA
+from hoomd_organics.library import OPLS_AA, BeadSpring
 
 ASSETS_DIR = os.path.join(os.path.dirname(__file__), "assets")
 
@@ -121,11 +121,6 @@ def _benzene_molecule(n_mols):
 
         return _benzene_molecule
 
-    @pytest.fixture()
-    def benzene_cg(self, benzene_molecule, benzene_smiles):
-        benzene_molecule.coarse_grain(beads={"A": benzene_smiles})
-        return benzene_molecule
-
     @pytest.fixture()
     def ethane_molecule(self, ethane_smiles):
         def _ethane_molecule(n_mols):
@@ -224,6 +219,18 @@ def benzene_system(self, benzene_mb):
         )
         return system
 
+    @pytest.fixture()
+    def benzene_cg_system(self, benzene_molecule, benzene_smiles):
+        benzene_mols = benzene_molecule(n_mols=10)
+        benzene_mols.coarse_grain(beads={"A": benzene_smiles})
+        system = Pack(
+            molecules=[benzene_mols],
+            density=0.01,
+            r_cut=2.5,
+            auto_scale=False,
+        )
+        return system
+
     @pytest.fixture()
     def polyethylene_system(self, polyethylene):
         polyethylene_mol = polyethylene(num_mols=5, lengths=5)
@@ -244,3 +251,13 @@ def benzene_simulation(self, benzene_system):
             forcefield=benzene_system.hoomd_forcefield,
         )
         return sim
+
+    @pytest.fixture()
+    def cg_single_bead_ff(self):
+        ff = BeadSpring(
+            r_cut=2.5,
+            beads={
+                "A": dict(epsilon=1.0, sigma=1.0),
+            },
+        )
+        return ff.hoomd_forcefield

From df941e1e0ea90a85e9c37d17547db2826f50e548 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:11:16 -0600
Subject: [PATCH 13/41] add cg system test

---
 hoomd_organics/tests/base/test_simulation.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 9f624161..4f594f46 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -11,9 +11,18 @@
 class TestSimulate(BaseTest):
     def test_initialize_from_system(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
+        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500)
         assert len(sim.forces) == len(benzene_system.hoomd_forcefield)
         assert sim.reference_values == benzene_system.reference_values
 
+    def test_initialize_from_system_no_ff(
+        self, benzene_cg_system, cg_single_bead_ff
+    ):
+        sim = Simulation.from_system(
+            benzene_cg_system, forcefield=cg_single_bead_ff
+        )
+        sim.run_NVT(kT=0.1, tau_kt=10, n_steps=500)
+
     def test_initialize_from_state(self, benzene_system):
         Simulation.from_state(
             initial_state=benzene_system.hoomd_snapshot,

From 1f03e36f00f411627f3bfed6075438df70de1753 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:14:06 -0600
Subject: [PATCH 14/41] update docstring, raise error if both box length and
 density are provided

---
 hoomd_organics/base/simulation.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index d8c83a1b..448fd225 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -93,7 +93,8 @@ def __init__(
 
     @classmethod
     def from_system(cls, system, **kwargs):
-        """Initialize a simulation from a system object."""
+        """Initialize a simulation from a `hoomd_organics.base.System`
+        object."""
 
         if system.hoomd_forcefield:
             return cls(
@@ -435,7 +436,11 @@ def run_update_volume(
         """
         if final_box_lengths is None and final_density is None:
             raise ValueError(
-                "Must provide either final_box_lengths or final_density"
+                "Must provide either `final_box_lengths` or `final_density`"
+            )
+        if final_box_lengths is not None and final_density is not None:
+            raise ValueError(
+                "Cannot provide both `final_box_lengths` and `final_density`."
             )
         if final_box_lengths is not None:
             final_box = hoomd.Box(
@@ -687,7 +692,7 @@ def _add_hoomd_writers(self):
         )
 
         logger = hoomd.logging.Logger(categories=["scalar", "string"])
-        logger.add(self, quantities=["timestep", "tps"])
+        logger.add(self, quantitietes=["timestep", "tps"])
         thermo_props = hoomd.md.compute.ThermodynamicQuantities(
             filter=self.integrate_group
         )

From 1c7f9258b83d0eb77f30fa83dfc667e95892ae9c Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:20:54 -0600
Subject: [PATCH 15/41] add unit test for sim init without a ff

---
 hoomd_organics/base/simulation.py            | 2 +-
 hoomd_organics/tests/base/test_simulation.py | 7 ++++++-
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 448fd225..dd8b5acd 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -103,7 +103,7 @@ def from_system(cls, system, **kwargs):
                 reference_values=system.reference_values,
                 **kwargs,
             )
-        elif kwargs["forcefield"]:
+        elif kwargs.get("forcefield", None):
             return cls(
                 initial_state=system.hoomd_snapshot,
                 reference_values=system.reference_values,
diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 4f594f46..88a86e89 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -3,6 +3,7 @@
 
 import hoomd
 import numpy as np
+import pytest
 
 from hoomd_organics import Simulation
 from hoomd_organics.tests import BaseTest
@@ -15,7 +16,7 @@ def test_initialize_from_system(self, benzene_system):
         assert len(sim.forces) == len(benzene_system.hoomd_forcefield)
         assert sim.reference_values == benzene_system.reference_values
 
-    def test_initialize_from_system_no_ff(
+    def test_initialize_from_system_separate_ff(
         self, benzene_cg_system, cg_single_bead_ff
     ):
         sim = Simulation.from_system(
@@ -23,6 +24,10 @@ def test_initialize_from_system_no_ff(
         )
         sim.run_NVT(kT=0.1, tau_kt=10, n_steps=500)
 
+    def test_initialize_from_system_missing_ff(self, benzene_cg_system):
+        with pytest.raises(ValueError):
+            Simulation.from_system(benzene_cg_system)
+
     def test_initialize_from_state(self, benzene_system):
         Simulation.from_state(
             initial_state=benzene_system.hoomd_snapshot,

From b94f85a413ea5db36523a386e04fbd7e6c38f4a5 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:40:21 -0600
Subject: [PATCH 16/41] update apply parameters

---
 hoomd_organics/base/simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index dd8b5acd..9c565d81 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -692,7 +692,7 @@ def _add_hoomd_writers(self):
         )
 
         logger = hoomd.logging.Logger(categories=["scalar", "string"])
-        logger.add(self, quantitietes=["timestep", "tps"])
+        logger.add(self, quantities=["timestep", "tps"])
         thermo_props = hoomd.md.compute.ThermodynamicQuantities(
             filter=self.integrate_group
         )

From 7f11cec562b118f838e75fa151c3fbeb4c637b27 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:41:40 -0600
Subject: [PATCH 17/41] fix cg sim test

---
 hoomd_organics/base/system.py     | 11 +++--------
 hoomd_organics/tests/base_test.py |  2 +-
 2 files changed, 4 insertions(+), 9 deletions(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 8e8bc88a..1d562e14 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -294,8 +294,8 @@ def _apply_forcefield(self):
             self.gmso_system,
             self._gmso_forcefields_dict,
             identify_connections=True,
-            use_molecule_info=True,
-            identify_connected_components=False,
+            speedup_by_moltag=True,
+            speedup_by_molgraph=False,
         )
         if self.remove_charges:
             for site in self.gmso_system.sites:
@@ -316,12 +316,7 @@ def _apply_forcefield(self):
 
         self._reference_values["energy"] = energy_scale * epsilons[0].unit_array
         self._reference_values["length"] = length_scale * sigmas[0].unit_array
-        if self.auto_scale:
-            self._reference_values["mass"] = mass_scale * masses[
-                0
-            ].unit_array.to("amu")
-        else:
-            self._reference_values["mass"] = mass_scale * u.g / u.mol
+        self._reference_values["mass"] = mass_scale * masses[0].unit_array
 
     def set_target_box(
         self, x_constraint=None, y_constraint=None, z_constraint=None
diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py
index 9642bd08..9322e62a 100644
--- a/hoomd_organics/tests/base_test.py
+++ b/hoomd_organics/tests/base_test.py
@@ -225,7 +225,7 @@ def benzene_cg_system(self, benzene_molecule, benzene_smiles):
         benzene_mols.coarse_grain(beads={"A": benzene_smiles})
         system = Pack(
             molecules=[benzene_mols],
-            density=0.01,
+            density=0.5,
             r_cut=2.5,
             auto_scale=False,
         )

From bc69f94066ab117c0a507cd455c0397000d3141e Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 15:45:34 -0600
Subject: [PATCH 18/41] change class name to

---
 hoomd_organics/base/simulation.py            | 2 +-
 hoomd_organics/tests/base/test_simulation.py | 8 +++++---
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 9c565d81..63c0e695 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -116,7 +116,7 @@ def from_system(cls, system, **kwargs):
             )
 
     @classmethod
-    def from_state(cls, initial_state, forcefield, **kwargs):
+    def from_snapshot_forces(cls, initial_state, forcefield, **kwargs):
         """Initialize a simulation from an initial state object and a
         list of HOOMD forces."""
         return cls(initial_state=initial_state, forcefield=forcefield, **kwargs)
diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 88a86e89..0ac8e404 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -29,14 +29,14 @@ def test_initialize_from_system_missing_ff(self, benzene_cg_system):
             Simulation.from_system(benzene_cg_system)
 
     def test_initialize_from_state(self, benzene_system):
-        Simulation.from_state(
+        Simulation.from_snapshot_forces(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
             reference_values=benzene_system.reference_values,
         )
 
     def test_no_reference_values(self, benzene_system):
-        sim = Simulation.from_state(
+        sim = Simulation.from_snapshot_forces(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
         )
@@ -202,6 +202,8 @@ def test_save_restart_gsd(self, benzene_system):
         sim.pickle_forcefield("forcefield.pickle")
         f = open("forcefield.pickle", "rb")
         hoomd_ff = pickle.load(f)
-        Simulation.from_state(initial_state="restart.gsd", forcefield=hoomd_ff)
+        Simulation.from_snapshot_forces(
+            initial_state="restart.gsd", forcefield=hoomd_ff
+        )
         os.remove("forcefield.pickle")
         os.remove("restart.gsd")

From b6d7c7f6fe25206194be576d1bde60a99f267cff Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 16:27:18 -0600
Subject: [PATCH 19/41] raise error if update vol has density but no ref units

---
 hoomd_organics/base/simulation.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 63c0e695..df828028 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -449,8 +449,17 @@ def run_update_volume(
                 Lz=final_box_lengths[2],
             )
         else:
-            # TODO: should we convert the final_density to g/cm^3?
-            L = calculate_box_length(self.mass, final_density)
+            if self.reference_values is None:
+                raise ReferenceUnitError(
+                    "Missing simulation units. Please "
+                    "provide units for mass, length, and"
+                    " energy."
+                )
+
+            density_quantity = u.unyt_quantity(final_density, u.g / u.cm**3)
+            L = calculate_box_length(self.mass, density_quantity)
+            # convert L from cm to reference units
+            L = (L.to(self.reference_length.unit) * self.reference_length).value
             final_box = hoomd.Box(Lx=L, Ly=L, Lz=L)
 
         resize_trigger = hoomd.trigger.Periodic(period)

From ba7c59c68cce4f4d6ba0d7359aa62c33fa8e97c3 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 16:27:55 -0600
Subject: [PATCH 20/41] update calculate box length to get unyt_quantity

---
 hoomd_organics/base/system.py |  7 ++++++-
 hoomd_organics/utils/utils.py | 11 +++++------
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 1d562e14..67e5bf5f 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -374,7 +374,12 @@ def _calculate_L(self, fixed_L=None):
             when solving for L
 
         """
-        return calculate_box_length(self.mass, self.density, fixed_L=fixed_L)
+        mass_quantity = u.unyt_quantity(self.mass, u.g / u.mol)
+        density_quantity = u.unyt_quantity(self.density, u.g / u.cm**3)
+        L = calculate_box_length(
+            mass_quantity, density_quantity, fixed_L=fixed_L
+        )
+        return L.to("nm").value
 
 
 class Pack(System):
diff --git a/hoomd_organics/utils/utils.py b/hoomd_organics/utils/utils.py
index fa65be37..193289e9 100644
--- a/hoomd_organics/utils/utils.py
+++ b/hoomd_organics/utils/utils.py
@@ -103,9 +103,9 @@ def calculate_box_length(mass, density, fixed_L=None):
 
     Parameters
     ----------
-    mass : float, required
+    mass : unyt.unyt_quantity, required
         Mass of the system
-    density : float, required
+    density : unyt.unyt_quantity, required
         Target density of the system
     fixed_L : np.array, optional, defualt=None
         Array of fixed box lengths to be accounted for
@@ -116,8 +116,8 @@ def calculate_box_length(mass, density, fixed_L=None):
     L : float
         Box edge length
     """
-    # Convert from amu to grams
-    M = mass * 1.66054e-24
+    # Convert mass to grams
+    M = mass.to("g")
     vol = M / density  # cm^3
     if fixed_L is None:
         L = vol ** (1 / 3)
@@ -125,6 +125,5 @@ def calculate_box_length(mass, density, fixed_L=None):
         L = vol / np.prod(fixed_L)
         if len(fixed_L) == 1:  # L is cm^2
             L = L ** (1 / 2)
-    # Convert from cm back to nm
-    L *= 1e7
+    # L is cm
     return L

From f7b9bded6059f9be4451886f7a5a4977906c42ff Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 16:28:13 -0600
Subject: [PATCH 21/41] use reduced box length for shrink

---
 hoomd_organics/tests/base/test_simulation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 0ac8e404..28ad6141 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -92,7 +92,7 @@ def test_update_volume(self, benzene_system):
             tau_kt=0.01,
             n_steps=500,
             period=1,
-            final_box_lengths=sim.box_lengths * 0.5,
+            final_box_lengths=sim.box_lengths_reduced * 0.5,
         )
 
     def test_update_volume_walls(self, benzene_system):
@@ -103,7 +103,7 @@ def test_update_volume_walls(self, benzene_system):
             tau_kt=0.01,
             n_steps=500,
             period=5,
-            final_box_lengths=sim.box_lengths * 0.5,
+            final_box_lengths=sim.box_lengths_reduced * 0.5,
         )
 
     def test_change_methods(self, benzene_system):

From b81046292cc8656f26309e3336bd095f5f96663f Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 16:59:09 -0600
Subject: [PATCH 22/41] install gmso from source

---
 environment-cpu.yml | 2 +-
 environment-dev.yml | 2 +-
 environment-gpu.yml | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/environment-cpu.yml b/environment-cpu.yml
index cd9b7c77..e0253a24 100644
--- a/environment-cpu.yml
+++ b/environment-cpu.yml
@@ -4,7 +4,6 @@ channels:
 dependencies:
   - foyer
   - freud
-  - gmso>=0.11.1
   - gsd<3.0
   - hoomd=3.11=*cpu*
   - mbuild
@@ -19,3 +18,4 @@ dependencies:
   - pip:
       - git+https://github.com/cmelab/cmeutils.git@master
       - git+https://github.com/cmelab/grits.git@main
+      - git+https://github.com/mosdef-hub/gmso@main
diff --git a/environment-dev.yml b/environment-dev.yml
index ff931fd0..e984684f 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -4,7 +4,6 @@ channels:
 dependencies:
   - foyer
   - freud
-  - gmso>=0.11.1
   - gsd<3.0
   - hoomd=3.11=*cpu*
   - mbuild
@@ -20,3 +19,4 @@ dependencies:
   - pip:
       - git+https://github.com/cmelab/cmeutils.git@master
       - git+https://github.com/cmelab/grits.git@main
+      - git+https://github.com/mosdef-hub/gmso@main
diff --git a/environment-gpu.yml b/environment-gpu.yml
index c3e25a2b..b7168d7f 100644
--- a/environment-gpu.yml
+++ b/environment-gpu.yml
@@ -4,7 +4,6 @@ channels:
 dependencies:
   - foyer
   - freud
-  - gmso>=0.11.1
   - gsd<3.0
   - hoomd=3.11=*gpu*
   - mbuild
@@ -19,3 +18,4 @@ dependencies:
   - pip:
       - git+https://github.com/cmelab/cmeutils.git@master
       - git+https://github.com/cmelab/grits.git@main
+      - git+https://github.com/mosdef-hub/gmso@main

From cfa8602a6da1963b47d5c50424c9b3c57765a2fd Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 17:00:06 -0600
Subject: [PATCH 23/41] check if ref length is available for target box

---
 hoomd_organics/base/system.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 67e5bf5f..36d78a5a 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -215,7 +215,10 @@ def hoomd_forcefield(self):
 
     @property
     def target_box(self):
-        return self._target_box / self.reference_length.value
+        if self.reference_length:
+            return self._target_box / self.reference_length.value
+        else:
+            return self._target_box
 
     def _remove_hydrogens(self):
         """Call this method to remove hydrogen atoms from the system.

From 50e73171bb3347de6fe4f9256362ad92e2f4b80e Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Tue, 12 Sep 2023 17:00:36 -0600
Subject: [PATCH 24/41] update unit test name

---
 hoomd_organics/tests/base/test_simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 28ad6141..7bd5d686 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -85,7 +85,7 @@ def test_displacement_cap(self, benzene_system):
         sim.run_displacement_cap(n_steps=500, maximum_displacement=1e-4)
         assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped)
 
-    def test_update_volume(self, benzene_system):
+    def test_update_volume_target_box(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_update_volume(
             kT=1.0,

From 19706618f505391c3460a34baf679128771cf781 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Wed, 13 Sep 2023 12:52:58 -0600
Subject: [PATCH 25/41] add gmso dependencies to env

---
 environment-cpu.yml | 6 +++++-
 environment-dev.yml | 6 +++++-
 environment-gpu.yml | 6 +++++-
 3 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/environment-cpu.yml b/environment-cpu.yml
index e0253a24..1aa085d4 100644
--- a/environment-cpu.yml
+++ b/environment-cpu.yml
@@ -15,7 +15,11 @@ dependencies:
   - pytest-cov
   - python=3.9
   - unyt
+  - sympy
+  - boltons
+  - symengine
+  - python-symengine
   - pip:
       - git+https://github.com/cmelab/cmeutils.git@master
       - git+https://github.com/cmelab/grits.git@main
-      - git+https://github.com/mosdef-hub/gmso@main
+      - git+https://github.com/mosdef-hub/gmso.git@main
diff --git a/environment-dev.yml b/environment-dev.yml
index e984684f..1ea2cc39 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -16,7 +16,11 @@ dependencies:
   - python=3.9
   - unyt
   - pre-commit
+  - sympy
+  - boltons
+  - symengine
+  - python-symengine
   - pip:
       - git+https://github.com/cmelab/cmeutils.git@master
       - git+https://github.com/cmelab/grits.git@main
-      - git+https://github.com/mosdef-hub/gmso@main
+      - git+https://github.com/mosdef-hub/gmso.git@main
diff --git a/environment-gpu.yml b/environment-gpu.yml
index b7168d7f..54a42eef 100644
--- a/environment-gpu.yml
+++ b/environment-gpu.yml
@@ -15,7 +15,11 @@ dependencies:
   - pytest-cov
   - python=3.9
   - unyt
+  - sympy
+  - boltons
+  - symengine
+  - python-symengine
   - pip:
       - git+https://github.com/cmelab/cmeutils.git@master
       - git+https://github.com/cmelab/grits.git@main
-      - git+https://github.com/mosdef-hub/gmso@main
+      - git+https://github.com/mosdef-hub/gmso.git@main

From 38d23e39f1141756aa17b3affa7e55c3adf81af0 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Wed, 13 Sep 2023 14:52:32 -0600
Subject: [PATCH 26/41] unit test for update volume with density

---
 hoomd_organics/base/simulation.py            |  4 +++-
 hoomd_organics/tests/base/test_simulation.py | 12 ++++++++++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index df828028..82dccc7b 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -459,7 +459,9 @@ def run_update_volume(
             density_quantity = u.unyt_quantity(final_density, u.g / u.cm**3)
             L = calculate_box_length(self.mass, density_quantity)
             # convert L from cm to reference units
-            L = (L.to(self.reference_length.unit) * self.reference_length).value
+            L = (
+                L.to(self.reference_length.units) / self.reference_length.value
+            ).value
             final_box = hoomd.Box(Lx=L, Ly=L, Lz=L)
 
         resize_trigger = hoomd.trigger.Periodic(period)
diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 7bd5d686..548cf408 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -4,6 +4,7 @@
 import hoomd
 import numpy as np
 import pytest
+import unyt as u
 
 from hoomd_organics import Simulation
 from hoomd_organics.tests import BaseTest
@@ -106,6 +107,17 @@ def test_update_volume_walls(self, benzene_system):
             final_box_lengths=sim.box_lengths_reduced * 0.5,
         )
 
+    def test_update_volume_density(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_update_volume(
+            kT=1.0, tau_kt=0.01, n_steps=500, period=1, final_density=0.1
+        )
+        assert np.isclose(
+            sim.density.to(u.g / u.cm**3).value,
+            (0.1 * (u.g / u.cm**3)).value,
+            atol=1e-4,
+        )
+
     def test_change_methods(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)

From d684c30ca0515a96a76a3c4ce36eceb551af4ac1 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Wed, 13 Sep 2023 14:52:59 -0600
Subject: [PATCH 27/41] increase system size in fixture

---
 hoomd_organics/tests/base_test.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py
index 9322e62a..8a69d85b 100644
--- a/hoomd_organics/tests/base_test.py
+++ b/hoomd_organics/tests/base_test.py
@@ -209,10 +209,10 @@ def __init__(self, lengths, num_mols, **kwargs):
 
     @pytest.fixture()
     def benzene_system(self, benzene_mb):
-        benzene = Molecule(num_mols=5, compound=benzene_mb)
+        benzene = Molecule(num_mols=20, compound=benzene_mb)
         system = Pack(
             molecules=[benzene],
-            density=0.5,
+            density=0.2,
             r_cut=2.5,
             force_field=OPLS_AA(),
             auto_scale=True,

From 91c3c03acfea27623bf08bc06090c84b64234e51 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Wed, 13 Sep 2023 15:11:36 -0600
Subject: [PATCH 28/41] handle cases where density is unit_quantity

---
 hoomd_organics/base/simulation.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 82dccc7b..0a4cdcbd 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -456,7 +456,12 @@ def run_update_volume(
                     " energy."
                 )
 
-            density_quantity = u.unyt_quantity(final_density, u.g / u.cm**3)
+            if isinstance(final_density, u.unyt_quantity):
+                density_quantity = final_density.to(u.g / u.cm**3)
+            else:
+                density_quantity = u.unyt_quantity(
+                    final_density, u.g / u.cm**3
+                )
             L = calculate_box_length(self.mass, density_quantity)
             # convert L from cm to reference units
             L = (

From a10bd7324a161af50b1772c9ed990f80e2c28682 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Wed, 13 Sep 2023 15:12:11 -0600
Subject: [PATCH 29/41] unit test for update volume by density factor

---
 hoomd_organics/tests/base/sim_data.txt       |  0
 hoomd_organics/tests/base/test_simulation.py | 20 ++++++++++++++++++++
 2 files changed, 20 insertions(+)
 create mode 100644 hoomd_organics/tests/base/sim_data.txt

diff --git a/hoomd_organics/tests/base/sim_data.txt b/hoomd_organics/tests/base/sim_data.txt
new file mode 100644
index 00000000..e69de29b
diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 548cf408..97373c98 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -1,3 +1,4 @@
+import copy
 import os
 import pickle
 
@@ -118,6 +119,25 @@ def test_update_volume_density(self, benzene_system):
             atol=1e-4,
         )
 
+    def test_update_volume_by_density_factor(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        init_density = copy.deepcopy(sim.density)
+        sim.run_update_volume(
+            kT=1.0,
+            tau_kt=0.01,
+            n_steps=500,
+            period=1,
+            final_density=sim.density * 5,
+        )
+        assert np.isclose(
+            sim.density.value, (init_density * 5).value, atol=1e-4
+        )
+
+    def test_update_volume_missing_values(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        with pytest.raises(ValueError):
+            sim.run_update_volume(kT=1.0, tau_kt=0.01, n_steps=500, period=1)
+
     def test_change_methods(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)

From 5420e337df95e52d314352715f17603fe9c2ccf1 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 12:26:06 -0600
Subject: [PATCH 30/41] fix bug is finding constraints that aren't none

---
 hoomd_organics/base/system.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 36d78a5a..1029a94e 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -344,11 +344,11 @@ def set_target_box(
             Lx = Ly = Lz = self._calculate_L()
         else:
             constraints = np.array([x_constraint, y_constraint, z_constraint])
-            fixed_L = constraints[np.where(constraints is not None)]
+            fixed_L = constraints[np.not_equal(constraints, None).nonzero()]
             # Conv from nm to cm for _calculate_L
             fixed_L *= 1e-7
             L = self._calculate_L(fixed_L=fixed_L)
-            constraints[np.where(constraints is None)] = L
+            constraints[np.equal(constraints, None).nonzero()] = L
             Lx, Ly, Lz = constraints
 
         self._target_box = np.array([Lx, Ly, Lz])

From 1d7b8367368f9d3b5fa27aa804ab9e11e3d7827d Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 12:33:40 -0600
Subject: [PATCH 31/41] remove mass conversion to g from calculate_box_length

---
 hoomd_organics/base/simulation.py | 3 ++-
 hoomd_organics/base/system.py     | 2 +-
 hoomd_organics/utils/utils.py     | 5 ++---
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 0a4cdcbd..933acb43 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -462,7 +462,8 @@ def run_update_volume(
                 density_quantity = u.unyt_quantity(
                     final_density, u.g / u.cm**3
                 )
-            L = calculate_box_length(self.mass, density_quantity)
+            mass_g = self.mass.to("g")
+            L = calculate_box_length(mass_g, density_quantity)
             # convert L from cm to reference units
             L = (
                 L.to(self.reference_length.units) / self.reference_length.value
diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index 1029a94e..b60474bf 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -377,7 +377,7 @@ def _calculate_L(self, fixed_L=None):
             when solving for L
 
         """
-        mass_quantity = u.unyt_quantity(self.mass, u.g / u.mol)
+        mass_quantity = u.unyt_quantity(self.mass, u.g / u.mol).to("g")
         density_quantity = u.unyt_quantity(self.density, u.g / u.cm**3)
         L = calculate_box_length(
             mass_quantity, density_quantity, fixed_L=fixed_L
diff --git a/hoomd_organics/utils/utils.py b/hoomd_organics/utils/utils.py
index 193289e9..6536c15e 100644
--- a/hoomd_organics/utils/utils.py
+++ b/hoomd_organics/utils/utils.py
@@ -116,9 +116,8 @@ def calculate_box_length(mass, density, fixed_L=None):
     L : float
         Box edge length
     """
-    # Convert mass to grams
-    M = mass.to("g")
-    vol = M / density  # cm^3
+
+    vol = mass / density  # cm^3
     if fixed_L is None:
         L = vol ** (1 / 3)
     else:

From 76bc42131b04dccd6356063711dcda4262c25173 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 12:39:08 -0600
Subject: [PATCH 32/41] make fixed_L variable a unyt_array

---
 hoomd_organics/base/system.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index b60474bf..cbec2fbc 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -379,6 +379,7 @@ def _calculate_L(self, fixed_L=None):
         """
         mass_quantity = u.unyt_quantity(self.mass, u.g / u.mol).to("g")
         density_quantity = u.unyt_quantity(self.density, u.g / u.cm**3)
+        fixed_L = u.unyt_array(fixed_L, u.cm)
         L = calculate_box_length(
             mass_quantity, density_quantity, fixed_L=fixed_L
         )

From d30d9b8c1fe7a0050c87bd42852f464ce06b85bb Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 12:49:08 -0600
Subject: [PATCH 33/41] unit tests for calculate_box_length

---
 hoomd_organics/tests/utils/test_utils.py | 26 +++++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/hoomd_organics/tests/utils/test_utils.py b/hoomd_organics/tests/utils/test_utils.py
index b132cda7..16d38462 100644
--- a/hoomd_organics/tests/utils/test_utils.py
+++ b/hoomd_organics/tests/utils/test_utils.py
@@ -1,7 +1,11 @@
 import pytest
 import unyt as u
 
-from hoomd_organics.utils import check_return_iterable, validate_ref_value
+from hoomd_organics.utils import (
+    calculate_box_length,
+    check_return_iterable,
+    validate_ref_value,
+)
 from hoomd_organics.utils.exceptions import ReferenceUnitError
 
 
@@ -38,3 +42,23 @@ def test_validate_ref_value(self):
 
         with pytest.raises(ValueError):
             validate_ref_value("test g", u.dimensions.mass)
+
+    def test_calculate_box_length(self):
+        mass = u.unyt_quantity(4.0, u.g)
+        density = u.unyt_quantity(0.5, u.g / u.cm**3)
+        box_length = calculate_box_length(mass, density)
+        assert box_length == 2.0 * u.cm
+
+    def test_calculate_box_length_fixed_l_1d(self):
+        mass = u.unyt_quantity(6.0, u.g)
+        density = u.unyt_quantity(0.5, u.g / u.cm**3)
+        fixed_L = u.unyt_quantity(3.0, u.cm)
+        box_length = calculate_box_length(mass, density, fixed_L=fixed_L)
+        assert box_length == 2.0 * u.cm
+
+    def test_calculate_box_length_fixed_l_2d(self):
+        mass = u.unyt_quantity(12.0, u.g)
+        density = u.unyt_quantity(0.5, u.g / u.cm**3)
+        fixed_L = u.unyt_array([3.0, 2.0], u.cm)
+        box_length = calculate_box_length(mass, density, fixed_L=fixed_L)
+        assert box_length == 4.0 * u.cm

From 4496eb5912f15cee2f97342e200ff2106249808a Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 12:49:31 -0600
Subject: [PATCH 34/41] fix the bug is getting fixed_L size

---
 hoomd_organics/utils/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/hoomd_organics/utils/utils.py b/hoomd_organics/utils/utils.py
index 6536c15e..75826518 100644
--- a/hoomd_organics/utils/utils.py
+++ b/hoomd_organics/utils/utils.py
@@ -122,7 +122,7 @@ def calculate_box_length(mass, density, fixed_L=None):
         L = vol ** (1 / 3)
     else:
         L = vol / np.prod(fixed_L)
-        if len(fixed_L) == 1:  # L is cm^2
+        if fixed_L.size == 1:  # L is cm^2
             L = L ** (1 / 2)
     # L is cm
     return L

From 071f62f60ffa4252866843d64ea314992c80bee3 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 14:03:53 -0600
Subject: [PATCH 35/41] fix the bug when checking fixed_L is none

---
 hoomd_organics/base/system.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py
index cbec2fbc..1e3af69d 100644
--- a/hoomd_organics/base/system.py
+++ b/hoomd_organics/base/system.py
@@ -379,7 +379,8 @@ def _calculate_L(self, fixed_L=None):
         """
         mass_quantity = u.unyt_quantity(self.mass, u.g / u.mol).to("g")
         density_quantity = u.unyt_quantity(self.density, u.g / u.cm**3)
-        fixed_L = u.unyt_array(fixed_L, u.cm)
+        if fixed_L is not None:
+            fixed_L = u.unyt_array(fixed_L, u.cm)
         L = calculate_box_length(
             mass_quantity, density_quantity, fixed_L=fixed_L
         )

From 2f90756c91f602eeff95accb05435939b9de68dd Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 15:17:52 -0600
Subject: [PATCH 36/41] add more unit tests

---
 hoomd_organics/base/simulation.py            |  2 +-
 hoomd_organics/tests/base/test_simulation.py | 27 ++++++++++++++++++++
 2 files changed, 28 insertions(+), 1 deletion(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index 933acb43..c7d6e8a0 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -449,7 +449,7 @@ def run_update_volume(
                 Lz=final_box_lengths[2],
             )
         else:
-            if self.reference_values is None:
+            if not self.reference_values:
                 raise ReferenceUnitError(
                     "Missing simulation units. Please "
                     "provide units for mass, length, and"
diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 97373c98..8643e330 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -9,6 +9,7 @@
 
 from hoomd_organics import Simulation
 from hoomd_organics.tests import BaseTest
+from hoomd_organics.utils.exceptions import ReferenceUnitError
 
 
 class TestSimulate(BaseTest):
@@ -138,6 +139,32 @@ def test_update_volume_missing_values(self, benzene_system):
         with pytest.raises(ValueError):
             sim.run_update_volume(kT=1.0, tau_kt=0.01, n_steps=500, period=1)
 
+    def test_update_volume_two_values(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        with pytest.raises(ValueError):
+            sim.run_update_volume(
+                kT=1.0,
+                tau_kt=0.01,
+                n_steps=500,
+                period=1,
+                final_box_lengths=sim.box_lengths_reduced * 0.5,
+                final_density=0.1,
+            )
+
+    def test_update_volume_with_density_no_ref_values(self, benzene_system):
+        sim = Simulation(
+            initial_state=benzene_system.hoomd_snapshot,
+            forcefield=benzene_system.hoomd_forcefield,
+        )
+        with pytest.raises(ReferenceUnitError):
+            sim.run_update_volume(
+                kT=1.0,
+                tau_kt=0.01,
+                n_steps=500,
+                period=1,
+                final_density=0.1,
+            )
+
     def test_change_methods(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)

From f8fb43b9bd6a8693b3799e70dc525ba4e0b7adf1 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 16:43:02 -0600
Subject: [PATCH 37/41] validate reference values in setters

---
 hoomd_organics/base/simulation.py | 63 ++++++++-----------------------
 1 file changed, 16 insertions(+), 47 deletions(-)

diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py
index c7d6e8a0..910b1dcb 100644
--- a/hoomd_organics/base/simulation.py
+++ b/hoomd_organics/base/simulation.py
@@ -7,7 +7,12 @@
 import numpy as np
 import unyt as u
 
-from hoomd_organics.utils import StdOutLogger, UpdateWalls, calculate_box_length
+from hoomd_organics.utils import (
+    StdOutLogger,
+    UpdateWalls,
+    calculate_box_length,
+    validate_ref_value,
+)
 from hoomd_organics.utils.exceptions import ReferenceUnitError
 
 
@@ -145,50 +150,19 @@ def reference_values(self):
         return self._reference_values
 
     @reference_length.setter
-    def reference_length(self, length, unit=None):
-        if isinstance(length, u.array.unyt_quantity):
-            self._reference_values["length"] = length
-        elif isinstance(unit, str) and (
-            isinstance(length, float) or isinstance(length, int)
-        ):
-            self._reference_values["length"] = length * getattr(u, unit)
-        else:
-            raise ReferenceUnitError(
-                f"Invalid reference length input.Please provide reference "
-                f"length (number) and unit (string) or pass length value as an "
-                f"{str(u.array.unyt_quantity)}."
-            )
+    def reference_length(self, length):
+        validated_length = validate_ref_value(length, u.dimensions.length)
+        self._reference_values["length"] = validated_length
 
     @reference_energy.setter
-    def reference_energy(self, energy, unit=None):
-        if isinstance(energy, u.array.unyt_quantity):
-            self._reference_values["energy"] = energy
-        elif isinstance(unit, str) and (
-            isinstance(energy, float) or isinstance(energy, int)
-        ):
-            self._reference_values["energy"] = energy * getattr(u, unit)
-        else:
-            raise ReferenceUnitError(
-                f"Invalid reference energy input.Please provide reference "
-                f"energy (number) and unit (string) or pass energy value as an "
-                f"{str(u.array.unyt_quantity)}."
-            )
+    def reference_energy(self, energy):
+        validated_energy = validate_ref_value(energy, u.dimensions.energy)
+        self._reference_values["energy"] = validated_energy
 
     @reference_mass.setter
-    def reference_mass(self, mass, unit=None):
-        if isinstance(mass, u.array.unyt_quantity):
-            self._reference_values["mass"] = mass
-        elif isinstance(unit, str) and (
-            isinstance(mass, float) or isinstance(mass, int)
-        ):
-            self._reference_values["mass"] = mass * getattr(u, unit)
-        else:
-            raise ReferenceUnitError(
-                f"Invalid reference mass input.Please provide reference "
-                f"mass (number) and "
-                f"unit (string) or pass mass value as an "
-                f"{str(u.array.unyt_quantity)}."
-            )
+    def reference_mass(self, mass):
+        validated_mass = validate_ref_value(mass, u.dimensions.mass)
+        self._reference_values["mass"] = validated_mass
 
     @reference_values.setter
     def reference_values(self, ref_value_dict):
@@ -196,12 +170,7 @@ def reference_values(self, ref_value_dict):
         for k in ref_keys:
             if k not in ref_value_dict.keys():
                 raise ValueError(f"Missing reference for {k}.")
-            if not isinstance(ref_value_dict[k], u.array.unyt_quantity):
-                raise ReferenceUnitError(
-                    f"{k} reference value must be of type "
-                    f"{str(u.array.unyt_quantity)}"
-                )
-        self._reference_values = ref_value_dict
+            self.__setattr__(f"reference_{k}", ref_value_dict[k])
 
     @property
     def box_lengths_reduced(self):

From e798c8e48cc8ba5dc9efa2b1155165af61e69980 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 16:44:11 -0600
Subject: [PATCH 38/41] test for box_length without ref values

---
 hoomd_organics/tests/base/test_system.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py
index eadd4ed7..af303128 100644
--- a/hoomd_organics/tests/base/test_system.py
+++ b/hoomd_organics/tests/base/test_system.py
@@ -216,6 +216,18 @@ def test_target_box(self, benzene_molecule):
             low_density_system.target_box > high_density_system.target_box
         )
 
+    def test_target_box_without_ref_values(self, benzene_molecule):
+        benzene_mol = benzene_molecule(n_mols=3)
+        system_1 = Pack(
+            molecules=[benzene_mol],
+            density=0.1,
+            r_cut=2.5,
+            auto_scale=True,
+        )
+        assert np.allclose(
+            system_1.target_box, [1.57289499, 1.57289499, 1.57289499], atol=1e-4
+        )
+
     def test_mass(self, pps_molecule):
         pps_mol = pps_molecule(n_mols=20)
         system = Pack(molecules=[pps_mol], density=1.0, r_cut=2.5)

From 665868ba9d5d63b7551c9534022c540cd0fa273b Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 17:28:37 -0600
Subject: [PATCH 39/41] remove unit test

---
 hoomd_organics/tests/base/test_system.py | 12 ------------
 1 file changed, 12 deletions(-)

diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py
index af303128..eadd4ed7 100644
--- a/hoomd_organics/tests/base/test_system.py
+++ b/hoomd_organics/tests/base/test_system.py
@@ -216,18 +216,6 @@ def test_target_box(self, benzene_molecule):
             low_density_system.target_box > high_density_system.target_box
         )
 
-    def test_target_box_without_ref_values(self, benzene_molecule):
-        benzene_mol = benzene_molecule(n_mols=3)
-        system_1 = Pack(
-            molecules=[benzene_mol],
-            density=0.1,
-            r_cut=2.5,
-            auto_scale=True,
-        )
-        assert np.allclose(
-            system_1.target_box, [1.57289499, 1.57289499, 1.57289499], atol=1e-4
-        )
-
     def test_mass(self, pps_molecule):
         pps_mol = pps_molecule(n_mols=20)
         system = Pack(molecules=[pps_mol], density=1.0, r_cut=2.5)

From d63a80f1d8c0af5ae3c7ee544b98bf960225e2f5 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Thu, 14 Sep 2023 17:57:22 -0600
Subject: [PATCH 40/41] add unit tests for ref values

---
 hoomd_organics/tests/base/test_simulation.py | 39 ++++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 8643e330..66479d01 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -57,6 +57,45 @@ def test_reference_values(self, benzene_system):
         assert np.isclose(float(sim.mass.value), benzene_system.mass, atol=1e-4)
         assert np.allclose(benzene_system.box.lengths, sim.box_lengths.value)
 
+    def test_set_ref_values(self, benzene_system):
+        sim = Simulation(
+            initial_state=benzene_system.hoomd_snapshot,
+            forcefield=benzene_system.hoomd_forcefield,
+        )
+        ref_value_dict = {
+            "length": 1 * u.angstrom,
+            "energy": 3.0 * u.kcal / u.mol,
+            "mass": 1.25 * u.Unit("amu"),
+        }
+        sim.reference_values = ref_value_dict
+        assert sim.reference_length == ref_value_dict["length"]
+        assert sim.reference_energy == ref_value_dict["energy"]
+        assert sim.reference_mass == ref_value_dict["mass"]
+
+    def test_set_ref_length(self, benzene_system):
+        sim = Simulation(
+            initial_state=benzene_system.hoomd_snapshot,
+            forcefield=benzene_system.hoomd_forcefield,
+        )
+        sim.reference_length = 1 * u.angstrom
+        assert sim.reference_length == 1 * u.angstrom
+
+    def test_set_ref_energy(self, benzene_system):
+        sim = Simulation(
+            initial_state=benzene_system.hoomd_snapshot,
+            forcefield=benzene_system.hoomd_forcefield,
+        )
+        sim.reference_energy = 3.0 * u.kcal / u.mol
+        assert sim.reference_energy == 3.0 * u.kcal / u.mol
+
+    def test_set_ref_mass(self, benzene_system):
+        sim = Simulation(
+            initial_state=benzene_system.hoomd_snapshot,
+            forcefield=benzene_system.hoomd_forcefield,
+        )
+        sim.reference_mass = 1.25 * u.amu
+        assert sim.reference_mass == 1.25 * u.amu
+
     def test_NVT(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500)

From 332b46a0ba170b1be96650664c067ef47ff3cbe6 Mon Sep 17 00:00:00 2001
From: marjanalbouye <marjan.albouye73@gmail.com>
Date: Fri, 15 Sep 2023 11:56:02 -0600
Subject: [PATCH 41/41] fix simulation volume update unit test

---
 hoomd_organics/tests/base/test_simulation.py | 27 ++++++++++----------
 hoomd_organics/tests/base_test.py            |  4 +++
 2 files changed, 17 insertions(+), 14 deletions(-)

diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py
index 66479d01..65fd40c2 100644
--- a/hoomd_organics/tests/base/test_simulation.py
+++ b/hoomd_organics/tests/base/test_simulation.py
@@ -9,7 +9,6 @@
 
 from hoomd_organics import Simulation
 from hoomd_organics.tests import BaseTest
-from hoomd_organics.utils.exceptions import ReferenceUnitError
 
 
 class TestSimulate(BaseTest):
@@ -190,19 +189,19 @@ def test_update_volume_two_values(self, benzene_system):
                 final_density=0.1,
             )
 
-    def test_update_volume_with_density_no_ref_values(self, benzene_system):
-        sim = Simulation(
-            initial_state=benzene_system.hoomd_snapshot,
-            forcefield=benzene_system.hoomd_forcefield,
-        )
-        with pytest.raises(ReferenceUnitError):
-            sim.run_update_volume(
-                kT=1.0,
-                tau_kt=0.01,
-                n_steps=500,
-                period=1,
-                final_density=0.1,
-            )
+    # def test_update_volume_with_density_no_ref_values(self, benzene_system):
+    #     sim_no_ref = Simulation(
+    #         initial_state=benzene_system.hoomd_snapshot,
+    #         forcefield=benzene_system.hoomd_forcefield,
+    #     )
+    #     with pytest.raises(ReferenceUnitError):
+    #         sim_no_ref.run_update_volume(
+    #             kT=1.0,
+    #             tau_kt=0.01,
+    #             n_steps=500,
+    #             period=1,
+    #             final_density=0.1,
+    #         )
 
     def test_change_methods(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py
index 8a69d85b..c6666dbe 100644
--- a/hoomd_organics/tests/base_test.py
+++ b/hoomd_organics/tests/base_test.py
@@ -12,6 +12,10 @@
 
 
 class BaseTest:
+    @pytest.fixture(autouse=True)
+    def initdir(self, tmpdir):
+        tmpdir.chdir()
+
     @pytest.fixture()
     def benzene_smiles(self):
         return "c1ccccc1"