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

PACKMOL gets grouchy if there are zero molecules of a species to insert #769

Open
rsdefever opened this issue Jul 23, 2020 · 1 comment
Open
Assignees

Comments

@rsdefever
Copy link
Member

Bug summary
If you try to insert zero molecules of a species, PACKMOL is grumpy. Example:

    # Use 1000 total ions
    n_ions = 1000
    x_li = job.sp.conc

    # Taken from Wang et al PIM results
    # at 1 bar and 1100 K
    licl_density=1431 #kg/m^3
    kcl_density=1497 #kg/m^3

    # Compute mixture density with ideal mixing
    mixed_density = x_li * licl_density + (1.0-x_li) * kcl_density

    # Compute the number of ions of each type
    n_cl = int(0.5 * n_ions)
    n_li = int(0.5 * x_li * n_ions)
    n_k = n_ions - n_cl - n_li

    # Create the system
    system = mbuild.fill_box(
        [k, li, cl],
         n_compounds=[n_k, n_li, n_cl],
         density=mixed_density
     )

Why might you want to do this? Imagine spanning a concentration range of 0.0 to 1.0 to, e.g., compute the enthalpy of mixing. Since PACKMOL won't let you specify zero molecules of a species, the above section of pretty code turns ugly with if and elif conditions to handle the 0.0 and 1.0 limits.

** Fix **

I propose we add a snippet of code to packing.py that checks if any species have zero molecules and then cleans that up before calling PACKMOL.

@rsdefever rsdefever self-assigned this Jul 23, 2020
@mattwthompson
Copy link
Member

Definitely a bug, and I'm surprised it didn't pop up in the work Ray and I did.

For whoever wants to tackle this, the patch would take lines like this and add something like

        for comp, m_compounds, rotate in zip(compound, n_compounds, fix_orientation):
            m_compounds = int(m_compounds)
            if not m_compounds:
                continue

The search string = int(m in packing.py finds the four instances that come to mind. It should be straightforward enough (except whenever that is said, it's wrong). Just make sure to add some unit tests that catch this behavior

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

2 participants