# Variational Monte Carlo: Harmonic Oscillator


```python
import numpy as np
import matplotlib.pyplot as plt

%config InlineBackend.figure_format = 'retina'
%matplotlib inline
try: plt.style.use('./notebook.mplstyle')
except: pass

red,blue,green = '#e85c47','#4173b2','#7dcca4'
```

## Hamiltonian

\begin{equation}
\hat{H} = -\frac{1}{2} \frac{d^2}{dx^2} + \frac{1}{2} x^2
\end{equation}

where $\hbar = \omega = m = 1$.

## Local Energy and Transition Probability

\begin{equation}
E_L^\alpha (x) = \alpha + x^2\left(\frac{1}{2} - 2\alpha^2\right)
\end{equation}

\begin{equation}
\frac{\pi(x^\prime)}{\pi(x)} = \mathrm{e}^{-2\alpha({x^\prime}^2 - x^2)}
\end{equation}


```python
def EL(x,α):
    return α + x**2*(0.5-2*α**2)

def transition_probability(x,x̄,α):
    return np.exp(-2*α*(x̄**2-x**2))
```


```python
def vmc(num_walkers,num_MC_steps,num_equil_steps,α,δ=1.0):
    
    # initilaize walkers
    walkers = -0.5 + np.random.rand(num_walkers)
    
    # initialize energy and number of accepted updates
    estimator = {'E':np.zeros(num_MC_steps-num_equil_steps)}
    num_accepted = 0
    
    for step in range(num_MC_steps):
        
        # generate new walker positions 
        new_walkers = np.random.normal(loc=walkers, scale=δ, size=num_walkers)
        
        # test new walkers
        for i in range(num_walkers):
            if np.random.random() < transition_probability(walkers[i],new_walkers[i],α):
                num_accepted += 1
                walkers[i] = new_walkers[i]
                
            # measure energy
            if step >= num_equil_steps:
                measure = step-num_equil_steps
                estimator['E'][measure] = EL(walkers[i],α)
                
    # output the acceptance ratio
    print('accept: %4.2f' % (num_accepted/(num_MC_steps*num_walkers)))
    
    return estimator
```

## Perform the VMC Simulation


```python
α = 0.45
num_walkers = 100
num_MC_steps = 20000
num_equil_steps = 5000

np.random.seed(1173)

estimator = vmc(num_walkers,num_MC_steps,num_equil_steps,α)
```

    accept: 0.62


## Compute the average energy and standard error


```python
from scipy.stats import sem
Ē,ΔĒ = np.average(estimator['E']),sem(estimator['E'])

print('Ē = %f ± %f' % (Ē,ΔĒ))
```

    Ē = 0.501036 ± 0.000566


## Due a brute-force minimizaiton search over the variational parameter $\alpha$


```python
Ēmin = []
ΔĒmin = []
α = np.array([0.45, 0.475, 0.5, 0.525, 0.55])
for cα in α: 
    estimator = vmc(num_walkers,num_MC_steps,num_equil_steps,cα)
    Ē,ΔĒ = np.average(estimator['E']),sem(estimator['E'])
    Ēmin.append(Ē)
    ΔĒmin.append(ΔĒ)
    print('%5.3f \t %7.5f ± %f' % (cα,Ē,ΔĒ))
```

    accept: 0.62
    0.450 	 0.50387 ± 0.000635
    accept: 0.62
    0.475 	 0.50012 ± 0.000291
    accept: 0.61
    0.500 	 0.50000 ± 0.000000
    accept: 0.60
    0.525 	 0.50151 ± 0.000280
    accept: 0.59
    0.550 	 0.50239 ± 0.000543


## Compare VMC with the exact variational energy

\begin{equation}
E_v = \frac{\alpha}{2} + \frac{1}{8\alpha}
\end{equation}


```python
cα = np.linspace(α[0],α[-1],1000)
plt.plot(cα,0.5*cα + 1/(8*cα), '-', linewidth=1, color=green, zorder=-10, 
         label=r'$\frac{\alpha}{2} + \frac{1}{8\alpha}$')
plt.errorbar(α,Ēmin,yerr=ΔĒmin, linestyle='None', marker='o', elinewidth=1.0, 
             markersize=6, markerfacecolor=blue, markeredgecolor=blue, ecolor=blue, label='VMC')
plt.xlabel(r'$\alpha$')
plt.ylabel('E');
plt.xlim(0.44,0.56)
plt.legend(loc='upper center')
```


```python

```

8e92c7671066b09d3d95e58a9965328e11f3a97c
ipynb
Jupyter Notebook
vmc_ho.ipynb
agdelma/qmc_ho
1a5875836b8ed69384c5af49d9649cf9945e827c
MIT
2018-01-11T22:11:04.000Z
2021-03-04T10:13:55.000Z
vmc_ho.ipynb
agdelma/qmc_ho
1a5875836b8ed69384c5af49d9649cf9945e827c
MIT
vmc_ho.ipynb
agdelma/qmc_ho
1a5875836b8ed69384c5af49d9649cf9945e827c
MIT
2019-07-21T08:55:14.000Z
2020-10-21T11:35:50.000Z
Qwen/Qwen-72B
1. YES
2. YES