Issue with `GAFFTemplateGenerator.generate_residue_template(residue_atoms=xxx)`

[ZhaoKe-BIT]

Apologies in advance if the phrasing sounds unusual—this issue is entirely translated by an LLM due to my limited English proficiency.

When using the GAFFTemplateGenerator.generate_residue_template function with the residue_atoms argument, I encountered an issue in line 745 of the code. The logic should iterate over residue_atoms when the argument is provided, and fall back to molecule.atoms when it is not. After implementing this change locally, my code worked as expected, so I believe this adjustment is reasonable.

[mattwthompson]

Do you have a minimal reproducing example?

[ZhaoKe-BIT]

Apologies for the delayed response—I’ve been occupied with other tasks recently. Below is my test code along with the PDB file I used:

from openff.toolkit import Molecule
from openmmforcefields.generators.template_generators import GAFFTemplateGenerator, SMIRNOFFTemplateGenerator
from openmm.app import ForceField, PDBFile
from io import StringIO

gys = r'c1cc(ccc1/C=C\2/C(=O)N(C(=N2)[C@H](CO)NC(=O)C)CC(=O)NC)O'
inner_atoms = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 20, 21, 22, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 39, 40, 45}
capped_gys = Molecule.from_smiles(gys, name='GYS')
gaff = GAFFTemplateGenerator()
ffxml_contents = gaff.generate_residue_template(capped_gys, [capped_gys.atom(i) for i in inner_atoms])
forcefield = ForceField(
    "amber/protein.ff14SB.xml",
    "amber/tip3p_standard.xml",
    "amber/tip3p_HFE_multivalent.xml",
    gaff.gaff_xml_filename
)
forcefield.loadFile(StringIO(ffxml_contents))
pdb = PDBFile('./2wur_watered.pdb')
system = forcefield.createSystem(pdb.topology)

Error Without the Fix:
When running the unmodified code, it throws the following error:

ValueError: No template found for residue 63 (GYS). The set of atoms is similar to [H][O][c]1[c]([H])[c]([H])[c](/[C]([H])=[C]2\[N]=[C]([C@@]([H])([N]([H])[C](=[O])[C]([H])([H])[H])[C]([H])([H])[O][H])[N]([C]([H])([H])[C](=[O])[N]([H])[C]([H])([H])[H])[C]2=[O])[c]([H])[c]1[H], but it is missing 12 atoms. For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#template

After applying the fix, the code runs successfully without errors.

2wur_watered.zip

Thank you for your check!

[epretti]

Right, it seems to only respect this option for bonds, but not for the atoms themselves. generate_residue_template() for SMIRNOFFTemplateGenerator and EspalomaTemplateGenerator appears to completely ignore the residue_atoms parameter: presumably this is a bug? When constructing a SmallMoleculeTemplateGenerator and providing molecules directly, this parameter never gets assigned, and this is the only path that the test cases cover, so it wasn't caught.

I will be out of office this upcoming week but can take a look at fixing this when I return.

[epretti]

It looks like you're trying to use openmmforcefields to generate parameters for a non-standard residue that is part of a protein chain, rather than a standalone small molecule ligand.

(@peastman / @mattwthompson / anyone else who may know) Is this something that the current template generators are supposed to support? From my limited understanding, the procedure for doing this, at least in the Amber case, might require assigning both Amber protein force field and GAFF parameters and is a bit more involved than what GAFFTemplateGenerator currently does. See an Amber tutorial on this, for example.

[mattwthompson]

I hesitated to respond to this originally because I'm not sure if there's something obvious that makes this work. I couldn't think of anything using this code that makes this use case simple.

I can, however say that in general makes little sense to use SMIRNOFFTemplateGenerator for anything larger than you'd run AM1-BCC on (a few hundred atoms). There are ways to use SMIRNOFF force fields with non-standard residues, but using different infrastructure - with more improvements in the pipeline.

[epretti]

Right, in that case, is there any legitimate use for the residue_atoms parameter of GAFFTemplateGenerator.generate_residue_template()? It seems like it is supposed to parameterize the whole molecule, remove atoms not selected, and change bonds between selected and non-selected atoms into external bonds. I was going to try to implement the same behavior in the other template generators, whose generate_residue_template() methods take this parameter but then silently ignore it contrary to what is documented. But the validity of the results that GAFFTemplateGenerator produces with this option enabled isn't clear to me. I don't know if the resulting charge assignments, for example, will make any sense, once atoms are deleted.

[epretti]

One path forward if we wanted to keep this behavior in GAFFTemplateGenerator and implement it in the other template generators would be to emulate what Amber's prepgen does in this case, and redistribute the charge of the excluded atoms among the remaining atoms. Some shared logic in SmallMoleculeTemplateGenerator could handle this and the problem described in #373 for all generators.

What concerns me is that "fixing" this feature may suggest to users that this is a recommended, standard automated way to parameterize a custom residue or other polymeric unit. To my understanding, it isn't: for GAFF, one ideally wants to assign certain parameters from the protein force field rather than GAFF; for SMIRNOFF, you pointed out the prototype for doing this; for Espaloma, I don't know.

[mattwthompson]

I can't speak for the intended GAFF behavior, but I think we'd be best off removing it from SMIRNOFFTemplateGenerator altogether. @j-wags is on leave until mid-June but I think he'd agree with me that we don't endorse this code path as the way to do that. I'd be surprised if it worked at all!

I looked at the history of the file and the feature is relatively untouched since its introduction many years ago. I suspect most people don't really know about it

https://github.com/openmm/openmmforcefields/blame/08e4c22dcf0f4bc0b318699d3db8ad91c8affba6/openmmforcefields/generators/template_generators.py#L618

[epretti]

I can't speak for the intended GAFF behavior, but I think we'd be best off removing it from SMIRNOFFTemplateGenerator altogether. @j-wags is on leave until mid-June but I think he'd agree with me that we don't endorse this code path as the way to do that. I'd be surprised if it worked at all!

For SMIRNOFF the feature won't work as this parameter seems to be silently ignored, so I agree that it can and should be removed; this change would break the API, but that should be OK since code using that parameter is broken as is. Espaloma is the same (it is just ignored). For GAFF, I will note to discuss with @tmarkland and @peastman to decide what to do.

@ZhaoKe-BIT: Making the change you did related to residue_atoms will make the code run, but you should be aware that the parameters and charges assigned might not be accurate, as this isn't a recommended method for parameterizing non-standard amino acid residues in proteins. The template generators in openmmforcefields are primarily designed for parameterizing small molecules, and you may want to use a different tool for what you are trying to do. For information on parameterizing non-standard residues with Amber, you may want to consult its manual or tutorials including but not limited to the tutorial already linked. If you do this, and eventually use LEaP to parameterize your system, OpenMM can read the resulting Amber PRMTOP file.