# Variational Principle  using Symbolic Mathematics in Python 

## 1. Introduction

The variational principle tells us that we can use a trial wavefunction to solve the Schrodinger equation using the following theorem:

$${{\int {{\Psi ^*}\hat H{\rm{ }}\Psi } d\tau } \over {\int {{\Psi ^*}\Psi } d\tau }} \ge {E_0}$$

We will use Sympy to solve the particle in a box problem by guessing a trial wavefunction using variational principle


```python
import sympy as sym
```

This exercise is a bit more self-guided than the other notebooks we have done. One of the most useful things you can do is **open last week's notebook to remember the commands in sympy**. Also, remember that google is your friend:

1. [Sympy tutorial](https://docs.sympy.org/latest/tutorial/index.html)
2. [Stack Overflow](https://stackoverflow.com/search?q=sympy+)
3. [Stack Exchange](https://stackexchange.com/)

## 2. Particle in a box

The wave function that we pick for a particle in a box needs to have the following properties

1. single valued
1. normalizable
1. function and its first derivative are continuous 
1. boundary condition that the wave function goes to zero at the ends of the box



Particle in a box: a is a classical particle, red is real part, blue is imaginary part.

This particle only expericnes kinetic energy between the box, so the Hamiltonian for this system is

$$\hat H = {{ - {\hbar ^2}} \over {2m}}{{{d^2}} \over {d{x^2}}} + \left\{ {\matrix{{V(x) = 0} & {0 < x < a}  \cr {V(x) = \infty } & {x < 0\text{ }{\rm{ or}}\;x > a}  \cr } } \right.$$

For our purposes, that means we can consider the Hamiltonian to be 

$$\hat H = {{ - {\hbar ^2}} \over {2m}}{{{d^2}} \over {d{x^2}}}$$

as long as we keep the limits of integration to be $(0,a)$

### 2.1 Trial Wave function

Although the particle in box has a well known solution

[https://en.wikipedia.org/wiki/Particle_in_a_box](https://en.wikipedia.org/wiki/Particle_in_a_box)

(or check your favorite pchem book)

We are going to guess a trial wave function:

$$\Phi (x) = \left( {{x \over a} - {{{x^3}} \over a}} \right) + \alpha \left( {{{{x^5}} \over {{a^5}}} - {1 \over 2}\left( {{{{x^7}} \over {{a^7}}} + {{{x^7}} \over {{a^7}}}} \right)} \right)$$

### 2.2 Exercise: Variational Theorem

We are going to follow the following plan:

1. Solve for the energy of the trial wave function above  

$${E_{trial}} = {{\int\limits_0^a {\Phi (x){{ - {\hbar ^2}} \over {2m}}{{{d^2}} \over {d{x^2}}}\Phi (x)dx} } \over {\int\limits_0^a {\Phi {{(x)}^2}dx} }}$$

Your answer will be a function of $ m,a,\text{and } \alpha$ We will use $\alpha$ as the parameter we vary to minimize the energy and make a new trial wave function.

2. Minimize the trial energy 
We will use a first derivative of the trial energy $${d \over {d\alpha }}{E_{trial}}(\alpha )$$ to find the value of $\alpha$ that gives you the lowest energy


3. Plot your new wavefunction compared to the ground state particle in a box: $${\psi _{true}}(x) = {\left( {{2 \over a}} \right)^{1/2}}\sin {{n\pi x} \over a}$$ Plot as a function of $x/a$ from $0$ to $1$. Assuming this has $m=m_e$, and $a=a_0$ use atomic (theorist) units to plot the function.

4. Compare your trial energy to the actual energy (using atomic units)
$${E_{true}}(n = 1) = {{{\hbar ^2}{\pi ^2}} \over {2m{a^2}}}$$



```python
import sympy as sym
a,x,m,alpha = sym.symbols('a,x,m,alpha')
sym.Rational(1/2)
sym.init_printing()
```


```python
phi = ((x/a)-(x**3/a**3)+ alpha*((x**5/a**5)-1/2*(x**9/a**9+x**9/a**9)))
num1 = (phi*sym.Rational(-1/2))
num2 = sym.diff(sym.diff(phi,x),x)
den1 = (phi**2)
```


```python
Numerator = sym.integrate((num1*num2),(x,0,a))
```


```python
Denominator = sym.integrate((den1),(x,0,a))
```

Your descriptions/explanations here


```python
Numerator/Denominator
```


```python
E_trial = (Numerator/Denominator).subs(a,1)
E_trial
```


```python
sym.solveset(sym.diff(E_trial,alpha))
```


```python
E_trial.subs(alpha, -0.257747652303185)
```


```python
E_trial.subs(alpha, -4.73345207773701)
```


```python
Newphi = phi.subs(alpha, -0.257747652303185)
Newphi
```


```python
Newphi2 = Newphi.subs(a,1)
Newphi2
```


```python
sym.plot(Newphi2,(x,0,1))
```


```python
sym.plot((2**1/2)*sym.sin(sym.pi*x),(x,0,1))
```


```python
sym.plot((2**1/2)*sym.sin(sym.pi*x),Newphi2,(x,0,1))
```


```python
E_true = sym.pi**2/(2*m)
E_true
```

### 2.3 Exercise: New trial wavefunction

Determine the minimum energy of the particle in a box using a new trial wavefunction $$x^\alpha(x-a)^\alpha$$

1. Find the minimum energy, $E_{trial}$
2. Plot the new trial wavefunction and compare it to the true solution and the wavefunction you found above
3. Compare you new energy to the trial energy you found above
4. Which wavefunction is better? How do you know?


```python
#Plug in 1 for a, after you solve for values, first find the first derivative of part 1 function in terms of alpha 
#After you solve both plug them into E trial and see which gives you the lower E value
```


```python
phhi = (x**alpha)*(x-a)**alpha
phhi
```


```python
num11 = (phhi*sym.Rational(1/2))
num21 = sym.diff(sym.diff(phhi,x),x)
den11 = (phhi**2)
Numerator2 = sym.integrate((num11*num21),(x,0,a))
Numerator2

```


```python
Denominator2 = sym.integrate((den11),(x,0,a))
Denominator2
```




$\displaystyle \frac{a a^{4 \alpha} e^{2 i \pi \alpha} \Gamma\left(2 \alpha + 1\right) {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)}}{\Gamma\left(2 \alpha + 2\right)}$




```python
E_trial2 = (Numerator2/Denominator2).subs(a,1)
E_trial2
```




$\displaystyle \frac{\alpha e^{- 2 i \pi \alpha} \Gamma\left(2 \alpha + 2\right) \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx}{2 \Gamma\left(2 \alpha + 1\right) {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)}}$




```python
sym.solveset(sym.diff(E_trial2,alpha))
```




$\displaystyle \left\{\alpha \mid \alpha \in \mathbb{C} \wedge \left(- 2 \alpha {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)} \operatorname{polygamma}{\left(0,2 \alpha + 1 \right)} \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx + 2 \alpha {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)} \operatorname{polygamma}{\left(0,2 \alpha + 2 \right)} \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx + \alpha {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)} \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 x^{2} - 4 x + 2 \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right) \log{\left(x \right)} + 2 \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right) \log{\left(x - 1 \right)} + 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx - 2 i \pi \alpha {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)} \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx - \alpha \frac{d}{d \alpha} {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)} \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx + {{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)} \int\limits_{0}^{1} \frac{x^{2 \alpha} \left(x - 1\right)^{2 \alpha} \left(4 \alpha x^{2} - 4 \alpha x + \alpha - 2 x^{2} + 2 x - 1\right)}{x^{2} \left(x - 1\right)^{2}}\, dx\right) \Gamma\left(2 \left(\alpha + 1\right)\right) = 0 \right\} \setminus \left\{\alpha \mid \alpha \in \mathbb{C} \wedge e^{2 i \pi \alpha} \Gamma\left(2 \alpha + 1\right) {{{}_{2}F_{1}\left(\begin{matrix} - 2 \alpha, 2 \alpha + 1 \\ 2 \alpha + 2 \end{matrix}\middle| {1} \right)}}^{2} = 0 \right\}$




```python

```

Your descriptions/explanations here

### 2.4 Exercise: Design your own wavefunction!

**Now you get to make your own wavefunction!**

The only guidance I would give you is that it make sense mathematically and that it include $\alpha$ so that you can minimize the energy.

Remember that $a$ and $x$ are both length units, and that trigonometric, logarithmic, and exponential functions are all unitless  


Using your new wavefunction:

1. Find the minimum energy, $E_{trial}$
2. Plot the new trial wavefunction and compare it to the true solution and the wavefunction you found above
3. Compare you new energy to the trial energy you found above
4. Which wavefunction is better? How do you know?


```python
import sympy as sym
a,x,alpha = sym.symbols('a,x,alpha')
sym.Rational(1/2)
sym.init_printing()
```


```python
patsywavefunction = (x**2)*(x-a)**alpha
patsywavefunction
```


```python
num1 = (patsywavefunction*sym.Rational(-1/2))
num2 = sym.diff(sym.diff(patsywavefunction,x),x)
den1 = (patsywavefunction**2)
```


```python
value = (num1*num2)
```


```python
Numerator = sym.integrate((value),(x,0,a))
```


```python
Denominator = sym.integrate((den1),(x,0,a))
```


```python
Numerator/Denominator
```




$\displaystyle \left(\frac{1.42081447963801 \alpha^{2}}{a} - \frac{1.68831168831169 \alpha}{a} + \frac{0.2 \left(- 10.0 \alpha - 3.0\right)}{a} + \frac{0.111111111111111 \left(- 10.0 \alpha^{2} + 36.0 \alpha\right)}{a} + \frac{1.0}{a}\right) \left(\begin{cases} \frac{1}{- \log{\left(- a \right)} - \infty \operatorname{sign}{\left(a^{4} \right)} - \frac{25}{12}} & \text{for}\: \alpha = - \frac{5}{2} \\\frac{1}{- 4 a \log{\left(- a \right)} - \frac{4 a}{3} + \infty \operatorname{sign}{\left(a^{4} \right)}} & \text{for}\: \alpha = -2 \\\frac{1}{- 6 a^{2} \log{\left(- a \right)} - 3 a^{2} - \infty \operatorname{sign}{\left(a^{4} \right)}} & \text{for}\: \alpha = - \frac{3}{2} \\\frac{1}{- 4 a^{3} \log{\left(- a \right)} - 4 a^{3} + \infty \operatorname{sign}{\left(a^{4} \right)}} & \text{for}\: \alpha = -1 \\\frac{1}{- a^{4} \log{\left(- a \right)} - \infty \operatorname{sign}{\left(a^{4} \right)}} & \text{for}\: \alpha = - \frac{1}{2} \\\frac{\left(- a\right)^{- 2 \alpha} \left(8 \alpha^{5} + 60 \alpha^{4} + 170 \alpha^{3} + 225 \alpha^{2} + 137 \alpha + 30\right)}{6 a^{5}} & \text{otherwise} \end{cases}\right)$




```python
E_trial = (Numerator/Denominator).subs(a,1)
E_trial
```




$\displaystyle \left(0.309703368526898 \alpha^{2} + 0.311688311688312 \alpha + 0.4\right) \left(\begin{cases} \frac{1}{-\infty - i \pi} & \text{for}\: \alpha = - \frac{5}{2} \\\frac{1}{\infty - 4 i \pi} & \text{for}\: \alpha = -2 \\\frac{1}{-\infty - 6 i \pi} & \text{for}\: \alpha = - \frac{3}{2} \\\frac{1}{\infty - 4 i \pi} & \text{for}\: \alpha = -1 \\\frac{1}{-\infty - i \pi} & \text{for}\: \alpha = - \frac{1}{2} \\\frac{\left(-1\right)^{- 2 \alpha} \left(8 \alpha^{5} + 60 \alpha^{4} + 170 \alpha^{3} + 225 \alpha^{2} + 137 \alpha + 30\right)}{6} & \text{otherwise} \end{cases}\right)$




```python
sym.solveset(sym.diff(E_trial,alpha))
```


```python
sym.solveset(sym.diff(E_trial,alpha))
```


```python

```

Your descriptions/explanations here

# Reading Homework

Read the following sections in Kramer

- 4.2.3 Born-Oppenheimer approximation
- 4.3.2 Secular equation
- All of 4.5

For each subsection
- write down the subchapter name
- what was the most important idea
- draw an idea digram of the main idea

**Make sure to upload this to the assignment repository**

Example idea diagram:



```python

```

1d0b71d6bf7fb8580f088582b4aaa81134134594
ipynb
Jupyter Notebook
variational-principle.ipynb
sju-chem264-2019/10-3-19-lecture-PatsyGriffin
4724db01ebeda7cfcdcedbafc9bf3931c9190eb2
MIT
variational-principle.ipynb
sju-chem264-2019/10-3-19-lecture-PatsyGriffin
4724db01ebeda7cfcdcedbafc9bf3931c9190eb2
MIT
variational-principle.ipynb
sju-chem264-2019/10-3-19-lecture-PatsyGriffin
4724db01ebeda7cfcdcedbafc9bf3931c9190eb2
MIT
Qwen/Qwen-72B
1. YES
2. YES