Draconian_Fire
Draconian_Fire

Reputation: 3

How to plot the Density of states using histogram with a curve that follows the midpoints of each peak in python?

image

I want to plot the density of states as shown in the figure above (bottom right) where the general shape of the histogram is shown in the curve. This curve follows the general shape formed by the peaks' midpoint. Unfortunately, I can only plot the DOS like the top right plot where the curve looks like the histogram itself. Is there a way to fit the curve with points corresponding to the midpoint of histogram's peaks?

The image is from these lecture notes. Eigenstates: [ 0.81 0.81 0.81 0.75 0.75 0.75 0.54 0.54 0.54 0.47 0.47 0.47 0.33 0.33 0.75 0.81 0.27 0.27 0.33 0.15 0.15 0.27 0.54 0.15 0.47 0.33 -0.06 -0.06 -0.06 0.27 0.15 -0.08 -0.06 -0.2 -0.2 -0.2 -1.03 -0.08 -0.08 -0.27 -0.2 -1.03 -1.03 -1.03 -0.27 -0.27 -0.93 -0.93 -0.36 -0.36 -0.4 -0.4 -0.4 -0.27 -0.93 -0.36 -0.36 -0.4 -0.93 -0.93 -1.03 -1.03 -0.93 -1.03 -1.03 -1.03 -1.03 -0.93 -0.08 -0.93 -1.03 -1.03 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75 -0.75]

My code:

bins = 100

plt.figure(figsize=[8,5])
plt.rcParams['axes.titlesize'] = 17
plt.rcParams.update({'font.size': 15})
dN = plt.hist(Eig_nn, bins)
plt.xlabel('Eigenvalues')
plt.ylabel('N')
plt.title('All eigenvalues')
plt.show()

# DOS

dE = (max(Eig_nn) - min(Eig_nn)) / bins
DOS = dN[0]/dE
E = np.linspace(-6,4,bins)

plt.figure(figsize=[8,5])
plt.rcParams['axes.titlesize'] = 17
plt.rcParams.update({'font.size': 15})
plt.plot(E, DOS)
plt.xlabel('E')
plt.ylabel('DOS')
plt.title('All eigenvalues')
plt.show()

My current attempts

Upvotes: 0

Views: 1129

Answers (1)

JohanC
JohanC

Reputation: 80329

I'm unsure how the demo plot is generated. Maybe it is just connecting the tops of the histogram with a smoothed curve?

It looks a bit like a kde. Such a curve can be generated by scipy's gaussian_kde. Its first parameter are the datapoints, and the second determines the width of each gaussian. The optimal width depends on the data and how you want to interpret them. A value around 0.1 seems adequate for this example. It's unclear how to best align the histogram with the kde, so the code below plots the kde using a secondary y-axis.

Some code to experiment with:

from matplotlib import pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde

Eig_nn = np.array([...])

bins = 100
plt.rcParams.update({'font.size': 15})
fig, ax1 = plt.subplots(figsize=[8,5])
bin_values, bin_lims, _patches  = plt.hist(Eig_nn, bins, color='dodgerblue')

x = np.linspace(bin_lims[0]-0.2, bin_lims[-1]+0.2, 500) # create an x-axis
kde = gaussian_kde(Eig_nn, 0.1)
ax2 = ax1.twinx()  # create a secondary axis for the kde
ax2.plot(x, kde(x), color='crimson')

ax1.set_xlim(x[0], x[-1])  # set strict limits
ax1.set_xlabel('Eigenvalues')
ax1.set_ylabel('N')
ax2.set_ylabel('kde')
ax2.set_ylim(ymin=0) # put the zero of the secondary y-axis at the bottom
plt.title('All eigenvalues', fontsize=20)
plt.show()

sample plot

Upvotes: 1

Related Questions