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