Conversion script for Crafting LIPID topologies

groITP2namdTOP.py

import pandas as pd
import numpy as np
import sys
def Get_topology(name_itp):
    itp = open(name_itp, 'r')
    itp_lines = []
    for line in itp:
        line.strip()
        line = line.split(';')[0]
        if line:itp_lines.append(line.strip())
    itp.close()
    
    imp_data = {
        'moleculetype': 0,
        'atoms': 0,
        'bonds': 0,
        'angles': 0,
        'prop_dih': 0,
        'pairs': 0,
        'impr_dih': 0,
    }
    for n in range(len(itp_lines)):
        if 'moleculetype' in itp_lines[n]: imp_data['moleculetype'] = n 
        elif 'atoms' in itp_lines[n]: imp_data['atoms'] =n
        elif 'bonds' in itp_lines[n]: imp_data['bonds'] =n
        elif 'angles' in itp_lines[n]: imp_data['angles'] =n
        elif ('dihedrals' in itp_lines[n]) and (imp_data['prop_dih']==0): imp_data['prop_dih'] =n
        elif 'pairs' in itp_lines[n]: imp_data['pairs'] =n
        elif ('dihedrals' in itp_lines[n]) and (imp_data['prop_dih']>0): imp_data['impr_dih'] =n
    topology = {
                'RESID': itp_lines[imp_data['moleculetype']+1].split()[0],
                'ATOMS': itp_lines[imp_data['atoms']+1:imp_data['bonds']],
                'BONDS': itp_lines[imp_data['bonds']+1:imp_data['angles']],            
                'IMPRP': itp_lines[imp_data['impr_dih']+1:],            
    }
    #######################################
#    topology['RESID'] = topology.RESID.split()[0]
    atom_list = [] 
    for atom in topology['ATOMS']:
        atom = atom.split()
        if len(atom) > 0: atom_list.append([atom[4],atom[1],float(atom[-1])])
    atom_list=pd.DataFrame(atom_list)
    atom_list.columns = ['ATOM','TYPE','CHARGE']
    atom_list.index = range(1,len(atom_list.CHARGE)+1)
    topology['ATOMS'] = atom_list
    num2type = atom_list['ATOM'].to_dict()  
#    print(num2type)
    #######################################
    imp_list = [] 
    for line in topology['IMPRP']:
        line = line.split()
        if len(line) > 0: imp_list.append(list(map(int,line[0:4])))
    imp_list=pd.DataFrame(imp_list)
    imp_list.columns = ['BJ','BI','BK','BL']
    imp_list['I_TYPE']=[num2type[i] for i in imp_list.BI]
    imp_list['J_TYPE']=[num2type[i] for i in imp_list.BJ]
    imp_list['K_TYPE']=[num2type[i] for i in imp_list.BK]
    imp_list['L_TYPE']=[num2type[i] for i in imp_list.BL]
    topology['IMPRP'] = imp_list
    #######################################
    bond_list = [] 
    for bond in topology['BONDS']:
        bond = bond.split()
        if len(bond) == 3: 
            bond = list(map(int,bond))
            bond_list.append(bond[0:2])
    bond_list = pd.DataFrame(bond_list)
    bond_list.columns = ['BI','BJ']
    bond_list['I_TYPE']=[num2type[i] for i in bond_list.BI]
    bond_list['J_TYPE']=[num2type[i] for i in bond_list.BJ]
    topology['BONDS'] = bond_list
    #######################################
    return(topology)
def write_NAMD_rtf(topology): 
    df_atom = topology['ATOMS']
    df_bond = topology['BONDS']
    print('AUTO ANGLES DIHE')
    print('RESI %s 0.000'%topology['RESID'])
    for ind,row in df_atom.iterrows(): 
        print('ATOM %5s %5s %10.5f'%(row.ATOM,row.TYPE,row.CHARGE))
    for ind,row in df_bond.iterrows():
        print('BOND %5s %5s'%(row.I_TYPE,row.J_TYPE,))
    for ind,row in topology['IMPRP'].iterrows():
        print('IMPR %5s %5s %5s %5s'%(row.I_TYPE,row.J_TYPE,row.K_TYPE,row.L_TYPE))
    print('PATCH FIRST NONE LAST NONE \nEND' )
    return None
itp_name = 'dppc.itp'
topol = Get_topology(itp_name)
write_NAMD_rtf(topol)