```python
from sympy import *
from sympy.abc import r,x,y,z
from scipy.integrate import quad, nquad
import matplotlib.pyplot as plt
%matplotlib inline
init_printing()
```

# Energy of the Hydrogen Atom
The variational principle states a trial wavefunction will have an energy greater than or equal to the ground state energy.
$$\frac{\int \psi H \psi}{ \int \psi^2}  \ge E_0$$

First consider the hydogen atom. Let us use a trial wavefunction that is not the exact ground state.


```python
beta = Symbol('beta')
R_T = exp(-r - beta*r*r)
R_T
```

The Hamiltonian for this system is
$$-\frac{1}{2} \nabla^2 - \frac{1}{r}$$
The first term is the kinetic energy of the electron, and the second term is the Coulomb attraction between the electron and proton.

The first step is to compute the derivative of the trial wavefunction in spherical coordinates


```python
def del_spherical(e, r):
    """Compute Laplacian for expression e with respect to symbol r.
        Currently works only with radial dependence"""
    t1 = r*r*diff(e, r)
    t2 = diff(t1, r)/(r*r)
    return simplify(t2)
```


```python
del1 = del_spherical(R_T, r)

```

Construct $\psi H \psi$


```python
H = -1/S(2) * R_T * del1  - R_T*R_T/r
```


```python
simplify(H)
```

The integration occurs in 3D over the electron coordinates.  Because the integrand only has a dependence on $r$, it can be converted to spherical coordinates, and reduced to a 1D integral over $r$. (There should be an additional factor of $4 \pi$, but it will cancel since it occurs in the numerator and denominator)


```python
h1 = simplify(r*r*H)
```

Substitute a concrete value for $\beta$.


```python
h2 = h1.subs(beta, 0.1)
```

Perform the integral


```python
num = integrate(h2, (r, 0, oo)).evalf()
num
```

Also construct and integrate the denominator (the normalization).


```python
norm1 = r*r*R_T*R_T
norm2 = norm1.subs(beta, 0.1)
norm3 = simplify(norm2)
```


```python
denom = integrate(norm3, (r, 0, oo)).evalf()
simplify(denom).evalf()
```


```python
E = num/denom
E
```

And, as expected, energy is greater than the exact ground state energy of -0.5 Hartree.

## Find the minimum energy
Collect all the steps for computing the energy into a single function.   Even though this particular integral could be done symbolically, use numerical integration instead.


```python
def compute_energy(R_T, beta_val):
    """Energy given a value for beta"""
    
    # Normalization integrand (denominator)
    norm1 = r*r*R_T*R_T
    norm2 = norm1.subs(beta, beta_val)
    norm3 = simplify(norm2)

    # Integrand for the numerator
    del1 = del_spherical(R_T, r)
    # Construct psi * H * psi
    H = -1/S(2) * R_T * del1  - R_T*R_T/r
    h1 = simplify(r*r*H)
    h2 = h1.subs(beta, beta_val)
    
    lim = 20.0
    
    denom_func = lambdify([r], norm3)
    denom_res = quad(denom_func, 0.0, lim)
    
    num_func = lambdify([r], h2)
    num_res = quad(num_func, 0.0, lim)

    e = num_res[0]/denom_res[0]

    return e
   
```

Now the energy can be computed vs. $\beta$, and we can find the minimum energy.  In this case, the minimum occurs at $\beta = 0$, which we know is the exact wavefunction for the hydrogen atom.


```python
energies = []
betas = []
for i in range(10):
    beta_val = i*.01
    e = compute_energy(R_T, beta_val)
    betas.append(beta_val)
    energies.append(e)

plt.plot(betas, energies)
```


```python

```

2166ce9847e281a2a483c8efa4b2e05cb3914019
ipynb
Jupyter Notebook
Variational/Variational_Hydrogen.ipynb
QMCPACK/qmc_algorithms
015fd1973e94f98662149418adc6b06dcd78946d
MIT
2018-02-06T06:15:19.000Z
2019-11-26T23:54:53.000Z
Variational/Variational_Hydrogen.ipynb
chrinide/qmc_algorithms
015fd1973e94f98662149418adc6b06dcd78946d
MIT
2017-03-23T17:17:04.000Z
2017-03-23T17:17:04.000Z
Variational/Variational_Hydrogen.ipynb
chrinide/qmc_algorithms
015fd1973e94f98662149418adc6b06dcd78946d
MIT
2016-06-30T21:29:32.000Z
2019-10-22T16:10:03.000Z
Qwen/Qwen-72B
1. YES
2. YES