hdiff output

r33473/rbody_grot.py 2017-11-14 13:30:13.828917412 +0000 r33472/rbody_grot.py 2017-11-14 13:30:14.268923265 +0000
  1: import sys  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/AMBER/BHmutation_steps/rbody_grot.py' in revision 33472
  2:  
  3: class protein(): 
  4:     def __init__(self, atom_dic, res_dic): 
  5:         self.atom_dic = atom_dic  # dictionary of all atoms with x,y,z coordinates and res id and name and atom name 
  6:         self.res_dic = res_dic  # dictionary containing a list of atom ids for all residues 
  7:  
  8:     def num_atoms(self):  # number of atoms 
  9:         return len(self.atom_dic) 
 10:  
 11:     def num_res(self):  # number of residues 
 12:         return len(self.res_dic) 
 13:  
 14:     # returns the name of a given residue 
 15:     def get_residue_name(self, residue): 
 16:         return atom_dic[res_dic[residue][0]][1] 
 17:  
 18:     # atom id from atom name given a residue id 
 19:     def get_atom_id_from_name(self, atom_name, res_num): 
 20:         for i in self.res_dic[res_num]: 
 21:             if self.atom_dic[i][0] == atom_name: 
 22:                 return i 
 23:  
 24:     # returns a list of atoms in a residue 
 25:     def get_atom_list(self, residue): 
 26:         all_atoms = self.atom_dic.keys() 
 27:         atomlist = [] 
 28:         for i in range(1, len(self.atom_dic) + 1): 
 29:             atom = self.atom_dic[i] 
 30:             if residue == atom[2]: 
 31:                 atomlist.append(all_atoms[i - 1]) 
 32:         return atomlist 
 33:  
 34: ################################################# 
 35: def reading_pdb(pdb): 
 36:     dic = {} 
 37:     atom_list = [] 
 38:     with open(pdb, 'r') as f_pdb: 
 39:         for line in f_pdb: 
 40:             if line.split()[0] == 'ATOM': 
 41:                 atom_list.append((line[12:16].strip(), line[17:20].strip(), int(line[22:26]))) 
 42:     for i in range(1, len(atom_list) + 1): 
 43:         dic[i] = atom_list[i - 1] 
 44:     return dic 
 45:  
 46: def reading_rbodyconfig(fname): 
 47:     dic = dict() 
 48:     with open(fname , "r") as f: 
 49:         for line in f: 
 50:             if line.split()[0] == 'GROUP': 
 51:                 dic[len(dic) + 1] = list() 
 52:             else: 
 53:                 dic[len(dic)].append(int(line.split()[0])) 
 54:     return dic 
 55:  
 56: def create_res_dic(atom_dic): 
 57:     res_dic = {} 
 58:     res_list = [] 
 59:     for i in atom_dic.keys(): 
 60:         res_list.append(atom_dic[i][2]) 
 61:     res_list = list(set(res_list)) 
 62:     for j in res_list: 
 63:         res_dic[j] = [] 
 64:     for k in atom_dic.keys(): 
 65:         res_dic[atom_dic[k][2]].append(k) 
 66:     return res_dic 
 67:  
 68: #################################################        
 69:  
 70: prob_uni = 0.025 
 71:  
 72: rbody_old = sys.argv[1] 
 73: shift = int(sys.argv[2]) 
 74: start = int(sys.argv[3]) 
 75: pdb_file = sys.argv[4] 
 76: atom_dic = reading_pdb(pdb_file) 
 77: res_dic = create_res_dic(atom_dic) 
 78: loaded_mol = protein(atom_dic, res_dic) 
 79:  
 80: amp_dic = {'ALA': (1.0,), 
 81:            'ASN': (0.5, 1.0), 
 82:            'ARG': (0.2, 0.3, 0.5), 
 83:            'ASP': (0.5, 1.0, 0.0), 
 84:            'CYS': (1.0,), 
 85:            'GLN': (0.3, 0.5, 1.0), 
 86:            'GLU': (0.3, 0.5, 1.0), 
 87:            'HIS': (0.3, 0.5), 
 88:            'HIE': (0.3, 0.5), 
 89:            'HIP': (0.3, 0.5), 
 90:            'HID': (0.3, 0.5), 
 91:            'ILE': (0.5, 1.0), 
 92:            'LEU': (0.5, 1.0), 
 93:            'LYS': (0.2, 0.3, 0.5), 
 94:            'MET': (0.5, 0.7), 
 95:            'PHE': (0.3, 0.5), 
 96:            'SER': (1.0,), 
 97:            'THR': (1.0,), 
 98:            'TRP': (0.3, 0.4), 
 99:            'TYR': (0.3, 0.5), 
100:            'VAL': (1.0,), 
101:            'AMB0': (0.1,)  # default for all non residue specific rotations 
102:           } 
103:  
104: #first read in old rbodyconfig and adjust it to give new output 
105: rb_dict = reading_rbodyconfig(rbody_old) 
106: used_atoms = list() 
107: for key in rb_dict.keys(): 
108:     if rb_dict[key][0] > start: 
109:         rb_dict[key] = [(atom + shift) for atom in rb_dict[key]] 
110:     used_atoms = used_atoms + rb_dict[key] 
111:  
112: rboutput = open("rbodyconfig.new" , "w") 
113: for key in rb_dict.keys(): 
114:     rboutput.write("GROUP " + str(len(rb_dict[key])) + "\n") 
115:     for atom in rb_dict[key]: 
116:         rboutput.write(str(atom) + "\n") 
117: rboutput.close() 
118:  
119: # number: [name, ax_atom1, ax_atom2, natom, amp, prob, [atom1, atom2, ...]] 
120: grot_dic = {} 
121: residues = res_dic.keys() 
122: for res in residues: 
123:     gr_res_name = loaded_mol.get_residue_name(res) 
124:     if gr_res_name not in ['PRO' , 'HYP']: 
125:         #N-CA rotation 
126:         prob = prob_uni 
127:         amp = 0.1  # use default amplitude 
128:  
129:         ax1 = loaded_mol.get_atom_id_from_name('N', res) 
130:         ax2 = loaded_mol.get_atom_id_from_name('CA', res) 
131:         gr_name = 'NCA' + str(res) 
132:         gr_atoms = [] 
133:         for gr_atom in loaded_mol.get_atom_list(res): 
134:             if atom_dic[gr_atom][0] not in ['N',  'C', 'CA', 'HA', 'O', 'OXT', 'H1', 'H2', 'H3']: 
135:                 gr_atoms.append(gr_atom) 
136:         inrbody = False 
137:         for atom in gr_atoms: 
138:             if atom in used_atoms: 
139:                 inrbody = True 
140:         if not(inrbody): 
141:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
142:  
143:     if gr_res_name not in ['PRO' , 'HYP']: 
144:         #C-CA rotation 
145:         prob = prob_uni 
146:         amp = 0.1  # use default amplitude 
147:         ax1 = loaded_mol.get_atom_id_from_name('C', res) 
148:         ax2 = loaded_mol.get_atom_id_from_name('CA', res) 
149:         gr_name = 'CCA' + str(res) 
150:         gr_atoms = [] 
151:         for gr_atom in loaded_mol.get_atom_list(res): 
152:             if atom_dic[gr_atom][0] not in ['N',  'C', 'CA', 'HA', 'O', 'OXT', 'H1', 'H2', 'H3']: 
153:                 gr_atoms.append(gr_atom) 
154:         inrbody = False 
155:         for atom in gr_atoms: 
156:             if atom in used_atoms: 
157:                 inrbody = True 
158:         if not(inrbody): 
159:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
160:  
161:     if gr_res_name not in ['PRO', 'HYP' , 'ALA', 'GLY']: 
162:         # CA-CB rotation 
163:         prob = prob_uni 
164:         try: 
165:             amp = amp_dic[gr_res_name][0]  # use default amplitude 
166:         except KeyError: 
167:             amp = 0.1 
168:         ax1 = loaded_mol.get_atom_id_from_name('CA', res) 
169:         ax2 = loaded_mol.get_atom_id_from_name('CB', res) 
170:         gr_name = 'CACB' + str(res) 
171:         gr_atoms = [] 
172:         for gr_atom in loaded_mol.get_atom_list(res): 
173:             if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'CB', 'O', 'OXT', 'H1', 'H2', 'H3']: 
174:                 gr_atoms.append(gr_atom) 
175:         inrbody = False 
176:         for atom in gr_atoms: 
177:             if atom in used_atoms: 
178:                 inrbody = True 
179:         if not(inrbody): 
180:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
181:  
182:     if gr_res_name not in ['PRO', 'HYP' , 'ILE', 'ALA', 'GLY', 'SER', 'CYS', 'THR', 'VAL']: 
183:         # CB-CG rotation 
184:         prob = prob_uni 
185:         try: 
186:             amp = amp_dic[gr_res_name][1]  # use default amplitude 
187:         except KeyError: 
188:             amp = 0.1 
189:         ax1 = loaded_mol.get_atom_id_from_name('CB', res) 
190:         ax2 = loaded_mol.get_atom_id_from_name('CG', res) 
191:         gr_name = 'CBCG' + str(res) 
192:         gr_atoms = [] 
193:         for gr_atom in loaded_mol.get_atom_list(res): 
194:             if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'CB', 'CG', 'HB', 'HB1', 'HB2', 
195:                                             'HB3', 'O', 'OXT', 'H1', 'H2', 'H3']: 
196:                 gr_atoms.append(gr_atom) 
197:         inrbody = False 
198:         for atom in gr_atoms: 
199:             if atom in used_atoms: 
200:                 inrbody = True 
201:         if not(inrbody): 
202:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
203:  
204:     if gr_res_name in ['ARG', 'LYS', 'GLU', 'GLN']: 
205:         # CG-CD rotation 
206:         prob = prob_uni 
207:         try: 
208:             amp = amp_dic[gr_res_name][2]  # use default amplitude 
209:         except KeyError: 
210:             amp = 0.1 
211:         ax1 = loaded_mol.get_atom_id_from_name('CG', res) 
212:         ax2 = loaded_mol.get_atom_id_from_name('CD', res) 
213:         gr_name = 'CGCD' + str(res) 
214:         gr_atoms = [] 
215:         atom_exc = ['N', 'H', 'C', 'CA', 'HA', 'CB', 'CG', 'HB', 'HB1', 'HB2', 'HB3', 'O', 'OXT', 'H1', 
216:                     'H2', 'H3', 'CD', 'HG1', 'HG2', 'HG11', 'HG12', 'HG13', 'HG21', 'HG22', 'HG23'] 
217:         for gr_atom in loaded_mol.get_atom_list(res): 
218:             if atom_dic[gr_atom][0] not in atom_exc: 
219:                 gr_atoms.append(gr_atom) 
220:                 inrbody = False 
221:         for atom in gr_atoms: 
222:             if atom in used_atoms: 
223:                 inrbody = True 
224:         if not(inrbody): 
225:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
226:  
227:    
228: # write output for group rotation files 
229: groups_file = open('atomgroups', 'w') 
230: for group in range(1, len(grot_dic.keys()) + 1): 
231:     string = 'GROUP ' + grot_dic[group][0] + ' ' + str(grot_dic[group][1]) + ' ' + str(grot_dic[group][2]) 
232:     string += ' ' + str(grot_dic[group][3]) + ' ' + str(grot_dic[group][4]) + ' ' + str(grot_dic[group][5]) + "\n" 
233:     for atom in grot_dic[group][6]: 
234:         string += str(atom) + "\n" 
235:     groups_file.write(string) 
236: groups_file.close() 
237:  
238:  
239:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0