## 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 |

- Also see Quantum Espresso Walkthrough: QE Units

## 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

- Wu, Zhongqing, and Renata M. Wentzcovitch. “Quasiharmonic thermal elasticity of crystals: an analytical approach.”
*Physical Review B*83.18 (2011): 184115. - 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.