scripts/make_gromos_rtp.py

   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
  22 from string import atoi, split
  23
  24 GROMOS_FILE   = 'mtb43a1.dat'
  25 ATOMTYPE_FILE = 'ffG43a1.atp'
  26
  27 START  = '# RNME\012'
  28 STOP   = '#RNMES\012'
  29
  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'
  35 NB     = '#  NB\012'
  36 NBA    = '# NBA\012'
  37 NIDA   = '# NIDA\012'
  38 NDA    = '# NDA\012'
  39
  40
  41 #### functions
  42 def 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
  51 def 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
  92 class 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
 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])
 239
 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)
 243
 244 start = 0; stop = 92
 245
 246 " The rtp header "
 247 print "[ bondedtypes ]"
 248 print "; bonds  angles  dihedrals  impropers"
 249 print "    2       2          1          2"
 250
 251
 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
 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])