{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bond Length vs Amplitude of Vibration\n",
"\n",
"As a first example, let's look at the hydrogen molecule."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Jupyter and Python\n",
"\n",
"This is a [jupyter](https://jupyter.org/) [ipython](https://ipython.org/) notebook. This allows us to mix text/latex descriptions with [python](https://www.python.org/) code.\n",
"\n",
"We'll be mainly relying on the following libraries:\n",
"- [numpy](https://www.numpy.org/)\n",
" - This allows us to easily work with arrays of data.\n",
"- [scipy](https://www.scipy.org/)\n",
" - This gives us many numerical analysis tools.\n",
"- [matplotlib](https://matplotlib.org/)\n",
" - This allows us to easily generate plots.\n",
"\n",
"You may not be familiar with python, but if you've used Matlab or Mathematica, you should find the process here somewhat similar, but with slightly different commands.\n",
"\n",
"Jupyter and Python are freely licensed, so if you want to install them on your own computer, you can either follow the above links or you can get a complete package from [anaconda](https://www.anaconda.com/download/) which will give you everything you need on Windows, Mac or Linux. There is extensive documentation and usage examples at the above links.\n",
"\n",
"If you want information on what a function call we're using is doing, you can enter `help(functionname)` in a cell and it will output a detailed description of the input, output and method used.\n",
"\n",
"To evaluate the code in a jupyter cell press `ctrl+enter`, or press `shift+enter` to evaluate the cell and automatically move to the next cell."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# First we tell Python to load in functionality from a set of python libraries for numerical\n",
"# analysis, plotting, and some scientific linear algebra, interpolation, optimization and\n",
"# integration libraries.\n",
"# You'll need to be sure to evaluate this cell or the rest of the code in this notebook won't\n",
"# work.\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.linalg\n",
"import scipy.interpolate\n",
"import scipy.optimize\n",
"import scipy.integrate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## H2 Bond Length\n",
"\n",
"To examine this, we've used DFT to calculate the energy as a function of the separation between two H atoms on a very dense grid."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Load the data from the file to an array.\n",
"# Note bond lengths are in Bohr and energies in Ry.\n",
"evb = np.loadtxt(\"etot_v_bl.dat\")\n",
"\n",
"# bond_lengths are in the first column\n",
"bond_lengths = evb[:, 0]\n",
"\n",
"# energies are in the second column\n",
"energies = evb[:, 1]\n",
"\n",
"# plot the energy vs bond length\n",
"plt.plot(bond_lengths, energies)\n",
"plt.xlabel(\"Bond Length (Bohr)\")\n",
"plt.ylabel(\"Energy (Ry)\")\n",
"# After setting up the plot with the various commands above, plt.show() will display it.\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Let's find the separation and vibration frequency to begin with.\n",
"# To find the minimum, we can interpolate with a spline and find where this is minimized\n",
"\n",
"evb_interp = scipy.interpolate.InterpolatedUnivariateSpline(bond_lengths, energies)\n",
"\n",
"# The interpolation behaves like a scalar function, and we can use minimize_scalar\n",
"# to find where this is minimized. This is our equilibrium separation.\n",
"# This function returns an object with various information about the minimum.\n",
"# We just need the actual value which is stored as \"x\", so we append \".x\" to the function call.\n",
"\n",
"evb_min = scipy.optimize.minimize_scalar(evb_interp).x\n",
"\n",
"print(\"Equilibrium Separation:\", evb_min, \"Bohr\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Now let's find the oscillation frequency.\n",
"# We can find this from the second derivative at the minimum point.\n",
"# Our interpolation has a \"derivative\" method, we can use this to\n",
"# find the second derivative at the separation we found above.\n",
"k = evb_interp.derivative(2)(evb_min)\n",
"\n",
"# Now we need to calculate the frequency and convert the units\n",
"# into something understandable. We have a spring constant\n",
"# and we know f = 1/2pi sqrt(k/m).\n",
"# Define a unit conversion factor to convert the result to THz\n",
"Ry_to_J = 2.1798741e-18\n",
"Bohr_to_m = 5.29177e-11\n",
"\n",
"amu_to_kg = 1.66054e-27\n",
"f_to_THz = np.sqrt((Ry_to_J / Bohr_to_m ** 2) / amu_to_kg) * 1e-12\n",
"\n",
"# Since the energy vs bond length was calculated with one atom\n",
"# fixed and the other moving, we should use the reduced mass here.\n",
"# This is m/2 for two equal masses.\n",
"reduced_mass_amu = 1.00794 / 2\n",
"\n",
"freq_THz = 1/(2*np.pi) * np.sqrt(k/reduced_mass_amu) * f_to_THz\n",
"print(\"Oscillation frequency: \", freq_THz, \"THz\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We've now found the separation and oscillation frequency for small vibrations. However, as we can see from the plot we produced earlier, the potential is not harmonic. This means that for larger amplitude oscillations, the average separation will be different from the minimum value. As the potential looks like it's tilted to the right, we might also expect the average separation to increase with increasing amplitude of vibration.\n",
"\n",
"Let's try to calculate the average separation as a function of the average energy of the molecule."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's now look at the equation of motion of a particle in this potential. With\n",
"\n",
"$$ F = -\\frac{dV}{dx} = m\\ddot{x} $$\n",
"\n",
"we can simulate how a particle will behave in this potential.\n",
"\n",
"We'll first write a python function that we perform leapfrog integration of the time evolution of a mass in this potential, returning the position and velocitiy at each timestep."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define a function that will calculate the value of position as a function\n",
"# of time by solving the differential equations for position and velocity\n",
"# using the leapfrog method.\n",
"\n",
"def get_x_v_t(potential_deriv, mass, x0, v0, dt, nt):\n",
" '''Leapfrog integration to find x(t) and v(t).'''\n",
" # Initialize empty arrays that will store x and v for each timestep.\n",
" x = np.empty(nt)\n",
" v = np.empty(nt)\n",
" # Add in the initial conditions.\n",
" x[0] = x0\n",
" # Put velocity out of sync by half a step (leapfrog method)\n",
" v[0] = v0 + 0.5 * dt * potential_deriv(x[0]) / mass\n",
"\n",
" # Use a for loop to iterate over timesteps, updating the velocity\n",
" # and position each time.\n",
" for it in range(1, nt):\n",
" v[it] = v[it-1] - dt * potential_deriv(x[it-1]) / mass\n",
" x[it] = x[it-1] + dt * v[it]\n",
" return x, v"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use this function to see how different amplitude vibrations behave."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# If we keep v0 set to zero, the amplitude will only depend\n",
"# on the starting position x0. So we can use some arbitrary units\n",
"# for our mass and t that will depend on each other, but will not\n",
"# change the behaviour with amplitude.\n",
"\n",
"dt = 0.1\n",
"tmax = 30\n",
"nt = int(tmax/dt)\n",
"\n",
"# Set up an array storing the t values (for plotting)\n",
"tvals = np.linspace(0, tmax, nt)\n",
"\n",
"# Get the first derivative of the potential\n",
"V_deriv = evb_interp.derivative(1)\n",
"\n",
"# And now look x(t) for three different amplitudes:\n",
"xvals1, vvals = get_x_v_t(V_deriv, 1, 1.6, 0.0, dt, nt)\n",
"xvals2, vvals = get_x_v_t(V_deriv, 1, 1.8, 0.0, dt, nt)\n",
"xvals3, vvals = get_x_v_t(V_deriv, 1, 2.0, 0.0, dt, nt)\n",
"\n",
"# Plot them all together so we can compare\n",
"plt.plot(tvals, xvals1, label=\"x0 = 1.6\")\n",
"plt.plot(tvals, xvals2, label=\"x0 = 1.8\")\n",
"plt.plot(tvals, xvals3, label=\"x0 = 2.0\")\n",
"plt.xlabel(\"time\")\n",
"plt.ylabel(\"Position (Bohr)\")\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, the larger the amplitude the smaller the frequency (i.e. the fewer total oscillations in the time simulated), confirming that this is an anharmonic oscillation. The first peak of the vibration is reached at roughly the same time for all starting amplitudes, while at longer time we can see the peaks for the higher amplitude vibrations come later.\n",
"\n",
"While you may have seen a model of the vibration of a diatomic molecule where a quadratic potential, depending on some spring constant was assumed, we now have the full potential from DFT instead, which turns out to be slightly different from a quadratic, so we get the above behaviour.\n",
"\n",
"What about the bond length though? You may be able to see that for the high amplitude vibration, the trajectory is \"wider\" at high separations, and \"narrower\" at low separations, which should mean we're spending more time at high separations.\n",
"\n",
"If we calculate a very long trajectory (relative to the oscillation period) we can obtain the average position directly."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Say we calculate a long trajectory and find the average atomic position:\n",
"xvals1, vvals = get_x_v_t(V_deriv, 1, 1.6, 0.0, dt, 10000)\n",
"xvals2, vvals = get_x_v_t(V_deriv, 1, 1.8, 0.0, dt, 10000)\n",
"xvals3, vvals = get_x_v_t(V_deriv, 1, 2.0, 0.0, dt, 10000)\n",
"\n",
"print(\"Average postion if the atom starts at 1.6 Bohr: \", xvals1.sum() / len(xvals1))\n",
"print(\"Average postion if the atom starts at 1.8 Bohr: \", xvals2.sum() / len(xvals2))\n",
"print(\"Average postion if the atom starts at 2.0 Bohr: \", xvals3.sum() / len(xvals3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An alternative way to find the average position for a given starting value would be to integrate over a half-period of oscillation as follows:\n",
"\n",
"$$\n",
"\\langle x \\rangle = \\frac{2}{\\tau} \\int_0^{\\tau/2} x(t) dt\n",
"$$\n",
"\n",
"where $\\tau$ is the period of oscillation. Say we just want to use our calculated potential, we can use $v(x) = dx/dt$, so that $dt = dx/v$ and obtain\n",
"\n",
"$$\n",
"\\langle x \\rangle = \\frac{2}{\\tau} \\int_{x_1}^{x_2}\\frac{x}{v}dx\n",
"$$\n",
"\n",
"where $x_1$ and $x_2$ are the limits of oscillation. This is a little tricky to simply plug into some numerical integration routine however, as $v(x_1)=v(x_2)=0$ so that the integral blows up at the limits."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}