Halide perovskite CsPbBr3 soft mode frequencies greatly underestimated

[kw600]

Dear developers,

I am doing the anharmonic calculations for the high temperature cubic structure of halide perovskite CsPbBr3, which is regarded as a strong anharmonic material. The experimental Tc (cubic to tetragonal) is around 410K. The soft modes of this system at high temperature locate at M(0.5,0.5,0) and R(0.5,0.5,0.5) points, which are commensutae with a 222 supercell.

I have done anharmonic phonons at 400K, 500K and 700K. The SSCHA matrix (positive definite) seems to be correct since the lowest mode at M and R points have a frequency around 30 inverse cm and its eigenvector is the same as the imaginary mode in harmonic level. After applying the bubble correction (using ensemble.get_free_energy_hessian), this mode becomes -20 inverse cm even at 700K. And it seems increasing temperature only slightly make it less negative and does not help make the mode stable. This is a bit werid since usually soft modes are quickly stabilized by temperature. So I try to set include_v4 to be True in the function of ensemble.get_free_energy_hessian. In this way, at 500K the mode at R point becomes positive and the mode at M point is -8 inverse cm, which makes the soft mode more stable but still not enough. I am sure the free energy hessian is fully converged in the order of 1000 configurations (kong_liu_ration > 0.8) and increasing energy cutoff, kpoints and SOC does not help.

Interestingly, there is a PRL work in 2022 (https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.129.185901) claimed the Tc using self-consistent phonon (cite SSCHA and other anharmonic codes) is 177K and the soft mode frequency at M point is greatly overestimated (18 cm -1) (Fig1c), which is in the opposite direction as my results (-20 cm-1). And they claim the overestimation is due to the lack of bubble self energy in the self-consistent phonon. As far as I know, SSCHA has the bubble correction in the free energy hessian and also has the four phonon scattering, which motivates me to use include_v4 in the enesemble.get_free_energy_hessian but still not stable.

Thank you in advance for any suggestions on this system.

Best wishes,
Kang

[mesonepigreco]

Hi,
So, we have done some work on metal halide perovskites (namely on CsSnI3, see, e.g., L. Monacelli and N. Marzari, Chemistry of Materials, 35, 1702, 2023).
We found that the inclusion of the fourth order is necessary in this system, otherwise the high symmetry structure is always unstable. However, the number of configurations required to converge the fourth order are higher than those required to converge the SCHA minimization, therefore it is possible that 1000 configurations is not enough.
To estimate the error on the frequency, you can split the ensemble, as indicated in the FAQ: How do I understand if the free energy hessian calculation is converged?

I can suggest you to downsample your ensemble from 250, 500, and 1000, and track how the frequency of the low energy mode changes to determine if you are converging or you need more configuration with the fourth order. Also, a bigger cell than the one employed could be needed.

There are also other possible explanations. The transition may depend on the XC employed, which metal halide perovskites are very sensitive to, especially in the low-energy mode (see, for example, W. Kaiser et al, J. Phys. Chem. Lett. 12, 11886, 2021). Alternatively, the mechanism of the phase transition could be of a different kind than a displacive, as it happens in BaTiO3 (see, e.g., L. Gigli et al, npj computational materials, 8, 209, 2022).

Hope this helps,
Lorenzo

[kw600]

Dear Lorenzo,

Thank you very much for your helpful suggestions.

For the convergence testing, is it the same as ensemble split if I use ensemble.load with the N_random=250, 500, 1000 configurations and then do the minimization and update the ensemble weight?

Another question is the inclusion of fourth order. When I only include third order, the ensemble.get_free_energy_hessian generates d3 and save it into the file d3_realspace_sym.npy, which is then used to construct FC3. With FC3, the tensor3 can be constructed to calculate CC.Spectral.get_static_correction_interpolated. I can see d4 is calculated when we set include_v4 to be True and we can also save d4 to a npy file. I am just wondering how can I do the CC.Spectral.get_static_correction_interpolated with the correction from d4? Currently, the FC3 is only constructed from d3.

Best wishes,
Kang

[mesonepigreco]

Hi, there is currently no way to perform the interpolation using the d4. The problem is that the interpolation of d4 is much more complex than the one of d3, and several subtle issues arise on its acoustic sum rule and symmetries. This is one of the directions we are actively exploring and doing research on.

The solution right now is to use the hessian matrix that comes out from get_free_energy_hessian as a reference for the phase transition. This is usually reasonable, as fourth-order anharmonicity converges quicker with the q-mesh, and it is what we did in the work on CsSnI3.

[kw600]

Dear Lorenzo,

Many thanks for telling me this. I will use get_free_energy_hessian and include d4 to calculate the anharmonic phonon. I am just wondering if you could tell me how many configurations you were using to make d4 converged. You previously mentioned 1000 configurations may be not sufficient to make d4 converged. The number of configurations may be different for various systems but it is good for me to know the order of magnitude given these are both hailde perovskite materials.

Best wishes,
Kang

[mesonepigreco]

In my case, varying the temperature, I spanned between 800 to 4000, so I would expect your case to be on the same order of magnitude.
Anyway, you can test the convergence a posteriori (it is a good practice anyway to get an idea on the stochastic error), as reported in this FAQ: How do I understand if the free energy Hessian calculation is converged?