#!/usr/bin/python
from openeye.oechem import *
import os, string, sys

try:
    ifs = oemolistream()
    ifs.open(sys.argv[1])
    ofs_1 = oemolostream()
    ofs_2 = oemolostream()
    ofs_1.open('b1'+sys.argv[1][string.find(sys.argv[1], '_'):string.find(sys.argv[1], '_out')]+'_ln_out.ism')
    ofs_2.open('b1'+sys.argv[1][string.find(sys.argv[1], '_'):string.find(sys.argv[1], '_out')]+'_lp_out.ism')
    date = sys.argv[1][string.find(sys.argv[1], '_'):string.find(sys.argv[1], '_out')]
except IndexError:
    print '\nUsage: b1_rxn_carbonyl.py molfile\n'
    raise SystemExit()

# Reaction 1: OH- nucleophilic attack at A-C(=OvS)-X (amide, ester etc., OvS exocyclic) (NON-AROMATIC)
# Former Carbonyl-oxygen becomes OH(SH)-group, attacking OH- becomes O-

# leaving group NOT protonated  [*:4] included, because leaving group needs to be recognized only once (if there are 2!)
mol_hyd = OEGraphMol()
OEParseSmiles(mol_hyd, "[OH]")
libgen1 = OELibraryGen("[#6,#15,#8&!H1,#7,#16&!H1;!-:4][C:1](=[O,S:3])[#8&!H1,#7,#16&!H1;!-:2].[O:100][H:101]>>[*:4][C:1]([*:3][H:101])([O-:100])[*:2]")
libgen1.SetExplicitHydrogens(True)
libgen1.SetValenceCorrection(False)

for mol in ifs.GetOEMols():
    libgen1.AddStartingMaterial(mol, 0)
    libgen1.AddStartingMaterial(mol_hyd, 1)
    mol.Clear()
    mol_hyd.Clear()

title_f = open('b1_titles.txt', 'w')
title_list = []
for products in libgen1.GetProducts():
    num = 1
    prot = 0
    title = products.GetTitle()+str(num)+'p'+str(prot)
    while title in title_list:
        num+=1
        title = products.GetTitle()+str(num)+'p'+str(prot)
    title_list.append(title)
    products.SetTitle(title)
    OEWriteMolecule(ofs_1, products)
    title_f.write(title+'\n')

ofs_1.close()
ifs.close()
ifs.open('b1'+sys.argv[1][string.find(sys.argv[1], '_'):string.find(sys.argv[1], '_out')]+'_ln_out.ism')

mol_prot = OEGraphMol()
OEParseSmiles(mol_prot, "[H+]")
# leaving group protonated 
libgen2 = OELibraryGen("[C:1]([O:3])([O,S;-:4])[#8&!H1,#7,#16&!H1;!-:2].[H+:101]>>[C:1]([*:3])([A:4])[*+:2][H+0:101]")
libgen2.SetExplicitHydrogens(True)
libgen2.SetValenceCorrection(False)

for mol in ifs.GetOEMols():
    libgen2.AddStartingMaterial(mol, 0)
    libgen2.AddStartingMaterial(mol_prot, 1)
    mol.Clear()
    mol_prot.Clear()


for products in libgen2.GetProducts():
    prot = 0
    title = products.GetTitle()[:string.find(products.GetTitle(), 'p0')]+'p'+str(prot)
    while title in title_list:
        prot+=1
        title = products.GetTitle()[:string.find(products.GetTitle(), 'p')]+'p'+str(prot)
    title_list.append(title)
    products.SetTitle(title)
    OEWriteMolecule(ofs_2, products)
    title_f.write(title+'\n')

ofs_2.close()
ifs.close()
title_f.close()