hdiff output

r31935/genrigid_input.py 2017-02-17 12:30:08.542891054 +0000 r31934/genrigid_input.py 2017-02-17 12:30:08.978896893 +0000
  1: # SCRIPT TO DEFINE THE RIGIDIFICATION FOR PROTEINS  1: # SCRIPT TO DEFINE THE RIGIDIFICATION FOR PROTEINS
  2: # by Konstantin Roeder  2: #by Konstantin Roeder
  3: # creates input files rbodyconfig and coordsinirigid and additionally a log file rigid.log  3: #creates input files rbodyconfig and coordsinirigid and additionally a log file rigid.log
  4:   4: 
  5: # ! /usr/bin/env python  5: #! /usr/bin/env python
  6:   6: 
  7: import sys  7: import sys
  8: import time  8: import time
  9: import math  9: import math
 10:  10: 
 11: if len(sys.argv) > 1: 11: if len(sys.argv) > 1:
 12:     print 'Script creates input files for RIGIDINIT keyword, needs inpcrd and pdb' + \ 12:     print 'Script creates input files for RIGIDINIT keyword, needs inpcrd and pdb'
 13:           'for RNA and DNA only post AMBER12 atom names' 
 14:     print 13:     print
 15:     pdb_inp = sys.argv[1] 14:     pdb_inp = sys.argv[1]
 16:     coords_inp = sys.argv[2] 15:     coords_inp = sys.argv[2]
 17:     # defines the tolerance for the checks whether rigid bodies are linear 16:     #defines the tolerance for the checks whether rigid bodies are linear
 18:     try: 17:     try:
 19:         tolerance = float(sys.argv[3]) 18:         tolerance = float(sys.argv[3])
 20:     except IndexError: 19:     except IndexError:
 21:         tolerance = 0.01 20:         tolerance = 0.01
 22:     # Pymol use; change to False if pymol is not installed/used 21:     #Pymol use; change to False if pymol is not installed/used
 23:     try: 22:     try:
 24:         if sys.argv[4] == 'pymol': 23:         if sys.argv[4] == 'pymol':
 25:             pymol_check = True 24:             pymol_check = True
 26:     except IndexError: 25:     except IndexError:
 27:         pymol_check = False 26:         pymol_check = False
 28:  27: 
 29: else: 28: else:
 30:     tolerance = 0.01 29:     tolerance = 0.01
 31:     pymol_check = False 30:     pymol_check = False
 32:     pdb_inp = raw_input('PDB file: ') 31:     pdb_inp = raw_input('PDB file: ')
 33:     coords_inp = raw_input('Coords input: ') 32:     coords_inp = raw_input('Coords input: ')
 34:  33: 
 35: inp_sys = raw_input('Use script for peptides (1), RNA (2), DNA (3), or mixed (4)? ') 
 36: try: 
 37:     if int(inp_sys) == 1: 
 38:         protein_t = True 
 39:         rna_t = False 
 40:         dna_t = False 
 41:     elif int(inp_sys) == 2: 
 42:         protein_t = False 
 43:         rna_t = True 
 44:         dna_t = False 
 45:     elif int(inp_sys) == 3: 
 46:         protein_t = False 
 47:         rna_t = False 
 48:         dna_t = True 
 49:     else: 
 50:         protein_t = True 
 51:         rna_t = True 
 52:         dna_t = True 
 53: except ValueError: 
 54:     print 'Assume RNA and protein in mixed system.' 
 55:     protein_t = True 
 56:     rna_t = True 
 57:     dna_t = True 
 58:  
 59: if pymol_check: 34: if pymol_check:
 60:     import __main__ 35:     import __main__
 61:  36: 
 62:     __main__.pymol_argv = ['pymol', '-qix']  # Pymol: quiet and no GUI(internal and external) 37:     __main__.pymol_argv = ['pymol', '-qix']  # Pymol: quiet and no GUI(internal and external)
 63:     import pymol 38:     import pymol
 64:  39: 
 65:  40: #class containing the main methods for the script, new functionality should be embedded here
 66: # class containing the main methods for the script, new functionality should be embedded here 
 67: class protein(): 41: class protein():
 68:     def __init__(self, atom_dic, res_dic): 42:     def __init__(self, atom_dic, res_dic):
 69:         self.atom_dic = atom_dic  # dictionary of all atoms with x,y,z coordinates and res id and name and atom name 43:         self.atom_dic = atom_dic  #dictionary of all atoms with x,y,z coordinates and res id and name and atom name
 70:         self.res_dic = res_dic  # dictionary containing a list of atom ids for all residues 44:         self.res_dic = res_dic    #dictionary containing a list of atom ids for all residues
 71:  45: 
 72:     def num_atoms(self):  # number of atoms 46:     def num_atoms(self): #number of atoms
 73:         return len(self.atom_dic) 47:         return len(self.atom_dic)
 74:  48: 
 75:     def num_res(self):  # number of residues 49:     def num_res(self): #number of residues
 76:         return len(self.res_dic) 50:         return len(self.res_dic)
 77:  51: 
 78:     # returns the name of a given residue 52:     #returns the name of a given residue
 79:     def get_residue_name(self, residue): 53:     def get_residue_name(self, residue):
 80:         return atom_dic[res_dic[residue][0]][1] 54:         return atom_dic[res_dic[residue][0]][1]
 81:  55: 
 82:     # prints the atoms in a specified residue(need number of residue) 56:     #prints the atoms in a specified residue(need number of residue)
 83:     # includes the atom number, atom name and coordinates 57:     #includes the atom number, atom name and coordinates
 84:     def get_atoms_in_res(self, residue): 58:     def get_atoms_in_res(self, residue):
 85:         all_atoms = self.atom_dic.keys() 59:         all_atoms = self.atom_dic.keys()
 86:         atoms_in_res = [] 60:         atoms_in_res = []
 87:         for i in range(1, len(self.atom_dic) + 1): 61:         for i in range(1, len(self.atom_dic) + 1):
 88:             atom = self.atom_dic[i] 62:             atom = self.atom_dic[i]
 89:             if residue == atom[2]: 63:             if residue == atom[2]:
 90:                 atoms_in_res.append(all_atoms[i - 1]) 64:                 atoms_in_res.append(all_atoms[i - 1])
 91:                 residuename = atom[1] 65:                 residuename = atom[1]
 92:         if atoms_in_res == []: 66:         if atoms_in_res == []:
 93:             return 'Residue not in molecule' 67:             return 'Residue not in molecule'
 94:         print 'Residue:', residue, residuename 68:         print 'Residue:', residue, residuename
 95:         print 69:         print
 96:         print 'Atoms:' 70:         print 'Atoms:'
 97:         for i in range(0, len(atoms_in_res)): 71:         for i in range(0, len(atoms_in_res)):
 98:             j = atoms_in_res[i] 72:             j = atoms_in_res[i]
 99:             atom = self.atom_dic[j] 73:             atom = self.atom_dic[j]
100:             print j, atom[0], atom[3], atom[4], atom[5] 74:             print j, atom[0], atom[3], atom[4], atom[5]
101:         return 75:         return
102:  76: 
103:     # returns a list of atoms in a residue 77:     #returns a list of atoms in a residue
104:     def get_atom_list(self, residue): 78:     def get_atom_list(self, residue):
105:         all_atoms = self.atom_dic.keys() 79:         all_atoms = self.atom_dic.keys()
106:         atomlist = [] 80:         atomlist = []
107:         for i in range(1, len(self.atom_dic) + 1): 81:         for i in range(1, len(self.atom_dic) + 1):
108:             atom = self.atom_dic[i] 82:             atom = self.atom_dic[i]
109:             if residue == atom[2]: 83:             if residue == atom[2]:
110:                 atomlist.append(all_atoms[i - 1]) 84:                 atomlist.append(all_atoms[i - 1])
111:         return atomlist 85:         return atomlist
112:  86: 
113:     # allows the input of a centre from the user, is complemented by the selection of the CoM as centre 87:     #allows the input of a centre from the user, is complemented by the selection of the CoM as centre
114:     def choose_centre(self): 88:     def choose_centre(self):
115:         check = raw_input('Display atoms in residue (y/n)?  ') 89:         check = raw_input('Display atoms in residue (y/n)?  ')
116:         if check == 'y': 90:         if check == 'y':
117:             while 1:  # acts as input mask --> once entered any number of residues can be viewed at 91:             while 1:  #acts as input mask --> once entered any number of residues can be viewed at
118:                 residue_number = raw_input('Residue number: (enter x to leave) ') 92:                 residue_number = raw_input('Residue number: (enter x to leave) ')
119:                 if residue_number == 'x': 93:                 if residue_number == 'x':
120:                     print 94:                     print
121:                     break 95:                     break
122:                 else: 96:                 else:
123:                     try: 97:                     try:
124:                         self.get_atoms_in_res(int(residue_number)) 98:                         self.get_atoms_in_res(int(residue_number))
125:                         print 99:                         print
126:                     except ValueError:100:                     except ValueError:
127:                         pass101:                         pass
128:         while 1:  # loop for choice of molecule, only accepts atoms in the molecule102:         while 1:  # loop for choice of molecule, only accepts atoms in the molecule
129:             centre_input = int(raw_input('Choose atom:  '))103:             centre_input = int(raw_input('Choose atom:  '))
130:             if centre_input > self.num_atoms() or centre_input < 1:104:             if centre_input > self.num_atoms() or centre_input < 1:
131:                 print 'Not in molecule.'105:                 print 'Not in molecule.'
132:             else:106:             else:
133:                 break107:                 break
134:         return centre_input108:         return centre_input
135: 109: 
136:     # transforms a given list of atoms into a list of residues110:     #transforms a given list of atoms into a list of residues
137:     # can be used to find the residues in a sphere from the output of the atoms in sphere method111:     #can be used to find the residues in a sphere from the output of the atoms in sphere method
138:     def transform_atomlist_to_reslist(self, atomlist):112:     def transform_atomlist_to_reslist(self, atomlist):
139:         reslist = []113:         reslist = []
140:         for i in atomlist:114:         for i in atomlist:
141:             atom_info = self.atom_dic[i]115:             atom_info = self.atom_dic[i]
142:             k = 0116:             k = 0
143:             if reslist != []:117:             if reslist != []:
144:                 for j in range(0, len(reslist)):118:                 for j in range(0, len(reslist)):
145:                     if reslist[j] == atom_info[2]:119:                     if reslist[j] == atom_info[2]:
146:                         k = 1120:                         k = 1
147:             if k == 0:121:             if k == 0:
149:             else:123:             else:
150:                 pass124:                 pass
151:         return reslist125:         return reslist
152: 126: 
153:     def transform_reslist_to_atomlist(self, reslist):127:     def transform_reslist_to_atomlist(self, reslist):
154:         atomlist = []128:         atomlist = []
155:         for i in reslist:129:         for i in reslist:
156:             atomlist += self.res_dic[i]130:             atomlist += self.res_dic[i]
157:         return atomlist131:         return atomlist
158: 132: 
159:     # get residue name from id133:     #get residue name from id
160:     def res_name_from_id(self, residue_id):134:     def res_name_from_id(self, residue_id):
161:         atom_id = res_dic[residue_id][0]135:         atom_id = res_dic[residue_id][0]
162:         return atom_dic[atom_id][1]136:         return atom_dic[atom_id][1]
163: 137: 
164:     # checks for unknown residues/ligands, add standard residues here if needed138:     #checks for unknown residues/ligands, add standard residues here if needed
165:     def get_ligands_unknown_res(self):139:     def get_ligands_unknown_res(self):
166:         list_res_names = ['ALA', 'ARG', 'ASN', 'ASP', 'ASX', 'CYS', 'CYX', 'GLN', 'GLU', 'GLY', 'GLX', 'HIS', 'HIE',140:         list_res_names = ['ALA', 'ARG', 'ASN', 'ASP', 'ASX', 'CYS', 'CYX', 'GLN', 'GLU', 'GLY', 'GLX', 'HIS', 'HIE',
167:                           'HID', 'HIP', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL',141:                           'HID', 'HIP', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL',
168:                           'NALA', 'NARG', 'NASN', 'NASP', 'NASX', 'NCYS', 'NCYX', 'NGLN', 'NGLU', 'NGLY', 'NGLX',142:                           'NALA', 'NARG', 'NASN', 'NASP', 'NASX', 'NCYS', 'NCYX', 'NGLN', 'NGLU', 'NGLY', 'NGLX',
169:                           'NHIS', 'NHIE', 'NHID', 'NHIP', 'NILE', 'NLEU', 'NLYS', 'NMET', 'NPHE', 'NPRO', 'NSER',143:                           'NHIS', 'NHIE', 'NHID', 'NHIP', 'NILE', 'NLEU', 'NLYS', 'NMET', 'NPHE', 'NPRO', 'NSER',
170:                           'NTHR', 'NTRP', 'NTYR', 'NVAL', 'CALA', 'CARG', 'CASN', 'CASP', 'CASX', 'CCYS', 'CCYX',144:                           'NTHR', 'NTRP', 'NTYR', 'NVAL', 'CALA', 'CARG', 'CASN', 'CASP', 'CASX', 'CCYS', 'CCYX',
171:                           'CGLN', 'CGLU', 'CGLY', 'CGLX', 'CHIS', 'CHIE', 'CHID', 'CHIP', 'CILE', 'CLEU', 'CLYS',145:                           'CGLN', 'CGLU', 'CGLY', 'CGLX', 'CHIS', 'CHIE', 'CHID', 'CHIP', 'CILE', 'CLEU', 'CLYS',
172:                           'CMET', 'CPHE', 'CPRO', 'CSER', 'CTHR', 'CTRP', 'CTYR', 'CVAL', 'RAN', 'RCN', 'RGN', 'RUN',146:                           'CMET', 'CPHE', 'CPRO', 'CSER', 'CTHR', 'CTRP', 'CTYR', 'CVAL']
173:                           'A', 'A3', 'A5', 'AN', 'C', 'C3', 'C5', 'CN', 'G', 'G3', 'G5', 'GN', 'U', 'U3', 'U5', 'UN', 
174:                           'OHE', 'RA', 'RA3', 'RA5', 'RC', 'RC3', 'RC5', 'RG', 'RG3', 'RG5', 'RU', 'RU3', 'RU5', 
175:                           'DA', 'DA5', 'DC', 'DC5', 'DG', 'DT', 'DT5', 'DA3', 'DC3', 'DG3', 'DT3', 
176:                           'DAN', 'DCN', 'DGN', 'DTN'] 
177:         ligand_dic = {}147:         ligand_dic = {}
178:         for i in self.res_dic.keys():148:         for i in self.res_dic.keys():
179:             if self.res_name_from_id(i) not in list_res_names:149:             if self.res_name_from_id(i) not in list_res_names:
180:                 ligand_dic[i] = self.res_dic[i]150:                 ligand_dic[i] = self.res_dic[i]
181:         return ligand_dic151:         return ligand_dic
182: 152: 
183:     # returns a list of atoms within a sphere of entered radius153:     #returns a list of atoms within a sphere of entered radius
184:     # needs an atom number as centre154:     #needs an atom number as centre
185:     def get_atoms_in_sphere(self, centre_atom, radius):155:     def get_atoms_in_sphere(self, centre_atom, radius):
186:         all_atoms = self.atom_dic.keys()156:         all_atoms = self.atom_dic.keys()
187:         atomlist = []157:         atomlist = []
188:         for i in range(1, len(self.atom_dic) + 1):158:         for i in range(1, len(self.atom_dic) + 1):
189:             atom = self.atom_dic[i]159:             atom = self.atom_dic[i]
190:             distance = (atom[3] - centre_atom[0]) ** 2 + (atom[4] - centre_atom[1]) ** 2 + (atom[5] - centre_atom[160:             distance = (atom[3] - centre_atom[0]) ** 2 + (atom[4] - centre_atom[1]) ** 2 + (atom[5] - centre_atom[
191:                 2]) ** 2161:                 2]) ** 2
192:             if distance <= radius:162:             if distance <= radius:
193:                 atomlist.append(all_atoms[i - 1])163:                 atomlist.append(all_atoms[i - 1])
194:             else:164:             else:
195:                 pass165:                 pass
196:         return atomlist166:         return atomlist
197: 167: 
198:     # returns the atoms within an ellipsoidal volume168:     #returns the atoms within an ellipsoidal volume
199:     def get_atoms_in_ellipsoid(self, centre_atom, a, b, c):169:     def get_atoms_in_ellipsoid(self, centre_atom, a, b, c):
200:         all_atoms = self.atom_dic.keys()170:         all_atoms = self.atom_dic.keys()
201:         atomlist = []171:         atomlist = []
202:         for i in range(1, len(self.atom_dic) + 1):172:         for i in range(1, len(self.atom_dic) + 1):
203:             atom = self.atom_dic[i]173:             atom = self.atom_dic[i]
204:             x1 = ((atom[3] - centre_atom[0]) / a) ** 2174:             x1 = ((atom[3] - centre_atom[0]) / a) ** 2
205:             x2 = ((atom[4] - centre_atom[1]) / b) ** 2175:             x2 = ((atom[4] - centre_atom[1]) / b) ** 2
206:             x3 = ((atom[5] - centre_atom[2]) / c) ** 2176:             x3 = ((atom[5] - centre_atom[2]) / c) ** 2
207:             distance = x1 + x2 + x3177:             distance = x1 + x2 + x3
208:             if distance <= 1:178:             if distance <= 1:
209:                 atomlist.append(all_atoms[i - 1])179:                 atomlist.append(all_atoms[i - 1])
210:             else:180:             else:
211:                 pass181:                 pass
212:         return atomlist182:         return atomlist
213: 183: 
214:     # returns a list of atoms in a residue184:     #returns a list of atoms in a residue
215:     def get_atom_list(self, residue):185:     def get_atom_list(self, residue):
216:         all_atoms = self.atom_dic.keys()186:         all_atoms = self.atom_dic.keys()
217:         atomlist = []187:         atomlist = []
218:         for i in range(1, len(self.atom_dic) + 1):188:         for i in range(1, len(self.atom_dic) + 1):
219:             atom = self.atom_dic[i]189:             atom = self.atom_dic[i]
220:             if residue == atom[2]:190:             if residue == atom[2]:
221:                 atomlist.append(all_atoms[i - 1])191:                 atomlist.append(all_atoms[i - 1])
222:         return atomlist192:         return atomlist
223: 193: 
224:     # 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 below194:     #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
225:     def res_mass(self, residue):195:     def res_mass(self, residue):
226:         atom_list = self.get_atom_list(residue)196:         atom_list = self.get_atom_list(residue)
227:         mass_dic = {'H': 1.007825, 'D': 2.014102, 'C': 12.0116, 'N': 14.00728, 'O': 15.99977, 'P': 30.973762,197:         mass_dic = {'H': 1.007825, 'D': 2.014102, 'C': 12.0116, 'N': 14.00728, 'O': 15.99977, 'S': 32.076}
228:                     'S': 32.076} 
229:         mass = 0198:         mass = 0
230:         for i in atom_list:199:         for i in atom_list:
231:             for j in mass_dic.keys():200:             for j in mass_dic.keys():
232:                 if j == self.atom_dic[i][0][0]:201:                 if j == self.atom_dic[i][0][0]:
233:                     mass += mass_dic[j]202:                     mass += mass_dic[j]
234:         return mass203:         return mass
235: 204: 
236:     # gives the coordinates for the mass_weighted centre of a residue205:     #gives the coordinates for the mass_weighted centre of a residue
237:     def mass_weighted_centre(self, residue):206:     def mass_weighted_centre(self, residue):
238:         atom_list = self.get_atom_list(residue)207:         atom_list = self.get_atom_list(residue)
239:         mass_dic = {'H': 1.007825, 'D': 2.014102, 'C': 12.0116, 'N': 14.00728, 'O': 15.99977, 'P': 30.973762,208:         mass_dic = {'H': 1.007825, 'D': 2.014102, 'C': 12.0116, 'N': 14.00728, 'O': 15.99977, 'S': 32.076}
240:                     'S': 32.076} 
241:         X = 0209:         X = 0
242:         Y = 0210:         Y = 0
243:         Z = 0211:         Z = 0
244:         mass_res = float(self.res_mass(residue))212:         mass_res = float(self.res_mass(residue))
245:         for i in atom_list:213:         for i in atom_list:
246:             for j in mass_dic.keys():214:             for j in mass_dic.keys():
247:                 if j == self.atom_dic[i][0][0]:215:                 if j == self.atom_dic[i][0][0]:
248:                     X += mass_dic[j] * self.atom_dic[i][3]216:                     X += mass_dic[j] * self.atom_dic[i][3]
249:         for i in atom_list:217:         for i in atom_list:
250:             for j in mass_dic.keys():218:             for j in mass_dic.keys():
253:         for i in atom_list:221:         for i in atom_list:
254:             for j in mass_dic.keys():222:             for j in mass_dic.keys():
255:                 if j == self.atom_dic[i][0][0]:223:                 if j == self.atom_dic[i][0][0]:
256:                     Z += mass_dic[j] * self.atom_dic[i][5]224:                     Z += mass_dic[j] * self.atom_dic[i][5]
257:         X = X / mass_res225:         X = X / mass_res
258:         Y = Y / mass_res226:         Y = Y / mass_res
259:         Z = Z / mass_res227:         Z = Z / mass_res
260:         print 'Mass weighted centre: ', (X, Y, Z)228:         print 'Mass weighted centre: ', (X, Y, Z)
261:         return (X, Y, Z)229:         return (X, Y, Z)
262: 230: 
263:     # returns all residues of a list of type 'name'231: 
 232:     #returns all residues of a list of type 'name'
264:     def get_all_res_name(self, res_input, name):233:     def get_all_res_name(self, res_input, name):
265:         res_list = []234:         res_list = []
266:         for i in res_input:235:         for i in res_input:
267:             if name == self.get_residue_name(i):236:             if name == self.get_residue_name(i):
268:                 res_list.append(i)237:                 res_list.append(i)
269:         return res_list238:         return res_list
270: 239: 
271:     # atom id from atom name given a residue id240:     #atom id from atom name given a residue id
272:     def get_atom_id_from_name(self, atom_name, res_num):241:     def get_atom_id_from_name(self, atom_name, res_num):
273:         for i in self.res_dic[res_num]:242:         for i in self.res_dic[res_num]:
274:             if self.atom_dic[i][0] == atom_name:243:             if self.atom_dic[i][0] == atom_name:
275:                 return i244:                 return i
276: 245: 
277:     # dot product of two vectors246:     #dot product of two vectors
278:     def dot_product(self, vector1, vector2):247:     def dot_product(self, vector1, vector2):
279:         return (vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2])248:         return (vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2])
280: 249: 
281:     # magnitude of a vector250:     #magnitude of a vector
282:     def magnitude(self, vector):251:     def magnitude(self, vector):
283:         return math.sqrt(vector[0] ** 2 + vector[1] ** 2 + vector[2] ** 2)252:         return math.sqrt(vector[0] ** 2 + vector[1] ** 2 + vector[2] ** 2)
284: 253: 
285:     # angle between three atoms254:     #angle between three atoms
286:     def get_angle(self, (atom1, atom2, atom3)):255:     def get_angle(self, (atom1, atom2, atom3)):
287:         vector1 = (self.atom_dic[atom2][3] - self.atom_dic[atom1][3], self.atom_dic[atom2][4] - self.atom_dic[atom1][4]256:         vector1 = (self.atom_dic[atom2][3] - self.atom_dic[atom1][3], self.atom_dic[atom2][4] - self.atom_dic[atom1][4]
288:                    , self.atom_dic[atom2][5] - self.atom_dic[atom1][5])257:                    , self.atom_dic[atom2][5] - self.atom_dic[atom1][5])
289:         vector2 = (self.atom_dic[atom3][3] - self.atom_dic[atom1][3], self.atom_dic[atom3][4] - self.atom_dic[atom1][4]258:         vector2 = (self.atom_dic[atom3][3] - self.atom_dic[atom1][3], self.atom_dic[atom3][4] - self.atom_dic[atom1][4]
290:                    , self.atom_dic[atom3][5] - self.atom_dic[atom1][5])259:                    , self.atom_dic[atom3][5] - self.atom_dic[atom1][5])
291:         dot_product = self.dot_product(vector1, vector2)260:         dot_product = self.dot_product(vector1, vector2)
292:         magnitudes = self.magnitude(vector1) * self.magnitude(vector2)261:         magnitudes = self.magnitude(vector1) * self.magnitude(vector2)
293:         angle = math.acos(dot_product / magnitudes)262:         angle = math.acos(dot_product / magnitudes)
294:         angle = math.degrees(angle)263:         angle = math.degrees(angle)
295:         if (angle > 0.00 and angle <= 180.00):264:         if (angle > 0.00 and angle <= 180.00):
296:             return angle265:             return angle
297:         else:266:         else:
298:             return angle - 180.00267:             return angle - 180.00
299: 268: 
300:     # checks a list of atoms whether it is linear, creates a triple list of all possible combinations and then computes angles269:     #checks a list of atoms whether it is linear, creates a triple list of all possible combinations and then computes angles
301:     def check_linear(self, atom_list, tolerance):270:     def check_linear(self, atom_list, tolerance):
302:         triple_list = []  # list of atom triples to check for linearity271:         triple_list = []  # list of atom triples to check for linearity
303:         for atom1 in atom_list:272:         for atom1 in atom_list:
304:             for atom2 in atom_list:273:             for atom2 in atom_list:
305:                 for atom3 in atom_list:274:                 for atom3 in atom_list:
306:                     if (atom1 != atom2) and (atom1 != atom3) and (atom2 != atom3):275:                     if (atom1 != atom2) and (atom1 != atom3) and (atom2 != atom3):
307:                         triple_list.append((atom1, atom2, atom3))276:                         triple_list.append((atom1, atom2, atom3))
308:                     else:277:                     else:
309:                         pass278:                         pass
310:         for triple in triple_list:279:         for triple in triple_list:
311:             if (self.get_angle(triple) > (0.00 + tolerance)) and (self.get_angle(triple) < (180.00 - tolerance)):280:             if (self.get_angle(triple) > (0.00 + tolerance)) and (self.get_angle(triple) < (180.00 - tolerance)):
312:                 return False281:                 return False
313:         return True282:         return True
314: 283: 
315:     # returns a dictionary including all the information needed for the localised rigidification scheme,284:     #returns a dictionary including all the information needed for the localised rigidification scheme,
316:     # if an extension is made to the default options for rigidification, this must change the local_scheme list!285:     #if an extension is made to the default options for rigidification, this must change the local_scheme list!
317:     # the input scheme must be a list of 0,1,... for the possible options within the rigidification for a given residue286:     #the input scheme must be a list of 0,1,... for the possible options within the rigidification for a given residue
318:     # any additional scheme can be added below that will be available on default if chosen287:     #any additional scheme can be added below that will be available on default if chosen
319:     def get_rigid_for_local(self, res_list, local_scheme, atoms_used):288:     def get_rigid_for_local(self, res_list, local_scheme, atoms_used):
320:         # local_scheme=['PRO','ARG',('HIS','HIE','HID','HIP'),'LYS','ASP','ASN','GLU','GLN','PHE','TYR','TRP']289:         #local_scheme=['PRO','ARG',('HIS','HIE','HID','HIP'),'LYS','ASP','ASN','GLU','GLN','PHE','TYR','TRP']
321:         groups_res = {}290:         groups_res = {}
322:         i = 1291:         i = 1
323:         # all following cases are for a single residue name with given patterns292:         #all following cases are for a single residue name with given patterns
324:         # when adding additional ones check, that atoms can only be chosen once as no atom can be part of more than one rigid body293:         #when adding additional ones check, that atoms can only be chosen once as no atom can be part of more than one rigid body
325:         if local_scheme[0] == 1:294:         if local_scheme[0] == 1:
326:             l_res = self.get_all_res_name(res_list, 'PRO')295:             l_res = self.get_all_res_name(res_list, 'PRO')
327:             for j in l_res:296:             for j in l_res:
328:                 atom_list = []297:                 atom_list = []
329:                 for k in ['C', 'O', 'N', 'CD', 'CG', 'CB', 'CA']:298:                 for k in ['C', 'O', 'N', 'CD', 'CG', 'CB', 'CA']:
330:                     atom_list.append(self.get_atom_id_from_name(k, int(j)))299:                     atom_list.append(self.get_atom_id_from_name(k, int(j)))
331:                 groups_res[i] = [j, 'PRO'] + atom_list300:                 groups_res[i] = [j, 'PRO'] + atom_list
332:                 atoms_used += atom_list301:                 atoms_used += atom_list
333:                 i += 1302:                 i += 1
334:             l_res = []303:             l_res = []
480:                 for k in ['CE2', 'CZ2', 'HZ2', 'CH2', 'HH2', 'CZ3', 'CE3', 'HE3', 'CD2']:449:                 for k in ['CE2', 'CZ2', 'HZ2', 'CH2', 'HH2', 'CZ3', 'CE3', 'HE3', 'CD2']:
481:                     atom_list2.append(self.get_atom_id_from_name(k, int(j)))450:                     atom_list2.append(self.get_atom_id_from_name(k, int(j)))
482:                 groups_res[i] = [j, 'TRP'] + atom_list1451:                 groups_res[i] = [j, 'TRP'] + atom_list1
483:                 groups_res[i + 1] = [j, 'TRP'] + atom_list2452:                 groups_res[i + 1] = [j, 'TRP'] + atom_list2
484:                 atoms_used += atom_list1453:                 atoms_used += atom_list1
485:                 atoms_used += atom_list2454:                 atoms_used += atom_list2
486:                 i += 2455:                 i += 2
487:             l_res = []456:             l_res = []
488:         return (groups_res, atoms_used)457:         return (groups_res, atoms_used)
489: 458: 
490:     def loc_NA(self, res_list, local_scheme, atoms_used): 
491:         # local_scheme=[G_large, G_small, A_large, A_small, C_large, C_small, T_large, T_small, U_large, U_small]) 
492:         groups_res = {} 
493:         i = 1 
494:         # all following cases are for a single residue name with given patterns 
495:         # when adding additional ones check, that atoms can only be chosen once as no atom can be part of more than one rigid body 
496:         if local_scheme[0]: 
497:             l_res = [] 
498:             for RA_name in ['G', 'G3', 'G5', 'GN', 'DG', 'DG3', 'DG5', 'DGN', 'RG', 'RG3', 'RG5', 'RGN']: 
499:                 l_res += self.get_all_res_name(res_list, RA_name) 
500:             for j in l_res: 
501:                 atom_list = [] 
502:                 for k in ['N1', 'N9', 'C8', 'N7', 'C5', 'C6', 'H8', 'O6', 'C2', 'N3', 'H1', 'N2', 'H21', 'H22']: 
503:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
504:                 groups_res[i] = [j, 'G'] + atom_list 
505:                 atoms_used += atom_list 
506:                 i += 1 
507:             l_res = [] 
508:         if local_scheme[1]: 
509:             l_res = [] 
510:             for RA_name in ['G', 'G3', 'G5', 'GN', 'DG', 'DG3', 'DG5', 'DGN', 'RG', 'RG3', 'RG5', 'RGN']: 
511:                 l_res += self.get_all_res_name(res_list, RA_name) 
512:             for j in l_res: 
513:                 atom_list = [] 
514:                 for k in ['N9', 'C8', 'N7', 'C5', 'C6', 'H8', 'C2', 'N3', 'H1', 'N2']: 
515:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
516:                 groups_res[i] = [j, 'G'] + atom_list 
517:                 atoms_used += atom_list 
518:                 i += 1 
519:             l_res = [] 
520:         if local_scheme[2]: 
521:             l_res = [] 
522:             for RA_name in ['A', 'A3', 'A5', 'AN', 'DA', 'DA3', 'DA5', 'DAN', 'RA', 'RA3', 'RA5', 'RAN']: 
523:                 l_res += self.get_all_res_name(res_list, RA_name) 
524:             for j in l_res: 
525:                 atom_list = [] 
526:                 for k in ['N1', 'C2', 'H2', 'N3', 'C4', 'C5', 'C6', 'N7', 'C8', 'N9', 'H8', 'N6', 'H61', 'H62']: 
527:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
528:                 groups_res[i] = [j, 'A'] + atom_list 
529:                 atoms_used += atom_list 
530:                 i += 1 
531:             l_res = [] 
532:         if local_scheme[3]: 
533:             l_res = [] 
534:             for RA_name in ['A', 'A3', 'A5', 'AN', 'DA', 'DA3', 'DA5', 'DAN', 'RA', 'RA3', 'RA5', 'RAN']: 
535:                 l_res += self.get_all_res_name(res_list, RA_name) 
536:             for j in l_res: 
537:                 atom_list = [] 
538:                 for k in ['N1', 'C2', 'H2', 'N3', 'C4', 'C5', 'C6', 'N7', 'C8', 'N9', 'H8', 'N6']: 
539:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
540:                 groups_res[i] = [j, 'A'] + atom_list 
541:                 atoms_used += atom_list 
542:                 i += 1 
543:             l_res = [] 
544:         if local_scheme[4]: 
545:             l_res = [] 
546:             for RA_name in ['C', 'C3', 'C5', 'CN', 'DC', 'DC3', 'DC5', 'DCN', 'RC', 'RC3', 'RC5', 'RCN']: 
547:                 l_res += self.get_all_res_name(res_list, RA_name) 
548:             for j in l_res: 
549:                 atom_list = [] 
550:                 for k in ['N1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H5', 'H6', 'O2', 'N4', 'H41', 'H42']: 
551:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
552:                 groups_res[i] = [j, 'C'] + atom_list 
553:                 atoms_used += atom_list 
554:                 i += 1 
555:             l_res = [] 
556:         if local_scheme[5]: 
557:             l_res = [] 
558:             for RA_name in ['C', 'C3', 'C5', 'CN', 'DC', 'DC3', 'DC5', 'DCN', 'RC', 'RC3', 'RC5', 'RCN']: 
559:                 l_res += self.get_all_res_name(res_list, RA_name) 
560:             for j in l_res: 
561:                 atom_list = [] 
562:                 for k in ['N1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H5', 'H6', 'N4']: 
563:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
564:                 groups_res[i] = [j, 'C'] + atom_list 
565:                 atoms_used += atom_list 
566:                 i += 1 
567:             l_res = [] 
568:         if local_scheme[6]: 
569:             l_res = [] 
570:             for RA_name in ['DT', 'DT3', 'DT5', 'DTN']: 
571:                 l_res += self.get_all_res_name(res_list, RA_name) 
572:             for j in l_res: 
573:                 atom_list = [] 
574:                 for k in ['N1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H5', 'H3', 'H6', 'O4', 'O2', 'C7', 'H71', 'H72', 'H73']: 
575:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
576:                 groups_res[i] = [j, 'T'] + atom_list 
577:                 atoms_used += atom_list 
578:                 i += 1 
579:             l_res = [] 
580:         if local_scheme[7]: 
581:             l_res = [] 
582:             for RA_name in ['DT', 'DT3', 'DT5', 'DTN']: 
583:                 l_res += self.get_all_res_name(res_list, RA_name) 
584:             for j in l_res: 
585:                 atom_list = [] 
586:                 for k in ['N1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H5', 'H6', 'C7', 'H71', 'H72', 'H73']: 
587:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
588:                 groups_res[i] = [j, 'T'] + atom_list 
589:                 atoms_used += atom_list 
590:                 i += 1 
591:             l_res = [] 
592:         if local_scheme[8]: 
593:             l_res = [] 
594:             for RA_name in ['U', 'U3', 'U5', 'UN', 'RU', 'RU3', 'RU5', 'RUN']: 
595:                 l_res += self.get_all_res_name(res_list, RA_name) 
596:             for j in l_res: 
597:                 atom_list = [] 
598:                 for k in ['N1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H5', 'H3', 'H6', 'O4', 'O2']: 
599:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
600:                 groups_res[i] = [j, 'U'] + atom_list 
601:                 atoms_used += atom_list 
602:                 i += 1 
603:             l_res = [] 
604:         if local_scheme[9]: 
605:             l_res = [] 
606:             for RA_name in ['U', 'U3', 'U5', 'UN', 'RU', 'RU3', 'RU5', 'RUN']: 
607:                 l_res += self.get_all_res_name(res_list, RA_name) 
608:             for j in l_res: 
609:                 atom_list = [] 
610:                 for k in ['N1', 'C2', 'N3', 'C4', 'C5', 'C6', 'H5', 'H6']: 
611:                     atom_list.append(self.get_atom_id_from_name(k, int(j))) 
612:                 groups_res[i] = [j, 'U'] + atom_list 
613:                 atoms_used += atom_list 
614:                 i += 1 
615:             l_res = [] 
616:         return (groups_res, atoms_used) 
617: 459: 
618: 460: #class to implement a live feed in pymol - not yet tested
619: # class to implement a live feed in pymol - not yet tested 
620: class PymolView():461: class PymolView():
621:     # initialises pymol and starts it without GUI, loads a given file and sets it to a white cartoon representation462:     #initialises pymol and starts it without GUI, loads a given file and sets it to a white cartoon representation
622:     def __init__(self, file_name):463:     def __init__(self, file_name):
623:         print 'Launching Pymol'464:         print 'Launching Pymol'
624:         pymol.finish_launching()465:         pymol.finish_launching()
625:         pymol.cmd.load(file_name)466:         pymol.cmd.load(file_name)
626:         print '%s loaded in Pymol' % file_name467:         print '%s loaded in Pymol' % file_name
627: 468: 
628:         # show molecule as cartoon and hide lines469:         #show molecule as cartoon and hide lines
629:         pymol.cmd.show('cartoon')470:         pymol.cmd.show('cartoon')
630:         pymol.cmd.hide('lines')471:         pymol.cmd.hide('lines')
631:         pymol.cmd.set('cartoon_color', 0)472:         pymol.cmd.set('cartoon_color', 0)
632: 473: 
633:     def set_colour_range(self, start, finish, colour):474:     def set_colour_range(self, start, finish, colour):
634:         # setting residue range to any colour475:         #setting residue range to any colour
635:         pymol.cmd.set('cartoon_color', colour, 'resi %s-%s' % (str(start), str(finish)))476:         pymol.cmd.set('cartoon_color', colour, 'resi %s-%s' % (str(start), str(finish)))
636: 477: 
637:     # reset all of the molecule to a white colour478:     #reset all of the molecule to a white colour
638:     # def reset_colour(self,length):479:     #def reset_colour(self,length):
639:     #    pymol.cmd.set('cartoon_color',0,'resi 1-%s'%str(length))480:     #    pymol.cmd.set('cartoon_color',0,'resi 1-%s'%str(length))
640: 481: 
641:     def set_colour_res(self, res, colour):482:     def set_colour_res(self, res, colour):
642:         pymol.cmd.set('cartoon_color', colour, 'resi %s' % str(res))483:         pymol.cmd.set('cartoon_color', colour, 'resi %s' % str(res))
643: 484: 
644:     def __del__(self):485:     def __del__(self):
645:         print 'Quit pymol'486:         print 'Quit pymol'
646:         pymol.cmd.quit()487:         pymol.cmd.quit()
647: 488: 
648:     def apply_colour_scheme(self, atomistic_list, local_list, colour_atomistic, colour_local):489:     def apply_colour_scheme(self, atomistic_list, local_list, colour_atomistic, colour_local):
649:         # self.reset_colour()490:         #self.reset_colour()
650:         for i in atomistic_list:491:         for i in atomistic_list:
651:             self.set_colour_res(i, colour_atomistic)492:             self.set_colour_res(i, colour_atomistic)
652:         for j in local_list:493:         for j in local_list:
653:             self.set_colour_res(j, colour_local)494:             self.set_colour_res(j, colour_local)
654: 495: 
655: 496: 
656: ###FUNCTION TO PARSE THE PDB AND COORDS.INPCRD FILE ###################################################################497: ###FUNCTION TO PARSE THE PDB AND COORDS.INPCRD FILE ###################################################################
657: def reading_pdb_inpcrd(pdb, inpcrd):498: def reading_pdb_inpcrd(pdb, inpcrd):
658:     dic = {}499:     dic = {}
659:     try:500:     try:
734: 575: 
735: 576: 
736: def removal_selection(frozen, local, atomistic, sel_local, sel_atomistic):577: def removal_selection(frozen, local, atomistic, sel_local, sel_atomistic):
737:     frozen = merge_lists(frozen, sel_local)578:     frozen = merge_lists(frozen, sel_local)
738:     frozen = merge_lists(frozen, sel_atomistic)579:     frozen = merge_lists(frozen, sel_atomistic)
739:     local = remove_selection_from_list(local, sel_local)580:     local = remove_selection_from_list(local, sel_local)
740:     atomistic = remove_selection_from_list(atomistic, sel_atomistic)581:     atomistic = remove_selection_from_list(atomistic, sel_atomistic)
741:     return (sorted(atomistic), sorted(local), sorted(frozen))582:     return (sorted(atomistic), sorted(local), sorted(frozen))
742: 583: 
743: 584: 
744: # FUNCTION TO SELECT RESIDUES ########################################################################################585: ###FUNCTION TO SELECT RESIDUES ########################################################################################
745: def choose_res(mol):586: def choose_res(mol):
746:     print '1 - choose sphere/ellipsoid of unfrozen atoms and add localised rigidification'587:     print '1 - choose sphere/ellipsoid of unfrozen atoms and add localised rigidification'
747:     print '2 - choose a range of residues or single residues to unfreeze or locally rigidify'588:     print '2 - choose a range of residues or single residues to unfreeze or locally rigidify'
748:     print589:     print
749:     log = []590:     log = []
750:     try:591:     try:
751:         method = int(raw_input('Method:  '))592:         method = int(raw_input('Method:  '))
752:     except IOError:593:     except IOError:
753:         method = 0594:         method = 0
754:     unfrozen_res = []595:     unfrozen_res = []
755:     local_res = []596:     local_res = []
756: 597: 
757:     if method == 0:598:     if method == 0:
758:         print 'Invalid input'599:         print 'Invalid input'
759: 600: 
760:     # ellipsoidal or spherical volume601:     # ellipsoidal or spherical volume
761:     elif method == 1:602:     elif method == 1:
762:         while 1:  # choose the geometry: currently ellipsoidal 'e' or spherical 's'603:         while 1:  #choose the geometry: currently ellipsoidal 'e' or spherical 's'
763:             log.append('Method: sphere / ellipsoid')604:             log.append('Method: sphere / ellipsoid')
764:             while 1:605:             while 1:
765:                 centre_method = raw_input('Mass weighted centre of residue (1) or atom in a selected residue (2)?  ')606:                 centre_method = raw_input('Mass weighted centre of residue (1) or atom in a selected residue (2)?  ')
766:                 try:607:                 try:
767:                     if int(centre_method) == 1:608:                     if int(centre_method) == 1:
768:                         residue = raw_input('Residue:  ')609:                         residue = raw_input('Residue:  ')
769:                         centre_xyz = mol.mass_weighted_centre(int(residue))610:                         centre_xyz = mol.mass_weighted_centre(int(residue))
770:                         log.append('Mass-weighted centre: Residue' + residue + ' ' + mol.get_residue_name(611:                         log.append('Mass-weighted centre: Residue' + residue + ' ' + mol.get_residue_name(
771:                             int(residue)) + '  ' + str(centre_xyz[0]) + '  ' + str(centre_xyz[1]) + '  ' + str(612:                             int(residue)) + '  ' + str(centre_xyz[0]) + '  ' + str(centre_xyz[1]) + '  ' + str(
772:                             centre_xyz[2]))613:                             centre_xyz[2]))
797:                     l_atom = mol.get_atoms_in_ellipsoid(centre_xyz, float(l_radii[0]), float(l_radii[1]),638:                     l_atom = mol.get_atoms_in_ellipsoid(centre_xyz, float(l_radii[0]), float(l_radii[1]),
798:                                                         float(l_radii[2]))639:                                                         float(l_radii[2]))
799:                     l_res = mol.transform_atomlist_to_reslist(l_atom)640:                     l_res = mol.transform_atomlist_to_reslist(l_atom)
800:                     unfrozen_res += l_res641:                     unfrozen_res += l_res
801:                     break642:                     break
802:                 else:643:                 else:
803:                     print 'Wrong input'644:                     print 'Wrong input'
804:             except ValueError:645:             except ValueError:
805:                 print 'Invalid input'646:                 print 'Invalid input'
806:         print647:         print
807:         while 1:  # geometry is same as for atomistic region, if a change is wanted employ the above648:         while 1:  #geometry is same as for atomistic region, if a change is wanted employ the above
808:             check = raw_input('Add localised rigidification (y/n)?  ')649:             check = raw_input('Add localised rigidification (y/n)?  ')
809:             if check == 'yes' or check == 'y':650:             if check == 'yes' or check == 'y':
810:                 log.append('Local rigidification added')651:                 log.append('Local rigidification added')
811:                 radius_input = raw_input(652:                 radius_input = raw_input(
812:                     'Radius of atomistic region (enter three values for an ellipsoid and one for a sphere):  ')653:                     'Radius of atomistic region (enter three values for an ellipsoid and one for a sphere):  ')
813:                 l_radii_outer = radius_input.split()654:                 l_radii_outer = radius_input.split()
814:                 try:655:                 try:
815:                     if len(l_radii_outer) == 1:656:                     if len(l_radii_outer) == 1:
816:                         log.append('Outer radius: ' + l_radii_outer[0])657:                         log.append('Outer radius: ' + l_radii_outer[0])
817:                         l_atom_out = mol.get_atoms_in_sphere(centre_xyz, float(l_radii_outer[0]))658:                         l_atom_out = mol.get_atoms_in_sphere(centre_xyz, float(l_radii_outer[0]))
832:                             if i not in l_res:673:                             if i not in l_res:
833:                                 local_res.append(i)674:                                 local_res.append(i)
834:                         break675:                         break
835:                     else:676:                     else:
836:                         print 'Wrong input'677:                         print 'Wrong input'
837:                 except ValueError:678:                 except ValueError:
838:                     print 'Invalid input'679:                     print 'Invalid input'
839:             else:680:             else:
840:                 break681:                 break
841: 682: 
842:     # rigidification using a range of residues, sanity check in function that chooses683:     #rigidification using a range of residues, sanity check in function that chooses
843:     elif method == 2:684:     elif method == 2:
844:         log.append('Method:  range selection')685:         log.append('Method:  range selection')
845:         check = raw_input('Define atomistic residues (y/n)?  ')686:         check = raw_input('Define atomistic residues (y/n)?  ')
846:         if check == 'yes' or check == 'y':687:         if check == 'yes' or check == 'y':
847:             res = raw_input(688:             res = raw_input(
848:                 '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 num2689:                 '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
849:             res = res.split()690:             res = res.split()
850:             i = 0691:             i = 0
851:             log_res = ''692:             log_res = ''
852:             while i < len(res):693:             while i < len(res):
853:                 if res[i] == 'r':  # if a range is indicated694:                 if res[i] == 'r':  #if a range is indicated
854:                     i += 1695:                     i += 1
855:                     start_range = int(res[i])  # finds starting residue696:                     start_range = int(res[i])  #finds starting residue
856:                     i += 1697:                     i += 1
857:                     end_range = int(res[i])  # finds final residue698:                     end_range = int(res[i])  #finds final residue
858:                     for j in range(start_range, end_range + 1):699:                     for j in range(start_range, end_range + 1):
859:                         unfrozen_res.append(j)700:                         unfrozen_res.append(j)
860:                         log_res += str(j) + ' '701:                         log_res += str(j) + ' '
861:                     i += 1702:                     i += 1
862:                 else:703:                 else:
863:                     try:704:                     try:
864:                         unfrozen_res.append(int(res[i]))705:                         unfrozen_res.append(int(res[i]))
865:                         log_res += res[i] + ' '706:                         log_res += res[i] + ' '
866:                         i += 1707:                         i += 1
867:                     except ValueError:708:                     except ValueError:
868:                         i += 1709:                         i += 1
869:             log.append('Residues entered as atomistically: ' + log_res)710:             log.append('Residues entered as atomistically: ' + log_res)
870:         check = raw_input('Locally rigidify (y/n)?  ')711:         check = raw_input('Locally rigidify (y/n)?  ')
871:         if check == 'yes' or check == 'y':712:         if check == 'yes' or check == 'y':
872:             res = raw_input(713:             res = raw_input(
873:                 '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 num2714:                 '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
874:             res = res.split()715:             res = res.split()
875:             i = 0716:             i = 0
876:             log_res = ''717:             log_res = ''
877:             while i < len(res):718:             while i < len(res):
878:                 if res[i] == 'r':  # if a range is indicated719:                 if res[i] == 'r':  #if a range is indicated
879:                     i += 1720:                     i += 1
880:                     start_range = int(res[i])  # finds starting residue721:                     start_range = int(res[i])  #finds starting residue
881:                     i += 1722:                     i += 1
882:                     end_range = int(res[i])  # finds final residue723:                     end_range = int(res[i])  #finds final residue
883:                     i += 1724:                     i += 1
884:                     for j in range(start_range, end_range + 1):725:                     for j in range(start_range, end_range + 1):
885:                         if j <= mol.num_res():726:                         if j <= mol.num_res():
886:                             local_res.append(j)727:                             local_res.append(j)
887:                             log_res += str(j) + ' '728:                             log_res += str(j) + ' '
888:                         else:729:                         else:
889:                             break730:                             break
890:                 else:731:                 else:
891:                     try:732:                     try:
892:                         if int(res[i]) <= mol.num_res():733:                         if int(res[i]) <= mol.num_res():
903: 744: 
904: 745: 
905: def remove_res(mol):746: def remove_res(mol):
906:     unfrozen_res = []747:     unfrozen_res = []
907:     local_res = []748:     local_res = []
908:     log = []749:     log = []
909:     log.append('Removal of selected residues')750:     log.append('Removal of selected residues')
910:     check = raw_input('Remove from the atomistic residues (y/n)?  ')751:     check = raw_input('Remove from the atomistic residues (y/n)?  ')
911:     if check == 'yes' or check == 'y':752:     if check == 'yes' or check == 'y':
912:         res = raw_input(753:         res = raw_input(
913:             '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 num2754:             '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
914:         res = res.split()755:         res = res.split()
915:         i = 0756:         i = 0
916:         log_res = ''757:         log_res = ''
917:         while i < len(res):758:         while i < len(res):
918:             if res[i] == 'r':  # if a range is indicated759:             if res[i] == 'r':  #if a range is indicated
919:                 i += 1760:                 i += 1
920:                 start_range = int(res[i])  # finds starting residue761:                 start_range = int(res[i])  #finds starting residue
921:                 i += 1762:                 i += 1
922:                 end_range = int(res[i])  # finds final residue763:                 end_range = int(res[i])  #finds final residue
923:                 for j in range(start_range, end_range + 1):764:                 for j in range(start_range, end_range + 1):
924:                     unfrozen_res.append(j)765:                     unfrozen_res.append(j)
925:                     log_res += str(j) + ' '766:                     log_res += str(j) + ' '
926:                 i += 1767:                 i += 1
927:             else:768:             else:
928:                 try:769:                 try:
929:                     unfrozen_res.append(int(res[i]))770:                     unfrozen_res.append(int(res[i]))
930:                     log_res += res[i] + ' '771:                     log_res += res[i] + ' '
931:                     i += 1772:                     i += 1
932:                 except ValueError:773:                 except ValueError:
933:                     i += 1774:                     i += 1
934:         log.append('Residues entered to be removed: ' + log_res + "\n")775:         log.append('Residues entered to be removed: ' + log_res + "\n")
935:     check = raw_input('Remove from the local selection (y/n)?  ')776:     check = raw_input('Remove from the local selection (y/n)?  ')
936:     if check == 'yes' or check == 'y':777:     if check == 'yes' or check == 'y':
937:         res = raw_input(778:         res = raw_input(
938:             '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 num2779:             '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
939:         res = res.split()780:         res = res.split()
940:         i = 0781:         i = 0
941:         log_res = ''782:         log_res = ''
942:         while i < len(res):783:         while i < len(res):
943:             if res[i] == 'r':  # if a range is indicated784:             if res[i] == 'r':  #if a range is indicated
944:                 i += 1785:                 i += 1
945:                 start_range = int(res[i])  # finds starting residue786:                 start_range = int(res[i])  #finds starting residue
946:                 i += 1787:                 i += 1
947:                 end_range = int(res[i])  # finds final residue788:                 end_range = int(res[i])  #finds final residue
948:                 i += 1789:                 i += 1
949:                 for j in range(start_range, end_range + 1):790:                 for j in range(start_range, end_range + 1):
950:                     if j <= mol.num_res():791:                     if j <= mol.num_res():
951:                         local_res.append(j)792:                         local_res.append(j)
952:                         log_res += str(j) + ' '793:                         log_res += str(j) + ' '
953:                     else:794:                     else:
954:                         break795:                         break
955:             else:796:             else:
956:                 try:797:                 try:
957:                     if int(res[i]) <= mol.num_res():798:                     if int(res[i]) <= mol.num_res():
958:                         local_res.append(int(res[i]))799:                         local_res.append(int(res[i]))
959:                         log_res += res[i] + ' '800:                         log_res += res[i] + ' '
960:                         i += 1801:                         i += 1
961:                 except ValueError:802:                 except ValueError:
962:                     i += 1803:                     i += 1
963:         log.append('Residues entered to be removed from local rigidification: ' + log_res + "\n")804:         log.append('Residues entered to be removed from local rigidification: ' + log_res + "\n")
964:     return (unfrozen_res, local_res, log)805:     return (unfrozen_res, local_res, log)
965: 806: 
966: 807: 
967: # FUNCTIONS FOR GROUPROTATION FILE ###################################################################################808: ###FUNCTIONS FOR GROUPROTATION FILE ###################################################################################
968: def get_prob_amb():809: def get_prob_amb():
969:     while 1:810:     while 1:
970:         inp_p = 'Probability: '811:         inp_p = 'Probability: '
971:         inp_a = 'Amplitude: '812:         inp_a = 'Amplitude: '
972:         try:813:         try:
973:             probability = float(inp_p)814:             probability = float(inp_p)
974:             amplitude = float(inp_a)815:             amplitude = float(inp_a)
975:             if (probability > 1.0) or (probability < 0.0):816:             if (probability > 1.0) or (probability < 0.0):
976:                 print 'Invalid probability'817:                 print 'Invalid probability'
977:             else:818:             else:
978:                 break819:                 break
979:         except ValueError:820:         except ValueError:
980:             print 'Both values need to be numbers'821:             print 'Both values need to be numbers'
981:     return (probability, amplitude)822:     return (probability,amplitude)
982: 823: 
983: 824: 
984: # VARIABLE DEFINITION ################################################################################################825: ###VARIABLE DEFINITION ################################################################################################
985: atom_dic = {}  # dictionary containing all information sorted by atom number:826: atom_dic = {}  #dictionary containing all information sorted by atom number:
986: # the information is in a tuple: (atomname, residue name, residue number, x, y, z)827:                #the information is in a tuple: (atomname, residue name, residue number, x, y, z)
987: res_dic = {}  # dictionary containing all atoms per residue sorted by residue number828: res_dic = {}  #dictionary containing all atoms per residue sorted by residue number
988: res_atomistic = []  # residue list that are treated fully atomistic829: res_atomistic = []  #residue list that are treated fully atomistic
989: res_local = []  # residue list that are treated locally rigid830: res_local = []  #residue list that are treated locally rigid
990: res_frozen = []  # residue list that are treated fully rigid831: res_frozen = []  #residue list that are treated fully rigid
991: atoms_used = []  # atoms used in rigid bodies, prevents double usage!832: atoms_used = []  #atoms used in rigid bodies, prevents double usage!
992: groups = {}  # dictionary of rigid groups833: groups = {}  #dictionary of rigid groups
993: groups_local = {}834: groups_local = {}
994: num_rbody = 0  # number of rigid bodies defined835: num_rbody = 0  #number of rigid bodies defined
995: 836: 
996: # PARSING OF INPUT FILES #############################################################################################837: ###PARSING OF INPUT FILES #############################################################################################
997: while 1:  # starts input selection838: while 1:  # starts input selection
998:     atom_dic = reading_pdb_inpcrd(pdb_inp, coords_inp)  # opens files and loads them into a dictionary839:     atom_dic = reading_pdb_inpcrd(pdb_inp, coords_inp)  # opens files and loads them into a dictionary
999:     res_dic = create_res_dic(atom_dic)840:     res_dic = create_res_dic(atom_dic)
1000:     if atom_dic != {}:  # if the dictionary contains data the data is read into the molecule class841:     if atom_dic != {}:  # if the dictionary contains data the data is read into the molecule class
1001:         loaded_mol = protein(atom_dic, res_dic)842:         loaded_mol = protein(atom_dic, res_dic)
1002:         print843:         print
1003:         print 'Molecule was loaded'844:         print 'Molecule was loaded'
1004:         print 'Number of residues: ', loaded_mol.num_res()845:         print 'Number of residues: ', loaded_mol.num_res()
1005:         print 'Number of atoms: ', loaded_mol.num_atoms()846:         print 'Number of atoms: ', loaded_mol.num_atoms()
1006:         print847:         print
1011:         pdb_inp = raw_input('PDB file: ')852:         pdb_inp = raw_input('PDB file: ')
1012:         coords_inp = raw_input('Coords input: ')853:         coords_inp = raw_input('Coords input: ')
1013: 854: 
1014: check_file = open('rigid.log', 'w')855: check_file = open('rigid.log', 'w')
1015: check_file.write('Input file for structure: ' + pdb_inp + "\n")856: check_file.write('Input file for structure: ' + pdb_inp + "\n")
1016: check_file.write('Coordinates file used: ' + coords_inp + "\n" + "\n")857: check_file.write('Coordinates file used: ' + coords_inp + "\n" + "\n")
1017: check_file.write('Number of residues:  ' + str(loaded_mol.num_res()) + "\n")858: check_file.write('Number of residues:  ' + str(loaded_mol.num_res()) + "\n")
1018: check_file.write('Number of atoms:  ' + str(loaded_mol.num_atoms()) + "\n" + "\n")859: check_file.write('Number of atoms:  ' + str(loaded_mol.num_atoms()) + "\n" + "\n")
1019: check_file.write('Time:  ' + time.strftime("%a, %d %b %Y %H:%M:%S") + "\n" + "\n")860: check_file.write('Time:  ' + time.strftime("%a, %d %b %Y %H:%M:%S") + "\n" + "\n")
1020: 861: 
1021: # CHECK FOR LIGANDS AND UNKNOWN RESIDUES #############################################################################862: ###CHECK FOR LIGANDS AND UNKNOWN RESIDUES #############################################################################
1022: check_file.write('Ligands and unknown residues:' + "\n")863: check_file.write('Ligands and unknown residues:' + "\n")
1023: ligand_dic = loaded_mol.get_ligands_unknown_res()864: ligand_dic = loaded_mol.get_ligands_unknown_res()
1024: if ligand_dic == {}:865: if ligand_dic == {}:
1025:     print 'There are no ligands or unknown residues in the molecule'866:     print 'There are no ligands or unknown residues in the molecule'
1026:     check_file.write('None' + "\n" + "\n")867:     check_file.write('None' + "\n" + "\n")
1027: else:868: else:
1028:     print 'There are ligands or unknown residues.'869:     print 'There are ligands or unknown residues.'
1029:     ligand_list = ligand_dic.keys()870:     ligand_list = ligand_dic.keys()
1030:     for res_id in ligand_dic.keys():871:     for res_id in ligand_dic.keys():
1031:         print str(res_id) + ' - ' + loaded_mol.res_name_from_id(res_id)872:         print str(res_id) + ' - ' + loaded_mol.res_name_from_id(res_id)
1043:                 try:884:                 try:
1044:                     res_selected = [int(s) for s in ligand_inp.split()]885:                     res_selected = [int(s) for s in ligand_inp.split()]
1045:                     new_ligand_list = []886:                     new_ligand_list = []
1046:                     string_ligands = ''887:                     string_ligands = ''
1047:                     for i in ligand_list:888:                     for i in ligand_list:
1048:                         if i not in res_selected:889:                         if i not in res_selected:
1049:                             new_ligand_list.append(i)890:                             new_ligand_list.append(i)
1050:                             string_ligands += str(i) + ' '891:                             string_ligands += str(i) + ' '
1051:                     res_atomistic += new_ligand_list892:                     res_atomistic += new_ligand_list
1052:                     print 'The following residues will be treated atomistically: ' + string_ligands893:                     print 'The following residues will be treated atomistically: ' + string_ligands
1053:                     check_file.write('The following residues will be treated atomistically: ' +894:                     check_file.write('The following residues will be treated atomistically: ' + string_ligands
1054:                                      string_ligands + "\n" + "\n")895:                                      + "\n" + "\n")
1055:                     break896:                     break
1056:                 except ValueError:897:                 except ValueError:
1057:                     print 'Please enter integer numbers.'898:                     print 'Please enter integer numbers.'
1058:     else:899:     else:
1059:         res_atomistic += ligand_list900:         res_atomistic += ligand_list
1060:         print 'The above listed residues/ligands will be treated fully atomistically.'901:         print 'The above listed residues/ligands will be treated fully atomistically.'
1061:         check_file.write('All are treated fully atomistically.' + "\n" + "\n")902:         check_file.write('All are treated fully atomistically.' + "\n" + "\n")
1062: 903: 
1063: # SELECTION OF RESIDUES FOR FULLY AND LOCALLY RIGID REGIONS ##########################################################904: ###SELECTION OF RESIDUES FOR FULLY AND LOCALLY RIGID REGIONS ##########################################################
1064: print905: print
1065: print906: print
1066: print 'Selection of atomistic and locally rigidified residues'907: print 'Selection of atomistic and locally rigidified residues'
1067: print908: print
1068: check_file.write('Selection of atomistic and locally rigidified residues' + "\n")909: check_file.write('Selection of atomistic and locally rigidified residues' + "\n")
1069: res_frozen = remove_selection_from_list(res_dic.keys(), res_atomistic)910: res_frozen = remove_selection_from_list(res_dic.keys(), res_atomistic)
1070: x = choose_res(loaded_mol)  # initiates the first round of selections911: x = choose_res(loaded_mol)  #initiates the first round of selections
1071: y = merge_new_and_old_res_lists(res_frozen, res_local, res_atomistic, x[1], x[0])912: y = merge_new_and_old_res_lists(res_frozen, res_local, res_atomistic, x[1], x[0])
1072: res_frozen = y[2]913: res_frozen = y[2]
1073: res_local = y[1]914: res_local = y[1]
1074: res_atomistic = y[0]915: res_atomistic = y[0]
1075: log_string = ''916: log_string = ''
1076: for strings in x[2]:917: for strings in x[2]:
1077:     log_string += strings + "\n"918:     log_string += strings + "\n"
1078: check_file.write(log_string + "\n" + "\n")919: check_file.write(log_string + "\n" + "\n")
1079: 920: 
1080: while 1:  # allows to choose more residues with different methods921: while 1:  #allows to choose more residues with different methods
1081:     print922:     print
1082:     print 'The following residues are not frozen (atomistic):'923:     print 'The following residues are not frozen (atomistic):'
1083:     print res_atomistic924:     print res_atomistic
1084:     print925:     print
1085:     print 'The following residues are selected to be locally rigidified:'926:     print 'The following residues are selected to be locally rigidified:'
1086:     print res_local927:     print res_local
1087:     print928:     print
1088:     print 'The following residues are selected to be fully rigidified:'929:     print 'The following residues are selected to be fully rigidified:'
1089:     print res_frozen930:     print res_frozen
1090:     print931:     print
1091:     check = raw_input('Enter more residues (1) or remove residues from previous selection (2)? Any other key to exit. ')932:     check = raw_input('Enter more residues (1) or remove residues from previous selection (2)? Any other key to exit. ')
1092:     try:933:     try:
1093:         value = int(934:         value = int(
1094:             check)  # to enter more residues the choose_res function is called again with the same options as before935:             check)  #to enter more residues the choose_res function is called again with the same options as before
1095:         if int(check) == 1:936:         if int(check) == 1:
1096:             x = choose_res(loaded_mol)937:             x = choose_res(loaded_mol)
1097:             y = merge_new_and_old_res_lists(res_frozen, res_local, res_atomistic, x[1], x[0])938:             y = merge_new_and_old_res_lists(res_frozen, res_local, res_atomistic, x[1], x[0])
1098:             res_frozen = y[2]939:             res_frozen = y[2]
1099:             res_local = y[1]940:             res_local = y[1]
1100:             res_atomistic = y[0]941:             res_atomistic = y[0]
1101:             log_string = ''942:             log_string = ''
1102:             for strings in x[2]:943:             for strings in x[2]:
1103:                 log_string += strings + "\n"944:                 log_string += strings + "\n"
1104:             check_file.write(log_string + "\n" + "\n")945:             check_file.write(log_string + "\n" + "\n")
1105:         elif int(check) == 2:946:         elif int(check) == 2:
1106:             x = remove_res(loaded_mol)947:             x = remove_res(loaded_mol)
1107:             y = removal_selection(res_frozen, res_local, res_atomistic, x[1],948:             y = removal_selection(res_frozen, res_local, res_atomistic, x[1],
1108:                                   x[0])  # the user enters a selection for removal949:                                   x[0])  #the user enters a selection for removal
1109:             res_frozen = y[2]950:             res_frozen = y[2]
1110:             res_local = y[1]951:             res_local = y[1]
1111:             res_atomistic = y[0]952:             res_atomistic = y[0]
1112:             check_file.write('Remove residues' + "\n" + 'Atomistic residues removed:' + "\n")953:             check_file.write('Remove residues' + "\n" + 'Atomistic residues removed:' + "\n")
1113:             string = ''954:             string = ''
1114:             for i in x[0]:955:             for i in x[0]:
1115:                 string = string + str(i) + '  '956:                 string = string + str(i) + '  '
1116:             check_file.write(string + "\n" + "\n")957:             check_file.write(string + "\n" + "\n")
1117:             check_file.write('Locally rigidified residues removed:' + "\n")958:             check_file.write('Locally rigidified residues removed:' + "\n")
1118:             string = ''959:             string = ''
1119:             for i in x[1]:960:             for i in x[1]:
1120:                 string = string + str(i) + '  '961:                 string = string + str(i) + '  '
1121:             check_file.write(string + "\n" + "\n")962:             check_file.write(string + "\n" + "\n")
1122:         else:963:         else:
1123:             break964:             break
1124: 965: 
1125:     except ValueError:966:     except ValueError:
1126:         break967:         break
1127: check_file.write('Selection completed' + "\n" + "\n")968: check_file.write('Selection completed' +  "\n" + "\n")
1128: string_atomistic = ''969: string_atomistic = ''
1129: for res in res_atomistic:970: for res in res_atomistic:
1130:     string_atomistic += str(res) + ', '971:     string_atomistic += str(res) + ', '
1131: string_local = ''972: string_local = ''
1132: for res in res_local:973: for res in res_local:
1133:     string_local += str(res) + ', '974:     string_local += str(res) + ', '
1134: check_file.write('Atomistic residues: ' + string_atomistic + "\n")975: check_file.write('Atomistic residues: ' + string_atomistic + "\n" )
1135: check_file.write('Locally rigidified residues: ' + string_local + "\n")976: check_file.write('Locally rigidified residues: ' + string_local + "\n" )
1136: 977: 
1137: # GROUPING OF THE RIGID BODIES ######################################################################################978: ###GROUPING OF THE RIGID BODIES #######################################################################################
1138: if res_frozen != []:979: if res_frozen != []:
1139:     for i in res_frozen:980:     for i in res_frozen:
1140:         atoms_used += res_dic[i]981:         atoms_used += res_dic[i]
1141:     print982:     print
1142:     print 'Grouping of the rigid body'983:     print 'Grouping of the rigid body'
1143:     print984:     print
1144:     print '1-Group all parts of the rigid system as one body'985:     print '1-Group all parts of the rigid system as one body'
1145:     print '2-Define groups'986:     print '2-Define groups'
1146:     print987:     print
1147:     while 1:988:     while 1:
1148:         check = raw_input('Grouping method:  ')989:         check = raw_input('Grouping method:  ')
1149:         try:990:         try:
1150:             value = int(check)991:             value = int(check)
1151:             if int(check) == 1:  # all residues are grouped as one rigid body992:             if int(check) == 1:  #all residues are grouped as one rigid body
1152:                 num_rbody += 1993:                 num_rbody += 1
1153:                 check_file.write('All frozen residues are rigidified as one body' + "\n" + "\n" + "\n")994:                 check_file.write('All frozen residues are rigidified as one body' + "\n" + "\n" + "\n")
1154:                 groups = {1: res_frozen}995:                 groups = {1: res_frozen}
1155:                 break996:                 break
1156:             elif int(check) == 2:  # residues are frozen in a number of rigid bodies997:             elif int(check) == 2:  #residues are frozen in a number of rigid bodies
1157:                 check_file.write('The frozen residues are rigidified as more than one body' + "\n" + "\n" + "\n")998:                 check_file.write('The frozen residues are rigidified as more than one body' + "\n" + "\n" + "\n")
1158:                 print 'The following residues are frozen:'999:                 print 'The following residues are frozen:'
1159:                 print res_frozen1000:                 print res_frozen
1160:                 print1001:                 print
1161:                 res_frozen_left = res_frozen1002:                 res_frozen_left = res_frozen
1162:                 num_groups = int(raw_input('How many groups will be defined?  '))  # number of rigid bodies is defined1003:                 num_groups = int(raw_input('How many groups will be defined?  '))  #number of rigid bodies is defined
1163:                 counter = 01004:                 counter = 0
1164:                 while counter < num_groups:1005:                 while counter < num_groups:
1165:                     counter += 11006:                     counter += 1
1166:                     num_rbody += 11007:                     num_rbody += 1
1167:                     group = raw_input('Residues for rigid body: ')1008:                     group = raw_input(
1168:                     # input allows to enter 'r num1 num2' being the range between num1 and num21009:                         'Residues for rigid body: ')  #input allows to enter 'r num1 num2' being the range between num1 and num2
1169:                     group = group.split()1010:                     group = group.split()
1170:                     i = 01011:                     i = 0
1171:                     group_res = []1012:                     group_res = []
1172:                     while i < len(group):1013:                     while i < len(group):
1173:                         if group[i] == 'r':  # if a range is indicated1014:                         if group[i] == 'r':  #if a range is indicated
1174:                             i += 11015:                             i += 1
1175:                             start_range = int(group[i])  # finds starting residue1016:                             start_range = int(group[i])  #finds starting residue
1176:                             i += 11017:                             i += 1
1177:                             end_range = int(group[i])  # finds final residue1018:                             end_range = int(group[i])  #finds final residue
1178:                             for j in range(start_range, end_range):1019:                             for j in range(start_range, end_range):
1179:                                 if j in res_frozen_left:  # appends list if residues have not be assigned earlier1020:                                 if j in res_frozen_left:  #appends list if residues have not be assigned earlier
1180:                                     group_res.append(j)1021:                                     group_res.append(j)
1181:                         else:1022:                         else:
1182:                             try:1023:                             try:
1183:                                 if int(group[i]) in res_frozen_left:1024:                                 if int(group[i]) in res_frozen_left:
1184:                                     group_res.append(int(group[i]))  # appends single residues entered1025:                                     group_res.append(int(group[i]))  #appends single residues entered
1185:                                 i += 11026:                                 i += 1
1186:                             except ValueError:1027:                             except ValueError:
1187:                                 i += 11028:                                 i += 1
1188:                     res_frozen_left = remove_selection_from_list(res_frozen_left,1029:                     res_frozen_left = remove_selection_from_list(res_frozen_left,
1189:                                                                  group_res)1030:                                                                  group_res)  #removes newly assigned residues from rest of residues
1190:                     # removes newly assigned residues from rest of residues1031:                     groups[num_rbody] = group_res  #defines a new entry in the dictionary that saves the groups
1191:                     groups[num_rbody] = group_res  # defines a new entry in the dictionary that saves the groups1032:                     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
1192:                     if res_frozen_left == []: 
1193:                         # automatically terminates when no residues are left to be assigned even if the number of rigid 
1194:                         # bodies entered is not reached 
1195:                         print 'All residues are assigned.'1033:                         print 'All residues are assigned.'
1196:                         break1034:                         break
1197:                     if counter == num_groups and res_frozen_left != []:1035:                     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
1198:                         # if all rigid bodies are assigned and there are still residues not assigned 
1199:                         # it deletes all groups and restarts the process 
1200:                         print 'Not all residues were assigned. Please start again.'1036:                         print 'Not all residues were assigned. Please start again.'
1201:                         counter = 01037:                         counter = 0
1202:                         groups = {}1038:                         groups = {}
1203:                         res_frozen_left = res_frozen1039:                         res_frozen_left = res_frozen
1204:                 break1040:                 break
1205:             else:1041:             else:
1206:                 pass1042:                 pass
1207:         except ValueError:1043:         except ValueError:
1208:             print 'Wrong input'1044:             print 'Wrong input'
1209:     print1045:     print
1210:     print 'The following rigid bodies have been specified:'  # prints all the entered rigid bodies from above1046:     print 'The following rigid bodies have been specified:'  #prints all the entered rigid bodies from above
1211:     for i in groups.keys():1047:     for i in groups.keys():
1212:         print 'Rigid body: ', i1048:         print 'Rigid body: ', i
1213:         print groups[i]1049:         print groups[i]
1214:         print1050:         print
1215: 1051: 
1216: else:1052: else:
1217:     groups = {}1053:     groups = {}
1218: 1054: 
1219: # LOCALLY RIGIDIFICATION, STANDARD GROUPS ###########################################################################1055: ### LOCALLY RIGIDIFICATION, STANDARD GROUPS ###########################################################################
1220: if res_local != []:1056: if res_local != []:
1221:     print 'Grouping of the local rigid bodies'1057:     print 'Grouping of the local rigid bodies'
1222:     print1058:     print
1223:     local_scheme = []1059:     local_scheme = []
1224:     if protein_t: 
1225:         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)? ') 
1226:         if proline_check == '1': 
1227:             check_file.write('Grouped proline as rigid body C-O-N-CD-CG-CB-CA' + "\n") 
1228:             local_scheme.append(1) 
1229:         elif proline_check == '2': 
1230:             check_file.write('Grouped proline as rigid body  N-CD-CG-CB-CA' + "\n") 
1231:             local_scheme.append(2) 
1232:         else: 
1233:             local_scheme.append(0) 
1234:  
1235:         arginine_check = raw_input('Group arginine as rigid body NE-HE-CZ-NH1-HH11-HH12-NH2-HH21-HH22 (y/n)? ') 
1236:         if arginine_check == 'y': 
1237:             check_file.write('Grouped arginine as rigid body NE-HE-CZ-NH1-HH11-HH12-NH2-HH21-HH22' + "\n") 
1238:             local_scheme.append(1) 
1239:         else: 
1240:             local_scheme.append(0) 
1241:  
1242:         histidine_check = raw_input('Group histidine (HIS, HIE, HID, HIP) as rigid body CG-ND1-CE1-NE2-CD2 (y/n)? ') 
1243:         if histidine_check == 'y': 
1244:             check_file.write('Grouped histidine (HIS, HIE, HID, HIP) as rigid body CG-ND1-CE1-NE2-CD2' + "\n") 
1245:             local_scheme.append(1) 
1246:             local_scheme.append(1) 
1247:             local_scheme.append(1) 
1248:             local_scheme.append(1) 
1249:         else: 
1250:             local_scheme.append(0) 
1251:             local_scheme.append(0) 
1252:             local_scheme.append(0) 
1253:             local_scheme.append(0) 
1254:  
1255:         lysine_check = raw_input('Group lysine as rigid body NZ-HZ1-HZ2-HZ3 (y/n)? ') 
1256:         if lysine_check == 'y': 
1257:             check_file.write('Grouped lysine as rigid body NZ-HZ1-HZ2-HZ3' + "\n") 
1258:             local_scheme.append(1) 
1259:         else: 
1260:             local_scheme.append(0) 
1261: 1060: 
1262:         aspartic_check = raw_input('Group aspartic acid as rigid body CB-CG-OD1-OD2 (y/n)? ')1061:     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)? ')
1263:         if aspartic_check == 'y':1062:     if proline_check == '1':
1264:             check_file.write('Grouped aspartic acid as rigid body CB-CG-OD1-OD2' + "\n")1063:         check_file.write('Grouped proline as rigid body C-O-N-CD-CG-CB-CA' + "\n")
1265:             local_scheme.append(1)1064:         local_scheme.append(1)
1266:         else:1065:     elif proline_check == '2':
1267:             local_scheme.append(0)1066:         check_file.write('Grouped proline as rigid body  N-CD-CG-CB-CA' + "\n")
1268: 1067:         local_scheme.append(2)
1269:         asparagine_check = raw_input('Group asparagine as rigid body CB-CG-OD1-ND2-HD21-HD22 (y/n)? ')1068:     else:
1270:         if asparagine_check == 'y':1069:         local_scheme.append(0)
1271:             check_file.write('Grouped asparagine as rigid body GB-CG-OD1-ND2-HD21-HD22' + "\n") 
1272:             local_scheme.append(1) 
1273:         else: 
1274:             local_scheme.append(0) 
1275:  
1276:         glutamic_check = raw_input('Group glutamic acid as rigid body CG-CD-OE1-OE2 (y/n)? ') 
1277:         if glutamic_check == 'y': 
1278:             check_file.write('Grouped glutamic acid as rigid body CG-CD-OE1-OE2' + "\n") 
1279:             local_scheme.append(1) 
1280:         else: 
1281:             local_scheme.append(0) 
1282:  
1283:         glutamine_check = raw_input('Group glutamine as rigid body CG-CD-OE1-NE2-HE21-HE22 (y/n)? ') 
1284:         if glutamine_check == 'y': 
1285:             check_file.write('Grouped glutamine as rigid body CG-CD-OE1-NE2-HE21-HE22' + "\n") 
1286:             local_scheme.append(1) 
1287:         else: 
1288:             local_scheme.append(0) 
1289:  
1290:         phenylalanine_check = raw_input( 
1291:             'Group phenylalanine as rigid body CG-CD1-HD1-CE1-HE1-CZ-HZ-CE2-HE2-CD2-HD2 (y/n)? ') 
1292:         if phenylalanine_check == 'y': 
1293:             check_file.write('Grouped phenylalanine as rigid body CG-CD1-HD1-CE1-HE1-CZ-HZ-CE2-HE2-CD2-HD2' + "\n") 
1294:             local_scheme.append(1) 
1295:         else: 
1296:             local_scheme.append(0) 
1297: 1070: 
1298:         tyrosine_check = raw_input('Group tyrosine as rigid body CG-CD1-HD1-CE1-HE1-CZ-OH-CE2-HE2-CD2-HD2 (y/n)? ')1071:     arginine_check = raw_input('Group arginine as rigid body NE-HE-CZ-NH1-HH11-HH12-NH2-HH21-HH22 (y/n)? ')
1299:         if tyrosine_check == 'y':1072:     if arginine_check == 'y':
1300:             check_file.write('Grouped tyrosine as rigid body CG-CD1-HD1-CE1-HE1-CZ-OH-CE2-HE2-CD2-HD2' + "\n")1073:         check_file.write('Grouped arginine as rigid body NE-HE-CZ-NH1-HH11-HH12-NH2-HH21-HH22' + "\n")
1301:             local_scheme.append(1)1074:         local_scheme.append(1)
1302:         else:1075:     else:
1303:             local_scheme.append(0)1076:         local_scheme.append(0)
1304: 1077: 
1305:         tryptophan_check = raw_input(1078:     histidine_check = raw_input('Group histidine (HIS, HIE, HID, HIP) as rigid body CG-ND1-CE1-NE2-CD2 (y/n)? ')
1306:             'Group tryptophan as either: 1 rigid body:(CG-CD1-HD1-NE1-HE1-CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-CE3-HE3-CD2) '1079:     if histidine_check == 'y':
1307:             'or 2 rigid bodies:(CG-CD1-HD1-NE1-HE1)-(CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-CE3-HE3-CD2)? ')1080:         check_file.write('Grouped histidine (HIS, HIE, HID, HIP) as rigid body CG-ND1-CE1-NE2-CD2' + "\n")
1308:         if tryptophan_check == '1':1081:         local_scheme.append(1)
1309:             check_file.write(1082:         local_scheme.append(1)
1310:                 'Grouped tryptophan as rigid body CG-CD1-HD1-NE1-HE1-CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-HZ3-CE3-HE3-CD2' +1083:         local_scheme.append(1)
1311:                 "\n" + "\n")1084:         local_scheme.append(1)
1312:             local_scheme.append(1)1085:     else:
1313:         elif tryptophan_check == '2':1086:         local_scheme.append(0)
1314:             check_file.write(1087:         local_scheme.append(0)
1315:                 'Grouped tryptophan as two rigid bodies CG-CD1-HD1-NE1-HE1 and' +1088:         local_scheme.append(0)
1316:                 ' CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-HZ3-CE3-HE3-CD2' + "\n" + "\n")1089:         local_scheme.append(0)
1317:             local_scheme.append(2)1090: 
1318:         else:1091:     lysine_check = raw_input('Group lysine as rigid body NZ-HZ1-HZ2-HZ3 (y/n)? ')
1319:             local_scheme.append(0)1092:     if lysine_check == 'y':
 1093:         check_file.write('Grouped lysine as rigid body NZ-HZ1-HZ2-HZ3' + "\n")
 1094:         local_scheme.append(1)
 1095:     else:
 1096:         local_scheme.append(0)
1320: 1097: 
1321:         loc_protein = loaded_mol.get_rigid_for_local(res_local, local_scheme, atoms_used)1098:     aspartic_check = raw_input('Group aspartic acid as rigid body CB-CG-OD1-OD2 (y/n)? ')
1322:         groups_local = loc_protein[0]1099:     if aspartic_check == 'y':
1323:         atoms_used = loc_protein[1]1100:         check_file.write('Grouped aspartic acid as rigid body CB-CG-OD1-OD2' + "\n")
1324:     print1101:         local_scheme.append(1)
 1102:     else:
 1103:         local_scheme.append(0)
1325: 1104: 
1326:     protein_res = ['ALA', 'ARG', 'ASN', 'ASP', 'ASX', 'CYS', 'CYX', 'GLN', 'GLU', 'GLY', 'GLX', 'HIS', 'HIE',1105:     asparagine_check = raw_input('Group asparagine as rigid body CB-CG-OD1-ND2-HD21-HD22 (y/n)? ')
1327:                    'HID', 'HIP', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL',1106:     if asparagine_check == 'y':
1328:                    'NALA', 'NARG', 'NASN', 'NASP', 'NASX', 'NCYS', 'NCYX', 'NGLN', 'NGLU', 'NGLY', 'NGLX',1107:         check_file.write('Grouped asparagine as rigid body GB-CG-OD1-ND2-HD21-HD22' + "\n")
1329:                    'NHIS', 'NHIE', 'NHID', 'NHIP', 'NILE', 'NLEU', 'NLYS', 'NMET', 'NPHE', 'NPRO', 'NSER',1108:         local_scheme.append(1)
1330:                    'NTHR', 'NTRP', 'NTYR', 'NVAL', 'CALA', 'CARG', 'CASN', 'CASP', 'CASX', 'CCYS', 'CCYX',1109:     else:
1331:                    'CGLN', 'CGLU', 'CGLY', 'CGLX', 'CHIS', 'CHIE', 'CHID', 'CHIP', 'CILE', 'CLEU', 'CLYS',1110:         local_scheme.append(0)
1332:                    'CMET', 'CPHE', 'CPRO', 'CSER', 'CTHR', 'CTRP', 'CTYR', 'CVAL'] 
1333:  
1334:     if protein_t: 
1335:         check_peptide_bonds = raw_input('Shall the peptide bonds be treated as rigid bodies (excl. PRO) (y/n)? ') 
1336:         if check_peptide_bonds == 'y' or 'yes': 
1337:             pairlist = []  # all peptide bonds are part of two residues 
1338:             for i in res_local: 
1339:                 if (i - 1, i) not in pairlist and (i - 1 > 0): 
1340:                     if loaded_mol.get_residue_name(i) != 'PRO': 
1341:                         if (loaded_mol.get_residue_name(i) in protein_res) and \ 
1342:                                 (loaded_mol.get_residue_name(i - 1) in protein_res): 
1343:                             pairlist.append((i - 1, i)) 
1344:                 if (i, i + 1) not in pairlist and (i + 1 <= loaded_mol.num_res()): 
1345:                     if loaded_mol.get_residue_name(i + 1) != 'PRO': 
1346:                         if (loaded_mol.get_residue_name(i) in protein_res) and \ 
1347:                                 (loaded_mol.get_residue_name(i + 1) in protein_res): 
1348:                             pairlist.append((i, i + 1)) 
1349:             for j in pairlist: 
1350:                 atomlist = [loaded_mol.get_atom_id_from_name('C', j[0]), loaded_mol.get_atom_id_from_name('O', j[0]), 
1351:                             loaded_mol.get_atom_id_from_name('N', j[1]), loaded_mol.get_atom_id_from_name('H', j[1])] 
1352:                 check_double = 0 
1353:                 for k in atomlist: 
1354:                     if k in atoms_used: 
1355:                         check_double = 1 
1356:                 if check_double == 0: 
1357:                     atoms_used += atomlist 
1358:                     check_file.write( 
1359:                         'Peptide bond rigidified for residues ' + str(j[0]) + ' and ' + str(j[1]) + "\n") 
1360:                     groups_local[len(groups_local) + 1] = [j, (loaded_mol.get_residue_name(j[0]), 
1361:                                                                loaded_mol.get_residue_name(j[1]))] + atomlist 
1362:  
1363:     if rna_t or dna_t: 
1364:         phosphate_check = raw_input("Group phosphate as (O3'-POO-O5') (1) or POO (2)? ") 
1365:         if phosphate_check == '1': 
1366:             NArigid_PO4_large = True 
1367:             NArigid_PO4_small = False 
1368:             check_file.write("Grouped phosphate group between nucleic acids as (O3'-POO-O5')" + "\n") 
1369:         elif phosphate_check == '2': 
1370:             NArigid_PO4_large = False 
1371:             NArigid_PO4_small = True 
1372:             check_file.write("Grouped phosphate group between nucleic acids as (-POO-)" + "\n") 
1373:         else: 
1374:             NArigid_PO4_large = False 
1375:             NArigid_PO4_small = False 
1376: 1111: 
1377:         sugar_check = raw_input("Rigidify sugar (y/n)? ")1112:     glutamic_check = raw_input('Group glutamic acid as rigid body CG-CD-OE1-OE2 (y/n)? ')
1378:         if sugar_check == 'y':1113:     if glutamic_check == 'y':
1379:             NArigid_sugar = True1114:         check_file.write('Grouped glutamic acid as rigid body CG-CD-OE1-OE2' + "\n")
1380:             check_file.write('Grouped sugar group in nucleic acids as rigid body' + "\n")1115:         local_scheme.append(1)
1381:         else:1116:     else:
1382:             NArigid_sugar = False1117:         local_scheme.append(0)
1383: 1118: 
1384:         guanine_check = raw_input('Group guanine base as one rigid body (1) or leave contacts free (2)? ')1119:     glutamine_check = raw_input('Group glutamine as rigid body CG-CD-OE1-NE2-HE21-HE22 (y/n)? ')
1385:         if guanine_check == '1':1120:     if glutamine_check == 'y':
1386:             NArigid_G_1 = True1121:         check_file.write('Grouped glutamine as rigid body CG-CD-OE1-NE2-HE21-HE22' + "\n")
1387:             NArigid_G_2 = False1122:         local_scheme.append(1)
1388:         elif guanine_check == '2':1123:     else:
1389:             NArigid_G_1 = False1124:         local_scheme.append(0)
1390:             NArigid_G_2 = True 
1391:         else: 
1392:             NArigid_G_1 = False 
1393:             NArigid_G_2 = False 
1394: 1125: 
1395:         adenine_check = raw_input('Group adenine base as one rigid body (1) or leave contacts free (2)? ')1126:     phenylalanine_check = raw_input('Group phenylalanine as rigid body CG-CD1-HD1-CE1-HE1-CZ-HZ-CE2-HE2-CD2-HD2 (y/n)? ')
1396:         if adenine_check == '1':1127:     if phenylalanine_check == 'y':
1397:             NArigid_A_1 = True1128:         check_file.write('Grouped phenylalanine as rigid body CG-CD1-HD1-CE1-HE1-CZ-HZ-CE2-HE2-CD2-HD2' + "\n")
1398:             NArigid_A_2 = False1129:         local_scheme.append(1)
1399:         elif adenine_check == '2':1130:     else:
1400:             NArigid_A_1 = False1131:         local_scheme.append(0)
1401:             NArigid_A_2 = True 
1402:         else: 
1403:             NArigid_A_1 = False 
1404:             NArigid_A_2 = False 
1405: 1132: 
1406:         cystosine_check = raw_input('Group cystosine base as one rigid body (1) or leave contacts free (2)? ')1133:     tyrosine_check = raw_input('Group tyrosine as rigid body CG-CD1-HD1-CE1-HE1-CZ-OH-CE2-HE2-CD2-HD2 (y/n)? ')
1407:         if cystosine_check == '1':1134:     if tyrosine_check == 'y':
1408:             NArigid_C_1 = True1135:         check_file.write('Grouped tyrosine as rigid body CG-CD1-HD1-CE1-HE1-CZ-OH-CE2-HE2-CD2-HD2' + "\n")
1409:             NArigid_C_2 = False1136:         local_scheme.append(1)
1410:         elif cystosine_check == '2':1137:     else:
1411:             NArigid_C_1 = False1138:         local_scheme.append(0)
1412:             NArigid_C_2 = True 
1413:         else: 
1414:             NArigid_C_1 = False 
1415:             NArigid_C_2 = False 
1416: 1139: 
1417:         NArigid_T_1 = False1140:     tryptophan_check = raw_input(
1418:         NArigid_T_2 = False1141:         'Group tryptophan as either: 1 rigid body:(CG-CD1-HD1-NE1-HE1-CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-CE3-HE3-CD2) '
1419:         NArigid_U_1 = False1142:         'or 2 rigid bodies:(CG-CD1-HD1-NE1-HE1)-(CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-CE3-HE3-CD2)? ')
1420:         NArigid_U_2 = False1143:     if tryptophan_check == '1':
1421:         if dna_t:1144:         check_file.write(
1422:             thymine_check = raw_input('Group thymine base as one rigid body (1) or leave contacts free (2)? ')1145:             'Grouped tryptophan as rigid body CG-CD1-HD1-NE1-HE1-CE2-CZ2-HZ2-CH2-HH2-CZ3-HZ3-HZ3-CE3-HE3-CD2' + "\n" + "\n")
1423:             if thymine_check == '1':1146:         local_scheme.append(1)
1424:                 NArigid_T_1 = True1147:     elif tryptophan_check == '2':
1425:                 NArigid_T_2 = False1148:         check_file.write(
1426:             elif thymine_check == '2':1149:             '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")
1427:                 NArigid_T_1 = False1150:         local_scheme.append(2)
1428:                 NArigid_T_2 = True1151:     else:
1429:             else:1152:         local_scheme.append(0)
1430:                 NArigid_T_1 = False 
1431:                 NArigid_T_2 = False 
1432: 1153: 
1433:         if rna_t:1154:     loc = loaded_mol.get_rigid_for_local(res_local, local_scheme, atoms_used)
1434:             uracil_check = raw_input('Group uracil base as one rigid body (1) or leave contacts free (2)? ')1155:     groups_local = loc[0]
1435:             if uracil_check == '1':1156:     atoms_used = loc[1]
1436:                 NArigid_U_1 = True1157:     print
1437:                 NArigid_U_2 = False 
1438:             elif uracil_check == '2': 
1439:                 NArigid_U_1 = False 
1440:                 NArigid_U_2 = True 
1441:             else: 
1442:                 NArigid_U_1 = False 
1443:                 NArigid_U_2 = False 
1444: 1158: 
1445:         gloc_NA, atoms_used = loaded_mol.loc_NA(res_local, 
1446:                                                 [NArigid_G_1, NArigid_G_2, NArigid_A_1, NArigid_A_2, NArigid_C_1, 
1447:                                                  NArigid_C_2, NArigid_T_1, NArigid_T_2, NArigid_U_1, NArigid_U_2], 
1448:                                                 atoms_used) 
1449:  
1450:         print gloc_NA 
1451:  
1452:         for key in gloc_NA.keys(): 
1453:             groups_local[len(groups_local) + 1] = gloc_NA[key]  # add NA groups to overall groups 
1454:  
1455:         # now do all sugars and phosphate if necessary 
1456:  
1457:         NA_res = ['RAN', 'RCN', 'RGN', 'RUN', 'A', 'A3', 'A5', 'AN', 'C', 'C3', 'C5', 'CN', 'G', 'G3', 'G5', 'GN', 'U', 
1458:                   'U3', 'U5', 'UN', 'OHE', 'RA', 'RA3', 'RA5', 'RC', 'RC3', 'RC5', 'RG', 'RG3', 'RG5', 'RU', 'RU3', 
1459:                   'RU5', 'DA', 'DA5', 'DC', 'DC5', 'DG', 'DG5', 'DT', 'DT5', 'DA3', 'DC3', 'DG3', 'DT3', 'DAN', 'DCN', 
1460:                   'DGN', 'DTN'] 
1461:  
1462:         if NArigid_PO4_large: 
1463:             pairlist = []  # all phosphates that are part of two NAs 
1464:             for i in res_local: 
1465:                 if loaded_mol.get_residue_name(i) in NA_res: 
1466:                     if (i - 1, i) not in pairlist and (i - 1 > 0): 
1467:                         pairlist.append((i - 1, i)) 
1468:                     if (i, i + 1) not in pairlist and (i + 1 <= loaded_mol.num_res()): 
1469:                         pairlist.append((i, i + 1)) 
1470:             for j in pairlist: 
1471:                 atomlist = [loaded_mol.get_atom_id_from_name("O3'", j[0]), loaded_mol.get_atom_id_from_name('P', j[1]), 
1472:                             loaded_mol.get_atom_id_from_name('OP1', j[1]), 
1473:                             loaded_mol.get_atom_id_from_name('OP2', j[1]), 
1474:                             loaded_mol.get_atom_id_from_name("O5'", j[1])] 
1475: 1159: 
 1160:     check_peptide_bonds = raw_input('Shall the peptide bonds be treated as rigid bodies (excl. PRO) (y/n)? ')
 1161:     if check_peptide_bonds == 'y' or 'yes':
 1162:         pairlist = []  #all peptide bonds are part of two residues
 1163:         for i in res_local:
 1164:             if (i - 1, i) not in pairlist and (i - 1 > 0):
 1165:                 if loaded_mol.get_residue_name(i) != 'PRO':
 1166:                     pairlist.append((i - 1, i))
 1167:             if (i, i + 1) not in pairlist and (i + 1 <= loaded_mol.num_res()):
 1168:                 if loaded_mol.get_residue_name(i + 1) != 'PRO':
 1169:                     pairlist.append((i, i + 1))
 1170:         for j in pairlist:
 1171:             atomlist = []
 1172:             atomlist.append(loaded_mol.get_atom_id_from_name('C', j[0]))
 1173:             atomlist.append(loaded_mol.get_atom_id_from_name('O', j[0]))
 1174:             atomlist.append(loaded_mol.get_atom_id_from_name('N', j[1]))
 1175:             atomlist.append(loaded_mol.get_atom_id_from_name('H', j[1]))
 1176:             check_double = 0
 1177:             for k in atomlist:
 1178:                 if k in atoms_used:
 1179:                     check_double = 1
 1180:             if check_double == 0:
1476:                 atoms_used += atomlist1181:                 atoms_used += atomlist
1477:                 check_file.write(1182:                 check_file.write(
1478:                     'Phosphate (all) rigidified for residues ' + str(j[0]) + ' and ' + str(j[1]) + "\n")1183:                     'Peptide bond rigidified for residues ' + str(j[0]) + ' and ' + str(j[1]) + "\n")
1479:                 groups_local[len(groups_local) + 1] = [j, (loaded_mol.get_residue_name(j[0]),1184:                 groups_local[len(groups_local) + 1] = [j, (loaded_mol.get_residue_name(j[0]),
1480:                                                            loaded_mol.get_residue_name(j[1]))] + atomlist1185:                                                            loaded_mol.get_residue_name(j[1]))] + atomlist
1481:         if NArigid_PO4_small: 
1482:             for i in res_local: 
1483:                 if loaded_mol.get_residue_name(i) in NA_res: 
1484:                     atomlist = [loaded_mol.get_atom_id_from_name('P', i), loaded_mol.get_atom_id_from_name('OP1', i), 
1485:                                 loaded_mol.get_atom_id_from_name('OP2', i)] 
1486:                     atoms_used += atomlist 
1487:                     check_file.write( 
1488:                         'Phosphate (only POO) rigidified for residue ' + str(i) + "\n") 
1489:                     groups_local[len(groups_local) + 1] = [i, loaded_mol.get_residue_name(i)] + atomlist 
1490:  
1491:         if NArigid_sugar: 
1492:             for i in res_local: 
1493:                 if loaded_mol.get_residue_name(i) in NA_res: 
1494:                     atomlist = [loaded_mol.get_atom_id_from_name("C5'", i), loaded_mol.get_atom_id_from_name("H5'", i), 
1495:                                 loaded_mol.get_atom_id_from_name("H5''", i), loaded_mol.get_atom_id_from_name("C4'", i), 
1496:                                 loaded_mol.get_atom_id_from_name("H4'", i), loaded_mol.get_atom_id_from_name("O4'", i), 
1497:                                 loaded_mol.get_atom_id_from_name("C3'", i), loaded_mol.get_atom_id_from_name("H3'", i), 
1498:                                 loaded_mol.get_atom_id_from_name("C2'", i), loaded_mol.get_atom_id_from_name("H2'", i), 
1499:                                 loaded_mol.get_atom_id_from_name("C1'", i), loaded_mol.get_atom_id_from_name("H1'", i)] 
1500:                     if loaded_mol.get_residue_name(i) in ['DA', 'DA5', 'DC', 'DC5', 'DG', 'DG5', 'DT', 'DT5', 'DA3', 
1501:                                                           'DC3', 'DG3', 'DT3', 'DAN', 'DCN','DGN', 'DTN']: 
1502:                         atomlist.append(loaded_mol.get_atom_id_from_name("H2''", i)) 
1503:                     else: 
1504:                         atomlist.append(loaded_mol.get_atom_id_from_name("O2'", i)) 
1505:                         atomlist.append(loaded_mol.get_atom_id_from_name("HO2'", i)) 
1506:                      
1507:                     atoms_used += atomlist 
1508:                     check_file.write( 
1509:                         'Sugar rigidified for residue ' + str(i) + "\n") 
1510:                     groups_local[len(groups_local) + 1] = [i, loaded_mol.get_residue_name(i)] + atomlist 
1511:  
1512:     check_file.write("\n")1186:     check_file.write("\n")
1513:     check_new_groups = raw_input('Enter additional user-defined rigid bodies (y/n)? ')1187:     check_new_groups = raw_input('Enter additional user-defined rigid bodies (y/n)? ')
1514:     if check_new_groups in ['yes', 'y']:1188:     if check_new_groups in ['yes', 'y']:
1515:         print '1 - select a residues by name'1189:         print '1 - select a residues by name'
1516:         print '2 - select one or more residues by id'1190:         print '2 - select one or more residues by id'
1517:         while 1:1191:         while 1:
1518:             while 1:1192:             while 1:
1519:                 local_method = raw_input('Method (1/2): ')1193:                 local_method = raw_input('Method (1/2): ')
1520:                 try:1194:                 try:
1521:                     if int(local_method) in [1, 2]:1195:                     if int(local_method) in [1, 2]:
1528:             if int(local_method) == 1:1202:             if int(local_method) == 1:
1529:                 while 1:1203:                 while 1:
1530:                     exit_1 = 01204:                     exit_1 = 0
1531:                     res_name = raw_input('Residue name (three letter code, n to exit): ')1205:                     res_name = raw_input('Residue name (three letter code, n to exit): ')
1532:                     if res_name == 'n':1206:                     if res_name == 'n':
1533:                         exit_1 = 11207:                         exit_1 = 1
1534:                         break1208:                         break
1535:                     try:1209:                     try:
1536:                         res_list = loaded_mol.get_all_res_name(res_local, res_name)1210:                         res_list = loaded_mol.get_all_res_name(res_local, res_name)
1537:                         if res_list != []:1211:                         if res_list != []:
1538:                             check_file.write('Additional rigid bodies for residues: ' + res_name + "\n")1212:                             check_file.write('Additional rigid bodies for residues: ' + res_name +  "\n")
1539:                             break1213:                             break
1540:                         else:1214:                         else:
1541:                             print 'None of the residues for local rigidification is of that type.'1215:                             print 'None of the residues for local rigidification is of that type.'
1542:                     except ValueError:1216:                     except ValueError:
1543:                         print 'Invalid input'1217:                         print 'Invalid input'
1544:                 if exit_1 == 0:1218:                 if exit_1 == 0:
1545:                     all_atoms_res_name = []1219:                     all_atoms_res_name = []
1546:                     atoms_used_name = []1220:                     atoms_used_name = []
1547:                     for i in res_list:1221:                     for i in res_list:
1548:                         all_atoms_res_name += res_dic[i]1222:                         all_atoms_res_name += res_dic[i]
1571:                             print 'Not all atoms entered are in the residue.'1245:                             print 'Not all atoms entered are in the residue.'
1572:                             check_file.write('For residue: ' + str(m) + ' - not all entered atoms in residue' + "\n")1246:                             check_file.write('For residue: ' + str(m) + ' - not all entered atoms in residue' + "\n")
1573:                             check_res = 11247:                             check_res = 1
1574:                         if len(atoms_rbody) < 3:1248:                         if len(atoms_rbody) < 3:
1575:                             print 'Rigid bodies need to have at least 3 atoms.'1249:                             print 'Rigid bodies need to have at least 3 atoms.'
1576:                             check_file.write('For residue: ' + str(m) + ' - not 3 or more atoms in residue' + "\n")1250:                             check_file.write('For residue: ' + str(m) + ' - not 3 or more atoms in residue' + "\n")
1577:                             check_res = 11251:                             check_res = 1
1578:                         for o in atoms_rbody:1252:                         for o in atoms_rbody:
1579:                             if o in atoms_used_name:1253:                             if o in atoms_used_name:
1580:                                 print 'Atoms cannot belong to two rigid bodies.'1254:                                 print 'Atoms cannot belong to two rigid bodies.'
1581:                                 check_file.write('For residue: ' + str(m) +1255:                                 check_file.write('For residue: ' + str(m) + ' - atoms cannot belong to two rigid bodies'
1582:                                                  ' - atoms cannot belong to two rigid bodies' + "\n" + "\n")1256:                                                  + "\n" + "\n")
1583:                                 check_res = 11257:                                 check_res = 1
1584:                                 break1258:                                 break
1585:                         if loaded_mol.check_linear(atoms_rbody, tolerance):1259:                         if loaded_mol.check_linear(atoms_rbody, tolerance):
1586:                             print 'Linear arrangement of atoms cannot be a rigid body.'1260:                             print 'Linear arrangement of atoms cannot be a rigid body.'
1587:                             check_file.write('For residue: ' + str(m) + ' - linear group entered' + "\n")1261:                             check_file.write('For residue: ' + str(m) + ' - linear group entered' + "\n")
1588:                             check_res = 11262:                             check_res = 1
1589:                         if check_res == 0:1263:                         if check_res == 0:
1590:                             groups_local[len(groups_local) + 1] = [m, loaded_mol.get_residue_name(m)] + atoms_rbody1264:                             groups_local[len(groups_local) + 1] = [m, loaded_mol.get_residue_name(m)] + atoms_rbody
1591:                             atoms_used += atoms_rbody1265:                             atoms_used += atoms_rbody
1592:                             check_file.write('For residue: ' + str(m) + ' new group: ' + atom_seq + "\n")1266:                             check_file.write('For residue: ' + str(m) + ' new group: ' + atom_seq + "\n")
1597:                     res_inp = raw_input('Residue numbers (n to exit): ')1271:                     res_inp = raw_input('Residue numbers (n to exit): ')
1598:                     if res_inp == 'n':1272:                     if res_inp == 'n':
1599:                         exit_2 = 11273:                         exit_2 = 1
1600:                         break1274:                         break
1601:                     try:1275:                     try:
1602:                         res_list_str = res_inp.split()1276:                         res_list_str = res_inp.split()
1603:                         res_list_int = []1277:                         res_list_int = []
1604:                         for res in res_list_str:1278:                         for res in res_list_str:
1605:                             res_list_int.append(int(res))1279:                             res_list_int.append(int(res))
1606:                         if res_list_int != []:1280:                         if res_list_int != []:
1607:                             check_file.write('Additional rigid bodies for residues: ' + res_inp + "\n")1281:                             check_file.write('Additional rigid bodies for residues: ' + res_inp  + "\n")
1608:                             break1282:                             break
1609:                         else:1283:                         else:
1610:                             print 'No residues entered'1284:                             print 'No residues entered'
1611:                     except ValueError:1285:                     except ValueError:
1612:                         print 'Invalid input'1286:                         print 'Invalid input'
1613:                 if exit_2 == 0:1287:                 if exit_2 == 0:
1614:                     all_atoms_res = []1288:                     all_atoms_res = []
1615:                     atoms_used_res = []1289:                     atoms_used_res = []
1616:                     for i in res_list_int:1290:                     for i in res_list_int:
1617:                         all_atoms_res += res_dic[i]1291:                         all_atoms_res += res_dic[i]
1620:                             atoms_used_res.append(j)1294:                             atoms_used_res.append(j)
1621:                     atoms_allowed = []1295:                     atoms_allowed = []
1622:                     for atom in all_atoms_res:1296:                     for atom in all_atoms_res:
1623:                         if atom not in atoms_used_res:1297:                         if atom not in atoms_used_res:
1624:                             atoms_allowed.append(atom)1298:                             atoms_allowed.append(atom)
1625:                     if atoms_used_res != []:1299:                     if atoms_used_res != []:
1626:                         print 'The following atoms are already in local rigid bodies: '1300:                         print 'The following atoms are already in local rigid bodies: '
1627:                         for k in atoms_used_res:1301:                         for k in atoms_used_res:
1628:                             print str(k) + ' ' + atom_dic[k][0] + ' ' + atom_dic[k][1] + ' ' + str(atom_dic[k][2])1302:                             print str(k) + ' ' + atom_dic[k][0] + ' ' + atom_dic[k][1] + ' ' + str(atom_dic[k][2])
1629:                     print1303:                     print
1630:                     print 'The residues %s contain the following atoms that are not assigned to rigid bodies:' \1304:                     print 'The residues %s contain the following atoms that are not assigned to rigid bodies:' % res_list_str
1631:                           % res_list_str 
1632:                     string_atoms = ''1305:                     string_atoms = ''
1633:                     for atom in atoms_allowed:1306:                     for atom in atoms_allowed:
1634:                         string_atoms += str(atom) + ' - ' + atom_dic[atom][0] + ', '1307:                         string_atoms += str(atom) + ' - ' + atom_dic[atom][0] + ', '
1635:                     print string_atoms1308:                     print string_atoms
1636:                     print1309:                     print
1637:                     atom_seq = raw_input('Enter a sequence of atom numbers separated by spaces: ')1310:                     atom_seq = raw_input('Enter a sequence of atom numbers separated by spaces: ')
1638:                     check_res = 01311:                     check_res = 0
1639:                     atoms_rbody = []1312:                     atoms_rbody = []
1640:                     check_rbody = 01313:                     check_rbody = 0
1641:                     for atom_str in atom_seq.split():1314:                     for atom_str in atom_seq.split():
1642:                         try:1315:                         try:
1643:                             if int(atom_str) in atoms_allowed:1316:                             if int(atom_str) in atoms_allowed:
1644:                                 atoms_rbody.append(int(atom_str))1317:                                 atoms_rbody.append(int(atom_str))
1645:                             elif int(atom_str) in atoms_used_res:1318:                             elif int(atom_str) in atoms_used_res:
1646:                                 print 'Atoms cannot belong to two rigid bodies.'1319:                                 print 'Atoms cannot belong to two rigid bodies.'
1647:                                 check_file.write('Atoms cannot belong to two rigid bodies' + "\n")1320:                                 check_file.write('Atoms cannot belong to two rigid bodies'  + "\n")
1648:                                 check_rbody = 11321:                                 check_rbody = 1
1649:                                 break1322:                                 break
1650:                         except ValueError:1323:                         except ValueError:
1651:                             pass1324:                             pass
1652:                     print len(atoms_rbody)1325:                     print len(atoms_rbody)
1653:                     print len(atom_seq.split())1326:                     print len(atom_seq.split())
1654:                     if len(atoms_rbody) != len(atom_seq.split()):1327:                     if len(atoms_rbody) != len(atom_seq.split()):
1655:                         print 'Not all atoms entered are in the residue.'1328:                         print 'Not all atoms entered are in the residue.'
1656:                         check_file.write('Not all entered atoms in residue' + "\n")1329:                         check_file.write('Not all entered atoms in residue' + "\n")
1657:                         check_rbody = 11330:                         check_rbody = 1
1665:                         check_rbody = 11338:                         check_rbody = 1
1666:                     if check_rbody == 0:1339:                     if check_rbody == 0:
1667:                         groups_local[len(groups_local) + 1] = ['', ''] + atoms_rbody1340:                         groups_local[len(groups_local) + 1] = ['', ''] + atoms_rbody
1668:                         atoms_used += atoms_rbody1341:                         atoms_used += atoms_rbody
1669:                         check_file.write('For residues: ' + res_inp + ' new group: ' + atom_seq + "\n")1342:                         check_file.write('For residues: ' + res_inp + ' new group: ' + atom_seq + "\n")
1670:                     check_file.write("\n")1343:                     check_file.write("\n")
1671:             check_more = raw_input('Enter more rigid bodies (y/n)? ')1344:             check_more = raw_input('Enter more rigid bodies (y/n)? ')
1672:             if check_more not in ['y', 'yes']:1345:             if check_more not in ['y', 'yes']:
1673:                 break1346:                 break
1674: 1347: 
1675: # WRITING OUTPUT FILES ##############################################################################################1348: ### WRITING OUTPUT FILES ##############################################################################################
1676: 1349: 
1677: print 'Output files will be written now.'1350: print 'Output files will be written now.'
1678: 1351: 
1679: # COORDSINIRIGID #1352: ### COORDSINIRIGID ###
1680: coords_out = open('coordsinirigid', 'w')1353: coords_out = open('coordsinirigid', 'w')
1681:  
1682: for i in range(1, len(atom_dic.keys()) + 1):1354: for i in range(1, len(atom_dic.keys()) + 1):
1683:     atom = str(atom_dic[i][3]) + '   ' + str(atom_dic[i][4]) + '   ' + str(atom_dic[i][5])1355:     atom = str(atom_dic[i][3]) + '   ' + str(atom_dic[i][4]) + '   ' + str(atom_dic[i][5])
1684:     coords_out.write(atom)1356:     coords_out.write(atom)
1685:     if i < len(atom_dic.keys()):1357:     if i < len(atom_dic.keys()):
1686:         coords_out.write("\n")1358:         coords_out.write("\n")
1687: coords_out.close()1359: coords_out.close()
1688: 1360: 
1689: # RBODYCONFIG FILE #1361: ### RBODYCONFIG FILE ###
1690: groups_file = open('rbodyconfig', 'w')1362: groups_file = open('rbodyconfig', 'w')
1691: natom_rbody = 01363: natom_rbody = 0
1692: for group in range(1, len(groups.keys()) + 1):1364: for group in range(1, len(groups.keys()) + 1):
1693:     string = ''1365:     string = ''
1694:     counter = 01366:     counter = 0
1695:     for res in groups[group]:1367:     for res in groups[group]:
1696:         for atom in res_dic[res]:1368:         for atom in res_dic[res]:
1697:             string += str(atom) + "\n"1369:             string += str(atom) + "\n"
1698:             counter += 11370:             counter += 1
1699:             natom_rbody += 11371:             natom_rbody += 1
1701: for group_local in range(1, len(groups_local.keys()) + 1):1373: for group_local in range(1, len(groups_local.keys()) + 1):
1702:     string = ''1374:     string = ''
1703:     counter = 01375:     counter = 0
1704:     for atom in groups_local[group_local][2:]:1376:     for atom in groups_local[group_local][2:]:
1705:         string += str(atom) + "\n"1377:         string += str(atom) + "\n"
1706:         counter += 11378:         counter += 1
1707:         natom_rbody += 11379:         natom_rbody += 1
1708:     groups_file.write('GROUP ' + str(counter) + "\n" + string)1380:     groups_file.write('GROUP ' + str(counter) + "\n" + string)
1709: groups_file.close()1381: groups_file.close()
1710: 1382: 
1711: # FINISH LOG FILE FOR RIGID BODIES #1383: ### FINISH LOG FILE FOR RIGID BODIES###
1712: check_file.write("\n" + "\n" + 'Number of rigid bodies: ' + str(len(groups.keys()) + len(groups_local.keys())))1384: check_file.write("\n" + "\n" + 'Number of rigid bodies: ' + str(len(groups.keys()) + len(groups_local.keys())))
1713: check_file.write(1385: check_file.write("\n" + 'Degrees of freedom: ' + str(3 * natom_rbody + 6* (len(groups.keys()) + len(groups_local.keys()))))
1714:     "\n" + 'Degrees of freedom: ' + str(3 * natom_rbody + 6 * (len(groups.keys()) + len(groups_local.keys()))))1386: 
 1387: 
1715: 1388: 
1716: # CREATING GROUP ROTATION FILES #####################################################################################1389: ### CREATING GROUP ROTATION FILES #####################################################################################
1717: print ''1390: print ''
1718: group_rot_t = raw_input('Create group rotation file (y/n)? ')1391: group_rot_t = raw_input('Create group rotation file (y/n)? ')
1719: if group_rot_t in ['y', 'yes']:1392: if group_rot_t in ['y', 'yes']:
1720:     check_file.write("\n" + "\n" + 'Create group rotation file' + "\n")1393:     check_file.write("\n" + "\n" + 'Create group rotation file' + "\n")
1721:     while 1:1394:     while 1:
1722:         prob_uni_inp = raw_input('Uniform selection probability (suggested: 0.025): ')1395:         prob_uni_inp = raw_input('Uniform selection probability (suggested: 0.025): ')
1723:         try:1396:         try:
1724:             prob_uni = float(prob_uni_inp)1397:             prob_uni = float(prob_uni_inp)
1725:             if prob_uni > 0.1:1398:             if prob_uni > 0.1:
1726:                 print 'This is a large probability (values recommended: 0.01 to 0.1).'1399:                 print 'This is a large probability (values recommended: 0.01 to 0.1).'
1727:                 prob_uni_check = raw_input(' Continue with the selected value (y/n)? ')1400:                 prob_uni_check = raw_input(' Continue with the selected value (y/n)? ')
1728:                 if prob_uni_check in ['y', 'yes']:1401:                 if prob_uni_check in ['y','yes']:
1729:                     break1402:                     break
1730:             elif prob_uni < 0.01:1403:             elif prob_uni < 0.01:
1731:                 print 'This is a small probability (values recommended: 0.01 to 0.1).'1404:                 print 'This is a small probability (values recommended: 0.01 to 0.1).'
1732:                 prob_uni_check = raw_input(' Continue with the selected value (y/n)? ')1405:                 prob_uni_check = raw_input(' Continue with the selected value (y/n)? ')
1733:                 if prob_uni_check in ['y', 'yes']:1406:                 if prob_uni_check in ['y','yes']:
1734:                     break1407:                     break
1735:             else:1408:             else:
1736:                 break1409:                 break
1737:         except ValueError:1410:         except ValueError:
1738:             print 'Input needs to be a float'1411:             print 'Input needs to be a float'
1739:     check_file.write('Selection probability: ' + str(prob_uni) + "\n")1412:     check_file.write('Selection probability: ' + str(prob_uni) + "\n")
1740:     while 1:1413:     while 1:
1741:         grouprot_reg_inp = raw_input(1414:         grouprot_reg_inp = raw_input('Create the group rotation file for the all atom region (1), the locally rigidified region (2), or both regions (3)? ')
1742:             'Create the group rotation file for the all atom region (1), the locally rigidified region (2),' + 
1743:             ' or both regions (3)? ') 
1744:         try:1415:         try:
1745:             group_reg = int(grouprot_reg_inp)1416:             group_reg = int(grouprot_reg_inp)
1746:             if group_reg == 1:1417:             if group_reg == 1:
1747:                 grot_res = res_atomistic1418:                 grot_res = res_atomistic
1748:                 check_file.write('Only all atom region' + "\n")1419:                 check_file.write('Only all atom region' + "\n")
1749:                 break1420:                 break
1750:             elif group_reg == 2:1421:             elif group_reg == 2:
1751:                 grot_res = res_local1422:                 grot_res = res_local
1752:                 check_file.write('Only local rigid region' + "\n")1423:                 check_file.write('Only local rigid region' + "\n")
1753:                 break1424:                 break
1754:             elif group_reg == 3:1425:             elif group_reg == 3:
1755:                 grot_res = res_atomistic + res_local1426:                 grot_res = res_atomistic + res_local
1756:                 check_file.write('Local rigid and all atom regions' + "\n")1427:                 check_file.write('Local rigid and all atom regions' + "\n")
1757:                 break1428:                 break
1758:             else:1429:             else:
1759:                 print 'Please choose one of the options.'1430:                 print 'Please choose one of the options.'
1760:         except ValueError:1431:         except ValueError:
1761:             print 'Please choose one of the options.'1432:             print 'Please choose one of the options.'
1762: 1433: 
1763:     group_rotations = [False, False, False, False, False]1434:     group_rotations = [False , False , False , False , False]
1764:     def_values = [True, True, True, True, True]1435:     def_values = [True , True , True , True ,True]
1765:     # default amplitudes (copied from original Fortran script)1436:     #default amplitudes (copied from original Fortran script)
1766:     amp_dic = {'ALA': (1.0,),1437:     amp_dic = {'ALA' : (1.0 ,),
1767:                'ASN': (0.5, 1.0),1438:                'ASN' : (0.5 , 1.0),
1768:                'ARG': (0.2, 0.3, 0.5),1439:                'ARG' : (0.2 , 0.3 , 0.5),
1769:                'ASP': (0.5, 1.0, 0.0),1440:                'ASP' : (0.5 , 1.0 , 0.0),
1770:                'CYS': (1.0,),1441:                'CYS' : (1.0 ,),
1771:                'GLN': (0.3, 0.5, 1.0),1442:                'GLN' : (0.3 , 0.5 , 1.0),
1772:                'GLU': (0.3, 0.5, 1.0),1443:                'GLU' : (0.3 , 0.5 , 1.0),
1773:                'HIS': (0.3, 0.5),1444:                'HIS' : (0.3 , 0.5),
1774:                'HIE': (0.3, 0.5),1445:                'HIE' : (0.3 , 0.5),
1775:                'HIP': (0.3, 0.5),1446:                'HIP' : (0.3 , 0.5),
1776:                'HID': (0.3, 0.5),1447:                'HIS' : (0.3 , 0.5),
1777:                'ILE': (0.5, 1.0),1448:                'HID' : (0.3 , 0.5),
1778:                'LEU': (0.5, 1.0),1449:                'ILE' : (0.5 , 1.0),
1779:                'LYS': (0.2, 0.3, 0.5),1450:                'LEU' : (0.5 , 1.0),
1780:                'MET': (0.5, 0.7),1451:                'LYS' : (0.2 , 0.3 , 0.5),
1781:                'PHE': (0.3, 0.5),1452:                'MET' : (0.5 , 0.7),
1782:                'SER': (1.0,),1453:                'PHE' : (0.3 , 0.5),
1783:                'THR': (1.0,),1454:                'SER' : (1.0 ,),
1784:                'TRP': (0.3, 0.4),1455:                'THR' : (1.0 ,),
1785:                'TYR': (0.3, 0.5),1456:                'TRP' : (0.3 , 0.4 ),
1786:                'VAL': (1.0,),1457:                'TYR' : (0.3 , 0.5 ),
1787:                'AMB0': (0.1,)  # default for all non residue specific rotations1458:                'VAL' : (1.0 ,),
 1459:                'AMB0': (0.1 ,)     #default for all non residue specific rotations
1788:                }1460:                }
1789:     group_rotations_NA = [False, False, False, False] 
1790:     def_values_NA = [True, True, True, True] 
1791:     # default amplitudes (copied from original Fortran script) 
1792:     amp_dic_NA = {'A': (0.2, 0.3, 0.3, 0.1), 
1793:                   'C': (0.2, 0.3, 0.3, 0.1), 
1794:                   'G': (0.2, 0.3, 0.3, 0.1), 
1795:                   'U': (0.2, 0.3, 0.3, 0.1), 
1796:                   'T': (0.2, 0.3, 0.3, 0.1) 
1797:                   } 
1798: 1461: 
1799:     # number: [name, ax_atom1, ax_atom2, natom, amp, prob, [atom1, atom2, ...]]1462:     #number: [name, ax_atom1, ax_atom2, natom, amp, prob, [atom1, atom2, ...]]
1800:     grot_dic = {}1463:     grot_dic = {}
1801: 1464: 
1802:     if protein_t:1465:     N_CA = raw_input('Group rotations of side chains about N-CA axis (not PRO) (y/n)? ')
1803:         N_CA = raw_input('Group rotations of side chains about N-CA axis (not PRO) (y/n)? ')1466:     if N_CA in ['y' , 'yes']:
1804:         if N_CA in ['y', 'yes']:1467:         group_rotations[0] = True
1805:             group_rotations[0] = True1468:         check_file.write('Group rotation for N - CA' + "\n")
1806:             check_file.write('Group rotation for N - CA' + "\n")1469:         N_CA_def = raw_input('Use default values (y/n)? ')
1807:             N_CA_def = raw_input('Use default values (y/n)? ')1470:         if N_CA_def not in ['y' , 'yes']:
1808:             if N_CA_def not in ['y', 'yes']:1471:             def_values[0] = False
1809:                 def_values[0] = False1472:     C_CA = raw_input('Group rotations of side chains about C-CA axis (not PRO) (y/n)? ')
1810:         C_CA = raw_input('Group rotations of side chains about C-CA axis (not PRO) (y/n)? ')1473:     if C_CA in ['y' , 'yes']:
1811:         if C_CA in ['y', 'yes']:1474:         check_file.write('Group rotation for C - CA' + "\n")
1812:             check_file.write('Group rotation for C - CA' + "\n")1475:         group_rotations[1] = True
1813:             group_rotations[1] = True1476:         C_CA_def = raw_input('Use default values (y/n)? ')
1814:             C_CA_def = raw_input('Use default values (y/n)? ')1477:         if C_CA_def not in ['y' , 'yes']:
1815:             if C_CA_def not in ['y', 'yes']:1478:             def_values[1] = False
1816:                 def_values[1] = False1479:     CA_CB = raw_input('Group rotations of side chains about CA-CB axis (not ALA, GLY, PRO) (y/n)? ')
1817:         CA_CB = raw_input('Group rotations of side chains about CA-CB axis (not ALA, GLY, PRO) (y/n)? ')1480:     if CA_CB in ['y' , 'yes']:
1818:         if CA_CB in ['y', 'yes']:1481:         check_file.write('Group rotation for CA - CB' + "\n")
1819:             check_file.write('Group rotation for CA - CB' + "\n")1482:         group_rotations[2] = True
1820:             group_rotations[2] = True1483:         CA_CB_def = raw_input('Use default values (y/n)? ')
1821:             CA_CB_def = raw_input('Use default values (y/n)? ')1484:         if CA_CB_def not in ['y' , 'yes']:
1822:             if CA_CB_def not in ['y', 'yes']:1485:             def_values[2] = False
1823:                 def_values[2] = False1486:     CB_CG = raw_input('Group rotations of side chains about CB-CG axis (not ALA, GLY, PRO, SER, CYS, THR, VAL) (y/n)? ')
1824:         CB_CG = raw_input(1487:     if CB_CG in ['y' , 'yes']:
1825:             'Group rotations of side chains about CB-CG axis (not ALA, GLY, PRO, SER, CYS, THR, VAL) (y/n)? ')1488:         check_file.write('Group rotation for CB - CG' + "\n")
1826:         if CB_CG in ['y', 'yes']:1489:         group_rotations[3] = True
1827:             check_file.write('Group rotation for CB - CG' + "\n")1490:         CB_CG_def = raw_input('Use default values (y/n)? ')
1828:             group_rotations[3] = True1491:         if CB_CG_def not in ['y' , 'yes']:
1829:             CB_CG_def = raw_input('Use default values (y/n)? ')1492:             def_values[3] = False
1830:             if CB_CG_def not in ['y', 'yes']:1493:     CG_CD = raw_input('Group rotations of side chains about CG-CD axis (only ARG, LYS, GLU, GLN) (y/n)? ')
1831:                 def_values[3] = False1494:     if CG_CD in ['y' , 'yes']:
1832:         CG_CD = raw_input('Group rotations of side chains about CG-CD axis (only ARG, LYS, GLU, GLN) (y/n)? ')1495:         check_file.write('Group rotation for CG - CD' + "\n")
1833:         if CG_CD in ['y', 'yes']:1496:         group_rotations[4] = True
1834:             check_file.write('Group rotation for CG - CD' + "\n")1497:         CG_CD_def = raw_input('Use default values (y/n)? ')
1835:             group_rotations[4] = True1498:         if CG_CD_def not in ['y' , 'yes']:
1836:             CG_CD_def = raw_input('Use default values (y/n)? ')1499:             def_values[4] = False
1837:             if CG_CD_def not in ['y', 'yes']:1500: 
1838:                 def_values[4] = False1501:     for res_gr in grot_res:
1839: 1502:         gr_res_name = loaded_mol.get_residue_name(res_gr)
1840:         for res_gr in grot_res:1503:         if group_rotations[0] and gr_res_name not in ['PRO']:
1841:             gr_res_name = loaded_mol.get_residue_name(res_gr)1504:             #N-CA rotation
1842:             if gr_res_name in protein_res:1505:             if def_values[0]:
1843:                 if group_rotations[0] and gr_res_name not in ['PRO']:1506:                 prob = prob_uni
1844:                     # N-CA rotation1507:                 amp = 0.1 #use default amplitude
1845:                     if def_values[0]:1508:             else:
1846:                         prob = prob_uni1509:                 prob,amp = get_prob_amb()
1847:                         amp = 0.1  # use default amplitude1510:             ax1 = loaded_mol.get_atom_id_from_name('N',res_gr)
1848:                     else:1511:             ax2 = loaded_mol.get_atom_id_from_name('CA',res_gr)
1849:                         prob, amp = get_prob_amb()1512:             gr_name = 'NCA'+str(res_gr)
1850:                     ax1 = loaded_mol.get_atom_id_from_name('N', res_gr)1513:             gr_atoms = []
1851:                     ax2 = loaded_mol.get_atom_id_from_name('CA', res_gr)1514:             for gr_atom in loaded_mol.get_atom_list(res_gr):
1852:                     gr_name = 'NCA' + str(res_gr)1515:                 if atom_dic[gr_atom][0] not in ['N','H','C','CA','HA','O','OXT','H1','H2','H3']:
1853:                     gr_atoms = []1516:                     gr_atoms.append(gr_atom)
1854:                     for gr_atom in loaded_mol.get_atom_list(res_gr):1517:             grot_dic[len(grot_dic)+1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]
1855:                         if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'O', 'OXT', 'H1', 'H2', 'H3']:1518: 
1856:                             gr_atoms.append(gr_atom)1519:         if group_rotations[1] and gr_res_name not in ['PRO']:
1857:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]1520:             #C-CA rotation
1858: 1521:             if def_values[1]:
1859:                 if group_rotations[1] and gr_res_name not in ['PRO']:1522:                 prob = prob_uni
1860:                     # C-CA rotation1523:                 amp = 0.1 #use default amplitude
1861:                     if def_values[1]:1524:             else:
1862:                         prob = prob_uni1525:                 prob,amp = get_prob_amb()
1863:                         amp = 0.1  # use default amplitude1526:             ax1 = loaded_mol.get_atom_id_from_name('C',res_gr)
1864:                     else:1527:             ax2 = loaded_mol.get_atom_id_from_name('CA',res_gr)
1865:                         prob, amp = get_prob_amb()1528:             gr_name ='CCA'+str(res_gr)
1866:                     ax1 = loaded_mol.get_atom_id_from_name('C', res_gr)1529:             gr_atoms = []
1867:                     ax2 = loaded_mol.get_atom_id_from_name('CA', res_gr)1530:             for gr_atom in loaded_mol.get_atom_list(res_gr):
1868:                     gr_name = 'CCA' + str(res_gr)1531:                 if atom_dic[gr_atom][0] not in ['N','H','C','CA','HA','O','OXT','H1','H2','H3']:
1869:                     gr_atoms = []1532:                     gr_atoms.append(gr_atom)
1870:                     for gr_atom in loaded_mol.get_atom_list(res_gr):1533:             grot_dic[len(grot_dic)+1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]
1871:                         if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'O', 'OXT', 'H1', 'H2', 'H3']:1534: 
1872:                             gr_atoms.append(gr_atom)1535:         if group_rotations[2] and gr_res_name not in ['PRO','ALA','GLY']:
1873:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]1536:             #CA-CB rotation
1874: 1537:             if def_values[2]:
1875:                 if group_rotations[2] and gr_res_name not in ['PRO', 'ALA', 'GLY']:1538:                 prob = prob_uni
1876:                     # CA-CB rotation1539:                 try:
1877:                     if def_values[2]:1540:                     amp = amp_dic[gr_res_name][0] #use default amplitude
1878:                         prob = prob_uni1541:                 except KeyError:
1879:                         try:1542:                     print 'No default data available for ' + gr_res_name
1880:                             amp = amp_dic[gr_res_name][0]  # use default amplitude1543:                     prob,amp = get_prob_amb()
1881:                         except KeyError:1544:             else:
1882:                             print 'No default data available for ' + gr_res_name1545:                 prob,amp = get_prob_amb()
1883:                             prob, amp = get_prob_amb()1546:             ax1 = loaded_mol.get_atom_id_from_name('CA',res_gr)
1884:                     else:1547:             ax2 = loaded_mol.get_atom_id_from_name('CB',res_gr)
1885:                         prob, amp = get_prob_amb()1548:             gr_name = 'CACB'+str(res_gr)
1886:                     ax1 = loaded_mol.get_atom_id_from_name('CA', res_gr)1549:             gr_atoms = []
1887:                     ax2 = loaded_mol.get_atom_id_from_name('CB', res_gr)1550:             for gr_atom in loaded_mol.get_atom_list(res_gr):
1888:                     gr_name = 'CACB' + str(res_gr)1551:                 if atom_dic[gr_atom][0] not in ['N','H','C','CA','HA','CB','O','OXT','H1','H2','H3']:
1889:                     gr_atoms = []1552:                     gr_atoms.append(gr_atom)
1890:                     for gr_atom in loaded_mol.get_atom_list(res_gr):1553:             grot_dic[len(grot_dic)+1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]
1891:                         if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'CB', 'O', 'OXT', 'H1', 'H2', 'H3']:1554: 
1892:                             gr_atoms.append(gr_atom)1555:         if group_rotations[3] and gr_res_name not in ['PRO','ALA','GLY','SER','CYS','THR','VAL']:
1893:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]1556:             #CB-CG rotation
1894: 1557:             if def_values[3]:
1895:                 if group_rotations[3] and gr_res_name not in ['PRO', 'ALA', 'GLY', 'SER', 'CYS', 'THR', 'VAL']:1558:                 prob = prob_uni
1896:                     # CB-CG rotation1559:                 try:
1897:                     if def_values[3]:1560:                     amp = amp_dic[gr_res_name][1] #use default amplitude
1898:                         prob = prob_uni1561:                 except KeyError:
1899:                         try:1562:                     print 'No default data available for ' + gr_res_name
1900:                             amp = amp_dic[gr_res_name][1]  # use default amplitude1563:                     prob,amp = get_prob_amb()
1901:                         except KeyError:1564:             else:
1902:                             print 'No default data available for ' + gr_res_name1565:                 prob,amp = get_prob_amb()
1903:                         prob, amp = get_prob_amb()1566:             ax1 = loaded_mol.get_atom_id_from_name('CB',res_gr)
1904:                     else:1567:             ax2 = loaded_mol.get_atom_id_from_name('CG',res_gr)
1905:                         prob, amp = get_prob_amb()1568:             gr_name = 'CBCG'+str(res_gr)
1906:                     ax1 = loaded_mol.get_atom_id_from_name('CB', res_gr)1569:             gr_atoms = []
1907:                     ax2 = loaded_mol.get_atom_id_from_name('CG', res_gr)1570:             for gr_atom in loaded_mol.get_atom_list(res_gr):
1908:                     gr_name = 'CBCG' + str(res_gr)1571:                 if atom_dic[gr_atom][0] not in ['N','H','C','CA','HA','CB','CG','HB','HB1','HB2','HB3','O','OXT','H1','H2','H3']:
1909:                     gr_atoms = []1572:                     gr_atoms.append(gr_atom)
1910:                     for gr_atom in loaded_mol.get_atom_list(res_gr):1573:             grot_dic[len(grot_dic)+1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]
1911:                         if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'CB', 'CG', 'HB', 'HB1', 'HB2',1574: 
1912:                                                         'HB3',1575:         if group_rotations[4] and gr_res_name in ['ARG','LYS','GLU','GLN']:
1913:                                                         'O', 'OXT', 'H1', 'H2', 'H3']:1576:             #CG-CD rotation
1914:                             gr_atoms.append(gr_atom)1577:             if def_values[4]:
1915:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]1578:                 prob = prob_uni
1916: 1579:                 try:
1917:                 if group_rotations[4] and gr_res_name in ['ARG', 'LYS', 'GLU', 'GLN']:1580:                     amp = amp_dic[gr_res_name][2] #use default amplitude
1918:                     # CG-CD rotation1581:                 except KeyError:
1919:                     if def_values[4]:1582:                     print 'No default data available for ' + gr_res_name
1920:                         prob = prob_uni1583:                     prob,amp = get_prob_amb()
1921:                         try:1584:             else:
1922:                             amp = amp_dic[gr_res_name][2]  # use default amplitude1585:                 prob,amp = get_prob_amb()
1923:                         except KeyError:1586:             ax1 = loaded_mol.get_atom_id_from_name('CG',res_gr)
1924:                             print 'No default data available for ' + gr_res_name1587:             ax2 = loaded_mol.get_atom_id_from_name('CD',res_gr)
1925:                             prob, amp = get_prob_amb()1588:             gr_name = 'CGCD'+str(res_gr)
1926:                     else:1589:             gr_atoms = []
1927:                         prob, amp = get_prob_amb()1590:             atom_exc = ['N','H','C','CA','HA','CB','CG','HB','HB1','HB2','HB3','O','OXT','H1','H2','H3','CD','HG1','HG2','HG11','HG12','HG13','HG21','HG22','HG23']
1928:                     ax1 = loaded_mol.get_atom_id_from_name('CG', res_gr)1591:             for gr_atom in loaded_mol.get_atom_list(res_gr):
1929:                     ax2 = loaded_mol.get_atom_id_from_name('CD', res_gr)1592:                 if atom_dic[gr_atom][0] not in atom_exc:
1930:                     gr_name = 'CGCD' + str(res_gr)1593:                     gr_atoms.append(gr_atom)
1931:                     gr_atoms = []1594:             grot_dic[len(grot_dic)+1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]
1932:                     atom_exc = ['N', 'H', 'C', 'CA', 'HA', 'CB', 'CG', 'HB', 'HB1', 'HB2', 'HB3', 'O', 'OXT', 'H1', 
1933:                                 'H2', 'H3', 'CD', 
1934:                                 'HG1', 'HG2', 'HG11', 'HG12', 'HG13', 'HG21', 'HG22', 'HG23'] 
1935:                     for gr_atom in loaded_mol.get_atom_list(res_gr): 
1936:                         if atom_dic[gr_atom][0] not in atom_exc: 
1937:                             gr_atoms.append(gr_atom) 
1938:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
1939:  
1940:     if rna_t or dna_t: 
1941:  
1942:         SugarBaseCheck = raw_input("Group rotations of base about C1'-N axis (y/n)? ") 
1943:         if SugarBaseCheck in ['y', 'yes']: 
1944:             check_file.write("Group rotation for C1' - N" + "\n") 
1945:             group_rotations_NA[0] = True 
1946:             CG_CD_def = raw_input('Use default values (y/n)? ') 
1947:             if CG_CD_def not in ['y', 'yes']: 
1948:                 def_values_NA[0] = False 
1949:  
1950:         RotBaseCheck = raw_input("Group rotations of base about intracyclic axis (y/n)? ") 
1951:         if RotBaseCheck in ['y', 'yes']: 
1952:             check_file.write("Group rotation for intracyclic axis in base" + "\n") 
1953:             group_rotations_NA[1] = True 
1954:             group_rotations_NA[2] = True 
1955:             CG_CD_def = raw_input('Use default values (y/n)? ') 
1956:             if CG_CD_def not in ['y', 'yes']: 
1957:                 def_values_NA[1] = False 
1958:                 def_values_NA[2] = False 
1959:  
1960:         SugarBackCheck = raw_input("Group rotations about C4'-C3' axis (y/n)? ") 
1961:         if SugarBackCheck in ['y', 'yes']: 
1962:             check_file.write("Group rotation for C4'-C3' " + "\n") 
1963:             group_rotations_NA[3] = True 
1964:             CG_CD_def = raw_input('Use default values (y/n)? ') 
1965:             if CG_CD_def not in ['y', 'yes']: 
1966:                 def_values_NA[3] = False 
1967:  
1968:         for res_gr in grot_res: 
1969:             gr_res_name = loaded_mol.get_residue_name(res_gr) 
1970:             if gr_res_name in ['RAN', 'A', 'A3', 'A5', 'AN', 'RA', 'RA3', 'RA5', 'DA', 'DA5', 'DA3', 'DAN']: 
1971:                 gr_res_name_dummy = 'A' 
1972:             elif gr_res_name in ['RCN', 'C', 'C3', 'C5', 'CN', 'RC', 'RC3', 'RC5', 'DC', 'DC5', 'DC3', 'DCN']: 
1973:                 gr_res_name_dummy = 'C' 
1974:             elif gr_res_name in ['RGN', 'G', 'G3', 'G5', 'GN', 'RG', 'RG3', 'RG5', 'DG', 'DG5', 'DG3', 'DGN']: 
1975:                 gr_res_name_dummy = 'G' 
1976:             elif gr_res_name in ['RCN', 'U', 'U3', 'U5', 'UN', 'RU', 'RU3', 'RU5']: 
1977:                 gr_res_name_dummy = 'U' 
1978:             elif gr_res_name in ['DT', 'DT5', 'DT3', 'DTN']: 
1979:                 gr_res_name_dummy = 'T' 
1980:  
1981:             if gr_res_name in NA_res: 
1982:  
1983:                 if group_rotations_NA[0]: 
1984:                     # C1'- N rotation 
1985:                     if def_values[1]: 
1986:                         try: 
1987:                             prob = prob_uni 
1988:                             amp = amp_dic_NA[gr_res_name_dummy][0]  # use default amplitude 
1989:                         except KeyError: 
1990:                             print 'No default data available for ' + gr_res_name 
1991:                             prob, amp = get_prob_amb() 
1992:                     else: 
1993:                         prob, amp = get_prob_amb() 
1994:                     if gr_res_name_dummy in ['A', 'G']: 
1995:                         ax1 = loaded_mol.get_atom_id_from_name('N9', res_gr) 
1996:                         gr_name = 'N9C1p' + str(res_gr) 
1997:                     elif gr_res_name_dummy in ['C', 'U', 'T']: 
1998:                         ax1 = loaded_mol.get_atom_id_from_name('N1', res_gr) 
1999:                         gr_name = 'N1C1p' + str(res_gr) 
2000:                     ax2 = loaded_mol.get_atom_id_from_name("C1'", res_gr) 
2001:                     gr_atoms = [] 
2002:                     for gr_atom in loaded_mol.get_atom_list(res_gr): 
2003:                         if atom_dic[gr_atom][0] not in ["P", "OP1", "OP2", "O5'", "C5'", "H5'", "H5''", "C4'", "H4'", 
2004:                                                         "O4'", "C1'", "H1'", "C3'", "H3'", "O3'", "C2'", "H2'", "O2'", 
2005:                                                         "HO2'", "HO3'", "HO5'", "H2''"]: 
2006:                             gr_atoms.append(gr_atom) 
2007:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
2008:  
2009:                 if group_rotations_NA[1]: 
2010:                     # intracyclic rotation 
2011:                     if def_values[1]: 
2012:                         try: 
2013:                             prob = prob_uni 
2014:                             amp = amp_dic_NA[gr_res_name_dummy][1]  # use default amplitude 
2015:                         except KeyError: 
2016:                             print 'No default data available for ' + gr_res_name 
2017:                             prob, amp = get_prob_amb() 
2018:                     else: 
2019:                         prob, amp = get_prob_amb() 
2020:                     if gr_res_name_dummy in ['A', 'G']: 
2021:                         ax1 = loaded_mol.get_atom_id_from_name('N9', res_gr) 
2022:                         ax2 = loaded_mol.get_atom_id_from_name('C6', res_gr) 
2023:                         gr_name = 'N9C6' + str(res_gr) 
2024:                     elif gr_res_name_dummy in ['C', 'U', 'T']: 
2025:                         ax1 = loaded_mol.get_atom_id_from_name('N1', res_gr) 
2026:                         ax2 = loaded_mol.get_atom_id_from_name('C4', res_gr) 
2027:                         gr_name = 'N1C4' + str(res_gr) 
2028:                     gr_atoms = [] 
2029:                     for gr_atom in loaded_mol.get_atom_list(res_gr): 
2030:                         if atom_dic[gr_atom][0] not in ["P", "OP1", "OP2", "O5'", "C5'", "H5'", "H5''", "C4'", "H4'", 
2031:                                                         "O4'", "C1'", "H1'", "C3'", "H3'", "O3'", "C2'", "H2'", "O2'", 
2032:                                                         "HO2'", "HO3'", "HO5'", "H2''"]: 
2033:                             gr_atoms.append(gr_atom) 
2034:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
2035:  
2036:                 if group_rotations_NA[2]: 
2037:                     # intracyclic rotation 2 
2038:                     if def_values[2]: 
2039:                         try: 
2040:                             prob = prob_uni 
2041:                             amp = amp_dic_NA[gr_res_name_dummy][2]  # use default amplitude 
2042:                         except KeyError: 
2043:                             print 'No default data available for ' + gr_res_name 
2044:                             prob, amp = get_prob_amb() 
2045:                     else: 
2046:                         prob, amp = get_prob_amb() 
2047:                     if gr_res_name_dummy in ['A', 'G']: 
2048:                         ax1 = loaded_mol.get_atom_id_from_name('N9', res_gr) 
2049:                         ax2 = loaded_mol.get_atom_id_from_name('N1', res_gr) 
2050:                         gr_name = 'N9N1' + str(res_gr) 
2051:                     elif gr_res_name_dummy in ['C', 'U', 'T']: 
2052:                         ax1 = loaded_mol.get_atom_id_from_name('N1', res_gr) 
2053:                         ax2 = loaded_mol.get_atom_id_from_name('N3', res_gr) 
2054:                         gr_name = 'N1N3' + str(res_gr) 
2055:                     gr_atoms = [] 
2056:                     for gr_atom in loaded_mol.get_atom_list(res_gr): 
2057:                         if atom_dic[gr_atom][0] not in ["P", "OP1", "OP2", "O5'", "C5'", "H5'", "H5''", "C4'", "H4'", 
2058:                                                         "O4'", "C1'", "H1'", "C3'", "H3'", "O3'", "C2'", "H2'", "O2'", 
2059:                                                         "HO2'", "HO3'", "HO5'", "H2''"]: 
2060:                             gr_atoms.append(gr_atom) 
2061:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
2062:  
2063:                 if group_rotations_NA[3]: 
2064:                     # sugarback bone rotation 
2065:                     if def_values[3]: 
2066:                         try: 
2067:                             prob = prob_uni 
2068:                             amp = amp_dic_NA[gr_res_name_dummy][3]  # use default amplitude 
2069:                         except KeyError: 
2070:                             print 'No default data available for ' + gr_res_name 
2071:                             prob, amp = get_prob_amb() 
2072:                     else: 
2073:                         prob, amp = get_prob_amb() 
2074:                     ax1 = loaded_mol.get_atom_id_from_name("C4'", res_gr) 
2075:                     ax2 = loaded_mol.get_atom_id_from_name("C3'", res_gr) 
2076:                     gr_name = 'N9N1' + str(res_gr) 
2077:  
2078:                     gr_atoms = [] 
2079:                     for gr_atom in loaded_mol.get_atom_list(res_gr): 
2080:                         if atom_dic[gr_atom][0] not in ["P", "OP1", "OP2", "O5'", "C5'", "H5'", "H5''", "C4'", 
2081:                                                         "C3'", "H3'", "O3'", "HO3'", "HO5'"]: 
2082:                             gr_atoms.append(gr_atom) 
2083:                     grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
2084: 1595: 
2085:     # user defined groups1596:     #user defined groups
2086:     gr_user_inp = raw_input('Add user defined rotational groups? ')1597:     gr_user_inp = raw_input('Add user defined rotational groups? ')
2087:     if gr_user_inp in ['y', 'yes']:1598:     if gr_user_inp in ['y','yes']:
2088:         gr_user_t = True1599:         gr_user_t = True
2089:     else:1600:     else:
2090:         gr_user_t = False1601:         gr_user_t = False
2091:     useri = 11602:     useri = 1
2092:     while gr_user_t:1603:     while gr_user_t:
2093:         print 'No tests for valid groups!'1604:         print 'No tests for valid groups!'
2094:         gr_name = 'USERD' + str(useri)1605:         gr_name = 'USERD' + str(useri)
2095:         try:1606:         try:
2096:             ax1 = int(raw_input('Atom number for axis1: '))1607:             ax1 = int(raw_input('Atom number for axis1: '))
2097:             ax2 = int(raw_input('Atom number for axis2: '))1608:             ax2 = int(raw_input('Atom number for axis2: '))
2098:             prob, amp = get_prob_amb()1609:             prob,amp = get_prob_amb()
2099:             gr_atoms = [int(x) for x in (raw_input('List of atoms in group separated by spaces: ')).split()]1610:             gr_atoms = [int(x) for x in split(raw_input('List of atoms in group separated by spaces: '))]
2100:         except:1611:         except:
2101:             print 'Invalid input'1612:             print 'Invalid input'
2102:         else:1613:         else:
2103:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]1614:             grot_dic[len(grot_dic)+1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms]
2104:         finally:1615:         finally:
2105:             gr_user_inp = raw_input('Add user defined rotational group? ')1616:             gr_user_inp = raw_input('Add user defined rotational group? ')
2106:             if gr_user_inp in ['y', 'yes']:1617:             if gr_user_inp in ['y','yes']:
2107:                 gr_user_t = True1618:                 gr_user_t = True
2108:             else:1619:             else:
2109:                 gr_user_t = False1620:                 gr_user_t = False
2110: 1621: 
2111:     # write output for group rotation files1622:     #write output for group rotation files
2112:     groups_file = open('atomgroups', 'w')1623:     groups_file = open('atomgroups', 'w')
2113:     for group in range(1, len(grot_dic.keys()) + 1):1624:     for group in range(1, len(grot_dic.keys()) + 1):
2114:         string = 'GROUP ' + grot_dic[group][0] + ' ' + str(grot_dic[group][1]) + ' ' + str(grot_dic[group][2])1625:         string = 'GROUP ' + grot_dic[group][0] + ' ' + str(grot_dic[group][1]) + ' ' + str(grot_dic[group][2])
2115:         string += ' ' + str(grot_dic[group][3]) + ' ' + str(grot_dic[group][4]) + ' ' + str(grot_dic[group][5]) + "\n"1626:         string += ' ' + str(grot_dic[group][3]) + ' ' + str(grot_dic[group][4]) + ' ' + str(grot_dic[group][5]) + "\n"
2116:         for atom in grot_dic[group][6]:1627:         for atom in grot_dic[group][6]:
2117:             string += str(atom) + "\n"1628:             string += str(atom) + "\n"
2118:         groups_file.write(string)1629:         groups_file.write(string)
2119: 1630: 
2120: check_file.close()1631: check_file.close()
2121: 1632: 
2122: # Writing in file parmed.py to exclude interactions ####1633: ### Writing in file parmed.py to exclude interactions ####
2123: parmed_in = open('parmed_in', "w")1634: parmed_in = open('parmed_in', "w")
2124: for group in range(1, len(groups.keys()) + 1):1635: for group in range(1, len(groups.keys()) + 1):
2125:     string = 'addExclusions '1636:     string = 'addExclusions '
2126:     string2 = '@'1637:     string2 = '@'
2127:     counter = 11638:     counter = 1
2128:     for res in groups[group]:1639:     for res in groups[group]:
2129:         for atom in res_dic[res]:1640:         for atom in res_dic[res]:
2130:             if counter == 1:1641:             if counter == 1:
2131:                 string2 += str(atom)1642:                 string2 += str(atom)
2132:                 counter = 21643:                 counter = 2
2133:             else:1644:             else:
2134:                 string2 += ',' + str(atom)1645:                 string2 += ',' + str(atom)
2135:     parmed_in.write(string + string2 + ' ' + string2 + "\n")1646:     parmed_in.write(string + string2 + ' ' +string2 +"\n")
2136: for group_local in range(1, len(groups_local.keys()) + 1):1647: for group_local in range(1, len(groups_local.keys()) + 1):
2137:     string = 'addExclusions '1648:     string = 'addExclusions '
2138:     string2 = '@'1649:     string2 = '@'
2139:     counter = 11650:     counter = 1
2140:     for atom in groups_local[group_local][2:]:1651:     for atom in groups_local[group_local][2:]:
2141:         if counter == 1:1652:         if counter == 1:
2142:             string2 += str(atom)1653:             string2 += str(atom)
2143:             counter = 21654:             counter = 2
2144:         else:1655:         else:
2145:             string2 += ',' + str(atom)1656:             string2 += ',' + str(atom)
2146:     parmed_in.write(string + string2 + ' ' + string2 + "\n")1657:     parmed_in.write(string + string2 + ' ' + string2 + "\n")
2147: parmed_in.write('parmout coords.prmtop.ni' + "\n" + 'go' + "\n")1658: parmed_in.write('parmout coords.prmtop.ni' + "\n" + 'go' + "\n")
2148: 1659: 
 1660: 
2149: print 'Selection process completed - change topology file to exclude interactions within the rigid bodies.'1661: print 'Selection process completed - change topology file to exclude interactions within the rigid bodies.'
2150: 1662: 
2151: if pymol_check:1663: if pymol_check:
2152:     mol_view = PymolView(pdb_inp)1664:     mol_view = PymolView(pdb_inp)
2153:     mol_view.apply_colour_scheme(res_atomistic, res_local, 2, 3)1665:     mol_view.apply_colour_scheme(res_atomistic, res_eical, 2, 3)
2154:     pymol.cmd.save('rigid_%s' % pdb_inp, format='png')1666:     pymol.cmd.save('rigid_%s' % pdb_inp, format='png')
2155:     time.sleep(10)1667:     time.sleep(10)
2156:     mol_view.__del__()1668:     mol_view.__del__()


r31935/README_genrigid 2017-02-17 12:30:08.314888000 +0000 r31934/README_genrigid 2017-02-17 12:30:08.762893999 +0000
  2:   2: 
  3: HOW TO USE IT  3: HOW TO USE IT
  4:   4: 
  5: 1. Either start genrigid_input.py or start with options python  genrigid_input.py [pdb name] [coords.inpcrd] [tolerance] [pymol].  5: 1. Either start genrigid_input.py or start with options python  genrigid_input.py [pdb name] [coords.inpcrd] [tolerance] [pymol].
  6:         [pdb name]        required file name  6:         [pdb name]        required file name
  7:         [coords.inpcrd]   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-  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  9:                           linear), default 0.01
 10:         [pymol]           requires pymol installed, 'pymol' loads a graphical representation of the set rigidification 10:         [pymol]           requires pymol installed, 'pymol' loads a graphical representation of the set rigidification
 11:  11: 
 12:    Then choose what kind of system you have: RNA, DNA, protein or a mixed system. (This refers to which regions you want to consider) 
 13:  
 14: 2. Check number of residues and number of atoms, if they are not what they should be, check the input files. 12: 2. Check number of residues and number of atoms, if they are not what they should be, check the input files.
 15:  13: 
 16: 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  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 
 17:    keys add these residues to the atomistic list. 15:    keys add these residues to the atomistic list.
 18:  16: 
 19: 4. Select and remove from selection following the script. The residue selections can be either one by one with spaces or using 'r integer1 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
 20:    integer2' to indicate the range from integer1 to integer2. If a residue is selected to be atomistic, selecting it again for the local rigid 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
 21:    bodies will not change the assignement (this change can only be done by removing it and then add it again). 19:    bodies will not change the assignement (this change can only be done by removing it and then add it again).
 22:  20: 
 23: 5. Grouping of the fully rigid region: enter residue by residue or ranges for option 2. 21: 5. Grouping of the fully rigid region: enter residue by residue or ranges for option 2.
 24:  22: 
 25: 6. Local rigid bodies: peptide bonds will be rigidified for all residues in the local rigid selection. If the neighbouring residue belongs to 23: 6. Local rigid bodies: peptide bonds will be rigidified for all residues in the local rigid selection. If the neighbouring residue belongs to
 26:    a fully rigid body, this peptide bond will not be used. If it belongs to an atomistic residue it will be used. 24:    a fully rigid body, this peptide bond will not be used. If it belongs to an atomistic residue it will be used.
 27:  25: 
 28: 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    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   
 29:    selection is ignored. 27:    selection is ignored.
 30:  28: 
 31: 8. After the rigid bodies are created you can choose to create a group rotation file. 
 32:  
 33: HOW TO CHANGE IT 29: HOW TO CHANGE IT
 34:  30: 
 35: 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. 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