Computational Chemistry 2025

Aims to be an 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

Introduction to Python: Three without given solutions.

Exercise 2: Little Potential

Task 1: Free range atoms

An (overly simplistic) model is to estimate the energy of each atom to be \(-\frac{1}{2}Z^{7/3}\). Write a python function get_energy(list_of_elements: list[str]) -> float which returns the total energy for elements until \(Z=10\). For carbon monoxide, the function would be called like get_energy(['C', 'O']).

Optional: for those who are a bit more fluent in python, write the function to instead accept stoichiometries as strings, such as "C6H6", i.e. get_energy(elements: str) -> float.

Task 2: Harmonic co-existence

In task 1, there were no bonds. A simple model of bonds is a harmonic potential. Consider carbon monoxide this time. Assume a potential where the bond \(b\) has an equilibrium bond length of \(b_0\) and the energy depending on the actual bond length is \(a(b-b_0)^2\) with some constant \(a>0\). Write a function co_energy(b: float, a: float, b_0: float) -> float. Choose some values for both \(a\) and \(b_0\) and plot the potential energy surface if the central carbon atom is moved along the axis of a linear carbon dioxide molecule. You may use matplotlib for plotting:

import matplotlib.pyplot as plt
positions = ...
energies = ...
plt.plot(positions, energies, label="energy")
plt.xlabel(...)
plt.ylabel(...)

Task 3: Make it double

What could possibly go wrong in this harmonic potential model for a) dimers as above and b) more complex molecules. Find at least three situations where the model must yield unphysical results for any real system.

Exercise 3: Force field

We build a small force field. If you are advanced in Python: solve the tasks in an object-oriented manner with the classes Atom, Bond, ClassicalPotential instead.

Task 1: Book keeping

For the atoms, we need to know element identities encoded as nuclear charges and their positions. Write a function which takes a molecule as string in the following XYZ format and extracts positions and nuclear charges, returning them as numpy arrays.

O        0.000      0.000      0.000
H       -0.770     -0.573      0.000
H        0.775     -0.566      0.000

Task 2: Connecting the dots

Write a function get_energy(molspec: str, bonds: list[tuple[int, int]]) which takes a molecule in the XYZ format (re-use the function from task 1) and a list of zero-based bond indices to evaluate the following potential for each bond of length \(b\): \(E(b) = Z_1Z_2(b-1)^2\) add the Coulomb repulsion for all pairwise atoms of distance \(r\): \(E_C = Z_1Z_2/r\). For water, one would call this as follows:

water = """O        0.000      0.000      0.000
H       -0.770     -0.573      0.000
H        0.775     -0.566      0.000"""
get_energy(water, [(0, 1), (0, 2)])

Task 3: So much potential

Think about the energy landscape defined by the potential in task 2 and predict the minimum energy geometry. Then visualise the potential energy landscape by plotting results of get_energy. If you scale the Coulomb repulsion by some prefactor, how does the result change? What is missing in this potential and how does it show?

Exercise 4: Function optimisation

Task 1: Step by step

One way to obtain the gradient of an arbitrary function \(f\) is called finite differences. This method evaluates \(f\) at two places with a very small displacement \(h\):

\[\left.\frac{\partial f}{\partial x}\right|_{x=x_0} \approx \frac{f(x_0+h)-f(x_0-h)}{2h}\]

Implement a function gradient_descent(h, step_scaling) which optimizes the function \(f(x)=x^2\) using gradient descent and gradients from finite differences. The starting point should be at \(x=2\), and reasonable settings are \(h=0.001\) and \(s=0.1\).

Optional: Write the function to accept an arbitrary function to optimize \(f\) and a starting point \(x_\textrm{i}\): gradient_descent(f, x_i) where \(h\) and step_scaling should be chosen automatically.

Task 2: Line rider

Use your code from task 1 and add a test whether the position found is an actual minimum. This is done by testing whether the second derivative at that point is positive. You can obtain the second derivative of a function with finite differences like this:

\[\left.\frac{\partial^2 f}{\partial x^2}\right|_{x=x_0} \approx \frac{f(x_0+h)-2f(x_0)+f(x_0-h)}{h^2}\]

Test your code on the function \(f\) from task 1 and on \(g(x)=x^3\), each time starting from \(x_\textrm{i}=1\). What do you observe?

Task 3: If you know, you know

Minimize the function \(\sin(x)/x\) for different starting points. What do you observe? Which consequence do you see for geometry optimisation of molecules?

Exercise 5: Clusters

Task 1: Build it

Implement the Lennard-Jones potential for a cluster of atoms, setting \(\sigma\) and \(\epsilon\) to 1. Make the function accept a 3 by \(n\) array of coordinates. Plot the potential to verify its shape. Where is the analytical minimum?

Task 2: Group it

Without running any code, name the minimum energy configurations of clusters with three and four atoms.

Then use scipy.optimize.minimize to check your results starting from some random coordinates. For comparison, the energies for the clusters are \(-3\) and \(-6\), respectively. Then use the same code to obtain the solution for five (energy \(-9.103852\)) and six atoms (energy \(-12.712062\)). Can you find the geometries?

Task 3: Fix it

Switch to scipy.optimize.basinhopping and try to find the minimum for six, ten (\(-28.422532\)) and 15 (\(-52.322627\)) atoms. Look at the solutions you know. Can you hypothesize how to improve the search efficiency without just making more steps? Hint: what property are solutions likely to have?

Exercise 6: Optimisation Strategies

Task 1: Initial Steps

Many optimizers require an initial geometry to start from. Propose (do not implement!) a procedure how one could find an initial geometry for a molecule if only the molecular graph was given. Make sure your method can deal with rings.

Task 2: Minima

Once we have one minimum geometry for a molecule, how could we either find the other thermally populated ones or could exclude that there are more? Describe one idea with and one without molecular dynamics.

Task 3: Transitions

Transition state geometries are saddle points and as such are not strictly downhill. How could we obtain such a geometry for a certain reaction if reactant and product are available? Hint: what if you went on a hike with the largest rubber band ever made?