$\newcommand{\AA}{\text{Å}}$

CRYSTALpytools.thermodynamics module

A post-processing module for DFT lattice dynamics by harmonic and quasiharmonic approximations (HA/QHA).

class Harmonic(temperature=[], pressure=[], filename=None, autocalc=True, use_old=False)

Bases: object

A class for harmonic phonon calclulations. It can be parameterized from a CRYSTAL output file, phonopy ouput file or by setting all the information (usually for QHA).

Parameters:

temperature (array|float) – Optional Temperatures where thermodynamic properties are computed. Unit: K
pressure (array|float) – Optional Pressures where thermodyanmic properties are calculated. Unit: GPa
filename (str | None) – Name of the printed-out file. If None, do not print out file.
autocalc (bool) – Automatically launch calculations.
use_old (bool) – Use the deprecated text fomat output, which cannot be used for restarting calculation.

Temperatures and pressures can also be defined by self.thermodynamics, whose entries always cover the entries here.

Phonon dispersions are forced to be summed if the automatic scheme (autocalc=True) is launched. To get verbose outputs, call self.thermodynamics() first and then call self.print_results().

Usage:

ha = Harmonic(temperature=[0, 100, 200, 300], pressure=[0.,])
ha.from_file('harmonic_phonon.out')

property edft: For compatibility

classmethod restart(filename)

Restart calculation from a dumped YAML file. autocalc is switched off. Instantiation will not update the dumped file, but thermodynamic calculations will.

Parameters:: filename (str) – YAML file
Returns:: obj (Harmonic)

from_file(filename, source='crystal', scale=1.0, scale_range=[], imaginary_tol=-0.0001, q_overlap_tol=0.0001, q_id=None, q_coord=None, read_eigvt=False, **kwargs)

Generate the Harominc object from an output file. Files generated by the following methods are supported (also the accepted entries for source):

'crystal': CRYSTAL harmonic phonon outputs (Gamma, dispersion).
'crystal-QHA': CRYSTAL QHA outputs.
'phonopy': Phonopy qpoints and mesh files.

Note

In rare cases, the user might want to collect data of a few specific q points. q_id and q_coord can be set but not simultaneously. If set, q_id takes priority and q_coord is ignored.

Parameters:

filename (str) – Name of the output file.
source (str) – The program used for the output file.
scale (float|array) – Scaling factors
scale_range (array) – nScale*2 array. The frequency range of the scaling factor. Including the minimum and excluding the maximum. Use empty list for a uniform scaling.
imaginary_tol (float) – The threshold of negative frequencies.
q_overlap_tol (float) – The threshold of overlapping points, defined as the 2nd norm of the difference of fractional q vectors
q_id (array) – Index (from 0) of q points to be read. 1*nqpoint.
q_coord (array) – Fractional coordinates of q points to read. nqpoint*3 list.
read_eigvt – Deprecated.
**kwargs – Other keywords passed to corresponding I/O functions. See below.
qha_index (int) – crystal-QHA. The index of calculation to read (from 0).
struc_yaml (str) – phonopy. ‘qpoints.yaml’ only.
u_0 (float) – phonopy. Static internal energy in kJ/mol.
irreps_yaml (array[str]) – phonopy. Read ‘irreps.yaml’ for mode symmetry symbols. The fractional coordinates in the file are used for q_coord if not specified. Phonopy ‘irreps.yaml’ only contains symmetry info of a q point. Use a list of filenames for multiple q points.

Returns:

self (Harmonic) – Attributes listed below.
structure (CStructure) – Extended Pymatgen structure class. For the finite displacement method, this is the cell before supercell expansion.
natom (int) – Number of atoms in the cell.
volume (float) – Volume of the cell. Unit: $\AA^{3}$
u_0 (float) – Internal energy. Unit: kJ/mol
nqpoint (int) – Number of q points
qpoint (list) – nQpoint*4 list. The first three elements are fractional coordinates in reciprocal space and the last is its weight, i.e., number of equivalent q points.
nmode (int) – Number of modes.
frequency (array[float]) – nQpoint*nMode array of frequencies in THz.
mode_symm (array[str]) – nQpoint*nMode array of the irreducible representations. In Mulliken symbols.
eigenvector (array[float]) – nQpoint*nMode*nAtom*3 array of phonon eigenvectors.

from_frequency(u_0, qpoint, frequency, eigenvector, symmetry, structure=None, natom=None, volume=None)

Generate a Harmonic object by specifying frequency and eigenvector. Usually called by other methods.

Parameters:

u_0 (float) – Static total energy in kJ/mol.
qpoint (array[float]) – nQpoint*4 array. The first three elements are fractional coordinate and the last is weight.
frequency (array[float]) – Array of frequencies. Unit: THz
eigenvector (array[float]) – Mass weighted and phased eigenvectors
1. (normalized to)
symmetry (array[str]) – Irreducible representations in Mulliken symbols.
structure (CStructure) – Extended Pymatgen structure class. For the finite displacement method, this is the cell before expansion.
natom (int) – Number of atoms in the reduced cell.
volume (float) – Volume of the reduced cell. Unit: $\AA^{3}$.

Note

The user should define structure. natom and volume are for developers only.

Returns:

self (Harmonic) – Attributes see from_file().

Raises:

Exception – If computational data is stored in the object.
ValueError – If neither of the 2 available options are defined.

thermodynamics(**kwargs)

Calculate the thermodynamic properties (zp_energy, u_vib, entropy, c_v and Gibbs and Helmholtz free energy) of the HA system.

Zero-point energy $E^{zp}$:

\[E^{zp}=\sum_{i,\mathbf{q}}\frac{1}{2}\hbar\omega_{i,\mathbf{q}}\]

Vibration contribution to internal energy $U^{vib}$:

\[U^{vib}\left(T\right)=E^{zp}+\sum_{i,\mathbf{q}} \frac{\hbar\omega_{i,\mathbf{q}}}{\exp{\left( \frac{\hbar\omega_{i,\mathbf{q}}}{k_{B}T} \right)}-1}\]

Entropy $S$:

\[S\left(T\right)=k_{B}\sum_{i,\mathbf{q}}\left\{ \frac{\hbar\omega_{i,\mathbf{q}}}{k_{B}T\left[ \exp{\left( \frac{\hbar\omega_{i,\mathbf{q}}}{k_{B}T} \right)}-1 \right]}-\ln{\left[ 1-\exp{\left( -\frac{\hbar\omega_{i,\mathbf{q}}}{k_{B}T} \right)} \right]} \right\}\]

Constant volume specific heat $C^{V}$:

\[C_{V}\left(T\right)=\sum_{i,\mathbf{q}} \frac{\left(\hbar\omega_{i,\mathbf{q}}\right)^{2}}{k_{B}T^{2}} \frac{\exp{ \left( \frac{\hbar\omega_{i,\mathbf{q}}}{k_{B}T} \right)} }{\left[ \exp{\left( \frac{\hbar\omega_{i,\mathbf{q}}}{k_{B}T} \right)-1} \right]^{2} }\]

The Helmholtz and Gibbs free energies are defined as:

\[ \begin{align}\begin{aligned}F(T) = U_{0} + F_{vib}(T) = U_{0} + U_{vib}(T) - TS(T)\\G(T, p) = F + pV\end{aligned}\end{align} \]

Parameters:

temperature (array|float) – Optional Unit: K
pressure (array|float) – Optional Unit: GPa
sumphonon (bool) – Deprecated. Not used.

Returns:

self (Harmonic) – Attributes listed below.
self.helmholtz (array) – 1*nTemperature. Unit: KJ/mol
self.gibbs (array) – nPressure*nTemperature. Unit: KJ/mol
self.zp_energy (float) – Zero-point energy. Unit: KJ/mol
self.u_vib (array) – Vibrational contribution to internal energy. 1*nTemperature. Unit: KJ/mol
self.entropy (array) – 1*nTemperature. Unit: J/mol*K
self.c_v (array) – Constant volume specific heat. 1*nTemperature. Unit: J/mol*K

Raises:

Exception – If temperature and pressure are defined neither here nor during initialization

mode_thermo(q, m, **kwargs)

Return to thermodynamic properties of the specified mode(s). Units are consistent with the thermodynamics() method.

Parameters:

q (int|array) – Indices of q point, from 0.
m (int|array) – Indices of mode, from 0.
**kwargs – See below.
temperature (array|float) – If not present, use self.temperature. It will not update the attribute.
pressure (array|float) – If not present, use self.pressure. It will not update the attribute.

Returns:

zp_energy (float)
u_vib (array) – 1*nTemperture
entropy (array) – 1*nTemperture
c_v (array) – 1*nTemperture
helmholtz (array) – 1*nTemperture, excluding DFT total energy of the system
gibbs (array) – nPressure*nTemperture, excluding DFT total energy of the system

write_HA_result(): Write output. Typically the code dumps data to file automatically.

from_phonopy(phono_yaml, struc_yaml=None, edft=None, scale=1.0, imaginary_tol=0.0001, q_overlap_tol=0.0001, q_id=None, q_coord=None): Deprecated. Call from_file() directly.

class Quasi_harmonic(temperature=[], pressure=[], filename=None, use_old=False)

Bases: object

Generate and rearrange harmonic phonons, store the fitted, volume-dependent QHA phonon information and obtain the QHA thermodynamic properties.

Parameters:

temperature (array|float) – Should be given in ascending order, or will be sorted. Unit: K
pressure (array|float) – Should be given in ascending order, or will be sorted. Unit: GPa
filename (str) – Name of the output file in YAML format for restarting calculation.
use_old (bool) – Use the deprecated text fomat output, which cannot be used for restarting calculation.

Temperatures and pressures can also be defined by self.thermodynamics, whose entries always cover the entries here.

Usage:

qha = Quasi_harmonic()
qha.from_QHA_file('qha_phonon.out')
qha.thermo_freq(eos_method='birch_murnaghan', temperature=[0, 100, 200, 300], pressure=[0., 0.5, 1.]):

from_file(*filename, source='crystal', mode_sort_tol=0.4, **kwargs)

Instatiation from output file(s). Files generated by the following methods are supported (also the accepted entries for source):

'crystal': CRYSTAL harmonic phonon outputs (Gamma, dispersion).
'crystal-QHA': CRYSTAL QHA outputs.
'phonopy': Phonopy qpoints and mesh files.

Instatiation from outputs.

Parameters:

*filename (str) – Phonon output filename(s). Extendable.
source (str) – The program used for the output file.
mode_sort_tol (float | None) – The threshold of close mode overlaps with eigenvectors if present. If none, do not sort modes and eigenvector is not read.
**kwargs – Optional keywords passed to Harmonic().from_file().
scale (float|array)
scale_range (array)
imaginary_tol (float)
q_overlap_tol (float)
q_id (array) – Applied to all HA calculations.
q_coord (array) – Applied to all HA calculations.
struc_yaml (list[str]) – phonopy. ‘qpoints.yaml’ only.
u_0 (list[float]) – phonopy. Static internal energy in kJ/mol.
irreps_yaml (array[str]) – phonopy. nCalc*nQpoint array. Read ‘irreps.yaml’ for mode symmetry symbols. The fractional coordinates in the file are used for q_coord if not specified. Phonopy ‘irreps.yaml’ only contains symmetry info of a q point. Use a list of filenames for multiple q points.

Returns:

self (Quasi_harmonic) – New Attributes listed below
ncalc (int) – Number of HA phonon calculations.
combined_volume (array[float]) – 1*nCalc volumes. Unit: $\AA^{3}$
combined_struc (list[CStructure]) – 1*nCalc HA structures.
combined_u0 (array[float]) – 1*nCalc static internal energies. Unit: kJ/mol
combined_freq (array[float]) – nQpoint*nMode*nCalc frequencies. Unit: THz
combined_symm (array[str]) – nQpoint*nMode*nCalc mode symmetry in Mulliken symbol.
nqpoint (int) – Number of q points.
qpoint (array[float]) – nQpoint*4 array. The first three elements are fractional coordinates, and the last is the weight.
nmode (int) – Number of modes.
natom (int) – Number of atoms.

property combined_edft: For compatibility with older versions

classmethod restart(filename)

Restart calculation from a dumped YAML file. Instantiation will not update the dumped file, but thermodynamic calculations will.

Parameters:: filename (str) – YAML file
Returns:: obj (Quasi_harmonic)

thermo_freq(eos_method='birch_murnaghan', poly_order=[2, 3], min_method='BFGS', volume_bound=None, mutewarning=False, **kwargs)

Obtain thermodynamic properties by explicitly fitting phonon frequencies as polynomial functions of volume. DFT total energies are fitted as a function of volume by equation of states (EOS).

The equilibrium volume is fitted by minimizing Gibbs free energy at constant temperature and pressure.

\[V(T,p)=\text{min}[G(V;T,p)]=\text{min}[E_{0}(V)+F_{vib}(V;T,p)+pV)]\]

Parameterized and tested algorithms for min_method:

BFGS(no boundary)
L-BFGS-B(with boundary)

Isothermal bulk modulus $K_{T}$ is obtained by a secondary fit of $F(V; T)$ with the same eos_method.

\[K_{T}(p;T) = V(p;T)\left(\frac{\partial^{2}F(V;T)}{\partial V^{2}}\right)_{T}\]

Note

For a good fitting of isothermal bulk modulus $K_{T}$, the number of pressures must >= 5.

The mode-specific Grüneisen parameters are obtained by deriving the the fitted frequency at given temperature, pressure and equilibrium volume:

\[\gamma_{\textbf{q}i}(T,p)=-\frac{V_{0}}{\omega_{\textbf{q}i}}\frac{\text{d}\omega_{\textbf{q}i}}{\text{d}V_{0}}\]

Then the macroscopic Grüneisen parameters and thermal expansions are obtained with the standard Grüneisen method. See the thermo_gruneisen() method below.

Parameters:

eos_method (str) – EOS used to fit DFT total energy and Helmholtz free energy (to get bulk modules). Valid entries are consistent with PyMatGen eos module.
poly_order (array[int]) – The order of polynomials used to fit frequency as the function of volumes.
min_method (string) – Minimisation algorithms, ‘BFGS’/’L-BFGS-B’. See above.
volume_bound (tuple-like) – volumes. Unit: $\AA^{3}$
mutewarning (bool) – Whether print out warning messages.
**kwargs – See below
temperature (array[float]) – Should be given in ascending order, or will be sorted. Unit: K
pressure (array[float]) – Should be given in ascending order, or will be sorted. Unit: GPa
order (int) – For DeltaFactor / Polynomial EOSs.
min_ndata_factor (int) – For Numerical EOS.
max_poly_order_factor (int) – For Numerical EOS.
min_poly_order_factor (int) – For Numerical EOS.

Returns:

self (Quasi_harmonic) – New Attributes listed below
self.temperature (array) – Unit: K
self.pressure (array) – Unit: GPa
self.volume (array) – nPressure*nTemperature, same below. Equilibrium volumes. Unit: $\AA^{3}$
self.helmholtz (array) – Helmholtz free energy. Unit: kJ/mol
self.gibbs (array) – Gibbs free energy. Unit: kJ/mol
self.entropy (array) – Entropy. Unit: $J.mol^{-1}.K^{-1}$
self.gruneisen (array) – npressure*ntemperature, same below. Macroscopic Grüneisen parameter.
self.alpha_v (array) – Thermal expansion coefficient by Grüneisen method.
self.c_v (array) – Constant volume specific heat. Unit: $J.mol^{-1}.K^{-1}$
self.c_p (array) – Constant pressure specific heat by Grüneisen method. Unit: $J.mol^{-1}.K^{-1}$
self.k_t (array) – Isothermal bulk modulus. Unit: GPa.
self.k_s (array) – Adiabatic bulk modulus by Grüneisen method. Unit: GPa.
self.eos_method (str) – Name of the EOS.
self.u0_fit (Pymatgen EOS) – EOS used to fit DFT energy.
self.eos (list[Pymatgen EOS]) – Pymatgen EOS objects. EOS used to fit Helmholtz free energy.
self.e0_eos_method – Deprecated. Synonym of ‘self.eos_method’.
self.e0_eos – Deprecated. Synonym of ‘self.u0_fit’.

Raises:

Exception – If temperature or pressure is defined neither here nor during initialization.

thermo_gruneisen(eos_method='birch_murnaghan', min_method='BFGS', volume_bound=None, mutewarning=False, **kwargs)

Obtain thermodynamic properties by Grüneisen parameters.

The pressure and temperature independent mode-specific Grüneisen parameters are obtained by:

\[\gamma_{\textbf{q}i}=-\frac{\ln{\omega}-\ln{\omega_{0}}}{\ln{V}-\ln{V_{0}}}\]

$V_{0}$ is the most compact structure sampled. Frequency at given volume can be expressed analytically by:

\[\omega_{\textbf{q}i}(V) = \omega_{\textbf{q}i}(V_{0})\left(\frac{V}{V_{0}}\right)^{-\gamma_{\textbf{q}i}}\]

Then the Gibbs free energy is minimized in the same way as thermo_freq().

The macroscopic Grüneisen parameter is defined as:

\[\gamma=\sum_{\textbf{q}i}\frac{\gamma_{\textbf{q}i}C_{V,\textbf{q}i}}{C_{V}}\]

Thermal expansion coefficient in Grüneisen model:

\[\alpha_{V}^{gru}=\frac{\gamma C_{V}}{K_{T}V}\]

Therefore this method does not require polynomial fittings of volume-temperature for expansion rate, adiabatic bulk modulus and constant pressure specific heat.

Parameters:

eos_method (str) –
EOS used to fit DFT total energy and Helmholtz free energy (to get bulk modules). Valid entries are consistent with PyMatGen eos module.
min_method (string) – Minimisation algorithms, ‘BFGS’/’L-BFGS-B’. See thermo_freq().
volume_bound (tuple-like) – volumes. Unit: $\AA^{3}$
mutewarning (bool) – Whether print out warning messages.
**kwargs – See below
temperature (array[float]) – Should be given in ascending order, or will be sorted. Unit: K
pressure (array[float]) – Should be given in ascending order, or will be sorted. Unit: GPa
order (int) – For DeltaFactor / Polynomial EOSs.
min_ndata_factor (int) – For Numerical EOS.
max_poly_order_factor (int) – For Numerical EOS.
min_poly_order_factor (int) – For Numerical EOS.

Returns:

self (Quasi_harmonic) – New Attributes listed below
self.temperature (array) – Unit: K
self.pressure (array) – Unit: GPa
self.volume (array) – nPressure*nTemperature, same below. Equilibrium volumes. Unit: $\AA^{3}$
self.helmholtz (array) – Helmholtz free energy. Unit: kJ/mol
self.gibbs (array) – Gibbs free energy. Unit: kJ/mol
self.entropy (array) – Entropy. Unit: $J.mol^{-1}.K^{-1}$
self.gruneisen (array) – npressure*ntemperature, same below. Macroscopic Grüneisen parameter.
self.alpha_v (array) – Thermal expansion coefficient by Grüneisen method.
self.c_v (array) – Constant volume specific heat. Unit: $J.mol^{-1}.K^{-1}$
self.c_p (array) – Constant pressure specific heat by Grüneisen method. Unit: $J.mol^{-1}.K^{-1}$
self.k_t (array) – Isothermal bulk modulus. Unit: GPa.
self.k_s (array) – Adiabatic bulk modulus by Grüneisen method. Unit: GPa.
self.eos_method (str) – Name of the EOS.
self.u0_fit (Pymatgen EOS) – EOS used to fit DFT energy.
self.eos (list[Pymatgen EOS]) – Pymatgen EOS objects. EOS used to fit Helmholtz free energy.
self.alpha_vgru – Deprecated. Synonym of ‘self.alpha_v’.
self.c_pgru – Deprecated. Synonym of ‘self.c_p’.
self.k_sgru – Deprecated. Synonym of ‘self.k_s’.

thermo_eos(eos_method='birch_murnaghan', poly_order=[3, 4], min_method='BFGS', volume_bound=None, mutewarning=False, **kwargs)

Obtain thermodynamic properties by fitting EOS, which is fitted by the Helmholtz free energies of sampled harmonic phonons. The explicit sorting and fitting of frequency-volume relationship is disabled.

The equilibrium volume is fitted by minimizing the difference between analytical pressure and the given pressure at constant temperature.

\[V(T,p_{0})=\text{min}\left[ -\left( \frac{\partial F(V)}{\partial V} \right)_{T} - p_{0} \right]^{2}\]

Parameterized and tested algorithms for min_method:

BFGS(no boundary)
L-BFGS-B(with boundary)

Entropy is obtained by taking the derivative of Gibbs free energy at constant pressure.

\[S=-\left(\frac{\partial G}{\partial T}\right)_{p}\]

Constant pressure specific heat is obtained by taking the second derivative of $G$.

\[C_{p}=-T\left(\frac{\partial^{2}G}{\partial T^{2}}\right)_{p}\]

$G(T)$ is fitted as polynomial specified by poly_order. The first order is forced to be 0 to ensure $S=0$ at 0 K.

Note

For a good fitting of $G(T)$, entropy and c_p, the number of temperatures must >= 5. poly_order must > 2.

Parameters:

eos_method (str) –
EOS used to fit DFT total energy and Helmholtz free energy (to get bulk modules). Valid entries are consistent with PyMatGen eos module.
poly_order (array[int]) – The order of polynomials used to fit Gibbs free energy as the function of volumes.
min_method (string, optional) – Minimisation algorithms, ‘BFGS’/’L-BFGS-B’. See above.
volume_bound (tuple-like) – volumes. Unit: $\AA^{3}$
mutewarning (bool) – Whether print out warning messages.
**kwargs – See below
temperature (array[float]) – Unit: K
pressure (array[float]) – Unit: GPa
order (int) – For DeltaFactor / Polynomial EOSs.
min_ndata_factor (int) – For Numerical EOS.
max_poly_order_factor (int) – For Numerical EOS.
min_poly_order_factor (int) – For Numerical EOS.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.method (str) – ‘thermo_eos’
self.temperature (array) – Should be given in ascending order, or will be sorted. Unit: K
self.pressure (array) – Should be given in ascending order, or will be sorted. Unit: GPa
self.volume (array) – nPressure*nTemperature, same below. Equilibrium volumes. Unit: $\AA^{3}$
self.helmholtz (array) – Helmholtz free energy. Unit: kJ/mol
self.gibbs (array) – Gibbs free energy. Unit: kJ/mol
self.entropy (array) – Entropy. Unit: $J.mol^{-1}.K^{-1}$
self.c_p (array) – Constant pressure specific heat. Unit: $J.mol^{-1}.K^{-1}$
self.k_t (array) – Isothermal bulk modulus. Unit: GPa.
self.eos (list[Pymatgen EOS]) – nTemperature*1 list of Pymatgen EOS objects. EOSs used to fit HA free energy at constant temperature.
self.eos_method (str) – The name of EOS used.
self.gibbs_fit (list[Polynomial]) – nPressure list of Polynomials of fitted Gibbs free energy, for entropy and c_p.
self.fe_eos – Deprecated. Synonym of ‘self.eos’.
self.fe_eos_method – Deprecated. Synonym of ‘self.fe_eos_method’.

Raises:

Exception – If the number of HA calculations is less than 4.
Exception – If temperature or pressure is defined neither here nor during initialization.

expansion_vol(poly_order=[2, 3], plot=True, fit_fig='expansion_fit.png')

Fit the thermal expansion curve as polynomials and get thermal expansion coefficients at equilibrium volumes.

The volumetric thermal expansion coefficient at constant pressure:

\[\alpha_{V}(T) = \frac{1}{V(T)}\left(\frac{\partial V(T)}{\partial T}\right)_{p}\]

thermo_freq() and thermo_gruneisen() methods have fitted thermal expansion with Grüneisen model. Calling this will overwrite the fittings.

Note

For a good fitting. poly_order must be no smaller than 2. The length of self.temperature must no smaller than poly_order+1.

Parameters:

poly_order (list[int]) – Order of polynomials.
plot (bool) – Plot V-T curves to examine the goodness of fitting. The order with the largest $R^{2}$ is automatically selected.
fit_fig (str) – File name for the fitting plot.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.vol_fit (list) – 1*nPressure. List of Numpy polynomial object, the fitted volume V(T).
self.alpha_v (array) – nPressure*nTemperature. Expansion coefficients at equilibrium volumes.

bulk_modulus(**kwargs)

Calculate isothermal and adiabatic bulk moduli at equilibrium volumes.

The following equation is used:

\[K_{S} = K_{T} + \frac{\alpha^{2}_{V}VTK^{2}_{T}}{C_{V}}\]

Note

With thermo_eos(), must call expansion_vol() first.
With thermo_freq() and thermo_gruneisen(), unless the
expansion_vol() has been called, it will not change the results.

Parameters:

**kwargs – Only functions as a container of deprecated keywords. Not used.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.k_t (array) – Isothermal bulk modulus. Unit: GPa.
self.k_s (array) – nPressure*nTemperature. Adiabatic bulk modulus. Unit: GPa.

specific_heat()

Calculate constant volume or pressure specific heat at equilibrium volumes.

The following equation is used:

\[C_{p} - C_{V} = \alpha_{V}^{2}K_{T}VT\]

Note

With thermo_eos(), must call expansion_vol() first.
With thermo_freq() and thermo_gruneisen(), unless the
expansion_vol() has been called, it will not change the results.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.c_v (array) – nPressure*nTemperature. Constant volume specific heat. Unit: J/mol/K
self.c_p (array) – nPressure*nTemperature. Constant pressure specific heat. Unit: J/mol/K.

lattice_poly(poly_order=[2, 3], interp=0)

Fit anisotropic expansions, i.e., lattice parameters, by polynomial fitting lattice parameters with repect to volume. ($V^{1/3}$ is used.) More robust than thermo_hess().

Lattice will be expanded and aligned to the standard conventional cell. Volume constrain is applied. Minimum HA references: poly_order+2.

Note

Use interp to linearly interpolate reference geometries and insert reference Gibbs free energies. Only available to objects generated by thermo_freq() and thermo_gruneisen(). interp must less than self.ncalc.

Currently only the fitted self.lattice can be read from dump file.

Parameters:

poly_order (list[int]) – The order of polynomials used to fit lattice parameters as function of volumes.
interp (int) – Number of interpolated HA calculations.

Returns:

self – New attribute listed below.
self.sg (int) – International space group number.
self.lattice (array) – nPress*nTempt*3*3 lattice matrix at given temperature and pressure. Standard conventional cell is used.
self.latt_params (array) – nPress*nTempt*nLatt minimal set of lattice parameters. The standard conventional cell is used.

lattice_hess(interp=0, **kwargs)

Fit anisotropic expansions, i.e., lattice parameters, by the second-order perturbation of Gibbs free energy $G$.

\[G(\mathbf{v})=G_{0}(\mathbf{v_{0}})+\Delta\mathbf{v}\mathbf{H}\Delta\mathbf{v}^{T}\]

$\mathbf{v}$ is the row vector of inequivalent lattice parameters. $\mathbf{v_0}$ is the row vector of equilibrium lattice parameters. $\Delta\mathbf{v}$ is the difference between $\mathbf{v_0}$ and $\mathbf{v}$. $\mathbf{H}$ is the Hessian matrix.

Note

Inspired by Phys. Rev. Mater., 2019, 3, 053605.

The root-mean-squared deviations (RMSD) between $\Delta\mathbf{v}\mathbf{H}\Delta\mathbf{v}^{T}$ and G-G_{0} is minimized at constant temperature and pressure. Lattice will be expanded and aligned to the standard conventional cell.

This method requires a larger number of HA calculations to ensure a successful fitting, which depends on the symmetry of the system. Volume constrain is applied to reduce one inequivalent lattice parameter, but its dimension in $\mathbf{H}$ is kept.

Minimum HA references:

Hexagonal, trigonal and tetragonal: 4
Orthorhombic: 8
Monoclinic: 13
Triclinic: 26

Note

Use interp to linearly interpolate reference geometries and insert reference Gibbs free energies. Only available to objects generated by thermo_freq() and thermo_gruneisen(). interp must less than self.ncalc.

Currently only the fitted self.lattice can be read from dump file.

Parameters:

interp (int) – Number of interpolated HA calculations.
**kwargs – Passed to scipy.optimize.minimize() Passed keywords listed below.
method (str) – Default ‘BFGS’.
jac (str) – Default ‘3-point’.
hess – Scipy default settings.
hessp – Scipy default settings.
bounds – Scipy default settings.
constraints – Scipy default settings.
tol (float) – Default 0.05.
options – Scipy default settings.

Returns:

self – New attribute listed below.
self.sg (int) – International space group number.
self.lattice (array) – nPress*nTempt*3*3 lattice matrix at given temperature and pressure. Standard conventional cell is used.
self.latt_params (array) – nPress*nTempt*nLatt minimal set of lattice parameters. The standard conventional cell is used.
self.latt_cart (array) – nPress*nTempt*3 projected lengths of cell on Cartesian coordinates. The standard conventional cell is used.

property latt_params: Generate nPress*nTempt*nLattice minimum set of lattice parameters from lattice matrices of the standard conventional cell. Callable only if the equilibrium lattice parametes are fitted.

property latt_cart: Generate nPress*nTempt*3 set of projected lengths of cell on Cartesian coordinates from lattice matrices of the standard conventional cell. Callable only if the equilibrium lattice parametes are fitted.

expansion_latt(poly_order=[2, 3], Cartesian=False, plot=True, fit_fig='expansion_latt_{}.png')

Fit the thermal expansion curve as polynomials and get linear thermal expansion coefficients.

The linear thermal expansion coefficient at constant pressure:

\[\alpha_{l}(T) = \frac{1}{l(T)}\left(\frac{\partial l(T)}{\partial T}\right)_{p}\]

Note

For a good fitting. poly_order must be no smaller than 2. The length of self.temperature must no smaller than poly_order+1.

Parameters:

poly_order (list[int]) – Order of polynomials.
Cartesian (bool) – Fit linear thermal expansion rate with reference to Cartesian coordinate. Useful for monoclinic and triclinic systems.
plot (bool) – Plot a-T curves to examine the goodness of fitting. The order with the largest $R^{2}$ is automatically selected.
fit_fig (str) – File name for the fitting plots. Must include the format string for string variables.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.latt_fit (list) – nPressure*nLattice. List of Numpy polynomial object, the fitted lattice parameters. nLattice=3 if Cartesian=True.
self.alpha_l (array) – nPressure*nTemperature*nLattice. Expansion coefficients. nLattice=3 if Cartesian=True.

eos_fit(volume, energy, method, **kwargs)

Fit energy-volume relationship by equation of states. Expected to be called internally.

Parameters:

volume (array[float]) – Unit: Angstrom^3
energy (array[float]) – Unit: kJ/mol
method (str) – Name of EoS used. Consistent with Pymatgen.
order (int) – For DeltaFactor / Polynomial methods.
min_ndata_factor (int) – For the NumericalEOS method.
max_poly_order_factor (int) – For the NumericalEOS method.
min_poly_order_factor (int) – For the NumericalEOS method.

Returns:

eos (Pymatgen EOS) – The fitted equation of state.

freq_polynomial_fit(order)

Fit phonon frequencies as polynomial functions of volumes. $\omega(\Delta V)$ is fitted. Expected to be called internally.

Parameters:

order (list[int] | array[int]) – The order of polynomials used.

Returns:

r2tot (array) – nQpoint*nMode array of R^2.
r2avg (array) – Averaged R^2 at each q point.
self – New attributes listed below.
fit_order (array) – Same as order.
freq_fit (list) – nQpoint*nMode list of fitted polynomials.

freq_gruneisen_fit()

Fit mode-specific Grüneisen parameters from sampled HA calculations. Expected to be called internally.

Returns:

rmsdtot (array) – nQpoint*nMode array of R^2.
rmsdavg (array) – Averaged R^2 at each q point.
self – New attributes listed below.
gru_fit (array) – nQpoint*nMode array of Grüneisen parameters.

_combine_data(ha_list, mode_sort_tol)

Combine the HA calculation data and rearrange it in the ascending order of volumes.

Parameters:

ha_list (list[Harmonic]) – List of harmonic objects.
mode_sort_tol (float | None)

Returns:

close_overlap (array[int]) – nQpoint*nMode_ref*nCalc*nMode_sort boolean array (only 1 and 0). Whether close overlap is identified between the n-1 (ref) and the n (sort) calculations. The first Calc is useless (always 0).
dot_pdt (array[float]) – nQpoint*(nCalc-1) array, dot products between n-1 and n calculations.
self – Attributes listed below.
combined_volume (array[float]) – 1*nCalc
combined_struc (list[CStructure]) – 1*nCalc
combined_u0 (array[float]) – 1*nCalc
combined_freq (array[float]) – nQpoint*nMode*nCalc
combined_symm (array[str]) – nQpoint*nMode*nCalc
nqpoint (int)
qpoint (array[float]) – nQpoint*4
nmode (int)
natom (int)

Raises:

Exception – If number of q points, modes or atoms are not consistent across the HA calculations.

static _phonon_continuity(freq, eigvt, symm=None, mode_sort_tol=0.4)

Rearrange phonon modes by their continuity. If the difference between the maximum scalar product of corresponding eigenvectors (normalized to 1) and scalar products of other modes is less than 0.4, warning is printed due to the potential overlap of modes. Adopted from CRYSTAL17.

Note

Erba, J. Chem. Phys., 2014, 141, 124115.

Parameters:

freq (array[float]) – Phonon frequencies. Unit: THz
eigvt (array[float]) – Eigenvectores normalized to 1
symm (array[str]) – Irreducible representations in Mulliken symbols.
mode_sort_tol (float) – The threshold of close mode overlaps.

Returns:

freq (array[float]) – Sorted phonon frequencies
eigvt (array[float]) – Sorted eigenvectores
close_overlap (array[float]) – ncalc*nmode*nmode boolean matrix (only 1 and 0). Whether close overlap is identified between the previous calculation (2nd dimension) and the current one (3rd).
dot_pdt (array[float]) – 1*(nCalc-1) array, dot products between n-1 and n calculations.

_clean_thermoprop(): Clean fitted data. Called by ‘thermo_*’ methods only.

_interp_HA(interp, thermo)

Interpolate harmonic phonons for lattice fitting. Interpolation is done by reducing the volume steplength by half, and add points closest to the boundary first.

Lattice will be expanded to conventional cell and aligned. $a$ will be aligned to $x$ and $b$ in $xOy$ plane.

Note

Only available to objects generated by thermo_freq() and thermo_gruneisen().

Parameters:

interp (int) – Number of interpolations. Must be smaller than self.ncalc.
thermo (bool) – Whether to run HA thermodynamics.

Returns:

sg (int) – International space group number.
vscale (float) – Volume ratio between conventional and primitive cell.
rot (array) – Rotation matrix from invernal convention to input / standard conventional cell.
lmx (array) – (nCalc+nInterp)*3*3 array of lattice matrices. The first nCalc elements are sampled lattices.
Gibbs (array) – (nCalc+nInterp)*nPress*nTempt array of Gibbs free energy in kJ/mol.

from_HA_files(*input_files, imaginary_tol=0.0001, q_overlap_tol=0.0001, mode_sort_tol=0.4): Deprecated. Call from_file() directly.

from_QHA_file(input_file, imaginary_tol=0.0001, q_overlap_tol=0.0001, mode_sort_tol=0.4): Deprecated. Call from_file() directly.

from_phonopy_files(phono_yaml, struc_yaml=None, edft=None, scale=1.0, imaginary_tol=0.0001, q_overlap_tol=0.0001, q_id=None, q_coord=None): Deprecated. Call from_file() directly.

expansion_lin(poly_order=[2, 3], interp=0): Deprecated and not compatibable. Only prints error message.

class Phonopy

Bases: object

Deprecated.

classmethod read_structure(file)

classmethod read_frequency(file, q_id=None, q_coord=None)

classmethod write_force_constants(hessfile='HESSFREQ.DAT', phonopyfile='FORCE_CONSTANTS')

Write Phonopy/VASP FORCE_CONSTANTS file by CRYSTAL HESSFREQ.DAT file.

For example, to convert the calculation ‘example’ with a 4*4*4 supercelland get phonon frequencies at Gamma point, use the following code:

>>> from CRYSTALpytools.thermodynamics import Phonopy
>>> Phonopy.write_force_constants(hessfile='example.HESSFREQ')
>>> phonopy --crystal --qpoints='0 0 0' -c example.out --dim='4 4 4' --readfc

Parameters:

hessfile (str) – The HESSFREQ.DAT file
phonopyfile (str) – The output name

class Output

Bases: object

Deprecated. Only prints out error messages.

classmethod write_HA_result()

classmethod write_QHA_combinedata()

classmethod write_QHA_sortphonon()

classmethod write_QHA_eosfit()

classmethod write_QHA_polyfit()

classmethod write_QHA_thermofreq()

classmethod write_QHA_thermogru()

classmethod write_QHA_thermoeos()

classmethod write_expansion_vol()

classmethod write_bulk_modulus()

classmethod write_specific_heat()

classmethod write_expansion_latt()