Statistics
| Revision:

root / PyOpenGL-Demo / dek / OglSurface / NumericPDB.py @ 1

History | View | Annotate | Download (4.5 kB)

1
# This is statement is required by the build system to query build info
2
if __name__ == '__build__':
3
        raise Exception
4

    
5
import sys
6
import os
7
import string
8
try:
9
        import numpy as Numeric
10
except ImportError, err:
11
        try: 
12
                import Numeric
13
        except ImportError, err:
14
                print "This demo requires the numpy or Numeric extension, sorry"
15
                import sys
16
                sys.exit()
17
import Geometry
18
import copy
19

    
20
GUAAtoms = [ "N2", "O6", "C6", "C5", "N7", "C8", "N9", "C4", "N3", "C2", "N1" ]
21

    
22
def rotmat(phi, theta, psi, tx=0, ty=0, tz=0):
23

    
24
        s1 = Numeric.sin(phi)
25
        s2 = Numeric.sin(theta)
26
        s3 = Numeric.sin(psi)
27
        c1 = Numeric.cos(phi)
28
        c2 = Numeric.cos(theta)
29
        c3 = Numeric.cos(psi)
30

    
31
        newmat = Numeric.array([
32
                [ c2*c3, s2*s1*c3 - c1*s3, s2*c1*c3 + s1*s3, 0],
33
                [ c2*s3, s2*s1*s3 + c1*c3, s2*c1*s3 - s1*c3, 0],
34
                [-s2, c2*s1, c2*c1, 0],
35
                [tx, ty, tz, 1]])
36

    
37
        return newmat
38

    
39
def matrix_apply(crd, mat):
40
        o = [crd[0], crd[1], crd[2], 1.0]
41
        n = [0.0, 0.0, 0.0, 0.0]
42

    
43
        for i in range(4):
44
                for j in range(4):
45
                        n[i] = n[i] + o[j] * mat[j][i]
46
                        
47
        new = [None, None, None]
48
        
49
        for i in range(3):
50
                new[i] = (n[i]/n[3])
51
        return new
52

    
53

    
54

    
55
class PDBRecord:
56
        def __init__(self, type=None, anum=None, atom=None, residue=None, chain=None, rnum=None):
57
                self.type = type
58
                self.anum = anum
59
                self.atom = atom
60
                self.residue = residue
61
                self.chain = chain
62
                self.rnum = rnum
63

    
64

    
65

    
66
class PDB:
67
        def __init__(self, filename = None, crds = None, records = None, connect = None):
68
                if crds == None:
69
                        crds = []
70
                if records == None:
71
                        records = []
72
                self.records = records
73
                self.crds = crds
74
                self.connect = connect
75

    
76
                if filename != None:
77
                        sys.stderr.write("Reading in new PDB %s\n" % filename)
78
                        self.Read(filename)
79

    
80

    
81
## note: we read the anum here, which doesn't necessarily correspond to the actual
82
## record number.
83
        def Read(self, filename):
84
                pdbfile = open(filename)
85
                sys.stderr.write("Opened '%s' for reading as PDB\n" % filename)
86

    
87
                while(1):
88
                        line = pdbfile.readline()
89

    
90
                        if line == '':
91
                                break
92

    
93
                        if line[0:4] == 'ATOM' or line[0:6] == 'HETATM':
94
                                type = string.strip(line[0:6])
95
                                anum = string.atoi(string.strip(line[7:11]))
96
                                atom = string.strip(line[12:17])
97
                                residue = string.strip(line[17:20])
98
                                chain = string.strip(line[21:22])
99
                                rnum = string.atoi(string.strip(line[23:26]))
100
                                x = string.atof(string.strip(line[31:38]))
101
                                y = string.atof(string.strip(line[39:46]))
102
                                z = string.atof(string.strip(line[47:54]))
103

    
104
                                self.records.append(PDBRecord(type, anum, atom, residue, chain, rnum))
105
                                self.crds.append((x,y,z))
106
                                
107
                
108
                self.crds = Numeric.array(self.crds)
109

    
110
        def Write(self, filename):
111
                self.pdbfd=open(filename,"w")
112
                        
113
                for i in range(len(self.records)):
114
                        record = self.records[i]
115
                        self.pdbfd.write("%-6s%5d %-4s%c%-4s%c%4d%c   %8.3f%8.3f%8.3f\n" % \
116
                                                         ( record.type, record.anum, record.atom,' ', record.residue,' ', record.rnum, ' ', 
117
                                                           self.crds[i][0], self.crds[i][1], self.crds[i][2]))
118

    
119
                self.pdbfd.write("TER\n")
120

    
121
                if self.connect != None:
122
                        for i in self.connect:
123
                                self.pdbfd.write("CONECT")
124
                                for j in i:
125
                                        self.pdbfd.write("%5d" % (j+1))
126
                                self.pdbfd.write("\n")
127

    
128
                self.pdbfd.write("END\n")
129
                self.pdbfd.close()
130

    
131

    
132
        def Rotate(self,alpha, beta, gamma, tx, ty, tz):
133
                r = rotmat(alpha, beta, gamma, tx, ty, tz)
134
                self.RotateMatrix(r)
135

    
136
        def RotateMatrix(self,r):
137
                for i in range(len(self.crds)):
138
                        self.crds[i] = matrix_apply(self.crds[i], r)
139

    
140
        def Center(self):
141
                center = Numeric.add.reduce(self.crds) / len(self.crds)
142
                self.crds = Numeric.subtract(self.crds, center)
143

    
144
        def Print(self):
145
                for i in self.records:
146
                        print i.type, i.anum, i.atom, i.residue, i.chain, i.rnum
147

    
148
        def ReturnAnum(self, atom, rnum):
149
                for i in range(len(self.records)):
150
                        record = self.records[i]
151
                        if record.atom == atom and record.rnum == rnum:
152
                                 return i
153
                sys.stderr.write("Unable to find atom '%s' in residue %d\n" % (atom, rnum))
154
                return None
155

    
156
        def CrdByName(self, atom, res):
157
                x = self.ReturnAnum(atom, res)
158
                if x == None:
159
                        return None
160
                return self.crds[x]
161

    
162
        def FixResNum(self):
163
                rnum = 1
164
                for i in range(len(self.records)):
165
                        newrnum = rnum
166
                        if i < len(self.records)-1 and \
167
                           self.records[i+1].rnum != self.records[i].rnum:
168
                                rnum = rnum + 1
169
                        self.records[i].rnum = newrnum
170
                        self.records[i].anum = i
171

    
172
        def Copy(self):
173
                crds = Numeric.array(self.crds)
174
                records = copy.copy(self.records)
175
                return PDB(None, crds, records)
176

    
177

    
178
        def Append(self, struct2):
179
                records = self.records + copy.copy(struct2.records)
180
                crds = Numeric.concatenate((self.crds, struct2.crds))
181

    
182
                return PDB(None, crds, records)
183