{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Predicting Molecular Properties - Theory and EDA"},{"metadata":{},"cell_type":"markdown","source":"# Table of Contents\n1. [Background](#Background)\n\n2. [Loading Modules and Data](#Loading-Modules-and-Data)\n\n3. [Exploratory Data Analysis](#Exploratory-Data-Analysis)"},{"metadata":{},"cell_type":"markdown","source":"### Background"},{"metadata":{},"cell_type":"markdown","source":"In this competition, we will develop an algorithm that can predict the magnetic interaction between two atoms in a molecule (i.e., the scalar coupling constant). "},{"metadata":{},"cell_type":"markdown","source":"A qualitative understanding of the scalar coupling constant and mechanics behind it may help guide in engineering features that succinctly capture the interactions between the atoms."},{"metadata":{},"cell_type":"markdown","source":"A series of snippets from a number of different sources lay out the qualitative characteristics of scalar coupling based on the theoretical background:\n\n> - Scalar couplings arise from spin-spin interactions that occur via bonding electrons. \n    - More specifically, it arises from the interaction of the nuclear magnetic moment with the electrons involved in the chemical bond. The nuclear spin polarization of one atom affects the polarization of the surrounding electrons. The electron polarization subsequently produces a change in the magnetic field that is sensed by the coupled spin.\n    - Consequently, they provide information on the chemical connectivity between atoms. \n    -  the sizes of three bond scalar couplings are sensitive to the electron distribution of the intervening bonds.\n    - The nomenclature that is used to describe the coupling is as follows: \"nJAB\" where n refers to the number of intervening bonds, and A and B identify the two coupled spins. \n[(source)](https://www.springer.com/gp/book/9781402034992)\n\n> As the\ncoupling is mediated by bonding electrons. it provides important information about the constitution of molecules in terms of\nthe connectivity of thc coupled nuclei. The size of the coupling\nconstant depends not only on the number of bonds that separate\nthe coupled nuclei, but also on the configuration of the electrons\nand their spatial arrangement. [(source)](https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.199516711)\n\n> The most important contribution to the scalar coupling is the Fermi contact term.  This term relies on the probability of finding an electron at the site of the two coupled nuclei.  With this term, it is therefore expected that s electrons will play a very significant role since these are the only electrons that do not have nodes at the nuclear sites. [(source)](https://bouman.chem.georgetown.edu/nmr/scalar/scalar.htm)\n\n> More on Fermi Contact: The scalar interaction arises between two different\nnuclear spins, I1 and I2, and is mediated by the\nelectrons surrounding these two spins. Through the\nFermi contact, the electrons are polarised in the\nopposite direction to the nucleus they are interacting\nwith. This polarisation in turn has an effect on the\nother electrons in close proximity, which in turn\naffects the neighbouring nuclei.[(source)](https://groups.chem.ubc.ca/straus/l3.pdf)\n\n> The Karplus equation, named after Martin Karplus, describes the correlation between $^3J$-coupling constants and dihedral torsion angles: \n\\begin{equation*} J(\\sigma) = Ccos(2\\sigma) + Bcos(\\sigma) + A\\end{equation*}\nwhere J is the $^3J$ coupling constant,  $\\phi$  is the dihedral angle, and $A$, $B$, and $C$ are empirically derived parameters whose values depend on the atoms and substituents involved.[(source)](https://en.wikipedia.org/wiki/Karplus_equation)\n\n> Atom-centered Symmetry Functions:  \nIn chemistry, the most common way of representing molecules is to use Cartesian coordinates. The problem of Cartesian coordinates is that if one rotates or translates a molecule, the coordinates change.  However, properties like the total energy of the molecule, the dipole moment or the atomisation energy remain unchanged.  Therefore, if you want to learn the properties of a molecule with a neural network, you want to represent the molecular structure in a way that doesn’t change when you rotate or translate the molecule.  \nBy doing this, you make the training process more efficient because the neural network doesn’t need to learn that many different Cartesian coordinates represent the same structure.  This is where Atom Centred Symmetry Functions (ACSF) come in handy. They are a way of representing a molecular structure which remains the same when the molecule is rotated or translated. [(source)](https://samabilino.wordpress.com/2018/07/07/atom-centred-symmetry-functions/)\n"},{"metadata":{},"cell_type":"markdown","source":"\n### Loading Modules and Data"},{"metadata":{},"cell_type":"markdown","source":"Load data analysis modules we may require:"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\n#Load data analysis/plotting modules\nfrom tqdm import tqdm_notebook as tqdm\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport pandas as pd\nimport numpy as np\nimport plotly\nimport sympy\nimport scipy\nimport numpy\nimport os\n\n\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"list the csv files for the competition:"},{"metadata":{"trusted":true},"cell_type":"code","source":"os.listdir('../input/champs-scalar-coupling/')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading all the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"PATH = \"../input/champs-scalar-coupling/\"\ntrain = pd.read_csv(PATH + \"train.csv\")\ntest = pd.read_csv(PATH + \"test.csv\")\nstruct = pd.read_csv(PATH + \"structures.csv\")\ndpm = pd.read_csv(PATH + \"dipole_moments.csv\")\nmst = pd.read_csv(PATH + \"magnetic_shielding_tensors.csv\")\nmlk = pd.read_csv(PATH + \"mulliken_charges.csv\")\npe = pd.read_csv(PATH + \"potential_energy.csv\")\nscc = pd.read_csv(PATH + \"scalar_coupling_contributions.csv\")\nsmpsub = pd.read_csv(PATH + \"sample_submission.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n### Exploratory Data Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"#preview train,test, and structure dataframes\ndisplay(train.head(), train.shape)\n\ndisplay(test.head(), test.shape)\n\ndisplay(struct.head(), struct.shape)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train.molecule_name.nunique(), \"molecules in train dataset.\")\nprint(test.molecule_name.nunique(), \"molecules in test dataset.\")\nprint(struct.molecule_name.nunique(), \"molecules in structures dataset.\")\n# structures is the sum of both train+test molecules ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nratio = (test.shape[0]/(train.shape[0]+test.shape[0]))\nprint('ratio:', 1-ratio) \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have roughly a 60/40 split for the train/test datasets.\nThe notation for the type of coupling follows the method mentioned in the background section.   \nLet's look at the distribution of types between the train/test dataset:"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train.type.unique()) # 8 different coupling types (between carbon, nitrogen, and hydrogen)\nprint(test.type.unique())\nprint(set(train.type.unique()) == set(test.type.unique())) #same types exist in both train/test datasets","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f,ax=plt.subplots(1,2,figsize=(15,5))\ntrain.type.value_counts().plot.bar(ax=ax[0], color = 'blue')\nax[0].set_title('Train: Number of Coupling Types')\ntest.type.value_counts().plot.bar(ax=ax[1], color = 'red')\nax[1].set_title('Test: Number of Coupling Types')\nplt.show()\nprint('train:')\ndisplay(train.type.str[0].value_counts(), train.type.str[-2:].value_counts()) # number of coupling interactions for intervening bonds and atom pairs.\nprint('test:')\ndisplay(test.type.str[0].value_counts(), test.type.str[-2:].value_counts())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The datasets have a suprisingly similar distribution of coupling types.  This is good, since it means we won't have to worry about oversampling certain types of coupling when we train our model on the training dataset.  \nBut, intuitively it seems that coupling constants with 3 intervening bonds would be the most challenging to predict (due to long-range interactions) and seem to be the most abundant.\n\nH-C interactions also make up a disproportionately large number of coupling compared to the others. ( H-C >> H-H > H-N)\n\nLet's see what the distribution of the number of atoms per molecule looks like:"},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndata_train = struct[struct.molecule_name.isin(train.molecule_name.unique())].molecule_name.value_counts()\nbins_train = data_train.nunique()\ndata_test = struct[struct.molecule_name.isin(test.molecule_name.unique())].molecule_name.value_counts()\nbins_test = data_test.nunique()\n\nf,ax=plt.subplots(1,2,figsize=(15,5))\ndata_train.plot.hist(ax=ax[0], color = 'blue', bins=bins_train)\nax[0].set_title('Distribution of # of atoms per Molecule in Train')\ndata_test.plot.hist(ax=ax[1], color = 'red', bins=bins_test)\nax[1].set_title('Distribution of # of atoms per Molecule in Test')\nplt.show()\ndisplay(data_train.describe())\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Once again, it seems that the train/test split were made to have identical distributions in molecule size and coupling type.  \n\nIt is likely that the small deviation from a 60/40 split is due to the need to balance the coupling type/molecule size between the train/test datasets.  Maybe we can take this into consideration when we create the validation dataset.  \n\nLet's look at the coupling contributions with respect to the resultant coupling constant:"},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.concat([scc, train.scalar_coupling_constant], axis=1).head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The Fermi contact term accounts for most of the coupling constant contribution (as was mentioned in the background section).  Predicting this term would be the primary objective of the model.\n"},{"metadata":{},"cell_type":"markdown","source":"A great molecule visualization function was written by [Mykola Zotko](https://www.kaggle.com/mykolazotko/3d-visualization-of-magnetic-interactions):\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nimport plotly.graph_objs as go\nfrom plotly.offline import init_notebook_mode, iplot\nfrom sympy.geometry import Point3D\n\n\n# initiate the plotly notebook mode\ninit_notebook_mode(connected=True)\n    \n\ndef plot_interactions(molecule_name, structures, train_df):\n    \"\"\"Creates a 3D plot of the molecule\"\"\"\n\n\n    atomic_radii = dict(C=0.77, F=0.71, H=0.38, N=0.75, O=0.73)  \n    cpk_colors = dict(C='black', F='green', H='white', N='blue', O='red')\n    \n    if molecule_name not in train_df.molecule_name.unique():\n        print(f'Molecule \"{molecule_name}\" is not in the training set!')\n        return\n    \n    molecule = structures[structures.molecule_name == molecule_name]\n    coordinates = molecule[['x', 'y', 'z']].values\n    x_coordinates = coordinates[:, 0]\n    y_coordinates = coordinates[:, 1]\n    z_coordinates = coordinates[:, 2]\n    elements = molecule.atom.tolist()\n    radii = [atomic_radii[element] for element in elements]\n    \n    data_train = train_df[train_df.molecule_name == molecule_name][['atom_index_0', 'atom_index_1', 'scalar_coupling_constant']]\n    interactions = data_train.groupby('atom_index_0')['atom_index_1'].apply(set).to_dict()\n    coupling_constants = data_train.set_index(['atom_index_0', 'atom_index_1']).round(2).to_dict()['scalar_coupling_constant']\n    \n    def get_bonds():\n        \"\"\"Generates a set of bonds from atomic cartesian coordinates\"\"\"\n        ids = np.arange(coordinates.shape[0])\n        bonds = dict()\n        coordinates_compare, radii_compare, ids_compare = coordinates, radii, ids\n        \n        for _ in range(len(ids)):\n            coordinates_compare = np.roll(coordinates_compare, -1, axis=0)\n            radii_compare = np.roll(radii_compare, -1, axis=0)\n            ids_compare = np.roll(ids_compare, -1, axis=0)\n            distances = np.linalg.norm(coordinates - coordinates_compare, axis=1)\n            bond_distances = (radii + radii_compare) * 1.3\n            mask = np.logical_and(distances > 0.1, distances <  bond_distances)\n            distances = distances.round(2)\n            new_bonds = {frozenset([i, j]): dist for i, j, dist in zip(ids[mask], ids_compare[mask], distances[mask])}\n            bonds.update(new_bonds)\n        return bonds      \n            \n    def atom_trace():\n        \"\"\"Creates an atom trace for the plot\"\"\"\n        colors = [cpk_colors[element] for element in elements]\n        markers = dict(color=colors, line=dict(color='lightgray', width=2), size=7, symbol='circle', opacity=0.8)\n        trace = go.Scatter3d(x=x_coordinates, y=y_coordinates, z=z_coordinates, mode='markers', marker=markers,\n                             text=elements, name='')\n        return trace\n\n    def bond_trace():\n        \"\"\"\"Creates a bond trace for the plot\"\"\"\n        trace = go.Scatter3d(x=[], y=[], z=[], hoverinfo='none', mode='lines',\n                             marker=dict(color='grey', size=7, opacity=1), line=dict(width=5))\n        for i, j in bonds.keys():\n            trace['x'] += (x_coordinates[i], x_coordinates[j], None)\n            trace['y'] += (y_coordinates[i], y_coordinates[j], None)\n            trace['z'] += (z_coordinates[i], z_coordinates[j], None)\n        return trace\n    \n    def interaction_trace(atom_id):\n        \"\"\"\"Creates an interaction trace for the plot\"\"\"\n        trace = go.Scatter3d(x=[], y=[], z=[], hoverinfo='none', mode='lines',\n                             marker=dict(color='pink', size=7, opacity=0.5),\n                            visible=False)\n        for i in interactions[atom_id]:\n            trace['x'] += (x_coordinates[atom_id], x_coordinates[i], None)\n            trace['y'] += (y_coordinates[atom_id], y_coordinates[i], None)\n            trace['z'] += (z_coordinates[atom_id], z_coordinates[i], None)\n        return trace\n    \n    bonds = get_bonds()\n    \n    zipped = zip(range(len(elements)), x_coordinates, y_coordinates, z_coordinates)\n    annotations_id = [dict(text=num, x=x, y=y, z=z, showarrow=False, yshift=15, font = dict(color = \"blue\"))\n                      for num, x, y, z in zipped]\n    \n    annotations_length = []\n    for (i, j), dist in bonds.items():\n        p_i, p_j = Point3D(coordinates[i]), Point3D(coordinates[j])\n        p = p_i.midpoint(p_j)\n        annotation = dict(text=dist, x=float(p.x), y=float(p.y), z=float(p.z), showarrow=False, yshift=10)\n        annotations_length.append(annotation)\n    \n    annotations_interaction = []\n    for k, v in interactions.items():\n        annotations = []\n        for i in v:\n            p_i, p_j = Point3D(coordinates[k]), Point3D(coordinates[i])\n            p = p_i.midpoint(p_j)\n            constant = coupling_constants[(k, i)]\n            annotation = dict(text=constant, x=float(p.x), y=float(p.y), z=float(p.z), showarrow=False, yshift=25,\n                              font = dict(color = \"hotpink\"))\n            annotations.append(annotation)\n        annotations_interaction.append(annotations)\n    \n    buttons = []\n    for num, i in enumerate(interactions.keys()):\n        mask = [False] * len(interactions)\n        mask[num] = True\n        button = dict(label=f'Atom {i}',\n                      method='update',\n                      args=[{'visible': [True] * 2 + mask},\n                            {'scene.annotations': annotations_id + annotations_length + annotations_interaction[num]}])\n        buttons.append(button)\n        \n    updatemenus = list([\n        dict(buttons = buttons,\n             direction = 'down',\n             xanchor = 'left',\n             yanchor = 'top'\n            )\n    ])\n    \n    data = [atom_trace(), bond_trace()]\n    \n    # add interaction traces\n    for num, i in enumerate(interactions.keys()):\n        trace = interaction_trace(i)\n        if num == 0:\n            trace.visible = True \n        data.append(trace)\n        \n    axis_params = dict(showgrid=False, showticklabels=False, zeroline=False, titlefont=dict(color='white'))\n    layout = dict(scene=dict(xaxis=axis_params, yaxis=axis_params, zaxis=axis_params,\n                             annotations=annotations_id + annotations_length + annotations_interaction[0]),\n                  margin=dict(r=0, l=0, b=0, t=0), showlegend=False, updatemenus=updatemenus)\n\n    fig = go.Figure(data=data, layout=layout)\n    iplot(fig)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_interactions('dsgdb9nsd_000001', struct, train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A good start for generating deeper features.  Let's compile a list of structural properties that can be engineered as features based on what we can intuit from the theory.  This [source](https://www.ucl.ac.uk/nmr/NMR_lecture_notes/L3_3_97_web.pdf) has a great summary of the variables to consider.  Influentual variables can be placed in three groups :\n\n### 1.  The hybridization of the atoms involved in the coupling.\n\nThe hybridization of the atoms along the coupling pathway has an effect on the scalar coupling constant (scc).  For instance, a 3JHH coupling would have a different scc between the molecule H-C-C-H and H-C=C-H, (where '-' and '=' represent a single and double bond, respectively).\n    \n### 2. Dihedral/bond angles of the atoms involved in the coupling.\n\nThe angles between the atoms involved in the coupling is also an important variable.  For 1J coupling, the angle would be 0 since we just have two atom bonded directly together.  For 2J coupling it would be the bond angle between the two atoms and for 3J coupling it would be the dihedral angle.  The dihedral angle is the angle between the atoms with respect to a common rotational axis (usually a pair of carbon atoms). Looking down the axis of the angle in question would look like this ([source](https://www.ucl.ac.uk/nmr/NMR_lecture_notes/L3_3_97_web.pdf) for figures below):\n\n<img src='https://i.imgur.com/RKcjhFZ.png'>    \n    \n\n### 3. Substituent effects due to bonded/adjacent chemical elements.\n\nThe molecular structures attached to the primary coupling pathway must also be considered.\n    \nFor 1J and 2J couplings, the electronegativity of $\\alpha$ and $\\beta$ substituents attached adjacent to the coupling also have an influence on the scc.  The $\\alpha$ substituent being the first atom bonded to the coupling and the $\\beta$ substituent being the second, we can define a third as a $\\gamma$ substituent but the effect of its electronegativity is negligible.  It would be useful to chiefly consider the electronegativity of the  $\\alpha$ substituent and average the electronegativities of all elements attached to it.\n\nFor 3J couplings, the electronegativity is still a factor but now it requires a deeper analysis of its location/angle along the 3J coupling pathway.  For example, in the case of 3JHH coupling, the position of the OH substituent plays a role([source](https://www.ucl.ac.uk/nmr/NMR_lecture_notes/L3_3_97_web.pdf) for figures below):\n        \n<img src=\"https://i.imgur.com/necGAcQ.png\">    \n\n\nThe presence of $\\pi$ (double) bonds adjacent to the coupling also affects the scc.  In this case, we need to take into account the angle between the $\\pi$ bond and the bond for the coupling pathway in question, where the effect is maximized if they're parallel ([source](http://www.sliderbase.com/spitem-114-3.html) for figures below).\n<table><tr>\n<td> <img src=\"https://i.imgur.com/mWYLFIQ.png\" alt=\"Drawing\" style=\"width: 300px;\"/> </td> \n<td> <img src=\"https://i.imgur.com/32thny7.png\" alt=\"Drawing\" style=\"width: 300px;\"/> </td>\n</tr></table> \n\nThe complexity of the problem becomes apparent when we layout all of the variables we need to consider for an accurate prediction of the scc.  Most of the features are going to rely on a structural approach to the problem; we need to define the primary coupling pathway and the substituents bonded/adjacent to it seperately.  \n\nThese substituents then need to be characterized by their overall electronegativity and distance/angle w.r.t to the coupling pathway.  A similar approach would be applied to $\\pi$ bonds(not on the coupling pathway) but instead of the electronegativity, we could consider the distance/angle. \n\n"},{"metadata":{},"cell_type":"markdown","source":"**Note:** There's a *dist* feature that seems to be a significant predictor of the coupling coefficient.  Though I haven't begun generating features and training a model based on my EDA, it would be interesting to consider why the *dist* features does so well as a means of predicting other features.  Such as:\n\n- the average distance between coupling atoms is likely to increase as a result of a longer coupling pathway, ergo more intervening bonds(1J, 2J, 3J).  \n\n\n- is it possible that ignoring substituent effects and other characteristics that could describe the fc contribution to the scc more accurately is creating a plateau for max lb scores?\n\n\n- what would happen if we include the thereoretically derived features along with dist, would the model become more accurate or would it suffer due to the dist feature?\n\n\n- In the case above, it may be wise to remove the dist feature if we're experimenting with new/different features (LOFO analysis?)."}],"metadata":{"kernelspec":{"display_name":"Python 3","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.6.8"}},"nbformat":4,"nbformat_minor":1}