10
$\begingroup$

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()
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$ Nov 16 '20 at 16:48
  • $\begingroup$ @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. $\endgroup$ Nov 17 '20 at 18:13
  • $\begingroup$ 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)? $\endgroup$ Aug 14 at 17:47