GIT repositories md_scripts / master python / translate_backwards_map.py
master

Tree @master (Download .tar.gz)

translate_backwards_map.py @masterraw · history · blame

#!/usr/bin/python
import sys
import os
import argparse

# default input/output files
#--------------------------
oldmapfile = 'chol.charmm36.map'
newmapfile = 'chol.opls.map'
itpfile = 'chol-opls.itp'
transfile = 'chol-charmm2opls.trans'


# commandline parser
#--------------------------------
parser = argparse.ArgumentParser()
parser.add_argument("-i", help='old map file')
parser.add_argument("-p", help='new topology file')
parser.add_argument("-o", help='output file')
parser.add_argument("-t", help='translation file')
args = parser.parse_args()

if args.i:
  oldmapfile = args.i
if args.p:
  itpfile = args.p
if args.o:
  newmapfile = args.o
if args.t:
  transfile = args.t

print(newmapfile)
if os.path.exists(newmapfile):
  print(newmapfile + ' already exists.')
  exit(0)

# load translation
#--------------------------
transstream = open(transfile, 'r')
newtoold = {}
oldtonew = {}

for line in transstream:

  if line.startswith(';'):
    continue

  columns = line.split()
  if len(columns) < 2:
    continue

  newtoold[columns[1]] = columns[0]
  oldtonew[columns[0]] = columns[1]


transstream.close()



# load atom order from itp
#--------------------------
itpstream = open(itpfile, 'r')
order = []

# skip stuf  before [ atoms ] section
while itpstream.readline().find('[ atoms ]') < 0:
  continue

# skip header
itpstream.readline()

# write [ atoms ] section to tempitp
while True:
  line = itpstream.readline()
  if line.startswith('['):
    break
  if line.startswith(';'):
    continue

  columns = line.split()
  if len(columns) < 5:
    continue

  atom = columns[4]
  order.append(atom)

itpstream.close()



# load cg-beads from old map file
#--------------------------
oldmapstream = open(oldmapfile, 'r')
beads = {}

# skip stuf  before [ atoms ] section
while oldmapstream.readline().find('[ atoms ]') < 0:
  continue

# save beads
while True:
  line = oldmapstream.readline()
  if line.startswith('['):
    break
  if line.startswith(';'):
    continue

  columns = line.split()
  if len(columns) < 3:
    continue

  atom = columns[1]
  beadset = columns[2:]
  beads[atom] = beadset

oldmapstream.close()




# write new file
#--------------------------
oldmapstream = open(oldmapfile, 'r')
newmapstream = open(newmapfile, 'w')

# Until [ atoms ] section, just copy everything
while True:
  line = oldmapstream.readline()
  newmapstream.write(line)
  if line.find('atoms') > 0:
    break
  if len(line) < 1:
    print '[ atoms ] section not found.'
    sys.exit(0)


# write [ atoms ] section
num = 1
for atom in order:
   
  if atom in newtoold:
    oldatom = newtoold[atom]
    beadstr = ' '.join(beads[oldatom])
  else:
    beadstr = ''

  atomstr = "{0:5}{1:6}".format(str(num), atom)
  newline = atomstr+ '   ' + beadstr + '\n'
  newmapstream.write(newline)
  num = num + 1


# find next section after [ atoms ]
while True:
  line = oldmapstream.readline()
  if line.startswith(';') or line.startswith('[') or line.isspace():
    newmapstream.write(line)
    break

# Switch atom names in the last sections
while True:

  line = oldmapstream.readline()

  if line.startswith(';') or line.startswith('[') or line.isspace():
    newmapstream.write(line)
    continue

  if len(line) < 1:
    print 'end of file'
    break

  newline=''
  oldatoms = line.split()
  for a in oldatoms:
    newline = newline + oldtonew[a] + ' '

  newmapstream.write(newline + '\n')


oldmapstream.close()
newmapstream.close()