Adding two terminal caps to a single residue

[ljmartin]

Hi PDBFixer, I'd like to model single-residue fragments in OpenMM, which have been extracted from a protein chain. Since they lack backbone bonds to other residues, they won't parameterize automatically in OpenMM. I thought to add terminal patches with pdbfixer, but only one gets generated:

import prody
name = prody.fetchPDB('3GBZ')
prody_pdb = prody.parsePDB(name)
prody_sel = prody_pdb.select(f'resnum {prody_pdb.getResnums()[0]}')
print(prody_sel.getResnames()[0])
prody.writePDB('ser.pdb', prody_sel)

from pdbfixer import PDBFixer
from openmm.app import PDBFile
fixer = PDBFixer(filename='ser.pdb')
fixer.findMissingResidues()
fixer.findMissingAtoms()

print(fixer.missingTerminals)

output:

{<Residue 0 (SER) of chain 0>: ['OXT']}

Indeed the result isn't recognised:

fixer.addMissingAtoms()
PDBFile.writeFile(fixer.topology, fixer.positions, open('output.pdb', 'w'))

from openmm.app import ForceField
pdbfile = PDBFile('output.pdb')
forcefield = ForceField('amber/protein.ff14SB.xml',)
system = forcefield.createSystem(pdbfile.topology)

returns No template found for residue 1 (SER)

If this is even sensible - Is there a way to add both terminal caps to a single residue? Thanks

[peastman]

In principle yes, you should be able to do that. The problem is that Amber doesn't support it. For each amino acid it includes three versions: C terminal, N terminal, and middle of chain. It doesn't define parameters for ones that are terminated at both ends.

I believe CHARMM36 does support isolated amino acids, so you might be able to use it.

missingTerminals only includes heavy atoms, not hydrogens. That's why you don't see the extra N-terminal hydrogen there.

[ljmartin]

Thank you very much, that works fine.

fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
PDBFile.writeFile(fixer.topology, fixer.positions, open('output.pdb', 'w'))

from openmm.app import ForceField
pdbfile = PDBFile('output.pdb')
forcefield = ForceField('charmm36.xml',)
system = forcefield.createSystem(pdbfile.topology)