Lattice Dynamics 1: Harmonic Approximation

The thermodynamics module is developed for both harmonic approximated (HA) and quasi-harmonic approximated (QHA) lattice dynamics. For simplicity, this part is focused only on HA phonons. For QHA, please refer to Lattice Dynamics 2: Quasi-Harmonic Approximation example book in the same catagory. For phonon band and density of states, please refer to the phonons module.

Instantiation

No data is read at the object instantiation step. Instead, important controlling parameters are defined.

Define the room temperature and pressure with temperature and pressure (optional) and specify whether to write data into an output file (filename) and whether to automatically run the calculation (autocalc)

[1]:

from CRYSTALpytools.thermodynamics import Harmonic

obj = Harmonic(filename='ha_paracetamolG.txt', autocalc=False,
               temperature=[298.15], pressure=[0.10132500E-3])

File I/O

Currently the user can read geometries, internal energies, phonon frequencies, eigenvectors from CRYSTAL or phonopy (not internal energies). CRYSTALpytools helps to calculate important thermodynamic properties at specified temperature and pressure.

The ‘from_file()’ method.

Please note that the get_phonon() method defined in crystal_io.Crystal_output is a general-propose method designed to be called (in most cases) internally. Users are always recommended to use the from_file() method after instantiation.

Call the thermodynamics() method after reading data. Temperature and pressure can also be specified here.

[2]:

obj.from_file('ha_paracetamolG.out')

print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode[0]))
print('The first non-translational mode frequnecy at Gamma: {:.4f} THz'.format(
    obj.frequency[0, 3]))

Number of q / k points: 1
Number of modes at Gamma: 240
The first non-translational mode frequnecy at Gamma: 1.0304 THz

This method also accepts outputs by ‘QHA’ keyword of CRYSTAL. qha_index should be specified to substract the desired data (starting from 0).

Get the equilibrium geometry of QHA calculation of Al\(_{2}\)O\(_{3}\).

[3]:

from CRYSTALpytools.thermodynamics import Harmonic

obj = Harmonic(filename=None, autocalc=False).from_file('qha_corundumG.out',
                                                        qha_index=3)
print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode[0]))
print('The first non-translational mode frequnecy at Gamma: {:.4f} THz'.format(
    obj.frequency[0, 3]))

Number of q / k points: 1
Number of modes at Gamma: 90
The first non-translational mode frequnecy at Gamma: 4.4021 THz

/home/huanyu/apps/anaconda3/envs/crystal_py3.9/lib/python3.9/site-packages/CRYSTALpytools/thermodynamics.py:379: UserWarning: QHA output found, reading the 3 HA calculation from file.
  warnings.warn("QHA output found, reading the {:d} HA calculation from file.".format(qha_index))

Also might be useful to read phonon dispersions in the first Brillouin Zone for a better description of thermodynamics.

[4]:

from CRYSTALpytools.thermodynamics import Harmonic

obj = Harmonic(filename=None, autocalc=False).from_file('ha_paracetamolDisp.out')
print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode[0]))
print('The first non-translational mode frequnecy at Gamma: {:.4f} THz'.format(
    obj.frequency[0, 3]))

/home/huanyu/apps/anaconda3/envs/crystal_py3.9/lib/python3.9/site-packages/CRYSTALpytools/thermodynamics.py:556: UserWarning: Negative frequencies detected.
  Calculated thermodynamics might be inaccurate. Negative frequencies will be substituted by 0.
  self = PhononBASE.clean_imaginary(self, threshold=imaginary_tol)

Number of q / k points: 216
Number of modes at Gamma: 240
The first non-translational mode frequnecy at Gamma: 1.2303 THz

The ‘from_phonopy()’ method

Output files of phonopy calculations in ‘yaml’ format are accepted. Use phono_yaml for frequency files (‘band.yaml’ or ‘qpoint.yaml’) and struc_yaml for geometry (‘phonopy.yaml’ or ‘phonopy_disp.yaml’).

DFT total energy is unknown, so it is important to specified that. A warning message is triggered in the following case.

NOTE

Phonopy is not in the dependency list of CRYSTALpytools.

[5]:

from CRYSTALpytools.thermodynamics import Harmonic

obj = Harmonic(filename=None, autocalc=False).from_phonopy(
    'ha_freq.yaml', 'ha_struc.yaml')

print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode[0]))
print('The first non-translational mode frequnecy at Gamma: {:.4f} THz'.format(
    obj.frequency[0, 3]))

/home/huanyu/apps/anaconda3/envs/crystal_py3.9/lib/python3.9/site-packages/CRYSTALpytools/thermodynamics.py:468: UserWarning: DFT energy is set to 0.
  warnings.warn('DFT energy is set to 0.')

Number of q / k points: 1
Number of modes at Gamma: 240
The first non-translational mode frequnecy at Gamma: 1.0301 THz

The ‘Phonopy.write_force_constants’ method

The numerical Hessian is printed out in a formatted ‘HESSFREQ.DAT’ file by CRYSTAL, which can be converted into a Phonopy FORCE_CONSTANTS file by the thermodynamics.Phonopy.write_force_constants() method and used by phonopy in future analysis.

NOTE

Phonopy is not in the dependency list of CRYSTALpytools.

[6]:

from CRYSTALpytools.thermodynamics import Phonopy

Phonopy.write_force_constants(hessfile='ha_paracetamolG.HESSFREQ',
                              phonopyfile='FORCE_CONSTANTS')
# call phonopy
! phonopy --crystal --qpoints='0 0 0' -c ha_paracetamolG.out --dim='1 1 1' --readfc > /dev/null 2>&1

Frequencies are computed by phonopy in the ‘qpoints.yaml’ file of the same directory.

Thermodynamics

Gamma point

Call the thermodynamics() method after reading data. Temperature and pressure can also be specified here if not specified at instantiation stage.

[7]:

from CRYSTALpytools.thermodynamics import Harmonic

obj = Harmonic(filename=None, autocalc=False, temperature=[298.15],
               pressure=[0.10132500E-3])
obj.from_file('ha_paracetamolG.out').thermodynamics()

print('DFT total energy (EL) = ', obj.edft, 'KJ/mol')
print('Helmholtz free energy (EL+E0+ET-TS) = ', obj.helmholtz[0, 0], ' KJ/mol')
print('Gibbs free energy (EL+E0+ET+pV-TS) = ', obj.gibbs[0, 0, 0], ' KJ/mol')
print('Zero-point energy (E0) = ', obj.zp_energy[0], ' KJ/mol')
print('Vibrational contribution to interla energy - E0 (ET) = ', obj.u_vib[0, 0] - obj.zp_energy[0], ' KJ/mol')
print('Entropy*Temperature (TS) = ', obj.entropy[0, 0] * 298.15, ' J/mol')
print('Heat capacity = ', obj.c_v[0, 0], ' J/mol*K')
print('Entropy = ', obj.entropy[0, 0], ' J/mol*K')

DFT total energy (EL) =  -5402456.254965135 KJ/mol
Helmholtz free energy (EL+E0+ET-TS) =  -5400800.797249377  KJ/mol
Gibbs free energy (EL+E0+ET+pV-TS) =  -5400800.7512229495  KJ/mol
Zero-point energy (E0) =  1759.100392885878  KJ/mol
Vibrational contribution to interla energy - E0 (ET) =  107.20364114600739  KJ/mol
Entropy*Temperature (TS) =  210846.31827405942  J/mol
Heat capacity =  665.5752603999974  J/mol*K
Entropy =  707.1820166830771  J/mol*K

References from ‘ha_paracetamolG.out’ (line 11703~11738) are attached below. Here are several findings:

Parameters from frequencies agree well with reference data. The choice of physical constants and different unit conversion coefficients might lead to the discrepancy, where the physical constants play a major role (See below).
The difference in DFT total energy lies in different coefficients for unit conversion. In CRYSTALpytools, constants from `scipy.constants <https://docs.scipy.org/doc/scipy/reference/constants.html>`__ is used.
The PV term is neglected at room condition. Helmholtz free energy is the focus of this tutorial since the Harmonic class was initially developed for QHA fittings. To calculate Gibbs free energy at HA level, use this equation: \(G(T,p) = F(T) + pV\)

*******************************************************************************

HARMONIC VIBRATIONAL CONTRIBUTIONS TO THERMODYNAMIC FUNCTIONS AT GIVEN
TEMPERATURE AND PRESSURE:

(EL = ELECTRONIC ENERGY
 E0 = ZERO-POINT ENERGY
 ET = THERMAL CONTRIBUTION TO THE VIBRATIONAL ENERGY
 PV = PRESSURE * VOLUME
 TS = TEMPERATURE * ENTROPY)

                         AU/CELL             EV/CELL                 KJ/MOL
EL            :   -2057.686915559598  -55992.507576455653    -5402456.23545757
E0            :       0.670005954735      18.231788914567        1759.10038625


*******************************************************************************

THERMODYNAMIC FUNCTIONS WITH VIBRATIONAL CONTRIBUTIONS

AT (T =  298.15 K, P =   0.10132500E+00 MPA):

                         AU/CELL             EV/CELL                 KJ/MOL
ET            :       0.040831798552       1.111089725297         107.20387199
PV            :       0.000017530544       0.000477030356           0.04602644
TS            :       0.080307302905       2.185272809177         210.84679406
ET+PV-TS      :      -0.039457973810      -1.073706053524        -103.59689564
EL+E0+ET+PV-TS:   -2057.056367578673  -55975.349493594607    -5400800.73196695

OTHER THERMODYNAMIC FUNCTIONS:

                     mHARTREE/(CELL*K)     mEV/(CELL*K)              J/(MOL*K)
ENTROPY       :       0.269352013769       7.329440916241         707.18361249
HEAT CAPACITY :       0.253504784422       6.898215882635         665.57671770

*******************************************************************************

Dispersion

Similarly, get thermodynamic properties from phonon dispersion calculations. Properties at q points are summed according to their weights. The output data is written in the output file ‘thermo_paracetamolDisp.txt’.

Note that negative frequencies are all in the translational modes. If more than 3 negative frequencies are found at the same q point, another warning message will be given.

[8]:

from CRYSTALpytools.thermodynamics import Harmonic
import numpy as np

obj = Harmonic(filename='ha_paracetamolDisp.txt', autocalc=False,
               temperature=np.linspace(0, 300, 11), pressure=np.linspace(0, 1, 5))
obj.from_file('ha_paracetamolDisp.out').thermodynamics()

[8]:

<CRYSTALpytools.thermodynamics.Harmonic at 0x7f078f3b9f40>

For more details, please refer to the API documentations.