Skipping unknown residues in the model builder

[BastiaanvHoorn]

Hi,

I am using the system modeller for dealing with proteins for the EVB workflow that I am implementing. I have my own forcefield parameters for the active site of the protein and the reactant, and I want to use OpenMM to automatically assign the parameters in the rest of the protein based on the topology. However, when calling createsystem with the topology, I get an error because OpenMM does not know what to do with the unknown residues taking part in the reaction. The solution that I have now is that I manually delete any unknown residues, call the buildsystem, and then add them again, as follows:

forcefield = mmapp.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

pdb = mmapp.PDBFile('protein.pdb')
topology = pdb.topology
positions = pdb.positions

modeller = mmapp.Modeller(topology, positions)
residues = [residue for residue in topology.residues()]
chains = [chain for chain in topology.chains()]

modeller.delete([chains[-1]])

residueTemplates = {res: temp.name for res, temp in zip(residues, templates)}
system = forcefield.createSystem(
    modeller.topology,
)

However, this won't work when one of these unknown residues is covalently bound to the backbone. In this case, I cannot delete the unknown residue without causing other problems. Ideally in my case, there would be some way to force OpenMM to ignore residues that it does not know what to do with, but still create a system based on the rest of the protein where the interactions of the unknown residue are simply missing. I already add and/or overwrite the parameters of the atoms parttaking in the reaction anyway. I have considered using openff interchange for this to register a small molecule forcefield generator based on the structures of the unknown residues. However, if at all possible, I would like to not have to rely on another dependency for this.

edit: after more reading, I found the page on creating residue generators. However, I am not sure how to create a generator that would create a simple empty residue template. I tried the following code, but this gives me a different error because it doesn't recognise the atom types

...
templates, residues = forcefield.generateTemplatesForUnmatchedResidues(topology)
for temp in templates:
    forcefield.registerResidueTemplate(temp)

residueTemplates = {res: temp.name for res, temp in zip(residues, templates)}
system = forcefield.createSystem(
    topology,
    residueTemplates=residueTemplates,
)

Attached is the PDB file I'm using.
protein.pdb.zip

How could I go about tackling this?

[peastman]

I'm not completely sure I understand what you're doing. The goal is to end up with a single system that contains a complete, fully parameterized model, right? If you're building it with ForceField, the straightforward solution is to create a single ForceField that can parameterize the entire system, both the standard and nonstandard residues. That means writing an XML file with the parameters for the nonstandard residues, and loading it along with the Amber files. The user guide has a chapter on how to do that.

[BastiaanvHoorn]

The goal is to parameterize the protein-part of any protein system possibly containing non-standard residues with minimal user-intervention. I have my own code that automatically parameterizes the rest of the system, including non-standard residues.

[peastman]

Can you make your code generate the XML file automatically? A template generator would also be an option. Either way, it will be simpler if you can have the ForceField parameterize everything at once.

[BastiaanvHoorn]

Letting the forefield parameterize everything at once isn't possible because I'm using a bunch of non-standard combinations of forces to create alchemical states along a reaction pathway. I went with the following code in the end. The main difficulty was that the documentation for registeratomtype and registerresiduetemplate was rather limited and I ended up digging in the source code to figure out how to use these function.

pdb_file = mmapp.PDBFile(self.pdb)
topology = pdb_file.getTopology()
forcefield = mmapp.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
templates, residues = forcefield.generateTemplatesForUnmatchedResidues(topology)

for t, template in enumerate(templates):
    for i, atom in enumerate(template.atoms):
        forcefield.registerAtomType({
            'name':
            f'evb_{atom.name}_{t}',
            'class':
            'evb_placeholder',
            'mass':
            atom.element.mass.value_in_unit(mmunit.dalton),
            'element':
            atom.element
        })
        template.atoms[i].type = f'evb_{atom.name}_{t}'
    forcefield.registerResidueTemplate(template)

nbforce = forcefield.getGenerators()[2]
nbforce.registerAtom({
    'class': 'evb_placeholder',
    'charge': 0.0,
    'epsilon': 0.0,
    'sigma': 1.0,
})

system = forcefield.createSystem(
    topology,
    removeCMMotion=True,
)