{"cells":[{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"library(data.table)\nlibrary(magrittr)\nlibrary(ggplot2)\nlibrary(gridExtra)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this report I demonstrate **step-by-step** how I merged the provided datasets into the training set using R data.table. <br>\nI took a lot of inspiration from the \"Deep Dive into Atoms and Molecules\" of Victor Basu.\n\nHere is an **overview of the final features** of my training dataset:\n\n* 1. Euclidean distance between xyz coordinates\n* 2. Molecule dipole moment along x-axis\n* 3. Molecule dipole moment along y-axis\n* 4. Molecule dipole moment along z-axis\n* 5. Radius of the atom with index zero\n* 6. Radius of the atom with index one\n* 7. The difference in electronegativity between the two atoms of the chemical bond\n* 8. Mulliken charge of the atom with index zero\n* 9. Mulliken charge of the atom with index one\n* 10. Potential energy of the molecule\n* 11. Fermi contact contribution to the scalar coupling constant\n* 12. Spin dipolar contribution to the scalar coupling constant\n* 13. Paramagnetic spin-orbit contribution to the scalar coupling constant\n* 14. diamagnetic spin-orbit contribution to the scalar coupling constant\n\n# 1. Loading and previewing the data"},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train <- fread(\n  \"../input/train.csv\", \n  col.names = c(\"id\", \"molecule.name\", \"atom.index.zero\", \n                 \"atom.index.one\", \"type\", \"scalar.coupling.constant\"))\n\ndata.test  <- fread(\n  \"../input/test.csv\",\n  col.names = c(\"id\", \"molecule.name\", \"atom.index.zero\", \n                \"atom.index.one\", \"type\"))\n\ndata.structures <- fread(\n  \"../input/structures.csv\",\n  col.names = c(\"molecule.name\", \"atom.index\", \n                \"atom.element\", \"x\", \"y\", \"z\"))\n\ndata.dipole <- fread(\n    \"../input/dipole_moments.csv\",  \n    col.names = c(\"molecule.name\", \"dipole.x\", \"dipole.y\", \"dipole.z\"))\n\ndata.mulliken <- fread(\n  \"../input/mulliken_charges.csv\", \n  col.names = c(\"molecule.name\", \"atom.index\", \n                \"mulliken.charge\"))\n\ndata.potential.energy <- fread(\n  \"../input/potential_energy.csv\",\n  col.names = c(\"molecule.name\", \"potential.energy\"))\n\ndata.scalar.coupling.contributions <- fread(\n  \"../input/scalar_coupling_contributions.csv\",\n  col.names = c(\"molecule.name\", \"atom.index.zero\", \"atom.index.one\",\n                \"type\", \"fermi.contact.contribution\", \n                \"spin.dipolar.contribution\", \"paramagnetic.so.contribution\",\n                \"diamagnetic.so.contribution\"))\n\ndata.magnetic.shielding <- fread(\n  \"../input/magnetic_shielding_tensors.csv\")\nnames(data.magnetic.shielding) <-  c(c(\"molecule.name\", \"atom.index\"), tolower(names(data.magnetic.shielding)[3:11]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Preview** the first 10 rows of the:\n\n* 1. training data\n* 2. test data\n* 3. tructures data"},{"metadata":{"trusted":true},"cell_type":"code","source":"sapply(list(data.train, data.test, data.structures), function(x) head(x, 10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As an example, the first molecule is methane with the chemical formula $CH_4$. <br>\nThere are **eight binding types and five elements** present in the molecules."},{"metadata":{"trusted":true},"cell_type":"code","source":"p1 <- ggplot(data.train, aes(type)) +\n  geom_bar(stat=\"count\", col = \"blue\") +\n  ggtitle(\"Counts per binding type in the training data\") +\n  theme_minimal()\n\np2 <- ggplot(data.structures, aes(atom.element)) +\n  geom_bar(stat=\"count\", col = \"blue\") +\n  ggtitle(\"Counts per atom element in the structures data\") +\n  theme_minimal()\n\ngridExtra::grid.arrange(p1, p2, nrow = 2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2. Feature creation\n\nIn this part the supplied data is used to create features. <br>\nSome interesting features are **graphically explored**.\n\n## 2.1 Euclidean distance\n\nThe **distance between two three-dimensional points** (also known as the euclidean distance) is calculated from their coordinates as follows:\n\n$$ distance = \\sqrt{(x_1-x_2)² + (y_1-y_2)² + (z_1-z_2)²}$$\n\nCalculate and add the Euclidean distance as a variable to the train data."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train[, euclidean.distance := (data.structures[\n  data.train, on = .(molecule.name, atom.index == atom.index.zero), .(x, y, z)] -\n    data.structures[\n      data.train, on = .(molecule.name, atom.index == atom.index.one),.(x, y, z)])^2 \n  %>% rowSums\n  %>% sqrt]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Mean euclidean distance per atom pair and binding type. <br>\nThe **1JHC, 1JHN and 2JHH seem to have a narrow distribution** of the Euclidean distance."},{"metadata":{"trusted":true},"cell_type":"code","source":"ggplot(data.train, aes(euclidean.distance, fill = type)) +\n  geom_density(alpha=0.5) +\n  scale_fill_brewer(palette = \"Dark2\") +\n  ggtitle(\"Density plots of the Euclidean distance per binding type\")  +\n  theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2.2 Dipole moments\n\nMerge dipole x, y, z data into the train data by reference."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train[,  c(\"dipole.x\", \"dipole.y\", \"dipole.z\") := data.dipole[\n  data.train, on = .(molecule.name), .(dipole.x, dipole.y, dipole.z)]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2.3 Atom physical properties\n\nAdding atomic radius, mass and electronegativity data for each atom. <br>\nThe values for the **radius, mass and electronegativity** were taken from Wikipedia."},{"metadata":{"trusted":true},"cell_type":"code","source":"atom.physical.properties <- data.table(\n  element = c(\"C\", \"H\", \"N\", \"O\", \"F\"),\n  mass = c(12.0107, 1.00784, 14.0067, 15.999, 18.998403),\n  radius = c(70, 53, 65, 60, 147),\n  electronegativity =c(2.55, 2.2, 3.04, 3.44, 3.98))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Merge mass and radius into the structures data by reference.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"data.structures[, c(\"mass\", \"radius\", \"electronegativity\") := atom.physical.properties[\n  data.structures, on = .(element == atom.element), .(mass, radius, electronegativity)]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add the atom radii of the two atoms into the train data. <br>\nNote that **atom.index.zero always refers to hydrogen**."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train[, radius.atom.zero := data.structures[\n  data.train, on = .(molecule.name, atom.index == atom.index.zero), .(radius)]]\n\ndata.train[, radius.atom.one := data.structures[\n  data.train, on = .(molecule.name, atom.index == atom.index.one), .(radius)]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Calculate **difference in electronegativty** of the two atoms in the bond. <br>\nThe difference in electronegativity is an indication for the **chemical bond strength**."},{"metadata":{},"cell_type":"markdown","source":"## 2.4 Mulliken charges\n\nMerge Mulliken charges data into the training data."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train[, mulliken.charge.atom.zero := data.mulliken[\n  data.train, on = .(molecule.name, atom.index == atom.index.zero),\n  .(mulliken.charge)]]\n\ndata.train[, mulliken.charge.atom.one := data.mulliken[\n  data.train, on = .(molecule.name, atom.index == atom.index.one),\n  .(mulliken.charge)]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A two dimensional density plot shows the **Mulliken charge versus atom index**."},{"metadata":{"trusted":true},"cell_type":"code","source":"ggplot(data.mulliken, aes(atom.index, mulliken.charge)) +\n  geom_density_2d() +\n  theme_minimal()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2.5 Potential energy\n\nMerge the potential energy into the training data."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train[, potential.energy := data.potential.energy[\n  data.train, on = .(molecule.name), potential.energy]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2.6 Scalar coupling contributions \n\nMerge the scalar contributions variables into the training data."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.train[, c(\"fermi.contact.contribution\", \"spin.dipolar.contribution\",\n               \"paramagnetic.so.contribution\",\"diamagnetic.so.contribution\") :=\n             data.scalar.coupling.contributions[\n               data.train, on = .(molecule.name, atom.index.one, atom.index.zero),\n               .(fermi.contact.contribution, spin.dipolar.contribution,\n                 paramagnetic.so.contribution, diamagnetic.so.contribution)]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Merge the scalar coupling constant from the training data into the scalar coupling contributions data."},{"metadata":{"trusted":true},"cell_type":"code","source":"data.scalar.coupling.contributions[, scalar.coupling.constant := data.train[\n  data.scalar.coupling.contributions, on = .(molecule.name, type, atom.index.one, atom.index.zero), .(scalar.coupling.constant)]]\n\nplot.data <- melt(data.scalar.coupling.contributions, \n                  id.vars = c(\"molecule.name\", \"atom.index.zero\", \"atom.index.one\", \"type\", \"scalar.coupling.constant\"),\n                  value.name = \"contribution.value\",\n                  variable.name = \"contribution.type\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A two dimensional density plot shows the various **contributions to the scalar coupling constant**. <br>\nOnly the **first one million observations** of each contribution type are taken into account."},{"metadata":{"trusted":true},"cell_type":"code","source":"ggplot(plot.data[, .SD[1:1000000], by = contribution.type ], aes(contribution.value, scalar.coupling.constant, col = contribution.type)) +\n  geom_density_2d() + \n  ggtitle(\"Individual contributions to the scalar coupling constant\") +\n  theme_minimal() +\n  scale_color_brewer(palette = \"Dark2\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. Summary of the data\n\nSummary statistics for the 14 features."},{"metadata":{"trusted":true},"cell_type":"code","source":"summary(data.train[, -(1:5)])","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}