2_charge_nano.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d3dd5b82-8c54-4b12-a2e6-5406974db3a7",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Charge distributions in copper nanostructures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edbab541-083b-471a-aa19-0684f7ec6e3b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "⚠️ Before you embark on this adventure, make sure you have worked your way through the preceding notebook `1_slice_of_numpy.ipynb` to learn the basics of NumPy.\n",
    "\n",
    "This notebook contains a series of problems to be solved by computing the charge distribution of a copper nanostructure.\n",
    "All required nanostructures can be loaded from `.txt` files, except for one which you are expected to create with NumPy code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7bb55b4-f315-4b4d-8324-ac5cd6c09e84",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Add additional imports here.\n",
    "import urllib\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d41e766c-026a-4c85-a300-cd27679cf0e6",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Configuration of your Notebook environment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ada3a19f-6a26-4a23-b612-6e981ad9721f",
   "metadata": {},
   "source": [
    "### Jupyter Lab and VSCode"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c04b8eb-eca0-4fbf-a8b9-f056555a9b01",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "You will need to install the following packages using [pip](https://pip.pypa.io/en/stable/installation/): [NumPy](https://numpy.org/), [SciPy](https://scipy.org/), [matplotlib](https://matplotlib.org/) and [ipympl](https://matplotlib.org/ipympl/).\n",
    "We assume that you have already done this before starting VSCode or Jupyter Lab. If not, execute the following cell and restart your notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1661572-0ee1-4281-a532-799ec5eb5192",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%pip install numpy scipy matplotlib ipympl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f681801-a219-41cb-b990-fdc85d798e8c",
   "metadata": {},
   "source": [
    "### Google Colab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "861f7de6-f1f6-4cfa-b548-d9e019ec862e",
   "metadata": {},
   "source": [
    "Google Colab comes with NumPy, SciPy and matplotlib packages preinstalled.\n",
    "However, it lacks support for interactive plots by default.\n",
    "This can be fixed with the code cell below.\n",
    "While it technically works with this fix,\n",
    "interactive plots on Google Colab still tend to be slow."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b150fdf",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The following code cell detects whether it is running on Google Colab or not, and executes the appropriate commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44dcd814-3ea7-4386-b41e-b3c6fcd63fe4",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    from google.colab import output\n",
    "\n",
    "    output.enable_custom_widget_manager()\n",
    "    %pip install -q ipympl\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f62b1ed9-b510-4cc8-bf9b-47ba4c8044d6",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### Any Jupyter Environment (VSCode, Jupyter or Google)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a2e072-dc91-4acd-b7da-3072b3584943",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The following command activates interactive plots (assuming `ipympl` was successfully installed)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "586d2740-7243-4c61-96a4-4d13eb2f4b2b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15438fe6-c580-43c1-a1d4-75c1483697df",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Download TXT files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4194bfe-b9ac-41ca-b251-99989e295414",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "You may already have the TXT files, in which case you can skip the download.\n",
    "For those using Google Colab, this download is definitely required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9bf5a63-19bc-47f1-b00f-97a436b577f7",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "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": "8800bd1d-ae1b-44b7-a463-70f483051f3d",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Visualization of the charge distribution of a nanostructure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc48cf4d-eb89-4c13-89c8-11af543d13f7",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This section introduces a convenient visualization of a charged nanostructure.\n",
    "Atoms are shown as circles and their color represents their charge.\n",
    "The first cell defines the function and the second cell shows how to use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb52ac12-6bc0-4e19-8584-ba012767d300",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_nano(label, atcoords, charges, size=100):\n",
    "    \"\"\"Plot the charge distribution of a nanocrystal.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    label\n",
    "        A unique string, has to be different for different plots.\n",
    "    atcoords\n",
    "        An array with Cartesian coordinates of the atoms, shape (natom, 3).\n",
    "    charges\n",
    "        An array with atomic charges, shape (natom,).\n",
    "    size\n",
    "        The size of the atomic spheres.\n",
    "    \"\"\"\n",
    "    # Matplotlib 3D scatter plot\n",
    "    natom = atcoords.shape[0]\n",
    "    num = f\"label_{natom}\"\n",
    "    plt.close(num)\n",
    "    fig = plt.figure(figsize=(6, 5), num=num)\n",
    "    ax = fig.add_subplot(projection=\"3d\")\n",
    "    scale = abs(charges).max()\n",
    "    sc = ax.scatter(\n",
    "        atcoords[:, 0],\n",
    "        atcoords[:, 1],\n",
    "        atcoords[:, 2],\n",
    "        c=charges,\n",
    "        s=size,\n",
    "        ec=\"k\",\n",
    "        cmap=\"coolwarm_r\",\n",
    "        depthshade=False,\n",
    "        vmin=-scale,\n",
    "        vmax=scale,\n",
    "    )\n",
    "    fig.colorbar(sc, label=\"Atomic charge [e]\")\n",
    "    qtot = charges.sum()\n",
    "    ax.set_title(f\"Cu$_{{{natom}}}$ nanostructure, $q_\\\\mathrm{{tot}}$ = {qtot:.2f} e\")\n",
    "    ax.set_xlabel(\"x [$a_0$]\")\n",
    "    ax.set_ylabel(\"y [$a_0$]\")\n",
    "    ax.set_zlabel(\"z [$a_0$]\")\n",
    "    ax.dist = 10\n",
    "    ax.set_box_aspect(\n",
    "        aspect=(np.ptp(atcoords[:, 0]), np.ptp(atcoords[:, 1]), np.ptp(atcoords[:, 2]))\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f57270fd-1cc4-4c58-a17f-7d5437824fd5",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The following is a small demo of the visualization function.\n",
    "It loads the nanocrystal from disk and visualizes it, using the x-coordinate as a surrogate for the partial charges.\n",
    "\n",
    "If all goes well, you should be able to rotate the 3D view by dragging the mouse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c60dd81c-eeea-4dee-9085-3907586521c3",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "atcoords = np.loadtxt(\"crystal.txt\")\n",
    "plot_nano(\"demo\", atcoords, charges=atcoords[:, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76cd351-9d5c-4551-9add-f0a54b21ae8a",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Challenge 1: the charge distribution of a nanocrystal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb424559-74f5-427d-b83e-522626801674",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Now it is up to you: you should be able to combine all the information from the theory slides and the previous notebook to construct the correct linear system of equations to calculate the charges.\n",
    "Do this calculation for a crystal with a net charge of $+10\\,\\mathrm{e}$.\n",
    "Once you have solved the equations, you can visualize the atomic charges using the `plot_nano` function.\n",
    "You should get a result similar to the \"goals\" slide from the theory part.\n",
    "\n",
    "You can check your answer by computing the potential energy of the optimal charge distribution. This should be $U_\\text{mol}\\approx-2.8571\\,\\mathrm{E_h}$.\n",
    "\n",
    "While solving and visualizing the charges, try to answer the following questions:\n",
    "\n",
    "- **Q1.**\n",
    "  What features of the charge distribution can you explain with what you have learned from the electrostatics chapters?\n",
    "- **Q2.**\n",
    "  Construct a boolean index array to select atoms with a negative $x$-coordinate.\n",
    "  Use this boolean array to visualize only one half of the crystal, so that you can look inside.\n",
    "  Can you explain the internal charge distribution?\n",
    "- **Q3.**\n",
    "  Compute the energy of the charged particle and compare it with the energy of a charged metal sphere of approximately the same radius.\n",
    "  Do you get consistent results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c282391b-dfa3-4a53-94c4-0ac016745123",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de78118c-8b19-4f79-9b9c-9101bb1fc79e",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The following cell contains a unit test to perform a few simple sanity checks to confirm that the implementation is correct."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c5472ec-e547-43bb-92fb-c24020450db8",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Challenge 2: A nanoscale capacitor consisting of two crystaline slabs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aa9fc8b-c716-4eda-9854-8c3c16302191",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The file `plates.txt` contains two parallel copper plates that can be used to model a nanoscale capacitor.\n",
    "Each plate consists of the same number of copper atoms.\n",
    "\n",
    "Load this file and compute a charge distribution when one plate is constrained to have a charge of $+10\\,e$ and the other plate is constrained to have a charge of $-10\\,e$.\n",
    "\n",
    "As you solve and visualize the charges, try to answer the following questions:\n",
    "\n",
    "- **Q4.**\n",
    "  How can you derive the capacitance? If implemented correctly, you should find a capacitance $C\\approx 6.2533\\,\\mathrm{e^2\\,E_h^{-1}}$.\n",
    "- **Q5.**\n",
    "  Try to estimate the capacitance by approximating this capacitor as an ideal flat plate capacitor.\n",
    "  Do you get consistent results?\n",
    "- **Q6 (difficult).**\n",
    "  What is the physical significance of the Lagrange multipliers in this case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff4405dd-4f4e-4883-b3c1-5d4392a09203",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccf9747b",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Optional extra challenges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c770557e",
   "metadata": {
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "For those with good programming skills, here are some additional challenges to keep you entertained. To get acceptable performance, you have to be careful how you vectorize the calculations with NumPy.\n",
    "\n",
    "- Visualize the electric field in the $(x, y)$-plane of the capacitor with a [stream plot](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.streamplot.html).\n",
    "  Do you see a region where the field is homogeneous?\n",
    "\n",
    "- Visualize the isopotential surfaces for the capacitor model.\n",
    "  You can do this using the [Marching Cubes](https://en.wikipedia.org/wiki/Marching_cubes) algorithm as implemented in [scikit-image](https://scikit-image.org/docs/stable/auto_examples/edges/plot_marching_cubes.html#marching-cubes).\n",
    "  Do you see flat surfaces between the plates?\n",
    "  "
   ]
  }
 ],
 "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
}