Cannot add missing atoms

[kntkb]

Problem

I'm having trouble adding missing atoms with PDBFixer.addMissingAtoms. It looks like the missing atoms for residue 16 (LYS) is not added properly. I've attached the error and script for reproducibility.

Error

Traceback (most recent call last):
  File "/Users/takabak/work/playground/openmm_ecosystem/pdbfixer/example.py", line 40, in <module>
    model.addSolvent(ff, model=water_model, padding=box_padding, ionicStrength=salt_conc)
  File "/Users/takabak/mambaforge/envs/espaloma/lib/python3.9/site-packages/openmm/app/modeller.py", line 483, in addSolvent
    system = forcefield.createSystem(self.topology)
  File "/Users/takabak/mambaforge/envs/espaloma/lib/python3.9/site-packages/openmm/app/forcefield.py", line 1212, in createSystem
    templateForResidue = self._matchAllResiduesToTemplates(data, topology, residueTemplates, ignoreExternalBonds)
  File "/Users/takabak/mambaforge/envs/espaloma/lib/python3.9/site-packages/openmm/app/forcefield.py", line 1427, in _matchAllResiduesToTemplates
    raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
ValueError: No template found for residue 16 (LYS).  The set of atoms is similar to ILE, but it is missing 1 hydrogen atoms.

Script

#!/usr/bin/env python
# coding: utf-8
import os, sys, math
import glob as glob
import numpy as np
import re
import warnings
import urllib
from openmm import *
from openmm.app import *
from openmm.unit import *
from pdbfixer import PDBFixer

# download pdb
urllib.request.urlretrieve('http://files.rcsb.org/download/4w52.pdb', '4w52.pdb')

# pdbfixer
pdbfile = "4w52.pdb"
fixer = PDBFixer(filename=pdbfile)
fixer.findMissingResidues()
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()
fixer.removeHeterogens(True)
fixer.findMissingAtoms()
fixer.addMissingHydrogens(7.0)  # default: 7
PDBFile.writeFile(fixer.topology, fixer.positions, open('t4-lysozyme-L99A.pdb', 'w'))

# define system
box_padding = 12.0 * angstrom
salt_conc = 0.15 * molar
nb_cutoff = 10 * angstrom
H_mass = 3.5 * amu #Might need to be tuned to 3.5 amu 
water_model='tip3p'

# solvate 
ff = ForceField('amber/RNA.OL3.xml', 'amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml')
model = Modeller(fixer.topology, fixer.positions)
model.addSolvent(ff, model=water_model, padding=box_padding, ionicStrength=salt_conc)

# save solvated system
topology = model.getTopology()
positions = model.getPositions()
PDBFile.writeFile(topology, positions, file=open('solvated.pdb', 'w'))

# create system
system = ff.createSystem(model.topology, nonbondedMethod=PME, nonbondedCutoff=nb_cutoff, constraints=HBonds, rigidWater=True, hydrogenMass=H_mass)

[peastman]

You're missing a function call. Immediately after the call to findMissingAtoms(), you need to add the line

fixer.addMissingAtoms()

Otherwise the missing atoms it has identified don't get added.

[kntkb]

Oh this is a careless mistake. How can I have missed that! Thank you.