Clarify default args for box in Compound.save()

[uppittu11]

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:

mbuild/mbuild/compound.py

Lines 1839 to 1842 in fe5292e

	box : mb.Box, optional, default=self.boundingbox (with buffer)
	Box information to be written to the output file. If 'None', a
	bounding box is used with 0.25nm buffers at each face to avoid
	overlapping atoms.

Where the actual box dims come from:

mbuild/mbuild/compound.py

Line 1929 in fe5292e

structure = self.to_parmed(box=box, residues=residues,

mbuild/mbuild/compound.py

Lines 2442 to 2454 in fe5292e

	if box is None:
	box = self.boundingbox
	box_vec_max = box.maxs.tolist()
	box_vec_min = box.mins.tolist()
	for dim, val in enumerate(self.periodicity):
	if val:
	box_vec_max[dim] = val
	box_vec_min[dim] = 0.0
	if not val:
	box_vec_max[dim] += 0.25
	box_vec_min[dim] -= 0.25
	box.mins = np.asarray(box_vec_min)
	box.maxs = np.asarray(box_vec_max)

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.

[mattwthompson]

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

mbuild/mbuild/compound.py

Lines 2327 to 2332 in fe5292e

	box : mb.Box, optional, default=self.boundingbox (with buffer)
	Box information to be used when converting to a `Structure`.
	If 'None', a bounding box is used with 0.25nm buffers at
	each face to avoid overlapping atoms, unless `self.periodicity`
	is not None, in which case those values are used for the
	box lengths.

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:

mbuild/mbuild/compound.py

Lines 759 to 761 in fe5292e

	if (inherit_periodicity and isinstance(new_child, Compound) and
	new_child.periodicity.any()):
	self.periodicity = new_child.periodicity

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.

[CalCraven]

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)