where ê are phonon eigenvectors, and the temperature dependence enters through phonon mode population n λ given by Bose-Einstein distribution. Once a temperaturedependent ensemble of configurations is sampled, we in effect use IFCs and not DFT to predict their forces. Specifically, by collecting the sampled configurations into a new displacement matrix A(T ) matching the clusters present in the IFCs, we use the anharmonic IFCs to expediently predict anharmonic forces in those configurations via F = A(T )Φ. We then obtain effective temperaturedependent harmonic IFCs, Φ TD 2 , that describe the computed forces. These amount to calculating Φ