Conversion file from mol2 atom types (generated from Chimera software) to Amber atom types

mol2_file_conversion.py

#!/usr/bin/env python
#
#
####   DATE        : 08/July/2018
####   Author      : Asif Iqbal
####   Script      : For changing mol2 variables
####   USE         : ./script.py filename.mol2
####   Edit mol2 file to remove empty line generated from chimera
#      coordinate counting should start form the 7th ROW!!!!!
#
#
import subprocess
from sys import argv
import linecache
from termcolor import colored, cprint
import numpy
import scipy, random
import math, time
import os,re, sys
import string
import fileinput
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.size
rank = comm.rank
status = MPI.Status()

#############################################################################

print "***************************************************************************"
print colored("  Script to change variable from mol2 file (generated with Chimera) with",'green')
print colored("  new variable (MCPB Amber code).",'green')
print
print colored("  Edit mol2 file to remove empty line generated from chimera.", 'yellow')
print colored("  Coordinate counting should start form the 7th ROW!!!!! ", 'yellow')
print
print colored("  The order of coordinates show be Fe - NN ... NN ... NN ...  ", 'yellow')
print colored("  Then all TFSI- and then all alkyl chains. This is necessary ", 'yellow')
print colored("  Since NN atoms have a cyclic period of i mod 6 == 0  ", 'yellow')
print "***************************************************************************"

##################################################################################
def loop(res_name, N_AT, atom_name, atom_type, no_atm):
	C=0
	pp=[]
	name = MPI.Get_processor_name()
	comm.send(None, dest=0, tag=0)
	task = comm.recv(source=0, tag=MPI.ANY_TAG, status=status)
	tag = status.Get_tag()
	for i in range(no_atm):
		source={}
		if atom_type[i] == 'nb':
			source['AT'] = atom_type[i]
			source['AN'] = atom_name[i]
			pp.append(source)
	for j in range(N_AT):
		pp[j]["AT"] =  'Y'+str(C+1)
		C+= 1
		if C % 6 == 0:
			C = 0
	
	for j in range(no_atm):
#		res_name[j] = 'POL'
		if atom_type[j] == 'n2':
			atom_type[j] = 'ne'
		
	for k in range(no_atm):
		if atom_type[k] == 'nb':
			atom_type[k] = pp[C]['AT']
			C+=1
			if C % N_AT == 0:
				break
		
	for t in range(no_atm):
		if atom_type[t] == 'sy':
			if t%2 != 0:
				atom_type[t] = 's6'
			else:
				atom_type[t] = 'sy'
    #        print atom_type[t]
    #print atom_type
	return atom_type

##################### Adding charges on the system ####################

def sys_charge(charge, atom_name, atom_type, no_atm):
    MMM=[]
    sum_chg=0
    chg_Fe=0
    counter=0
    for i in range(no_atm):
        source = {}
        source['AT_NUM'] = atom_name[i]
        source['AT_TYP'] = atom_type[i]
        source['AT_CHR'] = float(charge[i])
        MMM.append(source)
        if counter == no_atm:
            break

    for h in range(no_atm):
        sum_chg += MMM[h]['AT_CHR']
#        print MMM[h]
    print colored( "    Sum of the charge on the system is:  %6f  " %sum_chg, 'red')

    for h in range(no_atm):
        if MMM[h]['AT_TYP'] == 'FE':
            chg_Fe += MMM[h]['AT_CHR']
    print colored( "    Sum of the charge on the Fe ion is:  %6f  " %chg_Fe, 'red')

######################## Read each line from the file ##########################
def main_Engine():
    #script,filename = argv
    file = open(sys.argv[1], 'r')
    contents = file.readlines()
    file.close()
    #print contents

    ##################### Reading number of coordinates ######################

    #thirdline = linecache.getline('contents', 3)
    #bbb = thirdline.split( )
    bbb = contents[2]
    bbb = bbb.strip('\n')
    bbb = bbb.split(  )
#    print (bbb)
    no_atm = int(bbb[0])
    no_bds = int(bbb[1])
    print ("    Number of atoms in the file ", no_atm)
    print ("    Number of bonds in the file ", no_bds)

    ################ Variable for extacting data ######################

    atom_number    =[]
    atom_name    =[]
    res_number  =[]
    res_name    =[]
    x            =[]
    y            =[]
    z            =[]
    atom_type    =[]
    charge       =[]
    temp         =[]
    print colored('    Reading File ....', 'green')

    ######################### read from specific line #######################
    for l in range(6, no_atm+6):
		items = contents[l].split( )
		atom_number.append(int(items[0]))
		atom_name.append(items[1])
		x.append(float(items[2]))
		y.append(float(items[3]))
		z.append(float(items[4]))
		atom_type.append(items[5])
		res_number.append(int(items[6]))
		res_name.append(items[7])
		charge.append(float(items[8]))	
    #print (atom_type[0:10])
    N_AT = atom_type.count('nb')
    FE_AT = atom_type.count('FE')
    N2_AT = atom_type.count('n2')
    sy_AT = atom_type.count('sy')
	
    print ("    # of nb atomtypes corresponding to Ligands N atoms    : ", N_AT)
    print ("    # of FE atoms           : ", FE_AT, ': which means' , FE_AT, 'monomers units in a polymer')
    print ("    # of n2 atomtypes corresponding to TFSI- N atoms      : ", N2_AT)
    print ("    # of sy atomtypes corresponding to TFSI- S atoms      : ", sy_AT)

    ########## Here we replaces the atom types with MCPB variable for generating
    ########## Topology files. xLEaP will understand this variable.
    print
    loop(res_name, N_AT, atom_name, atom_type, no_atm)
    sys_charge(charge, atom_name, atom_type, no_atm)
    print colored('    Writing to a File ....', 'red', attrs=['reverse', 'blink'])
    print
    ################################# Writing TO A FILE ############################
	
    f1 = open('file_output.mol2','w')
    f1.write(''.join(contents[:6]))

    for p in range(0, no_atm):
        f1.write("%6d %5s %15f %10f %10f  %5s %6d %5s %10f\n" % (atom_number[p], atom_name[p], x[p], y[p], z[p], atom_type[p], res_number[p], res_name[p], charge[p]))
    f1.write(''.join(contents[no_atm+6:]))
#    subprocess.call(["sed -i -e 's/ar/1/g' file_output.mol2"], shell=True)
    f1.close()
  #################################################################################


if __name__ == '__main__':

    main_Engine()