week 5 updates

Author: bjmorgan

CH40208/geometry_optimisation/Lennard_Jones_optimisation.ipynb

@@ -5,34 +5,124 @@
    "id": "455b7ed0-22a2-4d49-852c-db66b191c027",
    "metadata": {},
    "source": [
-    "### Exercise: Geometry Optimisation of a Lennard-Jones Potential\n",
+    "## Synoptic Exercise: Geometry Optimisation of a Lennard-Jones Potential\n",
     "\n",
-    "The potential energy surface for a pair of “Lennard–Jones” is given by\n",
+    "The Lennard-Jones potential describes the interaction between a pair of non-bonded atoms:\n",
     "\n",
-    "$$U_\\mathrm{LJ} = \\frac{A}{r^{12}} - \\frac{B}{r^6}$$\n",
+    "$$U_\\mathrm{LJ}(r) = \\frac{A}{r^{12}} - \\frac{B}{r^6}$$\n",
     "\n",
-    "#### 1.\n",
+    "where $r$ is the interatomic separation. The repulsive term ($r^{-12}$) represents electron cloud overlap at short distances, whilst the attractive term ($r^{-6}$) represents van der Waals interactions.\n",
     "\n",
-    "Plot this function for $A$ = 1 × 10<sup>5</sup> eV Å<sup>12</sup> and $B$ = 40 eV Å<sup>6</sup> for the range $3.6 \\leq r \\leq 8.0$.\n",
+    "In this exercise, you will use both gradient descent and Newton-Raphson methods to find the equilibrium separation for a Lennard-Jones pair, compare their performance, and explore how the choice of starting position affects convergence.\n",
     "\n",
+    "We will use the parameters $A$ = 1 × 10<sup>5</sup> eV Å<sup>12</sup> and $B$ = 40 eV Å<sup>6</sup>.\n",
     "\n",
-    "#### 2.\n",
+    "The derivatives of the Lennard-Jones potential are:\n",
     "\n",
-    "The first derivative of this potential energy function is given by\n",
+    "**First derivative:**\n",
+    "$$U'_\\mathrm{LJ}(r) = -12 \\frac{A}{r^{13}} + 6\\frac{B}{r^7}$$\n",
     "\n",
-    "$$U^\\prime_\\mathrm{LJ} = -12 \\frac{A}{r^{13}} + 6\\frac{B}{r^7}$$\n",
+    "**Second derivative:**\n",
+    "$$U''_\\mathrm{LJ}(r) = 156\\frac{A}{r^{14}} - 42\\frac{B}{r^8}$$\n",
     "\n",
-    "Write code to find the lowest-energy interatomic separation for a pair of Lennard–Jones atoms (with the potential parameters given above) using gradient descent. Start with $r = 5.0$ and $\\alpha=100$. How does your optimisation code perform with different values of $\\alpha$ and different starting values of $r$ (try $r$ = 3.2 Å, $r$ = 4.4 Å, and $r$ = 6.0 Å)?\n",
+    "### Part A: Visualising the Potential Energy Surface\n",
     "\n",
-    "#### 3.\n",
+    "Before attempting optimisation, it is important to understand the shape of the potential energy surface.\n",
     "\n",
-    "Write code to perform geometry optimisation of your Lennard-Jones potential using the Newton-Raphson method.\n",
+    "1. Write a function `lennard_jones(r, A, B)` that calculates $U_\\mathrm{LJ}(r)$ for a given separation $r$ and parameters $A$ and $B$.\n",
     "\n",
-    "The second derivative of $U_\\mathrm{LJ}$ is\n",
+    "2. Create an array of $r$ values from 3.6 Å to 8.0 Å using `np.linspace()` with at least 200 points.\n",
     "\n",
-    "$$U^{\\prime\\prime} = 156\\frac{A}{r^{14}} - 42\\frac{B}{r^8}$$\n",
+    "3. Plot the potential energy as a function of $r$ for the given parameters. Label your axes appropriately and include a title.\n",
     "\n",
-    "Investigate how the Newton-Raphson method performs with starting points of $r$ = 3.2 Å, $r$ = 4.4 Å, and $r$ = 6.0 Å."
+    "4. By examining your plot:\n",
+    "   - Estimate the location of the minimum (the equilibrium bond length)\n",
+    "   - Identify the approximate energy at the minimum\n",
+    "   - Note where the potential energy crosses zero (where attractive and repulsive forces balance to give $U = 0$)\n",
+    "\n",
+    "### Part B: Gradient Descent Optimisation\n",
+    "\n",
+    "Now we will use gradient descent with an adaptive step size to find the minimum.\n",
+    "\n",
+    "5. Write a function `lennard_jones_gradient(r, A, B)` that calculates the first derivative $U'_\\mathrm{LJ}(r)$.\n",
+    "\n",
+    "6. Implement a gradient descent loop that:\n",
+    "   - Starts from $r = 5.0$ Å\n",
+    "   - Uses the update rule: $r_{i+1} = r_i - \\alpha U'(r_i)$\n",
+    "   - Uses $\\alpha = 100$ as the learning rate\n",
+    "   - Runs for a maximum of 50 iterations\n",
+    "   - Stops early (using `break`) when $|U'(r)| < 0.001$\n",
+    "   - Stores the position and gradient at each iteration in lists (for later analysis)\n",
+    "\n",
+    "7. After your loop completes, print:\n",
+    "   - The final predicted equilibrium separation\n",
+    "   - The number of iterations required\n",
+    "   - The final gradient value\n",
+    "\n",
+    "8. Create a plot showing how the position $r$ changes with iteration number. This visualises the convergence path.\n",
+    "\n",
+    "**Questions to consider:**\n",
+    "- How many iterations were required to converge?\n",
+    "- Does the convergence appear smooth or does the algorithm oscillate?\n",
+    "- Look at your convergence plot: does the step size decrease as you approach the minimum? (Remember: the step size is $\\alpha \\times |U'(r)|$, and $U'(r)$ should decrease as you approach the minimum.)\n",
+    "\n",
+    "### Part C: Exploring Different Starting Positions with Gradient Descent\n",
+    "\n",
+    "The choice of starting position can significantly affect optimisation performance.\n",
+    "\n",
+    "9. Test your gradient descent code with three different starting positions:\n",
+    "   - $r = 3.2$ Å (close to the minimum, approaching from below)\n",
+    "   - $r = 4.4$ Å (at the minimum, approximately)\n",
+    "   - $r = 6.0$ Å (far from the minimum, approaching from above)\n",
+    "\n",
+    "   For each starting position, record:\n",
+    "   - The final optimised position\n",
+    "   - The number of iterations required\n",
+    "   - Whether the algorithm converged successfully (i.e., did it find $|U'(r)| < 0.001$ within 50 iterations?)\n",
+    "\n",
+    "10. Summarise your results, showing the three starting positions and their outcomes.\n",
+    "\n",
+    "**Questions to consider:**\n",
+    "- Do all starting positions converge to the same minimum?\n",
+    "- Which starting position converges most quickly? Can you explain why based on your plot from Part A?\n",
+    "- Did any starting positions fail to converge within 50 iterations?\n",
+    "- How does the learning rate $\\alpha = 100$ perform for the different starting positions? Based on the gradient magnitudes at different $r$ values, does this $\\alpha$ seem appropriate everywhere, or only in certain regions?\n",
+    "\n",
+    "### Part D: Newton-Raphson Optimisation\n",
+    "\n",
+    "The Newton-Raphson method uses both first- and second-derivative information to converge more rapidly.\n",
+    "\n",
+    "11. Write a function `lennard_jones_second_derivative(r, A, B)` that calculates $U''_\\mathrm{LJ}(r)$.\n",
+    "\n",
+    "12. Implement a Newton-Raphson loop that:\n",
+    "    - Starts from $r = 5.0$ Å\n",
+    "    - Uses the update rule: $r_{i+1} = r_i - \\frac{U'(r_i)}{U''(r_i)}$\n",
+    "    - Runs for a maximum of 50 iterations\n",
+    "    - Stops early when $|U'(r)| < 0.001$\n",
+    "    - Stores the position and gradient at each iteration\n",
+    "\n",
+    "13. After your loop completes, print the final equilibrium separation and number of iterations required.\n",
+    "\n",
+    "14. Create a plot comparing the convergence paths of gradient descent (from Part B) and Newton-Raphson, both starting from $r = 5.0$ Å. Plot both curves on the same axes showing position vs. iteration number.\n",
+    "\n",
+    "**Questions to consider:**\n",
+    "- How does the convergence speed of Newton-Raphson compare to gradient descent?\n",
+    "- Why doesn't Newton-Raphson converge in a single step for the Lennard-Jones potential, unlike the harmonic potential?\n",
+    "- Which method required fewer function evaluations? (Remember: gradient descent needs one derivative per step, whilst Newton-Raphson needs two.)\n",
+    "\n",
+    "### Part E: Comparing Methods Across Different Starting Positions\n",
+    "\n",
+    "15. Test Newton-Raphson with the same three starting positions you used for gradient descent in Part C:\n",
+    "    - $r = 3.2$ Å\n",
+    "    - $r = 4.4$ Å\n",
+    "    - $r = 6.0$ Å\n",
+    "\n",
+    "    For each starting position, record the final position and the number of iterations required.\n",
+    "\n",
+    "**Questions to consider:**\n",
+    "- Which method is more robust to different starting positions?\n",
+    "- For which starting position(s) does Newton-Raphson show the greatest advantage?\n",
+    "- Based on your results, which method would you recommend for optimising a Lennard-Jones potential? What factors influenced your decision?"
    ]
   },
   {