1#!/usr/bin/env python2 2 3# usage: make_gromos_rtp.py > ffG43a1.rtp 4# this script tries to make a residue topology database for GROMACS from 5# the GROMOS version of this file. It needs ffG43a1.atp to fill in the 6# atom types for the atom integer codes. It converts until it finds the string 7# "#RNMES", which indicates the start of the solvent part 8# of mtb43a1.dat. Solvents are treated differently in GROMACS 9# The vaiables GROMOS_FILE and ATOMTYPE_FILE has to be specified as required. 10 11# The output file requires some manual modifications: 12# add ACE and NH2 13# add 4 dihedrals with parameters 0 0 2 in HEME 14# remove I in front of all atoms names in SO42- 15# make a copy of H2O called HOH 16 17# author: Klaas Dijkstra, Univ. Groningen 18# Oct. 2000 19# minor modifications by Berk Hess 20 21 22from string import atoi, split 23 24GROMOS_FILE = 'mtb43a1.dat' 25ATOMTYPE_FILE = 'ffG43a1.atp' 26 27START = '# RNME\012' 28STOP = '#RNMES\012' 29 30NATOM = '# NMAT,NLIN\012' 31EXCLUS = '#ATOM MAE MSAE\012' 32ATOM1 = '#ATOM ANM IACM MASS CGMICGM MAE MSAE\012' 33ATOM2 = '#ATOM ANM IACM MASS CGM ICGM MAE MSAE\012' 34ATOMtr = '#ATOM ANM IACM MASS CGMICGM\012' 35NB = '# NB\012' 36NBA = '# NBA\012' 37NIDA = '# NIDA\012' 38NDA = '# NDA\012' 39 40 41#### functions 42def translate(t): # translate number into atom code 43 if t == 0 : t = '-O' 44 elif t == -1 : t = '-C' 45 elif t == -2 : t = '-CA' 46 elif t == f.natom+1 : t = '+N' 47 else : t = f.atoms[t-1][1] 48 return t 49 50 51def findbonds(atomnumber): 52 atom = eval(f.atoms[atomnumber][0]) 53 54 onebond = [] 55 twobond = [] 56 ext = [] 57 bond = [] 58 ind = [] 59 excl = [] 60 61 62 for i in f.bonds: 63 a, b = i[0:2] 64 bond.append(eval(a),eval(b)) 65 bond.append(eval(b),eval(a)) 66 67 for i in bond: 68 if i[0] == atom: onebond.append(i[1]) 69 70 for i in bond: 71 for j in onebond: 72 if i[0] == j and atom < i[1]: twobond.append(i[1]) 73 74 ext = onebond 75 ext.extend(twobond) 76 77 for i in ext: 78 if i<atom: 79 ind.append(i) 80 for i in ind: 81 ext.remove(i) 82 83 ext.sort() 84 85 for i in ext: 86 excl.append(`i`) 87 88 return excl 89 90 91#### general class to unravel a gromos topology file 92class Cin: 93 def __init__(self, filename): 94 f=open(filename) 95 self.G = open(filename).readlines() 96 self.nlines = len(self.G) 97 98 def index(self): 99 " Find all indexes of lines with string '# RNME\012'and '#RNMES\012' " 100 self.ind = [] 101 for i in range(self.nlines): 102 if self.G[i] == STOP: 103 self.ind.append(i) 104 return 105 if self.G[i] == START: 106 self.ind.append(i) 107 108 def mkres(self): 109 " Returns list of residues " 110 self.all = [] 111 length = len(self.ind) 112 for i in range(length): 113 res = [] 114 start = self.ind[i] + 1 115 try: stop = self.ind[i+1] - 2 116 except: return 117 for j in range(start,stop): 118 res.append(self.G[j]) 119 self.all.append(res) 120 121 def gets(self,list,string): 122 " Returns linenumber of string in list" 123 for i in range(len(list)): 124 if list[i] == string: 125 return i 126 print >> sys.stderr, "Could not find string", string, "in list of length", len(list) 127 return -1 128 129#--------------------------# 130# unravel gromos list 131 132 def getRESname(self, res): 133 " Get the name of the current residue " 134 self.residue = split(res[0]) 135 136 def getNATOM(self, res): 137 " Get number of atoms, number of preceding excl. " 138 ind = self.gets(res, NATOM) 139 self.natom, self.nlin = split(res[ind+1]) 140 self.natom = atoi(self.natom) 141 self.nlin = atoi(self.nlin) 142 143 def getEXCLUS(self, res): 144 " Get preceding exclusions " 145 self.exclus = [] 146 ind = self.gets(res, EXCLUS) 147 for i in range(self.nlin): 148 self.exclus.append(split(res[ind+i+1])) 149 150 def getATOM(self, res): 151 " Get atoms" 152 self.atoms = [] 153 try: ind = self.gets(res, ATOM1) 154 except: ind = self.gets(res, ATOM2) 155 i = 0 156 cntr = 0 157 while cntr < self.natom - self.nlin: 158 i = i + 1 159 line = split(res[ind+i]) # get next line 160 noflo = (atoi(line[6])-1)/8 # if MAE > 8 cont. on next line 161 if noflo < 0: noflo = 0 162 for j in range(noflo): 163 i = i + 1 164 line1 = split(res[ind+i]) # overflow line 165 "print line1" 166 line = line + line1 167 self.atoms.append(line) 168 cntr = cntr + 1 169 170 def getATOMtr(self, res): 171 " Get trailing atoms" 172 self.atomtr = [] 173 try: ind = self.gets(res, ATOMtr) 174 except: return 175 for i in range(self.nlin): 176 self.atomtr.append(split(res[ind+i+1])) 177 self.atoms.append(split(res[ind+i+1])) 178 179 def getBOND(self, res): 180 " Get bonds" 181 self.bonds = [] 182 ind = self.gets(res, NB) 183 self.nb = atoi(split(res[ind+1])[0]) 184 j = 0 185 for i in range(self.nb): 186 line = split(res[ind+i+j+3]) 187 if line[0] == '#': 188 line = split(res[ind+i+j+4]) 189 j = j+1 190 self.bonds.append(line) 191 192 def getNBA(self, res): 193 " Get bond angles" 194 self.ba = [] 195 ind = self.gets(res, NBA) 196 self.nba = atoi(split(res[ind+1])[0]) 197 j = 0 198 for i in range(self.nba): 199 line = split(res[ind+i+j+3]) 200 if line[0] == '#': 201 line = split(res[ind+i+j+4]) 202 j = j + 1 203 self.ba.append(line) 204 205 def getNIDA(self, res): 206 " Get improper dihedrals" 207 self.ida = [] 208 ind = self.gets(res, NIDA) 209 self.nida = atoi(split(res[ind+1])[0]) 210 j = 0 211 for i in range(self.nida): 212 line = split(res[ind+i+j+3]) 213 if line[0] == '#': 214 line = split(res[ind+i+j+4]) 215 j = j + 1 216 self.ida.append(line) 217 218 def getNDA(self, res): 219 " Get dihedrals" 220 self.da = [] 221 ind = self.gets(res, NDA) 222 j = 0 223 self.nda = atoi(split(res[ind+1])[0]) 224 for i in range(self.nda): 225 line = split(res[ind+i+j+3]) 226 if line[0] == '#': 227 line = split(res[ind+i+j+4]) 228 j = j + 1 229 self.da.append(line) 230 231 232#-----------------------------# 233# main program 234 235typ = open(ATOMTYPE_FILE) # translate numbers to atoms 236typelines = typ.readlines() 237for i in range(len(typelines)): 238 typelines[i]=split(typelines[i]) 239 240f=Cin(GROMOS_FILE) # bind class instance 241f.index() # mark all residues (f.ind) 242f.mkres() # put all residues into list (f.all) 243 244start = 0; stop = 92 245 246" The rtp header " 247print "[ bondedtypes ]" 248print "; bonds angles dihedrals impropers" 249print " 2 2 1 2" 250 251 252for resnum in range(start,stop): # loop through all residues 253 f.getRESname(f.all[resnum]) # residue name 254 f.getNATOM (f.all[resnum]) # number of atoms 255 256 if f.nlin != 0: # 0 for a separate molecule 257 f.getEXCLUS(f.all[resnum]) # number of exclusions 258 259 f.getATOM (f.all[resnum]) # atoms => f.atoms 260 f.getATOMtr (f.all[resnum]) # trailing atoms => f.atomtr 261 f.getBOND (f.all[resnum]) # bonds => f.bonds 262 f.getNBA (f.all[resnum]) # bond angles => f.ba 263 f.getNIDA (f.all[resnum]) # improper dihedrals => f.ida 264 f.getNDA (f.all[resnum]) # dihedrals => f.da 265 266 267 268 # output to Gromacs format 269 #-------------------------# 270 271##### 272 # atoms 273 print "" 274 print "[",f.residue[0],"]" 275 print " [ atoms ]" 276 chargegroup = 0 277 for j in range(f.natom - f.nlin): 278 try: 279 atomtype = atoi(f.atoms [j][2]) - 1 280 atomfield = typelines[atomtype][0] 281 print "%5s %5s %11s %5s" % \ 282 (f.atoms [j][1],atomfield,f.atoms [j][4],chargegroup) 283 chargegroup = chargegroup + atoi(f.atoms[j][5]) 284 except: 285 print j 286 287##### 288 # trailing atoms 289 for j in range(f.nlin): 290 atomtype = atoi(f.atomtr [j][2]) - 1 291 atomfield = typelines[atomtype][0] 292 print "%5s %5s %11s %5s" % \ 293 (f.atomtr[j][1],atomfield,f.atomtr[j][4][:-2],chargegroup) 294 chargegroup = chargegroup + atoi(f.atomtr[j][5]) 295 296##### 297 # bonds 298 print " [ bonds ]" 299 for j in range(f.nb): 300 t1 = atoi(f.bonds [j][0]) 301 t2 = atoi(f.bonds [j][1]) 302 " Only special bonds go to 0 or less " 303 if t1 >= 1 and t2 >= 1: 304 t1 = translate(t1) 305 t2 = translate(t2) 306 print "%5s %5s gb_%-5s" % \ 307 (t1, t2, f.bonds[j][2]) 308 309##### 310 # exclusions 311 ne = 0 312 for j in range(f.natom - f.nlin): 313 aaa = findbonds(j) 314 bbb = f.atoms[j][7:] 315 for i in aaa: 316 if i in bbb: bbb.remove(i) 317 for i in bbb: 318 " Ignore special exclusions " 319 t1 = atoi(i) 320 if t1 >= 0: 321 t1 = translate(t1) 322 t2 = atoi(f.atoms[j][0]) 323 t2 = translate(t2) 324 if ne == 0: print " [ exclusions ]\n; ai aj" 325 print "%5s %5s" % (t2,t1) 326 ne = ne + 1 327 328##### 329 # angles 330 print " [ angles ]" 331 print "; ai aj ak gromos type" 332 for j in range(f.nba): 333 t1 = atoi(f.ba [j][0]) 334 t2 = atoi(f.ba [j][1]) 335 t3 = atoi(f.ba [j][2]) 336 if t1 >= -2 and t2 >= -2 and t3 >= -2: 337 t1 = translate(t1) 338 t2 = translate(t2) 339 t3 = translate(t3) 340 print "%5s %5s %5s ga_%-5s" % \ 341 (t1,t2,t3,f.ba[j][3]) 342 343##### 344 # improper dihedrals 345 print " [ impropers ]" 346 print "; ai aj ak al gromos type" 347 for j in range(f.nida): 348 t1 = atoi(f.ida [j][0]) 349 t2 = atoi(f.ida [j][1]) 350 t3 = atoi(f.ida [j][2]) 351 t4 = atoi(f.ida [j][3]) 352 if t1 >= -2 and t2 >= -2 and t3 >= -2 and t4 >= -2: 353 t1 = translate(t1) 354 t2 = translate(t2) 355 t3 = translate(t3) 356 t4 = translate(t4) 357 print "%5s %5s %5s %5s gi_%-5s" % \ 358 (t1,t2,t3,t4,f.ida[j][4]) 359 360##### 361 # dihedrals 362 print " [ dihedrals ]" 363 print "; ai aj ak al gromos type" 364 for j in range(f.nda): 365 t1 = atoi(f.da [j][0]) 366 t2 = atoi(f.da [j][1]) 367 t3 = atoi(f.da [j][2]) 368 t4 = atoi(f.da [j][3]) 369 if t1 >= -2 and t2 >= -2 and t3 >= -2 and t4 >= -2: 370 t1 = translate(t1) 371 t2 = translate(t2) 372 t3 = translate(t3) 373 t4 = translate(t4) 374 print "%5s %5s %5s %5s gd_%-5s" % \ 375 (t1,t2,t3,t4,f.da[j][4]) 376