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=None, autocalc=False,
temperature=[298.15], pressure=[0.])
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 phonon read methods defined in io modules are general-propose methods designed to be called (in most cases) internally. Users are always recommended to use the from_file() method after object instantiation.
[2]:
from CRYSTALpytools.thermodynamics import Harmonic
obj = Harmonic(filename=None, autocalc=False,
temperature=[298.15], pressure=[0.])
obj.from_file('ha_Fe2O3G.out')
print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode))
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: 30
The first non-translational mode frequnecy at Gamma: 5.9206 THz
This method also accepts outputs by ‘QHA’ keyword of CRYSTAL with source='crystal-QHA'. 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',
source='crystal-QHA',
qha_index=3)
print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode))
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
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_Fe2O3Disp.out')
print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode))
print('The first non-translational mode frequnecy at Gamma: {:.4f} THz'.format(
obj.frequency[0, 3]))
Number of q / k points: 27
Number of modes at Gamma: 30
The first non-translational mode frequnecy at Gamma: 5.9211 THz
Output files of phonopy calculations in ‘yaml’ format are accepted with source='phonopy'. Use phono_yaml for frequency files (‘mesh.yaml’, ‘band.yaml’ or ‘qpoints.yaml’). ‘qpoints.yaml’ file does not save structure information, so struc_yaml need to be defined 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.
When using ‘mesh.yaml’, the length unit is unknown, so the default unit of ‘angstrom’ is assumed. A warning messange 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_file('ha_Fe2O3mesh.yaml', source='phonopy')
print('Number of q / k points: {:d}'.format(obj.nqpoint))
print('Number of modes at Gamma: {:d}'.format(obj.nmode))
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:165: UserWarning: Static internal energy not available. Setting it to 0.
phonon = Phonon.from_file(filename, source=source,
/home/huanyu/apps/anaconda3/envs/crystal_py3.9/lib/python3.9/site-packages/CRYSTALpytools/phonons.py:457: UserWarning: Unknown length unit. 'angstrom' is assumed.
pobj = YAML.read(strucfile, phonon=filename)
Number of q / k points: 28
Number of modes at Gamma: 30
The first non-translational mode frequnecy at Gamma: 5.9193 THz
Restart calculation
Information can be retrieved from the dumped YAML file with the classmethod Harmonic.restart(). No modification of parameter (k points, scale factors etc.) is allowed.
Thermodynamic calculations will update the dumped file.
[6]:
from CRYSTALpytools.thermodynamics import Harmonic
obj = Harmonic(filename='ha_Fe2O3G_restart.yaml', autocalc=False,
temperature=[298.15], pressure=[0.])
obj.from_file('ha_Fe2O3G.out')
obj = Harmonic.restart('ha_Fe2O3G_restart.yaml')
obj.thermodynamics()
print('Zero-point Energy: {:.4f} kJ/mol'.format(obj.zp_energy))
print('Gibbs Free Energy: {:.4f} kJ/mol'.format(obj.gibbs[0, 0]))
/home/huanyu/apps/anaconda3/envs/crystal_py3.9/lib/python3.9/site-packages/CRYSTALpytools/thermodynamics.py:182: UserWarning: The existing HA file will be overwritten.
ThermoHA.write(self)
Zero-point Energy: 64.2232 kJ/mol
Gibbs Free Energy: -14456416.9077 kJ/mol
NOTE
The text output in < 2025 versions can be generated by setting use_old=True during initialization. But it cannot be used to restart calculations.
[7]:
from CRYSTALpytools.thermodynamics import Harmonic
obj = Harmonic(filename='ha_Fe2O3G_old.txt', use_old=True,
temperature=[298.15], pressure=[0.])
obj.from_file('ha_Fe2O3G.out')
/tmp/ipykernel_32673/1292865050.py:3: UserWarning: The text output is deprecated. Please use the dumping file in YAML format.
obj = Harmonic(filename='ha_Fe2O3G_old.txt', use_old=True,
[7]:
<CRYSTALpytools.thermodynamics.Harmonic at 0x7fdb565e9d00>
Thermodynamics
Gamma point
Call the thermodynamics() method after reading data. Temperature and pressure can also be specified here if not specified at instantiation stage.
[8]:
from CRYSTALpytools.thermodynamics import Harmonic
obj = Harmonic(filename=None, autocalc=False, temperature=[30, 300],
pressure=[0, 1])
obj.from_file('ha_Fe2O3Disp.out').thermodynamics()
print('30 K, 0 GPa')
print('DFT total energy (EL) = ', obj.u_0, 'KJ/mol')
print('Helmholtz free energy (EL+E0+ET-TS) = ', obj.helmholtz[0], ' KJ/mol')
print('Gibbs free energy (EL+E0+ET+pV-TS) = ', obj.gibbs[0, 0], ' KJ/mol')
print('Zero-point energy (E0) = ', obj.zp_energy, ' KJ/mol')
print('Vibrational contribution to interla energy - E0 (ET) = ', obj.u_vib[0] - obj.zp_energy, ' KJ/mol')
print('Entropy*Temperature (TS) = ', obj.entropy[0] * 30, ' J/mol')
print('Heat capacity = ', obj.c_v[0], ' J/mol*K')
print('Entropy = ', obj.entropy[0], ' J/mol*K')
30 K, 0 GPa
DFT total energy (EL) = -14456467.710887987 KJ/mol
Helmholtz free energy (EL+E0+ET-TS) = -14456400.140414735 KJ/mol
Gibbs free energy (EL+E0+ET+pV-TS) = -14456400.140414737 KJ/mol
Zero-point energy (E0) = 67.57226510101925 KJ/mol
Vibrational contribution to interla energy - E0 (ET) = 0.01025610055198456 KJ/mol
Entropy*Temperature (TS) = 12.047951449420218 J/mol
Heat capacity = 2.038146719204007 J/mol*K
Entropy = 0.4015983816473406 J/mol*K
References from ‘ha_Fe2O3Disp.out’ (30 K, 0 GPa) are attached below.
*******************************************************************************
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 : -5506.177755084074 -149830.713962144044 -14456467.65868747
E0 : 0.025736916490 0.700337102153 67.57226472
*******************************************************************************
THERMODYNAMIC FUNCTIONS WITH VIBRATIONAL CONTRIBUTIONS
AT (T = 30.00 K, P = 0.00000000E+00 MPA):
AU/CELL EV/CELL KJ/MOL
ET : 0.000003906371 0.000106297755 0.01025618
PV : 0.000000000000 0.000000000000 0.00000000
TS : 0.000004588856 0.000124869134 0.01204804
ET+PV-TS : -0.000000682486 -0.000018571378 -0.00179187
EL+E0+ET+PV-TS: -5506.152018850070 -149830.013643613260 -14456400.08821461
OTHER THERMODYNAMIC FUNCTIONS:
mHARTREE/(CELL*K) mEV/(CELL*K) J/(MOL*K)
ENTROPY : 0.000152961883 0.004162304453 0.40160137
HEAT CAPACITY : 0.000776294034 0.021124034580 2.03815970
*******************************************************************************
[9]:
from CRYSTALpytools.thermodynamics import Harmonic
obj = Harmonic(filename=None, autocalc=False, temperature=[30, 300],
pressure=[0, 1])
obj.from_file('ha_Fe2O3Disp.out').thermodynamics()
print('300 K, 1 GPa')
print('Helmholtz free energy (EL+E0+ET-TS) = ', obj.helmholtz[1], ' KJ/mol')
print('Gibbs free energy (EL+E0+ET+pV-TS) = ', obj.gibbs[1, 1], ' KJ/mol')
print('Vibrational contribution to interla energy - E0 (ET) = ', obj.u_vib[1] - obj.zp_energy, ' KJ/mol')
print('Entropy*Temperature (TS) = ', obj.entropy[1] * 300, ' J/mol')
print('Heat capacity = ', obj.c_v[1], ' J/mol*K')
print('Entropy = ', obj.entropy[1], ' J/mol*K')
300 K, 1 GPa
Helmholtz free energy (EL+E0+ET-TS) = -14456419.45244564 KJ/mol
Gibbs free energy (EL+E0+ET+pV-TS) = -14456357.394993309 KJ/mol
Vibrational contribution to interla energy - E0 (ET) = 28.954255811802184 KJ/mol
Entropy*Temperature (TS) = 48268.078570283185 J/mol
Heat capacity = 187.33142960669963 J/mol*K
Entropy = 160.89359523427729 J/mol*K
References from ‘ha_Fe2O3Disp.out’ (300 K, 1 GPa) are attached below.
*******************************************************************************
THERMODYNAMIC FUNCTIONS WITH VIBRATIONAL CONTRIBUTIONS
AT (T = 300.00 K, P = 0.10000000E+04 MPA):
AU/CELL EV/CELL KJ/MOL
ET : 0.011028118621 0.300090363974 28.95432136
PV : 0.023636440017 0.643180231504 62.05746452
TS : 0.018384386708 0.500264595280 48.26820050
ET+PV-TS : 0.016280171930 0.443006000198 42.74358538
EL+E0+ET+PV-TS: -5506.135737995654 -149829.570619041682 -14456357.34283737
OTHER THERMODYNAMIC FUNCTIONS:
mHARTREE/(CELL*K) mEV/(CELL*K) J/(MOL*K)
ENTROPY : 0.061281289026 1.667548650933 160.89400166
HEAT CAPACITY : 0.071350898293 1.941556649380 187.33175707
*******************************************************************************
Here are several findings:
Parameters from frequencies agree well with reference data. The choice of physical constants and unit conversion coefficients might lead to the discrepancy.
The difference in DFT total energy lies in different coefficients for unit conversion. In CRYSTALpytools, constants and coefficients from
`scipy.constants<https://docs.scipy.org/doc/scipy/reference/constants.html>`__ is used.
For more details, please refer to the API documentations.