{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Intro\nThis kernel is a reaction to the [\"Predicting Mulliken Charges With ACSF Descriptors\"](https://www.kaggle.com/borisdee/predicting-mulliken-charges-with-acsf-descriptors)-kernel from Boris. He introduced the very neat concept of ACSF descriptors and showed how it can be used to predict the Mulliken Charge. IMO unfortunately, he didn't want to share his calculated ACSF values. I'm creating this kernel to fill the gap.\n\nAfter some digging, I found the very nice \"DScribe\" package, including out-of-the-box ACSF functions. Let's go."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# import some stuff\n\nimport ase as ase\nimport pandas as pd\nimport numpy as np\nimport time, copy\nimport dscribe as ds\nfrom dscribe import descriptors","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will 'manually' calculate acsf for our beloved 'dsgdb9nsd_000001' molecule."},{"metadata":{"trusted":true},"cell_type":"code","source":"structure = pd.read_csv('../input/champs-scalar-coupling/structures.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# using rcut, g2 and g4 params from Boris\n\n#For all ACSF functions: rcut = 10.0\nrcut = 10.0\n\n#G2 - eta/Rs couples:\ng2_params = [[1, 2], [0.1, 2], [0.01, 2],\n           [1, 6], [0.1, 6], [0.01, 6]]\n\n#G4 - eta/ksi/lambda triplets:\ng4_params = [[1, 4,  1], [0.1, 4,  1], [0.01, 4,  1], \n           [1, 4, -1], [0.1, 4, -1], [0.01, 4, -1]]\n\ng3_params = None\ng5_params = None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# subselect structure\ntmp_structure = structure.loc[structure.molecule_name==\"dsgdb9nsd_000001\", : ].copy()\n\n# define acsf calculator\nspecies = tmp_structure.atom.unique() #array(['C', 'H'], dtype=object)\nacsf = ds.descriptors.ACSF(\n        species=species,\n        rcut=rcut,\n        g2_params=g2_params,\n        g3_params=g3_params,\n        g4_params=g4_params,\n        g5_params=g5_params,\n    )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create ase.Atoms from structure data\nmolecule_atoms = tmp_structure.loc[:, 'atom']\nmolecule_positions = tmp_structure.loc[:, ['x','y','z']]\n\nmolecule_system = ase.atoms.Atoms(symbols=molecule_atoms, positions=molecule_positions)\n\nprint(molecule_system)\nprint(molecule_system.get_atomic_numbers())\nprint(molecule_system.get_positions())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ok ready to create acsf features\nacsf_features = acsf.create(molecule_system, n_jobs=1) # structure of return is [[#acsf features] for each position in molecule_system]\nacsf_features[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One issue I have with this result is that features don't have any describing labels. Of course, I can't estimate what a g4 with params [x,y,z] would be, but this output alone doesn't describe anything. Unfortunately, I did not find any label_generator in their [git](https://github.com/SINGROUP/dscribe/blob/development/dscribe), so I tried to shadow their implementation and create my own set of labels."},{"metadata":{},"cell_type":"markdown","source":"# Connecting the dots\nNow that you know how to calculate acsf features yourself, we should put the pieces together and create acsf features for all molecules"},{"metadata":{"trusted":true},"cell_type":"code","source":"# some functions to help out\ndef create_feature_labels(species,rcut,\n                          g2_params=None,\n                          g3_params=None,\n                          g4_params=None,\n                          g5_params=None,\n                          transform_to_symbols=True):\n\n    #sub function to transform from atom numbers to chemical symbols\n    def get_atom_id(atom_nr, tranform_to_symbols):\n\n        if transform_to_symbols == True:\n            atom_id = nr_to_symbol[atom_nr]\n        else:\n            atom_id = atom_nr\n\n        return atom_id\n\n    feature_label = []\n\n    g_params={\n        'g1': [rcut],\n        'g2': g2_params,\n        'g3': g3_params,\n        'g4': g4_params,\n        'g5': g5_params\n    }\n\n\n    # create_atom_numbers -> symbol dict\n    tmp_system = ase.Atoms(species, [[0,0,0]]*len(species))\n\n    nr_to_symbol = {number:symbol for symbol, number in\n                    zip(tmp_system.get_chemical_symbols(),tmp_system.get_atomic_numbers())\n                    }\n\n\n    atomic_numbers = sorted(tmp_system.get_atomic_numbers())\n\n    for atom_nr in atomic_numbers:\n\n        atom_id = get_atom_id(atom_nr, transform_to_symbols)\n\n        for g in [\"g1\", \"g2\",\"g3\"]:\n\n            params = g_params[g]\n\n            if params is not None:\n\n                for para in params:\n\n                    feature_label.append(f'feat_acsf_{g}_{atom_id}_{para}')\n\n    for atom_nr in atomic_numbers:\n\n        atom_id = get_atom_id(atom_nr, transform_to_symbols)\n\n        for i in range(0, atom_nr+1):\n\n            if i in atomic_numbers:\n\n                atom_id_2 = get_atom_id(i, transform_to_symbols)\n\n                for g in [\"g4\",\"g5\"]:\n\n                    params = g_params[g]\n\n                    if params is not None:\n\n                        for para in params:\n\n                            feature_label.append(f'feat_acsf_{g}_{atom_id}_{atom_id_2}_{para}')\n\n    return feature_label\n\ndef calculate_symmetric_functions(df_structure, rcut, g2_params=None,\n                                  g3_params=None,\n                                  g4_params=None,\n                                  g5_params=None):\n\n    species = df_structure.atom.unique()\n\n    acsf = ds.descriptors.ACSF(\n        species=species,\n        rcut=rcut,\n        g2_params=g2_params,\n        g3_params=g3_params,\n        g4_params=g4_params,\n        g5_params=g5_params,\n    )\n\n    structure_molecules = df_structure.molecule_name.unique()\n\n    acsf_feature_labels = create_feature_labels(species=species,\n                                                rcut=rcut,\n                                                g2_params=g2_params,\n                                                g3_params=g3_params,\n                                                g4_params=g4_params,\n                                                g5_params=g5_params,\n                                                )\n\n    df_structure= df_structure.reindex(columns = df_structure.columns.tolist() + acsf_feature_labels)\n\n    df_structure = df_structure.sort_values(['molecule_name','atom_index'])\n\n    acsf_structure_chunks = calculate_acsf_in_chunks(structure_molecules, df_structure, acsf, acsf_feature_labels)\n\n    acsf_structure = pd.DataFrame().append(acsf_structure_chunks)\n\n    return acsf_structure\n\ndef calculate_acsf_in_chunks(structure_molecules, df_structure, acsf, acsf_feature_labels, step_size=2000):\n\n    mol_counter = 0\n    max_counter = len(structure_molecules)\n    all_chunks = []\n    tic = time.time()\n    while mol_counter*step_size < max_counter:\n\n        tmp_molecules = structure_molecules[mol_counter*step_size:(mol_counter+1)*step_size]\n\n        tmp_structure = df_structure.loc[df_structure.molecule_name.isin(tmp_molecules),:].copy()\n\n        tmp_results = calculate_acsf_multiple_molecules(tmp_molecules, tmp_structure, acsf, acsf_feature_labels)\n\n        all_chunks.append(tmp_results.copy())\n\n        print((mol_counter+1)*step_size, time.time()-tic)\n\n        mol_counter += 1\n\n    return all_chunks\n\n\ndef calculate_acsf_multiple_molecules(molecule_names, df_structure, acsf, acsf_feature_labels):\n\n    #acsf_feature_labels = [f'feat_acsf_{nr}' for nr in range(0, acsf.get_number_of_features())]\n    #df_molecules = df_structure.loc[df_structure.molecule_name.isin(molecule_names),:].copy()\n    counter = 0\n    tic = time.time()\n    for molecule_name in molecule_names:\n\n        df_molecule = df_structure.loc[df_structure.molecule_name == molecule_name,:]\n        acsf_values = calculate_acsf_single_molecule(df_molecule, acsf)\n\n\n        df_structure.loc[df_structure.molecule_name==molecule_name, acsf_feature_labels] = copy.copy(acsf_values)\n\n        counter += 1\n\n    #print(counter, time.time() - tic)\n\n    return df_structure\n\ndef calculate_acsf_single_molecule(df_molecule, acsf):\n\n    molecule_atoms = df_molecule.loc[:, 'atom']\n    molecule_positions = df_molecule.loc[:, ['x','y','z']]\n\n    molecule_system = ase.atoms.Atoms(symbols=molecule_atoms, positions=molecule_positions)\n\n    return acsf.create(molecule_system, n_jobs=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate acsf features with Boris parameter\n# this should take ~ 4 hours on kaggle\n# there is some issue with ase using 10 cores by default that I couldn't disable. \n# It should be possible to calculate way faster\n\n# I'm only using our beloved molecule to show how the output would look like. Remove the .loc condition if you want to recalc everything.\nacsf_structure = calculate_symmetric_functions(structure.loc[structure.molecule_name=='dsgdb9nsd_000001',:].copy(), rcut, \n                                                   g2_params=g2_params,\n                                                  g4_params=g4_params)\n\n\nacsf_structure.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check labels\nLet's calculate some g1 and g2s to check the labels"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Boris has a better way of doing this, but I'm trying to keep it simple here\ndef dist(coord_0, coord_1):    \n    return np.sqrt(np.sum((coord_0-coord_1)**2))\n\ndef fc(dist, rcut):\n    return  0.5*(np.cos(np.pi * dist / rcut)+1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# my beloved molecule\ntest_molecule = structure.loc[structure.molecule_name=='dsgdb9nsd_000001',:]\ncoord_c = test_molecule.loc[test_molecule.atom == 'C', ['x','y','z']].values[0]\n\n# G1 in regards to atoms of type H\ng1_H = 0\nfor coord_h in test_molecule.loc[test_molecule.atom == 'H', ['x','y','z']].values:\n    \n    dist_h_c = dist(coord_c, coord_h)\n    \n    if dist_h_c <= rcut:\n        g1_H += fc(dist_h_c, rcut)\n        \nprint(f'g1 value is {g1_H}, using rcut: {rcut}')\n\n\nfor para in g2_params:\n    eta= para[0]\n    rs = para[1]\n    g2_H = 0\n    for coord_h in test_molecule.loc[test_molecule.atom == 'H', ['x','y', 'z']].values:\n\n        dist_h_c = dist(coord_c, coord_h)\n\n        g2_H += np.exp(-eta*(dist_h_c-rs)**2) * fc(dist_h_c, rcut)\n    \n    print(f'g2 value is {g2_H}, using eta: {eta}, rs: {rs}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Compare the values with the labes above - looks good to me","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Compare the results\nNow we calculated a lot of numbers but do they have any meaning?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# load full dataset\n\nstructure_acsf = pd.read_csv('../input/molecules-structure-acsf/structure_with_acsf.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's compare feature number.\n# Interestingly Boris has twice as many features (!). I honestly have little to no clue about acsf but I'm going to trust the dscribe package on this one.\nfeature_columns = [col for col in structure_acsf.columns if col.startswith('feat_acsf')]\nlen_features = len(feature_columns)\nprint(f\"We have {len_features} feautres\")\nprint(f\"Boris announced ~ 250\")\nprint(\"Maybe he is using two sets of rcut (?)\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"structure_acsf.loc[structure_acsf.molecule_name == 'dsgdb9nsd_000001',:]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's check the numbers\n# A quick look revealed the following mappings:\n#BorisFeatNr -> FeatureLabel\n0 -> feat_acsf_g1_H_10.0\n1 -> feat_acsf_g2_H_[0.01, 2]\n\n4 -> feat_acsf_g2_H_[1, 2]\n5 -> feat_acsf_g2_H_[0.1, 2]\n6 -> feat_acsf_g2_H_[0.01, 2]\n\n12 -> feat_acsf_g2_H_[0.01, 6]\n18 -> feat_acsf_g2_C_[0.1, 2]\n25 -> feat_acsf_g2_C_[0.01, 6]\n\n# I couldn't map all of them but I think this is enough indication to give it a shot","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Words of warning and encouragement\nThis is my first kernel + public dataset so please let me know if I broke some rules or should do things differently.\nI tried to check everything (eg. created labels) but they still might be wrong. If you find inconsistencies or bugs please let me know."}],"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}