1_slice_of_numpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d640cadb-69d0-44f7-950f-fc22e5f45af6",
   "metadata": {},
   "source": [
    "# A first slice of NumPy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33bd6805-a238-45f6-8702-6ef5b8e38eac",
   "metadata": {},
   "source": [
    "This notebook walks you through just those aspects of NumPy that you need to complete the steps in the charge distribution notebook.\n",
    "If you want to get a deeper understanding of how NumPy works and how to use it effectively, the [NumPy tutorials](https://numpy.org/learn/) are an excellent source of learning material.\n",
    "\n",
    "Whenever you want to read a bit more about any of the functions shown below, you can look them up on the [NumPy API reference](https://numpy.org/doc/stable/reference/index.html).\n",
    "(Click on the looking glass and type the name of a NumPy function.)\n",
    "\n",
    "Each section ends with 1 or 2 challenges.\n",
    "These must be solved without using any loops (`for` or `while`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0edcebdc-6028-4d6d-89b6-b194c062f860",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to install all required Python packages for both notebooks.\n",
    "%pip install -q numpy scipy matplotlib ipympl\n",
    "# In VSCode, you will need to \"Restart\" the notebook\n",
    "# before the newly installed backages can be imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25b7c6f1-9625-4f41-984d-31a95344e1c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# All import statements\n",
    "import urllib\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f9f4f83-27e6-42f0-8461-94f7f67bb019",
   "metadata": {},
   "source": [
    "## Lists versus arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f290ed3-3793-4805-bb2c-14be194b4725",
   "metadata": {},
   "source": [
    "A Python list is a very flexible data type that can store a sequence of values of any type, including combinations of different types. It is possible to add or remove any element. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10a74f4-e2c2-4753-9a17-54c9e0c1f959",
   "metadata": {},
   "outputs": [],
   "source": [
    "elements = [1, \"egg\", [4, 5, \"foo\"]]\n",
    "print(elements[0])\n",
    "print(elements[1])\n",
    "del elements[1]\n",
    "elements.append(\"hello\")\n",
    "elements.insert(1, \"inserted\")\n",
    "print(elements)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d0bf724-f6fa-4b39-84c2-7cfea59127ab",
   "metadata": {},
   "source": [
    "(NumPy) arrays are much more rigid: the number of elements is fixed when the array is created, and for the sake of efficient computation, all elements must be of the same type. The following example creates an array from a list of three values and shows how to access and modify its elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75d4f899-dbdc-4071-b54f-4d49fd8b9793",
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.array([1, 2, 4.0])\n",
    "print(arr)\n",
    "print(arr[0])\n",
    "arr[0] = 5\n",
    "print(arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d140ad2d-9306-4675-a93d-3e87258878c3",
   "metadata": {},
   "source": [
    "Note that all array elements are floats, even though some of the values in the first line are integers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cc7cf50-9041-44e9-b312-a5c2101928f1",
   "metadata": {},
   "source": [
    "Unlike lists, arrays can have multiple indexes, e.g. for rows and columns, whereas a list is always a 1D sequence. You can represent a matrix as a list of lists, but that already involves multiple lists. Given such a list of lists of values, you can construct the corresponding array, as in the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c0c1bcd-d604-45c6-9b71-4ba8f6bb5846",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array([[1, 2, 3], [4, 5, 6]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3c7c3ca-2a50-4b1d-abb9-113fe0cac293",
   "metadata": {},
   "source": [
    "## Array creation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c738582-2164-40b2-b464-31fe11213978",
   "metadata": {},
   "source": [
    "There are many ways to create arrays without first putting numbers into lists.\n",
    "The cells below show some examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc80b4cb-a8db-46cd-be37-35c819a933c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros((3, 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c23aae4-3ff1-4c36-9b9f-2101126d9b3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ones((4, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5ee8d48-3407-416b-a6b9-187f8985f0be",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arange(40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "868d6ac5-9e97-418a-8cf0-5ba01ac7d733",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linspace(-1, 1, 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7745df7-5064-4c14-a4dc-97a2051bb73f",
   "metadata": {},
   "source": [
    "SciPy also has several functions for constructing arrays for special purposes.\n",
    "One of them is a function to construct what is called a distance matrix.\n",
    "This is 2D array of all possible distances between two sets of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d298b30c-2580-4e22-b296-808092ae67bf",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "points1 = np.array([[1, 2], [4, 5], [3, 1]])\n",
    "points2 = np.array([[0, 2], [5, 5], [2, 2]])\n",
    "sp.spatial.distance_matrix(points1, points2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "515da5c9-07fd-4286-b0f4-dceb5eb7def8",
   "metadata": {},
   "source": [
    "**Challenge 1:**\n",
    "Construct a square array of `natom` rows and columns filled with zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d565d25-0295-47b9-b75c-6f512bcb844f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here.\n",
    "natom = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "709ce810-ca16-4448-988a-552e1a6fc00b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Challenge 2:** Construct a distance matrix for all distances between points in an array. (What do you expect to find on the diagonal?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4d59eef-fff9-42de-ae5c-5bad7188b405",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here.\n",
    "points = np.array([[7, 1], [8, 5], [1, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1611ef41-8230-4cb3-a7cc-5eccd59ce41e",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Slicing parts from arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35069d0-b9de-4e2e-909c-26fc16804678",
   "metadata": {},
   "source": [
    "You can slice arrays, just like lists, but with more possibilities of how to slice exactly.\n",
    "All examples below perform slicing operations on the following array `a`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c53b330-f3af-4ef3-a662-60162680840b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f83b9d8-47d8-4f07-986e-55bd93d469f9",
   "metadata": {},
   "source": [
    "An array behaves like a list of rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "692e56d0-781f-4562-99f7-2ecf9a7a3ead",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a30762c-7a97-4b74-b165-6fffd94e9026",
   "metadata": {},
   "source": [
    "You can slice columns in a similar way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f67d170f-11aa-401b-bec9-41551498c02b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[:, :2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e0553e4-c0f0-4180-9cd1-937abee29e25",
   "metadata": {},
   "source": [
    "This means: take all rows (first index), and only the first two columns (second index)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42ebed1d-c440-4efc-a7a3-8ce22f58d8c5",
   "metadata": {},
   "source": [
    "Slicing is not only useful for taking parts of arrays, you can also use slicing to assign values.\n",
    "For example, the following fills the third column with `-1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c95ff6-6a3c-4bfa-bb88-c6d3ed5c8269",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[:, 2] = -1\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7224d13b-ef0a-467b-9ba1-7b3b8e2f0afc",
   "metadata": {},
   "source": [
    "Similarly, the first two elements of the second row can be assigned `-2` and `-3`, respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c4e44b3-d3b4-47ca-a053-0a0737d7768d",
   "metadata": {},
   "outputs": [],
   "source": [
    "a[1, :2] = np.array([-2, -3])\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c968c110-2e9f-4fa4-b5c1-f66346feb4d2",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Challenge 3:**\n",
    "Construct a square matrix with `natom+1` rows and columns. Fill the last row and column with ones, except for the last diagonal elements. Make sure that all other elements are zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "780b20ea-fcf6-4ade-a3c7-f7ffb2ed9335",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here\n",
    "natom = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aecc08a-ab72-48a9-990c-a1b18165a7d9",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Array operations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8cfdbf5-8b44-4357-9da5-dbe9a28d9dc6",
   "metadata": {},
   "source": [
    "NumPy arrays understand mathematical operations, just like ordinary numbers.\n",
    "By default, all operations are element-wise, safe for some exceptions.\n",
    "The main advantage of using array operations is that they are faster than the equivalent implementation using a Python loop.\n",
    "\n",
    "The following examples all perform operations on the two arrays `a` and `b`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c89cfcfd-d9bc-44a8-9353-e8f961099506",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(5)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91a47b71-94af-4a24-8801-cbad2a0b219d",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.arange(5, 10)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4f4bc0f-af0e-42ad-a682-1016e2e9a377",
   "metadata": {},
   "outputs": [],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0779ec4-db1d-412f-96ba-740f21cdbae7",
   "metadata": {},
   "outputs": [],
   "source": [
    "a * b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8f5b5e3-827f-4f73-87bf-09af704106f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b63f75ac-76b0-43e6-9070-7b64f02a0256",
   "metadata": {},
   "source": [
    "NumPy can also combine scalar and array quantities in mathematical operations.\n",
    "Here are a few examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffa26501-66c6-4337-bd9a-742d2120e61f",
   "metadata": {},
   "outputs": [],
   "source": [
    "a**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc329a44-c057-4024-8662-03a3e451d9de",
   "metadata": {},
   "outputs": [],
   "source": [
    "1 / a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68838ce5",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "b + 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c462000-4a3b-433c-a49d-12e186fe0eb8",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "You can tell NumPy to ignore specific floating point errors with [`np.errstate`](https://numpy.org/doc/stable/reference/generated/numpy.errstate.html), e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10a34e0e-947d-4314-ab4b-2490386e6920",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "with np.errstate(divide=\"ignore\"):\n",
    "    print(1 / a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12885c0b-9728-446f-9646-cfb0836dd2fd",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Challenge 4:**\n",
    "Compute the element-wise inverse of the distance matrix.\n",
    "How would this be useful in finding the charge distribution that minimizes the potential energy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d5b0f06-34ec-4582-af15-38f5921c495b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here\n",
    "points = np.array([[7, 1, 0], [8, 5, 1], [1, 2, 1], [6, 3, 7]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b595980a-62c6-414b-abe9-140a0b95d9c7",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Challenge 5:** Compute for the same points the square matrix with damped Coulomb interactions from the theory slides:\n",
    "\n",
    "$$\n",
    "A_{ij} = \\frac{1 - \\exp(-r_{ij}/\\sigma)}{r_{ij}}\n",
    "$$\n",
    "\n",
    "with $\\sigma=2.5$.\n",
    "You can ignore the division by zero on the diagonal for now.\n",
    "To improve numerical stability, you can use [`np.expm1`](https://numpy.org/doc/stable/reference/generated/numpy.expm1.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "015af27a-5eca-438d-956d-0945f18ed2c1",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee17ffbc-5603-49ae-a6af-1c716494b31a",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Array reductions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33fee14f-9362-40a3-89b9-77da1f2f5ce6",
   "metadata": {},
   "source": [
    "There are many more array operations than those discussed so far.\n",
    "Some of the more interesting ones are array reductions, e.g. to compute the sum, mean, standard deviation etc. of values in an array.\n",
    "Some examples are shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96903f97-7e9c-4593-9d42-b362e16adb15",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(5)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e9b567b-1006-442c-9d48-6f3b34109f8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "a.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eec8bc3-719e-4056-99e4-14cfdde3a51d",
   "metadata": {},
   "outputs": [],
   "source": [
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf6fdb61-2841-4817-8a6e-5f262f8c0d34",
   "metadata": {},
   "outputs": [],
   "source": [
    "points.mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36e8abac-2a26-4d1b-93f6-2800f18adce6",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Challenge 6:** consider an array with atomic charges `[0.9, 0.2, 0.3, -0.4]`. Compute the total charge using NumPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79fc309f-92e1-446b-ac5e-3306f0a9268a",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d3bbc9b-24a0-4e48-8f2f-22ef1f144315",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Boolean array indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ef1fba2-35ce-44cd-ac32-5632f8ef894a",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "In addition to the slicing seen above, there are a few other mechanisms for taking parts of an array to construct a new array.\n",
    "One of these is the boolean index array: an array of booleans that can be used to select specific rows or columns.\n",
    "This is best illustrated with an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91475407-b7af-438f-9024-25b5e3b10c5a",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "a = np.array([-4, 3, 0, 1, -3, 1, 3])\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22cc1d52-070a-4855-8997-0fab7fa722f0",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We can easily construct a boolean array that contains `True` for all strictly positive elements of `a`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df9617fe-d723-4b7e-9ccb-eddb8f172122",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "mask = a >= 0\n",
    "mask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "519f0199-7f1a-4a3a-baf6-1f355d6fa339",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We can now use such a boolean array to select the values of `a` that are strictly positive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42e97704-b0cd-4b5d-853b-0fa3a53d1bee",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "a[mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a73afffb-bf62-47f5-921a-b31106aa6010",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This feature was used to create the visualization of half of the nanocrystal: a mask is used to select all atoms with a negative $x$-coordinate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10d58b0a-4e77-44c2-8b0f-fa1df476831b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Loading and saving arrays from/to text files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86156b28-1530-49bf-9b90-ec502485d538",
   "metadata": {},
   "source": [
    "For this section, you need the files `crystal.txt` and `plates.txt`.\n",
    "If you don't have them yet, you can download them by executing the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01222810-515f-4b69-ab71-ba324638dd46",
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_txt():\n",
    "    filenames = [\"crystal.txt\", \"plates.txt\"]\n",
    "    for filename in filenames:\n",
    "        print(f\"Downloading {filename}\")\n",
    "        url = f\"https://raw.githubusercontent.com/molmod/chargedist/refs/heads/main/{filename}\"\n",
    "        urllib.request.urlretrieve(url, filename)\n",
    "\n",
    "\n",
    "download_txt()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d99c46fe-bb1a-4d6b-a93b-d9fbceb270eb",
   "metadata": {},
   "source": [
    "In the next notebook, you will get files with XYZ coordinates of copper atoms.\n",
    "These files can be loaded into a NumPy array with [`np.loadtxt`](https://numpy.org/doc/stable/reference/generated/numpy.loadtxt.html) as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28e04465-1259-42f4-9e10-d82d4a0ae6a8",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "atcoords = np.loadtxt(\"crystal.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "926eb598-ef88-41da-ac93-e23e8a238b56",
   "metadata": {},
   "source": [
    "To get a first impression of the data, you can print the shape and the data type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20be878e-d921-4878-9bb2-879b720bbae5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(atcoords.shape)\n",
    "print(atcoords.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fd2ed71-0662-49b0-8388-a6b6addfa73f",
   "metadata": {},
   "source": [
    "The file contains a few thousand atomic coordinates. The columns correspond to $x$, $y$ and $z$ coordinates and the data type is a 64-bit floating point number."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7cfb692-1356-417a-b26f-03822db16ea5",
   "metadata": {},
   "source": [
    "Whenever you want to save an array to a file, you can do so with [`np.savetxt`](https://numpy.org/doc/stable/reference/generated/numpy.savetxt.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeefd61e-4e31-4b79-9229-7efa81cdaaa6",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Challenge 7:**\n",
    "Load the atomic coordinates from the file `plates.txt` and compute the center of mass of the first half of the atoms in this file.\n",
    "\n",
    "You should find something like `[-20.0 0.0 0.0]`, but there should be small deviations due to rounding errors.\n",
    "\n",
    "Hints:\n",
    "\n",
    "- All rows contain positions of copper atoms, all having the same mass, so the center of mass is just the average of the positions.\n",
    "- You may want to look up the documentation of [`np.mean`](https://numpy.org/doc/stable/reference/generated/numpy.mean.html). It is possible to compute all relevant averages in one line.\n",
    "- You need to divide the length of the array by two, but slicing requires integer indexes, so you need to use integer division, e.g. `length // 2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac011294-38b6-4eaa-8495-752491c78bbe",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af23fa47-82b2-4c65-bfe3-65389142294b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "783c07df-af58-477f-88bf-85b1bc6930bf",
   "metadata": {},
   "source": [
    "The [`numpy.linalg`](https://numpy.org/doc/stable/reference/routines.linalg.html) subpackage provides many linear algebra operations.\n",
    "To calculate equilibrium charge distributions, we need to solve a linear system of equations.\n",
    "Let's demonstrate this by solving the following system,\n",
    "which is taken from the example on Kirchhoff's laws in section 26.3 of Giancoli:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "I_1 + I_2 - I_3 = 0\\,\\mathrm{A} \\\\\n",
    "30 I_1 + 41 I_3 = 45\\,\\mathrm{A} \\\\\n",
    "21 I_2 + 41 I_3 = 125\\,\\mathrm{A}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "This can be written in matrix notation as\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 & 1 & -1 \\\\\n",
    "30 & 0 & 41 \\\\\n",
    "0 & 21 & 41 \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "I_1 \\\\\n",
    "I_2 \\\\\n",
    "I_3\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "0\\,\\mathrm{A} \\\\\n",
    "45\\,\\mathrm{A} \\\\\n",
    "125\\,\\mathrm{A}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "To solve this equation with NumPy, we need to create the corresponding matrix and vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7b20fc7-5e30-4861-a2e4-9e742f56ccc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "coeffs = np.array([[1, 1, -1], [30, 0, 41], [0, 21, 41]])\n",
    "rhs = np.array([0, 45, 125])\n",
    "sol = np.linalg.solve(coeffs, rhs)\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ce080ca-3d90-42ce-9edb-42d374b85a64",
   "metadata": {},
   "source": [
    "The result is indeed consistent with the theory lecture on Kirchhoff's laws."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79c4c5ac-ca38-466f-91c8-081641edf8fd",
   "metadata": {},
   "source": [
    "You can check the result by computing the matrix-vector product of the coefficients with the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "586be312-1275-4d4a-9910-d37a9fc3e7d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.dot(coeffs, sol) - rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "398ec0cf-f96d-4370-933b-a540d6ca7251",
   "metadata": {},
   "source": [
    "The residual is not exactly zero due to round-off errors.\n",
    "The `dot` function can compute matrix-matrix, matrix-vector and vector-vector products.\n",
    "(This is not to be confused with the geometric vector product, which is implemented in [`np.cross`](https://numpy.org/doc/stable/reference/generated/numpy.cross.html).)\n",
    "Python also provides a special syntax to write `np.dot` more concisely:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8b05c01-19d5-4d44-9cf8-421b8d8a116f",
   "metadata": {},
   "outputs": [],
   "source": [
    "coeffs @ sol - rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af71b5c2-60f8-4bf0-b715-0cf5475eadb4",
   "metadata": {},
   "source": [
    "## Varia"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fcfd63b-45b0-4f02-bb2c-74aaf5ebc77d",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "There is one more NumPy function that you will need to construct the equations to solve the charge distribution.\n",
    "The `np.fill_diagonal` function replaces all the diagonal elements of an array with a new value.\n",
    "For example, the following constructs the 5x5 identity matrix multiplied by 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8d126c1-0eb9-4d5e-b2ee-c6d25df15c26",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.zeros((5, 5))\n",
    "np.fill_diagonal(t, 2)\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f89026-0c20-4d00-afd3-1c04e41973f1",
   "metadata": {},
   "source": [
    "Note that the array `t` is modified.\n",
    "The `fill_diagonal` function has no return value."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}