{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Synopsis\n\n![OpenBabel](http://openbabel.org/babel130.png)\n\nIn this kernel I will explore an alternative way of calculating local `mulliken charges`. \nYou may be interrested only if you wish to incorportate these features in your model.  \n=> The charges calculated with Open Babel are available as an autonomous dataset in [Open Babel Atom Charges](https://www.kaggle.com/asauve/open-babel-atom-charges)\n\nThe [Open Babel](http://openbabel.org/wiki/Main_Page) Package allow to\n* Read .xyz files (the ones provided in the structures directory)\n* Build bonding scheme (groovy baby!)\n* Compute partial charges, with several builtin methods\n\nHumm this last one is really interresting as it provides a quick way of computing atoms local charge with a method based on linear algebra. Hence it is fast and in the range of the allowed options for the competition.\n\nThe most promising method is Electronegativity Equalization Method  (EEM) which is describeed in this paper:  \n[High-quality and universal empirical atomic charges for chemoinformatics applications](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4667495/)\n\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Changelog\n\n* v7 : Include test molecules in a separate dataframe and updated the external dataset accordingly\n* v5 : building of output dataset `ob_charges` with train molecules\n* v4 : first public version\n"},{"metadata":{},"cell_type":"markdown","source":"# Definition of partial charge\n\nThere is no such thing as a space localized electron, they are rather something like a wave all around atomic nuclei. But there are some models for getting an approximation of the average charge at the location of atoms. This is this what is computed in this kernel under the name partial charge. \n\nOpen babel offers several methods to compute these values (you have to dig into the code to get them) which are all provided under the column with the matching name."},{"metadata":{},"cell_type":"markdown","source":"# Load Data"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nfrom scipy.stats.stats import pearsonr\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom tqdm import tqdm_notebook as tqdm\nimport seaborn as sns \nimport matplotlib.pyplot as plt\nsns.set()\nimport os","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load CSV data"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def load_dir_csv(directory, csv_files=None):\n    if csv_files is None:\n        csv_files = sorted( [ f for f in os.listdir(directory) if f.endswith(\".csv\") ])    \n    csv_vars  = [ filename[:-4] for filename in csv_files ]\n    gdict = globals()\n    for filename, var in zip( csv_files, csv_vars ):\n        print(f\"{var:32s} = pd.read_csv({directory}/{filename})\")\n        gdict[var] = pd.read_csv( f\"{directory}/{filename}\" )\n        print(f\"{'nb of cols ':32s} = \" + str(len(gdict[var])))\n        display(gdict[var].head())\n\nload_dir_csv(\"../input/\", \n             [\"train.csv\", \"test.csv\", \"structures.csv\", \"mulliken_charges.csv\"])\n","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# Install OpenBabel\n\nThe installation method was  demonstrated in: [Molecule with OpenBabel](https://www.kaggle.com/jmtest/molecule-with-openbabel)"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"!conda install -y -c openbabel openbabel \nimport openbabel as ob","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Mulliken Charges from train set"},{"metadata":{"trusted":true},"cell_type":"code","source":"# differentiate train and test set\ntrain_molecules = train.molecule_name.unique()\ntest_molecules  = test.molecule_name.unique()\n\n\nmulliken   = []\nmulliken_charges_idx = mulliken_charges.set_index(['molecule_name'])\n# ensure mulliken charges are in same order as for partial charges\nfor molecule_name in train_molecules:\n    mc  = mulliken_charges_idx.loc[molecule_name].sort_index()\n    mulliken.extend(mc.mulliken_charge.values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Compute atom partial charges with Open Babel "},{"metadata":{"trusted":true},"cell_type":"code","source":"##\n## Build molecules from files.xyz\n##\n\nobConversion = ob.OBConversion()\n#def read_ob_molecule(molecule_name, datadir=\"../input/champs-scalar-coupling/structures\"):\ndef read_ob_molecule(molecule_name, datadir=\"../input/structures\"):\n    mol = ob.OBMol()\n    path = f\"{datadir}/{molecule_name}.xyz\"\n    if not obConversion.ReadFile(mol, path):\n        raise FileNotFoundError(f\"Could not read molecule {path}\")\n    return mol\n    \n\nob_methods = [ \"eem\", \"mmff94\", \"gasteiger\", \"qeq\", \"qtpie\", \n               \"eem2015ha\", \"eem2015hm\", \"eem2015hn\", \"eem2015ba\", \"eem2015bm\", \"eem2015bn\" ]\n\nstructures_idx = structures.set_index( [\"molecule_name\"] )\ndef get_charges_df(molecule_names):\n    ob_methods_charges = [ [] for _ in ob_methods]\n    ob_molecule_name = []  # container for output  DF\n    ob_atom_index    = []  # container for output  DF\n    ob_error         = []\n    for molecule_name in molecule_names:\n        # fill data for output DF\n        ms = structures_idx.loc[molecule_name].sort_index()\n        natoms = len(ms)\n        ob_molecule_name.extend( [molecule_name] * natoms )\n        ob_atom_index.extend(    ms.atom_index.values )\n\n        # calculate open babel charge for each method\n        mol = read_ob_molecule(molecule_name)\n        assert( mol.NumAtoms() == natoms ) # consistency\n        error = 0\n        for method, charges in zip(ob_methods, ob_methods_charges):\n            ob_charge_model = ob.OBChargeModel.FindType(method)\n            if not ob_charge_model.ComputeCharges(mol):\n                error = 1\n            charges.extend( ob_charge_model.GetPartialCharges() )\n        ob_error.extend([error] * natoms)\n            \n    ob_charges = pd.DataFrame({\n        'molecule_name' : ob_molecule_name,\n        'atom_index'    : ob_atom_index}\n    )\n    for method, charges in zip(ob_methods, ob_methods_charges):\n        ob_charges[method] = charges\n    ob_charges[\"error\"] = ob_error\n    display(ob_charges.head())\n    return ob_charges","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time train_ob_charges = get_charges_df(train_molecules)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%time test_ob_charges = get_charges_df(test_molecules)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Compare Mulliken charges to Open Babel EEM"},{"metadata":{"trusted":true},"cell_type":"code","source":"# correlation plots\ncorrs = []\nfor method in ob_methods:\n    charges = train_ob_charges[method].values\n    fig = plt.figure()\n    ax = sns.scatterplot(mulliken, charges)\n    corr, pval = pearsonr(mulliken, charges)\n    corrs.append(corr)\n    title = f\"method = {method:10s}  corr = {corr:7.4f}\"\n    print(title)\n    plt.title(title)\n    plt.xlabel(\"Mulliken charge\")\n    plt.ylabel(f\"Open Babel {method}\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12,6))\ndata = pd.DataFrame( {'method':ob_methods, 'corr':[abs(c) for c in corrs]}).sort_values('corr', ascending=False)\nax = sns.barplot(data=data, x='corr', y='method', orient=\"h\", dodge=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Write Output dataframes"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_ob_charges.to_csv(\"train_ob_charges.csv\")\ntest_ob_charges.to_csv(\"test_ob_charges.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusion"},{"metadata":{},"cell_type":"markdown","source":"Open Babel local charge estimation is worth a try in our models!\n\nBest match is method **\"eem\"** and some of its variants!\n\nNote that some outliers in the scatter plots come from different convention for how ionization charge is located inside the molecule. For example COO- can have -1 charge set on one O, or let it be shared by both oxygens. Hence the fact that some models differ from the data in the table `mulliken_charges` does not mean that it is systematically worse for the end purpose to predict coupling.\n"}],"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}