Machine Learning for Materials and Chemistry

Aims to be a elective Master module, can be followed by Bachelor students. Developed at University of Kassel.

Slides

Exercises

Meant to both deepen understanding of the lecture content, but also to boost Python skills.

Exercise 1: Getting Started

Check out the Python course (Basics and NumPy until Array Generators) and do any three exercises which do not have solutions in the course and which you did not do before.

Exercise 2: Linear regression and model stability

Linear regression is a common technique used in statistics and machine learning to model the relationship between two variables. However, when using linear regression with a small number of data points, the model may be subject to noise and may not accurately capture the underlying pattern in the data. This can result in overfitting, underfitting, or inaccurate predictions.

As more data points are added, the model becomes more accurate and better able to capture the underlying pattern in the data. This is because the addition of more data points provides more information about the relationship between the variables and helps to reduce the impact of noise in the data.

The Theil-Sen regressor is a robust alternative to linear regression that is designed to overcome the issue of noise in small datasets. It uses a method called the "median slope" to estimate the slope of the line that best fits the data, rather than the mean slope used in linear regression. However, as the number of data points increases, the difference between the Theil-Sen regressor and linear regression tends to become smaller, and linear regression can provide accurate results when enough data is available.

We have used this particular regressor to improve Hammett models for reaction barriers.

Task 2.1: Implementing Linear Regression

Write a Python function randomized_points(n: int) -> tuple[np.ndarray, np.ndarray] that generates n random points (\(x\)) uniformly distributed over the interval (0, 1) and the corresponding \(y = ax + b + \epsilon\), where \(a\) is the slope (\(a=2\)), \(b\) is the intercept (\(b=1\)), and \(\epsilon\) is the error term that is drawn from a normal distribution with mean 0 and standard deviation 0.2. The function should return two arrays xs and ys, where xs contains the \(x\)-coordinates of the n data points and ys contains the corresponding \(y\)-coordinates.
Write a Python function linear_regression(xs: np.ndarray, ys: np.ndarray) -> tuple[float, float] that takes in the arrays xs and ys from the previous step and returns the slope and intercept of a linear regression line that fits the data points.

Task 2.2: Using Theil-Sen Regressor from Scipy

Write a Python function theil_sen(xs: np.ndarray, ys: np.ndarray) -> tuple[float, float] that performs linear regression using Theil-Sen Regressor from Scipy. The function should take in the arrays xs and ys from Task 1 and return the slope and intercept of the regression line.

Task 2.3: Assessing the Models

Write a Python function predict(xs: np.ndarray, slope: float, intercept: float) -> np.ndarray that takes in the array xs and the slope and intercept of a regression line and returns an array of predicted \(y\)-values for the corresponding \(x\)-values in xs.
Write a Python function calculate_prediction_error(n: int) -> Tuple[float, float] that generates n random points in the interval \((0, 1)\) and uses the predict function to estimate the corresponding \(y\)-values for each \(x\)-value using the slope and intercept of the regression line. The function should then calculate the prediction error by subtracting the true noise-free \(y\)-values from the estimated \(y\)-values. Finally, the function should return the mean and standard deviation of the prediction error.

Exercise 3: Turning a molecule into a mathematical object

When dealing with molecules in machine learning contexts, we need mathematical objects describing such molecules. The straightforward use of cartesian coordinates and nuclear charges (which enter the Hamiltonian) turns out to be disadvantageous. In this exercise, we explore two alternatives.

Normally, representations taylored to molecules such as the Coulomb Matrix perform better than one-hot encoding. However, sometimes meaningful categorial representations catch more of the Physics involved, such that one-hot-encoding leaves more advanced representations behind, e.g. in this work on learning reaction barriers.

Task 3.1: Implement the Coulomb Matrix representation

Write a python function coulomb_matrix(positions: np.ndarray, charges: np.ndarray) -> np.ndarray to calculate the coulomb matrix representation from the cartesian coordinates (xs, ys, zs) and the nuclear charges.

\[\begin{split}M_{ij} =\begin{cases}\tfrac{1}{2} Z_{i}^{2.4} & \text{if } i = j \\\frac{Z_{i}Z_{j}}{\| {\bf R}_{i} - {\bf R}_{j}\|} & \text{if } i \neq j\end{cases}\end{split}\]

You may use the following molecular geometries to test your function:

Water:

O 0.0 0.0 0.0
H -0.77 -0.573 0.0
H 0.775 -0.566 0.0

Caffeine:

C 3.184 1.228 0.0
N 2.239 0.145 0.001
C 2.488 -1.176 0.004
N 1.401 -1.918 0.003
C 0.383 -1.03 -0.002
C 0.869 0.27 -0.003
C 0.012 1.416 -0.002
O 0.357 2.585 -0.003
N -1.33 1.037 0.002
C -1.851 -0.253 0.001
O -3.047 -0.443 0.005
N -0.95 -1.296 0.003
C -1.43 -2.661 -0.009
C -2.327 2.096 0.0
H 3.838 1.161 0.871
H 3.794 1.2 -0.905
H 2.628 2.164 0.033
H 3.48 -1.583 0.007
H -0.746 -3.272 0.576
H -2.432 -2.674 0.417
H -1.468 -3.041 -1.031
H -2.942 2.026 0.897
H -2.971 1.993 -0.873
H -1.797 3.045 -0.025

Task 3.2: One-hot encoding

To learn categorial data, we need a representation that allows for regression to be done. One-hot encoding is such a representation, of which we will now consider a simplified, slightly overcomplete version for sake of simplicity of this exercise. In one-hot encoding, each category in the data is represented by a binary vector, with a length equal to the number of categories. Each vector has a value of 1 in the position corresponding to the category it represents, and a value of 0 in all other positions. For example, if we have three categories: alkanes, amines, and acids, then the one-hot encoding for alkanes would be [1,0,0], for amines [0,1,0], and for acids [0,0,1]. If there are multiple such categories, they will be appended: if we take the categories (alkanes, amines), then the one-hot encoding would be [1,0,0,0,1,0].

Write a function one_hot(categories: list[str], possible_categories: list[str]) -> list[list[int]] that converts a list of categories such as letters chosen from a set of possible categories into a one-hot encoded vector.

Examples:

one_hot(["A"], ["A", "B"]) -> [1,0]
one_hot(["A", "C"], ["A", "B", "C"]) -> [1,0,0,0,0,1]

Task 3.3: Invertible representations

Is the Coulomb matrix representation invertible, i.e. can you recover the molecule from the representation? Explain your answer.
Write a function sum_formula(coulomb_matrix: list[list[float]]) -> str that obtains the sum formula from the Coulomb matrix representation.

Exercise 4: Loss functions and their impact

Loss functions constitute the key metric that is minimized in order to improve a model. There are many different possible choices which impact the predictions of a model. This means that the loss function needs to be specified to fully qualify a model.

Different loss functions can yield highly disparate predictive power, exemplified in this work estimating bond dissociation energies.

Task 4.1: Loss functions

Let's fit a function of the form

\[y=a\exp(bx)\]

which can be logarithmized as

\[\log(y)=\log(a)+bx\]

Write loss functions for the squared error of the residuals in the original form and the logarithmized form. Use the following signature: original_loss(xs: np.ndarray, ys: np.ndarray, params: tuple[float]) -> float.

Task 4.2: Fit models

Write a function random_data(a: float, b: float, c: float) -> tuple(np.ndarray, np.ndarray) that returns the values 0-9 for xs and some exponential of the form

\[y=a\exp(bx)+c\epsilon\]

where epsilon is some elementwise normal distributed random noise. Fit the data using the loss functions of task 1 and scipy.optimize.minimize once for \(c=0\) and once for \(c\) being some non-zero value. Discuss the results.

Task 4.3: Motivation for loss functions

Consider the effects you saw in task 2. What could be possible motiviations to choose one loss function over another? Is there a best loss function? Which problem do you see for loss functions for multivariate models?

Exercise 5: Simple models

Most models for materials/chemistry applications build upon some description of a local environment. One way of doing so is to assign the same property to each similar environment. These kind of models can be refined by making the environment larger, e.g. by considering more neighbors. If the whole molecule is segmented into complete yet disjunct parts, this family of models is called a group contribution model/method. Often, models of increasing local environments yield acceptable results, but fail to converge, as typically many-body and long-range effects and interactions contribute significantly, especially for charged systems.

When considering the full environment, i.e. the whole molecule for small ones, it is possible to learn the 3D geometry from the graphs.

Task 5.1: Dressed atom from a database

The file bondorder.txt contains one line per molecule with three space-separated groups:

comma-separated elements
comma-separated bond order matrix as upper triangle (there are relevant numpy functions; note that bond orders of 1.5 are used to denote aromatic bonds)
energy in Hartrees

For water, the row would be

O,H,H 1,1,0 -74.9621236632301

Write a function that reads this file and estimates the energies using the dressed atom model, i.e. approximating the total energy as a sum of atom-wise contributions which only depend on the elements in question. For example, water has the sum formula of H2O, so the total energy of water can be approximated as twice the energy of a hydrogen plus once the energy of an oxygen. The atomic energies in this model are fitting parameters, which can be obtained via linear regression. Which relative error do you get for the energies?

Task 5.2: Bond counting

Now follow the same strategy as in task 1, but also consider all unique bonds, i.e. a CO bond, a CH bond, ... The total energy is now assumed to have one contribution from every element and another contribution from every bond. Which relative error do you get now?

Task 5.3: Generalisation

Give three different reasons (with examples) where this method must fail beyond the initial success you showed in task 1 and 2. Are these methods useful then? Discuss in the context of data standardisation.