Add ability to use number density when calculating target box lenghts for shrinking and NVT.

flowermd/base/simulation.py

@@ -651,7 +651,7 @@ def run_update_volume(
                     final_density, u.g / u.cm**3
                 )
             mass_g = self.mass.to("g")
-            L = calculate_box_length(mass_g, density_quantity)
+            L = calculate_box_length(mass=mass_g, density=density_quantity)
             # convert L from cm to reference units
             L = (
                 L.to(self.reference_length.units) / self.reference_length.value

flowermd/base/system.py

@@ -43,11 +43,6 @@ class System(ABC):
     ----------
     molecules : flowermd.Molecule or a list of Molecule objects, required
         The molecules to be placed in the system.
-    density : float, required
-        The desired density of the system (g/cm^3). Used to set the
-        target_box attribute. Can be useful when initializing
-        systems at low density and running a shrink simulation
-        to achieve a target density.
     base_units : dict, default {}
         Dictionary of base units to use for scaling.
         Dictionary keys are "length", "mass", and "energy". Values should be an
@@ -71,11 +66,9 @@ class System(ABC):
     def __init__(
         self,
         molecules,
-        density: float,
         base_units=dict(),
     ):
         self._molecules = check_return_iterable(molecules)
-        self.density = density
         self.all_molecules = []
         self.gmso_system = None
         self._reference_values = base_units
@@ -623,7 +616,7 @@ def apply_forcefield(
         self._hoomd_snapshot = self._create_hoomd_snapshot()
 
     def set_target_box(
-        self, x_constraint=None, y_constraint=None, z_constraint=None
+        self, density, x_constraint=None, y_constraint=None, z_constraint=None
     ):
         """Set the target box size of the system.
 
@@ -633,6 +626,8 @@ def set_target_box(
 
         Parameters
         ----------
+        density : float, required
+            Density used when calculating the lengths (g/cm^3).
         x_constraint : float, optional, defualt=None
             Fixes the box length (nm) along the x axis.
         y_constraint : float, optional, default=None
@@ -642,13 +637,13 @@ def set_target_box(
 
         """
         if not any([x_constraint, y_constraint, z_constraint]):
-            Lx = Ly = Lz = self._calculate_L()
+            Lx = Ly = Lz = self._calculate_L(density=density)
         else:
             constraints = np.array([x_constraint, y_constraint, z_constraint])
             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)
+            L = self._calculate_L(density=density, fixed_L=fixed_L)
             constraints[np.equal(constraints, None).nonzero()] = L
             Lx, Ly, Lz = constraints
 
@@ -663,7 +658,7 @@ def visualize(self):
                 "The initial configuraiton has not been created yet."
             )
 
-    def _calculate_L(self, fixed_L=None):
+    def _calculate_L(self, density, fixed_L=None):
         """Calculate the box length.
 
         Calculate the required box length(s) given the mass of a system and
@@ -680,11 +675,11 @@ 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)
+        density_quantity = u.unyt_quantity(density, u.g / u.cm**3)
         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
+            mass=mass_quantity, density=density_quantity, fixed_L=fixed_L
         )
         return L.to("nm").value
 
@@ -697,6 +692,11 @@ class Pack(System):
 
     Parameters
     ----------
+    density : float, required
+        The desired density of the system (g/cm^3). Used to set the
+        target_box attribute. Can be useful when initializing
+        systems at low density and running a shrink simulation
+        to achieve a target density.
     packing_expand_factor : int, default 5
         The factor by which to expand the box for packing.
     edge : float, default 0.2
@@ -731,17 +731,17 @@ def __init__(
         edge=0.2,
         overlap=0.2,
     ):
+        self.density = density * u.g / (u.cm ** (3))
         self.packing_expand_factor = packing_expand_factor
         self.edge = edge
         self.overlap = overlap
         super(Pack, self).__init__(
             molecules=molecules,
-            density=density,
             base_units=base_units,
         )
 
     def _build_system(self):
-        self.set_target_box()
+        self.set_target_box(density=self.density)
         system = mb.packing.fill_box(
             compound=self.all_molecules,
             n_compounds=[1 for i in self.all_molecules],
@@ -773,7 +773,6 @@ class Lattice(System):
     def __init__(
         self,
         molecules,
-        density: float,
         x: float,
         y: float,
         n: int,
@@ -786,7 +785,6 @@ def __init__(
         self.basis_vector = basis_vector
         super(Lattice, self).__init__(
             molecules=molecules,
-            density=density,
             base_units=base_units,
         )
 
@@ -820,7 +818,6 @@ def _build_system(self):
         x_len = bounding_box.lengths[0] + self.x
         y_len = bounding_box.lengths[1] + self.y
         z_len = bounding_box.lengths[2] + 0.2
-        self.set_target_box(x_constraint=x_len, y_constraint=y_len)
         # Center the lattice in its box
         system.box = mb.box.Box(np.array([x_len, y_len, z_len]))
         system.translate_to(

flowermd/library/surfaces.py

@@ -93,7 +93,6 @@ def __init__(
         surface_mol = Molecule(num_mols=1, compound=surface)
         super(Graphene, self).__init__(
             molecules=[surface_mol],
-            density=1.0,
             base_units=base_units,
         )
 

flowermd/tests/base/test_system.py

@@ -782,7 +782,6 @@ def test_lattice_polymer(self, polyethylene):
         polyethylene = polyethylene(lengths=2, num_mols=32)
         system = Lattice(
             molecules=[polyethylene],
-            density=1.0,
             x=1,
             y=1,
             n=4,
@@ -807,7 +806,6 @@ def test_lattice_molecule(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=32)
         system = Lattice(
             molecules=[benzene_mol],
-            density=1.0,
             x=1,
             y=1,
             n=4,

flowermd/tests/library/test_tensile.py

@@ -10,7 +10,6 @@ def test_tensile(self):
         pps = PPS(lengths=6, num_mols=32)
         system = Lattice(
             molecules=[pps],
-            density=1.0,
             x=1.2,
             y=1.2,
             n=4,

flowermd/tests/utils/test_utils.py

@@ -43,22 +43,44 @@ def test_validate_ref_value(self):
         with pytest.raises(ValueError):
             validate_ref_value("test g", u.dimensions.mass)
 
-    def test_calculate_box_length(self):
+    def test_calculate_box_length_mass_density(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)
+        box_length = calculate_box_length(density=density, mass=mass)
         assert box_length == 2.0 * u.cm
 
+    def test_calculate_box_length_number_density(self):
+        sigma = 1 * u.nm
+        n_beads = 100
+        density = 1 / sigma**3
+        box_length = calculate_box_length(density=density, n_beads=n_beads)
+        assert box_length == 100 ** (1 / 3) * u.nm
+
+    def test_calculate_box_length_bad_args(self):
+        mass_density = 1 * u.g / (u.cm**3)
+        number_density = 1 / (1 * u.nm) ** 3
+        invalid_density = 1 * u.J / u.kg
+        with pytest.raises(ValueError):
+            calculate_box_length(density=mass_density, n_beads=100)
+        with pytest.raises(ValueError):
+            calculate_box_length(density=number_density, mass=100)
+        with pytest.raises(ValueError):
+            calculate_box_length(density=invalid_density, mass=100)
+
     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)
+        box_length = calculate_box_length(
+            mass=mass, density=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)
+        box_length = calculate_box_length(
+            mass=mass, density=density, fixed_L=fixed_L
+        )
         assert box_length == 4.0 * u.cm

flowermd/utils/utils.py

@@ -96,21 +96,26 @@ def _is_float(num):
         )
 
 
-def calculate_box_length(mass, density, fixed_L=None):
+def calculate_box_length(density, mass=None, n_beads=None, 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.
+    #TODO: Add example of using box constraints
 
     Parameters
     ----------
-    mass : unyt.unyt_quantity, required
-        Mass of the system
     density : unyt.unyt_quantity, required
         Target density of the system
+    mass : unyt.unyt_quantity, required
+        Mass of the system.
+        Use when using mass density rather than number density.
+    n_beads : int, optional
+        Number of beads in the system.
+        Use when using number density rather than mass density.
     fixed_L : np.array, optional, defualt=None
         Array of fixed box lengths to be accounted for
         when solving for L
@@ -120,13 +125,33 @@ def calculate_box_length(mass, density, fixed_L=None):
     L : float
         Box edge length
     """
-
-    vol = mass / density  # cm^3
+    # Check units of density
+    mass_density = u.Unit("kg") / u.Unit("m**3")
+    number_density = u.Unit("m**-3")
+    if density.units.dimensions == mass_density.dimensions:
+        if not mass:
+            raise ValueError(
+                f"The given density has units of {mass_density.dimensions} "
+                "but the mass is not given."
+            )
+        vol = mass / density
+    elif density.units.dimensions == number_density.dimensions:
+        if not n_beads:
+            raise ValueError(
+                f"The given density has units of {number_density.dimensions} "
+                "but the number of beads is not given."
+            )
+        vol = n_beads / density
+    else:
+        raise ValueError(
+            f"Density units of {density.units.dimensions} were given. "
+            f"Only mass density ({mass_density.units.dimensions}) and "
+            f"number density ({number_density.dimensions}) are supported."
+        )
     if fixed_L is None:
         L = vol ** (1 / 3)
     else:
         L = vol / np.prod(fixed_L)
-        if fixed_L.size == 1:  # L is cm^2
+        if fixed_L.size == 1:  # L is units of area
             L = L ** (1 / 2)
-    # L is cm
     return L