Skip to content

Commit 04db40c

Browse files
bjmorganbjmorgan
committed
WIP
1 parent 48383ae commit 04db40c

3 files changed

Lines changed: 541 additions & 58 deletions

File tree

‎CH40208/geometry_optimisation/gradient_descent_method.ipynb‎

Lines changed: 90 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -9,26 +9,56 @@
99
"\n",
1010
"A more efficient approach is to use information about the slope of the potential energy surface to guide our search. \n",
1111
"\n",
12-
"The **gradient descent method** (also called the **steepest descent method**) works by analogy to releasing a ball on a hill and letting it roll to the bottom. At any point on our potential energy surface, the **gradient** tells us which direction is “uphill” and how steep the surface is at this point. With this information, we can step in the opposite direction (i.e., downhill), then recalculate the gradient at our new position, and repeat until we reach a point where the gradient is $\\sim0$.\n",
12+
"The **gradient descent method** (also called the **steepest descent method**) works by analogy to releasing a ball on a hill and letting it roll to the bottom. At any point on our potential energy surface, the **gradient** tells us which direction is \"uphill\" and how steep the surface is at this point. With this information, we can step in the opposite direction (i.e., downhill), then recalculate the gradient at our new position, and repeat until we reach a point where the gradient is $\\sim0$.\n",
1313
"\n",
1414
"The simplest implementation of this method is to move a fixed distance every step.\n",
15-
"#### Algorithm\n",
15+
"\n",
16+
"#### Algorithm (Fixed Step Size)\n",
1617
"\n",
1718
"1. Start at an initial guess position $r_0$\n",
18-
"2. Calculate the gradient $U^\\prime = \\frac{dU}{dr}$ at this point\n",
19-
"3. Move in the direction opposite to the gradient (i.e., downhill), by some amount $\\Delta r$.\n",
20-
"4. Repeat until the gradient is close to zero.\n",
19+
"2. Calculate the gradient $U^\\prime(r) = \\frac{\\mathrm{d}U}{\\mathrm{d}r}$ at this point\n",
20+
"3. Move in the direction opposite to the gradient (i.e., downhill):\n",
21+
" $$r_{i+1} = r_i - \\Delta r \\times \\mathrm{sign}(U^\\prime(r_i))$$\n",
22+
" where $\\Delta r$ is a fixed step size and $\\mathrm{sign}(U^\\prime)$ gives us the direction (+1 or -1)\n",
23+
"4. Repeat until the gradient is close to zero\n",
24+
"```{note}\n",
25+
"In molecular simulations, we often express the gradient of the potential energy in terms of the **force**. The force is defined as the negative gradient of the potential energy:\n",
26+
"\n",
27+
"$$F(r) = -\\frac{\\mathrm{d}U}{\\mathrm{d}r} = -U^\\prime(r)$$\n",
28+
"\n",
29+
"This relationship means that forces naturally point \"downhill\" towards lower energy, whilst the gradient points \"uphill\" towards higher energy. When particles move in the direction of the force, they move towards lower potential energy.\n",
30+
"```"
31+
]
32+
},
33+
{
34+
"cell_type": "markdown",
35+
"id": "cd899cfa-e2eb-4a33-b0c2-20860f6635cf",
36+
"metadata": {},
37+
"source": [
38+
"#### Exercise: Fixed Step Size Gradient Descent\n",
2139
"\n",
22-
"#### Exercise:\n",
2340
"Write a function to calculate the first derivative of a harmonic potential:\n",
2441
"\n",
25-
"$$U^\\prime = 2k(r - r_0)$$\n",
42+
"$$U^\\prime(r) = k(r - r_0)$$\n",
43+
"\n",
44+
"Using this function, write code to perform a gradient descent search to find the minimum of your harmonic potential energy surface. Use $r=1.0$ Å as your starting position.\n",
45+
"\n",
46+
"Your code should:\n",
47+
"- Use a `for` loop with a **maximum of 50 iterations** to prevent infinite loops\n",
48+
"- Store the position and gradient at each iteration in lists (for plotting later)\n",
49+
"- Update the position at each step using: $r_{i+1} = r_i - \\Delta r \\times \\mathrm{sign}(U^\\prime(r_i))$\n",
50+
"- Print the iteration number, position, and gradient at each step\n",
51+
"- Stop early (using `break`) when $|U^\\prime(r)| < 0.001$\n",
52+
"- Report whether the algorithm converged or hit the maximum iteration limit\n",
2653
"\n",
27-
"Using this function, write code to perform a gradient descent search, to find the minimum of your harmonic potential energy surface. Use $r=1.0$ Å as your starting position, and a step size $\\Delta r=0.1$ Å.\n",
54+
"**Part 1:** Start with $\\Delta r = 0.01$ Å. Does the algorithm converge? How many iterations does it require?\n",
2855
"\n",
29-
"_Note_: You should not need to run for more than 10 steps.\n",
56+
"**Part 2:** Now try $\\Delta r = 0.1$ Å to see if a larger step size converges faster. What behaviour do you observe? Does it still converge?\n",
3057
"\n",
31-
"What happens using this algorithm when you get close to the minimum? What happens if you decrease the step size to $\\Delta r=0.01$ Å?"
58+
"**Questions to consider:**\n",
59+
"- Why does the larger step size cause oscillation?\n",
60+
"- The minimum is at r = 0.74 Å and we start at r = 1.0 Å. Can you explain why Δr = 0.01 Å converges but Δr = 0.1 Å does not, based on this distance?\n",
61+
"- What would you predict for Δr = 0.05 Å?"
3262
]
3363
},
3464
{
@@ -37,40 +67,77 @@
3767
"metadata": {},
3868
"source": [
3969
"#### Rescaling the Step Size\n",
40-
"Using a fixed step-size makes our method very sensitive to the choice of step-size. Large step sizes will overshoot the minimum and then oscillate back and forth. Small step sizes will get closer to the minimum, but at the cost of needing to perform many more calculations.\n",
4170
"\n",
42-
"A common approach designed to address this problem is to rescale the step size for each iteration, based on how far we think we are from the minimum. A simple model is that we can expect the gradient to be steep if we are a long way from the minimum, but shallow if we are already close to the minimum, so we make our step-size proportional to the (negative) local gradient.\n",
71+
"Using a fixed step size makes our method very sensitive to the choice of step size. As you have seen, large step sizes overshoot the minimum and oscillate back and forth, whilst small step sizes converge more reliably but require many more iterations.\n",
72+
"\n",
73+
"A common approach designed to address this problem is to rescale the step size for each iteration based on how far we think we are from the minimum. A simple model assumes that the gradient will be steep if we are far from the minimum, but shallow if we are already close. Therefore, we make our step size proportional to the local gradient magnitude.\n",
74+
"\n",
75+
"Instead of moving a fixed distance $\\Delta r$, we take a step proportional to the gradient:\n",
4376
"\n",
4477
"$$\\Delta r_i = -\\alpha U^\\prime(r_i)$$\n",
4578
"\n",
46-
"or \n",
79+
"or equivalently, proportional to the force:\n",
4780
"\n",
4881
"$$\\Delta r_i = \\alpha F(r_i)$$\n",
4982
"\n",
50-
"where $F(r_i)$ is the **force** (i.e., the negative gradient of the energy) at $r_i$.\n",
83+
"where:\n",
84+
"- $\\alpha$ is a constant called the **learning rate** or **step size parameter**\n",
85+
"- $F(r_i) = -U^\\prime(r_i)$ is the force at position $r_i$\n",
5186
"\n",
52-
"#### Exercise:\n",
53-
"1. Write a new version of your steepest descent code that rescales the step to be proportional to the local force, with $\\alpha = 0.01$. You should write this as a loop that iteratively updates the current position, i.e.\n",
87+
"The update equation then becomes:\n",
5488
"\n",
55-
"$$r_{i+1} = r_i + \\alpha F(r_i).$$\n",
89+
"$$r_{i+1} = r_i - \\alpha U^\\prime(r_i)$$\n",
5690
"\n",
57-
"By combining a suitable `if` statement with `break`, have your code stop and report the predicted equilibrium bond-length when $U < \\left|0.001\\right|$.\n"
91+
"or equivalently:\n",
92+
"\n",
93+
"$$r_{i+1} = r_i + \\alpha F(r_i)$$\n",
94+
"\n",
95+
"With this adaptive approach:\n",
96+
"- When far from the minimum, $|U^\\prime|$ is large → we take large steps\n",
97+
"- When close to the minimum, $|U^\\prime|$ is small → we take small steps\n",
98+
"\n",
99+
"The parameter $\\alpha$ controls the overall \"aggressiveness\" of the optimisation. Too large and we overshoot; too small and we converge slowly."
58100
]
59101
},
60102
{
61103
"cell_type": "markdown",
62-
"id": "182260eb-a629-414d-9b1a-a4587ed2326e",
104+
"id": "e84df851-bdbd-48d5-99a1-da0df6ea6034",
63105
"metadata": {},
64106
"source": [
65-
"2. How does changing $\\alpha$ affect your rate of convergence? Experiment with larger and smaller values of $\\alpha$.\n",
107+
"#### Exercise: Adaptive Step Size Gradient Descent\n",
108+
"\n",
109+
"1. Write a new version of your gradient descent code that rescales the step to be proportional to the local force. You should write this as a loop that iteratively updates the current position:\n",
110+
"\n",
111+
"$$r_{i+1} = r_i + \\alpha F(r_i)$$\n",
112+
"\n",
113+
"where $F(r_i) = -U^\\prime(r_i) = -k(r_i - r_0)$.\n",
114+
"\n",
115+
"Your code should:\n",
116+
"- Start from $r = 1.0$ Å\n",
117+
"- Use a `for` loop with a **maximum of 50 iterations**\n",
118+
"- Update the position at each iteration\n",
119+
"- Print the iteration number, position, gradient, and step size taken at each iteration\n",
120+
"- Stop early (using `break`) when $|U^\\prime(r)| < 0.001$\n",
121+
"- Report the final predicted equilibrium bond length and number of iterations required\n",
122+
"\n",
123+
"**Part 1:** Start with $\\alpha = 0.01$. Does it converge? How many iterations does it take? Compare this to the fixed step size results.\n",
124+
"\n",
125+
"**Part 2:** Try $\\alpha = 0.001$. What happens to the convergence rate?\n",
126+
"\n",
127+
"**Part 3:** Try $\\alpha = 0.1$. Does the algorithm remain stable?\n",
66128
"\n",
67-
"_Note_: You might want to set a maximum number of iterations, e.g., 30."
129+
"**Questions to consider:**\n",
130+
"- How does the step size change as you approach the minimum? Watch the \"Step (Å)\" column in your output.\n",
131+
"- Why doesn't the adaptive method oscillate like the fixed step size method with $\\Delta r = 0.1$ Å did?\n",
132+
"- Compare $\\alpha = 0.01$ (adaptive) with $\\Delta r = 0.01$ Å (fixed): which converges faster and why?\n",
133+
"- What happens when $\\alpha$ is too large? Can you explain this behaviour?\n",
134+
"- Based on your experiments, what seems to be a good choice of $\\alpha$ for this problem?"
68135
]
69136
},
70137
{
71138
"cell_type": "code",
72139
"execution_count": null,
73-
"id": "5c131eda-e2e5-47c4-bfe2-92f9331b2df6",
140+
"id": "2100979f-a7ce-4818-b2fa-43bf4848966b",
74141
"metadata": {},
75142
"outputs": [],
76143
"source": []
@@ -92,7 +159,7 @@
92159
"name": "python",
93160
"nbconvert_exporter": "python",
94161
"pygments_lexer": "ipython3",
95-
"version": "3.11.10"
162+
"version": "3.12.3"
96163
}
97164
},
98165
"nbformat": 4,

‎CH40208/geometry_optimisation/scipy_optimize_minimize.ipynb‎

Lines changed: 31 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -7,14 +7,14 @@
77
"source": [
88
"## Minimisation with `scipy.optimize.minimize`\n",
99
"\n",
10-
"`scipy.optimize.minimize()` is a powerful tool for numerical optimization that finds the minimum of any scalar function (a function that returns a single numerical value). \n",
10+
"`scipy.optimize.minimize()` is a powerful tool for numerical optimisation that finds the minimum of any scalar function (a function that returns a single numerical value). Unlike the methods we have implemented so far, `scipy.optimize.minimize()` uses more sophisticated algorithms that automatically handle many of the challenges we have encountered, such as choosing appropriate step sizes and dealing with difficult potential energy surfaces.\n",
1111
"\n",
1212
"The following code demonstrates how `minimize()` can be used to find the minimum of a harmonic potential."
1313
]
1414
},
1515
{
1616
"cell_type": "code",
17-
"execution_count": 7,
17+
"execution_count": 1,
1818
"id": "6e7de093-14a9-4de8-894e-33c44dc039e9",
1919
"metadata": {},
2020
"outputs": [
@@ -33,7 +33,7 @@
3333
" njev: 3"
3434
]
3535
},
36-
"execution_count": 7,
36+
"execution_count": 1,
3737
"metadata": {},
3838
"output_type": "execute_result"
3939
}
@@ -42,22 +42,22 @@
4242
"from scipy.optimize import minimize\n",
4343
"\n",
4444
"def harmonic_potential(r, k, r0):\n",
45-
" \"\"\"Calculate the energy in eV for a harmonic potential U(r) = k/2 * (r - r0)^2\n",
45+
" \"\"\"Calculate the energy in eV for a harmonic potential U(r) = k/2 * (r - r0)^2.\n",
4646
" \n",
47-
" Parameters:\n",
48-
" r (float): Distance parameter to optimize\n",
49-
" k (float): Spring constant in eV/Ų\n",
50-
" r0 (float): Equilibrium distance in Å\n",
47+
" Args:\n",
48+
" r (float): Distance parameter to optimise.\n",
49+
" k (float): Spring constant in eV/Ų.\n",
50+
" r0 (float): Equilibrium distance in Å.\n",
5151
" \n",
5252
" Returns:\n",
53-
" float: Potential energy in eV\n",
53+
" float: Potential energy in eV.\n",
5454
" \"\"\"\n",
5555
" return 0.5 * k * (r - r0)**2\n",
5656
"\n",
57-
"# initial guess and parameters\n",
58-
"x0 = 1.0 # starting guess for r\n",
59-
"k = 36.0 # spring constant\n",
60-
"r0 = 0.74 # equilibrium distance\n",
57+
"# Initial guess and parameters\n",
58+
"x0 = 1.0 # Starting guess for r\n",
59+
"k = 36.0 # Spring constant\n",
60+
"r0 = 0.74 # Equilibrium distance\n",
6161
"\n",
6262
"result = minimize(harmonic_potential, x0=x0, args=(k, r0))\n",
6363
"\n",
@@ -72,48 +72,44 @@
7272
"Let us break this down step by step:\n",
7373
"\n",
7474
"First, we import the `minimize()` function from `scipy.optimize`.\n",
75-
"\n",
7675
"```python\n",
7776
"from scipy.optimize import minimize\n",
7877
"```\n",
7978
"\n",
8079
"Next, we define a function for our harmonic potential:\n",
81-
"\n",
8280
"```python\n",
8381
"def harmonic_potential(r, k, r0):\n",
84-
" \"\"\"Calculate the energy in eV for a harmonic potential U(r) = k/2 * (r - r0)^2\n",
82+
" \"\"\"Calculate the energy in eV for a harmonic potential U(r) = k/2 * (r - r0)^2.\n",
8583
" \n",
86-
" Parameters:\n",
87-
" r (float): Distance parameter to optimize\n",
88-
" k (float): Spring constant in eV/Ų\n",
89-
" r0 (float): Equilibrium distance in Å\n",
84+
" Args:\n",
85+
" r (float): Distance parameter to optimise.\n",
86+
" k (float): Spring constant in eV/Ų.\n",
87+
" r0 (float): Equilibrium distance in Å.\n",
9088
" \n",
9189
" Returns:\n",
92-
" float: Potential energy in eV\n",
90+
" float: Potential energy in eV.\n",
9391
" \"\"\"\n",
9492
" return 0.5 * k * (r - r0)**2\n",
9593
"```\n",
9694
"\n",
97-
"Note that we have ordered the three arguments to this function so that the parameter we want to minimise with respect to, `r`, is the first parameter. The second and third parameters, `k` and `r0` will be fixed for a given potential energy surface.\n",
98-
"\n",
99-
"Then, we define variables to store our initial guess for the equilibrium bond length and for the parameters that define our harmonic potential.\n",
95+
"Note that we have ordered the three arguments to this function so that the parameter we want to minimise with respect to, `r`, is the first parameter. The second and third parameters, `k` and `r0`, will be fixed for a given potential energy surface.\n",
10096
"\n",
97+
"Then, we define variables to store our initial guess for the equilibrium bond length and the parameters that define our harmonic potential.\n",
10198
"```python\n",
102-
"# initial guess and parameters\n",
103-
"x0 = 1.0 # starting guess for r\n",
104-
"k = 36.0 # spring constant\n",
105-
"r0 = 0.74 # equilibrium distance\n",
99+
"# Initial guess and parameters\n",
100+
"x0 = 1.0 # Starting guess for r\n",
101+
"k = 36.0 # Spring constant\n",
102+
"r0 = 0.74 # Equilibrium distance\n",
106103
"```\n",
107104
"\n",
108105
"Now we call the `minimize()` function:\n",
109-
"\n",
110106
"```python\n",
111107
"result = minimize(harmonic_potential, x0=x0, args=(k, r0))\n",
112108
"```\n",
113109
"\n",
114-
"We pass in three arguments.\n",
115-
"- The **first** argument is the name of the function that we want to minimise. In this case our function is called `harmonic_potential`\n",
116-
"- The **second** argument is our starting guess for the parameter we are minimising with respect to, so in this example, we pass our guess for our equilibrium bond length of 1 &angst;.\n",
110+
"We pass in three arguments:\n",
111+
"- The **first** argument is the name of the function that we want to minimise. In this case our function is called `harmonic_potential`.\n",
112+
"- The **second** argument is our starting guess for the parameter we are minimising with respect to, so in this example, we pass our guess for our equilibrium bond length of 1 Å.\n",
117113
"- The **third** argument is a **tuple** of parameters that will also get passed to the function we are minimising, that are fixed during the minimisation. In this case, we need to define the shape of our harmonic potential energy surface by providing values for `k` and `r0`. Because our `harmonic_potential()` function takes two additional parameters, the tuple we pass to `minimize()` needs to contain two parameters. The order of these parameters also needs to match the order of the corresponding arguments in the function being minimised; so our tuple is structured as `(k, r0)`."
118114
]
119115
},
@@ -171,9 +167,9 @@
171167
"\n",
172168
"The potential energy associated with changes in the Cl-C-C-Cl dihedral angle, $\\theta$ in 1,2-dichloroethane can be modelled as\n",
173169
"\n",
174-
"$$U = \\frac{1}{2}\\left\\{A_1\\left[1+\\cos\\theta\\right] + A_2\\left[1+\\cos2\\theta\\right] A_3\\left[1+\\cos3\\theta\\right]\\right\\}$$\n",
170+
"$$U = \\frac{1}{2}\\left\\{A_1\\left[1+\\cos\\theta\\right] + A_2\\left[1+\\cos2\\theta\\right]+ A_3\\left[1+\\cos3\\theta\\right]\\right\\}$$\n",
175171
"\n",
176-
"with $A_1$ = 55.229 kJ mol<sup>&minus;1</sup>, $A_2$ = 3.3.472 kJ mol<sup>&minus;1</sup>, and $A_3$ = -58.576 kJ mol<sup>&minus;1</sup>.\n",
172+
"with $A_1$ = 55.229 kJ mol<sup>&minus;1</sup>, $A_2$ = 3.3472 kJ mol<sup>&minus;1</sup>, and $A_3$ = -58.576 kJ mol<sup>&minus;1</sup>.\n",
177173
"\n",
178174
"- Plot this function for $-\\pi \\leq \\theta \\leq \\pi$.\n",
179175
"\n",
@@ -205,7 +201,7 @@
205201
"name": "python",
206202
"nbconvert_exporter": "python",
207203
"pygments_lexer": "ipython3",
208-
"version": "3.11.10"
204+
"version": "3.12.3"
209205
}
210206
},
211207
"nbformat": 4,

0 commit comments

Comments
 (0)