openmm / pdbfixer

PDBFixer fixes problems in PDB files
Other
443 stars 112 forks source link

Invalide chirality when fixing neural network generated structure #274

Open ryankzhu opened 1 year ago

ryankzhu commented 1 year ago

Hi. I am trying to use pdbfixer to add sidechain and hydrogen atoms to a neural network generated peptide backbone (which contains C, Ca, N, and O atoms only). The quality of the artificial conformation is poor: Screenshot from 2023-06-06 13-25-11 So I ran an energy minimization for the system after modelling the missing atoms and solvating it in water, which fixed the problematic bond lengths and atom clashing. However, there are some chiral inversions in the energy minimized structure. For example, the Ca of the proline: image and the Cb of the threonine sidechain: image I wonder is it possible to enforce the hydrogen atoms to be added on the 'correct' side? How can I get the correct stereoisomers?

peastman commented 1 year ago

What force field are you minimizing with? Does it include energy terms to enforce correct chirality?

ryankzhu commented 1 year ago

I am using ff14sb for the peptide. I am not familiar with the energy terms in the forcefield. Which forcefield am I supposed to use in this scenario? Here is my codes:

from pdbfixer import PDBFixer
from openmm.app import PDBFile
from openmm.unit import *
from openmm.app.modeller import Modeller
from openmm import *

fixer = PDBFixer('tmp.pdb')
fixer.findMissingResidues()
fixer.findMissingAtoms()
fixer.addMissingAtoms()
PDBFile.writeFile(fixer.topology, fixer.positions, open('tmp_fixed.pdb', 'w'))

pdb = PDBFile('tmp_fixed.pdb')
mdl = Modeller(pdb.topology, pdb.positions)
forcefield = ForceField('amber14/protein.ff14SB.xml', 'amber14/tip3p.xml')
mdl.addHydrogens(forcefield, pH=7.0)
mdl.addSolvent(forcefield, padding=1*nanometers)
PDBFile.writeFile(mdl.topology, mdl.positions, open('solvated.pdb', 'w'))

system = forcefield.createSystem(mdl.topology, 
                                 nonbondedMethod=PME,
                                 nonbondedCutoff=1*nanometer, 
                                 constraints=HBonds)
simulation = Simulation(mdl.topology, system, integrator)
simulation.context.setPositions(mdl.positions)
simulation.minimizeEnergy()

positions = simulation.context.getState(getPositions=True).getPositions()
PDBFile.writeFile(simulation.topology, positions, open('em.pdb', 'w'))
peastman commented 1 year ago

I don't believe Amber includes any terms to enforce chirality. As long as you're starting from a good conformation that generally isn't a problem. But with a very noisy starting structure it could easily end up with the wrong chirality, and the force field doesn't do anything to prevent it. From a quick look at the Amber manual, I see it includes a program makeCHIR_RST to generate chirality restraints in this situation.

Does anyone else have experience with the best way of handling this situation?