# CHEM 1000 - Spring 2022
Prof. Geoffrey Hutchison, University of Pittsburgh

## 6 Optimizing Functions (Examples)

Chapter 6 in [*Mathematical Methods for Chemists*](http://sites.bu.edu/straub/mathematical-methods-for-molecular-science/)

By the end of this session, you should be able to:
- Understand general approaches to optimize functions either for minima or maxima (i.e. "extrema")
- Understand how to use derivatives in multiple dimensions to categorize extrema
- Using `scipy.optimize` to do numeric optimization of complicated functions
  - (We'll have more examples for both algebra/calculus and numerical optimization in recitation)

### Quick Summary

In chemistry and physics, we often want to determine the maximum or minimum value of a function of one or many variables. Examples include characterizing the minima, maxima, and saddle points on a potential energy surface.

Optimizing functions in one dimension is pretty easy, if sometimes tedious.
1. Take derivatives and find where the first derivative is zero
2. Look at the second derivatives, to categorize as a minima / maxima / inflection point
3. Then compare values of the function at those points to see if it's a local minima / max or the global minima / max.

### Many Variables

Not surprisingly, we can use a similar technique in multiple dimensions.

If we have a function $f(x,y)$ in two dimensions, then to have an extrema:

$$
\frac{\partial f}{\partial x}=0 \quad \frac{\partial f}{\partial y}=0
$$

In other words, we need to see the partial derivative with respect to **both / all** variables be zero.

We can then categorize the type of minima / maxima with the [Hessian]. (Later, we will see that this is the *determinant* of the Hessian matrix, for when we have more than 2 variables.)

$$
D=\left(\frac{\partial^{2} f}{\partial x^{2}}\right)\left(\frac{\partial^{2} f}{\partial y^{2}}\right)-\left(\frac{\partial^{2} f}{\partial x \partial y}\right)\left(\frac{\partial^{2} f}{\partial y \partial x}\right)
$$


$$
\left.\frac{\partial^{2} f}{\partial x^{2}}\right|_{\left(x^{*}, y^{*}\right)}=\left\{\begin{array}{ll}
<0 & D>0 & \text { maximum } \\
>0 & D>0 & \text { minimum } \\
 & D < 0 & \text { saddle-point }
\end{array}\right.
$$

### Example:

Let's try the so-called "butterfly potential"

$$
V(x, y)=\left((x-y)^{2}-1\right)^{2}+10 x^{2} y^{2}
$$




```python
from sympy import init_session
init_session()
```


```python
V = ((x-y)**2-1)**2 + 10*x**2*y**2
V
```


```python
# okay, let's see where the partial derivative with x is zero...
diff(V, x)
```


```python
# now let's do y
diff(V, y)
```

Okay, we can't do this so easily - we have two equations and two unknowns.

Fortunately, sympy offers some nice features:
- Eq(*expression*, value) sets up an equation like $x = 0$
  - We can set our derivatives as equal to zero
- solve() can take multiple equations and multiple unknowns


```python
eq1 = Eq(diff(V, x), 0) # partial derivative with x == 0
eq2 = Eq(diff(V, y), 0) # partial derivative with y == 0
solve((eq1,eq2), (x, y)) # solve two equations for two unknowns
```


```python
# here's the two-variable Hessian test...
D = diff(V, x, 2)*diff(V, y, 2) - diff(V, x, y)*diff(V, y, x)
D # print it
```

Okay, I like when the computer does math for me...


```python
# we can use the .subs() method to substitute values (-1, 0)
D.subs([ (x, -1), (y, 0) ])
```


```python
# check (0,0)
D.subs([ (x, 0), (y, 0) ])
```

The origin is inconclusive - because it's something like a saddle point, but there are four wells around it.


```python
# now check (+1, 0)
D.subs([ (x, 1), (y, 0) ])
```

### Water

We looked at a one-dimensional Lennard-Jones potential.

Here, I've calculated the O-H bond stretch in water, using the RI-MP2 method and the cc-pVTZ basis set. (Don't worry, you don't need to know what that means, but we're aiming for a relatively accurate quantum chemical method.)



Most bonds fit well to a [Morse Potential](https://en.wikipedia.org/wiki/Morse_potential) like this. There are reasons to use other functional forms (e.g., speed), but this relatively simple form captures a lot of how real covalent bonds behave.

I also calculated the H-O-H bond angle bend using the same method. Fortunately, this (like many angle bends) fit well to a quadratic function:
          


The total potential energy will be:

$$
V = V_{bond} + V_{angle}
$$

In principle, there are "cross-terms" (e.g., stretch-bend) involving both $r$ and $\theta$ but it's always good to start with a simple model and get more complex once we understand that.


```python
# We'll declare new symbols for now
r, theta = symbols("r theta")
bond = 265.79*(1 - exp((-2.2534)*(r - 0.96242)))**2
angle = 167.16 - 3.217*theta + 0.01548*theta**2
V = bond + angle
```


```python
# let's solve for the derivative w.r.t. theta equal to zero
solve( diff(angle, theta) )
```


```python
diff(bond, r)
```

Oof, that's a hard one to solve. We know that the minimum value is around 0.96Å, so we can try this numerically:


```python
lengths = [0.93, 0.94, 0.95, 0.96, 0.97, 0.98]
for l in lengths:
    print(l, round(bond.subs(r, l), 4), round(dr.subs(r, l), 4))
```

Clearly 0.96Å is very close. Notice that the partial derivative goes from negative to positive between 0.96Å and 0.97Å.

We can interpolate:


```python
for digit in range(0, 10):
    l = 0.96 + digit/1000
    print(l, round(bond.subs(r, l), 4), round(dr.subs(r, l), 4))
```

I can go further, but ~0.9625Å seems like a fairly accurate estimate of the optimal bond length.

I picked this example because sometimes we want to use several methods (e.g., both calculus and numerical methods) to solve problems.

Now to be completely honest, a better description of the potential energy function:

$$
V(r, \theta) = V_{bond}(r) + V_{angle}(\theta) + V_{stretch-bend}(r, \theta)
$$

Usually those stretch-bend terms look something like:

$$
V_{\mathrm{str} / \mathrm{bend}}=k^{A B C}\left(\theta^{A B C}-\theta_{0}^{A B C}\right)\left[\left(R^{A B}-R_{0}^{A B}\right)+\left(R^{B C}-R_{0}^{B C}\right)\right]
$$

-------
This notebook is from Prof. Geoffrey Hutchison, University of Pittsburgh
https://github.com/ghutchis/chem1000

<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"></a>

134a8e2029bc0f03e38d94325c75f773cffe8d72
ipynb
Jupyter Notebook
lectures/06a-optimizing-examples1.ipynb
ghutchis/chem1000
07a7eac20cc04ee9a1bdb98339fbd5653a02a38d
CC-BY-4.0
2020-06-23T18:44:37.000Z
2022-03-14T10:13:05.000Z
lectures/06a-optimizing-examples1.ipynb
ghutchis/chem1000
07a7eac20cc04ee9a1bdb98339fbd5653a02a38d
CC-BY-4.0
lectures/06a-optimizing-examples1.ipynb
ghutchis/chem1000
07a7eac20cc04ee9a1bdb98339fbd5653a02a38d
CC-BY-4.0
2021-07-29T10:45:23.000Z
2021-10-16T09:51:00.000Z
Qwen/Qwen-72B
1. YES
2. YES