1 #!/usr/freeware/bin/python
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.
11 # The output file requires some manual modifications:
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
17 # author: Klaas Dijkstra, Univ. Groningen
19 # minor modifications by Berk Hess
22 from string import atoi, split
24 GROMOS_FILE = 'mtb43a1.dat'
25 ATOMTYPE_FILE = 'ffG43a1.atp'
30 NATOM = '# NMAT,NLIN\012'
31 EXCLUS = '#ATOM MAE MSAE\012'
32 ATOM1 = '#ATOM ANM IACM MASS CGMICGM MAE MSAE\012'
33 ATOM2 = '#ATOM ANM IACM MASS CGM ICGM MAE MSAE\012'
34 ATOMtr = '#ATOM ANM IACM MASS CGMICGM\012'
42 def translate(t): # translate number into atom code
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]
51 def findbonds(atomnumber):
52 atom = eval(f.atoms[atomnumber][0])
64 bond.append(eval(a),eval(b))
65 bond.append(eval(b),eval(a))
68 if i[0] == atom: onebond.append(i[1])
72 if i[0] == j and atom < i[1]: twobond.append(i[1])
91 #### general class to unravel a gromos topology file
93 def __init__(self, filename):
95 self.G = open(filename).readlines()
96 self.nlines = len(self.G)
99 " Find all indexes of lines with string '# RNME\012'and '#RNMES\012' "
101 for i in range(self.nlines):
102 if self.G[i] == STOP:
105 if self.G[i] == START:
109 " Returns list of residues "
111 length = len(self.ind)
112 for i in range(length):
114 start = self.ind[i] + 1
115 try: stop = self.ind[i+1] - 2
117 for j in range(start,stop):
118 res.append(self.G[j])
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:
126 print >> sys.stderr "Could not find string",string,"in list of length",len(list)
129 #--------------------------#
130 # unravel gromos list
132 def getRESname(self, res):
133 " Get the name of the current residue "
134 self.residue = split(res[0])
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)
143 def getEXCLUS(self, res):
144 " Get preceding exclusions "
146 ind = self.gets(res, EXCLUS)
147 for i in range(self.nlin):
148 self.exclus.append(split(res[ind+i+1]))
150 def getATOM(self, res):
153 try: ind = self.gets(res, ATOM1)
154 except: ind = self.gets(res, ATOM2)
157 while cntr < self.natom - self.nlin:
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):
164 line1 = split(res[ind+i]) # overflow line
167 self.atoms.append(line)
170 def getATOMtr(self, res):
171 " Get trailing atoms"
173 try: ind = self.gets(res, ATOMtr)
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]))
179 def getBOND(self, res):
182 ind = self.gets(res, NB)
183 self.nb = atoi(split(res[ind+1])[0])
185 for i in range(self.nb):
186 line = split(res[ind+i+j+3])
188 line = split(res[ind+i+j+4])
190 self.bonds.append(line)
192 def getNBA(self, res):
195 ind = self.gets(res, NBA)
196 self.nba = atoi(split(res[ind+1])[0])
198 for i in range(self.nba):
199 line = split(res[ind+i+j+3])
201 line = split(res[ind+i+j+4])
205 def getNIDA(self, res):
206 " Get improper dihedrals"
208 ind = self.gets(res, NIDA)
209 self.nida = atoi(split(res[ind+1])[0])
211 for i in range(self.nida):
212 line = split(res[ind+i+j+3])
214 line = split(res[ind+i+j+4])
216 self.ida.append(line)
218 def getNDA(self, res):
221 ind = self.gets(res, NDA)
223 self.nda = atoi(split(res[ind+1])[0])
224 for i in range(self.nda):
225 line = split(res[ind+i+j+3])
227 line = split(res[ind+i+j+4])
232 #-----------------------------#
235 typ = open(ATOMTYPE_FILE) # translate numbers to atoms
236 typelines = typ.readlines()
237 for i in range(len(typelines)):
238 typelines[i]=split(typelines[i])
240 f=Cin(GROMOS_FILE) # bind class instance
241 f.index() # mark all residues (f.ind)
242 f.mkres() # put all residues into list (f.all)
247 print "[ bondedtypes ]"
248 print "; bonds angles dihedrals impropers"
252 for 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
256 if f.nlin != 0: # 0 for a separate molecule
257 f.getEXCLUS(f.all[resnum]) # number of exclusions
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
268 # output to Gromacs format
269 #-------------------------#
274 print "[",f.residue[0],"]"
277 for j in range(f.natom - f.nlin):
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])
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])
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:
306 print "%5s %5s gb_%-5s" % \
307 (t1, t2, f.bonds[j][2])
312 for j in range(f.natom - f.nlin):
316 if i in bbb: bbb.remove(i)
318 " Ignore special exclusions "
322 t2 = atoi(f.atoms[j][0])
324 if ne == 0: print " [ exclusions ]\n; ai aj"
325 print "%5s %5s" % (t2,t1)
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:
340 print "%5s %5s %5s ga_%-5s" % \
341 (t1,t2,t3,f.ba[j][3])
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:
357 print "%5s %5s %5s %5s gi_%-5s" % \
358 (t1,t2,t3,t4,f.ida[j][4])
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:
374 print "%5s %5s %5s %5s gd_%-5s" % \
375 (t1,t2,t3,t4,f.da[j][4])