# Phonon dispersion calculation based on velocity autocorrelation function in reciprocal space?

I would like to ask about the method of calculating phonon dispersion relation. Up to my knowledge, there are 2 methods to calculate the phonon dispersion:

1. By diagonalizing the dynamics matrix. There are an open source code that implement this method in many molecular dynamics packages such as fix phonon in LAMMPS.

2. By Fourier transforming the velocity correlation function such that the velocity is in the reciprocal space. I want to apply this method. I found some papers that list the equations of calculation, however, the method is not clear. I wounder if some of you have experience about this topic .

One of the papers is: "Non-local modeling of epoxy using an atomistically-informed kernel"

They use the following equation to convert the velocity to the reciprocal space:

My question is about the position $r$ value, we know that the position is vector and has three components $(x, y, z)$, in the equation it is written but not clear if it is converted to a scalar value, and how? or they took the position as the same component of velocity. Another question is: as we will have three three correlation functions, one for each component $(x, y, z)$, Do they take their average or summation? So we will have one Fourier transform. If not we will have three three Fourier transforms and I don't know how to take the peaks from three Fouriers.

Please if any body have an idea about this topic, it is so helpful for me. Even if you know any code written for this purpose it will be helpful so I can compare it with my code to find the answers of my questions.

The impression I got form the equation is that you model your system as being composed of a collection of point particles. The index $i$ indexes this set of point particles (as opposed to indexing the components of a vector).

One could create a spatial velocity distribution by saying $\vec{v}(\vec{r})=\sum_i \vec{v}_i \delta(\vec{r} - \vec{r}_i)$, where $\vec{v}_i$ is the velocity of the $i$th particle and $\vec{r}_i$ is the position of the $i$th particle. Then the spatial fourier transform of this $\vec{v}(\vec{r})$ is $\vec{v}(\vec{k})$ or, as you have written it, $\vec{v}_\vec{k}$.

To get the dispersion relation, we put back in the time dependence of the velocities to get the function $\vec{v}_i(t)$. From this we form the velocity distribution $\vec{v}(\vec{r},t) = \sum_i \vec{v}_i(t) \delta(\vec{r} - \vec{r}_i)$. Now we can again do a spatial fourier transform to get $\vec{v}(\vec{k},t)$. We expect an oscillation to with wavevector $\vec{k}$ to oscillate at some frequency $\omega$, to find this frequency $\omega$, we take the temporal fourier transform of $\vec{v}(\vec{k},t)$ to get $\vec{v}(\vec{k},\omega)$. Now if we want to know the frequency corresponding to a given wavevector $\vec{k}_0$, then we just look at the frequency dependence of $\vec{v}(\vec{k}_0,\omega)$, and the $\omega$ for which this is non-zero give the frequency of the oscillation. So this procedure from getting an $\omega$ from a $k$ is a one way of getting the dispersion relation.

There is one more subtlety which is that $\vec{v}(\vec{k},\omega)$ has three components. Thus at a given $\vec{k}$ and $\omega$, the component of $\vec{v}(\vec{k},\omega)$ parallel to $\vec{k}$ may be non-zero but the perpendicular component may be zero. This tells you that longitudinal and transverse waves oscillate at different frequencies.

# Edit: How to get acoustic and optical bands.

I think the problem is when you do things this way, you don't put an upper limit on $\vec{k}$. Typically people restrict $\vec{k}$ to the first brillouin zone by means of defining a projection from the first brillouin zone by means of projection. Let me describe what I mean by an example.

Suppose we have a 1D system in 1D space. Suppose this system is a chain of ions whose species alternates between sodium and chlorine. Suppose the spacing between two adjacent ions is always $a$. Suppose the ions interact with a harmonic potential interaction: $U(x_\mathrm{Cl}, x_\mathrm{Na}) = \frac{1}{2}k(x_\mathrm{Cl}- x_\mathrm{Na}-a)^2$, where $k$ is a spring constant.

Now let us find the dispersion relation. As before, we look at velocity as a function of space and time. We do a spatial fourier transform to pick out velocity oscillations of a particular wavevector, and then we do a temporal fourier transform to get the frequency corresponding to each wavevector.

What do we expect to see? Let's start at $k=0$ and go up (the dispersion relation should have the symmetry $\omega(-k)=\omega(k)$). We expect $\omega(0) = 0$, and we expect that $\omega(k)$ should have a finite derivative as you increase $k$ from $0$. (This derivative is the speed of sound for this material.) Now we expect $\omega$ to keep increasing (albeit nonlinearly and even discontinuously at one point) until it gets to some maximum. This maximum occurs when the velocity oscillation of sodium is $pi$ out of phase with the velocity oscillation of chlorine. The half-wavelength of the velocity oscillation is $a$, so the wavelength must be $2a$ and the wavenumber must be $k_M=\pi/a$.

What happens if we increase $k$ past $k_M=\pi/a$? Well due to aliasing effects, a velocity oscillation with wave number $k_M + \Delta k$ gives the same particle velocities as a velocity oscillation of $k_M - \Delta k$ (since $\mathrm{Re}(\exp[i(k_Mna \pm \Delta k na)])=\mathrm{Re}(\exp[i(\pi n \pm \Delta k na)])=(-1)^n\mathrm{Re}(\exp[\pm i \Delta k na)])$, which is independent of the choice of plus vs minus). This tells us that $\omega(k_M + k) = \omega(k_M - k)$. Now we know what $\omega(k)$ is for $k< 2k_M$. To go beyond this we note that $\omega(k+2k_M) = \omega(k)$.

From the equations $\omega(-k)=\omega(k)$, $\omega(k_M + k) = \omega(k_M - k)$, and $\omega(k+2k_M) = \omega(k)$. It is clear that the represent the whole dispersion relation we need only consider $\omega(k)$ in the region $0<k<k_M$. This is one way of picking a "fundamental zone" of reciprocal space, and concentrating on that zone.

Now I still haven't gotten to acoustic and optical phonons. These come up when you pick the "fundamental zone" differently. What people usually done is to pick the fundamental zone to be half the size of the previously described one. That is, it is the region $0<k<k_M/2$. So that the information about the dispersion relation for $k>k_M/2$ doesn't get lost, we need to make two dispersion relations in the region $0<k<k_M/2$. One is the acoustic phonon dispersion relation defined by $\omega_a(k) = \omega(k)$ and the other is the optical phonon dispersion relation given by $\omega_o(k) = \omega(k_M-k)$, so that the optical phonon dispersion relation is the right half of the full dispersion relation folded over $k_M/2$.

So to get acoustic and optical phonons from the $\omega(k)$ in your example you need to find this second form of the "fundamental zone" in reciprocal space, which is conventionally called the "first brillouin zone", and you need to do the appropriate folding of $\omega(\vec{k})$ into this zone. This is probably a little trickier in 3d.

• Thank you for your answer. In fact I have applied the method you presented, and after calculating velocity in reciprocal space, I calculate velocity autocorrelation function and then Fourier transforming as you have mentioned.The correlation takes care of longitudinal and transverse differences. This method is giving me only the acoustic phonons and I need also the optical phonons, I dont know how! I am not physicist,I am computer science and need phonon dispersion to validate my simulation, so I try to find any reference explains a code to produce it – user35384 Dec 8 '13 at 17:44
• how many atoms do you have in your unit cell? – Brian Moths Dec 9 '13 at 3:43
• There are two carbon atoms in the unit cell. The material is one layer of graphene. The complete sheet of graphene used in my simulation is composed of 2508 carbon atoms. – user35384 Dec 9 '13 at 16:26
• @user35384 see my edit. Have you gotten optical phonons in the simple 1d lattice with two species? Also here is the obligatory wikipedia link to an article about the phonon dispersion relation – Brian Moths Dec 9 '13 at 18:28
• Thank you so much. Your explanation helped me alot. I will calculate the second form of fundamental zone, in my solution I have calculated the first one only. Also I was thinking that the maximum k is (2*pi/a) and I will change it now. About wikipedia method, it depeneds on the first approach (diagonalizing the dynamic matrix) because it depends on the constant forces that can be calculated from the potential function during the Molecular dynamics simulation. Thank you so much again :) – user35384 Dec 10 '13 at 6:02