diff --git a/doc/auxiliary-tools.md b/doc/auxiliary-tools.md new file mode 100644 index 00000000..140dbf02 --- /dev/null +++ b/doc/auxiliary-tools.md @@ -0,0 +1,368 @@ +(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 + +$$ +\kappa^\text{c}(\omega) = + \int^\omega_0 \frac{1}{N} \sum_\lambda +\kappa_\lambda \delta(\omega_\lambda - \omega') d\omega' +$$ + +$\kappa_\lambda$ is the contribution to $\kappa$ from the +phonon mode $\lambda$, which is defined as + +$$ +\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. + +```bash +% 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: + +```bash +% 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, + +```bash +% 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. + +```{image} Si-kaccum.png +:width: 50% +``` + +With $19\times 19\times 19$ mesh: + +```{image} Si-kaccum-m191919.png +:width: 25% +``` + +###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: + +```bash +% phono3py-kaccum --qe --pa="F" -c Si.in kappa-m191919.hdf5 +``` + +```{image} Si-kaccum-pwscf.png +:width: 25% +``` + +#### `--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 + +$$ +\mathbf{l}_\lambda = \mathbf{v}_\lambda \tau_\lambda. +$$ + +The MFP cumulative $\kappa^\text{c}(l)$ is given by + +$$ +\kappa^\text{c}(l) = + \int^l_0 \frac{1}{N} \sum_\lambda +\kappa_\lambda \delta(l_\lambda - l') dl' +$$ + +where $l_\lambda = |\mathbf{l}_\lambda|$ and +$\kappa_\lambda$ is the contribution to $\kappa$ from the +phonon mode $\lambda$ in the single-mode RTA, which is defined +as + +$$ +\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 +$19\times 19\times 19$ and $11\times 11\times 11$ sampling +meshes used for the lattice thermal conductivity calculation. They look +differently. Especially for the result of the $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 $19\times 19\times 19$ sampling +mesh. To show this type of plot, be careful about the sampling mesh +convergence. + +```{image} Si-kaccum-MFP.png +:width: 50% +``` + +(This plot is based on the `Si-PBEsol` example.) + +#### `--gv` + +Outer product of group velocities $\mathbf{v}_\lambda \otimes +\mathbf{v}_\lambda$ divided by primitive cell volume (in $\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` + +$\Gamma_\lambda(\omega_\lambda)$ (in THz). + +#### `--tau` + +Lifetime $\tau_\lambda = \frac{1}{2\Gamma_\lambda(\omega_\lambda)}$ (in ps) + +#### `--cv` + +Modal heat capacity $C_\lambda$ (in eV/K). + +#### `--gv-norm` + +Absolute value of group velocity $|\mathbf{v}_\lambda|$ (in +$\text{THz}\cdot\text{Angstrom}$). + +#### `--pqj` + +Averaged phonon-phonon interaction $P_{\mathbf{q}j}$ (in $\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: + +```bash +% 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 $19 \times 19 \times 19$ mesh and nbins=200 and the +`Si-PBEsol` example is used to generate the data. The colormap of +'jet' in matplotlib is used. + +```{image} Si-kdeplot.png +:width: 50% +``` + +### 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. + +```bash +% 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. + +```bash +% 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. + +```bash +% 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. + +```bash +% phono3py-kdeplot --cmap="OrRd" kappa-m111111.hdf5 +``` + +The following figures are those drawn with `jet`, `bwr`, +`seismic`, `gnuplot`, `hsv`, and `OrRd` colormaps. + + +```{image} Si-kdeplot-jet.png +:width: 25% +``` +```{image} Si-kdeplot-bwr.png +:width: 25% +``` +```{image} Si-kdeplot-seismic.png +:width: 25% +``` +```{image} Si-kdeplot-gnuplot.png +:width: 25% +``` +```{image} Si-kdeplot-hsv.png +:width: 25% +``` +```{image} Si-kdeplot-OrRd.png +:width: 25% +``` diff --git a/doc/changelog.md b/doc/changelog.md new file mode 100644 index 00000000..729ed487 --- /dev/null +++ b/doc/changelog.md @@ -0,0 +1,408 @@ +(changelog)= +# Change Log + +## Jul-22-2021: Version 2.0.0 + +This is a major version release. There are some backword-incompatible changes. + +1. Grid point indexing system to address grid points of q-points + is changed. +2. Most of integer array data type is changed to + `dtype='int_'` from `dtype='intc'`. + +To emulate the version 1.x behaviour in `phono3py` command, +try `--v1` option. To emurate the version 1.x behaviour in API, +specify `store_dense_gp_map=False` +and `store_dense_svecs=False` in instatiation of `Phono3py` class +or phon3py loader. + +## Mar-17-2021: Version 1.22.3 + +- Fix `--read-gamma` option to work. + +## Feb-21-2021: Version 1.22.2 + +- Fix PyPI source distribution package + +## Feb-21-2021: Version 1.22.1 + +- `phono3py` command didn't work. This was fixed. +- Fix behaviour when specifying `--thm` and `--sigma` simultaneously. + +## Jan-29-2021: Version 1.22.0 + +- Maintenance release to follow phonopy v2.9.0. + +## Sep-30-2020: Version 1.21.0 + +- Maintenance release to follow the change of phonopy at v2.8.1 +- Improvements of phono3py loader (`phono3py.load`), `phono3py-load` + command, API, and `phono3py_disp.yaml`. +- Harmonic phonon calculation on mesh was multithreaded. This is + effective when using very dense mesh with non-analytical term + correction (probably rare case). +- Real and imaginary parts of self energy and spectral function of + bubble diagram at API level + +## Mar-3-2020: Version 1.20.0 + +- `phono3py_disp.yaml` is made when creating displacements in + addition to `disp_fc3.yaml` and + `disp_fc2.yaml`. `phono3py_disp.yaml` will be used instead of + `disp_fc3.yaml` and `disp_fc2.yaml` in the future major release + (v2.0). + +## Mar-3-2020: Version 1.19.1 + +- Release for pypi and conda (atztogo channel) packagings + +## Mar-2-2020: Version 1.19.0 + +- Improvements of phono3py loader and API. +- Improvement of interfaces to calculators. Now it is expected to be + much easier to implement calculator interface if it exists in + phonopy. +- Fixed dependency to phonopy v2.6.0. + +## Dec-22-2019: Version 1.18.2 + +- Initial version of phono3py loader (`phono3py.load`) was + implemented. See docstring of `phono3py.load`. + +## Oct-17-2019: Version 1.18.1 + +- Fix of phono3py-kaccum to follow the latest phonopy. + +## Oct-9-2019: Version 1.18.0 + +- Maintenance release + +## Apr-18-2019: Version 1.17.0 + +- `--cfz` option was made to subtract residual forces. See + {ref}`cfz_option`. +- `--cutoff-pair` was made to override the cutoff pair distance + written in `disp_fc3.yaml` when using on calculating force + constants. This is useful when checking cutoff distance + dependency. So the use case of having fully computed `FORCES_FC3` + is assumed. +- TURBOMOLE interface is provided by Antti Karttunen + (`--turbomole`). +- Compatibility of `fc2.hdf5` and `force_constants.hdf5` was + improved for all calculators to store physical unit information in + the hdf5 file. See {ref}`file_format_compatibility`. + +## Mar-24-2019: Version 1.16.0 + +- Bug fixes and catching up the updates of phonopy. +- Most of hdf5 output files are compressed by `gzip` as + default. This compression can be set off, see + {ref}`hdf5_compression_option`. +- (Experimental) `phono3py` command accepts `phono3py.yaml` type + file as an input crystal structure by `-c` option. When `DIM` + and any structure file are not given, `phono3py_disp.yaml` + (primary) or `phono3py.yaml` (secondary) is searched in the current + directory. Then `phono3py.yaml` type file is used as the input. + By this, semi-automatic phono3py mode is invocked, which acts as + + 1. `supercell_matrix` corresponding to `DIM` in the + `phono3py.yaml` type file is used if it exists. + 2. `phonon_supercell_matrix` corresponding to `DIM_FC2` in the + `phono3py.yaml` type file is used if it exists. + 3. `primitive_matrix` in the `phono3py.yaml` type file + is used if it exists. Otherwise, set `PRIMITIVE_AXES = AUTO` + when `PRIMITIVE_AXES` is not given. + 4. NAC params are read (`NAC = .TRUE.`) if NAC params are + contained (primary) in the `phono3py.yaml` type file or if + `BORN` file exists in the current directory (secondary). + +## Nov-22-2018: version 1.14.3 + +- Update to work with phonopy v1.14.2. +- Ph-ph interaction can be read (`--read-pp`) and write + (`--write-pp`) in RTA thermal conductivity calculation, too. Mind + that the data stored are different with and without + `--full-pp`. Wihtout `--full-pp` the data are stored in + complicated way to save data side, so it is not considered readable + by usual users. + +## June-20-2018: version 1.13.3 + +- `--lw` (linewidth) option was removed. Use `--br` option and + find 2*gamma values as linewidths in `kappa-xxx.hdf5` file. +- Documentation of `--lbte` option is available at + {ref}`direct_solution`. +- This version is dependent on phonopy>=1.13.2. + +## May-17-2018: version 1.13.1 + +- Compact force constants are implemented (See {ref}`compact_fc_option`). + +## Mar-16-2018: version 1.12.9 + +- Definition of `mode_kappa` values in output hdf5 file is + changed. Previously they were divided by number of grid points, but + now not. Therefore users who compute `kappa` from `mode_kappa` + need to be careful about this change. This does not affect to + `phono3py-kaccum` results. + +## Feb-1-2018: version 1.12.7 + +- `--tsym` option is removed. Now with `--sym-fc3r` and + `--sym-fc2` options, + translational invariance symmetry is also applied. +- `--sym-fc` option is added. This is just an alias to specify both + `--sym-fc3r` and `--sym-fc2` together. +- Documentation on `--write-phonon` and `--read-phonon` options is + written. These options are used to save harmonic phonon infromation + on strage. + +## Nov-22-2017: version 1.12.5 + +- Bug fix of RTA thermal conductivity. This bug exists from version + 1.10.11.18 (git e40cd059). This bug exhibits when all the following + conditions are met: + + 1. RTA thermal conductivity calculation. + 2. Tetrahedron method for Brillouin zone integration is used. + 3. Number of triplets is smaller than number of bands at each grid point. + 4. Without using `--full-pp`. + + + (3) happens when the primitive cell is relatively large. Number of + triplets can be shown using `--stp` option. A race condition of + OpenMP multithreading is the source of the bug. Therefore, if it + occurs, the same calculation comes up with the different thermal + conductivity value in every run time, for whcih it behaves like + randomly. + +- RTA thermal conductivity with smearing method (`--sigma`) is made + to run with smaller memory consumption as similar as tetrahedron + method (`--thm`). + +## Nov-17-2017: version 1.12.3 + +- Command option parser of the phonopy tools is replaced from + `optparse` to `argparse`. +- The filenames used with these options were the positional arguments + previously. Now they are the command-line arguments, i.e., filenames + have to be put just after the option name like `-f vasprun.xml-001 + vasprun.xml-002 ...`. +- The names of auxiliary tools (`kdeplot` and `kaccum`) are + changed, for which the prefix phono3py- is added to the old names to + avoid accidental conflict with other script names already existing + under bin directory. +- {ref}`sigma_cutoff_option` option was created. + +## Jun-18-2017: version 1.11.13 + +- {ref}`normal_umklapp_option` option was made. +- Many minor updates: fixing bugs, improving usabilities. +- Improve of {ref}`auxiliary_tools_kaccum` and {ref}`auxiliary_tools_kdeplot`. + +## Mar-31-2017: version 1.11.11 + +- Abinit code interface is implemented and now under the testing. +- Reduction of memory usage in RTA thermal conductivity + calculation. This is especially effective for larger unit cell + case. Currently combinations with --full_pp, --write_gamma_detail, + and --simga(smearing method) are not supported for this. Performance + tuning is under going. In some case, computation can be slower than + the previous versions. + +## Feb-9-2017: version 1.11.9 + +- This version works coupled with phonopy-1.11.8 or later. +- CRYSTAL code interface is implemented by Antti J. Karttunen. + +## Dec-14-2016: version 1.11.7 + +- This is a maintenance release. This version must be used with + phonopy-1.11.6 or later. + +## Nov-27-2016: version 1.11.5 + +- `gaccum` is merged to `kaccum`. `gaccum` is removed. See + {ref}`auxiliary_tools_kaccum`. +- `kdeplot` is added. See {ref}`auxiliary_tools_kdeplot`. + +## Apr-24-2016: version 1.10.9 + +- Failure of writing `kappa-mxxx-gx.hdf5` was fixed. + +## Apr-16-2016: version 1.10.7 + +- API example is prepared and it is found in `Si` example. No + doucment yet. +- Si pwscf example was placed in `example-phono3py` directory. +- User interface bug fix. + +## Mar-15-2016: version 1.10.5 + +- Numbering way of phono3py version was just changed (No big updates + were made against previous version.) The number is given based on + the phonopy version. For example, the harmonic part of + phono3py-1.10.5 is based on the code close to phonopy-1.10.4. +- Python3 support +- For the RTA thermal conductivity calculation mode with using the + linear tetrahedron method, only necessary part of phonon-phonon + interaction strengh among phonons. This improves lifetime + calculation performance, but as the drawback, averaged ph-ph + interaction strength can not be given. See {ref}`full_pp_option`. +- Pwscf interface ({ref}`calculator_interfaces`) + +## Oct-10-2015: version 0.9.14 + +- Computational performance tuning for phonon-phonon interaction + strength calculation was made by Jonathan Skelton. Depending on + systems, but 10-20% performance improvement may be possible. +- `--stp` option is created to show numbers of q-point triplets to + be calculated. See {ref}`command_options`. +- `--write_gamma` and `--read_gamma` support using with `--bi` + option. Therefore a thermal conductivity calculation can be + distributed over band index, too. This may be useful for the system + whose unit cell is large. + +## Sep-26-2015: version 0.9.13 + +- Changed so that `--wgp` option writes `grid_address-mxxx.hdf5` + instead of `grid_address-mxxx.dat`. +- `--write_detailed_gamma` is implemented. See {ref}`command_options`. +- When running without setting `--thm` and `--sigma` options, + linear tetrahedron method corresponding to `--thm` is used as the + default behavior. +- `--ise` options is created. + +## Aug-12-2015: version 0.9.12 + +- Spglib update to version 1.8.2.1. +- Improve computational performance of `kaccum` and `gaccum`. + +## Jun-18-2015: version 0.9.10.1 + +- Bug fix of `gcaccum` + +## Jun-17-2015: version 0.9.10 + +- Fix bug in `kaccum`. When using with `--pa` option, irreducible + q-points were incorrectly indexed. +- `gaccum` is implemented. `gaccum` is very similar to `kaccum`, + but for $\Gamma_\lambda(\omega_\lambda)$. +- spglib update. + +## Changes in version 0.9.7 + +- The definition of MSPP is modified so as to be averaged ph-ph + interaction defined as $P_{\mathbf{q}j}$ in the arXiv + manuscript. The key in the kappa hdf5 file is changed from `mspp` + to `ave_pp`. The physical unit of $P_{\mathbf{q}j}$ is set + to $\text{eV}^2$. + +## Changes in version 0.9.6 + +- Silicon example is put in `example-phono3py` directory. +- Accumulated lattice thermal conductivity is calculated by `kaccum` + script. +- JDOS output format was changed. + +## Changes in version 0.9.5 + +- In `kappa-xxx.hdf5` file, `heat_capacity` format was changed + from `(irreducible q-point, temperature, phonon band)` to + `(temperature, irreducible q-point, phonon band)`. For `gamma`, + previous document was wrong in the array shape. It is + `(temperature, irreducible q-point, phonon band)` + + +## Changes in version 0.9.4 + +- The option of `--cutoff_mfp` is renamed to `--boundary_mfp` and + now it's on the document. +- Detailed contribution of `kappa` at each **q**-point and phonon + mode is output to .hdf5 with the keyword `mode_kappa`. + +## Changes in version 0.8.11 + +- A new option of `--cutoff_mfp` for including effective boundary + mean free path. +- The option name `--cutfc3` is changed to `--cutoff_fc3`. +- The option name `--cutpair` is changed to `--cutoff_pair`. +- A new option `--ga` is created. +- Fix spectrum plot of joint dos and imaginary part of self energy + +## Changes in version 0.8.10 + +- Different supercell size of fc2 from fc3 can be specified using + `--dim_fc2` option. +- `--isotope` option is implemented. This is used instead of + `--mass_variances` option without specifying the values. Mass + variance parameters are read from database. + +## Changes in version 0.8.2 + +- Phono3py python interface is rewritten and a lot of changes are + introduced. +- `FORCES_SECOND` and `FORCES_THIRD` are no more used. Instead just + one file of `FORCES_FC3` is used. Now `FORCES_FC3` is generated + by `--cf3` option and the backward compatibility is simple: `cat + FORCES_SECOND FORCES_THIRD > FORCES_FC3`. +- `--multiple_sigmas` is removed. The same behavior is achieved by + `--sigma`. + +## Changes in version 0.8.0 + +- `--q_direction` didn't work. Fix it. +- Implementation of tetrahedron method whcih is activated by + `--thm`. +- Grid addresses are written out by `--wgp` option. + +## Changes in version 0.7.6 + +- Cut-off distance for fc3 is implemented. This is activated by + `--cutfc3` option. FC3 elements where any atomic pair has larger + distance than cut-off distance are set zero. +- `--cutpair` works only when creating displacements. The cut-off + pair distance is written into `disp_fc3.yaml` and FC3 is created + from `FORCES_THIRD` with this information. Usually sets of pair + displacements are more redundant than that needed for creating fc3 + if index permutation symmetry is considered. Therefore using index + permutation symmetry, some elements of fc3 can be recovered even if + some of supercell force calculations are missing. In paticular, all + pair distances among triplet atoms are larger than cutoff pair + distance, any fc3 elements are not recovered, i.e., the element will + be zero. + +## Changes in version 0.7.2 + +- Default displacement distance is changed to 0.03. +- Files names of displacement supercells now have 5 digits numbering, + `POSCAR-xxxxx`. +- Cutoff distance between pair displacements is implemented. This is + triggered by `--cutpair` option. This option works only for + calculating atomic forces in supercells with configurations of pairs + of displacements. + +## Changes in version 0.7.1 + +- It is changed to sampling q-points in Brillouin zone. Previously + q-points are sampled in reciprocal primitive lattice. Usually this + change affects very little to the result. +- q-points of phonon triplets are more carefully sampled when a + q-point is on Brillouin zone boundary. Usually this + change affects very little to the result. +- Isotope effect to thermal conductivity is included. + +## Changes in version 0.6.0 + +- `disp.yaml` is renamed to `disp_fc3.yaml`. Old calculations with + `disp.yaml` can be used without any problem just by changing the + file name. +- Group velocity is calculated from analytical derivative of dynamical + matrix. +- Group velocities at degenerate phonon modes are better handled. + This improves the accuracy of group velocity and thus for thermal + conductivity. +- Re-implementation of third-order force constants calculation from + supercell forces, which makes the calculation much faster +- When any phonon of triplets can be on the Brillouin zone boundary, i.e., + when a mesh number is an even number, it is more carefully treated. diff --git a/doc/changelog.rst b/doc/changelog.rst deleted file mode 100644 index 86a060f4..00000000 --- a/doc/changelog.rst +++ /dev/null @@ -1,433 +0,0 @@ -.. _changelog: - -Change Log -========== - -Mar-17-2021: Version 1.22.3 ---------------------------- -- Fix ``--read-gamma`` option to work. - -Feb-21-2021: Version 1.22.2 ---------------------------- -- Fix PyPI source distribution package - -Feb-21-2021: Version 1.22.1 ---------------------------- -- ``phono3py`` command didn't work. This was fixed. -- Fix behaviour when specifying ``--thm`` and ``--sigma`` simultaneously. - -Jan-29-2021: Version 1.22.0 ---------------------------- -- Maintenance release to follow phonopy v2.9.0. - -Sep-30-2020: Version 1.21.0 ---------------------------- - -- Maintenance release to follow the change of phonopy at v2.8.1 -- Improvements of phono3py loader (``phono3py.load``), ``phono3py-load`` - command, API, and ``phono3py_disp.yaml``. -- Harmonic phonon calculation on mesh was multithreaded. This is - effective when using very dense mesh with non-analytical term - correction (probably rare case). -- Real and imaginary parts of self energy and spectral function of - bubble diagram at API level - -Mar-3-2020: Version 1.20.0 --------------------------- - -- ``phono3py_disp.yaml`` is made when creating displacements in - addition to ``disp_fc3.yaml`` and - ``disp_fc2.yaml``. ``phono3py_disp.yaml`` will be used instead of - ``disp_fc3.yaml`` and ``disp_fc2.yaml`` in the future major release - (v2.0). - -Mar-3-2020: Version 1.19.1 --------------------------- - -- Release for pypi and conda (atztogo channel) packagings - -Mar-2-2020: Version 1.19.0 --------------------------- - -- Improvements of phono3py loader and API. -- Improvement of interfaces to calculators. Now it is expected to be - much easier to implement calculator interface if it exists in - phonopy. -- Fixed dependency to phonopy v2.6.0. - -Dec-22-2019: Version 1.18.2 ---------------------------- - -- Initial version of phono3py loader (``phono3py.load``) was - implemented. See docstring of ``phono3py.load``. - -Oct-17-2019: Version 1.18.1 ---------------------------- - -- Fix of phono3py-kaccum to follow the latest phonopy. - -Oct-9-2019: Version 1.18.0 ---------------------------- - -- Maintenance release - -Apr-18-2019: Version 1.17.0 ---------------------------- -- ``--cfz`` option was made to subtract residual forces. See - :ref:`cfz_option`. -- ``--cutoff-pair`` was made to override the cutoff pair distance - written in ``disp_fc3.yaml`` when using on calculating force - constants. This is useful when checking cutoff distance - dependency. So the use case of having fully computed ``FORCES_FC3`` - is assumed. -- TURBOMOLE interface is provided by Antti Karttunen - (``--turbomole``). -- Compatibility of ``fc2.hdf5`` and ``force_constants.hdf5`` was - improved for all calculators to store physical unit information in - the hdf5 file. See :ref:`file_format_compatibility`. - -Mar-24-2019: Version 1.16.0 ---------------------------- -- Bug fixes and catching up the updates of phonopy. -- Most of hdf5 output files are compressed by ``gzip`` as - default. This compression can be set off, see - :ref:`hdf5_compression_option`. -- (Experimental) ``phono3py`` command accepts ``phono3py.yaml`` type - file as an input crystal structure by ``-c`` option. When ``DIM`` - and any structure file are not given, ``phono3py_disp.yaml`` - (primary) or ``phono3py.yaml`` (secondary) is searched in the current - directory. Then ``phono3py.yaml`` type file is used as the input. - By this, semi-automatic phono3py mode is invocked, which acts as - - (1) ``supercell_matrix`` corresponding to ``DIM`` in the - ``phono3py.yaml`` type file is used if it exists. - (2) ``phonon_supercell_matrix`` corresponding to ``DIM_FC2`` in the - ``phono3py.yaml`` type file is used if it exists. - (3) ``primitive_matrix`` in the ``phono3py.yaml`` type file - is used if it exists. Otherwise, set ``PRIMITIVE_AXES = AUTO`` - when ``PRIMITIVE_AXES`` is not given. - (4) NAC params are read (``NAC = .TRUE.``) if NAC params are - contained (primary) in the ``phono3py.yaml`` type file or if - ``BORN`` file exists in the current directory (secondary). - -Nov-22-2018: version 1.14.3 ----------------------------- -- Update to work with phonopy v1.14.2. -- Ph-ph interaction can be read (``--read-pp``) and write - (``--write-pp``) in RTA thermal conductivity calculation, too. Mind - that the data stored are different with and without - ``--full-pp``. Wihtout ``--full-pp`` the data are stored in - complicated way to save data side, so it is not considered readable - by usual users. - -June-20-2018: version 1.13.3 ----------------------------- - -- ``--lw`` (linewidth) option was removed. Use ``--br`` option and - find 2*gamma values as linewidths in ``kappa-xxx.hdf5`` file. -- Documentation of ``--lbte`` option is available at - :ref:`direct_solution`. -- This version is dependent on phonopy>=1.13.2. - -May-17-2018: version 1.13.1 ----------------------------- - -- Compact force constants are implemented (See :ref:`compact_fc_option`). - -Mar-16-2018: version 1.12.9 ----------------------------- - -- Definition of ``mode_kappa`` values in output hdf5 file is - changed. Previously they were divided by number of grid points, but - now not. Therefore users who compute ``kappa`` from ``mode_kappa`` - need to be careful about this change. This does not affect to - ``phono3py-kaccum`` results. - -Feb-1-2018: version 1.12.7 ----------------------------- - -- ``--tsym`` option is removed. Now with ``--sym-fc3r`` and - ``--sym-fc2`` options, - translational invariance symmetry is also applied. -- ``--sym-fc`` option is added. This is just an alias to specify both - ``--sym-fc3r`` and ``--sym-fc2`` together. -- Documentation on ``--write-phonon`` and ``--read-phonon`` options is - written. These options are used to save harmonic phonon infromation - on strage. - -Nov-22-2017: version 1.12.5 ------------------------------ - -- Bug fix of RTA thermal conductivity. This bug exists from version - 1.10.11.18 (git e40cd059). This bug exhibits when all the following - conditions are met: - - 1. RTA thermal conductivity calculation. - 2. Tetrahedron method for Brillouin zone integration is used. - 3. Number of triplets is smaller than number of bands at each grid point. - 4. Without using ``--full-pp``. - - - (3) happens when the primitive cell is relatively large. Number of - triplets can be shown using ``--stp`` option. A race condition of - OpenMP multithreading is the source of the bug. Therefore, if it - occurs, the same calculation comes up with the different thermal - conductivity value in every run time, for whcih it behaves like - randomly. - -- RTA thermal conductivity with smearing method (``--sigma``) is made - to run with smaller memory consumption as similar as tetrahedron - method (``--thm``). - -Nov-17-2017: version 1.12.3 ----------------------------- - -- Command option parser of the phonopy tools is replaced from - ``optparse`` to ``argparse``. -- The filenames used with these options were the positional arguments - previously. Now they are the command-line arguments, i.e., filenames - have to be put just after the option name like ``-f vasprun.xml-001 - vasprun.xml-002 ...``. -- The names of auxiliary tools (``kdeplot`` and ``kaccum``) are - changed, for which the prefix phono3py- is added to the old names to - avoid accidental conflict with other script names already existing - under bin directory. -- :ref:`sigma_cutoff_option` option was created. - -Jun-18-2017: version 1.11.13 ----------------------------- - -- :ref:`normal_umklapp_option` option was made. -- Many minor updates: fixing bugs, improving usabilities. -- Improve of :ref:`auxiliary_tools_kaccum` and :ref:`auxiliary_tools_kdeplot`. - -Mar-31-2017: version 1.11.11 ----------------------------- - -- Abinit code interface is implemented and now under the testing. -- Reduction of memory usage in RTA thermal conductivity - calculation. This is especially effective for larger unit cell - case. Currently combinations with --full_pp, --write_gamma_detail, - and --simga(smearing method) are not supported for this. Performance - tuning is under going. In some case, computation can be slower than - the previous versions. - -Feb-9-2017: version 1.11.9 ---------------------------- - -- This version works coupled with phonopy-1.11.8 or later. -- CRYSTAL code interface is implemented by Antti J. Karttunen. - -Dec-14-2016: version 1.11.7 ------------------------------- - -- This is a maintenance release. This version must be used with - phonopy-1.11.6 or later. - -Nov-27-2016: version 1.11.5 ------------------------------- - -- ``gaccum`` is merged to ``kaccum``. ``gaccum`` is removed. See - :ref:`auxiliary_tools_kaccum`. -- ``kdeplot`` is added. See :ref:`auxiliary_tools_kdeplot`. - -Apr-24-2016: version 1.10.9 ------------------------------- - -- Failure of writing ``kappa-mxxx-gx.hdf5`` was fixed. - -Apr-16-2016: version 1.10.7 ------------------------------- - -- API example is prepared and it is found in ``Si`` example. No - doucment yet. -- Si pwscf example was placed in ``example-phono3py`` directory. -- User interface bug fix. - -Mar-15-2016: version 1.10.5 ------------------------------- - -- Numbering way of phono3py version was just changed (No big updates - were made against previous version.) The number is given based on - the phonopy version. For example, the harmonic part of - phono3py-1.10.5 is based on the code close to phonopy-1.10.4. -- Python3 support -- For the RTA thermal conductivity calculation mode with using the - linear tetrahedron method, only necessary part of phonon-phonon - interaction strengh among phonons. This improves lifetime - calculation performance, but as the drawback, averaged ph-ph - interaction strength can not be given. See :ref:`full_pp_option`. -- Pwscf interface (:ref:`calculator_interfaces`) - -Oct-10-2015: version 0.9.14 ------------------------------- - -- Computational performance tuning for phonon-phonon interaction - strength calculation was made by Jonathan Skelton. Depending on - systems, but 10-20% performance improvement may be possible. -- ``--stp`` option is created to show numbers of q-point triplets to - be calculated. See :ref:`command_options`. -- ``--write_gamma`` and ``--read_gamma`` support using with ``--bi`` - option. Therefore a thermal conductivity calculation can be - distributed over band index, too. This may be useful for the system - whose unit cell is large. - -Sep-26-2015: version 0.9.13 ------------------------------- - -- Changed so that ``--wgp`` option writes ``grid_address-mxxx.hdf5`` - instead of ``grid_address-mxxx.dat``. -- ``--write_detailed_gamma`` is implemented. See :ref:`command_options`. -- When running without setting ``--thm`` and ``--sigma`` options, - linear tetrahedron method corresponding to ``--thm`` is used as the - default behavior. -- ``--ise`` options is created. - -Aug-12-2015: version 0.9.12 ------------------------------- - -- Spglib update to version 1.8.2.1. -- Improve computational performance of ``kaccum`` and ``gaccum``. - -Jun-18-2015: version 0.9.10.1 ------------------------------- - -- Bug fix of ``gcaccum`` - -Jun-17-2015: version 0.9.10 ----------------------------- - -- Fix bug in ``kaccum``. When using with ``--pa`` option, irreducible - q-points were incorrectly indexed. -- ``gaccum`` is implemented. ``gaccum`` is very similar to ``kaccum``, - but for :math:`\Gamma_\lambda(\omega_\lambda)`. -- spglib update. - -Changes in version 0.9.7 -------------------------- - -- The definition of MSPP is modified so as to be averaged ph-ph - interaction defined as :math:`P_{\mathbf{q}j}` in the arXiv - manuscript. The key in the kappa hdf5 file is changed from ``mspp`` - to ``ave_pp``. The physical unit of :math:`P_{\mathbf{q}j}` is set - to :math:`\text{eV}^2`. - -Changes in version 0.9.6 ------------------------- - -- Silicon example is put in ``example-phono3py`` directory. -- Accumulated lattice thermal conductivity is calculated by ``kaccum`` - script. -- JDOS output format was changed. - -Changes in version 0.9.5 -------------------------- - -- In ``kappa-xxx.hdf5`` file, ``heat_capacity`` format was changed - from ``(irreducible q-point, temperature, phonon band)`` to - ``(temperature, irreducible q-point, phonon band)``. For ``gamma``, - previous document was wrong in the array shape. It is - ``(temperature, irreducible q-point, phonon band)`` - - -Changes in version 0.9.4 ------------------------- - -- The option of ``--cutoff_mfp`` is renamed to ``--boundary_mfp`` and - now it's on the document. -- Detailed contribution of ``kappa`` at each **q**-point and phonon - mode is output to .hdf5 with the keyword ``mode_kappa``. - -Changes in version 0.8.11 -------------------------- - -- A new option of ``--cutoff_mfp`` for including effective boundary - mean free path. -- The option name ``--cutfc3`` is changed to ``--cutoff_fc3``. -- The option name ``--cutpair`` is changed to ``--cutoff_pair``. -- A new option ``--ga`` is created. -- Fix spectrum plot of joint dos and imaginary part of self energy - -Changes in version 0.8.10 -------------------------- - -- Different supercell size of fc2 from fc3 can be specified using - ``--dim_fc2`` option. -- ``--isotope`` option is implemented. This is used instead of - ``--mass_variances`` option without specifying the values. Mass - variance parameters are read from database. - -Changes in version 0.8.2 ------------------------- - -- Phono3py python interface is rewritten and a lot of changes are - introduced. -- ``FORCES_SECOND`` and ``FORCES_THIRD`` are no more used. Instead just - one file of ``FORCES_FC3`` is used. Now ``FORCES_FC3`` is generated - by ``--cf3`` option and the backward compatibility is simple: ``cat - FORCES_SECOND FORCES_THIRD > FORCES_FC3``. -- ``--multiple_sigmas`` is removed. The same behavior is achieved by - ``--sigma``. - -Changes in version 0.8.0 ------------------------- - -- ``--q_direction`` didn't work. Fix it. -- Implementation of tetrahedron method whcih is activated by - ``--thm``. -- Grid addresses are written out by ``--wgp`` option. - -Changes in version 0.7.6 ------------------------- - -- Cut-off distance for fc3 is implemented. This is activated by - ``--cutfc3`` option. FC3 elements where any atomic pair has larger - distance than cut-off distance are set zero. -- ``--cutpair`` works only when creating displacements. The cut-off - pair distance is written into ``disp_fc3.yaml`` and FC3 is created - from ``FORCES_THIRD`` with this information. Usually sets of pair - displacements are more redundant than that needed for creating fc3 - if index permutation symmetry is considered. Therefore using index - permutation symmetry, some elements of fc3 can be recovered even if - some of supercell force calculations are missing. In paticular, all - pair distances among triplet atoms are larger than cutoff pair - distance, any fc3 elements are not recovered, i.e., the element will - be zero. - -Changes in version 0.7.2 ------------------------- - -- Default displacement distance is changed to 0.03. -- Files names of displacement supercells now have 5 digits numbering, - ``POSCAR-xxxxx``. -- Cutoff distance between pair displacements is implemented. This is - triggered by ``--cutpair`` option. This option works only for - calculating atomic forces in supercells with configurations of pairs - of displacements. - -Changes in version 0.7.1 ------------------------- - -- It is changed to sampling q-points in Brillouin zone. Previously - q-points are sampled in reciprocal primitive lattice. Usually this - change affects very little to the result. -- q-points of phonon triplets are more carefully sampled when a - q-point is on Brillouin zone boundary. Usually this - change affects very little to the result. -- Isotope effect to thermal conductivity is included. - -Changes in version 0.6.0 ------------------------- - -- ``disp.yaml`` is renamed to ``disp_fc3.yaml``. Old calculations with - ``disp.yaml`` can be used without any problem just by changing the - file name. -- Group velocity is calculated from analytical derivative of dynamical - matrix. -- Group velocities at degenerate phonon modes are better handled. - This improves the accuracy of group velocity and thus for thermal - conductivity. -- Re-implementation of third-order force constants calculation from - supercell forces, which makes the calculation much faster -- When any phonon of triplets can be on the Brillouin zone boundary, i.e., - when a mesh number is an even number, it is more carefully treated. diff --git a/doc/citation.md b/doc/citation.md new file mode 100644 index 00000000..a7864ce3 --- /dev/null +++ b/doc/citation.md @@ -0,0 +1,50 @@ +# 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} + } + ``` \ No newline at end of file diff --git a/doc/command-options.md b/doc/command-options.md new file mode 100644 index 00000000..3692b084 --- /dev/null +++ b/doc/command-options.md @@ -0,0 +1,1041 @@ +(command_options)= +# Command options / Setting tags + +Phono3py is operated with command options or with a configuration file +that contains setting tags. In this page, the command options are +explained. Most of command options have their respective setting +tags. + +A configuration file with setting tags like phonopy can be used +instead of and together with the command options. The setting tags are +mostly equivalent to the respective most command options, but when +both are set simultaneously, the command options are preferred. An +example of configuration (e.g., saved in a file `setting.conf`) is +as follow: + +```bash +DIM = 2 2 2 +DIM_FC2 = 4 4 4 +PRIMITIVE_AXES = 0 1/2 1/2 1/2 0 1/2 1/2 1/2 0 +MESH = 11 11 11 +BTERTA = .TRUE. +NAC = .TRUE. +READ_FC2 = .TRUE. +READ_FC3 = .TRUE. +CELL_FILENAME = POSCAR-unitcell +``` + +where the setting tag names are case insensitive. This is run by + +```bash +% phono3py setting.conf [comannd options] +``` + +or + +```bash +% phono3py [comannd options] -- setting.conf +``` + +```{contents} +:depth: 2 +:local: +``` + +## Input cell file name + +### `-c` (`CELL_FILENAME`) + +This specifies input unit cell filename. + +```bash +% phono3py -c POSCAR-unitcell ... (many options) +``` + +## Calculator interface + + +### `--qe` (`CALCULATOR = QE`) + +Quantum espresso (pw) interface is invoked. +See the detail at {ref}`qe_interface`. + +### `--crystal` (`CALCULATOR = CRYSTAL`) + +CRYSTAL interface is invoked. +See the detail at {ref}`crystal_interface`. + +### `--turbomole` (`CALCULATOR = TURBOMOLE`) + +TURBOMOLE interface is invoked. +See the details at {ref}`turbomole_interface`. + +## Utilities to create default input files + +These options have no respective configuration file tags. + +### `--cf3` (command option only) + +This is used to create `FORCES_FC3` from `phono3py_disp.yaml` and +force calculator outputs containing forces in supercells. `phono3py_disp.yaml` +has to be located at the current directory. Calculator interface has +to be specified except for VASP (default) case. + +```bash +% phono3py --cf3 disp-{00001..00755}/vasprun.xml +``` + +```bash +% phono3py --cf3 supercell_out/disp-{00001..00111}/Si.out +``` + +(cf3_file_option)= +### `--cf3-file` (command option only) + +This is used to create `FORCES_FC3` from a text file containing a +list of calculator output file names. `phono3py_disp.yaml` has to be +located at the current directory. Calculator interface has to be +specified except for VASP (default) case. + +```bash +% phono3py --cf3-file file_list.dat +``` + +where `file_list.dat` contains file names that can be recognized +from the current directory and is expected to be like: + +```bash +disp-00001/vasprun.xml +disp-00002/vasprun.xml +disp-00003/vasprun.xml +disp-00004/vasprun.xml +... +``` + +The order of the file names is important. This option may be useful +to be used together with `--cutoff-pair` option. + +(cf2_option)= +### `--cf2` (command option only) + +This is used to create `FORCES_FC2` similarly to `--cf3` +option. `phono3py_disp.yaml` has to be located at the current +directory. This is optional. Calculator interface has to be specified +except for VASP (default) case. `FORCES_FC2` is necessary to run +with `--dim-fc2` option. + +```bash +% phono3py --cf2 disp_fc2-{00001..00002}/vasprun.xml +``` + +(cfz_option)= +### `--cfz` (command option only) + +This is used to create `FORCES_FC3` and `FORCES_FC2` subtracting +residual forces combined with `--cf3` and `--cf2`, +respectively. Calculator interface has to be specified except for VASP +(default) case. + +In the following example, it is supposed that +`disp3-00000/vasprun.xml` and `disp2-00000/vasprun.xml` contain +the forces of the perfect supercells. In ideal case, these forces are +zero, but often they are not. Here, this is called "residual +forces". Sometimes quality of force constants is improved in this way. + +```bash +% phono3py --cf3 disp3-{00001..01254}/vasprun.xml --cfz disp3-00000/vasprun.xml +% phono3py --cf2 disp2-{00001..00006}/vasprun.xml --cfz disp2-00000/vasprun.xml +``` + +(fs2f2_option)= +### `--fs2f2` or `--force-sets-to-forces-fc2` (command option only) + +`FORCES_FC2` is created from phonopy's `FORCE_SETS` file. +Necessary yaml lines for `phono3py_disp.yaml` is displayed as text. + +```bash +% phono3py --fs2f2 +``` + +(cfs_option)= +### `--cfs` or `--create-force-sets` (command option only) + +Phonopy's `FORCE_SETS` is created from +`FORCES_FC3` and `phono3py_disp.yaml`. + +```bash +% phono3py --cfs +``` + +In conjunction with {ref}`--dim-fc2 `, phonopy's +`FORCE_SETS` is created from `FORCES_FC2` and `phono3py_disp.yaml` +instead of `FORCES_FC3` and `phono3py_disp.yaml`. + +```bash +% phono3py --cfs --dim-fc2="x x x" +``` + +## Supercell and primitive cell + +(dim_option)= +### `--dim` (`DIM`) + +Supercell dimension is specified. See the +detail at http://phonopy.github.io/phonopy/setting-tags.html#dim. +When a proper `phono3py_disp.yaml` exists in the current directory, +this is unnecessary to be specified. + +(dim_fc2_option)= +### `--dim-fc2` (`DIM_FC2`) + +Supercell dimension for 2nd order force constants (for harmonic +phonons) is specified. This is optional. +When a proper `phono3py_disp.yaml` exists in the current directory, +this is unnecessary to be specified. + +A larger and different supercell size for 2nd order force constants +than that for 3rd order force constants can be specified with this +option. Often interaction between a pair of atoms has longer range in +real space than interaction among three atoms. Therefore to reduce +computational demand, choosing larger supercell size only for 2nd +order force constants may be a good idea. + +Using this option with `-d` option, the structure files +(e.g. `POSCAR_FC2-xxxxx` or equivalent files for the other +interfaces) and `phono3py_disp.yaml` are created. These are used to +calculate 2nd order force constants for the larger supercell size and +these force calculations have to be done in addition to the usual +force calculations for 3rd order force constants. + +```bash +% phono3py -d --dim="2 2 2" --dim-fc2="4 4 4" -c POSCAR-unitcell +``` + +After the force calculations, `--cf2` option is used to create +`FORCES_FC2`. + +```bash +% phono3py --cf2 disp-{001,002}/vasprun.xml +``` + +To calculate 2nd order force constants for the larger supercell size, +`FORCES_FC2` and `phono3py_disp.yaml` are necessary. Whenever running +phono3py for the larger 2nd order force constants, `--dim-fc2` +option has to be specified. `fc2.hdf5` created as a result of +running phono3py contains the 2nd order force constants with +larger supercell size. The filename is the same as that created in the +usual phono3py run without `--dim-fc2` option. + +```bash +% phono3py --dim="2 2 2" --dim_fc2="4 4 4" -c POSCAR-unitcell ... (many options) +``` + +(pa_option)= +### `--pa`, `--primitive-axes` (`PRIMITIVE_AXES`) + +Transformation matrix from a non-primitive cell to the primitive +cell. See phonopy `PRIMITIVE_AXES` tag (`--pa` option) at +http://phonopy.github.io/phonopy/setting-tags.html#primitive-axis. +When a proper `phono3py_disp.yaml` exists in the current directory, +this is unnecessary to be specified. + +## Displacement creation + +(create_displacements_option)= +### `-d` (`CREATE_DISPLACEMENTS = .TRUE.`) + +Supercell with displacements are created. Using with `--amplitude` +option, atomic displacement distances are controlled. With this +option, files for supercells with displacements and `phono3py_disp.yaml` +file are created. `--pa` should be specified if the input unit cell +structure is not a primitive cell, e.g., `--pa="F"` if the input +unit cell has F-centring. + +(amplitude_option)= +### `--amplitude` (`DISPLACEMENT_DISTANCE`) + +Atomic displacement distance is specified. This +value may be increased for the weak interaction systems and descreased +when the force calculator is numerically very accurate. + +The default value depends on calculator. See +{ref}`default_displacement_distance_for_calculator`. + +## Force constants + +(compact_fc_option)= +### `--cfc` or `--compact-fc` (`COMPACT_FC = .TRUE.`) + +When creating force constants from `FORCES_FC3` and/or +`FORCES_FC2`, force constants that use smaller data size are +created. The shape of the data array is `(num_patom, num_satom)` for +fc2 and `(num_patom, num_satom, num_satom)` for fc3, where +`num_patom` and `num_satom` are the numbers of atoms in primtive +cell and supercell. In the full size force constants case, +`num_patom` is replaced by `num_satom`. Therefore if the supercell +dimension is large, this reduction of data size becomes large. If the +input crystal structure has centring {ref}`--pa ` is +necessary to have smallest data size. In this case, `--pa` option +has to be specified on reading. Otherwise phono3py can recognize if +`fc2.hdf5` and `fc3.hdf5` are compact or full automatically. When +using with `--sym-fc`, the calculated results will become slightly +different due to imperfect symmetrization scheme that phono3py +employs. + +```bash +% phono3py --dim="2 2 2" --cfc --pa="F" -c POSCAR-unitcell +``` + +(symmetrization_option)= +### `--sym-fc` (`FC_SYMMETRY = .TRUE.`) + +Second- and third-order force constants are symmetrized. The index +exchange of real space force constantsand translational invariance +symmetry are applied in a simple way. This symmetrization just removes +drift force constants evenly from all elements and then applies +averaging index-exchange equivalent elements. Therefore the different +symmetries are not simultaneously enforced. For better symmetrization, +it is recommended to use an external force constants calculator like ALM. + +The symmetrizations for the second and third orders can be +independently applied by `--sym-fc2` (`SYMMETRIZE_FC2 = .TRUE.`) +and `--sym-fc3r` (`SYMMETRIZE_FC3 = .TRUE.`), , respectively. + +(cf3_option)= +### `--cutoff-fc3` or `--cutoff-fc3-distance` (`CUTOFF_FC3_DISTANCE`) + +This option is **not** used to reduce number of supercells with +displacements, but this option is used to set zero in elements of +given third-order force constants. The zero elements are selected by +the condition that any pair-distance of atoms in each atom triplet is +larger than the specified cut-off distance. + +If one wants to reduce number of supercells, the first choice is to +reduce the supercell size and the second choice is using +`--cutoff-pair` option. + +(cutoff_pair_option)= +### `--cutoff-pair` or `--cutoff-pair-distance` (`CUTOFF_PAIR_DISTANCE`) + +This option is only used together with `-d` option. + +A cutoff pair-distance in a supercell is used to reduce the number of +necessary supercells with displacements to obtain third order force +constants. As the drawback, a certain number of +third-order-force-constants elements are abandoned or computed with +less numerical accuracy. More details are found at {ref}`command_cutoff_pair`. + +## Reciprocal space sampling mesh and grid points, and band indices + +### `--mesh` (`MESH` or `MESH_NUMBERS`) + +Mesh sampling grids in reciprocal space are generated with the +specified numbers. This mesh is made along reciprocal axes and +is always Gamma-centered. Except for that this mesh is always +Gamma-centered, this works in the same way as written here, +https://phonopy.github.io/phonopy/setting-tags.html#mesh-mp-or-mesh-numbers. + +(gp_option)= +### `--gp` (`GRID_POINTS`) + +Grid points are specified by their unique indices, e.g., for selecting +the q-points where imaginary parts of self energees are +calculated. For thermal conductivity calculation, this can be +used to distribute its calculation over q-points (see +{ref}`workload_distribution`). + +Indices of grid points are specified by space or comma (`,`) +separated numbers. The mapping table between grid points to its +indices is obtained by running with `--loglevel=2` option. + +```bash +% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --fc3 --fc2 --br --write-gamma --gp="0 1 2 3 4 5" +``` + +where `--gp="0 1 2 3 4 5"` can be also written +`--gp="0,1,2,3,4,5"`. `--ga` option below can be used similarly +for the same purpose. + +(ga_option)= +### `--ga` (`GRID_ADDRESSES`) + +This is used to specify grid points like `--gp` option but in their +addresses represented by integer numbers. For example with +`--mesh="16 16 16"`, a q-point of (0.5, 0.5, 0.5) is given by +`--ga="8 8 8"`. The values have to be integers. If you want to +specify the point on a path, `--ga="0 0 0 1 1 1 2 2 2 3 3 3 ..."`, +where each three values are recogninzed as a grid point. The grid +points given by `--ga` option are translated to grid point indices +as given by `--gp` option, and the values given by `--ga` option +will not be shown in log files. + +(bi_option)= +### `--bi` (`BAND_INDICES`) + +Band indices are specified. The output file name will be, e.g., +`gammas-mxxx-gxx(-sx)-bx.dat` where `bxbx...` shows the band +indices used to be averaged. The calculated values at indices +separated by space are averaged, and those separated by comma are +separately calculated. + +```bash +% phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --gp="34" --bi="4 5, 6" +``` + +This may be also useful to distribute the computational demand +such like that the unit cell is large and the calculation of +phonon-phonon interaction is heavy. + +(wgp_option)= +### `--wgp` (command option only) + +Irreducible grid point indices and related information are written +into `ir_grid_points.yaml`. This information may be used when we +want to distribute thermal conductivity calculation into small pieces +or to find specific grid points to calculate imaginary part of self +energy, for which {ref}`--gp option ` can be used to +specify the grid point indices. + +`grid_address-mxxx.hdf5` is also written. This file contains all the +grid points and their grid addresses in integers. Q-points +corresponding to grid points are calculated divided these integers by +sampling mesh numbers for respective reciprocal axes. + +```bash +% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --wgp +``` + +(stp_option)= +### `--stp` (command option only) + +Numbers of q-point triplets to be calculated for irreducible grid +points for specified sampling mesh numbers are shown. This can be used +to estimate how large a calculation is. Only those for specific grid +points are shown by using with `--gp` or `--ga` option. + +```bash +% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="19 19 19" --stp --gp 20 +``` + +## Brillouin zone integration + +(thm_option)= +### `--thm` (`TETRAHEDRON = .TRUE.`) + +Tetrahedron method is used for calculation of imaginary part of self +energy. This is the default option. Therefore it is not necessary to +specify this unless both results by tetrahedron method and +smearing method in one time execution are expected. + +(sigma_option)= +### `--sigma` (`SIGMA`) + +$\sigma$ value of Gaussian function for smearing when +calculating imaginary part of self energy. See the detail at +{ref}`brillouinzone_sum`. + +Multiple $\sigma$ values are also specified by space separated +numerical values. This is used when we want to test several +$\sigma$ values simultaneously. + +(sigma_cutoff_option)= +### `--sigma-cutoff` (`SIGMA_CUTOFF_WIDTH`) + +The tails of the Gaussian functions that are used to replace delta +functions in the equation shown at {ref}`--full-pp ` +are cut with this option. The value is specified in number of standard +deviation. `--sigma-cutoff=5` gives the Gaussian functions to be cut +at $5\sigma$. Using this option scarifies the numerical +accuracy. So the number has to be carefully tested. But computation of +phonon-phonon interaction strength becomes much faster in exchange for +it. + +(full_pp_option)= +### `--full-pp` (`FULL_PP = .TRUE.`) + +For thermal conductivity calculation using the linear tetrahedron +method (from version 1.10.5) and smearing method with +`--simga-cutoff` (from version 1.12.3), only necessary elements +(i.e., that have non-zero delta functions) of phonon-phonon interaction strength, +$\bigl|\Phi_{-\lambda\lambda'\lambda''}\bigl|^2$, is calculated +due to delta functions in calculation of +$\Gamma_\lambda(\omega)$, + +$$ +\Gamma_\lambda(\omega) = \frac{18\pi}{\hbar^2} + \sum_{\lambda' \lambda''} + \bigl|\Phi_{-\lambda\lambda'\lambda''}\bigl|^2 + \left\{(n_{\lambda'}+ n_{\lambda''}+1) + \delta(\omega-\omega_{\lambda'}-\omega_{\lambda''}) \right. + + (n_{\lambda'}-n_{\lambda''}) + \left[\delta(\omega+\omega_{\lambda'}-\omega_{\lambda''}) +- \left. \delta(\omega-\omega_{\lambda'}+\omega_{\lambda''}) +\right]\right\}. +$$ + +But using this option, full elements of phonon-phonon interaction +strength are calculated and averaged phonon-phonon interaction +strength ($P_{\mathbf{q}j}$, see {ref}`--ave-pp +`) is also given and stored. + +## Method to solve BTE + +### `--br` (`BTERTA = .TRUE.`) + +Run calculation of lattice thermal conductivity tensor with the single +mode relaxation time approximation (RTA) and linearized phonon +Boltzmann equation. Without specifying `--gp` (or `--ga`) option, +all necessary phonon lifetime calculations for grid points are +sequentially executed and then thermal conductivity is calculated +under RTA. The thermal conductivity and many related properties are +written into `kappa-mxxx.hdf5`. + +With `--gp` (or `--ga`) option, +phonon lifetimes on the specified grid points are calculated. To save +the results, `--write-gamma` option has to be specified and the +physical properties belonging to the grid +points are written into `kappa-mxxx-gx(-sx).hdf5`. + +### `--lbte` (`LBTE = .TRUE.`) + +Run calculation of lattice thermal conductivity tensor with a direct +solution of linearized phonon Boltzmann equation. The basis usage of +this option is equivalent to that of `--br`. More detail is +documented at {ref}`direct_solution`. + +## Scattering + +### `--isotope` (`ISOTOPE =.TRUE.`) + +Phonon-isotope scattering is calculated based on the formula by +Shin-ichiro Tamura, Phys. Rev. B, **27**, 858 (1983). Mass variance +parameters are read from database of the natural abundance data for +elements, which refers Laeter *et al.*, Pure Appl. Chem., **75**, 683 +(2003). + +```bash +% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --isotope +``` + +### `--mass-variances` or `--mv` (`MASS_VARIANCES`) + +Mass variance parameters are specified by this option to include +phonon-isotope scattering effect in the same way as `--isotope` +option. For example of GaN, this may be set like `--mv="1.97e-4 +1.97e-4 0 0"`. The number of elements has to correspond to the number +of atoms in the primitive cell. + +Isotope effect to thermal conductivity may be checked first running +without isotope calculation: + +```bash +% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br +``` + +Then running with isotope calculation: + +```bash +% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --read-gamma --mv="1.97e-4 1.97e-4 0 0" +``` + +In the result hdf5 file, currently isotope scattering strength is not +written out, i.e., `gamma` is still imaginary part of self energy of +ph-ph scattering. + +### `--boundary-mfp`, `--bmfp` (`BOUNDARY_MFP`) + +A most simple phonon boundary scattering treatment is +included. $v_g/L$ is just used as the scattering rate, where +$v_g$ is the group velocity and $L$ is the boundary mean +free path. The value is given in micrometre. The default value, 1 +metre, is just used to avoid divergence of phonon lifetime and the +contribution to the thermal conducitivity is considered negligible. + +(ave_pp_option)= +### `--ave-pp` (`USE_AVE_PP = .TRUE.`) + +Averaged phonon-phonon interaction strength ($P_{\mathbf{q}j}=P_\lambda$) +is used to calculate imaginary part of self energy in thermal +conductivity calculation. $P_\lambda$ is defined as + +$$ +P_\lambda = \frac{1}{(3n_\text{a})^2}\sum_{\lambda' +\lambda''}|\Phi_{\lambda \lambda' \lambda''}|^2, +$$ + +where $n_\text{a}$ is the number of atoms in unit cell. This is +roughly constant with respect to the sampling mesh density for +converged $|\Phi_{\lambda \lambda' \lambda''}|^2$. Then for all +$\mathbf{q}',j',j''$, + +$$ +|\Phi_{\mathbf{q}j,\mathbf{q}'j',\mathbf{G-q-q'}j''}|^2 := +P_{\mathbf{q}j} / N, +$$ + +where $N$ is the number of grid points on the sampling +mesh. $\Phi_{\lambda \lambda' \lambda''} \equiv 0$ unless +$\mathbf{q} + \mathbf{q}' + \mathbf{q}'' = \mathbf{G}$. + +This option works only when `--read-gamma` +and `--br` options are activated where the averaged phonon-phonon +interaction that is read from `kappa-mxxx(-sx-sdx).hdf5` file is +used if it exists in the file. Therefore the averaged phonon-phonon +interaction has to be stored before using this option (see +{ref}`--full-pp `). The calculation result **overwrites** +`kappa-mxxx(-sx-sdx).hdf5` file. Therefore to use this option +together with `-o` option is strongly recommended. + +First, run full conductivity calculation, + +```bash +% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br +``` + +Then + +```bash +% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --read-gamma --ave-pp -o ave_pp +``` + +### `--const-ave-pp` (`CONST_AVE_PP = .TRUE.`) + +Averaged phonon-phonon interaction ($P_{\mathbf{q}j}$) is +replaced by this constant value and $|\Phi_{\lambda \lambda' +\lambda''}|^2$ are set as written in {ref}`--ave-pp ` for thermal +conductivity calculation. This option works only when `--br` options +are activated. Therefore third-order force constants are not necessary +to input. The physical unit of the value is $\text{eV}^2$. + +```bash +% phono3py --dim="3 3 2" -v --mesh="32 32 20" -c POSCAR-unitcell --br --const-ave-pp=1e-10 +``` + +(normal_umklapp_option)= +### `--nu` (`N_U = .TRUE.`) + +Integration over q-point triplets for the calculation of +$\Gamma_\lambda(\omega_\lambda)$ is made separately for normal +$\Gamma^\text{N}_\lambda(\omega_\lambda)$ and Umklapp +$\Gamma^\text{U}_\lambda(\omega_\lambda)$ processes. The sum of +them is usual $\Gamma_\lambda(\omega_\lambda) = +\Gamma^\text{N}_\lambda(\omega_\lambda) + +\Gamma^\text{U}_\lambda(\omega_\lambda)$ and this is used to calcualte +thermal conductivity in single-mode RTA. The separation, i.e., the +choice of G-vector, is made based on the first Brillouin zone. + +The data are stored in `kappa-mxxx(-gx-sx-sdx).hdf5` file and +accessed by `gamma_N` and `gamma_U` keys. The shape of the arrays +is the same as that of `gamma` (see +{ref}`kappa_hdf5_file_gamma`). An example (Si-PBEsol) is shown below: + +```bash +% phono3py --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="11 11 11" --fc3 --fc2 --br --nu +... +% ipython + +In [1]: import h5py + +In [2]: f = h5py.File("kappa-m111111.hdf5", 'r') + +In [3]: list(f) +Out[3]: +['frequency', + 'gamma', + 'gamma_N', + 'gamma_U', + 'group_velocity', + 'gv_by_gv', + 'heat_capacity', + 'kappa', + 'kappa_unit_conversion', + 'mesh', + 'mode_kappa', + 'qpoint', + 'temperature', + 'weight'] + +In [4]: f['gamma'].shape +Out[4]: (101, 56, 6) + +In [5]: f['gamma_N'].shape +Out[5]: (101, 56, 6) + +In [6]: f['gamma_U'].shape +Out[6]: (101, 56, 6) +``` + +## Temperature + +(ts_option)= +### `--ts` (`TEMPERATURES`): Temperatures + +Specific temperatures are specified by `--ts`. + +```bash +% phono3py --fc3 --fc2 --dim="2 2 2" -v --mesh="11 11 11" -c POSCAR-unitcell --br --ts="200 300 400" +``` + +### `--tmax`, `--tmin`, `--tstep` (`TMAX`, `TMIN`, `TSTEP`) + +Temperatures at equal interval are specified by `--tmax`, +`--tmin`, `--tstep`. See phonopy's document for the same tags at +http://phonopy.github.io/phonopy/setting-tags.html#tprop-tmin-tmax-tstep. + +```bash +% phono3py --fc3 --fc2 --dim="2 2 2" -v --mesh="11 11 11" -c POSCAR-unitcell --br --tmin=100 --tmax=1000 --tstep=50 +``` + +## Non-analytical term correction + +(nac_option)= +### `--nac` (`NAC = .TRUE.`) + +Non-analytical term correction for harmonic phonons. Like as phonopy, +`BORN` file has to be put on the same directory. Always the default +value of unit conversion factor is used even if it is written in the +first line of `BORN` file. + +### `--q-direction` (`Q_DIRECTION`) + +This is used with `--nac` to specify reciprocal-space direction +at $\mathbf{q}\rightarrow \mathbf{0}$. See the detail +at http://phonopy.github.io/phonopy/setting-tags.html#q-direction. + +(write_gamma_option)= +## Imaginary part of self energy + +(ise_option)= +### `--ise` (`IMAG_SELF_ENERGY = .TRUE.`) + +Imaginary part of self energy $\Gamma_\lambda(\omega)$ is +calculated with respect to $\omega$. The output is written to +`gammas-mxxx-gx(-sx)-tx-bx.dat` in THz (without $2\pi$) +with respect to frequency in THz (without $2\pi$). Frequency sampling +points can be specified by {ref}`freq_sampling_option`. + +```bash +% phono3py --fc3 --fc2 --dim="2 2 2" --mesh="16 16 16" -c POSCAR-unitcell --nac --q-direction="1 0 0" --gp=0 --ise --bi="4 5, 6" +``` + +## Joint density of states + +(jdos_option)= +### `--jdos` (`JOINT_DOS = .TRUE.`) + +Two classes of joint density of states (JDOS) are calculated. The +result is written into `jdos-mxxx-gx(-sx-sdx).dat` in +$\text{THz}^{-1}$ (without $(2\pi)^{-1}$) with +respect to frequency in THz (without $2\pi$). Frequency sampling +points can be specified by {ref}`freq_sampling_option`. + +The first column is the frequency, and the second and third columns +are the values given as follows, respectively, + +$$ +&D_2^{(1)}(\mathbf{q}, \omega) = \frac{1}{N} +\sum_{\lambda',\lambda''} \Delta(-\mathbf{q}+\mathbf{q}'+\mathbf{q}'') +\left[\delta(\omega+\omega_{\lambda'}-\omega_{\lambda''}) + +\delta(\omega-\omega_{\lambda'}+\omega_{\lambda''}) \right], \\ +&D_2^{(2)}(\mathbf{q}, \omega) = \frac{1}{N} +\sum_{\lambda',\lambda''} +\Delta(-\mathbf{q}+\mathbf{q}'+\mathbf{q}'') \delta(\omega-\omega_{\lambda'} +-\omega_{\lambda''}). +$$ + +```bash +% phono3py --fc2 --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="16 16 16" --jdos --ga="0 0 0 8 8 8" +``` + +When temperatures are specified, two classes of weighted JDOS are +calculated. The result is written into +`jdos-mxxx-gx(-sx)-txxx.dat` in $\text{THz}^{-1}$ (without +$(2\pi)^{-1}$) with respect to frequency in THz (without +$2\pi$). In the file name, `txxx` shows the temperature. The +first column is the frequency, and the second and third columns are +the values given as follows, respectively, + +$$ +&N_2^{(1)}(\mathbf{q}, \omega) = \frac{1}{N} +\sum_{\lambda'\lambda''} \Delta(-\mathbf{q}+\mathbf{q}'+\mathbf{q}'') +(n_{\lambda'} - n_{\lambda''}) [ \delta( \omega + \omega_{\lambda'} - +\omega_{\lambda''}) - \delta( \omega - \omega_{\lambda'} + +\omega_{\lambda''})], \\ +&N_2^{(2)}(\mathbf{q}, \omega) = \frac{1}{N} +\sum_{\lambda'\lambda''} \Delta(-\mathbf{q}+\mathbf{q}'+\mathbf{q}'') +(n_{\lambda'}+ n_{\lambda''}+1) \delta( \omega - \omega_{\lambda'} - +\omega_{\lambda''}). +$$ + +```bash +% phono3py --fc2 --dim="2 2 2" --pa="F" -c POSCAR-unitcell --mesh="16 16 16" --jdos --ga="0 0 0 8 8 8" --ts=300 +``` + +This is an example of `Si-PBEsol`. + +```{image} Si-JDOS.png +:width: 50% +``` + +## Sampling frequency for distribution functions + +(freq_sampling_option)= +### `--num-freq-points`, `--freq-pitch` (`NUM_FREQUENCY_POINTS`) + +For spectrum like calculations of imaginary part of self energy and +JDOS, number of sampling frequency points is controlled by +`--num-freq-points` or `--freq-pitch`. + +## Mode-Gruneisen parameter from 3rd order force constants + +### `--gruneisen` (`GRUNEISEN = .TRUE.`) + +Mode-Gruneisen-parameters are calculated from fc3. + +Mesh sampling mode: + +```bash +% phono3py --fc3 --fc2 --dim="2 2 2" -v --mesh="16 16 16" -c POSCAR-unitcell --nac --gruneisen +``` + +Band path mode: + +```bash +% phono3py --fc3 --fc2 --dim="2 2 2" -v -c POSCAR-unitcell --nac --gruneisen --band="0 0 0 0 0 1/2" +``` + +## File I/O + +### `--fc2` (`READ_FC2 = .TRUE.`) + +Read 2nd order force constants from `fc2.hdf5`. + +### `--fc3` (`READ_FC3 = .TRUE.`) + +Read 3rd order force constants from `fc3.hdf5`. + +### `--write-gamma` (`WRITE_GAMMA = .TRUE.`) + +Imaginary parts of self energy at harmonic phonon frequencies +$\Gamma_\lambda(\omega_\lambda)$ are written into file in hdf5 +format. The result is written into `kappa-mxxx-gx(-sx-sdx).hdf5` or +`kappa-mxxx-gx-bx(-sx-sdx).hdf5` with `--bi` option. With +`--sigma` and `--sigma-cutoff` options, `-sx` and `--sdx` are +inserted, respectively, in front of `.hdf5`. + +(read_gamma_option)= +### `--read-gamma` (`READ_GAMMA = .TRUE.`) + +Imaginary parts of self energy at harmonic phonon frequencies +$\Gamma_\lambda(\omega_\lambda)$ +are read from `kappa` file in hdf5 format. Initially the usual +result file of `kappa-mxxx(-sx-sdx).hdf5` is searched. Unless it is +found, it tries to read `kappa` file for each grid point, +`kappa-mxxx-gx(-sx-sdx).hdf5`. Then, similarly, +`kappa-mxxx-gx(-sx-sdx).hdf5` not found, +`kappa-mxxx-gx-bx(-sx-sdx).hdf5` files for band indices are searched. + +(write_detailed_gamma_option)= +### `--write-gamma-detail` (`WRITE_GAMMA_DETAIL = .TRUE.`) + +Each q-point triplet contribution to imaginary part of self energy is +written into `gamma_detail-mxxx-gx(-sx-sdx).hdf5` file. Be careful +that this is large data. + +In the output file in hdf5, following keys are used to extract the +detailed information. + +```{table} +| dataset | Array shape | +|-----------------------------|----------------------------------------------------------------------------------------------------------------------------------------| +| gamma_detail for `--ise` | (temperature, sampling frequency point, symmetry reduced set of triplets at a grid point, band1, band2, band3) in THz (without $2\pi$) | +| gamma_detail for `--br` | (temperature, symmetry reduced set of triplets at a grid point, band1, band2, band3) in THz (without $2\pi$) | +| mesh | Numbers of sampling mesh along reciprocal axes. | +| frequency_point for `--ise` | Sampling frequency points in THz (without $2\pi$), i.e., $\omega$ in $\Gamma_\lambda(\omega)$ | +| temperature | (temperature,), Temperatures in K | +| triplet | (symmetry reduced set of triplets at a grid point, 3), Triplets are given by the grid point indices (see below). | +| weight | (symmetry reduced set of triplets at a grid point,), Weight of each triplet to imaginary part of self energy | +``` + +Imaginary part of self energy (linewidth/2) is recovered by the +following script: + +```python +import h5py +import numpy as np + +gd = h5py.File("gamma_detail-mxxx-gx.hdf5") +temp_index = 30 # index of temperature +temperature = gd['temperature'][temp_index] +gamma_tp = gd['gamma_detail'][:].sum(axis=-1).sum(axis=-1) +weight = gd['weight'][:] +gamma = np.dot(weight, gamma_tp[temp_index]) +``` + +For example, for `--br`, this `gamma` gives +$\Gamma_\lambda(\omega_\lambda)$ of the band indices at the grid +point indicated by $\lambda$ at the temperature of index 30. If +any bands are degenerated, those `gamma` in +`kappa-mxxx-gx(-sx-sdx).hdf5` or `gamma-mxxx-gx(-sx-sdx).hdf5` +type file are averaged, but the `gamma` obtained here in this way +are not symmetrized. Apart from this symmetrization, the values must +be equivalent between them. + +To understand each contribution of triptle to imaginary part of self +energy, reading `phonon-mxxx.hdf5` is useful (see +{ref}`write_phonon_option`). For example, +phonon triplets of three phonon scatterings are obtained by + +```python +import h5py +import numpy as np + +gd = h5py.File("gamma_detail-mxxx-gx.hdf5", 'r') +ph = h5py.File("phonon-mxxx.hdf5", 'r') +gp1 = gd['grid_point'][()] +triplets = gd['triplet'][:] # Sets of (gp1, gp2, gp3) where gp1 is fixed +mesh = gd['mesh'][:] +grid_address = ph['grid_address'][:] +q_triplets = grid_address[triplets] / mesh.astype('double') +# Phonons of triplets[2] +phonon_tp = [(ph['frequency'][i], ph['eigenvector'][i]) for i in triplets[2]] +# Fractions of contributions of tripltes at this grid point and temperture index 30 +gamma_sum_over_bands = np.dot(weight, gd['gamma_detail'][30].sum(axis=-1).sum(axis=-1).sum(axis=-1)) +contrib_tp = [gd['gamma_detail'][30, i].sum() / gamma_sum_over_bands for i in range(len(weight))] +np.dot(weight, contrib_tp) # is one +``` + +(write_phonon_option)= +### `--write-phonon` (`WRITE_PHONON = .TRUE.`) + +Phonon frequencies, eigenvectors, and grid point addresses are stored +in `phonon-mxxx.hdf5` file. {ref}`--pa ` and {ref}`--nac +` may be required depending on calculation setting. + +```bash +% phono3py --fc2 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --write-phoonon +``` + +Contents of `phonon-mxxx.hdf5` are watched by: + +```bash +In [1]: import h5py + +In [2]: ph = h5py.File("phonon-m111111.hdf5", 'r') + +In [3]: list(ph) +Out[3]: ['eigenvector', 'frequency', 'grid_address', 'mesh'] + +In [4]: ph['mesh'][:] +Out[4]: array([11, 11, 11], dtype=int32) + +In [5]: ph['grid_address'].shape +Out[5]: (1367, 3) + +In [6]: ph['frequency'].shape +Out[6]: (1367, 6) + +In [7]: ph['eigenvector'].shape +Out[7]: (1367, 6, 6) +``` + +The first axis of `ph['grid_address']`, `ph['frequency']`, and +`ph['eigenvector']` corresponds to the number of q-points where +phonons are calculated. Here the number of phonons may not be equal to +product of mesh numbers ($1367 \neq 11^3$). This is because all +q-points on Brillouin zone boundary are included, i.e., even if +multiple q-points are translationally equivalent, those phonons are +stored separately though these phonons are physically equivalent +within the equations employed in phono3py. Here Brillouin zone is +defined by Wigner–Seitz cell of reciprocal primitive basis +vectors. This is convenient to categorize phonon triplets into Umklapp +and Normal scatterings based on the Brillouin zone. + +(read_phonon_option)= +### `--read-phonon` (`READ_PHONON = .TRUE.`) + +Phonon frequencies, eigenvectors, and grid point addresses are read +from `phonon-mxxx.hdf5` file and the calculation is continued using +these phonon values. This is useful when we want to use fixed phonon +eigenvectors that can be different for degenerate bands when using +different eigenvalue solvers or different CPU +architectures. {ref}`--pa ` and {ref}`--nac ` +may be required depending on calculation setting. + +```bash +% phono3py --fc2 --fc3 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --read-phoonon --br +``` + +(write_read_pp_option)= +### `--write-pp` (`WRITE_PP = .TRUE.`) and `--read-pp` (`READ_PP = .TRUE.`) + +Phonon-phonon (ph-ph) intraction strengths are written to and read +from `pp-mxxx-gx.hdf5`. This works only in the calculation of +lattice thermal conductivity, i.e., usable only with `--br` or +`--lbte`. The stored data are different with and without specifying +`--full-pp` option. In the former case, all the ph-ph interaction +strengths among considered phonon triplets are stored in a simple +manner, but in the later case, only necessary elements to calculate +collisions are stored in a complicated way. In the case of RTA +conductivity calculation, in writing and reading, ph-ph interaction +strength has to be stored in memory, so there is overhead in memory +than usual RTA calculation. + +```bash +% phono3py --fc2 --fc3 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --write-pp --br --gp=1 +``` + +```bash +% phono3py --fc2 --dim="2 2 2" --pa="F" --mesh="11 11 11" -c POSCAR-unitcell --nac --read-pp --br --gp=1 +``` + +(hdf5_compression_option)= +### `--hdf5-compression` (command option only) + +Most of phono3py HDF5 output file is compressed by default with the +`gzip` compression filter. To avoid compression, +`--hdf5-compression=None` has to be set. Other filters (`lzf` or +integer values of 0 to 9) may be used, see h5py +documentation +(http://docs.h5py.org/en/stable/high/dataset.html#filter-pipeline). + +(output_filename_option)= +### `-o` (command option only) + +This modifies default output file names to write. + +Using this option, output file names are slightly modified. For example, +with `-o iso`, a file name `kappa-m191919.hdf5` is changed +to `kappa-m191919.iso.hdf5`. + +This rule is applied to + +- `fc3.hdf5` +- `fc2.hdf5` +- `kappa-xxx.hdf5` +- `phonon-xxx.hdf5` +- `pp-xxx.hdf5` +- `gamma_detail-xxx.hdf5` (write only) + +(input_filename_option)= +### `-i` (command option only) + +This modifies default input file names to read. + +Using this option, input file names are slightly modified. For example, +specifying `-i iso --fc3`, a file name `fc3.iso.hdf5` is read +instead of `fc3.hdf5`. + +This rule is applied to + +- `fc3.hdf5` +- `fc2.hdf5` +- `kappa-xxx.hdf5` +- `phonon-xxx.hdf5` +- `pp-xxx.hdf5` + +### `--io` (command option only) + +This modifies default input and output file names. + +This is equivalent to setting `-i` and `-o` simultaneously. diff --git a/doc/crystal.md b/doc/crystal.md new file mode 100644 index 00000000..6cdd00c6 --- /dev/null +++ b/doc/crystal.md @@ -0,0 +1,95 @@ +(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: + +```bash +% 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) + + ```bash + % 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: + + ```bash + % runcry17 supercell-001.d12 + ``` + +3) Collect forces in `FORCES_FC3` and `FORCES_FC2`: + + ```bash + % 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`: + + ```bash + % 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: + + ```bash + % 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`! diff --git a/doc/cutoff-pair.md b/doc/cutoff-pair.md new file mode 100644 index 00000000..0e1e9ac5 --- /dev/null +++ b/doc/cutoff-pair.md @@ -0,0 +1,429 @@ +(command_cutoff_pair)= +# Force constants calculation with cutoff pair-distance + +Here the detail of the command option {ref}`--cutoff_pair ` 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, +$\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 +($\Delta\mathbf{r}_1$ and $\Delta\mathbf{r}_2$) but one of +them is included in a force given by the force calculator such as +$\mathbf{F}_3=-\frac{\partial V}{\partial\mathbf{r}_3}$. The +cutoff distance $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 $\{(\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 $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: + +```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 51. +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`). + +```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 ` may be recommended +when the number of force files is large. + +```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: + +```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. + +```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. + +```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. + +```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: + +```python +if data2['included']: +``` + +by + +```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"`. + + +## Tests + +### 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`. + +```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 + +```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 +``` \ No newline at end of file diff --git a/doc/direct-solution.md b/doc/direct-solution.md new file mode 100644 index 00000000..5e9537e2 --- /dev/null +++ b/doc/direct-solution.md @@ -0,0 +1,345 @@ +(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 `). + +```{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 `. An example to run with the direct solution of +LBTE for `example/Si-PBEsol` is as follows: + +```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 $(\text{number of +irreducible grid points} \times \text{number of bands} \times 3)^2$ +for the 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 ` together with +`--lbte` option, estimated memory space needed for storing collision +matrix is presented. An example for `example/Si-PBEsol` is as follows: + +``` +% 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 `, 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 +` 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 ` and {ref}`--read-gamma ` 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 ` when reading. + +The summary of the procedure is as follows: + +1. Running at each grid point with {ref}`--gp ` (or + {ref}`--ga `) 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`: + +```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`): + +``` +% 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: + +```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 `, the tetrahedron method or +smearing approach with {ref}`--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 ` and {ref}`--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 ` (or + {ref}`--ga `) 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 ` (or + {ref}`--ga `) 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` + +```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 + +```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`) + +```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 + +```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: + +```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 $10^{-4}$ and $10^{-16}$, which +is a good sign. The values whose absolute values are smaller than +$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. diff --git a/doc/examples.md b/doc/examples.md new file mode 100644 index 00000000..13d38d85 --- /dev/null +++ b/doc/examples.md @@ -0,0 +1,16 @@ +(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. diff --git a/doc/external-tools.md b/doc/external-tools.md new file mode 100644 index 00000000..173052c6 --- /dev/null +++ b/doc/external-tools.md @@ -0,0 +1,29 @@ +(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/ diff --git a/doc/hdf5_howto.md b/doc/hdf5_howto.md new file mode 100644 index 00000000..4041ab27 --- /dev/null +++ b/doc/hdf5_howto.md @@ -0,0 +1,337 @@ +(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. + +```bash +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, $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 ($\tau_\lambda=1/2\Gamma_\lambda(\omega_\lambda)$) may +be estimated from `gamma`. $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 +$\Gamma$ point. The below is the copy-and-paste from the +previous section to show how to obtain phonon lifetime in pico +second: + +```bash +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 $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, $\nabla_\mathbf{q}\omega_\lambda$. The +physical unit is $\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 + +$$ +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 (${}^*\mathbf{k}$): + +$$ +\sum_{\mathbf{q} \in {}^*\mathbf{k}} \kappa_{\mathbf{q}j}. +$$ + +The sum of this over ${}^*\mathbf{k}$ corresponding to +irreducible q-points divided by number of grid points gives +$\kappa$ ({ref}`output_kappa`), e.g.,: + +```python +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 ${}^*\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 +(${}^*\mathbf{k}$) for each irreducible q-point and phonon band +($j$): + +$$ +\sum_{\mathbf{q} \in {}^*\mathbf{k}} \mathbf{v}_{\mathbf{q}j} \otimes +\mathbf{v}_{\mathbf{q}j}. +$$ + +The physical unit is +$\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 $\text{eV}^2, +$P_{\mathbf{q}j}$: + +$$ +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 + +$$ +\kappa_\lambda = \frac{1}{V_0} C_\lambda \mathbf{v}_\lambda \otimes +\mathbf{v}_\lambda \tau_\lambda^{\mathrm{SMRT}}. +$$ + +For example, $\kappa_{\lambda,{xx}}$ is calculated by: + +```badh +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: + +```bash +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`, + +```bash +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 +``` diff --git a/doc/index.md b/doc/index.md new file mode 100644 index 00000000..f7d944a2 --- /dev/null +++ b/doc/index.md @@ -0,0 +1,85 @@ +# 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`) +- 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 `, + {ref}`QE (pw) `, {ref}`CRYSTAL `, + {ref}`TURBOMOLE `, and Abinit + are built in (see {ref}`calculator_interfaces`). +- API is prepared to operate phono3py from Python ([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`). + +```{image} Si-kaccum.png +:width: 20% +``` +```{image} Si-kaccum-MFP.png +:width: 20% +``` +```{image} Si-kdeplot.png +:width: 22% +``` + + +## Documentation + +```{toctree} +:maxdepth: 1 +install +examples +Interfaces to calculators (VASP, QE, CRYSTAL, Abinit, TURBOMOLE) +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/) diff --git a/doc/install.md b/doc/install.md new file mode 100644 index 00000000..5a7c908a --- /dev/null +++ b/doc/install.md @@ -0,0 +1,333 @@ +(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: + +```bash + % 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: + +```bash + % 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: + +```python + 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: + +```bash + % 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 by + +```bash + % 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 + +```bash + 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: + +```python + 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: + +```python + 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`): + +```bash + % 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. + +```bash + % 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: + + ```bash + % wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh -O ~/miniconda.sh + ``` + + For macOS: + + ```bash + % wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O ~/miniconda.sh + ``` + + The installation is made by + + ```bash + % 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 + + ```bash + % conda init + ``` + + `` 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. + + ```bash + % 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 + + ```bash + % 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 + + ```bash + export HDF5_USE_FILE_LOCKING=FALSE + ``` + + Install the latest phonopy and phono3py from github sources: + + ```bash + % 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. diff --git a/doc/interfaces.md b/doc/interfaces.md new file mode 100644 index 00000000..5f318776 --- /dev/null +++ b/doc/interfaces.md @@ -0,0 +1,74 @@ +(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. + +```{table} +| calculator | unit-cell | FORCES_FC3 | phono3py_disp.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_FC3`, `FORCES_FC2`, and `phono3py_disp.yaml` have the same physical units. + +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 $\text{eV}/\text{Angstrom}^2$ and +$\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:: + +```{table} +| calculator| filename | +|-----------|-------------| +| 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:: + +```{table} +| calculator | distance | +|------------|----------------| +| VASP | 0.03 Angstrom | +| QE | 0.06 au (bohr) | +| CRYSTAL | 0.03 Angstrom | +| Abinit | 0.06 au (bohr) | +| TURBOMOLE | 0.06 au (bohr) | +``` \ No newline at end of file diff --git a/doc/output-files.md b/doc/output-files.md new file mode 100644 index 00000000..153e182d --- /dev/null +++ b/doc/output-files.md @@ -0,0 +1,132 @@ +(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 +$\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 $\Phi_{\alpha\beta\gamma}(l\kappa, l'\kappa', +l''\kappa'')$. The first three `num_atom` are the atom indices in supercell +corresponding to $l\kappa$, $l'\kappa'$, +$l''\kappa''$, respectively. The last three elements are the Cartesian +coordinates corresponding to $\alpha$, $\beta$, +$\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 `. 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, + +$$ +(\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 +$\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 $\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 +($\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 +$\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` diff --git a/doc/qe.md b/doc/qe.md new file mode 100644 index 00000000..db34cce7 --- /dev/null +++ b/doc/qe.md @@ -0,0 +1,79 @@ +(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, + +```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 + + ```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: + + ```bash + % phono3py --cf3 disp-00001/Si-supercell.out disp-00002/Si-supercell.out ... + ``` + + or in recent bash or zsh: + + ```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: + + ```bash + % phono3py --sym-fc + ``` + + where `--sym-fc` symmetrizes fc3 and fc2. + +5) Calculate lattice thermal conductivity, e.g., by: + + ```bash + % phono3py --mesh="11 11 11" --fc3 --fc2 --br + ``` diff --git a/doc/tips.md b/doc/tips.md new file mode 100644 index 00000000..39760005 --- /dev/null +++ b/doc/tips.md @@ -0,0 +1,111 @@ +(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 `, default option) and a smearing method ({ref}`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 $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. + +```{image} Si-convergence.png +:width: 25% +``` + +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 +$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 `. 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 +$\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 ` 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 ` 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 ` 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. diff --git a/doc/turbomole.md b/doc/turbomole.md new file mode 100644 index 00000000..fe2af512 --- /dev/null +++ b/doc/turbomole.md @@ -0,0 +1,85 @@ +(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: + +```bash +% 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) + + ```bash + % 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: + + ```bash + % riper > supercell-00001.out + ``` + +3) Collect forces in `FORCES_FC3` and `FORCES_FC2`:: + + ```bash + % 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`: + + ```bash + % 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: + + ```bash + % 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`! diff --git a/doc/vasp.md b/doc/vasp.md new file mode 100644 index 00000000..30dcb0b0 --- /dev/null +++ b/doc/vasp.md @@ -0,0 +1,112 @@ +(vasp_interface)= +# VASP & phono3py calculation + +(vasp_workflow)= +## Workflow + +1. Create POSCARs with displacements + + This is the same way as usual phonopy: + + ```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: + + ```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. + +2. 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. + +3. Collect `vasprun.xml`'s + + When VASP is used as the force calculator, force sets to calculate + fc3 and fc2 are created as follows. + + ```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: + + ```bash + % phono3py --cf2 disp_fc2-{00001..00002}/vasprun.xml + ``` + `phono3py_displ.yaml` is necessary in this case and `FORCES_FC2` is + created. + +4. Create `fc2.hdf` and `fc3.hdf` + + ```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. + +5. Thermal conductivity calculation + + An example of thermal conductivity calculation is: + + ``` + % 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`: + + ``` + % 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: + + ``` + % phono3py --mesh="11 11 11" --br --write-gamma --gp="[grid ponit(s)]" + ``` + + After finishing all distributed calculations, run with + `--read-gamma` option: + + ``` + % 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. diff --git a/doc/workload-distribution.md b/doc/workload-distribution.md new file mode 100644 index 00000000..48224467 --- /dev/null +++ b/doc/workload-distribution.md @@ -0,0 +1,142 @@ +(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 `), over phonon bands ({ref}`--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 `. Once all data are +obtained, those data are collected by {ref}`--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: + +```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 ` + +```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`. + +```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: + +```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 ` and {ref}`-i ` +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: + +```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. + +```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`, + +```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, + +```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): + +```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 +``` \ No newline at end of file