Scripts/MMTK: MMTK2Molecule.py

File MMTK2Molecule.py, 4.9 KB (added by Conrad Huang, 16 years ago)

Conversion code from MMTK universe to Chimera molecule

Line 
1def convert(universe, defaultCS=None, usePDBnames=True):
2 """Accepts an MMTK universe and returns (Chimera molecule,
3 map from MMTK atom to Chimera atom) tuple"""
4
5 c = _Converter(universe)
6 return c.convert(defaultCS=defaultCS, usePDBnames=usePDBnames)
7
8class _Converter:
9
10 def __init__(self, universe):
11 self.universe = universe
12
13 def convert(self, defaultCS, usePDBnames):
14 # First create the reference molecule
15 self.usePDBnames = usePDBnames
16 from chimera import Molecule
17 m = Molecule()
18 if defaultCS is not None:
19 m.activeCoordSet = m.newCoordSet(defaultCS)
20 self._wholeSeq = 1
21 self._mmtk2chimera = {} # For updating coordinates
22 for obj in self.universe.objectList():
23 if not obj.is_chemical_object:
24 continue
25 # Add the atoms
26 try:
27 resList = obj.residues()
28 except:
29 # No residues in the object, so just
30 # fake one for the entire object
31 self._addWholeMolecule(m, obj)
32 else:
33 # Add residues one by one
34 self._addResidues(m, resList)
35 # Add the bonds
36 self._addBonds(m, obj)
37 return m, self._mmtk2chimera
38
39 def _addWholeMolecule(self, m, obj):
40 # This code is untested. It "should" work.
41 r = m.newResidue("UNK", "het", self._wholeSeq, " ")
42 self._wholeSeq += 1
43 self._addAtoms(m, r, obj)
44
45 def _addResidues(self, m, resList):
46 # This code is untested. It "should" work.
47 # This code has only been tested with "bala1" which happens
48 # to be a protein whose atom names map well to PDB
49 # conventions. Changes will be needed for chains whose
50 # residues are not amino acids or nucleic bases.
51 for mr in resList:
52 seq = mr.sequence_number
53 rtype = mr.symbol
54 if self.usePDBnames:
55 rtype = rtype.upper()[:3]
56 if mr.parent.is_chain:
57 chainId = mr.parent.name
58 if self.usePDBnames:
59 chainId = chainId.replace("chain", "")
60 chainId = chainId.upper()[:1]
61 else:
62 chainId = ' '
63 r = m.newResidue(rtype, chainId, seq, " ")
64 self._addAtoms(m, r, mr)
65
66 def _addAtoms(self, m, r, obj):
67 from chimera import Element, Coord
68 # The MMTK atom names look like "C_beta" while Chimera really
69 # expects something like "CB". So we pull the conversion
70 # map out from the residue type's "pdbmap" attribute.
71 # The data structures are a little convoluted, but constructed
72 # "mmtk2pdb" is a dictionary whose keys are MMTK names
73 # and whose values are PDB names. If an MMTK name
74 # is not found in mmtk2pdb, we just use the MMTK name and
75 # hope for the best.
76 mmtk2pdb = {}
77 try:
78 if not self.usePDBnames:
79 raise AttributeError("skip PDB name code")
80 for _, atomMap in obj.type.pdbmap:
81 for name, atom in atomMap.iteritems():
82 mname = obj.type.atoms[atom.number].name
83 mmtk2pdb[mname] = name
84 except AttributeError:
85 pass
86 for ma in obj.atomList():
87 aname = mmtk2pdb.get(ma.name, ma.name)
88 e = Element(ma.type.symbol)
89 a = m.newAtom(aname, e)
90 r.addAtom(a)
91 x, y, z = ma.position() * 10
92 # Multiply coordinates by 10 to convert from MMTK nm
93 # to Chimera Angstrom
94 a.setCoord(Coord(x, y, z))
95 # Cache MMTK->Chimera atom mapping. We can do this
96 # using an attribute in ma, but I prefer to leave
97 # MMTK instances pristine
98 self._mmtk2chimera[ma] = a
99
100 def _addBonds(self, m, obj):
101 for mbu in obj.bondedUnits():
102 try:
103 bonds = mbu.bonds
104 except AttributeError:
105 continue
106 for mb in bonds:
107 try:
108 a1 = self._mmtk2chimera[mb.a1]
109 a2 = self._mmtk2chimera[mb.a2]
110 except KeyError:
111 pass
112 else:
113 b = m.newBond(a1, a2)
114
115def updateChimera(thread, universe, atomMap, frames, callback=None):
116 "Update Chimera coordinates from MMTK coordinates periodically."
117 from chimera import triggers
118 updater = _ChimeraUpdater(thread, universe, atomMap, frames, callback)
119 h = triggers.addHandler("check for changes", updater.update, None)
120 updater.setHandler(h)
121
122class _ChimeraUpdater:
123
124 def __init__(self, thread, universe, atomMap, frames, callback):
125 self.thread = thread
126 self.universe = universe
127 self.atomMap = atomMap
128 self.frames = frames
129 self.callback = callback
130 self.handler = None
131 self.curFrame = 0
132
133 def setHandler(self, h):
134 self.handler = h
135
136 def update(self, triggerName, myData, triggerData):
137 self.curFrame += 1
138 if self.curFrame < self.frames:
139 return
140 from chimera import replyobj
141 self.curFrame = 0
142 self._getCoords()
143 replyobj.status("MMTK updated coordinates", log=True)
144 if not self.thread.is_alive():
145 if self.handler is not None:
146 from chimera import triggers
147 triggers.deleteHandler("check for changes",
148 self.handler)
149 self.handler = None
150 replyobj.status("MMTK operation finished", log=True)
151
152 def _getCoords(self):
153 from chimera import Coord
154 state = self.thread.copyState()
155 if state is None:
156 conf = self.universe.copyConfiguration()
157 else:
158 conf = state["configuration"]
159 for ma in self.universe.atomList():
160 x, y, z = conf[ma] * 10
161 self.atomMap[ma].setCoord(Coord(x, y, z))
162 if self.callback:
163 self.callback(self.universe, self.atomMap, state)