Issue 1080

.pre-commit-config.yaml

@@ -17,7 +17,7 @@ repos:
     - id: trailing-whitespace
       exclude: setup.cfg
 - repo: https://github.com/psf/black
-  rev: 22.12.0
+  rev: 23.1.0
   hooks:
     - id: black
       args: [--line-length=80]

docs/sphinxext/notebook_sphinxext.py

@@ -61,7 +61,6 @@ class NotebookDirective(Directive):
     final_argument_whitespace = True
 
     def run(self):
-
         # check if raw html is supported
         if not self.state.document.settings.raw_enabled:
             raise self.warning('"%s" directive disabled.' % self.name)

mbuild/compound.py

@@ -2427,7 +2427,6 @@ def _energy_minimize_openbabel(
                 # if an individual entry is a list, validate the input
                 if isinstance(fixed_temp, list):
                     if len(fixed_temp) == 2:
-
                         msg1 = (
                             "Expected tuple or list of length 3 to set"
                             "which dimensions to fix motion."

mbuild/formats/cassandramcf.py

@@ -165,7 +165,6 @@ def write_mcf(
 
     # Now we write the MCF file
     with open(filename, "w") as mcf_file:
-
         header = (
             "!***************************************"
             "****************************************\n"

mbuild/formats/lammpsdata.py

@@ -37,6 +37,7 @@ def write_lammpsdata(
     sigma_conversion_factor=None,
     epsilon_conversion_factor=None,
     mass_conversion_factor=None,
+    charge_conversion_factor=True,
     zero_dihedral_weighting_factor=False,
     moleculeID_offset=1,
 ):
@@ -52,51 +53,56 @@ def write_lammpsdata(
         ParmEd structure object
     filename : str
         Path of the output file
-    atom_style: str
+    atom_style: str, optional, default='full'
         Defines the style of atoms to be saved in a LAMMPS data file. The
         following atom styles are currently supported:
         'full', 'atomic', 'charge', 'molecular'
         see http://lammps.sandia.gov/doc/atom_style.html for more information
         on atom styles.
-    unit_style: str
-        Defines to unit style to be save in a LAMMPS data file.  Defaults to
-        'real' units. Current styles are supported: 'real', 'lj'
-        see https://lammps.sandia.gov/doc/99/units.html for more information
-        on unit styles
-    mins : list
-        minimum box dimension in x, y, z directions, nm
-    maxs : list
-        maximum box dimension in x, y, z directions, nm
-    pair_coeff_label : str
+    unit_style : str, optional, default='real'
+        Defines to unit style to be save in a LAMMPS data file. Current styles
+        are supported: 'real', 'lj', see lammps unit style documentation:
+        https://lammps.sandia.gov/doc/99/units.html for more information.
+    mins : list, optional, default=None
+        Minimum box dimension in x, y, z directions, nm
+    maxs : list, optional, default=None
+        Maximum box dimension in x, y, z directions, nm
+    pair_coeff_label : str, optional, default=None
         Provide a custom label to the pair_coeffs section in the lammps data
-        file. Defaults to None, which means a suitable default will be chosen.
-    detect_forcefield_style: boolean
+        file. A value of None means a suitable default will be chosen.
+    detect_forcefield_style : bool, optional, default=True
         If True, format lammpsdata parameters based on the contents of
         the parmed Structure
-    use_urey_bradleys: boolean
+    use_urey_bradleys : bool, optional, default=False
         If True, will treat angles as CHARMM-style angles with urey bradley
         terms while looking for `structure.urey_bradleys`
-    use_rb_torsions:
+    use_rb_torsions : bool, optional, default=True
         If True, will treat dihedrals OPLS-style torsions while looking for
         `structure.rb_torsions`
-    use_dihedrals:
+    use_dihedrals : bool, optional, default=False
         If True, will treat dihedrals as CHARMM-style dihedrals while looking
         for `structure.dihedrals`
-    zero_dihedral_weighting_factor:
+    zero_dihedral_weighting_factor : bool, optional, default=False
         If True, will set weighting parameter to zero in CHARMM-style dihedrals.
         This should be True if the CHARMM dihedral style is used in non-CHARMM forcefields.
-    sigma_conversion_factor: None, float
+    sigma_conversion_factor : float, optional, default=None
         If unit_style is set to 'lj', then sigma conversion factor is used to non-dimensionalize.
-        Assume to be in units of nm. Default is None. If None, will take the largest sigma value in
+        Assume to be in units of nm. If None, will take the largest sigma value in
         the structure.atoms.sigma values.
-    epsilon_conversion_factor: None, float
+    epsilon_conversion_factor : float, optional, default=None
         If unit_style is set to 'lj', then epsilon conversion factor is used to non-dimensionalize.
-        Assume to be in units of kCal/mol. Default is None. If None, will take the largest epsilon value in
+        Assume to be in units of kCal/mol. If None, will take the largest epsilon value in
         the structure.atoms.epsilon values.
-    mass_conversion_factor: None, float
+    mass_conversion_factor : float, optional, default=None
         If unit_style is set to 'lj', then mass conversion factor is used to non-dimensionalize.
-        Assume to be in units of amu. Default is None. If None, will take the largest mass value in
+        Assume to be in units of amu. If None, will take the largest mass value in
         the structure.atoms.mass values.
+    charge_conversion_factor : bool, optional, default=True
+        If unit_style is set to 'lj', then charge conversion factor may or may not be used to
+        non-dimensionalize. Assume to be in elementary charge units. If ``True``, the charges
+        are scaled by ``np.sqrt(4*np.pi()*eps_0*sigma_conversion_factor*epsilon_conversion_factor)``.
+        If ``False``, the charges are not scaled and the user must be wary to choose the dielectric
+        constant properly, which may be more convenient to implement an implicit solvent.
     moleculeID_offset : int , optional, default=1
         Since LAMMPS treats the MoleculeID as an additional set of information
         to identify what molecule an atom belongs to, this currently
@@ -287,16 +293,17 @@ def write_lammpsdata(
         )
         # Convert coordinates and charges to LJ units
         xyz = xyz / sigma_conversion_factor
-        charges = (charges * ELEM_TO_COUL) / np.sqrt(
-            4
-            * np.pi
-            * sigma_conversion_factor
-            * NM_TO_ANG**-1
-            * epsilon_conversion_factor
-            * KCAL_TO_KJ
-            * epsilon_0
-            * 10**-6
-        )
+        if charge_conversion_factor:
+            charges = (charges * ELEM_TO_COUL) / np.sqrt(
+                4
+                * np.pi
+                * sigma_conversion_factor
+                * NM_TO_ANG**-1
+                * epsilon_conversion_factor
+                * KCAL_TO_KJ
+                * epsilon_0
+                * 10**-6
+            )
         charges[np.isinf(charges)] = 0
     else:
         sigma_conversion_factor = 1
@@ -360,7 +367,6 @@ def write_lammpsdata(
 
     # Write lammps data file https://docs.lammps.org/2001/data_format.html
     with open(filename, "w") as data:
-
         data.write(f"{filename} - created by mBuild; units = {unit_style}\n")
         if unit_style == "lj":
             data.write("#Normalization factors: ")
@@ -938,7 +944,6 @@ def _get_dihedral_types(
             for dihedral in structure.rb_torsions
         ]
     elif use_dihedrals:
-
         dihedral_types = [
             unique_dihedral_types[dihedral_info]
             for dihedral_info in charmm_dihedrals

mbuild/lattice.py

@@ -523,7 +523,6 @@ def _from_lattice_vectors(self):
         return np.asarray([alpha, beta, gamma], dtype=np.float64)
 
     def _sanitize_populate_args(self, x, y, z):
-
         error_dict = {0: "X", 1: "Y", 2: "Z"}
 
         # Make sure neither x, y, z input is None

mbuild/lib/recipes/monolayer.py

@@ -71,7 +71,6 @@ def __init__(
             # Attach chains of each type to binding sites based on
             # respective fractions.
             for chain, fraction in zip(chains[:-1], fractions[:-1]):
-
                 # Create sub-pattern for this chain type
                 subpattern = deepcopy(pattern)
                 n_points = int(round(fraction * n_chains))

mbuild/lib/recipes/tiled_compound.py

@@ -86,7 +86,6 @@ def __init__(self, tile, n_tiles, name=None, **kwargs):
         # Bonds that were periodic in the original tile.
         periodic_bonds = set()
         for particle1, particle2 in tile.bonds():
-
             if np.linalg.norm(particle1.pos - particle2.pos) > dist_thresh:
                 periodic_bonds.add((particle1.index, particle2.index))
 
@@ -98,7 +97,10 @@ def __init__(self, tile, n_tiles, name=None, **kwargs):
         bonds_to_remove = set()
         bonds_to_add = set()
         for particle1, particle2 in self.bonds():
-            if (particle1.index, particle2.index,) not in periodic_bonds and (
+            if (
+                particle1.index,
+                particle2.index,
+            ) not in periodic_bonds and (
                 particle2.index,
                 particle1.index,
             ) not in periodic_bonds:

mbuild/packing.py

@@ -316,9 +316,13 @@ def fill_region(
     fix_orientation : bool or list of bools
         Specify that compounds should not be rotated when filling the box,
         default=False.
-    bounds : list-like of floats [minx, miny, minz, maxx, maxy, maxz], units nm, default=None
-        Bounding within box to pack compounds, if you want to pack within a bounding
+    bounds : list-like of list-likes of floats [[min_x, min_y, min_z, max_x, max_y, max_z], ...], units nm, default=None
+        Bounding(s) within box to pack compound(s). To pack within a bounding
         area that is not the full extent of the region, bounds are required.
+        Each item of `compound` must have its own bound specified. Use `None`
+        to indicate a given compound is not bounded, e.g.
+        `[ [0., 0., 1., 2., 2., 2.], None]` to bound only the first element
+        of `compound` and not the second.
     temp_file : str, default=None
         File name to write PACKMOL raw output to.
     update_port_locations : bool, default=False
@@ -354,6 +358,23 @@ def fill_region(
                 "`compound`, `n_compounds`, and `fix_orientation` must be of "
                 "equal length."
             )
+    if bounds is not None:
+        if not isinstance(bounds, (list)):
+            raise TypeError(
+                "`bounds` must be a list of one or more bounding boxes "
+                "and/or `None` for each item in `compound`."
+            )
+        if len(bounds) != len(n_compounds):
+            raise ValueError(
+                "`compound` and `bounds` must be of equal length. Use `None` "
+                "for non-bounded items in `compound`."
+            )
+        for bound in bounds:
+            if not isinstance(bound, (Box, list)):
+                raise ValueError(
+                    "Each bound in `bounds` must be `None`, `Box`, or a "
+                    "list of [min_x, min_y, min_z, max_x, max_y, max_z]."
+                )
 
     # See if region is a single region or list
     my_regions = []
@@ -373,7 +394,7 @@ def fill_region(
             f"expected a list of type: list or mbuild.Box, was provided {region} of type: {type(region)}"
         )
     container = []
-    if not bounds:
+    if bounds is None:
         bounds = []
     for bound, reg in zip_longest(bounds, my_regions, fillvalue=None):
         if bound is None:
@@ -720,9 +741,7 @@ def solvate(
             overlap, solvated_xyz.name, seed, sidemax * 10
         ) + PACKMOL_SOLUTE.format(solute_xyz.name, *center_solute.tolist())
 
-        for (solv, m_solvent, rotate) in zip(
-            solvent, n_solvent, fix_orientation
-        ):
+        for solv, m_solvent, rotate in zip(solvent, n_solvent, fix_orientation):
             m_solvent = int(m_solvent)
             solvent_xyz = _new_xyz_file()
             solvent_xyz_list.append(solvent_xyz)

mbuild/pattern.py

@@ -241,7 +241,6 @@ class Grid2DPattern(Pattern):
     """
 
     def __init__(self, n, m, **kwargs):
-
         points = np.zeros(shape=(n * m, 3), dtype=float)
         for i, j in product(range(n), range(m)):
             points[i * m + j, 0] = i / n

mbuild/tests/test_compound.py

@@ -1579,7 +1579,6 @@ def test_energy_minimize_fix_compounds(self, octane):
             octane.energy_minimize(fixed_compounds=[methyl_end0, (True)])
 
         with pytest.raises(Exception):
-
             octane.energy_minimize(
                 fixed_compounds=[methyl_end0, ("True", True, True)]
             )

mbuild/tests/test_packing.py

@@ -138,6 +138,33 @@ def test_fill_region_incorrect_type(self, ethane):
                 compound=[ethane], n_compounds=[2], region=box1, bounds=None
             )
 
+    def test_fill_region_bounds_not_list(self, ethane):
+        box1 = Box(lengths=[2, 2, 2], angles=[90.0, 90.0, 90.0])
+        with pytest.raises(TypeError):
+            mb.fill_region(
+                compound=[ethane], n_compounds=[2], region=box1, bounds=box1
+            )
+
+    def test_fill_region_incorrect_bounds_amount(self, ethane, h2o):
+        box1 = Box(lengths=[2, 2, 2], angles=[90.0, 90.0, 90.0])
+        with pytest.raises(ValueError):
+            mb.fill_region(
+                compound=[ethane, h2o],
+                n_compounds=[2, 2],
+                region=box1,
+                bounds=[box1],
+            )
+
+    def test_fill_region_incorrect_bounds_types(self, ethane, h2o):
+        box1 = Box(lengths=[2, 2, 2], angles=[90.0, 90.0, 90.0])
+        with pytest.raises(ValueError):
+            mb.fill_region(
+                compound=[ethane, h2o],
+                n_compounds=[2, 2],
+                region=box1,
+                bounds=[1.0, 1.0, 1.0],
+            )
+
     def test_box_no_bound(self, ethane):
         box1 = Box(lengths=[2, 2, 2], angles=[90.0, 90.0, 90.0])
         mb.fill_region(

setup.py

@@ -60,7 +60,6 @@ def compile_proto(protoc):
 
 
 if __name__ == "__main__":
-
     proto_procedure()
 
     setup(