Version v2.0.0 of ShengBTE has been released and can be obtained from the Downloads section. It is the largest release since v1.0.0: a thoroughly modernized version of the code, written in standard Fortran 2008 and built with CMake, with an optional test suite. It adds the contribution of coherences to the thermal conductivity, following the Wigner formulation of thermal transport (flag wigner), and custom isotopic compositions through a new ISOTOPES file, which replaces the masses and gfactors variables of CONTROL. Three-phonon processes are computed much faster and with much less memory, and several bugs have been fixed, including one that made the results of the iterative solver depend on the number of MPI processes.
The full list of changes, from the annotation of the v2.0.0 tag in the Bitbucket repository, follows.
New features:
Coherence contribution to the thermal conductivity: the new flag wigner adds the contribution of the coherences between phonon modes, from the Wigner formulation of thermal transport [M. Simoncelli, N. Marzari and F. Mauri, Nat. Phys. 15, 809 (2019); Phys. Rev. X 12, 041011 (2022)], written to BTE.KappaTensorVsT_COH. It describes the wave-like tunnelling of heat between phonon branches and matters in complex crystals with low or glass-like thermal conductivity. For La2Zr2O7 (22 atoms) on a 16x16x16 grid, the populations and coherences conductivities agree with Fig. 3 of the Phys. Rev. X paper within 1% and 2%, respectively.
The velocity operator is computed at the irreducible q points, split among the MPI processes, from the derivative of the dynamical matrix that both input formats already build, and transformed from the step-like phase convention used by ShengBTE to the smooth one, the only one that gives a conductivity independent of the choice of unit cell.
The linewidths are the total zeroth-order scattering rates. Pairs of modes closer in frequency than the new parameter wigner_degeneracy (1e-4 rad/ps by default) are treated as degenerate and left to the usual conductivity, as are modes with imaginary frequencies and the Gamma point.
Custom isotopic compositions: an optional ISOTOPES file, read when isotopes=.true. and autoisotopes=.false., gives the mass and concentration of each isotope of each element. The concentrations of each element are normalized; every element must appear, and unknown elements, malformed lines, negative concentrations and concentrations that add up to zero are reported as errors.
Bug fixes:
Multi-process iterative solver: after each iteration, only the first MPI process restored the symmetry of the solution, but all of them used it in the next iteration. Results therefore depended on the number of MPI processes: up to 6e-7 in the conductivity and 42% in some rates of BTE.w_final for Test-VASP. All processes now agree with a serial run to round-off.
Isotopic scattering: fixed a data race between OpenMP threads in the iterative solver, which could lose contributions and make results change from run to run.
Results that changed from run to run: the sums over q points of the conductivities, cumulative conductivities, heat capacity and total Grüneisen parameter added the contributions of the OpenMP threads in the order in which they finished, so their last digits could change between identical runs. They are now summed in a fixed order, and the results are the same for any number of OpenMP threads.
MPI reductions: the send and receive buffers of MPI_Reduce were aliased on non-root processes.
Dense q-point grids: fixed an integer overflow in the normalization of the phase space for grids with more than 46340 points.
Error handling: error paths no longer call MPI_Finalize and carry on; input errors stop all processes with a clear message.
Unknown elements: element names missing from the table of natural isotopes silently gave a zero mass and an undefined g factor; they are now reported as an error.
Small-grain conductivity: it was the only conductivity tensor written without being symmetrized. With degenerate modes its diagonal elements could differ by a few percent in a cubic crystal (2.4% for La2Zr2O7 on a 4x4x4 grid); it is now symmetrized like the others.
Progress of the matrix elements: the percentage was printed with backspaces and without a newline, so it never appeared when the output went to a file or through an MPI launcher.
Spurious floating-point note: runs stopped with onlyharmonic=.true. ended with "The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO", raised by the boundary rates of modes with zero velocity and, harmlessly, inside LAPACK's singular value decomposition. Neither raises the flag any more.
Improvements:
Much faster three-phonon matrix elements: for each incoming phonon and q', the contraction of the third-order force constants with the phase factors and its eigenvector is computed once and shared by all the allowed processes (a single process still uses the direct formula, which is cheaper). The matrix elements are computed 1.75 and 2.4 times faster for the InAs and GaAs test cases, and 127 times faster for La2Zr2O7 on a 4x4x4 grid; on its 16x16x16 grid (2.2e9 processes), they took 82 minutes with 32 MPI processes x 8 OpenMP threads (dual AMD EPYC 7713).
Faster search for three-phonon processes: the width of the Gaussian of a process is bounded by those of its second and third phonons, so for each second phonon only the bands of the third one with frequencies in a window around energy conservation are tested, found by binary search. The processes found and the results do not change; a run of La2Zr2O7 on an 8x8x8 grid with onlyharmonic=.true., which counts the processes and computes their weighted phase space, takes 39% less time.
Much less memory for the phase factors: they are stored per distinct lattice vector instead of per triplet of atoms and q point: 2.9 MB instead of 29.5 GB per MPI process for La2Zr2O7 on a 16x16x16 grid, which makes calculations with many MPI processes per node possible.
Lower peak memory: the third-order force constants, the Grüneisen parameters and the per-mode conductivity are no longer kept for the whole run.
Faster setup on dense grids: the irreducible q points are found in linear instead of quadratic time.
Diagnostics:
the version of ShengBTE, on the first line of the output;
a warning about modes with imaginary frequencies (other than the acoustic modes at Gamma), whose coherences are left out of the coherence conductivity;
a warning about overdamped modes (linewidth larger than frequency), for which the Wigner formulation is not valid;
the relative asymmetry removed by symmetrization from the small-grain, RTA and converged conductivities, as a measure of how well the grid and the eigenvectors respect the symmetry of the crystal;
the progress of the computation of the matrix elements, in complete lines every 10% with the elapsed time and an estimate of the remaining time, preceded by the minimum, mean and maximum number of processes per MPI process.
Build system:
CMake (3.20 or later) replaces the Makefile and arch.make, with automatic detection of MPI, LAPACK/BLAS, spglib and OpenMP; spglib can also be downloaded and built (-DSHENGBTE_FETCH_SPGLIB=ON).
Standard Fortran 2008, checked by the compiler, and the mpi_f08 MPI bindings.
Version: set only in the project() command of the top-level CMakeLists.txt.
Code structure:
Derived types instead of global variables: all data travel as arguments (crystal, q-point grid, symmetry, parameters, phonons, scattering processes, ...), and no module keeps global mutable state.
Stages and output: the main program is a short sequence of stages (shengbte_workflow), and every output file is written by one module (shengbte_output).
Pluggable elastic scattering: isotopic scattering has its own module and reaches the solver only through a generic description of elastic processes, so other models can be added without touching the solver.
No duplicated variants: absorption and emission processes, and bulk and nanowire solutions, share the same code.
Consistent naming and formatting: a shengbte_ module prefix and descriptive snake_case names, enforced with fprettify; see docs/STYLE.md.
Testing:
Opt-in test suite (-DSHENGBTE_BUILD_TESTS=ON, then ctest), with:
unit tests (test-drive), including the ISOTOPES file, the search for and the matrix elements of three-phonon processes, and the velocity operator and coherence conductivity (against finite differences of the square root of the dynamical matrix, a direct evaluation, the sum over the whole grid, and a primitive cell against a supercell);
tests that invalid CONTROL and ISOTOPES files are rejected with a clear message, and that the version is reported;
regression tests that compare each Test-* case, including the new Test-wigner, with its Reference/ directory, with tolerances that allow for a different machine or toolchain;
consistency tests across MPI x OpenMP layouts and between repeated runs;
a formatting check.
References: the Reference/ directories were generated with this version and checked against the outputs of the public version built with the same toolchain; the public version is no longer built by the test suite.
Documentation:
README: building with CMake, running the tests, the ISOTOPES file, the wigner flag and wigner_degeneracy parameter, and the new output file.
AUTHORS: the authors, with their e-mail addresses and the years of their contributions; the copyright notices of the sources now refer to "The ShengBTE Authors".
docs/STYLE.md: coding conventions and performance guidelines for contributors.
nthreads: now documented correctly. It only reports a value; the number of threads is set with OMP_NUM_THREADS.
Compatibility:
Inputs: masses and gfactors can no longer be set in CONTROL; files that still set them are rejected with a message pointing to ISOTOPES, and element names must be chemical symbols unless the composition is read from ISOTOPES. The flag wigner and the parameter wigner_degeneracy are new and optional. Everything else in CONTROL is unchanged.
Outputs: the names and formats of all output files are unchanged; BTE.KappaTensorVsT_COH is new, and the output now starts with the version and includes the diagnostic messages described above.
Results: the same as those of the public version up to round-off (below 1e-9 relative in the test cases, in quantities that depend on the three-phonon matrix elements), except for the small-grain conductivity, which is now symmetrized, and for iterative solutions with several MPI processes, which change slightly because of the bug fix above. They no longer depend on the number of OpenMP threads.
Requirements: a Fortran 2008 compiler, an MPI library providing the mpi_f08 module, and CMake 3.20 or later.
The latest release of the 1.x series, v1.5.1, will remain available from the Downloads section, and any older release can of course be obtained from the git repository, where each one is tagged (v0.9.0 to v1.5.1).