Update simulation interface and add density to update_volume

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
 
 
@@ -52,6 +52,7 @@ def __init__(
         self,
         initial_state,
         forcefield=None,
+        reference_values=dict(),
         r_cut=2.5,
         dt=0.0001,
         device=hoomd.device.auto_select(),
@@ -83,12 +84,30 @@ def __init__(
         self.integrator = None
         self._dt = dt
         self._reference_values = dict()
+        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:
@@ -379,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
@@ -395,20 +415,33 @@ 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:
+            # 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(
             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,

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,
@@ -62,7 +63,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 +209,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.value
+
     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 +349,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:
@@ -356,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
@@ -370,25 +375,29 @@ 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):
     """Uses PACKMOL via mbuild.packing.fill_box.
     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 +437,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,
         )

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,
+)

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