hdiff output

r30309/perm-pdb.py 2017-01-21 10:36:59.518250000 +0000 r30308/perm-pdb.py 2017-01-21 10:37:00.002250000 +0000
 98:     count=2      # number of permutable atoms 98:     count=2      # number of permutable atoms
 99:     groupcount=0 # number of permutable atoms in group 99:     groupcount=0 # number of permutable atoms in group
100:     swap=0       # number of other pairs of atoms that must swap if the first pair is permuted100:     swap=0       # number of other pairs of atoms that must swap if the first pair is permuted
101:     group2=[]101:     group2=[]
102:     swap2=0102:     swap2=0
103:     groupcount2=0103:     groupcount2=0
104: 104: 
105:     ###############105:     ###############
106:     ## amino acids106:     ## amino acids
107:     ###############107:     ###############
108:     if ATMlist[0].acidname in ['NGL', 'CGL', 'NAR' , 'CAR' , 'NVA' ,'CVA' , 'NAS' , 'CAS' , 'NLE' , 'CLE' , 'NLY' , 'CLY' , 'NCY' , 'CCY' , 'NPR' , 'CPR' , 'NTH' , 'CTH' , 'NAL' , 'CAL' , 'NHS' , 'CHS' , 'NHI' , 'CHI', 'NTR', 'CTR' , 'NSE' , 'CSE' , 'NIL', 'CIL' , 'NME' , 'CME' , 'NPH' , 'CPH' , 'NTY' , 'CTY' , 'CHY' , 'NHY']:108:     if ATMlist[0].acidname in ['NGL', 'CGL', 'NAR' , 'CAR' , 'NVA' ,'CVA' , 'NAS' , 'CAS' , 'NLE' , 'CLE' , 'NLY' , 'CLY' , 'NCY' , 'CCY' , 'NPR' , 'CPR' , 'NTH' , 'CTH' , 'NAL' , 'CAL' , 'NHS' , 'CHS' , 'NHI' , 'CHI', 'NTR', 'CTR' , 'NSE' , 'CSE' , 'NIL', 'CIL' , 'NME' , 'CME' , 'NPH' , 'CPH' , 'NTY' , 'CTY' ]:
109:        ATMlist[0].acidname = ATMlist[0].acidname_ter109:        ATMlist[0].acidname = ATMlist[0].acidname_ter
110:     if(ATMlist[0].acidname=='GLN'):110:     if(ATMlist[0].acidname=='GLN'):
111:       atnum.append(ATMlist[5].index)111:       atnum.append(ATMlist[5].index)
112:       atnum.append(ATMlist[6].index)112:       atnum.append(ATMlist[6].index)
113:       count2=2113:       count2=2
114:       atnum2.append(ATMlist[8].index)114:       atnum2.append(ATMlist[8].index)
115:       atnum2.append(ATMlist[9].index)115:       atnum2.append(ATMlist[9].index)
116:       if (amber):116:       if (amber):
117:          count3=2117:          count3=2
118:          atnum3.append(ATMlist[13].index)118:          atnum3.append(ATMlist[13].index)
990:          atnum3.append(ATMlist[9].index)990:          atnum3.append(ATMlist[9].index)
991:       elif (charmm):991:       elif (charmm):
992:          atnum2.append(ATMlist[7].index)992:          atnum2.append(ATMlist[7].index)
993:          atnum2.append(ATMlist[8].index)993:          atnum2.append(ATMlist[8].index)
994:          count3=2994:          count3=2
995:          atnum3.append(ATMlist[10].index)995:          atnum3.append(ATMlist[10].index)
996:          atnum3.append(ATMlist[11].index)996:          atnum3.append(ATMlist[11].index)
997:       count4=2997:       count4=2
998:       atnum4.append(ATMlist[13].index)998:       atnum4.append(ATMlist[13].index)
999:       atnum4.append(ATMlist[14].index)999:       atnum4.append(ATMlist[14].index)
1000:     elif(ATMlist[0].acidname=='HYP'): 
1001:       atnum.append(ATMlist[2].index) 
1002:       atnum.append(ATMlist[3].index) 
1003:       count2=2 
1004:       if (amber): 
1005:          atnum3.append(ATMlist[9].index) 
1006:          atnum3.append(ATMlist[10].index) 
1007:       elif (charmm): 
1008:          print 'Check permutations for CHYP 
1009:       #   atnum2.append(ATMlist[7].index) 
1010:       #   atnum2.append(ATMlist[8].index) 
1011:       #   count3=2 
1012:       #   atnum3.append(ATMlist[10].index) 
1013:       #   atnum3.append(ATMlist[11].index) 
1014:     elif(ATMlist[0].acidname=='CHYP'): 
1015:       atnum.append(ATMlist[2].index) 
1016:       atnum.append(ATMlist[3].index) 
1017:       count2=2 
1018:       if (amber): 
1019:          atnum3.append(ATMlist[9].index) 
1020:          atnum3.append(ATMlist[10].index) 
1021:       elif (charmm): 
1022:           print 'Check permutations for CHYP' 
1023:       #   atnum2.append(ATMlist[7].index) 
1024:       #   atnum2.append(ATMlist[8].index) 
1025:       #   count3=2 
1026:       #   atnum3.append(ATMlist[10].index) 
1027:       #   atnum3.append(ATMlist[11].index) 
1028:       count3=2 
1029:       atnum4.append(ATMlist[14].index) 
1030:       atnum4.append(ATMlist[15].index) 
1031:     elif(ATMlist[0].acidname=='ACE'):1000:     elif(ATMlist[0].acidname=='ACE'):
1032:       count=31001:       count=3
1033:       atnum.append(ATMlist[0].index)1002:       atnum.append(ATMlist[0].index)
1034:       atnum.append(ATMlist[2].index)1003:       atnum.append(ATMlist[2].index)
1035:       atnum.append(ATMlist[3].index)1004:       atnum.append(ATMlist[3].index)
1036:     elif(ATMlist[0].acidname=='NME'):1005:     elif(ATMlist[0].acidname=='NME'):
1037:       count=31006:       count=3
1038:       atnum.append(ATMlist[3].index)1007:       atnum.append(ATMlist[3].index)
1039:       atnum.append(ATMlist[4].index)1008:       atnum.append(ATMlist[4].index)
1040:       atnum.append(ATMlist[5].index)1009:       atnum.append(ATMlist[5].index)
1121:     elif(ATMlist[0].acidname=='DT5' or ATMlist[0].acidname=='DTN'):1090:     elif(ATMlist[0].acidname=='DT5' or ATMlist[0].acidname=='DTN'):
1122:       atnum.append(ATMlist[3].index)1091:       atnum.append(ATMlist[3].index)
1123:       atnum.append(ATMlist[4].index)1092:       atnum.append(ATMlist[4].index)
1124:       count2=31093:       count2=3
1125:       atnum2.append(ATMlist[15].index)1094:       atnum2.append(ATMlist[15].index)
1126:       atnum2.append(ATMlist[16].index)1095:       atnum2.append(ATMlist[16].index)
1127:       atnum2.append(ATMlist[17].index)1096:       atnum2.append(ATMlist[17].index)
1128:       count4=21097:       count4=2
1129:       atnum4.append(ATMlist[27].index)1098:       atnum4.append(ATMlist[27].index)
1130:       atnum4.append(ATMlist[28].index)1099:       atnum4.append(ATMlist[28].index)
1131:     elif(ATMlist[0].acidname=='RA' or ATMlist[0].acidname=='RA3' or ATMlist[0].acidname=='RC' or ATMlist[0].acidname=='RC3' or ATMlist[0].acidname in ['A','A3','C','C3']):1100:     elif(ATMlist[0].acidname=='RA' or ATMlist[0].acidname=='RA3' or ATMlist[0].acidname=='RC' or ATMlist[0].acidname=='RC3'):
1132:       atnum.append(ATMlist[1].index)1101:       atnum.append(ATMlist[1].index)
1133:       atnum.append(ATMlist[2].index)1102:       atnum.append(ATMlist[2].index)
1134:       count2=21103:       count2=2
1135:       atnum2.append(ATMlist[5].index)1104:       atnum2.append(ATMlist[5].index)
1136:       atnum2.append(ATMlist[6].index)1105:       atnum2.append(ATMlist[6].index)
1137:       count3=21106:       count3=2
1138:       atnum3.append(ATMlist[19].index)1107:       atnum3.append(ATMlist[19].index)
1139:       atnum3.append(ATMlist[20].index)1108:       atnum3.append(ATMlist[20].index)
1140:     elif(ATMlist[0].acidname=='RA5' or ATMlist[0].acidname=='RAN' or ATMlist[0].acidname=='RC5' or ATMlist[0].acidname=='RCN' or ATMlist[0].acidname in ['A5','AN','C5','CN']):1109:     elif(ATMlist[0].acidname=='RA5' or ATMlist[0].acidname=='RAN' or ATMlist[0].acidname=='RC5' or ATMlist[0].acidname=='RCN'):
1141:       atnum.append(ATMlist[3].index)1110:       atnum.append(ATMlist[3].index)
1142:       atnum.append(ATMlist[4].index)1111:       atnum.append(ATMlist[4].index)
1143:       count2=21112:       count2=2
1144:       atnum2.append(ATMlist[17].index)1113:       atnum2.append(ATMlist[17].index)
1145:       atnum2.append(ATMlist[18].index)1114:       atnum2.append(ATMlist[18].index)
1146:     elif(ATMlist[0].acidname=='RG' or ATMlist[0].acidname=='RG3' or ATMlist[0].acidname in ['G','G3']):1115:     elif(ATMlist[0].acidname=='RG' or ATMlist[0].acidname=='RG3'):
1147:       atnum.append(ATMlist[1].index)1116:       atnum.append(ATMlist[1].index)
1148:       atnum.append(ATMlist[2].index)1117:       atnum.append(ATMlist[2].index)
1149:       count2=21118:       count2=2
1150:       atnum2.append(ATMlist[5].index)1119:       atnum2.append(ATMlist[5].index)
1151:       atnum2.append(ATMlist[6].index)1120:       atnum2.append(ATMlist[6].index)
1152:       count3=21121:       count3=2
1153:       atnum3.append(ATMlist[23].index)1122:       atnum3.append(ATMlist[23].index)
1154:       atnum3.append(ATMlist[24].index)1123:       atnum3.append(ATMlist[24].index)
1155:     elif(ATMlist[0].acidname=='RG5' or ATMlist[0].acidname=='RGN' or ATMlist[0].acidname in ['G5','GN']):1124:     elif(ATMlist[0].acidname=='RG5' or ATMlist[0].acidname=='RGN'):
1156:       atnum.append(ATMlist[3].index)1125:       atnum.append(ATMlist[3].index)
1157:       atnum.append(ATMlist[4].index)1126:       atnum.append(ATMlist[4].index)
1158:       count2=21127:       count2=2
1159:       atnum2.append(ATMlist[21].index)1128:       atnum2.append(ATMlist[21].index)
1160:       atnum2.append(ATMlist[22].index)1129:       atnum2.append(ATMlist[22].index)
1161:     elif(ATMlist[0].acidname=='RU' or ATMlist[0].acidname=='RU3' or ATMlist[0].acidname in ['U','U3']):1130:     elif(ATMlist[0].acidname=='RU' or ATMlist[0].acidname=='RU3'):
1162:       atnum.append(ATMlist[1].index)1131:       atnum.append(ATMlist[1].index)
1163:       atnum.append(ATMlist[2].index)1132:       atnum.append(ATMlist[2].index)
1164:       count2=21133:       count2=2
1165:       atnum2.append(ATMlist[5].index)1134:       atnum2.append(ATMlist[5].index)
1166:       atnum2.append(ATMlist[6].index)1135:       atnum2.append(ATMlist[6].index)
1167:     elif((ATMlist[0].acidname=='RU5') or (ATMlist[0].acidname=='RUN') or ATMlist[0].acidname in ['U5','UN']):1136:     elif((ATMlist[0].acidname=='RU5') or (ATMlist[0].acidname=='RUN')):
1168:       atnum.append(ATMlist[3].index)1137:       atnum.append(ATMlist[3].index)
1169:       atnum.append(ATMlist[4].index)1138:       atnum.append(ATMlist[4].index)
1170:     else:1139:     else:
1171:       print 'Neither amino acid nor nucleic residue - please check residue %s' % ATMlist[0].acidname1140:       print 'Neither amino acid nor nucleic residue - please check residue %s' % ATMlist[0].acidname
1172: 1141: 
1173: 1142: 
1174:     if(els[0]!='END'):1143:     if(els[0]!='END'):
1175:       ATMlist=[]1144:       ATMlist=[]
1176:       prev=each[22:26].strip()1145:       prev=each[22:26].strip()
1177:       AM=readatom(each)1146:       AM=readatom(each)


r30309/perm-prmtop.ff14.py 2017-01-21 10:36:59.274250000 +0000 r30308/perm-prmtop.ff14.py 2017-01-21 10:36:59.770250000 +0000
  1: #!/usr/bin/env python  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/AMBER/symmetrise_prmtop/perm-prmtop.ff14.py' in revision 30308
  2:  
  3: import os 
  4: import os.path 
  5: import sys 
  6: import string 
  7:  
  8: ############################################################### 
  9: ## Adapted from original to account for ff14SB                #  
 10: ## by K. Roeder (kr366)                                       # 
 11: ##                                                            # 
 12: ## Edyta Malolepsza                                           # 
 13: ## David Wales' group, University of Cambridge                # 
 14: ## in case of problems please send email: em427@cam.ac.uk     # 
 15: ##                                                            # 
 16: ############################################################### 
 17: ##                                                            # 
 18: ## program finds in prmtop file from LEaP wrong defined order # 
 19: ## of atoms in IMPROPER, permutes appropriate atoms and write # 
 20: ## new prmtop file                                            # 
 21: ##                                                            # 
 22: ## how to use:                                                # 
 23: ## ./perm-top.py NAME_OF_OLD_PRMTOP NAME_OF_NEW_PRMTOP        # 
 24: ##                                                            # 
 25: ## IMPORTANT:                                                 # 
 26: ## 1. please change names of terminal amino acid residues     # 
 27: ##    according to warnings below                             # 
 28: ## 2. please change path to libraries                         # 
 29: ## 3. program changes the atom order ONLY for amino acid and  # 
 30: ##    nucleic residues                                        # 
 31: ##                                                            # 
 32: ############################################################### 
 33:  
 34: # khs26> changed the path to use the $AMBERHOME environment variable 
 35: amberhome = os.environ["AMBERHOME"] 
 36: path = os.path.join(amberhome, "dat/leap/lib") 
 37:  
 38: ######################### 
 39: ## some useful functions 
 40: ######################### 
 41:  
 42: def exchange_atoms(atom_type, a, aa, residue, dihedrals, currentAtomNumber): 
 43:   find_atom = a[aa.index(residue)].index(atom_type) 
 44:   atomNumber = find_atom+currentAtomNumber 
 45:   atomNumberIndex = atomNumber*3 
 46:   for j in range(len(dihedrals)): 
 47:     if (dihedrals[j][1]==str(atomNumberIndex)): 
 48:       d1 = dihedrals[j][0] 
 49:       d2 = dihedrals[j][1] 
 50:       dihedrals[j][0] = d2 
 51:       dihedrals[j][1] = d1 
 52:  
 53: def exchange_atoms_nt(atom_type, a, aa, residue, dihedrals): 
 54:   find_atom = a[aa.index(residue)].index(atom_type) 
 55:   for j in range(len(dihedrals)): 
 56:     if (dihedrals[j][1]==str(atomIndex[find_atom])): 
 57:       d1 = dihedrals[j][0] 
 58:       d2 = dihedrals[j][1] 
 59:       dihedrals[j][0] = d2 
 60:       dihedrals[j][1] = d1 
 61:  
 62: def exchange_atoms_arg(a, aa, residue, dihedrals, currentAtomNumber): 
 63:       ## IMPROPER responsible for trouble with NH2 group permutation:  
 64:       find_atom1 = a[aa.index(residue)].index('NE') 
 65:       atomNumber1 = find_atom1+currentAtomNumber 
 66:       atomNumberIndex1 = atomNumber1*3 
 67:       find_atom2 = a[aa.index(residue)].index('NH1') 
 68:       atomNumber2 = find_atom2+currentAtomNumber 
 69:       atomNumberIndex2 = atomNumber2*3 
 70:       find_atom3 = a[aa.index(residue)].index('CZ') 
 71:       atomNumber3 = find_atom3+currentAtomNumber 
 72:       atomNumberIndex3 = atomNumber3*3 
 73:       find_atom4 = a[aa.index(residue)].index('NH2') 
 74:       atomNumber4 = find_atom4+currentAtomNumber 
 75:       atomNumberIndex4 = atomNumber4*3 
 76:       for j in range(len(dihedrals)): 
 77:         if ((dihedrals[j][0]==str(atomNumberIndex1)) and (dihedrals[j][1]==str(atomNumberIndex2))): 
 78:           d0 = dihedrals[j][0] 
 79:           d1 = dihedrals[j][1] 
 80:           dihedrals[j][0] = d1 
 81:           dihedrals[j][1] = d0 
 82:  
 83: def exchange_atoms_ring1(a, aa, residue, dihedrals): 
 84:   find_atom1 = a[aa.index(residue)].index('CD1') 
 85:   atomNumber1 = find_atom1+currentAtomNumber 
 86:   atomNumberIndex1 = atomNumber1*3 
 87:   find_atom2 = a[aa.index(residue)].index('CD2') 
 88:   atomNumber2 = find_atom2+currentAtomNumber 
 89:   atomNumberIndex2 = atomNumber2*3 
 90:   for j in range(len(dihedrals)): 
 91:     if ((dihedrals[j][0]==str(atomNumberIndex1)) and (dihedrals[j][1]==str(atomNumberIndex2))): 
 92:       d0 = '-'+dihedrals[j][0] 
 93:       d1 = dihedrals[j][1] 
 94:       d3 = dihedrals[j][3][1:] 
 95:       dihedrals[j][0] = d1 
 96:       dihedrals[j][1] = d3 
 97:       dihedrals[j][3] = d0 
 98:  
 99: def exchange_atoms_ring2(a, aa, residue, dihedrals): 
100:   find_atom1 = a[aa.index(residue)].index('CG') 
101:   atomNumber1 = find_atom1+currentAtomNumber 
102:   atomNumberIndex1 = atomNumber1*3 
103:   find_atom2 = a[aa.index(residue)].index('CE2') 
104:   atomNumber2 = find_atom2+currentAtomNumber 
105:   atomNumberIndex2 = atomNumber2*3 
106:   find_atom3 = a[aa.index(residue)].index('CZ') 
107:   atomNumber3 = find_atom3+currentAtomNumber 
108:   atomNumberIndex3 = atomNumber3*3 
109:   find_atom4 = a[aa.index(residue)].index('CD2') 
110:   atomNumber4 = find_atom4+currentAtomNumber 
111:   atomNumberIndex4 = atomNumber4*3 
112:   find_atom5 = a[aa.index(residue)].index('CD1') 
113:   atomNumber5 = find_atom5+currentAtomNumber 
114:   atomNumberIndex5 = atomNumber5*3 
115:   find_atom6 = a[aa.index(residue)].index('CE1') 
116:   atomNumber6 = find_atom6+currentAtomNumber 
117:   atomNumberIndex6 = atomNumber6*3 
118: #  for j in range(len(dihedrals)): # this is ok 
119: #    if ((dihedrals[j][0]==str(atomNumberIndex1)) and (dihedrals[j][1]==str(atomNumberIndex2))): 
120:   for j in range(len(dihedrals)): 
121:     if ((dihedrals[j][0]==str(atomNumberIndex3)) and (dihedrals[j][1]==str(atomNumberIndex4))): 
122:       d1 = '-'+dihedrals[j][1] 
123:       d2 = dihedrals[j][3][1:] 
124:       dihedrals[j][1] = d2 
125:       dihedrals[j][3] = d1 
126:   for j in range(len(dihedrals)): 
127:     if ((dihedrals[j][0]==str(atomNumberIndex5)) and (dihedrals[j][1]==str(atomNumberIndex3))): 
128:       d1 = '-'+dihedrals[j][1] 
129:       d2 = dihedrals[j][3][1:] 
130:       dihedrals[j][1] = d2 
131:       dihedrals[j][3] = d1 
132:   for j in range(len(dihedrals)): 
133:     if ((dihedrals[j][0]==str(atomNumberIndex1)) and (dihedrals[j][1]==str(atomNumberIndex6))): 
134:       ## to compare IMPROPER before and after permutation 
135: ##test      a1 = (int(dihedrals[j][0])-currentAtomNumber)/3 
136: ##test      a2 = (int(dihedrals[j][1])-currentAtomNumber)/3 
137: ##test      a3 = (int(dihedrals[j][2][1:])-currentAtomNumber)/3 
138: ##test      a4 = (int(dihedrals[j][3][1:])-currentAtomNumber)/3 
139: ##test      print dihedrals[j], a[aa.index(residue)][a1], a[aa.index(residue)][a2], a[aa.index(residue)][a3], a[aa.index(residue)][a4] 
140:       d1 = '-'+dihedrals[j][0] 
141:       d2 = dihedrals[j][3][1:] 
142:       dihedrals[j][0] = d2 
143:       dihedrals[j][3] = d1 
144: ##test      a1 = (int(dihedrals[j][0])-currentAtomNumber)/3 
145: ##test      a2 = (int(dihedrals[j][1])-currentAtomNumber)/3 
146: ##test      a3 = (int(dihedrals[j][2][1:])-currentAtomNumber)/3 
147: ##test      a4 = (int(dihedrals[j][3][1:])-currentAtomNumber)/3 
148: ##test      print dihedrals[j], a[aa.index(residue)][a1], a[aa.index(residue)][a2], a[aa.index(residue)][a3], a[aa.index(residue)][a4] 
149:  
150: def exchange_atoms_ring3(a, aa, residue, dihedrals): 
151:   find_atom1 = a[aa.index(residue)].index('CE1') 
152:   atomNumber1 = find_atom1+currentAtomNumber 
153:   atomNumberIndex1 = atomNumber1*3 
154:   find_atom2 = a[aa.index(residue)].index('CE2') 
155:   atomNumber2 = find_atom2+currentAtomNumber 
156:   atomNumberIndex2 = atomNumber2*3 
157:   for j in range(len(dihedrals)): 
158:     if ((dihedrals[j][0]==str(atomNumberIndex1)) and (dihedrals[j][1]==str(atomNumberIndex2))): 
159:       d0 = '-'+dihedrals[j][0] 
160:       d1 = dihedrals[j][1] 
161:       d3 = dihedrals[j][3][1:] 
162:       dihedrals[j][0] = d1 
163:       dihedrals[j][1] = d3 
164:       dihedrals[j][3] = d0 
165:  
166: #################################### 
167: ## reading amino12.lib library 
168: #################################### 
169:  
170: print '\nDear user, please notice that only residues from the following libraries are taken into account:' 
171: print '  ions94.lib' 
172:  
173: print '  amino12.lib' 
174: aalib = open("%s/amino12.lib" % path).read() 
175: aa = string.split(aalib, "\n") 
176: q1 = aa.index("!!index array str") 
177: q2 = aa.index("!entry.ALA.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg") 
178:  
179: aaNames = [] # amino acid names 
180: aTypes = []  # atom types 
181: aNames = []  # atom names 
182:  
183: for i in range(q2-q1-1): 
184:   aaNames.append(aa[q1+1+i][2:5]) 
185:  
186: for i in range(len(aaNames)): 
187:   q1 = aa.index("!entry.%s.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg" % aaNames[i]) 
188:   q2 = aa.index("!entry.%s.unit.atomspertinfo table  str pname  str ptype  int ptypex  int pelmnt  dbl pchg" % aaNames[i]) 
189:   aT = [] 
190:   aN = [] 
191:   for j in range(q2-q1-1): 
192:     aT.append((string.split(aa[q1+1+j])[0]).replace('"','')) 
193:     aN.append((string.split(aa[q1+1+j])[1]).replace('"','')) 
194:   aTypes.append(aT) 
195:   aNames.append(aN) 
196:  
197: ###################################### 
198: ## reading aminont12.lib library 
199: ###################################### 
200:  
201: print '  aminont12.lib' 
202: aalib = open("%s/aminont12.lib" % path).read() 
203: aa = string.split(aalib, "\n") 
204: q1 = aa.index("!!index array str") 
205: q2 = aa.index("!entry.ACE.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg") 
206:  
207: aantNames = [] # N terminus amino acid names 
208: antTypes = []  # N terminus atom types 
209: antNames = []  # N terminus atom names 
210:  
211: for i in range(q2-q1-1): 
212:   aantNames.append((aa[q1+1+i].replace('"','')).replace(' ','')) 
213:  
214: for i in range(len(aantNames)): 
215:   q1 = aa.index("!entry.%s.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg" % aantNames[i]) 
216:   q2 = aa.index("!entry.%s.unit.atomspertinfo table  str pname  str ptype  int ptypex  int pelmnt  dbl pchg" % aantNames[i]) 
217:   aT = [] 
218:   aN = [] 
219:   for j in range(q2-q1-1): 
220:     aT.append((string.split(aa[q1+1+j])[0]).replace('"','')) 
221:     aN.append((string.split(aa[q1+1+j])[1]).replace('"','')) 
222:   antTypes.append(aT) 
223:   antNames.append(aN) 
224:  
225: ###################################### 
226: ## reading aminoct12.lib library 
227: ###################################### 
228:  
229: print '  aminoct12.lib' 
230: aalib = open("%s/aminoct12.lib" % path).read() 
231: aa = string.split(aalib, "\n") 
232: q1 = aa.index("!!index array str") 
233: q2 = aa.index("!entry.CALA.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg") 
234:  
235: aactNames = [] # C terminus amino acid names 
236: actTypes = []  # C terminus atom types 
237: actNames = []  # C terminus atom names 
238:  
239: for i in range(q2-q1-1): 
240:   aactNames.append((aa[q1+1+i].replace('"','')).replace(' ','')) 
241:  
242: for i in range(len(aactNames)): 
243:   q1 = aa.index("!entry.%s.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg" % aactNames[i]) 
244:   q2 = aa.index("!entry.%s.unit.atomspertinfo table  str pname  str ptype  int ptypex  int pelmnt  dbl pchg" % aactNames[i]) 
245:   aT = [] 
246:   aN = [] 
247:   for j in range(q2-q1-1): 
248:     aT.append((string.split(aa[q1+1+j])[0]).replace('"','')) 
249:     aN.append((string.split(aa[q1+1+j])[1]).replace('"','')) 
250:   actTypes.append(aT) 
251:   actNames.append(aN) 
252:  
253: ##################################### 
254: ## reading nucleic12.lib library 
255: ##################################### 
256:  
257: print '  nucleic12.lib' 
258: # check to see if the file is where expected or moved (as in AMBER14+) 
259: aalib = open("%s/nucleic12.lib" % path).read() 
260: aa = string.split(aalib, "\n") 
261: q1 = aa.index("!!index array str") 
262: q2 = aa.index("!entry.A.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg") 
263:  
264: nucNames = [] # nucleic names 
265: nucTypes = []  # nucleic atom types 
266: nucaNames = []  # nucleic atom names 
267:  
268: for i in range(q2-q1-1): 
269:   nucNames.append((aa[q1+1+i].replace('"','')).replace(' ','')) 
270:  
271: for i in range(len(nucNames)): 
272:   q1 = aa.index("!entry.%s.unit.atoms table  str name  str type  int typex  int resx  int flags  int seq  int elmnt  dbl chg" % nucNames[i]) 
273:   q2 = aa.index("!entry.%s.unit.atomspertinfo table  str pname  str ptype  int ptypex  int pelmnt  dbl pchg" % nucNames[i]) 
274:   aT = [] 
275:   aN = [] 
276:   for j in range(q2-q1-1): 
277:     aT.append((string.split(aa[q1+1+j])[0]).replace('"','')) 
278:     aN.append((string.split(aa[q1+1+j])[1]).replace('"','')) 
279:   nucTypes.append(aT) 
280:   nucaNames.append(aN) 
281:  
282: ################################# 
283: ## reading original prmtop file 
284: ################################# 
285:  
286: prmtop = open(sys.argv[1]).read() 
287: f = string.split(prmtop, "\n") 
288:  
289: q0 = f.index("%FLAG POINTERS                                                                  ") 
290: q1 = f.index("%FLAG ATOM_NAME                                                                 ") 
291: q2 = f.index("%FLAG CHARGE                                                                    ") 
292: q3 = f.index("%FLAG RESIDUE_LABEL                                                             ") 
293: q4 = f.index("%FLAG RESIDUE_POINTER                                                           ") 
294: q5 = f.index("%FLAG DIHEDRALS_INC_HYDROGEN                                                    ") 
295: q6 = f.index("%FLAG DIHEDRALS_WITHOUT_HYDROGEN                                                ") 
296: q7 = f.index("%FLAG EXCLUDED_ATOMS_LIST                                                       ") 
297:  
298: ## names of tables are related to names in prmtop file 
299:  
300: atomNumber = int(string.split(f[q0+2])[0]) 
301:  
302: atomName = [] 
303: residueLabel = [] 
304: dihedralsIncHydrogen = [] 
305: dihedralsWithoutHydrogen = [] 
306: atomIndex = [] 
307:  
308: an = 0 
309: line = 0 
310: while (an<atomNumber): 
311:   for j in range(20): 
312:     if (an<atomNumber): 
313:       an = an+1 
314:       atomName.append(f[q1+2+line][j*4:(j+1)*4].strip()) 
315:     else: 
316:       break 
317:   line = line+1 
318:  
319: for i in range(q4-q3-2): 
320:   for j in range((len(f[q3+2+i])+1)/4): 
321:     residueLabel.append(string.strip(f[q3+2+i][j*4:4*(j+1)])) 
322:  
323: caa = 0 
324: naa = 0 
325: for i in range(len(residueLabel)): 
326:   if (aactNames.count(residueLabel[i])>0): caa = caa+1 
327:   if (aantNames.count(residueLabel[i])>0): naa = naa+1 
328:  
329: if (caa==0): 
330:   print "\n-----------------------------------------------------------------------------" 
331:   print 'There is no C terminus amino acid in topology file!' 
332:   print 'If system does not contain amino acids - continue.' 
333:   print 'Otherwise please rename each C terminal residue by adding \'C\' to the name, e.g. ALA -> CALA' 
334:   print 'Remember to follow format of topology file!' 
335:   print "-----------------------------------------------------------------------------\n" 
336:   ## the only exception for C terminus amino acids is NME 
337:  
338: if (naa==0): 
339:   print "-----------------------------------------------------------------------------" 
340:   print 'There is no N terminus amino acid in topology file!' 
341:   print 'If system does not contain amino acids - continue.' 
342:   print 'Otherwise please rename each N terminal residue by adding \'N\' to the name, e.g. ALA -> NALA' 
343:   print 'Remember to follow format of topology file!' 
344:   print "-----------------------------------------------------------------------------" 
345:   ## the only exception for N terminus amino acids is ACE 
346:  
347: for i in range(q6-q5-2): 
348:   for j in range(len(string.split(f[q5+2+i]))/5): 
349:     dihedralsIncHydrogen.append(string.split(f[q5+2+i][j*40:40*(j+1)])) 
350:  
351: for i in range(q7-q6-2): 
352:   for j in range(len(string.split(f[q6+2+i]))/5): 
353:     dihedralsWithoutHydrogen.append(string.split(f[q6+2+i][j*40:40*(j+1)])) 
354:  
355: for i in range(len(atomName)): 
356:   atomIndex.append(i*3) 
357:  
358: ############################################################ 
359: ## groups of amino acids according to permutation behaviour 
360: ############################################################ 
361:  
362: ## group0 and group0n: nothing to do 
363: ## group1: problem only with C terminus COO- group: CA-O-C-OXT -> O-CA-C-OXT 
364:  
365: group0 = ['GLY','ALA','VAL','LEU','MET','ILE','SER','THR','CYS','LYS','HIP','HIE','HID'] 
366: group0n = ['NGLY','NALA','NVAL','NLEU','NMET','NILE','NSER','NTHR','NCYS','NLYS','NHIP','NHIE','NHID'] 
367: group1 = ['CGLY','CALA','CVAL','CLEU','CMET','CILE','CSER','CTHR','CCYS','CLYS','CPRO','CTRP','CHIP','CHIE','CHID','CHYP'] 
368:  
369: ################################################################# 
370: ## groups of nucleic residues according to permutation behaviour 
371: ################################################################# 
372:  
373: group2 = ['DA', 'DA3', 'DA5', 'DAN', 'A', 'A3', 'A5', 'AN'] 
374: group3 = ['DC', 'DC3', 'DC5', 'DCN', 'C', 'C3', 'C5', 'CN'] 
375: group4 = ['DG', 'DG3', 'DG5', 'DGN', 'G', 'G3', 'G5', 'GN'] 
376: group5 = ['DT', 'DT3', 'DT5', 'DTN', 'U', 'U3', 'U5', 'UN' , 'OHE' ] ## nothing to do 
377:  
378: ##################################################### 
379: ## groups of species without any IMPROPER to correct 
380: ##################################################### 
381:  
382: group7 = ['WAT', 'CIO', 'Cl-', 'Cs+', 'IB', 'K+', 'Li+', 'MG2', 'Na+', 'Rb+'] 
383:  
384: ## groupUknown: residue not counted in library 
385: groupUknown = [] 
386:  
387: for i in range(len(residueLabel)): 
388:   if ((aaNames.count(residueLabel[i])==0) and (aactNames.count(residueLabel[i])==0) and (aantNames.count(residueLabel[i])==0) and (nucNames.count(residueLabel[i])==0) and (group7.count(residueLabel[i])==0)): 
389:     groupUknown.append(residueLabel[i]) 
390:  
391: if (len(groupUknown)!=0):  
392:   print '\nThere are some residues missing in considered libraries:', groupUknown 
393:   print 'Program just skips them\n' 
394:  
395: currentAtomNumber = 0 
396:  
397: ###################################################################### 
398: ## main part - permutation of atom positions in appropriate IMPROPERs 
399: ###################################################################### 
400:  
401: for i in range(len(residueLabel)): 
402: #  print '----', i+1, '----', residueLabel[i] 
403:   if (group7.count(residueLabel[i])>0): continue 
404:  
405:   elif (groupUknown.count(residueLabel[i])>0): continue 
406:     ################################## 
407:     ## residue not counted in library 
408:     ################################## 
409:  
410:   elif (aantNames.count(residueLabel[i])>0): 
411:     ######################### 
412:     ## N terminus amino acid 
413:     ######################### 
414:  
415:     if (group0n.count(residueLabel[i])>0):  
416:       currentAtomNumber = currentAtomNumber+len(antTypes[aantNames.index(residueLabel[i])]) 
417:       continue 
418:     elif (residueLabel[i]=='NASN'): 
419:       exchange_atoms_nt('HD21', antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
420:     elif (residueLabel[i]=='NGLN'): 
421:       exchange_atoms_nt('HE21', antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
422:     elif (residueLabel[i]=='NARG'): 
423:       exchange_atoms_nt('HH11', antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
424:       exchange_atoms_nt('HH21', antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
425:       exchange_atoms_arg(antTypes, aantNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
426:     elif (residueLabel[i]=='NASP'): 
427:       exchange_atoms_nt('OD1', antTypes, aantNames, residueLabel[i], dihedralsWithoutHydrogen) 
428:     elif (residueLabel[i]=='NGLU'): 
429:       exchange_atoms_nt('OE1', antTypes, aantNames, residueLabel[i], dihedralsWithoutHydrogen) 
430:     elif (residueLabel[i]=='NPHE'): 
431:       exchange_atoms_ring1(antTypes, aantNames, residueLabel[i], dihedralsWithoutHydrogen) 
432:       exchange_atoms_ring2(antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
433:       exchange_atoms_ring3(antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
434:     elif (residueLabel[i]=='NTYR'): 
435:       exchange_atoms_ring1(antTypes, aantNames, residueLabel[i], dihedralsWithoutHydrogen) 
436:       exchange_atoms_ring2(antTypes, aantNames, residueLabel[i], dihedralsIncHydrogen) 
437:       exchange_atoms_ring3(antTypes, aantNames, residueLabel[i], dihedralsWithoutHydrogen) 
438: #        if (dihedralsWithoutHydrogen[j].count(str(atomNumberIndex1))>0): 
439: #          print '"""', dihedralsWithoutHydrogen[j], atomNumberIndex1, dihedralsWithoutHydrogen[j].count(str(atomNumberIndex1)) 
440: #          print '!!!', dihedralsWithoutHydrogen[j], atomNumberIndex1, dihedralsWithoutHydrogen[j].index(str(atomNumberIndex1)) 
441:  
442:     currentAtomNumber = currentAtomNumber+len(antTypes[aantNames.index(residueLabel[i])]) 
443:  
444:   elif (aactNames.count(residueLabel[i])>0): 
445:     ######################### 
446:     ## C terminus amino acid 
447:     ######################### 
448:  
449:     if (group1.count(residueLabel[i])>0): ## res belongs to group1 
450:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
451:     elif (residueLabel[i]=='CASP'): 
452:       exchange_atoms('OD1', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
453:     elif (residueLabel[i]=='CGLU'): 
454:       exchange_atoms('OE1', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
455:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
456:     elif (residueLabel[i]=='CGLN'): 
457:       exchange_atoms('HE21', actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
458:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
459:     elif (residueLabel[i]=='CASN'): 
460:       exchange_atoms('HD21', actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
461:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
462:     elif (residueLabel[i]=='CARG'): 
463:       exchange_atoms('HH11', actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
464:       exchange_atoms('HH21', actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
465:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
466:       exchange_atoms_arg(actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
467:     elif (residueLabel[i]=='CPHE'): 
468:       exchange_atoms_ring1(actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen) 
469:       exchange_atoms_ring2(actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen) 
470:       exchange_atoms_ring3(actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen) 
471:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
472:     elif (residueLabel[i]=='CTYR'): 
473:       exchange_atoms_ring1(actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen) 
474:       exchange_atoms_ring2(actTypes, aactNames, residueLabel[i], dihedralsIncHydrogen) 
475:       exchange_atoms('O', actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
476:       exchange_atoms_ring3(actTypes, aactNames, residueLabel[i], dihedralsWithoutHydrogen) 
477:  
478:     currentAtomNumber = currentAtomNumber+len(actTypes[aactNames.index(residueLabel[i])]) 
479:  
480:   elif (aaNames.count(residueLabel[i])>0): 
481:     ########################### 
482:     ## not terminal amino acid 
483:     ########################### 
484:  
485:     if (group0.count(residueLabel[i])>0):  
486:       currentAtomNumber = currentAtomNumber+len(aTypes[aaNames.index(residueLabel[i])]) 
487:       continue 
488:     elif (residueLabel[i]=='GLU'): 
489:       exchange_atoms('OE1', aTypes, aaNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
490:     elif (residueLabel[i]=='ASP'): 
491:       exchange_atoms('OD1', aTypes, aaNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
492:     elif (residueLabel[i]=='ASN'): 
493:       exchange_atoms('HD21', aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
494:     elif (residueLabel[i]=='GLN'): 
495:       exchange_atoms('HE21', aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
496:     elif (residueLabel[i]=='ARG'): 
497:       exchange_atoms('HH11', aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
498:       exchange_atoms('HH21', aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
499:       exchange_atoms_arg(aTypes, aaNames, residueLabel[i], dihedralsWithoutHydrogen, currentAtomNumber) 
500:     elif (residueLabel[i]=='PHE'): 
501:       exchange_atoms_ring1(aTypes, aaNames, residueLabel[i], dihedralsWithoutHydrogen) 
502:       exchange_atoms_ring2(aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen) 
503:       exchange_atoms_ring3(aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen) 
504:     elif (residueLabel[i]=='TYR'): 
505:       exchange_atoms_ring1(aTypes, aaNames, residueLabel[i], dihedralsWithoutHydrogen) 
506:       exchange_atoms_ring2(aTypes, aaNames, residueLabel[i], dihedralsIncHydrogen) 
507:       exchange_atoms_ring3(aTypes, aaNames, residueLabel[i], dihedralsWithoutHydrogen) 
508:  
509:     currentAtomNumber = currentAtomNumber+len(aTypes[aaNames.index(residueLabel[i])]) 
510:  
511:   elif (nucNames.count(residueLabel[i])>0): 
512:     ################### 
513:     ## nucleic residue 
514:     ################### 
515:  
516:     if (group5.count(residueLabel[i])>0): 
517:       currentAtomNumber = currentAtomNumber+len(nucTypes[nucNames.index(residueLabel[i])]) 
518:       continue 
519:     elif (group2.count(residueLabel[i])>0): 
520:       exchange_atoms('H61', nucTypes, nucNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
521:     elif (group3.count(residueLabel[i])>0): 
522:       exchange_atoms('H41', nucTypes, nucNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
523:     elif (group4.count(residueLabel[i])>0): 
524:       exchange_atoms('H21', nucTypes, nucNames, residueLabel[i], dihedralsIncHydrogen, currentAtomNumber) 
525:  
526:     currentAtomNumber = currentAtomNumber+len(nucTypes[nucNames.index(residueLabel[i])]) 
527:  
528:   else: 
529:     print 'Something strange happened... residue %d is neither in libraries and nor out of libraries' % residueLabel[i] 
530:     sys.exit() 
531:  
532: ################################## 
533: ## preparation of new prmtop file 
534: ################################## 
535:  
536: newprmtop = open(sys.argv[2],'w') 
537:  
538: for i in range(q5+2): 
539:   newprmtop.write("%s\n" % f[i]) 
540:  
541: an = 0 
542: line = 0 
543: while (an<len(dihedralsIncHydrogen)): 
544:   for j in range(2): 
545:     if (an<len(dihedralsIncHydrogen)): 
546:       for k in range(5): 
547:         newprmtop.write("%8s" % (dihedralsIncHydrogen[an][k])) 
548:       an = an+1 
549:     else: 
550:       break 
551:   newprmtop.write("\n") 
552:   line = line+1 
553:  
554: for i in range(2): 
555:   newprmtop.write("%s\n" % f[q6+i]) 
556:  
557: an = 0 
558: line = 0 
559: while (an<len(dihedralsWithoutHydrogen)): 
560:   for j in range(2): 
561:     if (an<len(dihedralsWithoutHydrogen)): 
562:       for k in range(5): 
563:         newprmtop.write("%8s" % (dihedralsWithoutHydrogen[an][k])) 
564:       an = an+1 
565:     else: 
566:       break 
567:   newprmtop.write("\n") 
568:   line = line+1 
569:  
570: for i in range(len(f)-q7-1): 
571:   newprmtop.write("%s\n" % f[q7+i]) 
572:  
573: newprmtop.close() 
574:  
575: print "\n-----------------------------------------------------------------" 
576: print 'If you added \'N\' and \'C\' to the names of terminal amino acids' 
577: print 'in the topology file, please remove them now to avoid problems' 
578: print 'with programs that produce PDB files' 
579: print "-----------------------------------------------------------------\n" 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0