Add Debye-Waller factors as an option of the phonon makers

[hrushikesh-s]

Summary

Include a summary of major changes in bullet points:

• Add the Debye-Waller factors as an option of the phonon makers. With create_debye_waller=True, a last job takes the thermal displacement matrix U of each site of the primitive cell at each temperature from the phonon output. pymatgen then computes the XRD, neutron and TEM patterns with the anisotropic Debye-Waller factors exp(-2π² gᵀUg) at each temperature, and once without them.
• New code: the compute_debye_waller job and the DebyeWallerDocument schema in common, and the create_debye_waller and debye_waller_kwargs fields of the phonopy BasePhononMaker. All phonon makers (VASP, aims, force fields, torchsim, pheasy and hiPhive) share its make, so they all get the option. The flow output is still the phonon document. The patterns are the output of the last job. A Debye-Waller subsection in the phonon part of docs/user/codes/vasp.md.
• U is taken from the phonon flow with create_thermal_displacements=True, computed as in Breaking: Leave out the acoustic modes at Gamma in the thermal displacement matrices #1590. The Debye-Waller job only computes the patterns. The q-point mesh, the temperatures and the frequency cutoff are set with kpoint_density_thermal_displacements, tmin/tmax/tstep_thermal_displacements and freq_min_thermal_displacements in generate_frequencies_eigenvectors_kwargs of the phonon maker, as for every phonon flow. If the mesh has imaginary modes, the phonon flow stores no U and the Debye-Waller job fails with a message that says so. With imaginary_modes_thermal_displacements="exclude" or "absolute" (Breaking: Leave out the acoustic modes at Gamma in the thermal displacement matrices #1590), U comes from the real modes, or also from the imaginary modes taken as real modes of the same absolute frequency, and the workflow runs. create_thermal_displacements is on by default for VASP and aims. For the other phonon makers, create_debye_waller=True also needs create_thermal_displacements=True, otherwise make raises an error.
• The patterns go to the data store. The document also stores the calculator settings and the TEM spots as typed rows (TEMSpot).
• Tests: an end-to-end EMT run for Cu checks that U equals the U of the phonon flow and that the XRD, ND and TEM intensities follow the Debye-Waller factors. A mocked VASP run for NaCl uses the NAC. Further tests check that the job gets the phonon output and the settings, and the error messages.
• This branch contains Breaking: Leave out the acoustic modes at Gamma in the thermal displacement matrices #1590, which should be merged first. The tutorial is in Add a tutorial for the Debye-Waller workflow #1592. Part of Tracking: atomate2 PRs from Anubhav Jain's group #1591.
• The physics and the first implementation are from Ruby Keesey (@rkeesey).

Additional dependencies introduced (if any)

• List all new dependencies needed and justify why. While adding dependencies that bring
  significantly useful functionality is perfectly fine, adding ones that add trivial
  functionality, e.g., to use one single easily implementable function, is frowned upon.
  Justify why that dependency is needed. Especially frowned upon are circular dependencies.

• This PR needs the anisotropic Debye-Waller factors of Add anisotropic Debye-Waller factors to XRDCalculator, NDCalculator and TEMCalculator pymatgen#4715 in the XRD, ND and TEM calculators. Until a pymatgen release has them, a debye-waller dependency group installs commit c6de1930 of that PR, since PyPI does not accept git URLs in extras.

• The strict extras pin pymatgen==2026.9.24, so CI installs the group in a second step. The new matrix entry test-forcefields (3.12, debye-waller) runs only the two Debye-Waller test files. It skips the pymatgen dispatch step, since pymatgen master would replace the PR commit there. The other jobs skip the tests that need the group.

TODO (if any)

If this is a work-in-progress, write something about what else needs to be done.

• Wait for Breaking: Leave out the acoustic modes at Gamma in the thermal displacement matrices #1590 and Add anisotropic Debye-Waller factors to XRDCalculator, NDCalculator and TEMCalculator pymatgen#4715. Once #4715 is in a pymatgen release, remove the debye-waller group, its CI entry, check_pymatgen and the skip conditions, and raise the pymatgen lower bound.

Checklist

Work-in-progress pull requests are encouraged, but please put [WIP] in the pull request
title.

Before a pull request can be merged, the following items must be checked:

• Code is in the standard Python style.
  The easiest way to handle this is to run the following in the correct sequence on
  your local machine. Start with running ruff and ruff format on your new code. This will
  automatically reformat your code to PEP8 conventions and fix many linting issues.
• Doc strings have been added in the Numpy docstring format.
  Run ruff on your code.
• Type annotations are highly encouraged. Run mypy to
  type check your code.
• Tests have been added for any new functionality or bug fixes.
• All linting and tests pass.

Note that the CI system will run all the above checks. But it will be much more
efficient if you already fix most errors prior to submitting the PR. It is highly
recommended that you use the pre-commit hook provided in the repository. Simply run
pre-commit install and a check will be run prior to allowing commits.

[JaGeo]

Could you check the frequency_cutoff and convergence with q points?

11 years ago, I did some relevant tests: https://doi.org/10.1039/c5ce01219h

EDIT:
I now reread the description that phonopy is not used in this context.
I would be interested in a comparison between this method and phonopy's with frequency_cutoff.

[JaGeo]

Also, there is a bit of a disconnect now. The phonon workflow already has an option to compute thermal displacement parameters and can save it as well

atomate2/src/atomate2/common/schemas/phonons.py

Line 137 in 60e6234

class ThermalDisplacementData(BaseModel):

[hrushikesh-s]

Could you check the frequency_cutoff and convergence with q points?

11 years ago, I did some relevant tests: https://doi.org/10.1039/c5ce01219h

EDIT: I now reread the description that phonopy is not used in this context. I would be interested in a comparison between this method and phonopy's with frequency_cutoff.

@JaGeo, yes, I am still running the benchmarking and comparison tests for this workflow from @rkeesey.
I will post the results here once they are ready and are properly checked from my end.

Also, there is a bit of a disconnect now. The phonon workflow already has an option to compute thermal displacement parameters and can save it as well

Yes, let me look into this as well...

[hrushikesh-s]

Could you check the frequency_cutoff and convergence with q points?

11 years ago, I did some relevant tests: https://doi.org/10.1039/c5ce01219h

EDIT: I now reread the description that phonopy is not used in this context. I would be interested in a comparison between this method and phonopy's with frequency_cutoff.

With a cutoff of 0 THz on the default q-point mesh, the acoustic modes at Gamma made U far too large or crashed phonopy, and U converges only as 1/N with the mesh number N. I fixed this for all phonon workflows in #1590, with a benchmark on 51 materials, and #1580 now takes U directly from the phonon flow, so there is no second implementation and no separate frequency_cutoff. Does this make sense?

[JaGeo]

This part I don't fully understand tbh

[hrushikesh-s]

An imaginary mode has no harmonic U, and leaving it out (“exclude”) lowers U. “absolute” keeps it as a real mode of the same absolute frequency, as an estimate of U. phonopy only sums modes above freq_min, so this block adds them in a second run, and the sign array undoes phonopy’s sign flip at T <= 1 K. I have added this to the comment. Does it make sense now?

[hrushikesh-s]

I have also added why we need this for the Debye-Waller workflow. Without U, a material with imaginary modes gets no Debye-Waller factors, and 18 of our 51 benchmark materials have imaginary modes at the default q-point density. The comment and the docstrings now say this.

[hrushikesh-s]

Could you check the frequency_cutoff and convergence with q points?

11 years ago, I did some relevant tests: https://doi.org/10.1039/c5ce01219h

EDIT: I now reread the description that phonopy is not used in this context. I would be interested in a comparison between this method and phonopy's with frequency_cutoff.

@JaGeo, I have now also compared this with phonopy’s frequency cutoff in two tests:

First --> the 51 materials of #1590. A cutoff of 0.01 THz gives the same U as leaving out the acoustic modes at Gamma. A cutoff of 0.5 THz makes U about 8% too small, and a denser mesh does not fix this.
frequency_cutoff_benchmark.zip

Second --> pentachloropyridine from your 2015 paper. I took your 100 K structure (Crystallography Open Database entry 7221873), kept its cell, and relaxed the atoms with MACE-OMAT-0-medium, which has no dispersion term. I computed the force constants on your 4x4x2 supercell, and U from 100 to 300 K on q-point meshes from 5x5x1 to 71x73x25, and on your 70x70x28 mesh. pentachloropyridine_comparison.zip

With #1590, U converges smoothly as the mesh gets denser. At 71x73x25, it is within 1% of U on your mesh without a cutoff.
With your 0.13 THz cutoff, U is about 6% smaller on every dense mesh, because the cutoff removes real modes.
At 100 K, U from #1590 agrees with your measured U_eq within 6% (mean |ΔU_ij| 0.0007 Ų). With the 0.13 THz cutoff, the mean |ΔU_ij| is 0.0010 Ų. Above 100 K, U is too small, as expected for one fixed cell and harmonic phonons.

Therefore, I have kept freq_min at 0 and left out only the acoustic modes at Gamma. Does this make sense?

[JaGeo]

@hrushikesh-s This looks already great. I would nevertheless recommend a larger benchmark once you make a production run!

Naive question: do we need a separate workflow or should this rather be another option of the phonon workflow? The only difference between the phonon workflow and this one is an additional postprocessing step and additional entries in the output document, right? I am worries about too many workflows and that they might be easily overlooked otherwise.

[hrushikesh-s]

@JaGeo, thanks, I will check with @rkeseey, and request her to run a larger benchmark with this workflow.

on your question on whether we need a seperate workflow -->
The Debye-Waller step only needs U, and U can come from any phonon flow (phonopy, pheasy, and hiPhive now, and later the finite-temperature phonons).

Should I drop the Debye-Waller makers and keep only the compute_debye_waller job, documented in the phonon section? Users would then add this job after their phonon flow themselves.

Or would you rather have it as an option of the phonon makers? Then the phonopy, pheasy, and hiPhive makers, and every new phonon maker, would each need the option in their own make(). Until a pymatgen release has materialsproject/pymatgen#4715, all phonon flows would also depend on it. Or is there something else that you suggest to do instead?

[JaGeo]

Or would you rather have it as an option of the phonon makers? Then the phonopy, pheasy, and hiPhive makers, and every new phonon maker, would each need the option in their own make(). Until a pymatgen release has materialsproject/pymatgen#4715, all phonon flows would also depend on it. Or is there something else that you suggest to do instead?

I had this in mind. However, the options would be part of the initialisation. I would anyhow recommend that we wait with merging until a new pymatgen release exists.

[hrushikesh-s]

Or would you rather have it as an option of the phonon makers? Then the phonopy, pheasy, and hiPhive makers, and every new phonon maker, would each need the option in their own make(). Until a pymatgen release has materialsproject/pymatgen#4715, all phonon flows would also depend on it. Or is there something else that you suggest to do instead?

I had this in mind. However, the options would be part of the initialisation. I would anyhow recommend that we wait with merging until a new pymatgen release exists.

I changed it as you suggested. The Debye-Waller factors are now an option of the phonon makers: PhononMaker(create_debye_waller=True, debye_waller_kwargs={...}). The separate Debye-Waller makers are removed. The flow output is still the phonon document.

[hrushikesh-s]

The tutorial in #1592 now also uses this option, and its rerun gives the same results for the six materials.

[rkeesey]

@hrushikesh-s Phonopy's default TDM matrices use Cartesian coordinates, and the CIF convention is only applied when specifically requested (the YAML files with the TDMs in the phonon database include both conventions). Would it be beneficial to implement a method of DWF calculation using Cartesian coordinates as well?

[hrushikesh-s]

@hrushikesh-s Phonopy's default TDM matrices use Cartesian coordinates, and the CIF convention is only applied when specifically requested (the YAML files with the TDMs in the phonon database include both conventions). Would it be beneficial to implement a method of DWF calculation using Cartesian coordinates as well?

@rkeesey, thanks! The Debye-Waller job already starts from phonopy’s Cartesian U, and its output stores both conventions. It converts U to the CIF convention only because pymatgen’s XRD, ND, and TEM calculators read the U11_cif to U12_cif site properties. Both are the same tensor in different bases, so the factors match exp(-2π² gᵀUg) from the Cartesian U to machine precision, which I checked on a triclinic cell.

[hrushikesh-s]

Or would you rather have it as an option of the phonon makers? Then the phonopy, pheasy, and hiPhive makers, and every new phonon maker, would each need the option in their own make(). Until a pymatgen release has materialsproject/pymatgen#4715, all phonon flows would also depend on it. Or is there something else that you suggest to do instead?

I had this in mind. However, the options would be part of the initialisation. I would anyhow recommend that we wait with merging until a new pymatgen release exists.

@JaGeo, materialsproject/pymatgen#4715 is now merged.