Euphonic: inelastic neutron scattering simulations from force constants and visualisation tools for phonon properties

Euphonic has been developed following software development best practice: continuous integration processes test the functionality and validate numerical results as features are added, and the Python API has been explicitly designed to make it easy to use Euphonic with other projects. Euphonic is open-source under the GNU General Public License v3 – the source code is is freely available on Github, and packages are distributed by PyPI and Conda-forge for the pip and conda package managers. Links to the source code and online documentation, including installation instructions and how to use Euphonic, are available from the DOI landing page https://doi.org/10.5286/SOFTWARE/EUPHONIC [20].

Throughout this paper a variety of simulated and experimental datasets will be used to compare Euphonic outputs to both other software tools and experimental data. The chosen materials are Nb, Al, Si, quartz and \ceLa2Zr2O7, and will be described and used in each section as appropriate. These datasets have been chosen to give a variety of unit cell sizes and crystal symmetries, and include both polar (quartz) and non-polar materials. Simulated force constants have been produced with either CASTEP or VASP [36] plus Phonopy in order to thoroughly test the features of Euphonic. Details of the force constant calculations for each material and where to obtain the results of the simulations and the neutron experimental data can be found in the supplementary material.

Vibrational excitations — phonons — in crystalline materials are described within the theory of lattice dynamics [14]. Displacement of a single atom in an infinite periodic crystal gives rise to forces on every other atom in the crystal, which in the harmonic approximation are linearly related to the displacement by the force constant matrix

where $\Phi$ is the force constant matrix, the unit cell containing the displaced atom is labelled $0$, $l$ runs over unit cells in the crystal, $\kappa$ runs over atoms in the unit cell, $\alpha$ runs over the Cartesian directions, $u_{\kappa{\alpha}l}$ is the displacement of atom $\kappa$ in cell $l$ in direction $\alpha$ from its equilibrium position, and $E$ is the total lattice energy. These force constants decay rapidly with distance according to a power-law $|R|^{-n}$ where $n=5-7$ in non-polar crystals. This allows for a truncation of $\Phi(0,l)$ at some radius $R_{c}$ beyond which the residual force constants may be neglected. $R_{c}\approx{8-10}$ Å  is sufficient for almost all practical cases.

Substituting a plane-wave solution into the equation of motion yields an eigenvalue equation

where $D_{\alpha{\alpha}^{\prime}}^{\kappa{\kappa}^{\prime}}$($\mathbf{q}$) is the dynamical matrix at wave vector $\mathbf{q}$, the eigenvalues ${\omega}_{\mathbf{q}\nu}^{2}$ are the square of the phonon frequencies of mode $\nu$ at $\mathbf{q}$, and the eigenvectors ${\mathbf{e}}_{\mathbf{q}\nu\kappa}$ are the phonon polarisation vectors at $\mathbf{q}$ of mode $\nu$ for atom $\kappa$. The dynamical matrix can be calculated as the mass-weighted Fourier transform of the force constant matrix

where $M_{\kappa}$ is the mass of atom $\kappa$, $\mathbf{R}_{l}$ is the vector from the origin to the $l^{\textrm{th}}$ unit cell, and $l$ includes unit cells with $|\mathbf{R}_{l}|<R_{c}$.

The coherent one-phonon scattering structure factors can be calculated at momentum transfer $\mathbf{Q}$ of scattered neutrons [15, 49]

where $\mathbf{q}$ is the reduced wave vector in the first Brillouin zone (i.e. $\mathbf{Q}=\mathbf{q}+\bm{\tau}$ where $\bm{\tau}$ is a reciprocal lattice vector), $b_{\kappa}$ is the coherent neutron scattering length of atom $\kappa$, $\mathbf{r}_{\kappa}$ is the vector from the origin to atom $\kappa$ within the unit cell, $M_{\kappa}$ is the mass of atom $\kappa$, ${\mathbf{e}}_{\mathbf{q}\nu\kappa}$ is the phonon polarisation vector at $\mathbf{q}$ of mode $\nu$ for atom $\kappa$ and ${\omega}_{\mathbf{q}\nu}$ is the phonon frequency of mode $\nu$ at $\mathbf{q}$. The term $\exp\left(-W_{\kappa}\right)$ is the anisotropic Debye–Waller factor for atom $\kappa$, where the exponents can be written

where the sum is over wave vectors and modes in the first Brillouin Zone (BZ), $N_{q^{\prime}}$ is the number of $\mathbf{q}$-points in the sum, $T$ is the temperature, $\alpha$ and $\beta$ run over the Cartesian directions, $\hbar$ is the reduced Planck constant and $k_{B}$ is the Boltzmann constant. The Debye–Waller factor has been written in this form to make it explicit that the expensive computation of the set of $3\times 3$ matrices $W_{\kappa}^{\alpha\beta}$ need only be performed once over an appropriately fine grid in the first Brillouin Zone, leaving computation of the Debye–Waller factor for an arbitrary $\mathbf{Q}$ as the fast evaluation of a quadratic form for each atom, Eq. 8.

From the one-phonon structure factors the neutron dynamical structure factor $S(\mathbf{Q},\omega)$ can be calculated

where the upper and lower signs refer to phonon creation and annihilation respectively and $n_{\mathbf{q}\nu}$ is the Bose population function

Finally, the neutron scattering cross-section per unit cell in term of $S(\mathbf{Q},\omega)$ is

where $k_{i}$ and $k_{f}$ are the incident and scattered neutron wave vectors.

In a typical Horace workflow, the scattering intensities at millions of $\mathbf{Q}$-points can be combined to produce multidimensional plots of measured data. Prior to the availability of Euphonic, the workflow for generating such plots was limited by the need to read pre-computed phonon frequencies and eigenvectors from files produced by other programs. The bottleneck here becomes the activity of reading and writing to disk. For a system of $N$ atoms per unit cell, each $\mathbf{Q}$-point has a set of eigenvectors whose storage requirements of $18N^{2}$ floating-point numbers can become impractically large. For example, a 22-atom unit cell with a modest 25,000 $\mathbf{Q}$-points would equate to a 5 GB text file in CASTEP .phonon format. In a typical cluster environment running the CASTEP phonons tool, more than 85% of run time is spent writing this file to disk (see Table 1 for specific $N^{2}$ timing examples). Furthermore, on a shared cluster resource this can be exacerbated by local network capacity and the activities of other users.

The requirement of efficiency for Euphonic has driven the decision to implement Fourier interpolation of phonon frequencies and eigenvectors directly from force constants. This allows for calculation of data for each $\mathbf{Q}$-point on demand, removing any need for file-based data transfer from other codes. Accordingly, the representation of force constants as a class and associated methods forms the core of the Euphonic API, illustrated in Fig. 2 and described for reference below.

ForceConstants objects can be instantiated from Python data objects such as Numpy [30] arrays, but would typically be created from the data outputs of an external modelling code – currently CASTEP .castep_bin output and Phonopy phonopy.yaml, FORCE_CONSTANTS and force_constants.hdf5 output are supported. Through Phonopy, Euphonic force constants can be obtained using a wide variety of atomistic codes such as VASP [36], Abinit [27] and Quantum Espresso [24], making Euphonic accessible to a large portion of the materials modelling community. From the force constants, phonon frequencies and eigenvectors may be calculated at arbitrary $\mathbf{Q}$-points using the methods described in Section 2.1.

QpointPhononModes represents phonon mode data: the $\mathbf{Q}$-points, phonon frequencies and eigenvectors. QpointPhononModes can also be instantiated from external data files (currently, these are CASTEP .phonon and Phonopy mesh/band/qpoints files). The Pint Python library [29] is extensively used in Euphonic to represent dimensioned data as a Quantity with both a magnitude and a unit. This makes the units explicit, and facilitates conversion to the end user’s preferred units, usually one of meV, cm^-1 or THz, in the case of phonon frequencies. No particular $\mathbf{Q}$-point sampling is enforced; while, for example, a Monkhorst–Pack mesh is recommended for traditional DOS plotting, phonon modes can also be calculated along a high-symmetry path or sampled at random points depending on the use case. From the phonon modes, Euphonic can compute quantities such as the mode-resolved structure factors, the Debye–Waller exponent and total, partial and neutron-weighted DOS.

StructureFactor is derived from the phonon modes. This object contains the neutron structure factors resolved by $\mathbf{Q}$-point and phonon mode index ($\nu$) and also the $\mathbf{Q}$-points and frequencies. This can be binned in energy to produce a two-dimensional (2D) $S(\mathbf{Q},\omega)$ plot, or averaged over the contained $\mathbf{Q}$-points to produce a one-dimensional (1D) $S(\omega)$ spectrum.

DebyeWaller represents the temperature dependent Debye–Waller factor; specifically, it contains the set of $3\times 3$ matrices, one per atom, $W_{\kappa}^{\alpha\beta}$ defined in Eq. 9. A DebyeWaller object is typically pre-computed over a uniform $\mathbf{Q}$-point mesh and then applied during a structure factor calculation to perform the fast calculation of the atomic Debye–Waller exponent $W_{\kappa}$ via Eq. 8 over whichever (potentially large) set of $\mathbf{Q}$-points are appropriate to the task at hand.

Spectrum2D, Spectrum1D and Spectrum1DCollection are classes representing generic spectrum objects that can be used for various purposes, such as representing the DOS with Spectrum1D, or a $S(\mathbf{Q},\omega)$ or $S(|\mathbf{Q}{}|,\omega)$ intensity map with Spectrum2D. Band structure and partial DOS data are represented with Spectrum1DCollection, which ensures consistent bins are used and allows individual lines to be tagged with metadata. The plotting tools in turn work with the generic spectrum objects to produce plots of phonon band structure, DOS, and intensity maps.

Crystal is a simple class containing the crystal structure information: the cell vectors, atom positions, species and masses. The above ForceConstants, QpointPhononModes, StructureFactor and DebyeWaller classes all contain an instance of this crystal class as an attribute, to ensure the data in each class remains complete and unambiguous.

A widely used software application for analysis and visualisation of multidimensional time-of-flight inelastic neutron scattering data from single crystal experiments is Horace [17]. In addition to handling and plotting the data, it also allows simulation and fitting of these datasets with user-created models of the scattering function. A Matlab add-on has been developed, Horace-Euphonic-Interface [19] which provides interface functions that allow Euphonic to be used to simulate datasets directly in Horace. The way this works is illustrated in Fig. 1. First a user sets up a model with Horace-Euphonic-Interface, giving the path to the force constants file or folder, and adding other optional parameters such as the sample temperature or the Debye–Waller grid size. The user then calls a Horace simulation function with the dataset to be simulated and the model they have just created. Horace will automatically provide the $\mathbf{Q}$-points to be simulated to Horace-Euphonic-Interface, which then calls Euphonic to calculate the structure factors and phonon frequencies at those $\mathbf{Q}$-points. Horace-Euphonic-Interface then converts the output from Euphonic into the required form for Horace to create the simulated dataset.

This is a significant improvement over previous workflows to simulate scattering from phonons, which required users to program their own functions to calculate the structure factors from ab initio calculations. These would not usually include Fourier interpolation of the force constants, and hence were restricted to the $\mathbf{Q}$-points output by the ab initio calculations. These bespoke scripts were also not typically optimised for fast computation. Horace-Euphonic-Interface has allowed users to easily simulate on exactly the same axes and in the same software as the experimental data, making fitting of scaling factors and quantitative comparisons much more convenient. The performance of Euphonic has also made application of the Monte Carlo based instrumental resolution convolution method [43] available in Horace tractable with phonons for the first time. Examples of data simulated and fitted with Euphonic and Horace are shown in Section 6.3. Horace-Euphonic-Interface is open source and distributed as a Matlab Toolbox file on Github and the Matlab File Exchange. Links to the source code and online documentation are available from the DOI landing page https://doi.org/10.5286/SOFTWARE/HORACEEUPHONICINTERFACE [19].

AbINS is a code which simulates powder-averaged INS spectra in an analytic incoherent approximation, and resides in the Mantid framework used for experimental data reduction and analysis [16, 3, 2]. AbINS has recently been updated to make use of Euphonic; as of Mantid version 6.3, it is possible to select force constants data in CASTEP or Phonopy format as an input file. End users do not need to know anything about Euphonic or change their workflow, except to ensure that force constants data is present in their .castep_bin or phonopy.yaml files. A single-parameter cutoff distance (as defined by [40], and hidden from the user interface) is used to determine a default Monkhorst–Pack mesh for the given crystal structure; eigenvalues and eigenvectors are computed using the Euphonic Python API and passed on to the usual incoherent inelastic scattering computation [16]. AbINS and Horace have different target users: with Euphonic as a common dependency, it becomes possible for these neutron-scattering simulation codes to develop overlapping feature sets while sharing implementation work.

All profiling in this section has been completed on the STFC Scientific Computing Department’s SCARF Cluster using the SCARF 18 hardware, which contains 2 Intel Gold 6126 processors per node (24 cores per node). Euphonic 0.6.1 and CASTEP 19.1 were used. Both the Euphonic C extension and the CASTEP phonons tool have been compiled with the Intel 2018.3 compiler, OpenMPI 3.1.1 and linked against Intel MKL 2018.3. The profiling results can be obtained at https://github.com/pace-neutrons/euphonic-performance.

This section illustrates some of the main features of Euphonic, using a variety of samples for both single crystals and powders, from a range of modelling codes. An overview of each of the command-line tools in Euphonic is shown in Table 5, with example figures referenced for each. Many of these tools allow sampling parameters (e.g. energy bins, broadening) and appearance options (e.g. unit conversions, axis labels and styling) to be easily specified via command-line arguments.

For custom plots it may be necessary to write a Python script. For example, Fig. 6 shows the neutron dynamical structure factor sampled along arbitrary $\mathbf{Q}$ directions. A sample script to achieve this kind of result is shown in Fig. 9.

Here, prior experimental measurements are compared with newly simulated spectra [18]. Measurements of powdered elemental samples of Al (at 5 K with 60 meV incident energy and the Gd monochromating chopper running at 200 Hz) and Si (at 300 K with 80 meV incident energy with the “sloppy”¹¹1The term “sloppy” here means the chopper has relative wide openings to maximize flux over a wide range of energies (particularly low energies ¡ 100meV), in contrast to other choppers which are optimized for sharp resolution at specific (high) energies but at the cost of lower transmission and flux. monochromating chopper at 250 Hz) were recorded on the MARI instrument at ISIS, the data were reduced using Mantid and binned to 2D $S(|\mathbf{Q}{}|,\omega)$ maps using MSLICE in Mantid (Fig. 10(a) and Fig. 11(a)). For the Al calculations, VASP 5.4.4 was used, and the force constants were obtained by finite displacements using Phonopy [36, 37, 50]. For Si, CASTEP was used, obtaining force constants by density-functional perturbation theory (DFPT) [13, 45]. For more details, see the supplementary material.

The resulting castep.bin and phonopy.yaml files were used with the euphonic-powder-map command-line tool to generate numerically sampled maps of $S_{\mathrm{coh}}(|\mathbf{Q}{}|{},\omega)$. (Details of the parameters are given in the supplementary material). These results are plotted in Fig. 10(b) and Fig. 11(b) and include the main inelastic scattering features in the positive region of the experimental measurements. It is easy to see in both the experimental and simulated data how spherically averaged periodic features collide to create broad regions of higher intensity. In the case of the 300 K Si data, both the experimental and simulated results have significant intensity in the negative energy transfer region, whereas for the 5 K Al data the intensity is suppressed in this region. This arises from the Bose population factor suppressing the low temperature scattering. The high intensities seen around zero energy transfer in the experimental measurements in Fig. 10(a) and Fig. 11(a) are due to scattering from Bragg peaks which is not currently modelled by Euphonic, hence they are absent from the corresponding simulations in Fig. 10(b). and Fig. 11(b). The high intensity feature in Fig. 11(a) at 80 meV is an artefact due to a later pulse of neutrons which passes through the chopper system, and is not due to scattering from the sample, so is not reproduced in the simulated data in Fig. 11(b).

Simulations created from pre-existing quartz force constants calculations have been compared with newly collected experimental data [18]. Original quartz force constants calculations performed with CASTEP 6.1 have been re-processed for this work with CASTEP 19.1; more details are in the supplementary material. For the experimental measurements, a large single crystal of natural quartz was aligned with the c-axis vertical on an Al plate, secured in place with Al wire. It was cooled to 10 K using a closed cycle refrigerator. The MERLIN ‘G’ chopper at 350 Hz was phased to select 45 meV incident energy neutrons. The sample was rotated over 180 ^∘ in 0.5 ^∘ steps. The data for each individual angle were reduced using Mantid and then combined using Horace.

Several cuts through experimental and simulated data for quartz are shown in Fig. 12. The neutron dynamical structure factor has been computed with Euphonic, with the instrumental resolution accounted for by the Monte Carlo resolution convolution method [43] implemented in Horace. In addition, 1D cuts are shown in Fig. 13, with the simulated scattering likewise convolved with the instrumental resolution function. The only adjustable parameter in the comparison is a global scaling factor which has been determined by setting the integrated areas of the experimental and simulated data in Fig. 13(a) between 12 and 37 meV to be equal, which has then been applied to all 1D and 2D cuts. Aside from this (arbitrary) choice of intensity scale, no adjustable parameters have been used in any part of the calculation. There are small disagreements, for example the phonon frequencies do not match up perfectly. This is to be expected, as the LDA tends to overestimate bond strengths (and thus phonon frequencies). In addition there are small differences in particular mode intensities, for example around 27 meV in Fig. 13(b) where the intensity is underestimated or at 23 meV in Fig. 13(c) where it is overestimated. Given the lack of adjustable parameters however the agreement is remarkably strong. Notwithstanding disagreements that arise from the ab initio code not fully capturing the physics of the material, Fig. 13(a-d) show the importance of accounting for instrumental resolution when comparing calculation with experimental data.