{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction\n\nSummary: efficient kernel using the `qml` package to generate eigenvalues of Coulomb matrix/atom-centered Coulomb matrix, and merge eigenvalues into the structures dataframe. \n\nThe Coulomb potential matrix is a symmetric square matrix $\\mathbf{C}_{ij} = (C_{ij})_{n\\times n}$ where $n$ is the total number of atoms in a molecule:\n$$\nC_{ij} =\n  \\begin{cases}\n     \\dfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \n     \\\\\n     \\dfrac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|}       & \\text{if } i \\neq j\n  \\end{cases}\n$$\nwhich is essentially the core of the [QM7 dataset](http://quantum-machine.org/datasets/) which are using a different set of molecules with the Champs competition.\n\n\nIn the formula above, \n- $i, j$: atom indices within a molecule.\n- $Z$: nuclear charge.\n- $\\mathbf{R}$: the coordinate in the Euclidean space.\n\nThe reason it is called \"potential\" is because its spacial gradient gives the Coulomb force (inverse square law).\n\nOne big challege to extract features from the Coulomb matrix is that its size varies from different molecules. So a simple idea is to use its eigenvalues (and/or eigenvectors).\n\nTo use the final `struct_eig` dataframe, you can `pd.merge` it with the original structure file by `molecule_name` and `atom_index`, or more simply, you can drop duplicated columns first then `pd.concat` the two dataframe along axis 1.\n\n\n### Reference\n* http://www.qmlcode.org/tutorial.html\n* Scirpus's kernel: https://www.kaggle.com/scirpus/coulomb-matrix-representation"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install qml","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import qml","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom numpy.linalg import eig, svd\nfrom sklearn.decomposition import PCA\nfrom tqdm import tqdm\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"folder = '../input/structures/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures = pd.read_csv('../input/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_molecule_names = structures['molecule_name'].unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Coulomb matrix playground\n\nHere I demostrate how the matrix is like using a medium-sized molecule (17 atoms)."},{"metadata":{"trusted":true},"cell_type":"code","source":"filenames = ['dsgdb9nsd_133883.xyz']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for filename in tqdm(filenames):\n    # Create the compound object mol from the file which happens to be methane\n    mol = qml.Compound(xyz=folder+filename)\n    mol.generate_coulomb_matrix(size=mol.natoms, sorting=\"unsorted\")\n    a = mol.representation","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## this is the structure of molecule dsgdb9nsd_133883\nstructures.loc[structures.molecule_name==filename[:-4]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(mol.coordinates, '\\n')\nprint(mol.atomtypes, '\\n')\nprint(mol.nuclear_charges,'\\n')\nprint(mol.name, '\\n')\nprint(mol.natoms, '\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`a` is a matrix represented using linear index, the follow routine converts it back to matrix. Since the matrix it represented is symmetric, we only need $n(n+1)/2 = 153$ entries to reprensent the Coulomb matrix."},{"metadata":{"trusted":true},"cell_type":"code","source":"c_mat1 = qml.representations.vector_to_matrix(a)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# this is 17x17 matrix, and it is symmetric\nprint(a.shape, c_mat1.shape, np.allclose(c_mat1, c_mat1.T))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pca = PCA(n_components = min(7,mol.natoms))\npca.fit(c_mat1) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(pca.mean_, '\\n\\n', pca.singular_values_)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Atom-centered Coulomb matrix\n\nThe atom-centered Coulomb matrix is similar to Coulomb matrix except it adds cut-off based on a center atom the user chooses. Now the entry of the matrix is:\n\n$$\nM_{ij}(k) =\n  \\begin{cases}\n     \\dfrac{1}{2} Z_{i}^{2.4} \\cdot f_{ik}^2 & \\text{if } i = j \\\\[5pt]\n     \\dfrac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} \\cdot f_{ik}f_{jk}f_{ij} \n     & \\text{if } i \\neq j,\n\\end{cases}\n$$\n\nwhere \n$$\nf_{ij} =\n  \\begin{cases}\n     1 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \\Delta r \\\\[5pt]\n     \\dfrac{1}{2} \\left(1 + \\cos\\Big(\\pi \\frac{\\|{\\bf R}_{i} - {\\bf R}_{j} \\|\n        - r + \\Delta r}{\\Delta r} \\Big)\\right)     \n        & \\text{if } r - \\Delta r < \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \\Delta r \n        \\\\[5pt]\n     0 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| > r\n  \\end{cases}\n$$\n\nis used as cut-offs similar to Atom-Centered Symmetric Functions (ACSF)."},{"metadata":{"trusted":true},"cell_type":"code","source":"for filename in tqdm(filenames):    \n    mol = qml.Compound(xyz=folder+filename)\n    mol.generate_atomic_coulomb_matrix(size=mol.natoms, \n                                       sorting=\"distance\", \n                                       central_cutoff=6.0, \n                                       central_decay=3, \n                                       interaction_cutoff=3.0, \n                                       interaction_decay=1)\n    a = mol.representation","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now `a` is a 2D array such that its `i`-th row represents `i`-th atom's atom-centered Coulomb matrix. Let us consider the first atom as follows."},{"metadata":{"trusted":true},"cell_type":"code","source":"c_mat2 = qml.representations.vector_to_matrix(a[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 17x17 matrix, and it is symmetric\nprint(a.shape, c_mat2.shape, np.allclose(c_mat2, c_mat2.T))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pca = PCA(n_components=min(7,mol.natoms))\npca.fit(c_mat2)  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(pca.mean_, '\\n\\n', pca.singular_values_)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Generation of eigenvalue related features\n\nThe Coulomb potential matrix $C$ can be used as a distance measure to build the graph Laplacian, which is quite simple:\n$$L=D-A,$$\nwhere where $D$ is the degree matrix and $A$ is the adjacency matrix of the graph. Intuitively speaking, if $C_{ij}$ is big, the interaction between two atoms is strong.\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"filenames = [i + '.xyz' for i in all_molecule_names]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_laplacian(A, tol=1e-10):\n    '''\n    input: square Coulomb matrix\n    '''\n    A = A + tol\n    L = np.exp(-1/A)\n    L[np.diag_indices_from(L)] = 0\n    G = np.diag((L*(L>tol)).sum(axis=1)) - L\n    \n    return G","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stats = []\nNUM_SINGULAR_VALUES = 5\nNUM_PCA_COMPONENTS = 7\nCUTOFF = 1e-8 # cut-off for zero interaction\nTOL = 1e-8 # tol for zero eigenvalue\n\nfor filename in tqdm(filenames):\n    \n    # Create the compound object mol from the file which happens to be methane\n    mol = qml.Compound(xyz=folder+filename)\n    natoms = mol.natoms\n    mol.generate_atomic_coulomb_matrix(size=natoms, \n                                       sorting='distance', \n                                       central_cutoff=6.0, \n                                       central_decay=3, \n                                       interaction_cutoff=3.0, \n                                       interaction_decay=1)\n    ac_c_mat = mol.representation # atom-centered Coulomb matrices collection\n    for i in range(natoms): # a loop for every atoms in this molecule\n        atomstats = {}\n        atomstats['molecule_name'] = filename[:-4]\n        atomstats['atom_index'] = i\n        \n        a = qml.representations.vector_to_matrix(ac_c_mat[i])\n        _, eigvals, _ = svd(a)\n        atomstats['eigv_min'] = eigvals[np.abs(eigvals)>TOL].min()\n        atomstats['eigv_max'] = eigvals.max()\n        atomstats['eigv_gap'] = atomstats['eigv_max'] - atomstats['eigv_min']\n        \n        L = get_laplacian(a, tol=CUTOFF)\n        _, eigvals, _ = svd(L)\n        \n        atomstats['fiedler_eig'] = eigvals[eigvals>TOL][-1]\n        atomstats['connectedness'] = (eigvals<TOL).sum()\n        \n        pca = PCA(n_components = min(NUM_PCA_COMPONENTS, natoms))\n        pca.fit(a)\n        sv = pca.singular_values_\n        atomstats['sv_min'] = sv[sv>TOL][-1]\n        atomstats['coulomb_mean'] = pca.mean_[0]\n        \n        if natoms < NUM_SINGULAR_VALUES: # if there are less than certain atoms/singular values\n            sv = np.r_[sv, np.zeros(NUM_SINGULAR_VALUES-natoms)] \n            \n        for k in range(NUM_SINGULAR_VALUES):\n            atomstats['sv_'+str(k)] = sv[k]\n        stats.append(atomstats)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"struct_eig = pd.DataFrame(stats)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"struct_eig.tail(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structures.tail(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"struct_eig.corrwith(struct_eig['eigv_min'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(struct_eig['eigv_min']);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(struct_eig['fiedler_eig']);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(struct_eig['coulomb_mean']);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(struct_eig['sv_1']);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"struct_eig.to_csv('struct_eigen.csv',index=False)","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.7.3"}},"nbformat":4,"nbformat_minor":1}