Different behavior in `fill_box` and `fill_region`. `fill_region` errors when using a list of Compounds
chrisjonesBSU opened this issue · 5 comments
Bug summary
the fill_region function in packing.py isn't working as expected, specifically when using a list of compounds and a list of n_compounds
Code to reproduce the behavior
import mbuild as mb
methane = mb.load("C", smiles=True)
ethane = mb.load("CC", smiles=True)
# This works with `fill_box`
box_compound = mb.fill_box(
compound=[methane, ethane], n_compounds=[10, 10], box=[2,2,2]
)
# This fails with `fill_region`
box = mb.box.Box(lengths=[2, 2, 2], angles=[90, 90, 90])
region_compound = mb.fill_region(
compound=[methane, ethane],
n_compounds=[10, 10],
region=box,
bounds = [[0, 0, 1, 2, 2, 2]]
)I get the following error when running the mb.fill_region portion:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[3], line 5
2 ethane = mb.load("CC", smiles=True)
4 box = mb.box.Box(lengths=[2, 2, 2], angles=[90, 90, 90])
----> 5 region_compound = mb.fill_region(
6 compound=[methane, ethane],
7 n_compounds=[10, 10],
8 region=box,
9 bounds = [[0, 0, 1, 2, 2, 2]]
10 )
File ~/miniconda3/envs/test-mbuild/lib/python3.11/site-packages/mbuild/packing.py:431, in fill_region(compound, n_compounds, region, overlap, bounds, seed, sidemax, edge, fix_orientation, temp_file, update_port_locations)
429 filled = Compound()
430 filled = _create_topology(filled, compound, n_compounds)
--> 431 filled.update_coordinates(
432 filled_xyz.name, update_port_locations=update_port_locations
433 )
434 finally:
435 for file_handle in compound_xyz_list:
File ~/miniconda3/envs/test-mbuild/lib/python3.11/site-packages/mbuild/compound.py:1818, in Compound.update_coordinates(self, filename, update_port_locations)
1816 self._update_port_locations(xyz_init)
1817 else:
-> 1818 self = conversion.load(filename, compound=self, coords_only=True)
File ~/miniconda3/envs/test-mbuild/lib/python3.11/site-packages/mbuild/conversion.py:127, in load(filename_or_object, relative_to_module, compound, coords_only, rigid, smiles, infer_hierarchy, backend, ignore_box_warn, **kwargs)
117 return load_smiles(
118 smiles_or_filename=filename_or_object,
119 compound=compound,
(...)
123 **kwargs,
124 )
125 # Last, if none of the above, load from file
126 else:
--> 127 return load_file(
128 filename=filename_or_object,
129 relative_to_module=relative_to_module,
130 compound=compound,
131 coords_only=coords_only,
132 rigid=rigid,
133 backend=backend,
134 infer_hierarchy=infer_hierarchy,
135 **kwargs,
136 )
File ~/miniconda3/envs/test-mbuild/lib/python3.11/site-packages/mbuild/conversion.py:412, in load_file(filename, relative_to_module, compound, coords_only, rigid, backend, infer_hierarchy, **kwargs)
410 tmp = read_xyz(filename)
411 if tmp.n_particles != compound.n_particles:
--> 412 raise ValueError(
413 "Number of atoms in {filename}"
414 "does not match {compound}".format(**locals())
415 )
416 ref_and_compound = zip(
417 tmp._particles(include_ports=False),
418 compound.particles(include_ports=False),
419 )
420 for ref_particle, particle in ref_and_compound:
ValueError: Number of atoms in /tmp/tmpn75wd68o.xyzdoes not match <Compound 130 particles, 110 bonds, non-periodic, id: 140559818198608>
fill_region does work when compound and n_compounds aren't passed in as lists larger than length 1:
methane = mb.load("C", smiles=True)
box = mb.box.Box(lengths=[2, 2, 2], angles=[90, 90, 90])
region_compound = mb.fill_region(
compound=[methane],
n_compounds=[20],
region=box,
bounds = [[0, 0, 1, 2, 2, 2]]
)
Software versions
- Which version of mBuild are you using? 0.15.1
- Which version of Python (
python --version)? 3.11 - Which operating system? Linux
Following!
Ok, turns out this wasn't really a bug in the functionality, but poor documentation and a missing check/validation step.
397 for comp, m_compounds, rotate, items_n in zip(
398 compound, n_compounds, fix_orientation, container
399 ):
In this for loop/zip, if container is only of length 1 while compound and n_compounds are of length loner than 1, it only goes through this for loop once (or as long as the smallest length item in zip)
I think we'll have to handle bounds in the same way we do fix_orientation
343 if not isinstance(fix_orientation, (list, set)):
344 fix_orientation = [fix_orientation] * len(compound)
In the example of my original comment, I only have 1 boundary condition, that I want to fill with both methane and ethane (i.e. fill the region above a surface), so I would have to do something like:
import mbuild as mb
methane = mb.load("C", smiles=True)
ethane = mb.load("CC", smiles=True)
pack_box = mb.box.Box(lengths=[2, 2, 2], angles=[90, 90, 90])
fill_bounds = [0, 0, 1, 2, 2, 2]
region_compound = mb.fill_region(
compound=[methane, ethane],
n_compounds=[10, 10],
region=pack_box,
bounds=[[fill_bounds, fill_bounds]]
)
I think it would make sense to be able to use a single boundary condition that can be used with any length of compound and n_compounds, while still having the functionality to define boundaries on a molecule-by-molecule basis, which is basically how fix_orientation works.
I'll work on a PR.
Although, we could keep it as is, and require that the user is explicit in defining packing boundaries (like my example above), then we just need a check to make sure bounds is the same length as compound and n_compounds. It is convenient to only have to define the boundary condition once if you want to pack it with a mixture, but also, requiring the user to be more explicit/thoughtful in what parameters they are passing in helps avoid confusion and assumptions that cause issues later.
Good find! Yeah, I agree a smitch more documentation + better validation (so more meaningful error) would be great!
Closed by #1088