In part based on Fortran code from Furio Ercolessi and the CHEM126 course of Kalju Khan.

```
In [1]:
```from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())

```
Out[1]:
```

In lab 1 we worked with the Lennard-Jones potential. Other simple potentials include the **harmonic** potential

with harmonic constant $k$ and equilibrium position $r_{\text{eq}}$, the **Kratzer** potential

with dissociation energy $D_0$, and the **Morse** potential

where $\alpha$ is the Morse parameter.

For example, carbon monoxide (CO) consists of one carbon atom (atomic weight 12.011) and one oxygen atom (atomic weight 15.999) separated by 1.1283 Angstroms. This can be modelled by

- a harmonic potential with $k = 2743.0$ and $r_{\text{eq}} = 1.1283$ Angstroms, or
- a Kratzer potential with $D_0 = 258.9$ kcal/mol and the same equilibrium distance, or
- a Morse potential with $\alpha = 2.302$ Angstroms${}^{-1}$, and the same dissociation energy and equilibrium distance.

Using standard minimization algorithms (e.g. `scipy.optimize.minimize`

) show that, for suitable initial guesses, all of the above potentials have minima at the expected location.

By plotting the potentials, understand why the success of the algorithms changes.

Briefly test how changing the potential changes the number of steps required by the algorithm.

Also compute the force $\frac{dV}{dr}$ at the minimum point for each potential to check that it vanishes as expected.

```
In [2]:
```%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
rcParams['figure.figsize'] = (12,6)
from scipy.optimize import minimize

```
In [2]:
```

*not* use periodic boundaries! Try using $\Delta t = 0.001$ up to $t=1$. Start from the equilibrium position and show that it does not evolve. The move the location of the oxygen atom slightly and see how it evolves.

```
In [2]:
```

The paper of Praprotnik, Janezic, and Mavri suggests a potential for a single water molecule (one oxygen atom and two hydrogen with atomic weights $1.008$) of

$$ \begin{equation} V = \sum_{l = 1}^2 D_0 \left[ 1 - e^{\alpha \Delta r_{OH_l}} \right]^2 + \frac{1}{2} k_{\theta} \Delta r_{HH}^2 + k_{r \theta} \Delta r_{HH} \left( \Delta r_{OH_1} + \Delta r_{OH_2} \right) + k_{rr} \Delta r_{OH_1} \Delta r_{OH_2}. \end{equation} $$The first term is the Morse potential again, and the other terms are the intra-molecule bond potentials. $\Delta r_{OH_l} = r_{OH_l} - r_{OH_{\text{eq}}}$ is the stretch in the distance between the oxygen atom and the $l^{\text{th}}$ hydrogen atom $r_{OH_l}$ and its equilibrium value, and $\Delta r_{HH} = r_{HH_l} - r_{HH_{\text{eq}}}$ is the stretch in the distance between the hydrogen atoms. The parameter values are

- $D_0 = 101.9188$ kcal/mol.
- $\alpha = 2.567$ Angstrom${}^{-1}$
- $k_{\theta} = 328.645606$ kcal/mol/Angstrom${}^2$
- $k_{r \theta} = -211.4672$ kcal/mol/Angstrom${}^2$
- $k_{rr} = 111.70765$ kcal/mol/Angstrom${}^2$
- $r_{OH_{\text{eq}}} = 1$ Angstrom
- $r_{HH_{\text{eq}}} = 1.633$ Angstrom

Fix the oxygen atom at the origin and start the hydrogen atoms at $(\pm 0.8, 0.6, 0)$. Use the minimization techniques to find the equilibrium position of the atoms, comparing against the known structure of water.

**NOTE**: if using `scipy`

's minimization routine, you may have to increase the tolerance as high as $10^{-4}$ to make it converge.

```
In [2]:
```

Using your integrator, evolve the molecule and see how it behaves.

```
In [2]:
```