File: MolecularSystem.py

package info (click to toggle)
mgltools-pyautodock 1.5.7-3
  • links: PTS, VCS
  • area: non-free
  • in suites: bullseye, buster, sid
  • size: 45,148 kB
  • sloc: python: 4,540; sh: 78; makefile: 13
file content (344 lines) | stat: -rw-r--r-- 13,068 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
############################################################################
#
# Authors: William Lindstrom, Ruth Huey
#
# Copyright: A. Olson TSRI 2004
#
#############################################################################

#
# $Id: MolecularSystem.py,v 1.9.6.1 2016/02/12 08:01:40 annao Exp $
#

from MolKit.molecule import AtomSet, Atom
from MolecularSystemAdapter import Adapters
import numpy
from math import sqrt
import types


class AbstractMolecularSystem:
    """This class defines the interface that a MolecularSystem must support
for use with a Scorer or other object(s). This class maintains a list of
sets of entities (like Atoms).

The get_supported_attributes method returns the list of attributes that
each Entity supports.

A distance matrix is maintained pairwise by set and pairwise by entities
within the set. NB: is actually the distance, **not** distance-squared.
The reason for this is that for clients (like WeightedMultiTerm scorers)
who must get_dist_mat repeatedly, forcing them to take the sqrt for each
term is a false economy.

An index is returned when adding a set of entities (by the add_entities method).
This index for an entity set can also be found using the get_entity_set_index method.
The molecular system attribute configuration is a 2-tuple of validated
entity_set_index. In its use, the Autodock scorer uses the first as the Receptor
and second as the Ligand.


This class is intended to abstract any number of underlying molecular
representations such as MolKit, mmLib, MMTK, or just a PDB file. The
underlying molecular representation is abstracted as a list of entities
which support the 'supported_attributes'. So, for instance, if your molecular
representation doesn't have charges, its abstraction as a specialized
implementation of this interface (ie the subclass of MolecularSystemBase)
may have to calculate per entity charges.
"""
    
    def __init__(self): pass
    def add_entities(self, entity_set, coords_access_func=None): pass
    def get_entity_set_index(self, entity_set): pass
    def set_coords(self, entity_set_index, coords): pass
    
    def get_dist_mat(self, ix, jx): pass
    def clear_dist_mat(self, ix): pass

    def get_supported_attributes(self): pass
    def check_required_attributes(self, entity_set_index, attribute_list): pass
    
    # @@ do we want to support the following interface??
    def get_coords(self, entity_set_index): pass
    
    # @@ or do we want to support the following interface??
    def get_entities(self, entity_set_num): pass

    def set_configuration(self, ix=None, jx=None ): pass
    

    
class MolecularSystem(AbstractMolecularSystem):
    """concrete subclass implementation of AbstractMolecularSystem

    """
    def __init__(self, e_set=None, coords_access_func=None):
        AbstractMolecularSystem.__init__(self)
        self.entity_sets = []
        self._dist_mats = {}
        # key is entity_set_index: values is dict
        # with key second entity_set_index, value distance mat
        # _dist_mat holds the square of distances between pairs of
        # entities, Distance matrices get computed only when needed
        self.configuration = (None, None)
        self.coords_access_func_list = []
        if e_set is not None:
            self.add_entities(e_set, coords_access_func)
            # single entity_set results in configuration of (0,0)
            # self.configuration = (0, None)


    def get_entity_set_index(self, entity_set):
        return self.entity_sets.index(entity_set)


    def _validate_entity_set_index(self, entity_set_ix):
        assert (entity_set_ix < len(self.entity_sets)), \
               "Invalid entity_set_index: %d; only %d entity_sets." % \
               (entity_set_ix, len(self.entity_sets) )

        
    def get_supported_attributes(self):
        raise NotImplementedError

    # for check_required_attributes
    # the hydrogen bonding term has the required attr of bonds
    # for subsetA, so when we implement this make sure that
    # 1. bonds is on the attribute_list when Hbonding is involved
    # 2. that subsetA does have bonds (or build them!)
    #
    # this also raises the issue that not all atom_sets need all attributes
    # to support a specific scorer term.
    #
    # also, for desolvation subsetA needs AtVols

    def check_required_attributes(self, entity_set_index, attribute_list):
        """Check that all the elements of the given entity_set have all
the attributes in the attribute_list.
        """
        # get_entities will validate entity_set_index
        entities = self.get_entities( entity_set_index)

        if type(entities) == types.GeneratorType:
            # @@ figure out how to check generator attributes !!
            return

        # check each attribute
        typeList = [types.ListType, types.DictType]#, types.InstanceType]
        from mglutil.util.misc import isInstance
        for attr in attribute_list:
            # get returns the subset of entities
            # for which the expression is true
            ents = entities.get( lambda x: hasattr(x, attr))
            if ents==None or len(ents)!=len(entities):
                raise AttributeError, \
                      "All Entities do not have required Entity attribute (%s) " % attr
            ent0 = entities[0]
            if type(getattr(ent0, attr)) in typeList \
              or isInstance(getattr(ent0, attr)) is True:
                if len(getattr(ent0, attr))==0:
                    raise ValueError, \
                          "Required Entity attribute (%s) has zero length" % attr


    def _compute_distance_matrix(self, ix, jx):
        """Compute the distance matrix for the given sets

This method should be called from within the class by get_dist_mat
which knows when it's really necessary.
        """
        # straight forward distance matrix calculation
        mat = []
        for c1 in self.get_coords(ix):
            l = []
            cx, cy, cz = c1
            mat.append(l)
            for c2 in self.get_coords(jx):
                c2x, c2y, c2z = c2
                d = cx-c2x, cy-c2y, cz-c2z
                l.append(sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2]) )
        return mat


    def check_distance_cutoff(self, ix, jx, cutoff):
        """check if a distance in the distance matrix for the specified sets
is below a user specified cutoff.

If the distance matrix does not exist if is computed up to the first encounter
of a distance below the cutoff. If such a distance is encountered it is
returned. If the matrix exists the minimum distance is returned if it is below
the cutoff. If distance are all above the cutoff, None is returned.
"""
        mat = self._dist_mats[ix][jx]
        if mat is None:
            mat = self._dist_mats[jx][ix]

        result = None
        if mat is None:
            mat = []
            for e1 in self.get_entities(ix):
                l = []
                cx, cy, cz = e1.coords
                mat.append(l)
                for e2 in self.get_entities(jx):
                    c2x, c2y, c2z = e2.coords
                    d = cx-c2x, cy-c2y, cz-c2z
                    dist = sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2])
                    if dist < cutoff:
                        return dist
                    l.append( dist)
            #if the whole dist matrix is ok, save it
            self._dist_mats[ix][jx] = mat
        else:
            mini = min( map( min, mat))
            if mini < cutoff:
                result = mini
        return result
    

    def get_dist_mat(self, ix, jx):
        """return the distance matrix for the given two sets
        """
        mat = self._dist_mats[ix][jx]
        tmat = self._dist_mats[jx][ix]
        
        if not mat:
            if tmat:
                # we just transpose
                #mat = numpy.swapaxes(numpy.array(tmat), 0, 1).tolist()
                mat = [] #eg: 3x4->4x3
                for j in xrange(len(tmat[0])):
                    mat.append([])
                    for i in xrange(len(tmat)):
                        mat[j].append(tmat[i][j])
            else:
                # we must compute
                mat = self._compute_distance_matrix(ix,jx)
            self._dist_mats[ix][jx] = mat
        return mat


    def clear_dist_mat(self, ix):
        """clear all distance-squared matrices for the given entity_set
        """
        
        # first clear all the pairs where entity_set is the first key
        for set in self._dist_mats[ix].keys():
            self._dist_mats[ix][set] = None

        # now clear all matrices where entity_set is the second key
        for set in self._dist_mats.keys():
            self._dist_mats[set][ix] = None


    def set_configuration(self, ix=None, jx=None ):
        """For now this is a 2-tuple of entity_set_index that describes
for a scorer which two entity sets are  being considered with respect to
one another.
        """
        if ix is not None:
            self._validate_entity_set_index(ix)
            self.configuration = tuple( [ix, self.configuration[1]] )
        if jx is not None:
            self._validate_entity_set_index(jx)
            self.configuration = tuple( [self.configuration[0], jx] )

    
    def add_entities(self, entity_set, coords_access_func=None):
        """add a MolKit.EntitySet to the MolecularSystem
        
Returns the entity_set_num for the new EntitySet. First two calls to
add_entities set the configuration to use those EntitySets.
        """
        entity_set_ix =  len(self.entity_sets)

        # initialize the dictionary of distance matricies for this set
        self._dist_mats[entity_set_ix] = {}

        # create an adapter for this entity_set ...
        adapter = Adapters[str(entity_set.__class__)](entity_set)
#        adapter = MolecularSystemAdapter( entity_set)
        # ... and add it to the entity_sets list
        self.entity_sets.append(adapter)
        
        # add the coords_access_func
        self.coords_access_func_list.append(coords_access_func)
        
        # create a new dict for every set vs. this one
        for ix, set in enumerate(self.entity_sets):
            self._dist_mats[ix][entity_set_ix] = None
            self._dist_mats[entity_set_ix][ix] = None

        # set the configuration (only on the first two calls)
        if entity_set_ix == 0:
            self.set_configuration( ix = entity_set_ix)
            # single entity_set results in configuration of (0,0)
            self.set_configuration( jx = entity_set_ix)
        elif entity_set_ix == 1:
            self.set_configuration( jx = entity_set_ix)
            
        return entity_set_ix


    def get_entities(self, entity_set_ix=0):
        """return the specified entity_set.

entity_set_ix is an integer that indexes into self.entity_sets"""
        # make sure the entity_set_index is valid
        self._validate_entity_set_index(entity_set_ix)
        return self.entity_sets[entity_set_ix].get_iterator()



    # @@ ISSUE: should the MolecularSystem maintain the coordinates
    # @@ as state. OR, should this actually change the values in whatever
    # @@ underlying Molecule representation that there is ??
    #
    # for now set_coords changes the MolKit values
    
    def set_coords(self, entity_set_ix, coords, coords_set_func=None):
        """provide new coords for the specified entitySet.
entity_set_ix is the index of an existing entitySet.
coords is a list (of proper length) of anything you want to use
to update coordinate with the coords_set_func you supply.
        """
        # validate input
        self._validate_entity_set_index(entity_set_ix)
        # assert <that coords is a list of coordinates>
        # check compatibility
        assert len(coords)==len(self.get_entities(entity_set_ix))

        if coords_set_func:
            for ex, e in enumerate(self.get_entities(entity_set_ix)):
                coords_set_func(e, coords[ex])
        else:
            #NB: this is MolKit AtomSet specific
            self.get_entities(entity_set_ix).updateCoords(coords)

        # clear distance matrix
        self.clear_dist_mat(entity_set_ix)


    def get_coords(self, entity_set_ix=0):
        """return the coordinates of the specified entity_set.

        DEPRECATED: use get_entities().coords instead.
        entity_set_ix is an integer that indexes into self.entity_sets
        """
        # raise DeprecationWarning, "use get_entities().coords instead"
        # validate input
        self._validate_entity_set_index(entity_set_ix)
        
        # get the coordinates accessor function
        accessor = self.coords_access_func_list[entity_set_ix]
        if accessor:
            coords = [accessor(e) for e in self.get_entities(entity_set_ix)]
        else:
            #in case getattr hasnot been overridden:
            coords = [e.coords for e in self.get_entities(entity_set_ix)] 
        return coords

# MolecularSystem

if __name__ == '__main__':
    print "unittests in Tests/test_MolecularSystem.py"