different harmonic phonon dispersions computed from kaldo and phonopy

[ZengZezhu]

Dear kaldo Team,

I am trying to use the kaldo (wonderful package to do thermal transport simulations) to compute the thermal conductivity of a disordered crystal.

What interesting is that I found the phonon frequency calculated based on the Kaldo and Phonopy has a large difference.
The maximum phonon frequency computed from phonopy is about 8THz, while kaldo only gives about 6THz. I am confused about the source of this difference.

I firstly did the 2nd and 3rd force constants calculations using lammps:

using MLP with lammps to compute the 2nd and 3rd IFCs, here is the input file:

units metal processors boundary p p p atom_style atomic atom_modify map yes box tilt large read_data lmp.data pair_style nep nep.txt pair_coeff * print 'calculating dynmat' dynamical_matrix all eskm 1e-06 file Dyn.form binary no print 'calculating third_order' third_order all eskm 1e-06 file THIRD binary no

then used the kaldo to compute the harmonic frequency and plot the phonon dispersion with this below code:

After doing these, I can successfully get the phonon dispersion and thermal conductivity. However, when I use the phonopy + the same MLP to compute the phonon dispersion, I found there is a large difference between two results from Phonopy and Kaldo, respectively. I attached the figures below for your review:

from kaldo:

from Phonopy:

I tried to understand this by doing some tests:

1. I changed the displacement magnitude of the finite difference method in lammps from 1e-6 to 1e-2 (because I know that Phonopy usually moves the atoms with a displacement of 1e-2), but the phonon dispersions are no distinct changes.

2. I understand that the k-paths are different from two packages to plot the phonon dispersion, while the maximum frequencies at Gamma point are also different from two packages, which is unacceptable.

   If it's possible, could you give me some advices to understand this? Many thanks in advance.

Cheers, zezhu

[niklundgren]

Hi Zezhu, If I'm reading your dispersion from kALDo correctly, it looks like the acoustic bands have zero frequency along the entire dispersion. At the very least they are extremely small. Here's my two lead suspicions-

1. Energy minimization I notice in your LAMMPS script, you did not do an energy minimization before you calculated the dynamical matrix. I would highly recommend optimizing again in LAMMPS and rerunning the dynamical matrix. Sometimes translating between two force calculators (LAMMPS vs GPUMD) introduces small errors and a minimization should take very little time. You can add it into your script right before the dynamical matrix calculation with the following: """ min_style cg minimize 0 1e-7 10000 20000 write_data replicated_atoms.lmp dynamical_matrix all eskm 1e-2 binary no file dyn.lmp """ Then you can use this "replicated_atoms.lmp" file to generate the "replicated_atoms.xyz"

2. Consistent Atom Ordering Another tricky thing I'm working on incorporating into our docs is that LAMMPS does not output atoms sorted by atom ID, however, this is a requirement for kALDo. If you converted LAMMPS output, you need to ensure you sorted the atoms correctly which can be done with the ASE python package if you'd like. Here's the way to read+write with the correct sorting """ from ase.io import read atoms = read('replicated_atoms.lmp', format='lammps-data', style='atomic', sort_by_ID=True, Z_of_type={...}) atoms.write('replicated_atoms.xyz') """ Make sure to edit the dictionary for the "Z_of_type" argument! https://wiki.fysik.dtu.dk/ase/_modules/ase/io/lammpsdata.html

If you think those two factors were taken care of correctly: Would you be willing to share the xyz file and dynamical matrix? I may have a better answer for you if I can confirm a few things myself.

[ZengZezhu]

Hi Nik,

Thanks so much for your information and suggestion!

1. I do agree with you that energy minimization is important before doing phonon calculations. While I use the same initial structure to do the calculations in Phonopy and lammps, and then plot the phonon dispersion to make a comparison. In this case may we guess that there is no distinct difference for two phonon dispersions? Actually I got my initial structure after doing the relaxation in ASE. I will follow your suggestions and do the relaxation in lammps, and test its phonon dispersion.

2. Thanks for the reminder. I do notice this issue as I had a similar experience when I used the Dynaphopy with Phonopy to get the phonon properties. Therefore I write my lammps data files using ASE (write_lammps_data function).

I will do the test you mentioned firstly, and if still does not work, I am very willing to send my all input files to you! Thanks a lot. I am also discussing with Zekun on the Wechat:-) He also gives me some hints to do more tests.

Cheers, zezhu

[ZengZezhu]

Hi Nik,

I maybe still a little bit confused but I have a naive idea now. I guess both of these two "phonon" dispersions are "correct". Maybe my understanding is wrong.

Today I just found a interesting result, if I use the 2nd and 3rd force constants generated from Phonopy and Phono3py within finite displacement method to compute the thermal conductivity of my disordered system with 112 atoms using the below command line:

phono3py --cell POSCAR --dim 1 1 1 --mesh 2 2 2 --fc3 --fc2 --tmin=100 --tmax=400 --tstep=100 --sym-fc --isotope --br --wigner

I found that the maximum vibrational frequency generated by the Phono3py is also around 6.4 THz. It is not 8.4 THz that I showed above using Phonopy API to plot! Which means the the results of maximum vibrational frequency computed from kaldo and phonopy could be consistent.

I guess when I used the Phonopy API with python to plot the dispersion, due to I manually set up the k-path, the Phonopy can recognize symmetry of my disordered structure and treat it as a crystal. While in the kaldo, it has no symmetry detected and it can be regarded as amorphous material. In addition, when I use the phono3py to compute the kappa, due to I did not explicitly give the k-path, the phono3py also does not recognize its symmetry.

In fact, for the disordered crystal, I guess there is no phonon dispersion? the only results we can access is the vibrational dos at Gamma point? If my understanding is right, then the "phonon" dispersion calculation based on Phonopy is also not meaningful.

Could you please send me your email? I will send my input files to you by email. Thanks.

Cheers, zezhu

[gbarbalinardo]

Hi @ZengZezhu. I'm just catching up on this interesting conversation. I agree that the dispersion relation makes less sense when you use QHGK because you are doing the calculation at gamma. Can you please send your files to giuseppe.barbalinardo@gmail.com or share a Google Drive link with me (same email)?

Cheers, Giuseppe