{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook generates the bond-topology of the structures and then calculates the shortest path beween the interaction partners. If the shortest path lenght is 2 bonds i.e. there is only one atom between the interaction partners, the angle is calculated. If there are 3 bonds in between, the dihedral-angle is calculated. The angle improves the model for 2J models, the dihedral-angle improves the 3J models."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"%pylab inline\n\nimport numpy as np\nimport os\nimport pandas as pd\nimport seaborn as sns\nfrom glob import glob\nfrom multiprocessing import Process, Pool, Manager, cpu_count\nfrom ipywidgets import IntProgress as Progress, Layout\nimport time\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\npd.options.display.max_columns = 100\n\nimport networkx as nx\nimport numpy as np","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"'''https://chem.libretexts.org/Ancillary_Materials/Reference/Reference_Tables/Atomic_and_Molecular_Properties/A3%3A_Covalent_Radii'''\nbond_lengths = {'H': np.array([0.32, np.NaN, np.NaN]),\n                'C': np.array([.75, .67, .60]),\n                'O': np.array([.63, .57, .54]),\n                'N': np.array([.71, .60, .54]),\n                'F': np.array([.64, .59, .53])}\nbond_lengths","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's take only the first 100k rows for demonstration\ntrain = pd.read_csv('../input/train.csv', nrows=100000)\ntest = pd.read_csv('../input/test.csv', nrows=100000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')\nstructures['possible_radii'] = [bond_lengths[a] for a in structures.atom]\ndisplay(structures.head())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_graph(structure, cutoff=0.2):\n    structure_x = get_structure_x(structure)\n    structure_x = structure_x[~((structure_x.atom_x == 'H') & \n                                (structure_x.atom_y == 'H'))]\n    g = nx.Graph()\n    positions = structure[['x', 'y', 'z']].values\n    g.add_nodes_from(structure['atom_index'], pos=positions)\n    \n    cols = ['atom_index_x', 'atom_index_y', 'distance', 'cutoff']\n    X = [ tuple(i) for i in structure_x[cols].values ]\n    \n    for i_a, i_b, d, c in X:\n        if c < cutoff:\n            g.add_edge(i_a, i_b, length=d)\n        elif nx.number_connected_components(g) == 1:\n            break\n        else:\n            g.add_edge(i_a, i_b, length=d)\n    return g\n\ndef get_structure(molecule_name):\n    return structures[structures['molecule_name'] == molecule_name]\n\ndef get_structure_x(structure):\n    positions = structure[['x', 'y', 'z']].values\n    structure_X = pd.merge(structure, structure, how='outer', on=['molecule_name'])\n    structure_X['distance'] = np.linalg.norm(structure_X[['x_x', 'y_x', 'z_x']].values - structure_X[['x_y', 'y_y', 'z_y']].values, axis=1)\n    structure_X = structure_X[(structure_X.atom_index_x > structure_X.atom_index_y)]\n    structure_X['cutoff'] = (structure_X.distance - (structure_X.possible_radii_x + structure_X.possible_radii_y)).apply(abs).apply(min)\n    structure_X.sort_values('cutoff', inplace=True)\n    return structure_X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_graph(g):\n    #labels = {i[0]: f'{i[0]} {i[1]}' for i in structure[['atom_index', 'atom']].values}\n    nx.draw_spring(g, with_labels=True)\n    show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_relations(structure, g, atom_index_0, atom_index_1):\n\n    shortest_path = nx.shortest_path(g, atom_index_0, atom_index_1)\n    shortest_path_atoms = ''.join([structure[structure.atom_index == i].atom.values[0] for i in shortest_path[1:-1]])\n    shortest_path_n_bonds = len(shortest_path)-1\n    \n    cos = None\n    dihe = None\n    \n    if shortest_path_n_bonds == 2:\n        x0 = structure[structure.atom_index == shortest_path[0]][['x', 'y', 'z']].values[0]\n        x1 = structure[structure.atom_index == shortest_path[1]][['x', 'y', 'z']].values[0]\n        x2 = structure[structure.atom_index == shortest_path[2]][['x', 'y', 'z']].values[0]\n        cos = cosinus(x0, x1, x2)\n        \n    if shortest_path_n_bonds == 3:\n        x0 = structure[structure.atom_index == shortest_path[0]][['x', 'y', 'z']].values[0]\n        x1 = structure[structure.atom_index == shortest_path[1]][['x', 'y', 'z']].values[0]\n        x2 = structure[structure.atom_index == shortest_path[2]][['x', 'y', 'z']].values[0]\n        x3 = structure[structure.atom_index == shortest_path[3]][['x', 'y', 'z']].values[0]\n        dihe = dihedral(x0, x1, x2, x3)\n        \n    results = {\n        'molecule_name': structure.molecule_name.values[0],\n        'atom_index_0': atom_index_0,\n        'atom_index_1': atom_index_1,\n        'shortest_path_atoms': shortest_path_atoms,\n        'shortest_path_n_bonds': shortest_path_n_bonds,\n        'cosinus': cos,\n        'dihedral': dihe\n               }\n    \n    return pd.DataFrame(results, index=[0])\n\n\ndef cosinus(x0, x1, x2):\n    e0 = (x0-x1)\n    e1 = (x2-x1)\n    e0 = (e0 / np.linalg.norm(e0))\n    e1 = (e1 / np.linalg.norm(e1))\n    cosinus = np.dot(e0, e1)\n    return np.round(cosinus, 5)\n\n\ndef dihedral(x0, x1, x2, x3):\n\n    b0 = -1.0 * (x1 - x0)\n    b1 = x2 - x1\n    b2 = x3 - x2\n\n    b0xb1 = np.cross(b0, b1)\n    b1xb2 = np.cross(b2, b1)\n\n    b0xb1_x_b1xb2 = np.cross(b0xb1, b1xb2)\n\n    y = np.dot(b0xb1_x_b1xb2, b1)*(1.0/np.linalg.norm(b1))\n    x = np.dot(b0xb1, b1xb2)\n    \n    grad = np.arctan2(y, x)\n    return grad","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structure = get_structure('dsgdb9nsd_003830')\ng = get_graph(structure)\ndisplay(get_relations(structure, g, 3, 6))\nplot_graph(g)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def process(args):\n    molecule_name = args['molecule_name']\n    structure = get_structure(molecule_name)\n    g = get_graph(structure)\n    ndxs = zip(args['atom_index_0s'], args['atom_index_1s'])\n    return pd.concat([get_relations(structure, g, i, j) for i,j in ndxs])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.concat([train, test])[['molecule_name', 'atom_index_0', 'atom_index_1']]\ndel train, test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import multiprocessing\nmultiprocessing.cpu_count()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pool = Pool(processes=4)\nm = Manager()\nq = m.Queue()\n\nargs = []\nstart = time.time()\n\nfor molecule_name, df in df.groupby('molecule_name'):\n    args.append({'molecule_name': molecule_name,\n                 'atom_index_0s': df.atom_index_0.values,\n                 'atom_index_1s': df.atom_index_1.values})\n\nresults = pool.map_async(process, args)\n\npool.close()\npool.join()\nend = time.time()\n\nprint(end - start)\nprint('Run time:', np.round((end - start) / 60 / 60, 2), 'h')\n\nresult = pd.concat(results.get())\nresult.to_csv('angles.csv', index=False)\n\nprint(len(result))\ndisplay(result)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/train.csv')\ncompare = pd.merge(train, result, on=['molecule_name', 'atom_index_0', 'atom_index_1'])\n\nfor t, df_tmp in compare.groupby('type'):\n    figure(figsize=(12,4))\n    subplot(1,3,1)\n    plt.scatter(df_tmp.scalar_coupling_constant, df_tmp.cosinus)\n    plt.title('cos')\n    subplot(1,3,2)\n    plt.scatter(df_tmp.scalar_coupling_constant, df_tmp.dihedral)\n    plt.xlabel('scalar_coupling_constant')\n    plt.title('dihedral')\n    ax = subplot(1,3,3)\n    df_tmp.shortest_path_n_bonds.value_counts().sort_index().plot.bar(ax=ax)\n    plt.yscale('log')\n    plt.title('shortest path n_bonds')\n    plt.suptitle(t, x=0.1, size=18)\n    plt.tight_layout()\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare['n_type'] = compare.type.apply(lambda x: x[0])\nfailed = compare[compare.n_type.astype(int) != compare.shortest_path_n_bonds.astype(int)]\nfailed","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"list(failed.molecule_name.drop_duplicates())","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}