The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects

In the simplest and most standard usage, the SSCHA free energy functional that is minimized depends on the centroid positions $\boldsymbol{\mathcal{R}}$ and the auxiliary FCs Φ as

$$\begin{equation}\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]={\left\langle K+V(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}-T{S}_{\text{ion}}\left[{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}\right].\end{equation} \tag{ 7 }$$

Here, we explicit that the entropy S_ion only accounts for ionic degrees of freedom (not electronic). After the SSCHA minimization, the final estimate of the equilibrium free energy is given by

$$\begin{equation}F=\underset{\boldsymbol{\mathcal{R}},\enspace \mathbf{\Phi }}{\mathrm{min}}\enspace \mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]=\mathcal{F}[{\boldsymbol{\mathcal{R}}}_{\text{eq}},{\mathbf{\Phi }}_{\text{eq}}].\end{equation} \tag{ 8 }$$

Therefore, the final result of a SSCHA free energy calculation is given in terms of the equilibrium configuration ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$, the free energy F, and the SSCHA auxiliary FCs Φ_eq. The final free energy accounts for quantum and thermal ionic fluctuations without approximating the BO energy surface: it is valid to study thermodynamic properties. The ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$ positions determine the average atomic positions also taken into account quantum/thermal fluctuations and anharmonicity. It is important to remark, however, that the Gaussian variance Φ has, in principle, no relation with the experimentally observed phonon frequencies, as it is just a variable parameterizing the density matrix. The relation of it with the physical phonon frequencies is discussed in section 4.3. For clarity, we report in table 1 a collection of symbols used in the equations, with their meanings and the reference on the equation where they have been introduced.

Table 1. Collection of some symbols frequently used in the main text. First column, the symbol used. Second column, a short description. Third column, equation or page-column of the first occurrence.

Symbol	Meaning	First use
R	Atomic position (canonical variable)	equation (3)
V( R )	Potential energy	equation (3)
$\boldsymbol{\mathcal{R}}$	Trial centroid positions (parameter)	equation (4)
u	Displacement from the average atomic position $\boldsymbol{\mathcal{R}}$	equation (15)
Φ	Trial harmonic matrix (parameter)	equation (4)
Ψ	Static displacement–displacement correlation matrix relative to Φ	equation (15)
${\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$	Trial harmonic Hamiltonian for given $\boldsymbol{\mathcal{R}}$, Φ	equation (4)
${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$	Density matrix of ${\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$ (trial density matrix)	Page 3-1
${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}(\boldsymbol{R})$	Gaussian positional probability density for ${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$	Page 3-1
$\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]$	SSCHA Helmholtz free energy functional	equation (7)
$\mathcal{G}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]$	SSCHA Gibbs free energy functional	equation (10)
${\mathbf{f}}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})$	Forces for the Hamiltonian ${\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$ acting on the ions when they are in the R positions	equation (12)
f(BO)( R )	Born–Oppenheimer forces acting on the ions when they are in the R positions	equation (12)
P (BO)( R )	Born–Oppenheimer stress tensor when the ions are in the R positions	equation (19)
${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$	2nd SSCHA force constants for a given $\boldsymbol{\mathcal{R}}$, it is the trial Φ that minimizes $\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]$	Page 12-2
${\boldsymbol{D}}_{\boldsymbol{\mathcal{R}}}$	${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$ divided by the square root of the masses	equation (52)
${(3){\mathbf{\Phi }}}_{\boldsymbol{\mathcal{R}}},\enspace {(4){\mathbf{\Phi }}}_{\boldsymbol{\mathcal{R}}}$	3rd and 4th order SSCHA force constants for a given $\boldsymbol{\mathcal{R}}$	equations (53) and (54)
${(3){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}},{(4){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}}$	${(3){\mathbf{\Phi }}}_{\boldsymbol{\mathcal{R}}},{(4){\mathbf{\Phi }}}_{\boldsymbol{\mathcal{R}}}$ divided by the square root of the masses	equations (53) and (54)
$F(\boldsymbol{\mathcal{R}})$	SSCHA positional Helmholtz free energy, given by $\mathcal{F}[\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}]$	equation (50)
${\boldsymbol{\mathcal{R}}}_{\text{eq}}$	SSCHA equilibrium centroids, trial centroids that minimizes $F(\boldsymbol{\mathcal{R}})$	equation (8)
Φeq	SSCHA harmonic matrix ${\mathbf{\Phi }}_{{\boldsymbol{\mathcal{R}}}_{\text{eq}}}$, the trial Φ at the minimum of the free energy functional	equation (8)
${\mathcal{H}}^{(\text{S})}$	SSCHA effective harmonic Hamiltonian, given by ${\mathcal{H}}_{{\boldsymbol{\mathcal{R}}}_{\text{eq}},{\mathbf{\Phi }}_{\text{eq}}}$	equation (60)
D (S)	Dynamical matrix of ${\mathcal{H}}^{(\text{S})}$, given by ${\boldsymbol{D}}_{{\boldsymbol{\mathcal{R}}}_{\text{eq}}}$	Page 14-1
${(3){\boldsymbol{D}}}_{\text{eq}},\enspace {(4){\boldsymbol{D}}}_{\text{eq}}$	Symbols indicating ${(3){\boldsymbol{D}}}_{{\boldsymbol{\mathcal{R}}}_{\text{eq}}}$ and ${(4){\boldsymbol{D}}}_{{\boldsymbol{\mathcal{R}}}_{\text{eq}}}$, respectively	equation (62)
D (F)	Positional Helmholtz free energy Hessian divided by the square root of the masses	Page 12-2
G (z)	One-phonon Green function	equation (67)
Π(0), Π(z)	Static and dynamic SSCHA self-energy	equations (62) and (68)
$(\text{B}){\mathbf{\Pi }}(0),\enspace (\text{B}){\mathbf{\Pi }}(z)$	Static and dynamic SSCHA bubble self-energy	equations (64) and (69)
σ( q , Ω)	Phonon spectral function (with the reciprocal lattice vector made explicit)	equation (70)
ωμ ( q )	Frequency of the (μ, q ) SSCHA auxiliary phonon	equation (16)
Ωμ ( q )	Frequency of the (μ, q ) static approximation phonon from D (F)	Page 15-1
Ωμ ( q ), Γμ ( q )	Frequency and linewidth of the (μ, q ) anharmonic phonon in the Lorentzian approximation	equation (81)

The SSCHA can also perform the free energy minimization at fixed pressure instead. In this case, the Gibbs–Bogoliubov inequality is satisfied by the Gibbs free energy G, defined as

$$\begin{equation}G=F+{P}^{{\ast}}{{\Omega}}_{\text{Vol}},\end{equation} \tag{ 9 }$$

where P* is the target pressure, Ω_Vol is the simulation box volume, and F is the Helmholtz free energy. In this case, the code optimizes

$$\begin{equation}G{\leqslant}\mathcal{G}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]=\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]+{P}^{{\ast}}{{\Omega}}_{\text{Vol}},\end{equation} \tag{ 10 }$$

which can be used, for instance, to estimate the structural changes imposed by pressure by fully accounting for fluctuations.

As made explicit in equation (7), only thermal effects on the ions are taken into account so far, whereas the electrons are considered at zero temperature. However, at very high temperatures the entropy associated to electrons may be important. Within the SSCHA, it is possible to explicitly include finite-temperature effects on the electrons too. The key is to replace in equation (7) the electronic ground state energy V( R ) with the finite-temperature electronic free energy F_el( R ) = E_el( R ) − TS_el( R ) (if electrons have finite temperature, in the adiabatic approximation forces and equilibrium position of the ions are ruled by the electronic free energy). In this case the SSCHA method minimizes the functional

$$\begin{align}\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]& ={\left\langle K+{F}_{\text{el}}(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}-T\enspace {S}_{\text{ion}}\left[{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}\right]\\ & ={\left\langle K+{E}_{\text{el}}(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}-TS\left[{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}\right],\end{align} \tag{ 11 }$$

where $S\left[{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}\right]={\left\langle {S}_{\text{el}}(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}+{S}_{\text{ion}}\left[{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}\right]$. The same trick can be applied to the Gibbs free energy minimization as well. Therefore, the SSCHA estimation of the system's entropy can also incorporate contributions from both electrons (averaged through the ionic distribution ${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$) and ions. In a DFT framework, for example, this simply comes down to including the electronic temperature in the energy/forces/stress calculations for the ensemble elements through the Fermi–Dirac occupation of the Kohn–Sham states [41].

In the SSCHA code, the minimization of $\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]$ is performed through a preconditioned gradient descent approach, which requires the calculation of the gradient of the free energy with respect to the centroid positions $\boldsymbol{\mathcal{R}}$ and the auxiliary FCs Φ. The partial derivatives are evaluated through the exact analytic formulas

$$\begin{equation}\frac{\partial \mathcal{F}}{\partial {\mathcal{R}}_{a}}=-{\left\langle {f}_{a}^{(\text{BO})}(\boldsymbol{R})-{f}_{a}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}\end{equation} \tag{ 12 }$$

and

$$\begin{align}\frac{\partial \mathcal{F}}{\partial {{\Phi}}_{cd}}& =\frac{1}{2}\sum\limits _{ab}\frac{\partial {{\Psi}}_{ab}}{\partial {{\Phi}}_{cd}}\left\langle \left({f}_{b}^{(\mathrm{B}\mathrm{O})}(\boldsymbol{R})-{f}_{b}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})\right)\right.\\ & \quad {\times}{\left.\sum\limits _{e}{{\Psi}}_{ae}^{-1}\left({R}_{e}-{\mathcal{R}}_{e}\right)\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}.\end{align} \tag{ 13 }$$

Here f^(BO)( R ) are the BO forces that act on the ions when they are in the R positions; ${\mathbf{f}}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})$ is the force given by the auxiliary harmonic Hamiltonian ${\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$,

$$\begin{equation}{f}_{a}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})=-\sum\limits _{b}\enspace {{\Phi}}_{ab}({R}_{b}-{\mathcal{R}}_{b});\end{equation} \tag{ 14 }$$

and Ψ is the displacement–displacement correlation matrix

$$\begin{equation}{{\Psi}}_{ab}={\left\langle {u}_{a}{u}_{b}\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}},\end{equation} \tag{ 15 }$$

where with $\boldsymbol{u}=\boldsymbol{R}-\boldsymbol{\mathcal{R}}$ we indicate the displacement from the average atomic position. Explicitly,

$$\begin{equation}{{\Psi}}_{ab}=\frac{1}{\sqrt{{M}_{a}{M}_{b}}}\sum\limits _{\mu }\frac{\hslash (2{n}_{\mu }+1)}{2{\omega }_{\mu }}{e}_{\mu }^{a}{e}_{\mu }^{b}.\end{equation} \tag{ 16 }$$

In equation (16), ω_μ and e _μ are the eigenvalues and eigenvectors of the mass rescaled auxiliary FCs ${{\Phi}}_{ab}/\sqrt{{M}_{a}{M}_{b}}$, and n_μ is the Bose–Einstein occupation number for the ω_μ frequency. We underline again here that ω_μ are not the phonon frequencies of the system, but just the frequencies of the auxiliary harmonic Hamiltonian ${\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$. In other words, they are only used to define the trial density matrix ${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$. We show how to compute the physical anharmonic phonon frequencies of the system in section 4.

It is convenient to give an explicit expression for the gradient of the free energy with respect to the auxiliary FCs in terms of the ω_μ eigenvalues and e _μ eigenvectors. As shown in reference [23] (see appendix B), the gradient can be rewritten as

$$\begin{align}\frac{\partial \mathcal{F}}{\partial {{\Phi}}_{cd}}& =\sum\limits _{ab}\frac{{\Lambda}{[0]}^{abcd}}{\sqrt{{M}_{a}{M}_{b}{M}_{c}{M}_{d}}}\left\langle \left({f}_{b}^{(\mathrm{B}\mathrm{O})}(\boldsymbol{R})-{f}_{b}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})\right)\right.\\ & \quad {\times}{\left.\sum\limits _{e}{{\Psi}}_{ae}^{-1}\left({R}_{e}-{\mathcal{R}}_{e}\right)\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}},\end{align} \tag{ 17 }$$

where

$$\begin{align}{\Lambda}{[0]}^{abcd}& =\sum\limits _{\mu \nu }\frac{\hslash }{4{\omega }_{\nu }{\omega }_{\mu }}{e}_{\nu }^{a}{e}_{\mu }^{b}{e}_{\nu }^{c}{e}_{\mu }^{d}\\ & \quad {\times}\begin{cases}\frac{\mathrm{d}{n}_{\mu }}{\mathrm{d}{\omega }_{\mu }}-\frac{2{n}_{\mu }+1}{2{\omega }_{\mu }},\qquad {\omega }_{\nu }={\omega }_{\mu }\\ \frac{{n}_{\mu }-{n}_{\nu }}{{\omega }_{\mu }-{\omega }_{n}u}-\frac{1+{n}_{\mu }+{n}_{\nu }}{{\omega }_{\mu }+{\omega }_{n}u},\quad {\omega }_{\nu }\ne {\omega }_{\mu }\end{cases}.\end{align} \tag{ 18 }$$

Here, ${n}_{\mu }=1/({\text{e}}^{\beta \hslash {\omega }_{\mu }}-1)$. The reason why we have introduced the Λ[0] tensor will be evident in section 4. Even if equation (17) looks different to the gradient introduced in the original SSCHA work in reference [22], it can be demonstrated that both expressions are equivalent by simply playing with the permutation symmetry of ${\left\langle \frac{{\partial }^{2}V}{\partial {R}_{a}\partial {R}_{b}}\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}$. However, equation (18) unambiguously determines the value taken by the Λ[0] tensor for the ω_ν = ω_μ case, while the gradient in reference [22] did not describe explicitly what to do in this degenerate limit.

At the end of the SSCHA optimization, apart from the temperature-dependent ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$ positions and the equilibrium auxiliary FC matrix Φ_eq, the code also calculates the anharmonic stress tensor P , which includes both quantum and thermal ionic fluctuations, as derivatives of the free energy with respect to a strain tensor ɛ :

$$\begin{align}{P}_{\alpha \beta }& =-{\left.\frac{1}{{{\Omega}}_{\text{Vol}}}\frac{\partial \mathcal{F}}{\partial {\varepsilon }_{\alpha \beta }}\right\vert }_{\boldsymbol{\varepsilon }=0}={\left\langle {P}_{\alpha \beta }^{(\text{BO})}(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}\\ & \quad -\frac{1}{2{{\Omega}}_{\text{Vol}}}\sum\limits _{s=1}^{{N}_{\text{a}}}{\left\langle {u}_{s}^{\alpha }{{f}^{(\text{BO})}}_{s}^{\beta }+{u}_{s}^{\beta }{{f}^{(\text{BO})}}_{s}^{\alpha }\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}.\end{align} \tag{ 19 }$$

Here, we have made explicit the atomic index s (lower index) and Cartesian α, β (upper index) of u and f^(BO). P ^(BO)( R ) is the BO stress tensor of the configuration with ions displaced in the R coordinates. This equation is slightly different from the stress tensor equation presented in [24]. The two equations coincide at equilibrium, but this is more general. The derivation of equation (19) is reported in appendix A. Thanks to the temperature-dependent stress, the SSCHA code can optimize also the lattice parameters and the volume. Thus, by relaxing the lattice at different temperatures, we get the thermal expansion straightforwardly.

Remarkably, the stress tensor of equation (19) can be computed with a single SSCHA minimization at fixed volume. This is a huge advantage with respect to the standard quasi-harmonic approximation, not only because it includes quantum and anharmonic effects, but also because it is computationally much more efficient. In fact, the quasiharmonic approximation requires performing harmonic phonon calculations at different volumes (and/or internal lattice positions) to estimate the minimum of the quasi-harmonic free energy with finite differences. This process is extremely cumbersome for crystals with few symmetries and lots of internal degrees of freedom in the structure.

In the current implementation of the code, the symmetries of the space group are imposed a posteriori on the gradients of equations (12) and (13), as well as on (19). This assures that the density matrix satisfies all the symmetries at each step of the minimization. Thus, during the geometry optimization, the system cannot lose any symmetry, though it can gain them. The symmetries are imposed following the methodology explained in appendix D, which is different to the method originally conceived [22]. The current SSCHA code can also work without imposing symmetries, allowing for symmetry loss, though the stochastic number of configurations needed to converge the minimization is larger (see section 3.2.1).

3.2.1. The stochastic sampling

The stochastic nature of the SSCHA comes from the Monte Carlo evaluation of the averages in equations (12), (13), and (19). A set of random ionic configurations are created in a chosen supercell according to the Gaussian ionic probability distribution

$$\begin{align}{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}(\boldsymbol{R})& =\sqrt{\mathrm{det}({\mathbf{\Psi }}^{-1}/2\pi )}\\ & \quad {\times}\enspace \mathrm{exp}\left[-\frac{1}{2}\sum\limits _{ab}({R}_{a}-{\mathcal{R}}_{a}){{{\Psi}}^{-1}}_{ab}({R}_{b}-{\mathcal{R}}_{b})\right].\end{align} \tag{ 20 }$$

The Monte Carlo average of a generic observable O( R ), function only of the ionic position R , is calculated then as weighted sum over the created ensemble:

$$\begin{equation}{\left\langle O(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}=\frac{1}{{\sum }_{j=1}^{{N}_{\text{c}}}\enspace {\rho }_{j}}\sum\limits _{j=1}^{{N}_{\text{c}}}\enspace {\rho }_{j}O({\boldsymbol{R}}_{j}).\end{equation} \tag{ 21 }$$

Here, N_c is the total number of configurations in the ensemble, while R _j is the jth ionic randomly displaced configuration. Each of the R _{j} configurations is generated according to the initial trial ionic distribution ${\tilde {\rho }}_{{\boldsymbol{\mathcal{R}}}^{(0)},{\mathbf{\Phi }}^{(0)}}(\boldsymbol{R})$ from which the minimization starts. To improve the stochastic accuracy, for each R _j configuration also − R _j is created, benefiting from ${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}(\boldsymbol{R})={\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}(-\boldsymbol{R})$ property of the Gaussian distribution.

The ρ_j weights are computed and updated along the free energy minimization as the values of $\boldsymbol{\mathcal{R}}$ and Φ change:

$$\begin{equation}{\rho }_{j}=\frac{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}\left({R}_{j}\right)}{{\tilde {\rho }}_{{\boldsymbol{\mathcal{R}}}^{(0)},{\mathbf{\Phi }}^{(0)}}\left({R}_{j}\right)}.\end{equation} \tag{ 22 }$$

At the beginning, when the ensemble has just been generated and $\boldsymbol{\mathcal{R}}={\boldsymbol{\mathcal{R}}}^{(0)}$ and Φ = Φ⁽⁰⁾, all values of ρ_j = 1. However, as the $\boldsymbol{\mathcal{R}}$ and Φ are updated during the minimization, the weights change. This reweighting technique is commonly used in Monte Carlo methods [42, 43] and takes the name of importance sampling. This allows avoiding generating a new ensemble and computing ab initio energies and forces at each step of the minimization, speeding up the SSCHA calculation.

3.2.2. Minimization algorithm

The minimization strategy implemented in the SSCHA code for the free energy is based on a preconditioned gradient descent. At each step, the $\boldsymbol{\mathcal{R}}$ and Φ are updated as

$$\begin{equation}{\mathbf{\Phi }}^{(n+1)}={\mathbf{\Phi }}^{(n)}-{\lambda }_{\mathbf{\Phi }}\sum\limits _{ab}{\left(\frac{{\partial }^{2}\mathcal{F}}{\partial \mathbf{\Phi }\partial {{\Phi}}_{ab}}\right)}^{-1}\frac{\partial \mathcal{F}}{\partial {{\Phi}}_{ab}}\end{equation} \tag{ 23 }$$

$$\begin{equation}{\boldsymbol{\mathcal{R}}}^{(n+1)}={\boldsymbol{\mathcal{R}}}^{(n)}-{\lambda }_{\boldsymbol{\mathcal{R}}}\sum\limits _{a}{\left(\frac{{\partial }^{2}\mathcal{F}}{\partial \boldsymbol{\mathcal{R}}\partial {\mathcal{R}}_{a}}\right)}^{-1}\frac{\partial \mathcal{F}}{\partial {\mathcal{R}}_{a}}.\end{equation} \tag{ 24 }$$

In a perfectly quadratic landscape, this algorithm assures the convergence in just one step if both λ_Φ and ${\lambda }_{\boldsymbol{\mathcal{R}}}$ are set equal to one. However, in order to avoid too big steps in the minimization, often it is more convenient to chose ${\lambda }_{\mathbf{\Phi }\vert \boldsymbol{\mathcal{R}}}{< }1$. This algorithm, with the Hessian matrix that multiplies the gradient, is the preconditioned steepest descent. If the preconditioning option is set to false, a standard steepest descent minimization is followed instead, with λ_Φ and ${\lambda }_{\boldsymbol{\mathcal{R}}}$ re-scaled to the maximum eigenvalue of the $\frac{{\partial }^{2}\mathcal{F}}{\partial {\mathbf{\Phi }}^{2}}$ and $\frac{{\partial }^{2}\mathcal{F}}{\partial {\boldsymbol{\mathcal{R}}}^{2}}$ Hessian matrices, respectively, in order to have a dimensional values independent on the system.

The preconditioning Hessian matrices that multiplies the gradients in equation (27) and (28) are approximated by the code. Following the procedure introduced in reference [24], we use the exact Hessian in the minimum of a perfectly harmonic oscillator with the same frequencies as the SCHA auxiliary Hamiltonian. In particular, they are:

$$\begin{equation}\frac{{\partial }^{2}\mathcal{F}}{\partial {{\Phi}}_{ab}\partial {{\Phi}}_{cd}}\approx \frac{1}{2}\frac{\partial {{\Psi}}_{ab}}{\partial {{\Phi}}_{cd}}\end{equation} \tag{ 25 }$$

and

$$\begin{equation}\frac{{\partial }^{2}\mathcal{F}}{\partial \boldsymbol{\mathcal{R}}\partial \boldsymbol{\mathcal{R}}}\approx \mathbf{\Phi }.\end{equation} \tag{ 26 }$$

Equation (25) is presented differently from the original work in which it was derived [24]. We prove in appendix B that they are exactly the same. Considering that the Hessian preconditioner cancels out the $\frac{1}{2}\frac{\partial {{\Psi}}_{ab}}{\partial {{\Phi}}_{cd}}$ term in equation (13), the resulting update of the variational parameters at each step in the minimization is performed as

$$\begin{align}{{\Phi}}_{ab}^{(n+1)}& ={{\Phi}}_{ab}^{(n)}+-{\lambda }_{\mathbf{\Phi }}\left\langle \left({f}_{b}^{(\text{BO})}(\boldsymbol{R})-{f}_{b}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})\right)\right.\\ & \quad {\times}{\left.\sum\limits _{c}\enspace {{{\Psi}}^{-1}}_{ac}\left({R}_{c}-{\mathcal{R}}_{c}\right)\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}\end{align} \tag{ 27 }$$

and

$$\begin{equation}{\mathcal{R}}_{a}^{(n+1)}={\mathcal{R}}_{a}^{(n)}+{\lambda }_{\boldsymbol{\mathcal{R}}}\sum\limits _{b}\enspace {{\Phi}}_{ab}^{-1}{\left\langle {f}_{b}^{(\text{BO})}(\boldsymbol{R})-{f}_{b}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}(\boldsymbol{R})\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}.\end{equation} \tag{ 28 }$$

This implementation is very efficient, especially for equation (27), as there is no need to calculate the ${{\Lambda}}_{\boldsymbol{\mathcal{R}}}[0]$ tensor. Therefore, computing directly equation (27) is much faster than calculating the gradient of equation (13).

The code allows the user to select a different minimization algorithm specifically for the minimization with respect to Φ: the root representation. Since the minimization with respect to Φ is the more challenging, this technique aims to further improving the Φ optimization. In particular, the gradient has a stochastic error and the minimization is performed with a finite step size. For these reasons, Φ could become non positive definite during the optimization (i.e. the dynamical matrix has imaginary frequencies). If this occurs, the minimization is halted raising an error, as the density matrix of equation (20) diverges. In such a case, the minimization must be manually restarted, either by taking a smaller step or by stopping the minimization before reaching imaginary frequencies (fixing the maximum number of steps). This kind of halts do not occur often when using the preconditioning in the minimization. However, they may be encountered if few configurations are generated for each ensemble or the starting dynamical matrix is very far from equilibrium.

To solve these problems, we implement the root representation, in which, instead of updating Φ as in equation (27), it updates a root of Φ:

$$\begin{equation}{\sqrt[n]{\mathbf{\Phi }}}^{(i+1)}={\sqrt[n]{\mathbf{\Phi }}}^{(i)}-{\lambda }_{\mathbf{\Phi }}{\boldsymbol{\mathcal{G}}}_{\boldsymbol{n}}.\end{equation} \tag{ 29 }$$

The updating direction ${\boldsymbol{\mathcal{G}}}_{\boldsymbol{n}}$ depends on the root order n:

$$\begin{equation}{\boldsymbol{\mathcal{G}}}_{\boldsymbol{2}}=\sqrt{\mathbf{\Phi }}\cdot \frac{\partial \mathcal{F}}{\partial \mathbf{\Phi }}+\frac{\partial \mathcal{F}}{\partial \mathbf{\Phi }}\cdot \sqrt{\mathbf{\Phi }},\end{equation} \tag{ 30 }$$

where with the ⋅ we indicate a matrix product. Similarly,

$$\begin{equation}{\boldsymbol{\mathcal{G}}}_{\boldsymbol{4}}=\sqrt[4]{\mathbf{\Phi }}\cdot {\boldsymbol{\mathcal{G}}}_{\boldsymbol{2}}+{\boldsymbol{\mathcal{G}}}_{\boldsymbol{2}}\cdot \sqrt[4]{\mathbf{\Phi }}.\end{equation} \tag{ 31 }$$

We select the positive definite root matrix. Indeed, after the step of equation (29) the original FC matrix is obtained as

$$\begin{equation}{\mathbf{\Phi }}^{(i+1)}={\left({\sqrt[n]{\mathbf{\Phi }}}^{(i+1)}\right)}^{n}.\end{equation} \tag{ 32 }$$

Thanks to the definition in equation (32), the dynamical matrix is always positive definite for any even value of n.

The root representation is independent of the preconditioning. With preconditioning, we replace the free energy gradient in equation (30) with the preconditioned direction in equation (27) (the gradient multiplied by the approximated Hessian). This is different from what was proposed in the original work [24], where the Hessian matrix was computed also for the $\sqrt{\mathbf{\Phi }}$ and $\sqrt[4]{\mathbf{\Phi }}$ cases. However, we noticed that in systems with many atoms, the Hessian matrix calculation becomes the bottleneck as it scales with ${N}_{\text{a}}^{6}$. The implementation here described allows for a much faster Φ update and avoids calculating the Hessian matrix. The drawback is that the optimization step is not as optimal as it would be if the proposal in reference [24] was followed. The code offers six combinations for the minimization procedure: no root, square root (n = 2), and fourth-root (n = 4), all of them with or without the preconditioned direction. The optimal minimization step is n = 1 with preconditioning. If the square root is employed (n = 2), it is preferable to use the preconditioning. If fourth-root is employed (n = 4), the best performances are without preconditioning.

3.2.3. The lattice geometry optimization

The lattice degrees of freedom { a _i } are relaxed only after the minimization of the free energy with respect to $\boldsymbol{\mathcal{R}}$ and Φ at a constant volume stops (see section 3.2.4 for a detailed description of the stopping criteria). For this reason, the lattice geometry optimization is an 'outer' optimization: at each step of the lattice geometry optimization, we perform a full free energy minimization with respect to the centroids $\boldsymbol{\mathcal{R}}$ and auxiliary FCs Φ. This means that each step of the lattice geometry optimization is performed with a different ensemble, whose configurations are all generated with the same lattice vectors.

To update the lattice, the code calculates the stress tensor with equation (19), and generates a strain for the lattice as

$$\begin{equation}{\varepsilon }_{\alpha \beta }={{\Omega}}_{\text{Vol}}\left({P}_{\alpha \beta }-{P}^{{\ast}}{\delta }_{\alpha \beta }\right),\end{equation} \tag{ 33 }$$

where P* is the target pressure of the relaxation and δ_αβ is the Kronecker delta. The lattice parameters { a _i } are updated as

$$\begin{equation}{{a}_{i}^{\prime }}_{\alpha }={{a}_{i}}_{\alpha }+{\lambda }_{\left\{{\boldsymbol{a}}_{i}\right\}}\sum\limits _{\beta }\enspace {\varepsilon }_{\alpha \beta }{{a}_{i}}_{\beta },\end{equation} \tag{ 34 }$$

where ${\lambda }_{\left\{{\boldsymbol{a}}_{i}\right\}}$ is the update step. Since each step requires a new ensemble, it is crucial to reduce the number of steps to reach convergence by properly picking the right value for ${\lambda }_{\left\{{\boldsymbol{a}}_{i}\right\}}$. In an isotropic material with a constant bulk modulus

$$\begin{equation}{B}_{0}={{\Omega}}_{\text{Vol}}{\left.\frac{{\mathrm{d}}^{2}{V}^{(\mathrm{B}\mathrm{O})}}{\mathrm{d}{{\Omega}}_{\text{Vol}}^{2}}\right\vert }_{\boldsymbol{R}=\boldsymbol{\mathcal{R}}},\end{equation} \tag{ 35 }$$

the optimal value of the step is

$$\begin{equation}{\lambda }_{\left\{{\boldsymbol{a}}_{i}\right\}}=\frac{1}{3{{\Omega}}_{\text{Vol}}{B}_{0}}.\end{equation} \tag{ 36 }$$

B₀ is an input parameter given in GPa units. Good values of B₀ may range from 10 GPa for crystals at ambient conditions, like ice, up to 800 GPa in systems at Mbar pressures (or for diamond). Remember that increasing the value of B₀ produces smaller steps in the cell parameters. The user can estimate the optimal value of B₀ to assure the fastest convergence by manually computing it from equation (35), by taking finite differences of the pressure obtained at two uniformly strained volumes, or by looking for the experimental value of similar compounds.

Alternatively to the fixed pressure optimization, it is also possible to perform the geometry lattice optimization at fixed volume. In this case P* is recomputed at each step so that $\mathrm{T}\mathrm{r}\left[\boldsymbol{\varepsilon }\right]=0$. In this case, the final lattice parameters are also rescaled so that the final volume matches the one before the step. Since this algorithm has one less degree of freedom than the fixed pressure one, it usually converges faster.

3.2.4. The code flowchart

To start a SSCHA simulation, we need a starting guess on the trial positive definite FCs matrix Φ⁽⁰⁾ and on the average atomic positions ${\boldsymbol{\mathcal{R}}}^{(0)}$. Even if in principle the starting point is arbitrary, the closer to the solution we begin, the faster the minimization converges. Thus, the coordinates at the minimum of the BO energy landscape and the harmonic FCs are usually good starting points, which can be obtained from any code that computes phonons. The supercell of the simulation is given by the dimension of the input FCs matrix, while the centroids R are defined in the unit cell (they satisfy translational symmetry). If the original dynamical matrix contains imaginary frequencies, it can be reverted to positive definite as

$$\begin{equation}{{\Phi}}_{ab}^{(0)}=\sqrt{{M}_{a}{M}_{b}}\sum\limits _{\mu }\vert {\omega }_{\mu }^{2}\vert {e}_{\mu }^{a}{e}_{\mu }^{b}.\end{equation} \tag{ 37 }$$

Then, the first random ensemble (that we call population in the SSCHA language) can be generated. For each configuration inside the population, its total BO energy as well as its classical atomic forces f^(BO) and stress tensor P ^(BO) must be computed. This is done with an external code, either manually (by computing externally the energies, forces, and stress tensor, and loading them back into the SSCHA code), or automatically (with an appropriate configuration discussed in section 5). Once the BO energies, forces, and stress tensor of all the configurations have been computed, the minimization starts. The gradients of the free energy are computed as described in section 3.2.2 and the minimization continues either until the stochastic sampling is not good or the algorithm converges.

If $\boldsymbol{\mathcal{R}}$ and Φ change a lot during the minimization of the free energy, the original ensemble no longer describes well the new probability distribution ${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}(\boldsymbol{R})$, and the stochastic error increases. This occurrence is automatically checked by the SSCHA code calculating the Kong–Liu [44, 45] effective sample size N_eff:

$$\begin{equation}{N}_{\text{eff}}=\frac{{\sum }_{j=1}^{{N}_{\text{c}}}\enspace {\rho }_{j}^{2}}{{\left({\sum }_{j=1}^{{N}_{\text{c}}}\enspace {\rho }_{j}\right)}^{2}}.\end{equation} \tag{ 38 }$$

We halt the minimization when the ratio between N_eff and the number of configurations N_c is lower than a parameter η defined by the user:

$$\begin{equation}\frac{{N}_{\text{eff}}}{{N}_{\text{c}}}{< }\eta .\end{equation} \tag{ 39 }$$

A standard value of η that ensures a correct minimization is 0.5, but it can be convenient to lower it a bit to accelerate convergence in the first steps.

The convergence, on the contrary, is achieved only if the two gradients with respect to $\boldsymbol{\mathcal{R}}$ and Φ are lower than a given threshold:

$$\begin{equation}\left\vert \frac{\partial \mathcal{F}}{\partial \mathbf{\Phi }}\right\vert {< }{\delta }_{\mathbf{\Phi }}\end{equation} \tag{ 40 }$$

$$\begin{equation}\left\vert \frac{\partial \mathcal{F}}{\partial \boldsymbol{\mathcal{R}}}\right\vert {< }{\delta }_{\boldsymbol{\mathcal{R}}}.\end{equation} \tag{ 41 }$$

The δ threshold is provided by the user and re-scaled at each step by the estimation of the stochastic error on the corresponding gradient (meaningful_factor). So, at each step, δ is

$$\begin{equation}{\delta }_{\mathbf{\Phi }}=\text{meaningful}_\text{factor}\cdot \left\vert {\Delta}\frac{\partial \mathcal{F}}{\partial \mathbf{\Phi }}\right\vert \end{equation} \tag{ 42 }$$

$$\begin{equation}{\delta }_{\boldsymbol{\mathcal{R}}}=\text{meaningful}_\text{factor}\cdot \left\vert {\Delta}\frac{\partial \mathcal{F}}{\partial \boldsymbol{\mathcal{R}}}\right\vert .\end{equation} \tag{ 43 }$$

In this way, the user-provided variable meaningful_factor is independent on the system size or the number of configurations used.

If the lattice parameters are free to move, then an additional condition must be fulfilled in order to end the minimization: each component of the strain per unit-cell volume Ω_Vol must be smaller than the stochastic error on the stress tensor:

$$\begin{equation}\frac{{\varepsilon }_{\alpha \beta }}{{{\Omega}}_{\text{Vol}}}{\leqslant}{\Delta}{P}_{\alpha \beta }.\end{equation} \tag{ 44 }$$

If equation (44) is not fulfilled, even if all the gradients are lower than the chosen threshold, the code generates a new ensemble and continues (until both conditions are satisfied).

At the end of the minimization, the output of the SSCHA minimization gives the total free energy (with stochastic error), the average equilibrium ionic positions ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$, the equilibrium auxiliary FC matrix Φ_eq, and the stress tensor P . All output quantities are temperature-dependent, and include quantum-thermal fluctuations and anharmonicity. A flowchart that represents the whole execution of a SSCHA run is presented in figure 2. If the SSCHA code is coupled with an ab initio total-energy engine, the most expensive calculation in the flowchart is by far the calculation of BO energy, forces, and stress tensors on the whole ensemble, which may contain up to several hundreds or thousands of configurations. For this reason, the pretty complex workflow we set up is aimed to pass by the calculation of a new ensemble as few times as possible. Most materials studied and presented in section 7 are converged within 3 populations, and the CPU time required to minimize each population is few minutes on a single CPU of modern laptops.

Zoom In Zoom Out Reset image size

The preconditioned gradient descent approach sketched above offers a very efficient implementation of the SSCHA theory, in which the anharmonic free energy is optimized by all degrees of freedom in the crystal structure, including internal coordinates as well as lattice vectors. If the centroid positions $\boldsymbol{\mathcal{R}}$ are kept fixed in the minimization, the SSCHA self-consistent equation

$$\begin{equation}{{\Phi}}_{ab}(\boldsymbol{\mathcal{R}})={\left\langle \frac{{\partial }^{2}V}{\partial {R}_{a}\partial {R}_{b}}\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }(\boldsymbol{\mathcal{R}})}}\end{equation} \tag{ 45 }$$

offers an alternative way of implementing the SSCHA theory (see reference [23] for a proof of equation (45)). It is important to underline the self-consistent condition required by the equation above, as the quantum statistical average is taken with a density matrix dependent on $\mathbf{\Phi }(\boldsymbol{\mathcal{R}})$, which must equal the result of the average. As in this approach the centroid positions are not optimized, the obtained auxiliary FC matrix depends parametrically on $\boldsymbol{\mathcal{R}}$.

The self-consistent equation can be implemented stochastically, following the procedure outlined in section 3.2.1. By using integration by parts [23], the right-hand-side of equation (45) can be rewritten in terms of forces and displacements. Thus, with the importance sampling technique and reweighting, the equation can be solved by calculating forces in supercells generated with the SSCHA density matrix. An equivalent approach [30] is to extract the auxiliary FCs by fitting the obtained forces in the supercells generated with the SSCHA density matrix to equation (14). This approach has been followed recently [30, 46, 47], where a least-squares technique is followed for the fitting.

The use of the self-consistent equation is valid, thus, only for fixed centroid positions. If the centroid positions want to be optimized as well within this approach, the self-consistent procedure should be repeated for different values of $\boldsymbol{\mathcal{R}}$, calculate the free energy for these positions, and see where its minimum is. Clearly this is a very cumbersome procedure unless centroid positions are fixed by symmetry. Moreover, solving the self-consistent equation fixing the centroid positions at the classical R ₀ positions, which it is usually the case [30, 46, 47], neglects all the effects of quantum and thermal fluctuations on the structure. Since within our approach based on the gradient descent we can optimize the free energy not only with respect to the auxiliary FCs but also all degrees of freedom in the crystal structure, the workflow outlined in section 3.2 provides a full picture of the effect of quantum/thermal fluctuations as well as anharmonicity on crystals, much more efficient than the approaches based on equation (45).

As shown in section 2, for a given temperature the free energy at equilibrium F of a system with Hamiltonian H is obtained by minimizing the density-matrix functional $\mathcal{F}[\tilde {\rho }]={\left\langle K+V(\boldsymbol{R})\right\rangle }_{\tilde {\rho }}-TS[\tilde {\rho }]$:

$$\begin{equation}F=\underset{\tilde {\rho }}{\mathrm{min}}\enspace \mathcal{F}[\tilde {\rho }]=\mathcal{F}[\rho ],\end{equation} \tag{ 46 }$$

where ρ is the equilibrium density matrix of the system obtained at the minimum. The average atomic positions at equilibrium are ${\left\langle \boldsymbol{R}\right\rangle }_{\rho }={\boldsymbol{\mathcal{R}}}_{\text{eq}}$. By minimizing the functional keeping fixed the average atomic positions in a generic configuration $\boldsymbol{\mathcal{R}}$, ${\left\langle \boldsymbol{R}\right\rangle }_{\tilde {\rho }}=\boldsymbol{\mathcal{R}}$, we define the positional free energy $F(\boldsymbol{\mathcal{R}})$:

$$\begin{equation}F(\boldsymbol{\mathcal{R}})=\underset{\begin{subarray}{c}\tilde {\rho }\\ {\left\langle \boldsymbol{R}\right\rangle }_{\tilde {\rho }}=\boldsymbol{\mathcal{R}}\end{subarray}}{\mathrm{min}}\mathcal{F}[\tilde {\rho }]=\mathcal{F}[{\rho }_{\boldsymbol{\mathcal{R}}}],\end{equation} \tag{ 47 }$$

where ${\rho }_{\boldsymbol{\mathcal{R}}}$ is the density matrix giving the constrained minimum for the considered average position $\boldsymbol{\mathcal{R}}$. Since

$$\begin{equation}F=F({\boldsymbol{\mathcal{R}}}_{\text{eq}})=\underset{\boldsymbol{\mathcal{R}}}{\mathrm{min}}\enspace F(\boldsymbol{\mathcal{R}}),\end{equation} \tag{ 48 }$$

$\boldsymbol{\mathcal{R}}$ and $F(\boldsymbol{\mathcal{R}})$ can be interpreted as a multidimensional order parameter and a thermodynamic potential, respectively, in the study of displacive phase transitions according to Landau's theory. Properly speaking, the 'order' parameter would be $\boldsymbol{\mathcal{R}}-{\boldsymbol{\mathcal{R}}}_{\text{hs}}$, where ${\boldsymbol{\mathcal{R}}}_{\text{hs}}$ is the average position of the atoms when the system is in the high-symmetry phase. Therefore, the knowledge of the positional free energy landscape as a function of external parameters, like temperature or pressure, gives crucial information about the structural stability and evolution of a system, as it allows to determine the (meta-)stable configurations corresponding to (local) minima of the positional free energy.

The Hessian of the positional free energy $F(\boldsymbol{\mathcal{R}})$ in the equilibrium configuration is the inverse response function to a static perturbation on the nuclei (i.e. the inverse of the static susceptibility). In presence of a second order phase transition the static response function diverges, which results in one or more eigenvalues of the positional free energy Hessian going to zero. This means that the occurrence of displacive second-order phase transitions, like CDW or ferroelectric transitions [41, 49–52], can be characterized by analyzing the evolution with temperature of the eigenvalues of the equilibrium positional free energy Hessian. Typically, in these cases at high temperature the minimum point of the free energy, i.e. the equilibrium configuration ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$, is a high-symmetry configuration ${\boldsymbol{\mathcal{R}}}_{\text{hs}}$. Therefore, at high temperature the free energy Hessian in ${\boldsymbol{\mathcal{R}}}_{\text{hs}}$ is positive definite, i.e. its eigenvalues are positive. As the temperature decreases, the minimum in ${\boldsymbol{\mathcal{R}}}_{\text{hs}}$ becomes less and less deep, until it becomes a saddle point at the transition temperature (i.e. at least one eigenvalue is zero), so that a second-order displacive phase transition occurs and the equilibrium configuration ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$ moves towards lower-symmetry configurations that reduce the free energy as the temperature decreases further (following the pattern indicated by the eigenvector of the vanishing eigenvalue). Using the same approach, it is possible to characterize second-order displacive phase transitions driven by other external parameters, like the pressure in high-pressure superconducting hydrides [21, 34, 36, 48, 53].

The role played by the eigenvalues and eigenvectors of the positional free energy Hessian in tracing the system's structural stability recalls the role played by the harmonic dynamical matrix in the standard harmonic approximation, but now including lattice thermal and quantum anharmonic effects in the dynamics of the nuclei. Therefore, the Hessian of the positional free energy, divided by the masses, ${D}_{ab}^{(\text{F})}={\left.{\partial }^{2}F/\partial {\mathcal{R}}^{a}\partial {\mathcal{R}}^{b}\right\vert }_{{\boldsymbol{\mathcal{R}}}_{\text{eq}}}/\sqrt{{M}_{a}{M}_{b}}$, can be considered a natural generalization of the harmonic dynamical matrix that, however, includes thermal and quantum effects.

What explained hitherto about the role played by the positional free energy and its Hessian is general. In particular, the evaluation of the positional free energy within the SSCHA is pretty straightforward. Indeed, the average position for a trial harmonic density matrix ${\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}$ coincides with the centroid parameter $\boldsymbol{\mathcal{R}}$,

$$\begin{equation}{\left\langle \boldsymbol{R}\right\rangle }_{{\tilde {\rho }}_{\boldsymbol{\mathcal{R}},\mathbf{\Phi }}}=\boldsymbol{\mathcal{R}}.\end{equation} \tag{ 49 }$$

Thus, within the SSCHA the positional free energy is obtained by minimizing the SSCHA free energy functional $\mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }]$ with respect to the trial quadratic amplitude Φ only:

$$\begin{equation}F(\boldsymbol{\mathcal{R}})=\underset{\mathbf{\Phi }}{\mathrm{min}}\enspace \mathcal{F}[\boldsymbol{\mathcal{R}},\mathbf{\Phi }].\end{equation} \tag{ 50 }$$

The auxiliary FCs that minimize equation (50) for a given $\boldsymbol{\mathcal{R}}$ position of the centroids will be labeled in the following as ${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$. Solving equation (50) allows to employ the SSCHA code to have direct access to $F(\boldsymbol{\mathcal{R}})$ for any $\boldsymbol{\mathcal{R}}$ and, in principle, to compute the Hessian by finite differences. However, as discussed above, such a finite-difference approach would be extremely expensive for two main reasons. First, it would require a large number of configurations in the ensemble to reduce the stochastic error and calculate the derivatives by finite differences. Second, because the large number of degrees of freedom in $\boldsymbol{\mathcal{R}}$ prevents any realistic finite-difference approach. Luckily the SSCHA code allows to avoid any cumbersome finite-difference approach by exploiting an analytic formula for the positional free energy Hessian.

Before describing the analytic formula, let us introduce a notation that will simplify the mathematical expressions. Given two tensors X and Y , with the single dot product X ⋅ Y we will indicate the contraction of the last index of X with the first index of Y , ∑_c X_...c Y_c.... Likewise, with the double-dot product X : Y we will indicate the contraction of the last two indices of X with the first two indices of Y , ∑_cd X_...cd Y_cd.... Moreover, any fourth-order tensor X_pqlm can be interpreted as a 'super' matrix X_AB , with the composite indices A = (pq) and B = (lm), and vice versa (through this correspondence we can define, for example, the inverse of a fourth-order tensor and the identity fourth-order tensor $\mathbb{1}$). Using this notation, we can express the positional free energy Hessian, $1/\sqrt{{M}_{a}{M}_{b}}\enspace \enspace {\partial }^{2}F/\partial {\mathcal{R}}^{a}\partial {\mathcal{R}}^{b}$, in component-free form as

$$\begin{align}& \frac{1}{\sqrt{\boldsymbol{M}}}\cdot \frac{{\partial }^{2}F}{\partial \boldsymbol{\mathcal{R}}\partial \boldsymbol{\mathcal{R}}}\cdot \frac{1}{\sqrt{\boldsymbol{M}}}\\ & ={\boldsymbol{D}}_{\boldsymbol{\mathcal{R}}}+{(3){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[0]\enspace \text{:}\enspace {\left[\mathbb{1}-{(4){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[0]\right]}^{-1}\enspace \text{:}\enspace {(3){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}},\end{align} \tag{ 51 }$$

where M_ab = δ_ab M_a is the mass matrix,

$$\begin{align}{({D}_{\boldsymbol{\mathcal{R}}})}_{ab}& =\frac{1}{\sqrt{{M}_{a}{M}_{b}}}{\left\langle \frac{{\partial }^{2}V}{\partial {R}^{a}\partial {R}^{b}}\right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}\\ & =\frac{{({{\Phi}}_{\boldsymbol{\mathcal{R}}})}_{ab}}{\sqrt{{M}_{a}{M}_{b}}},\end{align} \tag{ 52 }$$

$$\begin{align}{({(3){D}}_{\boldsymbol{\mathcal{R}}})}_{abc}& =\frac{1}{\sqrt{{M}_{a}{M}_{b}{M}_{c}}}{\left\langle \frac{{\partial }^{3}V}{\partial {R}^{a}\partial {R}^{b}\partial {R}^{c}}\right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}\\ & =\frac{{({(3){{\Phi}}}_{\boldsymbol{\mathcal{R}}})}_{abc}}{\sqrt{{M}_{a}{M}_{b}{M}_{c}}},\end{align} \tag{ 53 }$$

$$\begin{align}{({(4){D}}_{\boldsymbol{\mathcal{R}}})}_{abcd}& =\frac{1}{\sqrt{{M}_{a}{M}_{b}{M}_{c}{M}_{d}}}{\left\langle \frac{{\partial }^{4}V}{\partial {R}^{a}\partial {R}^{b}\partial {R}^{c}\partial {R}^{d}}\right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}\\ & =\frac{{({(4){{\Phi}}}_{\boldsymbol{\mathcal{R}}})}_{abcd}}{\sqrt{{M}_{a}{M}_{b}{M}_{c}{M}_{d}}},\end{align} \tag{ 54 }$$

and ${\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[0]$ is the z = 0 value of the fourth-order tensor ${\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[z]$, already introduced in equation (18). In the equations above the quantum statistical averages are taken with ${\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}$, which for a given $\boldsymbol{\mathcal{R}}$ position of the centroids is taken with the ${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$ auxiliary FCs that minimize the free energy. ${\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[z]$ is given in components by

$$\begin{equation}{\left({{\Lambda}}_{\boldsymbol{\mathcal{R}}}[z]\right)}^{abcd}=\sum\limits _{\mu \nu }\enspace \mathcal{F}(z,{\omega }_{\mu },{\omega }_{\nu })\enspace {e}_{\nu }^{a}{e}_{\mu }^{b}{e}_{\nu }^{c}{e}_{\mu }^{d},\end{equation} \tag{ 55 }$$

where ${\omega }_{\mu }^{2}$ and ${e}_{\nu }^{a}$ are the eigenvalues and eigenvectors of ${\boldsymbol{D}}_{\boldsymbol{\mathcal{R}}}$, and

$$\begin{align}\mathcal{F}(z,{\omega }_{\nu },{\omega }_{\mu })& =\frac{\hslash }{4{\omega }_{\mu }{\omega }_{\nu }}\left[\frac{({\omega }_{\mu }-{\omega }_{\nu })({n}_{\mu }-{n}_{\nu })}{{({\omega }_{\mu }-{\omega }_{\nu })}^{2}-{z}^{2}}+\right.\\ & \quad \left.-\frac{({\omega }_{\mu }+{\omega }_{\nu })(1+{n}_{\mu }+{n}_{\nu })}{{({\omega }_{\mu }+{\omega }_{\nu })}^{2}-{z}^{2}}\right].\end{align} \tag{ 56 }$$

The only difference between ${\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[0]$ and Λ[0] (introduced in equation (18)) is that in the former the eigenvalues and eigenvectors entering the equation are those associated to the ${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$ auxiliary FCs at the centroid positions $\boldsymbol{\mathcal{R}}$, while in the latter this is not necessarily the case. The subindex $\boldsymbol{\mathcal{R}}$ in the equations above precisely indicates that the averages are calculated with a density matrix defined by $\boldsymbol{\mathcal{R}}$ and ${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$ (after a full SCHA relaxation at fixed nuclei position $\boldsymbol{\mathcal{R}}$). We will refer to the ${(\text{n}){{\Phi}}}_{\boldsymbol{\mathcal{R}}}$ tensors as the nth-order SSCHA FCs. Note that for the second-order we drop the (2) upper index.

The SSCHA code computes the free energy positional Hessian through equation (51). At the end of a SSCHA free energy functional minimization, the SSCHA matrix equation (52), with its eigenvectors and eigenvalues, is available. Thus, ${\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[0]$ is readily computable and the only quantities that need some effort to be calculated are the averages of equations (53) and (54). The code computes them through these equivalent expressions (obtained by integrating by parts):

$$\begin{equation}{({(3){{\Phi}}}_{\boldsymbol{\mathcal{R}}})}_{abc}=-\sum\limits _{pq}\enspace {({{\Psi}}_{\boldsymbol{\mathcal{R}}}^{-1})}_{ap}\enspace {({{\Psi}}_{\boldsymbol{\mathcal{R}}}^{-1})}_{bq}{\left\langle \enspace {u}^{p}{u}^{q}{\mathbb{f}}_{\mathrm{c}}\enspace \right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}\end{equation} \tag{ 57a }$$

$$\begin{equation}{({(4){{\Phi}}}_{\boldsymbol{\mathcal{R}}})}_{abcd}=-\sum\limits _{pqr}\enspace {({{\Psi}}_{\boldsymbol{\mathcal{R}}}^{-1})}_{ap}\enspace {({{\Psi}}_{\boldsymbol{\mathcal{R}}}^{-1})}_{bq}\enspace {({{\Psi}}_{\boldsymbol{\mathcal{R}}}^{-1})}_{cr}{\left\langle \enspace {u}^{p}{u}^{q}{u}^{r}{\mathbb{f}}_{d}\enspace \right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}},\end{equation} \tag{ 57b }$$

where ${\mathbf{\Psi }}_{\boldsymbol{\mathcal{R}}}$ is the Ψ matrix with $\mathbf{\Phi }={\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$ and

$$\begin{equation}\mathbb{f}\left(\boldsymbol{R}\right)={\mathbf{f}}^{(\text{BO})}\left(\boldsymbol{R}\right)-{\left\langle {\mathbf{f}}^{(\text{BO})}\left(\boldsymbol{R}\right)\right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}-{\mathbf{f}}^{{\mathcal{H}}_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}\left(\boldsymbol{R}\right).\end{equation} \tag{ 58 }$$

These averages are computed employing the stochastic approach already described in section 3.2.1 (indeed, as explained in reference [23], the choice of equation (58), among other possible alternatives, aims at reducing the statistical noise). Note that if the calculation of the free energy Hessian is performed at ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$, ${\left\langle {\mathbf{f}}^{(\text{BO})}(\boldsymbol{R})\right\rangle }_{ {\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}$ vanishes.

In order to minimize the number of energy-force calculations needed, it is advisable to compute these averages using the same ensemble used to minimize the free energy functional and obtain ${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$ (at most adding new elements to reduce the statistical noise, if needed). Of course, since in this case the ensemble is not generated from ${\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}$, an importance sampling reweighting has to be employed in order to evaluate the averages ${\left\langle \enspace \cdot \enspace \right\rangle }_{{\rho }_{\boldsymbol{\mathcal{R}},{\mathbf{\Phi }}_{\boldsymbol{\mathcal{R}}}}}$. After computing the averages, the code symmetrizes the results with respect to the space group symmetries (including the lattice translation symmetries) and the index-permutation symmetry, following the approach described in appendix D.

In order to reduce the computational cost, the SSCHA code can also compute the free energy Hessian discarding the contribution coming from the higher-order terms of the geometric-series expansion in equation (51), i.e. discarding the terms coming from ${(4){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}}$. In many cases this approximation is extremely good, but it must be checked case by case. Within this so called 'bubble' approximation, the free energy Hessian becomes

$$\begin{equation}\frac{1}{\sqrt{\boldsymbol{M}}}\cdot \frac{{\partial }^{2}F}{\partial \boldsymbol{\mathcal{R}}\partial \boldsymbol{\mathcal{R}}}\cdot \frac{1}{\sqrt{\boldsymbol{M}}}\simeq {\boldsymbol{D}}_{\boldsymbol{\mathcal{R}}}+{(3){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\boldsymbol{\mathcal{R}}}[0]\enspace \text{:}\enspace {(3){\boldsymbol{D}}}_{\boldsymbol{\mathcal{R}}}.\end{equation} \tag{ 59 }$$

Using equation (51), or its approximated expression equation (59), the SSCHA code can compute the Hessian of the free energy at any $\boldsymbol{\mathcal{R}}$. However, as said, its most significant usage is in ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$, due to its relevance to characterize displacive second-order phase transitions. In this case, equation (51) can be written in a quite explanatory form. At the end of a full SSCHA minimization, the obtained ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$ and Φ_eq define the so-called SSCHA effective harmonic Hamiltonian

$$\begin{equation}{\mathcal{H}}^{(\text{S})}=K+\frac{1}{2}(\boldsymbol{R}-{\boldsymbol{\mathcal{R}}}_{\text{eq}})\cdot {\mathbf{\Phi }}_{\text{eq}}\cdot (\boldsymbol{R}-{\boldsymbol{\mathcal{R}}}_{\text{eq}}),\end{equation} \tag{ 60 }$$

which replaces the conventional harmonic Hamiltonian to define non-interacting bosonic quasiparticles as a basis to describe the collective vibrational excitations in presence of strong anharmonic effects. In terms of the dynamical matrix ${D}_{ab}^{(\text{S})}={({{\Phi}}_{\text{eq}})}_{ab}/\sqrt{{M}_{a}{M}_{b}}$ of the SSCHA Hamiltonian ${\mathcal{H}}^{(\text{S})}$, the anharmonic generalization of the dynamical matrix D ^(F) can be written as

$$\begin{equation}{\boldsymbol{D}}^{(\text{F})}={\boldsymbol{D}}^{(\text{S})}+\mathbf{\Pi }(0),\end{equation} \tag{ 61 }$$

where

$$\begin{equation}\mathbf{\Pi }(0)={(3){\boldsymbol{D}}}_{\text{eq}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\text{eq}}[0]\enspace \text{:}\enspace {\left[\mathbb{1}-{(4){\boldsymbol{D}}}_{\text{eq}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\text{eq}}[0]\right]}^{-1}\enspace \text{:}\enspace {(3){\boldsymbol{D}}}_{\text{eq}}\end{equation} \tag{ 62 }$$

is the static SSCHA self-energy (the reason behind the use of this name will be clear in the section 4.3). In particular, in the bubble approximation we have

$$\begin{equation}{\boldsymbol{D}}^{(\text{F})}={\boldsymbol{D}}^{(\text{S})}+(\text{B}){\mathbf{\Pi }}(0),\end{equation} \tag{ 63 }$$

where

$$\begin{equation}(\text{B}){\mathbf{\Pi }}(0)={(3){\boldsymbol{D}}}_{\text{eq}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\text{eq}}[0]\enspace \text{:}\enspace {(3){\boldsymbol{D}}}_{\text{eq}}\end{equation} \tag{ 64 }$$

is the so called 'bubble' static self-energy.

In conclusion, after the SSCHA minimization, the code allows to compute the high-order SSCHA FCs, equation (57), and the free energy Hessian dynamical matrix

(depending on whether the full or only the 'bubble' static self-energy is computed) on the q -points belonging to the reciprocal space grid commensurate with the real space supercell used to generate the ensemble. Here we are explicitly using the reciprocal-space formalism, i.e. we are Fourier transforming the quantities with respect to the lattice vector indices (see appendix E1 for more details). From the softening of the eigenvalues of D ^(F)( q ) as a function of external parameters (like temperature or pressure), it is possible to observe the occurrence of second order displacive phase transitions, characterize the distortion patterns and compute the critical value of the external parameters driving it. Examples of the employment of this method are given for H₃S in figure 3 of reference [48], with the softening of an optical mode driven by pressure release, and for SnSe in figure 2 of reference [50], with the softening of the distortion mode obtained by decreasing the temperature.

The SSCHA code also allows to compute the free energy Hessian dynamical matrix D ^(F)( q ) on any reciprocal space q -point, allowing to analyze the structural instabilities incommensurate with the used supercell. After the free energy evaluation and the subsequent free energy Hessian calculation, the real-space D ^(S)( l ₁, l ₂) and ${(3){\boldsymbol{D}}}_{\text{eq}}({\boldsymbol{l}}_{1},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ are available. Here, ${D}_{ab}^{(\text{S})}({\boldsymbol{l}}_{1},{\boldsymbol{l}}_{2})$ is the real space D ^(S) matrix in which we made explicit the dependence of the lattice vectors l ₁ and l ₂ that indentify the unit cells in which atom a and b are located, respectively. Using them, the code allows to compute the static bubble $(\text{B}){\mathbf{\Pi }}(\boldsymbol{q},0)$ in any q -point, through the formula

$$\begin{align}{(\text{B}){{\Pi}}}_{\mu \nu }(\boldsymbol{q},0)& =\frac{1}{{N}_{\boldsymbol{k}}}\enspace \sum\limits _{\begin{subarray}{c}{\boldsymbol{k}}_{1}{\boldsymbol{k}}_{2}\\ {\rho }_{1}{\rho }_{2}\end{subarray}}\enspace \sum\limits _{\boldsymbol{G}}\enspace {\delta }_{\boldsymbol{G},\boldsymbol{q}+{\boldsymbol{k}}_{1}+{\boldsymbol{k}}_{2}}\mathcal{F}(0,{\omega }_{{\rho }_{1}}({\boldsymbol{k}}_{1}),{\omega }_{{\rho }_{2}}({\boldsymbol{k}}_{2}))\\ & \quad {\times}{(3){D}}_{\mu {\rho }_{1}{\rho }_{2}}(-\boldsymbol{q},-{\boldsymbol{k}}_{1},-{\boldsymbol{k}}_{2})\enspace {(3){D}}_{{\rho }_{1}{\rho }_{2}\nu }({\boldsymbol{k}}_{1},{\boldsymbol{k}}_{2},\boldsymbol{q})\enspace .\end{align} \tag{ 66 }$$

This equation is equation (64) written in reciprocal space and SSCHA normal mode components, i.e. in components of the D ^(S)( q )'s eigenvector basis. In equation (66) the k _i sums are performed on a Brillouin zone (BZ) mesh of N _k points; ${(3){D}}_{\mu {\rho }_{1}{\rho }_{2}}({\boldsymbol{k}}_{1},{\boldsymbol{k}}_{2},\boldsymbol{q})$ are the SSCHA normal components of ${(3){\boldsymbol{D}}}_{\text{eq}}({\boldsymbol{k}}_{1},{\boldsymbol{k}}_{2},\boldsymbol{q})$, the Fourier transform of ${(3){\boldsymbol{D}}}_{\text{eq}}({\boldsymbol{l}}_{1},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ in ( k ₁, k ₂, q ); ω_ρ ( k ) are the frequencies of D ^(S)( k ); the function $\mathcal{F}$ is defined by equation (56); G are reciprocal lattice vectors; and ${\delta }_{\boldsymbol{G},\boldsymbol{q}+{\boldsymbol{k}}_{1}+{\boldsymbol{k}}_{2}}$ preserves crystal momentum. In this formula the q and the k _i 's are not confined to the grid commensurate with the supercell used in the SSCHA minimization, as long as the ranges of the real space D ^(S)( l ₁, l ₂) and $(3){\boldsymbol{D}}({\boldsymbol{l}}_{1},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ are smaller than the supercell size, so as to be legitimately Fourier interpolated on any reciprocal space points (more about the Fourier interpolation in appendix E). This allows to obtain two results at once. First, the k -mesh in equation (66) can be arbitrarily increased up to convergence, so as to reach the thermodynamic limit in the evaluation of the bubble static self-energy. Second, from $(\text{B}){\mathbf{\Pi }}(\boldsymbol{q},0)$ and D ^(S)( q ), through equation (65b) the code allows to compute the free energy Hessian dynamical matrix D ^(F)( q ) (useful to detect and characterize the system instabilities) in q points not necessarily commensurate with the supercell (at variance with what is obtained with the simple Hessian calculation). This can be used, for instance, to study incommensurate second-order displacive phase transitions. In particular, this is the correct way to compute the frequencies Ω_μ ( q ) along a reciprocal-space path, where ${{\Omega}}_{\mu }^{2}(\boldsymbol{q})$ are the eigenvalues of D ^(F)( q ). This is, for example, the procedure followed to compute the (static) SCHA phonon dispersions of NbS₂ shown in figures 2 and 3 of reference [51], and to compute the interpolation-based convergence analysis shown in figure 3 of reference [54] for TiSe₂ monolayer.

The anharmonic generalization of the harmonic dynamical matrix described in the previous sections is the starting point to build a quantum anharmonic ionic dynamical theory. As shown in references [23, 48], in the context of the SSCHA it is possible to formulate an ansatz to give the expression of the one-phonon Green function G (z) for the variable $\sqrt{{M}_{a}}({R}^{a}-{\mathcal{R}}_{\text{eq}}^{a})$. This ansatz has been rigorously proved within the time dependent self-consistent harmonic approximation [38, 39]. In this dynamical theory

$$\begin{equation}{\boldsymbol{G}}^{-1}(z)={z}^{2}\mathbb{1}-\left({\boldsymbol{D}}^{(\text{S})}+\mathbf{\Pi }(z)\right),\end{equation} \tag{ 67 }$$

where D ^(S) is the dynamical matrix of the SSCHA effective harmonic Hamiltonian ${\mathcal{H}}^{(\text{S})}$, and Π(z) is the SSCHA self-energy, in general given by

$$\begin{equation}\mathbf{\Pi }(z)={(3){\boldsymbol{D}}}_{\text{eq}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\text{eq}}(z)\enspace \text{:}\enspace {\left[\mathbb{1}-{(4){\boldsymbol{D}}}_{\text{eq}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\text{eq}}(z)\right]}^{-1}\enspace \text{:}\enspace \enspace {(3){\boldsymbol{D}}}_{\text{eq}},\end{equation} \tag{ 68 }$$

and in the bubble approximation by

$$\begin{equation}(\text{B}){\mathbf{\Pi }}(z)={(3){\boldsymbol{D}}}_{\text{eq}}\enspace \text{:}\enspace {\mathbf{\Lambda }}_{\text{eq}}(z)\enspace \text{:}\enspace {(3){\boldsymbol{D}}}_{\text{eq}}.\end{equation} \tag{ 69 }$$

In the equations above we use the 'eq' subindex to specify that the eigenvalues and eigenfunctions entering the equations are obtained from Φ_eq with the centroid positions at ${\boldsymbol{\mathcal{R}}}_{\text{eq}}$.

With the Green function we obtain the spectral function $\sigma ({\Omega})=-2\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}\left[\boldsymbol{G}({\Omega}+\mathrm{i}{0}^{+})\right]$, which provides the information that can be obtained with inelastic scattering experiments. Taking explicitly into account the lattice translational symmetry (i.e. Fourier transforming the quantities with respect to the lattice vector indices) we can write

$$\begin{equation}\sigma (\boldsymbol{q},{\Omega})=-\frac{{\Omega}}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}\left[\boldsymbol{G}(\boldsymbol{q},{\Omega}+\mathrm{i}{0}^{+})\right]\end{equation} \tag{ 70 }$$

or

$$\begin{align}\sigma (\boldsymbol{q},{\Omega})& =-\frac{{\Omega}}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}\left[{\left({\Omega}+\mathrm{i}{0}^{+}\right)}^{2}\mathbb{1}+\right.\\ & \quad {\left.-\left({\boldsymbol{D}}^{(\text{S})}(\boldsymbol{q})+\mathbf{\Pi }(\boldsymbol{q},{\Omega}+\mathrm{i}{0}^{+})\right)\right]}^{-1},\end{align} \tag{ 71 }$$

where the multiplicative factor Ω/2π has been included to have, for each q , a function that integrated on the real axis gives the total number of modes 3n_a (n_a is the number of atoms in the unit cell, that may be different from the total number of atoms in the supercell N_a).

In the so-called 'static approximation', we replace the full self-energy Π(z) with the static self-energy Π(0), where z is blocked at zero. In this case the spectral function is

$$\begin{align}(\text{stat}){\sigma }(\boldsymbol{q},{\Omega})& =-\frac{{\Omega}}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}\left[{\left({\Omega}+\mathrm{i}{0}^{+}\right)}^{2}\mathbb{1}-\left({\boldsymbol{D}}^{(\text{S})}(\boldsymbol{q})\right.\right.\\ & \quad {\left.\left.+\mathbf{\Pi }(\boldsymbol{q},0)\right)\right]}^{-1}\\ & =-\frac{{\Omega}}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\left[{\left({\Omega}+\mathrm{i}{0}^{+}\right)}^{2}\mathbb{1}-{\boldsymbol{D}}^{(\text{F})}(\boldsymbol{q})\right]}^{-1},\end{align} \tag{ 72 }$$

where in the last line we have used equation (65). Therefore,

$$\begin{equation}(\text{stat}){\sigma }(\boldsymbol{q},{\Omega})=\sum\limits _{\mu }{(\text{stat}){\sigma }}_{\mu }(\boldsymbol{q},{\Omega})\enspace \end{equation} \tag{ 73 }$$

with

$$\begin{equation}{(\text{stat}){\sigma }}_{\mu }(\boldsymbol{q},{\Omega})=\frac{1}{2}\left[\delta ({\Omega}-{{\Omega}}_{\mu }(\boldsymbol{q}))+\delta ({\Omega}+{{\Omega}}_{\mu }(\boldsymbol{q}))\right],\end{equation} \tag{ 74 }$$

where ${{\Omega}}_{\mu }^{2}(\boldsymbol{q})$ are the eigenvalues of the free energy Hessian matrix D ^(F)( q ). In other words, the spectral function in the static limit is formed with delta peaks at the eigenvalues of D ^(F)( q ).

In the current version, the SSCHA code computes the full dynamical SSCHA self-energy (z ≠ 0) only in the bubble approximation with the equation (see equations (64), (66), and (69))

$$\begin{align}{(\text{B}){{\Pi}}}_{\mu \nu }(\boldsymbol{q},{\Omega}+\mathrm{i}{\delta }_{\text{se}})& =\frac{1}{{N}_{\mathrm{c}}}\enspace \sum\limits _{\begin{subarray}{c}{\boldsymbol{k}}_{1}{\boldsymbol{k}}_{2}\\ {\rho }_{1}{\rho }_{2}\end{subarray}}\enspace \sum\limits _{\boldsymbol{G}}\enspace {\delta }_{\boldsymbol{G},\boldsymbol{q}+{\boldsymbol{k}}_{1}+{\boldsymbol{k}}_{2}}\\ & \quad {\times}\enspace \mathcal{F}({\Omega}+\mathrm{i}{\delta }_{\text{se}},{\omega }_{{\rho }_{1}}({\boldsymbol{k}}_{1}),{\omega }_{{\rho }_{2}}({\boldsymbol{k}}_{2}))\\ & \quad {\times}\enspace {(3){D}}_{\mu {\rho }_{1}{\rho }_{2}}(-\boldsymbol{q},-{\boldsymbol{k}}_{1},-{\boldsymbol{k}}_{2})\enspace \enspace {(3){D}}_{{\rho }_{1}{\rho }_{2}\nu }\\ & \quad {\times}({\boldsymbol{k}}_{1},{\boldsymbol{k}}_{2},\boldsymbol{q}),\end{align} \tag{ 75 }$$

where the summation k -grid can be arbitrarily fine as long as the interpolation of the third-order SSCHA FCs can be performed (as in the static case equation (66)), and δ_se is an arbitrary small, but finite, positive smearing value used to obtain converged results in the computation. In fact, the exact result corresponds to the limiting value obtained with an infinite k -grid and a zero δ_se smearing. In actual, finite-time calculations, the converged value of the dynamic self-energy is therefore estimated in this way. For a given summation k -grid, the corresponding self-energy converged value is estimated by analyzing the result given by equation (75) for smaller and smaller δ_se values (for the used k -grid, there will be a minimum value of δ_se under which the result shows numerical instability). This analysis is performed with finer and finer summation grids until the converged value in the thermodynamic limit is obtained. In principle, a dedicated convergence study of this kind has to be performed for all the specific observables of interest.

With the dynamical SSCHA bubble self-energy, the code allows to compute the spectral function by the equation (see equation (71))

$$\begin{align}\sigma (\boldsymbol{q},{\Omega})& =-\frac{{\Omega}}{\pi }\enspace \mathrm{I}\mathrm{m}\mathrm{T}\mathrm{r}\left[{\left({\Omega}+\mathrm{i}{\delta }_{\text{id}}\right)}^{2}\mathbb{1}+\right.\\ & {\left.\quad ,-\left({\boldsymbol{D}}^{(\text{S})}(\boldsymbol{q})+(\text{B}){\mathbf{\Pi }}(\boldsymbol{q},{\Omega}+\mathrm{i}{\delta }_{\text{se}})\right)\right]}^{-1},\end{align} \tag{ 76 }$$

with δ_id another arbitrary small, but finite, positive smearing value. The role of δ_id is significant when the imaginary part of the self-energy is small. A prominent example where this happens is when the spectral function is calculated in the static approximation, i.e. when the bubble self-energy is kept fixed at the static value ${(\text{B}){{\Pi}}}_{\mu \nu }(\boldsymbol{q},0)$ (see equation (72)). Indeed, in this case the self-energy is real (Hermitian) and the computed spectral function becomes

$$\begin{align}(\text{stat}){\sigma }(\boldsymbol{q},{\Omega})& =\sum\limits _{\mu }\frac{1}{2}\enspace \mathrm{I}\mathrm{m}\left[\frac{1}{\pi }\frac{1}{{\Omega}-{{\Omega}}_{\mu }(\boldsymbol{q})+\mathrm{i}{\delta }_{\text{id}}}\right.\\ & \left.\quad +\frac{1}{\pi }\frac{1}{{\Omega}+{{\Omega}}_{\mu }(\boldsymbol{q})+\mathrm{i}{\delta }_{\text{id}}}\right],\end{align} \tag{ 77 }$$

where ${{\Omega}}_{\mu }^{2}(\boldsymbol{q})$ are the eigenvalues of D ^(F)( q ) in the bubble approximation. Therefore, for the numerical computation of the static spectral function, a finite δ_id value is necessary to recover the analytical result, equation (74), but with smeared Dirac delta functions. Actually, this is not just an extreme example, since the code really gives the opportunity to compute the spectral function in the static approximation, replacing in equation (76) the full bubble self-energy with its static value computed through equation (66). This can be used to double-check that, as expected from equations (74) and (77), the obtained spectral function is given by spikes around the frequencies of the Hessian free energy matrix D ^(F) (computed in the bubble approximation). However, the role played by δ_id is not as critical as δ_se since it is not typically system-dependent and it does not require a convergence study: in the code its default value is automatically set depending on the spacing of the energy Ω-grid used to compute the spectral function.

Given a q , the calculation of the full spectral function σ( q , Ω) through equation (76) turns out to be quite a heavy task due to the inversion of a different 3n_a × 3n_a matrix for each Ω value. The code also allows to employ a much less computational demanding approach by discarding the off-diagonal elements of the computed dynamical self-energy in the SSCHA normal modes components (i.e. the components in the D ^(S)( q )'s eigenvector basis). Within this 'no mode-mixing' approximation, which usually proves to be extremely good, the SSCHA modes keep their individuality even after the renormalization due to anharmonic effects. Indeed in this case, as in the static approximation, equation (73), the total spectral function is given by the superposition of individual mode spectral functions:

$$\begin{equation}\sigma (\boldsymbol{q},{\Omega})=\sum\limits _{\mu }\enspace {\sigma }_{\mu }(\boldsymbol{q},{\Omega}),\end{equation} \tag{ 78 }$$

where now the ( q , μ)-mode spectral function σ_μ ( q , Ω) is computed with

$$\begin{align}{\sigma }_{\mu }(\boldsymbol{q},{\Omega})=& \enspace \frac{1}{2}\left[\frac{1}{\pi }\frac{-\mathrm{I}\mathrm{m}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})}{{[{\Omega}-\mathrm{R}\mathrm{e}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})]}^{2}+{[\mathrm{I}\mathrm{m}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})]}^{2}}\right.\\ & \left.+\frac{1}{\pi }\frac{\mathrm{I}\mathrm{m}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})}{{[{\Omega}+\mathrm{R}\mathrm{e}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})]}^{2}+{[\mathrm{I}\mathrm{m}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})]}^{2}}\right]\enspace \end{align} \tag{ 79 }$$

and

$$\begin{equation}{\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})=\sqrt{{\omega }_{\mu }^{2}(\boldsymbol{q})+{{\Pi}}_{\mu \mu }(\boldsymbol{q},{\Omega}+\mathrm{i}{\delta }_{\text{se}})}.\end{equation} \tag{ 80 }$$

Therefore, computing the spectral function in the 'no mode-mixing' approximation, by measuring the deviation of σ_μ ( q , Ω) from a Dirac delta function around Ω_μ ( q ), it is possible to asses the impact that anharmonicity has on the different SSCHA modes ( q , μ), separately.

The form of the ( q , μ)-mode spectral function σ_μ ( q , Ω) in equation (79) resembles a Lorentzian, but with frequency-dependent center and width, meaning that the actual form of the spectral function σ( q , Ω) can be quite different from the superposition of true Lorentzian functions. However, in some cases the σ_μ ( q , Ω) can be expressed with good approximation as a true Lorentzian with a certain half width at half maximum (HWHM) Γ_μ ( q ) and center Ω_μ ( q ),

$$\begin{align}{\sigma }_{\mu }(\boldsymbol{q},{\Omega})& =\frac{1}{2}\left[\frac{1}{\pi }\frac{{{\Gamma}}_{\mu }(\boldsymbol{q})}{{[{\Omega}-{{\Omega}}_{\mu }(\boldsymbol{q})]}^{2}+{[{{\Gamma}}_{\mu }(\boldsymbol{q})]}^{2}}\right.\\ & \quad +\left.\frac{1}{\pi }\frac{{{\Gamma}}_{\mu }(\boldsymbol{q})}{{[{\Omega}+{{\Omega}}_{\mu }(\boldsymbol{q})]}^{2}+{[{{\Gamma}}_{\mu }(\boldsymbol{q})]}^{2}}\right],\end{align} \tag{ 81 }$$

meaning that the quasiparticle picture is still valid, even after the inclusion of anharmonicity, with the (μ, q ) quasiparticle having frequency (energy) Ω_μ ( q ) and lifetime τ_μ ( q ) = 1/(2Γ_μ ( q )). The difference between the renormalized and the 'bare' SSCHA frequency, Δ_μ ( q ) = Ω_μ ( q ) − ω_μ ( q ), is called the frequency shift of the (μ, q ) mode.

The SSCHA code offers several tools to perform such a 'Lorentzian analysis'. In general, the best Lorentzian approximation is obtained with

$$\begin{equation}{{\Omega}}_{\mu }(\boldsymbol{q})=\mathrm{R}\mathrm{e}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{{\Omega}}_{\mu }(\boldsymbol{q}))\end{equation} \tag{ 82 }$$

$$\begin{equation}{{\Gamma}}_{\mu }(\boldsymbol{q})=-\mathrm{I}\mathrm{m}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{{\Omega}}_{\mu }(\boldsymbol{q})).\end{equation} \tag{ 83 }$$

Once the dynamical self-energy and the ${\mathcal{Z}}_{\mu }(\boldsymbol{q},{\Omega})$ are computed, the SSCHA code allows to compute the single-mode spectral functions in the Lorentzian approximation, estimating the frequency Ω_μ ( q ) and HWHMs Γ_μ ( q ) in different ways. One, optional, possibility is to solve self-consistently equation (82) to estimate Ω_μ ( q ), and then Γ_μ ( q ), through equation (83). However, by default, the 'one-shot' approximation is employed with

$$\begin{equation}{(\text{os}){{\Omega}}}_{\mu }(\boldsymbol{q})=\mathrm{R}\mathrm{e}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\omega }_{\mu }(\boldsymbol{q}))\end{equation} \tag{ 84 }$$

$$\begin{equation}{(\text{os}){{\Gamma}}}_{\mu }(\boldsymbol{q})=-\mathrm{I}\mathrm{m}\enspace {\mathcal{Z}}_{\mu }(\boldsymbol{q},{\omega }_{\mu }(\boldsymbol{q})).\end{equation} \tag{ 85 }$$

If the SSCHA self-energy Π is a (small) perturbation on the SSCHA free propagator (not meaning that we are in a perturbative regime with respect to the harmonic approximation), then perturbation theory can be employed to evaluate the spectral function. If we keep the first order in the self-consistent equations equation (85), we get:

$$\begin{equation}{(\text{pert}){{\Omega}}}_{\mu }(\boldsymbol{q})=\frac{1}{2{\omega }_{\mu }(\boldsymbol{q})}\enspace \mathrm{R}\mathrm{e}\enspace {{\Pi}}_{\mu \mu }(\boldsymbol{q},{\omega }_{\mu }(\boldsymbol{q}))\end{equation} \tag{ 86 }$$

$$\begin{equation}{(\text{pert}){{\Gamma}}}_{\mu }(\boldsymbol{q})=-\frac{1}{2{\omega }_{\mu }(\boldsymbol{q})}\enspace \mathrm{I}\mathrm{m}\enspace {{\Pi}}_{\mu \mu }(\boldsymbol{q},{\omega }_{\mu }(\boldsymbol{q})+\mathrm{i}{\delta }_{\text{se}}).\end{equation} \tag{ 87 }$$

This perturbative approach is also employed by the SSCHA code to evaluate the quasiparticles' energies and lifetimes. Examples of spectral function calculations done with equations (76), (78)–(80) and (81) can be found in figure 4 of reference [48]. In figure 5 of the same reference, the anharmonic phonon frequencies and linewidths along a path, computed using the Lorenztian approximation through equations (82) and (83), are shown. The spectral function computed with equation (79) along a path is shown with a colorplot in figure 3 of reference [49] for PbTe, and in figure 4 of reference [50] for SnSe.

In conclusion, with the SSCHA code we can calculate three frequencies for a mode ( q , μ): ω_μ ( q ), Ω_μ ( q ), and Ω_μ ( q ), which are the frequency of the SSCHA auxiliary boson, the frequency coming from the SSCHA free energy Hessian (i.e. from the static approximation), and the frequency of the SSCHA quasiparticle in the Lorentzian approximation. Only the last one is a true physical quantity as it can be measured in experiments. However, the static Ω_μ ( q ) is also a physical meaningful quantity, as its zero value corresponds to a structural instability driving a second-order phase transition along the pattern characterized by the mode ( q , μ). The SSCHA provides a specific physical meaning of each of these frequencies, in contrast to other approaches used to estimate anharmonic phonons, where no distinction is usually done.

Most of the program is written in Python with an object-oriented style. The system status (density matrix) is described by a class defined in CellConstructor (phonons), that contains all the information about the system, including lattice parameters, atomic positions, and the auxiliary force-constant matrix (plus eventual extra data, as effective charges used for post-processing purposes). Methods of this class allow the user to impose symmetries on the system, constrain the auxiliary force to be positive definite (equation (32)), extract auxiliary phonon frequencies and polarization vectors, or interpolate them to other points in the BZ.

All the calculations related to the SSCHA averages are performed by the ensemble class (inside Python-sscha). This class generates and stores all the randomly displaced ionic configurations, and can submit or load the results of the energy, forces, and stress tensors calculations. It also computes the quantities related to averages on the ensemble, as the free energy, the gradients, the stress tensor, and the free energy Hessian.

Finally there are other classes, which employ the ensemble and perform the minimization of the free energy, take care of communicating with a remote cluster to run the calculation of forces and energies (see next section), and manage the post-processing to compute the spectral function (the full description of them is provided within the documentation of the code).

Most of the code is written in Python, however, the heaviest CPU-intensive calculation is written in Fortran and interfaced with Python through the f2py utility provided by numpy [55]. In particular, the calculation of the free energy gradient, the free energy Hessian, the spectral functions, the interpolation, and the symmetrization are performed by a Fortran module compiled with the code. For this reason, in order to compile and use the code, a Fortran compiler as well as LAPACK and BLAS libraries are required.

The nature of the algorithm makes it very simple to exploit massive parallelization strategies available in high performance computing (HPC) facilities. In particular, the most expensive part of the code is the calculation of BO energies, forces, and stress tensors of the generated ionic configurations in each population (the red shaded cell in the code flowchart in figure 2). Each of these calculations is independent from the others, so they can be trivially run in parallel on different computing nodes. This is a huge advantage with respect to other methods based on AIMD or PIMD, which mimic a time evolution of the system and thus require to calculate atomic forces on one configuration after the other.

The SSCHA code does not include a particular engine for computing energies, forces and stresses, but relies on external software. For this reason, it is possible to exploit the efficient parallelization already implemented by the chosen software. For example, the widely used Quantum ESPRESSO package recently implemented also a hybrid parallelization that exploits together multi-threading (OpenMP), multiprocessing (MPI), and GPU (CUDA) parallelization [56]. In this way the SSCHA code stands on the shoulders of giants, exploiting the most efficient parallelization available today.

All other steps of the code are generally computationally very cheap compared to the energy and force calculation, especially when an ab initio approach is followed. The SSCHA minimization cannot exploit so well the possibilities offered by parallelization, since each step of the main cycle depends on the previous one. Most of the computations executed in the cycle are linear algebra calculations carried out with the numpy library [55], some of them speeded up with an explicit Fortran implementation. Thanks to the numpy implementation [57], if this library is correctly compiled, the linear algebra calculations will exploit multi-threading. For this reason, the best performances of the SSCHA are obtained by executing the ab initio calculations on an HPC facility, while the SSCHA minimization on a commercial workstation in which the minimization can take few seconds.

Post-processing calculations, like the free energy Hessian and the phonon dynamical spectral functions, may be executed with additional Python scripts after the end of the SSCHA run. The calculation of the free energy Hessian has been parallelized with OpenMP (multi-threading), while the calculation of the spectral functions, which may require a dense k-point grid for the interpolation, exploits multiprocessor parallelization through MPI (both mpi4py and pypar can be used [58–60]).

Since the best performances of the code are obtained by running it in different computers, we introduced three different execution modes: manual, automatic local, and automatic remote submission.

In the manual mode, the code stops after generating the ensemble and printing on files the structures of the randomly distributed ionic configurations. At this point the user must feed these structures to a total-energy engine, e.g. a DFT code, to calculate their BO energies, atomic forces, and stress tensors. The user should prepare later specific files with the output of these calculations. Then, the SSCHA code should be manually restarted; it reads the output of energies, atomic forces, and stress tensors and runs the minimization until the exit criteria is fulfilled. After, it is up to the user to decide whether to start a new population or not. In this sense, the manual mode does not require direct interaction between the SSCHA code and any other external software, and consequently this execution mode does not require the installation of the SSCHA code on an HPC facility.

In the local automatic mode that can be scripted in Python, the code has to be supplied with an interface to an external code that is able to compute atomic forces and energies. This can be done through the atomic simulation environment (ASE) library [61], which already implements interfaces with most common ab initio codes like Quantum ESPRESSO [62, 63], VASP [64], SIESTA [65], CP2K [66], and many more. Force-field codes like LAMMPS [67] may also be used. In this execution mode, the code will proceed automatically to perform the calculations locally, and the full flowchart in figure 2 is executed without requiring any direct interaction with the user. While this execution mode is very useful, as it does not waste human time to manually restart the code at each population, it requires the most expensive part of the code, the calculation of total energies and forces, to be executed on the same machine as the SSCHA algorithm. This has a drawback when running the whole process in an HPC facility: the overall cost in terms of hours and parallel resources that needs to be allocated for the ensemble computation could be very expensive, and the SSCHA code will not exploit this amount of resources during the minimization. For this reason, the automatic local mode is indicated only when the calculation of energies and forces is fast and the requested resources are not so expensive, for example when force fields are used, such as in the SnTe example provided below.

Lastly, a remote automatic mode is also implemented. In this case the software will submit the energy and force calculations into a server through a queue job manager, and retrieve the results when finished. This last mode is the most suited for standard calculations as it exploits the HPC parallelization when computing the total energies and forces of the configurations, but runs the SSCHA minimization on a local computer, which benefits from the high speed multi-thread processors of commercial workstations. Moreover, as in the manual mode, there is no need to install the SSCHA code on a HPC cluster. Thanks to the complete automatic workflow, the only effort required by the user is to setup the communication with the clusters, which is mostly system independent.

The package is distributed as a standard Python application, and can be installed with a setup.py script. Since parts of the code are written in Fortran and C, it requires the appropriate compilers with LAPACK and BLAS libraries to be installed. Part of the Fortran subroutines are modified versions of Quantum ESPRESSO subroutines from the PHonon package, especially those regarding the symmetries. Together with the github page, we also provide the stable release in the pip repository, to facilitate installation. A different setup.py script is provided to facilitate the installation of the package on clusters to fully exploit MPI parallelization for the post-processing. The code is documented with Sphinx. We release the package and the source code under the GPLv3 license.

Hydrogen is the lightest atom in the periodic table, and, consequently, it is subject to high amplitude fluctuations even at zero Kelvin. Hydrogen atoms thus sample the V( R ) potential far from its minima. Not surprisingly, it has been shown with the SSCHA that the phonons of many hydrogen-based compounds and hydrogen itself are characterized by a huge anharmonic renormalization, impossible to capture within perturbative approaches [21, 22, 34–36, 48, 53, 69, 70]. The anharmonic renormalization of phonons in these compounds has been crucial to explain the superconducting properties of many hydrogen-based superconductors. For instance, the anomalous inverse isotope effect on palladium hydrides, which makes the deuterium compound acquire a larger superconducting critical temperature T_c than the protium compound [71, 72], is a consequence of a huge anharmonic renormalization of the phonons [21]. Also, the experimentally found high-temperature superconductivity in H₃S around 200 K [73] and in LaH₁₀ around 250 K [74, 75] at high pressures can only be explained if phonon frequencies renormalized by anharmonicity are considered in the superconductivity equations [34, 36]. In figure 8 we show the huge anharmonic renormalization of the phonon frequencies for H₃S [34, 48]. Superconductivity in hydrogen compounds can be both largely suppressed but also enhanced by anharmonicity depending on the system [76].

Zoom In Zoom Out Reset image size

The quantum effects and anharmonicity that the SSCHA captures go beyond the renormalization of phonon frequencies. For crystals with Wyckoff positions not fixed by symmetry, quantum or thermal fluctuations may strongly modify the atomic positions, resulting in a structure with atoms far from the positions that minimize the V( R ) potential, occupying, instead, those that minimize the quantum $F(\boldsymbol{\mathcal{R}})$ free energy. The change in the structure can eventually be so large that changes the crystal symmetry. For instance, the experimental crystal structure of both H₃S and LaH₁₀ compounds is stable thanks to quantum effects in the pressure range where they highest superconducting critical temperatures have been experimentally observed [34, 36]. A large modification of the structure of molecular phases of hydrogen has also been predicted within the SSCHA, which is crucial to understand the experimental Raman and infrared spectra [35, 37]. The change in the crystal structure that the SSCHA captures goes beyond the internal degrees of freedom and can largely impact also the cell parameters. In figure 8 we illustrate the apparent difference between the structure found classically from the minimum of V( R ) and the one obtained from the quantum energy landscape for LaH₁₀.

It has been recently argued [36] that the large impact of quantum effects and anharmonicity on hydrogen-based compounds is precisely due to the large electron–phonon coupling of these compounds. This means that quantum effects will lower the pressure needed to synthesize these compounds with superconducting T_c's approaching room temperature. The SSCHA method will be of great importance in the quest of new high-T_c compounds at low pressures as it can be used for crystal structure predictions in the quantum energy landscape thanks to its capacity to relax crystal structures including quantum and anharmonic effects at any target pressure.

A CDW is a structural phase transition that induces a static modulation of the electronic density. CDW transitions are often second-order phase transitions in which the frequency of the phonon mode that drives the CDW instability rapidly softens as temperature is lowered and vanishes exactly at the CDW temperature T_cdw [77–79]. As the temperature dependence of phonon frequencies is a purely anharmonic property, the SSCHA has been used to predict from first principles T_cdw in several transition metal dichalcogenides (TMDs) both in the bulk and the monolayer [41, 51, 54, 80, 81].

The standard procedure in these calculations is to apply the SSCHA for the high-symmetry phase at different temperatures and calculate the spectra associated to the free energy Hessian D ^(F). These phonons represent the static limit of the physical phonons observed experimentally, which can be accessed with the SSCHA by calculating instead the spectral function as described in section 4. At the temperature at which D ^(F) develops a null eigenvalue, the high-symmetry structure is no longer a minimum of the free energy and the CDW distortion occurs leading the structure into a phase modulated by the wave vector at which the phonon collapse occurs. In figure 9 we show as an example the temperature dependence of the phonon spectra derived from the free energy Hessian in monolayer NbSe₂ and the consequent theoretical determination of the CDW temperature [41]. In most of the cases the calculation of D ^(F) within the 'bubble' approximation yields good results for the calculation of T_cdw, and setting ${(4){\boldsymbol{D}}}_{\text{eq}}=0$ in equation (62) seems in general a good approximation. However, converging T_cdw is rather sensitive to the SSCHA supercell and rather large supercells may be needed to converge the CDW transition temperatures [41, 54].

Zoom In Zoom Out Reset image size

The capacity of the SSCHA of predicting T_cdw purely ab initio without empirical fitting parameters offers a fantastic tool to determining the physics behind CDW transitions. The force calculations needed for the SSCHA variational minimization can be performed at different theoretical levels or at different thermodynamic conditions, disentangling the driving forces of the instability. For instance, calculations within the SSCHA have enlightened the sensitivity of CDW transitions in monolayer TMDs to strain [51] and doping [54], the difference (or similarities) between the CDW transitions in bulk and the corresponding two-dimensional structures [51, 54], as well as the importance of Van der Waals forces in the melting of CDW transitions [81]. Consequently the SSCHA program is expected to have a large impact on theoretical studies of CDW transitions for many type of materials, not just TMDs.

Phase transitions related to soft phonons are also very common in ferroelectric, thermoelectric, and other functional materials. The SSCHA is again a perfect method to study these phase transitions considering that many of the high-temperature phases of these compounds are not a minimum of the BO potential V( R ), but saddle points. Thus, it becomes imperative to adopt a non-perturbative treatment of anharmonicity in order to study their thermodynamic and transport properties such as the thermal conductivity. Whether these phase transitions are purely second-order or first order it is not always evident experimentally, unless a clear softening to zero of a phonon mode at the transition temperature is observed. The SSCHA can distinguish between continuous and discontinuous transitions as discussed in the practical example provided in section 6. For instance, in order to convincingly show that the transition between the high-temperature Cmcm phase of SnSe and the low-temperature Pnma is second order, at the temperature at which the free energy Hessian developed a negative eigenvalue a SSCHA relaxation was performed starting from the low-temperature phase. It was shown that the Pnma phase relaxed at this temperature into the Cmcm, showing that the Pnma phase is no longer a minimum of the free energy [50]. The SSCHA has also been used to study phase transitions in the similar SnS [52] and the ferroelectric SnTe [49].

Many of these semiconducting calchogenides are among the most efficient thermoelectric materials due to their very low thermal conductivity. The low value of the thermal conductivity of these materials is linked to the very large linewidth of its phonon modes. The anharmonic interaction is the main responsible for the large linewidths of the phonons and, consequently, their low lifetimes. Thanks to the strong anharmonic coupling, many of these compounds develop very anomalous spectral functions with satellite peaks and a clear departure from the Lorentzian-like behavior. Such anomalies can be very misleading for the interpretation of experiments, since the emergence of extra peaks can be misinterpreted with phase transitions. The SSCHA is a perfect method for capturing these subtleties as it provides the spectral function σ( q , Ω) without the Lorentzian approximation. It has been used to understand the complex σ( q , Ω) in PbTe, SnTe, SnSe, and SnS [49, 50, 52]. In figure 10 we show the spectral function calculated within the SSCHA for Cmcm SnSe at 800 K, where the σ_μ ( q , Ω) contribution of some particular modes is clearly anomalous and deviates from the standard Lorentzian picture.

Zoom In Zoom Out Reset image size

With the phonon frequencies and the phonon linewidths obtained with the SSCHA, transport properties such as the thermal conductivity can be calculated with an external code, for instance, within Boltzmann transport equations [82, 83]. It has been shown that employing the SSCHA phonon scattering tensor ${(3){{\Phi}}}_{\boldsymbol{\mathcal{R}}}$ in the thermal transport calculations leads to a very good agreement with experimental results, in contrast with what was obtained by employing the standard third-order derivatives of the BO total energy. The difference is that the former includes higher-order anharmonic terms coming from the average over the thermal ensemble (equation (53)). Both in good thermoelectric SnSe and SnS compounds, higher order terms captured by ${(3){{\Phi}}}_{\boldsymbol{\mathcal{R}}}$ reduce considerably the thermal conductivity, bringing it closer to the experimentally observed values [50, 52]. Therefore, the SSCHA also provides a fantastic platform to calculate the basic ingredients for transport properties when high-order terms of the BO potential are important both for the renormalization of the phonon frequencies and the anharmonic scattering terms. Due to the large effort devoted currently to the quest of more efficient thermoelectric materials, the SSCHA may become a reference method to understand the thermoelectric properties of materials and predict the efficiency of new promising compounds, which overcomes the limits of standard harmonic and perturbative approaches.

As mentioned above, beyond those examples listed above, the SSCHA code provides an efficient platform to calculate any property affected by ionic fluctuations, specially when it is affected by strong anharmonicity. For instance, it has been used to determine the muon implantation sites in metallic systems and to understand the effect of the large muon quantum fluctuations on the contact hyperfine field [84]. The SSCHA has also been employed to understand the thermal expansion and the behavior of low-energy acoustic modes of graphene [85], finally explaining the origin of the sound propagation and the non-diverging bending rigidity of graphene as well as any other strictly two-dimensional membrane. Many other exciting applications of the SSCHA code to interesting physical problems are expected in the coming years.

Besides space group symmetries, also the acoustic sum rule (ASR) must be imposed. The ASR is a condition that arises from the momentum conservation (the center of mass of the system is fixed). The energy must not change after a rigid translation of the whole system. This can be translated in a trivial condition for the FC matrix in the supercell:

$$\begin{equation}\sum\limits _{t}\enspace {{\Phi}}_{st}^{\alpha \beta }=\sum\limits _{s}\enspace {{\Phi}}_{st}^{\alpha \beta }=0.\end{equation} \tag{ D4 }$$

In general, the SSCHA gradient computed from a finite ensemble violates this condition due to the stochastic noise. We enforce the sum rule on the gradient at each step. As for the symmetries, also matrices that satisfy the ASR define a linear subspace. Thus we define the orthogonal projector operator that enforces the ASR as

$$\begin{equation}{\mathbf{\Phi }}^{(\mathrm{a}\mathrm{s}\mathrm{r})}=\boldsymbol{P}\mathbf{\Phi }{\boldsymbol{P}}^{{\dagger}}.\end{equation} \tag{ D5 }$$

The projection matrix in real space is

$$\begin{equation}{P}_{st}^{\alpha \beta }={\delta }_{st}{\delta }_{\alpha \beta }-\frac{{\delta }_{\alpha \beta }}{{n}_{\mathrm{a}}}\sum\limits _{u=1}^{{n}_{\mathrm{a}}}\enspace {\delta }_{tu}.\end{equation} \tag{ D6 }$$

This operation only affects the dynamical matrix at Γ. The same projector is employed to impose the ASR on the forces:

$$\begin{equation}{{f}^{(\mathrm{a}\mathrm{s}\mathrm{r})}}_{s}^{\alpha }={(\boldsymbol{P}f)}_{s}^{\alpha }={f}_{s}^{\alpha }-\frac{1}{{n}_{\mathrm{a}}}\sum\limits _{k=1}^{{n}_{\mathrm{a}}}\enspace {f}_{k}^{\alpha }\end{equation} \tag{ D7 }$$

Notably, it can be proved that this procedure does not spoil the symmetrization described above.

This ASR imposition procedure analytically cancels out the frequencies of acoustic modes at Γ and any rigid translation of the atomic positions, thus it is the most indicated for the SSCHA minimization. A different approach, implemented for the Fourier interpolation, is described in appendix E2. The latter affects not just phonons at Γ, and, thus, it is more suited for interpolating dynamical matrices close to the BZ center.

The SSCHA code is designed to be used with crystals, thus it takes advantage of lattice periodicity and Fourier transforms the relevant quantities with respect to the lattice vectors. That allows to make independent analysis for each q point in reciprocal space. When we need to stress this aspect, we will modify the notation adopted, partitioning the supercell atomic index into a unit-cell index plus a lattice index (s, l ), with s now ranging from 1 to n_a (the number of atoms in the unit cell), and l being a 3 dimensional integer vector assuming N_c total values (the number of unit cells forming the supercell). Thus, in this notation, in general we will have nth order tensors in a 3n_a dimensional space (indicated with bold symbols, in free-component notation), which in real space depend on n lattice-vector parameters, $(n){\boldsymbol{D}}({\boldsymbol{l}}_{1},\dots {\boldsymbol{l}}_{n})$ (to be precise, due to the translation symmetry, this real-space tensor actually depends only on n − 1 independent values l _i ). The reciprocal-space expression of such a tensor, $(n){\boldsymbol{D}}({\boldsymbol{q}}_{1},\dots ,{\boldsymbol{q}}_{n})$, is obtained through the Fourier transform

$$\begin{equation}(n){\boldsymbol{D}}({\boldsymbol{q}}_{1},\dots ,{\boldsymbol{q}}_{n})=\frac{1}{{N}_{\text{c}}}\sum\limits _{{\boldsymbol{l}}_{1}\dots {\boldsymbol{l}}_{n}}{\text{e}}^{\text{i}\sum\limits _{h}{\boldsymbol{q}}_{h}\cdot {\boldsymbol{l}}_{h}}\enspace (n){\boldsymbol{D}}({\boldsymbol{l}}_{1},\dots ,{\boldsymbol{l}}_{n}).\end{equation} \tag{ E1 }$$

Notice that, due to the lattice translation symmetry, D ^(S)( q ₁, ..., q _n ) is zero unless ∑_h q _h is a reciprocal lattice vector, thus we have again only n − 1 independent parameters q _i . In particular, after the calculation performed on a real-space supercell, for each q point of the commensurate grid of the reciprocal-space unit cell, the SSCHA code computes the Fourier transformed matrices D ^(S)(− q , q ), which we will shortly indicate as D ^(S)( q ), and the relative eigenvalues ω_μ ( q ) and eigenvectors e _μ ( q ). Similarly, the Hessian calculation provides the matrix D ^(F)(− q , q ), which we indicate as D ^(F)( q ), where (see equation (61))

$$\begin{equation}{\boldsymbol{D}}^{(\text{F})}(\boldsymbol{q})={\boldsymbol{D}}^{(\text{S})}(\boldsymbol{q})+\mathbf{\Pi }(\boldsymbol{q},0),\end{equation} \tag{ E2 }$$

and its eigenvalues Ω_μ ( q ) and eigenvectors f _μ ( q ).

The SSCHA code computes the FCs in real space supercells with periodic boundary conditions (PBCs). As shown in the previous section, a crucial feature of the SSCHA code is the use of the Fourier interpolation technique in order to extrapolate the results to the thermodynamic limit (infinite supercell results) without recurring to expensive large supercell calculations. In order to Fourier interpolate the computed FCs on arbitrary points of the reciprocal space, as a first thing it is necessary to reconstruct the real-space infinite-crystal FCs from them. Roughly speaking, this is done by removing the PBCs, i.e. superlattice equivalent atoms are not considered identical anymore, and assuming that only the FCs between atoms in the same supercell are different from zero. Of course, this gives correct results as long as the supercell used in the calculations is large enough to consider negligible the FCs between atoms at distances comparable with the distances between the periodic boundary replica. However, an intrinsic arbitrariness is present in this recipe, due to the fact that the supercell is not univocally defined and the choice of different supercells leads to different interpolation results (i.e. as long as the reciprocal-space point in which we are interpolating does not belong to the original commensurate grid, different—yet superlattice equivalent—lattice points give different contributions to the Fourier transform). This problem is solved by wisely selecting the supercell according to a prescription based on a physical principle: among equivalent superlattice points, the ones closest to each other must be selected. This procedure defines the so called 'centering' of the FCs and, as explained, it is a necessary step to be done before Fourier interpolating the real space FCs. The SSCHA code centers 2nd and 3rd order FCs (with a procedure that can be generalized to any nth-order FCs. In particular, the next release of the code will apply the same procedure to center and interpolate the 4th order FCs). Here we explictly describe the 3rd FCs centering algorithm [90].

The PBCs are defined on a superlattice ${\mathcal{R}}_{\text{lat}}^{(\text{S})}$ of the original lattice ${\mathcal{R}}_{\text{lat}}$. The lattice vectors set ${\mathcal{R}}_{\text{lat}}$ can be equivalently described as the superlattice ${\mathcal{R}}_{\text{lat}}^{(\text{S})}$ plus the basis given by the lattice vectors in a superlattice unit cell SC. In other words, a lattice vector $\boldsymbol{l}{\in \mathcal{R}}_{\text{lat}}$ identifies a set of superlattice-equivalent lattice vectors $\left\{\boldsymbol{l}+\boldsymbol{T}\enspace \;\text{with}\enspace \boldsymbol{T}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}\right\}$, and we have ${\mathcal{R}}_{\text{lat}}=\left\{\boldsymbol{l}+\boldsymbol{T}\enspace \text{with}\enspace \boldsymbol{l}\in SC\enspace ,\boldsymbol{T}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}\right\}$. Given three atoms s₁, s₂, s₃ in the unit cells 0, l ₂, l ₃, respectively (due to the lattice translation symmetry we can confine the first atom to the origin unit cell), they indentify a triangle with vertices in ${\boldsymbol{\tau }}_{{s}_{1}},{\boldsymbol{\tau }}_{{s}_{2}}+{\boldsymbol{l}}_{2},{\boldsymbol{\tau }}_{{s}_{3}}+{\boldsymbol{l}}_{3}$ (${\boldsymbol{\tau }}_{{s}_{i}}$ is the position vector of atom s_i in the original unit cell). For these three points we define the weight ${\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\boldsymbol{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ in this way: it is zero if there is at least another 'equivalent-vertices' triangle having as vertices points ${\boldsymbol{\tau }}_{{s}_{1}}$, ${\boldsymbol{\tau }}_{{s}_{2}}+{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2}$, ${\boldsymbol{\tau }}_{{s}_{3}}+{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3}$ with ${\boldsymbol{T}}_{2},{\boldsymbol{T}}_{3}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}$ (i.e. points that are superlattice-equivalent to ${\boldsymbol{\tau }}_{{s}_{1}},{\boldsymbol{\tau }}_{{s}_{2}}+{\boldsymbol{l}}_{2},{\boldsymbol{\tau }}_{{s}_{3}}+{\boldsymbol{l}}_{3}$) with smaller perimeter, otherwise it is the inverse of the number of equivalent-vertices triangles having the same (minimal) perimeter. In formulas, indicated with ${\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\boldsymbol{0},{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})$ the perimeter of the triangle with vertices ${\boldsymbol{\tau }}_{{s}_{1}}$, ${\boldsymbol{\tau }}_{{s}_{2}}+{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2}$, ${\boldsymbol{\tau }}_{{s}_{3}}+{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3}$, this amounts to

$$\begin{equation}{\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})=\left\{\begin{array}{ccc}& 0\quad \hspace{199.16928pt} & \text{if}\quad \begin{array}{c}\qquad \exists \enspace {\boldsymbol{T}}_{2},{\boldsymbol{T}}_{3}\in {\mathcal{R}}_{\text{lat}}^{(\text{S})}:\\ {\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3}){< }{\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\end{array}\\ & {\left[\begin{array}{c}\#({\boldsymbol{T}}_{2},{\boldsymbol{T}}_{3})\in {\mathcal{R}}_{\text{lat}}^{(\text{S})}:\\ {\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})={\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\end{array}\right]}^{-1}\quad & \text{if}\quad \begin{array}{c}\qquad \nexists \enspace {\boldsymbol{T}}_{2},{\boldsymbol{T}}_{3}\in {\mathcal{R}}_{\text{lat}}^{(\text{S})}:\\ {\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3}){< }{\mathcal{P}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\end{array}\end{array}\right.\end{equation} \tag{ E3 }$$

The weights ${\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ are pure geometrical factors, different from zero for 'compact' three-atom clusters, and they satisfy the normalization

$$\begin{equation}\begin{aligned}\sum\limits _{{\boldsymbol{T}}_{2},{\boldsymbol{T}}_{3}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}}\enspace {\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})=1\\ \begin{aligned}& \forall \enspace {\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}{\in \mathcal{R}}_{\text{lat}}\\ & \forall \enspace {s}_{1},{s}_{2},{s}_{3}\in \left\{1,\dots ,{n}_{\text{a}}\right\}.\end{aligned}\end{aligned}\end{equation} \tag{ E4 }$$

The weights are used to define the centering. Given a 3rd-order FCs, ${{\Phi}}_{{s}_{1},{s}_{2},{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$, its 'centered' version $(\text{cent}){{\Phi}}{}_{{s}_{1}{s}_{2}{s}_{3}}{}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ is given by

$$\begin{equation}(\text{cent}){{\Phi}}{}_{{s}_{1}{s}_{2}{s}_{3}}{}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})={{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\enspace {\times}\enspace {\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}),\end{equation} \tag{ E5 }$$

where we have separately indicated Cartesian (α_h ) and atomic (s_h ) indices. The idea behind this definition is pretty simple: once the PBCs are discarded, of the infinite set of superlattice-equivalent atoms only the 'closest one' are characterized by a FC different from zero. If there are several equivalent triplets at the minimal reciprocal distance, all of them are considered (to preserve the symmetry) and the FCs are consequently scaled (to avoid a wrong multiple counting effect). The centering definition has some degree of arbitrariness, though, due to the arbitrariness of the criterion employed to evaluate the 'size' of a three atoms cluster. We took the perimeter of the triangle, a criterion that is a direct generalization of the distance between atoms, which is the one used in the 2nd order FCs centering. More in general, for an n-atoms cluster this size measure is readily generalized as the sum of the distances between all the n(n − 1)/2 couples of atoms. However, even if in principle other choices could be done, this arbitrariness is immaterial as long as the supercell calculation is large enough (in the thermodynamic limit all the possible choices, if reasonable, are expected to be equivalent).

A delicate issue is associated with the centering: the spoiling of the ASR. For a nth-order FC, the ASR is

$$\begin{equation}\sum\limits _{{\boldsymbol{l}}_{i}}\enspace \sum\limits _{{s}_{i}}\enspace {{\Phi}}_{{s}_{1}\dots {s}_{i}\dots {s}_{n}}^{{\alpha }_{1}\dots {\alpha }_{i}\dots {\alpha }_{n}}({\boldsymbol{l}}_{1},\dots ,{\boldsymbol{l}}_{i},\dots ,{\boldsymbol{l}}_{n})=0\begin{array}{cccc}& \forall \enspace {\alpha }_{h}\in \left\{x,y,z\right\}& & \\ & \forall \enspace {s}_{h}\in \left\{1\dots {n}_{\text{a}}\right\}& & \text{with}\enspace h\ne i\\ & \forall \enspace {\boldsymbol{l}}_{h}\in {\mathcal{R}}_{\text{lat}}& & \text{with}\enspace h\ne i\end{array}\quad .\end{equation} \tag{ E6 }$$

The ASRis crucial, among other things, to have the correct acoustic phonon dispersion at and close to Γ. The FCs computed with SSCHA (in supercells with PBCs) fulfill the ASR but, in general, the centering spoils it (except for the n = 2 case). In fact, in general the centered FCs fulfill a 'weaker' version of the ASR in equation (E6), since only the simultaneous sum on n − 1 indices is zero:

$$\begin{equation}\sum\limits _{{\boldsymbol{l}}_{{i}_{1}},\dots ,{\boldsymbol{l}}_{{i}_{n-1}}}\enspace \sum\limits _{{s}_{{i}_{1}},\dots ,{s}_{{i}_{n-1}}}(\text{cent}){{\Phi}}{}_{{s}_{1}\dots {s}_{n}}{}^{{\alpha }_{1}\dots {\alpha }_{n}}({\boldsymbol{l}}_{1},\dots ,{\boldsymbol{l}}_{n})=0.\end{equation} \tag{ E7 }$$

In particular, this explains why the centering of 2nd-order FCs does not spoil the ASR (for n = 2 the weak ASR is nothing but the proper ASR, as the sum over n − 1 indices coincides with the sum over one index).

In order to see why this happens let us consider, as an example, the 3rd-order FCs case and the sum over the third index. It is:

$$\begin{equation}\sum\limits _{{s}_{3}}\enspace \sum\limits _{{\boldsymbol{l}}_{3}{\in \mathcal{R}}_{\text{lat}}}\enspace (\text{cent}){{\Phi}}{}_{{s}_{1}{s}_{2}{s}_{3}}{}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\end{equation} \tag{ E8 }$$

$$\begin{equation}\qquad =\sum\limits _{{s}_{3}}\enspace \sum\limits _{{\boldsymbol{l}}_{3}\in SC}\enspace \sum\limits _{{\boldsymbol{T}}_{3}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}}(\text{cent}){{\Phi}}{}_{{s}_{1}{s}_{2}{s}_{3}}{}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})\end{equation} \tag{ E9 }$$

$$\begin{equation}\begin{array}{lc}& \qquad =\sum\limits _{{s}_{3}}\enspace \sum\limits _{{\boldsymbol{l}}_{3}\in SC}\enspace \sum\limits _{{\boldsymbol{T}}_{3}\in {\mathcal{R}}_{\text{lat}}^{(\text{S})}}{{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})\enspace \enspace \\ & \quad \qquad {\times}\enspace {\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})\end{array}\end{equation} \tag{ E10 }$$

$$\begin{equation}\begin{array}{lc}& \qquad =\sum\limits _{{s}_{3}}\enspace \sum\limits _{{\boldsymbol{l}}_{3}\in SC}\enspace \sum\limits _{{\boldsymbol{T}}_{3}\in {\mathcal{R}}_{\text{lat}}^{(\text{S})}}{{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\enspace \enspace \\ & \qquad {\times}\enspace {\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})\end{array}\phantom{\rule{0ex}{0ex}}\end{equation} \tag{ E11 }$$

$$\begin{align}& \qquad =\sum\limits _{{s}_{3}}\enspace \sum\limits _{{\boldsymbol{l}}_{3}\in SC}\enspace {{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\\ & \quad \qquad {\times}\underset{\text{nonconstant}\;\text{w.r.t.}\;{s}_{3}\;\text{and}\;{\boldsymbol{l}}_{3}}{\underbrace{\sum\limits _{{\boldsymbol{T}}_{3}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}}\enspace {\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})}}.\end{align} \tag{ E12 }$$

Since the last factor, highlighted with a brace under, in general is not constant with respect to s₃ and l ₃, it cannot be factored out from the sums, so that the ASR for the original ${{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$

$$\begin{equation}\sum\limits _{{s}_{3}}\enspace \sum\limits _{{\boldsymbol{l}}_{3}\in SC}\enspace {{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})=0\end{equation} \tag{ E13 }$$

cannot be used to obtain the ASR for the centered $(\text{cent}){{\Phi}}{}_{{s}_{1}{s}_{2}{s}_{3}}{}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$. However, using equation (E4), with similar passages we can show that the sum over the last two indices is zero:

$$\begin{equation}\sum\limits _{{s}_{2},{s}_{3}}\enspace \enspace \sum\limits _{{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}{\in \mathcal{R}}_{\text{lat}}}\enspace (\text{cent}){{\Phi}}{}_{{s}_{1}{s}_{2}{s}_{3}}{}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\end{equation} \tag{ E14 }$$

$$\begin{align}& \qquad =\sum\limits _{{s}_{2},{s}_{3}}\enspace \enspace \sum\limits _{{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\in SC}\enspace {{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\\ & \quad \qquad {\times}\underset{=1}{\underbrace{\sum\limits _{{\boldsymbol{T}}_{2},{\boldsymbol{T}}_{3}{\in \mathcal{R}}_{\text{lat}}^{(\text{S})}}\enspace {\mathcal{W}}_{{s}_{1}{s}_{2}{s}_{3}}(\mathbf{0},{\boldsymbol{l}}_{2}+{\boldsymbol{T}}_{2},{\boldsymbol{l}}_{3}+{\boldsymbol{T}}_{3})}}\end{align} \tag{ E15 }$$

$$\begin{equation}\qquad =\sum\limits _{{s}_{2},{s}_{3}}\enspace \enspace \sum\limits _{{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\in SC}\enspace {{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\end{equation} \tag{ E16 }$$

$$\begin{equation}\qquad =0,\end{equation} \tag{ E17 }$$

thus the weak ASR is fulfilled.

The spoiling of the ASR after the centering dictates to impose it. In principle, there is not a unique way of doing it, as imposing the ASR on FCs simply consists in finding new FCs that fulfill the ASR and differ the least from the original FCs (according to some reasonable but arbitrary metric). In this release of the SSCHA code we impose the ASR by employing an iterative procedure, consisting of two steps [90]. First, the ASR is imposed on one index (the last one, for example). This spoils the permutation symmetry, which is consequently imposed. In general, the resulting permutation-symmetric FCs do not fulfill the ASR yet, thus this procedure is repeatedly applied until the permutation-symmetric FCs fulfill the ASR within a certain tolerance. The imposition of the permutation symmetry is a straightforward task. The ASR is imposed on the third index of a centered 3rd-order FCs by updating its values on the compact three-atom clusters that defined the centering (in order to preserve the short-sightedness of the centered FCs even after the ASR imposition). Given a centered ${{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$, the ${\tilde {{\Phi}}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ that fulfills the ASR on the third index is computed with

$$\begin{align}{\tilde {{\Phi}}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})& ={{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})-{\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\vert p)\\ & \quad {\times}\sum\limits _{{\bar{s}}_{3},{\bar{\boldsymbol{l}}}_{3}}\enspace {{\Phi}}_{{s}_{1}{s}_{2}{\bar{s}}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\bar{\boldsymbol{l}}}_{3}),\end{align} \tag{ E18 }$$

where ${\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\vert p)$ is the scaling factor

$$\begin{equation}{\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\vert p)=\left\{\right.\begin{aligned}& \frac{{\vert {{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})\vert }^{p}}{\sum\limits _{{\bar{s}}_{3},{\bar{\boldsymbol{l}}}_{3}}{\vert {{\Phi}}_{{s}_{1}{s}_{2}{\bar{s}}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\bar{\boldsymbol{l}}}_{3})\vert }^{p}}& & \quad \text{if}& & p=0\enspace \text{or}\enspace p{ >}0\;\text{and}\;\sum\limits _{{\bar{s}}_{3},{\bar{\boldsymbol{l}}}_{3}}\vert {{\Phi}}_{{s}_{1}{s}_{2}{\bar{s}}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\bar{\boldsymbol{l}}}_{3})\vert \ne 0\\ & 0& & \quad \text{if}& & p{ >}0\;\text{and}\;\sum\limits _{{\bar{s}}_{3},{\bar{\boldsymbol{l}}}_{3}}\vert {{\Phi}}_{{s}_{1}{s}_{2}{\bar{s}}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\bar{\boldsymbol{l}}}_{3})\vert =0\end{aligned},\end{equation} \tag{ E19 }$$

with p a non-negative real number which can be arbitrarily fixed to optimize the calculation performances (in the equation above the convention 0⁰ = 1 has been adopted). The ${\tilde {{\Phi}}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ defined through equation (E18) fulfills the ASR on the third index, since for any p ⩾ 0 the scaling factor fulfills the normalization condition

$$\begin{equation}\sum\limits _{{s}_{3},{\boldsymbol{l}}_{3}}\enspace {\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\vert p)=1.\end{equation} \tag{ E20 }$$

The value of p has effects on the way the different terms of ${{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3})$ are scaled. For p = 0 the scaling factor is a pure geometric quantity related to the three atoms clusters. Indeed, given s₁, s₂, l ₂, the scaling factor ${\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}\vert p=0)$ is fully determined (it is the same for all the α_h , s₃, l ₃) and, in particular, it does not depend on the FCs value. On the contrary, for p ≠ 0, given α_h , s₁, s₂, l ₂ we have

$$\begin{equation}\frac{{\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}^{\prime }}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}^{\prime }\vert p)}{{\mathcal{K}}_{{s}_{1}{s}_{2}{s}_{3}^{{\prime\prime}}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}^{{\prime\prime}}\vert p)}={\left\vert \frac{{{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}^{{\prime\prime}}}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}^{{\prime\prime}})}{{{\Phi}}_{{s}_{1}{s}_{2}{s}_{3}^{\prime }}^{{\alpha }_{1}{\alpha }_{2}{\alpha }_{3}}(\mathbf{0},{\boldsymbol{l}}_{2},{\boldsymbol{l}}_{3}^{\prime })}\right\vert }^{p}\end{equation} \tag{ E21 }$$

so that if p > 1 the scaling factor is higher (lower) for FCs have lower (higher) absolute value, otherwise the opposite.

In ionic crystals the nuclei displacement induces dipoles (proportional to the Born effective charge tensors), and this adds a dipole–dipole interaction term to the interatomic forces. This contribution, because of its long-range character (it goes as the inverse of the third power of the nuclei distances), is not suited to be Fourier interpolated and it is at the origin of the nonanalytic behavior of the dynamical matrix at Γ, with (in general anisotropic) LO–TO splitting of the phonon frequencies at BZ center. The long-range dipole–dipole contribution to the FCs can be calculated analytically since it is fully determined by the Born effective charges ${({Z}_{s}^{{\ast}})}^{\alpha \beta }$ (effective charge tensor of atom s) and the electronic dielectric permittivity tensor ${({{\epsilon}}^{\infty })}^{\alpha \beta }$, which can both be calculated from first principles. For a given q ∈ BZ, this dipole–dipole contribution is given by [91, 92]

$$\begin{equation}{\mathbf{\Phi }}_{st}^{(\text{dd})}(\boldsymbol{q})=\hat{\mathbf{\Phi }}{}_{st}{}^{(\text{dd})}(\boldsymbol{q})-{\delta }_{st}\sum\limits _{\bar{t}}\enspace \hat{\mathbf{\Phi }}{}_{s\bar{t}}{}^{(\text{dd})}(\boldsymbol{q}=\mathbf{0})\end{equation} \tag{ E22 }$$

with

$$\begin{align}\hat{\mathbf{\Phi }}{}_{st}{}^{(\text{dd})}(\boldsymbol{q})& =\frac{4\pi }{{{\Omega}}_{\text{Vol}}}{\sum\limits _{\boldsymbol{G}}}^{\prime }\frac{\left[(\boldsymbol{G}+\boldsymbol{q})\cdot {\boldsymbol{Z}}_{s}^{{\ast}}\right]\otimes \left[(\boldsymbol{G}+\boldsymbol{q})\cdot {\boldsymbol{Z}}_{t}^{{\ast}}\right]}{(\boldsymbol{G}+\boldsymbol{q})\cdot {\boldsymbol{{\epsilon}}}^{\infty }\cdot (\boldsymbol{G}+\boldsymbol{q})}\\ & \quad {\times}{\mathrm{e}}^{-\frac{(\boldsymbol{G}+\boldsymbol{q})\cdot {\boldsymbol{{\epsilon}}}^{\infty }\cdot (\boldsymbol{G}+\boldsymbol{q})}{4{\eta }^{2}}}\enspace {\text{e}}^{\text{i}(\boldsymbol{G}+\boldsymbol{q})\cdot ({\boldsymbol{\tau }}_{s}-{\boldsymbol{\tau }}_{t})},\end{align} \tag{ E23 }$$

where we have explicitly indicated only the atomic indices (i.e. we are using component-free notation for the Cartesian indices), η is a parameter whose value has to be large enough to allow to include only the reciprocal space terms in the Ewald sum, and ∑ _G ' is the sum over reciprocal lattice vectors such that G + q ≠ 0 (the sum includes as many G 's as it is necessary to reach the convergence for the considered η) [62].

Once ${\boldsymbol{Z}}_{s}^{{\ast}}$ and ^∞ are available, the problem caused to the Fourier interpolation by the long-range dipole–dipole interaction is thus bypassed in the SSCHA code in two steps. First, from the Φ( q ) calculated on a (coarse) grid of q point of the BZ, the corresponding dipole–dipole terms Φ^(dd)( q ) are subtracted and the resulting short range FCs is Fourier transformed to the real space. Subsequently, this real space short-range FCs, Φ^(sr)( l ), can be Fourier transformed back to any k ∈ BZ and the corresponding long-range dipole–dipole analytical contribution Φ^(dd)( k ) is added [62]:

$$\begin{equation}\begin{gathered}\mathbf{\Phi }(\boldsymbol{q})\hspace{2pt}\text{on}\enspace \mathrm{B}\mathrm{Z}\enspace \boldsymbol{q}-\text{grid}\end{gathered}\begin{subarray}{c}\text{Subtract}\;\text{dipole}-\text{dipole}\;\text{interaction}\;\text{terms}\\ {\mathbf{\Phi }}^{(dd)}(\boldsymbol{q})\\ +\\ \text{Fourier}\;\text{transform}\\ \text{to}\;\text{real}\;\text{space}\end{subarray}{\to }\hspace{2pt}{\mathbf{\Phi }}^{(\text{sr})}(\boldsymbol{l})\quad \begin{subarray}{c}\text{Fourier}\;\text{transform}\\ \text{back}\;\text{to}\enspace \enspace \boldsymbol{k}\enspace \enspace \in \enspace \enspace \mathrm{B}\mathrm{Z}\\ +\\ \text{Add}\;\text{dipole}-\text{dipole}\;\text{interaction}\;\text{term}\\ {\mathbf{\Phi }}^{\left(dd\right)}\left(\boldsymbol{k}\right)\end{subarray}{\to }\mathbf{\Phi }\left(\boldsymbol{k}\right)\end{equation} \tag{ E24 }$$

The dipole–dipole correction to the FCs given by equations (E22) and (E23) is nonanalytic at zone center and its q → 0 limit depends on the direction $\hat{\boldsymbol{q}}=\boldsymbol{q}/{\Vert}\boldsymbol{q}{\Vert}$ along which the limit is performed:

$$\begin{equation}\underset{\delta \to {0}^{+}}{\mathrm{lim}}\enspace {\mathbf{\Phi }}_{st}^{(\text{dd})}(\delta \hat{\boldsymbol{q}})={\mathbf{\Phi }}_{st}^{(\text{dd})}(\mathbf{0})+{\mathbf{\Phi }}_{st}^{(\text{dd}-\text{na})}(\hat{\boldsymbol{q}}),\end{equation} \tag{ E25 }$$

where

$$\begin{equation}{\mathbf{\Phi }}_{st}^{(\text{dd}-\text{na})}(\hat{\boldsymbol{q}})=\frac{4\pi }{{{\Omega}}_{\text{Vol}}}\frac{\left[\hat{\boldsymbol{q}}\cdot {\boldsymbol{Z}}_{s}^{{\ast}}\right]\otimes \left[\hat{\boldsymbol{q}}\cdot {\boldsymbol{Z}}_{t}^{{\ast}}\right]}{\hat{\boldsymbol{q}}\cdot {\boldsymbol{{\epsilon}}}^{\infty }\cdot \hat{\boldsymbol{q}}}\end{equation} \tag{ E26 }$$

is the nonanalytic zone-center correction term. When a phonon dispersion through Γ is calculated, the SSCHA code includes the nonanalytic correction term in the zone center, with the direction given by the followed path [62]. When the SSCHA code calculate the spectral properties (static or dynamic), it adds the nonanalytic correction term in the zone center dynamical matrix (necessary for the integral over the BZ) from a random direction.