Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Different behavior in fill_box and fill_region. fill_region errors when using a list of Compounds #1080

Closed
chrisjonesBSU opened this issue Jan 25, 2023 · 5 comments

Comments

@chrisjonesBSU
Copy link
Contributor

chrisjonesBSU commented Jan 25, 2023

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
@RainierBarrett
Copy link
Contributor

Following!

@chrisjonesBSU
Copy link
Contributor Author

chrisjonesBSU commented Jan 31, 2023

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.

@chrisjonesBSU
Copy link
Contributor Author

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.

@daico007
Copy link
Member

daico007 commented Feb 1, 2023

Good find! Yeah, I agree a smitch more documentation + better validation (so more meaningful error) would be great!

@RainierBarrett RainierBarrett mentioned this issue Feb 3, 2023
4 tasks
@chrisjonesBSU
Copy link
Contributor Author

Closed by #1088

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants