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

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

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} $$


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.


Regular output (elast_x00.dat)


  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.