hdiff output

r29059/genrigid_input.py 2017-01-21 10:36:19.622250000 +0000 r29058/genrigid_input.py 2017-01-21 10:36:20.158250000 +0000
  1: # SCRIPT TO DEFINE THE RIGIDIFICATION FOR PROTEINS 
  2: #by Konstantin Röder 
  3: #creates input files rbodyconfig and coordsinirigid and additionally a log file rigid.log 
  4:  
  5: #! /usr/bin/env python 
  6:  
  7: import sys 
  8: import time 
  9: import math 
 10:  
 11: if len(sys.argv) > 1: 
 12:     print 'Script creates input files for RIGIDINIT keyword, needs inpcrd and pdb' 
 13:     print 
 14:     pdb_inp = sys.argv[1] 
 15:     coords_inp = sys.argv[2] 
 16:     #defines the tolerance for the checks whether rigid bodies are linear 
 17:     try: 
 18:         tolerance = float(sys.argv[3]) 
 19:     except IndexError: 
 20:         tolerance = 0.01 
 21:     #Pymol use; change to False if pymol is not installed/used 
 22:     try: 
 23:         if sys.argv[4] == 'pymol': 
 24:             pymol_check = True 
 25:     except IndexError: 
 26:         pymol_check = False 
 27:  
 28: else: 
 29:     tolerance = 0.01 
 30:     pymol_check = False 
 31:     pdb_inp = raw_input('PDB file: ') 
 32:     coords_inp = raw_input('Coords input: ') 
 33:  
 34: if pymol_check: 
 35:     import __main__ 
 36:  
 37:     __main__.pymol_argv = ['pymol', '-qix']  # Pymol: quiet and no GUI(internal and external) 
 38:     import pymol 
 39:  
 40: #class containing the main methods for the script, new functionality should be embedded here 
 41: class protein(): 
 42:     def __init__(self, atom_dic, res_dic): 
 43:         self.atom_dic = atom_dic  #dictionary of all atoms with x,y,z coordinates and res id and name and atom name 
 44:         self.res_dic = res_dic    #dictionary containing a list of atom ids for all residues 
 45:  
 46:     def num_atoms(self): #number of atoms 
 47:         return len(self.atom_dic) 
 48:  
 49:     def num_res(self): #number of residues 
 50:         return len(self.res_dic) 
 51:  
 52:     #returns the name of a given residue 
 53:     def get_residue_name(self, residue): 
 54:         return atom_dic[res_dic[residue][0]][1] 
 55:  
 56:     #prints the atoms in a specified residue(need number of residue) 
 57:     #includes the atom number, atom name and coordinates 
 58:     def get_atoms_in_res(self, residue): 
 59:         all_atoms = self.atom_dic.keys() 
 60:         atoms_in_res = [] 
 61:         for i in range(1, len(self.atom_dic) + 1): 
 62:             atom = self.atom_dic[i] 
 63:             if residue == atom[2]: 
 64:                 atoms_in_res.append(all_atoms[i - 1]) 
 65:                 residuename = atom[1] 
 66:         if atoms_in_res == []: 
 67:             return 'Residue not in molecule' 
 68:         print 'Residue:', residue, residuename 
 69:         print 
 70:         print 'Atoms:' 
 71:         for i in range(0, len(atoms_in_res)): 
 72:             j = atoms_in_res[i] 
 73:             atom = self.atom_dic[j] 
 74:             print j, atom[0], atom[3], atom[4], atom[5] 
 75:         return 
 76:  
 77:     #returns a list of atoms in a residue 
 78:     def get_atom_list(self, residue): 
 79:         all_atoms = self.atom_dic.keys() 
 80:         atomlist = [] 
 81:         for i in range(1, len(self.atom_dic) + 1): 
 82:             atom = self.atom_dic[i] 
 83:             if residue == atom[2]: 
 84:                 atomlist.append(all_atoms[i - 1]) 
 85:         return atomlist 
 86:  
 87:     #allows the input of a centre from the user, is complemented by the selection of the CoM as centre 
 88:     def choose_centre(self): 
 89:         check = raw_input('Display atoms in residue (y/n)?  ') 
 90:         if check == 'y': 
 91:             while 1:  #acts as input mask --> once entered any number of residues can be viewed at 
 92:                 residue_number = raw_input('Residue number: (enter x to leave) ') 
 93:                 if residue_number == 'x': 
 94:                     print 
 95:                     break 
 96:                 else: 
 97:                     try: 
 98:                         self.get_atoms_in_res(int(residue_number)) 
 99:                         print 
100:                     except ValueError: 
101:                         pass 
102:         while 1:  # loop for choice of molecule, only accepts atoms in the molecule 
103:             centre_input = int(raw_input('Choose atom:  ')) 
104:             if centre_input > self.num_atoms() or centre_input < 1: 
105:                 print 'Not in molecule.' 
106:             else: 
107:                 break 
108:         return centre_input 
109:  
110:     #transforms a given list of atoms into a list of residues 
111:     #can be used to find the residues in a sphere from the output of the atoms in sphere method 
112:     def transform_atomlist_to_reslist(self, atomlist): 
113:         reslist = [] 
114:         for i in atomlist: 
115:             atom_info = self.atom_dic[i] 
116:             k = 0 
117:             if reslist != []: 
118:                 for j in range(0, len(reslist)): 
119:                     if reslist[j] == atom_info[2]: 
120:                         k = 1 
121:             if k == 0: 
122:                 reslist.append(atom_info[2]) 
123:             else: 
124:                 pass 
125:         return reslist 
126:  
127:     def transform_reslist_to_atomlist(self, reslist): 
128:         atomlist = [] 
129:         for i in reslist: 
130:             atomlist += self.res_dic[i] 
131:         return atomlist 
132:  
133:     #get residue name from id 
134:     def res_name_from_id(self, residue_id): 
135:         atom_id = res_dic[residue_id][0] 
136:         return atom_dic[atom_id][1] 
137:  
138:     #checks for unknown residues/ligands, add standard residues here if needed 
139:     def get_ligands_unknown_res(self): 
140:         list_res_names = ['ALA', 'ARG', 'ASN', 'ASP', 'ASX', 'CYS', 'CYX', 'GLN', 'GLU', 'GLY', 'GLX', 'HIS', 'HIE', 
141:                           'HID', 'HIP', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 
142:                           'NALA', 'NARG', 'NASN', 'NASP', 'NASX', 'NCYS', 'NCYX', 'NGLN', 'NGLU', 'NGLY', 'NGLX', 
143:                           'NHIS', 'NHIE', 'NHID', 'NHIP', 'NILE', 'NLEU', 'NLYS', 'NMET', 'NPHE', 'NPRO', 'NSER', 
144:                           'NTHR', 'NTRP', 'NTYR', 'NVAL', 'CALA', 'CARG', 'CASN', 'CASP', 'CASX', 'CCYS', 'CCYX', 
145:                           'CGLN', 'CGLU', 'CGLY', 'CGLX', 'CHIS', 'CHIE', 'CHID', 'CHIP', 'CILE', 'CLEU', 'CLYS', 
146:                           'CMET', 'CPHE', 'CPRO', 'CSER', 'CTHR', 'CTRP', 'CTYR', 'CVAL'] 
147:         ligand_dic = {} 
148:         for i in self.res_dic.keys(): 
149:             if self.res_name_from_id(i) not in list_res_names: 
150:                 ligand_dic[i] = self.res_dic[i] 
151:         return ligand_dic 
152:  
153:     #returns a list of atoms within a sphere of entered radius 
154:     #needs an atom number as centre 
155:     def get_atoms_in_sphere(self, centre_atom, radius): 
156:         all_atoms = self.atom_dic.keys() 
157:         atomlist = [] 
158:         for i in range(1, len(self.atom_dic) + 1): 
159:             atom = self.atom_dic[i] 
160:             distance = (atom[3] - centre_atom[0]) ** 2 + (atom[4] - centre_atom[1]) ** 2 + (atom[5] - centre_atom[ 
161:                 2]) ** 2 
162:             if distance <= radius: 
163:                 atomlist.append(all_atoms[i - 1]) 
164:             else: 
165:                 pass 
166:         return atomlist 
167:  
168:     #returns the atoms within an ellipsoidal volume 
169:     def get_atoms_in_ellipsoid(self, centre_atom, a, b, c): 
170:         all_atoms = self.atom_dic.keys() 
171:         atomlist = [] 
172:         for i in range(1, len(self.atom_dic) + 1): 
173:             atom = self.atom_dic[i] 
174:             x1 = ((atom[3] - centre_atom[0]) / a) ** 2 
175:             x2 = ((atom[4] - centre_atom[1]) / b) ** 2 
176:             x3 = ((atom[5] - centre_atom[2]) / c) ** 2 
177:             distance = x1 + x2 + x3 
178:             if distance <= 1: 
179:                 atomlist.append(all_atoms[i - 1]) 
180:             else: 
181:                 pass 
182:         return atomlist 
183:  
184:     #returns a list of atoms in a residue 
185:     def get_atom_list(self, residue): 
186:         all_atoms = self.atom_dic.keys() 
187:         atomlist = [] 
188:         for i in range(1, len(self.atom_dic) + 1): 
189:             atom = self.atom_dic[i] 
190:             if residue == atom[2]: 
191:                 atomlist.append(all_atoms[i - 1]) 
192:         return atomlist 
193:  
194:     #gives the mass of any residue containing only C, N, O, H, D and S, other masses can only be used if they are added to the dictionary mass_dic below 
195:     def res_mass(self, residue): 
196:         atom_list = self.get_atom_list(residue) 
197:         mass_dic = {'H': 1.007825, 'D': 2.014102, 'C': 12.0116, 'N': 14.00728, 'O': 15.99977, 'S': 32.076} 
198:         mass = 0 
199:         for i in atom_list: 
200:             for j in mass_dic.keys(): 
201:                 if j == self.atom_dic[i][0][0]: 
202:                     mass += mass_dic[j] 
203:         return mass 
204:  
205:     #gives the coordinates for the mass_weighted centre of a residue 
206:     def mass_weighted_centre(self, residue): 
207:         atom_list = self.get_atom_list(residue) 
208:         mass_dic = {'H': 1.007825, 'D': 2.014102, 'C': 12.0116, 'N': 14.00728, 'O': 15.99977, 'S': 32.076} 
209:         X = 0 
210:         Y = 0 
211:         Z = 0 
212:         mass_res = float(self.res_mass(residue)) 
213:         for i in atom_list: 
214:             for j in mass_dic.keys(): 
215:                 if j == self.atom_dic[i][0][0]: 
216:                     X += mass_dic[j] * self.atom_dic[i][3] 
217:         for i in atom_list: 
218:             for j in mass_dic.keys(): 
219:                 if j == self.atom_dic[i][0][0]: 
220:                     Y += mass_dic[j] * self.atom_dic[i][4] 
221:         for i in atom_list: 
222:             for j in mass_dic.keys(): 
223:                 if j == self.atom_dic[i][0][0]: 
224:                     Z += mass_dic[j] * self.atom_dic[i][5] 
225:         X = X / mass_res 
226:         Y = Y / mass_res 
227:         Z = Z / mass_res 
228:         print 'Mass weighted centre: ', (X, Y, Z) 
229:         return (X, Y, Z) 
230:  
231:  
232:     #returns all residues of a list of type 'name' 
233:     def get_all_res_name(self, res_input, name): 
234:         res_list = [] 
235:         for i in res_input: 
236:             if name == self.get_residue_name(i): 
237:                 res_list.append(i) 
238:         return res_list 
239:  
240:     #atom id from atom name given a residue id 
241:     def get_atom_id_from_name(self, atom_name, res_num): 
242:         for i in self.res_dic[res_num]: 
243:             if self.atom_dic[i][0] == atom_name: 
244:                 return i 
245:  
246:     #dot product of two vectors 
247:     def dot_product(self, vector1, vector2): 
248:         return (vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2]) 
249:  
250:     #magnitude of a vector 
251:     def magnitude(self, vector): 
252:         return math.sqrt(vector[0] ** 2 + vector[1] ** 2 + vector[2] ** 2) 
253:  
254:     #angle between three atoms 
255:     def get_angle(self, (atom1, atom2, atom3)): 
256:         vector1 = (self.atom_dic[atom2][3] - self.atom_dic[atom1][3], self.atom_dic[atom2][4] - self.atom_dic[atom1][4] 
257:                    , self.atom_dic[atom2][5] - self.atom_dic[atom1][5]) 
258:         vector2 = (self.atom_dic[atom3][3] - self.atom_dic[atom1][3], self.atom_dic[atom3][4] - self.atom_dic[atom1][4] 
259:                    , self.atom_dic[atom3][5] - self.atom_dic[atom1][5]) 
260:         dot_product = self.dot_product(vector1, vector2) 
261:         magnitudes = self.magnitude(vector1) * self.magnitude(vector2) 
262:         angle = math.acos(dot_product / magnitudes) 
263:         angle = math.degrees(angle) 
264:         if (angle > 0.00 and angle <= 180.00): 
265:             return angle 
266:         else: 
267:             return angle - 180.00 
268:  
269:     #checks a list of atoms whether it is linear, creates a triple list of all possible combinations and then computes angles 
270:     def check_linear(self, atom_list, tolerance): 
271:         triple_list = []  # list of atom triples to check for linearity 
272:         for atom1 in atom_list: 
273:             for atom2 in atom_list: 
274:                 for atom3 in atom_list: 
275:                     if (atom1 != atom2) and (atom1 != atom3) and (atom2 != atom3): 
276:                         triple_list.append((atom1, atom2, atom3)) 
277:                     else: 
278:                         pass 
279:         for triple in triple_list: 
280:             if (self.get_angle(triple) > (0.00 + tolerance)) and (self.get_angle(triple) < (180.00 - tolerance)): 
281:                 return False 
282:         return True 
283:  
284:     #returns a dictionary including all the information needed for the localised rigidification scheme, 
285:     #if an extension is made to the default options for rigidification, this must change the local_scheme list! 
286:     #the input scheme must be a list of 0,1,... for the possible options within the rigidification for a given residue 
287:     #any additional scheme can be added below that will be available on default if chosen 
288:     def get_rigid_for_local(self, res_list, local_scheme, atoms_used): 
289:         #local_scheme=['PRO','ARG',('HIS','HIE','HID','HIP'),'LYS','ASP','ASN','GLU','GLN','PHE','TYR','TRP'] 
290:         groups_res = {} 
291:         i = 1 
292:         #all following cases are for a single residue name with given patterns 
293:         #when adding additional ones check, that atoms can only be chosen once as no atom can be part of more than one rigid body 
294:         if local_scheme[0] == 1: 
295:             l_res = self.get_all_res_name(res_list, 'PRO') 
296:             for j in l_res: 
297:                 atom_list = [] 
298:                 for k in ['C', 'O', 'N', 'CD', 'CG', 'CB', 'CA']: 
299:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
300:                 groups_res[i] = [j, 'PRO'] + atom_list 
301:                 atoms_used += atom_list 
302:                 i += 1 
303:             l_res = [] 
304:         if local_scheme[0] == 2: 
305:             l_res = self.get_all_res_name(res_list, 'PRO') 
306:             for j in l_res: 
307:                 atom_list = [] 
308:                 for k in ['N', 'CD', 'CG', 'CB', 'CA']: 
309:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
310:                 groups_res[i] = [j, 'PRO'] + atom_list 
311:                 atoms_used += atom_list 
312:                 i += 1 
313:             l_res = [] 
314:         if local_scheme[1] == 1: 
315:             l_res = self.get_all_res_name(res_list, 'ARG') 
316:             for j in l_res: 
317:                 atom_list = [] 
318:                 for k in ['NE', 'HE', 'CZ', 'NH1', 'HH11', 'HH12', 'NH2', 'HH21', 'HH22']: 
319:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
320:                 groups_res[i] = [j, 'ARG'] + atom_list 
321:                 atoms_used += atom_list 
322:                 i += 1 
323:             l_res = [] 
324:         if local_scheme[2] == 1: 
325:             l_res = self.get_all_res_name(res_list, 'HIS') 
326:             for j in l_res: 
327:                 atom_list = [] 
328:                 for k in ['CG', 'ND1', 'CE1', 'NE2', 'CD2']: 
329:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
330:                 groups_res[i] = [j, 'HIS'] + atom_list 
331:                 atoms_used += atom_list 
332:                 i += 1 
333:             l_res = [] 
334:             l_res = self.get_all_res_name(res_list, 'HIE') 
335:             for j in l_res: 
336:                 atom_list = [] 
337:                 for k in ['CG', 'ND1', 'CE1', 'NE2', 'CD2']: 
338:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
339:                 groups_res[i] = [j, 'HIE'] + atom_list 
340:                 atoms_used += atom_list 
341:                 i += 1 
342:             l_res = [] 
343:             l_res = self.get_all_res_name(res_list, 'HID') 
344:             for j in l_res: 
345:                 atom_list = [] 
346:                 for k in ['CG', 'ND1', 'CE1', 'NE2', 'CD2']: 
347:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
348:                 groups_res[i] = [j, 'HID'] + atom_list 
349:                 atoms_used += atom_list 
350:                 i += 1 
351:             l_res = [] 
352:             l_res = self.get_all_res_name(res_list, 'HIP') 
353:             for j in l_res: 
354:                 atom_list = [] 
355:                 for k in ['CG', 'ND1', 'CE1', 'NE2', 'CD2']: 
356:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
357:                 groups_res[i] = [j, 'HIP'] + atom_list 
358:                 atoms_used += atom_list 
359:                 i += 1 
360:             l_res = [] 
361:         if local_scheme[3] == 1: 
362:             l_res = self.get_all_res_name(res_list, 'LYS') 
363:             for j in l_res: 
364:                 atom_list = [] 
365:                 for k in ['NZ', 'HZ1', 'HZ2', 'HZ3']: 
366:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
367:                 groups_res[i] = [j, 'LYS'] + atom_list 
368:                 atoms_used += atom_list 
369:                 i += 1 
370:             l_res = [] 
371:         if local_scheme[4] == 1: 
372:             l_res = self.get_all_res_name(res_list, 'ASP') 
373:             for j in l_res: 
374:                 atom_list = [] 
375:                 for k in ['CB', 'CG', 'OD1', 'OD2']: 
376:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
377:                 groups_res[i] = [j, 'ASP'] + atom_list 
378:                 atoms_used += atom_list 
379:                 i += 1 
380:             l_res = [] 
381:         if local_scheme[5] == 1: 
382:             l_res = self.get_all_res_name(res_list, 'ASN') 
383:             for j in l_res: 
384:                 atom_list = [] 
385:                 for k in ['CB', 'CG', 'OD1', 'ND2', 'HD21', 'HD22']: 
386:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
387:                 groups_res[i] = [j, 'ASN'] + atom_list 
388:                 atoms_used += atom_list 
389:                 i += 1 
390:             l_res = [] 
391:         if local_scheme[6] == 1: 
392:             l_res = self.get_all_res_name(res_list, 'GLU') 
393:             for j in l_res: 
394:                 atom_list = [] 
395:                 for k in ['CG', 'CD', 'OE1', 'OE2']: 
396:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
397:                 groups_res[i] = [j, 'GLU'] + atom_list 
398:                 atoms_used += atom_list 
399:                 i += 1 
400:             l_res = [] 
401:         if local_scheme[7] == 1: 
402:             l_res = self.get_all_res_name(res_list, 'GLN') 
403:             for j in l_res: 
404:                 atom_list = [] 
405:                 for k in ['CG', 'CD', 'OE1', 'NE2', 'HE21', 'HE22']: 
406:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
407:                 groups_res[i] = [j, 'GLN'] + atom_list 
408:                 atoms_used += atom_list 
409:                 i += 1 
410:             l_res = [] 
411:         if local_scheme[8] == 1: 
412:             l_res = self.get_all_res_name(res_list, 'PHE') 
413:             for j in l_res: 
414:                 atom_list = [] 
415:                 for k in ['CG', 'CD1', 'HD1', 'CE1', 'HE1', 'CZ', 'HZ', 'CE2', 'HE2', 'CD2', 'HD2']: 
416:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
417:                 groups_res[i] = [j, 'PHE'] + atom_list 
418:                 atoms_used += atom_list 
419:                 i += 1 
420:             l_res = [] 
421:         if local_scheme[9] == 1: 
422:             l_res = self.get_all_res_name(res_list, 'TYR') 
423:             for j in l_res: 
424:                 atom_list = [] 
425:                 for k in ['CG', 'CD1', 'HD1', 'CE1', 'HE1', 'CZ', 'OH', 'CE2', 'HE2', 'CD2', 'HD2']: 
426:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
427:                 groups_res[i] = [j, 'TYR'] + atom_list 
428:                 atoms_used += atom_list 
429:                 i += 1 
430:             l_res = [] 
431:         if local_scheme[10] == 1: 
432:             l_res = self.get_all_res_name(res_list, 'TRP') 
433:             for j in l_res: 
434:                 atom_list = [] 
435:                 for k in ['CG', 'CD1', 'HD1', 'NE1', 'HE1', 'CE2', 'CZ2', 'HZ2', 'CH2', 'HH2', 'CZ3', 'CE3', 'HE3', 
436:                           'CD2']: 
437:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
438:                 groups_res[i] = [j, 'TRP'] + atom_list 
439:                 atoms_used += atom_list 
440:                 i += 1 
441:             l_res = [] 
442:         if local_scheme[10] == 2: 
443:             l_res = self.get_all_res_name(res_list, 'TRP') 
444:             for j in l_res: 
445:                 atom_list1 = [] 
446:                 atom_list2 = [] 
447:                 for k in ['CG', 'CD1', 'HD1', 'NE1', 'HE1']: 
448:                     atom_list1.append(self.get_atom_id_from_name(k, int(j))) 
449:                 for k in ['CE2', 'CZ2', 'HZ2', 'CH2', 'HH2', 'CZ3', 'CE3', 'HE3', 'CD2']: 
450:                     atom_list2.append(self.get_atom_id_from_name(k, int(j))) 
451:                 groups_res[i] = [j, 'TRP'] + atom_list1 
452:                 groups_res[i + 1] = [j, 'TRP'] + atom_list2 
453:                 atoms_used += atom_list1 
454:                 atoms_used += atom_list2 
455:                 i += 2 
456:             l_res = [] 
457:         return (groups_res, atoms_used) 
458:  
459:  
460: #class to implement a live feed in pymol - not yet tested 
461: class PymolView(): 
462:     #initialises pymol and starts it without GUI, loads a given file and sets it to a white cartoon representation 
463:     def __init__(self, file_name): 
464:         print 'Launching Pymol' 
465:         pymol.finish_launching() 
466:         pymol.cmd.load(file_name) 
467:         print '%s loaded in Pymol' % file_name 
468:  
469:         #show molecule as cartoon and hide lines 
470:         pymol.cmd.show('cartoon') 
471:         pymol.cmd.hide('lines') 
472:         pymol.cmd.set('cartoon_color', 0) 
473:  
474:     def set_colour_range(self, start, finish, colour): 
475:         #setting residue range to any colour 
476:         pymol.cmd.set('cartoon_color', colour, 'resi %s-%s' % (str(start), str(finish))) 
477:  
478:     #reset all of the molecule to a white colour 
479:     #def reset_colour(self,length): 
480:     #    pymol.cmd.set('cartoon_color',0,'resi 1-%s'%str(length)) 
481:  
482:     def set_colour_res(self, res, colour): 
483:         pymol.cmd.set('cartoon_color', colour, 'resi %s' % str(res)) 
484:  
485:     def __del__(self): 
486:         print 'Quit pymol' 
487:         pymol.cmd.quit() 
488:  
489:     def apply_colour_scheme(self, atomistic_list, local_list, colour_atomistic, colour_local): 
490:         #self.reset_colour() 
491:         for i in atomistic_list: 
492:             self.set_colour_res(i, colour_atomistic) 
493:         for j in local_list: 
494:             self.set_colour_res(j, colour_local) 
495:  
496:  
497: ###FUNCTION TO PARSE THE PDB AND COORDS.INPCRD FILE ################################################################### 
498: def reading_pdb_inpcrd(pdb, inpcrd): 
499:     dic = {} 
500:     try: 
501:         atom_list = [] 
502:         with open(pdb, 'r') as f_pdb: 
503:             for line in f_pdb: 
504:                 if line.split()[0] == 'ATOM': 
505:                     atom_list.append((line[12:16].strip(), line[17:20].strip(), int(line[22:26]))) 
506:         with open(inpcrd, 'r') as f_crd: 
507:             for _ in xrange(2): 
508:                 next(f_crd) 
509:             i = 1 
510:             for line in f_crd: 
511:                 atom_line = line.split() 
512:                 try: 
513:                     atom_list[2 * i - 2] = atom_list[2 * i - 2] + ( 
514:                         float(atom_line[0]), float(atom_line[1]), float(atom_line[2])) 
515:                 except IndexError: 
516:                     pass 
517:                 try: 
518:                     atom_list[2 * i - 1] = atom_list[2 * i - 1] + ( 
519:                         float(atom_line[3]), float(atom_line[4]), float(atom_line[5])) 
520:                     i = i + 1 
521:                 except IndexError: 
522:                     pass 
523:         for i in range(1, len(atom_list) + 1): 
524:             dic[i] = atom_list[i - 1] 
525:         return dic 
526:     except IOError: 
527:         return dic 
528:  
529:  
530: def create_res_dic(atom_dic): 
531:     res_dic = {} 
532:     res_list = [] 
533:     for i in atom_dic.keys(): 
534:         res_list.append(atom_dic[i][2]) 
535:     res_list = list(set(res_list)) 
536:     for j in res_list: 
537:         res_dic[j] = [] 
538:     for k in atom_dic.keys(): 
539:         res_dic[atom_dic[k][2]].append(k) 
540:     return res_dic 
541:  
542:  
543: ###FUNCTIONS HELPING TO MANAGE THE LISTS AND DICTIONARY ############################################################### 
544: def merge_lists(list_a, list_b): 
545:     return sorted(list(set(list_a + list_b))) 
546:  
547:  
548: def remove_selection_from_list(list_res, selection): 
549:     new_res_list = [] 
550:     for i in list_res: 
551:         if i not in selection: 
552:             new_res_list.append(i) 
553:     return sorted(new_res_list) 
554:  
555:  
556: def check_res_exists(res_list, selection): 
557:     new_selection = [] 
558:     for i in selection: 
559:         if i in res_list: 
560:             new_selection.append(i) 
561:     return new_selection 
562:  
563:  
564: def merge_new_and_old_res_lists(frozen, local, atomistic, sel_local, sel_atomistic): 
565:     new_sel_atomistic = check_res_exists(frozen, sel_atomistic) 
566:     new_sel_atomistic_from_local = check_res_exists(local, sel_atomistic) 
567:     frozen = remove_selection_from_list(frozen, new_sel_atomistic) 
568:     atomistic = merge_lists(atomistic, new_sel_atomistic) 
569:     atomistic = merge_lists(atomistic, new_sel_atomistic_from_local) 
570:     local = remove_selection_from_list(local, new_sel_atomistic_from_local) 
571:     new_sel_local = check_res_exists(frozen, sel_local) 
572:     local = merge_lists(local, new_sel_local) 
573:     frozen = remove_selection_from_list(frozen, new_sel_local) 
574:     return (sorted(atomistic), sorted(local), sorted(frozen)) 
575:  
576:  
577: def removal_selection(frozen, local, atomistic, sel_local, sel_atomistic): 
578:     frozen = merge_lists(frozen, sel_local) 
579:     frozen = merge_lists(frozen, sel_atomistic) 
580:     local = remove_selection_from_list(local, sel_local) 
581:     atomistic = remove_selection_from_list(atomistic, sel_atomistic) 
582:     return (sorted(atomistic), sorted(local), sorted(frozen)) 
583:  
584:  
585: ###FUNCTION TO SELECT RESIDUES ######################################################################################## 
586: def choose_res(mol): 
587:     print '1 - choose sphere/ellipsoid of unfrozen atoms and add localised rigidification' 
588:     print '2 - choose a range of residues or single residues to unfreeze or locally rigidify' 
589:     print 
590:     log = [] 
591:     try: 
592:         method = int(raw_input('Method:  ')) 
593:     except IOError: 
594:         method = 0 
595:     unfrozen_res = [] 
596:     local_res = [] 
597:  
598:     if method == 0: 
599:         print 'Invalid input' 
600:  
601:     # ellipsoidal or spherical volume 
602:     elif method == 1: 
603:         while 1:  #choose the geometry: currently ellipsoidal 'e' or spherical 's' 
604:             log.append('Method: sphere / ellipsoid') 
605:             while 1: 
606:                 centre_method = raw_input('Mass weighted centre of residue (1) or atom in a selected residue (2)?  ') 
607:                 try: 
608:                     if int(centre_method) == 1: 
609:                         residue = raw_input('Residue:  ') 
610:                         centre_xyz = mol.mass_weighted_centre(int(residue)) 
611:                         log.append('Mass-weighted centre: Residue' + residue + ' ' + mol.get_residue_name( 
612:                             int(residue)) + '  ' + str(centre_xyz[0]) + '  ' + str(centre_xyz[1]) + '  ' + str( 
613:                             centre_xyz[2])) 
614:                         break 
615:                     elif int(centre_method) == 2: 
616:                         centre = mol.choose_centre() 
617:                         centre_xyz = (mol.atom_dic[centre][3], mol.atom_dic[centre][4], mol.atom_dic[centre][5]) 
618:                         log.append('Centre: Atom  ' + str(centre) + '  ' + str(centre_xyz[0]) + '  ' + str( 
619:                             centre_xyz[1]) + '  ' + str(centre_xyz[2])) 
620:                         break 
621:                     else: 
622:                         print 'Invalid choice' 
623:                 except ValueError: 
624:                     print 'Invalid input' 
625:  
626:             radius_input = raw_input( 
627:                 'Radius of atomistic region (enter three values for an ellipsoid and one for a sphere):  ') 
628:             l_radii = radius_input.split() 
629:             try: 
630:                 if len(l_radii) == 1: 
631:                     log.append('Radius of atomistic region: ' + l_radii[0]) 
632:                     l_atom = mol.get_atoms_in_sphere(centre_xyz, float(l_radii[0])) 
633:                     l_res = mol.transform_atomlist_to_reslist(l_atom) 
634:                     unfrozen_res += l_res 
635:                     break 
636:                 elif len(l_radii) == 3: 
637:                     log.append('Radii for the atomistic region: ' + l_radii[0] + '  ' + l_radii[1] + '  ' + l_radii[2]) 
638:                     l_atom = mol.get_atoms_in_ellipsoid(centre_xyz, float(l_radii[0]), float(l_radii[1]), 
639:                                                         float(l_radii[2])) 
640:                     l_res = mol.transform_atomlist_to_reslist(l_atom) 
641:                     unfrozen_res += l_res 
642:                     break 
643:                 else: 
644:                     print 'Wrong input' 
645:             except ValueError: 
646:                 print 'Invalid input' 
647:         print 
648:         while 1:  #geometry is same as for atomistic region, if a change is wanted employ the above 
649:             check = raw_input('Add localised rigidification (y/n)?  ') 
650:             if check == 'yes' or check == 'y': 
651:                 log.append('Local rigidification added') 
652:                 radius_input = raw_input( 
653:                     'Radius of atomistic region (enter three values for an ellipsoid and one for a sphere):  ') 
654:                 l_radii_outer = radius_input.split() 
655:                 try: 
656:                     if len(l_radii_outer) == 1: 
657:                         log.append('Outer radius: ' + l_radii_outer[0]) 
658:                         l_atom_out = mol.get_atoms_in_sphere(centre_xyz, float(l_radii_outer[0])) 
659:                         l_res_out = mol.transform_atomlist_to_reslist(l_atom_out) 
660:                         for i in l_res_out: 
661:                             if i not in l_res: 
662:                                 local_res.append(i) 
663:                         break 
664:                     elif len(l_radii_outer) == 3: 
665:                         log.append( 
666:                             'Radii for locally rigid region: ' + l_radii_outer[0] + '  ' + l_radii_outer[1] + '  ' + 
667:                             l_radii_outer[2]) 
668:                         l_atom_out = mol.get_atoms_in_ellipsoid(centre_xyz, float(l_radii_outer[0]), 
669:                                                                 float(l_radii_outer[1]), 
670:                                                                 float(l_radii_outer[2])) 
671:                         l_res_out = mol.transform_atomlist_to_reslist(l_atom_out) 
672:                         for i in l_res_out: 
673:                             if i not in l_res: 
674:                                 local_res.append(i) 
675:                         break 
676:                     else: 
677:                         print 'Wrong input' 
678:                 except ValueError: 
679:                     print 'Invalid input' 
680:             else: 
681:                 break 
682:  
683:     #rigidification using a range of residues, sanity check in function that chooses 
684:     elif method == 2: 
685:         log.append('Method:  range selection') 
686:         check = raw_input('Define atomistic residues (y/n)?  ') 
687:         if check == 'yes' or check == 'y': 
688:             res = raw_input( 
689:                 'Enter residues separated by spaces, use r followed by two residues to define a range: ')  #input allows to enter 'r num1 num2' being the range between num1 and num2 
690:             res = res.split() 
691:             i = 0 
692:             log_res = '' 
693:             while i < len(res): 
694:                 if res[i] == 'r':  #if a range is indicated 
695:                     i += 1 
696:                     start_range = int(res[i])  #finds starting residue 
697:                     i += 1 
698:                     end_range = int(res[i])  #finds final residue 
699:                     for j in range(start_range, end_range + 1): 
700:                         unfrozen_res.append(j) 
701:                         log_res += str(j) + ' ' 
702:                     i += 1 
703:                 else: 
704:                     try: 
705:                         unfrozen_res.append(int(res[i])) 
706:                         log_res += res[i] + ' ' 
707:                         i += 1 
708:                     except ValueError: 
709:                         i += 1 
710:             log.append('Residues entered as atomistically: ' + log_res) 
711:         check = raw_input('Locally rigidify (y/n)?  ') 
712:         if check == 'yes' or check == 'y': 
713:             res = raw_input( 
714:                 'Enter residues separated by spaces, use r followed by two residues to define a range: ')  #input allows to enter 'r num1 num2' being the range between num1 and num2 
715:             res = res.split() 
716:             i = 0 
717:             log_res = '' 
718:             while i < len(res): 
719:                 if res[i] == 'r':  #if a range is indicated 
720:                     i += 1 
721:                     start_range = int(res[i])  #finds starting residue 
722:                     i += 1 
723:                     end_range = int(res[i])  #finds final residue 
724:                     i += 1 
725:                     for j in range(start_range, end_range + 1): 
726:                         if j <= mol.num_res(): 
727:                             local_res.append(j) 
728:                             log_res += str(j) + ' ' 
729:                         else: 
730:                             break 
731:                 else: 
732:                     try: 
733:                         if int(res[i]) <= mol.num_res(): 
734:                             local_res.append(int(res[i])) 
735:                             log_res += res[i] + ' ' 
736:                             i += 1 
737:                     except ValueError: 
738:                         i += 1 
739:             log.append('Residues entered for local rigidification: ' + log_res) 
740:  
741:     else: 
742:         return 'Wrong input' 
743:     return (unfrozen_res, local_res, log) 
744:  
745:  
746: def remove_res(mol): 
747:     unfrozen_res = [] 
748:     local_res = [] 
749:     log = [] 
750:     log.append('Removal of selected residues') 
751:     check = raw_input('Remove from the atomistic residues (y/n)?  ') 
752:     if check == 'yes' or check == 'y': 
753:         res = raw_input( 
754:             'Enter residues separated by spaces, use r followed by two residues to define a range: ')  #input allows to enter 'r num1 num2' being the range between num1 and num2 
755:         res = res.split() 
756:         i = 0 
757:         log_res = '' 
758:         while i < len(res): 
759:             if res[i] == 'r':  #if a range is indicated 
760:                 i += 1 
761:                 start_range = int(res[i])  #finds starting residue 
762:                 i += 1 
763:                 end_range = int(res[i])  #finds final residue 
764:                 for j in range(start_range, end_range + 1): 
765:                     unfrozen_res.append(j) 
766:                     log_res += str(j) + ' ' 
767:                 i += 1 
768:             else: 
769:                 try: 
770:                     unfrozen_res.append(int(res[i])) 
771:                     log_res += res[i] + ' ' 
772:                     i += 1 
773:                 except ValueError: 
774:                     i += 1 
775:         log.append('Residues entered to be removed: ' + log_res + "\n") 
776:     check = raw_input('Remove from the local selection (y/n)?  ') 
777:     if check == 'yes' or check == 'y': 
778:         res = raw_input( 
779:             'Enter residues separated by spaces, use r followed by two residues to define a range: ')  #input allows to enter 'r num1 num2' being the range between num1 and num2 
780:         res = res.split() 
781:         i = 0 
782:         log_res = '' 
783:         while i < len(res): 
784:             if res[i] == 'r':  #if a range is indicated 
785:                 i += 1 
786:                 start_range = int(res[i])  #finds starting residue 
787:                 i += 1 
788:                 end_range = int(res[i])  #finds final residue 
789:                 i += 1 
790:                 for j in range(start_range, end_range + 1): 
791:                     if j <= mol.num_res(): 
792:                         local_res.append(j) 
793:                         log_res += str(j) + ' ' 
794:                     else: 
795:                         break 
796:             else: 
797:                 try: 
798:                     if int(res[i]) <= mol.num_res(): 
799:                         local_res.append(int(res[i])) 
800:                         log_res += res[i] + ' ' 
801:                         i += 1 
802:                 except ValueError: 
803:                     i += 1 
804:         log.append('Residues entered to be removed from local rigidification: ' + log_res + "\n") 
805:     return (unfrozen_res, local_res, log) 
806:  
807:  
808: ###VARIABLE DEFINITION ################################################################################################ 
809: atom_dic = {}  #dictionary containing all information sorted by atom number: 
810:                #the information is in a tuple: (atomname, residue name, residue number, x, y, z) 
811: res_dic = {}  #dictionary containing all atoms per residue sorted by residue number 
812: res_atomistic = []  #residue list that are treated fully atomistic 
813: res_local = []  #residue list that are treated locally rigid 
814: res_frozen = []  #residue list that are treated fully rigid 
815: atoms_used = []  #atoms used in rigid bodies, prevents double usage! 
816: groups = {}  #dictionary of rigid groups 
817: groups_local = {} 
818: num_rbody = 0  #number of rigid bodies defined 
819:  
820: ###PARSING OF INPUT FILES ############################################################################################# 
821: while 1:  # starts input selection 
822:     atom_dic = reading_pdb_inpcrd(pdb_inp, coords_inp)  # opens files and loads them into a dictionary 
823:     res_dic = create_res_dic(atom_dic) 
824:     if atom_dic != {}:  # if the dictionary contains data the data is read into the molecule class 
825:         loaded_mol = protein(atom_dic, res_dic) 
826:         print 
827:         print 'Molecule was loaded' 
828:         print 'Number of residues: ', loaded_mol.num_res() 
829:         print 'Number of atoms: ', loaded_mol.num_atoms() 
830:         print 
831:         break 
832:     else: 
833:         print 'Files could not be found or they are empty - please try again' 
834:         print 
835:         pdb_inp = raw_input('PDB file: ') 
836:         coords_inp = raw_input('Coords input: ') 
837:  
838: check_file = open('rigid.log', 'w') 
839: check_file.write('Input file for structure: ' + pdb_inp + "\n") 
840: check_file.write('Coordinates file used: ' + coords_inp + "\n" + "\n") 
841: check_file.write('Number of residues:  ' + str(loaded_mol.num_res()) + "\n") 
842: check_file.write('Number of atoms:  ' + str(loaded_mol.num_atoms()) + "\n" + "\n") 
843: check_file.write('Time:  ' + time.strftime("%a, %d %b %Y %H:%M:%S") + "\n" + "\n") 
844:  
845: ###CHECK FOR LIGANDS AND UNKNOWN RESIDUES ############################################################################# 
846: check_file.write('Ligands and unknown residues:' + "\n") 
847: ligand_dic = loaded_mol.get_ligands_unknown_res() 
848: if ligand_dic == {}: 
849:     print 'There are no ligands or unknown residues in the molecule' 
850:     check_file.write('None' + "\n" + "\n") 
851: else: 
852:     print 'There are ligands or unknown residues.' 
853:     ligand_list = ligand_dic.keys() 
854:     for res_id in ligand_dic.keys(): 
855:         print str(res_id) + ' - ' + loaded_mol.res_name_from_id(res_id) 
856:         check_file.write(str(res_id) + ' - ' + loaded_mol.res_name_from_id(res_id) + "\n") 
857:     check_ligand = raw_input('Shall some or all of the above be treated as normal residues, all other residues will' 
858:                              ' be treated as fully atomistic? (y/n) ') 
859:     if check_ligand == ('y' or 'yes'): 
860:         while 1: 
861:             ligand_inp = raw_input('Enter the residue numbers separated by spaces (enter a for all residues): ') 
862:             if ligand_inp == 'a': 
863:                 print 'All residues/ligands are treated normally.' 
864:                 check_file.write('All are treated normally.' + "\n" + "\n") 
865:                 break 
866:             else: 
867:                 try: 
868:                     res_selected = [int(s) for s in ligand_inp.split()] 
869:                     new_ligand_list = [] 
870:                     string_ligands = '' 
871:                     for i in ligand_list: 
872:                         if i not in res_selected: 
873:                             new_ligand_list.append(i) 
874:                             string_ligands += str(i) + ' ' 
875:                     res_atomistic += new_ligand_list 
876:                     print 'The following residues will be treated atomistically: ' + string_ligands 
877:                     check_file.write('The following residues will be treated atomistically: ' + string_ligands 
878:                                      + "\n" + "\n") 
879:                     break 
880:                 except ValueError: 
881:                     print 'Please enter integer numbers.' 
882:     else: 
883:         res_atomistic += ligand_list 
884:         print 'The above listed residues/ligands will be treated fully atomistically.' 
885:         check_file.write('All are treated fully atomistically.' + "\n" + "\n") 
886:  
887: ###SELECTION OF RESIDUES FOR FULLY AND LOCALLY RIGID REGIONS ########################################################## 
888: print 
889: print 
890: print 'Selection of atomistic and locally rigidified residues' 
891: print 
892: check_file.write('Selection of atomistic and locally rigidified residues' + "\n") 
893: res_frozen = remove_selection_from_list(res_dic.keys(), res_atomistic) 
894: x = choose_res(loaded_mol)  #initiates the first round of selections 
895: y = merge_new_and_old_res_lists(res_frozen, res_local, res_atomistic, x[1], x[0]) 
896: res_frozen = y[2] 
897: res_local = y[1] 
898: res_atomistic = y[0] 
899: log_string = '' 
900: for strings in x[2]: 
901:     log_string += strings + "\n" 
902: check_file.write(log_string + "\n" + "\n") 
903:  
904: while 1:  #allows to choose more residues with different methods 
905:     print 
906:     print 'The following residues are not frozen (atomistic):' 
907:     print res_atomistic 
908:     print 
909:     print 'The following residues are selected to be locally rigidified:' 
910:     print res_local 
911:     print 
912:     print 'The following residues are selected to be fully rigidified:' 
913:     print res_frozen 
914:     print 
915:     check = raw_input('Enter more residues (1) or remove residues from previous selection (2)? Any other key to exit. ') 
916:     try: 
917:         value = int( 
918:             check)  #to enter more residues the choose_res function is called again with the same options as before 
919:         if int(check) == 1: 
920:             x = choose_res(loaded_mol) 
921:             y = merge_new_and_old_res_lists(res_frozen, res_local, res_atomistic, x[1], x[0]) 
922:             res_frozen = y[2] 
923:             res_local = y[1] 
924:             res_atomistic = y[0] 
925:             log_string = '' 
926:             for strings in x[2]: 
927:                 log_string += strings + "\n" 
928:             check_file.write(log_string + "\n" + "\n") 
929:         elif int(check) == 2: 
930:             x = remove_res(loaded_mol) 
931:             y = removal_selection(res_frozen, res_local, res_atomistic, x[1], 
932:                                   x[0])  #the user enters a selection for removal 
933:             res_frozen = y[2] 
934:             res_local = y[1] 
935:             res_atomistic = y[0] 
936:             check_file.write('Remove residues' + "\n" + 'Atomistic residues removed:' + "\n") 
937:             string = '' 
938:             for i in x[0]: 
939:                 string = string + str(i) + '  ' 
940:             check_file.write(string + "\n" + "\n") 
941:             check_file.write('Locally rigidified residues removed:' + "\n") 
942:             string = '' 
943:             for i in x[1]: 
944:                 string = string + str(i) + '  ' 
945:             check_file.write(string + "\n" + "\n") 
946:         else: 
947:             break 
948:  
949:     except ValueError: 
950:         break 
951: check_file.write('Selection completed' +  "\n" + "\n") 
952: string_atomistic = '' 
953: for res in res_atomistic: 
954:     string_atomistic += str(res) + ', ' 
955: string_local = '' 
956: for res in res_local: 
957:     string_local += str(res) + ', ' 
958: check_file.write('Atomistic residues: ' + string_atomistic + "\n" ) 
959: check_file.write('Locally rigidified residues: ' + string_local + "\n" ) 
960:  
961: ###GROUPING OF THE RIGID BODIES ####################################################################################### 
962: if res_frozen != []: 
963:     for i in res_frozen: 
964:         atoms_used += res_dic[i] 
965:     print 
966:     print 'Grouping of the rigid body' 
967:     print 
968:     print '1-Group all parts of the rigid system as one body' 
969:     print '2-Define groups' 
970:     print 
971:     while 1: 
972:         check = raw_input('Grouping method:  ') 
973:         try: 
974:             value = int(check) 
975:             if int(check) == 1:  #all residues are grouped as one rigid body 
976:                 num_rbody += 1 
977:                 check_file.write('All frozen residues are rigidified as one body' + "\n" + "\n" + "\n") 
978:                 groups = {1: res_frozen} 
979:                 break 
980:             elif int(check) == 2:  #residues are frozen in a number of rigid bodies 
981:                 check_file.write('The frozen residues are rigidified as more than one body' + "\n" + "\n" + "\n") 
982:                 print 'The following residues are frozen:' 
983:                 print res_frozen 
984:                 print 
985:                 res_frozen_left = res_frozen 
986:                 num_groups = int(raw_input('How many groups will be defined?  '))  #number of rigid bodies is defined 
987:                 counter = 0 
988:                 while counter < num_groups: 
989:                     counter += 1 
990:                     num_rbody += 1 
991:                     group = raw_input( 
992:                         'Residues for rigid body: ')  #input allows to enter 'r num1 num2' being the range between num1 and num2 
993:                     group = group.split() 
994:                     i = 0 
995:                     group_res = [] 
996:                     while i < len(group): 
997:                         if group[i] == 'r':  #if a range is indicated 
998:                             i += 1 
999:                             start_range = int(group[i])  #finds starting residue 
1000:                             i += 1 
1001:                             end_range = int(group[i])  #finds final residue 
1002:                             for j in range(start_range, end_range): 
1003:                                 if j in res_frozen_left:  #appends list if residues have not be assigned earlier 
1004:                                     group_res.append(j) 
1005:                         else: 
1006:                             try: 
1007:                                 if int(group[i]) in res_frozen_left: 
1008:                                     group_res.append(int(group[i]))  #appends single residues entered 
1009:                                 i += 1 
1010:                             except ValueError: 
1011:                                 i += 1 
1012:                     res_frozen_left = remove_selection_from_list(res_frozen_left, 
1013:                                                                  group_res)  #removes newly assigned residues from rest of residues 
1014:                     groups[num_rbody] = group_res  #defines a new entry in the dictionary that saves the groups 
1015:                     if res_frozen_left == []:  #automatically terminates when no residues are left to be assigned even if the number of rigid bodies entered is not reached 
1016:                         print 'All residues are assigned.' 
1017:                         break 
1018:                     if counter == num_groups and res_frozen_left != []:  #if all rigid bodies are assigned and there are still residues not assigned it deletes all groups and restarts the process 
1019:                         print 'Not all residues were assigned. Please start again.' 
1020:                         counter = 0 
1021:                         groups = {} 
1022:                         res_frozen_left = res_frozen 
1023:                 break 
1024:             else: 
1025:                 pass 
1026:         except ValueError: 
1027:             print 'Wrong input' 
1028:     print 
1029:     print 'The following rigid bodies have been specified:'  #prints all the entered rigid bodies from above 
1030:     for i in groups.keys(): 
1031:         print 'Rigid body: ', i 
1032:         print groups[i] 
1033:         print 
1034:  
1035: else: 
1036:     groups = {} 
1037:  
1038: ### LOCALLY RIGIDIFICATION, STANDARD GROUPS ########################################################################### 
1039: if res_local != []: 
1040:     print 'Grouping of the local rigid bodies' 
1041:     print 
1042:     local_scheme = [] 
1043:  
1044:     proline_check = raw_input('Group proline as rigid body C-O-N-CD-CG-CB-CA(1) or as  N-CD-CG-CB-CA(2)? ') 
1045:     if proline_check == '1': 
1046:         check_file.write('Grouped proline as rigid body C-O-N-CD-CG-CB-CA' + "\n") 
1047:         local_scheme.append(1) 
1048:     elif proline_check == '2': 
1049:         check_file.write('Grouped proline as rigid body  N-CD-CG-CB-CA' + "\n") 
1050:         local_scheme.append(2) 
1051:     else: 
1052:         local_scheme.append(0) 
1053:  
1054:     arginine_check = raw_input('Group arginine as rigid body NE-HE-CZ-NH1-HH11-HH12-NH2-HH21-HH22 (y/n)? ') 
1055:     if arginine_check == 'y': 
1056:         check_file.write('Grouped arginine as rigid body NE-HE-CZ-NH1-HH11-HH12-NH2-HH21-HH22' + "\n") 
1057:         local_scheme.append(1) 
1058:     else: 
1059:         local_scheme.append(0) 
1060:  
1061:     histidine_check = raw_input('Group histidine (HIS, HIE, HID, HIP) as rigid body CG-ND1-CE1-NE2-CD2 (y/n)? ') 
1062:     if histidine_check == 'y': 
1063:         check_file.write('Grouped histidine (HIS, HIE, HID, HIP) as rigid body CG-ND1-CE1-NE2-CD2' + "\n") 
1064:         local_scheme.append(1) 
1065:         local_scheme.append(1) 
1066:         local_scheme.append(1) 
1067:         local_scheme.append(1) 
1068:     else: 
1069:         local_scheme.append(0) 
1070:         local_scheme.append(0) 
1071:         local_scheme.append(0) 
1072:         local_scheme.append(0) 
1073:  
1074:     lysine_check = raw_input('Group lysine as rigid body NZ-HZ1-HZ2-HZ3 (y/n)? ') 
1075:     if lysine_check == 'y': 
1076:         check_file.write('Grouped lysine as rigid body NZ-HZ1-HZ2-HZ3' + "\n") 
1077:         local_scheme.append(1) 
1078:     else: 
1079:         local_scheme.append(0) 
1080:  
1081:     aspartic_check = raw_input('Group aspartic acid as rigid body CB-CG-OD1-OD2 (y/n)? ') 
1082:     if aspartic_check == 'y': 
1083:         check_file.write('Grouped aspartic acid as rigid body CB-CG-OD1-OD2' + "\n") 
1084:         local_scheme.append(1) 
1085:     else: 
1086:         local_scheme.append(0) 
1087:  
1088:     asparagine_check = raw_input('Group asparagine as rigid body CB-CG-OD1-ND2-HD21-HD22 (y/n)? ') 
1089:     if asparagine_check == 'y': 
1090:         check_file.write('Grouped asparagine as rigid body GB-CG-OD1-ND2-HD21-HD22' + "\n") 
1091:         local_scheme.append(1) 
1092:     else: 
1093:         local_scheme.append(0) 
1094:  
1095:     glutamic_check = raw_input('Group glutamic acid as rigid body CG-CD-OE1-OE2 (y/n)? ') 
1096:     if glutamic_check == 'y': 
1097:         check_file.write('Grouped glutamic acid as rigid body CG-CD-OE1-OE2' + "\n") 
1098:         local_scheme.append(1) 
1099:     else: 
1100:         local_scheme.append(0) 
1101:  
1102:     glutamine_check = raw_input('Group glutamine as rigid body CG-CD-OE1-NE2-HE21-HE22 (y/n)? ') 
1103:     if glutamine_check == 'y': 
1104:         check_file.write('Grouped glutamine as rigid body CG-CD-OE1-NE2-HE21-HE22' + "\n") 
1105:         local_scheme.append(1) 
1106:     else: 
1107:         local_scheme.append(0) 
1108:  
1109:     phenylalanine_check = raw_input('Group phenylalanine as rigid body CG-CD1-HD1-CE1-HE1-CZ-HZ-CE2-HE2-CD2-HD2 (y/n)? ') 
1110:     if phenylalanine_check == 'y': 
1111:         check_file.write('Grouped phenylalanine as rigid body CG-CD1-HD1-CE1-HE1-CZ-HZ-CE2-HE2-CD2-HD2' + "\n") 
1112:         local_scheme.append(1) 
1113:     else: 
1114:         local_scheme.append(0) 
1115:  
1116:     tyrosine_check = raw_input('Group tyrosine as rigid body CG-CD1-HD1-CE1-HE1-CZ-OH-CE2-HE2-CD2-HD2 (y/n)? ') 
1117:     if tyrosine_check == 'y': 
1118:         check_file.write('Grouped tyrosine as rigid body CG-CD1-HD1-CE1-HE1-CZ-OH-CE2-HE2-CD2-HD2' + "\n") 
1119:         local_scheme.append(1) 
1120:     else: 
1121:         local_scheme.append(0) 
1122:  
1123:     tryptophan_check = raw_input( 
1124:         'Group tryptophan as either: 1 rigid body:(CG-CD1-HD1-NE1-HE1-CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-CE3-HE3-CD2) ' 
1125:         'or 2 rigid bodies:(CG-CD1-HD1-NE1-HE1)-(CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-CE3-HE3-CD2)? ') 
1126:     if tryptophan_check == '1': 
1127:         check_file.write( 
1128:             'Grouped tryptophan as rigid body CG-CD1-HD1-NE1-HE1-CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-HZ3-CE3-HE3-CD2' + "\n" + "\n") 
1129:         local_scheme.append(1) 
1130:     elif tryptophan_check == '2': 
1131:         check_file.write( 
1132:             'Grouped tryptophan as two rigid bodies CG-CD1-HD1-NE1-HE1 and CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-HZ3-CE3-HE3-CD2' + "\n" + "\n") 
1133:         local_scheme.append(2) 
1134:     else: 
1135:         local_scheme.append(0) 
1136:  
1137:     loc = loaded_mol.get_rigid_for_local(res_local, local_scheme, atoms_used) 
1138:     groups_local = loc[0] 
1139:     atoms_used = loc[1] 
1140:     print 
1141:  
1142:  
1143:     check_peptide_bonds = raw_input('Shall the peptide bonds be treated as rigid bodies (excl. PRO) (y/n)? ') 
1144:     if check_peptide_bonds == 'y' or 'yes': 
1145:         pairlist = []  #all peptide bonds are part of two residues 
1146:         for i in res_local: 
1147:             if (i - 1, i) not in pairlist and (i - 1 > 0): 
1148:                 if loaded_mol.get_residue_name(i) != 'PRO': 
1149:                     pairlist.append((i - 1, i)) 
1150:             if (i, i + 1) not in pairlist and (i + 1 <= loaded_mol.num_res()): 
1151:                 if loaded_mol.get_residue_name(i + 1) != 'PRO': 
1152:                     pairlist.append((i, i + 1)) 
1153:         for j in pairlist: 
1154:             atomlist = [] 
1155:             atomlist.append(loaded_mol.get_atom_id_from_name('C', j[0])) 
1156:             atomlist.append(loaded_mol.get_atom_id_from_name('O', j[0])) 
1157:             atomlist.append(loaded_mol.get_atom_id_from_name('N', j[1])) 
1158:             atomlist.append(loaded_mol.get_atom_id_from_name('H', j[1])) 
1159:             check_double = 0 
1160:             for k in atomlist: 
1161:                 if k in atoms_used: 
1162:                     check_double = 1 
1163:             if check_double == 0: 
1164:                 atoms_used += atomlist 
1165:                 check_file.write( 
1166:                     'Peptide bond rigidified for residues ' + str(j[0]) + ' and ' + str(j[1]) + "\n") 
1167:                 groups_local[len(groups_local) + 1] = [j, (loaded_mol.get_residue_name(j[0]), 
1168:                                                            loaded_mol.get_residue_name(j[1]))] + atomlist 
1169:     check_file.write("\n") 
1170:     check_new_groups = raw_input('Enter additional user-defined rigid bodies (y/n)? ') 
1171:     if check_new_groups in ['yes', 'y']: 
1172:         print '1 - select a residues by name' 
1173:         print '2 - select one or more residues by id' 
1174:         while 1: 
1175:             while 1: 
1176:                 local_method = raw_input('Method (1/2): ') 
1177:                 try: 
1178:                     if int(local_method) in [1, 2]: 
1179:                         break 
1180:                     else: 
1181:                         print 'Invalid choice' 
1182:                 except ValueError: 
1183:                     print 'Invalid input' 
1184:  
1185:             if int(local_method) == 1: 
1186:                 while 1: 
1187:                     exit_1 = 0 
1188:                     res_name = raw_input('Residue name (three letter code, n to exit): ') 
1189:                     if res_name == 'n': 
1190:                         exit_1 = 1 
1191:                         break 
1192:                     try: 
1193:                         res_list = loaded_mol.get_all_res_name(res_local, res_name) 
1194:                         if res_list != []: 
1195:                             check_file.write('Additional rigid bodies for residues: ' + res_name +  "\n") 
1196:                             break 
1197:                         else: 
1198:                             print 'None of the residues for local rigidification is of that type.' 
1199:                     except ValueError: 
1200:                         print 'Invalid input' 
1201:                 if exit_1 == 0: 
1202:                     all_atoms_res_name = [] 
1203:                     atoms_used_name = [] 
1204:                     for i in res_list: 
1205:                         all_atoms_res_name += res_dic[i] 
1206:                     for j in all_atoms_res_name: 
1207:                         if j in atoms_used: 
1208:                             atoms_used_name.append(j) 
1209:                     if atoms_used_name != []: 
1210:                         print 'The following atoms are already in local rigid bodies: ' 
1211:                         for k in atoms_used_name: 
1212:                             print str(k) + ' ' + atom_dic[k][0] + ' ' + atom_dic[k][1] + ' ' + str(atom_dic[k][2]) 
1213:                     print 
1214:                     print 'The residues %s contain the following atoms:' % res_name 
1215:                     atoms_in_res = loaded_mol.get_atom_list(res_list[0]) 
1216:                     for l in atoms_in_res: 
1217:                         print atom_dic[l][0] 
1218:                     print 
1219:                     atom_seq = raw_input('Enter a sequence of atom names separated by spaces: ') 
1220:                     for m in res_list: 
1221:                         check_res = 0 
1222:                         atomlist = res_dic[m] 
1223:                         atoms_rbody = [] 
1224:                         for n in atomlist: 
1225:                             if atom_dic[n][0] in atom_seq.split(): 
1226:                                 atoms_rbody.append(n) 
1227:                         if len(atoms_rbody) != len(atom_seq.split()): 
1228:                             print 'Not all atoms entered are in the residue.' 
1229:                             check_file.write('For residue: ' + str(m) + ' - not all entered atoms in residue' + "\n") 
1230:                             check_res = 1 
1231:                         if len(atoms_rbody) < 3: 
1232:                             print 'Rigid bodies need to have at least 3 atoms.' 
1233:                             check_file.write('For residue: ' + str(m) + ' - not 3 or more atoms in residue' + "\n") 
1234:                             check_res = 1 
1235:                         for o in atoms_rbody: 
1236:                             if o in atoms_used_name: 
1237:                                 print 'Atoms cannot belong to two rigid bodies.' 
1238:                                 check_file.write('For residue: ' + str(m) + ' - atoms cannot belong to two rigid bodies' 
1239:                                                  + "\n" + "\n") 
1240:                                 check_res = 1 
1241:                                 break 
1242:                         if loaded_mol.check_linear(atoms_rbody, tolerance): 
1243:                             print 'Linear arrangement of atoms cannot be a rigid body.' 
1244:                             check_file.write('For residue: ' + str(m) + ' - linear group entered' + "\n") 
1245:                             check_res = 1 
1246:                         if check_res == 0: 
1247:                             groups_local[len(groups_local) + 1] = [m, loaded_mol.get_residue_name(m)] + atoms_rbody 
1248:                             atoms_used += atoms_rbody 
1249:                             check_file.write('For residue: ' + str(m) + ' new group: ' + atom_seq + "\n") 
1250:                     check_file.write("\n") 
1251:             elif int(local_method) == 2: 
1252:                 while 1: 
1253:                     exit_2 = 0 
1254:                     res_inp = raw_input('Residue numbers (n to exit): ') 
1255:                     if res_inp == 'n': 
1256:                         exit_2 = 1 
1257:                         break 
1258:                     try: 
1259:                         res_list_str = res_inp.split() 
1260:                         res_list_int = [] 
1261:                         for res in res_list_str: 
1262:                             res_list_int.append(int(res)) 
1263:                         if res_list_int != []: 
1264:                             check_file.write('Additional rigid bodies for residues: ' + res_inp  + "\n") 
1265:                             break 
1266:                         else: 
1267:                             print 'No residues entered' 
1268:                     except ValueError: 
1269:                         print 'Invalid input' 
1270:                 if exit_2 == 0: 
1271:                     all_atoms_res = [] 
1272:                     atoms_used_res = [] 
1273:                     for i in res_list_int: 
1274:                         all_atoms_res += res_dic[i] 
1275:                     for j in all_atoms_res: 
1276:                         if j in atoms_used: 
1277:                             atoms_used_res.append(j) 
1278:                     atoms_allowed = [] 
1279:                     for atom in all_atoms_res: 
1280:                         if atom not in atoms_used_res: 
1281:                             atoms_allowed.append(atom) 
1282:                     if atoms_used_res != []: 
1283:                         print 'The following atoms are already in local rigid bodies: ' 
1284:                         for k in atoms_used_res: 
1285:                             print str(k) + ' ' + atom_dic[k][0] + ' ' + atom_dic[k][1] + ' ' + str(atom_dic[k][2]) 
1286:                     print 
1287:                     print 'The residues %s contain the following atoms that are not assigned to rigid bodies:' % res_list_str 
1288:                     string_atoms = '' 
1289:                     for atom in atoms_allowed: 
1290:                         string_atoms += str(atom) + ' - ' + atom_dic[atom][0] + ', ' 
1291:                     print string_atoms 
1292:                     print 
1293:                     atom_seq = raw_input('Enter a sequence of atom numbers separated by spaces: ') 
1294:                     check_res = 0 
1295:                     atoms_rbody = [] 
1296:                     check_rbody = 0 
1297:                     for atom_str in atom_seq.split(): 
1298:                         try: 
1299:                             if int(atom_str) in atoms_allowed: 
1300:                                 atoms_rbody.append(int(atom_str)) 
1301:                             elif int(atom_str) in atoms_used_res: 
1302:                                 print 'Atoms cannot belong to two rigid bodies.' 
1303:                                 check_file.write('Atoms cannot belong to two rigid bodies'  + "\n") 
1304:                                 check_rbody = 1 
1305:                                 break 
1306:                         except ValueError: 
1307:                             pass 
1308:                     print len(atoms_rbody) 
1309:                     print len(atom_seq.split()) 
1310:                     if len(atoms_rbody) != len(atom_seq.split()): 
1311:                         print 'Not all atoms entered are in the residue.' 
1312:                         check_file.write('Not all entered atoms in residue' + "\n") 
1313:                         check_rbody = 1 
1314:                     if len(atoms_rbody) < 3: 
1315:                         print 'Rigid bodies need to have at least 3 atoms.' 
1316:                         check_file.write('Not 3 or more atoms in residue' + "\n") 
1317:                         check_rbody = 1 
1318:                     if loaded_mol.check_linear(atoms_rbody, tolerance): 
1319:                         print 'Linear arrangement of atoms cannot be a rigid body.' 
1320:                         check_file.write('Linear group entered' + "\n") 
1321:                         check_rbody = 1 
1322:                     if check_rbody == 0: 
1323:                         groups_local[len(groups_local) + 1] = ['', ''] + atoms_rbody 
1324:                         atoms_used += atoms_rbody 
1325:                         check_file.write('For residues: ' + res_inp + ' new group: ' + atom_seq + "\n") 
1326:                     check_file.write("\n") 
1327:             check_more = raw_input('Enter more rigid bodies (y/n)? ') 
1328:             if check_more not in ['y', 'yes']: 
1329:                 break 
1330:  
1331: ### WRITING OUTPUT FILES ############################################################################################## 
1332:  
1333: print 'Output files will be written now.' 
1334:  
1335: ### COORDSINIRIGID ### 
1336: coords_out = open('coordsinirigid', 'w') 
1337: for i in range(1, len(atom_dic.keys()) + 1): 
1338:     atom = str(atom_dic[i][3]) + '   ' + str(atom_dic[i][4]) + '   ' + str(atom_dic[i][5]) 
1339:     coords_out.write(atom) 
1340:     if i < len(atom_dic.keys()): 
1341:         coords_out.write("\n") 
1342: coords_out.close() 
1343:  
1344: ### RBODYCONFIG FILE ### 
1345: groups_file = open('rbodyconfig', 'w') 
1346: for group in range(1, len(groups.keys()) + 1): 
1347:     string = '' 
1348:     counter = 0 
1349:     for res in groups[group]: 
1350:         for atom in res_dic[res]: 
1351:             string += str(atom) + "\n" 
1352:             counter += 1 
1353:     groups_file.write('GROUP ' + str(counter) + "\n" + string) 
1354: for group_local in range(1, len(groups_local.keys()) + 1): 
1355:     string = '' 
1356:     counter = 0 
1357:     for atom in groups_local[group_local][2:]: 
1358:         string += str(atom) + "\n" 
1359:         counter += 1 
1360:     groups_file.write('GROUP ' + str(counter) + "\n" + string) 
1361: groups_file.close() 
1362:  
1363: ### FINISH LOG FILE ### 
1364: check_file.close() 
1365:  
1366: print 'Selection process completed - change topology file to exclude interactions within the rigid bodies.' 
1367:  
1368: if pymol_check: 
1369:     mol_view = PymolView(pdb_inp) 
1370:     mol_view.apply_colour_scheme(res_atomistic, res_local, 2, 3) 
1371:     pymol.cmd.save('rigid_%s' % pdb_inp, format='png') 
1372:     time.sleep(10) 
1373:     mol_view.__del__() 
1374:   1: 
   2: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/AMBER/rigidbody/genrigid_input.py' in revision 29058


r29059/README_genrigid 2017-01-21 10:36:19.398250000 +0000 r29058/README_genrigid 2017-01-21 10:36:19.910250000 +0000
  1: #Script by Konstantin Roeder kr366  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/AMBER/rigidbody/README_genrigid' in revision 29058
  2:  
  3: HOW TO USE IT 
  4:  
  5: 1. Either start genrigid_input.py or start with options genrigid_input.py [pdb name] [coords.inpcrd] [tolerance] [pymol]. 
  6:         [pdb name]        required file name 
  7:         [coords.inpcrd]   required file name 
  8:         [tolerance]       tolerance for linearity check (angles need to be between 0 + tolerance and 180 - tolerance to be counted as non- 
  9:                           linear), default 0.01 
 10:         [pymol]           requires pymol installed, 'pymol' loads a graphical representation of the set rigidification 
 11:  
 12: 2. Check number of residues and number of atoms, if they are not what they should be, check the input files. 
 13:  
 14: 3. If the script detects unknown residues, they will be highlighted next. Only if yes or y is entered residues can be set as normal, all other  
 15:    keys add these residues to the atomistic list. 
 16:  
 17: 4. Select and remove from selection following the script. The residue selections can be either one by one with spaces or using 'r integer1 
 18:    integer2' to indicate the range from integer1 to integer2. If a residue is selected to be atomistic, selecting it again for the local rigid 
 19:    bodies will not change the assignement (this change can only be done by removing it and then add it again). 
 20:  
 21: 5. Grouping of the fully rigid region: enter residue by residue or ranges for option 2. 
 22:  
 23: 6. Local rigid bodies: peptide bonds will be rigidified for all residues in the local rigid selection. If the neighbouring residue belongs to 
 24:    a fully rigid body, this peptide bond will not be used. If it belongs to an atomistic residue it will be used. 
 25:  
 26: 7. Atoms entered for self-made rigid bodies are checked for linearity, double usage of atoms and existence. If one of the checks fails, the    
 27:    selection is ignored. 
 28:  
 29: HOW TO CHANGE IT 
 30:  
 31: The script relies on two dictionaries: res_dic and atom_dic. All information is stored in them, and they are read into the protein class. Furthermore the rigid bodies are stored in groups (for the fully rigid region, residue ids are stored) and groups_local (residue names, ids and atom ids are stored). The list atoms_used contains all atoms that are already assigned to rigid bodies and needs to be updated for every new rigid body defined. 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0