week 5 updates

Author: bjmorgan

CH40208/geometry_optimisation/analytical_solution_for_a_harmonic_potential.md

@@ -26,4 +26,4 @@ Your function should take three arguments as input; $r$, $k$, and $r_0$; and ret
 
 2. Plot this function for H<sub>2</sub> ($r_0$ = 0.74 &#8491;, $k$ = 36.0 eV &angst;<sup>&minus;2</sup>) for $0.38\leq r \leq 1.1$.
 
-3. Use your function to calculate the potential energy at $r=r_0$. Add a point to your plot at $(U(r_0), r_0)$ and confirm visually that this is the minimum of the potential energy surface.
+3. Use your function to calculate the potential energy at $r=r_0$. Add a point to your plot at $(r_0, U(r_0))$ and confirm visually that this is the minimum of the potential energy surface.

CH40208/geometry_optimisation/geometry_optimisation.md

@@ -23,12 +23,12 @@ Let us start with the simplest possible case: that of a diatomic molecule. For a
 
 In this simple case, our potential energy surface is just a one-dimensional function $U(r)$ that gives us the potential energy for any given internuclear separation $r$. The equilibrium geometry corresponds to the value of $r$ where $U(r)$ reaches its minimum.
 
-This concept of a potential energy surface generalises to more complex systems: for any molecular or crystal system, we can define a function $U(\left\{r\right\})$ that gives us the potential energy for any set of atomic positions $\left\{r\right\}$. Our challenge then becomes:
+This concept of a potential energy surface generalises to more complex systems. For any molecular or crystal system, we can define a function $U(\left\{r\right\})$ that gives us the potential energy for any set of atomic positions $\left\{r\right\}$. Our challenge then becomes:
 Given the ability to calculate $U(\left\{r\right\})$ for any set of atomic positions $\left\{r\right\}$, how can we algorithmically find the set of positions that minimises $U(\left\{r\right\})$?
 This is fundamentally a mathematical optimisation problem. However, it has several characteristics that make it particularly challenging:
 1. The function $U(\{r\})$ is typically not known analytically — we can only evaluate it at specified points.
 2. For all but the simplest systems, we are dealing with many dimensions ($3N-6$ internal coordinates for a molecule with $N$ atoms).
 3. The potential energy surface often has multiple local minima, making it crucial to distinguish between local and global minimum energy structures.
 
-For now we will only worry about point 1., and will return to points 2. and 3. later.
+For now, we will only focus on point 1 and return to points 2 and 3 later.
 

CH40208/geometry_optimisation/newton_raphson_method.ipynb

@@ -74,13 +74,33 @@
    "source": [
     "### Exercise\n",
     "\n",
+    "In this exercise you will explore how the Newton-Raphson method performs on a harmonic potential, where we know the analytical solution is $r = r_0 = 0.74$ Å.\n",
+    "\n",
+    "**Part 1: Define the second derivative**\n",
+    "\n",
     "1. Write a function to calculate the second derivative of a harmonic potential:\n",
     "\n",
-    "$$U^{\\prime\\prime} = k$$\n",
+    "$$U^{\\prime\\prime}(r) = k$$\n",
+    "\n",
+    "Note that for a harmonic potential, the second derivative is simply the force constant $k$ and does not depend on $r$.\n",
+    "\n",
+    "**Part 2: Demonstrate one-step convergence**\n",
+    "\n",
+    "2. Write code to perform Newton-Raphson optimisation on your harmonic potential starting from $r = 1.0$ Å. Your code should:\n",
+    "   - Calculate $U'(r)$ and $U''(r)$ at the starting position\n",
+    "   - Apply the Newton-Raphson equation once: $r_1 = r_0 - \\frac{U'(r_0)}{U''(r_0)}$\n",
+    "   - Print both the starting position and the position after one step\n",
+    "   - Compare the result to the analytical solution ($r_0 = 0.74$ Å)\n",
     "\n",
-    "2. Show that the Newton-Raphson method finds the optimal bond length for your harmonic potential in one step, irrespective of your starting point.\n",
+    "3. Repeat this calculation for at least three different starting positions:\n",
+    "   - $r = 0.5$ Å (starting below the minimum)\n",
+    "   - $r = 1.0$ Å (starting above the minimum)  \n",
+    "   - $r = 1.5$ Å (starting well above the minimum)\n",
     "\n",
-    "*Hint: Consider why the Newton-Raphson method works perfectly for quadratic functions.*"
+    "**Questions to consider:**\n",
+    "- Does each starting point converge to exactly $r = 0.74$ Å in one step?\n",
+    "- Why does the Newton-Raphson method work perfectly for a harmonic potential? (Hint: Consider what happens when you substitute the derivatives of a quadratic function into the Newton-Raphson equation. The method assumes the function is locally quadratic—what happens when the function is *exactly* quadratic everywhere?)\n",
+    "- Will this one-step convergence property hold for other potential energy functions like the Lennard-Jones potential? Why or why not?"
    ]
   },
   {

CH40208/geometry_optimisation/scipy_optimize_minimize.ipynb

@@ -247,6 +247,14 @@
     "\n",
     "12. How many distinct minima does your scatter plot reveal? (Look for horizontal bands of points.)\n",
     "\n",
+    "```{note}\n",
+    "You may notice points clustered near both $\\theta \\approx +3.14$ (near $+\\pi$) and $\\theta \\approx -3.14$ (near $-\\pi$). Before counting these as separate minima, remember that dihedral angles are periodic: $\\theta = +\\pi$ and $\\theta = -\\pi$ represent the same molecular conformation (both are 180°). \n",
+    "\n",
+    "Check your potential energy plot from Part A: you should see that $U(+\\pi) = U(-\\pi)$. These two sets of points represent convergence to the same minimum approached from different directions, not two different minima.\n",
+    "\n",
+    "This is a consequence of how the optimiser works: it finds a nearby minimum by following gradients, and doesn't \"wrap around\" the periodic boundary. Starting from $\\theta = +2$ radians naturally leads to the minimum expressed as $\\theta \\approx +\\pi$, whilst starting from $\\theta = -2$ radians leads to the same minimum expressed as $\\theta \\approx -\\pi$.\n",
+    "```\n",
+    "\n",
     "13. For each minimum, estimate the range of starting angles that lead to that minimum. These ranges are called the \"basin of attraction\" for each minimum.\n",
     "\n",
     "14. Based on your results, if you wanted to find the global minimum (the lowest energy structure), what strategy would you need to use? (Hint: look back at your potential energy plot from Part A to see which minimum has the lowest energy.)\n"