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

rst files were transformed to md files.

This commit is contained in:
Atsushi Togo 2021-07-22 21:17:12 +09:00
parent ec1c657649
commit 4e30cf8e7c
19 changed files with 1 additions and 4132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

410
doc/auxiliary-tools.rst
View File

@ -1,410 +0,0 @@
.. _auxiliary_tools:
Auxiliary tools
===============
.. contents::
:depth: 3
:local:
.. _auxiliary_tools_kaccum:
``phono3py-kaccum``
--------------------
Cumulative physical properties with respect to frequency or mean free
path are calculated using this command.
For example, cumulative thermal conductivity is defined by
.. math::
\kappa^\text{c}(\omega) =
\int^\omega_0 \frac{1}{N} \sum_\lambda
\kappa_\lambda \delta(\omega_\lambda - \omega') d\omega'
:math:`\kappa_\lambda` is the contribution to :math:`\kappa` from the
phonon mode :math:`\lambda`, which is defined as
.. math::
\kappa_\lambda = \frac{1}{V_0}
C_\lambda \mathbf{v}_\lambda \otimes \mathbf{v}_\lambda
\tau_\lambda.
(The notations are found in http://arxiv.org/abs/1501.00691.)
How to use ``phono3py-kaccum``
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Let's computer lattice thermal conductivity of Si using the ``Si-PBEsol``
example found in the example directory.
.. code-block:: shell
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="11 11 11" --sym-fc --br
Then using the output file, ``kappa-m111111.hdf5``, run
``phono3py-kaccum`` as follows:
.. code-block:: shell
% phono3py-kaccum --pa="F" -c POSCAR-unitcell kappa-m111111.hdf5 |tee kaccum.dat
Here ``--pa`` is optional. The definition of ``--pa`` option is same
as :ref:`pa_option`. ``POSCAR-unitcell`` is the unit cell filename
that is specified with ``-c`` option. ``kappa-m111111.hdf5`` is the
output file of thermal conductivity calculation, which is passed to
``phono3py-kaccum`` as the first argument.
The format of the output is as follows: The first column gives
frequency in THz, and the second to seventh columns give the
cumulative lattice thermal conductivity of 6 elements, xx, yy, zz, yz,
xz, xy. The eighth to 13th columns give the derivatives. There are
sets of frequencies, which are separated by blank lines. Each set is
for a temperature. There are the groups corresponding to the number of
temperatures calculated.
To plot the output by gnuplot at temperature index 30 that may
correspond to 300 K,
.. code-block:: shell
% gnuplot
...
gnuplot> p "kaccum.dat" i 30 u 1:2 w l, "kaccum.dat" i 30 u 1:8 w l
The plot like below is displayed.
.. |Si-kaccum| image:: Si-kaccum.png
:width: 50%
|Si-kaccum|
With :math:`19\times 19\times 19` mesh:
.. |Si-kaccum-m191919| image:: Si-kaccum-m191919.png
:width: 25%
|Si-kaccum-m191919|
General options
~~~~~~~~~~~~~~~~
``--pa``
^^^^^^^^^
See :ref:`pa_option`.
``-c``
^^^^^^^
Unit cell filename is specified with this option, e.g., ``-c
POSCAR-unitcell``.
``--qe``
^^^^^^^^
Let ``phono3py-kaccum`` read a QE (pw) unit cell file with ``-c``
option, for example:
.. code-block:: bash
phono3py-kaccum --qe --pa="F" -c Si.in kappa-m191919.hdf5
.. |ipwscf| image:: Si-kaccum-pwscf.png
:width: 25%
|ipwscf|
``--crystal``
^^^^^^^^^^^^^
Analogous to ``--qe``, but to be used with the CRYSTAL interface.
``--turbomole``
^^^^^^^^^^^^^^^
Analogous to ``--qe``, but to be used with the TURBOMOLE interface
``--temperature``
^^^^^^^^^^^^^^^^^^
Pick up one temperature point. For example, ``--temperature=300`` for
300 K, which works only if thermal conductivity is calculated at
temperatures including 300 K.
``--nsp``
^^^^^^^^^^
Number of points to be sampled in the x-axis.
Options for tensor properties
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For cummulative thermal conductivity, the last value is given as the
thermal conductivity in W/mK. For the other properties, the last value
is effectively the sum of values on all mesh grids divided by number
of mesh grids. This is understood as normalized for one primitive
cell. Before version 1.11.13.1, the last value for gv_by_gv (--gv
option) was further divided by the primitive cell volume.
Number of columns of output data is 13 as explained above. With
``--average`` and ``--trace`` options, number of columns of output
data becomes 3.
``--mfp``
^^^^^^^^^^
Mean free path (MFP) is used instead of frequency for the x-axis. MFP
is defined in the single-mode RTA by a vector
.. math::
\mathbf{l}_\lambda = \mathbf{v}_\lambda \tau_\lambda.
The MFP cumulative :math:`\kappa^\text{c}(l)` is given by
.. math::
\kappa^\text{c}(l) =
\int^l_0 \frac{1}{N} \sum_\lambda
\kappa_\lambda \delta(l_\lambda - l') dl'
where :math:`l_\lambda = |\mathbf{l}_\lambda|` and
:math:`\kappa_\lambda` is the contribution to :math:`\kappa` from the
phonon mode :math:`\lambda` in the single-mode RTA, which is defined
as
.. math::
\kappa_\lambda = \frac{1}{V_0} C_\lambda \mathbf{v}_\lambda \otimes
\mathbf{v}_\lambda \tau_\lambda = \frac{1}{V_0} C_\lambda
\mathbf{v}_\lambda \otimes \mathbf{l}_\lambda.
The physical unit of MFP is Angstrom.
The figure below shows the results of Si example with the
:math:`19\times 19\times 19` and :math:`11\times 11\times 11` sampling
meshes used for the lattice thermal conductivity calculation. They look
differently. Especially for the result of the :math:`11\times 11\times
11` sampling mesh, the MFP seems converging but we can see it's not
true to look at that of the :math:`19\times 19\times 19` sampling
mesh. To show this type of plot, be careful about the sampling mesh
convergence.
.. |iMFP| image:: Si-kaccum-MFP.png
:width: 50%
|iMFP|
(This plot is based on the ``Si-PBEsol`` example.)
``--gv``
^^^^^^^^
Outer product of group velocities :math:`\mathbf{v}_\lambda \otimes
\mathbf{v}_\lambda` divided by primitive cell volume (in :math:`\text{THz}^2 /
\text{Angstrom}`)
``--average``
^^^^^^^^^^^^^
Output the traces of the tensors divided by 3 rather than the unique
elements.
``--trace``
^^^^^^^^^^^
Output the traces of the tensors rather than the unique elements.
Options for scalar properties
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For the following properties, those values are normalized by the
number of full grid points. This is understood as normalized for one
primitive cell.
Number of columns of output data is three,
frequency, cumulative property, and derivative of cumulative property
such like DOS.
``--gamma``
^^^^^^^^^^^
:math:`\Gamma_\lambda(\omega_\lambda)` (in THz)
``--tau``
^^^^^^^^^
Lifetime :math:`\tau_\lambda = \frac{1}{2\Gamma_\lambda(\omega_\lambda)}` (in ps)
``--cv``
^^^^^^^^
Modal heat capacity :math:`C_\lambda` (in eV/K)
``--gv-norm``
^^^^^^^^^^^^^
Absolute value of group velocity :math:`|\mathbf{v}_\lambda|` (in
:math:`\text{THz}\cdot\text{Angstrom}`)
``--pqj``
^^^^^^^^^
Averaged phonon-phonon interaction :math:`P_{\mathbf{q}j}` (in :math:`\text{eV}^2`)
``--dos``
^^^^^^^^^
Constant value of 1. This results in phonon DOS.
.. _auxiliary_tools_kdeplot:
``phono3py-kdeplot``
---------------------
**This script is under the development and may contain bugs.** But a
feature is briefly introduced below since it may be useful. Scipy is
needed to use this script.
This script draws density of phonon modes in the frequency-lifetime
plane. Its density is estimated using Gaussian-KDE using `scipy
<https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html>`_.
Then (frequency, lifetime)-data points are superimposed on the density
plot.
``phono3py-kdeplot`` reads a result of the thermal conductivity
calculation as the first argument:
.. code-block:: shell
% phono3py-kdeplot kappa-m111111.hdf5
This calculation takes some time from minutes to hours depending on
mesh numbers and nbins. Therefore it is recommended to start with
smaller mesh and gradually to increase mesh numbers and nbins up to
satisfaction.
After finishing the calculation, the plot is saved in
``lifetime.png``. The black dots show original data points. The
drawing area is automatically set to make the look good, where its
higher lifetime side is not drawn if all density beyond a lifetime
value is smaller than some ratio (see
:ref:`kdeplot_density_ratio`) of the maximum density.
The following plot is drawn with a 19x19x19 mesh and nbins=200 and the
``Si-PBEsol`` example is used to generate the data. The colormap of
'jet' in matplotlib is used.
.. |ikde| image:: Si-kdeplot.png
:width: 50%
|ikde|
Options
~~~~~~~
``--temperature``
^^^^^^^^^^^^^^^^^^
Pick up one temperature point. For example, ``--temperature=300`` for
300 K, which works only if thermal conductivity is calculated at
temperatures including 300 K.
Without specifying this option, the 31st temperature index is
chosen. This often corresponds to 300 K if phono3py ran without
setting temperature range and step.
``--nbins``
^^^^^^^^^^^^
This option controls the resolution of the density plot. The default
value is 100. With larger nbins, the resolution of the plot becomes
better, but the computation will take more time.
.. code-block:: shell
% phono3py-kdeplot --nbins=200 kappa-m111111.hdf5
``--cutoff``, ``--fmax``
^^^^^^^^^^^^^^^^^^^^^^^^^^
The option ``--cutoff`` (``--fmax``) sets the maximum value of
lifetime (frequency) to be included as data points **before**
Gaussian-KDE. Normally increasing this value from the chosen value
without specifying this option does nothing since automatic control of
drawing area cuts high lifetime (frequency) side if the density is low.
``--xmax`` and ``--ymax``
^^^^^^^^^^^^^^^^^^^^^^^^^^
Maximum values of drawing region of phonon frequency (x-axis) and
lifetime (y-axis) are specified by ``--xmax`` and ``--ymax``,
respectively.
``--ymax`` switches off automatic determination of maximum value
of drawing region along y-axis, therefore as a side effect, the
computation will be roughly twice faster.
.. code-block:: shell
% phono3py-kdeplot --ymax=60 kappa-m111111.hdf5
``--zmax``
^^^^^^^^^^^
Maximum value of the density is specified with this option. The color
along colorbar saturates by choosing a smaller value than the maximum value
of density in the data.
.. _kdeplot_density_ratio:
``--dr``, ``--density-ratio``
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The density threshold is specified by the ratio with respect to
maximum density. Normally smaller value results in larger drawing
region. The default value is 0.1. When ``--ymax`` is specified
together, this option is ignored.
.. code-block:: shell
% phono3py-kdeplot --dr=0.01 kappa-m111111.hdf5
``--cmap``
^^^^^^^^^^^
Color map to be used for the density plot. It's given by the name
presented at the matplotlib documentation,
https://matplotlib.org/users/colormaps.html. The default colormap
depends on your matplotlib environment.
.. code-block:: shell
% phono3py-kdeplot --cmap="OrRd" kappa-m111111.hdf5
The following figures are those drawn with ``jet``, ``bwr``,
``seismic``, ``gnuplot``, ``hsv``, and ``OrRd`` colormaps.
.. |ikde-jet| image:: Si-kdeplot-jet.png
:width: 25%
.. |ikde-bwr| image:: Si-kdeplot-bwr.png
:width: 25%
.. |ikde-seismic| image:: Si-kdeplot-seismic.png
:width: 25%
.. |ikde-gnuplot| image:: Si-kdeplot-gnuplot.png
:width: 25%
.. |ikde-hsv| image:: Si-kdeplot-hsv.png
:width: 25%
.. |ikde-OrRd| image:: Si-kdeplot-OrRd.png
:width: 25%
|ikde-jet| |ikde-bwr| |ikde-seismic| |ikde-gnuplot| |ikde-hsv| |ikde-OrRd|

54
doc/citation.rst
View File

@ -1,54 +0,0 @@
How to cite phono3py
====================
Citation of phono3py
--------------------
If you have used phono3py, please cite the following article:
- "Distributions of phonon lifetimes in Brillouin zones",
Atsushi Togo, Laurent Chaput, and Isao Tanaka, Phys. Rev. B, **91**, 094306 (2015)
::
@article{phono3py,
title = {Distributions of phonon lifetimes in Brillouin zones},
author = {Togo, Atsushi and Chaput, Laurent and Tanaka, Isao},
journal = {Phys. Rev. B},
volume = {91},
issue = {9},
pages = {094306},
numpages = {31},
year = {2015},
month = {Mar},
publisher = {American Physical Society},
doi = {10.1103/PhysRevB.91.094306},
}
.. _citation_direct_solution_lbte:
Citation of direct solution of LBTE (``--lbte``)
-------------------------------------------------
If you have used direct solution of LBTE in phono3py, please cite the
following article:
- "Direct Solution to the Linearized Phonon Boltzmann Equation"
Laurent Chaput, Phys. Rev. Lett., **110**, 265506 (2013)
::
@article{PhysRevLett.110.265506,
title = {Direct Solution to the Linearized Phonon Boltzmann Equation},
author = {Chaput, Laurent},
journal = {Phys. Rev. Lett.},
volume = {110},
issue = {26},
pages = {265506},
numpages = {5},
year = {2013},
month = {Jun},
publisher = {American Physical Society},
doi = {10.1103/PhysRevLett.110.265506},
url = {https://link.aps.org/doi/10.1103/PhysRevLett.110.265506}
}

1174
doc/command-options.rst
View File

File diff suppressed because it is too large Load Diff

1
doc/conf.py
View File

@ -29,6 +29,7 @@ import os
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
# extensions = ['sphinx.ext.imgmath']
extensions = ['sphinx.ext.mathjax', 'myst_parser']
myst_enable_extensions = ["linkify", "dollarmath", "amsmath"]
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']

90
doc/crystal.rst
View File

@ -1,90 +0,0 @@
.. _crystal_interface:
CRYSTAL & phono3py calculation
==============================
CRYSTAL program package has a robust built-in phonon calculation
workflow for harmonic phonon properties. However, combining CRYSTAL
with Phono3py enables convenient access to anharmonic phonon properties.
An example for CRYSTAL is found in the ``example/Si-CRYSTAL`` directory.
To invoke the CRYSTAL interface, ``--crystal`` option has to be always
specified::
% phono3py --crystal [options] [arguments]
When the file name of the unit cell is different from the default one
(see :ref:`default_unit_cell_file_name_for_calculator`), ``-c`` option
is used to specify the file name. CRYSTAL unit cell file parser used in
phono3py is the same as that in phonopy. It reads a limited number of
keywords that are documented in the phonopy web site
(http://phonopy.github.io/phonopy/crystal.html#crystal-interface).
.. _crystal_workflow:
Workflow
---------
In this example (Si-CRYSTAL), the CRYSTAL output file is crystal.o.
This is the default file name for the CRYSTAL interface,
so the -c crystal.o parameter is not needed.
1) Create supercells with displacements
(4x4x4 for 2nd order FC, 2x2x2 for 3rd order FC)
::
% phono3py --crystal --dim="2 2 2" --dim-fc2="4 4 4" -d
57 supercell files (``supercell-xxx.d12/.ext``) for the third order
force constants are created. In addition, one supercell file
(``supercell_fc2-00001.d12/.ext``) is created for the second order
force constants.
2) To make valid CRYSTAL input files, there are two possible options:
a) Manually: modify the generated supercell-xxx.d12 and supercell_fc2-xxxxx.d12
files by replacing the line ``Insert basis sets and parameters here`` with the
basis set and computational parameters.
b) Recommended option: before generating the supercells, include files named
``TEMPLATE`` and ``TEMPLATE3`` in the current directory. These files should
contain the basis sets and computational parameters for CRYSTAL (see the example).
When phono3py finds these files it automatically generates complete
CRYSTAL input files in the step 1.
Note that supercells with displacements must not be relaxed in the
force calculations, because atomic forces induced by a small atomic
displacement are what we need for phonon calculation. To get accurate
forces, TOLDEE parameter should be 10 or higher. Phono3py includes this
parameter and the necessary GRADCAL keyword automatically in the inputs.
Then, CRYSTAL supercell calculations are executed to obtain forces on
atoms, e.g., as follows::
% runcry17 supercell-001.d12
3) Collect forces in ``FORCES_FC3`` and ``FORCES_FC2``::
% phono3py --crystal --cf3 supercell-*o
% phono3py --crystal --cf2 supercell_fc2-*o
``disp_fc3.yaml`` and ``disp_fc2.yaml`` are used to create ``FORCES_FC3`` and
``FORCES_FC2``, therefore they must exist in current directory.
4) Calculate 3rd and 2nd order force constants in files ``fc3.hdf5`` and ``fc2.hdf5``::
% phono3py --crystal --dim="2 2 2" --dim-fc2="4 4 4" --sym-fc
``--sym-fc`` is used to symmetrize second- and third-order force constants.
5) Thermal conductivity calculation::
% phono3py --crystal --fc3 --fc2 --dim="2 2 2" --dim-fc2="4 4 4" --mesh="20 20 20" --br
``--br`` invokes the Relaxation Time Approximation.
Add ``--isotope`` for isotope scattering.
Check the effect of ``--nac`` for polar systems.
Carefully test the convergence with respect to ``--mesh``!

443
doc/cutoff-pair.rst
View File

@ -1,443 +0,0 @@
.. _command_cutoff_pair:
Force constants calculation with cutoff pair-distance
=====================================================
Here the detail of the command option :ref:`--cutoff_pair
<cutoff_pair_option>` is explained.
.. contents::
:depth: 2
:local:
What is cutff pair-distance in phono3py?
----------------------------------------
Using ``--cutoff-pair`` option, number of supercells with
displacements to be calculated is reduced. But of course this
sacrifices the accuracy of third-order force constants (fc3).
In phono3py, to obtain supercell-fc3,
:math:`\Phi_{\alpha\beta\gamma}(jl, j'l', j''l'')`, forces in many
supercells having different pairs of displaced atoms are computed
using some force-calculator such as ab-initio code. In the phono3py
default behaviour, full elements of supercell-fc3 are computed. In
this case, though depending on the number of atoms in the supercell
and the crystal symmetry, the number of atomic-pair configuration can
be huge and beyond our computational resource.
Sometimes we may expect that interaction range of fc3 among triplets
of atoms is shorter than chosen supercell size. If it is the case, we
may be allowed to omit computing some elements of supercell-fc3. This
is what achieved by ``--cutoff-pair`` option.
A supercell-fc3 element is specified by three atomic
displacements. Two of three are finitely displaced
(:math:`\Delta\mathbf{r}_1` and :math:`\Delta\mathbf{r}_2`) but one of
them is included in a force given by the force calculator such as
:math:`\mathbf{F}_3=-\frac{\partial V}{\partial\mathbf{r}_3}`. The
cutoff distance :math:`d_\text{cut}` is defined as the upper bound of
the distance between the atoms 1 and 2. By this, the set of atomic
pairs :math:`\{(\mathbf{r}_1,\mathbf{r}_2)| |\mathbf{r}_1 -
\mathbf{r}_2| < d_\text{cut}\}` is selected for the supercell-fc3
calculation. By this, when three distances among the three atoms of
triplets are all larger than :math:`d_\text{cut}`, those fc3 elements
can not be obtained and so they are simply set zero.
Usage
------
Creating supercells with displacements
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
``--cutoff-pair`` option is employed when creating supercells with
displacements, therefore this option must be used with ``-d`` option
when running phono3py, for the Si-PBEsol example:
.. code-block:: bash
% phono3py --cutoff-pair=5 -d --dim="2 2 2" -c POSCAR-unitcell
...
Unit cell was read from "POSCAR-unitcell".
Displacement distance: 0.03
Number of displacements: 111
Cutoff distance for displacements: 5.0
Number of displacement supercell files created: 51
Displacement dataset was written in "phono3py_disp.yaml".
...
% ls POSCAR-0*
POSCAR-00001 POSCAR-00032 POSCAR-00043 POSCAR-00080 POSCAR-00097
POSCAR-00002 POSCAR-00033 POSCAR-00070 POSCAR-00081 POSCAR-00098
POSCAR-00003 POSCAR-00034 POSCAR-00071 POSCAR-00082 POSCAR-00099
POSCAR-00016 POSCAR-00035 POSCAR-00072 POSCAR-00083 POSCAR-00100
POSCAR-00017 POSCAR-00036 POSCAR-00073 POSCAR-00084 POSCAR-00101
POSCAR-00018 POSCAR-00037 POSCAR-00074 POSCAR-00085 POSCAR-00102
POSCAR-00019 POSCAR-00038 POSCAR-00075 POSCAR-00086 POSCAR-00103
POSCAR-00024 POSCAR-00039 POSCAR-00076 POSCAR-00087
POSCAR-00025 POSCAR-00040 POSCAR-00077 POSCAR-00088
POSCAR-00026 POSCAR-00041 POSCAR-00078 POSCAR-00089
POSCAR-00027 POSCAR-00042 POSCAR-00079 POSCAR-00096
% ls POSCAR-0*|wc -l
51
``Number of displacements: 111`` shows the number of supercells with
displacements when this is run without ``--cutoff-pair``
option. ``Number of displacement supercell files created: 51`` gives
the contracted number of supercells with displacements by
``--cutoff-pair`` option. There number of ``POSCAR-0xxxx`` files is found
1. At this step, a special ``phono3py_disp.yaml`` is created. This
contains information on this contraction and used in the other
calculation step, therefore this file must be kept carefully.
Supercell files
^^^^^^^^^^^^^^^
``POSCAR-xxxxx`` (in the other calculator interface, the prefix of the
filename is different) are not generated if distance between a pair of
atoms to be displaced is larger than the specified cutoff pair
distance. The indexing number (``xxxxx``) corresponds to that of the
case without setting this option, i.e., the same ``POSCAR-xxxxx``
files are created for the same configurations of pairs of
displacements but ``POSCAR-xxxxx`` files not being included are not
generated. The reason of this indexing is that it can be useful when
changing the cutoff-pair-distance.
Special ``phono3py_disp.yaml``
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Using ``--cutoff-pair`` option together with ``-d`` option, a special
``phono3py_disp.yaml`` is created. This contains information on distances
between displaced atomic-pairs and whether those pairs are to be
computed or not. This special ``phono3py_disp.yaml`` is necessary to create
fc3, therefore be careful not to overwrite it by running the option
``-d`` without ``--cutoff-pair`` or with different ``--cutoff-pair``
with different value.
Making ``FORCES_FC3``
~~~~~~~~~~~~~~~~~~~~~~
To create ``FORCES_FC3``, only output files of the supercells created
using ``--cutoff-pair`` option are passed to ``phono3py`` as the
arguments. The special ``phono3py_disp.yaml`` file is necessary to be
located at current directory.
An example is shown below for the Si example. Here, it is supposed
that forces are calculated using VASP in ``disp-xxxxx``
directories. After running force calculations, there should be the
output file containing forces in each directory (for VASP
``vasprun.xml``).
.. code-block:: bash
% phono3py --cf3 disp-{00001,00002,00003,00016,00017,00018,00019,00024,00025,00026,00027,00032,00033,00034,00035,00036,00037,00038,00039,00040,00041,00042,00043,00070,00071,00072,00073,00074,00075,00076,00077,00078,00079,00080,00081,00082,00083,00084,00085,00086,00087,00088,00089,00096,00097,00098,00099,00100,00101,00102,00103}/vasprun.xml
...
Displacement dataset is read from phono3py_disp.yaml.
counter (file index): 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
FORCES_FC3 has been created.
Using :ref:`--cf3-file option <cf3_file_option>` may be recommended
when the number of force files is large.
.. code-block:: bash
% for i in `ls POSCAR-0*|sed s/POSCAR-//`;do echo disp-$i/vasprun.xml;done > file_list.dat
% phono3py --cf3-file file_list.dat
Using a python script, ``phono3py_disp.yaml`` is easily parsed. So
it is also easy to create the file list by a python
script:
.. code-block:: python
#!/usr/bin/env python
from phono3py.interface.phono3py_yaml import Phono3pyYaml
p3yaml = Phono3pyYaml()
p3yaml.read("phono3py_disp.yaml")
dds = p3yaml.dataset
file_name_tmpl = "disp-%05d/vasprun.xml"
count = 1
for d1 in dds['first_atoms']:
print(file_name_tmpl % count)
count += 1
for d1 in dds['first_atoms']:
for d2 in d1['second_atoms']:
if d2['included']:
print(file_name_tmpl % count)
count += 1
Running phonon-phonon interaction calculation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
To create fc3, ``--cutoff-pair`` option is not necessary but the
special ``phono3py_disp.yaml`` is required.
.. code-block:: bash
% phono3py
...
----------------------------- General settings -----------------------------
Run mode: None
HDF5 data compression filter: gzip
Crystal structure was read from "phono3py_disp.yaml".
Supercell (dim): [2 2 2]
Primitive matrix:
[0. 0.5 0.5]
[0.5 0. 0.5]
[0.5 0.5 0. ]
Spacegroup: Fd-3m (227)
Use -v option to watch primitive cell, unit cell, and supercell structures.
----------------------------- Force constants ------------------------------
Imposing translational and index exchange symmetry to fc2: False
Imposing translational and index exchange symmetry to fc3: False
Imposing symmetry of index exchange to fc3 in reciprocal space: False
Displacement dataset for fc3 was read from "phono3py_disp.yaml".
Sets of supercell forces were read from "FORCES_FC3".
Computing fc3[ 1, x, x ] using numpy.linalg.pinv with a displacement:
[ 0.0300 0.0000 0.0000]
Expanding fc3.
Cutting-off fc3 (cut-off distance: 5.000000)
Building atom mapping table...
Creating contracted fc3...
Writing fc3 to "fc3.hdf5".
Max drift of fc3: 0.047748 (yxz) 0.047748 (xyz) 0.047748 (xzy)
Displacement dataset for fc2 was read from "phono3py_disp.yaml".
Sets of supercell forces were read from "FORCES_FC3".
Writing fc2 to "fc2.hdf5".
Max drift of fc2: -0.000001 (zz) -0.000001 (zz)
----------- None of ph-ph interaction calculation was performed. -----------
Summary of calculation was written in "phono3py.yaml".
...
Once ``fc3.hdf5`` and ``fc2.hdf5`` are created, ``--cutoff-pair``
option and the special ``phono3py_disp.yaml`` are not needed anymore.
.. code-block:: bash
% phono3py --mesh="11 11 11" --fc3 --fc2 --br
...
300.0 108.051 108.051 108.051 0.000 0.000 -0.000
...
A script extract supercell IDs from ``phono3py_disp.yaml``
----------------------------------------------------------
The following script is an example to collect supercell IDs
with the cutoff-pair distance, for which ``included: true`` is used
to hook them. ``duplicates`` in ``phono3py_disp.yaml`` gives the pairs of
exactly same supercells having different IDs. Therefore one of each
pair is necessary to calculate. As a ratio, the number is not many,
but if we want to save very much the computational demand, it is good
to consider.
.. code-block:: python
#!/usr/bin/env python
from phono3py.interface.phono3py_yaml import Phono3pyYaml
p3yaml = Phono3pyYaml()
p3yaml.read("phono3py_disp.yaml")
data = p3yaml.dataset
disp_ids = []
for data1 in data['first_atoms']:
disp_ids.append(data1['id'])
for data1 in data['first_atoms']:
for data2 in data1['second_atoms']:
if data2['included']:
disp_ids.append(data2['id'])
# To remove duplicates
# duplicates = dict(data['duplicates'])
# disp_ids_nodup = [i for i in disp_ids if i not in duplicates]
print(" ".join(["%05d" % i for i in disp_ids]))
Even for the case that ``phono3py_disp.yaml`` was created without
``--cutoff-pair`` option, if we replace the line in the above script:
.. code-block:: python
if data2['included']:
by
.. code-block:: python
if data2['distance'] < 5.0: # 5.0 is cutoff-pair distance
we can find the supercell IDs almost equivalent to those obtained
above for ``--cutoff-pair 5.0``.
The Test
---------
Si-PBE
~~~~~~
For testing, thermal conductivities with respect to ``--cutoff-pair``
values are calculated as follows. Note that if ``FORCES_FC3`` for full
fc3 elements exists, the same ``FORCES_FC3`` file can be used for
generating contracted fc3 for each special ``phono3py_disp.yaml``.
.. code-block:: bash
% egrep '^\s+pair_distance' phono3py_disp.yaml|awk '{print $2}'|sort|uniq
0.000000
2.366961
3.865232
4.532386
5.466263
5.956722
6.694777
7.100884
7.730463
9.467845
% cp phono3py_disp.yaml phono3py_disp.orig.yaml
% for i in {2..10};do d=`grep pair_distance phono3py_disp.orig.yaml|awk '{print $2}'|sort|uniq|sed "${i}q;d"`; d=$((d+0.1)); phono3py --cutoff-pair=$d -o $i -d --pa="F" --dim="2 2 2" -c POSCAR-unitcell; mv phono3py_disp.yaml phono3py_disp.$i.yaml; done
% ls phono3py_disp.*.yaml
% ls phono3py_disp.*.yaml
phono3py_disp.10.yaml phono3py_disp.5.yaml phono3py_disp.9.yaml
phono3py_disp.2.yaml phono3py_disp.6.yaml phono3py_disp.orig.yaml
phono3py_disp.3.yaml phono3py_disp.7.yaml
phono3py_disp.4.yaml phono3py_disp.8.yaml
% for i in {2..10};do grep number_of_pairs_in_cutoff phono3py_disp.$i.yaml;done
number_of_pairs_in_cutoff: 10
number_of_pairs_in_cutoff: 30
number_of_pairs_in_cutoff: 50
number_of_pairs_in_cutoff: 55
number_of_pairs_in_cutoff: 75
number_of_pairs_in_cutoff: 95
number_of_pairs_in_cutoff: 103
number_of_pairs_in_cutoff: 108
number_of_pairs_in_cutoff: 110
% for i in {2..10};do cp phono3py_disp.$i.yaml phono3py_disp.yaml; phono3py --mesh="11 11 11" --br|tee std.$i.out;done
% for i in {2..10};do egrep '^\s+300' std.$i.out;done
300.0 123.606 123.606 123.606 -0.000 -0.000 0.000
300.0 118.617 118.617 118.617 -0.000 -0.000 0.000
300.0 118.818 118.818 118.818 -0.000 -0.000 0.000
300.0 118.879 118.879 118.879 -0.000 -0.000 0.000
300.0 119.468 119.468 119.468 -0.000 -0.000 0.000
300.0 119.489 119.489 119.489 -0.000 -0.000 0.000
300.0 119.501 119.501 119.501 -0.000 -0.000 0.000
300.0 119.483 119.483 119.483 -0.000 -0.000 0.000
300.0 119.481 119.481 119.481 -0.000 -0.000 0.000
% for i in {2..10};do cp phono3py_disp.$i.yaml phono3py_disp.yaml; phono3py --sym-fc --mesh="11 11 11" --br|tee std.sym-$i.out;done
% for i in {2..10};do egrep '^\s+300' std.sym-$i.out;done
300.0 124.650 124.650 124.650 -0.000 -0.000 0.000
300.0 119.765 119.765 119.765 -0.000 -0.000 0.000
300.0 118.847 118.847 118.847 -0.000 -0.000 0.000
300.0 118.782 118.782 118.782 -0.000 -0.000 0.000
300.0 119.471 119.471 119.471 -0.000 -0.000 0.000
300.0 119.366 119.366 119.366 -0.000 -0.000 0.000
300.0 119.350 119.350 119.350 -0.000 -0.000 0.000
300.0 119.339 119.339 119.339 -0.000 -0.000 0.000
300.0 119.337 119.337 119.337 -0.000 -0.000 0.000
AlN-LDA
~~~~~~~
.. code-block:: bash
% egrep '^\s+pair_distance' phono3py_disp.yaml|awk '{print $2}'|sort|uniq
0.000000
1.889907
1.901086
3.069402
3.076914
3.111000
3.640065
3.645881
4.370303
4.375582
4.743307
4.743308
4.788360
4.978000
5.364501
5.388410
5.672503
5.713938
5.870162
6.205027
6.469591
7.335901
% cp phono3py_disp.yaml phono3py_disp.orig.yaml
% for i in {2..21};do d=`grep pair_distance phono3py_disp.orig.yaml|awk '{print $2}'|sort|uniq|sed "${i}q;d"`; d=$((d+0.0001)); phono3py --cutoff-pair=$d -o $i -d --dim="3 3 2" -c POSCAR-unitcell; mv phono3py_disp.yaml phono3py_disp.$i.yaml; done
% for i in {2..21};do grep number_of_pairs_in_cutoff phono3py_disp.$i.yaml;done
number_of_pairs_in_cutoff: 72
number_of_pairs_in_cutoff: 92
number_of_pairs_in_cutoff: 196
number_of_pairs_in_cutoff: 216
number_of_pairs_in_cutoff: 312
number_of_pairs_in_cutoff: 364
number_of_pairs_in_cutoff: 460
number_of_pairs_in_cutoff: 564
number_of_pairs_in_cutoff: 660
number_of_pairs_in_cutoff: 712
number_of_pairs_in_cutoff: 764
number_of_pairs_in_cutoff: 784
number_of_pairs_in_cutoff: 888
number_of_pairs_in_cutoff: 928
number_of_pairs_in_cutoff: 980
number_of_pairs_in_cutoff: 1020
number_of_pairs_in_cutoff: 1116
number_of_pairs_in_cutoff: 1156
number_of_pairs_in_cutoff: 1208
number_of_pairs_in_cutoff: 1248
% for i in {2..21};do cp phono3py_disp.$i.yaml phono3py_disp.yaml; phono3py --mesh="13 13 9" --br --nac --io $i|tee std.$i.out; done
% for i in {2..21};do egrep '^\s+300\.0' std.$i.out;done
300.0 205.550 205.550 193.665 -0.000 -0.000 -0.000
300.0 218.963 218.963 204.942 -0.000 -0.000 -0.000
300.0 213.624 213.624 193.863 -0.000 -0.000 -0.000
300.0 219.932 219.932 199.819 -0.000 -0.000 -0.000
300.0 235.516 235.516 218.843 -0.000 -0.000 -0.000
300.0 234.750 234.750 217.384 -0.000 -0.000 -0.000
300.0 234.355 234.355 218.030 -0.000 -0.000 -0.000
300.0 235.381 235.381 218.609 -0.000 -0.000 -0.000
300.0 235.996 235.996 219.785 -0.000 -0.000 -0.000
300.0 236.220 236.220 219.867 -0.000 -0.000 -0.000
300.0 236.161 236.161 219.298 -0.000 -0.000 -0.000
300.0 236.096 236.096 219.313 -0.000 -0.000 -0.000
300.0 234.602 234.602 217.064 -0.000 -0.000 -0.000
300.0 235.914 235.914 218.689 -0.000 -0.000 -0.000
300.0 235.049 235.049 217.935 -0.000 -0.000 -0.000
300.0 235.877 235.877 219.065 -0.000 -0.000 -0.000
300.0 236.133 236.133 219.364 -0.000 -0.000 -0.000
300.0 236.207 236.207 219.595 -0.000 -0.000 -0.000
300.0 236.035 236.035 219.463 -0.000 -0.000 -0.000
300.0 236.104 236.104 219.348 -0.000 -0.000 -0.000
% for i in {2..21};do cp phono3py_disp.$i.yaml phono3py_disp.yaml; phono3py --mesh="13 13 9" --br --nac --io $i|tee std.$i.out; done|tee std.sym-$i.out; done
% for i in {2..21};do egrep '^\s+300\.0' std.sym-$i.out;done
300.0 232.964 232.964 216.333 0.000 -0.000 -0.000
300.0 235.442 235.442 219.602 0.000 -0.000 -0.000
300.0 235.521 235.521 217.767 0.000 -0.000 -0.000
300.0 235.581 235.581 217.687 0.000 -0.000 -0.000
300.0 236.837 236.837 219.933 0.000 -0.000 -0.000
300.0 236.020 236.020 219.324 0.000 -0.000 -0.000
300.0 235.482 235.482 218.633 0.000 -0.000 -0.000
300.0 236.313 236.313 219.677 0.000 -0.000 -0.000
300.0 236.308 236.308 219.955 0.000 -0.000 -0.000
300.0 236.074 236.074 219.882 0.000 -0.000 -0.000
300.0 235.520 235.520 219.450 0.000 -0.000 -0.000
300.0 235.769 235.769 219.562 0.000 -0.000 -0.000
300.0 235.441 235.441 219.168 0.000 -0.000 -0.000
300.0 235.892 235.892 219.590 0.000 -0.000 -0.000
300.0 235.509 235.509 219.167 0.000 -0.000 -0.000
300.0 235.646 235.646 219.521 0.000 -0.000 -0.000
300.0 235.783 235.783 219.311 0.000 -0.000 -0.000
300.0 235.887 235.887 219.301 0.000 -0.000 -0.000
300.0 235.642 235.642 219.348 0.000 -0.000 -0.000
300.0 235.728 235.728 219.102 0.000 -0.000 -0.000

363
doc/direct-solution.rst
View File

@ -1,363 +0,0 @@
.. _direct_solution:
Direct solution of linearized phonon Boltzmann equation
========================================================
This page explains how to use the direct solution of LBTE by
`L. Chaput, Phys. Rev. Lett. 110, 265506 (2013)
<https://doi.org/10.1103/PhysRevLett.110.265506>`_ (:ref:`citation <citation_direct_solution_lbte>`).
.. contents::
:depth: 2
:local:
How to use
----------
As written in the following sections, this calculation requires large
memory space. When running multiple temperature points, simply the
memory space needed is multiplied by the number of the temperature
points. Therefore it is normally recommended to specify :ref:`--ts
option <ts_option>`. An example to run with the direct solution of
LBTE for ``example/Si-PBEsol`` is as follows:
.. code-block:: bash
% phono3py --dim="2 2 2" --sym-fc -c POSCAR-unitcell
...
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="11 11 11" --fc3 --fc2 --lbte --ts=300
...
=================== End of collection of collisions ===================
- Averaging collision matrix elements by phonon degeneracy [0.031s]
- Making collision matrix symmetric (built-in) [0.000s]
----------- Thermal conductivity (W/m-k) with tetrahedron method -----------
Diagonalizing by lapacke dsyev... [0.111s]
Calculating pseudo-inv with cutoff=1.0e-08 (np.dot) [0.001s]
# T(K) xx yy zz yz xz xy
300.0 111.588 111.588 111.588 0.000 0.000 -0.000
(RTA) 109.009 109.009 109.009 0.000 0.000 -0.000
----------------------------------------------------------------------------
Thermal conductivity and related properties were written into
"kappa-m111111.hdf5".
Eigenvalues of collision matrix were written into "coleigs-m111111.hdf5"
...
Memory usage
-------------
The direct solution of LBTE needs diagonalization of a large collision
matrix, which requires large memory space. This is the largest
limitation of using this method. The memory size needed for one
collision matrix at a temperature point is :math:`(\text{number of
irreducible grid points} \times \text{number of bands} \times 3)^2`
for the symmetrized collision matrix.
..
and :math:`(\text{number of grid
points} \times \text{number of bands})^2` for the non-symmetrized
collision matrix.
These collision matrices contain real values and are supposed to be
64bit float symmetric matrices. During the diagonalization of each
collision matrix with LAPACK ``dsyev`` solver, around 1.2 times more
memory space is consumed in total.
When phono3py runs with :ref:`--wgp option <wgp_option>` together with
``--lbte`` option, estimated memory space needed for storing collision
matrix is presented. An example for ``example/Si-PBEsol`` is as follows:
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="40 40 40" --wgp --lbte
...
Memory requirements:
- Piece of collision matrix at each grid point, temp and sigma: 0.00 Gb
- Full collision matrix at each temp and sigma: 7.15 Gb
- Phonons: 0.04 Gb
- Grid point information: 0.00 Gb
- Phonon properties: 0.14 Gb
...
With :ref:`--stp option <stp_option>`, estimated
memory space needed for ph-ph interaction strengths is shown.
Work load distribution
-----------------------
The other difficulty compared with RTA is the workload
distribution. Currently there are two ways to distribute the
calculation: (1) Collision matrix is divided and the pieces are
distributed into computing nodes. (2) Ph-ph interaction strengths at
grid points are distributed into computing nodes. These two can not be
mixed, so one of them has to be chosen. In either case, the
distribution is done simply running a set of phono3py calculations
over grid points and optionally band indices. The data computed on
each computing node are stored in an hdf5 file. Increasing the
calculation size, e.g., larger mesh numbers or larger number of atoms
in the primitive cell, large files are created.
.. _distribution_colmat:
Distribution of collision matrix
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A full collision matrix is divided into pieces at grid points of
irreducible part of Brillouin zone. Each piece is calculated
independently from the other pieces. After finishing the calculations
of these pieces, the full collision matrix is diagonzalized to obtain
the thermal conductivity.
File size of Each piece of the collision matrix can be
large. Therefore it is recommended to use :ref:`--ts option
<ts_option>` to limit the number of temperature points, e.g.,
``--ts="100 200 300 400 500``, depending on the memory size installed
on each computing node. To write them into files,
``--write-collision`` option must be specified, and to read them from
files, ``--read-collision`` option is used. These are similarly used
as :ref:`--write-gamma <write_gamma_option>` and :ref:`--read-gamma
<read_gamma_option>` options for RTA calculation as shown in
:ref:`workload_distribution`.
``--read-collision`` option collects the pieces and make one full
collision matrix, then starts to diagonalize it. This option requires
one argument to specify an index to read the collision matrix at one
temperature point, e.g., the collision matrix at 200K is read with
``--read-collision=1`` for the (pieces of) collision matrices created
with ``--ts="100 200 300 400 500"`` (corresponding to 0, 1, 2, 3,
4). The temperature (e.g. 200K) is also read from the file, so it is
unnecessary to specify :ref:`--ts option <ts_option>` when reading.
The summary of the procedure is as follows:
1. Running at each grid point with :ref:`--gp <gp_option>` (or
:ref:`--ga <ga_option>`) option and
saving the piece of the collision matrix to an hdf5 file with
``--write-collision`` option. It is probably OK to calculate and
store the pieces of the collision matrices at multiple temperatures
though it depends on memory size of the computer node. This
calculation has to be done at all irreducible grid points.
2. Collecting and creating all necessary pieces of the collision
matrix with ``--read-collision=num`` (``num``: index of
temperature). By this one full collision matrix at the selected
temperature is created and then diagonalized. An option ``-o num``
may be used together with ``--read-collision`` to distinguish the
file names of the results at different temperatures.
Examples of command options are shown below using ``Si-PBE`` example.
Irreducible grid point indices are obtained by :ref:`--wgp option<wgp_option>`:
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --lbte --wgp
and the information is given in ``ir_grid_points.yaml``. For
distribution of collision matrix calculation
(see also :ref:`workload_distribution`):
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --lbte --ts=300 --write-collision --gp="grid_point_numbers..."
To collect distributed pieces of the collision matrix:
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --lbte --ts=300 --read-collision=0
Distribution of phonon-phonon interaction strengths
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The distribution of pieces of collision matrix is straightforward and
is recommended to use if the number of temperature points is
small. However increasing data file size, network communication
becomes to require long time to send the files from a master node to
computation nodes. In this case, the distribution over ph-ph
interaction strengths can be another choice. Since, without using
:ref:`--full-pp option <full_pp_option>`, the tetrahedron method or
smearing approach with :ref:`--sigma-cutoff option
<sigma_cutoff_option>` results in the sparse ph-ph interaction
strength data array, i.e., most of the elements are zero, the data
size can be reduced by only storing non-zero elements. Not like the
collision matrix, the ph-ph interaction strengths in phono3py are
independent from temperature though it is not the case if the force
constants provided are temperature dependent. Once
stored, they are used to create the collision matrices at
temperatures. Using ``--write-pp`` and ``--read-pp``, they are written
into and read from hdf5 files at grid points.
It is also recommended to use :ref:`--write-phonon option
<write_phonon_option>` and :ref:`--read-phonon option
<read_phonon_option>` to use identical phonon eigenvectors among the
distributed nodes.
The summary of the procedure is as follows:
1. Running at each grid point with :ref:`--gp <gp_option>` (or
:ref:`--ga <ga_option>`) option and storing the ph-ph interaction
strengths to an hdf5 file with ``--write-pp`` option. This calculation
has to be done at all irreducible grid points.
2. Running with ``--read-pp`` option and without :ref:`--gp <gp_option>` (or
:ref:`--ga <ga_option>`) option. By this one full collision matrix at the
selected temperature is created and then diagonalized. An option
``-o num`` may be used together with ``--read-collision`` to
distinguish the file names of the results at different
temperatures.
Examples of command options are shown below using ``Si-PBE`` example.
Irreducible grid point indices are obtained by :ref:`--wgp option<wgp_option>`
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --lbte --wgp
and the grid point information is provided in
``ir_grid_points.yaml``. All phonons on mesh grid points are saved
by
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc2 --write-phonon
For distribution of ph-ph interaction strength calculation (see also
:ref:`workload_distribution`)
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --lbte --ts=300 --write-pp --gp="grid_point_numbers..." --read-phonon
Here one temperature has to be specified but any one of temperatures
is OK since ph-ph interaction strength computed here is assumed to be
temperature independent. Then the computed ph-ph interaction strengths
are read and used to compute collision matrix and lattice thermal
conductivity at a temperature by
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --lbte --ts=300 --read-pp --read-phonon
This last command is repeated at different temperatures to obtain the
properties at multiple temperatures.
.. _diagonzalization_solver:
Cutoff parameter of pseudo inversion
-------------------------------------
To achieve a pseudo inversion, a cutoff parameter is used to find null
space, i.e., to select the nearly zero eigenvalues. The default cutoff
value is ``1e-8``, and this hopefully works in many cases. But if a
collision matrix is numerically not very accurate, we may have to
carefully choose the value by ``--pinv-cutoff`` option. It is safer to
plot the absolute values of eigenvalues in log scale to see if there
is clear gap between non-zero eigenvalue and nearly-zero eigenvalues.
After running the direct solution of LBTE, ``coleigs-mxxx.hdf5`` is
created. This contains the eigenvalues of the collision matrix (either
symmetrized or non-symmetrized). The eigenvalues are plotted using
``phono3py-coleigplot`` in the phono3py package:
.. code-block:: bash
phono3py-coleigplot coleigs-mxxx.hdf5
It is assumed that only one set of eigenvalues at a temperature point
is contained.
.. figure:: Si-coleigplot.png
:width: 50%
:name: coleigplot
Eigenvalues are plotted in log scale (Si-PBEsol exmaple with
15x15x15 mesh). The number in x-axis is just the index where each
eigenvalue is stored. Normally the eigenvalues are stored ascending
order. The bule points show the positive values, and
the red points show the negative values as positive values
(absolute values) to be able to plot in log scale. In this plot, we
can see the gap between :math:`10^{-4}` and :math:`10^{-16}`, which
is a good sign. The values whose absolute values are smaller than
:math:`10^{-8}` are treated as 0 and those solutions are considered
as null spaces.
Installation of diagonalization solvers with multithreaded BLAS
---------------------------------------------------------------
Multithreaded BLAS is recommended to use for the calculation of the
direct solution of LBTE since the diagonalization of the collision
matrix is computationally demanding. A few examples of how to install
multithreded BLAS libraries are presented below.
MKL linked scipy
^^^^^^^^^^^^^^^^^^
Scipy (also numpy) has an interface to LAPACK dsyev
(``scipy.linalg.lapack.dsyev``). An MKL LAPACK linked scipy (also
numpy) gives very good computing performance and is easily obtained
using the anaconda package manager. In this choice, usual installation
of LAPACKE is necessary for running ``dgesvd`` and ``zheev``. When
using anaconda, installing OpenBLAS is the easiest way to do. See
:ref:`install_openblas_lapacke`
OpenBLAS or MKL via LAPACKE
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Using LAPACKE via python C-API is implemented. By this, phono3py can
use LAPACK dsyev. This uses smaller memory space than using MKL linked
scipy. Practically there are two choices, OpenBLAS and MKL. For MKL,
proper installatin of the MKL package is necessary. The MKL
library installed obtained from anaconda can not be used.
OpenBLAS
~~~~~~~~~
Use of OpenBLAS is an easy choice if the anaconda package is used.
See :ref:`install_openblas_lapacke`.
MKL
~~~~
The BLAS multithread performance may be better in that in MKL. Using
MKL-LAPACK via MKL-LAPACKE via python C-API is also implemented if the
link is succeeded. See :ref:`install_mkl_lapacke`.
Solver choice for diagonalization
---------------------------------
For larger systems, diagonalization of collision matrix takes longest
time and requires large memory space. Phono3py relies on LAPACK for
the diagonalization and so the performance is dependent on the choice
of the diagonalization solver.
Using multithreaded BLAS with many-core computing node, computing time
may be well reduced and the calculation can finish in a realistic
time. Currently scipy, numpy and LAPACKE can be used as the LAPACK
wrapper in phono3py. Scipy and numpy distributed by anaconda are MKL
linked, therefore MKL multithread BLAS is used through
them. Multithreaded OpenBLAS is installed by conda and can be used via
LAPACKE in phono3py. MKL LAPACK and BLAS are also able to be used via
LAPACKE in phono3py with appropriate setting in ``setup.py``.
Using ``--pinv-solver=[number]``, one of the following solver is
chosen:
1. Lapacke ``dsyev``: Smaller memory consumption than ``dsyevd``, but
slower. This is the default solver when MKL LAPACKE is integrated or
scipy is not installed.
2. Lapacke ``dsyevd``: Larger memory consumption than ``dsyev``, but
faster. This is not recommended because sometimes a wrong result is
obtained.
3. Numpy's ``dsyevd`` (``linalg.eigh``). This is not recommended
because sometimes a wrong result is obtained.
4. Scipy's ``dsyev``: This is the default solver when scipy is
installed and MKL LAPACKE is not integrated.
5. Scipy's ``dsyevd``. This is not recommended because sometimes a
wrong result is obtained.
The solver choices other than ``--pinv-solver=1`` and
``--pinv-solver=4`` are dangerous and not recommend. They exist just
for the tests.

18
doc/examples.rst
View File

@ -1,18 +0,0 @@
.. _examples_link:
Examples
========
.. contents::
:depth: 2
:local:
Example files are found at
https://github.com/phonopy/phono3py/tree/master/example. The same are
found in the example directory of the phono3py package downloaded at
https://github.com/phonopy/phono3py/archive/master.zip.
Examples of silicon with VASP and Pwscf as calculators are given in
the ``example`` directory. An example using phono3py API is found in
the ``example/Si-PBEsol/Si.py`` but the API document has not yet written.

33
doc/external-tools.rst
View File

@ -1,33 +0,0 @@
.. _external_tools:
External tools
===============
Here external tools related to phono3py but supported by the groups out
of the phono3py project are listed.
Each of the tools is not supported by the phono3py project because of
the difficulties of the maintenance and the test by main developers
of phono3py under current style of phono3py development. However
useful tools should be known. If developers want to use here to notify
their tools, please contact via the phonopy mailing list.
Phono3py-Power-Tools
---------------------
Phono3py-Power-Tools is a collection of codes and scripts from the Skelton group for Phono(3)py "power users", and includes additional tools for data analysis and plotting.
https://github.com/skelton-group/Phono3py-Power-Tools
Kinetic Collective Model
-------------------------
The Kinetic Collective Model (KCM), developed by the nanoTransport
group of the Universitat Autònoma de Barcelona, is a generalization of
the Guyer and Krumhansl solution of the phonon Boltzmann Transport
Equation. KCM allows to compute the thermal conductivity of
semiconductors in a fast and low memory way from first principles
calculations.
https://physta.github.io/

357
doc/hdf5_howto.rst
View File

@ -1,357 +0,0 @@
.. _hdf5_howto:
How to read the results stored in hdf5 files
=============================================
.. contents::
:depth: 3
:local:
How to use HDF5 python library
-------------------------------
It is assumed that ``python-h5py`` is installed on the computer you
interactively use. In the following, how to see the contents of
``.hdf5`` files in the interactive mode of Python. The basic usage of
reading ``.hdf5`` files using ``h5py`` is found at `here
<http://docs.h5py.org/en/latest/high/dataset.html#reading-writing-data>`_.
Usually for running interactive python, ``ipython`` is recommended to
use but not the plain python. In the following example, an MgO result
of thermal conductivity calculation is loaded and thermal conductivity
tensor at 300 K is watched.
::
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5")
In [3]: list(f)
Out[3]:
['frequency',
'gamma',
'group_velocity',
'gv_by_gv',
'heat_capacity',
'kappa',
'kappa_unit_conversion',
'mesh',
'mode_kappa',
'qpoint',
'temperature',
'weight']
In [4]: f['kappa'].shape
Out[4]: (101, 6)
In [5]: f['kappa'][:]
Out[5]:
array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 2.11702476e+05, 2.11702476e+05, 2.11702476e+05,
6.64531043e-13, 6.92618921e-13, -1.34727352e-12],
[ 3.85304024e+04, 3.85304024e+04, 3.85304024e+04,
3.52531412e-13, 3.72706406e-13, -7.07290889e-13],
...,
[ 2.95769356e+01, 2.95769356e+01, 2.95769356e+01,
3.01803322e-16, 3.21661793e-16, -6.05271364e-16],
[ 2.92709650e+01, 2.92709650e+01, 2.92709650e+01,
2.98674274e-16, 3.18330655e-16, -5.98999091e-16],
[ 2.89713297e+01, 2.89713297e+01, 2.89713297e+01,
2.95610215e-16, 3.15068595e-16, -5.92857003e-16]])
In [6]: f['temperature'][:]
Out[6]:
array([ 0., 10., 20., 30., 40., 50., 60., 70.,
80., 90., 100., 110., 120., 130., 140., 150.,
160., 170., 180., 190., 200., 210., 220., 230.,
240., 250., 260., 270., 280., 290., 300., 310.,
320., 330., 340., 350., 360., 370., 380., 390.,
400., 410., 420., 430., 440., 450., 460., 470.,
480., 490., 500., 510., 520., 530., 540., 550.,
560., 570., 580., 590., 600., 610., 620., 630.,
640., 650., 660., 670., 680., 690., 700., 710.,
720., 730., 740., 750., 760., 770., 780., 790.,
800., 810., 820., 830., 840., 850., 860., 870.,
880., 890., 900., 910., 920., 930., 940., 950.,
960., 970., 980., 990., 1000.])
In [7]: f['kappa'][30]
Out[7]:
array([ 1.09089896e+02, 1.09089896e+02, 1.09089896e+02,
1.12480528e-15, 1.19318349e-15, -2.25126057e-15])
In [8]: f['mode_kappa'][30, :, :, :].sum(axis=0).sum(axis=0) / weight.sum()
Out[8]:
array([ 1.09089896e+02, 1.09089896e+02, 1.09089896e+02,
1.12480528e-15, 1.19318349e-15, -2.25126057e-15])
In [9]: g = f['gamma'][30]
In [10]: import numpy as np
In [11]: g = np.where(g > 0, g, -1)
In [12]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
.. _kappa_hdf5_file:
Details of ``kappa-*.hdf5`` file
---------------------------------
Files name, e.g. ``kappa-m323220.hdf5``, is determined by some
specific options. ``mxxx``, show the numbers of sampling
mesh. ``sxxx`` and ``gxxx`` appear optionally. ``sxxx`` gives the
smearing width in the smearing method for Brillouin zone integration
for phonon lifetime, and ``gxxx`` denotes the grid number. Using the
command option of ``-o``, the file name can be modified slightly. For
example ``-o nac`` gives ``kappa-m323220.nac.hdf5`` to
memorize the option ``--nac`` was used.
Currently ``kappa-*.hdf5`` file (not for the specific grid points)
contains the properties shown below.
mesh
^^^^^
(Versions 1.10.11 or later)
The numbers of mesh points for reciprocal space sampling along
reciprocal axes, :math:`a^*, b^*, c^*`
frequency
^^^^^^^^^^
Phonon frequencies. The physical unit is THz, where THz
is in the ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band).
.. _kappa_hdf5_file_gamma:
gamma
^^^^^^
Imaginary part of self energy. The physical unit is THz, where THz
is in the ordinal frequency not the angular frequency.
The array shape for all grid-points (irreducible q-points) is
(temperature, irreducible q-point, phonon band).
The array shape for a specific grid-point is
(temperature, phonon band).
Phonon lifetime (:math:`\tau_\lambda=1/2\Gamma_\lambda(\omega_\lambda)`) may
be estimated from ``gamma``. :math:`2\pi` has to be multiplied with
``gamma`` values in the hdf5 file to convert the unit of ordinal
frequency to angular frequency. Zeros in ``gamma`` values mean that
those elements were not calculated such as for three acoustic modes at
:math:`\Gamma` point. The below is the copy-and-paste from the
previous section to show how to obtain phonon lifetime in pico
second::
In [8]: g = f['gamma'][30]
In [9]: import numpy as np
In [10]: g = np.where(g > 0, g, -1)
In [11]: lifetime = np.where(g > 0, 1.0 / (2 * 2 * np.pi * g), 0)
gamma_isotope
^^^^^^^^^^^^^^
Isotope scattering of :math:`1/2\tau^\mathrm{iso}_\lambda`.
The physical unit is same as that of gamma.
The array shape is same as that of frequency.
group_velocity
^^^^^^^^^^^^^^^
Phonon group velocity, :math:`\nabla_\mathbf{q}\omega_\lambda`. The
physical unit is :math:`\text{THz}\cdot\text{Angstrom}`, where THz
is in the ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band, 3 = Cartesian coordinates).
heat_capacity
^^^^^^^^^^^^^^
Mode-heat-capacity defined by
.. math::
C_\lambda = k_\mathrm{B}
\left(\frac{\hbar\omega_\lambda}{k_\mathrm{B} T} \right)^2
\frac{\exp(\hbar\omega_\lambda/k_\mathrm{B}
T)}{[\exp(\hbar\omega_\lambda/k_\mathrm{B} T)-1]^2}.
The physical unit is eV/K.
The array shape is (temperature, irreducible q-point, phonon band).
.. _output_kappa:
kappa
^^^^^^
Thermal conductivity tensor. The physical unit is W/m-K.
The array shape is (temperature, 6 = (xx, yy, zz, yz, xz, xy)).
.. _output_mode_kappa:
mode-kappa
^^^^^^^^^^^
Thermal conductivity tensors at k-stars (:math:`{}^*\mathbf{k}`):
.. math::
\sum_{\mathbf{q} \in {}^*\mathbf{k}} \kappa_{\mathbf{q}j}.
The sum of this over :math:`{}^*\mathbf{k}` corresponding to
irreducible q-points divided by number of grid points gives
:math:`\kappa` (:ref:`output_kappa`), e.g.,::
kappa_xx_at_index_30 = mode_kappa[30, :, :, 0].sum()/ weight.sum()
Be careful that until version 1.12.7, mode-kappa values were divided
by number of grid points.
The physical unit is W/m-K. Each tensor element is the sum of tensor
elements on the members of :math:`{}^*\mathbf{k}`, i.e., symmetrically
equivalent q-points by crystallographic point group and time reversal
symmetry.
The array shape is (temperature, irreducible q-point, phonon band, 6 =
(xx, yy, zz, yz, xz, xy)).
gv_by_gv
^^^^^^^^^
Outer products of group velocities for k-stars
(:math:`{}^*\mathbf{k}`) for each irreducible q-point and phonon band
(:math:`j`):
.. math::
\sum_{\mathbf{q} \in {}^*\mathbf{k}} \mathbf{v}_{\mathbf{q}j} \otimes
\mathbf{v}_{\mathbf{q}j}.
The physical unit is
:math:`\text{THz}^2\cdot\text{Angstrom}^2`, where THz is in the
ordinal frequency not the angular frequency.
The array shape is (irreducible q-point, phonon band, 6 = (xx, yy, zz,
yz, xz, xy)).
q-point
^^^^^^^^
Irreducible q-points in reduced coordinates.
The array shape is (irreducible q-point, 3 = reduced
coordinates in reciprocal space).
temperature
^^^^^^^^^^^^
Temperatures where thermal conductivities are calculated. The physical
unit is K.
weight
^^^^^^^
Weights corresponding to irreducible q-points. Sum of weights equals to
the number of mesh grid points.
ave_pp
^^^^^^^
Averaged phonon-phonon interaction in :math:`\text{eV}^2`,
:math:`P_{\mathbf{q}j}`:
.. math::
P_{\mathbf{q}j} = \frac{1}{(3n_\mathrm{a})^2} \sum_{\lambda'\lambda''}
|\Phi_{\lambda\lambda'\lambda''}|^2.
This is not going to be calculated in the RTA thermal coductivity
calculation mode by default. To calculate this, ``--full_pp`` option
has to be specified (see :ref:`full_pp_option`).
kappa_unit_conversion
^^^^^^^^^^^^^^^^^^^^^^
This is used to convert the physical unit of lattice thermal
conductivity made of ``heat_capacity``, ``group_velocity``, and
``gamma``, to W/m-K. In the single mode relaxation time (SMRT) method,
a mode contribution to the lattice thermal conductivity is given by
.. math::
\kappa_\lambda = \frac{1}{V_0} C_\lambda \mathbf{v}_\lambda \otimes
\mathbf{v}_\lambda \tau_\lambda^{\mathrm{SMRT}}.
For example, :math:`\kappa_{\lambda,{xx}}` is calculated by::
In [1]: import h5py
In [2]: f = h5py.File("kappa-m111111.hdf5")
In [3]: kappa_unit_conversion = f['kappa_unit_conversion'][()]
In [4]: weight = f['weight'][:]
In [5]: heat_capacity = f['heat_capacity'][:]
In [6]: gv_by_gv = f['gv_by_gv'][:]
In [7]: gamma = f['gamma'][:]
In [8]: kappa_unit_conversion * heat_capacity[30, 2, 0] * gv_by_gv[2, 0] / (2 * gamma[30, 2, 0])
Out[8]:
array([ 1.02050241e+03, 1.02050241e+03, 1.02050241e+03,
4.40486382e-15, 0.00000000e+00, -4.40486382e-15])
In [9]: f['mode_kappa'][30, 2, 0]
Out[9]:
array([ 1.02050201e+03, 1.02050201e+03, 1.02050201e+03,
4.40486209e-15, 0.00000000e+00, -4.40486209e-15])
How to know grid point index number corresponding to grid address
-----------------------------------------------------------------
Runngin with ``--write-gamma``, hdf5 files are written out with file names
such as ``kappa-m202020-g4448.hdf5``. You may want to know the grid point
index number with given grid address. This is done as follows:
.. code-block:: python
In [1]: import phono3py
In [2]: ph3 = phono3py.load("phono3py_disp.yaml")
In [3]: ph3.mesh_numbers = [20, 20, 20]
In [4]: ph3.grid.get_indices_from_addresses([0, 10, 10])
Out[4]: 4448
This index number is different between phono3py version 1.x and 2.x.
To get the number corresponding to the phono3py version 1.x,
``store_dense_gp_map=False`` should be specified in ``phono3py.load``,
.. code-block:: python
In [5]: ph3 = phono3py.load("phono3py_disp.yaml", store_dense_gp_map=False)
In [6]: ph3.mesh_numbers = [20, 20, 20]
In [7]: ph3.grid.get_indices_from_addresses([0, 10, 10])
Out[7]: 4200

93
doc/index.rst
View File

@ -1,93 +0,0 @@
=======================================================
Welcome to phono3py
=======================================================
This software calculates phonon-phonon interaction and related
properties using the supercell approach. For example, the following
physical properties are obtained:
- Lattice thermal conductivity (RTA and :ref:`direct solution of LBTE
<direct_solution>`)
- Phonon lifetime/linewidth
- Imaginary part of self energy
- Joint density of states (JDOS) and weighted-JDOS
Some papers that may introduce phono3py well:
- Theoretical background is summarized in this paper:
http://dx.doi.org/10.1103/PhysRevB.91.094306 (arxiv
http://arxiv.org/abs/1501.00691).
- Introduction to phono3py application:
https://doi.org/10.1103/PhysRevB.97.224306 (open access).
The source code is found at https://github.com/phonopy/phono3py
(BSD-3-Clause). The code is written in Python extended with C and
written as:
- Works at least on Linux systems and MacOS easily.
- Each calculation is distributed over CPU-cores by OpenMP.
- Phonon lifetime (or ph-ph collision) calculations of respective
phonon modes can be executed as independent calculations.
- Thermal conductivity calculations are highly efficiently
distributed over nodes (see :ref:`workload_distribution`).
- User interfaces for :ref:`VASP <vasp_interface>`,
:ref:`QE (pw) <qe_interface>`, :ref:`CRYSTAL <crystal_interface>`,
:ref:`TURBOMOLE <turbomole_interface>`, and Abinit
are built in (see :ref:`calculator_interfaces`).
- API is prepared to operate from Python script (`example
<https://github.com/phonopy/phono3py/blob/master/example/Si-PBEsol/Si.py>`_)
Some tools to analyze the calculated results are prepared (see
:ref:`auxiliary_tools`).
.. |i0| image:: Si-kaccum.png
:width: 20%
.. |iMFP| image:: Si-kaccum-MFP.png
:width: 20%
.. |ikde| image:: Si-kdeplot.png
:width: 22%
|i0| |iMFP| |ikde|
Documentation
=============
.. toctree::
:maxdepth: 1
install
examples
Interfaces to calculators (VASP, QE, CRYSTAL, Abinit, TURBOMOLE) <interfaces>
command-options
output-files
hdf5_howto
auxiliary-tools
direct-solution
workload-distribution
cutoff-pair
external-tools
tips
citation
changelog
Mailing list
============
For questions, bug reports, and comments, please visit following
mailing list:
https://lists.sourceforge.net/lists/listinfo/phonopy-users
Message body including attached files has to be smaller than 300 KB.
License
=======
BSD-3-Clause (New BSD)
Contact
=======
* Author: `Atsushi Togo <http://atztogo.github.io/>`_

323
doc/install.rst
View File

@ -1,323 +0,0 @@
.. _install:
Installation
=============
The detailed installation processes for different environments
are described below. The easiest installation with a good computation
performance is achieved by using the phono3py conda package (see
:ref:`install_an_example`).
.. contents::
:depth: 3
:local:
Installation of phonopy before the installation of phono3py is
required. See how to install phonopy at
https://phonopy.github.io/phonopy/install.html. Phono3py relies on
phonopy, so please use the latest release of phonopy when installing
phono3py.
Installation using conda
-----------------------------
Using conda is the easiest way for installation of phono3py if you are
using x86-64 linux system or macOS. These packages are made and
maintained by Jan Janssen. The installation is simply done by::
% conda install -c conda-forge phono3py
All dependent packages should be installed.
Installation using pip is not recommended
-----------------------------------------
PyPI packages are prepared for phonopy and phono3py releases. When
installing with PyPI, ``setup.py`` is executed locally to compile the
part of the code written in C, so a few libraries such as
lapacke must exist in the system. Those necessary libraries are
explained in the next section.
Installation from source code
------------------------------
When installing phono3py using ``setup.py`` from the source code, a
few libraries are required before running ``setup.py`` script.
For phono3py, OpenMP library is necessary for the multithreding
support. In additon, BLAS, LAPACK, and LAPACKE are also needed. These
packages are probably installed using the package manager for each OS
or conda.
When using gcc to compile phono3py, ``libgomp1`` is necessary to
enable OpenMP multithreading support. This library is probably
installed already in your system. If you don't have it and you use
Ubuntu linux, it is installed by::
% sudo apt-get install libgomp1
.. _install_lapacke:
Installation of LAPACKE
~~~~~~~~~~~~~~~~~~~~~~~~
LAPACK library is used in a few parts of the code to diagonalize
matrices. LAPACK*E* is the C-wrapper of LAPACK and LAPACK relies on
BLAS. Both single-thread or multithread BLAS can be
used in phono3py. In the following, multiple different ways of
installation of LAPACKE are explained.
.. _install_mkl_lapacke:
MKL LAPACKE (with multithread BLAS)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Phono3py can be compiled with MKL for using LAPACKE. If ``setup.py``
finds the file named ``setup_mkl.py``, the contents of ``setup_mkl.py`` is read
and those are included in the compilation setting. For example, the
following setting prepared as ``setup_mkl.py`` seems working on Ubuntu 16.04
system::
intel_root = "/opt/intel/composer_xe_2015.7.235"
mkl_root = "%s/mkl" % intel_root
compiler_root = "%s/compiler" % intel_root
mkl_extra_link_args_lapacke = ['-L%s/lib/intel64' % mkl_root,
'-lmkl_rt']
mkl_extra_link_args_lapacke += ['-L%s/lib/intel64' % compiler_root,
'-lsvml',
'-liomp5',
'-limf',
'-lpthread']
mkl_include_dirs_lapacke = ["%s/include" % mkl_root]
This setting considers to use ``icc`` but it may be compiled with
``gcc``. With ``gcc``, the compiler related setting shown above (i.e.,
around ``compiler_root``) is unnecessary. To achieve this
installation, not only the MKL library but also the header file are
necessary. The libraries are linked dynamically, so in most of the
cases, ``LD_LIBRARY_PATH`` environment variable has to be correctly
specified to let phono3py find those libraries.
.. _install_openblas_lapacke:
OpenBLAS provided by conda (with multithread BLAS)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The installtion of LAPACKE is easy by conda. It is::
% conda install -c conda-forge openblas libgfortran
The recent change of openblas package provided from anaconda makes to
install nomkl, i.e., numpy and scipy with Intel MKL cannot be used
together with openblas. At this moment, this is avoided to install
openblas from conda-forge channel. If the python libraries are not yet
installed::
% conda install -c conda-forge numpy scipy h5py pyyaml matplotlib
Note that using hdf5 files on NFS mouted file system, you may have to
disable file locking by setting
::
export HDF5_USE_FILE_LOCKING=FALSE
This openblas package contains BLAS, LAPACK, and LAPACKE. When this
``libopenblas`` is linked and the ``else`` statement of the C macro
definition section in ``setup.py`` is executed, the following macro
are activated::
if use_setuptools:
extra_compile_args += ['-DMULTITHREADED_BLAS']
else:
define_macros += [('MULTITHREADED_BLAS', None)]
Libraries or headers are not found at the build by ``setup.py``, the
following setting may be of the help::
extra_link_args_lapacke += ['-lopenblas', '-lgfortran']
include_dirs_lapacke += [
os.path.join(os.environ['CONDA_PREFIX'], 'include'), ]
Netlib LAPACKE provided by Ubuntu package manager (with single-thread BLAS)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In the versions of Ubuntu-12.10 or later, LAPACKE
(http://www.netlib.org/lapack/lapacke.html) can be installed from the
package manager (``liblapacke`` and ``liblapacke-dev``)::
% sudo apt-get install liblapack-dev liblapacke-dev
Compiling Netlib LAPACKE
^^^^^^^^^^^^^^^^^^^^^^^^^
The compilation procedure is found at the LAPACKE web site. After
creating the LAPACKE library, ``liblapacke.a`` (or the dynamic link
library), ``setup.py`` must be properly modified to link it. As an
example, the procedure of compiling LAPACKE is shown below.
::
% tar xvfz lapack-3.6.0.tgz
% cd lapack-3.6.0
% cp make.inc.example make.inc
% make lapackelib
BLAS, LAPACK, and LAPACKE, these all may have to be compiled
with ``-fPIC`` option to use it with python.
Building using setup.py
~~~~~~~~~~~~~~~~~~~~~~~~
If package installation is not possible or you want to compile with
special compiler or special options, phono3py is built using
setup.py. In this case, manual modification of ``setup.py`` may be
needed.
1. Download the latest source code at
https://pypi.python.org/pypi/phono3py
2. and extract it::
% tar xvfz phono3py-1.11.13.39.tar.gz
% cd phono3py-1.11.13.39
The other option is using git to clone the phonopy repository from github::
% git clone https://github.com/phonopy/phono3py.git
% cd phono3py
2. Set up C-libraries for python C-API and python codes. This can be
done as follows:
Run ``setup.py`` script via pip::
% pip install -e .
3. Set :envvar:`$PATH` and :envvar:`$PYTHONPATH`
``PATH`` and ``PYTHONPATH`` are set in the same way as phonopy, see
https://phonopy.github.io/phonopy/install.html#building-using-setup-py.
.. _install_an_example:
Installation instruction of latest development version of phono3py
------------------------------------------------------------------
When using conda, ``PYTHONPATH`` should not be set if possible because
potentially wrong python libraries can be imported.
This installation instruction supposes linux x86-64 environment.
1. Download miniconda package
Miniconda is downloaded at https://conda.io/miniconda.html.
For usual 64-bit Linux system::
% wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh -O ~/miniconda.sh
For macOS::
% wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O ~/miniconda.sh
The installation is made by
::
% bash ~/miniconda.sh -b -p $HOME/miniconda
% export PATH="$HOME/miniconda/bin:$PATH"
The detailed installation instruction is found at
https://conda.io/projects/conda/en/latest/user-guide/install/index.html. To
use the conda-forge channel as the default one, ``~/.condarc`` may be
written following the conda-forge official documentation
(https://conda-forge.org/docs/user/tipsandtricks.html).
2. Initialization of conda and setup of conda environment
::
% conda init <your_shell>
``<your_shell>`` is often ``bash`` but may be something else. It is
important that after running ``conda init``, your shell is needed
to be closed and restarted. See more information by ``conda init
--help``.
Then conda allows to make conda installation isolated by using conda's
virtual environment.
::
% conda create -n phono3py -c conda-forge python=3.8
% conda activate phono3py
Use of this is strongly recommended, otherwise proper settings of
``CONDA_PREFIX``, ``C_INCLUDE_PATH``, and ``LD_LIBRARY_PATH`` will
be necessary.
2. Install necessary conda packages for phono3py
::
% conda install -c conda-forge numpy scipy h5py pyyaml matplotlib-base compilers "libblas=*=*mkl" spglib mkl-include
Note that using hdf5 files on NFS mouted file system, you may have to
disable file locking by setting
::
export HDF5_USE_FILE_LOCKING=FALSE
Install the latest phonopy and phono3py from github sources::
% mkdir dev
% cd dev
% git clone https://github.com/phonopy/phonopy.git
% git clone https://github.com/phonopy/phono3py.git
% cd phonopy
% git checkout develop
% python setup.py build
% pip install -e .
% cd ../phono3py
% git checkout develop
% python setup.py build
% pip install -e .
The conda packages dependency can often change and this recipe may
not work properly. So if you find this instruction doesn't work, it
is very appreciated if letting us know it in the phonopy mailing
list.
Multithreading and its controlling by C macro
----------------------------------------------
Phono3py uses multithreading concurrency in two ways. One is that
written in the code with OpenMP ``parallel for``. The other is
achieved by using multithreaded BLAS. The BLAS multithreading is
depending on which BLAS library is chosen by users and the number of
threads to be used may be controlled by the library's environment
variables (e.g., ``OPENBLAS_NUM_THREADS`` or ``MKL_NUM_THREADS``). In
the phono3py C code, these two are written in a nested way, but of
course the nested use of multiple multithreadings has to be
avoided. The outer loop of the nesting is done by the OpenMP
``parallel for`` code. The inner loop calls LAPACKE functions and then
the LAPACKE functions call the BLAS routines. If both of the inner and
outer multithreadings can be activated, the inner multithreading must
be deactivated at the compilation time. This is achieved by setting
the C macro ``MULTITHREADED_BLAS``, which can be written in
``setup.py``. Deactivating the multithreading of BLAS using the
environment variables is not recommended because it is used in the
non-nested parts of the code and these multithreadings are
unnecessary to be deactivated.
Trouble shooting
-----------------
1. Phonopy version should be the latest to use the latest phono3py.
2. There are other pitfalls, see
https://phonopy.github.io/phonopy/install.html#trouble-shooting.

75
doc/interfaces.rst
View File

@ -1,75 +0,0 @@
.. _calculator_interfaces:
Interfaces to calculators
==========================
Currently the built-in interfaces for VASP, QUANTUM ESPRESSO (QE),
CRYSTAL, Abinit, and TURBOMOLE are prepared. VASP is the default interface and no
special option is necessary to invoke it, but for the other
interfaces, each special option has to be specified, e.g. ``--qe``,
``--crystal``, ``--abinit``, or ``--turbomole``
.. toctree::
:maxdepth: 1
vasp
qe
crystal
turbomole
Calculator specific behaviors
------------------------------
Physical unit
^^^^^^^^^^^^^^
The interfaces for VASP, QE (pw), CRYSTAL, Abinit, and TURBOMOLE are
built in to the phono3py command.
For each calculator, each physical unit system is used. The physical
unit systems used for the calculators are summarized below.
::
| unit-cell FORCES_FC3 disp_fc3.yaml
-----------------------------------------------
VASP | Angstrom eV/Angstrom Angstrom
QE (pw) | au (bohr) Ry/au au
CRYSTAL | Angstrom eV/Angstrom Angstrom
Abinit | au (bohr) eV/Angstrom au
TURBOMOLE | au (bohr) hartree/au au
``FORCES_FC2`` and ``disp_fc2.yaml`` have the same physical units as
``FORCES_FC3`` and ``disp_fc3.yaml``, respectively.
Always (irrespective of calculator interface) the physical units of
2nd and 3rd order force constants that are to be stored in
``fc2.hdf5`` and ``fc3.hdf5`` are :math:`\text{eV}/\text{Angstrom}^2` and
:math:`\text{eV}/\text{Angstrom}^3`, respectively.
.. _default_unit_cell_file_name_for_calculator:
Default unit cell file name
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Default unit cell file names are also changed according to the calculators::
VASP | POSCAR
QE | unitcell.in
CRYSTAL | crystal.o
Abinit | unitcell.in
TURBOMOLE | control
.. _default_displacement_distance_for_calculator:
Default displacement distance created
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Default displacement distances created by ``-d`` option without
``--amplitude`` option are respectively as follows::
VASP | 0.03 Angstrom
QE | 0.06 au (bohr)
CRYSTAL | 0.03 Angstrom
Abinit | 0.06 au (bohr)
TURBOMOLE | 0.06 au (bohr)

145
doc/output-files.rst
View File

@ -1,145 +0,0 @@
.. _output_files:
Output files
============
.. contents::
:depth: 3
:local:
The calculation results are written into files. Mostly the data are
stored in HDF5 format, therefore how to read the data
from HDF5 files is also shown.
Intermediate text files
------------------------
The following files are not compatible with phonopy. But phonopy's
``FORCE_SETS`` file can be created using phono3py command options from
the following files. See the detail at :ref:`file_format_compatibility`.
``phono3py_disp.yaml``
^^^^^^^^^^^^^^^^^^^^^^
This is created with ``-d`` option. See :ref:`create_displacements_option`.
``FORCES_FC3``
^^^^^^^^^^^^^^^
This is created with ``--cf3`` option. See :ref:`cf3_option`.
``FORCES_FC2``
^^^^^^^^^^^^^^^
This is created with ``--cf2`` option. See :ref:`cf2_option` and
:ref:`dim_fc2_option`.
HDF5 files
-------------
``kappa-*.hdf5``
^^^^^^^^^^^^^^^^^
See the detail at :ref:`kappa_hdf5_file`.
.. _fc3_hdf5_file:
``fc3.hdf5``
^^^^^^^^^^^^^
Third order force constants (in real space) are stored in
:math:`\mathrm{eV}/\text{Angstrom}^3`.
In phono3py, this is stored in the numpy array ``dtype='double'`` and
``order='C'`` in the shape of::
(num_atom, num_atom, num_atom, 3, 3, 3)
against :math:`\Phi_{\alpha\beta\gamma}(l\kappa, l'\kappa',
l''\kappa'')`. The first three ``num_atom`` are the atom indices in supercell
corresponding to :math:`l\kappa`, :math:`l'\kappa'`,
:math:`l''\kappa''`, respectively. The last three elements are the Cartesian
coordinates corresponding to :math:`\alpha`, :math:`\beta`,
:math:`\gamma`, respectively.
If you want to import a supercell structure and its fc3, you may
suffer from matching its atom index between the supercell and an
expected unit cell. This may be easily dealt with by letting phono3py
see your supercell as the unit cell (e.g., ``POSCAR``,
``unitcell.in``, etc) and find the unit (primitive) cell using
:ref:`--pa option <pa_option>`. For example, let us assume your
supercell is the 2x2x2 multiples of your unit cell that has no
centring, then your ``--pa`` setting will be ``1/2 0 0 0 1/2 0 0 1/2
0``. If your unit cell is a conventional unit cell and has a centring,
e.g., the face centring,
.. math::
(\mathbf{a}_\text{p}, \mathbf{b}_\text{p}, \mathbf{c}_\text{p}) =
(\mathbf{a}_\text{s}, \mathbf{b}_\text{s}, \mathbf{c}_\text{s})
\begin{pmatrix}
\frac{{1}}{2} & 0 & 0 \\
0 & \frac{{1}}{2} & 0 \\
0 & 0 & \frac{{1}}{2}
\end{pmatrix}
\begin{pmatrix}
0 & \frac{{1}}{2} & \frac{{1}}{2} \\
\frac{{1}}{2} & 0 & \frac{{1}}{2} \\
\frac{{1}}{2} & \frac{{1}}{2} & 0
\end{pmatrix} =
(\mathbf{a}_\text{s}, \mathbf{b}_\text{s}, \mathbf{c}_\text{s})
\begin{pmatrix}
0 & \frac{{1}}{4} & \frac{{1}}{4} \\
\frac{{1}}{4} & 0 & \frac{{1}}{4} \\
\frac{{1}}{4} & \frac{{1}}{4} & 0
\end{pmatrix}.
So what you have to set is ``--pa="0 1/4 1/4 1/4 0 1/4 1/4 1/4 0"``.
.. _fc2_hdf5_file:
``fc2.hdf5``
^^^^^^^^^^^^^
Second order force constants are stored in
:math:`\mathrm{eV}/\text{Angstrom}^2`.
In phono3py, this is stored in the numpy array ``dtype='double'`` and
``order='C'`` in the shape of::
(num_atom, num_atom, 3, 3)
against :math:`\Phi_{\alpha\beta}(l\kappa, l'\kappa')`. More detail is
similar to the case for :ref:`fc3_hdf5_file`.
``gamma-*.hdf5``
^^^^^^^^^^^^^^^^^
Imaginary parts of self energies at harmonic phonon frequencies
(:math:`\Gamma_\lambda(\omega_\lambda)` = half linewidths) are stored in
THz. See :ref:`write_gamma_option`.
``gamma_detail-*.hdf5``
^^^^^^^^^^^^^^^^^^^^^^^^
Q-point triplet contributions to imaginary parts of self energies at
phonon frequencies (half linewidths) are stored in THz. See
:ref:`write_detailed_gamma_option`.
Simple text file
-----------------
``gammas-*.dat``
^^^^^^^^^^^^^^^^^
Imaginary parts of self energies with respect to frequency
:math:`\Gamma_\lambda(\omega)` are stored in THz. See :ref:`ise_option`.
``jdos-*.dat``
^^^^^^^^^^^^^^^
Joint densities of states are stored in Thz. See :ref:`jdos_option`.
``linewidth-*.dat``
^^^^^^^^^^^^^^^^^^^^

83
doc/qe.rst
View File

@ -1,83 +0,0 @@
.. _qe_interface:
Quantum ESPRESSO (pw) & phono3py calculation
============================================
Quantum espresso package itself has a set of the force constants
calculation environment based on DFPT. But the document here explains how
to calculate phonon-phonon interaction and related properties using
phono3py, i.e., using the finite displacement and supercell approach.
An example for QE (pw) is found in the ``example-phono3py/Si-QE`` directory.
Unless a proper ``phono3py_disp.yaml`` containing calculator information,
to invoke the QE (pw) interface, ``--qe`` option has to be specified,
.. code-block:: bash
% phono3py --qe [options] [arguments]
When the file name of the unit cell is different from the default one
(see :ref:`default_unit_cell_file_name_for_calculator`), ``-c`` option
is used to specify the file name. QE (pw) unit cell file parser used in
phono3py is the same as that in phonopy. It can read
only limited number of keywords that are shown in the phonopy web site
(http://phonopy.github.io/phonopy/qe.html#qe-interface).
.. _qe_workflow:
Workflow
---------
1. Create supercells with displacements
.. code-block:: bash
% phono3py --qe -d --dim="2 2 2" --pa="F" -c Si.in
In this example, probably 111 different supercells with
displacements are created. Supercell files (``supercell-xxx.in``)
are created but they contain only the crystal
structures. Calculation setting has to be added before running the
calculation. In this step, the option ``--qe`` is necessary.
2. Run QE (pw) for supercell force calculations
Let's assume that the calculations have been made in ``disp-xxx``
directories with the file names of ``Si-supercell.in``. Then after
finishing those calculations, ``Si-supercell.out`` may be created
in each directory.
3. Collect forces
``FORCES_FC3`` is obtained with ``--cf3`` options collecting the
forces on atoms in QE (pw) calculation results:
.. code-block:: bash
% phono3py --cf3 disp-00001/Si-supercell.out disp-00002/Si-supercell.out ...
or in recent bash or zsh:
.. code-block:: bash
% phono3py --cf3 disp-{00001..00111}/Si-supercell.out
``phono3py_disp.yaml`` is used to create ``FORCES_FC3``, therefore it
must exist in current directory.
4) Calculate 3rd and 2nd order force constants
``fc3.hdf5`` and ``fc2.hdf5`` files are created by:
.. code-block:: bash
% phono3py --sym-fc
where ``--sym-fc`` symmetrizes fc3 and fc2.
5) Calculate lattice thermal conductivity, e.g., by:
.. code-block:: bash
% phono3py --mesh="11 11 11" --fc3 --fc2 --br

125
doc/tips.rst
View File

@ -1,125 +0,0 @@
.. _tips:
Tips
=====
.. contents::
:depth: 2
:local:
.. _brillouinzone_sum:
Brillouin zone summation
-------------------------
Brillouin zone (BZ) summations appear at different two points in
phonon lifetime calculation. First it is used for the Fourier
transform of force constants, and second to obtain imaginary part of
phonon-self-energy. For the summation, usually uniform sampling meshes
are employed. To obtain more accurate result, it is always better to
use denser meshes. But the denser mesh requires more computational
demand.
The second BZ summation contains delta functions. In
phono3py calculation, a linear tetrahedron method (:ref:`thm option
<thm_option>`, default option) and a smearing method (:ref:`sigma
option <sigma_option>`) can be used for this BZ
integration. In most cases, the tetrahedron method is better. Especially in high
thermal conductivity materials, the smearing method results in
underestimation of lattice thermal conductivity.
The figure below shows Si thermal conductivity convergence with
respect to number of mesh points along an axis from n=19 to 65. This
is calculated with RTA and the linear tetrahedron method. Within the
methods and phono3py implementation, it is converging at around n=55,
however this computational demand is not trivial. As far as observing
this result, an extrapolation to :math:`1/n \rightarrow 0` seems not a
good idea, since it gives overestimation in the case of this Si
example. This plot tells that we have to decide how much value is
acceptable as lattice thermal conductivity value. Therefore it
important to describe the number of sampling mesh and method of BZ
integration to let other people reproduce the computational results.
.. |isiconv| image:: Si-convergence.png
:width: 25%
|isiconv|
In case the smearing method is necessary to use, the convergence of
q-point mesh together with smearing width has to be checked
carefully. Since smearing parameter is used to approximate delta
functions, small ``sigma`` value is better to describe the detailed
structure of three-phonon-space, however it requires a denser sampling
mesh to converge the result. To check the convergence with respect to
the ``sigma`` value, multiple sigma values can be set. This can be
computationally efficient, since it is avoided to re-calculate
phonon-phonon interaction strength for different ``sigma`` values in
this case. Convergence with respect to the sampling mesh and smearing
parameter strongly depends on materials. For Si example, a
:math:`20\times 20\times 20` sampling mesh (or 8000 reducible sampling
points) and 0.1 THz smearing value for reciprocal of the volume of an
atom may be a good starting choice. The tetrahedron method requires no
such parameter as the smearing width.
Importance of numerical quality of force constants
---------------------------------------------------
Third-order force constants (fc3) are much weaker to numerical noise
of a force calculator than second-order force constants
(fc2). Therefore supercell force calculations have to be done very
carefully.
Numerical quality of forces given by force calculators is the most
important factor for the numerical quality of lattice thermal
conductivity calculation. We may be able to apply symmetry constraints
to force constants, however even if force constants fulfill those
symmetries, the numerical quality of force constants is not guaranteed
since elements of force constants just suffice the symmetries but most
of those intensities are not constrained.
It is important to use the best possible force calculator in the
possibly best way. The knowledge of the force calculator from the
theory and method to the practical usage is required to obtain
good results of lattice thermal conductivity calculation.
In the following, a few things that may be good to know are
presented.
A practical way to check lattice thermal conductivity result
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Some feeling whether our calculation result is OK or not may be
obtained by comparing lattice thermal conductivities calculated with
and without :ref:`symmetrizations of force constants
<symmetrization_option>`. If they are enough different, e.g., more
than twice different, it is better to re-consider about the force
calculation. In the case of DFT calculations, the choice of input
settings such as k-point sampling mesh, plane-wave energy cutoff, and
exchange-correlational potential, etc, should be reconsidered.
Displacement distance of atoms
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The phono3py default displacement distance is 0.03
:math:`\text{Angstrom}`. In some cases, accurate result may not be obtained
due to the numerical noise of the force calculator. Usually increasing
the displacement distance by the :ref:`amplitude option
<amplitude_option>` reduces the numerical noise, but as its drawback
higher order anharmonicity is involved (renormalized) into fc3 and fc2.
.. _file_format_compatibility:
File format compatibility with phonopy
---------------------------------------
- ``FORCES_FC3`` and ``FORCES_FC2`` are not
compatible with phonopy's ``FORCE_SETS``.
- ``FORCE_SETS`` can be created using :ref:`--cfs <cfs_option>` from
``FORCES_FC3`` and ``phono3py_disp.yaml`` or ``FORCES_FC2`` and
``phono3py_disp.yaml`` (needs to specify ``--dim-fc2``).
- ``FORCES_FC2`` can be created using :ref:`--fs2f2
<fs2f2_option>` from ``FORCE_SETS``.
- ``fc2.hdf5`` can be used in phonopy in the ``hdf5`` mode when it is
renamed to ``force_constants.hdf5``. In the previous combinations of
phonopy and phono3py, depending on the physical unit of force
constants of calculators, the direct compatibility is not guranteed.

80
doc/turbomole.rst
View File

@ -1,80 +0,0 @@
.. _turbomole_interface:
TURBOMOLE & phono3py calculation
================================
The riper module of TURBOMOLE can be used to study periodic structures.
An example for TURBOMOLE is found in the ``example/Si-TURBOMOLE`` directory.
To invoke the TURBOMOLE interface, ``--turbomole`` option has to be always
specified::
% phono3py --turbomole [options] [arguments]
When the file name of the unit cell is different from the default one
(see :ref:`default_unit_cell_file_name_for_calculator`), ``-c`` option
is used to specify the file name. TURBOMOLE unit cell file parser used in
phono3py is the same as that in phonopy. It reads a limited number of
keywords that are documented in the phonopy web site
(http://phonopy.github.io/phonopy/turbomole.html#turbomole-interface).
.. _turbomole_workflow:
Workflow
---------
In the example Si-TURBOMOLE, the TURBOMOLE input file is ``control``.
This is the default file name for the TURBOMOLE interface,
so the ``-c control`` parameter is not needed.
1) Create supercells with displacements (2x2x2 conventional cell for
3rd order FC and 3x3x3 conventional cell for 2nd order FC)
::
% phono3py --turbomole --dim="2 2 2" --dim-fc2="3 3 3" -d
111 supercell directories (``supercell-00xxx``) for the third order
force constants are created. In addition, one supercell directory
(``supercell_fc2-00001``) is created for the second order
force constants.
2) Complete TURBOMOLE inputs need to be prepared manually in the subdirectories.
Note that supercells with displacements must not be relaxed in the
force calculations, because atomic forces induced by a small atomic
displacement are what we need for the phonon calculation. To get accurate
forces, $scfconv should be 10. Phono3py includes this data group automatically
in the ``control`` file. You also need to choose a k-point mesh for the force
calculations. TURBOMOLE data group $riper may need to be adjusted to improve
SCF convergence (see example files in subdirectory supercell-00001 for
further details)
Then, TURBOMOLE supercell calculations are executed to obtain forces on
atoms, e.g., as follows::
% riper > supercell-00001.out
3) Collect forces in ``FORCES_FC3`` and ``FORCES_FC2``::
% phono3py --turbomole --cf3 supercell-*
% phono3py --turbomole --cf2 supercell_fc2-*
``disp_fc3.yaml`` and ``disp_fc2.yaml`` are used to create ``FORCES_FC3`` and
``FORCES_FC2``, therefore they must exist in current directory. The Si-TURBOMOLE
example contains pre-calculated force files.
4) Calculate 3rd and 2nd order force constants in files ``fc3.hdf5`` and ``fc2.hdf5``::
% phono3py --turbomole --dim="2 2 2" --dim-fc2="3 3 3" --sym-fc
``--sym-fc`` is used to symmetrize second- and third-order force constants.
5) Thermal conductivity calculation::
% phono3py --turbomole --primitive-axis="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --fc3 --fc2 --dim="2 2 2" --dim-fc2="3 3 3" --mesh="20 20 20" --br
``--primitive-axis`` is used to get the results for the primitive 2-atom cell
``--br`` invokes the Relaxation Time Approximation.
Carefully test the convergence with respect to ``--mesh``!

117
doc/vasp.rst
View File

@ -1,117 +0,0 @@
.. _vasp_interface:
VASP & phono3py calculation
=============================
.. _vasp_workflow:
Workflow
---------
1. Create POSCARs with displacements
This is the same way as usual phonopy:
.. code-block:: bash
% phono3py -d --dim="2 2 2" --pa="F" -c POSCAR-unitcell
``phono3py_disp.yaml`` and ``POSCAR-xxxxx`` files are created.
If you want to use larger supercell size for
second-order force constants (fc2) calculation than that
for third-order force constants (fc3) calculation:
.. code-block:: bash
% phono3py -d --dim-fc2="4 4 4" --dim="2 2 2" --pa="F" -c POSCAR-unitcell
In this case, ``POSCAR_FC2-xxxxx`` files are also created.
1. Run VASP for supercell force calculations
To calculate forces on atoms in supercells, ``POSCAR-xxxxx`` (and
``POSCAR_FC2-xxxxx`` if they exist) are used as VASP (or any force
calculator) calculations.
It is supposed that each force calculation is executed under the
directory named ``disp-xxxxx`` (and ``disp_fc2-xxxxx``), where
``xxxxx`` is sequential number.
2. Collect ``vasprun.xml``'s
When VASP is used as the force calculator, force sets to calculate
fc3 and fc2 are created as follows.
.. code-block:: bash
% phono3py --cf3 disp-{00001..00755}/vasprun.xml
where 0755 is an example of the index of the last displacement
supercell. To perform this collection, ``phono3py_disp.yaml`` created at
step 1 is required. Then ``FORCES_FC3`` is created.
When you use larger supercell for fc2 calculation:
.. code-block:: bash
% phono3py --cf2 disp_fc2-{00001..00002}/vasprun.xml
``phono3py_displ.yaml`` is necessary in this case and ``FORCES_FC2`` is
created.
3. Create fc2.hdf and fc3.hdf
.. code-block:: bash
% phono3py --sym-fc
``--sym-fc`` symmetrizes fc3 and fc2. ``fc2.hdf5`` and ``fc3.hdf5``
are created from ``FORCES_FC3`` (and
optionally ``FORCES_FC2``) and ``phono3py_disp.yaml``. This step is
not mandatory, but you can avoid calculating fc2 and fc3 at every
run time when reading force constants from these files with
``--fc3`` and ``--fc2`` options.
4. Thermal conductivity calculation
An example of thermal conductivity calculation is:
.. code-block:: bash
% phono3py --mesh="11 11 11" --br
This calculation may take very long time. ``--thm`` invokes a
tetrahedron method for Brillouin zone integration for phonon
lifetime calculation, which is the default option. Instead,
``--sigma`` option can be used with the smearing widths.
In this command, phonon lifetimes at many grid points are
calculated in series. The phonon lifetime calculation at each grid
point can be separately calculated since they
are independent and no communication is necessary at the
computation. The procedure is as follows:
First run the same command with the addition option of ``--wgp``:
.. code-block:: bash
% phono3py --fc3 --fc2 --mesh="11 11 11" --br --wgp
``ir_grid_points.yaml`` is obtained. In this file, irreducible
q-points are shown. Then distribute calculations of phonon
lifetimes on grid points with ``--write-gamma`` option by:
.. code-block:: bash
% phono3py --mesh="11 11 11" --br --write-gamma --gp="[grid ponit(s)]"
After finishing all distributed calculations, run with
``--read-gamma`` option:
.. code-block:: bash
% phono3py --fc3 --fc2 --mesh="11 11 11" --br --read-gamma
Once this calculation runs without problem, separately calculated
hdf5 files on grid points are no more necessary and may be deleted.

149
doc/workload-distribution.rst
View File

@ -1,149 +0,0 @@
.. _workload_distribution:
Workload distribution
======================
Workload of thermal conductivity calculation can be distributed into
computer nodes. The distribution over q-point grid points (:ref:`--gp
option <gp_option>`), over phonon bands (:ref:`--bi option
<bi_option>`), and over both of them are supported. Unless necessary,
the distribution over bands is not recommended since it has some
amount of overhead in the part of Fourier transformation of force
constants. Therefore the distribution over grid-points is explained
below. However since the distribution over bands works quite similarly as
that over q-points, the usage can be easily guessed.
On each computer node, pieces of lattice thermal conductivity
calculation are executed. The resulting data for each grid point are
stored in its ``kappa-mxxx-gx.hdf5`` file on each node by setting
:ref:`--write_gamma option <write_gamma_option>`. Once all data are
obtained, those data are collected by :ref:`--read_gamma option
<read_gamma_option>` and the lattice thermal conductivity is obtained.
.. contents::
:depth: 2
:local:
How to do it
------------
The following example is executed in the ``Si-PBE`` example.
To avoid re-calculating fc3 and fc2, ``fc3.hdf5`` and ``fc2.hdf5`` are
created on a single node:
.. code-block:: bash
% phono3py --dim="2 2 2" --sym-fc -c POSCAR-unitcell
The indices of the irreducible grid-points neccesarry to specify
``--ga`` option are found by :ref:`--wgp option <wgp_option>`
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --wgp
and they are stored in ``ir_grid_points.yaml``.
.. code-block:: bash
% egrep '^- grid_point:' ir_grid_points.yaml|awk '{printf("%d,",$3)}'
0,1,2,3,4,5,6,7,8,9,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,60,61,62,63,64,65,66,67,68,69,70,71,72,73,80,81,82,83,84,85,86,87,88,89,90,91,100,101,102,103,104,105,106,107,108,109,120,121,122,123,124,125,126,127,140,141,142,143,144,145,160,161,162,163,180,181,402,403,404,405,406,407,408,409,422,423,424,425,426,427,428,429,430,431,432,433,434,435,442,443,444,445,446,447,448,449,450,451,452,453,462,463,464,465,466,467,468,469,470,471,482,483,484,485,486,487,488,489,502,503,504,505,506,507,522,523,524,525,542,543,804,805,806,807,808,809,824,825,826,827,828,829,830,831,832,833,844,845,846,847,848,849,850,851,864,865,866,867,868,869,884,885,886,887,904,905,1206,1207,1208,1209,1226,1227,1228,1229,1230,1231,1246,1247,1248,1249,1266,1267,1608,1609,1628,1629,
The calculated data on all the grid points shown above as indices are
necessary to obtain lattice thermal conductivity. To distribute
computational demands into computer nodes, a set of the grid-point
indices are chosen and executed as follows:
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --gp="0,1,2,3,4,5,6,7,8,9,20,21,22,23,24,25" --write-gamma
Then many ``kappa-m191919-gx.hdf5`` files are generated. These file
names should not be altered because in reading the data by phono3py,
those file names are supposed to be so, though there is a little
freedom to arrange those file names, for which see :ref:`-o
<output_filename_option>` and :ref:`-i <input_filename_option>`
options. After completing calculations for all irreducible grid-point
indices, the RTA thermal conductivity is computed by another run in a
short time from the stored data:
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --read-gamma
A convenient script
--------------------
The following short script may be useful to splitting all irreducible
grid-point indices into of reasonable number of sets of grid-point
indices for workload distribution.
.. code-block:: python
#!/usr/bin/env python
import sys
import yaml
if len(sys.argv) > 1:
num = int(sys.argv[1])
else:
num = 1
with open("ir_grid_points.yaml") as f:
data = yaml.load(f)
gps = [gp['grid_point'] for gp in data['ir_grid_points']]
gp_lists = [[] for i in range(num)]
for i, gp in enumerate(gps):
gp_lists[i % num].append(gp)
for gp_set in gp_lists:
print(",".join(["%d" % gp for gp in gp_set]))
Supposed that this script is saved as ``divide_gps.py``,
.. code-block:: bash
% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --wgp
...
% python divide_gps.py 20
0,30,52,82,120,402,434,468,524,844,1206
1,31,53,83,121,403,435,469,525,845,1207
2,32,54,84,122,404,442,470,542,846,1208
3,33,55,85,123,405,443,471,543,847,1209
4,34,60,86,124,406,444,482,804,848,1226
5,35,61,87,125,407,445,483,805,849,1227
6,36,62,88,126,408,446,484,806,850,1228
7,37,63,89,127,409,447,485,807,851,1229
8,40,64,90,140,422,448,486,808,864,1230
9,41,65,91,141,423,449,487,809,865,1231
20,42,66,100,142,424,450,488,824,866,1246
21,43,67,101,143,425,451,489,825,867,1247
22,44,68,102,144,426,452,502,826,868,1248
23,45,69,103,145,427,453,503,827,869,1249
24,46,70,104,160,428,462,504,828,884,1266
25,47,71,105,161,429,463,505,829,885,1267
26,48,72,106,162,430,464,506,830,886,1608
27,49,73,107,163,431,465,507,831,887,1609
28,50,80,108,180,432,466,522,832,904,1628
29,51,81,109,181,433,467,523,833,905,1629
For example distributing into 20 computer nodes using a queueing
system,
.. code-block:: bash
% j=1; for i in `python divide_gps.py 20`;do echo $i; sed -e s/gps/$i/g -e s/num/$j/g job.sh|qsub; j=$((j+1)); done
with ``job.sh`` (here for grid-engine):
.. code-block:: bash
#$ -S /bin/zsh
#$ -cwd
#$ -N phono3py-num
#$ -pe mpi* 16
#$ -e err-phono3py-num.log
#$ -o std-phono3py-num.log
phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --gp="gps" --write-gamma