```
from sympy import *
import numpy as np
import math
```

<h2>Strutinsky Energy Averaging Method</h2>

Define a callable class for computing the smooth part of the density of states $\tilde{g}_m(E)$ with curvature
correction of order $2M$.


```
class SmoothDOS(object):

    def __init__(self, energies):
        self.energies = energies

    def gaussian_smoothing_func(self, M, x):
        return (mpmath.laguerre(M,0.5,x**2)*exp(-x**2)/sqrt(pi)).evalf()

    def __call__(self, e, M=3, gamma=1.0):
        return sum(self.gaussian_smoothing_func(M, (e-ei)/gamma)/gamma for ei in self.energies)
```


```
def avgf(M, x):
    return (mpmath.laguerre(M,0.5,x**2)*exp(-x**2)/sqrt(pi)).evalf()
```


```
def smooth_dos(gamma, M, e, energies):
    return sum(avgf(M, (e-ei)/gamma)/gamma for ei in energies)
```

<h2>1D Simple Harmonic Oscillator DOS</h2>


```
def sho_smooth_dos(e):
    """Compute the exact smooth DOS."""
    return 1
```


```
def sho_spectrum(nmax):
    """Compute the first nmax energies of the 1D SHO."""
    return [n+0.5 for n in range(nmax)]
```


```
sho_evalues = np.linspace(0.0,20,100)
```


```
sho_dos = SmoothDOS(sho_spectrum(30))
```


```
sho_exact_dos = [sho_smooth_dos(e) for e in sho_evalues]
```


```
sho_approx_dos = [sho_dos(e,M=3,gamma=10.0) for e in sho_evalues]
```


```
plot(sho_evalues, sho_exact_dos, label="exact")
plot(sho_evalues, sho_approx_dos, label="approx")
title("Smooth part of the DOS for the 1D SHO")
xlabel("Energy"); ylabel("$g(E)$")
legend(loc=4)
```

Note how the exact smoothed DOS and the Strutinsky approximation agree once the energy is away from
0. In general, these are edge effects and are also present at the right limit if you don't use enough energy
levels.  Here we have used 30, so the right side doesn't show these artifacts.

<h2>1D Particle in a BOX (PIAB) DOS</h2>


```
def piab_smooth_dos(e):
    """Compute the exact smooth DOS."""
    return 1.0/(2.0*math.sqrt(e*math.pi**2/2))
```


```
def piab_spectrum(nmax):
    """Compute the first nmax energies of the 1D PIAB."""
    return [0.5*math.pi**2*(n+1)**2 for n in range(nmax)]
```


```
piab_evalues = linspace(1.0,1000.0,100)
```


```
piab_dos = SmoothDOS(piab_spectrum(20))
```


```
piab_exact_dos = [piab_smooth_dos(e) for e in piab_evalues]
```


```
piab_approx_dos = [piab_dos(e, M=4, gamma=150.0) for e in piab_evalues]
```


```
plot(piab_evalues, piab_exact_dos, label="exact")
plot(piab_evalues, piab_approx_dos, label="approx")
title("Smooth part of the DOS for the 1D PIAB")
xlabel("Energy"); ylabel("$g(E)$")
legend(loc=1)
```


```

```

f62e593d182a5923e1884f315a5d4e3347c93448
ipynb
Jupyter Notebook
docs/examples/notebooks/smooth_dos.ipynb
tinyclues/ipython
71e32606b0242772b81c9be0d40751ba47d95f2c
BSD-3-Clause-Clear
2016-05-26T10:57:18.000Z
2016-05-26T10:57:18.000Z
docs/examples/notebooks/smooth_dos.ipynb
adgaudio/ipython
a924f50c0f7b84127391f1c396326258c2b303e2
BSD-3-Clause-Clear
docs/examples/notebooks/smooth_dos.ipynb
adgaudio/ipython
a924f50c0f7b84127391f1c396326258c2b303e2
BSD-3-Clause-Clear
Qwen/Qwen-72B
1. YES
2. YES