Phonon band structure using Phoebe

[wugaxp]

Dear developers,

Thank you very much for the great code!

I'm trying to use phoebe to estimate the phonon transport behavior. Currently I'm following the instruction on https://phoebe.readthedocs.io/en/develop/tutorials/mlPhononTransport.html and try Si as a first example. With the python scripts, I can generate the IFCs as fc2.hdf5 and fc3.hdf5 files. However, when I use phoebe to generate the phonon band structure, incorrect results are obtained. It seems that some phonon branches show imaginal frequencies.

The input files and the obtained phonon dispersion can be found in
Phoebe.zip

Could you please advise what I can done wrong? Thanks.

Looking forward to your suggestions.

Regards,
Gang

[jcoulter12]

Hi Gang,

We're glad you're enjoying Phoebe, and thanks for reaching out about this!

Your Phoebe input file (phband.phoebe) looks correct for the phonon bands app calculation, so I don't think you have any issue with the actual Phoebe calculation. The problem is likely coming from the force constants.

If your potential isn't of a good quality, the force constants could come out quite wrong, and then of course you could see an incorrect phonon dispersion. Classical Tersoff potentials may not work well enough for accurate phonon calculations, though I see here the regular Tersoff Si potential gets the phonons quite wrong, but they are stable.

Some suggestions to diagnose the issue:

• You can check if this is a "workflow" issue with the tutorial or your ifc.py script, versus an issue with the potential by modifying lines 48-63 of your ifc.py script to instead use DFT to calculate forces for each generated supercell. (As in, changing the atoms.calc line of the ASE calculator from a lammps call to a DFT code call).
• An alternative test to the above -- do you find that using a different potential drastically changes these results? That could be a sign if the problem is potential related or not.
• Another possible issue -- did you check the cutoff pair distance makes sense? To check that this is or isn't an issue, simply remove the cutoff_pair_distance argument.
  ph3.generate_displacements(cutoff_pair_distance=3.0)

Let me know what you learn after checking a few of these, and we can think about what else to try.
Best,
Jenny

[wugaxp]

Dear Jenny,

Thank you very much for your reply!
Actually I have benchmarked the performance of the forcefield. I used another code called phonolammps. With the same potential, it can produce reasonable phonon dispersion.

As you can see, no imaginal branches can be observed.

I have also tried to increase the cutoff_pair_distance to 6.0 Ang, although I believe that it only affect the interaction cutoff for 3rd order FCs. Also I tried to use a larger supercell (phonon_supercell_matrix=[8, 8, 8]) to accommodate the large cutoff.

Unfortunately, the phonon dispersion still looks weird.

Meanwhile, I notice that Prof. Togo has released a project cogue on github, in which he has a module for the calculation of 2nd and 3rd order IFCs (https://github.com/atztogo/cogue/blob/develop/cogue/task/phonon_fc3.py, lines 278-301). It is interesting to see that, although the displacements are generated using phono3py, the 2nd order IFCs seem to be calculated using phonopy. I'm not sure if the error comes from the different in phono3py and phonopy.

Regards,
Gang

[wugaxp]

I tried another empirical potential for Si, the Stillinger-Weber potential.
With PhonoLammps, the phonon spetrum is quite consistent with the one from Tersoff potential:

With IFCs.py+phoebe, the strange imaginal frequencies become more severe.

[YuuuXie]

Hi Gang, have you ever tried using phonopy instead of phono3py to generate 2nd FC and maybe use that in Phoebe to get phonon dispersion? I had the example that phonopy and phono3py generate different displacements for the FC2 calculation with the same system and unit cell, so maybe that could cause an issue

[jcoulter12]

Also, it would be great to see a phonon dispersion on the full band path of silicon -- I'd like to see how the dispersion matches up to expected results.

[jcoulter12]

Hi Gang,

Thank you for this detailed, thorough report of the issue -- it's very helpful to see this. I agree, it's definitely not the potential.

@YuuuXie, who replied in a comment and is the author of the tutorial, had some thoughts on what might cause the issue. Indeed, it may be some difference in the behavior of phonopy and phono3py.

I'm going to do some tests for myself and will get back to you soon with either more questions or a fix for this.

Thanks again for reporting the issue, we really appreciate it.
Jenny

[wugaxp]

Dear @YuuuXie and @jcoulter12,

Thank you for the suggestion. I tried to use phonopy instead of phono3py to generate the fc2.hdf5 file, and I was surprisingly to see that the phonon dispersion becomes correct!

Here I attach the fc2.hdf5 and disp.yaml files from phono3py and phonopy. However I cannot observe the difference which may cause the strange problem.
Si-fc2.zip

Have you ever tried to validate the third-order IFCs? For example, have you compared the relaxation time from phono3py+phoebe and Alamode/ShengBTE before? I prefer phoebe because the whole flow chart can be run with python scripts, the other two cannot (or I did not figure out the correct way). Do you have any suggestion on the interface between phoebe and Alamode/ShengBTE?

Regards,
Gang

[wugaxp]

Finally, I figure out the problem. When we deal with the IFCs in Phoebe, we are required to provide the disp.yaml file. If the supercell_matrix and the phonon_supercell_matrix parameters are set to different values, the disp.yaml file for the 2nd order should not use the file generated by phono3py, which is in a different supercell. Instead, if we feed in a yaml file generated by Phonopy, we can obtain the correct dispersion. Otherwise we must unify supercell_matrix and the phonon_supercell_matrix for a consistent disp.yaml file.
Other parameters are not so critical. I have tried a supercell of 2x2x2 for 2nd and 3rd order IFCs, and a cutoff_pair_distance of 3.0 Ang, the resulting phonon spectrum is still acceptable without imaginal frequency. I will also test the phonon lifetime and kappa for the convergence later.

Thank you very much for your help!

[jcoulter12]

Hi Gang,

Thanks to your helpful comments, I was able to pin down the issue! We really appreciate that you reported this.

The reason it works when the fc2 and fc3 supercells are the same is because the write_fc2_to_hdf5 block of code needs to be changed to:

write_fc2_to_hdf5(
  ph3.fc2,
  p2s_map=ph3.phonon_primitive.p2s_map,  # changed from ph3.primitive.p2s_map
  physical_unit="eV/angstrom^2",
)

We were previously writing the primitive cell to supercell mapping for the fc3s to fc2.hdf5. Changing that line to ph3.phonon_primitive.p2s_map in the write_fc2_to_hdf5 function instead should fix the problem -- please let me know ASAP if it doesn't. I'll also update the tutorial accordingly!

Additionally, on your question about validation vs other codes -- yes, we've checked that our lifetime and transport properties come out equivalently to ShengBTE and Phono3py in the past. An example of a comparison for GaN lifetimes between Phoebe (purple points) and Phono3py (black points) is attached (note, they aren't exactly the same in large part because they're done on different meshes, to make it easier to see that the Phono3py points sit on top of the Phoebe ones. Some additional small things, like ASR choice/smearing aren't identical in this particular comparison example, either).

I'm going to also double check that the fc2.hdf5 and fc3.hdf5 files generated by this script produce a thermal conductivity which matches the command line use of Phono3py (as in the Phonon Transport tutorial on the Phoebe website). I feel certain this will be fine, but I'll do it out of an abundance of caution.

Thank you again for reporting this in such a helpful and detailed manner, and do let us know if you need any other help using Phoebe.
Best,
Jenny

[jcoulter12]

Hi Gang,

Thanks also for noticing this -- fortunately, this is much simpler to figure out. :)

Just add
linewidths = linewidths[calcIndex].flatten()
To where tau and energies are flatted in the "noWindow" function at the top of the script. I'll also push this change to the github repository.

I simply never run calculations with noWindow, so I've never noticed this. In phonons, the window doesn't make a different at high temperatures because no phonon states are frozen out, but if you run at low temperatures or ever run an electron calculation, please be sure to add one, as it can reduce the calculation cost significantly.

Thanks!
Jenny

[wugaxp]

Hi Jenny,

Similar error appears when I run the projection of relaxation time to bands.
python ~/Install/phoebe/scripts/plotScripts/tau_path.py rta_ph_relaxation_times.json phonon_bands.json 0
produces

Calculation Temperature:  10.0 Calculation Chemical Potential: 0.0
Traceback (most recent call last):
  File "/home/wug/Install/phoebe/scripts/plotScripts/tau_path.py", line 184, in <module>
    punchPlotTau(plotFileName, tau, points, pathTicks, pathLabels)
  File "/home/wug/Install/phoebe/scripts/plotScripts/tau_path.py", line 28, in punchPlotTau
    plt.plot(x, y, label="band #{}".format(ib+1))
  File "/opt/intel/oneapi/intelpython/latest/lib/python3.9/site-packages/matplotlib/pyplot.py", line 3019, in plot
    return gca().plot(
  File "/opt/intel/oneapi/intelpython/latest/lib/python3.9/site-packages/matplotlib/axes/_axes.py", line 1605, in plot
    lines = [*self._get_lines(*args, data=data, **kwargs)]
  File "/opt/intel/oneapi/intelpython/latest/lib/python3.9/site-packages/matplotlib/axes/_base.py", line 315, in __call__
    yield from self._plot_args(this, kwargs)
  File "/opt/intel/oneapi/intelpython/latest/lib/python3.9/site-packages/matplotlib/axes/_base.py", line 501, in _plot_args
    raise ValueError(f"x and y must have same first dimension, but "
ValueError: x and y must have same first dimension, but have shapes (351,) and (1000,)

And running
python ~/Install/phoebe/scripts/plotScripts/tau_path_color.py rta_ph_relaxation_times.json phonon_bands.json 0
leads to

Traceback (most recent call last):
  File "/home/wug/Install/phoebe/scripts/plotScripts/tau_path_color.py", line 88, in <module>
    cf = ax.scatter(allPoints, (energies-mu), marker='o',
ValueError: operands could not be broadcast together with shapes (351,6) (15,)

projtau.zip

Regards,
Gang

[jcoulter12]

Hi Gang,

I think it's likely you're incorrectly applying this plot script. This is a script for plotting a band structure with lifetimes/linewidths projected onto it, for which you need the output of this tutorial.

If you want to make a lifetimes vs energy plot like the one I posted above, you need to use either tau.py or tau_bandResolved.py.

Let me know if that fixes it,
Jenny

[wugaxp]

Hi Jenny,

Yes, it works! I'm sorry for not reading the tutorial carefully -_-b
Both tau.py and tau_bandResolved.py can be run smoothly without any problem.

Regards,
Gang

[jcoulter12]

No worries -- happens to the best of us.
It's also an indication I should make the instructions for the plot scripts a bit more detailed, so it's also helpful to know the docs are confusing on this.

[wugaxp]

Hi Jenny,

In addition to the above question, I'm still struggling in two problems.

First is about the parallel computation of the forces. I don't want to take the mpi/omp parallelism of LAMMPS. Instead, I'd like to scatter the calculation of the displaced configurations to different subprocesses, which can achieve ideal parallelism in principle. Meanwhile, I like your use of tqdm. Now I'm trying to use the multiprocessing module in pathos lib. Do you have any suggestion on the better realization?

Second is about the compilation of Phoebe on GPGPU. I prefer Kokkos acceleration, however the compilation fails when I use the Nvidia hpcsdk + intel MKL. Can you suggest the best compilation environment? I will try first and post more details if I still fail.

Regards,
Gang

[jcoulter12]

Hi Gang,

I don't often try to parallelize python this way, so I'm going to think for a moment about the first question and get back to you.

On the second question, it's a bit hard to know what went wrong without seeing the compilation error. Definitely also make sure you use the appropriate -DKokkos_ARCH flag for your architecture. I see no reason to expect compilation difficulties with NVIDIA hpcsdk + intel MKL, however.

Thanks,
Jenny

[wugaxp]

Hi Jenny,

I tried the compilation again. My system is CentOS8, and I have Nvidia hpcsdk 22.5 and intel oneapi installed.
Before installation, I source the following files:

source /opt/intel/oneapi/intelpython/python3.9/env/vars.sh
source /opt/intel/oneapi/mkl/latest/env/vars.sh

CUDA_HOME is also set correctly.

Then I run following line under a build_gpu folder in phoebe source code

CC=mpicc CXX=mpicxx FC=mpif90 cmake .. -DHDF5_AVAIL=ON -DMPI_AVAIL=ON -DOMP_AVAIL=ON -DKokkos_ENABLE_CUDA=ON -DKokkos_ARCH_TURING75=ON -DCMAKE_CXX_STANDARD_LIBRARIES="-L/usr/lib64/openmpi/lib/" -DCMAKE_CXX_FLAGS="-I/usr/include/openmpi-x86_64"

Since I have two Nvidia Geforce RTX 2060 installed, I set cuda arch as TURING75.

However, cmake reports errors:

-- Configuring done
CMake Error at CMakeLists.txt:61 (add_executable):
  Target "runTests" links to target "CUDA::cusolver" but the target was not
  found.  Perhaps a find_package() call is missing for an IMPORTED target, or
  an ALIAS target is missing?

CMake Error at CMakeLists.txt:60 (add_executable):
  Target "phoebe" links to target "CUDA::cusolver" but the target was not
  found.  Perhaps a find_package() call is missing for an IMPORTED target, or
  an ALIAS target is missing?

-- Generating done

Previously I can generate the Makefile though, however I forget the setup... If I'm right, last time I failed when compiling files with HDF5 support. I will try again later.

Regards,
Gang

[anjohan]

Hi Gang,

1. Parallelization: Since you only need to do a handful of LAMMPS calculation for relatively small systems, this should not take very long. But if you are familiar with MPI, mpi4py provides an easy solution.

2. Kokkos: I have never tried to use NVIDIA HPC SDK, but in theory, it should work. I always just use GCC + CUDA. The specific complaint is about the cusolver library, which we use to diagonalize many small matrices in parallel. Do you have it installed?

Best wishes,
Anders

[wugaxp]

Hi Anders,

Thank you for the reply. I was busy these days, so I will try your suggestions and get back again.

1. Yes, I'm also considering the parallelization through something like mpi4py. Although the system itself may not be very large, the number of calculations can be huge if the cutoff for higher-order FC is not small. For example, for a (3x3x3) supercell of monoclinic ZrO2 and a cutoff of 3.0 Ang, more than 30000 configurations must be done. Thus it is helpful to parallel the computations for the configurations, especially when the symmetry is low. Previously I preferred parallelization through python multiprocessing because I need to use tqdm, which is not easy to be used with mpi4py.
2. Thank you for the suggestion. Actually I raised the question because I need to know the exact environment used to do the compilation. Nvidia HPC SDK does include the CUDA lib, but I will try the GCC+ openmpi + CUDA instead.

Best wishes,
Gang