{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Molecules EDA"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":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)\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport seaborn as sns\nsns.set()\n\nimport os\nprint(os.listdir(\"../input\"))\n\nimport plotly.offline as py\npy.init_notebook_mode(connected=True)\nimport plotly.graph_objs as go","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(\"../input/train.csv\")\ntest = pd.read_csv(\"../input/test.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## What is given by the molecule structures?"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"structures = pd.read_csv(\"../input/structures.csv\")\nstructures.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### How are atom types distributed in the space?"},{"metadata":{},"cell_type":"markdown","source":"Let's take a look at the spacial distribution of the atoms given simple scatter plots:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"M = 8000\nfig, ax = plt.subplots(1,3,figsize=(20,5))\n\ncolors = [\"darkred\", \"dodgerblue\", \"mediumseagreen\", \"gold\", \"purple\"]\natoms = structures.atom.unique()\n\nfor n in range(len(atoms)):\n\n    ax[0].scatter(structures.loc[structures.atom==atoms[n]].x.values[0:M],\n                  structures.loc[structures.atom==atoms[n]].y.values[0:M],\n                  color=colors[n], s=2, alpha=0.5, label=atoms[n])\n    ax[0].legend()\n    ax[0].set_xlabel(\"x\")\n    ax[0].set_ylabel(\"y\")\n    \n    ax[1].scatter(structures.loc[structures.atom==atoms[n]].x.values[0:M],\n                  structures.loc[structures.atom==atoms[n]].z.values[0:M],\n                  color=colors[n], s=2, alpha=0.5, label=atoms[n])\n    ax[1].legend()\n    ax[1].set_xlabel(\"x\")\n    ax[1].set_ylabel(\"z\")\n    \n    ax[2].scatter(structures.loc[structures.atom==atoms[n]].y.values[0:M],\n                  structures.loc[structures.atom==atoms[n]].z.values[0:M],\n                  color=colors[n], s=2, alpha=0.5, label=atoms[n])\n    ax[2].legend()\n    ax[2].set_xlabel(\"y\")\n    ax[2].set_ylabel(\"z\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Is there some kind of order how the atoms were added to the structures file?\n\nLet's color by the index:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"M = 200000\nfig, ax = plt.subplots(1,3,figsize=(20,5))\n\nax[0].scatter(structures.x.values[0:M],\n              structures.y.values[0:M],\n              c = structures.index.values[0:M],\n              s=2, alpha=0.5, cmap=\"magma\")\nax[0].set_xlabel(\"x\")\nax[0].set_xlabel(\"y\");\n\nax[1].scatter(structures.x.values[0:M],\n              structures.z.values[0:M],\n              c = structures.index.values[0:M],\n              s=2, alpha=0.5, cmap=\"magma\")\nax[1].set_xlabel(\"x\")\nax[1].set_xlabel(\"z\");\n\nax[2].scatter(structures.y.values[0:M],\n              structures.z.values[0:M],\n              c = structures.index.values[0:M],\n              s=2, alpha=0.5, cmap=\"magma\")\nax[2].set_xlabel(\"y\")\nax[2].set_xlabel(\"z\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This looks very interesting. It seems like there is some order. What do you think?"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"N = 100000\n\ntrace1 = go.Scatter3d(\n    x=structures.x.values[0:N], \n    y=structures.y.values[0:N],\n    z=structures.z.values[0:N],\n    mode='markers',\n    marker=dict(\n        color=structures.index.values[0:N],\n        colorscale = \"YlGnBu\",\n        opacity=0.3,\n        size=1\n    )\n)\n\nfigure_data = [trace1]\nlayout = go.Layout(\n    title = 'Coloring by index',\n    scene = dict(\n        xaxis = dict(title=\"x\"),\n        yaxis = dict(title=\"y\"),\n        zaxis = dict(title=\"z\"),\n    ),\n    margin=dict(\n        l=0,\n        r=0,\n        b=0,\n        t=0\n    ),\n    showlegend=True\n)\n\nfig = go.Figure(data=figure_data, layout=layout)\npy.iplot(fig, filename='simple-3d-scatter')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Were larger molecules added later to the file?"},{"metadata":{"trusted":true},"cell_type":"code","source":"molecules_size = structures.groupby(\"molecule_name\").atom.size()\nmean_size = molecules_size.rolling(window=500).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nplt.plot(molecules_size.values, '+')\nplt.plot(np.arange(250, len(molecules_size)-249), mean_size.dropna().values, '-')\nplt.xlabel(\"Molecule name number\")\nplt.ylabel(\"Number of atoms\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* It looks a bit like zick-zack. The sharp shifts at some points look very strange."},{"metadata":{},"cell_type":"markdown","source":"### Does the way atom types are ordered has something to do with the molecule size?"},{"metadata":{"trusted":true},"cell_type":"code","source":"structures[\"atom_count\"] = structures.groupby(\"molecule_name\").atom.transform(\"size\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"to_show = [\"atom_index\", \"atom_count\"]\nmy_choice = to_show[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"M = 10000\nfig, ax = plt.subplots(5,3, figsize=(20,25))\n\nfor n in range(len(atoms)):\n    ax[n,0].scatter(structures.loc[structures.atom==atoms[n]].x.values[0:M],\n                    structures.loc[structures.atom==atoms[n]].y.values[0:M],\n                    c=structures.loc[structures.atom==atoms[n]][my_choice].values[0:M],\n                    s=5, alpha=0.5, cmap=\"inferno\")\n    ax[n,0].set_title(atoms[n])\n    ax[n,0].set_xlabel(\"x\")\n    ax[n,0].set_ylabel(\"y\")\n    \n    ax[n,1].scatter(structures.loc[structures.atom==atoms[n]].x.values[0:M],\n                    structures.loc[structures.atom==atoms[n]].z.values[0:M],\n                    c=structures.loc[structures.atom==atoms[n]][my_choice].values[0:M],\n                    s=5, alpha=0.5, cmap=\"inferno\")\n    ax[n,1].set_title(atoms[n])\n    ax[n,1].set_xlabel(\"x\")\n    ax[n,1].set_ylabel(\"z\")\n    \n    ax[n,2].scatter(structures.loc[structures.atom==atoms[n]].y.values[0:M],\n                    structures.loc[structures.atom==atoms[n]].z.values[0:M],\n                    c=structures.loc[structures.atom==atoms[n]][my_choice].values[0:M],\n                    s=5, alpha=0.5, cmap=\"inferno\")\n    ax[n,2].set_title(atoms[n])\n    ax[n,2].set_xlabel(\"y\")\n    ax[n,2].set_ylabel(\"z\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Yes, the structure has a dependency on the molecule size, but it makes sense that this is not true for all points. Larger molecules have atoms close to the center 0 as well. Interestingly this is more true for C and H atoms.\n* If we color by the atom index of molecules this comes even more clear. Low index numbers are more centered whereas higher indices live more outside of the \"sphere\".\n* **Consequently the order we have found has definitely something to do with the size of the molecules**. The next question is now, why are molecules of similar or same size orientated in different angles? We can clearly see a symmetric order with respect to angles as well."},{"metadata":{},"cell_type":"markdown","source":"### What about the mulliken charges?"},{"metadata":{"trusted":true},"cell_type":"code","source":"mulliken = pd.read_csv(\"../input/mulliken_charges.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"train_structures = structures.loc[structures.molecule_name.isin(train.molecule_name)].copy()\ntrain_structures = train_structures.set_index([\"molecule_name\", \"atom_index\"])\nmulliken = mulliken.set_index([\"molecule_name\", \"atom_index\"])\nmulliken[\"atom\"] = train_structures.atom\n\nplt.figure(figsize=(20,5))\nsns.violinplot(x=mulliken.atom, y=mulliken.mulliken_charge)\nplt.title(\"How are the mulliken charges distributed given the atom type?\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For H, N, O and F atoms the charges have more sharp distributions with multiple, but clear modes in contrast to the C atom type. "},{"metadata":{"trusted":true},"cell_type":"code","source":"train_structures[\"mulliken\"] = mulliken.mulliken_charge\ntrain_structures.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"M = 100000\nfig, ax = plt.subplots(5,1,figsize=(20,100))\n\nfor n in range(len(atoms)):\n    try:\n        ax[n].scatter(train_structures.loc[train_structures.atom==atoms[n]].x.values[0:M],\n                      train_structures.loc[train_structures.atom==atoms[n]].z.values[0:M], s=1,alpha=0.5,\n                      c=train_structures.loc[train_structures.atom==atoms[n]].mulliken.values[0:M],\n                      cmap=\"coolwarm\")\n    except:\n        ax[n].scatter(train_structures.loc[train_structures.atom==atoms[n]].x.values,\n                      train_structures.loc[train_structures.atom==atoms[n]].z.values, s=1, alpha=0.5,\n                      c=train_structures.loc[train_structures.atom==atoms[n]].mulliken.values,\n                      cmap=\"coolwarm\")\n    ax[n].set_title(\"Mulliken charges of \" + str(atoms[n]))\n    ax[n].set_xlabel(\"x\")\n    ax[n].set_ylabel(\"z\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that there is a clear radial dependency of the charge as well as a dependency with specific accumulations (angles). "},{"metadata":{},"cell_type":"markdown","source":"### What about dipole moments?"},{"metadata":{},"cell_type":"markdown","source":"The dipole moment tells us something about the way charges are distributed within a molecule. Some couplings like FH causes electron to move closer to one atom of the pair than to the otherone. In this case a small separation of charges has happened which will result in a non-zero dipole moment. [Like to know more?](https://en.wikipedia.org/wiki/Dipole#Molecular_dipoles)"},{"metadata":{"trusted":true},"cell_type":"code","source":"dipoles = pd.read_csv(\"../input/dipole_moments.csv\")\ndipoles.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check if we have this information given for test as well:"},{"metadata":{"trusted":true},"cell_type":"code","source":"test.molecule_name.isin(dipoles.molecule_name).sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"No!\n\nHence it's an information that we can use to create new features that we can somehow map on test. \n\n### How do the dipoles look like?\n\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(dipoles.X, label=\"X\")\nsns.distplot(dipoles.Y, label=\"Y\")\nsns.distplot(dipoles.Z, label=\"Z\")\nplt.legend();\nplt.title(\"Coordinate specific dipole moment distribution\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In contrast to X the coordinates Y and Z show a bimodality. Besides that we can see extreme outliers. Let's check the absolute sum of dipole moments:"},{"metadata":{"trusted":true},"cell_type":"code","source":"dipoles[\"total\"] = np.abs(dipoles.X) + np.abs(dipoles.Y) + np.abs(dipoles.Z)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nsns.distplot(dipoles.total);\nplt.title(\"How is the absolute sum of X, Y, Z dipole moments distributed?\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There seem to be a few molecules that have almost no separation of charges within the molecule. The extreme outliers indicate some exotics with very high charge separation and consequently a very high dipole moment."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(3,1,figsize=(20,60))\nax[0].scatter(dipoles.X.values, dipoles.Y.values, s=0.5, alpha=0.1);\nax[0].set_xlabel(\"X\")\nax[0].set_ylabel(\"Y\")\nax[1].scatter(dipoles.X.values, dipoles.Z.values, s=0.5, alpha=0.1);\nax[1].set_xlabel(\"X\")\nax[1].set_ylabel(\"Z\")\nax[2].scatter(dipoles.Y.values, dipoles.Z.values, s=0.5, alpha=0.1);\nax[2].set_xlabel(\"Y\")\nax[2].set_ylabel(\"Z\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Potential energy"},{"metadata":{"trusted":true},"cell_type":"code","source":"energy = pd.read_csv(\"../input/potential_energy.csv\")\nenergy.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"atom_counts = structures[structures.molecule_name.isin(dipoles.molecule_name)].groupby(\"molecule_name\").atom.size()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max_mulliken = mulliken.groupby(\"molecule_name\").mulliken_charge.max()\nmin_mulliken = mulliken.groupby(\"molecule_name\").mulliken_charge.min()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,5))\nax[0].scatter(energy.potential_energy.values, dipoles.total.values, s=1, c=atom_counts.values, cmap=\"YlGnBu\");\nax[0].set_xlabel(\"Potential energy\")\nax[0].set_ylabel(\"Total absolute dipole moment\")\nax[0].set_title(\"Colored by atom counts of a molecule\");\n\nax[1].scatter(energy.potential_energy.values, max_mulliken.values, c=atom_counts.values, cmap=\"YlGnBu\", s=1);\nax[1].set_xlabel(\"Potential energy\")\nax[1].set_ylabel(\"Max mulliken charge\");\n\nax[2].scatter(energy.potential_energy.values, min_mulliken.values, c=atom_counts.values, cmap=\"YlGnBu\", s=1);\nax[2].set_xlabel(\"Potential energy\")\nax[2].set_ylabel(\"Min mulliken charge\");","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}