#!/usr/bin/env python

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Descriptors as Des
from rdkit.Chem import Draw
from rdkit.Chem import RWMol
from random import randint
import topomod
import filtermod as fm

def COsingle_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if ((bond.GetBondType() is Chem.BondType.SINGLE and
             bond.GetBeginAtom().GetAtomicNum() == 6 and
             bond.GetEndAtom().GetAtomicNum() == 8) or 
             (bond.GetBondType() is Chem.BondType.SINGLE and
             bond.GetBeginAtom().GetAtomicNum() == 8 and
             bond.GetEndAtom().GetAtomicNum() == 6)):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break
            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions

def OOsingle_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if (bond.GetBondType() is Chem.BondType.SINGLE and
            bond.GetBeginAtom().GetAtomicNum() == 8 and
            bond.GetEndAtom().GetAtomicNum() == 8):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions


def CCsingle_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if (bond.GetBondType() is Chem.BondType.SINGLE and
            bond.GetBeginAtom().GetAtomicNum() == 6 and
            bond.GetEndAtom().GetAtomicNum() == 6):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions

def HHsingle_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if (bond.GetBondType() is Chem.BondType.SINGLE and
            bond.GetBeginAtom().GetAtomicNum() == 1 and
            bond.GetEndAtom().GetAtomicNum() == 1):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions


def COdouble_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if ((bond.GetBondType() is Chem.BondType.DOUBLE and
             bond.GetBeginAtom().GetAtomicNum() == 6 and
             bond.GetEndAtom().GetAtomicNum() == 8) or
            (bond.GetBondType() is Chem.BondType.DOUBLE and
             bond.GetBeginAtom().GetAtomicNum() == 8 and
             bond.GetEndAtom().GetAtomicNum() == 6)):
            
            if(bond.GetBeginAtom().GetIsAromatic() or bond.GetEndAtom().GetIsAromatic()):
                continue

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(temp))
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions


def OOdouble_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if (bond.GetBondType() is Chem.BondType.DOUBLE and
            bond.GetBeginAtom().GetAtomicNum() == 8 and
            bond.GetEndAtom().GetAtomicNum() == 8):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions


def CCdouble_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if (bond.GetBondType() is Chem.BondType.DOUBLE and
            bond.GetBeginAtom().GetAtomicNum() == 6 and
            bond.GetEndAtom().GetAtomicNum() == 6):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions

def CCtriple_dis(ed):
    reactions = []
    for bond in ed.GetBonds():
        at1 = bond.GetBeginAtom().GetIdx()
        at2 = bond.GetEndAtom().GetIdx()
        if (bond.GetBondType() is Chem.BondType.TRIPLE and
            bond.GetBeginAtom().GetAtomicNum() == 6 and
            bond.GetEndAtom().GetAtomicNum() == 6):

            temp = RWMol(ed)
            temp.RemoveBond(at1,at2)
            err = Chem.SanitizeMol(temp,catchErrors=True)
            if err != 0:
                return(['','',''])
            #print(Chem.MolToSmiles(temp))
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at1)
            topomod.remove_H(temp,at2)
            topomod.remove_H(temp,at2)
            topomod.remove_H(temp,at2)
            # check for fragmentation
            frags = Chem.GetMolFrags(temp,asMols=True, sanitizeFrags=False)
            if len(frags) == 1:
                pro1 = frags[0]
                pro2 = Chem.MolFromSmiles('')
            elif len(frags) == 2:
                pro1 = frags[0]
                pro2 = frags[1]
            # check for neighboring radicals to modify bondorder
            for mol in [pro1,pro2]:
                for atom in mol.GetAtoms():
                    if atom.GetNumRadicalElectrons() > 0:
                        for neigh in atom.GetNeighbors():
                            if neigh.GetNumRadicalElectrons() > 0:
                                #mol = RWMol(mol)
                                #print("!!   !!",Chem.MolToSmiles(mol))
                                at1 = atom.GetIdx()
                                at2 = neigh.GetIdx()
                                if (mol.GetBondBetweenAtoms(at1,at2).GetBondType() != Chem.BondType.TRIPLE):
                                    topomod.add_H(mol,at1)
                                    topomod.add_H(mol,at2)
                                    topomod.increment_bond(mol,mol.GetBondBetweenAtoms(at1,at2).GetIdx())
                                    err = Chem.SanitizeMol(mol,catchErrors=True)
                                    if err != 0:
                                        return(['','',''])
                                    #print("!!   !!",Chem.MolToSmiles(mol))
                                    break

            duplicate = False
            for prere in reactions:
                if (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro2)):
                    duplicate = True
                elif (Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro2) and
                    Chem.MolToSmiles(prere[2]) == Chem.MolToSmiles(pro1)):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1),Chem.MolToSmiles(pro2))
    return reactions


def CHdis_modBO(mol):
    reactions = []
    ed = RWMol(mol)
    Chem.Kekulize(ed)
    for atom in ed.GetAtoms():
        if atom.GetAtomicNum() == 6:
            if atom.GetTotalNumHs() > 0:
                at1 = atom.GetIdx()
                radneighbor = False
                for neighbor in atom.GetNeighbors():
                    if neighbor.GetNumRadicalElectrons() > 0:
                        radneighbor = True
                        at2 = neighbor.GetIdx()
                if not radneighbor:
                    # This reaction is hardcoded, because the CO molecule is weird in SMILES
                    # update: this is obsolete, since we now use [C]=O for CO
                    if Chem.MolToSmiles(ed) == '[CH]=O':
                        pro1 = Chem.MolFromSmiles('[C-]#[O+]')
                        pro2 = Chem.MolFromSmiles('[H]')
                        # print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
                    else:
                        continue
                elif ed.GetBondBetweenAtoms(at1,at2).GetBondType() is Chem.BondType.TRIPLE:
                    continue 
                else:
                    pro1 = RWMol(ed)
                    #topomod.remove_H(pro1,at1)
                    topomod.add_H(pro1,at2)
                    topomod.increment_bond(pro1,ed.GetBondBetweenAtoms(at1,at2).GetIdx())
                    pro2 = Chem.MolFromSmiles('[H]')
                    err = Chem.SanitizeMol(ed,catchErrors=True)
                    if err != 0:
                        return(['','',''])
                    err = Chem.SanitizeMol(pro1,catchErrors=True)
                    if err != 0:
                        return(['','',''])
                duplicate = False
                for prere in reactions:
                    if Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1):
                        duplicate = True
                if not duplicate:
                    reactions.append([ed,pro1,pro2])
                    #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions

def CHdis_ring(mol):
    reactions = []
    ed = RWMol(mol)
    if Chem.GetSSSR(ed) > 0:
        return reactions
    for atom in ed.GetAtoms():
        rad = False
        at1 = atom.GetIdx()
        if atom.GetAtomicNum() == 6:
            if atom.GetTotalNumHs() > 0:
                for atom2 in ed.GetAtoms():
                    at2 = atom2.GetIdx()
                    if ((atom2.GetNumRadicalElectrons() > 0) and
                        (at2 != at1)): 
                        notNeigh = True
                        for neighbor in atom.GetNeighbors():
                            if neighbor.GetIdx() == at2:
                                notNeigh = False                            
                        if notNeigh:
                            rad = True
                            break
        if not rad:
            continue
        else:
            pro1 = RWMol(ed)
            #topomod.remove_H(pro1,at1)
            #topomod.add_H(pro1,at2)
            #topomod.increment_bond(pro1,ed.GetBondBetweenAtoms(at1,at2).GetIdx())
            #print(at1,at2)
            topomod.remove_radical(pro1,at2)
            topomod.add_bond(pro1,at1,at2)
            pro2 = Chem.MolFromSmiles('[H]')
            err = Chem.SanitizeMol(ed,catchErrors=True)
            if err != 0:
                return(['','',''])
            err = Chem.SanitizeMol(pro1,catchErrors=True)
            if err != 0:
                return(['','',''])
            duplicate = False
            for prere in reactions:
                if Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions


def OHdis_modBO(mol):
    reactions = []
    ed = RWMol(mol)
    Chem.Kekulize(ed)
    for atom in ed.GetAtoms():
        if atom.GetAtomicNum() == 8:
            if atom.GetTotalNumHs() > 0:
                at1 = atom.GetIdx()
                radneighbor = False
                for neighbor in atom.GetNeighbors():
                    if neighbor.GetNumRadicalElectrons() > 0:
                        radneighbor = True
                        at2 = neighbor.GetIdx()
                if not radneighbor:
                    continue
                elif ed.GetBondBetweenAtoms(at1,at2).GetBondType() is Chem.BondType.TRIPLE:
                    continue
                else:
                    pro1 = RWMol(ed)
                    #topomod.remove_H(pro1,at1)
                    topomod.add_H(pro1,at2)
                    topomod.increment_bond(pro1,ed.GetBondBetweenAtoms(at1,at2).GetIdx())
                    pro2 = Chem.MolFromSmiles('[H]')
                    err = Chem.SanitizeMol(ed,catchErrors=True)
                    if err != 0:
                        return(['','',''])
                    err = Chem.SanitizeMol(pro1,catchErrors=True)
                    if err != 0:
                        return(['','',''])
                duplicate = False
                for prere in reactions:
                    if Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1):
                        duplicate = True
                if not duplicate:
                    reactions.append([ed,pro1,pro2])
                    #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions

def OHdis_ring(mol):
    reactions = []
    ed = RWMol(mol)
    if Chem.GetSSSR(ed) > 0:
        return reactions
    for atom in ed.GetAtoms():
        rad = False
        at1 = atom.GetIdx()
        if atom.GetAtomicNum() == 8:
            if atom.GetTotalNumHs() > 0:
                for atom2 in ed.GetAtoms():
                    at2 = atom2.GetIdx()
                    if ((atom2.GetNumRadicalElectrons() > 0) and
                        (at2 != at1)):
                        notNeigh = True
                        for neighbor in atom.GetNeighbors():
                            if neighbor.GetIdx() == at2:
                                notNeigh = False
                        if notNeigh:
                            rad = True
                            break
        if not rad:
            continue
        else:
            pro1 = RWMol(ed)
            #topomod.remove_H(pro1,at1)
            #topomod.add_H(pro1,at2)
            #topomod.increment_bond(pro1,ed.GetBondBetweenAtoms(at1,at2).GetIdx())
            #print(at1,at2)
            topomod.remove_radical(pro1,at2)
            topomod.add_bond(pro1,at1,at2)
            pro2 = Chem.MolFromSmiles('[H]')
            err = Chem.SanitizeMol(ed,catchErrors=True)
            if err != 0:
                return(['','',''])
            err = Chem.SanitizeMol(pro1,catchErrors=True)
            if err != 0:
                return(['','',''])
            duplicate = False
            for prere in reactions:
                if Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1):
                    duplicate = True
            if not duplicate:
                reactions.append([ed,pro1,pro2])
                #print(Chem.MolToSmiles(ed),Chem.MolToSmiles(pro1))
    return reactions


def CHdis(ed):
    reactions = []
    for atom in ed.GetAtoms():
        if atom.GetAtomicNum() == 6:
            if atom.GetTotalNumHs() > 0:
                at = atom.GetIdx()
                pro1 = RWMol(ed)
                topomod.remove_H(pro1,at)
                pro2 = Chem.MolFromSmiles('[H]')
                duplicate = False
                for prere in reactions:
                    if Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1):
                        duplicate = True
                if not duplicate:
                    reactions.append([ed,pro1,pro2])
    return reactions 


def OHdis(ed):
    reactions = []
    for atom in ed.GetAtoms():
        if atom.GetAtomicNum() == 8:
            if atom.GetTotalNumHs() > 0:
                at = atom.GetIdx()
                pro1 = RWMol(ed)
                topomod.remove_H(pro1,at)
                pro2 = Chem.MolFromSmiles('[H]')
                duplicate = False
                for prere in reactions:
                    if Chem.MolToSmiles(prere[1]) == Chem.MolToSmiles(pro1):
                        duplicate = True
                if not duplicate:
                    reactions.append([ed,pro1,pro2])
    return reactions


def checkcomp(ed,pro):
    AtomTypes = [1,6,8]
    EdTypes = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    ProTypes = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]    

    for educt in ed:
        for atom in educt.GetAtoms():
           EdTypes[1] += atom.GetTotalNumHs()
           for atype in AtomTypes:
               if atom.GetAtomicNum() == atype:
                   EdTypes[atype] += 1

#   print(EdTypes)
    for product in pro:
        for atom in product.GetAtoms():
           ProTypes[1] += atom.GetTotalNumHs()
           for atype in AtomTypes:
               if atom.GetAtomicNum() == atype:
                   ProTypes[atype] += 1

#   print(ProTypes)

    check = True
    for i in range(len(EdTypes)):
        if EdTypes[i] != ProTypes[i]:
            check = False
#   print(check)
    return check




# check changes in the backbone
def classify_reactions(edlist,prolist,Verbose):
    EdBonds = []
    Unmatched = []
    nbonds = 0
    for ed in edlist:
        nbonds += ed.GetNumBonds()
    unmatch = 0
    for ed in edlist:
        for bond in ed.GetBonds():
            EdBonds.append([bond.GetBeginAtom().GetAtomicNum(),
                            bond.GetEndAtom().GetAtomicNum(),
                            bond.GetBondType()])

    for pro in prolist:
        for bond in pro.GetBonds():
            Found = False
            for i in range(len(EdBonds)):
                if((bond.GetBeginAtom().GetAtomicNum() == EdBonds[i][0]) and
                   (bond.GetEndAtom().GetAtomicNum() == EdBonds[i][1]) and
                   (bond.GetBondType() == EdBonds[i][2])):
                    del EdBonds[i]
                    Found = True
                    break
                elif ((bond.GetBeginAtom().GetAtomicNum() == EdBonds[i][1]) and
                   (bond.GetEndAtom().GetAtomicNum() == EdBonds[i][0]) and
                   (bond.GetBondType() == EdBonds[i][2])):
                    del EdBonds[i]
                    Found = True
                    break
            if not Found:
                unmatch += 1
                Unmatched.append([bond.GetBeginAtom().GetAtomicNum(),
                                  bond.GetEndAtom().GetAtomicNum(),
                                  bond.GetBondType()])

    UnmatchedEductBonds = len(EdBonds)
    UnmatchedProductBonds  = unmatch
#   print("    Unmatched bonds in educt   : %i" % len(EdBonds))
#   print("    Unmatched bonds in products: %i" % unmatch)

# Looking for H-H bonds
    prodHHsing = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 1 and
           Unmatched[i][1] == 1 and
           Unmatched[i][2] == Chem.BondType.SINGLE):
            prodHHsing += 1
    edHHsing = 0
    for i in range(len(EdBonds)):
        if(EdBonds[i][0] == 1 and
           EdBonds[i][1] == 1 and
           EdBonds[i][2] == Chem.BondType.SINGLE):
            edHHsing += 1
#   print("    Unmatched H-H in educt     : %i" % edHHsing)
#   print("    Unmatched H-H in product   : %i" % prodHHsing)


# Looking for C-C bonds
    prodCCsing = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 6 and 
           Unmatched[i][1] == 6 and 
           Unmatched[i][2] == Chem.BondType.SINGLE):
            prodCCsing += 1
    edCCsing = 0
    for i in range(len(EdBonds)):
        if(EdBonds[i][0] == 6 and
           EdBonds[i][1] == 6 and
           EdBonds[i][2] == Chem.BondType.SINGLE):
            edCCsing += 1
#   print("    Unmatched C-C in educt     : %i" % edCCsing)
#   print("    Unmatched C-C in product   : %i" % prodCCsing)

# Looking for C=C bonds
    prodCCdoub = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 6 and 
           Unmatched[i][1] == 6 and 
           Unmatched[i][2] == Chem.BondType.DOUBLE):
            prodCCdoub += 1
    edCCdoub = 0
    for i in range(len(EdBonds)):
        if(EdBonds[i][0] == 6 and
           EdBonds[i][1] == 6 and
           EdBonds[i][2] == Chem.BondType.DOUBLE):
            edCCdoub += 1
#   print("    Unmatched C=C in educt     : %i" % edCCdoub)
#   print("    Unmatched C=C in product   : %i" % prodCCdoub)

# Looking for C#C bonds
    prodCCtrip = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 6 and
           Unmatched[i][1] == 6 and 
           Unmatched[i][2] == Chem.BondType.TRIPLE):
            prodCCtrip += 1
    edCCtrip = 0
    for i in range(len(EdBonds)):
        if(EdBonds[i][0] == 6 and
           EdBonds[i][1] == 6 and 
           EdBonds[i][2] == Chem.BondType.TRIPLE):
            edCCtrip += 1
#   print("    Unmatched C#C in educt     : %i" % edCCtrip)
#   print("    Unmatched C#C in product   : %i" % prodCCtrip)

# Looking for O-O bonds
    prodOOsing = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 8 and 
           Unmatched[i][1] == 8 and 
           Unmatched[i][2] == Chem.BondType.SINGLE):
            prodOOsing += 1
    edOOsing = 0
    for i in range(len(EdBonds)):
        if(EdBonds[i][0] == 8 and
           EdBonds[i][1] == 8 and
           EdBonds[i][2] == Chem.BondType.SINGLE):
            edOOsing += 1
#   print("    Unmatched O-O in educt     : %i" % edOOsing)
#   print("    Unmatched O-O in product   : %i" % prodOOsing)

# Looking for O=O bonds
    prodOOdoub = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 8 and
           Unmatched[i][1] == 8 and 
           Unmatched[i][2] == Chem.BondType.DOUBLE):
            prodOOdoub += 1
    edOOdoub = 0
    for i in range(len(EdBonds)):
        if(EdBonds[i][0] == 8 and
           EdBonds[i][1] == 8 and 
           EdBonds[i][2] == Chem.BondType.DOUBLE):
            edOOdoub += 1
#   print("    Unmatched O=O in educt     : %i" % edOOdoub)
#   print("    Unmatched O=O in product   : %i" % prodOOdoub)


# Looking for C-O bonds
    prodCOsing = 0
    for i in range(len(Unmatched)):
        if(((Unmatched[i][0] == 6 and
             Unmatched[i][1] == 8) or
            (Unmatched[i][0] == 8 and
             Unmatched[i][1] == 6)) and
             Unmatched[i][2] == Chem.BondType.SINGLE):
            prodCOsing += 1
    edCOsing = 0
    for i in range(len(EdBonds)):
        if(((EdBonds[i][0] == 6 and
             EdBonds[i][1] == 8) or
            (EdBonds[i][0] == 8 and
             EdBonds[i][1] == 6)) and
             EdBonds[i][2] == Chem.BondType.SINGLE):
            edCOsing += 1
#   print("    Unmatched C-O in educt     : %i" % edCOsing)
#   print("    Unmatched C-O in product   : %i" % prodCOsing)


# Looking for C=O bonds
    prodCOdoub = 0
    for i in range(len(Unmatched)):
        if(((Unmatched[i][0] == 6 and
             Unmatched[i][1] == 8) or
            (Unmatched[i][0] == 8 and
             Unmatched[i][1] == 6)) and
             Unmatched[i][2] == Chem.BondType.DOUBLE):
            prodCOdoub += 1
    edCOdoub = 0
    for i in range(len(EdBonds)):
        if(((EdBonds[i][0] == 6 and
             EdBonds[i][1] == 8) or
            (EdBonds[i][0] == 8 and
             EdBonds[i][1] == 6)) and
             EdBonds[i][2] == Chem.BondType.DOUBLE):
            edCOdoub += 1
#   print("    Unmatched C=O in educt     : %i" % edCOdoub)
#   print("    Unmatched C=O in product   : %i" % prodCOdoub)


# Looking for C#O bonds
    prodCOtrip = 0
    for i in range(len(Unmatched)):
        if(((Unmatched[i][0] == 6 and
             Unmatched[i][1] == 8) or
            (Unmatched[i][0] == 8 and
             Unmatched[i][1] == 6)) and
             Unmatched[i][2] == Chem.BondType.TRIPLE):
            prodCOtrip += 1
    edCOtrip = 0
    for i in range(len(EdBonds)):
        if(((EdBonds[i][0] == 6 and
             EdBonds[i][1] == 8) or
            (EdBonds[i][0] == 8 and
             EdBonds[i][1] == 6)) and
             EdBonds[i][2] == Chem.BondType.TRIPLE):
            edCOtrip += 1
#   print("    Unmatched C#O in educt     : %i" % edCOtrip)
#   print("    Unmatched C#O in product   : %i" % prodCOtrip)


# check changes in hydrogenation
    EdAtoms = []
    Unmatched = []
    nAtoms = 0
    for ed in edlist:
        nAtoms += ed.GetNumHeavyAtoms()
    unmatch = 0
    for ed in edlist:
        for atom in ed.GetAtoms():
            if atom.GetAtomicNum() == 1:
                continue
            EdAtoms.append([atom.GetAtomicNum(),
                            atom.GetTotalNumHs()])
#   print(EdAtoms)

    for pro in prolist:
        for atom in pro.GetAtoms():
            Found = False
            if atom.GetAtomicNum() == 1:
                continue
            for i in range(len(EdAtoms)):
                if ((atom.GetAtomicNum()  == EdAtoms[i][0]) and
                    (atom.GetTotalNumHs() == EdAtoms[i][1])):
                    del EdAtoms[i]
                    Found = True
                    break
            if not Found:
                unmatch += 1
                Unmatched.append([atom.GetAtomicNum(),atom.GetTotalNumHs()])

#   print(EdAtoms)
#   print(Unmatched)

#   print("    Unmatched hydro in products: %i" % unmatch)
#   print("    Unmatched hydro in educt   : %i" % len(EdAtoms))
    diffCH = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 6):
            diffCH += Unmatched[i][1]
    for i in range(len(EdAtoms)):
        if(EdAtoms[i][0] == 6):
            diffCH -= EdAtoms[i][1]
#   print("    Difference in CH bonds     : %i" % diffCH)

    diffOH = 0
    for i in range(len(Unmatched)):
        if(Unmatched[i][0] == 8):
            diffOH += Unmatched[i][1]
    for i in range(len(EdAtoms)):
        if(EdAtoms[i][0] == 8):
            diffOH -= EdAtoms[i][1]
#   print("    Difference in OH bonds     : %i" % diffOH)


    if Verbose:
        print("    -------------------------------")
        print("    Unmatched bonds in educt   : %i" % len(EdBonds))
        print("    Unmatched bonds in products: %i" % unmatch)
        print("    Unmatched H-H in educt     : %i" % edHHsing)
        print("    Unmatched H-H in product   : %i" % prodHHsing)
        print("    Unmatched C-C in educt     : %i" % edCCsing)
        print("    Unmatched C-C in product   : %i" % prodCCsing)
        print("    Unmatched C=C in educt     : %i" % edCCdoub)
        print("    Unmatched C=C in product   : %i" % prodCCdoub)
        print("    Unmatched C#C in educt     : %i" % edCCtrip)
        print("    Unmatched C#C in product   : %i" % prodCCtrip)
        print("    Unmatched O-O in educt     : %i" % edOOsing)
        print("    Unmatched O-O in product   : %i" % prodOOsing)
        print("    Unmatched O=O in educt     : %i" % edOOdoub)
        print("    Unmatched O=O in product   : %i" % prodOOdoub)
        print("    Unmatched C-O in educt     : %i" % edCOsing)
        print("    Unmatched C-O in product   : %i" % prodCOsing)
        print("    Unmatched C=O in educt     : %i" % edCOdoub)
        print("    Unmatched C=O in product   : %i" % prodCOdoub)
        print("    Unmatched C#O in educt     : %i" % edCOtrip)
        print("    Unmatched C#O in product   : %i" % prodCOtrip)    
        print("    Unmatched hydro in products: %i" % unmatch)
        print("    Unmatched hydro in educt   : %i" % len(EdAtoms))
        print("    Difference in CH bonds     : %i" % diffCH)
        print("    Difference in OH bonds     : %i" % diffOH)

    # Determine Reaction Type
    REtype = "Non Elementary"
    if(len(edlist)==1 and len(prolist)==2):
    # Type A -> B + C
        if(UnmatchedProductBonds == 0 and UnmatchedEductBonds == 0):
        # H Abstractions
            if (diffCH == -1 and diffOH == 0):
                REtype = "CH-Abstraction"
            elif(diffCH == 0 and diffOH == -1):
                REtype = "OH-Abstraction"
        elif(UnmatchedProductBonds == 0 and UnmatchedEductBonds == 1):
        # Bond Cleavage
            if (edCCsing==1):
                REtype = "C-C-Cleavage"
            elif (edOOsing==1):
                REtype = "O-O-Cleavage"
            elif (edHHsing==1):
                REtype = "H-H-Cleavage"
            elif (edCOsing==1):
                REtype = "C-O-Cleavage"
            elif (edCCdoub==1):
                REtype = "C=C-Cleavage"
            elif (edOOdoub==1):
                REtype = "O=O-Cleavage"
            elif (edCOdoub==1):
                REtype = "C=O-Cleavage"
            elif (edCCtrip==1):
                REtype = "C#C-Cleavage"
            elif (edCOtrip==1):
                REtype = "C#O-Cleavage"
        elif(UnmatchedProductBonds == 1 and UnmatchedEductBonds == 1):
        # Radical at double bond to triple bond via H-abstraction
            if (edCCdoub==1 and prodCCtrip==1 and diffCH == -1 and diffOH == 0):
                REtype = "CH-Abstraction (= > #)"
            if (edCOdoub==1 and prodCOtrip==1 and diffCH == -1 and diffOH == 0):
                REtype = "CH-Abstraction (= > #)"
        # Radical at single bond to double bond via H-abstraction
            if (edCCsing==1 and prodCCdoub==1 and diffCH == -1 and diffOH == 0):
                REtype = "CH-Abstraction (- > =)"
            if (edCOsing==1 and prodCOdoub==1 and diffCH == -1 and diffOH == 0):
                REtype = "CH-Abstraction (- > =)"
            if (edCOsing==1 and prodCOdoub==1 and diffCH == 0 and diffOH == -1):
                REtype = "OH-Abstraction (- > =)"
        elif(UnmatchedProductBonds == 2 and UnmatchedEductBonds == 1 and prodHHsing == 1):
        # H2 abstraction with bond transformation 
            REtype = "H2-Abstraction"
        elif(UnmatchedProductBonds == 1 and UnmatchedEductBonds == 1 and prodHHsing == 1):
        # H2 abstraction without bond transformation
            REtype = "H2-Abstraction"
        elif(UnmatchedProductBonds == 1 and UnmatchedEductBonds == 0 and prodHHsing == 1):
        # H2 abstraction from monoatomic backbone (e.g. H2O)
            REtype = "H2-Abstraction"
    if(len(edlist)==2 and len(prolist)==2):
    # Type A + B -> C + D
        if(UnmatchedProductBonds == 0 and UnmatchedEductBonds == 0):
            pass
        elif((UnmatchedProductBonds == 0 and UnmatchedEductBonds == 1) or
             (UnmatchedProductBonds == 1 and UnmatchedEductBonds == 0)):
        # bond break due to hydrogenation
            pass
        elif(UnmatchedProductBonds == 1 and UnmatchedEductBonds == 1):
            pass

#   print("    Reaction Type:   %s" % REtype)
    return REtype


