<jupyter_start><jupyter_text>How to detect resonance structures?<jupyter_code>import pandas as pd
from openeye import oechem, oequacpac
def get_labeled_mol(smiles, label='heavy'):
    """
    returns an OEMol with heavy atoms labeled with a specific indice
    """
    mol = oechem.OEMol()
    if not oechem.OESmilesToMol(mol, smiles):
        print("Couldn't parse smiles (%s) returning None" % smiles)
        return None
    
    for idx, a in enumerate(mol.GetAtomIter(oechem.OEIsHeavy())):
        a.SetData('heavy', idx)
    
    return mol


def find_tautomer(mol1, new_taut):
    """
    returns a tautomer of mol1 that has the same cannonical smile as smiles2(new_taut)
    """
    taut_opts = oequacpac.OETautomerOptions()
    taut_opts.SetCarbonHybridization(True)
    taut_opts.SetLevel(5)
    taut_opts.SetMaxSearchTime(60.0)
    taut_opts.SetRankTautomers(True)
    taut_opts.SetSaveStereo(True) 

    mol = oechem.OEMol(mol1)
    temp_mol = oechem.OEMol()
    oechem.OESmilesToMol(temp_mol, new_taut)
    smiles = oechem.OEMolToSmiles(temp_mol)

    if smiles == oechem.OEMolToSmiles(mol):
        print("Initial mol is the same as the tautomer: {} and {}".format(microstate_ID1, microstate_ID2))
        return mol
    
    for taut in oequacpac.OEEnumerateTautomers(mol, taut_opts):
        s = oechem.OEMolToSmiles(taut)
        if s == smiles:
            return oechem.OEMol(taut)
        
    return None


def is_resonance(mol1, mol2, label='heavy'):
    """
    Returns True is the two molecules are resonance structures of each other
    
    This assumes that each heavy atom in mol1 and mol2 have been mapped correct so 
    identical heavy atoms have the same label 'heavy atom': # in the atom data
    """
    H_counts = dict()
    for a in mol1.GetAtomIter(oechem.OEIsHeavy()):
        H_counts[a.GetData(label)] = a.GetTotalHCount()

    for a in mol2.GetAtomIter(oechem.OEIsHeavy()):
        if a.GetTotalHCount() != H_counts[a.GetData(label)]:
            return False
    
    return True


def test_resonance(smiles1, smiles2):
    """
    Returns True if two smiles are resonance structures of the same microstate.
    Returns False if two smiles are determined not to be resonance structures, 
    because their heavy atom hydrogen counts don't match. 
    Returns None if it is undetermined if these two states are resonance 
    structures or not.
    """
    mol1 = get_labeled_mol(smiles1)
    mol2 = find_tautomer(mol1, smiles2)
    if mol2 is not None:
        res = is_resonance(mol1, mol2)
        #print("These are resonance structures: %s" % res)
        return res  
    else:
        #print("Could not find tautomer with smiles (%s)" % smiles2)
        return None
    
def find_resonance_pairs(df):
    # Check for resonance structure pairs in corrected list
    resonance_pairs_list = []

    # Pairwise comparison of microstates
    for i, row in enumerate(df.iterrows()):
        smiles1 = row[1][1]
        microstate_ID1 = row[1][0]
        # print(microstate_ID1, ":", smiles1)

        for j, row in enumerate(df.iterrows()):
            if j<= i:
                continue
            else:
                smiles2 = row[1][1]
                microstate_ID2 = row[1][0]

                #print("Comparing {} and {}...".format(microstate_ID1, microstate_ID2))
                
                #if total charge is the same test resonance
                if(compare_total_charge(smiles1, smiles2)):
                    #print("Total charge is the same. Must compare if these are resonance structures...")
                    if(test_resonance(smiles1, smiles2)): # if resonance structures
                        # print("{} and {} are resonance structures of the same microstate.".format(microstate_ID1, microstate_ID2))
                        resonance_pairs_list.append((microstate_ID1, microstate_ID2))
                else:
                    continue
    
    return resonance_pairs_list
    
    
    
def get_total_charge(smiles):
    mol = oechem.OEMol()
    oechem.OESmilesToMol(mol,smiles)
    
    total_charge = 0
    for a in mol.GetAtomIter():
        total_charge += a.GetFormalCharge()
    return total_charge

def compare_total_charge(smiles1, smiles2):
    total_charge1 = get_total_charge(smiles1)
    total_charge2 = get_total_charge(smiles2)
    if total_charge1 == total_charge2:
        return True
    else:
        return False
# These two are resonance structures of the same microstate:
res1 = 'c1cc2=[O+]C3=C(NCCCC3=c2cc1O)[O-]'
res2 = 'c1cc2c(cc1O)c3c(o2)C(=O)NCCC3'
test_resonance(res1, res2)<jupyter_output><empty_output><jupyter_text>Compare microstates of SM03 pairwise to eliminate resonance structures<jupyter_code>mol_ID = "SM03"
path_to_input_microstates = "corrections_for_v1_3_1/"
input_file_name = path_to_input_microstates + mol_ID +"_correction.csv"

df_microstates = pd.read_csv(input_file_name)
df_microstates
# Check if there is any deprecated microstate

correction = df_microstates["correction"]
deprecated_boolean = correction.isin(["deprecated"])

deprecated_label = False
for b in deprecated_boolean:
    if b == False:
        continue
    if b == True:
        print("Deprecated microstate found.")
        deprecated_label = True

if(deprecated_label):
    df_corrected = df_microstates.loc[df_microstates["correction"] != "deprecated"]
else:
    df_corrected = df_microstates
    
df_corrected
resonance_pairs_list = []

# Pairwise comparison of microstates
for i, row in enumerate(df_corrected.iterrows()):
    smiles1 = row[1][1]
    microstate_ID1 = row[1][0]
    # print(microstate_ID1, ":", smiles1)
    
    for j, row in enumerate(df_corrected.iterrows()):
        if j<= i:
            continue
        else:
            smiles2 = row[1][1]
            microstate_ID2 = row[1][0]
            
            #print("Comparing {} and {}...".format(microstate_ID1, microstate_ID2))
        
            if(test_resonance(smiles1, smiles2)): # if resonance structures
                # print("{} and {} are resonance structures of the same microstate.".format(microstate_ID1, microstate_ID2))
                resonance_pairs_list.append((microstate_ID1, microstate_ID2))

print("\n")
print(mol_ID, "resonance structures:")
for pair in resonance_pairs_list:
    print(pair)<jupyter_output>SM03 resonance structures:<jupyter_text>Compare microstates of all small molecules pairwise to eliminate resonance structuresAll small molecules except SM18 is scanned. SM18 has 74 microstates so I need a quicker method for that.<jupyter_code>path_to_input_microstates = "corrections_for_v1_3_1/"

for m in range(24):
    mol_ID = "SM"+str(m+1).zfill(2)
    
    # Skip SM18
    if mol_ID == "SM18":
        continue
        
    print(mol_ID, "...")
    
    input_file_name = path_to_input_microstates + mol_ID +"_correction.csv"
    df_microstates = pd.read_csv(input_file_name)
    
    
    # Check if there is any deprecated microstate and create corrected microstates list accordingly.

    correction = df_microstates["correction"]
    deprecated_boolean = correction.isin(["deprecated"])

    deprecated_label = False
    for b in deprecated_boolean:
        if b == False:
            continue
        if b == True:
            print("Deprecated microstate found.")
            deprecated_label = True

    if(deprecated_label):
        df_corrected = df_microstates.loc[df_microstates["correction"] != "deprecated"]
    else:
        df_corrected = df_microstates
        
    # Check for resonance structure pairs in corrected list
    
    resonance_pairs_list = []

    # Pairwise comparison of microstates
    for i, row in enumerate(df_corrected.iterrows()):
        smiles1 = row[1][1]
        microstate_ID1 = row[1][0]
        # print(microstate_ID1, ":", smiles1)

        for j, row in enumerate(df_corrected.iterrows()):
            if j<= i:
                continue
            else:
                smiles2 = row[1][1]
                microstate_ID2 = row[1][0]

                #print("Comparing {} and {}...".format(microstate_ID1, microstate_ID2))
                
                #if total charge is the same test resonance
                if(compare_total_charge(smiles1, smiles2)):
                    #print("Total charge is the same. Must compare if these are resonance structures...")
                    if(test_resonance(smiles1, smiles2)): # if resonance structures
                        # print("{} and {} are resonance structures of the same microstate.".format(microstate_ID1, microstate_ID2))
                        resonance_pairs_list.append((microstate_ID1, microstate_ID2))
                else:
                    continue

    print(mol_ID, "resonance structures:")
    for pair in resonance_pairs_list:
        print(pair)
    print("\n")<jupyter_output>SM01 ...
Deprecated microstate found.
Deprecated microstate found.
SM01 resonance structures:


SM02 ...
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
SM02 resonance structures:


SM03 ...
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
SM03 resonance structures:


SM04 ...
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
SM04 resonance structures:


SM05 ...
Deprecated microstate found.
Deprecated microstate found.
SM05 resonance structures:


SM06 ...
SM06 resonance structures:


SM07 ...
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
Deprecated microstate found.
SM07 resonance structures:


SM08 ...
[...]<jupyter_text>Detecting resonance structures of SM18Only structures with the same total charge could possibly be resonance structures.I will split SM18 microstates accordint to total charge sets and than do pairwise comparison.<jupyter_code>mol_ID = "SM18"
path_to_input_microstates = "corrections_for_v1_3_1/"
input_file_name = path_to_input_microstates + mol_ID +"_correction.csv"

df_microstates = pd.read_csv(input_file_name)

# Check if there is any deprecated microstate

correction = df_microstates["correction"]
deprecated_boolean = correction.isin(["deprecated"])

deprecated_label = False
for b in deprecated_boolean:
    if b == False:
        continue
    if b == True:
        print("Deprecated microstate found.")
        deprecated_label = True

if(deprecated_label):
    df_corrected = df_microstates.loc[df_microstates["correction"] != "deprecated"]
    df_corrected = df_corrected.reset_index(drop=True)
else:
    df_corrected = df_microstates
    
df_corrected["total charge"] = None
for i,row in enumerate(df_corrected.iterrows()):
    smiles = row[1][1]
    microstate_ID = row[1][0]
    df_corrected.loc[i,"total charge"] =get_total_charge(smiles)
    
df_corrected
max_charge = df_corrected["total charge"].max()
min_charge = df_corrected["total charge"].min()
print("Max total charge: ", max_charge)
print("Min total charge: ", min_charge)
df_SM18_minus2 = df_corrected.loc[df_corrected["total charge"] == -2]
print("Number of microstates: ", df_SM18_minus2.shape[0])

caught_resonance_pairs = find_resonance_pairs(df_SM18_minus2)
for pair in caught_resonance_pairs:
    print(pair)
    print("\n")
df_SM18_minus1 = df_corrected.loc[df_corrected["total charge"] == -1]
print("Number of microstates: ", df_SM18_minus1.shape[0])

caught_resonance_pairs = find_resonance_pairs(df_SM18_minus1)
for pair in caught_resonance_pairs:
    print(pair)
    print("\n")
df_SM18_0 = df_corrected.loc[df_corrected["total charge"] == 0]
print("Number of microstates: ", df_SM18_0.shape[0])

caught_resonance_pairs = find_resonance_pairs(df_SM18_0)
for pair in caught_resonance_pairs:
    print(pair)
df_SM18_plus1 = df_corrected.loc[df_corrected["total charge"] == 1]
print("Number of microstates: ", df_SM18_plus1.shape[0])

caught_resonance_pairs = find_resonance_pairs(df_SM18_plus1)
for pair in caught_resonance_pairs:
    print(pair)
df_SM18_plus2 = df_corrected.loc[df_corrected["total charge"] == 2]
print("Number of microstates: ", df_SM18_plus2.shape[0])

caught_resonance_pairs = find_resonance_pairs(df_SM18_plus2)
for pair in caught_resonance_pairs:
    print(pair)<jupyter_output>Number of microstates:  3
{'fraction_non_alphanumeric': 0.09229779707348448, 'fraction_numerical': 0.019858498150828106, 'mean_word_length': 3.2775103163686383, 'pattern_counts': {'":': 4, '<': 17, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '299154', 'license': 'MIT', 'n_tokens_mistral': 4190, 'n_tokens_neox': 3842, 'n_words': 977, 'path': 'pKa_microstate_enumeration/corrected_microstates_v1_3_2/detecting_resonance_structures.ipynb', 'repo_name': 'choderalab/sampl6-logD-compound-selection', 'size': '12438'}