I'm trying to model a zinc binding protein, where the zinc is
coordinated by two cysteines, a histidine and a glutamic acid. The bonds
are present in the CONECT records of the template.
Note that Modeller ignores the CONECT records.
I added upper bound restraints to restrain all distances between the
zinc and the four coordinating atoms, and used
"rsr.excluded_pairs.append" to exclude the bonds between the zinc and
the four coordinating atoms from the nonbonded list.
This sounds reasonable, although note that Modeller automatically builds
distance restraints between ligands/metal ions and nearby protein atoms
to maintain coordination (see
http://salilab.org/modeller/9v7/manual/node65.html), so these extra
restraints you have added may conflict. That might be the cause of your
"moving away" problem. You might want to try letting Modeller do its
default behavior here, as that might get better results. Alternatively,
if you want to define the restraints yourself, you could try removing
the default Modeller restraints by disabling the nonstd_restraints()
method, e.g. with something like:
class MyModel(automodel):
def nonstd_restraints(self, aln):
pass
a = MyModel(env, alnfile=..., knowns=..., sequence=...)
What is the best way to prevent this? Is there an easy way to define the
cysteines as unprotoneated?
Modeller's default topology includes no hydrogens anyway, but if you
want to change any residue, you would need to apply a patch by calling
model.patch() from within the special_patches() automodel method. See
http://salilab.org/modeller/9v7/manual/node24.html for an example. A
CHARMM patch can add or remove atoms from residues or generate new links
between pairs of residues. Refer to the CHARMM documentation for more
details.