smooth_dos.ipynb
201 lines
| 49.5 KiB
| text/plain
|
TextLexer
In [1]:
from sympy import *
import numpy as np
import math
Strutinsky Energy Averaging Method
Define a callable class for computing the smooth part of the density of states $\tilde{g}_m(E)$ with curvature correction of order $2M$.
In [2]:
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)
In [3]:
def avgf(M, x):
return (mpmath.laguerre(M,0.5,x**2)*exp(-x**2)/sqrt(pi)).evalf()
In [4]:
def smooth_dos(gamma, M, e, energies):
return sum(avgf(M, (e-ei)/gamma)/gamma for ei in energies)
1D Simple Harmonic Oscillator DOS
In [5]:
def sho_smooth_dos(e):
"""Compute the exact smooth DOS."""
return 1
In [6]:
def sho_spectrum(nmax):
"""Compute the first nmax energies of the 1D SHO."""
return [n+0.5 for n in range(nmax)]
In [7]:
sho_evalues = np.linspace(0.0,20,100)
In [8]:
sho_dos = SmoothDOS(sho_spectrum(30))
In [10]:
sho_exact_dos = [sho_smooth_dos(e) for e in sho_evalues]
In [23]:
sho_approx_dos = [sho_dos(e,M=3,gamma=10.0) for e in sho_evalues]
In [24]:
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)
Out[24]:
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.
1D Particle in a BOX (PIAB) DOS
In [16]:
def piab_smooth_dos(e):
"""Compute the exact smooth DOS."""
return 1.0/(2.0*math.sqrt(e*math.pi**2/2))
In [17]:
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)]
In [18]:
piab_evalues = linspace(1.0,1000.0,100)
In [19]:
piab_dos = SmoothDOS(piab_spectrum(20))
In [20]:
piab_exact_dos = [piab_smooth_dos(e) for e in piab_evalues]
In [21]:
piab_approx_dos = [piab_dos(e, M=4, gamma=150.0) for e in piab_evalues]
In [22]:
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)
Out[22]:
In [ ]: