comparing_smiles/comparing_smiles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4d115db",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import swifter\n",
    "from rdkit import Chem, RDLogger\n",
    "from rdkit.Chem import rdchem\n",
    "from rdkit.Chem.Draw import IPythonConsole\n",
    "from rdkit.Chem.MolStandardize import rdMolStandardize\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "RDLogger.DisableLog(\"rdApp.info\")\n",
    "\n",
    "\n",
    "from importlib import reload\n",
    "\n",
    "import compare_smiles_code\n",
    "import pandas as pd\n",
    "\n",
    "reload(compare_smiles_code)\n",
    "\n",
    "from compare_smiles_code import (\n",
    "    check_smiles,\n",
    "    check_smiles_fragments_resonance,\n",
    "    check_smiles_tmc_resonance,\n",
    "    count,\n",
    "    get_agreement_between_multiple_smiles,\n",
    "    get_agreement_between_smiles,\n",
    "    process_paralell,\n",
    ")\n",
    "from hydrogen_parser import *\n",
    "\n",
    "import xyz2mol_tm.NBO_to_smiles.mol_utils as mol_utils\n",
    "from xyz2mol_tm.NBO_to_smiles.mol_utils import strip_stereo\n",
    "\n",
    "IPythonConsole.molSize = 450, 450\n",
    "IPythonConsole.drawOptions.addAtomIndices = False\n",
    "IPythonConsole.drawOptions.useSVG = True\n",
    "IPythonConsole.drawOptions.minFontSize = 24\n",
    "# IPythonConsole.drawOptions.legendFraction = 0.2\n",
    "IPythonConsole.molSize = 500, 500\n",
    "\n",
    "\n",
    "def print_to_file(line, file_name):\n",
    "    with open(f\"{file_name}\", \"a\") as f:\n",
    "        f.write(line + \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d2e29ce",
   "metadata": {},
   "source": [
    "# Comparing the SMILES obtained with the 3 methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cb8f90b",
   "metadata": {},
   "source": [
    "### Load the smiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2407f52",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_frame = pd.read_csv(\"../SMILES_csvs/tmqmg_smiles.csv\")\n",
    "\n",
    "\n",
    "# Strip chiral tags:\n",
    "combined_frame[\"smiles_NBO_DFT_xyz\"] = combined_frame[\"smiles_NBO_DFT_xyz\"].swifter.apply(mol_utils.strip_stereo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72b0e8f6",
   "metadata": {},
   "source": [
    "### Instantiate the masks checking for invalid smiles values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42c64a51",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_csd = combined_frame[\"smiles_CSD_fix\"].apply(check_smiles)\n",
    "mask_tmqmg = combined_frame[\"smiles_NBO_DFT_xyz\"].apply(check_smiles)\n",
    "mask_huckel = combined_frame[\"smiles_huckel_DFT_xyz\"].apply(check_smiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5207185",
   "metadata": {},
   "source": [
    "### How many smiles for each method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7efe8ad",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Aligning the print statements\n",
    "print(f\"CSD smiles:     {len(combined_frame[mask_csd]):>5} / {len(combined_frame):>5}\")\n",
    "print(f\"NBO smiles:     {len(combined_frame[mask_tmqmg]):>5} / {len(combined_frame):>5}\")\n",
    "print(f\"Huckel smiles:  {len(combined_frame[mask_huckel]):>5} / {len(combined_frame):>5}\")\n",
    "print(f\"3 valid SMILES: {len(combined_frame[mask_csd & mask_tmqmg & mask_huckel]):>5} / {len(combined_frame):>5}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87ac9b36",
   "metadata": {},
   "outputs": [],
   "source": [
    "not_valid_smiles = combined_frame[~(mask_csd | mask_tmqmg | mask_huckel)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fade6cdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"{len(not_valid_smiles)} does not have a valid SMILES for either approach\")\n",
    "print(\"We remove them from the dataframe\")\n",
    "combined_frame.drop(not_valid_smiles.index.values, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "896c8602",
   "metadata": {},
   "source": [
    "### SMILES with incorrect formulas (missing hydrogens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4efa3077",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "mask_csd = combined_frame[\"smiles_CSD_fix\"].apply(check_smiles)\n",
    "mask_tmqmg = combined_frame[\"smiles_NBO_DFT_xyz\"].apply(check_smiles)\n",
    "mask_huckel = combined_frame[\"smiles_huckel_DFT_xyz\"].apply(check_smiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb238b3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = combined_frame[(mask_csd | mask_tmqmg | mask_huckel)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8235f1cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_frame.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a08132d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ground truth number of hydrogens from the CSD API formula\n",
    "df[\"number_hydrogen_csd_api\"] = df[\"formula_heaviest_fragment\"].apply(parse_formula, task=\"hydrogen\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d850d75",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Step 2: Apply the function only to the rows where the mask is True\n",
    "df[\"formula_smiles_huckel_DFT_xyz\"] = np.nan  # Initialize with NaN\n",
    "df.loc[mask_huckel, \"formula_smiles_huckel_DFT_xyz\"] = df.loc[mask_huckel, \"smiles_huckel_DFT_xyz\"].swifter.apply(get_formula_from_smiles)\n",
    "\n",
    "df[\"formula_csd_smiles\"] = np.nan  # Initialize with NaN\n",
    "df.loc[mask_csd, \"formula_csd_smiles\"] = df.loc[mask_csd, \"smiles_CSD_fix\"].swifter.apply(get_formula_from_smiles)\n",
    "\n",
    "df[\"formula_nbo_smiles\"] = np.nan  # Initialize with NaN\n",
    "df.loc[mask_tmqmg, \"formula_nbo_smiles\"] = df.loc[mask_tmqmg, \"smiles_NBO_DFT_xyz\"].swifter.apply(get_formula_from_smiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35b69e65",
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"number_hydrogen_huckel_formula\"] = df.swifter.apply(\n",
    "    lambda x: parse_formula(\n",
    "        x.formula_smiles_huckel_DFT_xyz,\n",
    "        task=\"hydrogen\",\n",
    "    ),\n",
    "    axis=1,\n",
    ")\n",
    "df[\"number_hydrogen_csd_formula\"] = df.swifter.apply(\n",
    "    lambda x: parse_formula(\n",
    "        x.formula_csd_smiles,\n",
    "        task=\"hydrogen\",\n",
    "    ),\n",
    "    axis=1,\n",
    ")\n",
    "df[\"number_hydrogen_nbo_formula\"] = df.swifter.apply(\n",
    "    lambda x: parse_formula(\n",
    "        x.formula_nbo_smiles,\n",
    "        task=\"hydrogen\",\n",
    "    ),\n",
    "    axis=1,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c63d997",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "\n",
    "csd_arr = df[mask_csd][\"number_hydrogen_csd_formula\"]\n",
    "ground_truth = df[\"number_hydrogen_csd_api\"]\n",
    "nbo_arr = df[mask_tmqmg][\"number_hydrogen_nbo_formula\"]\n",
    "huckel_arr = df[mask_huckel][\"number_hydrogen_huckel_formula\"]\n",
    "\n",
    "csd_comp = csd_arr != ground_truth[mask_csd]\n",
    "nbo_comp = nbo_arr != ground_truth[mask_tmqmg]\n",
    "huckel_comp = huckel_arr != ground_truth[mask_huckel]\n",
    "\n",
    "print(f\"Number of missing hydrogens for each method:\\n\")\n",
    "print(f\"CSD {sum(csd_comp)}\")\n",
    "print(f\"NBO: {sum(nbo_comp)}\")\n",
    "print(f\"Huckel: {sum(huckel_comp)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da816f61",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_wrong_formula = df[(csd_comp & nbo_comp & huckel_comp)]\n",
    "print(f\"{len(all_wrong_formula)} Entries have wrong or missing formula for all methods. We drop them\")\n",
    "df.drop(all_wrong_formula.index.values, inplace=True)\n",
    "print(len(df))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1ab2008",
   "metadata": {},
   "source": [
    "To consistently compare SMILES, we remove all rows where one of the three methods have a SMILES with incorrect number of hydrogens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c2ecb98",
   "metadata": {},
   "outputs": [],
   "source": [
    "one_wrong_formula = df[(csd_comp | nbo_comp | huckel_comp)]\n",
    "print(f\"{len(one_wrong_formula)} Entries have wrong or missing formula for at least one of the methods. We drop them\")\n",
    "df.drop(one_wrong_formula.index.values, inplace=True)\n",
    "print(\"Remaining entries in the dataset:\", len(df))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baafdbea",
   "metadata": {},
   "source": [
    "## Comparing SMILES directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f30570b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_csd = combined_frame[\"smiles_CSD_fix\"].apply(check_smiles)\n",
    "mask_tmqmg = combined_frame[\"smiles_NBO_DFT_xyz\"].apply(check_smiles)\n",
    "mask_huckel = combined_frame[\"smiles_huckel_DFT_xyz\"].apply(check_smiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c621d045",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_have_smiles = combined_frame[(mask_csd & mask_tmqmg) & (mask_huckel)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55fd777e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "r = all_have_smiles.swifter.apply(\n",
    "    lambda x: get_agreement_between_multiple_smiles(\n",
    "        x.smiles_CSD_fix,\n",
    "        x.smiles_NBO_DFT_xyz,\n",
    "        x.smiles_huckel_DFT_xyz,\n",
    "    ),\n",
    "    axis=1,\n",
    ")\n",
    "\n",
    "print_to_file(f\"Comparing the SMILES of all three methods: {sum(r)}/{len(all_have_smiles)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "167a5858",
   "metadata": {},
   "outputs": [],
   "source": [
    "comp1 = combined_frame[mask_csd & mask_tmqmg]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf383872",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = comp1.swifter.apply(\n",
    "    lambda x: get_agreement_between_smiles(\n",
    "        x.smiles_CSD_fix,\n",
    "        x.smiles_NBO_DFT_xyz,\n",
    "    ),\n",
    "    axis=1,\n",
    ")\n",
    "print_to_file(f\"Comparing smiles directly CSD/NBO: {sum(r)}/{len(comp1)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d9c45ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "comp2 = combined_frame[mask_csd & mask_huckel]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95141d64",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = comp2.swifter.apply(\n",
    "    lambda x: get_agreement_between_smiles(x.smiles_CSD_fix, x.smiles_huckel_DFT_xyz),\n",
    "    axis=1,\n",
    ")\n",
    "print_to_file(f\"Comparing smiles directly CSD/Huckel: {sum(r)}/{len(comp2)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2586221",
   "metadata": {},
   "outputs": [],
   "source": [
    "comp3 = combined_frame[mask_tmqmg & mask_huckel]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d7d44dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "r = comp3.swifter.apply(\n",
    "    lambda x: get_agreement_between_smiles(x.smiles_NBO_DFT_xyz, x.smiles_huckel_DFT_xyz),\n",
    "    axis=1,\n",
    ")\n",
    "print_to_file(f\"Comparing smiles directly NBO/Huckel: {sum(r)}/{len(comp3)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96781cfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40b88058",
   "metadata": {},
   "source": [
    "## Comparing SMILES with resonance forms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef7e493",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [\n",
    "    (sm1, sm2)\n",
    "    for sm1, sm2 in zip(\n",
    "        comp1.smiles_CSD_fix.to_list(),\n",
    "        comp1.smiles_NBO_DFT_xyz.to_list(),\n",
    "    )\n",
    "]\n",
    "results = process_paralell(check_smiles_tmc_resonance, arguments)\n",
    "print_to_file(f\"Comparing smiles with resonance CSD/NBO: {count(results)}/{len(results)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1090ca9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [(sm1, sm2) for sm1, sm2 in zip(comp2.smiles_CSD_fix.to_list(), comp2.smiles_huckel_DFT_xyz.to_list())]\n",
    "results = process_paralell(check_smiles_tmc_resonance, arguments)\n",
    "print_to_file(f\"Comparing smiles with resonance CSD/Huckel: {count(results)}/{len(results)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c1ce728",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [(sm1, sm2) for sm1, sm2 in zip(comp3.smiles_NBO_DFT_xyz.to_list(), comp3.smiles_huckel_DFT_xyz.to_list())]\n",
    "results = process_paralell(check_smiles_tmc_resonance, arguments)\n",
    "print_to_file(f\"Comparing smiles with resonance NBO/Huckel: {count(results)}/{len(results)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2ca2b79",
   "metadata": {},
   "source": [
    "## Compare SMILES with metal disconnector and resonance forms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da870522",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [\n",
    "    (sm1, sm2)\n",
    "    for sm1, sm2 in zip(\n",
    "        comp1.smiles_CSD_fix.to_list(),\n",
    "        comp1.smiles_NBO_DFT_xyz.to_list(),\n",
    "    )\n",
    "]\n",
    "\n",
    "results = process_paralell(check_smiles_fragments_resonance, arguments)\n",
    "print_to_file(f\"Comparing smiles with disconnected resonances CSD/NBO: {count(results)}/{len(results)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea7ca02",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [(sm1, sm2) for sm1, sm2 in zip(comp2.smiles_CSD_fix.to_list(), comp2.smiles_huckel_DFT_xyz.to_list())]\n",
    "results = process_paralell(check_smiles_fragments_resonance, arguments)\n",
    "print_to_file(f\"Comparing smiles with disconnected resonances CSD/Huckel: {count(results)}/{len(results)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50693d16",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [(sm1, sm2) for sm1, sm2 in zip(comp3.smiles_NBO_DFT_xyz.to_list(), comp3.smiles_huckel_DFT_xyz.to_list())]\n",
    "results = process_paralell(check_smiles_fragments_resonance, arguments)\n",
    "print_to_file(f\"Comparing smiles with disconnected resonances NBO/Huckel: {count(results)}/{len(results)}\", \"comparison.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1978b3fa",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "tmc_smiles_surt",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}