scnc
/
phonopy
mirror of https://github.com/phonopy/phonopy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Refactoring QHA

This commit is contained in:
Atsushi Togo 2018-07-07 12:01:38 +09:00
parent a2ea5c4605
commit 2befc4b00b
4 changed files with 41 additions and 170 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
doc/qha.rst
View File

@ -18,7 +18,7 @@ follows.
To watch selected plots::
phonopy-qha --sparse=50 -p e-v.dat thermal_properties.yaml-{-{5..1},{0..5}}
phonopy-qha -p e-v.dat thermal_properties.yaml-{-{5..1},{0..5}}
.. figure:: Si-QHA.png
@ -78,24 +78,6 @@ to be lower than the maximum temperature calculated in
``thermal_properties.yaml`` to let at least two temperature points
fewer. The default value is ``--tmax=1000``.
``-p``
~~~~~~~
The fitting results, volume-temperature relation, and thermal expansion
coefficient are plotted on the display.
``--sparse``
~~~~~~~~~~~~~~
This is used with ``-s`` or ``-p`` to thin out the number of plots of
the fitting results at temperatures. When ``--sparse=10``, 1/10 is
only plotted.
``-s``
~~~~~~~
The calculated values are written into files.
``--pressure``
~~~~~~~~~~~~~~~~
@ -109,7 +91,7 @@ bulk modulus obtained with this option than 0 GPa is incorrect.
Fitting volume-energy data to an EOS, and show bulk
modulus (without considering phonons). This is made by::
phonopy-qha -b e-v.dat
% phonopy-qha -b e-v.dat
``--eos``
~~~~~~~~~~~
@ -119,7 +101,26 @@ EOS is chosen among ``vinet``, ``birch_murnaghan``, and
::
phonopy-qha --eos='birch_murnaghan' -b e-v.dat
% phonopy-qha --eos='birch_murnaghan' -b e-v.dat
``-p``
~~~~~~~
The fitting results, volume-temperature relation, and thermal expansion
coefficient are plotted on the display.
``-s``
~~~~~~~
The calculated values are written into files.
``--sparse``
~~~~~~~~~~~~~~
This is used with ``-s`` or ``-p`` to thin out the number of plots of
the fitting results at temperatures. For example with ``--sparse=10``,
1/10 is only plotted.
.. _phonopy_qha_output_files:
@ -133,27 +134,24 @@ are used.
- Bulk modulus :math:`B_T` (GPa) vs :math:`T` (``bulk_modulus-temperature.*``)
- Gibbs free energy :math:`G` (eV) vs :math:`T` (``gibbs-temperature.*``)
- Volume change with respect to the volume at 300 K vs :math:`T`
(``volume_expansion.*``)
- Heat capacity at constant pressure :math:`C_p` (J/K/mol) vs
:math:`T` derived by :math:`-T\frac{\partial^2 G}{\partial T^2}`
(``Cp-temperature.*``)
:math:`T` computed by second order numerical differentiation of
:math:`-T\frac{\partial^2 G}{\partial T^2}` (``Cp-temperature.*``)
- Heat capacity at constant puressure :math:`C_p` (J/K/mol) vs
:math:`T` by polynomial fittings of :math:`C_V` and :math:`S`
(``Cp-temperature_polyfit.*``)
:math:`T` computed by polynomial fittings of :math:`C_V(V)`
(``Cv-volume.dat``) and :math:`S(V)` (``entropy-volume.dat``) for
:math:`\partial S/\partial V` (``dsdv-temperature.dat``) and
numerical differentiation of :math:`\partial V/\partial T`, e.g., see
Eq.(5) of PRB **81**, 17430 by Togo *et al.* (``Cp-temperature_polyfit.*``)
- Volumetric thermal expansion coefficient :math:`\beta` vs :math:`T`
computed by numerical differentiation (``thermal_expansion.*``)
- Volume vs :math:`T` (``volume-temperature.*``)
- Thermodynamics Grüneisen parameter :math:`\gamma = V\beta B_T/C_V`
(no unit) vs :math:`T` (``gruneisen-temperature.dat``)
- Helmholtz free energy (eV) vs volume
(``helmholtz-volume.*``). When ``--pressure`` option is specified,
energy offset of :math:`pV` is added. See also the following section
(:ref:`theory_of_qha`).
- Volume vs :math:`T` (``volume-temperature.*``)
- Volumetric thermal expansion coefficient :math:`\beta` vs :math:`T`
(``thermal_expansion.*``)
- Thermodynamics Grüneisen parameter :math:`\gamma = V\beta B_T/C_V`
(no unit) vs :math:`T` (``gruneisen-temperature.dat``)
``Cv-volume.dat``, ``entropy-volume.dat``,
and ``dsdv-temperature.dat`` (:math:`dS/dV`) are the data internally
used.
.. _theory_of_qha:

4
example/Si-QHA/README
View File

@ -3,10 +3,10 @@ The POSCAR's are the conventional unit cells used to calculate thermal propertie
The following command demonstrates the quasi-harmonic approximation calculation.
% phonopy-qha e-v.dat thermal_properties.yaml-{-{5..1},{0..5}} --sparse=50
% phonopy-qha e-v.dat thermal_properties.yaml-{-{5..1},{0..5}} --sparse=10
For more plots in pdf,
% phonopy-qha -s e-v.dat thermal_properties.yaml-{-{5..1},{0..5}} --sparse=50
% phonopy-qha -s e-v.dat thermal_properties.yaml-{-{5..1},{0..5}} --sparse=10
Before running Si.py, FORCE_SETs have to be created from vasprun.xml's, that are in the compressed file, with disp.yaml that is contained in this directory.

87
phonopy/qha/core.py
View File

@ -113,7 +113,6 @@ class QHA(object):
self._free_energies = None
self._thermal_expansions = None
self._volume_expansions = None
self._cp_numerical = None
self._volume_entropy_parameters = None
self._volume_cv_parameters = None
@ -174,7 +173,6 @@ class QHA(object):
# method is used. Therefore number of temperature points are needed
# than self._max_t_index that nearly equals to the temparature point
# we expect.
self._set_volume_expansion()
self._set_thermal_expansion() # len = len(t) - 1
self._set_heat_capacity_P_numerical() # len = len(t) - 2
self._set_heat_capacity_P_polyfit()
@ -342,58 +340,6 @@ class QHA(object):
self._thermal_expansions[i]))
w.close()
def get_volume_expansion(self):
return self._volume_expansions[:self._max_t_index]
def plot_volume_expansion(self, exp_data=None, symbol='o', _plt=None):
if self._temperatures[self._max_t_index] > 300:
if plt is None:
import matplotlib.pyplot as _plt
self._plot_volume_expansion(_plt,
exp_data=exp_data,
symbol=symbol)
return plt
else:
self._plot_volume_expansion(plt,
exp_data=exp_data,
symbol=symbol)
else:
return None
def plot_pdf_volume_expansion(self,
exp_data=None,
symbol='o',
filename='volume_expansion.pdf'):
if self._temperatures[self._max_t_index] > 300:
import matplotlib.pyplot as plt
plt.rcParams['backend'] = 'PDF'
plt.rcParams['pdf.fonttype'] = 3
plt.rcParams['font.family'] = 'serif'
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['figure.subplot.left'] = 0.15
plt.rcParams['figure.subplot.bottom'] = 0.15
plt.rcParams['figure.figsize'] = 8, 6
fig, ax = plt.subplots()
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_tick_params(which='both', direction='in')
ax.yaxis.set_tick_params(which='both', direction='in')
self._plot_volume_expansion(plt,
exp_data=exp_data,
symbol=symbol)
plt.savefig(filename)
plt.close()
def write_volume_expansion(self, filename='volume_expansion.dat'):
if self._temperatures[self._max_t_index] > 300:
w = open(filename, 'w')
for i in range(self._max_t_index):
w.write("%20.15f %25.15f\n" % (self._temperatures[i],
self._volume_expansions[i]))
w.close()
def get_gibbs_temperature(self):
return self._equiv_energies[:self._max_t_index]
@ -724,27 +670,6 @@ class QHA(object):
plt.xlabel(xlabel)
plt.ylabel(ylabel)
def _plot_volume_expansion(
self,
plt,
exp_data=None,
symbol='o',
xlabel='Temperature (K)',
ylabel=r'Volume expansion $\Delta L/L_0 \, (L=V^{\,1/3})$'):
plt.plot(self._temperatures[:self._max_t_index],
self._volume_expansions[:self._max_t_index],
'r-')
if exp_data:
plt.plot(exp_data[0],
(exp_data[1] / exp_data[1][0]) ** (1.0 / 3) - 1,
symbol)
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.xlim(self._temperatures[0],
self._temperatures[self._max_t_index])
def _plot_gibbs_temperature(self,
plt,
xlabel='Temperature (K)',
@ -831,18 +756,6 @@ class QHA(object):
self._thermal_expansions = beta
def _set_volume_expansion(self):
if self._temperatures[self._max_t_index] > 300:
l = np.array(self._equiv_volumes) ** (1.0 / 3)
for i in range(self._max_t_index):
t = self._temperatures[i]
if (abs(t - 300) <
(self._temperatures[1] - self._temperatures[0]) / 10):
l_0 = (self._equiv_volumes[i]) ** (1.0 / 3)
break
self._volume_expansions = l / l_0 - 1
def _set_heat_capacity_P_numerical(self):
cp = []
g = np.array(self._equiv_energies) * EvTokJmol * 1000

50
scripts/phonopy-qha
View File

@ -39,8 +39,7 @@ import sys
from phonopy.units import *
from phonopy.interface.vasp import read_vasp, write_vasp
from phonopy import PhonopyQHA
from phonopy.file_IO import (read_thermal_properties_yaml, read_v_e,
read_cp, read_ve)
from phonopy.file_IO import read_thermal_properties_yaml, read_v_e
def create_unitcells(unitcell_filename,
@ -74,17 +73,11 @@ def get_options():
is_bulk_modulus_only=False,
eos="vinet",
energy_shift=None,
symbol='o',
cp_filename=None,
ve_filename=None,
thin_number=10,
t_max=1000.0)
parser.add_argument(
"-b", dest="is_bulk_modulus_only", action="store_true",
help="Just show Bulk modulus from v-e data")
parser.add_argument(
"--cp", "--heat_capacity", dest="cp_filename", metavar="FILE",
help="Experimental data for heat capacity Cv")
parser.add_argument(
"--cu", "--create_unitcells", dest="is_create_unitcells",
action="store_true",
@ -115,9 +108,6 @@ def get_options():
parser.add_argument(
"-s", "--save", dest="is_graph_save", action="store_true",
help="Save plot data in pdf")
parser.add_argument(
"--symbol", dest="symbol",
help="Symbol used to plot experiment points")
parser.add_argument(
"--sparse", dest="thin_number", type=int,
help=("Thin out the F-V plots of temperature. The value is "
@ -125,9 +115,6 @@ def get_options():
parser.add_argument(
"--tmax", dest="t_max", type=float,
help="Maximum calculated temperature")
parser.add_argument(
"--ve", "--volume_expansion", dest="ve_filename", metavar="FILE",
help="Experimental data for volume expansion")
parser.add_argument(
"--vf", "--volume_factor", dest="volume_factor", type=float,
help="Conversion factor of volume to A^3")
@ -248,36 +235,18 @@ def main(args):
# - Volume vs Helmholtz free energy
# - Volume vs Temperature
# - Thermal expansion coefficient
if args.ve_filename is None:
exp_data = None
else:
exp_data = parse_ve(args.ve_filename)
phonopy_qha.plot_qha(thin_number=args.thin_number,
volume_temp_exp=exp_data).show()
phonopy_qha.plot_qha(thin_number=args.thin_number).show()
if args.is_graph_save:
# Volume vs Helmholts free energy
phonopy_qha.plot_pdf_helmholtz_volume(thin_number=args.thin_number)
# Volume vs Temperature
if args.ve_filename is None:
exp_data = None
else:
exp_data = parse_ve(args.ve_filename)
phonopy_qha.plot_pdf_volume_temperature(exp_data=exp_data)
phonopy_qha.plot_pdf_volume_temperature()
# Thermal expansion coefficient
phonopy_qha.plot_pdf_thermal_expansion()
# Volume expansion
if args.ve_filename is None:
exp_data = None
else:
exp_data = parse_ve(args.ve_filename)
phonopy_qha.plot_pdf_volume_expansion(exp_data=exp_data,
symbol=args.symbol)
# G vs Temperature
phonopy_qha.plot_pdf_gibbs_temperature()
@ -285,18 +254,10 @@ def main(args):
phonopy_qha.plot_pdf_bulk_modulus_temperature()
# C_P vs Temperature
if args.cp_filename is None:
exp_data = None
else:
exp_data = read_cp(args.cp_filename)
phonopy_qha.plot_pdf_heat_capacity_P_numerical(exp_data=exp_data)
phonopy_qha.plot_pdf_heat_capacity_P_numerical()
# C_P vs Temperature (poly fit)
if args.cp_filename is None:
exp_data = None
else:
exp_data = read_cp(args.cp_filename)
phonopy_qha.plot_pdf_heat_capacity_P_polyfit(exp_data=exp_data)
phonopy_qha.plot_pdf_heat_capacity_P_polyfit()
# Gruneisen parameter vs Temperature
phonopy_qha.plot_pdf_gruneisen_temperature()
@ -304,7 +265,6 @@ def main(args):
phonopy_qha.write_helmholtz_volume()
phonopy_qha.write_volume_temperature()
phonopy_qha.write_thermal_expansion()
phonopy_qha.write_volume_expansion()
phonopy_qha.write_gibbs_temperature()
phonopy_qha.write_bulk_modulus_temperature()
phonopy_qha.write_heat_capacity_P_numerical()