# Phonon density of states as a function of frequency vs energy [closed]

Lattice dynamics codes like Phonopy output density of states, $$g_{(\omega)}$$, as a function of frequency. Since $$\int g_{(\omega)}d\omega$$ is the number of states between $$\omega$$ and $$d\omega$$, can the same be extended to $$\epsilon$$, where $$\epsilon = 2\pi\omega$$?

To put it briefly: Are $$g_{(\omega_i)}$$ and $$g_{(\epsilon_i)}$$ the same for $$\epsilon_i = 2\pi\omega_i$$?

Intuitively, a phonon vibrating at $$\omega$$ will have an energy of $$\epsilon = 2\pi\omega$$, so the number of states between $$\omega$$ and $$d\omega$$ should be the same as that between $$\epsilon$$ and $$d\epsilon$$, where $$\epsilon = 2\pi\omega$$ and $$d\epsilon = 2\pi d\omega$$.

One way I see this intuition being wrong is $$d\epsilon \ne 2\pi d\omega$$ but $$\epsilon = 2\pi\omega$$.

All of this sprang up into my head while I was trying to calculate harmonic and quasi-harmonic properties from the Phonopy output.

I'm confused (which might be apparent from the overuse of the equation "$$\epsilon = 2\pi\omega$$" in the post). Will be glad to get some insight. Also, please feel free to rewrite the post in a way that those redundant equations can be removed.

Edit: I understand $$d\epsilon \ne 2\pi d\omega$$ sounds silly. To make it sound less silly, here are phonopy output files for Ni$$_3$$Al. Can you try and find the zero-point energy? It doesn't match with what the code outputs. Here is what I tried:

import numpy as np
import math

def main():

kB = 1.380649 * (10**-23) #Boltzmann Cconstant (J/K)

hbar = 1.054571817 * (10**-34) #Planck's constant (hbar in J.s)

NA = 6.02214 * (10**23) #Avogadro's number (per mole)

num_atoms = 4 #Number of atoms in the primitive cell that Phonopy finds, equal to natoms in thermal_properties.yaml

with open('freq.dat') as f: #freq.dat = cat mesh.yaml | grep frequency

lines = f.readlines() #Store each line in the list lines

A = np.empty(len(lines), dtype = float) #Create an empty array A to store the contents of the file

for i in range(len(lines)):
A[i] = float(lines[i].split()[1])

for i in range(len(A)):
A[i] = A[i] * 2 * math.pi * (10**12) #Convert frequency (THz) to angular frequency (radians/s)
A[i] = (A[i] * hbar) / 1000 #Convert angular frequency (radians/s) to energy (KJ)

#Calculate the Zero Point Energy (ZPE) (KJ/mol)
ZPE = 0

for i in range(len(A)):
ZPE = ZPE + (A[i] / 2) #(A[i] / 2) = hbar x omega / 2
ZPE = (ZPE * NA) / (num_atoms)
print("ZPE = {0}".format(ZPE))

if __name__ == '__main__':
main()

• I don't know this package, so I can't answer your actual question, but I can confirm that there is no theoretical obstacle that would prevent you from using a $g(\epsilon)$ instead of $g(\omega)$. Both are commonly done. It may be the case that there is something in the Phonopy package that does not allow for this. One suggestion: have you double check that you changed your limits of integration accordingly? Nov 16 '20 at 16:48
• @taciteloquence, yes, I changed integration limits accordingly. In fact, it's something about the phonopy output that I seem to misunderstand. The calculation for the zero-point energy, which doesn't require integration since it is the sum of $\bar{h}\omega/2$ over all frequencies of vibration at each q-point for all atoms, also doesn't match with the output from the code. Nov 17 '20 at 18:13
• When I click on your link: gofile.io/d/mY8qtd ... it says "file not found". Can you upload it here in a folder called 3744 (since that's the number in the URL for this question)? Aug 14 at 17:47