## Correction on Equations

In the paper Wu (2011) [@wu2011quasiharmonic] , such a term is identified in Equation (15.2) as

\begin{aligned} c_{iijj}^{th} = \frac{k_BT}{V} \sum_{qm} [ &-Q^2_{qm}\frac{e^{Q_{qm}}}{(e^{Q_{qm}} - 1)^2} \gamma^{ii}_{qm}\gamma^{jj}_{qm} \ &+ \frac{e^{Q_{qm}}}{(e^{Q_{qm}} - 1)}(\gamma^{ii}_{qm}\gamma^{jj}_{qm}-\frac{\partial \gamma^{ii}_{qm}}{\partial e_{jj}} + \gamma^{ii}_{qm}\delta_{ii})] \end{aligned}\tag{15.2}

It should be

\begin{aligned} c_{iijj}^{th} = \frac{k_BT}{V} \sum_{qm} [ &-Q^2_{qm}\frac{Q_{qm}}{(e^{Q_{qm}} - 1)^2} \gamma^{ii}_{qm}\gamma^{jj}_{qm} \ &+ \frac{Q_{qm}}{(e^{Q_{qm}} - 1)}(\gamma^{ii}_{qm}\gamma^{jj}_{qm}-\frac{\partial \gamma^{ii}_{qm}}{\partial e_{jj}} + \gamma^{ii}_{qm}\delta_{ii})] \end{aligned} \tag{15.2}

## Notes on Input

### Input for input01

#### Explaination of variables

Variable Meaning Usage notes
nv number of PVEs
nq number of $q$-vectors
np number of normal modes
nm number of formula units
na number of atoms

#### Extra notes

• 3 Lines before reading nv, nq, np, nm, na.
• No empty line before and after weight.
• The frequencies in thermo.f is $f$, in QHA is $\omega$ (see unit conversion section). QE’s output is $\omega$.

### Input for elast.dat

#### Explaination of variables

Variable Meaning Usage notes
V_0 calibration volume used when fitting $P(V)$, added by my self
NforE number of volumes used when reading static $c_{ij}$
cellmass mass of unit cell used when calculating sound velocities

### Parameters Specified in the Source Code thermo2.1.0.f

#### Explaination of variables

Variable Meaning Usage notes
mxdt dimension maximum of index for temperature nt < mxdt
mxdv dimension maximum of index for volume ntv < mxdv
mxdq dimension maximum of index for $q$-points nq < mxdq
mxdf dimension maximum of index for mode frequencies np < mxdf
• Remember also to change mxdq1, mxdf1, mxdv1 in subroutine readenergy.
• Make sure the demension for VforE, PforE, freq and cc0 is correct at the beginning of the code
• Make sure the dimension for freq is correct at the beginning of readenergy.
• Make sure the dimension for freq is correct at the beginning of modegamma.
• Make sure the dimension for VforE, PforE is correct at the beginning of PbyIP.

## Tidbits in Calculation

### EoS fitting

Equation of state (EoS) is calculated by fitting a Helmholtz free energy ($F$)-strain ($x_j$) relationship. The strain is calculated from $V$ as follow:

$$x_j = \frac{1}{2}\left[\left(\frac{V_\text{ref}}{V_j}\right)^\frac{2}{3}-1\right]$$

Or convert back:

$$V_j = V_{\text{ref}},(1 + 2x_j)^{-\frac{2}{3}}$$

Then fit free energy $F_j$ with polynomial function $F(x_j)$ (see function BMF).

Pressure is calculated with

$$P(T,V_j) = -\left.\frac{\partial F(T,V)}{\partial V}\right|_{V=V_j}$$

### modegamma

You need to keep an eye on the $\omega(V)$, $\gamma(V)$ fitting, it is very error-prone. Crossing in $\omega(V)$ is disasterous

### Unit conversion constants

Variable Expression Value Example
bol or b 0.00008617 $0.00008617$ Boltzmann’s constant $k$ = $8.617 \times 10 ^ {-5}$ eV / K
a0 0.529177d0 $0.529177$ Bohr radius $a_0$ 1 Bohr = 0.529177 Å
ry 13.6058d0 $13.6058$ Rydberg Constant1 Rydberg = 13.6058 eV
ev 1.602d-12 $1.602 \times 10^{-12}$ (???)
ev_A3 160.218 $160.218$ 1 eV / Å$^3$ = 160.218 GPa
pcv (a0**3.0)/ry/1.602d2 $0.000067985$ 1 GPa = $6.79855 \times 10^{-5}$ Ry / bohr$^3$
h 0.000123985 0.000123985 $hc$ = 0.000123985 eV⋅cm

## Notes on Compiling

This program is most compatible with ifort from Intel® C++ Compiler. Its compatibility with gfortran from GCC, the GNU Compiler Collection, is not warrented.

One might also use a input file with a lot of $q$-points, and a lot of volumes. In this case, you should really work on supercomputers with very large RAMs—there are large arrays like freq(mxdfil,mxdv,mxdq,mxdf).

The best friend you have while debugging this program, is WRITE(*,*). However, there is tools like gdb and can work with VSCode that is available, but you might not find it as helpful. Because of this, you might end up needing a machine with large hard drive also.

## Output

### Regular output (elast_x00.dat)

• $P$
• $V$
• $c_{11}^S(T,P)$, $c_{22}^S(T,P)$, $c_{33}^S(T,P)$
• $c_{12}^S(T,P)$, $c_{13}^S(T,P)$, $c_{23}^S(T,P)$
• $c_{44}^S(T,P)$, $c_{55}^S(T,P)$, $c_{66}^S(T,P)$
• $\beta_\text{V}$, $\beta_\text{R}$, $\beta_\text{VRH}$: bulk modulus
• $G_\text{V}$, $G_\text{R}$, $G_\text{VRH}$: shear modulus
• $v_p$, $v_s$: sound velocities
• $\beta^T(T,P)$, $\beta^S(T,P)$

## References

1. Wu, Zhongqing, and Renata M. Wentzcovitch. “Quasiharmonic thermal elasticity of crystals: an analytical approach.” Physical Review B 83.18 (2011): 184115.
2. Zou, Fan, et al. “An extended semianalytical approach for thermoelasticity of monoclinic crystals: application to diopside.” Journal of Geophysical Research: Solid Earth 123.9 (2018): 7629-7643.