openforcefield / openff-interchange

A project (and object) for storing, manipulating, and converting molecular mechanics data.
https://docs.openforcefield.org/projects/interchange
MIT License
71 stars 23 forks source link

Interchange subtracts tolerance when packing boxes, whereas Evaluator didn't. #1106

Open lilyminium opened 1 week ago

lilyminium commented 1 week ago

Description

Not so much a bug report as an FYI for discussion -- there are still differences between Packmol behaviour between Evaluator and Interchange, and I believe one of them is down to the line below. Evaluator's input file does not seem to have subtracted the tolerance:

Reproduction

Please include a minimally reproducing example of this bug.

https://github.com/openforcefield/openff-interchange/blob/8d2cee156c8a6763fde0e1938809a96b688325da/openff/interchange/components/_packmol.py#L477

Output

Using a mixture of triethanolamine and water, modern Evaluator builds the following input file:

Evaluator packmol_input.txt

tolerance 2.000000
filetype pdb
output packmol_output.pdb

structure 0.pdb
  number 490
  inside box 0. 0. 0. 57.63472087723553 57.63472087723553 57.63472087723553
end structure

structure 1.pdb
  number 510
  inside box 0. 0. 0. 57.63472087723553 57.63472087723553 57.63472087723553
end structure

Interchange 0.4.0:

tolerance 2.000000
filetype pdb
output packmol_output.pdb

structure _PACKING_MOLECULE0.pdb
  number 490
  inside box 0. 0. 0. 55.63472087723553 55.63472087723553 55.63472087723553
end structure

structure _PACKING_MOLECULE1.pdb
  number 510
  inside box 0. 0. 0. 55.63472087723553 55.63472087723553 55.63472087723553
end structure

And the box size from density from both returns ~57.63472087723553:

Evaluator:

from openff.evaluator.utils.packmol import _approximate_box_size_by_density
from openff.toolkit.topology.molecule import Molecule, unit
from openff.evaluator import unit as eval_unit
from openff.evaluator.utils.openmm import openmm_quantity_to_pint
import openff.evaluator
import openff.toolkit

smiles = "OCCN(CCO)CCO"

comp0 = Molecule.from_smiles(smiles)
comp1 = Molecule.from_smiles("O")
copies = [490, 510]
mass_density = 0.95 * eval_unit.grams / eval_unit.milliliters
box_aspect_ratio = [1, 1, 1]

box = _approximate_box_size_by_density(
    [comp0, comp1],
    copies,
    mass_density,
    box_aspect_ratio,
)

#   [57.63472087723553 57.63472087723553 57.63472087723553]

Interchange:

from openff.interchange.components._packmol import _box_from_density
from openff.toolkit.topology.molecule import Molecule, unit
import openff.interchange
import openff.toolkit
import numpy as np

box = _box_from_density(
    [comp0, comp1],
    copies,
    mass_density,
    np.identity(3),
)

# [[57.63472087723553 0.0 0.0]
# [0.0 57.63472087723553 0.0]
# [0.0 0.0 57.63472087723553]]

In a fun twist, I'm not getting Evaluator's input box to actually pack for me either in this separate workflow (it has worked fine in an actual run of Evaluator....) so there's also other stuff going on here, but thought I'd raise it.

Software versions

mattwthompson commented 1 week ago

cc @Yoshanuikabundi - I'm happy for this change to be included if you two agree it makes sense

lilyminium commented 1 week ago

I have a weak preference for putting this back in, mostly because I'm not sure it makes sense to subtract the tolerance from the initial box. In the "working" Evaluator output, the box does use the tolerance and goes up to 59 A. It's only a weak preference because you can get around it by either changing the density or scaleup factor (although the scaleup isn't currently a publically available kwarg).

lilyminium commented 1 week ago

FYI, in using only 1 conformer in #1107 worked for all the test cases I was looking at without removing this additional subtraction of tolerance, so my "weak" preference here is definitely very weak.