mbuild
mbuild copied to clipboard
Clarify default args for box in Compound.save()
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.
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.
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)