\(\newcommand{\AA}{\text{Å}}\)

CRYSTALpytools.thermodynamics module

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

class Mode(rank=0, frequency=[], volume=[], eigenvector=[], symmetry='')

Bases: object

Defines a vibrational mode and do analysis per mode.

Parameters:

rank (int) – The rank of the mode object, from 1.
frequency (array[float]) – Frequencies of the mode (Ncalc*1). Unit: THz. Note: NOT angular frequency, which is frequency*2pi.
volume (array[float]) – Lattice volumes of harmonic calculations (Ncalc*1). Unit: Angstrom^3
eigenvector (array[float]) – Corresponding normalized eigenvectors (Ncalc*Natom*3).
symmetry (str) – Irreducible representations in Mulliken symbols.

get_zp_energy()

Get the zero-point energy of a single mode. ncalc = 1 only.

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

Returns:: self.zp_energy (float) – Zero-point energy. Unit: KJ/mol

get_u_vib(temperature)

Get the vibration contribution to internal energy (including zero-point energy) of a single mode. ncalc = 1 only.

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

Parameters:: temperature (float) – Temperature where the quantity is computed. Unit: K
Returns:: self.u_vib (float) – Vibration contribution to internal energy. Unit: KJ/mol

get_entropy(temperature)

Get the entropy of a single mode. ncalc = 1 only.

\[S_{i,\mathbf{q}}\left(T\right)=k_{B}\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\}\]

Parameters:: temperature (float) – Unit: K
Returns:: self.entropy (float) – Entropy. Unit: J/mol*K

get_c_v(temperature)

Get the constant volume specific heat of a single mode. ncalc = 1 only.

\[C^{V}_{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} }\]

Parameters:: temperature (float) – Unit: K
Returns:: self.c_v (float) – Constant volume specific heat. Unit: J/mol*K

polynomial_fit(order)

Fit phonon frequency as the polynomial function of volume. ncalc > 1 only.

Note

To improve the accuracy of fittings, \(\Delta\omega(\Delta V)\) is fitted as a polynomial function without the constant term.

\(\Delta V=V-V_{min}\) is used so HA phonons of the most compact structure is kept. See ‘FIXINDEX’ keyword in CRYSTAL manual for further information.

Parameters:

order (array[int]) – Orders of polynomials.

Returns:

self.poly_fit (Dict[int, NumPy Polynomial]) – Key - orders of power, Value - fitted NumPy polynomials
self.poly_fit_rsquare (Dict[int, float]) – Key - orders of power, Value - goodness of fittings, characterized by R^2.

get_gruneisen(order, volume)

Return to mode Grüneisen parameter. ncalc > 1 only.

\[\gamma = -\frac{V}{\omega(V)}\frac{\partial\omega}{\partial V}\]

Note

order is used only for the dict key. In principle, order=1 for the Grüneisen model (see QHA.thermo_gruneisen()). Use other order numbers for frequencies fitted by QHA.thermo_freq().

Parameters:

order (int | list) – See polynomial_fit
volume (float | array) – Typically the equilibrium volume

Returns:

self.gruneisen (dict) – Key, order; Value, Gruneisen parameter

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

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.

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')

from_file(output_name, read_eigvt=False, imaginary_tol=-0.0001, q_overlap_tol=0.0001, qha_index=0)

Generate the Harominc object from a HA output file. Imaginary modes and overlapped q points are forced to be cleaned.

Parameters:

output_name (str) – Name of the output file.
read_eigvt (bool) – Whether to read eigenvectors from output.
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
qha_index (int) – If the input is a QHA file, specify the index of calculation to read (from 0).

Returns:

self (Harmonic) – Attributes listed below.
structure (CStructure) – Extended Pymatgen structure class. The supercell expanded by keyword ‘SCELPHONO’ is reduced.
natom (int) – Number of atoms in the reduced cell.
volume (float) – Volume of the reduced cell. Unit: \(\AA^{3}\)
edft (float) – Internal energy
nqpoint (int) – Number of q points
qpoint (list) – nQpoint*2 list. The first element is 1*3 array of fractional coordinates in reciprocal space and the second is its weight.
nmode (array[int]) – Number of modes at every q point.
mode (list[Mode]) – List of mode objects at all the qpoints.

Raises:

ValueError – If a QHA output file is read.

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)

Build a Harmonic object from Phonopy ‘band.yaml’ or ‘qpoints.yaml’ file.

Note

Mode symmetry and eigenvector are currently not available.

Note

q_id and q_coord should not be set simultaneously. If set, q_id takes priority and q_coord is ignored. If both are none, all the points will be read.

Parameters:

phono_yaml (str) – Phonopy band.yaml or qpoint.yaml file
struc_yaml (str) – Phonopy phonopy.yaml or phonopy_disp.yaml file. Needed only if a qpoint.yaml file is read.
edft (float) – DFT energy. Unit: kJ/mol
scale (float) – Scaling factor of phonon frequency.
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 (list[int]) – Specify the id (from 0) of q points to be read. nqpoint*1 list.
q_coord (list[list]) – Specify the coordinates of q points to be read. nqpoint*3 list.

Returns:

self (Harmonic) – Attributes listed below.
structure (CStructure) – Extended Pymatgen structure class. The supercell expanded by keyword ‘SCELPHONO’ is reduced.
natom (int) – Number of atoms in the reduced cell.
volume (float) – Volume of the reduced cell. Unit: \(\AA^{3}\).
edft (float) – Internal energy.
nqpoint (int) – Number of q points.
qpoint (list) – nQpoint*2 list. The first element is 1*3 array of fractional coordinates in reciprocal space and the second is its weight.
nmode (array[int]) – Number of modes at every q point.
mode (list[Mode]) – List of mode objects at all the qpoints.

Raises:

Exception – If the length unit in yaml file is neither ‘au’ nor ‘angstrom’.
Exception – If q point is not found.

from_frequency(edft, qpoint, frequency, eigenvector, symmetry, structure=None, natom=None, volume=None, imaginary_tol=-0.0001, q_overlap_tol=0.0001, ignore_natom=False)

Generate a Harmonic object by specifying frequency and eigenvector. Imaginary modes and overlapped q points are forced to be cleaned.

Parameters:

edft (float) – Electron total energy in kJ/mol.
qpoint (list[list[array[float], float]]) – nQpoint*2 list of: 1*3 array of fractional coordinate and weight of qpoint (maximum value normalized to 1).
frequency (array[float]) – Array of frequencies. Unit: THz
eigenvector (array[float]) – Normalized eigenvectors to 1.
symmetry (array[str]) – Irreducible representations in Mulliken symbols.
structure (CStructure) – Extended Pymatgen structure class. The supercell expanded by keyword ‘SCELPHONO’ is reduced.
natom (int) – Number of atoms in the reduced cell.
volume (float) – Volume of the reduced cell. Unit: \(\AA^{3}\)
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
ignore_natom (bool) – Developer only.

Note

The user should define either structure or natom + volume.

Returns:

self (Harmonic) – Attributes see from_file() and from_phonopy.

Raises:

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

thermodynamics(sumphonon=True, **kwargs)

Calculate the thermodynamic properties (zp_energy, u_vib, entropy, c_v and Gibbs and Helmholtz free energy) of the HA system at all qpoints and the whole temperature/pressure range.

Other parameters are the sum of corresponding attributes of all the Mode objects. The Helmholtz and Gibbs free energies are defined as:

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

Parameters:

temperature (array|float) – Optional Unit: K
pressure (array|float) – Optional Unit: GPa
sumphonon (bool) – Whether to sum up the phonon contributions across the sampled q points and take weighted-average.

Returns:

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

Note

If sumphonon = True, nqpoint = 1 but its dimension in attribute is kept for consistency.

Raises:: AttributeError – If temperature and pressure are defined neither here nor during initialization

write_HA_result()

_thermo_gibbs(temperature, pressure): Used only for QHA gibbs free energy minimization with Quasi_harmonic.thermo_freq. Not for general use.

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

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) – Unit: K
pressure (array|float) – Unit: GPa
filename (str) – Name of the output file. Valid if write_out = True.

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_HA_files(*input_files, imaginary_tol=-0.0001, q_overlap_tol=0.0001, mode_sort_tol=0.4)

Read data from individual HA calculation outputs. Imaginary modes and overlapped q points are forced to be cleaned.

Parameters:

*input_files (str) – Harmonic phonon output filenames. Extendable.
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
mode_sort_tol (float | None) – The threshold of close mode overlaps. If none, do not sort modes and eigenvector is not read.

Returns:

self (Quasi_harmonic) – New Attributes listed below
self.ncalc (int) – Number of HA phonon calculations.
self.combined_phonon (list[Harmonic]) – List of Harmonic objects.
self.combined_volume (list[float]) – Volumes. Unit: Angstrom^3
self.combined_edft (list[float]) – DFT total energies. Unit: KJ/mol
self.combined_mode (list[Mode]) – List of mode objects.

from_QHA_file(input_file, imaginary_tol=-0.0001, q_overlap_tol=0.0001, mode_sort_tol=0.4)

Read data from a single QHA calculation at Gamma point. Imaginary modes and overlapped q points are forced to be cleaned.

Parameters:

input_file (str) – Only 1 QHA file is permitted.
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
mode_sort_tol (float | None) – The threshold of close mode overlaps. If none, do not sort modes.

Returned attributes are consistent with Quasi_harmonic.from_HA_files.

Raises:: ValueError – If multiple files are defined.

from_phonopy_files(*phono_yaml, struc_yaml=None, edft=None, imaginary_tol=-0.0001, q_overlap_tol=0.0001, q_id=None, q_coord=None)

Build a QHA object from Phonopy ‘band.yaml’ or ‘qpoints.yaml’ file.

Note

Mode symmetry and eigenvector are currently not available.

Note

q_id and q_coord should be set once and are applicable to all the yaml files.

Parameters:

phono_yaml (str) – ncalc*1 list of Phonopy band.yaml or qpoint.yaml files
struc_yaml (list[str]) – ncalc*1 list of Phonopy phonopy.yaml or phonopy_disp.yaml files. Needed only if a qpoint.yaml file is read.
edft (list[float]) – ncalc*1 list / array of DFT energies.
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 (list[int]) – See Harmonic.from_phonopy.
q_coord (list[list]) – See Harmonic.from_phonopy.

Returned attributes are consistent with Quasi_harmonic.from_HA_files.

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)]\]

Parameters:

eos_method (str) – EOS used to fit DFT total energy and Helmholtz free energy (to get bulk modules).
poly_order (array[int]) – The order of polynomials used to fit frequency as the function of volumes.
min_method (string, optional) – Minimisation algorithms.
volume_bound (tuple-like) – volumes. Unit: Angstrom^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.

Note

Valid entries of eos_method are consistent with PyMatGen eos module.
Parameterized and tested algorithms for min_method:
- BFGS(no boundary)
- L-BFGS-B(with boundary)

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.c_v (array) – Constant volume specific heat. Unit: \(J.mol^{-1}.K^{-1}\)
self.e0_eos (Pymatgen EOS) – Pymatgen EOS object. EOS used to fit DFT energy.
self.e0_eos_method (str) – Name of the EOS.

Raises:

ValueError – 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)

Grüneisen parameters and related properties. 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}\]

Note

The Grüneisen model is used to fit frequencies, equivalent to using self.thermo_freq(poly_order=1).

For arguments, see self.thermo_freq.

Returns:

self (Quasi_harmonic) – New attributes listed below. Other attributes are the same as self.thermo_freq.
self.gruneisen (array) – npressure*ntemperature, same below. Macroscopic Grüneisen parameter. Temperature should > 0.
self.alpha_vgru (array) – Thermal expansion coefficient by Grüneisen method.
self.c_pgru (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_sgru (array) – Adiabatic bulk modulus by Grüneisen method. Unit: GPa.

thermo_eos(eos_method='birch_murnaghan', poly_order=[2, 3], 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.

Entropy is obtained by taking the derivation 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}\]

Note

poly_order should >= 2.

For arguments, see self.thermo_freq.

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.c_p (array) – Constant pressure specific heat. Unit: \(J.mol^{-1}.K^{-1}\)
self.fe_eos (list[Pymatgen EOS]) – nTemperature*1 list of Pymatgen EOS objects. EOSs used to fit HA free energy at constant temperature.
self.fe_eos_method (str) – The name of EOS used.

Raises:

Exception – If the number of HA calculations is less than 4.
ValueError – 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 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}\]

For thermal expansion with Grüneisen model, please refer to the thermo_gruneisen() method.

Parameters:

poly_order (list[int]) – method = ‘polynomial’, order of polynomials.
plot (bool) – Plot V-T curves to examine the goodness of fitting. An interactive window will pump out to let user to specify the optimial fitting.
fit_fig (str) – File name for fittings. A temperal figure is printed to help the user choose the optimal fitting.

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(adiabatic=True, **kwargs)

Calculate isothermal and adiabatic bulk moduli at equilibrium volumes.

The following equations are used:

\[ \begin{align}\begin{aligned}K_{T}(p;T) = V(p;T)\left(\frac{\partial^{2}F(V;T)}{\partial V^{2}}\right)_{T}\\K_{S} = K_{T} + \frac{\alpha^{2}_{V}VTK^{2}_{T}}{C_{V}}\end{aligned}\end{align} \]

To get \(K_{T}\), Helmholtz free energy is fit as isothermal EOSs. For self.thermo_eos(), that means doing nothing; For self.thermo_freq(), EOS fitting is required, whose form is the same as EOS used for \(E_{0}\).

Parameters:

adiabatic (bool) – Whether to fit adiabatic bulk modulus. Thermal expansion coefficient needed.
**kwargs – See below.
order – (int): To restore EOS.
min_ndata_factor – (int): To restore EOS.
max_poly_order_factor – (int): To restore EOS.
min_poly_order_factor – (int): To restore EOS.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.k_t (array) – nPressure*nTemperature, same below. Isothermal bulk modulus. Unit: GPa.
self.k_s (array) – 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\]

Returns:

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

Note

This method fits self.c_p by self.c_v when thermo_freq and thermo_gruneisen was used. self.c_v is obtained by when thermo_eos is used.

expansion_lin(poly_order=[2, 3], interp=None)

Fit linear expansions of lattice parameters by the 2-order Taylor expansion.

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

\(G\) is Gibbs free energy. \(\mathbf{p}\) is the vector of lattice parameters. \(\Delta\mathbf{p}\) means the difference between the fitted and equilibrium lattice parameters. \(\mathbf{H}\) is the Hessian of \(G\) and displacements along lattice parameters.

The RMS deviations (RMSD) of the following equation is minimized at constant temperature and pressure. But deviations from equilibrium volume might occur. RMSD of Gibbs free energy is available in output file only.

\[\mathbf{p_{0}} = \min\left\{\Delta\mathbf{p}^{T}\mathbf{H}\Delta\mathbf{p} - [G(\mathbf{p})-G_{0}(T,p)]\right\}\]

Note

Raimbault, V. Athavale and M. Rossi, Phys. Rev. Materials, 2019, 3, 053605.

This method requires a larger number of HA calculations to ensure a small RMSD. Typically the number of HA calculations should follow the equation below, otherwise the warning massage is given.

\[n_{HA} \geq n_{latt} + \sum_{i=1}^{n_{latt}}i\]

\(n_{latt}\) is the lenth of the minimial set of lattice parameters. The optimized lattice parameters at DFT level are used for fitting.

To avoid warnings, use interp to linearly interpolate the lattice parameters and to get thermodynamic properties from the fitted QHA object.

Parameters:

poly_order (list[int]) – Order of polynomials used to fit the linear expansion coefficients. The optimal fit across the sampled temperature and pressure range of a certain lattice parameter is automatically chosen based on \(R^{2}\).
interp (int) – Number of interpolated geometries. All the HA geometries are used besides the interpolated ones.

Returns:

self (Quasi_harmonic) – New attributes listed below
self.lattice (array) – nPressure*nTemperature*nLattice. The equilibrium values of minimal set of lattice parameters.
self.latt_fit (list) – nPressure*nLattice. Numpy polynomial object, the fitted a(v). Linear part only.
self.alpha_latt (array) – nPressure*nTemperature*nLattice. Linear expansion coefficients. Linear part only.

_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:

combined_phonon (list[Harmonic])
combined_volume (array[float])
combined_edft (array[float])
combined_mode (list[Mode])

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[bool]) – ncalc*nmode*nmode. Whether close overlap is identified between the previous calculation (2nd dimension) and the current one (3rd).

eos_fit(volume, energy, method, write_out=True, **kwargs)

Fit energy-volume relationship by equation of states.

Parameters:

volume (array[float]) – Unit: Angstrom^3
energy (array[float]) – Unit: kJ/mol
method (str) – Name of EoS used. Consistent with Pymatgen.
write_out (bool) – Whether to print EOS information.
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.
eos_method (string) – Name of the fitted equation of state

freq_polynomial_fit(order)

Fit phonon frequencies as polynomial functions of volumes.

Parameters:: order (list[int] | array[int]) – The order of polynomials used.
Returns:: self.fit_order (int) – The optimal order of polynomial fit.

Please also refer to self.poly_fit and self.poly_fit_rsquare attributes of Mode class.

_get_harmonic_phonon(volume)

Get numerical phonon frequencies from fitted analytical expressions and generate harmonic phonon objects. Not a standalone method.

Parameters:: volume (float) – Unit: Angstrom^3
Returns:: ha (Harmonic) – Harmonic phonon object with numerical data.
Raises:: Exception – If frequency is not fitted as function of volumes.

_minimize_gibbs(volume, temperature, pressure)

Get Gibbs free energy from the Harmonic phonon object. Used only for minimizing \(G(V;T, p)\) by SciPy.

Parameters:

volume (float)
temperature (float)
pressure (float)

Returns:

ha.gibbs (float) – Gibbs free energy. Unit: KJ/mol

static _poly_no_cst(param, x, y): Define a polynomial \(\Delta f(\Delta x)\) without constant term. Orders low to high. For SciPy. Functions of vectors not supported.

_clean_attr(): When temperature / pressure are changed, thermodynamic attributes are removed to keep consistency.

static _minimize_latt(x, fe_eq, latt_ref, fe_ref): Minimize the RMSD between \(\mathbf{p}^{T}\mathbf{Hp}\) and the difference of Gibbs free energy. For Scipy. For fitting lattice parameters in self.expansion_vol.

class Phonopy

Bases: object

The convertor between Phonopy and CRYSTALpytools file formats

classmethod read_structure(file)

Read geometry from Phonopy band.yaml or phonopy.yaml or phonopy_disp.yaml files.

Parameters:: file (str) – Phonopy yaml file
Returns:: struc (Pymatgen Structure)
Raises:: Exception – If the length unit in yaml file is neither ‘au’ nor ‘angstrom’.

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

Read phonon frequency from Phonopy band.yaml or qpoints.yaml files. Frequency units must be THz (default of Phonopy).

Parameters:

file (str) – Phonopy yaml file
q_id (list[int]) – Specify the id (from 0) of q points to be read. nqpoint*1 list.
q_coord (list[list]) – Specify the coordinates of q points to be read. nqpoint*3 list.

q_id and q_coord should not be set simultaneously. If set, q_id takes priority and q_coord is ignored. If both are none, all the points will be read.

Returns:

qpoint (list) – natom*2 list. 1st element: 3*1 array. Fractional coordinates of q points; 2nd element: float. Weight
frequency (array) – nqpint*nmode array. Phonon frequency in THz.

Raises:

Exception – If (some of) q point is not found.

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

Deal with output data file

classmethod write_HA_result(ha)

Write harmonic phonon information.

Parameters:: ha (Harmonic) – CRYSTALpytools.thermodynamic.Harmonic object

classmethod write_QHA_combinedata(qha)

Write QHA combined phonon information.

Parameters:: qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object

classmethod write_QHA_sortphonon(qha, close_overlap)

Write QHA phonon sort information.

Parameters:

qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.
close_overlap (array[int]) – nqpoint*nmode_ref*ncalc*nmode_sort. Number of close overlaps.

classmethod write_QHA_eosfit(qha, eos, method)

Write QHA phonon eos fit information.

Parameters:

qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.
order (list[int]) – Orders of polynomials
method (str) – Name of EoS used.

classmethod write_QHA_polyfit(qha, order, rsquare_q)

Write QHA phonon polynomial fit information.

Parameters:

qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.
order (list[int]) – List of polynomial orders.
rsquare_q (array) – Nqpoint*Norder list. Overall goodness at a q point.

classmethod write_QHA_thermofreq(qha, min_method, volume_bound)

Write QHA thermodynamics information (frequency fitting).

Parameters:: qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.

classmethod write_QHA_thermogru(qha)

Write QHA thermodynamics information (Grüneisen fitting).

Parameters:: qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.

classmethod write_QHA_thermoeos(qha)

Write QHA thermodynamics information (EOS fitting).

Parameters:: qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.

classmethod write_expansion_vol(qha, fit_order, fit_rs)

Write volumetric thermal expansions.

Parameters:

qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.
fit_order (int) – The order of polynomial used for fitting.
fit_rs (list[float]) – R^2 of fitting. nPressure*1 list.

classmethod write_bulk_modulus(qha, adiabatic)

Write bulk moduli.

Parameters:

qha (Quasi_harmonic) –

CRYSTALpytools.thermodynamic.Quasi_harmonic: object.
adiabatic (bool): Whether the adiabatic bulk modulus \(K_{S}\): is fitted.

classmethod write_specific_heat(qha)

Write bulk moduli.

Parameters:: qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.

classmethod write_expansion_latt(qha, e_err, fit_order, r_square)

Write linear expansion information.

Parameters:

qha (Quasi_harmonic) – CRYSTALpytools.thermodynamic.Quasi_harmonic object.
e_err (array) – RMS deviation of Gibbs free energy at T and p.
fit_order (array) – The order of polynomials used for fitting lattice parameter.
r_square (array) – R^2 of fitting. nLatt*nPress