{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Workshop General Info\n",
"## Connecting to Notebooks\n",
"\n",
"You can run this notebook locally on your HetSys laptop or in any other environment with Python, SciPy, NumPy and ASE present. \n",
"\n",
"\n",
"## Submitting the Assignment\n",
"\n",
"Before you turn this problem in, make sure everything runs as expected. First, **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart) and then **run all cells** (in the menubar, select Cell$\\rightarrow$Run All). For cells that invoke long DFT runs, you may wish to leave the notebook in a state where the calculations are not executed again but rather read in from the previous runs. If you do that, please ensure the code you submit matches that which generated the output from the DFT run. There are no marks for \"what my code used to do before I changed it\".\n",
"\n",
"Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\". Answer any text-based questions and double-check that you have given any variables you create the names you are asked to give them.\n",
"\n",
"A submission that passes all the tests will pass overall, but there are often marks for hidden tests you cannot see.\n",
"\n",
"I have to ask that people submit their work through two channels: sorry for the inconvenience but there are separate routes for auto-grading (via nbgrader) and eventual feedback and incorporation into module grades.\n",
"\n",
"For submission via the nbgrader interface: once finished, go back to the \"Assignments\" tab and click submit. You can do this again if you are not happy with your answers. This submits a version of the notebook with all outputs cleared.\n",
"\n",
"For the sake of auditability and providing a \"paper trail\" for external examiners, we also need a finished, exported html version submitted to moodle. Please choose \"Download as\" from the file menu and select \".html\" format, then upload that via Moodle.\n",
"\n",
"Please enter your name and SCRTP username in the box below in case files go astray. Feedback will be released through the nbgrader interface: you can go back to the Assignments tab again to download it. It will also be uploaded to Moodle.\n",
"\n",
"\n",
"## Remote Access\n",
"\n",
"From outside the Warwick network, the range of ports on which we run these notebooks is blocked. However, you can still get in if you establish an \"ssh tunnel\" if you know what port your jupyter notebook is running on (it tells you when you launch it - the default is 8888 but it will increment if that is already in use).\n",
"\n",
"For example, to map port localhost:8888 on stan4 to localhost:8888 on your own machine you would do:\n",
"\n",
" ssh -L localhost:8888:localhost:8888 stan4.csc.warwick.ac.uk\n",
"\n",
"Then you would point your browser at\n",
"\n",
" https://localhost:8888\n",
" \n",
"to launch the notebook \n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"NAME = \"\"\n",
"SCRTP_USERNAME = \"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "890cfe975efade216af3e9b045f98cf7",
"grade": false,
"grade_id": "cell-ce8e1032d27449af",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"***\n",
"***\n",
"\n",
"# HetSys Python Induction & Introduction to Jupyter Notebooks\n",
"\n",
"## Scope of this Workshop\n",
"\n",
"This Workshop will introduce you to using Jupyter Notebooks, by showing you how to use some simple (even trivial) python commands inside notebooks. This is intended to act as a refresher for those who have not used python for a while or have recently learned it from scratch. Those who have extensive Python experience will likely find it trivial. We are basically checking you understand how to take an equation and turn it into some simple python code.\n",
"\n",
"The workshop script is implemented as a Jupyter notebook, which is the same format we will use throughout the PX911 module. These notebooks are based on the python language. These nbgrader notebooks consist of text-only cells with instructions (white background), editable cells (grey background) where you can adjust the code as required to carry out the instructions, and non-editable cells (also grey) which contain tests that will check the result of your code. Some of these tests are hidden so you cannot just copy in the answer from the test.\n",
"\n",
"##### This workshop is not for credit! It is intended to ensure everyone is able to get off to a good start with the notebooks in the PX911 module."
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "24f8f079b4447a9a9e49796ef57d1c70",
"grade": false,
"grade_id": "cell-cd6ce0875213cfe1",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## What is Jupyter?\n",
"\n",
"Project Jupyter exists to develop open-source software, open-standards, and services for interactive computing across dozens of programming languages.\n",
"\n",
"The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text. Uses include: data cleaning and transformation, numerical simulation, statistical modeling, data visualization, machine learning, and much more.\n",
"\n",
"\n",
"https://jupyter.org/"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "f054bc94ea3d6186bb4bc089ee901f9b",
"grade": false,
"grade_id": "cell-4db2cc2455f1761a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"***\n",
"\n",
"## Importing Modules"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "391efc02e92720ca8b4cc37d2db2c4a3",
"grade": false,
"grade_id": "cell-3fc35daabab1406c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Our first exercise is to calculate a numerical value for the Bohr radius $a_0$, the expectation value of the distance between the electron and the nucleus in a hydrogen atom in its ground state, from fundamental constants.\n",
"\n",
"We can use the formula: $a_0 = \\displaystyle \\frac{4 \\pi \\epsilon_0 \\hbar^2}{m_e e^2}$ which is encountered in the derivation of hydrogen-like orbitals (ignoring relativistic effects).\n",
"\n",
"To use this we need numerical values for all the constants. The constants module from the SciPy library has appropriate definitions, so we will \"import\" it. The list of the available constants is at this link:\n",
"\n",
"https://docs.scipy.org/doc/scipy/reference/constants.html\n",
"\n",
"The constants are defined in SI units. That is straightfoward but is often not desirable, as we will see later.\n",
"\n",
"NB: I don't want to see any cut-and-paste defining of the values of physical constants in a notebook like this! That is almost never something you should be doing."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "d9d5a994c0ffc28a5f16ab9ba6b802cd",
"grade": false,
"grade_id": "cell-842a757d61391d15",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# import the SciPy constants module, and give it the short-name spc\n",
"import scipy.constants as spc\n",
"\n",
"# This module contains many useful scientific constants, such as epsilon_0, which you can now access as spc.epsilon_0\n",
"# \n",
"# Write a line of code that puts your result in a variable called 'a0'\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()\n",
"\n",
"# If you have set the variable, it will be printed here (otherwise this will cause an error)\n",
"print(f'a0 = {a0} m')\n",
"\n",
"# Note I used an f-string for printing the result - a handy modern way of printing formatted data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "0db7ee7d4e0b3e8d295934e35801cf42",
"grade": true,
"grade_id": "cell-3968b17fc5016c71",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"# Now we test your answer against a reference value, imported from a different library,\n",
"# namely the Atomic Simulation Environment, which we use heavily in PX911.\n",
"# ASE's internal units are Angstrom, so we multiply by 10^-10 to get an answer in meters\n",
"from ase.units import Bohr\n",
"\n",
"assert(a0 - Bohr*10**-10 < 1e-10)"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "61d000ad1c61b33c2c1f40e4fc924dd5",
"grade": false,
"grade_id": "cell-280a6de315a5b569",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Now we will calculate the \"Rydberg unit of energy\", 1 $\\text{Ry}$, which is equal to the binding energy of an electron in a hydrogen atom.\n",
"\n",
"An appropriate formula for the Rydberg is $\\text{Ry} = \\displaystyle \\frac{m_e e^4}{8\\epsilon_0^2 h^2}$\n",
"\n",
"Some electronic structure codes use the Rydberg as their internal unit of energy. Others use the Hartree, where $1\\text{ Ha} = 2\\text{Ry}$.\n",
"\n",
"The ASE uses the electron volt. $1\\text{ Ry} \\simeq13.606\\text{ eV}$ and $1\\text{ Ha} \\simeq 27.211\\text{ eV}$ are useful conversion factors to recall."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "ac759cbf383c3dec9082cea7f47d8aa8",
"grade": false,
"grade_id": "cell-749cb75578c3f226",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Calculate the Rydberg unit of energy and put it in the variable 'Ry'\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()\n",
"\n",
"# If you have set the variable, it will be printed here (otherwise this will cause an error)\n",
"print(f'Ry = {Ry} J')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "b0e60785eac4ae64157248f0d2fcf677",
"grade": true,
"grade_id": "cell-9404ea26ab4a6c94",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"# In this cell we test your answer expressed in electron volts, by comparing to ASE.\n",
"# However, this time the test is \"hidden\" (because I am mean).\n",
"# In my notebooks I will always flag up the existence of a hidden test, but I won't\n",
"# tell you what it is as that would defeat the point. These hidden tests are designed to\n",
"# make you think a bit more carefully and check your work rather than just \"code to pass\n",
"# the test\"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "ace8007e9a8d2f0626850bd3435691f4",
"grade": false,
"grade_id": "cell-a5ff93eeea012a2b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"***\n",
"\n",
"## Lists and Dictionaries"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we'll have a look at some python datatypes. First, the list - an ordered collection of elements (of any type) indexed by integers.\n",
"\n",
"We can define an empty list by assigning a variable to empty square brackets: []\n",
"```\n",
"my_list = []\n",
"```\n",
"\n",
"We can add a new entry to a list with ```append()```\n",
"```\n",
"my_list.append(25)\n",
"```\n",
"\n",
"There is loads more that is worth knowing about lists, particularly \"list comprehensions\" but I don't want to get bogged down. For more info see here: https://docs.python.org/3/tutorial/datastructures.html\n",
"\n",
"As an exercise, I want you to make a list called ```symbols``` containing the chemical symbols for the first 8 elements of the periodic table."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "404da00f9f196754e66b0adeb90ddfd3",
"grade": true,
"grade_id": "cell-045eea0ecf55c539",
"locked": false,
"points": 1,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# I've defined an empty list for you called \"symbols\" and added the first element\n",
"symbols = []\n",
"symbols.append('H')\n",
"\n",
"# You need to add the other 7\n",
"# Aim for conciseness of the solution\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "54f17c6270ad71713429f78ff77c8461",
"grade": true,
"grade_id": "cell-db80360cfea7f789",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"from ase.data import chemical_symbols\n",
"assert(symbols == ['H','He','Li','Be','B','C','N','O'])"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "ad23cadf2fb3ef41b4c06fda927279ea",
"grade": false,
"grade_id": "cell-8bd9a6fb1984ca6f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Another very useful and flexible built-in Python datatype is the \"dictionary\".\n",
"\n",
"Dictionaries can be created by placing a comma-separated list of `key: value` pairs within braces. Empty braces create an empty dictionary:\n",
"```\n",
"my_empty_dict = {}\n",
"my_initialised_dict = {\"hello\": \"goodbye\", 123: \"fish\"}\n",
"print(my_initialised_dict[\"hello\"]\n",
"```\n",
"which would print \"goodbye\".\n",
"\n",
"Dictionaries are particularly useful for defining \"lookup tables\" where the \"key\" that returns a particular value is not necessarily an integer: it could be a string, or any other value which is \"hashable\".\n",
"\n",
"In the box below, please create a dictionary which, if supplied with a \"key\" which is the single-letter string label for an angular momentum subshell, can supply the corresponding angular momentum quantum number $l$. I would like it to work from \"s\" up to \"g\", ie from l=0 to l=4.\n",
"\n",
"https://en.wikipedia.org/wiki/Azimuthal_quantum_number - see the table."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "9d7028e26cfabdd8fb14353750580c26",
"grade": false,
"grade_id": "cell-8fefeab26129067e",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Please store your dictionary as l_values\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "1b280274b8320f89a4ed96ad4bfcfa79",
"grade": true,
"grade_id": "cell-5254e602576046bd",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert(list(l_values.keys()) == [\"s\",\"p\",\"d\",\"f\",\"g\"])\n",
"assert(list(l_values.values()) == [0,1,2,3,4])"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "219d6485ceeb85e31e7da383e850c475",
"grade": false,
"grade_id": "cell-b2e9c406108e21a6",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"***\n",
"\n",
"## NumPy"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "9252fe1c1dcc416ac8c88e6d63315f98",
"grade": false,
"grade_id": "cell-9b950cc2a1cd0093",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"There is a huge amount to be said about NumPy, a widely-used python package adding support for large, multi-dimensional arrays and matrices. Crucially, it is the way Python (normally a relatively slow, interpreted language) can be made to perform in certain tasks nearly as efficiently as a compiled program in a low-level language like C or Fortran. I cannot do it justice in a few boxes here, but we will just aim to give you some familiarity with its routines when using NumPy arrays for holding positions of atoms, and \n",
"\n",
"A much \n",
"https://numpy.org/devdocs/user/quickstart.html\n",
"\n",
"In the box below, we first create an empty 3x3 array of floating-point values. Then we put some numbers in that represent coordinates of the atoms in a water molecule (in Angstrom units) w\n",
"\n",
"Your job is to shift them 3A along the z-axis by creating a vector of length 3 along z and adding it to each column of the position vector."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "3ac375c202f0cf31163dcca721448f80",
"grade": false,
"grade_id": "cell-a91ee18116c473e4",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Import numpy and give it a short name \"np\"\n",
"import numpy as np\n",
"\n",
"# Here we create an empty numpy array of 3 column vectors of length 3.\n",
"water_pos = np.zeros([3,3])\n",
"\n",
"# Set the O-H bond length (same as ASE's pre-made model)\n",
"rOH = 0.96856502\n",
"# Set the H-O-H bond angle (same as ASE's pre-made model)\n",
"theta = np.radians(104/2)\n",
"\n",
"# Put the oxygen atom at the origin\n",
"water_pos[0] = [0,0,0]\n",
"# Hydrogen atoms at appropriate locations in y-z plane\n",
"water_pos[1] = [0, rOH*np.sin(theta),-rOH*np.cos(theta)]\n",
"water_pos[2] = [0,-rOH*np.sin(theta),-rOH*np.cos(theta)]\n",
"\n",
"# Now translate these vectors by 3 along the z direction\n",
"# You could use a for loop for this for conciseness.\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()\n",
"\n",
"print(f'Positions of the atoms after translation:\\n {water_pos}')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "9df894df0cff096741c2b58ab9f3cd5f",
"grade": true,
"grade_id": "cell-f4fb4e0f29184d45",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"# Once again I will import a result from ASE to test your solution\n",
"# ASE's pre-made water molecule is at a slight offset from z=0 so we correct that first and shift by 3 along Z.\n",
"from ase.build import molecule\n",
"m = molecule('H2O')\n",
"# Put oxygen atom at z=3 (it starts at z=0.119262)\n",
"m.translate([0,0,-0.119262 + 3])\n",
"assert((water_pos-m.positions < 1e-6).all())"
]
},
{
"cell_type": "markdown",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "0ff368ae90aac64b3522cd5129f33d19",
"grade": false,
"grade_id": "cell-f8fe9fb98d42c788",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We are going to make a matrix representing the Hamiltonian matrix of a two-level system, with some small component of off-diagonal coupling. You need to find a numpy routine that will give you the eigenstates of this Hamiltonian.\n",
"\n",
"The premise is that we have two \"sites\" in our system, each represented with an orthogonal basis function $f_i$ for $i=1,2$. The first has a site energy $\\sigma_1$ and the second has site energy $\\sigma_2$. They are \"coupled\" by a coupling term $\\delta$. The Hamiltonian matrix therefore looks like this:\n",
"\n",
"$ \\qquad \\qquad \\mathbf{H} = \n",
"\\begin{pmatrix}\n",
"\\sigma_1 & \\delta\\\\\n",
"\\delta & \\sigma_2\n",
"\\end{pmatrix}$\n",
"\n",
"The box below sets this up and puts in some arbitrary values. You need to find the numpy routine that gives you the *eig*envalues $\\epsilon_n$ and *eig*envectors $c^n_i$ and puts the results in ```eigval``` and ```eigvec```. Try looking in the \"linalg\" module.\n",
"\n",
"These are the solutions of the secular equation:\n",
"\n",
"$ \\qquad \\qquad \\mathbf{H}.\\mathbf{c}_n = \\epsilon_n \\mathbf{c}_n$ for $n=1,2$\n",
"\n",
"and the corresponding wavefunctions would be $\\psi_n (\\mathbf{r}) = \\sum_i c^n_i f_i(\\mathbf{r})$. This represents just about the simplest possible meaningful quantum mechanics calculation. Real electronic structure calculations just involve scaling things up!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "26261f35137477d6b6ea6385dd8f79df",
"grade": false,
"grade_id": "cell-a0668b76eda59efa",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"\n",
"# On-site and coupling terms of the Hamiltonian matrix\n",
"sigma1 = 0.5\n",
"sigma2 = 0.55\n",
"delta = 0.02\n",
"\n",
"# Create a 2x2 matrix\n",
"ham = np.array([[sigma1,delta],[delta,sigma2]])\n",
"\n",
"# Find the eigenvalues\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()\n",
"\n",
"print(f'The eigenvalues are {eigval[0]} and {eigval[1]}')\n",
"\n",
"print(f'The eigenvectors are {eigvec[0]} and {eigvec[1]}')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "fe693120700f90bac5b77b375fe30791",
"grade": true,
"grade_id": "cell-0eed48d11b7b8ef0",
"locked": true,
"points": 1,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"# Check against correct values\n",
"eigval_sol = np.array([0.49298437881283574, 0.5570156211871643])\n",
"eigvec_sol = np.array([[-0.94362832, -0.33100694],[ 0.33100694, -0.94362832]])\n",
"tol = 1e-8\n",
"assert(((eigval-eigval_sol)