Seeking guidance on sampling strongly anharmonic materials

[Jarvis-Towne]

Hello, I am using TDEP to calculate temperature-dependent phonon spectra for strongly anharmonic materials, specifically PbTe. My current results deviate significantly from literature/experimental values, and I suspect my methodology, particularly the sampling, is flawed.

Here are my current steps:

1. AIMD Simulation: NVT ensemble, 600K and 1000K, 50 ps for PbTe.
2. Data Extraction: process_outcar_5.3.py OUTCAR
3. Force Constant Fitting: extract_forceconstants -rc2 5.5 -rc3 4 > 1.log
4. Phonon Calculation: phonon_dispersion_relations --gruneisen --dos --qpoint_grid 20 20 20 --sigma 5 > 2.log

I have the following questions:

1. R-squared for Force Constants: The R-squared values for both second and third-order force constants are exceptionally low (e.g., 0.002 at 1000K). Does such a low R-squared definitively indicate unreliable results? What is generally considered an acceptable R-squared for force constant fitting in TDEP?
2. Sampling Strategy: Based on discussions in the TDEP issues, I suspect my AIMD sampling method is inadequate for strongly anharmonic materials. What is a more robust and recommended sampling strategy for such systems?
3. System Size and Data Requirements: If I were to restart, what are the recommended supercell sizes (number of atoms) and AIMD simulation durations (amount of data/snapshots) for strongly anharmonic materials like PbTe?
4. Reference Sampling Method: I have attached a PDF (TDEP-ref.pdf) detailing a sampling method from a published work. Is this method generally considered reasonable for TDEP, and can I strictly follow it?

I am also attaching input files and output files below.

Thanks,
Jarvis

PbTe-test.zip

[OrbitalC]

Hello Jarvis,

Looking quickly at your input, I'm noticing that you are using the conventional unitcell instead of the primitive one in your infile.ucposcar
You should really use the primitive one, since it has a lower number of atoms -> higher number of symmetries -> less parameters -> better convergence
For PbTe - if you are doing the rocksalt structure - you should have only two atoms in the unitcell

To answer your questions:

1. The R-squared question is a bit delicate because it's really dependent on the system being studied. In a way, TDEP is not really fitting in a conventional way, because the phonons obtained have a physical meaning/are physical observables (check our papers for more info: https://doi.org/10.1063/5.0174255 and https://doi.org/10.1103/PhysRevB.111.094306) So a low R-squared has a physical meaning. Though, the R-squared you're showing at 1000K is particularly low. But it might be because you are particularly close to the melting temperature where the phonon picture doesn't really mean anything anymore. Also I'm more used to look at the anharmonicity measure (https://doi.org/10.1103/PhysRevMaterials.4.083809 for a reference of what it means) and something close to 1.0 as what you have is not unexpected. It simply means that you have a very anharmonic material.
2. AIMD should be be completely fine for anharmonic systems and I'm not seeing any problem in your process. Even the time scale for the sampling should be enough (if you remove the burn-up phase of the AIMD to be sure to have a canonical averaging). I would expect your problem to come from the size of the supercells and the use of non-primitive unitcell.
3. The size of the supercell is a property to converge as it directly relates to the maximum cutoff available for a simulation, which means that it is system dependent. For metals it often converges around 100 atoms and around 200300 atoms for semiconductor but again, it is system dependent. Fortunately, increasing the number of atoms in the supercell increase the data per configurations, so higher number of atoms also mean faster convergence with respect to the number of configurations (as long as these configurations are a good sample of the canonical ensemble)
4. The method is perfectly valid. The idea is to reduce the cost of the sampling by using Gamma-point (or reduced k-point) calculations to obtain the canonical ensemble and then extract configurations to then compute with a tight convergence on the DFT. So there is a bit of an approximation (the canonical ensemble is approximated) but it should work. Though you have to ensure that the approximated ensemble is not too far from the true one. In other words, you would have to check that Gamma-point only forces (for atoms with displacements) are close to full k-point forces for the same configurations. should work really fine for semi-conductor, but for metals it might still be necessary to have some k-point for the sampling.

Best,
Aloïs

[Jarvis-Towne]

Thanks for the prompt reply! I will proceed with further tests based on your suggestions.