Merge pull request #112 from chrisjonesBSU/number-density

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

.gitignore

@@ -2,6 +2,7 @@
 __pycache__/
 *.py[cod]
 *$py.class
+*.swp
 
 # C extensions
 *.so

.pre-commit-config.yaml

@@ -43,5 +43,5 @@ repos:
     rev: '6.3.0'
     hooks:
       - id: pydocstyle
-        exclude: ^(flowermd/tests/|flowermd/utils/|setup.py|flowermd/__version__.py|docs/)
+        exclude: ^(flowermd/tests/|flowermd/internal/|flowermd/utils|setup.py|flowermd/__version__.py|docs/)
         args: [ --convention=numpy ]

docs/source/index.rst

@@ -52,6 +52,7 @@ Resources
    base
    modules
    library
+   utils
 
 
 

docs/source/utils.rst

@@ -0,0 +1,18 @@
+Utils
+=======
+
+
+.. py:currentmodule:: flowermd.utils
+
+Utils
+-----
+
+.. automodule:: flowermd.utils.utils
+   :members:
+
+
+Base Types
+----------
+
+.. automodule:: flowermd.utils.base_types
+   :members:

flowermd/base/molecule.py

@@ -14,9 +14,9 @@
 from mbuild.lib.recipes import Polymer as mbPolymer
 
 from flowermd.base import BaseHOOMDForcefield, BaseXMLForcefield
-from flowermd.utils import check_return_iterable
-from flowermd.utils.exceptions import ForceFieldError, MoleculeLoadError
-from flowermd.utils.ff_utils import _validate_hoomd_ff
+from flowermd.internal import check_return_iterable
+from flowermd.internal.exceptions import ForceFieldError, MoleculeLoadError
+from flowermd.internal.ff_utils import _validate_hoomd_ff
 
 
 class Molecule:

flowermd/base/simulation.py

@@ -9,14 +9,9 @@
 import numpy as np
 import unyt as u
 
-from flowermd.utils import (
-    HOOMDThermostats,
-    StdOutLogger,
-    UpdateWalls,
-    calculate_box_length,
-    validate_ref_value,
-)
-from flowermd.utils.exceptions import ReferenceUnitError
+from flowermd.internal import validate_ref_value
+from flowermd.utils.actions import StdOutLogger, UpdateWalls
+from flowermd.utils.base_types import HOOMDThermostats
 
 
 class Simulation(hoomd.simulation.Simulation):
@@ -586,24 +581,27 @@ def remove_walls(self, wall_axis):
 
     def run_update_volume(
         self,
+        final_box_lengths,
         n_steps,
         period,
         kT,
         tau_kt,
-        final_box_lengths=None,
-        final_density=None,
         thermalize_particles=True,
         write_at_start=True,
     ):
         """Run an NVT simulation while shrinking or expanding simulation box.
 
         The simulation box is updated using `hoomd.update.BoxResize` and the
-        final box lengths are set to `final_box_lengths` or inferred from
-        `final_density` if `final_box_lengths` is not provided.
+        final box lengths are set to `final_box_lengths`.
 
+        See `flowermd.utils.get_target_volume_mass_density` and
+        `flowermd.utils.get_target_volume_number_density` which are
+        helper functions that can be used to get `final_box_lengths`.
 
         Parameters
         ----------
+        final_box_lengths : np.ndarray or unyt.array.unyt_array, shape=(3,), required # noqa: E501
+            The final box edge lengths in (x, y, z) order.
         n_steps : int, required
             Number of steps to run during volume update.
         period : int, required
@@ -612,52 +610,54 @@ def run_update_volume(
             The temperature to use during volume update.
         tau_kt : float, required
             Thermostat coupling period (in simulation time units).
-        final_box_lengths : np.ndarray, shape=(3,), dtype=float, default None
-            The final box edge lengths in (x, y, z) order.
-        final_density : float, default None
-            The final density of the simulation in g/cm^3.
         write_at_start : bool, default True
             When set to True, triggers writers that evaluate to True
             for the initial step to execute before the next simulation
             time step.
 
-        """
-        if final_box_lengths is None and final_density is None:
-            raise ValueError(
-                "Must provide either `final_box_lengths` or `final_density`"
+        Examples
+        --------
+        In this example, a low density system is initialized with `Pack`
+        and a box matching a density of 1.1 g/cm^3 is passed into
+        `final_box_lengths`.
+
+        ::
+
+            import unyt
+            from flowermd.base import Pack, Simulation
+            from flowermd.library import PPS, OPLS_AA_PPS
+
+            pps_mols = PPS(num_mols=20, lengths=15)
+            pps_system = Pack(
+                molecules=[pps_mols],
+                force_field=OPLS_AA_PPS(),
+                r_cut=2.5,
+                density=0.5,
+                auto_scale=True,
+                scale_charges=True
             )
-        if final_box_lengths is not None and final_density is not None:
-            raise ValueError(
-                "Cannot provide both `final_box_lengths` and `final_density`."
+            sim = Simulation(
+                initial_state=pps_system.hoomd_snapshot,
+                forcefield=pps_system.hoomd_forcefield
             )
-        if final_box_lengths is not None:
-            final_box = hoomd.Box(
-                Lx=final_box_lengths[0],
-                Ly=final_box_lengths[1],
-                Lz=final_box_lengths[2],
+            target_box = flowermd.utils.get_target_box_mass_density(
+                density=1.1 * unyt.g/unyt.cm**3, mass=sim.mass.to("g")
+            )
+            sim.run_update_volume(
+                n_steps=1e4, kT=1.0, tau_kt=1.0, final_box_lengths=target_box
             )
-        else:
-            if not self.reference_values:
-                raise ReferenceUnitError(
-                    "Missing simulation units. Please "
-                    "provide units for mass, length, and"
-                    " energy."
-                )
-
-            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
-                )
-            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
-            ).value
-            final_box = hoomd.Box(Lx=L, Ly=L, Lz=L)
 
+        """
+        if self.reference_length and hasattr(final_box_lengths, "to"):
+            ref_unit = self.reference_length.units
+            final_box_lengths = final_box_lengths.to(ref_unit)
+            final_box_lengths /= self.reference_length
+
+        final_box = hoomd.Box(
+            Lx=final_box_lengths[0],
+            Ly=final_box_lengths[1],
+            Lz=final_box_lengths[2],
+        )
         resize_trigger = hoomd.trigger.Periodic(period)
         box_ramp = hoomd.variant.Ramp(
             A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps)