The Stochastic Self-Consistent Harmonic Approximation: Calculating Vibrational Properties of Materials with Full Quantum and Anharmonic Effects

The efficient and accurate calculation of how ionic quantum and thermal fluctuations impact the free energy of a crystal, its atomic structure, and phonon spectrum is one of the main challenges of solid state physics, especially when strong anharmonicy invalidates any perturbative approach. To tackle this problem, we present the implementation on a modular Python code of the stochastic self-consistent harmonic approximation method. This technique rigorously describes the full thermodyamics of crystals accounting for nuclear quantum and thermal anharmonic fluctuations. The approach requires the evaluation of the Born-Oppenheimer energy, as well as its derivatives with respect to ionic positions (forces) and cell parameters (stress tensor) in supercells, which can be provided, for instance, by first principles density-functional-theory codes. The method performs crystal geometry relaxation on the quantum free energy landscape, optimizing the free energy with respect to all degrees of freedom of the crystal structure. It can be used to determine the phase diagram of any crystal at finite temperature. It enables the calculation of phase boundaries for both first-order and second-order phase transitions from the Hessian of the free energy. Finally, the code can also compute the anharmonic phonon spectra, including the phonon linewidths, as well as phonon spectral functions. We review the theoretical framework of the stochastic self-consistent harmonic approximation and its dynamical extension, making particular emphasis on the physical interpretation of the variables present in the theory that can enlighten the comparison with any other anharmonic theory. A modular and flexible Python environment is used for the implementation, which allows for a clean interaction with other packages. We briefly present a toy-model calculation to illustrate the potential of the code.

The efficient and accurate calculation of how ionic quantum and thermal fluctuations impact the free energy of a crystal, its atomic structure, and phonon spectrum is one of the main challenges of solid state physics, especially when strong anharmonicy invalidates any perturbative approach. To tackle this problem, we present the implementation on a modular Python code of the stochastic self-consistent harmonic approximation method. This technique rigorously describes the full thermodyamics of crystals accounting for nuclear quantum and thermal anharmonic fluctuations. The approach requires the evaluation of the Born-Oppenheimer energy, as well as its derivatives with respect to ionic positions (forces) and cell parameters (stress tensor) in supercells, which can be provided, for instance, by first principles density-functional-theory codes. The method performs crystal geometry relaxation on the quantum free energy landscape, optimizing the free energy with respect to all degrees of freedom of the crystal structure. It can be used to determine the phase diagram of any crystal at finite temperature. It enables the calculation of phase boundaries for both first-order and second-order phase transitions from the Hessian of the free energy. Finally, the code can also compute the anharmonic phonon spectra, including the phonon linewidths, as well as phonon spectral functions. We review the theoretical framework of the stochastic self-consistent harmonic approximation and its dynamical extension, making particular emphasis on the physical interpretation of the variables present in the theory that can enlighten the comparison with any other anharmonic theory. A modular and flexible Python environment is used for the implementation, which allows for a clean interaction with other packages. We briefly present a toy-model calculation to illustrate the potential of the code. Several applications of the method in superconducting hydrides, charge-density-wave materials, and thermoelectric compounds are also reviewed.

I. INTRODUCTION
Ions fluctuate at any temperature in matter, also at zero kelvin due to the quantum zero-point motion. Even if the energy of ionic fluctuations is considerably smaller than the electronic one, many physical and chemical properties of materials and molecules cannot be understood without considering ionic vibrations. Since ionic vibrations are excited at much lower temperatures than electrons, ionic fluctuations are mainly responsible for the temperature dependence of thermodynamic properties of materials. They also determine heat and electrical transport through the electron-phonon and/or phononphonon interactions, as well as spectroscopic signatures detected in infrared, Raman, and inelastic x-ray or neu-tron scattering experiments. The large computational power available today has paved the way to material design and characterization, but advanced and reliable methods that accurately calculate vibrational properties of materials in the limit of strong quantum anharmonicity and that are easily interfaced with modern ab initio codes are required for accurately describing materials' properties in silico.
Since electrons are faster than ions, the ionic motion is assumed to be described by the Born-Oppenheimer potential V (R), which, at an ionic configuration R, is given by the electronic ground state energy. In the standard harmonic approximation V (R) is Taylor-expanded up to second-order around the R 0 ionic positions that minimize V (R). The resulting Hamiltonian is exactly diagonalizable in terms of phonons, the quanta of vibrations. Harmonic phonons are well-defined quasiparticles whith an infinite lifetime, which energies do not depend on temperature. These two features are intrinsic failures of this approximation: phonons acquire a finite lifetime due to their anharmonic interaction with arXiv:2103.03973v1 [cond-mat.mtrl-sci] 5 Mar 2021 other phonons (also because of other types of interactions such as the electron-phonon coupling), and phonon energies do depend on temperature experimentally. When higher-order anharmonic terms are small compared to harmonic ones, anharmonicity can be treated within perturbation theory [1][2][3]. Even if within perturbative approaches phonons' temperature dependence and lifetimes can be understood, whenever anharmonic terms of the V (R) potential are similar or larger than the harmonic terms in the range sampled by the ionic fluctuations, perturbative approaches collapse and are not valid [4]. This is often the case when light ions are present, as well as when the system is close to melting or a displacive phase transition, such as a ferroeletric or charge-density wave (CDW) instability.
In order to calculate from first principles vibrational properties of solids beyond perturbation theory and overcome these difficulties, several methods have been developed in the last years . Many of them are based on extracting renormalized phonon frequencies from ab initio molecular dynamics (AIMD) through velocity autocorrelation functions [6][7][8][9] or by extracting effective force constants from the AIMD trajectory [10][11][12]. In order to include quantum effects on the ionic motion, which are neglected on AIMD, the AIMD trajectory may be substituted by a path-integral molecular dynamics (PIMD) one [13]. Other methods are based on variational principles [16,18,[22][23][24]31], which are mainly inspired on the self-consistent harmonic approximation [32] or vibrational self-consistent field [33] theories, and yield free energies and/or phonon frequencies corrected by anharmonicity non-perturbatively.
Even if these methods have often successfully incorporated the effect of anharmonicity beyond perturbation theory in different materials, they usually lack a consistent procedure that prevents them from capturing properly both quantum effects and anharmonicity in the compound. For instance, many of them simply correct the free energy and/or the phonon frequencies assuming that the ions remain fixed at the R 0 classical positions. However, as it has been shown recently in several compounds [34][35][36][37], the ionic positions can be strongly altered by quntum effects and anharmonicity even at zero Kelvin. The structural changes are important for both internal degrees of freedom (the Wyckoff positions), and the lattice parameters themselves. Moreover, in many of the aforementioned methods, it is not clear what the meaning of the renormalized phonon frequencies is, i.e., whether they are auxiliary phonon frequencies intrinsic to the devised theoretical framework or if they really represent the physical vibrational excitations probed experimentally.
The stochastic self-consistent harmonic approximation (SSCHA) [22][23][24] is a unique method that provides a full and complete way of incorporating ionic quantum and anharmonic effects on materials' properties without approximating the V (R) potential. The SSCHA is defined from a rigorous variational method that directly yields the anharmonic free energy. It can optimize completely the crystal structure, including both internal and lattice degrees of freedom, accounting for the quantum nature of the ions at any target pressure or temperature. It computes thermal expansion even in highly anharmonic crystals. Furthermore, the SSCHA provides a well-defined approach to estimate at which thermodynamic conditions displacive second-order phase transitions occur. This is particularly challenging in ab initio molecular dynamics simulations, both for the dynamical slowing down that may hamper the thermalization close to the critical point, and for the difficulties in resolving the two distinct phases that continuously transform one into the other. Also, the rigorous theoretical approach of the SSCHA yields a clear distinction between auxiliary phonons of the theory and the phonon spectra probed experimentally, which can be accessed from a rigorous dynamical extension of the theory [23,38,39]. Lastly, the code provides non-perturbative third-and fourth-order phonon-phonon scattering matrices that can be fed in any external thermal transport code to compute thermal conductivity and lattice transport properties. Here, we present an implementation of the full SSCHA theory in a modular Python software that can be easily and efficiently interfaced with any total-energy-force engine, e.g., density-functional-theory (DFT) first-principles codes. This paper is organized to introduce the reader to the SSCHA algorithm and to review the recent developments in the SSCHA theory that lead to the SSCHA code, following the typical usage of the final user. In Sec. II we give a simple overview of the method, presenting a simple picture of how it works with a model calculation on a highly anharmonic system with one particle in one dimension. Then, we review the full theory of the SSCHA in details, starting from the free energy calculation and structure optimization in Sec. III. Then, we describe, in Sec. IV, the post-processing features of the code, which include calculations of the free energy Hessian for secondorder phase transitions, as well as phonon spectral function and linewidth calculations. Each section is introduced by an overview of the theory to understand what the code is doing and, then, reports the details of the implementation, together with a guide for setting up a typical run. In Sec. V, specific details of the Python code are provided, including the different execution modes and installation tips. As a showcase of the SSCHA, we provide a simple example in a thermoelectric material in Sec. VI (SnTe), where we fully characterize the thermodynamics of the phase transition between the high-symmetry and low-symmetry phases. This is also a guide on how to correctly analyze the output of the SSCHA calculation and the physical interpretation of the different frequencies.
In Sec. VII, we review some important results obtained so far with the SSCHA code. Finally, in Sec. VIII, we summarize the main conclusions.

II. THE VARIATIONAL FREE ENERGY
The SSCHA is a theory that aims at describing the thermodynamics of a crystal, fully accounting for quantum, thermal, and anharmonic effects of nuclei within the Born-Oppenheimer approximation. The basis of all equilibrium thermodynamics is that a system in equilibrium at fixed volume, temperature, and number of particles is at the minimum of the free energy. The free energy is expressed by the sum of the internal energy E, which includes the energy of the interaction between the particles (kinetic and potential), and the product between the temperature T and the entropy S, which accounts for "disorder" and is related to the number of microstates corresponding to the same macrostate of the system: In a classical picture, the free energy can be thus expressed in terms of the microscopic states of the system, which are determined by the classical probability distribution of atoms ρ cla (R). We remind that R is a vector of coordinates of all atoms in the system (we will use bold symbols to denote vectors and tensors in component free notation). The same holds for a quantum system, but we need to account also for quantum interference. This is achieved by calculating the free energy with the manybody density matrix. As the system at equilibrium is at the minimum of the free energy, the Gibbs-Bogoliubov variational principle [40] states that between all possible trial density matricesρ, the true free energy of the system is reached at the minimum of the functional F[ρ]: where is the total energy (K is the kinetic energy operator and V (R) the potential energy), and and S[ρ] the entropy calculated with the trial density matrix. · ρ = T r[ρ ·] indicates the quantum average of the the operator · taken withρ. If we pick any trial density matrixρ, F[ρ] is an upper bound of the true free energy of the system. The SSCHA follows this principle: we optimize a trial density matrix ρ to minimize the free energy functional F[ρ] of Eq. (2). Performing the optimization on any possible trial density matrix is, however, an unfeasible task due its many-body character that hinders an efficient parametrization. This is true also for a classical system: no exact parametrization of ρ cla (R) can be obtained in a computer with a finite memory.
The SSCHA solves the problem by imposing a constraint on the density matrix. In particular, the quantum probability distribution function that the SSCHA density matrix defines,ρ R,Φ (R) = R|ρ R,Φ |R , is a Gaussian.ρ R,Φ (R) is the quantum analogue of ρ cla (R), and determines the probability to find the atoms in the configuration R. The trial SSCHA density matrixρ R,Φ is uniquely identified by the average atomic positions (centroids) R and the quantum-thermal fluctuations around them Φ (we have explicitly expressed the dependence of ρ on R and Φ by adding them as subindexes), just like any Gaussian is defined by the average and mean square displacements. Within the SSCHA, we optimize R and Φ to minimize the free energy of the system. In this way, we compress the memory requested to storeρ R,Φ , as R depends only on 3N a numbers (the coordinates of the atoms), while the fluctuations Φ are encoded in a symmetric square real matrix of 3N a × 3N a . N a is the total number of atoms in the system. The free parameters in R and Φ can be further reduced by exploiting translation and point group symmetries of the crystal, resulting in an efficient and compact representation of the density matrixρ R,Φ .
The "harmonic" in the SSCHA name comes from the fact that any Gaussian density matrix that describes a physical system is the equilibrium solution of a particular harmonic Hamiltonian. Therefore, there is a one-to-one mapping between the trial density matrixρ R,Φ and an auxiliary trial harmonic Hamiltonian H R,Φ : Here, R is a real vector and Φ a real matrix that parametrize the trial Hamiltonian, while K and R are quantum operators that measure the kinetic energy and the position of the state. For simplicity, unless otherwise specified, all indices a, b, etc. run over both atomic and Cartesian coordinates from 1 to 3N a . Let us note here, that, inspired by the harmonic shape of H R,Φ , we will also refer to Φ as the auxiliary force constants. This mapping with a harmonic Hamiltonian is very useful, as both K ρ R,Φ and S[ρ R,Φ ] become simply the kinetic energy and entropy of the auxiliary harmonic system H R,Φ , which are analytic functions of Φ. Hence, the only quantity that we really need to compute is the average over the interacting Born-Oppenheimer potential The potential V (R) is the Born-Oppenheimer energy landscape, and can be easily computed ab initio by any DFT code (or by any energy and force engine). The SSCHA algorithm starts with an initial guess on R and Φ, and proceeds as follows: • Use the trial Gaussian probability distribution functionρ R,Φ (R) to extract an ensemble of random nuclear configurations in a supercell.
• For each nuclear configuration in the ensemble, compute total energies and forces with an external code, either ab initio or via a force field.
• Use total energy and forces on the ensemble to compute the free energy functional and its derivatives with respect to the free parameters of our distribution R and Φ.
• Update R and Φ to minimize the free energy.
These steps are repeated until the minimum of the free energy is found.
To illustrate better the philosophy of the method, we report in Figure 1 a simple application of the SSCHA to a one particle in one dimension at T = 0 K. In panel a, we plot the very anharmonic "Born-Oppenheimer" (BO) energy landscape V (R) of our one-dimensional particle (of mass of an electron). In Hartree atomic units it is given by We first study the classical harmonic result, obtained by Taylor-expanding the potential in Eq. (6) to second order around the minimum R 0 . Then, we use the harmonic solution to build our initial guess for the SSCHA density matrixρ R,Φ and update the parameters (Φ and R) until we reach the minimum of the free energy. In Figure 1(a), we compare the average atomic position and equilibrium free energy obtained with the harmonic approximation, with the SCHA, and the result we obtained with the exact diagonalization of the potential. While the harmonic result clearly overestimates the energy and yields an average atomic position far from the exact result, the SSCHA energy and average position are very close to the exact solution. In Figure 1(b), we report the probability distribution functions of the particle for the different approximations compared with the exact result. By definition, both the harmonic and SSCHA results have Gaussian probability distributions. However, while the harmonic solution is centered in the minimum of the BO energy landscape (and the width is fixed by the harmonic frequencies), the SSCHA distribution is optimized to minimize the free energy. Notice how, even if the exact equilibrium distribution deviates from the Gaussian line-shape, the SSCHA energy and average nuclear position match almost perfectly the exact solution as stated above. The very good result on the free energy reflects that the SSCHA error is variational: the free energy of the exact density matrix is the minimum. This means that the free energy is stationary around the exact solution, assuring that even an approximate density matrix (like the SSCHA solution) describes very well the exact free energy. This is an excellent feature of the SCHA, as the free energy and its derivatives fully characterize thermodynamic properties. Even if this simple calculation is performed at T = 0 K, the SSCHA can simulate any finite temperature by mixing quantum and thermal fluctuations on the nuclear distribution. The previously outlined straightforward implementation of the SSCHA becomes too cumbersome on a real system composed of many particles, especially if ab initio methods are used to extract V (R). The reason is that at any minimization step we need to calculate total energies and forces for many ionic configurations with displaced atoms in a supercell. The bottleneck is the computational cost of the force engine adopted. In the next sections of the paper we will show how the number of force calculations can be minimized and how these issues can be overcome by the code implementation proposed here. The resulting SSCHA code is very efficient, and, in most of the core cases, much faster than standard AIMD, with the advantage of fully accounting for the quantum nature of nuclei.

III. STRUCTURE RELAXATION AND FREE ENERGY MINIMIZATION
In this section we explain the simplest and most common use of the code: the calculation of the free energy and the optimization of a structure by fully accounting for temperature and quantum effects. This enables the simulation of finite temperature and pressure phasediagrams (with first order boundaries), as well as the calculation of the lattice thermal expansion. We start by briefly reviewing the theory of the SSCHA method. Then, we will explain the details of the implementation, giving tips on how to run a simulation.

A. The SSCHA free energy minimization
In the simplest and most standard usage, the SSCHA free energy functional that is minimized depends on the centroid positions R and the auxiliary force constants Φ as Here, we explicit that the entropy S ion only accounts for ionic degrees of freedom (not electronic). After the SS-CHA minimization, the final estimate of the equilibrium free energy is given by Therefore, the final result of a SSCHA free energy calculation is in general given in terms of the equilibrium configuration R eq , the free energy F , and the SSCHA auxiliary force constants Φ eq . The final free energy accounts for quantum and thermal ionic fluctuations without approximating the BO energy surface, valid to study thermodynamic properties, and the R eq positions determine the most probable atomic positions also taken into account quantum/thermal flluctuations 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 parametrizing the density matrix. The relation Figure 1. Illustration of the SSCHA method to a one dimensional particle problem at T = 0 K. Panel a: The one dimensional Born-Oppenheimer energy landscape V (R) as a function of the particle position R. The points represent the solution of the Harmonic approximation, the SSCHA, and the exact solution. The y coordinate of the points are the quantum total energy (including the zero-point motion), while the x axis coordinate is the average position of the particle. The SSCHA outperforms the harmonic approximation and it is very close to the exact solution. Panel b: Representation of the nuclear quantum distribution functions in the different approaches. The arrows point the average position of the particle in each distribution. Both harmonic and the SCHA are Gaussians, while the exact solution is more complex. The harmonic solution is centered around the minimum of the energy landscape R0, while the SSCHA centroid position R and width are optimized to satisfy the least energy principle. The average position in the exact case is however obtained as R ρ exact .
of it with the physical phonon frequencies is discussed in Sec. IV C. 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 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 which can be used, for instance, to estimate the structural changes imposed by pressure by fully accounting for fluctuations. As made explicit in Eq. (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 Eq. (7) the electronic ground state energy V (R) with the finite-temperature electronic free energy F el (R) = E el (R) − T S 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 SS-CHA method minimizes the functional 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 distributioñ ρ R,Φ ) 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].

B. The implementation of the free energy minimization
In the SSCHA code, the minimization of F[R, Φ] 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 R and the auxiliary force constants Φ. The partial derivatives are evaluated through the exact analytic formulas ΦR, ΦR 3rd and 4th order SSCHA force constants for a given R Eqs. (53), (54) Dynamical matrix of H (S) , given by D Req Pag. 14-1

Symbols indicating
Here f (BO) (R) are the Born-Oppenheimer forces that act on the ions when they are in the R positions; f H R,Φ (R) is the force given by the auxiliary harmonic Hamiltonian and Ψ is the displacement-displacement correlation matrix where with u = R − R we indicate the displacement from the average atomic position. Explicitly, In Eq. (16), ω µ and e µ are the eigenvalues and eigenvectors of the mass rescaled auxiliary force constants Φ ab / √ 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 H R,Φ . In other words, they are only used to define the trial density matrixρ R,Φ . We show how to compute the physical anharmonic phonon frequencies of the system in Sec. IV.
It is convenient to give an explicit expression for the gradient of the free energy with respect to the auxiliary force constants in terms of the ω µ eigenvalues and e µ eigenvectors. As shown in Ref. [23] (see Appendix B), the gradient can be rewritten as where Here, n µ = 1/(e β ωµ −1). The reason why we have introduced the Λ[0] tensor will be evident in Sec. IV. Even if Eq. (17) looks different to the gradient introduced in the original SSCHA work in Ref. [22], it can be demonstrated that both expressions are equivalent by simply playing with the permutation symmetry of . However, Eq. (18) unambiguously determines the value taken by the Λ[0] tensor for the ω ν = ω µ case, while the gradient in Ref. [22] did not describe explicitly what to do in this degenerate limit.
At the end of the SSCHA optimization, apart from the temperature-dependent R eq positions and the equilibrium auxiliary force constant 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 ε: 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 Born-Oppenheimer 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 Eq. (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 Eq. (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 Eqs. (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 Sec. III B 1).

The stochastic sampling
The stochastic nature of the SSCHA comes from the Monte Carlo evaluation of the averages in Eqs. (12), (13), and (19). A set of random ionic configurations are created in a chosen supercell according to the Gaussian ionic probability distributioñ 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: Here, N c is the total number of configurations in the ensemble, while R {j} is the j-th ionic randomly displaced configuration. Each of the R {j} configurations is generated according to the initial trial ionic distributioñ ρ R (0) ,Φ (0) (R) from which the minimization starts. To improve the stochastic accuracy, for each R {j} configuration also −R {j} is created, benefiting fromρ R,Φ (R) = ρ R,Φ (−R) property of the Gaussian distribution. The ρ j weights are computed and updated along the free energy minimization as the values of R and Φ change: At the beginning, when the ensemble has just been generated and R = R (0) and Φ = Φ (0) , all values of ρ j = 1. However, as the 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.

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 R and Φ are updated as In a perfectly quadratic landscape, this algorithm assures the convergence in just one step if both λ Φ and λ R are set equal to one. However, in order to avoid too big steps in the minimization, often it is more convenient to chose λ Φ|R < 1. This algorithm, with Hessian matrices 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 λ R re-scaled to the maximum eigenvalue of the ∂ 2 F ∂Φ 2 and ∂ 2 F ∂R 2 Hessian matrices, respectively, in order to have adimensional values independent on the system.
The preconditioning Hessian matrices that multiplies the gradients in Eq. (27) and Eq. (28) are approximated by the code. Following the procedure introduced in Ref. [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: and Eq. (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 1 2 ∂Ψ ab ∂Φ cd term in Eq. (13), the resulting update of the variational parameters at each step in the minimization is performed as and 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 Eq. (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 Eq. (27), it updates a root of Φ: The updating direction G n depends depends on the root order n: where with the · we indicate a matrix product. Similarly, We select the positive definite root matrix. Indeed, after the step of Eq. (29) the original force constant matrix is obtained as Thanks to the definition in Eq. (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 Eq. (30) with the preconditioned direction in Eq. (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 √ Φ and 4 √ Φ cases. However, we noticed that in systems with many atoms, the Hessian matrix calculation becomes the bottleneck as it scales with N 6 a . 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 Ref. [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.

The lattice geometry optimization
The lattice degrees of freedom {a i } are relaxed only after the minimization of the free energy with respect to R and Φ at a constant volume stops (see Sec. III B 4 for a detailed description of the stopping criteria). For this reason, the lattice geometry optimization is an "outer" optimization: at each step the lattice geometry optimization, we perform a full free energy minimization with respect to the centroids R and auxiliary force constants Φ. 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 Eq. (19), and generates a strain for the lattice as where P * is the target pressure of the relaxation and δ αβ is the Kronecker delta. The lattice parameters {a i } are updated as where λ {ai} 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 λ {ai} . In an isotropic material with a constant bulk modulus the optimal value of the step is B 0 is an input parameter given in GPa units. Good values of B 0 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 0 produces smaller steps in the cell parameters. The user can estimate the optimal value of B 0 to assure the fastest convergence by manually computing it from Eq. (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 Tr [ε] = 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.

The code flowchart
To start a SSCHA simulation, we need a starting guess on the trial positive definite force constants matrix Φ (0) and on the average atomic positions 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 Born-Oppenheimer energy landscape and the harmonic force constants 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 force constants 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 Then, the first random ensemble (that we call population in the SSCHA language) can be generated. For each configuration inside the population, its total Born-Oppenheimer 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 Sec. V). Once the Born-Oppenheimer 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 Sec. III B 2 and the minimization continues either until the stochastic sampling is not good or the algorithm converges. If R and Φ change a lot during the minimization of the free energy, the original ensemble no longer describes well the new probability distributionρ R,Φ (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 ef f : We halt the minimization when the ratio between N ef f and the number of configurations N c is lower than a parameter η defined by the user: 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 R and Φ are lower than a given threshold: 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 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: If Eq. (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 R eq , the equilibrium auxiliary force constant 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 Born-Oppenheimer 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 Sec. VII are converged within 3 populations, and the CPU time required to minimize each population is few minutes on a single CPU of modern laptops. 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 R are kept fixed in the minimization, the SSCHA self-consistent equation offers an alternative way of implementing the SSCHA theory (see Ref. [23] for a proof of Eq. (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 Φ(R), which must equal the result of the average. As in this approach the centroid positions are not optimized, the obtained auxiliary force constant matrix depends parametrically on R.
The self-consistent equation can be implemented stochastically, following the procedure outlined in Sec. III B 1. By using integration by parts [23], the righthand-side of Eq. (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 SS-CHA density matrix. An equivalent approach [30] is to extract the auxiliary force constants by fitting the obtained forces in the supercells generated with the SS-CHA density matrix to Eq. (14). This approach has Get the update direction for the matrix Root algorithm ?
Yes No Figure 2. Flowchart of the SSCHA code. The most time consuming part of the diagram is the ab initio calculation of the Born-Oppenheimer forces, energies, and stress tensors for all the configurations inside the ensemble, and it is shaded in red. All the other steps usually take few seconds when executed on a standard workstation, even in systems that contain several hundreds of atoms.
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 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 0 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 force constants but also all degrees of freedom in the crystal structure, the workflow outlined in Sec. III B provides a full picture of the effect of quantum/thermal fluctuations as well as anharmonic-ity on crystals, much more efficient than the approaches based on Eq. (45).

IV. POST-MINIMIZATION TOOLS: POSITIONAL FREE ENERGY HESSIAN, PHONON SPECTRAL FUNCTIONS, AND PHONON LINEWIDTHS
In the previous section we described how to compute the free energy of a system and fully optimize its structure by taking into account the anharmonicity that arises from both thermal and quantum fluctuations. After the free energy functional minimization, additional information can be extracted from the results obtained, namely, the second derivative (Hessian) of the positional free energy with respect to the centroids, the anharmonic phonon spectral functions, and the anharmonic frequency linewidths and shifts. In the next subsections we will ex-plain why these quantities are of physical interest and what is the strategy adopted by the code to computing them. The theory here reviewed was introduced in Ref. [23] and extensively applied for the first time in ab initio calculations to H 3 S in Ref. [48].

A. Positional free energy Hessian
As shown in Sec. II, for a given temperature the free energy at equilibrium F of a system with Hamiltonian H is obtained by minimizing the density-matrix functional where ρ is the equilibrium density matrix of the system obtained at the minimum. The average atomic positions at equilibrium are R ρ = R eq . By minimizing the functional keeping fixed the average atomic positions in a generic configuration R, R ρ = R, we define the positional free energy F (R): where ρ R is the density matrix giving the constrained minimum for the considered average position R. Since R and F (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 R − R hs , where R hs is the average position of the atoms when the system is in the highsymmetry 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 (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][50][51][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 R eq , is a high-symmetry configuration R hs . Therefore, at high temperature the free energy Hessian in R hs is positive definite, i.e. its eigenvalues are positive. As the temperature decreases, the minimum in R 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 R 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, 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ρ R,Φ coincides with the centroid parameter R, Thus, within the SSCHA the positional free energy is obtained by minimizing the SSCHA free energy functional F[R, Φ] with respect to the trial quadratic amplitude Φ only: The auxiliary force constants that minimize Eq. (50) for a given R position of the centroids will be labeled in the following as Φ R . Solving Eq. (50) allows to employ the SSCHA code to have direct access to F (R) for any 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 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 1). Using this notation, we can express the positional free energy Hes- where M ab = δ ab M a is the mass matrix, and Λ R [0] is the z = 0 value of the fourth-order tensor Λ R [z], already introduced in Eq. (18). In the equations above the quantum statistical averages are taken with ρ R,Φ R , which for a given R position of the centroids is taken with the Φ R auxiliary force constants that minimize the free energy. Λ R [z] is given in components by where ω 2 µ and e a ν are the eigenvalues and eigenvectors of D R , and The only difference between Λ R [0] and Λ[0] (introduced in Eq. (18)) is that in the former the eigenvalues and eigenvectors entering the equation are those associated to the Φ R auxiliary force constants at the centroid positions R, while in the latter this is not necessarily the case. The subindex R in the equations above precisely indicates that the averages are calculated with a density matrix defined by R and Φ R (after a full SCHA relaxation at fixed nuclei position R). We will refer to the (n) Φ R tensors as the nth-order SSCHA force constants (FCs). Note that for the second-order we drop the (2) upper index.
The SSCHA code computes the free energy positional Hessian through Eq. (51). At the end of a SSCHA free energy functional minimization, the SSCHA matrix Eq. (52), with its eigenvectors and eigenvalues, is available. Thus, Λ R [0] is readily computable and the only quantities that need some effort to be calculated are the averages of Eqs. (53) and (54). The code computes them through these equivalent expressions (obtained by integrating by parts): These averages are computed employing the stochastic approach already described in Sec. III B 1 (indeed, as explained in Ref. [23], the choice of Eq. (58), among other possible alternatives, aims at reducing the statistical noise). Note that if the calculation of the free energy Hessian is performed at R eq , f (BO) (R) ρ R,Φ 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 Φ 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 Φ R , an importance sampling reweighting has to be employed in order to evaluate the averages · ρ R,Φ 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 Eq. (51), i.e. discarding the terms coming from (4) D 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 Using Eq. (51), or its approximated expression Eq. (59), the SSCHA code can compute the Hessian of the free energy at any R. However, as said, its most significant usage is in R eq , due to its relevance to characterize displacive second-order phase transitions. In this case, Eq. (51) can be written in a quite explanatory form. At the end of a full SSCHA minimization, the obtained R eq and Φ eq define the so-called SSCHA effective harmonic Hamiltonian 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 (S) ab = (Φ eq ) ab / √ M a M b of the SSCHA Hamiltonian H (S) , the anharmonic generalization of the dynamical matrix D (F) can be written as where D eq (62) is the static SSCHA self-energy (the reason behind the use of this name will be clear in the Sec. IV C). In particular, in the bubble approximation we have where (B) D eq : Λ eq [0] : is the so called "bubble" static self-energy.
In conclusion, after the SSCHA minimization, the code allows to compute the high-order SSCHA force constants, Eqs. (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 E 1 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 possbile 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 employ of this method are given for H 3 S in Fig. 3 of Ref. [48], with the softening of an optical mode driven by pressure release, and for SnSe in Fig. 2 of Ref. [54], with the softening of the distortion mode obtained by decreasing the temperature.
B. Static bubble self-energy calculation: improved free energy Hessian calculation 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 1 , l 2 ) and (3) D eq (l 1 , l 2 , l 3 ) are available. Here, D (S) ab (l 1 , l 2 ) is the real space D (S) matrix in which we made explicit the dependence of the lattice vectors l 1 and l 2 that identify the unit cells in which atom a and b are located, respectively. Using them, the code allows to compute the static bubble  D(l 1 , l 2 , 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 Eq. (66) can be arbitrarly increased up to convergence, so as to reach the thermodynamic limit in the evaluation of the bubble static self-energy. Second, from (B) Π(q, 0) and D (S) (q), through Eq. (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 Ω 2 µ (q) are the eigenvalues of D (F) (q). This is, for example, the procedure followed to compute the (static) SCHA phonon dispersions of NbS 2 shown in Figs. 2, 3 of Ref. [51], and to compute the interpolation-based convergence analysis shown in Fig. 3 of Ref. [55] for TiSe 2 monolayer.
C. Dynamic bubble self-energy calculation: spectral functions, phonon linewidth and shift 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 Refs. [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 √ M a (R a −R a eq ). This ansatz has been rigorously proved within the Time Dependent Self-Consistent Harmonic Approximation (TD-SCHA) [38,39]. In this dynamical theory where D (S) is the dynamical matrix of the SSCHA effective harmonic Hamiltonian H (S) , and Π(z) is the SSCHA self-energy, in general given by D eq : Λ eq (z) : 1 − D eq : Λ eq (z) and in the bubble approximation by D eq : Λ eq (z) : 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 R eq . With the Green function we obtain the spectral function σ(Ω) = −2 ImTr [G(Ω + i0 + )], 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 or 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 where in the last line we have used Eq. (65). Therefore, where Ω 2 µ (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 Eqs. (64), (66), and (69))) 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 Eq. (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 Eq. (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 Eq. (71)) 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 (B) Π µν (q, 0) (see Eq. (72)). Indeed, in this case the self-energy is real (Hermitian) and the computed spectral function becomes where Ω 2 µ (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, Eq. (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 Eq. (76) the full bubble selfenergy with its static value computed through Eq. (66). This can be used to double-check that, as expected from Eqs. (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 Eq. (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, Eq. (73), the total spectral function is given by the superposition of individual mode spectral functions: where now the (q, µ)-mode spectral function σ µ (q, Ω) is computed with and Z µ (q, Ω) = ω 2 µ (q) + Π µµ (q, Ω + iδ se ) .
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 Eq. (79) resembles a Lorentzian, but with frequencydependent 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), 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 Γ µ (q) = −ImZ µ (q, Ω µ (q)).
Once the dynamical self-energy and the Z µ (q, Ω) are computed, the SSCHA code allows to compute the singlemode spectral functions in the Lorentzian approximation, estimating the frequency Ω µ (q) and HWHMs Γ µ (q) in different ways. One, optional, possibility is to solve selfconsistently Eq. (82) to estimate Ω µ (q), and then Γ µ (q), through Eq. (83). However, by default, the "one-shot" approximation is employed with 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 Eq. (85), we get: This perturbative approach is also employed by the SS-CHA code to evaluate the quasiparticles' energies and lifetimes. Examples of spectral function calculations done with Eq. (76), Eqs. (78)- (80) and Eq. (81) can be found in Fig. 4 of Ref. [48]. In Fig. 5 of the same reference, the anaharmonic phonon frequencies and linewidths along a path, computed using the Lorenztian approximation through Eqs. (82) and (83), are shown. The spectral function computed with Eq. (79) along a path is shown with a colorpolot in Fig. 3 of Ref. [49] for PbTe, and in Fig. 4 of Ref. [54] 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.

V. THE PYTHON CODE
Two different Python libraries are provided with the SSCHA code: CellConstructor and Python-sscha. The latter is the library that performs the SSCHA minimization itself, while the former is a library that deals with the dynamical matrix, the crystal structure, the symmetrization, and performs the calculation of phonon spectral functions and linewidths as a post-processing tool.
The SSCHA code allows to set up the calculations with a simple Python script. In the standard calculation, the script loads the starting dynamical matrices; sets up the ensemble and the parameters for the SSCHA run; performs the calculations of the Born-Oppenheimer energies, forces, and stress tensors on the configurations in the ensemble by calling to a external total-energy-force engine; and starts the minimization of the free energy. A simple input script that performs all these steps requires less than 20 lines. Examples are provided within the code, as well as step-by-step tutorials to perform a full SSCHA calculation starting just with the structure in a cif file. Python scripting the SSCHA run makes it versatile, as it can be interfaced with other Python libraries to facilitate the analysis of the results. As an alternative, it is also possible to write an input file and run the SSCHA code as a stand-alone command-line program.

A. Code structure
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 (Eq. 32), extract auxiliary phonon frequencies and polarization vectors, or interpolate them to other points in the Brillouin zone.
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 [56]. 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.

B. Parallelization
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 Born-Oppenheimer 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 this calculation 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 multithreading (OpenMP), multiprocessing (MPI), and GPU (CUDA) parallelization [57]. 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 [56], some of them speeded up with an explicit Fortran implementation. Thanks to the numpy implementation [58], 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 [59][60][61]).

C. Execution modes
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 totalenergy engine, e.g. a DFT code, to calculate their Born-Oppenheimer 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 [62], which already implements interfaces with most common ab initio codes like Quantum ESPRESSO [63,64], VASP [65], SIESTA [66], CP2K [67], and many more. Force-field codes like LAMMPS [68] 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.

D. Distribution
The packages 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.

VI. A MODEL CALCULATION ON TIN TELLURIDE
To display the potentiality of the code, we provide an example calculation on a SnTe toy model force field, where the lattice has been artificially stretched to enhance the anharmonicity. More details on this force field can be found in Ref. [23]. We provide it as a separate package under GPL license.
SnTe, as other ferroelectric materials [50,52], undergoes a displacive phase-transition, where a phonon mode at Γ softens with temperature lowering and provokes a cell distortion from the high-temperature high-symmetry F m3m phase to the low-temperature R3m phase. The toy model is able to reproduce this behavior, although it does not pretend to accurately describe the real SnTe transition and it is just provided as an artificial example.
The system has an incipient ferroelectric instability, marked by the negative curvature of the energy in the high-symmetry position. This means that an optical vibrational mode has an imaginary frequency at Γ within the harmonic approximation. The BO energy of the toy model as a function of the atomic displacements projected onto the eigenvecotrs of the imaginary mode (the order parameter ∆) is reported in Figure 3(a), where it is clear that the high-symmetry F m3m is not at the minimum of V (R). However, as extensively discussed above, the stability of a structure is determined by the temperature-dependent free energy, and not the BO potential. Notably, E is not the energy profile reported in Figure 3(a), as it also includes the vibrational contribution to the energy. For this reason, the energy profile (and the harmonic approximation) does not correctly describe even the behavior at T = 0, where there is no entropy contribution. Since entropy usually is higher in high-symmetric positions (∆ = 0), the F m3m high symmetry phase will become progressively more stable as temperature increases.
In Figure 4 we show the evolution of the free energy, its gradient, and the frequencies of the auxiliary force constants during a typical SSCHA minimization at T = 250 K for this system. We start the minimization from the harmonic solution of the high-symmetry phase F m3m, which has imaginary frequencies. Since the system is strongly anharmonic, the starting solution is very far from the solution. To approach the minimum quickly and with low computational effort, we start the minimization with a small number of configurations (here 50). Here, we need only three ensembles to converge to the minimum, as a zero gradient is obtained with a reasonably large value of η. To improve the quality of the calculation (and decrease the stochastic error) we further run two more populations with 100 and 200 configurations, both converging in one population. Both the gradient and the free energy rapidly decrease and the result converges. An extra population with 1000 configurations is included to see that the result is converged. Figure 4(c) presents the evolution of the auxiliary phonon frequencies associated to the auxiliary force constant matrix. The small change in these frequencies when the number of configurations is increased means that a small number of configurations is sufficient to have a good estimate of the auxiliary frequencies. Indeed, a good check for a well-converged result is to verify that these frequencies are stationary and do not change more in the minimization. The SSCHA code prints in output, if requested, this information at each run. We provide the code scripts that produce this kind of graphs from the raw data generated by the code, which facilitates the user to control if the minimization is working correctly.
In Figure 3(b) we report the order parameter obtained at the end of the SSCHA minimizations at different temperatures. The starting structure at low temperatures is the low-symmetry R3m. When temperature is increased, the structure obtained at the previous lower temperature is used as input. At low temperatures the output structure reamins the R3m, with ∆ = 0, but at T = 205 K, the low-symmetry phase jumps into the high-symmetry phase, marking a first-order phase transition. We can confirm it is a first-order phase transition as we can start cooling down from the high-symmetry phase (without constraining the new symmetries acquired) and the system remains stable up to T = 160 K, when it transforms back to the low-symmetry phase. This is the hysteresis  cycle of the material. We can further analyze the thermodynamic properties. The SSCHA provides also the free energies of the two phases. We compare them in Figure 3(c). As clearly shown, the low-symmetry phase is more stable up to 200 K, so that the phase diagram in this model is formed by the R3m phase below 200 K and the F m3m above. We can also see whether the F m3m becomes dynamically unstable by calculating the hessian matrix of the free energy. We plot the second derivative of the free energy with respect to the order parameter in Figure 3(d). The free energy curvature becomes negative below T = 154 K. This is a threshold below which the F m3m phase is no longer stable, and cannot exist or be observed. Consequently, it coincides with the lower bound of the hysteresis cycle. On the other side, looking at the free energy Hessian of the R3m phase, we see that the frequency of the mode along the order parameter softens to zero at T = 210 K, marking an upper bound to the stability of the low symmetry phase. Interestingly, while in the high symmetry phase the "bubble approximation" (Eq. 59) is very accurate, the correct estimation of the free energy Hessian in the R3m phase requires the full expression of the Hessian (Eq. (51)).
This example shows that the SSCHA can fully characterize a complex first-order phase transition, and thanks to the possibility of exploiting symmetries, we can even study a phase that is dynamically unstable, i.e., the F m3m below the critical point. We can do simulations directly in the high-symmetry phase, with a considerable gain in the computational cost, and spot instabilities by the Hessian matrix calculation, as in Figure 3(d).
However, we can do even more: finite temperature structure search. To investigate whether the R3m is the  actual ground state within the toy model or a lower symmetry phase is energetically favored, we calculate the Hessian also in the R3m phase. We find that in the whole region of the simulation, the R3m phase is dynamically unstable and the system wants to break the symmetry once again. To find the real ground state, we release all the symmetry constrains in our simulation and perform a full relaxation with the SSCHA at T = 100 K. We discovered a new phase of Cc symmetry defined in a 1×2×1 supercell of the original cubic cell. In Figure 5 we report the final phase diagram for the SnTe toy model. The new Cc phase is found to be the ground state up to 250 K, where the cubic F m3m becomes again energetically fa-vorable. The Cc continues to exists until 280 K, where it transforms into the F m3m phase. We want to remark that the particular temperature, phase transitions, as well as the real existence of phase Cc are just features of the toy model and do not pretend to represent the physics of this system. The real SnTe has a ferroelectric R3m ground state at low temperatures and the phase transition to the F m3m phase is of second-order type (the order parameter does not jump, and the free energy lines of the two phases touch when the F m3m mode becomes imaginary). This system has been already studied with the SSCHA in Ref. [49] with ab initio energies and forces. However, even if it is just R3m Cc Figure 5. Full phase diagram of the SnTe toy model. In dashed lines we report the unstable phases (whose free energy Hessian has an imaginary mode).
a toy model, this example shows how the code can reach high-symmetry phases in an unsupervised way starting from the low-temperature structure. For this reason, the SSCHA is also attractive for a structure search perspective: it is able to perform a search of saddle-point structures in the classical Born-Oppenheimer energy landscape that become the ground state due to ionic quantum and/or thermal fluctuations. This can be a great advantage to find saddle-point structures in complex systems with many atoms in the unit cell, such as in molecular crystals, where symmetry constrains may be inefficient [69].
The other post-processing utility the code provides is the calculation of spectral functions and dynamical phonon spectra. We remark that the auxiliary phonons, i.e. the eigenvalues of D (S) , are just an auxiliary quantity used to define the density matrix. For this reason, they only describe quantum flucuations around the centroid positions. The eigenvalues of the Hessian matrix D (F) , instead, are the response to a static external perturbation, and describe the stability of the structure with respect to a spontaneous symmetry breaking. Last, physical phonons, those observed by experimental probes like vibrational spectroscopy and inelastic scattering, must be computed from the dynamical interacting Green function. While all these definitions of phonon frequencies coincide in perfectly harmonic crystals, when anharmonicity is involved, they can differ significantly. The SSCHA code offers a tool to easily compute the dynamical Green functions as a post-processing utility as discussed in Sec. IV C.
In Figure 6, we plot the phonon spectrum, computed as the spectral function obtained from the dynamical Green functions. As anticipated, the peaks of the spectral function in Figure 6 do not coincide with the dispersion obtained from the auxiliary dynamical matrix D (S) , and show a rather anomalous behavior. It is worth mentioning that effective charges are considered in the calcula-tion of the spectral functions. The effective charges are considered following the procedure outlined in Appendix E 3. In Fig. 7 we illustrate the convergence for a phonon linewidth 2Γ µ (q) with respect to the δ se parameter and the k-mesh summation grid in Eq. (75).
All the data of this simulation has been obtained in less than one hour, using a single processor on a laptop, proving the high-efficiency of the SSCHA package, which is beyond standard molecular dynamics software. We provide in the additional materials the Python scripts to run and analyze all the simulations here reported in this example.

VII. APPLICATIONS OF THE SSCHA METHOD
In order to illustrate some physical problems and materials that have already been efficiently tackled with the SSCHA, in this section we briefly overview some of the systems studied with this method. One should not consider that the applications are limited to these examples, the SSCHA provides a general utility to treat accurately and efficiently all materials where ionic vibrations play a crucial role both in the thermodynamic and transport properties.

A. Hydrogen-based compounds
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, 70, 71]. 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 [72,73], is a consequence of a huge anharmonic renormalization of the phonons [21]. Also, the experimentally found high-temperature superconductivity in H 3 S around 200 K [74] and in LaH 10 around 250 K [75,76] at high pressures can only be explained if phonon frequencies renormalized by anharmonicity are considered in the superconductivity equations [34,36]. In Fig. 8 we show the huge anharmonic renormalization of the phonon frequencies for H 3 S [34,48]. Superconductivity in hydrogen compounds can be both largely suppressed but also enhanced by anharmonicity depending on the system [77].  The quantum effects and anharmonicity that the SS-CHA 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 (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 3 S and LaH 10 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 Fig. 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 10 .
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. The SSCHA auxiliary phonon frequencies are given, those obtained diagonalizing Φeq, together with the phonon frequencies obtained from the spectal function in the Lorentzian approximation. The linewidth obtained in the latter case is also given. The harmonic phonons are also shown (bottom panel). We also provide the structure of Im3m H3S. Data taken from Ref. [48]. (b) Crystal structures for LaH10 obtained from the minimum of the classical Born-Oppenheimer energy landscape, C2, and from the SSCHA quantum energy landscape, F m3m.

B. Charge density wave materials
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 [78][79][80]. 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,55,81,82].
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 Sec. IV. 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 Fig. 9 we show as an example the temperature dependence of the phonon spectra derived from the free energy Hessian in monolayer NbSe 2 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) D eq = 0 in Eq. (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 temepratures [41,55].
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 [55], the difference (or similarities) between the CDW transitions in bulk and the corresponding two-dimensional structures [51,55], as well as the importance of Van der Waals forces in the melting of CDW transitions [82]. Consequently the SS-CHA program is expected to have a large impact on theoretical studies of CDW transitions for many type of materials, not just TMDs.

C. Phase transitions, spectral functions, and thermal conductivity in semiconducting materials
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 (a) (b)  the high-temperature phases of these compounds are not a minimum of the Born-Oppenheimer 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 Sec. VI. For instance, in order to convincingly show that the transition between the high-temperature Cmcm phase of SnSe and the lowtemperature P nma 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 P nma phase relaxed at this temperature into the Cmcm, showing that the P nma 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 Fig. 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.
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 [83,84]. It has been shown that employing the SSCHA phonon scattering tensor (3) Φ R in the thermal transport calculations leads to a very good agreement with experimental results, in contrast with what obtained by employing the standard third-order derivatives of the Born-Oppenheimer total energy. The difference is that the former includes higher-order anharmonic terms coming from the average over the thermal ensemble (Eq. 53). Both in good thermoelectric SnSe and SnS compounds, higher order terms captured by (3) Φ 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 Born-Oppenheimer 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 thermolectric 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.

D. Other type of materials
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 [85]. The SSCHA has also been employed to understand the thermal expansion and the behavior of low-energy acoustic modes of graphene [86], 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.

VIII. CONCLUSIONS
We present here the implementation of the SSCHA theory into an efficient modular Python code, which can be run in conjunction with other Python modules and interfaced with HPC clusters for the Born-Oppenheimer total energy, force, and stress tensor calculations needed.
The SSCHA provides an efficient way of calculating the effect of ionic quantum and/or thermal fluctuations on the free energy, as well as their impact on the atomic positions. It is a unique feature of the SSCHA to optimize the atomic positions, including the lattice degrees of freedom, by considering quantum and finite temperature fluctuations and without any approximation on the Born-Oppenheimer energy landscape. As a postprocessing, it calculates the free energy Hessian, which allows to infer the thermodynamic conditions at which second-order phase transitions occur. Furthermore, it enables the evaluation the interacting phonon spectral functions, predicting the outcome of most common experimental techniques (IR and Raman spectroscopies, x-ray and neutron inelastic scattering). It can also extract the phonon linewidths from the Lorentzian approximation of the spectral functions, which can be later interfaced with any code that calculates the thermal conductivity.
In conclusion the SSCHA code provides a complete and efficient software for studying vibrational properties of materials, particularly suitable to study systems with prominent quantum and/or thermal fluctuations that are thus largely affected by anharmonicity, which can be applied to study many relevant problems in physics, chemistry, and material science.
Let's note that after this change of variables. Note that in Eq. (A1) we explicitly indicate the dependence of the operator O on the lattice parameters {a i }. Thus, in that equation, centroid positions R refer only to the internal degrees of freedom of the crystal structure. When calculating the stress tensor from Eq. (19) we are deriving the free energy functional in the minimum of the SSCHA free energy with respect to the auxiliary force constants Φ for given centroid positions R. Thus, the stress tensor should be calculated considering the derivatives (A6) The final equation of the strain should however be calculated for the minimum of the free energy with respect to the centroid positions, in which case the second addend above vanishes. Therefore we will just give the expression for the equilibrium situation: This expression coincides with the one usually employed to compute the stress tensor from the BO energy surface, but with the V (R) potential substituted by the anharonic free energy. Let us write the free energy at the minimum as where Φ R is the dynamical matrix that minimizes F fixing the average atomic positions.
The first term in Eq. (A8) does not give any contribution to the derivative (as it depends on R through Φ R , which minimizes already the free energy). Therefore, the only therm that survives in the stress tensor is Joining Eq. (A10) with Eq. (A2) we can compute the derivative of an average in the SSCHA ensemble with respect to the strain: Replacing O by the BO energy landscape V (R) we get The term with the harmonic potential V can be derived writing its explicit dependence on the strain tensor ε: (A15) From which we obtain Combining Eqs. (A13) and (A16) with the definition of the stress tensor, it is trivial to get Eq. (19).

Appendix B: Gradient equation
The gradient equation presented here in Eq. (13) can be obtained starting from the equation obtained in Ref. [23]. By using the result proved in the same reference, we have Anologously, Substituting Eq. (B3) and (B4) into Eq. (B1) we get Eq. (13). In Ref. [23] it was also shown that where A is a symmetric matrix. This also proves Eq. .

Appendix C: The Hessian preconditioner
In Ref. [24] it was shown that However, due to the relationship in Eq. (B5), we can see that we can effectively extend this equality to With the latter result, it is trivial to see how the preconditioned gradient that is used along the minimization can be written as in Eq. (27).

Appendix D: Symmetries
The original algorithm proposed to account for symmetry in Ref. [22] was based on the Gram-Schmidt orthonormalization of the symmetry generators. This algorithm is suited for systems with a reduced number of atoms in the unit cell, but scales as n 6 a , with n a the number of atoms in the unit cell. This symmetrization procedure becomes thus a real bottleneck of the SSCHA code for systems with more than 30 atoms in the unit cell. In the version of the code we describe here, the orthonormalized generators are not calculated and, instead, the starting dynamical matrix and the gradient are directly symmetrized. The symmetrization of the dynamical matrix, or its gradient, is made in q space, which allows for a very fast implementation even for big supercells.
The code enforces all the symmetries in the auxiliary force constant matrix as where S i are the 3×3 point group matrices of the space group, N S the number of symmetries of the crystal, TŜ(q) are unitary matrices that represent the S symmetry in the q point. These matrices are reported in Refs. [87][88][89]. To find the symmetries given the structure, we wrapped into the SSCHA the symmetry module of Quantum ESPRESSO [63,64]. This operation is performed also on the gradient of the dynamical matrix each time it is computed. Since the dynamical matrices satisfying the symmetries define a linear subspace, if both the gradient and the original dynamical matrix belong to this subspace, any linear combination of them will also satisfy the symmetry constrains. Thus, it is necessary to symmetrize the dynamical matrix once at the beginning, and then apply the symmetry constrains only to the gradient to preserve the symmetries in the whole simulation.
The symmetry module from Quantum ESPRESSO only recognizes symmetries when the unit cell is the primitive one. Sometimes, it could be convenient to choose a different unit cell. Therefore, we also interfaced the SSCHA code with the spglib package [90] to improve the identification of symmetries. Instead of working in the unit cell in q space, spglib provides the symmetry operations in real space. In this case, the SSCHA code divides the symmetry matrices identified by spglib into pure translations and point group operations. Then, symmetries are enforced in real space by first imposing pure translations, followed by point group operations as Then, the permutation symmetry on the indices is imposed. Finally, the code transforms back the real space dynamical matrix (or the gradient) in q space. This operation takes more time than the symmetrization in q space, as it is performed in the supercell. However, due to its simplicity and to avoid the cumbersome q-space symmetrization of higher-order force constants, the same supercell approach is used to symmetrize the third-and fourth-order force constant matrices, namely (3) Φ and (4) Φ introduced in Sec. III B. Symmetries are also enforced on the average positions of the nuclei (and the forces). After computing the SCHA forces on atom t along direction α, f α t , we impose symmetry as where S −1 (t) is the atom in which the S symmetry maps the t one to. In this way, the forces are correctly directed only along the Wyckoff coordinates, and the atomic positions relax subsequently keeping the correct Wyckoff positions.

Acoustic sum rule on the auxiliary force constants
Besides space group symmetries, also the acoustic sum rule (ASR) must be imposed. The acoustic sum rule 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 force constant matrix in the supercell: 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 acoustic sum rule define a linear subspace. Thus we define the orthogonal projector operator that enforces the acoustic sum rule as The projection matrix in real space is This operation only affects the dynamical matrix at Γ. The same projector is employed to impose the ASR on the forces: 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 E 2. The latter affects not just phonons at Γ, and, thus, it is more suited for interpolating dynamical matrices close to the Brillouin zone center.
Appendix E: Reciprocal space formalism and Fourier interpolation

Reciprocal space formalism
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 numer of unit cells forming the supercell). Thus, in this notation, in general we will have n-th 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,   Notice that, due to the lattice translation symmetry, D (S) (q 1 , . . . , 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 Eq. (61)) and its eigenvalues Ω µ (q) and eigenvectors f µ (q).

Fourier interpolation: centering and acoustic sum rule
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 arbtrariness 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 equivalentlattice 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 [91].
The weights are used to define the centering. Given a 3rd-order FCs, Φ α1α2α3 s1,s2,s3 (0, l 2 , l 3 ), its "centered" version (cent) Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) is given by 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 force constant 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 force constants 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 acoustic sum rule (ASR). For a nth-order FC, the acoustic sum rule is li si Φ α1...αi...αn s1...si...sn (l 1 , . . . , l i , . . . , l n ) = 0 The ASR comes is 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 acoustic sum rule 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 Eq. (E6), since only the simultaneous sum on n − 1 indices is zero: 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: Since the last factor, highlighted with a brace under, in general is not constant with respect to s 3 and l 3 , it cannot be factored out from the sums, so that the ASR for the original Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) s3 l3∈SC Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) = 0 (E13) cannot be used to obtain the ASR for the centered (cent) Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ). However, using Eq. (E4), with similar passages we can show that the sum over the last two indices is zero: thus the weak ASR is fulfilled. The spoling 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 fullfil 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 [91]. 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-sightdness of the centered FCs even after the ASR imposition). Given a centered Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ), the Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) that fulfills the ASR on the third index is computed with Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) = Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) − K α1α2α3 s1s2s3 (0, l 2 , l 3 |p) × s3,l3 where K α1α2α3 s1s2s3 (0, l 2 , l 3 |p) is the scaling factor , (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 0 = 1 has been adopted). The Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) defined through Eq. (E18) fulfills the ASR on the third index, since for any p ≥ 0 the scaling factor fulfills the normalization condition s3,l3 K α1α2α3 s1s2s3 (0, l 2 , l 3 |p) = 1 .
The value of p has effects on the way the different terms of Φ α1α2α3 s1s2s3 (0, l 2 , l 3 ) are scaled. For p = 0 the scaling factor is a pure geometric quantity related to the three atoms clusters. Indeed, given s 1 , s 2 , l 2 , the scaling factor K α1α2α3 s1s2s3 (0, l 2 , l 3 |p = 0) is fully determined (it is the same for all the α h , s 3 , l 3 ) and, in particular, it does not depend on the FCs value. On the contrary, for p = 0, given α h , s 1 , s 2 , l 2 we have so that if p > 1 the scaling factor is higher (lower) for FCs have lower (higher) absolute value, otherwise the opposite.

Effective charges
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 force constants can be calculated analytically since it is fully determined by the Born effective charges (Z * s ) αβ (effective charge tensor of atom s) and the electronic dielectric permittivity tensor ( ∞ ) αβ , which can both be calculated from first principles. For a given q ∈ BZ, this dipole-dipole contribution is given by [92,93] 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 η) [63]. Once Z * s 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 Brillouin zone, 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 [63]: The dipole-dipole correction to the FCs given by Eqs. (E22), (E23) is nonanalytic at zone center and its q → 0 limit depends on the direction q = q/ q along which the limit is performed: where 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 [63]. 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.