# vim: set expandtab ts=4 sw=4: def simpleNormalizeName(name): return name.strip() class Unimplemented(ValueError): pass ExtraSteps = 2 # number of extra steps to add to MMTK minimization # Execute two extra steps because # step 0 only calculates energy and not gradient # step 1 calculates gradient but positions are same as 0 # step 2 is the first step where atoms move def mmtkIdent(obj): import chimera if isinstance(obj, chimera.Atom): ident = obj.name.replace('.', '_') if ident[0].isdigit(): ident = '_' + ident if len(obj.residue.atomsMap[obj.name]) > 1: # non-unique atom names ident += '_' + str(obj.residue.atomsMap[obj.name].index(obj)) return ident ident = obj.oslIdent().replace('#', '_0_').replace(':', '_1_').replace( '.', '__') if isinstance(obj, chimera.Residue): # so that residue type can be deduced when read back in, # as well as the ordering of the residues... ident = obj.type + '_' + str(obj.molecule.residues.index(obj)) \ + "_" + ident return ident class MMTKinter: def __init__(self, mols, exclres=set(), nogui=False, addhyd=True, memorize=False, prep=True, ljOptions=None, esOptions=None, callback=None): # MMTK lengths are in nm while Chimera ones are in Angstroms # so we need a scale factor when converting if not mols: raise ValueError("No molecules specified") self.tempDir = None self.molId = 0 self.ljOptions = ljOptions self.esOptions = esOptions self.callback = callback self.mols = mols self.exclres = exclres self._getParameters(mols, nogui, addhyd, memorize, prep) def _finishInit(self): timestamp("_finishInit") self.molecules = [] self.atomMap = {} try: for m in self.mols: self.molecules.extend(self._makeMmtkMolecules(m)) self._makeUniverse() if self.callback: self.callback(self) del self.callback finally: self._removeTempDir() def _makeUniverse(self): import os.path parmDir = os.path.dirname(__file__) timestamp("_makeUniverse") from MMTK import InfiniteUniverse from MMTK.ForceFields.Amber import AmberData from MMTK.ForceFields.Amber.AmberForceField import readAmber99 from MMTK.ForceFields.MMForceField import MMForceField # # GAFF uses lower case atom types to distinguish # from Amber atom types. MMTK, however, normalizes # all atom types to upper case. So we hack MMTK # and temporarily replace _normalizeName function # with ours while reading our parameter files. # (We have to reread the parameter file each time # because we potentially have different frcmod files # for the different universes.) # saveNormalizeName = AmberData._normalizeName AmberData._normalizeName = simpleNormalizeName modFiles = set([ str(m.frcmod) for m in self.molecules if m.frcmod is not None]) from AmberInfo import amberHome paramDir = os.path.join(amberHome, "dat", "leap", "parm") parameters = readAmber99(os.path.join(paramDir, "gaff.dat"), modFiles) self._mergeAmber99(parameters) bondedScaleFactor = 1.0 ff = MMForceField("Amber99/GAFF", parameters, self.ljOptions, self.esOptions, bondedScaleFactor) AmberData._normalizeName = saveNormalizeName timestamp("Made forcefield") self.universe = InfiniteUniverse(ff) timestamp("Made universe") for mm in self.molecules: self.universe.addObject(mm) timestamp("Added model %s" % mm.name) timestamp("end _makeUniverse") def _getParameters(self, mols, nogui, addhyd, memorize, prep): if not prep: self.addedAtoms = [] self._finishInit() return timestamp("_getParameters") import DockPrep import chimera self.originalAtoms = set([]) for m in mols: self.originalAtoms.update(m.atoms) from AddCharge import defaultChargeModel kw = { "doneCB": self._finishDockPrep, "gaffType": True, "chargeModel": defaultChargeModel } if nogui or chimera.nogui: if not addhyd: kw["addHFunc"] = None DockPrep.prep(mols, nogui=nogui, **kw) else: d = DockPrep.memoryPrep("Minimize", memorize, mols, nogui=nogui, **kw) if d: d.addHydVar.set(addhyd) d.writeMol2Var.set(False) def _finishDockPrep(self): timestamp("end _getParameters") self.addedAtoms = sum([[a for a in m.atoms if not a in self.originalAtoms] for m in self.mols], []) # If atom was added to a selected atom, then we add it # into the current selection as well. from chimera import selection attached = attachedAtoms(self.addedAtoms, selection.currentAtoms()) if len(attached) > 0: selection.addCurrent(attached) del self.originalAtoms self._finishInit() def _mergeAmber99(self, parm): "Merge parameters appropriate to charge model into our parameters" from MMTK.ForceFields.Amber.AmberForceField import readAmber99 chargeModels = set([getattr(m, 'chargeModel', None) for m in self.mols]) if None in chargeModels: if len(chargeModels) == 1: from AddCharge import defaultChargeModel chargeModels = set([defaultChargeModel]) else: chargeModels.remove(None) if len(chargeModels) > 1: from chimera import LimitationError raise LimitationError("Molecules have differing charge models (%s);" " don't know what parameter set to use." % ", ".join(list(chargeModels))) mainFile, modFiles = chargeModelToFiles(chargeModels.pop()) from chimera import replyobj import os.path if modFiles: replyobj.info("Using main parameter file %s modified by %s\n" % (os.path.basename(mainFile), ", ".join([os.path.basename(mf) for mf in modFiles]))) else: replyobj.info("Using main parameter file %s with no modifications\n" % os.path.basename(mainFile)) parm99 = readAmber99(main_file=mainFile, mod_files=[open(mf, "r") for mf in modFiles]) parm.atom_types.update(parm99.atom_types) parm.bonds.update(parm99.bonds) parm.bond_angles.update(parm99.bond_angles) parm.dihedrals.update(parm99.dihedrals) parm.dihedrals_2.update(parm99.dihedrals_2) parm.impropers.update(parm99.impropers) parm.impropers_1.update(parm99.impropers_1) parm.impropers_2.update(parm99.impropers_2) parm.hbonds.update(parm99.hbonds) parm.lj_equivalent.update(parm99.lj_equivalent) for name, ljpar_set99 in parm99.ljpar_sets.iteritems(): try: ljpar_set = parm.ljpar_sets[name] except KeyError: parm.ljpar_sets[name] = ljpar_set99 else: if ljpar_set.type != ljpar_set99.type: print "incompatible ljpar_set" print " GAFF type:", ljpar_set.type print " AMBER99 type:", ljpar_set99.type ljpar_set.entries.update(ljpar_set99.entries) def _makeMmtkMolecules(self, m): timestamp("_makeMmtkMolecules %s" % m.name) mols = MMTKChimeraModel(m, self.molId, self.exclres, self).retrieveMolecules() self.molId += 1 timestamp("end _makeMmtkMolecules") return mols def setFixed(self, atoms): for ma in self.universe.atomList(): ma.fixed = False for a in atoms: if a in self.atomMap: ma = self.atomMap[a] ma.fixed = True # Fix added atoms that are connected to fixed atoms. for a in attachedAtoms(self.addedAtoms, atoms): if a in self.atomMap: ma = self.atomMap[a] ma.fixed = True def loadMMTKCoordinates(self): "Load MMTK coordinates from Chimera models" import chimera from Scientific.Geometry import Vector from MMTK import Units s = Units.Ang for ma in self.universe.atomList(): try: ca = ma.getAtomProperty(ma, "chimera_atom") except AttributeError: pass else: c = ca.coord() p = Vector(c[0] * s, c[1] * s, c[2] * s) ma.setPosition(p) def saveMMTKCoordinates(self): "Save MMTK coordinates into Chimera models" import chimera from chimera import Coord from MMTK import Units s = Units.Ang sum = 0.0 count = 0 for ma in self.universe.atomList(): if ma.fixed: continue ca = ma.getAtomProperty(ma, "chimera_atom") p = ma.position() c = Coord(p[0] / s, p[1] / s, p[2] / s) dsq = (c - ca.coord()).sqlength() ca.setCoord(c) #print "%.6f" % dsq sum += dsq count += 1 import math if count > 0: print "Updated", count, "atoms. RMSD: %.6f" \ % math.sqrt(sum / count) else: print "No atoms updated." def minimize(self, nsteps, stepsize=0.02, cgsteps=0, cgstepsize=0.02, interval=None, action=None, **kw): timestamp("_minimize") if not interval: actions = [] else: import sys from MMTK.Trajectory import LogOutput report = interval + ExtraSteps actions = [ LogOutput(sys.stdout, ["energy"], report, None, report) ] from MMTK import Units if action is None or not interval: interval = None from MMTK.ForceFields.Amber import AmberData saveNormalizeName = AmberData._normalizeName AmberData._normalizeName = simpleNormalizeName from chimera import replyobj try: msg = "Initial energy: %f kJ/mol" % self.universe.energy() except KeyError, e: msg = ("Chimera/MMTK cannot minimize structure, " "probably because there is no parameter " "for '%s'" % e) replyobj.error(msg) return else: replyobj.info(msg) if nsteps: replyobj.status("starting %d steps of steepest descent" % nsteps) kw["step_size"] = stepsize * Units.Ang from MMTK.Minimization import SteepestDescentMinimizer minimizer = SteepestDescentMinimizer(self.universe, actions=actions, **kw) self._doMinimize(minimizer, "steepest descent", nsteps, interval, action) if cgsteps: replyobj.status("starting %d steps of conjugate gradient" % cgsteps) kw["step_size"] = cgstepsize * Units.Ang from MMTK.Minimization import ConjugateGradientMinimizer minimizer = ConjugateGradientMinimizer(self.universe, actions=actions, **kw) self._doMinimize(minimizer, "conjugate gradient", cgsteps, interval, action) AmberData._normalizeName = saveNormalizeName timestamp("end _minimize") def _doMinimize(self, minimizer, name, steps, interval, action): from chimera import replyobj remaining = steps while remaining > 0: timestamp(" minimize interval") if interval is None: realSteps = remaining else: realSteps = min(remaining, interval) minimizer(steps=realSteps + ExtraSteps) remaining -= realSteps if action is not None: action(self) timestamp(" finished %d steps" % realSteps) msg = "Finished %d of %d %s minimization steps" % ( steps - remaining, steps, name) replyobj.status(msg) replyobj.info(msg) replyobj.info("\n") def getTempDir(self): if self.tempDir: return self.tempDir from tempfile import mkdtemp self.tempDir = mkdtemp() #self.tempDir = "." return self.tempDir def _removeTempDir(self): if not self.tempDir: return if True: import os, os.path for filename in os.listdir(self.tempDir): os.remove(os.path.join(self.tempDir, filename)) os.rmdir(self.tempDir) else: print "Did not clean up temp dir", self.tempDir self.tempDir = None def attachedAtoms(atoms, baseAtoms): bset = set(baseAtoms) attached = [] for a in atoms: for b in a.bonds: oa = b.otherAtom(a) if oa in bset: attached.append(a) break return attached from MMTK.MoleculeFactory import MoleculeFactory, AtomTemplate, GroupTemplate class ChimeraGroupTemplate(GroupTemplate): """GroupTemplate that handles containing actual Groups/Atoms""" def __init__(self, *args, **kw): GroupTemplate.__init__(self, *args, **kw) def atomNameToPath(self, atom_name): atom_name = atom_name.split('.') object = self from MMTK.ChemicalObjects import Group, Atom try: for path_element in atom_name: if isinstance(object, Group): for subobj in object.atomList(): if subobj.name == path_element or subobj.chimera_atom.name == path_element: object = subobj if isinstance(object, Atom): return object break else: break else: object = object.children[object.names[path_element]] if not isinstance(object, (Atom, AtomTemplate)): raise ValueError("no atom " + '.'.join(atom_name)) except KeyError: raise ValueError("no atom " + '.'.join(atom_name)) return atom_name class ChimeraMoleculeFactory(MoleculeFactory): """MoleculeFactory that supports containing actual Groups/Atoms""" def createGroup(self, name): """ Create a new (initially empty) group object. """ if self.groups.has_key(name): raise ValueError("redefinition of group " + name) self.groups[name] = ChimeraGroupTemplate(name) def makeChemicalObjects(self, template, top_level): from MMTK import Bonds, ChemicalObjects self.groups[template.name].locked = True if top_level: if template.attributes.has_key('sequence'): object = ChemicalObjects.ChainMolecule(None) else: object = ChemicalObjects.Molecule(None) else: object = ChemicalObjects.Group(None) object.atoms = [] object.bonds = Bonds.BondList([]) object.groups = [] object.type = self.groups[template.name] object.parent = None child_objects = [] for child in template.children: if isinstance(child, GroupTemplate): group = self.makeChemicalObjects(child, False) object.groups.append(group) object.atoms.extend(group.atoms) object.bonds.extend(group.bonds) group.parent = object child_objects.append(group) elif not isinstance(child, AtomTemplate): # actual Molecule or Group, returned by # _findStandardResidue object.groups.append(child) object.atoms.extend(child.atoms) object.bonds.extend(child.bonds) child.parent = object child_objects.append(child) else: try: atom = ChemicalObjects.Atom(child.element) except IOError: from chimera import LimitationError raise LimitationError("Atom type \"%s\" " "is not supported by MMTK" % child.element) object.atoms.append(atom) atom.parent = object child_objects.append(atom) a = self.atomMap[child] del self.atomMap[child] self.atomMap[a] = atom for name, index in template.names.items(): setattr(object, name, child_objects[index]) child_objects[index].name = name for name, value in template.attributes.items(): path = name.split('.') setattr(self.namePath(object, path[:-1]), path[-1], value) for atom1, atom2 in template.bonds: if not isinstance(atom1, ChemicalObjects.Atom): atom1 = self.namePath(object, atom1) if not isinstance(atom2, ChemicalObjects.Atom): atom2 = self.namePath(object, atom2) object.bonds.append(Bonds.Bond((atom1, atom2))) for name, vector in template.positions.items(): path = name.split('.') self.namePath(object, path).setPosition(vector) return object class MMTKChimeraModel(object): def __init__(self, m, ident, exclres, owner): from MMTK import Bonds self.chimeraMolecule = m molResGroups = connectedResidues(m, exclres) self.atomMap = owner.atomMap self.mf = ChimeraMoleculeFactory() self.mf.atomMap = self.atomMap rtype2frcmod = {} self.frcmods = frcmods = [] res2grp = {} for i, mrg in enumerate(molResGroups): self.mf.createGroup(str(i)) molGroup = self.mf.groups[str(i)] needParmchk = {} for r in mrg: res2grp[r] = molGroup v = self._findStandardResidue(r) if v is None: self._makeNonStandardResidue(str(i), r) needParmchk[r.type] = r else: self._makeStandardResidue(molGroup, r, *v) if needParmchk: frcmod = None # are all needed residue types already in some # frcmod file we've computed? for rtype in needParmchk.keys(): if rtype in rtype2frcmod: if frcmod is None: frcmod = rtype2frcmod[rtype] elif frcmod != rtype2frcmod[rtype]: frcmod = None break else: break if frcmod is None: frcmod = self._runParmchk("%s-%d" % (ident, i), owner.getTempDir(), needParmchk) for rtype in needParmchk.keys(): rtype2frcmod[rtype] = frcmod frcmods.append(frcmod) else: frcmods.append(None) from chimera import Bond for b in m.bonds: if b.display == Bond.Never: continue a0, a1 = b.atoms if a0.residue is a1.residue: continue if a0.residue in exclres or a1.residue in exclres: continue res2grp[a0.residue].addBond( mmtkIdent(a0.residue) + "." + mmtkIdent(a0), mmtkIdent(a1.residue) + "." + mmtkIdent(a1)) def retrieveMolecules(self): mols = [] for i, frcmod in enumerate(self.frcmods): mol = self.mf.retrieveMolecule(str(i)) mol.type.frcmod = frcmod mols.append(mol) return mols def _findStandardResidue(self, r): try: amberName = r.amberName except AttributeError: pass else: from AddCharge import unitedAtomChargeModels if r.molecule.chargeModel in unitedAtomChargeModels: resMaps = [ UnitedAtomResidueNameMap, ResidueNameMap ] else: resMaps = [ ResidueNameMap ] for rm in resMaps: try: blueprint = rm[amberName] except KeyError: pass else: from MMTK.ChemicalObjects import Group return Group(blueprint), amberName try: blueprint = MoleculeNameMap[r.type] except KeyError: return None else: from MMTK.ChemicalObjects import Molecule return Molecule(blueprint), r.type def _makeStandardResidue(self, molGroup, r, mg, blueprint): chimera2mmtk = self._mapStandardAtoms(r, mg) if (len(chimera2mmtk) == len(mg.atomList()) and len(chimera2mmtk) == len(r.atoms)): self._addStandardResidue(molGroup, mg, chimera2mmtk) return # Some atoms were not used. # If residue is a histidine, try other protonation forms. resType = blueprint[-3:] hisList = [ "HIE", "HID", "HIP", "HIS" ] if resType in hisList: prefix = blueprint[:-3] from MMTK.ChemicalObjects import Group from chimera import LimitationError for hisType in hisList: if hisType == resType: continue newType = prefix + hisType mg = Group(ResidueNameMap[newType]) try: chimera2mmtk = self._mapStandardAtoms(r, mg) except LimitationError: # Must have hit a Chimera atom with # no corresponding MMTK atom when # using the wrong blueprint continue if (len(chimera2mmtk) == len(mg.atomList()) and len(chimera2mmtk) == len(r.atoms)): break from chimera import replyobj replyobj.warning("histidine %s reclassified " "from AMBER type %s to %s\n" % (r.oslIdent(), blueprint, newType)) self._addStandardResidue(molGroup, mg, chimera2mmtk) return # There is an irreconcilable atom mismatch. allMMTKAtoms = set(mg.atomList()) allChimeraAtoms = set(r.atoms) usedMMTKAtoms = set(chimera2mmtk.itervalues()) usedChimeraAtoms = set(chimera2mmtk.iterkeys()) extra = [ a.name for a in allChimeraAtoms - usedChimeraAtoms ] if extra: extraAtoms = " has extra " + makeAtomList(extra) else: extraAtoms = None missing = [ ma.name for ma in allMMTKAtoms - usedMMTKAtoms ] if missing: missingAtoms = " is missing " + makeAtomList(missing) else: missingAtoms = None if extra and missing: msg = extraAtoms + " and" + missingAtoms elif extra: msg = extraAtoms else: msg = missingAtoms from chimera import LimitationError raise LimitationError("Residue %s (%s/%s) %s" % (r.oslIdent(), r.type, mg.name, msg)) def _mapStandardAtoms(self, r, mg): pdbmap = {} altmap = {} self._getMaps(mg, pdbmap, altmap) try: subgroups = mg.groups except AttributeError: pass else: for subg in mg.groups: self._getMaps(subg, pdbmap, altmap) used = {} chimera2mmtk = {} for a in r.atoms: aname = self.getMMTKname(a.name, r.type, pdbmap, altmap) ma = mg.getAtom(aname) if ma in used: from chimera import LimitationError raise LimitationError("Residue %s (%s/%s) should have either atom %s " "or %s, but not both" % (r.oslIdent(), r.type, mg.name, used[ma], a.name)) chimera2mmtk[a] = ma used[ma] = a.name return chimera2mmtk def getMMTKname(self, aname, rtype, pdbmap, altmap): # If the name does not start with an alphabetic # character, try rotating it until it does. This works # for hydrogens in amino acids. letters = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ" n = aname while n[0] not in letters: n = n[1:] + n[0] n = altmap.get(n, n) if pdbmap.has_key(n): return pdbmap[n] # Nope. If the name contains a ', try replacing it with # a *. This works for hydrogens in nucleotides. n = aname.replace("'", '*') n = altmap.get(n, n) if pdbmap.has_key(n): return pdbmap[n] #print pdbmap.keys() #print altmap.keys() from chimera import LimitationError raise LimitationError("No MMTK name for atom \"%s\" " "in standard residue \"%s\"" % (aname, rtype)) def _getMaps(self, obj, pdbmap, altmap): try: pm = obj.pdbmap except AttributeError: pass else: for item in pm: for name, ref in item[1].iteritems(): atom = obj.getAtom(ref) pdbmap[name] = atom try: am = obj.pdb_alternative except AttributeError: pass else: altmap.update(am) def _addStandardResidue(self, molGroup, mg, chimera2mmtk): for a, ma in chimera2mmtk.iteritems(): properties = { "amber_charge": a.charge, "amber_atom_type": a.gaffType, "chimera_atom": a, } ma.addProperties(properties) self.atomMap[a] = ma molGroup.addSubgroup(mmtkIdent(a.residue), mg) def _makeNonStandardResidue(self, molGroupName, r): #print "makeNonstandardResidue", r.oslIdent(), r.type import chimera atoms = [] c2m = {} residueBonds = set([]) rIdent = mmtkIdent(r) self.mf.createGroup(rIdent) rg = self.mf.groups[rIdent] self.mf.addSubgroup(molGroupName, rIdent, rIdent) for a in r.atoms: aIdent = mmtkIdent(a) rg.addAtom(aIdent, a.element.name) ma = rg.children[-1] try: charge = a.charge except AttributeError: from chimera import LimitationError raise LimitationError("Unable to find " "partial charge for %s" % a.oslIdent()) try: gaffType = a.gaffType except AttributeError: from chimera import LimitationError raise LimitationError("Element %s (atom %s)" " is not currently supported [no GAFF type]" % (a.element.name, a.oslIdent())) for attrName, val in [ ("amber_charge", charge), ("amber_atom_type", gaffType), ("chimera_atom", a)]: rg.setAttribute(aIdent + "." + attrName, val) c2m[a] = ma for b in a.bonds: if b.display == chimera.Bond.Never: continue other = b.otherAtom(a) if other.residue == r: residueBonds.add(b) # yes, map is reversed for templates, # corrected when actual objects are made self.atomMap[ma] = a for b in residueBonds: a0, a1 = b.atoms rg.addBond(mmtkIdent(a0), mmtkIdent(a1)) if len(r.atoms) == 1: atom = r.atoms[0] if len(atom.neighbors) == 0: gaffType = atom.gaffType if gaffType.upper() == gaffType or not gaffType.isalnum(): # don't parmchk known ions return def _runParmchk(self, uniquifier, tempDir, needParmchk): import os, os.path, sys from subprocess import Popen, STDOUT, PIPE from chimera import replyobj parmDir = os.path.dirname(__file__) if not parmDir: parmDir = os.getcwd() parmchkIn = os.path.join(tempDir, "parmchk.in.%s" % uniquifier) self._writeParmchk(parmchkIn, needParmchk) frcmod = os.path.join(tempDir, "frcmod.%s" % uniquifier) chimeraRoot = os.environ["CHIMERA"] from AmberInfo import amberBin, amberHome command = [ os.path.join(amberBin, "parmchk"), "-i", parmchkIn, "-f", "mol2", "-o", frcmod, "-p", os.path.join(amberHome, "dat", "leap", "parm", "gaff.dat") ] replyobj.status("Running PARMCHK for %s" % self.chimeraMolecule.name, log=True) replyobj.info("command: %s\n" % " ".join(command)) p = Popen(command, stdin=PIPE, stdout=PIPE, stderr=STDOUT, cwd=tempDir, shell=False, env={"AMBERHOME": amberHome}, bufsize=1).stdout while True: line = p.readline() if not line: break replyobj.info("(parmchk) %s\n" % line.rstrip()) if not os.path.exists(frcmod): from chimera import LimitationError raise LimitationError("Unable to compute partial " "charges: PARMCHK failed.\n" "Check reply log for details.\n") return None replyobj.status("Finished PARMCHK for %s" % self.chimeraMolecule.name, log=True) return frcmod def _writeParmchk(self, filename, needParmchk): # generate Mol2 input file for parmchk import WriteMol2 from chimera.selection import ItemizedSelection from chimera import Bond rSet = set(needParmchk.values()) npSet = set(needParmchk.values()) for b in self.chimeraMolecule.bonds: if b.display == Bond.Never: continue a0, a1 = b.atoms r0 = a0.residue r1 = a1.residue if (r0 in npSet and r1 not in npSet): rSet.add(r1) elif (r0 not in npSet and r1 in npSet): rSet.add(r0) sel = ItemizedSelection() sel.add(rSet) WriteMol2.writeMol2(sel, filename, gaffType=True, temporary=True) class MMTKChimeraGroupType: # This class is needed to supply the necessary attributes # so that MMTK description() method will complete successfully. def __init__(self): self.name = "ChimeraGroup" self.groups = [] ChimeraGroupType = MMTKChimeraGroupType() from MMTK import Biopolymers class MMTKChimeraGroup(Biopolymers.Residue): def __init__(self, r, atomMap): import chimera from MMTK.ChemicalObjects import Atom from MMTK.Bonds import Bond self.type = ChimeraGroupType # see class above self.parent = None self.name = r.oslIdent() atoms = [] c2m = {} residueBonds = set([]) for a in r.atoms: try: ma = Atom(a.element.name) except IOError: from chimera import LimitationError raise LimitationError("Atom type \"%s\" " "is not supported by MMTK" % a.element.name) ma.parent = self try: charge = a.charge except AttributeError: from chimera import LimitationError raise LimitationError("Unable to find " "partial charge for %s" % a.oslIdent()) try: gaffType = a.gaffType except AttributeError: from chimera import LimitationError raise LimitationError("Element %s (atom " "%s) is not currently supported" % (a.element.name, a.oslIdent())) properties = { "amber_charge": charge, "amber_atom_type": gaffType, "chimera_atom": a, } ma.addProperties(properties) atoms.append(ma) c2m[a] = ma for b in a.bonds: if b.display == chimera.Bond.Never: continue other = b.otherAtom(a) if other.residue == r: residueBonds.add(b) atomMap[a] = ma bonds = [] for b in residueBonds: a0, a1 = b.atoms ma0 = c2m[a0] ma1 = c2m[a1] mb = Bond((ma0, ma1)) mb.parent = self bonds.append(mb) self.atoms = atoms self.bonds = bonds def _dumpChain(obj): while obj is not None: print "object", obj _dumpDataAttributes(obj) obj = getattr(obj, "parent", None) def _dumpDataAttributes(mg): for attr in dir(mg): if attr.startswith("__"): continue val = getattr(mg, attr) if callable(val): continue print " ", attr, val def connectedResidues(mol, excludedResidues): connectedSets = [] seenResidues = set() for res in mol.residues: if res in seenResidues or res in excludedResidues: continue roots = set([res.molecule.rootForAtom(a, True) for a in res.atoms]) connected = set([a.residue for rt in roots for a in mol.traverseAtoms(rt)]) connectedSets.append(connected) seenResidues.update(connected) return connectedSets def makeAtomList(names): # Assume names is not empty if len(names) == 1: return "atom %s" % names[0] else: return "atoms %s and %s" % (", ".join(names[:-1]), names[-1]) def chargeModelToFiles(chargeModel): expectedPrefix = "AMBER ff" from chimera import LimitationError if not chargeModel.startswith(expectedPrefix): raise LimitationError("Don't know how to find parameter files for charge model '%s'" % chargeModel) ffName = chargeModel[len(expectedPrefix):] if not ffName[:2].isdigit(): raise LimitationError("Charge model name not of the expected form (starting with" " %s[year])" % expectedPrefix) from AmberInfo import amberHome import os.path paramDir = os.path.join(amberHome, "dat", "leap", "parm") year = int(ffName[:2]) mainFile = os.path.join(paramDir, "parm" + ffName + ".dat") tries = 0 while not os.path.exists(mainFile): mainFile = os.path.join(paramDir, "parm%02d.dat" % year) year -= 1 if year < 0: year += 100 tries += 1 if tries > 100: raise AssertionError("No parm files in %s?!?" % paramDir) modFiles = [os.path.join(paramDir, "heme-iron.frcmod")] baseFF = os.path.join(paramDir, "frcmod.ff" + ffName[:2]) if os.path.exists(baseFF): baseSuffixes = ["ff" + ffName[:2]] else: baseSuffixes = [] for suffix in baseSuffixes + ["ff"+ffName, "parm"+ffName[2:], "ionsjc_tip3p"]: modFile = os.path.join(paramDir, "frcmod." + suffix) if os.path.exists(modFile): modFiles.append(modFile) return mainFile, modFiles def timestamp(s): pass #import time #print "%s: %s" % (time.ctime(time.time()), s) ResidueNameMap = { "ACE": "ace_beginning_nt", "NME": "methyl", "ALA": "alanine", "ARG": "arginine", "ASP": "aspartic_acid", "ASN": "asparagine", "CYS": "cysteine", "CYX": "cystine_ss", "GLU": "glutamic_acid", "GLN": "glutamine", "GLY": "glycine", "HIS": "histidine", "HID": "histidine_deltah", "HIE": "histidine_epsilonh", "HIP": "histidine_plus", "ILE": "isoleucine", "LEU": "leucine", "LYS": "lysine", "MET": "methionine", "PHE": "phenylalanine", "PRO": "proline", "SER": "serine", "NA": "sodium", "THR": "threonine", "TRP": "tryptophan", "TYR": "tyrosine", "VAL": "valine", "CALA": "alanine_ct", "CARG": "arginine_ct", "CASP": "aspartic_acid_ct", "CASN": "asparagine_ct", "CCYS": "cysteine_ct", "CCYX": "cystine_ss_ct", "CGLU": "glutamic_acid_ct", "CGLN": "glutamine_ct", "CGLY": "glycine_ct", "CHIS": "histidine_ct", "CHID": "histidine_deltah_ct", "CHIE": "histidine_epsilonh_ct", "CHIP": "histidine_plus_ct", "CILE": "isoleucine_ct", "CLEU": "leucine_ct", "CLYS": "lysine_ct", "CMET": "methionine_ct", "CPHE": "phenylalanine_ct", "CPRO": "proline_ct", "CSER": "serine_ct", "CNA": "sodium_ct", "CTHR": "threonine_ct", "CTRP": "tryptophan_ct", "CTYR": "tyrosine_ct", "CVAL": "valine_ct", "NALA": "alanine_nt", "NARG": "arginine_nt", "NASP": "aspartic_acid_nt", "NASN": "asparagine_nt", "NCYS": "cysteine_nt", "NCYX": "cystine_ss_nt", "NGLU": "glutamic_acid_nt", "NGLN": "glutamine_nt", "NGLY": "glycine_nt", "NHIS": "histidine_nt", "NHID": "histidine_deltah_nt", "NHIE": "histidine_epsilonh_nt", "NHIP": "histidine_plus_nt", "NILE": "isoleucine_nt", "NLEU": "leucine_nt", "NLYS": "lysine_nt", "NMET": "methionine_nt", "NPHE": "phenylalanine_nt", "NPRO": "proline_nt", "NSER": "serine_nt", "NNA": "sodium_nt", "NTHR": "threonine_nt", "NTRP": "tryptophan_nt", "NTYR": "tyrosine_nt", "NVAL": "valine_nt", "A": "adenine", "C": "cytosine", "G": "guanine", "T": "thymine", "U": "uracil", "DA": "d-adenosine", "DC": "d-cytosine", "DG": "d-guanosine", "DT": "d-thymine", "RA": "r-adenosine", "RC": "r-cytosine", "RG": "r-guanosine", "RU": "r-uracil", "DA3": "d-adenosine_3ter", "DC3": "d-cytosine_3ter", "DG3": "d-guanosine_3ter", "DT3": "d-thymine_3ter", "RA3": "r-adenosine_3ter", "RC3": "r-cytosine_3ter", "RG3": "r-guanosine_3ter", "RU3": "r-uracil_3ter", "DA5": "d-adenosine_5ter", "DC5": "d-cytosine_5ter", "DG5": "d-guanosine_5ter", "DT5": "d-thymine_5ter", "RA5": "r-adenosine_5ter", "RC5": "r-cytosine_5ter", "RG5": "r-guanosine_5ter", "RU5": "r-uracil_5ter", "DAN": "d-adenosine_5ter_3ter", "DCN": "d-cytosine_5ter_3ter", "DGN": "d-guanosine_5ter_3ter", "DTN": "d-thymine_5ter_3ter", "RAN": "r-adenosine_5ter_3ter", "RCN": "r-cytosine_5ter_3ter", "RGN": "r-guanosine_5ter_3ter", "RUN": "r-uracil_5ter_3ter", } UnitedAtomResidueNameMap = { "ALA": "alanine_uni", "ARG": "arginine_uni", "ASP": "aspartic_acid_uni", "ASN": "asparagine_uni", "CYS": "cysteine_uni", "CYX": "cystine_ss_uni", "GLU": "glutamic_acid_uni", "GLN": "glutamine_uni", "GLY": "glycine_uni", "HIS": "histidine_uni", "HID": "histidine_deltah_uni", "HIE": "histidine_epsilonh_uni", "HIP": "histidine_plus_uni", "ILE": "isoleucine_uni", "LEU": "leucine_uni", "LYS": "lysine_uni", "MET": "methionine_uni", "PHE": "phenylalanine_uni", "PRO": "proline_uni", "SER": "serine_uni", "NA": "sodium_uni", "THR": "threonine_uni", "TRP": "tryptophan_uni", "TYR": "tyrosine_uni", "VAL": "valine_uni", "CALA": "alanine_ct_uni", "CARG": "arginine_ct_uni", "CASP": "aspartic_acid_ct_uni", "CASN": "asparagine_ct_uni", "CCYS": "cysteine_ct_uni", "CCYX": "cystine_ss_ct_uni", "CGLU": "glutamic_acid_ct_uni", "CGLN": "glutamine_ct_uni", "CGLY": "glycine_ct_uni", "CHIS": "histidine_ct_uni", "CHID": "histidine_deltah_ct_uni", "CHIE": "histidine_epsilonh_ct_uni", "CHIP": "histidine_plus_ct_uni", "CILE": "isoleucine_ct_uni", "CLEU": "leucine_ct_uni", "CLYS": "lysine_ct_uni", "CMET": "methionine_ct_uni", "CPHE": "phenylalanine_ct_uni", "CPRO": "proline_ct_uni", "CSER": "serine_ct_uni", "CNA": "sodium_ct_uni", "CTHR": "threonine_ct_uni", "CTRP": "tryptophan_ct_uni", "CTYR": "tyrosine_ct_uni", "CVAL": "valine_ct_uni", "NALA": "alanine_nt_uni", "NARG": "arginine_nt_uni", "NASP": "aspartic_acid_nt_uni", "NASN": "asparagine_nt_uni", "NCYS": "cysteine_nt_uni", "NCYX": "cystine_ss_nt_uni", "NGLU": "glutamic_acid_nt_uni", "NGLN": "glutamine_nt_uni", "NGLY": "glycine_nt_uni", "NHIS": "histidine_nt_uni", "NHID": "histidine_deltah_nt_uni", "NHIE": "histidine_epsilonh_nt_uni", "NHIP": "histidine_plus_nt_uni", "NILE": "isoleucine_nt_uni", "NLEU": "leucine_nt_uni", "NLYS": "lysine_nt_uni", "NMET": "methionine_nt_uni", "NPHE": "phenylalanine_nt_uni", "NPRO": "proline_nt_uni", "NSER": "serine_nt_uni", "NNA": "sodium_nt_uni", "NTHR": "threonine_nt_uni", "NTRP": "tryptophan_nt_uni", "NTYR": "tyrosine_nt_uni", "NVAL": "valine_nt_uni", } MoleculeNameMap = { "HOH": "water", "WAT": "water", } if __name__ == "__main__" or __name__ == "chimeraOpenSandbox": def minimize(mi): def update(mi): import chimera mi.saveMMTKCoordinates() chimera.runCommand("wait") mi.loadMMTKCoordinates() mi.minimize(nsteps=100, interval=10, action=update) def test(): import chimera # Standard residues only #mList = chimera.openModels.open("testdata/small2.pdb") #mList = chimera.openModels.open("testdata/1gcn.pdb") # Non-standard residues only #mList = chimera.openModels.open("testdata/gdp.pdb") # Both standard and non-standard residues mList = chimera.openModels.open("testdata/3fx2.pdb") mi = MMTKinter(mList, callback=minimize) test()