Specify individual axes to constrain within PACKMOL in packing.py

[chrisjonesBSU]

PR Summary:

Addresses issue #725 and allows user to fix rotations about axes individually rather than all or none.

• Changes the fix_orientation parameter from True/False to a tuple where True/False is given for each axis (x, y, z). e.g. fix_orientation = (False, False, False)

• Added a new function packmol_constrain that builds up the PACKMOL_CONSTRAIN string variable depending on the contents of fix_orientation

I added doc strings to packmol_constrain, but haven't yet made the doc string updates to fill_box that reflect this change. I also have not created the necessary unit tests, or looked into what unit test already exists fill_box and fix_orientation

Ultimately, this gives PACKMOL and therefore the user a little more flexibility. Specifically, it worked for me when packing a layer of organic molecules above a crystal layer where I wanted to keep the molecules parallel with the crystal, but allow them to make planar rotations during packing.

PR Checklist

---

• Includes appropriate unit test(s)
• [] Appropriate docstring(s) are added/updated
• Code is (approximately) PEP8 compliant
• Issue(s) raised/addressed?

[chrisjonesBSU]

Also, I think it is worth mentioning that after talking with Jenny a bit about this change, she offered an alternative approach that greatly reduced the amount of if/else statements. I decided to keep my original changes because its a bit easier to read what is happening, but it is probably worth considering a more efficient approach given below:

    constrain = 'constrain_rotation {} 0. 0.\n'
    dirs = ['x','y','z']
    format_list = [constrain.format(xyz) if tf else "" for xyz,tf in zip(dirs,fix_orientation)]
    return "".join(format_list)

[justinGilmer]

This is a great addition! Thanks @chrisjonesBSU ! The next steps are the unit tests. Im not entirely sure the best way to validate this. Aside from asserting that the expected strings are written out for a variety of cases, and packmol successfully fills the box. A minimal example of a system with this applied might be useful to post here too. But that would be more for GitHub viewing rather than something we can test.

I also included a possible suggestion for the logic to reduce some of the needed conditionals, but we can choose however we want to do that as this is fleshed out with tests.

[justinGilmer]

One other option that might be aligned with jenny's suggestion

final_constraints = list()
for i,val in enumerate(fix_orientation):
  if val:
    final_constraints.append(constraints[i])
  else:
    final_constraints.append("")
PACKMOL_CONSTRAIN = CONSTRAIN.format(final_constraints**)

Yours is more readable, but this removes all of the conditionals from line 74 onwards i think

[rmatsum836]

Right now test_no_rotate and test_rotate in test_packing.py are failing.

[chrisjonesBSU]

@justinGilmer Since this involves a slight API change, should I add behaviour to handle the original case of fix_orientation not being iterable? If the user passes fix_orientation=True then it will default to True for all directions (and same with False). @rmatsum836 suggested we could do this and include a warning message about the change; just curious on what others think.

[rsdefever]

Since this involves a slight API change, should I add behaviour to handle the original case of fix_orientation not being iterable? If the user passes fix_orientation=True then it will default to True for all directions (and same with False)

I think this should still be supported. If the user passes True or False it should apply to all compounds and all rotation directions.

Also, double check that your code will give the expected behavior for multiple compounds. There is an isinstance(fix_orientation, (list, set)) line that I think might have been broken by your changes. I could be mistaken though.

[lgtm-com]

This pull request introduces 2 alerts and fixes 1 when merging 7f068d8 into 4b60ced - view on LGTM.com

new alerts:

• 1 for Syntax error
• 1 for 'import *' may pollute namespace

fixed alerts:

• 1 for Variable defined multiple times

[lgtm-com]

This pull request introduces 1 alert when merging 41669bc into 4b60ced - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[codecov]

Codecov Report

Merging #734 (c363358) into master (16c39b0) will decrease coverage by 0.20%.
The diff coverage is 70.21%.

@@            Coverage Diff             @@
##           master     #734      +/-   ##
==========================================
- Coverage   88.82%   88.62%   -0.21%     
==========================================
  Files          59       59              
  Lines        5156     5197      +41     
==========================================
+ Hits         4580     4606      +26     
- Misses        576      591      +15     

Impacted Files	Coverage Δ
mbuild/packing.py	90.85% <70.21%> (-3.92%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 16c39b0...c363358. Read the comment docs.

[lgtm-com]

This pull request introduces 1 alert when merging 033b576 into 2e00474 - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[lgtm-com]

This pull request introduces 1 alert when merging 3ba2b80 into 2e00474 - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[chrisjonesBSU]

I think this is ready. Per a conversation with @daico007 I left the original test_rotate and test_no_rotate which will still test for fix_orientation = (False, False, False) or (True, True, True).

The new unit test specify_axis tests for the 6 other possible combinations of fix_orientation. It's important to note that they only test that the correct PACKMOL constraint string is generated rather than going as far as building a small system using PACKMOL and checking that the atomic positions are consistent with the constraints generated.

[rmatsum836]

This looks great, I have just a couple small notes shown below. Thanks for the contribution!

[lgtm-com]

This pull request introduces 4 alerts when merging 475d17f into 2e1c1c9 - view on LGTM.com

new alerts:

• 4 for Except block handles 'BaseException'

[daico007]

LGTM!

[chrisjonesBSU]

Since the packmol_constrain function is strictly a behind-the-scenes function, should I rename it to _packmol_constrain and move it below the user-facing functions?

[lgtm-com]

This pull request introduces 4 alerts when merging 3b6f19a into 368236d - view on LGTM.com

new alerts:

• 4 for Except block handles 'BaseException'

[uppittu11]

This is awesome! Thanks for adding this @chrisjonesBSU. Just had a few small changes.

[uppittu11]

Ok I found a couple more cases that are unexpected. It's a little unclear what the accepted types and shapes of the fix_orientation arg are. It would be good if you add some more type checking and clarifying what inputs are allowed.

Here is a gist with some commands I expected to work: https://gist.github.com/uppittu11/f2acb70e12bff6ae249775b7bd24e6f8

[lgtm-com]

This pull request introduces 4 alerts when merging 85cf398 into 8a28926 - view on LGTM.com

new alerts:

• 4 for Except block handles 'BaseException'

[chrisjonesBSU]

Ok I found a couple more cases that are unexpected. It's a little unclear what the accepted types and shapes of the fix_orientation arg are. It would be good if you add some more type checking and clarifying what inputs are allowed.

Here is a gist with some commands I expected to work: https://gist.github.com/uppittu11/f2acb70e12bff6ae249775b7bd24e6f8

Thanks for the input! I pushed a change that fixed the typos and small errors, I'll work on addressing this next.

[chrisjonesBSU]

Ok I found a couple more cases that are unexpected. It's a little unclear what the accepted types and shapes of the fix_orientation arg are. It would be good if you add some more type checking and clarifying what inputs are allowed.

Here is a gist with some commands I expected to work: https://gist.github.com/uppittu11/f2acb70e12bff6ae249775b7bd24e6f8

With case 1, I agree there is an issue with passing fix_orientation as a list of bools of length 3.

I would expect fix_orientation = [True, True, True], fix_orientation = (True, True, True) and fix_orientation = True all to behave the same regardless of the length of compound

This most recent commit includes a fix to your scenario in case 1. If fix_orientation is a list of booleans of length 3, then it applies those constraints to the [x,y,z] axes for everything in compound

However, it gets tricky with something like fix_orientation = [False, True, True] and compound is a list of length 3. It’s impossible to know if the user wants fix_orientation = [[False, False, False], [True, True, True], [True, True, True]] (which would be the current behavior of mbuild), or if they want fix_orientation = [False, True, True]*3 (which is the current behavior via the changes in the most recent commit). This is only an issue when the length of compound is 3. I feel like we need to decide what the default behavior will be in this case, and if it is not what the user intended, leave it to them to read the doc strings and pass in a more detailed representation of what they want for fix_orientation. Given the changes of this PR, I feel like fix_orientation = [False, True, True] should be treated as a single argument passed rather than a list of 3 separate arguments.

Case 2 is not properly addressed in the doc strings, but is in the warning message that pops up when a single bool is passed. The thinking here was that the old behavior of passing a single bool to specify all 3 axes would eventually be deprecated, but it was never decided if that would be the case or not.. I changed the doc strings to clarify what is allowed and what the resulting behavior is.

Case 3, I personally wouldn't expect a numpy array of bools to work. To me, it is an awkward way of passing in this kind of parameter. Errors also result if compound is an array of mb.Compounds or even if n_compounds is an array of ints.

n = 10
ethanes = [mb.load('CC', smiles=True) for _ in range(n)]
n_compounds = np.full((n, 1), 100)
system = mb.fill_box(compound = ethanes,
                    n_compounds=n_compounds,
                    box=mb.Box(lengths=[5, 5, 5]),
                    fix_orientation=False)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
--> 117     arg_count = 3 - [n_compounds, box, density].count(None)
    118     if arg_count != 2:
    119         msg = ("Exactly 2 of `n_compounds`, `box`, and `density` "

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

And

n=10
ethane = mb.load('CC', smiles=True)
compounds = np.full((n, 1), ethane)
n_compounds = [100]*n
system = mb.fill_box(compound = compounds,
                    n_compounds=n_compounds,
                    box=mb.Box(lengths=[5, 5, 5]),
                    fix_orientation=False)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
    141         if len(compound) != len(n_compounds):
    142             msg = ("`compound` and `n_compounds` must be of equal length.")
--> 143             raise ValueError(msg)
    144
    145     if compound is not None:

ValueError: `compound` and `n_compounds` must be of equal length.

Perhaps this case could be addressed in a separate issue/PR that both condenses and makes more robust all of the type checking being done in the 4 different packing functions. I know that Ryan started something like this in a PR that was eventually dropped. He created a single function to handle all of the type checking that was being called in each of the separate packing functions. He did mention to me that I include that approach in this PR, but I think it would make more sense as a separate PR.

[lgtm-com]

This pull request introduces 4 alerts when merging d63f967 into ad61f7c - view on LGTM.com

new alerts:

• 4 for Except block handles 'BaseException'

[lgtm-com]

This pull request introduces 4 alerts when merging 052be47 into f2cd2a8 - view on LGTM.com

new alerts:

• 4 for Except block handles 'BaseException'

[lgtm-com]

This pull request introduces 1 alert and fixes 1 when merging 59cf662 into 636f9b1 - view on LGTM.com

new alerts:

• 1 for Syntax error

fixed alerts:

• 1 for Variable defined multiple times

[jennyfothergill]

tests are failing. Let's check if these two changes fix them and I'll review again c:

[jennyfothergill]

looks like the tests are passing! 👍 we can increase codecoverage by asserting that the TypeError is raised if fix_orientation is passed in as a noniterable. Then I think this is ready to go :)

[CalCraven]

How are we looking on this PR? Need any final help/input to gett it merged?

[chrisjonesBSU]

It's probably a bit outdated; I should have time over the next couple weeks to take another look at it. Ultimately, I still think this feature should be added, but there are some quirks with how it's handled.

[CalCraven]

Sounds good. Let me know if you need any help trying anything out. Seems like a lot of issues users have with mBuild are trying to do something complicated with PackMol, enough so that someone suggested we look into methods to replace packmol with a more accessible and efficient code. At any rate, having better support for packmol features is really good for mBuild.