Python_AddHydrogensToProtein

How can I add missing hydrogens to a protein?

BALL offers a fragment database BALL::FragmentDatabase to add missing information like hydrogens and bonds to a protein and to normalize the atom names.

In addition BALL offers a Minimizer class, that can be restricted by a selection, e.g. all hydrogens, and optimizes atom positions with respect to a ForceField.

import sys
from  BALL import *

# get the first system
system = getSystems()[0]

# print the number of atoms read from the file
print "read ", system.countAtoms(), " atoms."

# now we open a fragment database
print "reading fragment DB..." 
fragment_db = FragmentDB("")

# and normalize the atom names, i.e. we convert different
# naming standards to the PDB naming scheme - just in case!
print "normalizing names..." 
system.apply(fragment_db.normalize_names)

# now we add any missing hydrogens to the residues
# the data on the hydrogen positions stems from the
# fragment database. However the hydrogen positions 
# created in this way are only good estimates
print "creating missing atoms..." 
system.apply(fragment_db.add_hydrogens)
print "added ", fragment_db.add_hydrogens.getNumberOfInsertedAtoms(), " atoms" 

# now we create the bonds between the atoms (PDB files hardly
# ever contain a complete set of CONECT records)                                         
print "building bonds..."
system.apply(fragment_db.build_bonds)

# now we check whether the model we built is consistent
# The ResidueChecker checks for charges, bond lengths,
# and missing atoms
print "checking the built model..." 
checker = ResidueChecker(fragment_db)
system.apply(checker)

# now we create an AMBER force field 
print "setting up force field..." 
FF= AmberFF()

# we then select all hydrogens (element(H))
# using a specialized processor (Selector)
system.deselect()
FF.setup(system)
selector = Selector("element(H)")
system.apply(selector)

# just for curiosity: check how many atoms we are going
# to optimize
print "optimizing ", FF.getNumberOfMovableAtoms(), " out of ", system.countAtoms(), " atoms"    

# now we create a minimizer object that uses a conjugate 
# gradient algorithm to optimize the atom positions
minimizer = ConjugateGradientMinimizer()
initial_energy = FF.updateEnergy()
print "initial energy: ", initial_energy , " kJ/mol" 

# initialize the minimizer and perform (up to)
# 50 optimization steps
minimizer.setup(FF)
minimizer.setEnergyOutputFrequency(1)
minimizer.minimize(50)

# calculate the terminal energy and print it
terminal_energy = FF.getEnergy()
print "energy before/after minimization: ", initial_energy , "/", terminal_energy, " kJ/mol" 

# and trigger an update of the BALLView scene
getMainControl().update(system)

# done