{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport math\nimport h5py\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\nfrom tensorflow.python.framework import ops\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom pandas.plotting import scatter_matrix\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"#Load data\ntrain = pd.read_csv(\"../input/train.csv\")\ntest = pd.read_csv(\"../input/test.csv\")\n\n#Load sub data\nsub = pd.read_csv(\"../input/sample_submission.csv\")\nstructures = pd.read_csv(\"../input/structures.csv\")\nscalar_coupling_contributions = pd.read_csv(\"../input/scalar_coupling_contributions.csv\")\nmagnetic_shielding_tesors = pd.read_csv(\"../input/magnetic_shielding_tensors.csv\")\ndipole_moments = pd.read_csv(\"../input/dipole_moments.csv\")\npotential_energy = pd.read_csv(\"../input/potential_energy.csv\")\nmulliken_charges = pd.read_csv(\"../input/mulliken_charges.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train.shape)\nprint(test.shape)\nprint(sub.shape)\nprint(structures.shape)\nprint(scalar_coupling_contributions.shape)\nprint(magnetic_shielding_tesors.shape)\nprint(dipole_moments.shape)\nprint(potential_energy.shape)\nprint(mulliken_charges.shape)\n\n#check number of unique data \nprint(train[\"molecule_name\"].nunique())\nprint(test[\"molecule_name\"].nunique())\nprint(structures[\"atom\"].nunique())\nprint(train[\"type\"].nunique())","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":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","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Merge Datasets"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_structure = map_atom_info(train, 0)\ntrain_structure = map_atom_info(train_structure, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_structure = map_atom_info(test, 0)\ntest_structure = map_atom_info(test_structure, 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Computer Distance of atoms\nThe distance of atoms is sqrt((x_0 - x_1)^2 + (y_0 - y_1)^2 + (z_0 - z_1)^2), so I use numpy.linalg.norm."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_p_0 = train_structure[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train_structure[['x_1', 'y_1', 'z_1']].values\ntrain_structure['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_p_0 = test_structure[['x_0', 'y_0', 'z_0']].values\ntest_p_1 = test_structure[['x_1', 'y_1', 'z_1']].values\ntest_structure['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Academically, atomic number become large, scalar coupling constant also become large.\nSo, I want to add atomic number information in dataset."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_structure[\"atomic_number_0\"]=train_structure[\"atom_0\"].replace(\"H\", 1).replace(\"C\", 6).replace(\"N\", 7).replace(\"O\", 8).replace(\"F\", 9)\ntrain_structure[\"atomic_number_1\"]=train_structure[\"atom_1\"].replace(\"H\", 1).replace(\"C\", 6).replace(\"N\", 7).replace(\"O\", 8).replace(\"F\", 9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_structure[\"atomic_number_0\"]=test_structure[\"atom_0\"].replace(\"H\", 1).replace(\"C\", 6).replace(\"N\", 7).replace(\"O\", 8).replace(\"F\", 9)\ntest_structure[\"atomic_number_1\"]=test_structure[\"atom_1\"].replace(\"H\", 1).replace(\"C\", 6).replace(\"N\", 7).replace(\"O\", 8).replace(\"F\", 9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In addition, maybe coupling number have effect for scalar coupling constant, let's add that."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_structure[\"coupling_number_0\"]=train_structure[\"atom_0\"].replace(\"H\", 1).replace(\"C\", 4).replace(\"N\", 3).replace(\"O\", 2).replace(\"F\", 1)\ntrain_structure[\"coupling_number_1\"]=train_structure[\"atom_1\"].replace(\"H\", 1).replace(\"C\", 4).replace(\"N\", 3).replace(\"O\", 2).replace(\"F\", 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_structure[\"coupling_number_0\"]=test_structure[\"atom_0\"].replace(\"H\", 1).replace(\"C\", 4).replace(\"N\", 3).replace(\"O\", 2).replace(\"F\", 1)\ntest_structure[\"coupling_number_1\"]=test_structure[\"atom_1\"].replace(\"H\", 1).replace(\"C\", 4).replace(\"N\", 3).replace(\"O\", 2).replace(\"F\", 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_structure.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## New features using RDKIT "},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install --quiet cupy-cuda100\n!pip install --quiet chainer-chemistry\n!pip install --quiet chaineripy\n!conda install -y --quiet -c rdkit rdkit","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check correctly installed, and modules can be imported.\nimport chainer\nimport chainer_chemistry\nimport chaineripy\nimport cupy\nimport rdkit\n\nprint('chainer version: ', chainer.__version__)\nprint('cupy version: ', cupy.__version__)\nprint('chainer-chemistry version: ', chainer_chemistry.__version__)\nprint('rdkit version:', rdkit.__version__)\n\nfrom contextlib import contextmanager\nimport gc\nfrom pathlib import Path\nfrom time import time, perf_counter\n\nimport seaborn as sns\nfrom rdkit import Chem\nfrom chainer_chemistry.datasets.numpy_tuple_dataset import NumpyTupleDataset","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def timer(name):\n    t0 = perf_counter()\n    yield\n    t1 = perf_counter()\n    print('[{}] done in {:.3f} s'.format(name, t1-t0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\nCopied from\nhttps://github.com/jensengroup/xyz2mol/blob/master/xyz2mol.py\n\nModified `chiral_stereo_check` method for this task's purpose.\n\"\"\"\n##\n# Written by Jan H. Jensen based on this paper Yeonjoon Kim and Woo Youn Kim\n# \"Universal Structure Conversion Method for Organic Molecules: From Atomic Connectivity\n# to Three-Dimensional Geometry\" Bull. Korean Chem. Soc. 2015, Vol. 36, 1769-1777 DOI: 10.1002/bkcs.10334\n#\nfrom rdkit.Chem import AllChem\nimport itertools\nfrom rdkit.Chem import rdmolops\nfrom collections import defaultdict\nimport copy\nimport networkx as nx  # uncomment if you don't want to use \"quick\"/install networkx\n\nglobal __ATOM_LIST__\n__ATOM_LIST__ = [x.strip() for x in ['h ', 'he', \\\n                                     'li', 'be', 'b ', 'c ', 'n ', 'o ', 'f ', 'ne', \\\n                                     'na', 'mg', 'al', 'si', 'p ', 's ', 'cl', 'ar', \\\n                                     'k ', 'ca', 'sc', 'ti', 'v ', 'cr', 'mn', 'fe', 'co', 'ni', 'cu', \\\n                                     'zn', 'ga', 'ge', 'as', 'se', 'br', 'kr', \\\n                                     'rb', 'sr', 'y ', 'zr', 'nb', 'mo', 'tc', 'ru', 'rh', 'pd', 'ag', \\\n                                     'cd', 'in', 'sn', 'sb', 'te', 'i ', 'xe', \\\n                                     'cs', 'ba', 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu', 'gd', 'tb', 'dy', \\\n                                     'ho', 'er', 'tm', 'yb', 'lu', 'hf', 'ta', 'w ', 're', 'os', 'ir', 'pt', \\\n                                     'au', 'hg', 'tl', 'pb', 'bi', 'po', 'at', 'rn', \\\n                                     'fr', 'ra', 'ac', 'th', 'pa', 'u ', 'np', 'pu']]\n\n\ndef get_atom(atom):\n    global __ATOM_LIST__\n    atom = atom.lower()\n    return __ATOM_LIST__.index(atom) + 1\n\n\ndef getUA(maxValence_list, valence_list):\n    UA = []\n    DU = []\n    for i, (maxValence, valence) in enumerate(zip(maxValence_list, valence_list)):\n        if maxValence - valence > 0:\n            UA.append(i)\n            DU.append(maxValence - valence)\n    return UA, DU\n\n\ndef get_BO(AC, UA, DU, valences, UA_pairs, quick):\n    BO = AC.copy()\n    DU_save = []\n\n    while DU_save != DU:\n        for i, j in UA_pairs:\n            BO[i, j] += 1\n            BO[j, i] += 1\n\n        BO_valence = list(BO.sum(axis=1))\n        DU_save = copy.copy(DU)\n        UA, DU = getUA(valences, BO_valence)\n        UA_pairs = get_UA_pairs(UA, AC, quick)[0]\n\n    return BO\n\n\ndef valences_not_too_large(BO, valences):\n    number_of_bonds_list = BO.sum(axis=1)\n    for valence, number_of_bonds in zip(valences, number_of_bonds_list):\n        if number_of_bonds > valence:\n            return False\n\n    return True\n\n\ndef BO_is_OK(BO, AC, charge, DU, atomic_valence_electrons, atomicNumList, charged_fragments):\n    Q = 0  # total charge\n    q_list = []\n    if charged_fragments:\n        BO_valences = list(BO.sum(axis=1))\n        for i, atom in enumerate(atomicNumList):\n            q = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])\n            Q += q\n            if atom == 6:\n                number_of_single_bonds_to_C = list(BO[i, :]).count(1)\n                if number_of_single_bonds_to_C == 2 and BO_valences[i] == 2:\n                    Q += 1\n                    q = 2\n                if number_of_single_bonds_to_C == 3 and Q + 1 < charge:\n                    Q += 2\n                    q = 1\n\n            if q != 0:\n                q_list.append(q)\n\n    if (BO - AC).sum() == sum(DU) and charge == Q and len(q_list) <= abs(charge):\n        return True\n    else:\n        return False\n\n\ndef get_atomic_charge(atom, atomic_valence_electrons, BO_valence):\n    if atom == 1:\n        charge = 1 - BO_valence\n    elif atom == 5:\n        charge = 3 - BO_valence\n    elif atom == 15 and BO_valence == 5:\n        charge = 0\n    elif atom == 16 and BO_valence == 6:\n        charge = 0\n    else:\n        charge = atomic_valence_electrons - 8 + BO_valence\n\n    return charge\n\n\ndef clean_charges(mol):\n    # this hack should not be needed any more but is kept just in case\n    #\n\n    rxn_smarts = ['[N+:1]=[*:2]-[C-:3]>>[N+0:1]-[*:2]=[C-0:3]',\n                  '[N+:1]=[*:2]-[O-:3]>>[N+0:1]-[*:2]=[O-0:3]',\n                  '[N+:1]=[*:2]-[*:3]=[*:4]-[O-:5]>>[N+0:1]-[*:2]=[*:3]-[*:4]=[O-0:5]',\n                  '[#8:1]=[#6:2]([!-:6])[*:3]=[*:4][#6-:5]>>[*-:1][*:2]([*:6])=[*:3][*:4]=[*+0:5]',\n                  '[O:1]=[c:2][c-:3]>>[*-:1][*:2][*+0:3]',\n                  '[O:1]=[C:2][C-:3]>>[*-:1][*:2]=[*+0:3]']\n\n    fragments = Chem.GetMolFrags(mol, asMols=True, sanitizeFrags=False)\n\n    for i, fragment in enumerate(fragments):\n        for smarts in rxn_smarts:\n            patt = Chem.MolFromSmarts(smarts.split(\">>\")[0])\n            while fragment.HasSubstructMatch(patt):\n                rxn = AllChem.ReactionFromSmarts(smarts)\n                ps = rxn.RunReactants((fragment,))\n                fragment = ps[0][0]\n        if i == 0:\n            mol = fragment\n        else:\n            mol = Chem.CombineMols(mol, fragment)\n\n    return mol\n\n\ndef BO2mol(mol, BO_matrix, atomicNumList, atomic_valence_electrons, mol_charge, charged_fragments):\n    # based on code written by Paolo Toscani\n\n    l = len(BO_matrix)\n    l2 = len(atomicNumList)\n    BO_valences = list(BO_matrix.sum(axis=1))\n\n    if (l != l2):\n        raise RuntimeError('sizes of adjMat ({0:d}) and atomicNumList '\n                           '{1:d} differ'.format(l, l2))\n\n    rwMol = Chem.RWMol(mol)\n\n    bondTypeDict = {\n        1: Chem.BondType.SINGLE,\n        2: Chem.BondType.DOUBLE,\n        3: Chem.BondType.TRIPLE\n    }\n\n    for i in range(l):\n        for j in range(i + 1, l):\n            bo = int(round(BO_matrix[i, j]))\n            if (bo == 0):\n                continue\n            bt = bondTypeDict.get(bo, Chem.BondType.SINGLE)\n            rwMol.AddBond(i, j, bt)\n    mol = rwMol.GetMol()\n\n    if charged_fragments:\n        mol = set_atomic_charges(mol, atomicNumList, atomic_valence_electrons, BO_valences, BO_matrix, mol_charge)\n    else:\n        mol = set_atomic_radicals(mol, atomicNumList, atomic_valence_electrons, BO_valences)\n\n    return mol\n\n\ndef set_atomic_charges(mol, atomicNumList, atomic_valence_electrons, BO_valences, BO_matrix, mol_charge):\n    q = 0\n    for i, atom in enumerate(atomicNumList):\n        a = mol.GetAtomWithIdx(i)\n        charge = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])\n        q += charge\n        if atom == 6:\n            number_of_single_bonds_to_C = list(BO_matrix[i, :]).count(1)\n            if number_of_single_bonds_to_C == 2 and BO_valences[i] == 2:\n                q += 1\n                charge = 0\n            if number_of_single_bonds_to_C == 3 and q + 1 < mol_charge:\n                q += 2\n                charge = 1\n\n        if (abs(charge) > 0):\n            a.SetFormalCharge(int(charge))\n\n    # shouldn't be needed anymore bit is kept just in case\n    # mol = clean_charges(mol)\n\n    return mol\n\n\ndef set_atomic_radicals(mol, atomicNumList, atomic_valence_electrons, BO_valences):\n    # The number of radical electrons = absolute atomic charge\n    for i, atom in enumerate(atomicNumList):\n        a = mol.GetAtomWithIdx(i)\n        charge = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])\n\n        if (abs(charge) > 0):\n            a.SetNumRadicalElectrons(abs(int(charge)))\n\n    return mol\n\n\ndef get_bonds(UA, AC):\n    bonds = []\n\n    for k, i in enumerate(UA):\n        for j in UA[k + 1:]:\n            if AC[i, j] == 1:\n                bonds.append(tuple(sorted([i, j])))\n\n    return bonds\n\n\ndef get_UA_pairs(UA, AC, quick):\n    bonds = get_bonds(UA, AC)\n    if len(bonds) == 0:\n        return [()]\n\n    if quick:\n        G = nx.Graph()\n        G.add_edges_from(bonds)\n        UA_pairs = [list(nx.max_weight_matching(G))]\n        return UA_pairs\n\n    max_atoms_in_combo = 0\n    UA_pairs = [()]\n    for combo in list(itertools.combinations(bonds, int(len(UA) / 2))):\n        flat_list = [item for sublist in combo for item in sublist]\n        atoms_in_combo = len(set(flat_list))\n        if atoms_in_combo > max_atoms_in_combo:\n            max_atoms_in_combo = atoms_in_combo\n            UA_pairs = [combo]\n        #           if quick and max_atoms_in_combo == 2*int(len(UA)/2):\n        #               return UA_pairs\n        elif atoms_in_combo == max_atoms_in_combo:\n            UA_pairs.append(combo)\n\n    return UA_pairs\n\n\ndef AC2BO(AC, atomicNumList, charge, charged_fragments, quick):\n    # TODO\n    atomic_valence = defaultdict(list)\n    atomic_valence[1] = [1]\n    atomic_valence[6] = [4]\n    atomic_valence[7] = [4, 3]\n    atomic_valence[8] = [2, 1]\n    atomic_valence[9] = [1]\n    atomic_valence[14] = [4]\n    atomic_valence[15] = [5, 4, 3]\n    atomic_valence[16] = [6, 4, 2]\n    atomic_valence[17] = [1]\n    atomic_valence[32] = [4]\n    atomic_valence[35] = [1]\n    atomic_valence[53] = [1]\n\n    atomic_valence_electrons = {}\n    atomic_valence_electrons[1] = 1\n    atomic_valence_electrons[6] = 4\n    atomic_valence_electrons[7] = 5\n    atomic_valence_electrons[8] = 6\n    atomic_valence_electrons[9] = 7\n    atomic_valence_electrons[14] = 4\n    atomic_valence_electrons[15] = 5\n    atomic_valence_electrons[16] = 6\n    atomic_valence_electrons[17] = 7\n    atomic_valence_electrons[32] = 4\n    atomic_valence_electrons[35] = 7\n    atomic_valence_electrons[53] = 7\n\n    # make a list of valences, e.g. for CO: [[4],[2,1]]\n    valences_list_of_lists = []\n    for atomicNum in atomicNumList:\n        valences_list_of_lists.append(atomic_valence[atomicNum])\n\n    # convert [[4],[2,1]] to [[4,2],[4,1]]\n    valences_list = list(itertools.product(*valences_list_of_lists))\n\n    best_BO = AC.copy()\n\n    # implemenation of algorithm shown in Figure 2\n    # UA: unsaturated atoms\n    # DU: degree of unsaturation (u matrix in Figure)\n    # best_BO: Bcurr in Figure\n    #\n\n    for valences in valences_list:\n        AC_valence = list(AC.sum(axis=1))\n        UA, DU_from_AC = getUA(valences, AC_valence)\n\n        if len(UA) == 0 and BO_is_OK(AC, AC, charge, DU_from_AC, atomic_valence_electrons, atomicNumList,\n                                     charged_fragments):\n            return AC, atomic_valence_electrons\n\n        UA_pairs_list = get_UA_pairs(UA, AC, quick)\n        for UA_pairs in UA_pairs_list:\n            BO = get_BO(AC, UA, DU_from_AC, valences, UA_pairs, quick)\n            if BO_is_OK(BO, AC, charge, DU_from_AC, atomic_valence_electrons, atomicNumList, charged_fragments):\n                return BO, atomic_valence_electrons\n\n            elif BO.sum() >= best_BO.sum() and valences_not_too_large(BO, valences):\n                best_BO = BO.copy()\n\n    return best_BO, atomic_valence_electrons\n\n\ndef AC2mol(mol, AC, atomicNumList, charge, charged_fragments, quick):\n    # convert AC matrix to bond order (BO) matrix\n    BO, atomic_valence_electrons = AC2BO(AC, atomicNumList, charge, charged_fragments, quick)\n\n    # add BO connectivity and charge info to mol object\n    mol = BO2mol(mol, BO, atomicNumList, atomic_valence_electrons, charge, charged_fragments)\n\n    return mol\n\n\ndef get_proto_mol(atomicNumList):\n    mol = Chem.MolFromSmarts(\"[#\" + str(atomicNumList[0]) + \"]\")\n    rwMol = Chem.RWMol(mol)\n    for i in range(1, len(atomicNumList)):\n        a = Chem.Atom(atomicNumList[i])\n        rwMol.AddAtom(a)\n\n    mol = rwMol.GetMol()\n\n    return mol\n\n\ndef get_atomicNumList(atomic_symbols):\n    atomicNumList = []\n    for symbol in atomic_symbols:\n        atomicNumList.append(get_atom(symbol))\n    return atomicNumList\n\n\ndef read_xyz_file(filename):\n    atomic_symbols = []\n    xyz_coordinates = []\n\n    with open(filename, \"r\") as file:\n        for line_number, line in enumerate(file):\n            if line_number == 0:\n                num_atoms = int(line)\n            elif line_number == 1:\n                if \"charge=\" in line:\n                    charge = int(line.split(\"=\")[1])\n                else:\n                    charge = 0\n            else:\n                atomic_symbol, x, y, z = line.split()\n                atomic_symbols.append(atomic_symbol)\n                xyz_coordinates.append([float(x), float(y), float(z)])\n\n    atomicNumList = get_atomicNumList(atomic_symbols)\n\n    return atomicNumList, charge, xyz_coordinates\n\n\ndef xyz2AC(atomicNumList, xyz):\n    import numpy as np\n    mol = get_proto_mol(atomicNumList)\n\n    conf = Chem.Conformer(mol.GetNumAtoms())\n    for i in range(mol.GetNumAtoms()):\n        conf.SetAtomPosition(i, (xyz[i][0], xyz[i][1], xyz[i][2]))\n    mol.AddConformer(conf)\n\n    dMat = Chem.Get3DDistanceMatrix(mol)\n    pt = Chem.GetPeriodicTable()\n\n    num_atoms = len(atomicNumList)\n    AC = np.zeros((num_atoms, num_atoms)).astype(int)\n\n    for i in range(num_atoms):\n        a_i = mol.GetAtomWithIdx(i)\n        Rcov_i = pt.GetRcovalent(a_i.GetAtomicNum()) * 1.30\n        for j in range(i + 1, num_atoms):\n            a_j = mol.GetAtomWithIdx(j)\n            Rcov_j = pt.GetRcovalent(a_j.GetAtomicNum()) * 1.30\n            if dMat[i, j] <= Rcov_i + Rcov_j:\n                AC[i, j] = 1\n                AC[j, i] = 1\n\n    return AC, mol\n\n\ndef chiral_stereo_check(mol):\n    # Chem.SanitizeMol(mol)\n    num_error = Chem.SanitizeMol(mol, Chem.SANITIZE_ALL ^ Chem.SANITIZE_PROPERTIES, catchErrors=True)\n    if num_error != 0:\n        print('error id', num_error)\n\n    Chem.DetectBondStereochemistry(mol, -1)\n    Chem.AssignStereochemistry(mol, flagPossibleStereoCenters=True, force=True)\n    Chem.AssignAtomChiralTagsFromStructure(mol, -1)\n\n    return mol\n\n\ndef xyz2mol(atomicNumList, charge, xyz_coordinates, charged_fragments, quick):\n    # Get atom connectivity (AC) matrix, list of atomic numbers, molecular charge,\n    # and mol object with no connectivity information\n    AC, mol = xyz2AC(atomicNumList, xyz_coordinates)\n\n    # Convert AC to bond order matrix and add connectivity and charge info to mol object\n    new_mol = AC2mol(mol, AC, atomicNumList, charge, charged_fragments, quick)\n\n    # Check for stereocenters and chiral centers\n    new_mol = chiral_stereo_check(new_mol)\n\n    return new_mol","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def mol_from_xyz(filepath, to_canonical=False):\n    charged_fragments = True  # alternatively radicals are made\n\n    # quick is faster for large systems but requires networkx\n    # if you don't want to install networkx set quick=False and\n    # uncomment 'import networkx as nx' at the top of the file\n    quick = True\n\n    atomicNumList, charge, xyz_coordinates = read_xyz_file(filepath)\n    # print('atomicNumList', atomicNumList, 'charge', charge)\n    mol = xyz2mol(atomicNumList, charge, xyz_coordinates, charged_fragments, quick)\n\n    # Canonical hack\n    if to_canonical:\n        smiles = Chem.MolToSmiles(mol, isomericSmiles=True)\n        mol = Chem.MolFromSmiles(smiles)\n    return mol","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# This script is referred from http://rdkit.blogspot.jp/2015/02/new-drawing-code.html\n# and http://cheminformist.itmol.com/TEST/wp-content/uploads/2015/07/rdkit_moldraw2d_2.html\nfrom __future__ import print_function\nfrom rdkit import Chem\nfrom rdkit.Chem.Draw import IPythonConsole\nfrom IPython.display import SVG\n\nfrom rdkit.Chem import rdDepictor\nfrom rdkit.Chem.Draw import rdMolDraw2D\ndef moltosvg(mol,molSize=(450,150),kekulize=True):\n    mc = Chem.Mol(mol.ToBinary())\n    if kekulize:\n        try:\n            Chem.Kekulize(mc)\n        except:\n            mc = Chem.Mol(mol.ToBinary())\n    if not mc.GetNumConformers():\n        rdDepictor.Compute2DCoords(mc)\n    drawer = rdMolDraw2D.MolDraw2DSVG(molSize[0],molSize[1])\n    drawer.DrawMolecule(mc)\n    drawer.FinishDrawing()\n    svg = drawer.GetDrawingText()\n    return svg\n\ndef render_svg(svg):\n    # It seems that the svg renderer used doesn't quite hit the spec.\n    # Here are some fixes to make it work in the notebook, although I think\n    # the underlying issue needs to be resolved at the generation step\n    return SVG(svg.replace('svg:',''))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mol_names = train['molecule_name'].unique()\ntest_mol_names = test['molecule_name'].unique()\nmol_names = np.hstack((train_mol_names, test_mol_names))\nprint(mol_names)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#This for loop need a few minutes calculation\ninput_dir = Path('../input')\nmol_data=[]\nfor mol_id in mol_names:\n    filepath = input_dir/f'structures/{mol_id}.xyz'\n    mol_data.append(mol_from_xyz(filepath, to_canonical=True))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mol =pd.DataFrame({\"molecule_name\":mol_names, \"mol_data\":mol_data})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"smiles=[]\nfor i in range(mol.shape[0]):\n    smiles.append(Chem.MolToSmiles(mol.mol_data[i], isomericSmiles=False))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mol =pd.DataFrame({\"molecule_name\":mol_names, \"mol_data\":mol_data, \"smiles\":smiles})\n#It looks success to convert xyz files\nmol.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Figerprints"},{"metadata":{"trusted":true},"cell_type":"code","source":"from rdkit import rdBase, Chem, DataStructs\nprint(rdBase.rdkitVersion) # 2017.09.1\nfrom rdkit.Avalon import pyAvalonTools\nfrom rdkit.Chem import AllChem\nfrom rdkit.Chem.Fingerprints import FingerprintMols\nfrom rdkit.Chem.AtomPairs import Pairs, Torsions\nfrom rdkit.Chem import MACCSkeys\nfrom rdkit.Chem import Descriptors\nfrom rdkit.ML.Descriptors import MoleculeDescriptors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fps1 = [ MACCSkeys.GenMACCSKeys(mol1).ToBitString() for mol1 in mol.mol_data]\nfps2 = [ list(map(int,list(fps))) for fps in fps1]\nfps3 = np.array(fps2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fps3.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#change data type ndarray to dataframe\nfingerprints = pd.DataFrame(fps3)\nfingerprints.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#merge dataset\nmol = mol.join(fingerprints)\nmol.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Compute descriptors"},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_descriptors(mols, names=None, ipc_avg=False):\n    if names is None:\n        names = [d[0] for d in Descriptors._descList]\n    calc = MoleculeDescriptors.MolecularDescriptorCalculator(names)\n    descs = [calc.CalcDescriptors(mol) for mol in mols]\n    descs = pd.DataFrame(descs, columns=names)\n    if 'Ipc' in names and ipc_avg:\n        descs['Ipc'] = [Descriptors.Ipc(mol, avg=True) for mol in mols]      \n    return descs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"descriptors = calculate_descriptors(mol.mol_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Merge dataset\nmol = mol.join(descriptors)\nmol.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mol.to_csv('mol_dataset.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Evaluate 166 parameters - PCA"},{"metadata":{"trusted":true},"cell_type":"code","source":"def main():\n    mols = mol.mol_data\n\n    fps = [MACCSkeys.GenMACCSKeys(mol1) for mol1 in mols]\n    fpMtx = np.array([fp2arr(fp) for fp in fps])\n\n    plotPCA(fpMtx)\n\n# convert rdkit fingerprint to numpy array\ndef fp2arr(fp):\n    arr = np.zeros((1,))\n    DataStructs.ConvertToNumpyArray(fp, arr)\n    return arr\n\n\n# plot each compound using PCA\ndef plotPCA(fpMtx):\n    from sklearn.decomposition import PCA\n\n    pca = PCA(n_components=2)\n    res = pca.fit_transform(fpMtx)\n    # extract each component\n    pc = res.T\n\n    cum_cr = sum(pca.explained_variance_ratio_)\n    print(\"cumulative contribution ratio=%.2f\" % cum_cr)\n\n    plt.figure()\n    plt.scatter(pc[0], pc[1], marker=\".\")\n    plt.xlabel(\"PC1\")\n    plt.ylabel(\"PC2\")\n    plt.legend()\n    plt.savefig(\"pca_plot.png\")\n\nif __name__ == \"__main__\":\n    main()","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}