{"cells":[{"metadata":{},"cell_type":"markdown","source":"In a molecule, the coupling constants will depend, on a first order, on the geometry of the molecules: once the positions of the nuclei are fixed in space, \nthe electronic orbitals can be calculated resolving the shrodinger equation in the scope of the Born Oppenheiner approximation. From the potential energy surface obtained that way, coupling constants can be derived. Hence, as a first guess, I thought that an approach that would capture the global view of the geometry of the molecule could provide good results. In order to simplify the problem, I thought of reducing the complexity of the problem by integrating over a dimension of space, thus reducing the 3D geometrical structure of the molecule to a 2D image. Doing so, the problem could then be treated as a regression task on images. \n\nIn this notebook, I describe the way how such images are obtained.\n"},{"metadata":{},"cell_type":"markdown","source":"## 1. Description of the method"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install chart_studio ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Considering a vector $\\vec{d}$ along the molecular bond, we define a coordinate system aligned with $\\vec{d}$:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\n\nd = np.array([1, 2, 1])\nprint(f'with d = {d}')\n\n#_______________________\nnorm = np.linalg.norm(d)\nd = d / norm\n#______________________\nif d[2] != 0:\n    u = np.ones(3)\n    u[2] = -(d[0]*u[0] + d[1]*u[1]) / d[2]\n    norm = np.linalg.norm(u)\n    u = u / norm\nelif d[0] != 0:\n    u = np.ones(3)\n    u[0] = - d[1]*u[1] / d[0]\n    norm = np.linalg.norm(u)\n    u = u / norm\nelse:\n    u = np.ones(3)\n    u[1] = 0\n    norm = np.linalg.norm(u)\n    u = u / norm\n#____________________________________\nv = np.cross(d, u)\nnorm = np.linalg.norm(v)\nv = v / norm\nprint('the coordinate system will be defined through the vectors:'\n    '\\n', d,'\\n', u,'\\n', v)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using this coordinate system, we define 2 distances. The first, **R** corresponds to the distance along $\\vec{d}$ from the centroid of the molecular bond. The scond, **r**, defines the radius of a circle in a plane perpendicular to $\\vec{d}$. Defining a set of **(r, R)** points:"},{"metadata":{"trusted":true},"cell_type":"code","source":"values = [\n    [1, 5],\n    [2, 5],\n    [3, 5],\n    [1, 7],\n    [2, 7],\n    [3, 7]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"new_data = []\nfor r, R in values:\n    center_circle = R * d  # origine du cercle\n    #______________________________________________________________\n    val = []\n    for theta in np.linspace(0, 2 * np.pi):\n        xp = r * np.cos(theta)\n        yp = r * np.sin(theta)\n        val.append([xp, yp])\n    #________________________________\n    new_pt = []\n    for pt in val:\n        new_pt.append(center_circle + (pt[0] * u + pt[1] * v ))\n    new_data.append(np.array(new_pt))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We obtain the following set of points in space:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"from plotly.offline import init_notebook_mode, iplot\nimport pandas as pd\nimport chart_studio.plotly as py\nimport plotly.graph_objects as go\n\ndata = []\nclusters = []\ncolors = ['rgb(228,26,28)','rgb(55,126,184)','rgb(77,175,74)',\n          'rgb(155,126,184)','rgb(177,175,74)', 'rgb(77,175,174)',]\n\nfor i, pt in enumerate(new_data):\n    x = pt[:, 0]\n    y = pt[:, 1]\n    z = pt[:, 2]\n    trace = dict(\n            name = f'r={values[i][0]} R={values[i][1]}',\n            x = x, y = y, z = z,\n            type = \"scatter3d\",    \n            mode = 'markers',\n            marker = dict( size=3, color=colors[i], line=dict(width=0) ) )\n    data.append( trace )\n \n\nlayout = dict(\n    width=800,\n    height=550,\n    autosize=False,\n    title='Cone',\n    scene=dict(\n        xaxis=dict(\n            gridcolor='rgb(255, 255, 255)',\n            zerolinecolor='rgb(255, 255, 255)',\n            showbackground=True,\n            backgroundcolor='rgb(230, 230,230)'\n        ),\n        yaxis=dict(\n            gridcolor='rgb(255, 255, 255)',\n            zerolinecolor='rgb(255, 255, 255)',\n            showbackground=True,\n            backgroundcolor='rgb(230, 230,230)'\n        ),\n        zaxis=dict(\n            gridcolor='rgb(255, 255, 255)',\n            zerolinecolor='rgb(255, 255, 255)',\n            showbackground=True,\n            backgroundcolor='rgb(230, 230,230)'\n        ),\n        aspectratio = dict( x=1, y=1, z=1 ),\n        aspectmode = 'data'        \n    ),\n)\n\nfig = dict(data=data, layout=layout)\nurl = iplot(fig, filename='pandas-3d-iris', validate=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For every position defined that way, we calculate the density of atoms through\n$ \n \\sum_i 1/d_i^2\n$\n    where the sum if over all the atoms of a single type in the molecule. At this stage, in order to create the 2D images, we integrate over the circumference of the circles the density calculated above. At the end of the process, we thus have 5 *pseudo images* whose pixels are in the **(r, R)** coordinate system. "},{"metadata":{},"cell_type":"markdown","source":"## 2. Application to a chunk of 1JHN bonds \nBelow, I use the scheme defined above to some 1JHN bonds."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport pickle\nimport time\n\nstructures = pd.read_csv('../input/champs-scalar-coupling/structures.csv')\ntrain = pd.read_csv('../input/champs-scalar-coupling/test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"symetry = '1JHN'\ntrain = train[train['type'] == symetry]\nprint(f'selecting {train.shape[0]} bonds from {symetry}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def map_atom_info(df, atom_idx):\n    df = pd.merge(df, structures, how = 'left',\n                  left_on  = ['molecule_name', f'atom_index_{atom_idx}'],\n                  right_on = ['molecule_name',  'atom_index'])\n    \n    df = df.drop('atom_index', axis=1)\n    df = df.rename(columns={'atom': f'atom_{atom_idx}',\n                            'x': f'x_{atom_idx}',\n                            'y': f'y_{atom_idx}',\n                            'z': f'z_{atom_idx}'})\n    return df\n\ntrain = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntrain_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\nd_bond = train_p_0 - train_p_1\nnbonds = len(d_bond)\n\ntrain_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\ncenter = (train_p_0 + train_p_1) / 2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"we define a function that performs the integration along the circles:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def density_integration(df, nb_angles):\n    density_integrale = []\n    for icouple in range(df.shape[0]): \n\n        step_length = length[icouple]\n        p = df[icouple, :, :].copy()   \n        density = np.zeros(nb_angles)\n\n        for iatom in range(len(atoms)):\n            val = p - atoms[iatom]\n            dist = np.linalg.norm(val, axis=1)\n            density += 1 / dist**2\n\n        integrale = density.sum() * step_length\n        density_integrale.append(integrale)\n    return density_integrale","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We define images of size 35 x 35:"},{"metadata":{"trusted":true},"cell_type":"code","source":"nb_r = 35\nnb_R = 35\nnb_angles = 45\n\nvalues = [(r, R) for r in np.linspace(5, 0.01, nb_r) for R in np.linspace(-5, 5, nb_R)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"start_time = time.time()\n\nimage_file = 'images.pkl'\nbond_file = 'bond.pkl'\n\nbond_list = []\nimage_set = []\n\nids = train['id'].values\nfor ibond in range(nbonds): \n    #__________________\n    d = d_bond[ibond]\n    orig = center[ibond]\n\n    molec = train.iloc[ibond]['molecule_name']\n    bond_list.append(molec+'_'+str(ids[ibond]))\n\n    norm = np.linalg.norm(d)\n    d = d / norm\n    #______________________\n    if d[2] != 0:\n        u = np.ones(3)\n        u[2] = -(d[0]*u[0] + d[1]*u[1]) / d[2]\n        norm = np.linalg.norm(u)\n        u = u / norm\n    elif d[0] != 0:\n        u = np.ones(3)\n        u[0] = - d[1]*u[1] / d[0]\n        norm = np.linalg.norm(u)\n        u = u / norm\n    else:\n        u = np.ones(3)\n        u[1] = 0\n        norm = np.linalg.norm(u)\n        u = u / norm\n    #____________________________________\n    v = np.cross(d, u)\n    norm = np.linalg.norm(v)\n    v = v / norm\n    #_______________________________\n    new_data = []\n    length = []\n    for r, R in values:\n        center_circle = R * d  + orig \n        #______________________________________________________________\n        val = []\n        for theta in np.linspace(0, 2 * np.pi, nb_angles):\n            xp = r * np.cos(theta)\n            yp = r * np.sin(theta)\n            val.append([xp, yp])\n    #         print(theta)\n        #________________________________\n        new_pt = []\n        for pt in val:\n            new_pt.append(center_circle + (pt[0] * u + pt[1] * v ))\n        new_data.append(np.array(new_pt))\n\n        length.append(r * (2 * np.pi) / (nb_angles - 1))\n    #__________________________________\n    new_data = np.array(new_data)\n    #__________________________________\n    image = []\n    for atom in ['C', 'N', 'O', 'H', 'F']:\n        atoms = structures[(structures['molecule_name'] == molec) &\n                           (structures['atom'] == atom)][['x', 'y', 'z']].values\n\n        data = density_integration(new_data, nb_angles)\n        image.append(np.array(data).reshape(nb_r, nb_R))\n\n    image_set.append(image)\n    \n    # limiting the number of images \n    if ibond > 20: break\n    if (time.time() - start_time) / 3600 > 8: break\n\nwith open(image_file, 'wb') as fp:\n    pickle.dump(image_set, fp)\n\nwith open(bond_file, 'wb') as fp:\n    pickle.dump(bond_list, fp)\n        \ntime.time() - start_time","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":false},"cell_type":"markdown","source":"Below, we see how the 5 images looks like for a sample of bonds: "},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n%matplotlib inline\n    \nfor i, num_image in enumerate([1, 5, 7, 8]):\n    \n    image = image_set[num_image]\n    \n    fig = plt.figure(figsize=(16, 3))\n\n    clip_val = 255\n\n    ax1 = fig.add_subplot(151)\n    ax1.imshow(image[0], cmap=plt.cm.BrBG, interpolation='nearest', origin='lower')\n    plt.imshow(np.clip(image[0], 0, clip_val) / clip_val)\n    if i == 0: ax1.set_title('C')\n\n    ax1 = fig.add_subplot(152)\n    ax1.imshow(image[1], cmap=plt.cm.BrBG, interpolation='nearest', origin='lower')\n    plt.imshow(np.clip(image[1], 0, clip_val)/ clip_val)\n    if i == 0: ax1.set_title('N')\n\n    ax1 = fig.add_subplot(153)\n    ax1.imshow(image[2], cmap=plt.cm.BrBG, interpolation='nearest', origin='lower')\n    plt.imshow(np.clip(image[2], 0, clip_val)/ clip_val)\n    if i == 0: ax1.set_title('O')\n\n    ax1 = fig.add_subplot(154)\n    ax1.imshow(image[3], cmap=plt.cm.BrBG, interpolation='nearest', origin='lower')\n    plt.imshow(np.clip(image[3], 0, clip_val)/ clip_val)\n    if i == 0: ax1.set_title('H')\n    \n    ax1 = fig.add_subplot(155)\n    ax1.imshow(image[3], cmap=plt.cm.BrBG, interpolation='nearest', origin='lower')\n    plt.imshow(np.clip(image[4], 0, clip_val)/ clip_val)\n    if i == 0: ax1.set_title('F')\n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The next steps then consists in introducing these images in a CNN in order to perform the regression task. An example for the 1JHN symetry is given in [this kernel](https://www.kaggle.com/fabiendaniel/cnn-on-1jhn-2d-images)"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}