Packmol wrapper chokes on amino acids because of MDTraj-written PDB files

[lilyminium]

Describe the bug

If a property includes a molecule that is a recognisable amino acid, this will get named as such in Packmol. When that happens, the parameter assignment step fails at the Topology.from_openmm call. This is in the Equilibration workflow. Highly related to #691.

To Reproduce

Try to equilibrate the property in the attached CSV with this script.

issue-699.csv

Output

JSON file:

0cfeef75228241b0bfe095a1798139dd/dens_2059314443383903998_assign_parameters/dens_2059314443383903998|assign_parameters_output.json:{"message": ["Traceback (most recent call last):\n", "  File \"/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/evaluator/workflow/protocols.py\", line 1201, in _execute_protocol\n    protocol.execute(directory, available_resources)\n", "  File \"/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/evaluator/workflow/protocols.py\", line 688, in execute\n    self._execute(directory, available_resources)\n", "  File \"/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/evaluator/protocols/forcefield.py\", line 563, in _execute\n    topology = Topology.from_openmm(\n               ^^^^^^^^^^^^^^^^^^^^^\n", "  File \"/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/utilities/utilities.py\", line 81, in wrapper\n    return function(*args, **kwargs)\n           ^^^^^^^^^^^^^^^^^^^^^^^^^\n", "  File \"/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/toolkit/topology/topology.py\", line 1495, in from_openmm\n    raise ValueError(msg)\n", "ValueError: No match found for molecule C. This would be a very unusual molecule to try and parameterize, and it is likely that the data source it was read from does not contain connectivity information. If this molecule is coming from PDB, please ensure that the file contains CONECT records. The PDB format documentation (https://www.wwpdb.org/documentation/file-format-content/format33/sect10.html) states \"CONECT records are mandatory for HET groups (excluding water) and for other bonds not specified in the standard residue connectivity table.\"\n"], "protocol_id": "dens_2059314443383903998|assign_parameters", "@type": "openff.evaluator.workflow.exceptions.WorkflowException"}

Printed with formatting:

Traceback (most recent call last):
  File "/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/evaluator/workflow/protocols.py", line 1201, in _execute_protocol
    protocol.execute(directory, available_resources)
  File "/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/evaluator/workflow/protocols.py", line 688, in execute
    self._execute(directory, available_resources)
  File "/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/evaluator/protocols/forcefield.py", line 563, in _execute
    topology = Topology.from_openmm(
               ^^^^^^^^^^^^^^^^^^^^^
  File "/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/utilities/utilities.py", line 81, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/homezvol3/lilyw7/miniforge3/envs/evaluator-050/lib/python3.11/site-packages/openff/toolkit/topology/topology.py", line 1495, in from_openmm
    raise ValueError(msg)
ValueError: No match found for molecule C. This would be a very unusual molecule to try and parameterize, and it is likely that the data source it was read from does not contain connectivity information. If this molecule is coming from PDB, please ensure that the file contains CONECT records. The PDB format documentation (https://www.wwpdb.org/documentation/file-format-content/format33/sect10.html) states "CONECT records are mandatory for HET groups (excluding water) and for other bonds not specified in the standard residue connectivity table."

Computing environment (please complete the following information):

• Operating system
• Output of running conda list

Additional context

Culprit PDB (output.pdb from Packmol):

issue-699.pdb

[mattwthompson]

I agree with the assessment; alanine here isn't getting CONECT records:

$ grep TER issue-699.pdb
TER    1015      ALA A  78
TER   12080      AVL B 922
$ grep CONECT issue-699.pdb | head -n 10
CONECT 1016 1017 1021
CONECT 1017 1016 1018 1022 1023
CONECT 1018 1017 1019 1024 1025
CONECT 1019 1018 1020 1026 1027
CONECT 1020 1019
CONECT 1021 1016
CONECT 1022 1017
CONECT 1023 1017
CONECT 1024 1018
CONECT 1025 1018

I haven't yet figured out if it's a problem with Packmol, with how we're preparing molecules for it, or processing molecules after it. I suspect Evaluator has mostly not been used at all with amino acids?

---

I tried to shift this to a more targeted reproduction:

from openff.evaluator.datasets import PhysicalPropertyDataSet
import pandas
import openmm.app

from openff.toolkit import Topology, Molecule


unique_molecules = [
    Molecule.from_smiles(component.smiles)
    for component in [
        *PhysicalPropertyDataSet.from_pandas(
            pandas.read_csv("issue-699.csv")
        ).substances
    ][0].components
]

openmm_topology = openmm.app.PDBFile("issue-699.pdb").topology

try:
    Topology.from_openmm(
        openmm_topology,
        unique_molecules=unique_molecules,
    )
except ValueError as error:
    assert "No match found for molecule C" in str(error)

(Below the next divider is me going down the wrong path)

---

Topology.from_pdb fails, however, because the alanine doesn't have bonds:

In [5]: [*openmm_topology.bonds()][0]
Out[5]: Bond(<Atom 1014 (O1) of chain 1 residue 78 (AVL)>, <Atom 1015 (C1) of chain 1 residue 78 (AVL)>)

The atom ordering of the alanine in the PDB file doesn't match what Molecule.from_smiles guesses:

In [2]: [atom.element.symbol for atom in next(openmm_topology.residues()).atoms()]
Out[2]: ['C', 'C', 'H', 'N', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H']

In [3]: [atom.symbol for atom in unique_molecules[0].atoms]
Out[3]: ['C', 'C', 'N', 'C', 'O', 'O', 'H', 'H', 'H', 'H', 'H', 'H', 'H']

but that shouldn't matter since it should be doing isomorphism checks.

Still, this is a little surprising - I would guess that the atom ordering sent to Packmol (and its wrapper) are generated from the same code?

---

[mattwthompson]

I can also rule out RDKitToolkitWrapper not writing bonds:

In [3]:
   ...: for index, molecule in enumerate(unique_molecules):
   ...:     molecule.to_file(
   ...:         f"{index}.pdb", file_format="pdb",
   ...:         toolkit_registry=RDKitToolkitWrapper(),
   ...:     )
   ...:

In [4]: !cat 0.pdb 1.pdb
HETATM    1  C1  UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    2  C2  UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    3  N1  UNL     1       0.000   0.000   0.000  1.00  0.00           N
HETATM    4  C3  UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    5  O1  UNL     1       0.000   0.000   0.000  1.00  0.00           O
HETATM    6  O2  UNL     1       0.000   0.000   0.000  1.00  0.00           O
HETATM    7  H1  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM    8  H2  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM    9  H3  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   10  H4  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   11  H5  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   12  H6  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   13  H7  UNL     1       0.000   0.000   0.000  1.00  0.00           H
CONECT    1    2    7    8    9
CONECT    2    3    4   10
CONECT    3   11   12
CONECT    4    5    5    6
CONECT    6   13
END
HETATM    1  O1  UNL     1       0.000   0.000   0.000  1.00  0.00           O
HETATM    2  C1  UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    3  C2  UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    4  C3  UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    5 CL1  UNL     1       0.000   0.000   0.000  1.00  0.00          CL
HETATM    6  H1  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM    7  H2  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM    8  H3  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM    9  H4  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   10  H5  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   11  H6  UNL     1       0.000   0.000   0.000  1.00  0.00           H
HETATM   12  H7  UNL     1       0.000   0.000   0.000  1.00  0.00           H
CONECT    1    2    6
CONECT    2    3    7    8
CONECT    3    4    9   10
CONECT    4    5   11   12
END

[mattwthompson]

I had a closer look at how amino acid residue names are handled in Evaluator's Packmol wrapper.

Residues are used for molecule bookkeepping, which necessitates generating a random residue name on-the-fly for each residue. When this happens to generate a amino acid code¹ a new residue name is generated until it doesn't look like an amino acid code. So if I'm working with hexanol in water and something's randomly given GLY, the dice is rolled again so that it isn't.

When a molecule does look like an amino acid, its residue name is intentionally set to be the matching amino acid code. So when a solute is an amino acid, it is kept to look like one. Is this the behavior we want? It's unchanged for some 5 years and I don't know if amino acids are in existing training sets.

Footnotes

1. Or one of something else that has a pre-defined meaning, in a small pre-defined set ↩

[mattwthompson]

Evaluator uses MDTraj to write PDB files, which lack CONECT records of amino acids, so Packmol writes a PDB missing the same CONECT records, breaking Topology.from_openmm.

For the short version, just trust me on the following details¹ (each of which I have reproductions/lengthier explanations for if sought)

• Packmol will write (or not write) CONECT records based on the contents of the inputs
• OpenMM does its best to read what's in a PDB file, but won't fill in bonds when CONECT records are missing²
• The error bubbling up from openff.toolkit.Topology.from_openmm is because of missing CONECT records
  • The toolkit will then complain³ inside of Topology.from_pdb
• For molecules that don't look at all like amino acids, etc., MDTraj will write bonds as CONECT records
• Evaluator internally standardizes SMILES, but this doesn't make a difference here

In [1]: from openff.evaluator.utils.packmol import _create_trajectory

In [2]: from openff.toolkit import Molecule

In [3]: ala_raw_openff = Molecule.from_smiles("C[C@@H](C(=O)O)N")

In [4]: _create_trajectory(ala_raw_openff).save("x.pdb")

In [5]: cat x.pdb
ATOM      1  C1  ALA     1      -1.019   0.656  -0.494  1.00  0.00           C
ATOM      2  C2  ALA     1      -0.153  -0.474  -0.029  1.00  0.00           C
ATOM      3  C3  ALA     1       1.282  -0.046  -0.022  1.00  0.00           C
ATOM      4  O1  ALA     1       1.953  -0.154   1.056  1.00  0.00           O
ATOM      5  O2  ALA     1       1.884   0.463  -1.169  1.00  0.00           O
ATOM      6  N1  ALA     1      -0.512  -0.954   1.271  1.00  0.00           N
ATOM      7  H1  ALA     1      -0.632   1.632  -0.147  1.00  0.00           H
ATOM      8  H2  ALA     1      -2.042   0.576  -0.050  1.00  0.00           H
ATOM      9  H3  ALA     1      -1.135   0.627  -1.600  1.00  0.00           H
ATOM     10  H4  ALA     1      -0.238  -1.278  -0.791  1.00  0.00           H
ATOM     11  H5  ALA     1       2.761   0.990  -1.125  1.00  0.00           H
ATOM     12  H6  ALA     1      -1.009  -0.236   1.865  1.00  0.00           H
ATOM     13  H7  ALA     1      -1.139  -1.803   1.235  1.00  0.00           H
TER      14      ALA     1
END

There's a lot that happens before and after this, but not all prior steps are important and the subsequent steps make up much of the traceback.

Footnotes

1. I've only observed it for amino acids, but other standard-ish-looking things might be subject as well. ↩

2. I don't know if this is always true, but it's true enough in the files I've been looking at ↩

3. At least for now - I have not run this through Pablo ↩