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

Journal of Physics: Condensed Matter TOPICAL REVIEW • OPEN ACCESS The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects To cite this article: Lorenzo Monacelli et al 2021 J. Phys.: Condens. Matter 33 363001 View the article online for updates and enhancements. You may also like Structure of disordered materials under ambient to extreme conditions revealed by synchrotron x-ray diffraction techniques at SPring-8—recent instrumentation and synergic collaboration with modelling and topological analyses Koji Ohara, Yohei Onodera, Motohiko Murakami et al. - Particle-laden fluid/fluid interfaces: physico-chemical foundations Eduardo Guzmán, Irene Abelenda-Núñez, Armando Maestro et al. - Surface-assisted fabrication of lowdimensional carbon-based nanoarchitectures Dong Han and Junfa Zhu - This content was downloaded from IP address 161.111.10.229 on 27/01/2022 at 08:31 Journal of Physics: Condensed Matter J. Phys.: Condens. Matter 33 (2021) 363001 (34pp) https://doi.org/10.1088/1361-648X/ac066b Topical Review The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects Lorenzo Monacelli1,∗, Raffaello Bianco2,∗, Marco Cherubini1,3, Matteo Calandra4,5,IonErrea 2,6,7,∗and Francesco Mauri1,∗ 1Dipartimento di Fisica, Universit` a di Roma Sapienza, Piazzale Aldo Moro 5, 00185 Roma, Italy 2Centro de Física de Materiales (CSIC-UPV/EHU), Manuel de Lardizabal pasealekua 5, 20018 Donostia/San Sebastián, Spain 3Center for Life NanoScience, Istituto Italiano di Tecnologia, Viale ReginaElena 291, 00161 Rome, Italy 4Sorbonne Universit´ e, CNRS, Institut des Nanosciences de Paris, UMR7588, F-75252 Paris, France 5Dipartimento di Fisica, Universitá di Trento, Via Sommarive 14, 38123 Povo, Italy 6Fisika Aplikatua Saila, Gipuzkoako Ingeniaritza Eskola, University of the Basque Country (UPV/EHU), Europa Plaza 1, 20018 Donostia/San Sebastián, Spain 7Donostia International Physics Center (DIPC), Manuel Lardizabal pasealekua 4, 20018 Donostia/San Sebastián, Spain E-mail: [email protected],[email protected],ion.[email protected] and [email protected] Received 9 March 2021, revised 13 May 2021 Accepted for publication 28 May 2021 Published 13 July 2021 Abstract 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 (SSCHA) method. This technique rigorously describes the full thermodynamics 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 ∗Authors to whom any correspondence should be addressed. Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1361-648X/21/363001+34$33.00 1 ©2021 The Author(s). Published by IOP Publishing Ltd Printed in the UK J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review 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 SSCHA and its dynamical extension, making particular emphasis on the physical inter pretation 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. Keywords: anharmonicity, stochastic self-consistent harmonic approximation, computational methods, ionic fluctuations, quantum effects, first-principles methods (Some figures may appear in colour only in the online journal) 1. 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 phonon–phonon interactions, as well as spectroscopic signatures detected in infrared, Raman, and inelastic x-ray or neutron 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 (BO) potential V(R), which, at an ionic configuration R,isgiven by the electronic ground state energy. In the standard harmonic approximation V(R) is Taylor-expanded up to secondorder around the R0ionic positions that minimize V(R). The resulting Hamiltonian is exactly diagonalizable in terms of phonons,the quantaof vibrations.Harmonicphononsare welldefined quasiparticles with an infinite lifetime, whose 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 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–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 displacivephase transition, such as a ferroelectric or charge-density wave (CDW) instability. In orderto calculate from first principles vibrational properties of solids beyond perturbation theory and overcome these difficulties, several methods have been developed in the last years [5–30]. Many of them are based on extracting renormalized phonon frequencies from ab initio molecular dynamics (AIMD) through velocity autocorrelation functions [6–9]or by extracting effective force constants (FCs) from the AIMD trajectory [10–12]. In order to include quantum effects on the ionicmotion, which are neglectedon 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–24,31], which are mainly inspired on the self-consistent harmonic approximation [32]orvibrational 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, manyof them simply correctthe free energy and/orthe phonon frequencies assuming that the ions remain fixed at the R0classical positions. However, as it has been shown recently in several compounds [34–37], the ionic positions can be strongly altered by quantum 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 renormalizedphononfrequenciesis, i.e., whetherthey are auxiliary phonon frequencies intrinsic to the devised theoretical framework or if they really represent the physical vibrational excitations probed experimentally. 2 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review The stochastic self-consistent harmonic approximation (SSCHA) [22–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 canoptimizecompletelythecrystalstructure,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 AIMD 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-perturbativethirdand fourth-orderphonon–phononscattering 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 section 2we 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 section 3. Then, we describe, in section 4, the post-processingfeatures of the code, which include calculations of the free energy Hessian for second-order 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 section 5, 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 section 6(SnTe), where we fully characterize the thermodynamics of the phase transition between the highsymmetry 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 section 7, we review some important results obtained so far with the SSCHA code. Finally, in section 8, we summarize the main conclusions. 2. The variational free energy The SSCHA is a theory that aims at describing the thermodynamics of a crystal, fully accountingfor quantum, thermal, and anharmonic effects of nuclei within the BO approximation. Thebasis ofallequilibriumthermodynamicsis 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 Tand the entropy S, which accounts for ‘disorder’ and is related to the number of microstates corresponding to the same macrostate of the system: F=E−TS.(1) 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 Ris 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 many-body 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[˜ρ]: F[˜ρ]=E[˜ρ]−TS[˜ρ]⩾F,(2) where E[˜ρ]=K+V(R)˜ρ(3) is the total energy (Kis the kinetic energy operator and V(R) the potential energy), and S[˜ρ] the entropy calculated with the trial density matrix. ·˜ρ=Tr[˜ρ·] indicates the quantum average of the operator ·taken with ˜ρ. If we pick any trial density matrix ˜ρ,F[˜ρ] is an upperbound 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 equation (2). Performing the optimization on any possible trial density matrix is, however, an unfeasible task due its many-body character that hinders an efficient parameterization. This is true also for a classical system: no exact parameterization 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) Rand the quantum-thermal fluctuations around them Φ(we have explicitly expressed the dependence of ˜ρon Rand Φby adding them as subindexes), just like any 3 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Figure 1. Illustration of the SSCHA method to a one dimensional particle problem at T=0 K. Panel (a) the one dimensional BO 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 ycoordinate of the points are the quantum total energy (including the zero-point motion), while the xaxis 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 Rand width are optimized to satisfy the least energy principle. The average position in the exact case is however obtained as Rρexact . Gaussian is defined by the average and mean square displacements. Within the SSCHA, we optimize Rand Φto minimize the free energy of the system. In this way, we compress the memory requested to store ˜ρR,Φ,asRdepends only on 3Nanumbers (the coordinates of the atoms), while the fluctuations Φare encoded in a symmetric square real matrix of 3Na×3Na.Nais the total number of atoms in the system. The free parameters in Rand Φ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 HR,Φ: HR,Φ=K+1 2 ab (Ra−Ra)Φab(Rb−Rb).(4) Here, Ris a real vector and Φa real matrix that parameterize the trial Hamiltonian, while Kand Rare 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 3Na. Let us note here, that, inspired by the harmonicshape of HR,Φ, we will also refer to Φas the auxiliary FCs. This mapping with a harmonic Hamiltonian is very useful, as both K˜ρR,ΦandS[˜ρR,Φ] become simply the kinetic energy and entropy of the auxiliary harmonic system HR,Φ,which are analytic functions of Φ. Hence, the only quantity that we really need to compute is the average over the interacting BO potential V(R)˜ρ=dRV(R)˜ρ(R).(5) The potential V(R) is the BO 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 Rand Φ, 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 Rand Φ. •Update Rand Φ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 1a 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 V(R)=3R4+1 2R3−3R2.(6) 4 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review We first study the classical harmonic result, obtained by Taylor-expanding the potential in equation (6) to second order around the minimum R0. 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 veryclose 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 lineshape, 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 energyis 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. The only quantity in the SSCHA encoding the role of electrons is the total energy landscape V(R), stochastically sampled on random configurations and computed with an external code. For this reason, within the SSCHA, we can describe electronically excited states by replacing the ground state energy in V(R) with an excited state obtained, for example, fixing the electronic occupation to force an electron–hole excitation. 3. 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 phase-diagrams (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. 3.1. 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 Rand the auxiliary FCs Φas F[R,Φ]=K+V(R)˜ρR,Φ−TSion ˜ρR,Φ.(7) Here, we explicit that the entropy Sion only accounts for ionic degrees of freedom (not electronic). After the SSCHA minimization, the final estimate of the equilibrium free energy is given by F=min R,ΦF[R,Φ]=F[Req,Φeq].(8) Therefore, the final result of a SSCHA free energy calculation is given in terms of the equilibrium configuration Req, the free energy F, and the SSCHA auxiliary FCs Φeq.The final free energy accounts for quantum and thermal ionic fluctuations without approximating the BO energy surface: it is valid to study thermodynamic properties. The Req positions determinethe average atomic positions also taken into account quantum/thermal fluctuations and anharmonicity. It is important to remark, however, that the Gaussian variance Φhas, in principle,norelationwith the experimentallyobservedphonon frequencies, as it is just a variable parameterizing the density matrix. The relation of it with the physical phononfrequencies is discussed in section 4.3. For clarity, we report in table 1a collection of symbols used in the equations, with their meanings and the reference on the equation where they have been introduced. 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,definedas G=F+P∗ΩVol ,(9) where P∗is the target pressure, ΩVol is the simulation box volume, and Fis the Helmholtz free energy. In this case, the code optimizes G⩽G[R,Φ]=F[R,Φ]+P∗ΩVol , (10) which can be used, for instance, to estimate the structural changes imposed by pressure by fully accounting for fluctuations. 5 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Table 1. Collection of some symbols frequently used in the main text. First column, the symbol used. Second column, a short description. Third column, equation or page-column of the first occurrence. Symbol Meaning First use RAtomic position (canonical variable) equation (3) V(R) Potential energy equation (3) RTrial centroid positions (parameter) equation (4) uDisplacement from the average atomic position Requation (15) ΦTrial harmonic matrix (parameter) equation (4) ΨStatic displacement–displacement correlation matrix relative to Φequation (15) HR,ΦTrial harmonic Hamiltonian for given R,Φequation (4) ˜ρR,ΦDensity matrix of HR,Φ(trial density matrix) Page 3-1 ˜ρR,Φ(R) Gaussian positional probability density for ˜ρR,ΦPage 3-1 F[R,Φ] SSCHA Helmholtz free energy functional equation (7) G[R,Φ] SSCHA Gibbs free energy functional equation (10) fHR,Φ(R) Forces for the Hamiltonian HR,Φacting on the ions when they are in the Rpositions equation (12) f(BO)(R) Born–Oppenheimer forces acting on the ions when they are in the Rpositions equation (12) P(BO)(R) Born–Oppenheimer stress tensor when the ions are in the Rpositions equation (19) ΦR2nd SSCHA force constants for a given R, it is the trial Φthat minimizes F[R,Φ] Page 12-2 DRΦRdivided by the square root of the masses equation (52) (3) ΦR,(4) ΦR3rd and 4th order SSCHA force constants for a given Requations (53)and(54) (3) DR,(4) DR (3) ΦR,(4) ΦRdivided by the square root of the masses equations (53)and(54) F(R) SSCHA positional Helmholtz free energy, given by F[R,ΦR] equation (50) Req SSCHA equilibrium centroids, trial centroids that minimizes F(R) equation (8) Φeq SSCHA harmonic matrix ΦReq , the trial Φat the minimum of the free energy functional equation (8) H(S) SSCHA effective harmonic Hamiltonian, given by HReq,Φeq equation (60) D(S) Dynamical matrix of H(S),givenbyDReq Page 14-1 (3) Deq,(4) Deq Symbols indicating (3) DReq and (4) DReq , respectively equation (62) D(F) Positional Helmholtz free energy Hessian divided by the square root of the masses Page 12-2 G(z) One-phonon Green function equation (67) Π(0),Π(z) Static and dynamic SSCHA self-energy equations (62)and(68) (B) Π(0), (B) Π(z) Static and dynamic SSCHA bubble self-energy equations (64)and(69) σ(q,Ω) Phonon spectral function (with the reciprocal lattice vector made explicit) equation (70) ωμ(q) Frequency of the (μ,q) SSCHA auxiliary phonon equation (16) Ωμ(q) Frequency of the (μ,q) static approximation phonon from D(F) Page 15-1 Ωμ(q),Γμ(q) Frequency and linewidth of the (μ,q) anharmonic phonon in the Lorentzian approximation equation (81) As made explicit in equation (7), only thermal effects on the ions are taken into account so far, whereas the electrons are considered at zero temperature. However, at very high temperatures the entropy associated to electrons may be important. Within the SSCHA, it is possible to explicitly include finite-temperature effects on the electrons too. The key is to replace in equation (7) the electronic ground state energy V(R) with the finite-temperature electronic free energy Fel(R)= Eel(R)−TSel(R) (if electrons have finite temperature, in the adiabatic approximation forces and equilibrium position of the ions are ruled by the electronic free energy). In this case the SSCHA method minimizes the functional F[R,Φ]=K+Fel(R)˜ρR,Φ−TS ion ˜ρR,Φ =K+Eel(R)˜ρR,Φ−TS ˜ρR,Φ, (11) where S˜ρR,Φ=Sel(R)˜ρR,Φ+Sion ˜ρR,Φ.Thesame trick can be applied to the Gibbs free energy minimization as well. Therefore, the SSCHA estimation of the system’s entropy can also incorporate contributions from both electrons (averaged through the ionic distribution ˜ρ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]. 3.2. The implementation of the free energy minimization In the SSCHA code, the minimization of F[R,Φ]isperformed through a preconditioned gradient descent approach, whichrequires thecalculation of thegradient ofthe free energy with respect to the centroid positions Rand the auxiliary FCs Φ. The partial derivatives are evaluated through the exact analytic formulas ∂F ∂Ra =−f(BO) a(R)−fHR,Φ a(R)˜ρR,Φ (12) 6 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review and ∂F ∂Φcd =1 2 ab ∂Ψab ∂Φcd f(BO) b(R)−fHR,Φ b(R) × e Ψ−1 ae (Re−Re)˜ρR,Φ .(13) Here f(BO)(R) are the BO forces that act on the ions when they are in the Rpositions; fHR,Φ(R) is the force given by the auxiliary harmonic Hamiltonian HR,Φ, fHR,Φ a(R)=− b Φab(Rb−Rb); (14) and Ψis the displacement–displacement correlation matrix Ψab =uaub˜ρR,Φ, (15) where with u=R−Rwe indicate the displacement from the average atomic position. Explicitly, Ψab =1 √MaMb μ (2nμ+1) 2ωμ ea μeb μ.(16) In equation (16), ωμand eμare the eigenvalues and eigenvectors of the mass rescaled auxiliary FCs Φab/√MaMb,andnμ 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 HR,Φ. In other words, they are only used to define the trial density matrix ˜ρR,Φ.Weshowhowto compute the physical anharmonic phonon frequencies of the system in section 4. It is convenient to give an explicit expression for the gradient of the free energy with respect to the auxiliary FCs in terms of the ωμeigenvalues and eμeigenvectors. As shown in reference [23] (see appendix B), the gradient can be rewritten as ∂F ∂Φcd = ab Λ[0]abcd √MaMbMcMdf(BO) b(R)−fHR,Φ b(R) × e Ψ−1 ae (Re−Re)˜ρR,Φ , (17) where Λ[0]abcd = μν  4ωνωμ ea νeb μec νed μ ×⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ dnμ dωμ−2nμ+1 2ωμ ,ων=ωμ nμ−nν ωμ−ωnu−1+nμ+nν ωμ+ωnu,ων=ωμ . (18) Here, nμ=1/(eβωμ−1). The reason why we have introduced the Λ[0] tensor will be evident in section 4.Evenif equation (17) looks different to the gradient introduced in the originalSSCHA workin reference[22], it can be demonstrated that both expressionsareequivalent by simply playing with the permutation symmetry of ∂2V ∂Ra∂Rb˜ρR,Φ . However, equation (18) unambiguously determines the value taken by the Λ[0] tensor for the ων=ωμcase, while the gradient in reference [22] did not describe explicitly what to do in this degenerate limit. At the end of the SSCHA optimization, apart from the temperature-dependentReq positionsand the equilibriumauxiliary FC matrix Φeq, the code also calculates the anharmonic stress tensor P, whichincludes both quantumand thermalionic fluctuations, as derivatives of the free energy with respect to a strain tensor ε: Pαβ =−1 ΩVol ∂F ∂εαβ ε=0 =P(BO) αβ (R)˜ρR,Φ −1 2ΩVol Na  s=1uα sf(BO)β s+uβ sf(BO)α s˜ρR,Φ .(19) Here, we have made explicit the atomic index s(lower index) and Cartesian α,β(upper index) of uand f(BO).P(BO)(R)is the BO stress tensor of the configuration with ions displaced in the Rcoordinates. This equation is slightly different from the stress tensor equation presented in [24]. The two equations coincide at equilibrium, but this is more general. The derivation of equation (19) is reported in appendix A. Thanks to the temperature-dependent stress, the SSCHA code can optimize also the lattice parameters and the volume. Thus, by relaxing the lattice at different temperatures, we get the thermal expansion straightforwardly. Remarkably, the stress tensor of equation (19) can be computed with a single SSCHA minimization at fixed volume. This is a huge advantage with respect to the standard quasiharmonicapproximation,notonlybecauseitincludesquantum and anharmonic effects, but also because it is computationally much more efficient. In fact, the quasiharmonicapproximation requires performing harmonic phonon calculations at different volumes (and/or internal lattice positions) to estimate the minimum of the quasi-harmonic free energy with finite differences. This process is extremely cumbersome for crystals with few symmetries and lots of internal degrees of freedom in the structure. In the current implementation of the code, the symmetries of the space group are imposed a posteriori on the gradients of equations (12)and(13), as well as on (19). This assures that the density matrix satisfies all the symmetries at each step of the minimization. Thus, during the geometry optimization, the system cannotlose any symmetry,thoughit can gain them.The symmetries are imposed followingthe methodologyexplained in appendix D, which is different to the method originallyconceived [22]. The current SSCHA code can also work without imposing symmetries, allowing for symmetry loss, though the stochastic number of configurations needed to converge the minimization is larger (see section 3.2.1). 3.2.1. The stochastic sampling. The stochastic nature of the SSCHA comes from the Monte Carlo evaluation of the 7 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review averages in equations (12), (13), and (19). A set of random ionic configurations are created in a chosen supercell according to the Gaussian ionic probability distribution ˜ρR,Φ(R)=det(Ψ−1/2π) ×exp −1 2 ab (Ra−Ra)Ψ−1ab(Rb−Rb). (20) The Monte Carlo average of a generic observable O(R), function only of the ionic position R, is calculated then as weighted sum over the created ensemble: O(R)˜ρR,Φ=1 Nc j=1ρj Nc  j=1 ρjO(Rj).(21) Here, Ncis the total number of configurations in the ensemble, while Rjis the jth ionic randomly displaced configuration. Each of the R{j}configurations is generated according to the initial trial ionic distribution ˜ρR(0),Φ(0) (R) from which the minimization starts. To improve the stochastic accuracy, for each Rjconfiguration also −Rjis created, benefiting from ˜ρR,Φ(R)=˜ρR,Φ(−R) property of the Gaussian distribution. The ρjweights are computed and updated along the free energy minimization as the values of Rand Φchange: ρj=˜ρR,ΦRj ˜ρR(0),Φ(0) Rj.(22) At the beginning, when the ensemble has just been generated and R=R(0) and Φ=Φ(0), all values of ρj=1. However, as the Rand Φare updated during the minimization, the weights change. This reweighting technique is commonly used in Monte Carlo methods [42,43] and takes the name of importance sampling. This allows avoiding generating a new ensemble and computing ab initio energies and forces at each step of the minimization, speeding up the SSCHA calculation. 3.2.2. Minimization algorithm. The minimization strategy implemented in the SSCHA code for the free energy is based on a preconditioned gradient descent. At each step, the Rand Φare updated as Φ(n+1) =Φ(n)−λΦ ab ∂2F ∂Φ∂Φab −1∂F ∂Φab (23) R(n+1) =R(n)−λR a∂2F ∂R∂Ra−1∂F ∂Ra .(24) In a perfectly quadratic landscape, this algorithm assures the convergence in just one step if both λΦand λRare 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 the Hessian matrix that multiplies the gradient, is the preconditioned steepest descent. If the preconditioning option is set to false, a standard steepest descent minimization is followed instead, with λΦand λRre-scaled to the maximum eigenvalue of the ∂2F ∂Φ2and ∂2F ∂R2Hessian matrices, respectively, in order to havea dimensionalvalues independenton the system. The preconditioning Hessian matrices that multiplies the gradients in equation (27)and(28) are approximated by the code.Followingthe procedureintroducedin reference[24],we use the exact Hessian in the minimum of a perfectly harmonic oscillator with the same frequencies as the SCHA auxiliary Hamiltonian. In particular, they are: ∂2F ∂Φab∂Φcd ≈1 2 ∂Ψab ∂Φcd (25) and ∂2F ∂R∂R≈Φ.(26) Equation (25) is presented differently from the original work in which it was derived [24]. We prove in appendix Bthat they are exactly the same. Considering that the Hessian preconditioner cancels out the 1 2 ∂Ψab ∂Φcd term in equation (13), the resulting update of the variational parameters at each step in the minimization is performed as Φ(n+1) ab =Φ (n) ab +−λΦf(BO) b(R)−fHR,Φ b(R) × c Ψ−1ac (Rc−Rc)˜ρR,Φ (27) and R(n+1) a=R(n) a+λR b Φ−1 ab f(BO) b(R)−fHR,Φ b(R)˜ρR,Φ . (28) This implementation is very efficient, especially for equation (27), as there is no need to calculate the ΛR[0] tensor. Therefore, computing directly equation (27)ismuchfasterthan calculating the gradient of equation (13). The code allows the user to select a different minimization algorithm specifically for the minimization with respect to Φ: the root representation. Since the minimization with respect to Φis the more challenging, this technique aims to further improving the Φoptimization. In particular, the gradient has a stochastic error and the minimization is performed with a finite step size. For these reasons, Φcould become non positive definite during the optimization (i.e. the dynamical matrix has imaginary frequencies). If this occurs, the minimization is halted raising an error, as the density matrix of equation (20) diverges. In such a case, the minimization must be manually restarted, either by taking a smaller step or by stopping the minimization before reaching imaginary frequencies (fixing the maximum number of steps). This kind of halts do not occur often when using the preconditioning in the minimization.However, they may be encountered if few configurations are generatedfor each ensembleor the starting dynamical matrix is very far from equilibrium. To solve these problems, we implement the root representation, in which, instead of updating Φas in equation (27), it 8 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review or σ(q,Ω)=−Ω πIm Tr Ω+i0+2𝟙+ −D(S)(q)+Π(q,Ω+i0+)−1, (71) where the multiplicative factor Ω/2πhas been included to have, foreach q, a functionthat integratedon the real axis gives the total numberof modes3na(nais the number of atoms in the unit cell, that may be different from the total number of atoms in the supercell Na). In the so-called ‘static approximation’, we replace the full self-energy Π(z) with the static self-energy Π(0), where zis blocked at zero. In this case the spectral function is (stat) σ(q,Ω)=−Ω πIm Tr Ω+i0+2𝟙−D(S)(q) +Π(q,0)) ]−1 =−Ω πIm TrΩ+i0+2𝟙−D(F)(q) −1, (72) whereinthelastlinewehaveusedequation(65). Therefore, (stat) σ(q,Ω)= μ (stat) σμ(q,Ω) (73) with (stat) σμ(q,Ω)=1 2δ(Ω−Ωμ(q)) +δ(Ω+Ω μ(q)), (74) 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 equations (64), (66), and (69)) (B) Πμν(q,Ω+iδse)=1 Nc k1k2 ρ1ρ2 G δG,q+k1+k2 ×F(Ω+iδse,ωρ1(k1), ωρ2(k2)) ×(3) Dμρ1ρ2(−q,−k1,−k2)(3) Dρ1ρ2ν ×(k1,k2,q), (75) where the summation k-grid can be arbitrarily fine as long as the interpolation of the third-order SSCHA FCs can be performed (as in the static case equation (66)), and δse is an arbitrary small, but finite, positive smearing value used to obtain converged results in the computation. In fact, the exact result corresponds to the limiting value obtained with an infinite kgrid and a zero δse smearing. In actual, finite-time calculations, the converged value of the dynamic self-energy is therefore estimated in this way. For a given summation k-grid, the corresponding self-energy converged value is estimated by analyzing the result given by equation (75) for smaller and smaller δse values (for the used k-grid, there will be a minimum value of δse under which the result shows numerical instability). This analysis is performed with finer and finer summation grids until the converged value in the thermodynamic limit is obtained. In principle, a dedicated convergence study of this kind has to be performed for all the specific observables of interest. With the dynamical SSCHA bubble self-energy, the code allows to compute the spectral function by the equation (see equation (71)) σ(q,Ω)=−Ω πImTr (Ω+iδid)2𝟙+ −D(S)(q)+ (B) Π(q,Ω+iδse)−1 , (76) with δid another arbitrary small, but finite, positive smearing value. The role of δid is significant when the imaginary part of the self-energy is small. A prominent example where this happens is when the spectral function is calculated in the static approximation, i.e. when the bubble self-energy is kept fixed at the static value (B) Πμν(q, 0) (see equation (72)). Indeed, in this case the self-energy is real (Hermitian) and the computed spectral function becomes (stat) σ(q,Ω)= μ 1 2Im 1 π 1 Ω−Ωμ(q)+iδid +1 π 1 Ω+Ω μ(q)+iδid , (77) 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, equation (74), but with smeared Dirac delta functions. Actually, this is not just an extreme example, since the code really gives the opportunity to compute the spectral function in the static approximation, replacing in equation (76) the full bubble self-energy with its static value computed through equation (66). This can be used to double-check that, as expected from equations (74)and(77), the obtainedspectral function is given by spikes aroundthe frequencies of the Hessian free energy matrix D(F) (computed in the bubble approximation). However, the role played by δid is not as critical as δse since it is not typically system-dependent and it does not require a convergence study: in the code its default value is automatically set depending on the spacing of the energy Ω-grid used to compute the spectral function. Given a q, the calculation of the full spectral function σ(q,Ω) through equation (76) turns out to be quite a heavy task due to the inversion of a different 3na×3namatrix for each Ωvalue. The code also allows to employ a much less computational demanding approach by discarding the offdiagonal elements of the computed dynamical self-energy in the SSCHA normal modes components(i.e. the componentsin 15 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review the D(S)(q)’s eigenvector basis). Within this ‘no mode-mixing’ approximation,which usually provesto be extremelygood, the SSCHA modes keep their individuality even after the renormalization due to anharmonic effects. Indeed in this case, as in the static approximation,equation (73), the total spectral function is given by the superposition of individual mode spectral functions: σ(q,Ω)= μ σμ(q,Ω), (78) where now the (q,μ)-mode spectral function σμ(q,Ω)is computed with σμ(q,Ω)=1 21 π−Im Zμ(q,Ω) [Ω−ReZμ(q,Ω)]2+[Im Zμ(q,Ω)]2 +1 π Im Zμ(q,Ω) [Ω+Re Zμ(q,Ω)]2+[Im Zμ(q,Ω)]2 (79) and Zμ(q,Ω)=ω2 μ(q)+Π μμ(q,Ω+iδse).(80) Therefore, computing the spectral function in the ‘no mode-mixing’ approximation, by measuring the deviation of σμ(q,Ω) from a Dirac delta function around Ωμ(q), it is possible to asses the impact that anharmonicity has on the different SSCHA modes (q,μ), separately. The form of the (q,μ)-mode spectral function σμ(q,Ω) in equation (79) resembles a Lorentzian, but with 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), σμ(q,Ω)=1 21 π Γμ(q) [Ω−Ωμ(q)]2+[Γμ(q)]2 +1 π Γμ(q) [Ω+Ω μ(q)]2+[Γμ(q)]2, (81) meaning that the quasiparticle picture is still valid, even after the inclusion of anharmonicity, with the (μ,q) quasiparticle having frequency (energy) Ωμ(q) and lifetime τμ(q)= 1/(2Γμ(q)). The difference between the renormalized and the ‘bare’ SSCHA frequency, Δμ(q)=Ω μ(q)−ωμ(q), is called the frequency shift of the (μ,q) mode. The SSCHA code offers several tools to perform such a ‘Lorentzian analysis’. In general, the best Lorentzian approximation is obtained with Ωμ(q)=Re Zμ(q,Ωμ(q)) (82) Γμ(q)=−Im Zμ(q,Ωμ(q)).(83) Once the dynamical self-energy and the Zμ(q,Ω)arecomputed, the SSCHA code allows to compute the single-mode spectral functions in the Lorentzian approximation, estimating the frequency Ωμ(q) and HWHMs Γμ(q) in different ways. One, optional, possibility is to solve self-consistently equation (82) to estimate Ωμ(q), and then Γμ(q), through equation (83). However, by default, the ‘one-shot’ approximation is employed with (os) Ωμ(q)=Re Zμ(q,ωμ(q)) (84) (os) Γμ(q)=−Im Zμ(q,ωμ(q)).(85) If the SSCHA self-energy Πis a (small) perturbation on the SSCHA free propagator (not meaning that we are in a perturbative regime with respect to the harmonic approximation), then perturbationtheory can be employed to evaluate the spectral function. If we keep the first order in the self-consistent equations equation (85), we get: (pert) Ωμ(q)=1 2ωμ(q)ReΠμμ(q,ωμ(q)) (86) (pert) Γμ(q)=−1 2ωμ(q)Im Πμμ(q,ωμ(q)+iδse).(87) This perturbative approach is also employed by the SSCHA code to evaluate the quasiparticles’ energies and lifetimes. Examples of spectral function calculations done with equations (76), (78)–(80)and(81) can be found in figure 4 of reference [48]. In figure 5of the same reference, the anharmonic phonon frequencies and linewidths along a path, computed using the Lorenztian approximation through equations (82)and(83), are shown. The spectral function computedwith equation (79) along a path is shown with a colorplot in figure 3 of reference [49] for PbTe, and in figure 4 of reference [50] for SnSe. In conclusion, with the SSCHA code we can calculate three frequencies for a mode (q,μ): ωμ(q), Ωμ(q), and Ωμ(q), which are the frequencyof the SSCHA auxiliaryboson, 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 otherapproachesused to estimate anharmonicphonons,where no distinction is usually done. 5. 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 16 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review loads the starting dynamical matrices; sets up the ensemble and the parameters for the SSCHA run; performs the calculations of the BO energies, forces, and stress tensors on the configurations in the ensemble by calling to a external totalenergy-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-steptutorials 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 librariesto 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. The code is hosted on a GitHub repository at the date of publication, located at https://github.com/SSCHAcode.News about future releases, tutorials and documentation, are hosted in the website www.sscha.eu. 5.1. Code structure Most of the program is written in Python with an objectoriented 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 postprocessing purposes). Methods of this class allow the user to impose symmetries on the system, constrain the auxiliary force to be positivedefinite (equation (32)), extractauxiliary phonon frequencies and polarization vectors, or interpolate them to other points in the BZ. All the calculations related to the SSCHA averages are performed by the ensemble class (inside Python-sscha). This class generates and stores all the randomly displaced ionic configurations, and can submit or load the results of the energy, forces, and stress tensors calculations. It also computes the quantities related to averages on the ensemble, as the free energy, the gradients, the stress tensor, and the free energy Hessian. Finally there are other classes, which employ the ensemble and perform the minimization of the free energy, take care of communicating with a remote cluster to run the calculation of forces and energies (see next section), and manage the post-processing to compute the spectral function (the full description of them is provided within the documentation of the code). Most of the code is written in Python, however, the heaviest CPU-intensive calculation is written in Fortran and interfaced with Python through the f2py utility provided by numpy [55]. In particular, the calculation of the free energy gradient, the free energy Hessian, the spectral functions, the interpolation, and the symmetrization are performed by a Fortran module compiled with the code. For this reason, in order to compile and use the code, a Fortran compiler as well as LAPACK and BLAS libraries are required. 5.2. 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 BO energies, forces, and stress tensors of the generated ionic configurations in each population (the red shaded cell in the code flowchart in figure 2). Each of these calculations is independent from the others, so they can be trivially run in parallel on different computing nodes. This is a huge advantage with respect to other methods based on AIMD or PIMD, which mimic a time evolution of the system and thus require to calculate atomic forces on one configuration after the other. The SSCHA code does not include a particular engine for computing energies, forces and stresses, but relies on external software. For this reason, it is possible to exploit the efficient parallelization already implemented by the chosen software. For example, the widely used Quantum ESPRESSO package recentlyimplementedalso a hybrid parallelizationthatexploits together multi-threading (OpenMP), multiprocessing (MPI), and GPU (CUDA) parallelization [56]. In this way the SSCHA code stands on the shoulders of giants, exploiting the most efficient parallelization available today. All other steps of the code are generally computationally verycheap comparedto the energyand forcecalculation, especially when an ab initio approach is followed. The SSCHA minimization cannot exploit so well the possibilities offered by parallelization, since each step of the main cycle depends on the previous one. Most of the computations executed in the cycle are linear algebra calculations carried out with the numpy library [55], some of them speeded up with an explicit Fortran implementation. Thanks to the numpy implementation [57], if this library is correctly compiled, the linear algebra calculations will exploit multi-threading. For this reason, the best performances of the SSCHA are obtained by executing the ab initio calculations on an HPC facility, while the SSCHA minimization on a commercial workstation in which the minimization can take few seconds. Post-processing calculations, like the free energy Hessian and the phonon dynamicalspectral functions, may be executed with additional Python scripts after the end of the SSCHA run. The calculation of the free energy Hessian has been parallelized with OpenMP (multi-threading), while the calculation of the spectral functions, which may require a dense k-point grid for the interpolation, exploits multiprocessor parallelization through MPI (both mpi4py and pypar can be used [58–60]). 5.3. 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 17 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review feed these structures to a total-energy engine, e.g. a DFT code, to calculate their BO energies, atomic forces, and stress tensors. The user should preparelater 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 throughthe atomic simulation environment(ASE) library [61], which already implements interfaces with most common ab initio codes like Quantum ESPRESSO [62,63], VASP [64], SIESTA [65], CP2K [66], and many more. Forcefield codes like LAMMPS [67] may also be used. In this execution mode, the code will proceed automatically to perform the calculations locally, and the full flowchart in figure 2is 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 multithread 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. 5.4. Distribution The package is distributed as a standard Python application, and can be installed with a setup.py script. Since parts of the code are written in Fortran and C, it requires the appropriate compilers with LAPACK and BLAS libraries to be installed. Part of the Fortran subroutines are modified versions of Quantum ESPRESSO subroutines from the PHonon package, especially those regardingthe 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. 6. 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 reference [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 Fm¯ 3mphase to the low-temperature R3mphase. 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 energyof the toy modelas a functionof the atomicdisplacements projected onto the eigenvectors of the imaginary mode (the order parameter Δ) is reported in figure 3(a), where it is clear that the high-symmetry Fm¯ 3mis not at the minimum of V(R). However,as extensivelydiscussed above,the stability of a structure is determined by the temperature-dependentfree energy, F=E−TS, (88) and not the BO potential. Notably, Eis 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 eventhe behavior at T=0, where there is no entropycontribution. Since entropy usually is higher in high-symmetric positions (Δ=0), the Fm¯ 3mhigh symmetry phase will become progressively more stable as temperature increases. In figure 4we show the evolution of the free energy, its gradient, and the frequencies of the auxiliary FCs 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 Fm¯ 3m, 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). Figure 4(d) reports the stochastic condition to stop the minimization (and extract a new ensemble), as defined in equation (39) (we chose η=0.4). Here, we need only three ensembles 18 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Figure 3. (a) Born–Oppenheimer energy as a function of the order parameter of the ferroelectric phase transition of SnTe obtained with a model force field. (b) Hysteresis cycle between the ferroelectric phase R3mand the cubic paraelectric phase Fm¯ 3m. In both heating and cooling we constrained the SSCHA simulation to the R3msymmetry, which is a subgroup of Fm¯ 3m. (c) Free energy of the two phases. The high-symmetry phase becomes more stable around 200 K, slightly before the R3mfalls into the high-symmetry phase in the heating cycle. The dashed line indicates a dynamical instability. (d) The free energy curvature around the order parameter in the high-symmetry Fm¯ 3m phase. Positive values mean (meta)stability; negative values indicate a dynamical instability. The transition occurs at T=154 K, and it coincides with the lower bound for the Fm¯ 3min the hysteresis cycle. 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 FC 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 previouslower temperatureis used as input. At low temperatures the output structure remains the R3m, with Δ=0, but at T=205 K, the low-symmetry phase jumps into the highsymmetryphase,markinga first-orderphasetransition.We can confirm it is a first-order phase transition as we can start cooling down fromthe high-symmetryphase (withoutconstraining the new symmetries acquired) and the system remains stable up to T=160 K, when it transformsback 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 lowsymmetry phase is more stable up to 200 K, so that the phase diagram in this model is formed by the R3mphase below 200 K and the Fm¯ 3mabove. We can also see whether the Fm¯ 3mbecomes dynamically unstable by calculating the 19 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Figure 4. Convergence of the SSCHA minimization in the SnTe system in a 2 ×2×2 supercell at T=250 K. (a) Evolution of the free energy per unit cell during the minimization. The width of the line denotes the stochastic error. (b) Evolution of the modulus of the gradient of the auxiliary force constants. (c) Evolution of the frequencies of the auxiliary force constants. (d) The Kong–Liu effective sample size ratio (we used 0.4 as stochastic criterion for restarting) during the minimization. Vertical dotted lines indicate the new population after the simulation was out of the stochastic criterion. Vertical solid lines indicate a new population after the simulation converged by increasing the number of configurations to improve the accuracy. The simulation is performed starting with 50 configurations and increasing after successful convergence to 100 and 200. The overall number of total energy and forces calculations required to converge this example is 450. One last step is shown for demonstrative purposes, where we increased the number of configurations to 1000 to show how well the result converges already with few configurations (in particular the auxiliary frequencies of panel (c)). 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 Fm¯ 3mphase 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 R3mphase, 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’ (equation (59)) is very accurate, the correct estimation of the free energyHessian in the R3mphaserequires the full expression of the Hessian (equation (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 Fm¯ 3mbelow the critical point. We can do simulations directly in the highsymmetry 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 evenmore: finite temperature structure search. To investigate whether the R3mis 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 R3mphase 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 discovereda new phase of Cc symmetrydefinedina1×2×1 supercell of the original cubic cell. In figure 5we report the final phase diagram for the SnTe toy model. The new Cc phase 20 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review 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). is found to be the ground state up to 250 K, where the cubic Fm¯ 3mbecomes again energetically favorable. The Cc continues to exists until 280 K, where it transforms into the Fm¯ 3m 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 Fm¯ 3mphase is of second-order type (the order parameter does not jump, and the free energy lines of the two phases touch when the Fm¯ 3mmode becomes imaginary). This system has been already studied with the SSCHA in reference [49] with ab initio energies and forces. However, even if it is just 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 saddlepoint structures in the classical BO 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 complexsystems with many atoms in the unit cell, such as in molecular crystals, where symmetry constrains may be inefficient [68]. 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 fluctuations 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 postprocessing utility as discussed in section 4.3. 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 6do 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 calculation of the spectral functions. The effective charges are considered following the procedure outlined in appendix E3. In figure 7we illustrate the convergence for a phonon linewidth 2Γμ(q) with respect to the δse parameter and the k-meshsummationgridinequation(75). All the data of this simulation has been obtained in less than 1 h, 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. 7. Applications of the SSCHA method In orderto illustrate some physicalproblems 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. 21 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Figure 6. Phonon spectrum of SnTe at T=280 K. Left panel: spectral function at Γ. The two main peaks are the LO–TO splitting. Right panel: the full spectral function along a path in the Brillouin zone. The red dashed line is the dispersion of the auxiliary phonons (eigenvalues of D(S)). Figure 7. Convergence study of the linewidth (full width at half maximum) of the SnTe highest optical phonon frequency in Γ.For increasing size of the used k-mesh summation grid, the linewidth as a function of the smearing parameter δse is studied (see section 4.3 for details about these quantities). The result shows that the converged value of the linewidth (15 cm−1) is obtained with a 30 ×30 ×30 k-mesh summation grid, at least, and smearing δse around 0.8 cm−1. 7.1. 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 anharmonicrenormalization,impossible to capturewithin perturbative approaches [21,22,34–36,48,53,69,70]. The anharmonic renormalization of phonons in these compounds has been crucial to explain the superconducting properties of many hydrogen-based superconductors. For instance, the anomalousinverseisotopeeffect onpalladiumhydrides,which makes the deuterium compound acquire a larger superconducting critical temperature Tcthan the protium compound [71,72], is a consequence of a huge anharmonic renormalization of the phonons [21]. Also, the experimentally found hightemperature superconductivity in H3S around 200 K [73]and in LaH10 around 250 K [74,75] at high pressures can only be explained if phonon frequencies renormalized by anharmonicity are considered in the superconductivity equations [34,36]. In figure 8we show the huge anharmonic renormalization of the phonon frequencies for H3S[34,48]. Superconductivityin hydrogen compounds can be both largely suppressed but also enhanced by anharmonicity depending on the system [76]. The quantum effects and anharmonicity that the SSCHA captures go beyond the renormalization of phonon frequencies. For crystals with Wyckoff positions not fixed by symmetry, quantum or thermal fluctuations may strongly modify the atomic positions, resulting in a structure with atoms far from the positions that minimize the V(R) potential, occupying, instead, those that minimize the quantum F(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 H3SandLaH 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 experimentalRaman and infrared spectra [35,37]. The change in the crystal structure that the SSCHA captures goes beyond 22 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Figure 8. (a) Anharmonic phonon spectra obtained with the SSCHA for H3S at 158 GPa in the Im¯ 3mphase (top panel). 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 Im¯ 3mH3S. Data taken from reference [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, Fm¯ 3m. the internal degrees of freedom and can largely impact also the cell parameters.In figure8we illustrate the apparentdifference between the structure found classically from the minimum of V(R) and the one obtained from the quantum energy landscape for LaH10. 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 Tc’s approaching room temperature. The SSCHA method will be of great importance in the quest of new high-Tccompounds 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. 7.2. 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 Tcdw [77–79]. As the temperature dependence of phonon frequencies is a purely anharmonic property, the SSCHA has been used to predict from first principles Tcdw in several transition metal dichalcogenides(TMDs) both in the bulk and the monolayer [41,51,54,80,81]. The standard procedure in these calculations is to apply the SSCHA for the high-symmetryphase at different temperatures and calculate the spectra associated to the free energy Hessian D(F). These phonons represent the static limit of the physical phonons observed experimentally, which can be accessed with the SSCHA by calculating instead the spectral function as described in section 4. At the temperature at which D(F) develops a null eigenvalue,the high-symmetrystructureis 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 phononcollapse occurs. In figure9we show as an example the temperature dependenceof the phonon spectra derived from the free energy Hessian in monolayer NbSe2and the consequenttheoretical determinationof 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 Tcdw, and setting (4) Deq =0inequation(62) seems in general a good approximation.However, convergingTcdw is rather sensitive to the SSCHA supercell and rather large supercells may be needed to converge the CDW transition temperatures [41,54]. The capacity of the SSCHA of predicting Tcdw purely ab initio without empirical fitting parameters offers a fantastic tool to determining the physics behind CDW transitions. The forcecalculationsneededfortheSSCHA variational 23 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Figure 9. (a) Phonon spectra of monolayer NbSe2derived from the free energy Hessian as a function of temperature. (b) Squared phonon frequency of the lowest energy mode at q=1/3ΓMas a function of temperature and the determination of the CDW temperature. Data taken from reference [41]. 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 monolayerTMDs to strain [51] and doping [54], the difference (or similarities) between the CDW transitionsin bulkandthe correspondingtwo-dimensionalstructures [51,54], as well as the importance of Van der Waals forces in the melting of CDW transitions [81]. Consequently the SSCHA program is expected to have a large impact on theoretical studies of CDW transitions for many type of materials, not just TMDs. 7.3. Phase transitions, spectral functions, and thermal conductivity in semiconducting materials Phase transitions related to soft phononsare also very common in ferroelectric, thermoelectric, and other functional materials. The SSCHA is again a perfect method to study these phase transitions considering that many of the high-temperature phases of these compounds are not a minimum of the BO potential V(R), but saddle points. Thus, it becomes imperative to adopt a non-perturbative treatment of anharmonicity in order to study their thermodynamic and transport properties such as the thermal conductivity. Whether these phase transitions are purely second-order or first order it is not always evident experimentally, unless a clear softening to zero of a phonon mode at the transition temperature is observed. The SSCHA can distinguish between continuous and discontinuous transitions as discussed in the practical example provided in section 6. For instance, in order to convincingly show that the transition between the high-temperature Cmcm phase of SnSe and the low-temperature Pnma is second order, at the temperatureat which the free energyHessian developeda negative eigenvalue a SSCHA relaxation was performed starting from the low-temperature phase. It was shown that the Pnma phase relaxed at this temperature into the Cmcm, showing that the Pnma phaseisnolongera minimumofthefreeenergy[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 emergenceof extra peaks can be misinterpreted with phase transitions. TheSSCHA is a perfect methodfor capturing these subtleties as it provides the spectral function σ(q,Ω) without the Lorentzian approximation. It has been used to understand thecomplexσ(q,Ω) in PbTe, SnTe, SnSe, and SnS [49,50,52]. In figure 10 we show the spectral function calculated within the SSCHA for Cmcm SnSe at 800 K, where the σμ(q,Ω) contribution of some particular modes is clearly anomalous and deviates from the standard Lorentzian picture. With the phonon frequencies and the phonon linewidths obtained with the SSCHA, transport properties such as the thermal conductivity can be calculated with an external code, for instance, within Boltzmann transport equations [82,83]. It has been shown that employing the SSCHA phonon scattering tensor (3) ΦRin the thermal transport calculations leads to a very good agreement with experimental results, in contrast with what was obtained by employing the standard third-order derivatives of the BO total energy. The difference is that the former includes higher-order anharmonic terms coming from the average over the thermal ensemble (equation(53)). Both in 24 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review Kα1α2α3 s1s2s3(0,l2,l3|p)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Φα1α2α3 s1s2s3(0,l2,l3)p  s3,l3Φα1α2α3 s1s2s3(0,l2,l3) pif p=0orp>0and s3,l3Φα1α2α3 s1s2s3(0,l2,l3)=0 0ifp>0and s3,l3Φα1α2α3 s1s2s3(0,l2,l3)=0 , (E19) with pa non-negative real number which can be arbitrarily fixed to optimize the calculation performances (in the equation above the convention 00=1 has been adopted). The ' Φα1α2α3 s1s2s3(0,l2,l3) defined through equation (E18) fulfills the ASR on the third index, since for any p⩾0 the scaling factor fulfills the normalization condition  s3,l3Kα1α2α3 s1s2s3(0,l2,l3|p)=1.(E20) The value of phas effects on the way the different terms of Φα1α2α3 s1s2s3(0,l2,l3) are scaled. For p=0 the scaling factor is a pure geometric quantity related to the three atoms clusters. Indeed, given s1,s2,l2, the scaling factor Kα1α2α3 s1s2s3(0,l2,l3|p= 0) is fully determined(it is the same for all the αh,s3,l3) and,in particular, it doesnotdependonthe FCs value.On the contrary, for p=0, given αh,s1,s2,l2we have Kα1α2α3 s1s2s 3(0,l2,l 3|p) Kα1α2α3 s1s2s 3(0,l2,l 3|p)= Φα1α2α3 s1s2s 3(0,l2,l 3) Φα1α2α3 s1s2s 3(0,l2,l 3) p (E21) so that if p>1 the scaling factor is higher(lower) forFCs have lower (higher) absolute value, otherwise the opposite. E.3. 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 FCs 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 [91,92] Φ(dd) st (q)=( Φ(dd) st (q)−δst t( Φ(dd) st(q=0) (E22) with ( Φ(dd) st (q)=4π ΩVol  G (G+q)·Z∗ s⊗(G+q)·Z∗ t (G+q)·∞·(G+q) ×e−(G+q)·∞·(G+q) 4η2ei(G+q)·(τs−τt), (E23) where we have explicitly indicated only the atomic indices (i.e. we are using component-free notation for the Cartesian indices), ηis a parameter whose value has to be large enough to allow to includeonlythe reciprocalspace termsinthe Ewald sum, and  Gis the sum over reciprocal lattice vectors such thatG+q=0(thesum includesas manyG’s as it is necessary to reach the convergence for the considered η)[62]. Once Z∗ sand ∞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 qpoint of the BZ, the corresponding dipole–dipole terms Φ(dd)(q) are subtracted and the resulting short range FCs is Fourier transformed to the real space. Subsequently, this real space short-range FCs, Φ(sr)(l), can be Fourier transformed back to any k∈BZ and the corresponding long-range dipole–dipole analytical contribution Φ(dd)(k)is added [62]: Φ(q)onBZq-grid Subtract dipole-dipole interaction terms Φ(dd)(q) + Fourier transform to real space −−−−−−−−−−−−−−−−−−−−→Φ(sr)(l) Fourier transform back to k∈BZ + Add dipole-dipole interaction term Φ(dd)(k) −−−−−−−−−−−−−−−−−−→Φ(k) (E24) The dipole–dipole correction to the FCs given by equations (E22)and(E23) is nonanalytic at zone center and its q→0 limit depends on the direction (q=q/qalong which the limit is performed: lim δ→0+ Φ(dd) st (δ(q)=Φ(dd) st (0)+Φ(dd−na) st ((q), (E25) 31 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review where Φ(dd−na) st ((q)=4π ΩVol (q·Z∗ s⊗(q·Z∗ t (q·∞·(q(E26) is the nonanalyticzone-center correctionterm. When a phonon dispersion through Γis calculated, the SSCHA code includes the nonanalytic correction term in the zone center, with the direction given by the followed path [62]. When the SSCHA code calculate the spectral properties (static or dynamic), it adds the nonanalytic correction term in the zone center dynamical matrix (necessary for the integral over the BZ) from a random direction. ORCID iDs LorenzoMonacelli https://orcid.org/0000-0002-6381-3741 Ion Errea https://orcid.org/0000-0002-5719-6580 References [1] Maradudin A A and Fein A E 1962 Scattering of neutrons by an anharmonic crystal Phys. Rev. 128 2589–608 [2] Cowley R A 1968 Anharmonic crystals Rep. Prog. Phys. 31 123 [3] Calandra M, Lazzeri M and Mauri F 2007 Anharmonic and nonadiabatic effects in MgB2: implications for the isotope effect and interpretation of Raman spectra Physica C456 38–44 [4] Errea I 2016 Approaching the strongly anharmonic limit with ab initio calculations of materials’ vibrational properties—a colloquium Eur. Phys. J. B89 237 [5] Car R and Parrinello M 1985 Unified approach for molecular dynamics and density-functional theory Phys.Rev.Lett.55 2471–4 [6] Wang C Z, Chan C T and Ho K M 1990 Tight-binding molecular-dynamics study of phonon anharmonic effects in silicon and diamond Phys. Rev. B42 11276–83 [7] Ljungberg M P and Íñiguez J 2013 Temperature-dependent classical phonons from efficient nondynamical simulations Phys. Rev. Lett. 110 105503 [8] Magd˘ au I B and Ackland G J 2013 Identification of highpressure phases III and IV in hydrogen: simulating Raman spectra using molecular dynamics Phys. Rev. B87 174110 [9] Zhang D-B, Sun T and Wentzcovitch R M 2014 Phonon quasiparticles and anharmonic free energy in complex systems Phys. Rev. Lett. 112 058501 [10] Hellman O, Abrikosov I A and Simak S I 2011 Lattice dynamics of anharmonic solids from first principles Phys. Rev. B84 180301 [11] Hellman O, Steneteg P, Abrikosov I A and Simak S I 2013 Temperature dependent effective potential method for accurate free energy calculations of solids Phys. Rev. B87 104111 [12] Hellman O and Abrikosov I A 2013 Temperature-dependent effective third-order interatomic force constants from first principles Phys. Rev. B88 144301 [13] Ceperley D M 1995 Path integrals in the theory of condensed helium Rev. Mod. Phys. 67 279–355 [14] Zhong W, Vanderbilt D and Rabe K M 1995 First-principles theory of ferroelectric phase transitions for perovskites: the case of BaTiO3Phys. Rev. B52 6301–12 [15] Wojdea J C, Hermet P, Ljungberg M P, Ghosez P and Aiguez J 2013 First-principles model potentials for lattice-dynamical studies: general methodology and example of application to ferroic perovskite oxides J. Phys.: Condens. Matter 25 305401 [16] Errea I, Rousseau B and Bergara A 2011 Anharmonic stabilization of the high-pressure simple cubic phase of calcium Phys. Rev. Lett. 106 165501 [17] Zhou F, Nielson W, Xia Y and Ozoli V 2014 Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations Phys.Rev.Lett.113 185501 [18] Tadano T and Tsuneyuki S 2015 Self-consistent phonon calculations of lattice dynamical properties in cubic SrTiO3with first-principles anharmonic force constants Phys. Rev. B92 054301 [19] Tadano T, Gohda Y and Tsuneyuki S 2014 Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations J. Phys.: Condens. Matter 26 225402 [20] Azadi S, Monserrat B, Foulkes W M C and Needs R J 2014 Dissociation of high-pressure solid molecular hydrogen: a quantum Monte Carlo and anharmonic vibrational study Phys.Rev.Lett.112 165501 [21] Errea I, Calandra M and Mauri F 2013 First-principles theory of anharmonicity and the inverse isotope effect in superconducting palladium-hydride compounds Phys.Rev.Lett.111 177002 [22] Errea I, Calandra M and Mauri F 2014 Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: application to platinum and palladium hydrides Phys. Rev. B89 064302 [23] Bianco R, Errea I, Paulatto L, Calandra M and Mauri F 2017 Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: theory and stochastic implementation Phys. Rev. B96 014111 [24] Monacelli L, Errea I, Calandra M and Mauri F 2018 Pressure and stress tensor of complex anharmonic crystals within the stochastic self-consistent harmonic approximation Phys. Rev. B98 024106 [25] Souvatzis P, Eriksson O, Katsnelson M I and Rudin S P 2008 Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory Phys. Rev. Lett. 100 095901 [26] Tang X, Li C W and Fultz B 2010 Anharmonicity-induced phonon broadening in aluminium at high temperatures Phys. Rev. B82 184301 [27] Antolin N, Restrepo O D and Windl W 2012 Fast free-energy calculations for unstable high-temperature phases Phys. Rev. B86 054119 [28] Parlinski K 2018 Ab initio determination of anharmonic phonon peaks Phys. Rev. B98 054305 [29] Eriksson F, Fransson E and Erhart P 2019 The Hiphive Package for the extraction of high-order force constants by machine learning Adv. Theory Simul. 21800184 [30] van Roekeghem A, Carrete J and Mingo N 2020 Quantum selfconsistent ab initio lattice dynamics [31] Monserrat B, Drummond N D and Needs R J 2013 Anharmonic vibrational properties in periodic systems: energy, electron–phonon coupling, and stress Phys. Rev. B87 144302 [32] Hooton D J 1955 LI. A new treatment of anharmonicity in lattice thermodynamics: I London, Edinburgh Dublin Phil. Mag. J. Sci. 46 422–32 [33] Bowman J M 1978 Self-consistent field energies and wavefunctions for coupled oscillators J. Chem. Phys. 68 608 [34] Errea I et al 2016 Quantum hydrogen-bond symmetrization in the superconducting hydrogen sulfide system Nature 532 81–4 [35] Borinaga M, Riego P, Leonardo A, Calandra M, Mauri F, Bergara A and Errea I 2016 Anharmonic enhancement of superconductivity in metallic molecular Cmca—4 hydrogen at high pressure: a first-principles study J. Phys.: Condens. Matter 28 494001 32 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review [36] Errea I et al 2020 Quantum crystal structure in the 250 K superconducting lanthanum hydride Nature 578 66–9 [37] Monacelli L, Errea I, Calandra M and Mauri F 2020 Black metal hydrogen above 360 GPa driven by proton quantum fluctuations Nat. Phys. 17 63 [38] Monacelli L and Mauri F 2020 Time-dependent self-consistent harmonic approximation: Anharmonic nuclear quantum dynamics and time correlation functions Phys. Rev. B103 104305 [39] Lihm J-M and Park C-H 2020 Gaussian time-dependent variational principle for the finite-temperature anharmonic lattice dynamics (29 October 2020) (arXiv:2010.15725) [40] Huang G Q, Chen L F, Liu M and Xing D Y 2004 Electronic structure and electron–phonon interaction in the ternary silicides malsi (m=ca, sr, and ba) Phys. Rev. B69 064509 [41] Bianco R, Monacelli L, Calandra M, Mauri F and Errea I 2020 Weak dimensionality dependence and dominant role of ionic fluctuations in the charge-density-wave transition of NbSe2 Phys. Rev. Lett. 125 106101 [42] Neal R M 2001 Stat. Comput. 11 125–39 [43] Miotto M and Monacelli L 2018 Entropy evaluation sheds light on ecosystem complexity Phys. Rev. E98 042402 [44] Kong A, Liu J S and Wong W H 1994 Sequential imputations and Bayesian missing data problems J. Am. Stat. Assoc. 89 278–88 [45] Liu R H et al 2009 A large iron isotope effect in SmFeAsO1−xfx and Ba1−xkxFe2As2Nature 459 64–7 [46] Shulumba N, Hellman O and Minnich A J 2017 Intrinsic localized mode and low thermal conductivity of PbSe Phys. Rev. B95 014302 [47] Shulumba N, Hellman O and Minnich A J 2017 Lattice thermal conductivity of polyethylene molecular crystals from firstprinciples including nuclear quantum effects Phys.Rev.Lett. 119 185901 [48] Bianco R, Errea I, Calandra M and Mauri F 2018 Highpressure phase diagram of hydrogen and deuterium sulfides from first principles: structural and vibrational properties including quantum and anharmonic effects Phys. Rev. B97 214101 [49] Ribeiro G A S, Paulatto L, Bianco R, Errea I, Mauri F and Calandra M 2018 Strong anharmonicity in the phonon spectra of PbTe and SnTe from first principles Phys. Rev. B97 014306 [50] Aseginolaza U, Bianco R, Monacelli L, Paulatto L, Calandra M, Mauri F, Bergara A and Errea I 2019 Phonon collapse and second-order phase transition in thermoelectric SnSe Phys. Rev. Lett. 122 075901 [51] Bianco R, Errea I, Monacelli L, Calandra M and Mauri F 2019 Quantum enhancement of charge density wave in NbS2in the two-dimensional limit Nano Lett. 19 3098–103 PMID: 30932501 [52] Aseginolaza U, Bianco R, Monacelli L, Paulatto L, Calandra M, Mauri F, Bergara A and Errea I 2019 Strong anharmonicity and high thermoelectric efficiency in high-temperature SNS from first principles Phys. Rev. B100 214307 [53] Errea I et al 2015 High-pressure hydrogen sulfide from first principles: a strongly anharmonic phonon-mediated superconductor Phys. Rev. Lett. 114 157004 [54] Zhou J S, Monacelli L, Bianco R, Errea I, Mauri F and Calandra M 2020 Anharmonicity and doping melt the charge density wave in single-layer TiSe2Nano Lett. 20 4809–15 PMID: 32496779 [55] Harris C R et al 2020 Array programming with numpy Nature 585 357–62 [56] Giannozzi P et al 2020 Quantum ESPRESSO toward the exascale J. Chem. Phys. 152 154105 [57] van der Walt S, Colbert S C and Varoquaux G 2011 The numpy array: a structure for efficient numerical computation Comput. Sci. Eng. 13 22–30 [58] Dalcín L, Paz R and Storti M 2005 MPI for Python J. Parallel Distrib. Comput. 65 1108–15 [59] Dalcín L, Paz R, Storti M and D’Elía J 2008 MPI for Python: performance improvements and MPI-2 extensions J. Parallel Distrib. Comput. 68 655–62 [60] Dalcin L D, Paz R R, Kler P A and Cosimo A 2011 Parallel distributed computing using Python Adv. Water Resour. 34 1124–39 [61] Hjorth Larsen A et al 2017 The atomic simulation environment—a Python library for working with atoms J. Phys.: Condens. Matter 29 273002 [62] Giannozzi P et al 2009 QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials J. Phys.: Condens. Matter 21 395502 [63] Giannozzi P et al 2017 Advanced capabilities for materials modelling with quantum ESPRESSO J. Phys.: Condens. Matter 29 465901 [64] Kresse G and Furthmüller J 1996 Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set Phys. Rev. B54 11169–86 [65] Soler J M, Artacho E, Gale J D, García A, Junquera J, Ordej´ on P and Sánchez-Portal D 2002 The SIESTA method for ab initio order-nmaterials simulation J. Phys.: Condens. Matter 14 2745–79 [66] Kühne T D et al 2020 CP2k: an electronic structure and molecular dynamics software package—quickstep: efficient and accurate electronic structure calculations J. Chem. Phys. 152 194103 [67] Plimpton S 1995 Fast parallel algorithms for short-range molecular dynamics J. Comput. Phys. 117 1–19 [68] Monserrat B, Drummond N D, Dalladay-Simpson P, Howie R T, Ríos P L, Gregoryanz E, Pickard C J and Needs R J 2018 Structure and metallicity of phase V of hydrogen Phys. Rev. Lett. 120 255701 [69] Borinaga M, Errea I, Calandra M, Mauri F and Bergara A 2016 Anharmonic effects in atomic hydrogen: superconductivity and lattice dynamical stability Phys. Rev. B93 174308 [70] Biswas S, Errea I, Calandra M, Mauri F and Scandolo S 2019 Ab initio study of the LiH phase diagram at extreme pressures and temperatures Phys. Rev. B99 024108 [71] Stritzker B and Buckel W 1972 Superconductivity in the palladium–hydrogen and the palladium–deuterium systems Z. Phys. 257 1–8 [72] Schirber J E and Northrup C J M 1974 Concentration dependence of the superconducting transition temperature in pdhx and pddxPhys. Rev. B10 3818–20 [73] Drozdov A P, Eremets M I, Troyan I A, Ksenofontov V and Shylin S I 2015 Conventional superconductivity at 203 K at high pressures in the sulfur hydride system Nature 525 73–6 [74] Somayazulu M, Ahart M, Mishra A K, Geballe Z M, Baldini M, Meng Y, Struzhkin V V and Hemley R J 2019 Evidence for superconductivity above 260 K in lanthanum superhydride at megabar pressures Phys.Rev.Lett.122 027001 [75] Drozdov A P et al 2019 Superconductivity at 250 K in lanthanum hydride under high pressures Nature 569 528–31 [76] Pickard C J, Errea I and Eremets M I 2020 Superconducting hydrides under pressure Annu. Rev. Condens. Matter Phys. 11 57–76 [77] Weber F, Rosenkranz S, Castellan J-P, Osborn R, Karapetrov G, Hott R, Heid R, Bohnen K-P and Alatas A 2011 Electron–phonon coupling and the soft phonon mode in TiSe2 Phys.Rev.Lett.107 266401 [78] Weber F et al 2011 Extended phonon collapse and the origin of the charge-density wave in 2H-NbSe2Phys.Rev.Lett.107 107403 33 J. Phys.: Condens. Matter 33 (2021) 363001 Topical Review [79] Leroux M et al 2015 Strong anharmonicity induces quantum melting of charge density wave in 2H-NbSe2under pressure Phys. Rev. B92 140303 [80] Sky Zhou J, Bianco R, Monacelli L, Errea I, Mauri F and Calandra M 2020 Theory of the thickness dependence of the charge density wave transition in 1 t-TiTe22D Mater. 7 045032 [81] Diego J et al 2021 Van der Waals driven anharmonic melting of the 3d charge density wave in VSe2Nat. Commun. 12 598 [82] Fugallo G, Lazzeri M, Paulatto L and Mauri F 2013 Ab initio variational approach for evaluating lattice thermal conductivity Phys. Rev. B88 045430 [83] Li W, Carrete J, Katcho N A and Mingo N 2014 ShengBTE: a solver of the Boltzmann transport equation for phonons Comput. Phys. Commun. 185 1747–58 [84] Onuorah I J, Bonf` a P, Renzi R D, Monacelli L, Mauri F, Calandra M and Errea I 2019 Quantum effects in muon spin spectroscopy within the stochastic self-consistent harmonic approximation Phys. Rev. Mater. 3073804 [85] Aseginolaza U, Cea T, Bianco R, Monacelli L, Calandra M, Bergara A, Mauri F and Errea I 2020 Bending rigidity and sound propagation in graphene [86] Maradudin A A and Vosko S H 1968 Symmetry properties of the normal vibrations of a crystal Rev. Mod. Phys. 40 1–37 [87] Warren J L 1968 Further considerations on the symmetry properties of the normal vibrations of a crystal Rev. Mod. Phys. 40 38–76 [88] Hendrikse Z W, Elout M O and Maaskant W J A 1995 Computation of the independent elements of the dynamical matrix Comput. Phys. Commun. 86 297–311 [89] Togo A and Tanaka I 2018 Spglib: a software library for crystal symmetry search (5 August 2018) (arXiv:1808.01590) [90] Paulatto L, Mauri F and Lazzeri M 2013 Anharmonic properties from a generalized third-order ab initio approach: theory and applications to graphite and graphene Phys. Rev. B87 214303 [91] Gonze X and Lee C 1997 Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory Phys. Rev. B55 10355–68 [92] Giannozzi P, de Gironcoli S, Pavone P and Baroni S 1991 Ab initio calculation of phonon dispersions in semiconductors Phys. Rev. B43 7231–42 34