Mistakes in the doc of hdf5 files
Closed this issue · 13 comments
On the webpage about hdf5 files (How to read the results stored in hdf5 files), the physical unit of group velocity seems to be given wrong. And the unit "THz" here ought to be angular instead of ordinal. The same is true for gv_by_gv.
I used phono3py for graphene recently and found that I have to regard the group_velocity as angular to get the correct sound speed, just as below in python code:
sqrt(group_velocity[5,1,:]**2.sum())
# [5,1,:] ensures the point is near Gamma point(151*151 q-mesh) and in TA branch
# the code above gives 137.17 (THz*Ang), close to the value 12.9 km/s from current \
# research (https://doi.org/10.1016/j.carbon.2019.04.006) only if we consider it angularBesides, I can deduce the kappa_unit_conversion correctly only when I regard the THz in gv_by_gv angular. So I guess the instructions on the web page may be wrong. Or if it is my mistake, I hope you can point it out.
Have you compared the slope of phonon band structure and the group velocity in the same calculation? If they are inconsistent by 2pi, what you wrote it true.
Yes, I've calculated the group_velocity[5,1] from the slope method, and it turned out consistent with the existing one in hdf5. And this consistency is based on the fact that I consider the "frequency" in hdf5 ordinal just as said on the webpage, while the "group velocity" angular, contrary to the webpage.
Therefore, maybe you could correct the description on the webpage to avoid misunderstanding.
With the slope method, which physical unit have you used for the reciprocal space distance?
I used the same unit as the reciprocal lattice vector showed in the yaml file.
If what you wrote is correct, what is wrong is not the documentation, but the code.
I have checked the value, but I get what as I expect. In fact, in the documentation 2pi is omitted as a convention (phonopy/phonopy#142.) Maybe this is the reason of our disagreement. Otherwise can you give a python script to reproduce it? In addition, could you use the latest released versions and provide the calculation condition such as phonopy and phono3py versions.?
The code works fine here too. The fact that the slope method gave consistent value has proved the group_velocity is angular. As I have already converted the unit of frequency from ordinal to angular (that means a 2pi is already multiplied), the result of group velocity ought to be angular. So the problem here is not about code, but about the wrong "ordinal" description of group_velocity.
However, the "angular" and "ordinal" descriptions indeed differ by only one factor 2pi. If 2pi is usually omitted as you said, I guess the difference between "angular" and "ordinal" is negligible too?
We are not discussing the same point. I need to reproduce what is happening to your system. For this, I want to have a python script or a very solid way to reproduce your values.
I guess the difference between "angular" and "ordinal" is negligible too?
No.
The code of group velocity comparison is pasted below:
import h5py
import numpy as np
f = h5py.File("./kappa-m1511511.hdf5")
gv = f['group_velocity'][:] # group_velocity
fre = f['frequency'][:] # frequency(ordinal)
qp = f['qpoint'][:] # q points(reduced)
b1 = 2*np.pi*np.array([0.405181460517092,0.233931589981552]) # reciprocal lattice vector 1
b2 = 2*np.pi*np.array([0.00000000,0.467863179963104]) # reciprocal lattice vector 2
qp_2d = qp[:,0:2] # 2-D material, only consider in-plane coordinates
qp_2d_cart = (np.kron(qp_2d[:,0],b1)+np.kron(qp_2d[:,1],b2)).reshape(qp_2d.shape) # convert reduced coordinates to cartesian coordinates
gv_y = (fre[80,1]-fre[5,1])/(qp_2d_cart[80,1]-qp_2d_cart[5,1]) # group velocity y component(ordinal THz*A)
gv_x = (fre[79,1]-fre[6,1])/(qp_2d_cart[79,0]-qp_2d_cart[6,0]) # group velocity x component(ordinal THz*A)
gv_slope = np.array([gv_x,gv_y])
print("if consider group velocity ordinal as in doc\n then diff vector (gv_hdf-gv_slope):"+str(gv[5,1,:2]-gv_slope)) # consider group velocity ordinal just as said in doc
print("if consider group velocity angular\n then diff-vector (gv_hdf-gv_slope):"+str(gv[5,1,:2]/2/np.pi-gv_slope)) # consider group velocity angularThis code gives the following output:
if consider group velocity ordinal as in doc
then diff vector (gv_hdf-gv_slope):[100.98258118 56.12816668]
if consider group velocity angular
then diff vector (gv_hdf-gv_slope):[ 1.09656098 -1.54104526]
The index 5,6,79,80 are adjacent points in horizontal and vertical directions. Here it's apparently not a rigorous calculation, but just an approximation, which has clearly shown that the "THz" in group velocity is angular, not ordinal as in the doc.
The Cartesian coordinates of qpoints 5,6,79,80 are attached here:
5 [0.08429901 0.04867005]
6 [0.10115882 0.05840406]
79 [0.06743921 0.05840406]
80 [0.08429901 0.06813807]
I appreciate if you could write a bit more comprehensible code... Please discuss with your colleagues or supervisors how it can be more comprehensible when explaining.
But you multiply 2pi to reciprocal lattice vector. I guess you should not.
Sorry for my lack of experience in coding.
In the yaml file, there is a comment after the reciprocal lattice vectors saying that they are without 2pi. Since there is 2pi in the physical definition of the reciprocal lattice vectors, I multiplied them with 2pi. If I shouldn't do that, does that mean the reciprocal lattice vectors in yaml file have already been multiplied with 2pi?
2pi is a killer always. On wave vector, you may have a look at https://en.wikipedia.org/wiki/Wave_vector.
You are right. I did make a mistake about the 2pi in the previous Python code.
From the wiki, the crystallography wave vector without 2pi and ordinal frequency are bound together. And in that way, the group velocity calculated from the combination is the same as that from the angular frequency and traditional physics wave vector. That exactly means the group velocity is angular.
That exactly means the group velocity is angular.
No. See the definition, https://en.wikipedia.org/wiki/Group_velocity. Please discuss with your supervisors or colleagues about this. I will close this issue. Please reopen this if you need.