{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Synopsis\n\nThe purpose of this kernel to provide a simple *annotated* visualization baseline code for molecules.\n\nAdditional clickable 3D molecule views have been added. They use the python `ase` library which can be installed with pip, see [How To: Easy Visualization of Molecules](https://www.kaggle.com/borisdee/how-to-easy-visualization-of-molecules/comments). \n\n## Changelog\n\n* v7: Ionized atoms groups support with `NH3+` and 1.5 bond in `COO-` groups.\n* v4: Bonds have been added using [Dataset with number of bonds between atoms](https://www.kaggle.com/asauve/dataset-with-number-of-bonds-between-atoms)\n\n# General information\n\n![](https://storage.googleapis.com/kaggle-competitions/kaggle/14313/logos/thumb76_76.png?t=2019-05-16-16-56-19)\n\nThis kernel uses data from [Predicting Molecular Properties](https://www.kaggle.com/c/champs-scalar-coupling) which is intended to predict interactions between atoms in the domain of Nuclar Magnetic Resonnance (NMR). More precisely it is the scalar coupling constant between atoms which is to be predicted.\n\nAs this challenge is based uppon molecular topological properties, it can be useful to have an appropriate way of representing molecules.\nThese visualizations can then be used to infer useful hints to understand the coupling properties, engineer appropriate features and debug\nprediction failures.\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Load data"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nfrom ase import Atoms\nimport ase.visualize  # clickable 3D molecule viewer    # pip install ase\n\nfrom mpl_toolkits.mplot3d import Axes3D\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set()\n\nimport os\n!ls -al --color ../input","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_dir_csv(directory):\n    csv_files = sorted( [ f for f in os.listdir(directory) if f.endswith(\".csv\") ])    \n    csv_vars  = [ filename[:-4] for filename in csv_files ]\n    gdict = globals()\n    for filename, var in zip( csv_files, csv_vars ):\n        print(f\"{var:32s} = pd.read_csv({directory}/{filename})\")\n        gdict[var] = pd.read_csv( f\"{directory}/{filename}\" )\n        print(f\"{'nb of rows ':32s} = \" + str(len(gdict[var])))\n        display(gdict[var].head())\n\nload_dir_csv(\"../input/champs-scalar-coupling\")\nload_dir_csv(\"../input/predicting-molecular-properties-bonds/\")\n                       ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Structure of molecules\n\n## Molecule geometry\n\nTo begin with, the location of atoms can be found in `structures.csv`. The columns of interest are\n* atom : a letter for the atom type\n* x, y, z: 3D coordinrates of each atom"},{"metadata":{},"cell_type":"markdown","source":"## Atoms types"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.atom.unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok there are only four different types of atoms. \nLet use the standard [CPK](https://fr.wikipedia.org/wiki/Code_de_couleurs_CPK) color code for their representation:  \n\n<img src=\"https://images-na.ssl-images-amazon.com/images/I/610rjiH9f5L._SL1500_.jpg\" width=\"50%\" />"},{"metadata":{},"cell_type":"markdown","source":"## Coupling types"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.type.unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are only 8 different possible couplings types. The coupling type starts with a number `1`, `2` or `3` and is relative to the number of bonds involved.\nIt will be plotted with\n* black (transparent) line for 1 bond coupling\n* green dashed (transparent) line for two bonds coupling\n* red dotted (transparent) line for three bond coupling\n"},{"metadata":{},"cell_type":"markdown","source":"## Scalar coupling constant\n\nthis is the value which has to be predicted, it can be found in `train.sc`\n* coupling is shown as a transparent connection between coupled atoms, which thickness correspond to the scalar value\n\n## Bonds\n\nA special dataset has been calculated from molecule topology and most probable covalent bondings in [Dataset with number of bonds between atoms](https://www.kaggle.com/asauve/dataset-with-number-of-bonds-between-atoms).\n\nBonds are shown as\n* thick black lines for 1-bonds\n* thick*1.5 dark green for 1.5-bonds found in COO- groups\n* thick*2 green lines for 2-bonds\n* thick*3 red lines for 3-bonds\n\n## Visualization code"},{"metadata":{},"cell_type":"markdown","source":"\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"def view3d_molecule(name, xsize=\"200px\", ysize=\"200px\"):\n    \"\"\"Mouse clickeble 3D view\"\"\"\n    m = structures[structures.molecule_name == name]\n    positions = m[['x','y','z']].values\n    v = ase.visualize.view(Atoms(positions=positions, symbols=m.atom.values), \n                           viewer=\"x3d\") \n    return v\n\ncpk = { \n    'C': (\"black\", 2),\n    'H': (\"white\", 1),\n    'O': (\"red\",   2),\n    'N': (\"dodgerblue\", 2),\n    'F': (\"green\", 2) }\n\nbond_colors = {'1.0':'black', '1.5':'darkgreen', '2.0':'green', '3.0':'red'}\n\ndef bond_type_to_pair(bond_type):\n    return bond_type[3:]\ndef bond_type_to_n(bond_type):\n    return bond_type[0:3]\n\ndef plot_molecule(name, ax=None, bonds=None, charges=None, elev=0, azim=-60):\n    \"\"\"bonds = if provided add bonds display from the bond table dataset in https://www.kaggle.com/asauve/predicting-molecular-properties-bonds\n    elev = 3D elevation angle [degree] for the molecule view\n    azim = 3D azimut angle [degree]\n    \"\"\"\n    if not ax:\n        fig = plt.figure()\n        ax = fig.add_subplot(111, projection='3d')\n    if (elev != 0) or (azim != -60):\n        ax.view_init(elev=elev, azim=azim)\n    \n    # atoms location\n    m = structures[structures.molecule_name == name].sort_values(by='atom_index')\n    if (charges is not None):\n        charges = charges[charges.molecule_name == name].sort_values(by='atom_index')\n        if len(charges) != len(m):\n            print(f\"Warning bad charges data for molecule {name}\")\n    \n    # formula\n    acount = {a : 0 for a in cpk}\n    for a in m.atom:\n        acount[a] += 1\n    formula = \"\"\n    for a in acount:\n        if acount[a] == 1:\n            formula += a\n        elif acount[a] > 1:\n            formula += \"%s_{%d}\" % (a, acount[a])\n\n    ax.set_title(f'{name} ${formula}$')\n    \n    # display couplings (coupling is not molecular bonds!)\n    couples = train[train.molecule_name == name][['atom_index_0', 'atom_index_1', 'type', 'scalar_coupling_constant']]\n    for c in couples.itertuples():\n        m0 = m[m.atom_index == c.atom_index_0]\n        m1 = m[m.atom_index == c.atom_index_1]\n        ax.plot([float(m0.x), float(m1.x)],[float(m0.y), float(m1.y)],[float(m0.z), float(m1.z)],\n               linestyle = ['', '-', '--', 'dotted'][int(c.type[0])],\n               color     = ['', 'black', 'green', 'red' ][int(c.type[0])],\n               linewidth = abs(float(c.scalar_coupling_constant))/5,\n               alpha=0.2), \n    if bonds is not None:\n        for b in bonds[bonds.molecule_name == name].itertuples():\n            m0 = m[m.atom_index == b.atom_index_0]\n            m1 = m[m.atom_index == b.atom_index_1]\n            ax.plot([float(m0.x), float(m1.x)],[float(m0.y), float(m1.y)],[float(m0.z), float(m1.z)], 'black', \n                    linewidth=2*float(b.nbond),\n                    color=bond_colors[bond_type_to_n(b.bond_type)])\n            \n    # display atoms\n    ax.scatter(m.x, m.y, m.z, c=[cpk[a][0] for a in m.atom], s=[100*cpk[a][1] for a in m.atom], edgecolor='black')\n        \n    # display atom index and charges\n    for row in m.itertuples():\n        x = float(row.x) - 0.15 if row.x > ax.get_xlim()[0] + 0.15 else row.x\n        y = float(row.y) - 0.15 if row.y > ax.get_ylim()[0] + 0.15 else row.y\n        z = float(row.z) - 0.15 if row.z > ax.get_zlim()[0] + 0.15 else row.z\n        ax.text(x, y, z, str(row.atom_index), color='darkviolet')\n        if charges is not None:\n            ch = float(charges[charges.atom_index == row.atom_index].charge)\n            if ch != 0:\n                x = float(row.x) + 0.15 if row.x < ax.get_xlim()[1] - 0.15 else row.x\n                y = float(row.y) + 0.15 if row.y > ax.get_ylim()[1] - 0.15 else row.y\n                z = float(row.z) + 0.15 if row.z < ax.get_zlim()[1] - 0.15 else row.z\n                ax.text(x, y, z, f\"{ch:+.1f}\", color='orangered' if ch > 0 else 'blue',\n                        bbox=dict(boxstyle='round', facecolor='white', alpha=0.5, \n                                  edgecolor='black'))\n                \nax = plot_molecule(\"dsgdb9nsd_000007\", bonds=train_bonds)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This example shows the three kind of *scalar coupling*.   \nThe 1 bond coupling having the largest coupling constant (visible as the thickness of the bond)."},{"metadata":{"trusted":true},"cell_type":"code","source":"view3d_molecule(\"dsgdb9nsd_000007\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is the very same molecule but with the non-annotated `ase`  3D  view.   \nBoth visualizatons are complementary for complex molecules."},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{},"cell_type":"markdown","source":"# Visualization of the 20 first molecules\n\nNow let show the 20 first molecules. Most of them are common.\n\n* dsgdb9nsd_000001: $CH_4$ Methane \n* dsgdb9nsd_000002: $NH_3$ Ammonia\n* dsgdb9nsd_000003: $H_2O$ Water :-) \n* dsgdb9nsd_000005: $HCN$ Hydrogen cyanide (aka the infamous zyklon B) contains a triple `3CC` bond and exhibit a very high scalar coupling for the `H` on the other side\n* dsgdb9nsd_000007: $C_2H_6$ Ethane \n* dsgdb9nsd_000008: $CH-3OH$ Methanol\n* dsgdb9nsd_000009: $C_3H_4$ Propyne contains a triple `3CC` bond and has also a very high scalar coupling for the `H` on the other side\n* dsgdb9nsd_000010: $CH_3CN$ Acetonitrile\n* dsgdb9nsd_000011: $CH_3CHO$ Ethanal or Acetaldehyde\n* dsgdb9nsd_000012: $HCONH_2$  Formamide or Methanamid\n* dsgdb9nsd_000013: $C_3H_8$ Propane\n* dsgdb9nsd_000014: $C_2H_5OH$ Alcool :-) \n* dsgdb9nsd_000017: $C_2H_4O$ Ethylene oxide or oxyrane contains a CCO cycle\n* dsgdb9nsd_000018: $C_3H_6O$ Aceton\n* dsgdb9nsd_000019: $CH_3CONH_2$ Acetamide or Ethanamide\n* dsgdb9nsd_000021: $C_4H_{10}$ Butane\n* dsgdb9nsd_000023: $C4_H2$ Diacetylene or Butadiyne has a very uncommon linear structure with two triple `3CC` bonds and very high scalar coupling\n* dsgdb9nsd_000026: $HC_2CHO$ Propynal contains also a triple bond and a very high scalar coupling\n* dsgdb9nsd_000027: $HCOCN$ Formyl cyanide has an uncomon triple bond `3CN` and has also a particularly high 2-bond scalar coupling `2JHC`\n* dsgdb9nsd_000028: $OCHCHO$ Glyoxal\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"nrow = 5\nncol = 4\nfig = plt.figure(figsize=(20, 20))\nmolecules = train.molecule_name.unique()\nfor i in range(nrow*ncol):\n    ax = fig.add_subplot(nrow, ncol, 1+i, projection='3d')\n    plot_molecule(molecules[i], ax=ax, bonds=train_bonds)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The butadyne - see bellow - is an interresting case because of its inline structure and high coupling. The triple bonds effect are visible because `C` atoms are closer in `3CC` bonds."},{"metadata":{"trusted":true},"cell_type":"code","source":"view3d_molecule(\"dsgdb9nsd_000023\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Ionized molecules\n\nA few (about 3% in each set) of total molecules are ionized, here are some examples of these beasts."},{"metadata":{"trusted":true},"cell_type":"code","source":"ionized = train_charges[(train_charges.charge != 0)].molecule_name.unique()\n\n# filter out molecules with failed bonding\nerrors  = train_bonds[train_bonds.error == 1].molecule_name.unique()\nerrors  = {e:1 for e in errors} # convert to dict for fast lookup\nionized = [ name for name in ionized if not name in errors]\n\nnrow = 3\nncol = 3\nfig = plt.figure(figsize=(18, 18))\nfor i in range(nrow*ncol):\n    ax = fig.add_subplot(nrow, ncol, 1+i, projection='3d')\n    plot_molecule(ionized[i], ax=ax, bonds=train_bonds, charges=train_charges)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Currently the two ionized groups supported are `COO-` and `NH3+`. The negative charge, at this time is supposed to be shared by both atoms of the `COO-` group.\n\nOne of these molecule is shown bellow for appreciationg the absence of `H` atom on the `COO-` group.\nThere is one extra `H` on the `NH3+` group that is further away than the other two `H` (supposedly building regular bonds)."},{"metadata":{"trusted":true},"cell_type":"code","source":"view3d_molecule(\"dsgdb9nsd_076394\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# What to do next ?\n\nThe first run is really instructive and shows for example that inline HCCH structures with  `3CC` triple bonds produce a strong scalar value. \nHence a potentially interesting directions would be to\n* study the relation of the scalar constant with the number of bonds on the nearby atoms\n* build space oriented features\n* study cycles\n\nHave fun!\n"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}