mbuild icon indicating copy to clipboard operation
mbuild copied to clipboard

Clarify default args for box in Compound.save()

Open uppittu11 opened this issue 5 years ago • 2 comments

Bug summary

Docstring says that the box arg defaults to self.boundingbox, but defaults to periodicity IRL. We should either use the bounding box or track the periodicity better when adding compounds.

Docstring: https://github.com/mosdef-hub/mbuild/blob/fe5292e45343888ce223a9b78e16593c9ddc1317/mbuild/compound.py#L1839-L1842

Where the actual box dims come from: https://github.com/mosdef-hub/mbuild/blob/fe5292e45343888ce223a9b78e16593c9ddc1317/mbuild/compound.py#L1929

https://github.com/mosdef-hub/mbuild/blob/fe5292e45343888ce223a9b78e16593c9ddc1317/mbuild/compound.py#L2442-L2454

Code to reproduce the behavior

import mbuild

propane = mbuild.lib.recipes.Alkane(3)

box1 = mbuild.Box(mins=[0, 0, 0], maxs=[2, 2, 2])
box2 = mbuild.Box(mins=[0, 0, 3], maxs=[2, 2, 5])

filled_box1 = mbuild.fill_box(propane, 100, box=box1)
filled_box2 = mbuild.fill_box(propane, 100, box=box2)

system = mbuild.Compound()
system.add(filled_box1)
system.add(filled_box2)

print("System bounding box: ", system.boundingbox.lengths)

system.save("file.gro", residues=['Alkane'], overwrite=True)

import mdtraj
system_reloaded = mdtraj.load("file.gro")
print("Actual box lengths: ", system_reloaded.unitcell_lengths)

Output:

System bounding box:  [1.801899  1.801009  4.8007786]
Actual box lengths:  [[2. 2. 2.]]

I'd expect the actual box lengths to be equal to the bounding box, but that's not the case.

uppittu11 avatar Jan 29 '20 16:01 uppittu11

It looks like there's a discrepancy between the box docstrings for save and to_parmed (here):

https://github.com/mosdef-hub/mbuild/blob/fe5292e45343888ce223a9b78e16593c9ddc1317/mbuild/compound.py#L2327-L2332

I think this is the right behavior (bounding boxes are bad and periodicity should be used if it exists). Exactly how the periodicity got set is more opaque, it looks like the packing functions set the periodicity of the returned compounds to the dimensions of the boxes you passed in, and adding compounds into that container just took in the values of the first thing added. But after adding a second compound, the periodicity wasn't updated. It looks like there's an argument to add that allows you to set the periodicity to what's being added, but not intelligently update it after the fact:

https://github.com/mosdef-hub/mbuild/blob/fe5292e45343888ce223a9b78e16593c9ddc1317/mbuild/compound.py#L759-L761

You'd be able to use either box1 or box2, but not an updated version of each. So the root of the problem here, I think, is that the periodicity is garbage after adding things that cover different regions of space. Should add do some cleaning up to check that every particle is inside the self.periodicity? I think the "check for periodicity, if it's not there then fall back to boundingbox with a buffer" is generally the right behavior.

mattwthompson avatar Jan 29 '20 16:01 mattwthompson

This is still an issue, and leads to some ambiguity in the writers. This line shows that even if a box is not passed, the box of the compound can still be grabbed from compound.box. This should be a quick fix to the docs.

    if box is None:
        if compound.box is not None:
            box = deepcopy(compound.box)
        else:
            box = compound.get_boundingbox()
            # Pad by an extra 0.5 nm (0.25 on each side) from bounding box
            box = Box(lengths=np.array(box.lengths) + 0.5, angles=box.angles)

CalCraven avatar Aug 10 '22 21:08 CalCraven