For large-scale constitutive strength models the shear modulus is typically assumed to be linearly dependent on temperature. However, for materials compressed beyond the Hugoniot or in regimes where there is very little experimental data, accurate and validated models must be used. To this end, we present here a new methodology that fully accounts for electron- and ion-thermal contributions to the elastic moduli over broad ranges of temperature (<20,000 K) and pressure (<10 Mbar). In this approach, the full potential linear muffin-tin orbital (FP-LMTO) method for the cold and electron-thermal contributions is closely coupled with ion-thermal contributions. For the latter two separate approaches are used. In one approach, the quasi-harmonic, ion-thermal contribution is obtained through a Brillouin zone sum of strain derivatives of the phonons, and in the other a full anharmonic ion-thermal contribution is obtained directly through Monte Carlo (MC) canonical distribution averages of strain derivatives on the multi-ion potential itself.